require(stepPlr)   ###Implements logistic regression with quadratic penalization
require(glmnet)    ###Elastic net penalty for glmmodels 
require(rda)       ###Regularized discriminant analysis
require(ggplot2)

SRBCTxtrain <- read.table("http://www-stat-class.stanford.edu/~tibs/ElemStatLearn/datasets/khan.xtrain")
SRBCTytrain <- read.table("http://www-stat-class.stanford.edu/~tibs/ElemStatLearn/datasets/khan.ytrain")

SRBCTxtest <- read.table("http://www-stat-class.stanford.edu/~tibs/ElemStatLearn/datasets/khan.xtest")
SRBCTytest <- read.table("http://www-stat-class.stanford.edu/~tibs/ElemStatLearn/datasets/khan.ytest")

###Elastic net -- or with alpha=1 the default this is Lasso -- using the multinomial model

SRBCTglmnet <- glmnet(t(SRBCTxtrain),as.factor(t(SRBCTytrain)),family="multinomial")
par(mfcol=c(2,2))
plot(SRBCTglmnet)

###Computing the predictions on the training and test data and
###plotting the curves

SRBCTpredictTrain <- predict(SRBCTglmnet,t(SRBCTxtrain),type="class")
TrainErrors <- apply(apply(SRBCTpredictTrain,2,"!=", SRBCTytrain),2,sum,na.rm=TRUE)
qplot(SRBCTglmnet$lambda,TrainErrors,geom="line")

SRBCTpredict <- predict(SRBCTglmnet,t(SRBCTxtest),type="class")
TestErrors <- apply(apply(SRBCTpredict,2,"!=", SRBCTytest),2,sum,na.rm=TRUE)
qplot(SRBCTglmnet$lambda,TestErrors,geom="line")

###Nicer plot
Errors <- melt(data.frame(lambda=SRBCTglmnet$lambda,TrainErrors = TrainErrors,TestErrors = TestErrors),id="lambda")
qplot(lambda,value,data=Errors,colour=variable,geom="line")


###Extracting the non-null coefficients for s=0.1

SRBCTcoef <- coef(SRBCTglmnet,s=0.1)
SRBCTcoef <- as.data.frame(lapply(SRBCTcoef,as.matrix))
SRBCTcoef[apply(abs(SRBCTcoef),1,sum) > 0,]

###Regularized discriminant analysis

SRBCTrda <- rda(as.matrix(SRBCTxtrain),as.vector(unlist(SRBCTytrain)))