##sparse logistic regression

rm(list=ls())

library(glmnet)
set.seed(95373)


##make data
n = 200
d = 500
beta = c(2,2,2,2,2,rep(0,d-5))
beta = matrix(beta,ncol=1)
X = matrix(rnorm(n*d),n,d)
mu = X %*% beta
p = exp(mu)/(1+exp(mu))
Y = rbinom(n,size=1,prob=p)


## Run the logistic lasso

pdf("Sparse_Logistic1.pdf")

##cross validation

fit.cv = cv.glmnet(X,Y,family="binomial")
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=4,ylim=c(.5,max(up)),
             xlab="",ylab="",xaxt="n",cex.axis=1.3)
polygon(c(L,rev(L)),c(low,rev(up)),col="plum3",border=NA)
lines(L,cv,col="purple",lwd=4)
lines(L,cv+se,col="purple",lty=2,lwd=3)
lines(L,cv-se,col="purple",lty="dashed",lwd=3)
I = which.min(cv)
segments(L[I],.5-.003,L[I],min(low),col="red",lty=2,lwd=4)
title(ylab="Prediction Error",line=2.7, cex.lab=1.9)
title(xlab=expression(log(lambda)),line=1.3, cex.lab=1.8,adj=.9)
axis(side = 1, at = round(L[I],2),cex.axis=1.8,cex.lab=2,col="red")

dev.off()



pdf("Sparse_Logistic2.pdf")

##plot nonzero betas
b = as.matrix(coef(fit.cv))
L = fit.cv$lambda.min
i = which(fit.cv$lambda == L)
plot(b,type="h",lwd=2,col="purple",cex.axis=1.3,xlab="",ylab="")
abline(h=0, col="red",lwd=2)
title(ylab=expression(beta),line=2.4,cex.lab=1.9)
title(xlab="Index",line=3, cex.lab=1.8)
dev.off()