rm(list=ls())

library(glmnet)

d = read.csv("riboflavin.csv")
d = t(d)
names = d[1,]
d = d[-1,]
row.names(d) = NULL
d = as.numeric(d)
d = matrix(d,ncol=4089)
Y = d[,1]
X = d[,-1]
fit = glmnet(x=X,y=Y,alpha=0)
beta = fit$beta



pdf("RiboRidge1.pdf")

# Using glmnet function for ridge regression
fit <- glmnet(x=X,y=Y, alpha = 0)

beta = fit$beta
lambda = fit$lambda
matplot(log(fit$lambda),t(beta),type="l",lty=1,xlab="",ylab="",lwd=2)
title(ylab=expression(beta),line=2.2, cex.lab=1.8)
title(xlab=expression(log(lambda)),line=2.8, cex.lab=1.8)

dev.off()

pdf("RiboRidge2.pdf")

###find 20 genes that have non-zero betas and plot them
d = ncol(X)
k = length(lambda)
enter = rep(0,d)
for(i in 1:d){
     b = abs(beta[i,])
     if(sum(b) > 0)enter[i] = min((1:k)[b != 0])
     }
I = (1:d)[enter > 0]
I = I[1:20]
matplot(log(fit$lambda),t(beta[I,]),type="l",lty=1,xlab="",ylab="",lwd=2)
title(ylab=expression(beta),line=2.2, cex.lab=1.8))
title(xlab=expression(log(lambda),line=2.7, cex.lab=1.8)

dev.off()


pdf("RiboRidge3.pdf")

fit.cv = cv.glmnet(X,Y,alpha=0)
##cross validation
lambda = fit.cv$lambda
L = log(lambda)
cv = fit.cv$cvm
se = fit.cv$cvsd
low = cv-se
up = cv+se

plot(L,cv,type="l",col="purple",lwd=3,ylim=c(min(low),max(up)),xlab="",ylab="",xaxt="n")
polygon(c(L,rev(L)),c(low,rev(up)),col="plum3",border=NA)
lines(L,cv,col="purple",lwd=3)
lines(L,cv+se,col="purple",lty="dashed",lwd=3)
lines(L,cv-se,col="purple",lty="dashed",lwd=3)
I = which.min(cv)
abline(v=L[I],col="blue",lty=4,lwd=3)
#segments(L[I],0,L[I],0.35,col="navy",lty=2,lwd=4)
title(ylab="Prediction Error",line=1.8, cex.lab=1.8)
title(xlab=expression(log(lambda)),line=2.6, cex.lab=1.8)
axis(side = 1, at = L[I],cex.lab=3,col.lab="blue")
dev.off()