## Prostate data
require(glmnet)
prostate <- read.delim("../../data/prostate.txt")
prostate <- prostate[,-ncol(prostate)]
X <- model.matrix(lpsa~0+.,prostate)
y <- prostate$lpsa
fit <- glmnet(X,y)
plot(fit)
cvfit <- cv.glmnet(X,y)
plot(cvfit)
coef(fit,s=cvfit$lambda.min)

## Slide 14
lam <- c(0,exp(seq(log(0.01),log(1e8),len=201)))
fit <- lm.ridge(lpsa~.,prostate,lambda=lam)
l <- apply(coef(fit)[,-1]^2,1,sum)
s <- l/max(l)
matplot(c(s,0),rbind(coef(fit)[,-1],0),type="l",lwd=2,lty=1,col="slateblue",xlab="",ylab="Coefficients",xlim=c(0,1.15),ylim=c(-0.1,0.8),main="Ridge")
text(1.1,coef(fit)[1,-1],labels=colnames(X),cex=0.8)
mtext(expression(sum(beta[j]^2)/max(sum(beta[j]^2))),1,line=3)

cvfit <- cv.glmnet(X,y)
fit <- glmnet(X,y,lambda.min=0,nlambda=501)
l <- apply(abs(coef(fit)[-1,]),2,sum)
s <- l/max(l)
matplot(s,t(as.matrix(coef(fit)[-1,])),type="l",lwd=2,lty=1,col="slateblue",xlab="",ylab="Coefficients",xlim=c(0,1.15),ylim=c(-0.1,0.8),main="Lasso")
text(1.1,coef(fit)[-1,length(fit$lambda)],labels=colnames(X),cex=0.8)
mtext(expression(sum(abs(beta[j]))/max(sum(abs(beta[j])))),1,line=3)

## Selection of lambda
b <- coef(fit,cvfit$lambda.min)[-1]
abline(v=sum(abs(b))/max(l),col="gray80")
RSS <- apply(predict(fit,newx=X)-y,2,crossprod)
n <- length(y)
AIC <- n*log(RSS) + 2*fit$df
BIC <- n*log(RSS) + log(n)*fit$df
GCV <- RSS/n/(1-fit$df/n)^2
abline(v=s[which.min(AIC)],col=rgb(1,.5,.5))
abline(v=s[which.min(BIC)],col=rgb(.5,1,.5))