## Context and Data

To build LASSO models in `tidymodels`, first load the package and set the seed for the random number generator to ensure reproducible results:

```{r}
library(dplyr)
library(readr)
library(broom)
library(ggplot2)
library(tidymodels) 
tidymodels_prefer() # Resolves conflicts, prefers tidymodel functions
set.seed(___) # Pick your favorite number to fill in the parentheses
```

Then adapt the following code:

```{r}
# Lasso Model Spec
lm_lasso_spec <- 
  linear_reg() %>%
  set_args(mixture = 1, penalty = 0) %>% ## mixture = 1 indicates Lasso, we'll choose penalty later
  set_engine(engine = 'glmnet') %>% # glmnet does regularization (LASSO, ridge, elastic net)
  set_mode('regression') 

# Recipe with standardization (!)
data_rec <- recipe( ___ ~ ___ , data = ___) %>%
    step_nzv(all_predictors()) %>% # removes variables with the same value
    step_novel(all_nominal_predictors()) %>% # important if you have rare categorical variables 
    step_normalize(all_numeric_predictors()) %>%  # important standardization step for LASSO
    step_dummy(all_nominal_predictors())  # creates indicator variables for categorical variables

# Workflow (Recipe + Model)
lasso_wf <- workflow() %>% 
  add_recipe(data_rec) %>%
  add_model(lm_lasso_spec)

# Fit Model
lasso_fit <- lasso_wf %>% 
  fit(data = ___) # Fit to data
```

<br><br><br>

```{r}
set.seed(123)

cars2018 <- read_csv("https://raw.githubusercontent.com/juliasilge/supervised-ML-case-studies-course/master/data/cars2018.csv")

head(cars2018)

# Cleaning
cars2018 <- cars2018 %>%
    select(-model_index)
```


## Exercise 1: A least squares model

Let’s start by building an ordinary (not penalized) least squares model to review important concepts. We’ll fit a model to predict fuel efficiency measured in miles per gallon (mpg) with all possible predictors.

```{r}
lm_spec <-
    linear_reg() %>% 
    set_engine(engine = 'lm') %>% 
    set_mode('regression')

full_rec <- recipe(mpg ~ ., data = cars2018) %>%
    update_role(model, new_role = 'ID') %>% # we want to keep the name of the car model but not as a predictor or outcome
    step_nzv(all_predictors()) %>% # removes variables with the same value
    step_normalize(all_numeric_predictors()) %>% # important standardization step for LASSO
    step_dummy(all_nominal_predictors())  # creates indicator variables for categorical variables

full_lm_wf <- workflow() %>%
    add_recipe(full_rec) %>%
    add_model(lm_spec)
    
full_model <- fit(full_lm_wf, data = cars2018) 

full_model %>% tidy()
```

a. Use `tidymodels` to perform 10-fold cross-validation to estimate test MAE for this model.

```{r}
cars_cv <- vfold_cv(cars2018, v = ???)

full_mod_cv = fit_resamples(????,
  resamples = ????, 
  metrics = metric_set(????)
)

CV_maes = full_mod_cv %>% unnest(.metrics) %>% select(.estimate) 
mean(CV_maes$.estimate)
```

b. How do you think the estimated test error would change with fewer predictors?

c. This model fit with ordinary least squares corresponds to a special case of penalized least squares. What is the value of $\lambda$ in this special case?

d. As $\lambda$ increases, what would you expect to happen to the number of predictors that remain in the model?    


<br><br>


### Exercise 2: Fitting a LASSO model in `tidymodels` {-}

a. Adapt our general LASSO code to fit a set of LASSO models with the following parameters:

- Use 10-fold CV.
- Use mean absolute error (MAE) to select a final model.
- Select the simplest model for which the metric is within one standard error of the best metric.
- Use a sequence of 30 $\lambda$ values from 0.001 to 1.

Before running the code, enter `install.packages("glmnet")` in the Console.

Save the CV-fit models from `tune_grid()` as `tune_output`.


```{r}
# Fit LASSO models for a grid of lambda values
# Tune and fit a LASSO model to the data (with CV)
set.seed(2023)

# Create CV folds
data_cv10 <- vfold_cv(??, v = 10)

# Lasso Model Spec with tune
lm_lasso_spec_tune <- 
  linear_reg() %>%
  set_args(mixture = 1, penalty = tune()) %>% ## mixture = 1 indicates Lasso
  set_engine(engine = 'glmnet') %>% #note we are using a different engine
  set_mode('regression') 

# Workflow (Recipe + Model)
lasso_wf_tune <- workflow() %>% 
  add_recipe(full_rec) %>% # recipe defined above
  add_model(lm_lasso_spec_tune) 

# Tune Model (trying a variety of values of Lambda penalty)
penalty_grid <- grid_regular(
  penalty(range = c(??, ??)), #log10 transformed (e.g., put '-3' for 0.001) 
  levels = ??)

tune_output <- tune_grid( # new function for tuning parameters
  lasso_wf_tune, # workflow
  resamples = data_cv10, # cv folds
  metrics = metric_set(rmse, mae),
  grid = penalty_grid # penalty grid defined above
)
```


b. Let’s visualize the model evaluation metrics from tuning. We can use `autoplot()`.

```{r}
# Visualize Model Evaluation Metrics from Tuning
autoplot(tune_output) + theme_classic()
```

c. Inspect the shape of the plot. The errors go down at the very beginning then start going back up. Why do you think this happens? (This is an example of a very important idea that we’ll see shortly: the bias-variance tradeoff.)
  

d. Next, we need to choose the lambda that leads to the best model. We can choose the lambda penalty value that leads to the lowest cross-validated MAE or we can take into account the variation of the cross-validated MAE and choose the largest lambda penalty value that is within 1 standard error of the lowest cross-validated MAE. How might the models that result from these two penalties differ?  
  
```{r}
best_penalty <- select_best(tune_output, metric = 'mae') # choose penalty value based on lowest cross-validated MAE
best_penalty

best_se_penalty <- select_by_one_std_err(tune_output, metric = 'mae', desc(penalty)) # choose largest penalty value within 1 se of the lowest cross-validated MAE
best_se_penalty
```

e. Now check your understanding by fitting both “final” models and comparing the coefficients. How are these two models different?

```{r}
# Fit Final Model

final_wf <- finalize_workflow(lasso_wf_tune, best_penalty) # incorporates penalty value to workflow
final_wf_se <- finalize_workflow(lasso_wf_tune, best_se_penalty) # incorporates penalty value to workflow

final_fit <- fit(final_wf, data = cars2018)
final_fit_se <- fit(final_wf_se, data = cars2018)

tidy(final_fit)
tidy(final_fit_se)
```

<br><br>

### Exercise 3: Examining output: plot of coefficient paths {-}

Once we’ve used cross validation, a useful plot allows us to examine coefficient paths resulting from the final fitted LASSO models: coefficient estimates as a function of $\lambda$.

Before running the code, run install.packages("stringr") and install.packages("purrr") in the Console.

```{r}
glmnet_output <- final_fit_se %>% extract_fit_parsnip() %>% pluck('fit') # get the original glmnet output

lambdas <- glmnet_output$lambda
coefs_lambdas <- 
  coefficients(glmnet_output, s = lambdas )  %>% 
  as.matrix() %>%  
  t() %>% 
  as.data.frame() %>% 
  mutate(lambda = lambdas ) %>% 
  select(lambda, everything(), -`(Intercept)`) %>% 
  pivot_longer(cols = -lambda, 
               names_to = "term", 
               values_to = "coef") %>%
  mutate(var = purrr::map_chr(stringr::str_split(term,"_"),~.[1]))

coefs_lambdas %>%
  ggplot(aes(x = lambda, y = coef, group = term, color = var)) +
  geom_line() +
  geom_vline(xintercept = best_se_penalty %>% pull(penalty), linetype = 'dashed') + 
  theme_classic() + 
  theme(legend.position = "bottom", legend.text=element_text(size=8))
```

There's a lot of information in this plot!

- Each colored line corresponds to a different predictor.
- The x-axis reflects the range of different $\lambda$ values.
- At each $\lambda$, the y-axis reflects the coefficient estimates for the predictors in the corresponding LASSO model.
- The vertical dashed line shows where the best penalty value (using the SE method) based on cross-validated MAE.

<br>

a. Why do all of the lines head toward y = 0 on the far right of the plot?

b. What variables seem to be more “important” or “persistent” (persistently present in the model) variable? Does this make sense in context?


<br><br>


### Exercise 4: Examining and evaluating the best LASSO model. {-}

a. Take a look at the predictors and coefficients for the “best” LASSO model. Are the predictors that remain in the model sensible? Do the coefficient signs make sense?

```{r}
# Obtain the predictors and coefficients of the "best" model
# Filter out the coefficient are 0
final_fit_se %>% tidy() %>% filter(estimate != 0)
```

b. Evaluate the best LASSO model:

- Contextually interpret (with units) the cross-validated MAE error for the best model by inspecting `tune_output %>% collect_metrics() %>% filter(penalty == (best_se_penalty %>% pull(penalty)))`.
- Make residual plots for the model by creating a dataset called `lasso_mod_out` which contains the original data as well as predicted values and residuals (`.pred` and `resid`).   
```{r}
lasso_mod_out <- final_fit_se %>%
    predict(new_data = cars2018) %>%
    bind_cols(cars2018) %>%
    mutate(resid = mpg - .pred)
```

**Note:** If you're curious about making plots that show both test error *estimates* and their *uncertainty*, look at Digging Deeper.


<br><br>


### Digging deeper {-}

We used the plot of coefficient paths to evaluate the variable importance of our predictors. The code below does this systematically for each predictor so that we don’t have to eyeball. Step through and work out what each part is doing. It may help to look up function documentation with ?function_name in the Console.

```{r}
glmnet_output <- final_fit_se %>% extract_fit_engine()
    
# Create a boolean matrix (predictors x lambdas) of variable exclusion
bool_predictor_exclude <- glmnet_output$beta==0

# Loop over each variable
var_imp <- sapply(seq_len(nrow(bool_predictor_exclude)), function(row) {
    this_coeff_path <- bool_predictor_exclude[row,]
    if(sum(this_coeff_path) == ncol(bool_predictor_exclude)){ return(0)}else{
    return(ncol(bool_predictor_exclude) - which.min(this_coeff_path) + 1)}
})

# Create a dataset of this information and sort
var_imp_data <- tibble(
    var_name = rownames(bool_predictor_exclude),
    var_imp = var_imp
)
var_imp_data %>% arrange(desc(var_imp))
```

If you want more practice, the `Hitters` data in the `ISLR` package (be sure to to install and load) contains the salaries and performance measures for 322 Major League Baseball players. Use LASSO to determine the “best” predictive model of `Salary`.
