--- title: "LASSO: Shrinkage / Regularization (Notes)" subtitle: "Stat 253" author: "Your Name" format: html: toc: true toc-depth: 2 embed-resources: true --- ```{r include = FALSE} knitr::opts_chunk$set( collapse = TRUE, warning = FALSE, message = FALSE, fig.height = 2.75, fig.width = 4.25, fig.env='figure', fig.pos = 'h', fig.align = 'center') ``` # Learning Goals {.unnumbered .smaller} - Explain how ordinary and penalized least squares are similar and different with regard to (1) the form of the objective function and (2) the goal of variable selection - Explain why variable scaling is important for the performance of shrinkage methods - Explain how the lambda tuning parameter affects model performance and how this is related to overfitting - Describe how output from LASSO models can give a measure of *variable importance* \ # Notes: LASSO {.unnumbered .smaller} ## Context {.unnumbered .smaller} ![](https://mac-stat.github.io/STAT253/images/MLdiagram2.jpg){width=100%} - **world = supervised learning** We want to model some output variable $y$ using a set of *potential* predictors ($x_1, x_2, ..., x_p$). - **task = regression** $y$ is quantitative - **model = linear regression** We'll assume that the relationship between $y$ and ($x_1, x_2, ..., x_p$) can be represented by $$y = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + ... + \beta_p x_p + \varepsilon$$ - [**estimation algorithm = LASSO**]{style="background-color:lightgray;"} (instead of least squares)
## Least Absolute Shrinkage and Selection Operator {.unnumbered .smaller} **LASSO:** **L**east **A**bsolute **S**hrinkage and **S**election **O**perator Dates back to 1996, proposed by Robert Tibshirani (one of the authors of ISLR) > Robert Tibshirani, Regression Shrinkage and Selection Via the Lasso, Journal of the Royal Statistical Society: Series B (Methodological), Volume 58, Issue 1, January 1996, Pages 267–288, [https://doi.org/10.1111/j.2517-6161.1996.tb02080.x](https://doi.org/10.1111/j.2517-6161.1996.tb02080.x) \ ## Goal: Model Selection {.unnumbered .smaller} Use the LASSO algorithm to help us *regularize* and *select* the "best" predictors $x$ to use in a **predictive** linear regression model of $y$: $$y = \hat{\beta}_0 + \hat{\beta}_1 x_1 + \cdots + \hat{\beta}_p x_p + \varepsilon$$ \ ## Big Idea {.unnumbered .smaller} - *Penalize* a predictor for adding complexity to the model (by penalizing its coefficient). - Track whether the predictor's *contribution* to the model (lowering RSS) is enough to offset this penalty. \ ## Algorithm Criterion {.unnumbered .smaller} Identify the model coefficients $\hat{\beta}_1, \hat{\beta}_2, ... \hat{\beta}_p$ that *minimize* the **penalized residual sum of squares**: $$ \begin{aligned} RSS + \lambda \sum_{j=1}^p \vert \hat{\beta}_j\vert &= \sum_{i=1}^n (y_i - \hat{y}_i)^2 + \lambda \sum_{j=1}^p \vert \hat{\beta}_j\vert \\ & = \sum_{i=1}^n (y_i - \hat{\beta}_0 - \hat{\beta}_1 x_{1i} - \cdots - \hat{\beta}_p x_{pi})^2 + \lambda \sum_{j=1}^p \vert \hat{\beta}_j\vert \end{aligned} $$ where - residual sum of squares (RSS) measures the overall model prediction error - the penalty term measures the overall size of the model coefficients - $\lambda \ge 0$ ("lambda") is a tuning parameter :::{.callout-note} This problem is equivalent to finding the model coefficients $\hat{\beta}_1, \hat{\beta}_2, ... \hat{\beta}_p$ that minimize the residual sum of squares *subject to the constraint* that $\sum_{j=1}^p \vert \beta_j \vert \le t$ where $t$ is a *budget* for how large $\sum_{j=1}^p \vert \beta_j$ can be. - *Large budgets*, $t$, equate with *small* penalties, $\lambda$. - *Small budgets*, $t$, equate with *large* penalties, $\lambda$. ::: \ \ # Small Group Discussion {.unnumbered .smaller} Discuss basic understanding from the video to help each other clear up concepts. ## Questions {.unnumbered .smaller} **1: LASSO vs other algorithms for building linear regression models** a. LASSO vs least squares - What's one advantage of LASSO vs least squares? - Which algorithm(s) require us (or R) to *scale* the predictors? b. What is one advantage of LASSO vs backward stepwise selection? \ **2: LASSO tuning** We have to pick a $\lambda$ penalty tuning parameter for our LASSO model. What's the impact of $\lambda$? a. When $\lambda$ is 0, ... b. As $\lambda$ increases, the predictor coefficient estimates .... c. **Goldilocks problem**: If $\lambda$ is too big, .... If $\lambda$ is too small, ... d. To decide between a LASSO that uses $\lambda = 0.01$ vs $\lambda = 0.1$ (for example), we can .... \ ______________________________ **COMMENT: Picking $\lambda$** We *cannot* know the "best" value for $\lambda$ in advance. This varies from analysis to analysis. We must try a reasonable *range* of possible values for $\lambda$. This also varies from analysis to analysis. In general, we have to use *trial-and-error* to identify a range that is... - wide enough that it doesn't miss the best values for $\lambda$ - narrow enough that it focuses on reasonable values for $\lambda$ \ \ # Exercises {-} ## Instructions {.unnumbered .smaller} - Open the QMD fo today and scroll down to the *Exercises* - Work on implementing LASSO to familiar data - Become familiar with the new code **structures**: - instead of `fit_resamples` to run CV, we'll use `tune_grid` to tune the algorithm with CV - new engine: `set_engine('glmnet')` - in general: focus on the concepts over the R code - As always: - Be kind to yourself/each other - Collaborate - Ask me questions as I move around the room \ ## Questions {.unnumbered .smaller} We'll use the LASSO algorithm to help us build a good *predictive* model of `height` using the collection of 12 possible predictors in the `humans` dataset: ```{r 06-read-data} #| eval: true #| echo: true # Load packages library(tidyverse) library(tidymodels) # Resolves package conflicts by preferring tidymodels functions tidymodels_prefer() # Load data humans <- read.csv("https://mac-stat.github.io/data/bodyfat2.csv") ``` \ \ Let's implement the LASSO. We'll pause to examine the code. The R code notes section, below, and R code tutorial video (see [Schedule](../logistics/schedule.html)) provide more detail. ```{r} #| eval: true #| echo: true # STEP 1: LASSO algorithm / model specification # NOTE: we're using a new engine now: glmnet instead of lm lasso_spec <- linear_reg() %>% set_mode("regression") %>% set_engine("glmnet") %>% set_args(mixture = 1, penalty = tune()) ``` ```{r} #| eval: true #| echo: true # STEP 2: variable recipe # NOTE: "y ~ ." is shorthand for "y as a function of all other variables" # NOTE: We'll discuss step_dummy() next class. variable_recipe <- recipe(height ~ ., data = humans) %>% step_dummy(all_nominal_predictors()) ``` ```{r} #| eval: true #| echo: true # STEP 3: workflow specification (model + recipe) lasso_workflow <- workflow() %>% add_recipe(variable_recipe) %>% add_model(lasso_spec) ``` ```{r} #| eval: true #| echo: true # STEP 4: Estimate 50 LASSO models using # lambda values on a "grid" or range from 10^(-5) to 10^(-0.1). # Calculate the CV MAE for each of the 50 models. # NOTE: we use tune_grid instead of fit_resamples to run CV # NOTE: I usually start with a range from 10^(-5) to 10^1 and tweak through trial-and-error. set.seed(253) lasso_models <- lasso_workflow %>% tune_grid( grid = grid_regular(penalty(range = c(-5, -0.1)), levels = 50), resamples = vfold_cv(humans, v = 10), metrics = metric_set(mae) ) ``` \ 1. **Examining the impact of $\lambda$** Let's compare the CV MAEs (y-axis) for our **50** LASSO models which used 50 different $\lambda$ values (x-axis): ```{r} #| eval: true #| code-fold: true autoplot(lasso_models) + scale_x_continuous() + xlab(expression(lambda)) ``` a. We told R to use a range of $\lambda$ from -5 to -0.1 on the log10 scale. Calculate this range on the non-log scale and confirm that it matches the x-axis. ```{r} 10^(-5) 10^(-0.1) ``` b. Explain why this plot displays the "Goldilocks" problem of tuning $\lambda$. \ 2. **Picking a $\lambda$ value** a. In the plot above, roughly which value of the $\lambda$ penalty parameter produces the *smallest* CV MAE? Check your approximation: ```{r} # get "best" lambda best_penalty <- lasso_models %>% select_best(metric = "mae") # print out result best_penalty ``` b. Suppose we prefer a **parsimonious** model. The plot below adds error bars to the CV MAE *estimates* of prediction error (+/- one standard error). Any model with a CV MAE that falls within another model's error bars is not significantly better or worse at prediction: ```{r} #| eval: true #| code-fold: true # with error bars # NOTE: we start with the same code as above (lines 1--3), # then add error bars (`geom_errorbar`) autoplot(lasso_models) + scale_x_continuous() + xlab(expression(lambda)) + geom_errorbar(data = collect_metrics(lasso_models), aes(x = penalty, ymin = mean - std_err, ymax = mean + std_err), alpha = 0.5) ``` Use this to approximate the largest $\lambda$, thus the most simple LASSO model, that produces a CV MAE that's within 1 standard error of the best model (thus is not significantly worse). Check your approximation: ```{r} #| eval: true # get "parsimonous" lambda parsimonious_penalty <- lasso_models %>% select_by_one_std_err(metric = "mae", desc(penalty)) ``` ```{r} # print out selected lambda parsimonious_penalty ``` c. Moving forward, we'll use the parsimonious LASSO model. Simply report the tuning parameter $\lambda$ here. Just as a radio show needs to tell its audience where to tune the radio dial, it's important to explicitly report $\lambda$ so that we and others can reproduce the model! \ > **PAUSE:** Picking a range to try for $\lambda$ > > The range of values we tried for $\lambda$ had the following nice properties. > If it didn't, we should adjust our range (make it narrower or wider). > > - Our range was wide enough. > We observed the goldilocks effect. Further, the "best" and "parsimonious" $\lambda$ values were *not* at the edges of the range, suggesting there aren't better $\lambda$ values outside our range. > > - Our range was narrow enough. > We didn't observe any loooooong flat lines in CV MAE, thus we narrowed in on the $\lambda$ values where the "action is happening", i.e. where changing $\lambda$ impacts the model. \ 3. **Finalizing our LASSO model** Let's *finalize* our parsimonious LASSO model: ```{r} #| eval: true # fit final LASSO model # using selected lambda (parameters = parsimonious_penalty) # and full dataset (data = humans) final_lasso <- lasso_workflow %>% finalize_workflow(parameters = parsimonious_penalty) %>% fit(data = humans) ``` ```{r} # look at coefficients final_lasso %>% tidy() ``` a. How many and which predictors were *kept* in this model? b. How do these compare to the 5-predictor model identified using the backward stepwise selection algorithm with this subset of data: weight, abdomen, thigh, neck, chest c. Through shrinkage, the LASSO coefficients lose some contextual meaning, so we typically shouldn't interpret them. Why don't we care?! THINK: What is the goal of LASSO modeling? d. The LASSO `tidy()` summary doesn't report p-values for testing the "significance" of our predictors. Why don't we care? (Name two reasons.) \ 4. **LASSO vs LASSO** a. Our parsimonious LASSO selected only 5 of the 12 possible predictors. Out of curiosity, how many predictors would have remained if we had used the `best_penalty` value for $\lambda$? ```{r eval = FALSE} lasso_workflow %>% finalize_workflow(parameters = ___) %>% fit(data = humans) %>% tidy() ``` ```{r} ``` b. Based on this example, do you think LASSO is a greedy algorithm? Are you "stuck" with your past locally optimal choices? Compare the predictors in this larger model with those in the smaller, parsimonious model. \ 5. **LASSO vs least squares** Let's compare our `final_lasso` model to the least squares model using all predictors: ```{r} #| eval: true # Build the least squares model using recipes and workflows lm_spec <- linear_reg() %>% set_mode("regression") %>% set_engine("lm") ls_workflow <- workflow() %>% add_model(lm_spec) %>% add_recipe(variable_recipe) ls_model <- ls_workflow %>% fit(data = humans) ``` ```{r} # examine coefficients ls_model %>% tidy() ``` ```{r} # get 10-fold CV MAE set.seed(253) ls_workflow %>% fit_resamples( resamples = vfold_cv(humans,v = 10), metrics = metric_set(mae) ) %>% collect_metrics() ``` a. Our `final_lasso` has 5 predictors and a CV MAE of 1.9 (calculated above). The `ls_model` has 12 predictors and a CV MAE of 1.8 (confirm). Comment. b. Use both `final_lasso` and `ls_model` to predict the height of the new patient below. How do these compare? Does this *add to* or *calm* any fears you might have had about shrinking coefficients?! ```{r} new_patient <- data.frame(age = 50, weight = 200, neck = 40, chest = 115, abdomen = 105, hip = 100, thigh = 60, knee = 38, ankle = 23, biceps = 32, forearm = 29, wrist = 19) ``` ```{r} #| eval: false # LS prediction ___ %>% predict(new_data = ___) ``` ```{r} ``` ```{r} #| eval: false # LASSO prediction ___ %>% predict(new_data = ___) ``` ```{r} ``` c. Which final model would you choose, the LASSO or least squares? \ **6. Visualizing LASSO shrinkage** Finally, let's zoom back out and compare the coefficients for all 50 LASSO models: ```{r} #| fig-width: 8 #| fig-height: 6 #| eval: true # Get output for each LASSO model all_lassos <- final_lasso %>% extract_fit_parsnip() %>% pluck("fit") # Plot coefficient paths as a function of lambda plot(all_lassos, xvar = "lambda", label = TRUE, col = rainbow(20)) # Codebook for which variables the numbers correspond to rownames(all_lassos$beta) ``` There's a *lot* of information in this plot! - **lines** = each line represents a different predictor. The small number to the left of each line indicates the predictor by its order in the `rownames()` list. Click "Zoom" to zoom in. - **x-axis** = our range of $\lambda$ values, on the log scale - **y-axis** = coefficient values at the corresponding $\lambda$ - **numbers above the plot** = how many predictors remain in the model with the corresponding $\lambda$ We'll process this information in the next 2 exercises. If you're curious, here is some code to recreate this plot using `ggplot`: ```{r} #| code-fold: true lasso_coefs <- all_lassos$beta %>% as.matrix() %>% t() %>% as.data.frame() %>% mutate(lambda = all_lassos$lambda ) %>% pivot_longer(cols = -lambda, names_to = "term", values_to = "coef") lasso_coefs %>% filter(coef != 0, lambda > 1) lasso_coefs %>% ggplot(aes(x = lambda, y = coef, color = term)) + geom_line() + geom_text(data = lasso_coefs %>% filter(lambda == min(lambda)), aes(x = 0.001, label = term), size = 2) + scale_x_log10(limits = c(-2, 2)) + guides(color = FALSE) ``` \ 7. **plot: examining specific predictors** Answer the following questions for predictor 7. a. Which predictor is this? b. Approximate the coefficient in the LASSO with $log(\lambda) \approx -5$. c. At what $log(\lambda)$ does the coefficient start to significantly shrink? d. At what $log(\lambda)$ does the predictor get dropped from the model? \ 8. **plot: big picture** a. How does this plot reflect the LASSO shrinkage phenomenon? b. What is one of the most "important" or "persistent" predictors? c. What is one of the least persistent predictors? d. Our parsimonious LASSO model had 5 predictors. How many predictors would remain if we had *minimized* the CV MAE using $\lambda \approx 0.0126$ ($log(\lambda) = -4.4$)? \ 9. **REVIEW: Model evaluation** Let's finalize our LASSO analysis. Just as in least squares, it's important to evaluate a LASSO model before applying it. We've already examined whether our LASSO model produces accurate predictions. Use a **residual plot** to determine if this model is *wrong*. NOTE: `augment()` gives predictions, but not residuals :/. You'll need to calculate them. ```{r} # Note what augment() gives us final_lasso %>% augment(new_data = humans) %>% names() ``` ```{r} #| eval: false # Now calculate and plot the residuals final_lasso %>% augment(new_data = humans) %>% mutate(.resid = ___) %>% ggplot(aes(x = ___, y = ___)) + geom_point() + geom_hline(yintercept = 0) ``` ```{r} ``` \ 10. **OPTIONAL: Practice on another dataset**. The `Hitters` data in the `ISLR` package contains the salaries and performance measures for 322 Major League Baseball players. Use LASSO to determine the "best" predictive model of player `Salary`. ```{r} # Load the data library(ISLR) data(Hitters) Hitters <- Hitters %>% filter(!is.na(Salary)) # IN THE CONSOLE (not in the QMD): Examine codebook #?Hitters ``` ```{r} ``` \ 11. **Reflection** This is the end of the (short!) Unit 2 on "Regression: Model Selection"! Let's reflect on the technical content of this unit: - What was the main motivation / goal behind this unit? - For each of the following algorithms, describe the steps, pros, cons, and comparisons to least squares: - best subset selection - backward stepwise selection - LASSO - In your own words, define the following: parsimonious models, greedy algorithms, Goldilocks problem. - Review the new tidymodels syntax from this unit. Identify key themes and patterns. \ \ # Deeper Learning (optional) {.unnumbered .smaller} If you're curious, consider how the LASSO connects to some other algorithms we won't cover in STAT 253. ## Ridge Regression {.unnumbered .smaller} LASSO isn't the only shrinkage & regularization algorithm. An alternative is **ridge regression**. This algorithm also seeks to build a (predictive) linear regression model of $y$: $$y = \hat{\beta}_0 + \hat{\beta}_1 x_1 + \cdots + \hat{\beta}_p x_p + \varepsilon$$ It also does so by selecting coefficients which **minimize a penalized residual sum of squares**. HOWEVER, the ridge regression penalty is based upon the sum of *squared* coefficients instead of the sum of *absolute* coefficients: $$RSS + \lambda \sum_{j=1}^p \hat{\beta}_j^2$$ This penalty regularizes / shrinks the coefficients. BUT, unlike the LASSO, ridge regression does NOT shrink coefficients to 0, thus cannot be used for variable selection. (Check out the ISLR text for a more rigorous, geometric explanation for why LASSO often shrinks coefficients to 0.) ## Elastic Net {.unnumbered .smaller} The **elastic net** is yet another shrinkage & regularization algorithm. It *combines* the penalties used in LASSO and ridge regression. This algorithm seeks to build a (predictive) linear regression model of $y$: $$y = \hat{\beta}_0 + \hat{\beta}_1 x_1 + \cdots + \hat{\beta}_p x_p + \varepsilon$$ by selecting coefficients which minimize the following penalized residual sum of squares: $$RSS + \lambda_1 \sum_{j=1}^p \vert \hat{\beta}_j \vert + \lambda_2 \sum_{j=1}^p \hat{\beta}_j^2$$ **NOTE:** - Elastic net depends upon **two** tuning parameters, $\lambda_1$ and $\lambda_2$, thus is more complicated than the LASSO. - In cases when we have a group of correlated predictors, LASSO tends to select only one of these predictors. The elastic net does not. ## Bayesian Connections {.unnumbered .smaller} IF you have taken or will take *Bayesian statistics* (STAT 454), we can also write the LASSO as a Bayesian model. Specifically, LASSO estimates are equivalent to the posterior mode estimates of $\beta_j$. $$\begin{split} Y_i | \beta_0,...\beta_k & \sim N(\beta_0 + \beta_1X_1 + \cdots + \beta_p X_p, \sigma^2) \\ \beta_j & \sim \text{ Laplace (double-exponential)}(0, f(\lambda)) \\ \sigma^2 & \sim \text{ some prior} \\ \end{split}$$ \ \ # Notes: R Code {.unnumbered .smaller} Suppose we want to build a model of response variable `y` using *all* possible predictors in our `sample_data`. ```{r eval = FALSE} # Load packages library(tidymodels) ``` ------------------------------------ **Build the model for a range of tuning parameters** ```{r eval = FALSE} # STEP 1: LASSO model specification lasso_spec <- linear_reg() %>% set_mode("regression") %>% set_engine("glmnet") %>% set_args(mixture = 1, penalty = tune()) ``` STEP 1 notes: - We use the `glmnet`, not `lm`, engine to build the LASSO. - The `glmnet` engine requires us to specify some arguments (`set_args`): - `mixture = 1` indicates LASSO. Changing this would run a different regularization algorithm (see Deeper Learning, above). - `penalty = tune()` indicates that we don't (yet) know an appropriate $\lambda$ penalty term. We need to tune it. ```{r eval = FALSE} # STEP 2: variable recipe # (You can add pre-processing steps. We will discuss step_dummy() in the next class.) variable_recipe <- recipe(y ~ ., data = sample_data) %>% step_dummy(all_nominal_predictors()) ``` ```{r eval = FALSE} # STEP 3: workflow specification (model + recipe) lasso_workflow <- workflow() %>% add_recipe(variable_recipe) %>% add_model(lasso_spec) ``` ```{r eval = FALSE} # STEP 4: Estimate multiple LASSO models using a range of possible lambda values set.seed(___) lasso_models <- lasso_workflow %>% tune_grid( grid = grid_regular(penalty(range = c(___, ___)), levels = ___), resamples = vfold_cv(sample_data, v = ___), metrics = metric_set(mae) ) ``` STEP 4 notes: - Since the CV process is random, we need to `set.seed(___)`. - We use `tune_grid()` instead of `fit()` since we have to build multiple LASSO models, each using a different tuning parameter. - `grid` specifies the values of tuning parameter $\lambda$ that we want to try. - `penalty(range = c(___, ___))` specifies a range of $\lambda$ values we want to try, on the log10 scale. You might start with `c(-5, 1)`, hence $\lambda$ from 0.00001 ($10^{-5}$) to 10 ($10^1$), and adjust from there. - `levels` is the *number* of $\lambda$ values to try in that range, thus how many LASSO models to build. - `resamples` and `metrics` indicate that we want to calculate a CV MAE for each LASSO model. \ **Tuning $\lambda$** ```{r eval = FALSE} # Calculate CV MAE for each LASSO model lasso_models %>% collect_metrics() # Plot CV MAE (y-axis) for the LASSO model from each lambda (x-axis) autoplot(lasso_models) + scale_x_log10() # plot lambda on log10 scale autoplot(lasso_models) + scale_x_continuous() + # plot lambda on original scale xlab(expression(lambda)) # CV MAE plot with error bars (+/- 1 standard error) # With error bars autoplot(lasso_models) + scale_x_continuous() + xlab(expression(lambda)) + geom_errorbar(data = collect_metrics(lasso_models), aes(x = penalty, ymin = mean - std_err, ymax = mean + std_err), alpha = 0.5) # Identify lambda which produced the lowest ("best") CV MAE best_penalty <- lasso_models %>% select_best(metric = "mae") best_penalty # Identify the largest lambda (hence simplest LASSO) for which the CV MAE is # larger but "roughly as good" (within one standard error of the lowest) parsimonious_penalty <- lasso_models %>% select_by_one_std_err(metric = "mae", desc(penalty)) parsimonious_penalty ``` \ **Finalizing the "best" LASSO model** ```{r eval = FALSE} # parameters = final lambda value (best_penalty or parsimonious_penalty) final_lasso_model <- lasso_workflow %>% finalize_workflow(parameters = ___) %>% fit(data = sample_data) # Check it out final_lasso_model %>% tidy() ``` \ **Using the final LASSO model to make predictions** ```{r eval = FALSE} # new data = some data frame with observations on each predictor final_lasso_model %>% predict(new_data = ___) ``` \ **Visualizing shrinkage: comparing LASSO coefficients under each $\lambda$** ```{r eval = FALSE} # Get output for each LASSO all_lassos <- final_lasso_model %>% extract_fit_parsnip() %>% pluck("fit") # Plot coefficient paths as a function of lambda plot(all_lassos, xvar = "lambda", label = TRUE, col = rainbow(20)) # Codebook for which variables the numbers correspond to rownames(all_lassos$beta) # e.g., What are variables 2 and 4? rownames(all_lassos$beta)[c(2,4)] ``` \ \ # Done! - Render your notes. - Check the solutions in the course website (Solution drop downs). - If you finish all that during class, start your homework!