Supervised ML: Lasso saves the job of OLS

Goals:

• Illustrate how Lasso works

• Illustrate that Lasso does not overfit

Data generating process

We use the same DGP as in the notebook Overfitting of OLS and value of training vs. test sample:

• \(p\geq10\) independent and standard normal covariates: \(X \sim N(0,I_p)\), where \(I_p\) is the \(p\)-dimensional identity matrix

• The outcome model is \(Y = \underbrace{\beta_0 + \beta_1 X_1 + ... + \beta_{10} X_{10}}_{\text{conditional expectation function }m(X)} + e\), where \(X_j\) is the \(j\)-th column of \(X\) and \(e \sim N(0,1)\)

• We consider the following parameters \(\beta_0 = 0\), \(\beta_1 = 1\), \(\beta_2 = 0.9\), …, \(\beta_9 = 0.2\), \(\beta_{10} = 0.1\) such that the first 10 variables have a decreasing impact on the outcome and any further variables are just irrelevant noise with no true predictive power

We set \(p=99\) and draw a training sample with \(N_{tr}=100\) and a test sample with \(N_{te}=10,000\)

# Load the packages required for later
if (!require("tidyverse")) install.packages("tidyverse", dependencies = TRUE); library(tidyverse)
if (!require("glmnet")) install.packages("glmnet", dependencies = TRUE); library(glmnet)
if (!require("hdm")) install.packages("hdm", dependencies = TRUE); library(hdm)
if (!require("plasso")) install.packages("plasso", dependencies = TRUE); library(plasso)

set.seed(1234) # For replicability

# Define the important parameters
n_tr = 100
n_te = 10000
p = 99
beta = c(0,seq(1,0.1,-0.1),rep(0,p-10))

# Combine constant and randomly drawn covariates
x_tr = cbind(rep(1,n_tr),matrix(rnorm(n_tr*p),ncol=p))
x_te = cbind(rep(1,n_te),matrix(rnorm(n_te*p),ncol=p))

# Create the CEF using matrix multiplication for compactness
cfe_tr = x_tr %*% beta
cfe_te = x_te %*% beta

# Create the "observed" outcomes by adding noise
y_tr = cfe_tr + rnorm(n_tr,0,1)
y_te = cfe_te + rnorm(n_te,0,1)

Lasso at work

We run Lasso with all 99 regressors. The glmnet command calculates the trajectory from an empty model to a nearly unpenalized model (read from right to left):

# As most of the commands that we use from now on, 
# glmnet takes inputs in matrix form and not as formulas
lasso = glmnet(x_tr,y_tr)
plot(lasso, xvar = "lambda")

In the beginning the coefficients of the really relevant variables are build up and at some point the excessive overfitting starts.

Cross-validation helps to figure out the sweet spot that ensures a good out-of-sample performance. The cv.glmnet command uses by default 10-fold cross-validation leading to the following result:

cv_lasso = cv.glmnet(x_tr,y_tr)
plot(cv_lasso)

We observe that the cross-validated MSE decreases as the relevant variables are selected and build up. At some point the reduced penalization leads to selection of noise variables and the cross-validated MSE deteriorated again. We select the penalty term where the curve indicates the lowest MSE (alternative is the 1SE rule indicated by the right dashed line, but we don’t go into the details here).

Post-Lasso at work

Cross-validated

Cross-validation for Post-Lasso is implemented in the plasso package:

# Provide the cov matrix w/o the constant x_tr[,-1]
# Increasing lambda.min.ratio ensures that Lasso does not overfit too heavily and reduces running time
post_lasso = plasso(x_tr[,-1],y_tr,lambda.min.ratio=0.01)
plot(post_lasso,xvar = "lambda")

While Lasso builds up coefficients gradually, Post-Lasso gives the full OLS coefficient as soon as the variable is selected \(\Rightarrow\) the coefficients paths of Post-Lasso are usually not smooth.

This is also observable in the cross-validation curve:

cv_plasso = cv.plasso(x_tr,y_tr,lambda.min.ratio=0.01)
plot(cv_plasso, legend_pos="bottomleft")

The comparison of Lasso and Post-Lasso cross-validation is instructive to understand some differences:

• The Post-Lasso cross-validation curve is more bumpy and has flat regions compared to the Lasso curve. Post-Lasso gives the full OLS coefficient as soon as the variable is selected. This explains the flat regions of the Post-Lasso curve where no additional variable is selected and thus MSE stays constant.

• Post-Lasso usually produces sparser cross-validated models because the unshrunken coefficients deliver more explanatory power with a smaller set of variables.

Fast implementation with hdm

While the cross-validated Post-Lasso is instructive, it is not fast and the rlasso function of the hdm package implements a faster way to choose the penalty parameter in a data-driven way:

post_hdm = rlasso(x_tr,y_tr)
summary(post_hdm)

Call:
rlasso.default(x = x_tr, y = y_tr)

Post-Lasso Estimation:  TRUE 

Total number of variables: 100
Number of selected variables: 6 

Residuals: 
     Min       1Q   Median       3Q      Max 
-2.78672 -0.76531 -0.04826  0.63651  3.37613 

            Estimate
(Intercept)   -0.097
V1             0.000
V2             0.994
V3             0.887
V4             0.833
V5             0.781
V6             0.692
V7             0.499
V8             0.000
V9             0.000
V10            0.000
V11            0.000
V12            0.000
V13            0.000
V14            0.000
V15            0.000
V16            0.000
V17            0.000
V18            0.000
V19            0.000
V20            0.000
V21            0.000
V22            0.000
V23            0.000
V24            0.000
V25            0.000
V26            0.000
V27            0.000
V28            0.000
V29            0.000
V30            0.000
V31            0.000
V32            0.000
V33            0.000
V34            0.000
V35            0.000
V36            0.000
V37            0.000
V38            0.000
V39            0.000
V40            0.000
V41            0.000
V42            0.000
V43            0.000
V44            0.000
V45            0.000
V46            0.000
V47            0.000
V48            0.000
V49            0.000
V50            0.000
V51            0.000
V52            0.000
V53            0.000
V54            0.000
V55            0.000
V56            0.000
V57            0.000
V58            0.000
V59            0.000
V60            0.000
V61            0.000
V62            0.000
V63            0.000
V64            0.000
V65            0.000
V66            0.000
V67            0.000
V68            0.000
V69            0.000
V70            0.000
V71            0.000
V72            0.000
V73            0.000
V74            0.000
V75            0.000
V76            0.000
V77            0.000
V78            0.000
V79            0.000
V80            0.000
V81            0.000
V82            0.000
V83            0.000
V84            0.000
V85            0.000
V86            0.000
V87            0.000
V88            0.000
V89            0.000
V90            0.000
V91            0.000
V92            0.000
V93            0.000
V94            0.000
V95            0.000
V96            0.000
V97            0.000
V98            0.000
V99            0.000
V100           0.000

Residual standard error: 1.061
Multiple R-squared:  0.7991
Adjusted R-squared:  0.7861
Joint significance test:
 the sup score statistic for joint significance test is 11.64 with a p-value of 0.004

OLS vs. (Post-)Lasso

We replicate the analysis of Notebook Overfitting of OLS and value of training vs. test sample and gradually add covariates to check how out-of-sample performance develops.

# Container of the results
results_ols = results_lasso = results_plasso = results_rlasso = matrix(NA,p-1,4)
colnames(results_ols) = colnames(results_lasso) = 
                        colnames(results_plasso) = 
                        colnames(results_rlasso) = c("Obs MSE train","Obs MSE test",
                                                    "Oracle MSE train","Oracle MSE test")

# Loop that gradually adds variables (start with 2, otherwise glmnet crashes)
for (i in 2:p) {
  # OLS
  temp_ols = lm(y_tr ~ x_tr[,2:(i+1)])
  temp_yhat_tr = predict(temp_ols)
  temp_yhat_te = x_te[,1:(i+1)] %*% temp_ols$coefficients
  # Calculate the observable MSEs in training and test sample
  results_ols[i-1,1] = mean((y_tr - temp_yhat_tr)^2) # in-sample MSE
  results_ols[i-1,2] = mean((y_te - temp_yhat_te)^2) # out-of-sample MSE
  # Calculate the oracle MSEs that are only observables b/c we know the true CEF
  results_ols[i-1,3] = var(y_tr - cfe_tr) + mean((cfe_tr - temp_yhat_tr)^2)
  results_ols[i-1,4] = var(y_te - cfe_te) + mean((cfe_te - temp_yhat_te)^2)
  
  # Lasso
  temp_lasso = cv.glmnet(x_tr[,2:(i+1)],y_tr)
  temp_yhat_tr = predict(temp_lasso,newx=x_tr[,2:(i+1)])
  temp_yhat_te = predict(temp_lasso,newx=x_te[,2:(i+1)])
  # Calculate the observable MSEs in training and test sample
  results_lasso[i-1,1] = mean((y_tr - temp_yhat_tr)^2) # in-sample MSE
  results_lasso[i-1,2] = mean((y_te - temp_yhat_te)^2) # out-of-sample MSE
  # Calculate the oracle MSEs that are only observables b/c we know the true CEF
  results_lasso[i-1,3] = var(y_tr - cfe_tr) + mean((cfe_tr - temp_yhat_tr)^2)
  results_lasso[i-1,4] = var(y_te - cfe_te) + mean((cfe_te - temp_yhat_te)^2)
  
  # plasso
  temp_plasso = cv.plasso(x_tr[,2:(i+1)],y_tr)
  temp_yhat_tr = predict(temp_plasso,newx=x_tr[,2:(i+1)])$plasso
  temp_yhat_te = predict(temp_plasso,newx=x_te[,2:(i+1)])$plasso
  # Calculate the observable MSEs in training and test sample
  results_plasso[i-1,1] = mean((y_tr - temp_yhat_tr)^2) # in-sample MSE
  results_plasso[i-1,2] = mean((y_te - temp_yhat_te)^2) # out-of-sample MSE
  # Calculate the oracle MSEs that are only observables b/c we know the true CEF
  results_plasso[i-1,3] = var(y_tr - cfe_tr) + mean((cfe_tr - temp_yhat_tr)^2)
  results_plasso[i-1,4] = var(y_te - cfe_te) + mean((cfe_te - temp_yhat_te)^2)
  
  # rlasso
  temp_rlasso = rlasso(x_tr[,2:(i+1)],y_tr)
  temp_yhat_tr = predict(temp_rlasso,newdata=x_tr[,2:(i+1)])
  temp_yhat_te = predict(temp_rlasso,newdata=x_te[,2:(i+1)])
  # Calculate the observable MSEs in training and test sample
  results_rlasso[i-1,1] = mean((y_tr - temp_yhat_tr)^2) # in-sample MSE
  results_rlasso[i-1,2] = mean((y_te - temp_yhat_te)^2) # out-of-sample MSE
  # Calculate the oracle MSEs that are only observables b/c we know the true CEF
  results_rlasso[i-1,3] = var(y_tr - cfe_tr) + mean((cfe_tr - temp_yhat_tr)^2)
  results_rlasso[i-1,4] = var(y_te - cfe_te) + mean((cfe_te - temp_yhat_te)^2)
}

Again, OLS explodes but the Lassos show a stable prediction performance regardless of the number of noise variables.

df = data.frame(Estimator = c(rep("OLS",p-1),rep("Lasso",p-1),rep("Post-Lasso CV",p-1),rep("Post-Lasso hdm",p-1)),
                Number.of.variables = c(2:p,2:p,2:p,2:p),
                Obs.MSE.test = c(results_ols[,2], results_lasso[,2], results_plasso[,2], results_rlasso[,2]),
                Oracle.MSE.test = c(results_ols[,4], results_lasso[,4], results_plasso[,4],results_plasso[,4]))

ggplot(subset(df), aes(x=Number.of.variables,y=Obs.MSE.test,colour=Estimator)) + geom_line(size=1)

Warnung: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
Please use `linewidth` instead.

ggplot(subset(df,df$Number.of.variables<70), aes(x=Number.of.variables,y=Obs.MSE.test,colour=Estimator))     + geom_line(size=1)  + geom_hline(yintercept = 0)

It seems that Post-Lasso has a slight advantage over plain Lasso. However, this is DGP dependent and could be flipped in other settings.

Take-away

• The penalization of coefficients ensures good out-of-sample performance \(\Rightarrow\) Lasso works

Suggestions to play with the toy model

Feel free to play around with the code. This is useful to sharpen and challenge your understanding of the methods. Think about the consequences of a modifications before you run it and check whether the results are in line with your expectation. Some suggestions:

• Modify DGP (change betas, change level of noise, introduce correlation between covariates, …)

• Modify seed

• Change training and test sample size

• Investigate the 1SE rule for cross-validated Lasso