Lasso and Ridge Regression
An implementation in R markdown
Paul Jozefek
2020-06-23
Shrinkage Methods
1
To find the best model we can use subset selection 1 or Model contains a subset of
predictors
shrinkage methods. 2
2
Model contains all p predcitors
Subset Selection - an iterative process that
evaluates each predictor for inclusion.
Shrinkage Methods - a technique that reduces
the coefficient estimates towards zero.
Ridge Regression Basics
Least squares fits a model by minimizing the sum of
squared residuals.
2
n p
RS S = ∑ (yi − β0 − ∑ βj x ij )
i=1 j=1
3 λ ≥ 0 is a tuning parameter and the
Ridge Regression is similar, but it includes another term. 3
second term, λ ∑ p
j=1
β
2
j
, is called a
2
shrinkage penalty
n p p
2
∑ (yi − β0 − ∑ βj x ij ) + λ∑β
j
i=1 j=1 j=1
2
= RS S + λ ∑ β
j
j=1
In order to minimze this equation β , … , β should be close
1 p
to zero and so it shrinks the coefficients. The tuning
4 λ = 0 will produce the least squares
parameter, λ , controls the impact. 4
estimate. λ → ∞ will produce
coefficients near zero
Ridge regression will produce a different set of coefficient
R
estimates, β^ , for each value of λ .
λ
Importance of Scaled Data
Least squares coefficients are scale equivalent.
5 Column of x values for the j
Multiplying X 5 by a constant c will change β^ 6 by 1/c.
th
j j
predictor
Therefore, X ^
jβ j will remain the same.
6 β^ value for the j th
predictor
With ridge regression the coefficient estimates can change
7 For example, dollars in thousand
significantly when multiplying a coefficient by a constant. 7
R
vs. dollars in millions
^
will depend on λ and the scaling of the j predictor.
X j β j,λ
th
To account for this the predictors must be standardized:
x ij
~
x ij = −−−−−−−−−−−−−−
1 n
2
√ ∑ (x ij − x̄ j )
n i=1
The demoninator is the estimated standard deviation of the
j predictor.
th
All of the standardized predictors will have a
standard deviation of one.
The final model will not depend on the scale of
the predictors.
Ridge Regression vs. Least Squares
Model selection often involves a bias-variance trade-
off. 8
8
8 For more information on bias-
As λ increases the model flexibility decreases, variance see my report on Resampling
Resampling
Resampling
which will reduce the variance. However, at Methods
Methods
Methods
some point the shrikage of the coefficients will
cause them to be significantly underestimated,
9
which will increase the bias. 9 The test mean squared error
(MSE) is a function of the
When the relationship between the response and variance plus the squared bias
the predictors is close to linear least squares
10
estimates will have low bias but high variance. 10 A small change in the test
data will have a large change in
When the number of variables p is almost as the coefficient estimates
large as the number of observations n the least
squares estimates will be extremely variable.
When p > n least squares can not produce a
unique solution.
Ridge regression works best when least squares
11
estimates have high variance. 11 Trading off a small increase
in bias for a large decrease in
It also has computational advatages. Subset variance
selection evaluates 2 models, whereas for a fixed
p
λ ridge regression only fits a single model.
The Lasso
Ridge regression does have some disadvantages.
Unlike subset selection, ridge regression includes
all p predictors.
The penalty term will shrink all of the
coefficients towards zero, but none of them will
be exactly zero.
Such a large model often makes interpretation
difficult.
The lasso helps overcome these problems. It is similar to
12 |β | instead of β
ridge regression, but the penalty term is slightly different. 12
2
j
j
2
n p p
∑ (yi − β0 − ∑ βj x ij ) + λ ∑ |βj |
i=1 j=1 j=1
= RS S + λ ∑ |βj |
j=1
Like ridge regression it shrinks the coefficients towards zero.
However, the lasso allows some of the coefficients to be
exactly zero.
The Lasso vs. Ridge Regression
Models produced by the lasso are easier to
13
interpret. 13 Fewer non-zero coefficients
The lasso will perform better when a relatively
small number of predictors have substantial
coefficients.
Ridge regression will outperform when the
dependent variable is a function of many
predictors.
Ridge regression tends to shrink the coefficients
by a similar proportion while the lasso tends to
14
shrink the coefficients by a similar amount. 14 The lasso shrinks small
coefficients all the way to zero
Both the lasso and ridge regression can generate
more accurate predictions when the least squares
estimates have very high variance.
Selecting the Tuning Parameter
We can select a tuning parameter by using cross
validation.
Choose a grid of λ values and compute the cross-
validation error for each.
Select a tuning parameter value for which the
cross-validation error is smallest.
Re-fit the model using all of the observations
and the selected value of λ .
Ridge Regression Example
We can use the Hitters dataset to predict a baseball
player’s Salary using ridge regression and the lasso. It is
important to first remove missing values.
> library(ISLR) #load package
> names(Hitters) #column names
[1] "AtBat" "Hits" "HmRun" "Runs" "RBI"
[7] "Years" "CAtBat" "CHits" "CHmRun" "CRuns
[13] "CWalks" "League" "Division" "PutOuts" "Assis
[19] "Salary" "NewLeague"
> dim(Hitters) #Rows and Columns
[1] 322 20
> #Number of NA values for salary
> sum([Link](Hitters$Salary))
[1] 59
We can see that there are 59 missing values for Salary.
They can be removed with the [Link] function.
> Hitters=[Link](Hitters) #Remove NAs
> dim(Hitters) #Rows and Columns
[1] 263 20
> sum([Link](Hitters)) #Number of NA values
[1] 0
Next, we must load the glmnet package. The glmnet()
function is used with an x matrix as well as a y vector. The
[Link]() function can create x by producing a matrix
of the 19 predictors. It also transforms any qualitative
15 Glmnet() can only take numerical
variables into dummy variables. 15
inputs
> #Create matrix of 19 predictors
> x=[Link](Salary~.,Hitters)[,-1]
> #Create vector of y values
> y=Hitters$Salary
Once we have the values for x and y we can create a grid of
lambda values. Using the seq function I have created a grid
of 100 values from λ = 10 to λ = 10 . The glmnet()
10 −2
function can then be used to perform ridge regression on the
6
16 If no λ values are specified it will
selected λ values. 16 By default the variables are
automatically select a range
automatically standardized so that they are on the same
scale. This can be changed with stadardize=FALSE.
> library(glmnet) #load package
> #create grid of 100 values
> grid=10^seq(10,-2,length=100)
> #alpha = 0 for ridge regression
> #alpha = 1 for lasso
> [Link]=glmnet(x,y,alpha=0,lambda=grid)
For each value of λ there is a vector of coefficients stored in
17 One for each predictor plus the
a matrix. In this case it has 20 rows 17 and 100 columns 18.
intercept
They can be accessed by coef().
18
One for each value of λ
> #coefficient matrix size
> dim(coef([Link]))
[1] 20 100
The coefficients should be smaller when a large value of λ is
used. The 50 λ value is 11498 and the sum of squared
th
19 ∑
coefficients 19 is 6.36
p
2
β
j=1 j
> #50th lambda value
> [Link]$lambda[50]
[1] 11497.57
> #coefficients for 50th lambda
> coef([Link])[,50]
(Intercept) AtBat Hits HmRun
407.356050200 0.036957182 0.138180344 0.524629976 0
RBI Walks Years CAtBat
0.239841459 0.289618741 1.107702929 0.003131815 0
CHmRun CRuns CRBI CWalks
0.087545670 0.023379882 0.024138320 0.025015421 0
DivisionW PutOuts Assists Errors
-6.215440973 0.016482577 0.002612988 -0.020502690 0
> #sum of squared coefficients
> sqrt(sum(coef([Link])[-1,50]^2))
[1] 6.360612
The 60 th
λ value is 705 and the sum of squared coefficients
is 57.1
> #60th lambda value
> [Link]$lambda[60]
[1] 705.4802
> #coefficients for 60th lambda
> coef([Link])[,60]
(Intercept) AtBat Hits HmRun
54.32519950 0.11211115 0.65622409 1.17980910 0.937
Walks Years CAtBat CHits C
1.31987948 2.59640425 0.01083413 0.04674557 0.337
CRBI CWalks LeagueN DivisionW Pu
0.09780402 0.07189612 13.68370191 -54.65877750 0.118
Errors NewLeagueN
-0.70358655 8.61181213
> #sum of squared coefficients
> sqrt(sum(coef([Link])[-1,60]^2))
[1] 57.11001
To generate the coefficients for a specific value of λ , in this
case 50, the predict() function can be used.
> # Coefficients for lambda 50
> predict([Link],s=50,
+ type="coefficients")[1:20,]
(Intercept) AtBat Hits HmRun
4.876610e+01 -3.580999e-01 1.969359e+00 -1.278248e+00 1
RBI Walks Years CAtBat
8.038292e-01 2.716186e+00 -6.218319e+00 5.447837e-03 1
CHmRun CRuns CRBI CWalks
6.244860e-01 2.214985e-01 2.186914e-01 -1.500245e-01 4
DivisionW PutOuts Assists Errors
-1.182011e+02 2.502322e-01 1.215665e-01 -3.278600e+00 -9
In order to estimate the test error we need to split the data
into test and training sets.
> # Make replicable
> [Link](1)
> #random sample of half
> train=sample(1:nrow(x),nrow(x)/2)
> #make the values negative
> test=(-train)
> #remove training values
> [Link]=y[test]
We can use fit a model on the training data and use the
predict() function with a selected λ value to calculate the
MSE for the test set. With λ = 4 the M S E = 101037.
> #ridge regression on training data
> [Link]=glmnet(x[train,],y[train],alpha=0,
+ lambda=grid,thresh = 1e-12)
> #test set predictions for lambda = 4
> [Link]=predict([Link],s=4,newx=x[test,])
> #calculate MSE
> mean(([Link])^2)
[1] 101036.8
20 For the training set
We can use the mean y value 20 as our prediction to
calculate the MSE for an intercept only model. The same
21 Shrinks coefficients to near 0
result can be achieved by using a very large λ value. 21
> #Use mean for prediction
> mean((mean(y[train])-[Link])^2)
[1] 193253.1
> #test set predictions for lambda = 1e10
> [Link]=predict([Link],s=1e10,
+ newx=x[test,])
> #calculate MSE
> mean(([Link])^2)
[1] 193253.1
A λ value of 0 will yield the least squares results.
> #Set lambda=0 for least squares
> [Link]=predict([Link],s=0, newx=x[test,],
+ exact=TRUE,x=x[train,],
+ y=y[train])
> #MSE for least squares
> mean(([Link] -[Link])^2)
[1] 114783.1
> #Linear regression
> lm(y~x, subset=train)
Call:
lm(formula = y ~ x, subset = train)
Coefficients:
(Intercept) xAtBat xHits xHmRun x
299.42849 -2.54027 8.36682 11.64512 -9.0
xWalks xYears xCAtBat xCHits xCH
9.23440 -22.93673 -0.18154 -0.11598 -1.3
xCRBI xCWalks xLeagueN xDivisionW xPut
0.07536 -1.07841 59.76065 -98.86233 0.3
xErrors xNewLeagueN
-0.64207 -0.67442
> #Compare to lambda=0
> predict([Link],s=0,exact=TRUE,
+ type="coefficients",
+ x=x[train,],y=y[train])[1:20,]
(Intercept) AtBat Hits HmRun
299.42883596 -2.54014665 8.36611719 11.64400720 -9.098
Walks Years CAtBat CHits C
9.23403909 -22.93584442 -0.18160843 -0.11561496 -1.338
CRBI CWalks LeagueN DivisionW Pu
0.07511771 -1.07828647 59.76529059 -98.85996590 0.340
Errors NewLeagueN
-0.64205839 -0.67606314
We can see that ridge regression with λ = 4 outperformed
the intercept only model or the least squares model.
Rather than selecting our own λ value, we can use cross-
validation to find the best result. In this case it was a λ of
212.
> #Make replicable
> [Link](1)
> #10 fold cross-validation
> [Link]=[Link](x[train,], y[train], alpha=0,
+ nfolds = 10)
> #plot results
> plot([Link])
> #lambda for min MSE
> bestlam=[Link]$[Link]
> bestlam
[1] 211.7416
The best MSE was 96016, which outperformed least squares
and the intercept only model.
> #test set predictions for lambda = 212
> [Link]=predict([Link],s=bestlam,
+ newx=x[test,])
> #calculate MSE
> mean(([Link])^2)
[1] 96015.51
We can then fit the model on the full dataset and use the
optimal λ value to generate the coefficients.
> #generate full model
> out=glmnet(x,y,alpha=0)
> #coefficients for best model
> predict(out,type = "coefficients",
+ s=bestlam)[1:20,]
(Intercept) AtBat Hits HmRun
9.88487157 0.03143991 1.00882875 0.13927624 1.113
Walks Years CAtBat CHits C
1.80410229 0.13074383 0.01113978 0.06489843 0.451
CRBI CWalks LeagueN DivisionW Pu
0.13737712 0.02908572 27.18227527 -91.63411282 0.191
Errors NewLeagueN
-1.81244470 7.21208394
Lasso Example
We can perform the same functions for the lasso by
changing the function input to alpha=1. We can see from
the plot that some coefficients can be zero depending on the
λ value.
> #generate lasso model
> [Link]=glmnet(x[train,],y[train],alpha=1,
+ lambda=grid)
> #plot coefficients
> plot([Link])
Cross-validation can be used to find the best value of λ and
the minimum M S E . In this case the best M S E was
100743. It is slightly higher that the best M S E for ridge
regression, which was 96016, but it outperformed least
squares.
> #make replicable
> [Link](1)
> # 10 fold CV
> [Link]=[Link](x[train,],y[train],alpha=1,
+ nfolds=10)
> #lambda for min MSE
> bestlam=[Link]$[Link]
> #test predictions with best lambda
> [Link]=predict([Link],s=bestlam,
+ newx=x[test,])
> #MSE
> mean(([Link])^2)
[1] 100743.4
However, the lasso model is much more sparse than it was
with ridge regression. It contains only 7 variables, while
ridge regesssion contains all 19.
> #full model
> out=glmnet(x,y,alpha=1,lambda=grid)
> #predictions with best lambda
> [Link]=predict(out,type="coefficients",
+ s=bestlam)[1:20,]
> #non-zero coefficients
> [Link][[Link]!=0]
(Intercept) Hits Walks CRuns
18.5394844 1.8735390 2.2178444 0.2071252 0.41
DivisionW PutOuts
-103.4845458 0.2204284