################################################### ### chunk number 1: pcdemo1 ################################################### #line 32 "Lect15c.Rnw" library(MASS) x<-mvrnorm(250,mu=c(0,0),Sigma=matrix(c(1.5,1,1,1),2,2)) plot(x) prx<-prcomp(x) print(prx$sdev^2/sum(prx$sdev^2)) abline(0,prx$rot[2,1]/prx$rot[1,1],col=2) abline(0,prx$rot[2,2]/prx$rot[1,2],col=3) ### can also use ### svd(t(x)%*%x) or svd(x) or eigen(t(x)%*%x) ################################################### ### chunk number 2: pcdemo2 ################################################### #line 75 "Lect15c.Rnw" SA<-read.table('SA.dat',header=T) SA$famhist<-SA$famhist-1 ii<-sample(seq(1,dim(SA)[1]),200) mm<-lm(ldl~age+sbp+adiposity+obesity+typea+alcohol+alcind+tobacco+as.factor(tobind)+as.factor(chd)+as.factor(famhist),data=SA,subset=ii) print(summary(mm)) selmm<-step(mm,trace=F) print(summary(selmm)) ################################################### ### chunk number 3: pcdemo3 ################################################### #line 85 "Lect15c.Rnw" pmm<-predict(mm,newdata=SA[-ii,-12]) pselmm<-predict(selmm,newdata=SA[-ii,-12]) print(sum((SA$ldl[-ii]-pmm)^2)/length(pmm)) print(sum((SA$ldl[-ii]-pselmm)^2)/length(pselmm)) ################################################### ### chunk number 4: pcdemo4 ################################################### #line 94 "Lect15c.Rnw" library(elasticnet) #standardize SA2<-SA standardize<-function(x) {x<-(x-mean(x))/sd(x)} SA2<-apply(SA2,2,standardize) SA2<-as.data.frame(SA2) names(SA2)<-names(SA) ################################################### ### chunk number 5: pcdemo5 ################################################### #line 105 "Lect15c.Rnw" mm<-lm(ldl~age+sbp+adiposity+obesity+typea+alcohol+alcind+tobacco+as.factor(tobind)+as.factor(chd)+as.factor(famhist),data=SA2,subset=ii) print(summary(mm)) ################################################### ### chunk number 6: pcdemo6 ################################################### #line 111 "Lect15c.Rnw" newx<-spca(SA2[,-12],K=3,para=c(0,0,0),type=c("predictor")) print(newx) SA3<-as.matrix(SA2[,-12])%*%newx$load SA3<-data.frame(cbind(SA3,SA2[,12])) names(SA3)<-c("pc1","pc2","pc3","ldl") mm<-lm(ldl~pc1+pc2+pc3,data=SA3,subset=ii) print(summary(mm)) selmm<-step(mm,trace=F) print(summary(selmm)) ################################################### ### chunk number 7: pcdemo7 ################################################### #line 117 "Lect15c.Rnw" newx<-spca(SA2[,-12],K=11,para=rep(0,11),type=c("predictor")) SA3<-SA2 SA3[,-12]<-as.matrix(SA2[,-12])%*%newx$load names(SA3)<-c("pc1","pc2","pc3","pc4","pc5","pc6","pc7","pc8","pc9","pc10","pc11","ldl") mm<-lm(ldl~pc1+pc2+pc3+pc4+pc5+pc6+pc7+pc8+pc9+pc10+pc11,data=SA3,subset=ii) print(summary(mm)) selmm<-step(mm,trace=F) print(summary(selmm)) ################################################### ### chunk number 8: pcdemo8 ################################################### #line 129 "Lect15c.Rnw" newx<-spca(SA2[,-12],K=11,para=rep(3,11),type=c("predictor")) print(newx) ################################################### ### chunk number 9: pcdemo9 ################################################### #line 135 "Lect15c.Rnw" SA3<-SA2 SA3[,-12]<-as.matrix(SA2[,-12])%*%newx$load names(SA3)<-c("pc1","pc2","pc3","pc4","pc5","pc6","pc7","pc8","pc9","pc10","pc11","ldl") mm<-lm(ldl~pc1+pc2+pc3+pc4+pc5+pc6+pc7+pc8+pc9+pc10+pc11,data=SA3,subset=ii) print(summary(mm)) selmm<-step(mm,trace=F) print(summary(selmm)) ################################################### ### chunk number 10: ridgedemo ################################################### #line 169 "Lect15c.Rnw" xte<-mvrnorm(2500,mu=c(0,0),Sigma=rbind(c(1,.95),c(.95,1))) lambda<-seq(0,1,by=.001) coefmat<-matrix(0,length(lambda),3) x<-mvrnorm(25,mu=c(0,0),Sigma=rbind(c(1,.95),c(.95,1))) y<-2+3*x[,1]+2*x[,2]+rnorm(25) mm<-lm.ridge(y~x,lambda=seq(0,1,by=.001)) coefmat<-coef(mm) ppm<-cbind(rep(1,dim(xte)[1]),xte)%*%t(coefmat) ynew<-2+3*xte[,1]+2*xte[,2]+rnorm(2500) pmse<-(apply((ynew-ppm)^2/dim(ppm)[1],2,sum)) plot(seq(0,1,by=.001),pmse,xlab="r",ylab="pMSE",type='l') abline(h=pmse[1],lty=2) ################################################### ### chunk number 11: penreg ################################################### #line 229 "Lect15c.Rnw" library(ellipse) plot(ellipse(.5*rbind(c(2,-.9),c(-.9,1)),centre=c(2,2)),type='l',xlim=c(-2,6),ylim=c(-2,6),xlab="beta1",ylab="beta2") points(2,2,pch=3) lines(ellipse(.25*rbind(c(1,0),c(0,1))),col=2,lwd=2) abline(h=0) abline(v=0) lines(ellipse(.4*rbind(c(2,-.9),c(-.9,1)),centre=c(2,2))) lines(ellipse(.2*rbind(c(2,-.9),c(-.9,1)),centre=c(2,2))) lines(ellipse(.1*rbind(c(2,-.9),c(-.9,1)),centre=c(2,2))) lines(ellipse(.75*rbind(c(2,-.9),c(-.9,1)),centre=c(2,2))) lines(ellipse(.5*rbind(c(1,0),c(0,1))),col=2,lwd=2) lines(ellipse(.75*rbind(c(1,0),c(0,1))),col=2,lwd=2) ################################################### ### chunk number 12: penreg2 ################################################### #line 255 "Lect15c.Rnw" plot(ellipse(.5*rbind(c(2,-.9),c(-.9,1)),centre=c(2,2)),type='l',xlim=c(-2,6),ylim=c(-2,6),xlab="beta1",ylab="beta2") points(2,2,pch=3) abline(h=0) abline(v=0) lines(ellipse(.4*rbind(c(2,-.9),c(-.9,1)),centre=c(2,2))) lines(ellipse(.2*rbind(c(2,-.9),c(-.9,1)),centre=c(2,2))) lines(ellipse(.1*rbind(c(2,-.9),c(-.9,1)),centre=c(2,2))) lines(ellipse(.75*rbind(c(2,-.9),c(-.9,1)),centre=c(2,2))) lines(c(0,1),c(1,0),col=2,lwd=2) lines(c(0,1),c(-1,0),col=2,lwd=2) lines(c(-1,0),c(0,-1),col=2,lwd=2) lines(c(-1,0),c(0,1),col=2,lwd=2) lines(c(0,2),c(2,0),col=2,lwd=2) lines(c(0,2),c(-2,0),col=2,lwd=2) lines(c(-2,0),c(0,-2),col=2,lwd=2) lines(c(-2,0),c(0,2),col=2,lwd=2) lines(c(0,3),c(3,0),col=2,lwd=2) lines(c(0,3),c(-3,0),col=2,lwd=2) lines(c(-3,0),c(0,-3),col=2,lwd=2) lines(c(-3,0),c(0,3),col=2,lwd=2) ################################################### ### chunk number 13: penreg3 ################################################### #line 293 "Lect15c.Rnw" library(ellipse) plot(ellipse(.5*rbind(c(2,-.9),c(-.9,1)),centre=c(2,2)),type='l',xlim=c(-2,6),ylim=c(-2,6),xlab="beta1",ylab="beta2") points(2,2,pch=3) abline(h=0) abline(v=0) lines(ellipse(.4*rbind(c(2,-.9),c(-.9,1)),centre=c(2,2))) lines(ellipse(.2*rbind(c(2,-.9),c(-.9,1)),centre=c(2,2))) lines(ellipse(.1*rbind(c(2,-.9),c(-.9,1)),centre=c(2,2))) lines(ellipse(.75*rbind(c(2,-.9),c(-.9,1)),centre=c(2,2))) lines(c(-1,1),c(-1,-1),col=2,lwd=2) lines(c(1,1),c(-1,1),col=2,lwd=2) lines(c(-1,1),c(1,1),col=2,lwd=2) lines(c(-1,-1),c(-1,1),col=2,lwd=2) lines(c(-2,2),c(-2,-2),col=2,lwd=2) lines(c(2,2),c(-2,2),col=2,lwd=2) lines(c(-2,2),c(2,2),col=2,lwd=2) lines(c(-2,-2),c(-2,2),col=2,lwd=2) lines(c(-1.5,1.5),c(-1.5,-1.5),col=2,lwd=2) lines(c(1.5,1.5),c(-1.5,1.5),col=2,lwd=2) lines(c(-1.5,1.5),c(1.5,1.5),col=2,lwd=2) lines(c(-1.5,-1.5),c(-1.5,1.5),col=2,lwd=2) ################################################### ### chunk number 14: betafcn ################################################### #line 328 "Lect15c.Rnw" plot(c(-1,1),c(-1,1),lwd=2,xlab="beta",ylab="reg-beta",type='l') lines(c(-1,1),.8*c(-1,1),col=2,lwd=2) bb<-seq(-1,1,by=.01) bb2<-bb bb2[abs(bb)<.25]<-0 bb2[bb>=.25]<-bb[bb>=.25]-.25 bb2[bb<=(-.25)]<-bb[bb<=(-.25)]+.25 lines(seq(-1,1,by=.01),bb2,lwd=2,col=3) abline(h=0) abline(v=0) ################################################### ### chunk number 15: demo15a ################################################### #line 353 "Lect15c.Rnw" library(lars) xx<-as.matrix(cbind(rep(1,dim(SA)[1]),as.matrix(SA[,-12]))) colnames(xx)<-c("int",names(SA)[-12]) yy<-SA[,12] ll<-lars(xx[ii,],yy[ii],intercept=F) print(summary(ll)) print(ll) ################################################### ### chunk number 16: demo15b ################################################### #line 365 "Lect15c.Rnw" plot(ll) ################################################### ### chunk number 17: demo15c ################################################### #line 380 "Lect15c.Rnw" selll<-ll$beta[which.min(ll$Cp),] pll<-xx[-ii,]%*%selll # print(sum(SA$ldl[-ii]-xx[-ii,]%*%ll$beta[dim(ll$beta)[1],])^2/length(pll)) print(sum((SA$ldl[-ii]-pll)^2)/length(pll)) ################################################### ### chunk number 18: demo15d ################################################### #line 392 "Lect15c.Rnw" library(elasticnet) xx<-as.matrix(cbind(rep(1,dim(SA)[1]),as.matrix(SA[,-12]))) colnames(xx)<-c("int",names(SA)[-12]) yy<-SA[,12] ll<-enet(xx[ii,],yy[ii],intercept=F) print(summary(ll)) print(ll) ################################################### ### chunk number 19: demo15e ################################################### #line 401 "Lect15c.Rnw" plot(ll) ################################################### ### chunk number 20: demo15f ################################################### #line 416 "Lect15c.Rnw" selll<-ll$beta[which.min(ll$Cp),] pll<-xx[-ii,]%*%selll # print(sum(SA$ldl[-ii]-xx[-ii,]%*%ll$beta[dim(ll$beta)[1],])^2/length(pll)) print(sum((SA$ldl[-ii]-pll)^2)/length(pll)) ################################################### ### chunk number 21: demo15g ################################################### #line 426 "Lect15c.Rnw" ll<-enet(xx[ii,],yy[ii],intercept=F,lambda=.4) print(summary(ll)) print(ll) ################################################### ### chunk number 22: demo15h ################################################### #line 431 "Lect15c.Rnw" plot(ll) ################################################### ### chunk number 23: demo15i ################################################### #line 445 "Lect15c.Rnw" selll<-ll$beta[which.min(ll$Cp),] pll<-xx[-ii,]%*%selll # print(sum(SA$ldl[-ii]-xx[-ii,]%*%ll$beta[dim(ll$beta)[1],])^2/length(pll)) print(sum((SA$ldl[-ii]-pll)^2)/length(pll))