R PROGRAMMING QUESTION

- Below I have code. For each double hashtag (##) can you comment on the code below it (Describe what is happening in the code below it next to each ##)

- Run the code and compare the confidence intervals

- I have to submit the confidence intervals, comments, and the code with the ## filled out

Leaps.then.press.plot.2<-function(xmat0,yvec,xpred,ncheck=20)

{

#

#input quadratic matrix with less than 30 columns eg. the result of x.auto2a<-matrix.2ndorder.make(xmat[,-7],F)

#also, no need for plotting, just pull out best, xpred is one of the row vectors from x.auto2a, but all terms with weight are divided by 2

#

##

leaps.str<-leaps(xmat,yvec)

##

z1<-leaps.str$Cp-leaps.str$size

##

o1<-order(z1)

matwhich<-(leaps.str$which[o1,])[1:ncheck,]

MPSEvec<-NULL

##

for(i in 1:ncheck){

ls.str0<-regpluspress(xmat[,matwhich[i,]],yvec)

##

parvec<-matwhich[i,]

npar<-sum(parvec)

## (WHY npar+1)

MPSE<-ls.str0$press/(length(yvec)-(npar+1))

MPSEvec<-c(MPSEvec,MPSE)

}

##

I1<-(MPSEvec==min(MPSEvec))

##

i<-c(1:ncheck)[I1]

##

xmat.out<-xmat[,matwhich[i,]]

##

xpred.out<-xpred[matwhich[i,]]

##

list(xmatout=xmat.out,yvec=yvec,xpredout=xpred.out)

}

Bootreg<-function(xmat,yvec,xpred,nboot=10000,alpha=0.05)

{

##

lstr0<-leaps.then.press.plot2(xmat,yvec,xpred)

xmat0<-lstr0$xmat.out

yvec0<-lstr0$yvec

xpred0<-lsstr0$xpredout

##

rprd.list<-regpred(xpred0,xmat0,yvec0)

ypred0<-rprd.list$pred

sdpred0<-rprd.list$sd

df<-rprd.list$df

##

bootvec<-NULL

nobs<-length(yvec0)

for(i in 1:nboot){

##

vboot<-sample(c(1:nobs),replace=T)

xmatb<-xmat0[vboot,]

yvecb<-yvec0[vboot]

##

lstrb<-leaps.then.press.plot2(xmatb,yvecb,xpred)

##

xmatb0<-lstrb$xmat.out

yvecb0<-lstrb$yvec

xpredb0<-lsstrb$xpredout

##

rprd.list<-regpred(xpred0,xmat0,yvec0)

ypredb<-rprd.list$pred

sdpredb<-rprd.list$sd

dfb<-rprd.list$df

##

bootvec<-c(bootvec,(ypredb-ypred0)/sdpredb)

}

##

lq<-quantile(bootvec,alpha/2)

uq<-quantile(bootvec,1-alpha/2)

##

LB<-ypred0-(sdpred0)*uq

UB<-ypred0-(sdpred0)*lq

##

NLB<-ypred0-(sdpred0)*qt(1-alpha/2,df0)

NUB<-ypred0+(sdpred0)*qt(1-alpha/2,df0)

list(bootstrap.confidence.interval=c(LB,UB),normal.confidence.interval=c(NLB,NUB))

}

> regpred<-

function(xpred,xmat,y){

##

ls.str<-lsfit(xmat,y)

#calculate prediction

ypred<-ls.str$coef%*%c(1,xpred)

#use ls.diag to extract covariance matrix

ycov<-ls.diag(ls.str)$cov.unscaled

#use ls.diag to extract std deviation

std.dev<-ls.diag(ls.str)$std.dev

#variance of data around line

v1<-std.dev^2

#variance of prediction

vpred<-v1*c(1,xpred)%*%ycov%*%c(1,xpred)

df=length(y)-length(diag(ycov))

list(pred=ypred,sd=sqrt(vpred),df=df)

}