library(mfp) #for body fat example
library(leaps)

##### 
#lasso: (don't show this? keep for an exercise?)
rm(list=ls())
data(bodyfat)
bodyfat2 <- bodyfat[,-c(1,2,4)]

n <- nrow(bodyfat2)
K <- 10
#create the vector that gives index of the fold (the last fold may have more/less obs)
I <- rep(K,n) #
for (k in 1:(K-1)) I[sample((1:n)[I==K],n/K)] <- k
table(I)

library(glmnet)
lasso1 <- glmnet(y=c(bodyfat2[,1]),x=as.matrix(bodyfat2[,-1]),family='gaussian',alpha=1)
lambdaseq <- lasso1$lambda

# I use the function cv.glmnet here to check if my results below are consistent with those obtained using the package
cvglmnet1 <- cv.glmnet(y=c(bodyfat2[,1]),x=as.matrix(bodyfat2[,-1]),family='gaussian',alpha=1,lambda=lambdaseq,foldid=I)
cvglmnet1$cvm
cvglmnet1$cvsd

yhat <- matrix(ncol=length(lambdaseq),nrow=n)
for (k in 1:K){
	bodyfat2k <- bodyfat2[!(I==k),]
	fitk <- glmnet(y=c(bodyfat2k[,1]),x=as.matrix(bodyfat2k[,-1]),family='gaussian',alpha=1,lambda=lambdaseq)
	yhat[(I==k),] <- predict.glmnet(fitk,newx=as.matrix(bodyfat2[(I==k),-1]))
}
Errip <- (yhat-matrix(bodyfat2$siri,ncol=length(lambdaseq),nrow=n))^2
cv <- apply(Errip,2,mean)
se <- apply(Errip,2, function(x) sd(unlist(lapply(split(x,I),mean))) )/sqrt(K) #I need to block the data again
cv
se #my SE are a little different but this seems ok
plot(lambdaseq,cv,ylim=range(c(cv+se,cv-se)),log='x')
segments(x0=lambdaseq,y0=cv-se,y1=cv+se,col='grey')
points(lambdaseq[which.min(cv)],cv[which.min(cv)],pch=20)
abline(h=(cv+se)[which.min(cv)],lty='dotted')
max(lambdaseq[which(cv<((cv+se)[which.min(cv)]))]) #my best lambda using one-SE rule