---
title: "Supervised Learning"
subtitle: "Regularization"
author: "Ron Yurko"
date: "06/23/2020"
---

```{r setup, echo = FALSE, warning = FALSE, message = FALSE}
library(tidyverse)
knitr::opts_chunk$set(fig.retina = 3, warning = FALSE, message = FALSE)
```


## Yesterday...

We talked about variable selection in the linear model context:

$$
\hat{Y} = \hat{\beta}_0 + \hat{\beta}_1 X_1 + \cdots + \hat{\beta}_p X_p
$$

Why would we attempt to select a __subset__ of the $p$ variables?


--
- *To improve prediction accuracy* 
  - Eliminating uninformative predictors can lead to lower variance in the test-set MSE, at the expense of a slight increase in bias


--
- *To improve model interpretability*
  - Eliminating uninformative predictors is obviously a good thing when your goal is to tell the story of how your predictors are associated with your response.


--
__Remember the bias-variance tradeoff__:


$$
{\rm MSE} = {\rm (Bias)}^2 + {\rm Variance}
$$

- Introduce bias but decrease variance to improve predictions

---

## Shrinkage methods: Ridge regression

__Remember__: linear regression estimates coefficients by minimizing:

 $$\text{RSS} = \sum_i^n \big(Y_i - \beta_0 -  \sum_j^p \hat{\beta}_j X_{ij}\big)^2$$


--
__Ridge regression__ introduces a __shrinkage penalty__ $\lambda \geq 0$ by minimizing:

$$\sum_i^n \big(Y_i - \beta_0 -  \sum_j^p \beta_j X_{ij}\big)^2 + \lambda \sum_j^p \beta_j^2 = \text{RSS} + \lambda \sum_j^p \beta_j^2$$

--
- as $\lambda$ increases $\Rightarrow$ flexibility of models decreases

  - __increases bias, but decreases variance__
 

--
- for fixed value of $\lambda$, ridge regression fits only a single model

  - need to use cross-validation to __tune__ $\lambda$

---

## Shrinkage methods: Ridge regression

For example: note how the magnitude of the coefficient for `Income` trends as $\lambda \rightarrow \infty$

```{r out.width='40%', echo = FALSE, fig.align='center'}
knitr::include_graphics("http://www.stat.cmu.edu/~pfreeman/Ridge.png")
```

--

The coefficient __shrinks towards zero__, but never actually reaches it

- `Income` is always a variable in the learned model, regardless of the value of $\lambda$


---

## Shrinkage methods: Lasso regression

Ridge regression __keeps all variables__ 

But we may believe there is a __sparse__ solution


--
__Lasso__ enables variable selection with $\lambda$ by minimizing:

$$\sum_i^n \big(Y_i - \beta_0 -  \sum_j^p \beta_j X_{ij}\big)^2 + \lambda \sum_j^p\vert  \beta_j \vert = \text{RSS} + \lambda \sum_j^p \vert \beta_j \vert$$

- Lasso uses an $\ell_1$ ("ell 1") penalty


--
- as $\lambda$ increases $\Rightarrow$ flexibility of models decreases

  - __increases bias, but decreases variance__
  
--

- Can handle the $p > n$ case, i.e. more variables than observations!
 

---

## Shrinkage methods: Lasso regression

Lasso regression __performs variable selection__ yielding __sparse__ models

```{r out.width='40%', echo = FALSE, fig.align='center'}
knitr::include_graphics("http://www.stat.cmu.edu/~pfreeman/Lasso.png")
```


--

The coefficient shrinks towards and __eventually equals zero__ at $\lambda ~ 1000$

- if the optimum value of $\lambda$ is larger, then `Income` would NOT be included in the learned model

---

## Which do we use?


--

```{r out.width='80%', echo = FALSE, fig.align='center'}
knitr::include_graphics("https://i.gifer.com/Xnj.gif")
```


---

## Best of both worlds? Elastic net

$$\sum_{i}^{n}\left(Y_{i}-\beta_{0}-\sum_{j}^{p} \beta_{j} X_{i j}\right)^{2}+\lambda\left[(1-\alpha)\|\beta\|_{2}^{2} / 2+\alpha\|\beta\|_{1} \right]$$
- $\vert \vert \beta \vert \vert_1$ is the $\ell_1$ norm: $\vert \vert \beta \vert \vert_1 = \sum_j^p \vert \beta_j \vert$

- $\vert \vert \beta \vert \vert_2$ is the $\ell_2$, Euclidean, norm: $\vert \vert \beta \vert \vert_2 = \sqrt{\sum_j^p \beta_j^2}$


--

- Ridge penalty: $\lambda \cdot (1 - \alpha) / 2$

- Lass penalty: $\lambda \cdot \alpha$


--
- $\alpha$ controls the __mixing__ between the two types, ranges from 0 to 1

  - $\alpha = 1$ returns lasso
  
  - $\alpha = 0$ return ridge

---

## Caveats to consider...

- For either ridge, lasso, or elastic net: __you should standardize your data__

- Common convention: within each column, compute then subtract off the sample mean, and compute the divide off the sample standard deviation:

$$\tilde{X}_{ij} = \frac{X_{ij} - \bar{X}_j}{s_{X,j}}$$

- [`glmnet`](https://web.stanford.edu/~hastie/glmnet/glmnet_alpha.html) package does this by default and reports coefficients on the original scale


--
- $\lambda$ and $\alpha$ are __tuning parameters__

- Have to select appropriate values based on test data / cross-validation

- When using `glmnet`, the `cv.glmnet()` function will perform the cross-validation for you

---

## Regularization example data

Will use the `Hitters` dataset from the [`ISLR`](https://cran.r-project.org/web/packages/ISLR/index.html) package

```{r}
library(ISLR)
data("Hitters")
Hitters <- as_tibble(Hitters) %>%
  filter(!is.na(Salary))
Hitters
```

---

## Introduction to `glmnet`

We will use the [`glmnet`](https://web.stanford.edu/~hastie/glmnet/glmnet_alpha.html#lin) package for ridge, lasso, and elastic net

```{r load-glmnet}
library(glmnet)
```


--
__Important__: `glmnet` does NOT use formula but instead `x` matrix of predictors and `y` vector of response

- use the `model.matrix()` function (which converts factors to 0-1 dummy variables!)

```{r model-matrix}
model_x <- model.matrix(Salary ~ ., Hitters)[, -1]
model_y <- Hitters$Salary
```


--
What do the initial regression coefficients look like?

```{r init-lm}
init_reg_fit <- lm(Salary ~ ., Hitters)
coef(init_reg_fit)
```

---

## Ridge regression example

We use cross-validation to select $\lambda$ with `cv.glmnet()` which uses 10-folds by default

- specify ridge regression with `alpha = 0`

```{r ridge-ex, fig.align='center', fig.height=5}
fit_ridge_cv <- cv.glmnet(model_x, model_y, alpha = 0)
plot(fit_ridge_cv)
```


---

## Ridge regression example


.pull-left[
Coefficients using the __1 standard error rule__ $\lambda$

```{r}
coef(fit_ridge_cv)
```

]
.pull-right[
Or specify $\lambda$ with minimum cross-validation error

```{r}
coef(fit_ridge_cv, fit_ridge_cv$lambda.min)

```


]


---

## Lasso regression example

Similar to ridge by specify `alpha = 1`:

```{r lasso-ex, fig.align='center', fig.height=5}
fit_lasso_cv <- cv.glmnet(model_x, model_y, alpha = 1)
plot(fit_lasso_cv)
```


---

## Lasso regression example

.pull-left[
Coefficients using the __1 standard error rule__ $\lambda$

```{r}
coef(fit_lasso_cv)
```

]
.pull-right[
Or specify $\lambda$ with minimum cross-validation error

```{r}
coef(fit_lasso_cv, fit_lasso_cv$lambda.min)

```


]

---

## Elastic net example

Need to tune both $\lambda$ and $\alpha$ - can do so manually with our own folds

```{r}
set.seed(2020)
fold_id <- sample(rep(1:10, length.out = nrow(model_x)))
```

Then use cross-validation with these folds for different candidate `alpha` values:

```{r}
cv_en_25 <- cv.glmnet(model_x, model_y, foldid = fold_id, alpha = .25)
cv_en_50 <- cv.glmnet(model_x, model_y, foldid = fold_id, alpha = .5)
cv_ridge <- cv.glmnet(model_x, model_y, foldid = fold_id, alpha = 0)
cv_lasso <- cv.glmnet(model_x, model_y, foldid = fold_id, alpha = 1)
```

Can see which one had the lowest CV error among its candidate $\lambda$ values:

```{r}
which.min(c(min(cv_en_25$cvm), min(cv_en_50$cvm), min(cv_ridge$cvm), min(cv_lasso$cvm)))
```



