Question

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)

}

Answer 1

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)

}