##################################################
# this script was used to produce Tables 4 and 5 
# and boxplots in the paper "The predictive Lasso"
##################################################

# small-p case

n = 500
ntest = n
nrep = 500						# number of replications	

p = 8
b0 = 2						# intercept
t.beta1 = as.matrix(c(3,1.5,0,0,.5,0.5,0,0))
true.beta = c(b0,t.beta1)
trueflag = abs(true.beta)>0
Sig = diag(p);
for (ii in 1:p)
	for (jj in 1:p)
		if (ii!=jj) Sig[ii,jj] = 0.5^(abs(ii-jj))

PPSlasso = PPSalasso = PPSnplasso = PPSplasso = rep(0,nrep)
MSElasso = MSEalasso = MSEnplasso = MSEplasso = rep(0,nrep)
NZElasso = NZEalasso = NZEnplasso = NZEplasso = rep(0,nrep)
CRPSlasso = CRPSalasso = CRPSnplasso = CRPSplasso = rep(0,nrep)
CFRlasso = CFRalasso = CFRnplasso = CFRplasso = 0

N = n+ntest

for (irep in 1:nrep) {

x = mvrnorm (N, rep(0,p),Sig)
y = rbinom (N, 1, invlogit(b0+x%*%t.beta1))
xtest = x[(ntest+1):N,]
ytest = y[(ntest+1):N]
x = x[1:ntest,]
y = y[1:ntest]

# handle Lasso

beta.lasso = lasso_logistic(x,y)
xb = cbind(as.matrix(rep(1,ntest)),xtest)%*%beta.lasso
PPSlasso[irep] = (-ytest%*%xb+ sum(log(1+exp(xb))))/ntest
MSElasso[irep] = sum((true.beta-beta.lasso)^2)
NZElasso[irep] = p+1-sum(abs(beta.lasso)>0)
CFRlasso = CFRlasso+as.numeric(sum((abs(beta.lasso)>0)==trueflag)==p+1)
muhat = exp(xb)/(1+exp(xb))
CRPSlasso[irep] = (sum((muhat[ytest==0])^2)+sum((1-muhat[ytest==1])^2))/ntest

# handle aLasso

beta.alasso = alasso_logistic(x,y)
xb = cbind(as.matrix(rep(1,ntest)),xtest)%*%beta.alasso
PPSalasso[irep] = (-ytest%*%xb+ sum(log(1+exp(xb))))/ntest
MSEalasso[irep] = sum((true.beta-beta.alasso)^2)
NZEalasso[irep] = p+1-sum(abs(beta.alasso)>0)
CFRalasso = CFRalasso+as.numeric(sum((abs(beta.alasso)>0)==trueflag)==p+1)
muhat = exp(xb)/(1+exp(xb))
CRPSalasso[irep] = (sum((muhat[ytest==0])^2)+sum((1-muhat[ytest==1])^2))/ntest

# handle npLasso

beta.nplasso = nplasso_logistic(x,y)
xb = cbind(as.matrix(rep(1,ntest)),xtest)%*%beta.nplasso
PPSnplasso[irep] = (-ytest%*%xb+ sum(log(1+exp(xb))))/ntest
MSEnplasso[irep] = sum((true.beta-beta.nplasso)^2)
NZEnplasso[irep] = p+1-sum(abs(beta.nplasso)>0)
CFRnplasso = CFRnplasso+as.numeric(sum((abs(beta.nplasso)>0)==trueflag)==p+1)
muhat = exp(xb)/(1+exp(xb))
CRPSnplasso[irep] = (sum((muhat[ytest==0])^2)+sum((1-muhat[ytest==1])^2))/ntest

# handle pLasso

beta.plasso = plasso_logistic(x,y)
xb = cbind(as.matrix(rep(1,ntest)),xtest)%*%beta.plasso
PPSplasso[irep] = (-ytest%*%xb+ sum(log(1+exp(xb))))/ntest
MSEplasso[irep] = sum((true.beta-beta.plasso)^2)
NZEplasso[irep] = p+1-sum(abs(beta.plasso)>0)
CFRplasso = CFRplasso+as.numeric(sum((abs(beta.plasso)>0)==trueflag)==p+1)
muhat = exp(xb)/(1+exp(xb))
CRPSplasso[irep] = (sum((muhat[ytest==0])^2)+sum((1-muhat[ytest==1])^2))/ntest

}
cbind(mean(PPSlasso),sd(PPSlasso),mean(PPSalasso),sd(PPSalasso),
	mean(PPSnplasso),sd(PPSnplasso),mean(PPSplasso),sd(PPSplasso))
cbind(mean(MSElasso),sd(MSElasso),mean(MSEalasso),sd(MSEalasso),
	mean(MSEnplasso),sd(MSEnplasso),mean(MSEplasso),sd(MSEplasso))
cbind(mean(NZElasso),sd(NZElasso),mean(NZEalasso),sd(NZEalasso),
	mean(NZEnplasso),sd(NZEnplasso),mean(NZEplasso),sd(NZEplasso))
cbind(mean(CRPSlasso),sd(CRPSlasso),mean(CRPSalasso),sd(CRPSalasso),
	mean(CRPSnplasso),sd(CRPSnplasso),mean(CRPSplasso),sd(CRPSplasso))
cbind(CFRlasso/nrep,CFRalasso/nrep,CFRnplasso/nrep,CFRplasso/nrep)

# boxplots

par(mfrow=c(2,2))
Lasso = PPSlasso; aLasso = PPSalasso; npLasso = PPSnplasso; pLasso = PPSplasso
PPS = cbind(Lasso,aLasso,npLasso,pLasso)
boxplot(PPS)
title('PPS')

Lasso = MSElasso; aLasso = MSEalasso; npLasso = MSEnplasso; pLasso = MSEplasso
MSE = cbind(Lasso,aLasso,npLasso,pLasso)
boxplot(MSE)
title('MSE')

Lasso = CRPSlasso; aLasso = CRPSalasso; npLasso = CRPSnplasso; pLasso = CRPSplasso
CRPS = cbind(Lasso,aLasso,npLasso,pLasso)
boxplot(CRPS)
title('CRPS')

Lasso = NZElasso; aLasso = NZEalasso; npLasso = NZEnplasso; pLasso = NZEplasso
NZE = cbind(Lasso,aLasso,npLasso,pLasso)
boxplot(NZE)
title('NZE')




#======================= large p case ==========================')
# large-p case

n = 200
ntest = n
nrep = 200

p = 100
b0 = 2
t.beta1 = as.matrix(c(3,1.5,.5,0.5,0,rep(0,p-5)))
true.beta = c(b0,t.beta1)
trueflag = abs(true.beta)>0
Sig = diag(p);
for (ii in 1:p)
	for (jj in 1:p)
		if (ii!=jj) Sig[ii,jj] = 0.5^(abs(ii-jj))


PPSlasso = PPSalasso = PPSnplasso = PPSplasso = rep(0,nrep)
MSElasso = MSEalasso = MSEnplasso = MSEplasso = rep(0,nrep)
NZElasso = NZEalasso = NZEnplasso = NZEplasso = rep(0,nrep)
CRPSlasso = CRPSalasso = CRPSnplasso = CRPSplasso = rep(0,nrep)
CFRlasso = CFRalasso = CFRnplasso = CFRplasso = 0

N = n+ntest

for (irep in 1:nrep) {

x = mvrnorm (N, rep(0,p),Sig)
y = rbinom (N, 1, invlogit(b0+x%*%t.beta1))
xtest = x[(ntest+1):N,]
ytest = y[(ntest+1):N]
x = x[1:ntest,]
y = y[1:ntest]

# handle Lasso

beta.lasso = lasso_logistic(x,y)
xb = cbind(as.matrix(rep(1,ntest)),xtest)%*%beta.lasso
PPSlasso[irep] = (-ytest%*%xb+ sum(log(1+exp(xb))))/ntest
MSElasso[irep] = sum((true.beta-beta.lasso)^2)
NZElasso[irep] = p+1-sum(abs(beta.lasso)>0)
CFRlasso = CFRlasso+as.numeric(sum((abs(beta.lasso)>0)==trueflag)==p+1)
muhat = exp(xb)/(1+exp(xb))
CRPSlasso[irep] = (sum((muhat[ytest==0])^2)+sum((1-muhat[ytest==1])^2))/ntest

# handle aLasso

beta.alasso = alasso_logistic(x,y)
xb = cbind(as.matrix(rep(1,ntest)),xtest)%*%beta.alasso
PPSalasso[irep] = (-ytest%*%xb+ sum(log(1+exp(xb))))/ntest
MSEalasso[irep] = sum((true.beta-beta.alasso)^2)
NZEalasso[irep] = p+1-sum(abs(beta.alasso)>0)
CFRalasso = CFRalasso+as.numeric(sum((abs(beta.alasso)>0)==trueflag)==p+1)
muhat = exp(xb)/(1+exp(xb))
CRPSalasso[irep] = (sum((muhat[ytest==0])^2)+sum((1-muhat[ytest==1])^2))/ntest

# handle npLasso

beta.nplasso = nplasso_logistic(x,y)
xb = cbind(as.matrix(rep(1,ntest)),xtest)%*%beta.nplasso
PPSnplasso[irep] = (-ytest%*%xb+ sum(log(1+exp(xb))))/ntest
MSEnplasso[irep] = sum((true.beta-beta.nplasso)^2)
NZEnplasso[irep] = p+1-sum(abs(beta.nplasso)>0)
CFRnplasso = CFRnplasso+as.numeric(sum((abs(beta.nplasso)>0)==trueflag)==p+1)
muhat = exp(xb)/(1+exp(xb))
CRPSnplasso[irep] = (sum((muhat[ytest==0])^2)+sum((1-muhat[ytest==1])^2))/ntest

# handle pLasso

beta.plasso = plasso_logistic(x,y)
xb = cbind(as.matrix(rep(1,ntest)),xtest)%*%beta.plasso
PPSplasso[irep] = (-ytest%*%xb+ sum(log(1+exp(xb))))/ntest
MSEplasso[irep] = sum((true.beta-beta.plasso)^2)
NZEplasso[irep] = p+1-sum(abs(beta.plasso)>0)
CFRplasso = CFRplasso+as.numeric(sum((abs(beta.plasso)>0)==trueflag)==p+1)
muhat = exp(xb)/(1+exp(xb))
CRPSplasso[irep] = (sum((muhat[ytest==0])^2)+sum((1-muhat[ytest==1])^2))/ntest

}
cbind(mean(PPSlasso),sd(PPSlasso),mean(PPSalasso),sd(PPSalasso),
	mean(PPSnplasso),sd(PPSnplasso),mean(PPSplasso),sd(PPSplasso))
cbind(mean(MSElasso),sd(MSElasso),mean(MSEalasso),sd(MSEalasso),
	mean(MSEnplasso),sd(MSEnplasso),mean(MSEplasso),sd(MSEplasso))
cbind(mean(NZElasso),sd(NZElasso),mean(NZEalasso),sd(NZEalasso),
	mean(NZEnplasso),sd(NZEnplasso),mean(NZEplasso),sd(NZEplasso))
cbind(mean(CRPSlasso),sd(CRPSlasso),mean(CRPSalasso),sd(CRPSalasso),
	mean(CRPSnplasso),sd(CRPSnplasso),mean(CRPSplasso),sd(CRPSplasso))
cbind(CFRlasso/nrep,CFRalasso/nrep,CFRnplasso/nrep,CFRplasso/nrep)

# boxplots

par(mfrow=c(2,2))
Lasso = PPSlasso; aLasso = PPSalasso; npLasso = PPSnplasso; pLasso = PPSplasso
PPS = cbind(Lasso,aLasso,npLasso,pLasso)
boxplot(PPS)
title('PPS')

Lasso = MSElasso; aLasso = MSEalasso; npLasso = MSEnplasso; pLasso = MSEplasso
MSE = cbind(Lasso,aLasso,npLasso,pLasso)
boxplot(MSE)
title('MSE')

Lasso = CRPSlasso; aLasso = CRPSalasso; npLasso = CRPSnplasso; pLasso = CRPSplasso
CRPS = cbind(Lasso,aLasso,npLasso,pLasso)
boxplot(CRPS)
title('CRPS')

Lasso = NZElasso; aLasso = NZEalasso; npLasso = NZEnplasso; pLasso = NZEplasso
NZE = cbind(Lasso,aLasso,npLasso,pLasso)
boxplot(NZE)
title('NZE')
