Last class

Knowledge check from last week

1) Which model is higher bias: A complex model or a simpler one?

2) Why do we split our data into training and testing datasets?

3) How do we compare models with continuous outcomes?

How forward stepwise selection works: Example from last class

1) Start with a null model (no covariates)

• Your best guess will be the average of the outcome in the training dataset!

2) Test out all models with one covariate, and select the best one:

• E.g. logins∼female$logins\sim female$, logins∼succession$logins\sim succession$, logins∼age$logins\sim age$, ...
• logins∼succession$logins\sim succession$ is the best one (according to RMSE)

3) Test out all models with two covariates, but that have succession$succession$!

• E.g. logins∼succession+female$logins\sim succession+female$, logins∼succession+age$logins\sim succession+age$, logins∼succession+city$logins\sim succession+city$, ...

4) You will end up with k$k$ possible models (k: total number of predictors).

• Choose the best one, depending on the RMSE.

Today: Continuing our journey

What is shrinkage?

• We reviewed the stepwise procedure: Subsetting model selection approach.

  • Select k$k$ out of p$p$ total predictors

• Shrinkage (a.k.a Regularization): Fitting a model with all p$p$ predictors, but introducing bias (i.e. shrinking coefficients towards 0) for improvement in variance.

  • Ridge regression

  • Lasso regression

Ridge Regression: An example

• Predict spending based on frequency of visits to a website

Ordinary Least Squares

• In an OLS: Minimize sum of squared-errors, i.e. minβ∑ni=1(spendi−(β0+β1freqi))2${min}_{\beta }{\sum }_{i=1}^n({\text{spend}}_i-({\beta }_0+{\beta }_1{\text{freq}}_i){)}^2$

What about fit?

• Does the OLS fit the testing data well?

Ridge Regression

• Let's shrink the coefficients!: Ridge Regression

Ridge Regression: What does it do?

• Ridge regression introduces bias to reduce variance in the testing data set.

• In a simple regression (i.e. one regressor/covariate):

minβn∑i=1(yi−β0−xiβ1)2⏟OLS

Ridge Regression: What does it do?

• Ridge regression introduces bias to reduce variance in the testing data set.

• In a simple regression (i.e. one regressor/covariate):

minβn∑i=1(yi−β0−xiβ1)2⏟OLS+λ⋅β21⏟RidgePenalty

• λ is the penalty factor → indicates how much we want to shrink the coefficients.

RMSE on the testing dataset: OLS

RMSE=√144∑i=1(spendi−(132.5−16.25⋅freqi))2=28.36

RMSE on the testing dataset: Ridge Regression

RMSE=√144∑i=1(spendi−(119.5−11.25⋅freqi))2=12.13

Ridge Regression in general

• For regressions that include more than one regressor:

minβn∑i=1(yi−p∑k=0xiβk)2⏟OLS+λ⋅p∑k=1β2k⏟RidgePenalty

• In our previous example, if we had two regressors, female and freq:

minβn∑i=1(spendi−β0−β1femalei−β2freqi)2+λ⋅(β21+β22)

• Because the ridge penalty includes the β's coefficients, scale matters:

  • Standardize variables (you will do that as an option in your code)

How do we choose λ?

Cross-validation!

1) Choose a grid of λ values

• The grid you choose will be context dependent (play around with it!)

2) Compute cross-validation error (e.g. RMSE) for each

3) Choose the smallest one.

λ vs RMSE?

λ vs RMSE? A zoom

How do we do this in R?

library(caret)
set.seed(100)
hbo = read.csv("https://raw.githubusercontent.com/maibennett/sta235/main/exampleSite/content/Classes/Week11/1_Shrinkage/data/hbomax.csv")
lambda_seq = seq(0, 20, length = 500)
ridge = train(logins ~ . - unsubscribe - id, 
            data = train.data, 
            method = "glmnet",
            preProcess = "scale",
            trControl = trainControl("cv", number = 10),
            tuneGrid = expand.grid(alpha = 0, 
                         lambda = lambda_seq)
  )
plot(ridge)

How do we do this in R?

library(caret) 
set.seed(100)
hbo = read.csv("https://raw.githubusercontent.com/maibennett/sta235/main/exampleSite/content/Classes/Week11/1_Shrinkage/data/hbomax.csv")
lambda_seq = seq(0, 20, length = 500)
ridge = train(logins ~ . - unsubscribe - id, 
            data = train.data, 
            method = "glmnet",
            preProcess = "scale",
            trControl = trainControl("cv", number = 10),
            tuneGrid = expand.grid(alpha = 0, 
                         lambda = lambda_seq)
  )
plot(ridge)

How do we do this in R?

library(caret) 
set.seed(100)
hbo = read.csv("https://raw.githubusercontent.com/maibennett/sta235/main/exampleSite/content/Classes/Week11/1_Shrinkage/data/hbomax.csv")
lambda_seq = seq(0, 20, length = 500)
ridge = train(logins ~ . - unsubscribe - id, 
            data = train.data, 
            method = "glmnet",
            preProcess = "scale",
            trControl = trainControl("cv", number = 10),
            tuneGrid = expand.grid(alpha = 0, 
                         lambda = lambda_seq)
  )
plot(ridge)

How do we do this in R?

library(caret) 
set.seed(100)
hbo = read.csv("https://raw.githubusercontent.com/maibennett/sta235/main/exampleSite/content/Classes/Week11/1_Shrinkage/data/hbomax.csv")
lambda_seq = seq(0, 20, length = 500)
ridge = train(logins ~ . - unsubscribe - id,
            data = train.data,
            method = "glmnet",
            preProcess = "scale",
            trControl = trainControl("cv", number = 10),
            tuneGrid = expand.grid(alpha = 0,
                         lambda = lambda_seq)
  )
plot(ridge)

How do we do this in R?

library(caret) 
set.seed(100)
hbo = read.csv("https://raw.githubusercontent.com/maibennett/sta235/main/exampleSite/content/Classes/Week11/1_Shrinkage/data/hbomax.csv")
lambda_seq = seq(0, 20, length = 500)
ridge = train(logins ~ . - unsubscribe - id, 
            data = train.data, 
            method = "glmnet",
            preProcess = "scale",
            trControl = trainControl("cv", number = 10),
            tuneGrid = expand.grid(alpha = 0, 
                         lambda = lambda_seq)
  )
plot(ridge)

How do we do this in R?

lm1 = lm(logins ~ succession + city,
          data = train.data)
coef(lm1)

## (Intercept)  succession        city 
##    7.035888   -6.306371    2.570454

rmse(lm1, test.data)

## [1] 2.089868

coef(ridge$finalModel, ridge$bestTune$lambda)

## 5 x 1 sparse Matrix of class "dgCMatrix"
##                       s1
## (Intercept)  6.564243424
## female       0.002726465
## city         0.824387472
## age          0.046468790
## succession  -2.639308962

rmse(ridge, test.data)

## [1] 2.097452

Lasso regression

• Very similar to ridge regression, except it changes the penalty term:

minβn∑i=1(yi−p∑k=0xiβk)2⏟OLS+λ⋅p∑k=1|βk|⏟LassoPenalty

• In our previous example:

minβn∑i=1(spendi−β0−β1femalei−β2freqi)2+λ⋅(|β1|+|β2|)

• Lasso regression is also called l1 regularization:

||β||1=p∑k=1|β|

Ridge vs Lasso

And how do we do Lasso in R?

library(caret) 
set.seed(100)
hbo = read.csv("https://raw.githubusercontent.com/maibennett/sta235/main/exampleSite/content/Classes/Week11/1_Shrinkage/data/hbomax.csv")
lambda_seq = seq(0, 20, length = 500)
lasso = train(logins ~ . - unsubscribe - id, data = train.data, 
            method = "glmnet",
            preProcess = "scale",
            trControl = trainControl("cv", number = 10),
            tuneGrid = expand.grid(alpha = 1,
                         lambda = lambda_seq)
  )
plot(lasso)

And how do we do Lasso in R?

coef(ridge$finalModel, ridge$bestTune$lambda)

## 5 x 1 sparse Matrix of class "dgCMatrix"
##                       s1
## (Intercept)  6.564243424
## female       0.002726465
## city         0.824387472
## age          0.046468790
## succession  -2.639308962

rmse(ridge, test.data)

## [1] 2.097452

coef(lasso$finalModel, lasso$bestTune$lambda)

## 5 x 1 sparse Matrix of class "dgCMatrix"
##                      s1
## (Intercept)  6.84122778
## female       .         
## city         0.87982819
## age          0.03099797
## succession  -2.83492585

rmse(lasso, test.data)

## [1] 2.09171

A note on binary outcomes

• If we are predicting binary outcomes, RMSE would not be an appropriate measure anymore!

  • We will use accuracy instead: The proportion (%) of correctly classified observations.

• For example:

set.seed(100)
lasso = train(factor(unsubscribe) ~ . - id, data = train.data,
            method = "glmnet", preProcess = "scale",
            trControl = trainControl("cv", number = 10),
            tuneGrid = expand.grid(alpha = 1, lambda = lambda_seq))
pred.values = lasso %>% predict(test.data)
mean(pred.values == test.data$unsubscribe)

## [1] 0.736

A note on binary outcomes

• If we are predicting binary outcomes, RMSE would not be an appropriate measure anymore!

  • We will use accuracy instead: The proportion (%) of correctly classified observations.

• For example:

set.seed(100)
lasso = train(factor(unsubscribe) ~ . - id, data = train.data,
            method = "glmnet", preProcess = "scale",
            trControl = trainControl("cv", number = 10),
            tuneGrid = expand.grid(alpha = 1, lambda = lambda_seq))
pred.values = lasso %>% predict(test.data)
mean(pred.values == test.data$unsubscribe)

## [1] 0.736

A note on binary outcomes

• If we are predicting binary outcomes, RMSE would not be an appropriate measure anymore!

  • We will use accuracy instead: The proportion (%) of correctly classified observations.

• For example:

set.seed(100)
lasso = train(factor(unsubscribe) ~ . - id, data = train.data,
            method = "glmnet", preProcess = "scale",
            trControl = trainControl("cv", number = 10),
            tuneGrid = expand.grid(alpha = 1, lambda = lambda_seq))
pred.values = lasso %>% predict(test.data)
mean(pred.values == test.data$unsubscribe)

## [1] 0.736

Main takeway points

References

• James, G. et al. (2021). "Introduction to Statistical Learning with Applications in R". Springer. Chapter 6.

• STDHA. (2018). "Penalized Regression Essentials: Ridge, Lasso & Elastic Net"