PH240D_hw2

Question 3. Elastic Net: Simulation Study
Part A: Simulation study.

#Create covariance matrix gamma
J<-10; rho<-0.5
element<-function(i,j){rho^abs(i-j)}
gamma<-outer(1:J,1:J,FUN=function(i,j) element(i,j))

#Create vector of beta coefficients
beta<-c(seq(from=-J/2+1,0,by=1),seq(from=0,to=J/2-1,by=1))/10

#Simulate the learning set
library(MASS)
sigma<-2
learn_x<-mvrnorm(n=100,mu=rep(0,10),Sigma=gamma)
learn_y<-sapply(learn_x%*%beta,function(x) rnorm(n=1,mean=x,sd=sigma))

#Here we provide a graphical summary of the true covariance matrix and the generated learning set covariance matrix
library(ggplot2)
library(reshape2)
ggplot(data = melt(gamma[,nrow(gamma):1]), aes(x=Var1, y=Var2, fill=value))+geom_raster()+guides(fill=FALSE)+theme(axis.text=element_blank(),axis.ticks=element_blank(),panel.background=element_blank(),panel.grid.minor=element_blank(),panel.grid.major=element_blank(),axis.title=element_blank(),title=element_text(size=15))+labs(title="Heatmap of Covariance Matrix Gamma")

ggplot(data = melt(cor(learn_x)[,nrow(gamma):1]), aes(x=Var1, y=Var2, fill=value)) + geom_raster()+guides(fill=FALSE)+theme(axis.text=element_blank(),axis.ticks=element_blank(),panel.background=element_blank(),panel.grid.minor=element_blank(),panel.grid.major=element_blank(),axis.title=element_blank(),title=element_text(size=15))+labs(title="Heatmap of Learning Set Covariance Matrix")

#numeric summary of learning set and boxplot to compare distributions
summary(learn_x)

##        V1                 V2                 V3          
##  Min.   :-2.63292   Min.   :-2.23400   Min.   :-2.22181  
##  1st Qu.:-0.69155   1st Qu.:-0.49738   1st Qu.:-0.68790  
##  Median : 0.11672   Median :-0.01375   Median :-0.07168  
##  Mean   : 0.04911   Mean   : 0.04634   Mean   :-0.11148  
##  3rd Qu.: 0.72366   3rd Qu.: 0.57262   3rd Qu.: 0.59712  
##  Max.   : 2.67611   Max.   : 2.95631   Max.   : 1.71558  
##        V4                 V5                 V6           
##  Min.   :-2.99094   Min.   :-3.19108   Min.   :-3.034298  
##  1st Qu.:-0.69310   1st Qu.:-0.55360   1st Qu.:-0.652532  
##  Median : 0.00577   Median : 0.06287   Median : 0.004186  
##  Mean   :-0.03481   Mean   : 0.03848   Mean   :-0.083192  
##  3rd Qu.: 0.63363   3rd Qu.: 0.78209   3rd Qu.: 0.505846  
##  Max.   : 2.34318   Max.   : 2.04255   Max.   : 2.591399  
##        V7                 V8                V9          
##  Min.   :-2.52713   Min.   :-2.8125   Min.   :-2.48339  
##  1st Qu.:-0.63728   1st Qu.:-0.9295   1st Qu.:-0.76183  
##  Median : 0.08059   Median :-0.1253   Median : 0.18988  
##  Mean   : 0.01647   Mean   :-0.1302   Mean   :-0.03533  
##  3rd Qu.: 0.64371   3rd Qu.: 0.6412   3rd Qu.: 0.77769  
##  Max.   : 2.48396   Max.   : 2.3929   Max.   : 1.90354  
##       V10          
##  Min.   :-2.11472  
##  1st Qu.:-0.62930  
##  Median :-0.03642  
##  Mean   :-0.01710  
##  3rd Qu.: 0.60993  
##  Max.   : 2.27727

learn_set<-cbind(learn_x,learn_y)
colnames(learn_set)[1:10]<-paste("x",1:10,sep="")
ggplot(data=melt(learn_set),aes(x=Var2,y=value,colour=Var2))+geom_boxplot(fill="steelblue")+guides(colour=FALSE)+labs(title="Boxplot of Learning Data",x="Learning Covariates and Response",y="Value")+theme(axis.title=element_text(size=12),title=element_text(size=15),axis.text=element_text(size=15))+scale_x_discrete(labels=c(paste("X",1:10,sep=""),"Y"))

#Simulate the testing set:
test_x<-mvrnorm(n=1000,mu=rep(0,10),Sigma=gamma)
test_y<-sapply(test_x%*%beta,function(x) rnorm(n=1,mean=x,sd=sigma))

Part B: Elastic net regression on learning set.

###Create glmnet fit objects for Ridge, LASSO, and Elastic Net###
library(glmnet)

## Loading required package: Matrix
## Loading required package: foreach
## Loaded glmnet 2.0-2

lambda_ridge<-0:100/100
fit_ridge<-glmnet(x=learn_x,y=learn_y,family="gaussian",alpha=0,lambda=lambda_ridge)

lambda_lasso<-0:100/(2*100)
fit_lasso<-glmnet(x=learn_x,y=learn_y,family="gaussian",alpha=1,lambda=lambda_lasso)

lambda_enet<-0:100*(3/(4*100))
fit_enet<-glmnet(x=learn_x,y=learn_y,family="gaussian",alpha=1/3,lambda=lambda_enet)

###Plot effective degrees of freedom vs. lambda for each estimator###
svd_x<-svd(learn_x)$d #get singular values of x

df_ridge<-sapply(0:100,function(x) sum(svd_x^2/((svd_x^2)+x)))
df_ridge_data<-data.frame(cbind(lambda=rev(lambda_ridge),df=rev(df_ridge)))

df_data<-data.frame(cbind(lambda=0:100,df_ridge=rev(df_ridge),df_lasso=fit_lasso$df,df_enet=fit_enet$df))
df_data<-melt(df_data,id.vars="lambda")

ggplot(data=df_data,aes(x=rev(lambda),y=value,colour=variable))+geom_line(size=1.25)+labs(title="Effective Degrees of Freedom",x="Lambda",y="Effective Degrees of Freedom")+scale_colour_discrete(name="Method",labels=c("Ridge","LASSO","Elastic Net"))+theme(axis.title=element_text(size=12),title=element_text(size=15),legend.title=element_text(size=11))

###Plot estimated regression coefficients vs. lambda for each estimator###
estcoef_ridge<-fit_ridge$beta[,101:1]
estcoef_ridge_data<-melt(as.matrix(estcoef_ridge))
estcoef_ridge_data$Var2<-rep(lambda_ridge,each=10)
ggplot(data=estcoef_ridge_data,aes(x=Var2,y=value,colour=Var1))+geom_line(size=1.25)+labs(title="Ridge Regression Solution Path",x="Lambda",y="Estimated Regression Coefficient")+theme(axis.title=element_text(size=12),title=element_text(size=15),legend.title=element_text(size=11))+scale_colour_discrete(name="Coefficient",labels=paste("Beta",1:10,sep=""))

estcoef_lasso<-fit_lasso$beta[,101:1]
estcoef_lasso_data<-melt(as.matrix(estcoef_lasso))
estcoef_lasso_data$Var2<-rep(lambda_lasso,each=10)
ggplot(data=estcoef_lasso_data,aes(x=Var2,y=value,colour=Var1))+geom_line(size=1.25)+labs(title="LASSO Solution Path",x="Lambda",y="Estimated Regression Coefficient")+theme(axis.title=element_text(size=12),title=element_text(size=15),legend.title=element_text(size=11))+scale_colour_discrete(name="Coefficient",labels=paste("Beta",1:10,sep=""))

estcoef_enet<-fit_enet$beta[,101:1]
estcoef_enet_data<-melt(as.matrix(estcoef_enet))
estcoef_enet_data$Var2<-rep(lambda_enet,each=10)
ggplot(data=estcoef_enet_data,aes(x=Var2,y=value,colour=Var1))+geom_line(size=1.25)+labs(title="Elastic Net Solution Path",x="Lambda",y="Estimated Regression Coefficient")+theme(axis.title=element_text(size=12),title=element_text(size=15),legend.title=element_text(size=11))+scale_colour_discrete(name="Coefficient",labels=paste("Beta",1:10,sep=""))

###Plot MSE vs. lambda for each estimator###
mse<-function(i,pred){sum((learn_y - pred[,i])^2)/100}

pred_ridge<-learn_x %*% fit_ridge$beta
mse_ridge<-sapply(seq_along(lambda_ridge),mse,pred=pred_ridge)

pred_lasso<-learn_x %*% fit_lasso$beta
mse_lasso<-sapply(seq_along(lambda_lasso),mse,pred=pred_lasso)

pred_enet<-learn_x %*% fit_enet$beta
mse_enet<-sapply(seq_along(lambda_enet),mse,pred=pred_enet)

mse_data<-data.frame(cbind(lambda=100:0,mse_ridge=mse_ridge,mse_lasso=mse_lasso,mse_enet=mse_enet))
mse_data<-melt(mse_data,id.vars="lambda")
ggplot(data=mse_data,aes(x=lambda,y=value,colour=variable))+geom_line(size=1.25)+labs(title="Learning Set MSE",x="Lambda",y="MSE")+scale_colour_discrete(name="Method",labels=c("Ridge","LASSO","Elastic Net"))+theme(axis.title=element_text(size=12),title=element_text(size=15),legend.title=element_text(size=11),legend.position="bottom")

Part C: Performance assessment on test set.

t_mse<-function(i,pred){sum((test_y - pred[,i])^2)/1000} #function to compute MSE for single value of lambda
#ridge MSE
pred_ridge_t<-test_x %*% fit_ridge$beta
mse_ridge_t<-sapply(seq_along(lambda_ridge),t_mse,pred=pred_ridge_t)

#lasso MSE
pred_lasso_t<-test_x %*% fit_lasso$beta
mse_lasso_t<-sapply(seq_along(lambda_lasso),t_mse,pred=pred_lasso_t)

#elastic net MSE
pred_enet_t<-test_x %*% fit_enet$beta
mse_enet_t<-sapply(seq_along(lambda_enet),t_mse,pred=pred_enet_t)

mse_data_t<-data.frame(cbind(lambda=100:0,mse_ridge_t=mse_ridge_t,mse_lasso_t=mse_lasso_t,mse_enet_t=mse_enet_t))
mse_data_t<-melt(mse_data_t,id.vars="lambda")
ggplot(data=mse_data_t,aes(x=lambda,y=value,colour=variable))+geom_line(size=1.25)+labs(title="Test Set MSE",x="Lambda",y="MSE")+scale_colour_discrete(name="Method",labels=c("Ridge","LASSO","Elastic Net"))+theme(axis.title=element_text(size=12),title=element_text(size=15),legend.title=element_text(size=11),legend.position="bottom")

#optimal estimators for which test set MSE is lowest
min_data_trio<-melt(cbind(ridge=fit_ridge$beta[,which.min(mse_ridge_t)],lass=fit_lasso$beta[,which.min(mse_lasso_t)],enet=fit_enet$beta[,which.min(mse_enet_t)]))
ggplot()+geom_bar(data=min_data_trio,aes(x=Var1,fill=Var2,y=value),position="dodge",stat="identity",colour="steelblue")+labs(title="Coefficient Estimates at Optimal Lambda Value",x="Coefficients",y="Estimated Coefficient Value")+theme(axis.title=element_text(size=12),title=element_text(size=15),legend.title=element_text(size=10),legend.position="bottom")+scale_fill_discrete(name=element_blank(),labels=c("Ridge","LASSO","Elastic Net"))+scale_x_discrete(labels=paste("Beta",1:10,sep=""))+geom_hline(yintercept=0,size=1,colour="steelblue")

Part D: Ridge regression: Bias, variance, and mean squared error of estimated regression coefficients.

#create simulation; we want 1000 iterations of ridge regression estimators
sim_beta<-function(){
  tmp_x<-mvrnorm(n=100,mu=rep(0,10),Sigma=gamma)
  tmp_y<-sapply(tmp_x%*%beta,function(x) rnorm(n=1,mean=x,sd=sigma))
  tmp_ridge<-glmnet(x=tmp_x,y=tmp_y,family="gaussian",alpha=0,lambda=lambda_ridge)
  return(tmp_ridge$beta)
}
sim_data<-replicate(1000,sim_beta())

#calculate bias of ridge regression estimator
bias_fun<-function(x,i){as.matrix(sim_data[[x]])[,i]-beta}
bias_ridge<-matrix(rowMeans(sapply(1:1000,function(x) bias_fun(x,1:101))),nrow=10,ncol=101,byrow=FALSE)
bias_ridge<-data.frame(cbind(t(bias_ridge),lambda=100:0))
bias_ridge_melt<-melt(bias_ridge,id.vars="lambda")
ggplot(data=bias_ridge_melt,aes(x=lambda,y=value,colour=variable))+geom_line(size=1.25)+labs(title="Ridge Estimator Bias (1000 Simulations)",x="Lambda",y="Bias of Estimator")+scale_colour_discrete(name="Coefficients")+theme(axis.title=element_text(size=12),title=element_text(size=15),legend.title=element_text(size=10))

#calculate variance of ridge regression estimator
var_ridge<-sapply(1:10,function(j) apply(sapply(1:1000,function(x) sim_data[[x]][j,]),1,var))
var_ridge<-data.frame(cbind(var_ridge,lambda=100:0))
var_ridge_melt<-melt(var_ridge,id.vars="lambda")
ggplot(data=var_ridge_melt,aes(x=lambda,y=value,colour=variable))+geom_line(size=1.25)+labs(title="Ridge Estimator Variance (1000 Simulations)",x="Lambda",y="Variance of Estimator")+scale_colour_discrete(name="Coefficients")+theme(axis.title=element_text(size=12),title=element_text(size=15),legend.title=element_text(size=10))

#calculate MSE of ridge regression estimator
#function to return mse for single simulation
sim_mse<-function(){
  tmp_x<-mvrnorm(n=100,mu=rep(0,10),Sigma=gamma)
  tmp_y<-sapply(tmp_x%*%beta,function(x) rnorm(n=1,mean=x,sd=sigma))
  tmp_ridge<-glmnet(x=tmp_x,y=tmp_y,family="gaussian",alpha=0,lambda=lambda_ridge)
  tmp_test_x<-mvrnorm(n=100,mu=rep(0,10),Sigma=gamma)
  tmp_test_y<-sapply(tmp_test_x%*%beta,function(x) rnorm(n=1,mean=x,sd=sigma))
  tmp_pred<-tmp_test_x %*% tmp_ridge$beta
  mse_calc<-function(i){sum((tmp_test_y - tmp_pred[,i])^2)/100}
  tmp_mse<-sapply(seq_along(lambda_ridge),mse_calc)
  return(tmp_mse)
}
mse_ridge_sim<-replicate(100,sim_mse())
mse_ridge_sim<-data.frame(cbind(mse_ridge_sim,lambda=100:0))
mse_ridge_sim_melt<-melt(mse_ridge_sim,id.vars="lambda")
ggplot(data=mse_ridge_sim_melt,aes(x=lambda,y=value,colour=variable))+geom_line(size=1.25)+labs(title="Ridge Estimator MSE (100 Simulations)",x="Lambda",y="MSE")+scale_color_discrete(guide=FALSE)+theme(axis.title=element_text(size=12),title=element_text(size=15))

ggplot(data=mse_ridge_sim_melt,aes(x=lambda,y=value,colour=variable))+geom_smooth(aes(group=1),size=1.25,method="loess")+labs(title="Ridge Estimator MSE (Smoothed over 100 Simulations)",x="Lambda",y="MSE")+theme(axis.title=element_text(size=12),title=element_text(size=15))

Part E: LASSO regression: Bias, variance, and mean squared error of estimated regression coefficients.

#create simulation; we want 1000 iterations of LASSO estimators
simL_beta<-function(){
  tmp_x<-mvrnorm(n=100,mu=rep(0,10),Sigma=gamma)
  tmp_y<-sapply(tmp_x%*%beta,function(x) rnorm(n=1,mean=x,sd=sigma))
  tmp_lasso<-glmnet(x=tmp_x,y=tmp_y,family="gaussian",alpha=1,lambda=lambda_lasso)
  return(tmp_lasso$beta)
}
simL_data<-replicate(1000,simL_beta())

#calculate bias of LASSO estimator
biasL_fun<-function(x,i){as.matrix(simL_data[[x]])[,i]-beta}
bias_lasso<-matrix(rowMeans(sapply(1:1000,function(x) biasL_fun(x,1:101))),nrow=10,ncol=101,byrow=FALSE)
bias_lasso<-data.frame(cbind(t(bias_lasso),lambda=100:0))
bias_lasso_melt<-melt(bias_lasso,id.vars="lambda")
ggplot(data=bias_lasso_melt,aes(x=lambda,y=value,colour=variable))+geom_line(size=1.25)+labs(title="LASSO Estimator Bias (1000 Simulations)",x="Lambda",y="Bias of Estimator")+scale_colour_discrete(name="Coefficients")+theme(axis.title=element_text(size=12),title=element_text(size=15),legend.title=element_text(size=10))

#calculate variance of LASSO estimator
var_lasso<-sapply(1:10,function(j) apply(sapply(1:1000,function(x) simL_data[[x]][j,]),1,var))
var_lasso<-data.frame(cbind(var_lasso,lambda=100:0))
var_lasso_melt<-melt(var_lasso,id.vars="lambda")
ggplot(data=var_lasso_melt,aes(x=lambda,y=value,colour=variable))+geom_line(size=1.25)+labs(title="LASSO Estimator Variance (1000 Simulations)",x="Lambda",y="Variance of Estimator")+scale_colour_discrete(name="Coefficients")+theme(axis.title=element_text(size=12),title=element_text(size=15),legend.title=element_text(size=10))

#calculate MSE of LASSO estimator
#function to return mse for single simulation
simL_mse<-function(){
  tmp_x<-mvrnorm(n=100,mu=rep(0,10),Sigma=gamma)
  tmp_y<-sapply(tmp_x%*%beta,function(x) rnorm(n=1,mean=x,sd=sigma))
  tmp_lasso<-glmnet(x=tmp_x,y=tmp_y,family="gaussian",alpha=1,lambda=lambda_lasso)
  tmp_test_x<-mvrnorm(n=100,mu=rep(0,10),Sigma=gamma)
  tmp_test_y<-sapply(tmp_test_x%*%beta,function(x) rnorm(n=1,mean=x,sd=sigma))
  tmp_pred<-tmp_test_x %*% tmp_lasso$beta
  mse_calc<-function(i){sum((tmp_test_y - tmp_pred[,i])^2)/100}
  tmp_mse<-sapply(seq_along(lambda_lasso),mse_calc)
  return(tmp_mse)
}
mse_lasso_sim<-replicate(100,simL_mse())
mse_lasso_sim<-data.frame(cbind(mse_lasso_sim,lambda=100:0))
mse_lasso_sim_melt<-melt(mse_lasso_sim,id.vars="lambda")
ggplot(data=mse_lasso_sim_melt,aes(x=lambda,y=value,colour=variable))+geom_line(size=1.25)+labs(title="LASSO MSE (100 Simulations)",x="Lambda",y="MSE")+scale_color_discrete(guide=FALSE)+theme(axis.title=element_text(size=15),title=element_text(size=15))

ggplot(data=mse_lasso_sim_melt,aes(x=lambda,y=value,colour=variable))+geom_smooth(aes(group=1),size=1.25,method="loess")+labs(title="LASSO Estimator MSE (Smoothed over 100 Simulations)",x="Lambda",y="MSE")+theme(axis.title=element_text(size=12),title=element_text(size=15))

Optimal estimator values

r_plot<-ggplot(data=mse_ridge_sim_melt,aes(x=lambda,y=value,colour=variable))+geom_smooth(aes(group=1),size=1.25,method="loess")+labs(title="Ridge Estimator MSE (Smoothed over 100 Simulations)",x="Lambda",y="MSE")+theme(axis.title=element_text(size=15),title=element_text(size=20))

loess_fitr<-ggplot_build(r_plot)$data[[1]][,1:2]
min_ridge<-floor(loess_fitr[,1][which.min(loess_fitr[,2])])

l_plot<-ggplot(data=mse_lasso_sim_melt,aes(x=lambda,y=value,colour=variable))+geom_smooth(aes(group=1),size=1.25,method="loess")+labs(title="LASSO Estimator MSE (Smoothed over 100 Simulations)",x="Lambda",y="MSE")+theme(axis.title=element_text(size=15),title=element_text(size=20))

loess_fitl<-ggplot_build(l_plot)$data[[1]][,1:2]
min_lasso<-floor(loess_fitl[,1][which.min(loess_fitl[,2])])

min_data<-melt(cbind(ridge=fit_ridge$beta[,101-min_ridge],lasso=fit_lasso$beta[,101-min_lasso]))
ggplot()+geom_bar(data=min_data,aes(x=Var1,fill=Var2,y=value),position="dodge",stat="identity",colour="steelblue")+labs(title="Coefficient Estimates at Optimal Lambda Value",x="Coefficients",y="Estimated Coefficient Value")+theme(axis.title=element_text(size=12),title=element_text(size=15),legend.title=element_text(size=10),legend.position="bottom")+scale_fill_discrete(name=element_blank(),labels=c("Ridge","LASSO"))+scale_x_discrete(labels=paste("Beta",1:10,sep=""))+geom_hline(yintercept=0,size=1,colour="steelblue")