###################################################
### chunk number 1: l6-fig4plot
###################################################
#line 44 "Lecture7.Rnw"
n<-50; noisesd<-2
x<-sort(rnorm(n))
K<-50
mma<-rep(0,2) ## creating a vector to store the model coefficients in
## The for-loop runs everything inside the curly brackets K times
for (kk in (1:K)) {
 y<-2*x-.5*x^2+rnorm(n,sd=noisesd)
 if (kk==1 ) {plot(x, y)} ## plot the data once to get the axes fixed
 mm<-lm(y~x)
 lines(x,mm$fit,col=3,lty=2) ## plotting the different fites in col=3 (green).
 mma<-mma+(1/K)*mm$coef } ## storing the coefficients
lines(x,2*x-.5*x^2,lwd=2) ## true model
lines(x,mma[1]+mma[2]*x,lwd=2,col=2) ## average model


###################################################
### chunk number 2: l6-fig5plot
###################################################
#line 59 "Lecture7.Rnw"
mma<-rep(0,3)
for (kk in (1:K)) {
 y2<-2*x-.5*x^2+rnorm(n,sd=noisesd)
 if (kk==1) {plot(x, y2,lwd=2)}
 mm<-lm(y2~x+I(x^2))
 lines(x,mm$fit,col=3,lty=2)
 mma<-mma+(1/K)*mm$coef}
lines(x,2*x-.5*x^2,lwd=2)
lines(x,mma[1]+mma[2]*x+mma[3]*x^2,lwd=2,col=2)


###################################################
### chunk number 3: l6-fig6plot
###################################################
#line 70 "Lecture7.Rnw"
mma<-rep(0,length(x))
for (kk in (1:K)) {
 y3<-2*x-.5*x^2+rnorm(n,sd=noisesd)
 if (kk==1) {plot(x, y3) }
 mm3<-loess(y3~x,degree=0,span=.1)
 lines(x,mm3$fit,col=3,lty=2)
 mma<-mma+(1/K)*mm3$fit }
lines(x,2*x-.5*x^2,lwd=2)
lines(x,mma,lwd=2,col=2)


###################################################
### chunk number 4: fig6plot
###################################################
#line 81 "Lecture7.Rnw"
mma<-rep(0,length(x))
for (kk in (1:K)) {
 y3<-2*x-.5*x^2+rnorm(n,sd=noisesd)
 if (kk==1) {plot(x, y3) }
 mm3<-loess(y3~x,degree=0,span=.1)
 lines(x,mm3$fit,col=3,lty=2)
 mma<-mma+(1/K)*mm3$fit }
lines(x,2*x-.5*x^2,lwd=2)
lines(x,mma,lwd=2,col=2)


###################################################
### chunk number 5: ldl1
###################################################
#line 129 "Lecture7.Rnw"
SA<-data.frame(read.table('SA.dat',sep='\t',header=T)) ## read in the data


###################################################
### chunk number 6: l6-rssfig
###################################################
#line 132 "Lecture7.Rnw"
library(leaps)
frac<-.5 ## creating a training data set of proportion frac
ii<-sample(seq(1,dim(SA)[1]),round(dim(SA)[1]*frac)) ## data points to use for training
yy<-SA[ii,12] ## training data
xx<-as.matrix(SA[ii,-c(12)])
yyt<-SA[-ii,12] ## test data
xxt<-as.matrix(SA[-ii,-c(12)])
###
rleaps<-regsubsets(xx,yy,int=T,nbest=250,nvmax=dim(SA)[2],really.big=T,method=c("ex")) ## all subset models
cleaps<-summary(rleaps,matrix=T) ## True/False matrix. The r-th is a True/False statement about which
## variables are included in model r.
tt<-apply(cleaps$which,1,sum) ## size of each model in the matrix
mses<-cleaps$rss/length(yy) ## corresponding MSEs
###
plot(tt,mses,xlab="number of variables",ylab="MSE",main="MSE for all subset models")
tmin<-min(tt)
tmax<-max(tt)
tsec<-seq(tmin,tmax)
msevec<-rep(0,length(tsec))
for (tk in 1:length(tsec)) {
 msevec[tk]<-min(mses[tt==tsec[tk]])} ## the best model for each size
lines(tsec,msevec,lwd=2,col=2) ## a line connecting the best models.


###################################################
### chunk number 7: l6-msefig
###################################################
#line 156 "Lecture7.Rnw"
plot(tsec,msevec,xlab="number of variables",ylab="MSE",main="MSE for best model of each size",type='b',col=4,lwd=2)


###################################################
### chunk number 8: l6-pmsefig
###################################################
#line 193 "Lecture7.Rnw"
pmses<-mses
for (ta in (2:dim(cleaps$which)[1])) {
 mmr<-lm(yy~xx[,cleaps$which[ta,-1]==T])
 PEcp<-sum((yyt-cbind(rep(1,dim(xxt)[1]),xxt[,cleaps$which[ta,-1]==T])%*%mmr$coef)^2)/length(yyt)
 pmses[ta]<-PEcp }
pmses[1]<-max(pmses)
pmsevec<-rep(0,length(tsec))
for (tk in 1:length(tsec)) {
 pmsevec[tk]<-min(pmses[tt==tsec[tk]])}
plot(tsec,pmsevec,xlab="number of variables",ylab="pMSE",main="prediction MSE", type='b',lwd=2,col=2)


###################################################
### chunk number 9: l6-whichwin
###################################################
#line 219 "Lecture7.Rnw"
ptmin<-which.min(pmses)
pmod<-cleaps$which[ptmin,]
winsize<-sum(pmod)
mtmin<-which.min(mses[tt==sum(pmod)])
mod<-(cleaps$which[tt==sum(pmod),])[mtmin,]


###################################################
### chunk number 10: l6-wintable
###################################################
#line 231 "Lecture7.Rnw"
library(xtable)
x<-rbind(as.character(pmod),as.character(mod))
dimnames(x) <- list(c("pMSE win", "MSE of same size"), c("Intercept",names(SA)[-12]))
x1<-x[,1:6]
x2<-x[,7:12]
dimnames(x1)<-list(c("pMSE win", "MSE of same size"),c("Intercept",names(SA)[1:5]))
dimnames(x2)<-list(c("pMSE win", "MSE of same size"),c(names(SA))[6:11])
xtable(rbind(x1), digits=c(0, rep(3,6)))
xtable(rbind(x2), digits=c(0, rep(3,6)),caption="Selected variables using pMSE",label="table:pmse")


###################################################
### chunk number 11: cars1
###################################################
#line 246 "Lecture7.Rnw"
cars<-data.frame(read.table('cars.dat',sep='\t',header=T)) ## read in the data
print(dim(cars)) ## dimensions
print(names(cars)) ## variable names


###################################################
### chunk number 12: cars1b
###################################################
#line 256 "Lecture7.Rnw"
ntrain<-60
ii<-sample(seq(1,dim(cars)[1]),ntrain)
# ntrain cars = training set
cars.train<-cars[ii,]
row.names(cars.train)<-seq(1,ntrain)
cars.test<-cars[-ii,]
row.names(cars.test)<-seq(1,dim(cars.test)[1])


###################################################
### chunk number 13: cars2a
###################################################
#line 267 "Lecture7.Rnw"
par(mfrow=c(2,2))
plot(cars.train$mid,cars.train$ci,main="citympg on price")
plot(cars.train$hw,cars.train$ci,main="citympg on hwmpg")
plot(cars.train$le,cars.train$ci,main="citympg on length")
plot(cars.train$engsize,cars.train$ci,main="citympg on enginesize")
p<-locator()


###################################################
### chunk number 14: cars2b
###################################################
#line 275 "Lecture7.Rnw"
par(mfrow=c(2,2))
plot(cars.train$le,cars.train$wi,main="width on length")
plot(cars.train$wi,cars.train$we,main="weight on width")
plot(cars.train$le,cars.train$wh,main="wheelbase on length")
plot(cars.train$lu,cars.train$rea,main="rearroom on luggage room")
p<-locator()


###################################################
### chunk number 15: cars3a
###################################################
#line 308 "Lecture7.Rnw"
par(mfrow=c(2,2))
plot(cars.train$ci,cars.train$mid,main="price on citymg")
plot(cars.train$we,cars.train$mid,main="price on weight")
plot(cars.train$ho,cars.train$mid,main="price on horsepower")
plot(cars.train$man,cars.train$mid,main="price on transmission")
p<-locator()


###################################################
### chunk number 16: cars3b
###################################################
#line 327 "Lecture7.Rnw"
par(mfrow=c(2,2))
plot(1/cars.train$ci,log(cars.train$mid),main="log(price) on 1/citympg")
plot(cars.train$we,log(cars.train$mid),main="log(price) on weight")
plot(cars.train$ho,log(cars.train$mid),main="log(price) on horsepower")
plot(cars.train$man,log(cars.train$mid),main="log(price) on transmission")
p<-locator()


###################################################
### chunk number 17: cars4
###################################################
#line 349 "Lecture7.Rnw"
cars.train[,1]<-log(cars.train$mid)
cars.train[,2]<-1/cars.train$ci
cars.train[,3]<-1/cars.train$hw
cars.test[,1]<-log(cars.test$mid)
cars.test[,2]<-1/cars.test$ci
cars.test[,3]<-1/cars.test$hw
names(cars.train)[c(1,2,3)]<-c("log.mid.price","city.gpm","hw.gpm")
names(cars.test)[c(1,2,3)]<-c("log.mid.price","city.gpm","hw.gpm")


###################################################
### chunk number 18: cars5
###################################################
#line 361 "Lecture7.Rnw"
mm<-lm(log.mid.price~city.gpm+hw.gpm+airbagstd+cylnbr+engsize+horsepwr+rpm.at.max+engrev.high+mantrans.op+fueltank+passengers+length+width+uturn+rearroom+luggage+weight+domestic,data=cars.train)
print(ms<-summary(mm))


###################################################
### chunk number 19: cars5b
###################################################
#line 365 "Lecture7.Rnw"
par(mfrow=c(2,2))
plot(mm)
p<-locator()


###################################################
### chunk number 20: cars6
###################################################
#line 385 "Lecture7.Rnw"
par(mfrow=c(1,1))
plot(mm$fit,mm$res,xlab="fitted values",ylab="residuals")
abline(h=0)
id1<-identify(mm$fit,mm$res,pos=T)


###################################################
### chunk number 21: cars6b
###################################################
#line 391 "Lecture7.Rnw"
plot(mm$fit,mm$res,xlab="fitted values",ylab="residuals")
abline(h=0)
if (max(id1$ind)!=-Inf) {text(mm$fit[id1$ind],mm$res[id1$ind],id1$ind,pos=id1$pos)}
p<-locator()


###################################################
### chunk number 22: cars7
###################################################
#line 408 "Lecture7.Rnw"
qq<-seq(.5/ntrain,(ntrain-.5)/ntrain,length=ntrain)
normq<-qnorm(p=qq)
rsort<-sort(rstandard(mm))
rlist<-sort.list(rstandard(mm))
plot(normq,rsort,xlab="Theoretical quantiles",ylab="Standardized residuals")
qr<-quantile(rstandard(mm))
qn<-quantile(qnorm(p=qq))
b<-(qr[4]-qr[2])/(qn[4]-qn[2])
a<-qr[4]-b*qn[4]
abline(a,b)
id2<-identify(normq,sort(rstandard(mm)),label=rlist,pos=T)


###################################################
### chunk number 23: cars7b
###################################################
#line 421 "Lecture7.Rnw"
plot(normq,rsort,xlab="Theoretical quantiles",ylab="Standardized residuals")
qr<-quantile(rstandard(mm))
qn<-quantile(qnorm(p=qq))
b<-(qr[4]-qr[2])/(qn[4]-qn[2])
a<-qr[4]-b*qn[4]
abline(a,b)
if (max(id2$ind)!=-Inf) {text(normq[id2$ind],rsort[id2$ind],rlist[id2$ind],pos=id2$pos)}
p<-locator()


###################################################
### chunk number 24: cars8
###################################################
#line 442 "Lecture7.Rnw"
plot(mm$fit,abs(rstandard(mm)),xlab="fitted values", ylab="|standardized residuals|")
id3<-identify(mm$fit,abs(rstandard(mm)),pos=T)


###################################################
### chunk number 25: cars8b
###################################################
#line 446 "Lecture7.Rnw"
plot(mm$fit,abs(rstandard(mm)),xlab="fitted values", ylab="|standardized residuals|")
if (max(id3$ind)!=-Inf) {text(mm$fit[id3$ind],abs(rstandard(mm))[id3$ind],id3$ind,pos=id3$pos)}
p<-locator()


###################################################
### chunk number 26: cars9
###################################################
#line 461 "Lecture7.Rnw"
lm1<-lm.influence(mm)
cooksd<-cooks.distance(mm)
plot(cooksd,main="Cook's Distance",type="h")
abline(h=qf(.95,1,mm$df),lty=2)
idc<-identify(cooksd,pos=T)


###################################################
### chunk number 27: cars9b
###################################################
#line 468 "Lecture7.Rnw"
plot(cooksd,main="Cook's Distance",type="h")
abline(h=qf(.95,1,mm$df),lty=2)
if (max(idc$ind)!=-Inf) {text(idc$ind,cooksd[idc$ind],idc$ind,pos=idc$pos) }


###################################################
### chunk number 28: cars10
###################################################
#line 483 "Lecture7.Rnw"
plot(lm1$hat,main="Leverage")
idlev<-identify(lm1$hat,pos=T)


###################################################
### chunk number 29: cars10b
###################################################
#line 487 "Lecture7.Rnw"
plot(lm1$hat,main="Leverage")
abline(h=3*length(mm$coef)/dim(cars.train)[1])
if (max(idlev$ind)!=-Inf) {text(idlev$ind,lm1$hat[idlev$ind],idlev$ind,pos=idlev$pos)}


###################################################
### chunk number 30: cars11
###################################################
#line 494 "Lecture7.Rnw"
plot(lm1$coeff[,4],main="change in slope 4")
id4<-identify(lm1$coeff[,4],pos=T)


###################################################
### chunk number 31: cars11b
###################################################
#line 498 "Lecture7.Rnw"
plot(lm1$coeff[,4],main="change in slope 4")
if (max(id4$ind)!=-Inf) {text(id4$ind,lm1$coeff[id4$ind,4],id4$ind,pos=id4$pos)}


###################################################
### chunk number 32: cars12
###################################################
#line 502 "Lecture7.Rnw"
plot(lm1$coeff[,7],main="change in slope 7")
id7<-identify(lm1$coeff[,7],pos=T)


###################################################
### chunk number 33: cars12b
###################################################
#line 506 "Lecture7.Rnw"
plot(lm1$coeff[,7],main="change in slope 7")
if (max(id7$ind)!=-Inf) {text(id7$ind,lm1$coeff[id7$ind,7],id7$ind,pos=id7$pos)}


###################################################
### chunk number 34: cars13
###################################################
#line 510 "Lecture7.Rnw"
plot(lm1$sig,main="change in sigma")
ids<-identify(lm1$sig,pos=T)


###################################################
### chunk number 35: cars13b
###################################################
#line 514 "Lecture7.Rnw"
plot(lm1$sig,main="change in sigma")
if (max(ids$ind)!=-Inf) {text(lm1$sig[ids$ind],ids$ind,pos=ids$pos)}


###################################################
### chunk number 36: cars14
###################################################
#line 533 "Lecture7.Rnw"
indvec<-sort(c(id1$ind,rlist[id2$ind],id3$ind,idlev$ind,id4$ind,id7$ind,ids$ind))
print(table(indvec)) ## how many of each?
maxid<-max(table(indvec))
indout<-unique(indvec)[table(indvec)==max(table(indvec))]


###################################################
### chunk number 37: cars15
###################################################
#line 546 "Lecture7.Rnw"
## regression without the identified outlier
mmb<-lm(log.mid.price~city.gpm+hw.gpm+airbagstd+cylnbr+engsize+horsepwr+rpm.at.max+engrev.high+mantrans.op+fueltank+passengers+length+width+uturn+rearroom+luggage+weight+domestic,data=cars.train,subset=-indout)
print(summary(mmb))


###################################################
### chunk number 38: cars15b
###################################################
#line 551 "Lecture7.Rnw"
par(mfrow=c(2,2))
plot(mmb)
p<-locator()


###################################################
### chunk number 39: cars16
###################################################
#line 570 "Lecture7.Rnw"
selectstep<-step(mmb,trace=F)


###################################################
### chunk number 40: cars17
###################################################
#line 575 "Lecture7.Rnw"
print(summary(selectstep))


###################################################
### chunk number 41: cars18
###################################################
#line 579 "Lecture7.Rnw"
predval<-predict(selectstep,newdata=cars.test)
prederror<-sum((cars.test[,1]-predval)^2)


###################################################
### chunk number 42: cars19
###################################################
#line 586 "Lecture7.Rnw"
# compare refit
mmtest<-lm(log.mid.price~city.gpm+hw.gpm+airbagstd+cylnbr+engsize+horsepwr+rpm.at.max+engrev.high+mantrans.op+fueltank+passengers+length+width+uturn+rearroom+luggage+weight+domestic,data=cars.test)
print(selecttest<-step(mmtest,trace=F))
fiterror<-sum(summary(selecttest)$res^2)
p<-locator()


###################################################
### chunk number 43: cars20
###################################################
#line 599 "Lecture7.Rnw"
yy<-cars.train[,1] ## training data
xx<-as.matrix(cars.train[,-1])
yyt<-cars.test[,1] ## test data
xxt<-as.matrix(cars.test[,-1])
###
rleaps<-regsubsets(xx,yy,int=T,nbest=250,nvmax=dim(cars)[2],really.big=T,method=c("ex")) ## all subset models
cleaps<-summary(rleaps,matrix=T) ## True/False matrix. The r-th is a True/False statement about which
## variables are included in model r.
tt<-apply(cleaps$which,1,sum) ## size of each model in the matrix
mses<-cleaps$rss/length(yy) ## corresponding MSEs


###################################################
### chunk number 44: cars20b
###################################################
#line 613 "Lecture7.Rnw"
plot(tt,mses,xlab="number of variables",ylab="MSE",main="MSE for all subset models")
tmin<-min(tt)
tmax<-max(tt)
tsec<-seq(tmin,tmax)
msevec<-rep(0,length(tsec))
for (tk in 1:length(tsec)) {
 msevec[tk]<-min(mses[tt==tsec[tk]])} ## the best model for each size
lines(tsec,msevec,lwd=2,col=2) ## a line connecting the best models.
p<-locator()


###################################################
### chunk number 45: cars20c
###################################################
#line 624 "Lecture7.Rnw"
plot(tsec,msevec,xlab="number of variables",ylab="MSE",main="MSE for best model of each size",type='b',col=4,lwd=2)
p<-locator()


###################################################
### chunk number 46: cars21
###################################################
#line 642 "Lecture7.Rnw"
pmses<-mses
for (ta in (2:dim(cleaps$which)[1])) {
 mmr<-lm(yy~xx[,cleaps$which[ta,-1]==T])
 PEcp<-sum((yyt-cbind(rep(1,dim(xxt)[1]),xxt[,cleaps$which[ta,-1]==T])%*%mmr$coef)^2)/length(yyt)
 pmses[ta]<-PEcp }
pmses[1]<-max(pmses)
pmsevec<-rep(0,length(tsec))
for (tk in 1:length(tsec)) {
 pmsevec[tk]<-min(pmses[tt==tsec[tk]])}


###################################################
### chunk number 47: cars21b
###################################################
#line 655 "Lecture7.Rnw"
plot(tsec,pmsevec,xlab="number of variables",ylab="pMSE",main="prediction MSE", type='b',lwd=2,col=2)


###################################################
### chunk number 48: cars22
###################################################
#line 671 "Lecture7.Rnw"
ptmin<-which.min(pmses)
pmod<-cleaps$which[ptmin,]
winsize<-sum(pmod)
mtmin<-which.min(mses[tt==sum(pmod)])
mod<-(cleaps$which[tt==sum(pmod),])[mtmin,]
print(names(mod[mod==T])[-1])
print(names(selectstep$model)[-1])