The R codes for the baseball example and the data are taken from the ISLR book by James, Witten, Hastie and Tibshirani.

This is a low-dimensional example, which means \(p < n\). We will see how LASSO and Ridge compares with OLS on this data-set.

We will see a high-dimensional example after this. Of course, there isn’t a unique OLS for high-dimensional data.

Major League Baseball Data from the 1986 and 1987 seasons.

library(ISLR)
#fix(Hitters)
names(Hitters)

##  [1] "AtBat"     "Hits"      "HmRun"     "Runs"      "RBI"       "Walks"    
##  [7] "Years"     "CAtBat"    "CHits"     "CHmRun"    "CRuns"     "CRBI"     
## [13] "CWalks"    "League"    "Division"  "PutOuts"   "Assists"   "Errors"   
## [19] "Salary"    "NewLeague"

dim(Hitters)

## [1] 322  20

Salary is missing for some players. The na.omit() function removes all of the rows that have missing values in any variable.

sum(is.na(Hitters$Salary))

## [1] 59

Hitters=na.omit(Hitters)
dim(Hitters)

## [1] 263  20

• Ridge, Lasso, Elastic Net are all based on R package glmnet.

• Elastic Net Penalty: \[ P_{\lambda,\alpha}(\beta) = \lambda \left \{ \frac{(1-\alpha)}{2} {\lVert \beta \rVert}_2^2 + \alpha {\lVert \beta \rVert}_1 \right \} \]

• \(\alpha = 0\) is Ridge, \(\alpha = 1\) is Lasso, anything in between (\(0 < \alpha < 1\) is called Elastic Net.

• \(\lambda = 0\) is ordinary least squares.

• ‘Elastic Net’ is a mix of Lasso and Ridge - it works better when the predictor variables are ‘grouped’.

• All the three methods can be fit using the glmnet package.

x=model.matrix(Salary~.,Hitters)[,-1]
y=Hitters$Salary
dim(x)

## [1] 263  19

length(y)

## [1] 263

• We’ll specify \(\alpha = 0\) and calculate the ridge solution for \(100\) different values of \(\lambda\), chosen by us.

library(glmnet)
grid=10^seq(10,-2,length=100)
ridge.mod=glmnet(x,y,alpha=0,lambda=grid)
dim(coef(ridge.mod))

## [1]  20 100

• We will use cross-validation to choose the best \(\lambda\) that minimizes the test error.

• We have discussed LOO-CV and \(k\)-fold CV in class. In general, LOO-CV has higher variance but \(k\)-fold is more stable.

• We’ll first use a single validation set and then show CV.

• Also, note that if you are splitting randomly, you should fix the seed to reproduce your results.

• Now split the samples into a training set and a test set in order to estimate the test error of ridge regression.

• Two ways to split:

  1. Produce a random vector of TRUE, FALSE elements and select the observations corresponding to TRUE for the training data.

  2. Randomly choose a subset of numbers between 1 and n - these can then be used as the indices for the training observations.

set.seed(1)
train=sample(1:nrow(x), nrow(x)/2)
test=(-train)
y.test=y[test]
dim(x[train,])

## [1] 131  19

dim(x[test,])

## [1] 132  19

• Next we fit a ridge regression model on the training set, and evaluate its MSE on the test set, using \(\lambda = 4\).
• The predict() function : here we get predictions for a test set, by replacing type="coefficients" with the newx argument.
• Specify s = 4 to fix \(\lambda = 4\).

ridge.mod=glmnet(x[train,],y[train],alpha=0,lambda=grid, thresh=1e-12)
ridge.pred=predict(ridge.mod,s=4,newx=x[test,])
mean((ridge.pred-y.test)^2)

## [1] 142199.2

The test MSE is 1.421991510^{5}

If we had instead simply fit a model with just an intercept, we would have predicted each test observation using the mean of the training observations:

mean((mean(y[train])-y.test)^2)

## [1] 224669.9

• We could also get the same result by fitting a ridge regression model with a very large value of \(\lambda\). Note that s=1e10 means \(\lambda = 10^{10}\).

ridge.pred=predict(ridge.mod,s=1e10,newx=x[test,])
mean((ridge.pred-y.test)^2)

## [1] 224669.8

• So fitting a ridge regression model with \(\lambda = 4\) leads to a much lower test MSE than fitting a model with just an intercept.

• Is it better than just performing least squares regression?
• Ridge regression reduces to OLS for \(\lambda = 0\).

ridge.pred=predict(ridge.mod,s=0,newx=x[test,])
mean((ridge.pred-y.test)^2)

## [1] 167789.8

• Ridge Regression MSE was 101036. OLS is worse.

• In general, instead of arbitrarily choosing \(\lambda = 4\), it would be better to use cross-validation to choose the tuning parameter.

• We can do this using the built-in cross-validation function, cv.glmnet().

• By default, the function cv.glmnet() performs ten-fold cross-validation, though this can be changed using the argument nfolds.

set.seed(1) 
cv.out=cv.glmnet(x[train,],y[train],alpha=0)
plot(cv.out)

bestlam=cv.out$lambda.min
bestlam

## [1] 326.0828

The value of \(\lambda\) that yields lowest CV error is 326.0827885, very different from 4.

We should calculate the test error at this “best” \(\lambda\) and check if it’s improving the fit.

bestlam=cv.out$lambda.min
bestlam

## [1] 326.0828

ridge.pred=predict(ridge.mod,s=bestlam,newx=x[test,])
mean((ridge.pred-y.test)^2)

## [1] 139856.6

• It does reduce the test MSE.

• Recall: Ridge, Lasso, Elastic Net are all based on R package glmnet.

• Elastic Net Penalty: \[ P_{\lambda,\alpha}(\beta) = \lambda \left \{ \frac{(1-\alpha)}{2} {\lVert \beta \rVert}_2^2 + \alpha {\lVert \beta \rVert}_1 \right \} \]

• \(\alpha = 0\) is Ridge, \(\alpha = 1\) is Lasso, anything in between is Elastic Net.

• All we need to do is change alpha = 1 in glmnet function call.

lasso.mod =glmnet (x[train ,], y[train], alpha = 1, lambda = grid)
plot(lasso.mod)

## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to unique
## 'x' values

set.seed (1)
cv.out =cv.glmnet( x[train ,], y[train], alpha = 1)
plot(cv.out)

bestlam =cv.out$lambda.min
lasso.pred=predict(lasso.mod, s = bestlam , newx = x[test,])
mean((lasso.pred - y.test)^2)

## [1] 143673.6

This is substantially lower than the test set MSE of the null model and of least squares, and very similar to the test MSE of ridge regression with \(\lambda\) chosen by cross-validation.

Here we see that 12 of the 19 coefficient estimates are exactly zero. So the lasso model with \(\lambda\) chosen by cross-validation contains only seven variables.

lasso.coef[lasso.coef !=0]

##   (Intercept)         AtBat          Hits         Walks         Years 
##    1.27479059   -0.05497143    2.18034583    2.29192406   -0.33806109 
##        CHmRun         CRuns          CRBI       LeagueN     DivisionW 
##    0.02825013    0.21628385    0.41712537   20.28615023 -116.16755870 
##       PutOuts        Errors 
##    0.23752385   -0.85629148

• Generate \(X\), \(y\) with \(n = 30, p = 60\) with a sparse \(\beta\) vector with 5 non-zero elements.
• Generate \(X\) matrix as \(n \times p\) matrix with all entries from a \(N(0,1)\) distribution.
• Take \(\beta = (3,\ldots,3, 0, \ldots, 0)\), and a small error distribution \(\epsilon \sim N(0,0.1)\).

library(Matrix)
set.seed(12345)
n = 30;p = 60
p1 = 5
beta <- c(rep(3,p1),rep(0,p-p1))
x=scale(matrix(rnorm(n*p),n,p))
eps=rnorm(n,mean=0,sd = 0.1)
fx = x %*% beta
y=drop(fx+eps)

• Our goal is to see if Lasso can recover the 5 non-zero \(\beta\)’s correctly.

• The OLS Solution will not work here, as expected.

xtx = t(x)%*%x
y1 <- t(x) %*% y 
betahat.1 <- solve(xtx,y1)

## Error in solve.default(xtx, y1): system is computationally singular: reciprocal condition number = 8.24087e-19

plot(betahat,col=2,ylim=c(0,3), pch = 1)
points(beta,col= 3, pch = 2)
legend(75,2,c("G-inv","True"),col=c(2,3),pch=c(1,2))

grid=10^seq(10,-2,length=100)
lasso.mod =glmnet(x,y,alpha =1,lambda =grid)
plot(lasso.mod)

## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to unique
## 'x' values

plot(lasso.coef[-1],col="red", pch = 1, ylim = c(-0.5,3.5))
points(beta,col= rgb(0,0,1,0.5),pch = 16)
legend(75,2,c("Lasso","True"),col=c("red","blue"),pch=c(17,16))