regression using R

final

ReferenceforAssignment.pdf

Home >Computer Science homework help >regression using R

Review: training and test error

We can evaluate the prediction error of our model either on the

training dataset (known as training error) or on a dataset not used

for training the model (known as test error).

When training (fitting) a supervised learning model, the algorithm’s

goal is to minimize training error.

However, what we really want to do is to minimize the test error, or

the prediction error on data we have not observed yet.

Overfitting is an issue, because we can begin modeling noise in the

training data rather than only the underlying function.

This can lead to low training error, but high test error.

Review: illustration of training vs. test error

Below is prediction error as a function of polynomial degree 𝑀 for both training and test datasets:

Training error always decreases with increasing model complexity

(flexibility), whereas test error has a U-shaped curve.

Bishop (2006) Pattern Recognition and Machine Learning. Springer.

Prediction

error

How do we

identify 𝑴 such that the inferred

model has good

test error?

How do we minimize test error?

We have two main goals:

1) Model selection – comparing the performance of different

models to select the best one

2) Model assessment – estimating test error of the selected model

If we have a lot of data, then the optimal approach for both of these is

to divide the dataset uniformly at random into three parts:

Training set

Validation set

Test set

Training, validation, and test sets

Training set: used for fitting the models

Validation set: used to estimate prediction error for model selection

Test set: used for assessing prediction (test) error on new data

The test set should only be used at the end of the analysis for model

assessment, and not for choosing the best model.

There is no general rule for how to split the dataset, but a typical split

is 50% for training, 25% for validation, and 25% for test.

Hastie et al (2009) The Elements of Statistical Learning: Data Mining, Inference, and Prediction. Springer, 2nd Ed.

Do we need to split data into three parts?

For many datasets it may not be possible.

That is, there may not be enough data to create a separate test

dataset without decreasing the sample size such that it is deleterious

to model selection.

Model selection is critical, and often our primary concern is training a

method that will perform as well as possible on unseen data.

We will therefore consider methods that split the dataset into only

training and validation sets.

In these methods, the validation set is used for estimating test error.

The validation (hold-out) set approach

The validation (or hold-out) set approach involves dividing the training

dataset uniformly at random into roughly two equally-sized parts: the

training set and the validation set (known as the hold-out set).

The model is rained on the traning se

The model is trained on the training set (blue), and the test error of

the trained model is estimated from the validation set (orange).

Models of different complexity can be trained and their test error

evaluated, and then the best-performing model can be selected.

James et al (2017) An Introduction to Statistical Learning: with Applications in R. Springer.

Validation set approach on polynomial regression

Validation (test) error for a single split into training and validation sets:

James et al (2017) An Introduction to Statistical Learning: with Applications in R. Springer.

Validation set approach on polynomial regression

Validation (test) error for different splits of the data:

Test error is highly variable across splits, and depends on which

observations are included in the training vs. validation sets.

This variation occurs because the validation set approach trains the

model on only half of the data (small sample).

Methods tend to perform worse when trained on fewer observations,

suggesting validation error overestimates test error on entire dataset. James et al (2017) An Introduction to Statistical Learning: with Applications in R. Springer.

Overcoming pitfalls of the validation set approach

We need the learned model to be more robust to variability in the

choice of training and validation set.

We also want the trained model to be based on a large sample size.

The solution to this problem is cross-validation (CV), which involves

repeatedly splitting the data into training and test sets, and then

assessing mean performance across iterations.

There are tow major CV approaches:

1) Leave-one-out cross-validation (LOOCV)

2) 𝑘-fold cross validation

Leave-one-out cross-validation (LOOCV)

James et al (2017) An Introduction to Statistical Learning: with Applications in R. Springer.

Leave-one-out cross-validation (𝐋𝐎𝐎𝐂𝐕) splits the data into a training set of 𝑛 − 1 observations and a validation set of the remaining 1 observation:

Leave-one-out cross-validation (LOOCV)

We start by removing training observation (𝑥1,𝑦1), and using only observations {(𝑥2,𝑦2),(𝑥3,𝑦3),…,(𝑥𝑛,𝑦𝑛)} for training.

The model is fit to the training data, and a prediction is made for

training observation (𝑥1,𝑦1).

This validation (test) error of training observation (𝑥1,𝑦1) is:

MSE1 = 𝑦1 − ො𝑦1 2

This estimate of the MSE of the method will be highly variable

because it is based on only a single observation.

Leave-one-out cross-validation (LOOCV)

We repeat this process for each of the 𝑛 training observations, such that we remove training observation (𝑥𝑖,𝑦𝑖), and use observations { 𝑥1,𝑦1 ,…, 𝑥𝑖−1,𝑦𝑖−1 , (𝑥𝑖+1,𝑦𝑖+1),…,(𝑥𝑛,𝑦𝑛)} for training.

The model is fit to the training data, and a prediction is made for

training observation (𝑥𝑖,𝑦𝑖).

The validation (test) error of training observation (𝑥𝑖,𝑦𝑖) is:

MSE𝑖 = 𝑦𝑖 − ො𝑦𝑖 2

When we have computed MSE for all 𝑛 training observations, we compute the LOOCV validation (test) error as:

CV(𝑛) = MSE1 + MSE2 + ⋯+ MSE𝑛

𝑛

Advantages of LOOCV over validation set

LOOCV has substantially less bias because the model is repeatedly

trained on a large sample set of size 𝑛 − 1 training observations, whereas the validation set approach trains models with a sample size

of approximately 𝑛/2 training observations.

Because of this, LOOCV tends to not overestimate the test error.

Moreover, LOOCV always yields the same test error and inferred

model, compared to variation attributed to different splits of the data

with the validation set approach.

However, LOOCV is computationally expensive, because the model

needs to be trained 𝑛 times.

LOOCV approach on polynomial regression

James et al (2017) An Introduction to Statistical Learning: with Applications in R. Springer.

𝑘-fold CV

James et al (2017) An Introduction to Statistical Learning: with Applications in R. Springer.

𝒌-fold cross-validation (𝒌−fold 𝐂𝐕) randomly splits data into 𝑘 “folds”, with a training set of 𝑘 − 1 folds and a test set of the remaining 1 fold.

𝑘-fold cross-validation (CV)

This process is repeated for each set (just like with LOOCV):

When set 𝑖 is held out for validation, the remaining 𝑘 − 1 sets (sets 1,…,𝑖 − 1,𝑖 + 1,…,𝑘) are used for training.

The MSE across the observations in the 𝑖th set is used to measure the validation (test) error, and is denoted by MSE𝑖, and the overall 𝑘- fold CV test error is computed as

CV(𝑘) = MSE1 + MSE2 + ⋯+ MSE𝑘

𝑘

Computational savings with 𝑘-fold CV

LOOCV is a special case of 𝑘-fold CV with 𝑘 = 𝑛.

Though LOOCV is superior to 𝑘-fold CV (𝑘 < 𝑛), it is also computationally expensive to train the model 𝑛 times.

In practice, we usually use 𝑘 = 5 or 𝑘 = 10 for 𝑘-fold CV.

This makes 𝑘-fold CV (𝑘 < 𝑛) much less computationally intensive than LOOCV.

In such cases, we only have to train 5 or 10 models, rather than 𝑛, which could lead to enormous savings for methods that are

computationally intensive to fit.

LOOCV vs. 10-fold CV on polynomial regression

James et al (2017) An Introduction to Statistical Learning: with Applications in R. Springer.

The overall performance is similar between LOOCV and 10-fold CV:

Though there is some variability across the 10-fold CV estimates, it is much lower than the variability from the validation set approach.

Non-computational advantages of 𝑘-fold CV (𝑘 < 𝑛)

Because LOOCV uses a greater number of observations in each

training set, its model estimates will be based on a larger sample

size, and should exhibit smaller bias relative to 𝑘-fold CV (𝑘 < 𝑛).

However, LOOCV has a higher variance than 𝑘-fold CV (𝑘 < 𝑛).

LOOCV is essentially averaging the output of 𝑛 fitted models, each of which is highly positively correlated with the others (as almost all

training observations overlap).

𝑘-fold CV is averaging the output of 𝑘 fitted models, each of which is partially correlated with the others.

Non-computational advantages of 𝑘-fold CV (𝑘 < 𝑛)

The mean of many highly correlated quantities has higher variance

than the mean of many quantities that are less correlated.

Because of this, the test error estimates resulting from 𝑘-fold CV (𝑘 < 𝑁) tend to be smaller than those from LOOCV.

Therefore, there is a bias-variance trade-off depending on the chosen

value of 𝑘 in 𝑘-fold CV.

For these reasons, 5- and 10-fold CV is often performed, as they

have been demonstrated empirically to provided estimates of the test

error that do not experience high bias or high variance.

Feature selection

A related problem to model selection is feature selection, which

involves determining which features are associated with the response

so that we can fit a model with only those features.

In feature selection, we seek to minimize test error, while also

minimizing the number of features used to make predictions.

After feature selection, chosen features would hopefully be the

important ones for making predictions about the response.

Irrelevant features not associated with the response are discarded.

Redundant (correlated) features are also discarded.

Why perform feature selection?

Feature selection can enhance interpretability of the fitted model

because fewer features are associated with the response.

The reduction of model parameters may lead to better generalization

or test error due to a reduction in overfitting.

Including fewer features can reduce the time it takes to train models,

as training time increases with increased dimensionality (more

complexity/ features).

Subset selection methods

We begin by considering subset selection methods.

These methods try to identify a subset of the 𝑝 features that will lead to reduced prediction error.

We will consider three approaches:

1) Best subset selection – compare all possible subsets of features

2) Forward selection – start with null model (no features) and add

3) Backward selection – start with full model (all features) and subtract

Best subset selection

1. Let ℳ0 denote the null model that contains no features, which is the model with only parameter 𝛽0.

2. For each 𝑘 = 1,2,…,𝑝:

(a) Fit all models of exactly 𝑘 features.

(b) Pick the best (smallest RSS) among these models, and call it

ℳ𝑘.

3. Select a single best model from {ℳ0,ℳ1,…ℳ𝑝} based on model

selection (e.g., cross-validation).

Adapted from James et al (2017) An Introduction to Statistical Learning: with Applications in R. Springer.

This is the ideal approach because we are comparing every

possible subset of features!

Illustration of best subset selection

RSS for dataset with 𝑝 = 8 features:

Choose model with lowest RSS among all sets of models.

The best model contains all 𝑝 = 8 features in the training set.

Hastie et al (2009) The Elements of Statistical Learning: Data Mining, Inference, and Prediction. Springer, 2nd Ed.

Problem with best subset selection

In best subset selection, we must perform least squares regression

for all subsets with 0,1,2,…,𝑝 features.

The set of all subsets is known as the power set, which is

exponentially increasing in size of the original set.

Because we have 𝑝 features, the total number of fitted models is 2𝑝, which is enormous, and computationally infeasible for anything but a

small number of features (maximum 𝑝 between 30 and 40 based on an efficient algorithm of Furnival and Wilson 1974).

For example, even when 𝑝 = 20 (small number of features for modern datasets), we must consider 220 = 1,048,576 models!

Furnival, Wilson (1974) Technometrics 16:499-511.

How can we overcome this computational burden?

Best subset selection is simply not computationally feasible for many

contemporary datasets.

The other two approaches that we will discuss examine a far more

restricted space compared to best subset selection:

2) Forward selection – start with null model (no features) and add

3) Backward selection – start with full model (all features) and subtract

Forward selection

1. Let ℳ0 denote the null model that contains no features, which is the model with only parameter 𝛽0.

2. For each 𝑘 = 0,1,…,𝑝 − 1:

(a) Fit all 𝑝 − 𝑘 models that augment the features in ℳ𝑘 with one additional feature.

(b) Pick the best (smallest RSS) among these 𝑝 − 𝑘 models, and call it ℳ𝑘+1.

3. Select a single best model from {ℳ0,ℳ1,…ℳ𝑝} based on model

selection (e.g., cross validation).

Adapted from James et al (2017) An Introduction to Statistical Learning: with Applications in R. Springer.

Backward selection

1. Let ℳ𝑝 denote the full model that contains all 𝑝 features.

2. For each 𝑘 = 𝑝,𝑝 − 1,…,1:

(a) Fit all 𝑘 models that have all but one of the features in ℳ𝑘, for a total of 𝑘 − 1 features.

(b) Pick the best (smallest RSS) among these 𝑘 models, and call it ℳ𝑘−1.

3. Select a single best model from {ℳ0,ℳ1,…ℳ𝑝} based on model

selection (e.g., cross validation).

Adapted from James et al (2017) An Introduction to Statistical Learning: with Applications in R. Springer.

Number of fitted models with forward & backward selection

The number of fitted models for either forward selection or backward

selection is:

1 + 𝑝 𝑝 + 1

Thus, they are equivalent in computational efficiency.

This number is also much smaller than the 2𝑝 models required by best subset selection.

For example, with 𝑝 = 20, best subset selection requires 1,048,576 fitted models, while forward and backward selection require only 211.

How to choose between forward & backward selection?

First, it is important to point out that both methods may not yield the

optimal solution, which can only be found by best subset selection.

The optimal solutions for forward and backward stepwise selection

also may differ.

Backward subset selection requires 𝑝 ≤ 𝑛, so that the full model can be fit without issue.

Forward subset selection can be used even when 𝑝 > 𝑛, as the algorithm can just stop at the subset of size 𝑘 = 𝑛 features.

However, forward selection is “greedy” and might include features

early that later become less important.

Can we get closer to best subset selection?

A hybrid approach between forward and backward stepwise selection

can be taken to get a better approximation to best subset selection.

With such an approach, features are added iteratively through

forward stepwise selection.

After each new feature is added, the method can remove any feature

from the model that no longer improves model fit.

Therefore, features added earlier in the process may be removed if

they are no longer important conditional on the set of added features

that were added later in the process.

This alleviates the “greedy” issue of forward selection.

Comparison of subset selection approaches

𝑁 = 30 observations with 𝑝 = 31 standard normal features with pairwise correlations of 0.85:

Adapted from Hastie et al (2009) The Elements of Statistical Learning: Data Mining, Inference, and Prediction. Springer, 2nd Ed.

Shrinkage methods as alternative to subset selection

Subset selection methods retain a subset of the features and discard

the remainder, often leading to a more interpretable model.

However, this discrete process of retaining or discarding features

often leads to high variance, and therefore does not reduce the

prediction error as much as hoped from the full model.

An alternative approach is to fit a model with all features, then shrink

estimated parameters toward zero relative to least squares estimates.

This can reduce variance of estimates, as well as potentially remove

(estimate parameters to be zero) unimportant features.

This approach is called regularization because it constrains, or

regularizes, parameter estimates.

Regularization techniques

Recall that in least squares, we estimate parameters 𝛽0 …𝛽𝑝 by

minimizing the cost function 𝐽 𝛽 :

𝐽 𝛽 = RSS = ෍

𝑖=1

𝑛

(𝑦𝑖 − ො𝑦𝑖) 2

Regularization shrinks estimates by minimizing a penalized variation.

We will discuss three common techniques:

1) Ridge regression

2) Lasso

3) Elastic net

Regularization techniques

Recall that in least squares, we estimate parameters 𝛽0 …𝛽𝑝 by

minimizing the cost function 𝐽 𝛽 :

𝐽 𝛽 = RSS = ෍

𝑖=1

𝑛

(𝑦𝑖 − ො𝑦𝑖) 2

Regularization shrinks estimates by minimizing a penalized variation.

We will discuss three common techniques:

1) Ridge regression

2) Lasso

3) Elastic net

Ridge regression

Recall that in least squares, we estimate parameters 𝛽0 …𝛽𝑝 by

minimizing the cost function 𝐽 𝛽 :

𝐽 𝛽 = RSS = ෍

𝑖=1

𝑛

(𝑦𝑖 − ො𝑦𝑖) 2

Ridge regression instead minimizes:

𝐽 𝛽,𝜆 = RSS + 𝜆 ෍

𝑗=1

𝑝

𝛽𝑗 2 ,

where the second term is a shrinkage (L2) penalty, and 𝜆 ≥ 0 is a tuning parameter that determines its influence on minimization.

How ridge regression works

Ridge regression minimizes:

𝐽 𝛽,𝜆 = RSS + 𝜆 ෍

𝑗=1

𝑝

𝛽𝑗 2 ,

where the second term is a shrinkage (L2) penalty, and 𝜆 ≥ 0 is a tuning parameter that determines its influence on minimization.

The shrinkage penalty is smallest when 𝛽1 …𝛽𝑝 are close to zero, thus

shrinking their estimates close to zero.

The shrinkage penalty has no influence when 𝜆 = 0 (least squares regression), and increasing influence as 𝜆 → ∞.

Choosing the tuning parameter 𝜆.

The original formulation of least squares regression minimizes the cost

function 𝐽(𝛽) = RSS, leading to a single set of parameter estimates.

However, ridge regression leads to different sets of estimates across

values of the tuning parameter 𝜆.

We choose the best set of parameter estimates under ridge regression

through cross-validation.

That is, we identify መ𝜆 through cross-validation that leads to the best መ𝛽.

We then estimate parameters based on this tuning parameter መ𝜆.

Example behavior of parameter estimates with increasing 𝜆

Parameter estimates (coefficients) for predicting credit card balance

from income, credit limit, credit rating, and student status as a

function of 𝜆 :

James et al (2017) An Introduction to Statistical Learning: with Applications in R. Springer.

Tuning parameter creates bias-variance tradeoff

The tuning parameter 𝜆 modifies model complexity.

That is, as the tuning parameter increases, model complexity

decreases, and we have fewer effective parameters.

This leads to a greater bias but a decreased variance.

The optimal 𝜆 value is chosen such that there is an optimal tradeoff between bias and variance, yielding improved test error.

Choosing the optimal value of the tuning parameter 𝜆

Depicted are MSEs (purple), squared bias (black), variance (green),

and irreducible error (dotted black) for different values of 𝜆:

James et al (2017) An Introduction to Statistical Learning: with Applications in R. Springer.

Optimal MSE (purple X)

occurs at intermediate

value of 𝜆

Advantages of ridge regression over other approaches

Least squares estimates have increasing variance as 𝑝 → 𝑛, and do not even have unique solutions when 𝑝 > 𝑛.

Ridge regression can still perform well when 𝑝 > 𝑛, trading small increases in bias for potentially larger decreases in variance.

Moreover, we must fit 2𝑝 models for best subset selection, and 1 + 𝑝(𝑝 + 1)/2 models for forward and backward selection.

In contrast, ridge regression only requires us to fit a single model for

any value of 𝜆, which can be done quickly.

But there is one obvious disadvantage to ridge regression relative to

all subset selection methods—it includes all 𝑝 features in the model.

Regularization techniques

Recall that in least squares, we estimate parameters 𝛽0 …𝛽𝑝 by

minimizing the cost function 𝐽 𝛽 :

𝐽 𝛽 = RSS = ෍

𝑖=1

𝑛

(𝑦𝑖 − ො𝑦𝑖) 2

Regularization shrinks estimates by minimizing a penalized variation.

We will discuss three common techniques:

1) Ridge regression

2) Lasso

3) Elastic net

Lasso is very similar in form to ridge regression

Ridge regression minimizes:

𝐽 𝛽,𝜆 = RSS + 𝜆 ෍

𝑗=1

𝑝

𝛽𝑗 2 ,

where the second term is a shrinkage (L2) penalty, and 𝜆 ≥ 0 is a tuning parameter that determines its influence on minimization.

Lasso instead minimizes:

𝐽 𝛽,𝜆 = RSS + 𝜆 ෍

𝑗=1

𝑝

𝛽𝑗 ,

where the second term is a shrinkage (L1) penalty, and 𝜆 ≥ 0 is a tuning parameter that determines its influence on minimization.

Lasso also works similarly to ridge regression

Lasso minimizes:

𝐽 𝛽,𝜆 = RSS + 𝜆 ෍

𝑗=1

𝑝

𝛽𝑗 ,

where the second term is a shrinkage (L1) penalty, and 𝜆 ≥ 0 is a tuning parameter that is selected by cross-validation.

As with ridge regression, the shrinkage penalty is smallest when 𝛽1 …𝛽𝑝 are close to zero, and has increasing influence as 𝜆 → ∞.

But unlike ridge regression, lasso can perform feature selection by

forcing some estimates to be exactly zero when is 𝜆 sufficiently large.

Tuning parameter creates bias-variance tradeoff

The tuning parameter 𝜆 modifies model complexity.

That is, as the tuning parameter increases, model complexity

decreases, and we have fewer effective parameters.

This leads to a greater bias but a decreased variance.

In addition, the decreased variance is associated with some features

becoming discarded as the tuning parameter increases in value.

The optimal 𝜆 value is chosen such that there is an optimal tradeoff between bias and variance, yielding improved test error.

Depicted are MSEs (purple), squared bias (black), and variance

(green) for different values of 𝜆:

Tuning parameter makes bias-variance tradeoff

James et al (2017) An Introduction to Statistical Learning: with Applications in R. Springer.

Optimal MSE (purple X)

again occurs at

intermediate value of 𝜆

Ridge regression vs. lasso: Which is better?

Some of the same advantages of ridge regression also apply to lasso

(perform well when 𝑝 > 𝑛 unlike least squares, only need to fit a single model for each 𝜆 unlike subset selection methods).

Though ridge regression and lasso may outcompete these other

methods, neither one is always better than the other.

Ridge regression will generally perform better when the response is a

function of many features (because it must include all features).

Lasso will generally perform better when the response is a function of

few features (because it can perform feature selection).

However, we do not typically know how many features are related to the

response a priori for real data.

Comparison of approaches for prostate cancer prediction

Parameter estimates and test errors for least squares (LS), best

subset selection, ridge regression, and lasso:

Adapted from Hastie et al (2009) The Elements of Statistical Learning: Data Mining, Inference, and Prediction. Springer, 2nd Ed.

Comparison of approaches for prostate cancer prediction

Parameter estimates and test errors for least squares (LS), best

subset selection, ridge regression, and lasso:

Which approach is best? What does this say about the data?

Adapted from Hastie et al (2009) The Elements of Statistical Learning: Data Mining, Inference, and Prediction. Springer, 2nd Ed.

What if the best solution lies between ridge and lasso?

The best solution may not be the one determined by ridge regression,

or by lasso, but rather may be somewhere in between the two.

In such cases, combining ridge regression and lasso may increase

predictive ability.

This solution, termed elastic net, was proposed such that both ridge

regression and lasso are nested models within it.

Regularization techniques

Recall that in least squares, we estimate parameters 𝛽0 …𝛽𝑝 by

minimizing the cost function 𝐽 𝛽 :

𝐽 𝛽 = RSS = ෍

𝑖=1

𝑛

(𝑦𝑖 − ො𝑦𝑖) 2

Regularization shrinks estimates by minimizing a penalized variation.

We will discuss three common techniques:

1) Ridge regression

2) Lasso

3) Elastic net

Elastic net combines ridge regression and lasso

Ridge regression minimizes:

𝐽 𝛽,𝜆 = RSS + 𝜆 ෍

𝑗=1

𝑝

𝛽𝑗 2

Lasso minimizes:

𝐽 𝛽,𝜆 = RSS + 𝜆 ෍

𝑗=1

𝑝

𝛽𝑗

Elastic net minimizes:

𝐽 𝛽,𝜆,𝛼 = RSS + 𝜆 1 − 𝛼 ෍

𝑗=1

𝑝

𝛽𝑗 2 + 𝛼 ෍

𝑗=1

𝑝

𝛽𝑗 ,

where the second term is a shrinkage penalty, 𝜆 ≥ 0 is a tuning parameter that determines its influence on minimization, and 0 ≤ 𝛼 ≤ 1 is another tuning parameter that determines the influence of ridge

regression and lasso on the estimated parameters.

Elastic net works similarly to ridge regression and lasso

Elastic net minimizes:

𝐽 𝛽,𝜆,𝛼 = RSS + 𝜆 1 − 𝛼 ෍

𝑗=1

𝑝

𝛽𝑗 2 + 𝛼 ෍

𝑗=1

𝑝

𝛽𝑗 ,

where the second term is a shrinkage penalty, and 𝜆 ≥ 0 and 0 ≤ 𝛼 ≤ 1 are tuning parameters that are selected by cross-validation.

As with ridge regression and lasso, the shrinkage penalty is smallest

when 𝛽1 …𝛽𝑝 are close to zero, and has increasing influence as 𝜆 → ∞.

The tuning parameter 𝛼 determines the influence of ridge regression and lasso on the estimated parameters.

We obtain ridge regression when 𝛼 = 0, and lasso when 𝛼 = 1.

Comparison of approaches for prostate cancer prediction

Parameter predictions and test errors for ordinary least squares

(OLS), ridge regression, lasso, naïve elastic net, and elastic net

(improved performance, citation below):

Note that 𝑠 is the fraction of the L1 norm 𝛽𝑗 , which is in the L1

penalty term in lasso.

Zou and Hastie (2005) J. R. Statist. Soc. B. 67:301-320.

Cross-validated ridge and lasso regression in R

R of course has the machinery for us to easily perform ridge and

lasso regression together with cross-validation to determine the

appropriate tuning parameter 𝜆 (i.e., model complexity).

We first need to install and load the library for the glmnet package,

which provides an implementation of elastic net (and therefore lasso

and ridge) regression together with cross-validation.

To install and load glmnet, use the command:

> install.packages("glmnet")

> library(glmnet)

We will apply these methods to the Advertising dataset from our

lecture on linear regression.

Review: Advertising dataset

We will next consider the Advertising dataset, which can be

downloaded from https://www.statlearning.com/resources-first-edition.

This dataset is tidied into a data frame with 𝑛 = 200 market observations (rows) on 𝑝 = 4 features (columns):

sales Sales (thousands of units)

TV TV advertising budget (thousands of dollars)

radio Radio advertising budget (thousands of dollars)

newspaper Newspaper advertising budget (thousands of dollars)

Our goal is to predict sales from the other three features.

https://www.statlearning.com/resources-first-edition

Creating the response and input feature matrices

The glmnet package requires that response and feature datasets be

provided separately and as matrices instead of data frames.

We can easily convert a data frame to a matrix with the as.matrix() function in R.

Using this function, we can obtain the response (𝑌) and feature (𝑋) matrices for use in the glmnet package:

> Y <- Advertising %>% select(sales) %>% as.matrix()

> X <- Advertising %>% select(TV, radio, newspaper) %>%

as.matrix()

The response and input feature matrices

We can use the head() function to peek at the first few rows of the

response and feature matrices:

> head(Y)

sales

[1,] 22.1

[2,] 10.4

[3,] 9.3

[4,] 18.5

[5,] 12.9

[6,] 7.2

> head(X)

TV radio newspaper

[1,] 230.1 37.8 69.2

[2,] 44.5 39.3 45.1

[3,] 17.2 45.9 69.3

[4,] 151.5 41.3 58.5

[5,] 180.8 10.8 58.4

[6,] 8.7 48.9 75.0

Specifying the set of tuning parameter (𝜆) values

Our goal is to identify the “best” 𝜆 and corresponding model through cross-validation.

Thus, we need to ensure that we explore a range of 𝜆 values that is wide enough to capture all important levels of model complexity.

We will consider 100 𝜆 values that span a range from 0.001 to 1000 on a logarithmic scale.

To do this, we will use the seq() function together with some

transformations:

> lambdas <- 10^seq(-3, 3, length.out = 100)

Obtaining ridge parameter estimates as a function of 𝜆

We can examine the set of all ridge regression parameter estimates as a function of the 𝜆 values in lambdas.

Specifically, we want to use ridge regression (𝛼 = 0) to estimate the four model parameters corresponding to the intercept and coefficients for TV, radio, and newspaper for each of the 100 𝜆 values.

To do this, we use the glmnet() function and set the argument

alpha (𝛼) to zero and provide the argument lambda (𝜆) our set of tuning parameter values lambdas:

> ridge.fit <- glmnet(X, Y, alpha = 0, lambda = lambdas)

Plotting parameter estimates as a function of 𝜆

The package glmnet comes with its own useful plotting functions.

We can view how parameter estimates shrink as 𝜆 increases:

> plot(ridge.fit, xvar = "lambda")

> legend("topright", lwd=1, col = 1:3, legend = colnames(X))

Using 𝑘-fold CV to fit a ridge regression model

We next wish to perform 𝑘 = 10-fold cross-validation as a function of the set of 𝜆 values in lambdas.

Specifically, we want to use ridge regression (𝛼 = 0) to calculate the 10-fold cross-validation error rate for each of the 100 𝜆 values.

To do this, we use the cv.glmnet() function and set the argument

alpha (𝛼) to zero, provide the argument lambda (𝜆) our set of tuning parameter values lambdas, and set the nfolds (𝑘) argument to 10:

> ridge.cv <- cv.glmnet(X, Y, alpha = 0, lambda = lambdas,

nfolds = 10)

Plotting the cross-validation error rate as a function of 𝜆

We can view the cross-validation error rate for ridge regression as a

function of 𝜆 with the command:

> plot(ridge.cv)

Numbers at the top show the number of selected features, which take

values of 0,1,…,𝑝 (ridge regression always retains all 𝑝 features).

Estimating the best model (smallest CV error rate)

We now would like to find the model with the smallest CV error rate.

To do so, we must first select the 𝜆 with the minimum CV error rate across the range of values that we explored.

Given our cross-validated object ridge.cv, we can access the 𝜆

with the minimum CV error rate using the command ridge.cv$lambda.min.

We can then find the best-fit ridge regression model by typing:

> ridge.best <- glmnet(X, Y, alpha = 0,

lambda = ridge.cv$lambda.min)

Obtaining parameter estimates

We can view the parameter estimates using the coef() function:

> coef(ridge.best)

4 x 1 sparse Matrix of class "dgCMatrix"

(Intercept) 2.9906421059

TV 0.0455213515

radio 0.1874124145

newspaper -0.0007094929

This is almost identical to the least squares regression estimates:

> lm.fit

Call:

lm(formula = sales ~ TV + radio + newspaper, data = Advertising)

Coefficients:

(Intercept) TV radio newspaper

2.938889 0.045765 0.188530 -0.001037

Making predictions with the fitted model

Given our fitted model ridge.best, we can now make predictions

from new data with the predict() function in glmnet.

We can make predictions from the best-fit ridge regression model as:

> predict(ridge.best, as.matrix(tibble(TV = 1000,

radio = 2000, newspaper = 3000)))

[1,] 421.2083

This is almost identical to our prediction with least squares:

> predict(lm.fit, tibble(TV = 1000, radio = 2000,

newspaper = 3000))

422.6511

Obtaining lasso parameter estimates as a function of 𝜆

The procedure for performing lasso regression in glmnet is akin to

ridge, except we set the argument alpha (𝛼) to one instead of zero:

> lasso.fit <- glmnet(X, Y, alpha = 1, lambda = lambdas)

> plot(lasso.fit, xvar = "lambda")

> legend("topright", lwd=1, col = 1:3, legend = colnames(X))

Using 𝑘-fold CV to fit a lasso regression model

We next wish to perform 𝑘 = 10-fold cross-validation as a function of the set of 𝜆 values in lambdas.

Specifically, we want to use lasso (𝛼 = 1) to calculate the 10-fold cross-validation error rate for each of the 100 𝜆 values.

To do this, we use the cv.glmnet() function and set the argument

alpha (𝛼) to one, provide the argument lambda (𝜆) our set of tuning parameter values lambdas, and set the nfolds (𝑘) argument to 10:

> lasso.cv <- cv.glmnet(X, Y, alpha = 1, lambda = lambdas,

nfolds = 10)

Plotting the cross-validation error rate as a function of 𝜆

We can view the cross-validation error rate of lasso as a function of 𝜆 with the command:

> plot(lasso.cv)

Recall that numbers at the top show the number of selected features,

which take values of 0,1,…,𝑝.

Estimating the best model (smallest CV error rate)

We now would like to find the model with the smallest CV error rate.

To do so, we must first select the 𝜆 with the minimum CV error rate across the range of values that we explored.

Given our cross-validated object ridge.cv, we can access the 𝜆

with the minimum CV error rate using the command ridge.cv$lambda.min.

We can then find the best-fit lasso regression model by typing:

> lasso.best <- glmnet(X, Y, alpha = 1,

lambda = lasso.cv$lambda.min)

Obtaining parameter estimates

We can view the parameter estimates using the coef() function:

> coef(lasso.best)

4 x 1 sparse Matrix of class "dgCMatrix"

(Intercept) 3.07623499

TV 0.04520385

radio 0.18480819

newspaper .

These are similar to ridge regression with newspaper removed:

> coef(ridge.best)

4 x 1 sparse Matrix of class "dgCMatrix"

(Intercept) 2.9906421059

TV 0.0455213515

radio 0.1874124145

newspaper -0.0007094929

Making predictions with the fitted model

Given our fitted model lasso.best, we can now make predictions

from new data with the predict() function in glmnet.

We can make predictions from the best-fit lasso regression model as:

> predict(lasso.best, as.matrix(tibble(TV = 1000,

radio = 2000, newspaper = 3000)))

[1,] 417.8965

This is similar to our prediction from the best-fit ridge regression model:

> predict(ridge.best, as.matrix(tibble(TV = 1000,

radio = 2000, newspaper = 3000)))

[1,] 421.2083

Application of regularization to classification problems

We can also apply penalization to logistic regression!

The ridge-penalized logistic regression cost function is:

𝐽 𝛽,𝜆 = −log𝐿(𝛽) + 𝜆 ෍

𝑗=1

𝑝

𝛽𝑗 2

The lasso-penalized logistic regression cost function is:

𝐽 𝛽,𝜆 = −log𝐿(𝛽) + 𝜆 ෍

𝑗=1

𝑝

𝛽𝑗

The elastic net-penalized logistic regression cost function is:

𝐽 𝛽,𝜆,𝛼 = −log𝐿(𝛽) + 𝜆 1 − 𝛼 ෍

𝑗=1

𝑝

𝛽𝑗 2 + 𝛼 ෍

𝑗=1

𝑝

𝛽𝑗

Cross-validated ridge and lasso logistic regression in R

Recall that R has the machinery for us to easily perform ridge and

lasso logistic regression together with cross-validation to determine

the appropriate tuning parameter 𝜆 (i.e., model complexity) using the elastic net model in the glmnet package.

To load glmnet, use the command:

> library(glmnet)

We will apply these methods to the Default dataset from our lecture

on classification with logistic regression.

Review: the Default dataset

The Default dataset is located in the ISLR package.

Simulated observations from 10,000 customers on four features:

default (Yes = defaulted, No = did not default)

student (Yes = student, No = not student)

balance (average balance after monthly payment)

income (annual income)

Our goal is to predict default from the other three features.

Cross-validated ridge and lasso logistic regression in R

Recall that to perform logistic regression, we coded the two default

classes (Yes and No) as 1 and 0:

> Default.recoded <- Default %>%

mutate(default = ifelse(default == "Yes", 1, 0),

student = ifelse(student == "Yes", 1, 0))

We will use these recoded data again for consistency, so that we can compare our results here to those obtained earlier.

Creating the response and input feature matrices

Recall that glmnet requires that response and feature datasets be

provided separately and as matrices instead of data frames.

We can easily convert a data frame to a matrix with the as.matrix() function in R.

Using this function, we can obtain the response (𝑌) and feature (𝑋) matrices for use in the glmnet package:

> Y <- Default.recoded %>% select(default) %>% as.matrix()

> X <- Default.recoded %>%

select(balance, income, student) %>%

as.matrix()

The response and input feature matrices

We can use the head() function to peek at the first few rows of the

response and feature matrices:

> head(Y)

default

[1,] 0

[2,] 0

[3,] 0

[4,] 0

[5,] 0

[6,] 0

> head(X)

balance income student

[1,] 729.5265 44361.625 0

[2,] 817.1804 12106.135 1

[3,] 1073.5492 31767.139 0

[4,] 529.2506 35704.494 0

[5,] 785.6559 38463.496 0

[6,] 919.5885 7491.559 1

Specifying the set of tuning parameter (𝜆) values

Recall that our goal is to identify the “best” 𝜆 and corresponding model through cross-validation.

Thus, we need to ensure that we explore a range of 𝜆 values that is wide enough to capture all important levels of model complexity.

We will consider 100 lambda values that span a range from 10−6 to 106 on a logarithmic scale.

To do this, we will use the seq() function together with some

transformations:

> lambdas <- 10^seq(-6, 6, length.out = 100)

Obtaining ridge parameter estimates as a function of 𝜆

We can examine the set of all ridge regression parameter estimates as a function of the 𝜆 values in lambdas.

Specifically, we want to use ridge regression (𝛼 = 0) to estimate the four model parameters corresponding to the intercept and coefficients for balance, income, and student for each of the 100 𝜆 values.

To do this, we use the glmnet() function and set the argument

family to "binomial" to indicate logistic regression, the argument

alpha (𝛼) to zero for ridge regression, and provide the argument lambda (𝜆) our set of tuning parameter values lambdas:

> ridge.fit <- glmnet(X, Y, family = "binomial", alpha = 0,

lambda = lambdas)

Plotting parameter estimates as a function of 𝜆

Recall that glmnet comes with its own useful plotting functions.

Let us view how parameter estimates shrink as 𝜆 increases:

> plot(ridge.fit, xvar = "lambda")

> legend("bottomright", lwd=1, col=1:3, legend = colnames(X))

Using 𝑘-fold CV to fit a ridge logistic regression model

We next wish to perform 𝑘 = 10-fold cross-validation as a function of the set of 𝜆 values in lambdas.

Specifically, we want to use ridge regression (𝛼 = 0) to calculate the 10-fold cross-validation error rate for each of the 100 𝜆 values.

To do this, we use the cv.glmnet() function and set the argument

family to "binomial" to indicate logistic regression, the argument

alpha (𝛼) to zero for ridge regression, provide the argument lambda (𝜆) our set of tuning parameter values lambdas, and set the nfolds

(𝑘) argument equal to 10:

> ridge.cv <- cv.glmnet(X, Y, family = "binomial", alpha = 0,

lambda = lambdas, nfolds = 10)

Plotting the cross-validation error rate as a function of 𝜆

We can view the cross-validation error rate for ridge logistic

regression as a function of 𝜆 with the command:

> plot(ridge.cv)

Numbers at top indicate the number of selected features, which take

values of 0,1,…,𝑝 (ridge regression always retrains all 𝑝 features).

Estimating the best model (simplest with smallest CV error)

We know how to find the model with the smallest CV error rate.

However, sometimes it is better to obtain the 𝜆 of the simplest model that is within one standard error of (i.e., not statistically different from)

the best model.

Given our cross-validated object ridge.cv, we can access this 𝜆 using the command ridge.cv$lambda.1se.

We can then find the best-fit, yet simple, ridge logistic regression

model by typing:

> ridge.1se <- glmnet(X, Y, family = "binomial", alpha = 0,

lambda = ridge.cv$lambda.1se)

Obtaining parameter estimates

We can view the parameter estimates using the coef() function:

> coef(ridge.1se)

4 x 1 sparse Matrix of class "dgCMatrix"

(Intercept) -8.128411e+00

balance 3.844218e-03

income 4.968183e-06

student -2.279153e-01

These are similar to our logistic regression estimates:

> glm.mult.fit

Call: glm(formula = default ~ balance + income + student, family = "binomial",

data = Default.recoded)

Coefficients:

(Intercept) balance income student

-1.087e+01 5.737e-03 3.033e-06 -6.468e-01

Degrees of Freedom: 9999 Total (i.e. Null); 9996 Residual

Null Deviance: 2921

Residual Deviance: 1572 AIC: 1580

Predicting probabilities with the fitted model

Given our fitted model ridge.1se, we can now make predictions

from new data with the predict() function in glmnet.

We can predict probabilities of default with:

> predict(ridge.1se, as.matrix(tibble(balance = 1000,

income = 45000, student = 1)), type = "response",

s = ridge.cv$lambda.1se)

[1,] 0.01353914

This is about an order of magnitude larger than for logistic regression:

> predict(glm.mult.fit, tibble(balance = 1000,

income = 45000, student = 1), type = "response")

0.003530389

Predicting classes with the fitted model

We can also predict classes instead of probabilities by setting the

argument type to "class" instead of "response".

We can therefore find the predicted class as:

> predict(ridge.1se, as.matrix(tibble(balance = 1000,

income = 45000, student = 1)), type = "class",

s = ridge.cv$lambda.1se)

[1,] "0"

This is sensible as the predicted default probability associated with

this class < 0.5.

Adding predicted probabilities and classes to the data frame

We may also wish to add predicted class probabilities and classes to

our original data frame:

> estProbs <- predict(ridge.1se, X, type = "response",

s = ridge.cv$lambda.1se)

> estClass <- predict(ridge.1se, X, type = "class",

s = ridge.cv$lambda.1se)

> Default.recoded.withProbs <- Default.recoded %>%

mutate(probs = c(estProbs), pred = c(estClass))

We can peek at this modified data frame using the command:

> head(Default.recoded.withProbs)

default student balance income probs pred

1 0 0 729.5265 44361.625 0.006038429 0

2 0 1 817.1804 12106.135 0.005738740 0

3 0 0 1073.5492 31767.139 0.020966741 0

4 0 0 529.2506 35704.494 0.002687456 0

5 0 0 785.6559 38463.496 0.007267246 0

6 0 1 919.5885 7491.559 0.008293071 0

Obtaining a confusion matrix and accuracy

A confusion matrix for our ridge logistic regression model fit to training

data is obtained by typing:

> Default.recoded.withProbs %>%

select(default, pred) %>%

table()

pred

default 0 1

0 9658 9

1 284 49

Our training accuracy can be obtained by typing:

> Default.recoded.withProbs %>%

summarize(accuracy = mean(pred == default))

accuracy

1 0.9707

Obtaining lasso parameter estimates as a function of 𝜆

The procedure for performing lasso regression in glmnet is akin to

ridge, except we set the argument alpha (𝛼) to one instead of zero:

> lasso.fit <- glmnet(X, Y, family = "binomial", alpha = 1,

lambda = lambdas)

> plot(lasso.fit, xvar = "lambda")

> legend("bottomright", lwd=1, col=1:3, legend = colnames(X))

Using 𝑘-fold CV to fit a lasso logistic regression model

We next wish to perform 𝑘 = 10-fold cross-validation as a function of the set of 𝜆 values in lambdas.

Specifically, we want to use lasso (𝛼 = 1) to calculate the 10-fold cross-validation error rate for each of the 100 𝜆 values.

To do this, we use the cv.glmnet() function and set the argument

family to "binomial" to indicate logistic regression, the argument

alpha (𝛼) to one for lasso regression, provide the argument lambda (𝜆) our set of tuning parameter values lambdas, and set the nfolds

(𝑘) argument equal to 10.

> lasso.cv <- cv.glmnet(X, Y, family = "binomial", alpha = 1,

lambda = lambdas, nfolds = 10)

Plotting the cross-validation error rate as a function of 𝜆

We can view the cross-validation error rate as a function of 𝜆 with the command:

> plot(lasso.cv)

Recall that numbers at the top show the number of selected features,

which take values of 0,1,…,𝑝.

Estimating the best model (simplest with smallest CV error)

Recall that it is sometimes better to obtain the 𝜆 of the simplest model that is within one standard error of (i.e., not statistically different from)

the best model.

Given our cross-validated object lasso.cv, we can access this 𝜆 using the command lasso.cv$lambda.1se.

We can then find this best-fit, yet simple, lasso logistic regression

model by typing:

> lasso.1se <- glmnet(X, Y, family = "binomial", alpha = 1,

lambda = lasso.cv$lambda.1se)

Obtaining parameter estimates

We can view the parameter estimates using the coef() function:

> coef(ridge.1se)

4 x 1 sparse Matrix of class "dgCMatrix"

(Intercept) -8.128411e+00

balance 3.844218e-03

income 4.968183e-06

student -2.279153e-01

These are close to our logistic regression estimates:

> glm.mult.fit

Call: glm(formula = default ~ balance + income + student, family = "binomial",

data = Default.recoded)

Coefficients:

(Intercept) balance income student

-1.087e+01 5.737e-03 3.033e-06 -6.468e-01

Degrees of Freedom: 9999 Total (i.e. Null); 9996 Residual

Null Deviance: 2921

Residual Deviance: 1572 AIC: 1580

Predicting class probabilities and predictions

Let us also predict class probabilities and labels, and add these

predictions to our original data frame in one step:

> estProbs <- predict(lasso.1se, X, type = "response",

s = lasso.cv$lambda.1se)

> estClass <- predict(lasso.1se, X, type = "class",

s = lasso.cv$lambda.1se)

> Default.recoded.withProbs <- Default.recoded %>%

mutate(probs = c(estProbs), pred = c(estClass))

We can peek at this modified data frame using the command:

> head(Default.recoded.withProbs)

default student balance income probs pred

1 0 0 729.5265 44361.625 0.004059856 0

2 0 1 817.1804 12106.135 0.005832021 0

3 0 0 1073.5492 31767.139 0.016725900 0

4 0 0 529.2506 35704.494 0.001771389 0

5 0 0 785.6559 38463.496 0.005120073 0

6 0 1 919.5885 7491.559 0.008895485 0

Obtaining a confusion matrix and accuracy

A confusion matrix for our lasso logistic regression model fit to

training data is obtained by typing:

> Default.recoded.withProbs %>%

select(default, pred) %>%

table()

pred

default 0 1

0 9651 16

1 273 60

Our training accuracy can be obtained by typing:

> Default.recoded.withProbs %>%

summarize(accuracy = mean(pred == default))

accuracy

1 0.9711

Logistic regression with 𝐾 > 2 classes

Recall that logistic regression models the probability that 𝑌 is a particular class 𝑗 given features 𝑋 = (𝑋1,𝑋2,…,𝑋𝑝) with the logistic

function:

𝑃 𝑌 = 𝑗|𝑋 = 𝑒𝛽0+𝛽1𝑋1+𝛽2𝑋2+⋯+𝛽𝑝𝑋𝑝

1 + 𝑒𝛽0+𝛽1𝑋1+𝛽2𝑋2+⋯+𝛽𝑝𝑋𝑝

Thus far, we have only considered 𝐾 = 2 classes.

Logistic regression with 𝐾 > 2 classes is multinomial regression.

We can implement a penalized version of multinomial regression using the glmnet() function, except we choose the family to be

"multinomial" instead of "binomial".

Review: the mpg data frame

234 observations on 11 features for 38 models of cars:

manufacturer Car manufacturer

model Model name

displ Engine displacement (in liters)

year Year of manufacture

cyl Number of cylinders

trans Type of transmission

drv f=front-wheel, r=rear wheel, 4=four wheel drive

cty City miles per gallon

hwy Highway miles per gallon

fl Fuel type

class Type of car

We will perform penalized logistic (multinomial) regression to predict drv (𝐾 = 3) from 𝑝 = 3 quantitative features (displ, cty, and hwy).

R example of multinomial regression classifier

Let’s first convert drv values of f, r, and 4 to classes 1, 2, and 3:

> mpg.recoded <- mpg %>%

mutate(drv = ifelse(drv == "f", 1,

ifelse(drv == "r", 2, 3)))

> mpg.recoded

# A tibble: 234 x 11

manufacturer model displ year cyl trans drv cty hwy fl class

1 audi a4 1.8 1999 4 auto(l5) 1 18 29 p compact

2 audi a4 1.8 1999 4 manual(m5) 1 21 29 p compact

3 audi a4 2 2008 4 manual(m6) 1 20 31 p compact

4 audi a4 2 2008 4 auto(av) 1 21 30 p compact

5 audi a4 2.8 1999 6 auto(l5) 1 16 26 p compact

6 audi a4 2.8 1999 6 manual(m5) 1 18 26 p compact

7 audi a4 3.1 2008 6 auto(av) 1 18 27 p compact

8 audi a4 quattro 1.8 1999 4 manual(m5) 3 18 26 p compact

9 audi a4 quattro 1.8 1999 4 auto(l5) 3 16 25 p compact

10 audi a4 quattro 2 2008 4 manual(m6) 3 20 28 p compact

# ... with 224 more rows

Creating the response and input feature matrices

Recall that glmnet requires that response and feature datasets be

provided separately and as matrices instead of data frames.

We will again use the as.matrix() function to convert data frames

to matrices.

Using this function, we obtain the response (𝑌) and feature (𝑋) matrices as:

> Y <- mpg.recoded %>% select(drv) %>% as.matrix()

> X <- mpg.recoded %>% select(displ, cty, hwy) %>% as.matrix()

The response and input feature matrices

We next use the head() function to peek at the first few rows of the

response and feature matrices:

> head(Y)

drv

[1,] 1

[2,] 1

[3,] 1

[4,] 1

[5,] 1

[6,] 1

> head(X)

displ cty hwy

[1,] 1.8 18 29

[2,] 1.8 21 29

[3,] 2.0 20 31

[4,] 2.0 21 30

[5,] 2.8 16 26

[6,] 2.8 18 26

Using penalized multinomial regression to predict drv

We will use penalized multinomial regression to develop a classifier to

predict whether a car uses front-wheel, rear-wheel, or four-wheel

drive based on its engine size and city and highway mileage.

We will fit this penalized multinomial model with a lasso penalty using

the command:

> lasso.fit <- glmnet(X, Y, family = "multinomial",

alpha = 1, lambda = lambdas)

If we wanted to perform ridge regression, then we would have instead set alpha = 0.

Obtaining lasso parameter estimates as a function of 𝜆

> plot(lasso.fit, xvar = "lambda")

> legend("bottomright", lwd=1, col=1:3, legend = colnames(X))

Obtaining lasso parameter estimates as a function of 𝜆

> plot(lasso.fit, xvar = "lambda")

> legend("bottomright", lwd=1, col=1:3, legend = colnames(X))

Using 𝑘-fold CV to fit a lasso logistic regression model

We next wish to perform 𝑘 = 10-fold cross-validation as a function of the set of 𝜆 values in lambdas.

Specifically, we want to use lasso (𝛼 = 1) to calculate the 10-fold cross-validation error rate for each of the 100 𝜆 values.

To do this, we use the cv.glmnet() function and set the argument

family to "multinomial" to indicate multinomial regression, the

argument alpha (𝛼) to one for lasso regression, provide the argument lambda (𝜆) our set of tuning parameter values lambdas, and set the nfolds (𝑘) argument equal to 10:

> lasso.cv <- cv.glmnet(X, Y, family = "multinomial",

alpha = 1, lambda = lambdas,

nfolds = 10)

Plotting the cross-validation error rate as a function of 𝜆

We can view the cross-validation error rate as a function of 𝜆 with the command:

> plot(lasso.cv)

Recall that numbers at the top show the number of selected features,

which take values of 0,1,…,𝑝.

Estimating the best model (simplest with smallest CV error)

Recall that it is sometimes it is better to obtain the 𝜆 of the simplest model that is within one standard error of (i.e., not statistically

different from) the best model.

Given our cross-validated object lasso.cv, we can access this 𝜆 using the command lasso.cv$lambda.1se.

We can then find this best-fit, yet simple, lasso multinomial regression

model by typing:

> lasso.1se <- glmnet(X, Y, family = "multinomial",

alpha = 1,

lambda = lasso.cv$lambda.1se)

Obtaining parameter estimates

We can view the parameter estimates using the coef() function:

> coef(lasso.1se) $`1`

4 x 1 sparse Matrix of class "dgCMatrix"

-2.2069847

displ .

cty .

hwy 0.1255963

$`2`

4 x 1 sparse Matrix of class "dgCMatrix"

-4.463812

displ 1.081043

cty .

hwy .

$`3`

4 x 1 sparse Matrix of class "dgCMatrix"

6.670797

displ .

cty .

hwy -0.249793

Estimating class probabilities

We may wish to know the probabilities of each of the 𝐾 = 3 classes, which we can estimate as:

> estProbs <- predict(lasso.1se, X, type = "response",

s = lasso.cv$lambda.1se)

We can peek at this set of estimated probabilities with the command:

> as_tibble(estProbs)

[1] # A tibble: 234 x 3

`1.1` `2.1` `3.1`

1 0.867 0.0166 0.116

2 0.867 0.0166 0.116

3 0.924 0.0171 0.0585

4 0.898 0.0189 0.0828

5 0.668 0.0551 0.277

6 0.668 0.0551 0.277

7 0.722 0.0726 0.205

8 0.694 0.0194 0.287

9 0.612 0.0194 0.369

10 0.818 0.0221 0.160

# ... with 224 more rows

Predicting classes and adding them to our data frame

Next, we predict classes and add them to our original data frame:

> estClass <- predict(lasso.1se, X, type = "class",

s = lasso.cv$lambda.1se)

> mpg.recoded.withPreds <- mpg.recoded %>%

mutate(pred = c(estClass))

We can peek at this modified data frame with the command:

> mpg.recoded.withPreds

# A tibble: 234 x 12

manufacturer model displ year cyl trans drv cty hwy fl class pred

1 audi a4 1.8 1999 4 auto(l5) 1 18 29 p compact 1

2 audi a4 1.8 1999 4 manual(m5) 1 21 29 p compact 1

3 audi a4 2 2008 4 manual(m6) 1 20 31 p compact 1

4 audi a4 2 2008 4 auto(av) 1 21 30 p compact 1

5 audi a4 2.8 1999 6 auto(l5) 1 16 26 p compact 1

6 audi a4 2.8 1999 6 manual(m5) 1 18 26 p compact 1

7 audi a4 3.1 2008 6 auto(av) 1 18 27 p compact 1

8 audi a4 quattro 1.8 1999 4 manual(m5) 3 18 26 p compact 1

9 audi a4 quattro 1.8 1999 4 auto(l5) 3 16 25 p compact 1

10 audi a4 quattro 2 2008 4 manual(m6) 3 20 28 p compact 1

# ... with 224 more rows

Obtaining a confusion matrix and accuracy

A confusion matrix for our lasso multinomial regression model fit to

training data is obtained by typing:

> mpg.recoded.withPreds %>%

select(drv, pred) %>%

table()

pred

drv 1 2 3

1 97 1 8

2 4 5 16

3 23 1 79

Our training accuracy can be obtained by typing:

> mpg.recoded.withPreds %>%

summarize(accuracy = mean(pred == drv))

# A tibble: 1 x 1

accuracy

<dbl>

1 0.774

Summary of this lecture

We assess test accuracy of models with cross-validation, which involves

randomly and repeatedly splitting the data into training and test sets,

and then computing mean performance across iterations.

There are two major types of cross-validation: 𝑘-fold and LOOCV.

In feature selection, we try to minimize test error while also minimizing

the number of features used to make predictions.

We discussed three important algorithms for “traditional” feature

selection: best subset, forward, and background selection.

We also discussed three regularization approaches that “shrink”

parameter estimates toward zero: ridge, lasso, and elastic net.