--- title: "LASSO: Shrinkage / Regularization" subtitle: "Notes and in-class exercises" format: html: embed-resources: true toc: 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') ``` You can download the .qmd file for this activity [here](../activity_templates/L06-lasso.qmd) and open in R-studio. The rendered version is posted in the [course website](https://mutasim221b.github.io/Mac-STAT-253-Sp-26/) (Activities tab). I often experiment with the class activities (and see it in live!) and make updates, but I always post the final version before class starts. To be sure you have the most up-to-date copy, please download it once you’ve settled in before class begins. # 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} ![](../images/MLdiagram.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 \vert$ can be. - *Large budgets*, $t$, equate with *small* penalties, $\lambda$. - *Small budgets*, $t$, equate with *large* penalties, $\lambda$. ::: \ \ # Team 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** Q. LASSO vs least squares - What's one advantage of LASSO vs least squares? - Which algorithm(s) require us (or R) to *scale* the predictors? A. LASSO vs least squares - LASSO helps with model selection, i.e. kicks some predictors out of the model, and preventing overfitting. - LASSO requires that the predictors be scaled, but R will do this for us! Least squares does not require us to scale predictors. (As we saw in HW1, scaling predictors doesn't affect our predictions or model quality if use least squares.) Q. What is one advantage of LASSO vs backward stepwise selection? A. LASSO isn’t greedy- LASSO doesn’t overestimate the significance of the predictors it retains (its variable selection isn’t based on p-values). \ **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, LASSO is equivalent to least squares. b. As $\lambda$ increases, the predictor coefficient estimates shrink toward or to 0. c. If $\lambda$ is too big: all predictors are kicked out of the model, and we're left with just an intercept. If $\lambda$ is too small: too few predictors are kicked out, hence the model is complicated and maybe overfit. d. To decide between a LASSO model that uses $\lambda = 0.01$ vs $\lambda = 0.1$, for example, we can compare CV MAE of the LASSO models that use each of these values of $\lambda$ and pick the one with the smaller CV MAE. (Unless the MAEs are similar, in which case we may prefer the model with the simpler model, i.e. larger $\lambda$, even if it has a slightly higher MAE.) \ ______________________________ **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") head(humans) ``` \ \ Let's implement the LASSO. We'll pause to examine the code. The R code notes section, below, and R code tutorial video (see course Schedule) 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" 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$**- Team 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$. > The Goldilocks problem is an idea, not a formal statistical theorem. It comes from the fairy tale Goldilocks and the Three Bears, and in statistics / data science it means: Choosing a value that is neither too small nor too large, but “just right.” CV MAE is large when $\lambda$ is either too small or too big. \ 2. **Picking a $\lambda$ value**- Team 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 (i.e. as simple as possible, but still explains the data well). The plot below adds error bars to the CV MAE *estimates* of prediction error (+/- one standard error). (The CV MAE results are *estimates* because they would change if we had different training data or if we had set a different random seed which would result in different folds.) 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! > 0.2 \ > **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** - Group 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**- Group 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} lasso_workflow %>% finalize_workflow(parameters = best_penalty) %>% fit(data = humans) %>% tidy() %>% filter(estimate != 0) # add this line to print only non-zero coefficients ``` 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. > ankle is not in this larger model but it is in the more parsimonious model. Therefore, the algorithm can’t be greedy. If it were greedy, then ankle would get removed for all smaller models. \ 5. **LASSO vs least squares**- Group 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. > LASSO model is much simpler, and has only slightly worse predictions (on the scale of inches). 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} # LS prediction ls_model %>% predict(new_data = new_patient) ``` ```{r} #| eval: false # LASSO prediction ___ %>% predict(new_data = ___) ``` ```{r} # LASSO prediction final_lasso %>% predict(new_data = new_patient) ``` > They’re very similar! Shrinking coefficients doesn’t mean our predictions are odd. c. Which final model would you choose, the LASSO or least squares? > I’d choose LASSO since it seems to yield comparably accurate predictions but is much simpler than the full least squares model. \ **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 right 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? > thigh b. Approximate the coefficient in the LASSO with $log(\lambda) \approx -5$. > very roughly -0.32 c. At what $log(\lambda)$ does the coefficient start to significantly shrink? > roughly -3.5 d. At what $log(\lambda)$ does the predictor get dropped from the model? > roughly -1.8 \ 8. **plot: big picture** a. How does this plot reflect the LASSO shrinkage phenomenon? > coefficients are shrinking toward or to 0 as $\lambda$ increases (ie $-\log(\lambda)$ decreases) b. What is one of the most "important" or "persistent" predictors? > weight (variable 2), knee (variable 8) c. What is one of the least persistent predictors? > lots of options here. look for the lines that drop to 0 sooner. 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$)? > 11 \ 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_smooth() + geom_hline(yintercept = 0) ``` ```{r} # Now calculate and plot the residuals final_lasso %>% augment(new_data = humans) %>% mutate(.resid = height - .pred) %>% # resid = observed - predicted ggplot(aes(x = .pred, y = .resid)) + geom_point() + geom_smooth() + geom_hline(yintercept = 0) + geom_smooth(se = FALSE) ``` \ 10. **Additional 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 # Show CV results for the chosen lambda # (Can interchange parsimonious_penalty below for best_penalty) lasso_models %>% collect_metrics() %>% filter(penalty==parsimonious_penalty$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)] ``` \ \

# Solutions {.unnumbered .smaller} ## Exercises {.unnumbered .smaller} 1. **Examining the impact of $\lambda$**
Solution a. Yep it matches. ```{r soln-1a} #| eval: true #| echo: true 10^(-5) 10^(-0.1) ``` b. CV MAE is large when $\lambda$ is either too small or too big.

2. **Picking a $\lambda$ value**
Solution a. $\lambda$ close to 0, specifically: ```{r} #| eval: true #| echo: true # plot autoplot(lasso_models) + scale_x_continuous() + xlab(expression(lambda)) # get "best" lambda best_penalty <- lasso_models %>% select_best(metric = "mae") # print out selected lambda best_penalty ``` b. $\lambda$ around 0.2 ```{r} #| eval: true #| echo: true # plot 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) # get "parsimonious" lambda parsimonious_penalty <- lasso_models %>% select_by_one_std_err(metric = "mae", desc(penalty)) # print out selected lambda parsimonious_penalty ``` c. 0.2

3. **Finalizing our LASSO model**
Solution ```{r} #| eval: true #| echo: true # fit final model with selected lambda (parameters = parsimonious_penalty) # note we are fitting to the entire dataset (data = humans) final_lasso <- lasso_workflow %>% finalize_workflow(parameters = parsimonious_penalty) %>% fit(data = humans) # look at coefficient estimates final_lasso %>% tidy() %>% filter(estimate != 0) # optional: focus only on non-zero coefficient estimates ``` a. 5: `weight`, `abdomen`, `thigh`, `knee`, `ankle` b. 3 of the predictors are the same (`weight`, `abdomen`, `thigh`). Interestingly, LASSO includes `knee` (which is the most highly correlated with height in this dataset) but the model chosen by the backward stepwise selection algorithm did not. ```{r} #| eval: true #| echo: true cor(humans)[,'height'] %>% sort() ``` If you're curious, here's some code for implementing backward stepwise section: ```{r backward-stepwise} #| code-fold: true #| eval: true # setup predictors <- c("age", "abdomen", "forearm", "thigh", "neck", "chest", "wrist", "biceps", "hip", "ankle", "weight", "knee") p <- length(predictors) kick_out <- rep(0, p) cvs <- rep(0, p) # model spec lm_spec <- linear_reg() %>% set_mode("regression") %>% set_engine("lm") # loop through predictors for(i in 1:12){ # fit model my_model <- lm_spec %>% fit(as.formula(paste("height ~ ", paste(predictors, collapse = "+"))), data = humans) %>% tidy() %>% filter(term != "(Intercept)") %>% arrange(desc(p.value)) # use 10-fold CV to get MAE for model set.seed(253) cv_process <- lm_spec %>% fit_resamples( as.formula(paste("height ~ ", paste(predictors, collapse = "+"))), resamples = vfold_cv(humans, v = 10), metrics = metric_set(mae) ) %>% collect_metrics() # get name of worst variable (biggest p-value) worst <- as.data.frame(my_model)[1,1] kick_out[i] <- worst # get rid of worst variable from predictor list predictors <- predictors[predictors != worst] # save CV MAE for this model cvs[i] <- as.data.frame(cv_process)$mean } kick_out # final 5 variables: chest, neck, thigh, abdomen, weight cvs # 5-variable model has MAE of 1.337 ``` c. We're using this model to give good predictions, not to explore / make inferences about relationships. d. The remaining predictors are those that have good predictive power in this linear regression model (thus we get conclusions like a hypothesis test without doing a test). Also, our goal is to build a good predictive model, not to do inference.

4. **LASSO vs LASSO**
Solution a. This would have 11 predictors. ```{r eval=TRUE, echo=TRUE} lasso_workflow %>% finalize_workflow(parameters = best_penalty) %>% fit(data = humans) %>% tidy() %>% filter(estimate != 0) # add this line to print only non-zero coefficients ``` b. ankle is not in this larger model but it is in the more parsimonious model. Therefore, the algorithm can't be greedy. If it were greedy, then ankle would get removed for all smaller models.

5. **LASSO vs least squares**
Solution ```{r eval=TRUE, echo=TRUE} # Build the LS model lm_spec <- linear_reg() %>% set_mode("regression") %>% set_engine("lm") # NOTE: we created variable_recipe above # here's what that looked like: ## variable_recipe <- recipe(height ~ ., data = humans) %>% ## step_dummy(all_nominal_predictors()) ls_workflow <- workflow() %>% add_model(lm_spec) %>% add_recipe(variable_recipe) ls_model <- ls_workflow %>% fit(data = humans) # examine coefficients ls_model %>% tidy() # 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. LASSO model is much simpler, and has only slightly worse predictions (on the scale of inches). b. They're very similar! Shrinking coefficients doesn't mean our predictions are odd. ```{r eval=TRUE, echo=TRUE} 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) # LS prediction ls_model %>% predict(new_data = new_patient) # LASSO prediction final_lasso %>% predict(new_data = new_patient) ``` c. I'd choose LASSO since it seems to yield comparably accurate predictions but is much simpler than the full least squares model.

6. **Visualizing LASSO shrinkage.**
Solution ```{r} #| 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) ``` `ggplot` version: ```{r} #| eval: 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 = "none") ```

7. **plot: examining specific predictors**
Solution ```{r} #| fig-width: 8 #| fig-height: 6 #| eval: false #| echo: false # 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) ``` a. thigh b. very roughly -0.32 c. roughly -3.5 d. roughly -1.8

8. **plot: big picture**
Solution a. coefficients are shrinking toward or to 0 as $\lambda$ increases (ie $-\log(\lambda)$ decreases) b. weight (variable 2), knee (variable 8) c. lots of options here. look for the lines that drop to 0 sooner. d. 11

9. **REVIEW: Model evaluation**
Solution ```{r} #| eval: true #| echo: true # use augment to get predictions (.pred) # NOTE: we don't get the residuals (.resid) automatically final_lasso %>% augment(new_data = humans) %>% names() # Now calculate and plot the residuals final_lasso %>% augment(new_data = humans) %>% mutate(.resid = height - .pred) %>% # resid = observed - predicted ggplot(aes(x = .pred, y = .resid)) + geom_point() + geom_smooth() + geom_hline(yintercept = 0) + geom_smooth(se = FALSE) ``` Do you notice any concerning patterns?

10. **OPTIONAL: Practice on another dataset.**
Solution Solutions to this exercise are not provided. Stop by office hours to discuss!

11. **Reflection.**
Solution Solutions to this exercise are not provided. Stop by office hours to discuss!

\ \ # Wrapping Up {.unnumbered .smaller} - [ ] Finish the activity, check the solutions, and reach out with questions. - [ ] Review the R code reference section at the end of today's notes and some optional R code tutorial videos posted for today (see [schedule](../logistics/schedule.qmd)). - [ ] If you're curious, there's also an optional "deeper learning" section that presents two other shrinkage algorithms we won't cover in this course - [ ] **This is the end of Unit 2!** Take some time to pause, review the material we've covered so far, organize your notes, etc. - [ ] **No checkpoint for next class!** See the [schedule](../logistics/schedule.qmd) for optional readings. - [ ] **Finish HW2** (due Tuesday). You have everything you need to complete this assignment after today's class.