---
title: "More Tidymodels"
subtitle: "Lecture 24"
author: "Dr. Colin Rundel"
footer: "Sta 323 - Spring 2025"
format:
  revealjs:
    theme: slides.scss
    transition: fade
    slide-number: true
    self-contained: true
execute: 
  echo: true
  message: true
---

```{r setup}
#| message: False
#| warning: False
#| include: False

knitr::opts_chunk$set(
  fig.align = "center", fig.retina = 2, dpi = 150,
  out.width = "100%"
)

library(tidymodels)
library(tidyverse)
library(patchwork)

ggplot2::theme_set(ggplot2::theme_bw())

options(width=85)
```

#

![](imgs/hex_tidymodels.png){fig-align="center" width="50%"}


## Hotels Data

Original data from [Antonio, Almeida, and Nunes (2019)](https://doi.org/10.1016/j.dib.2018.11.126), see data dictionary [here](https://github.com/rfordatascience/tidytuesday/tree/master/data/2020/2020-02-11#data-dictionary)

```{r}
hotels = read_csv(
  'https://tidymodels.org/start/case-study/hotels.csv'
) |>
  mutate(
    across(where(is.character), as.factor)
  )
```

::: {.aside}
This version of the data is slightly modified from the original data - see this [gist](https://gist.github.com/topepo/05a74916c343e57a71c51d6bc32a21ce) for the cleanup details
:::

## The data

```{r}
#| include: false
options(width=95)
```

::: {.small}
```{r}
glimpse(hotels)
```
:::

```{r}
#| include: false
options(width=60)
```

## The model

Our goal is to develop a predictive model that is able to predict whether a booking will include children or not based on the other characteristics of the booking.

. . .

```{r}
hotels |>
  count(children) |>
  mutate(prop = n/sum(n))
```


## Stratifying the test/train split

:::: {.columns .small}
::: {.column width='50%'}
```{r}
set.seed(123)

splits = initial_split(
  hotels, strata = children
)

hotel_train = training(splits)
hotel_test = testing(splits)
```

```{r}
dim(hotel_train)
dim(hotel_test)
```
:::

::: {.column width='50%' .fragment}
```{r}
hotel_train |>
  count(children) |>
  mutate(prop = n/sum(n))

hotel_test |>
  count(children) |>
  mutate(prop = n/sum(n))
```
:::
::::

## Logistic Regression model

::: {.xsmall}
```{r}
show_engines("logistic_reg")
```
:::

. . .

::: {.xsmall}
```{r}
lr_model = logistic_reg() |>
  set_engine("glm")
```

```{r}
translate(lr_model)
```
:::

## Recipe

::: {.small}
```{r}
holidays = c("AllSouls", "AshWednesday", "ChristmasEve", "Easter", 
              "ChristmasDay", "GoodFriday", "NewYearsDay", "PalmSunday")

lr_recipe = recipe(children ~ ., data = hotel_train) |> 
  step_date(arrival_date) |> 
  step_holiday(arrival_date, holidays = holidays) |> 
  step_rm(arrival_date) |> 
  step_rm(country) |>
  step_dummy(all_nominal_predictors()) |> 
  step_zv(all_predictors())
```
:::

##

```{r include=FALSE}
options(width=100)
options(tibble.width = 100)
```

::: {.small}
```{r}
lr_recipe |>
  prep() |>
  bake(new_data = hotel_train)
```
:::

```{r include=FALSE}
options(width=60)
options(tibble.width = 60)
```

## Workflow

::: {.small}
```{r}
( lr_work = workflow() |>
    add_model(lr_model) |>
    add_recipe(lr_recipe) 
)
```
:::

## Fit

::: {.small}
```{r}
( lr_fit = lr_work |>
    fit(data = hotel_train) )
```
:::

## Tidy

::: {.small}
```{r}
lr_fit |>
  broom::tidy()
```
:::


## Logistic regression predictions

:::: {.columns .xsmall}
::: {.column width='50%'}
```{r}
( lr_train_perf = lr_fit |>
    augment(new_data = hotel_train) |>
    select(children, starts_with(".pred")) )
```
:::

::: {.column width='50%' .fragment}
```{r}
( lr_test_perf = lr_fit |>
    augment(new_data = hotel_test) |>
    select(children, starts_with(".pred")) )
```
:::
::::


## Performance metrics (within-sample)

:::: {.columns .xsmall}
::: {.column width='50%'}
```{r}
conf_mat(lr_train_perf, children, .pred_class)

accuracy(lr_train_perf, children, .pred_class)

precision(lr_train_perf, children, .pred_class)
```
:::

::: {.column width='50%' .fragment}
```{r}
yardstick::roc_curve(lr_train_perf, children, .pred_children) |>
  autoplot()

roc_auc(lr_train_perf, children, .pred_children)
```
:::
::::

## Performance metrics (out-of-sample)

:::: {.columns .xsmall}
::: {.column width='50%'}
```{r}
conf_mat(lr_test_perf, children, .pred_class)

accuracy(lr_test_perf, children, .pred_class)

precision(lr_test_perf, children, .pred_class)
```
:::

::: {.column width='50%'}
```{r}
yardstick::roc_curve(lr_test_perf, children, .pred_children) |>
  autoplot()

roc_auc(lr_test_perf, children, .pred_children)
```
:::
::::


## Combining ROC curves

::: {.small}
```{r}
#| out-width: 75%
bind_rows(
  lr_train_perf |> mutate(name = "logistic - train"),
  lr_test_perf |> mutate(name = "logistic - test")
) |>
  group_by(name) |>
  yardstick::roc_curve(children, .pred_children) |>
  autoplot()
```
:::




# Lasso 

## Lasso Model

For this we will be using the `glmnet` package which supports fitting lasso, ridge and elastic net models. 

$$
\min _{\beta_0, \beta} \frac{1}{N} \sum_{i=1}^N w_i l\left(y_i, \beta_0+\beta^T x_i\right)+\lambda\left[(1-\alpha)\|\beta\|{ }_2^2 / 2+\alpha\|\beta\|_1\right],
$$

. . .


::: {.small}
```{r}
lasso_model = logistic_reg(penalty = tune(), mixture = 1) |>
  set_engine("glmnet")
```
:::

* `mixture` ($\alpha$) determines the type of model fit 
  
  * `1` for Lasso, 
  
  * `0` for Ridge,
  
  * other for elastic net.

* `penalty` ($\lambda$) is the penalty term for coefficient size.


##

::: {.small}
```{r}
lasso_model |> 
  hardhat::extract_parameter_set_dials() |>
  print()
```
:::

. . .

::: {.small}
```{r}
lasso_model |>
  translate()
```
:::



## Lasso Recipe

Lasso (and Ridge) models are sensitive to the scale of the model features, and so a standard approach is to normalize all features before fitting the model.

::: {.xsmall}
```{r}
lasso_recipe = lr_recipe |>
  step_normalize(all_predictors())
```
:::

. . .

::: {.xsmall}
```{r}
lasso_recipe |>
  prep() |>
  bake(new_data = hotel_train)
```
:::

## Lasso workflow

::: {.xsmall}
```{r}
( lasso_work = workflow() |>
    add_model(lasso_model) |>
    add_recipe(lasso_recipe)
)
```
:::


## v-folds for hyperparameter tuning

::: {.small}
```{r}
( hotel_vf = rsample::vfold_cv(hotel_train, v=5, strata = children) )
```
:::

## grid search

::: {.small}
```{r}
#| code-line-numbers: "|1|3|4-6|7|8"
( lasso_grid = lasso_work |>
    tune_grid(
      hotel_vf,
      grid = tibble(
        penalty = 10^seq(-4, -1, length.out = 10)
      ),
      control = control_grid(save_pred = TRUE),
      metrics = metric_set(roc_auc)
    )
)
```
:::

## Results

```{r include=FALSE}
options(width=50)
options(tibble.width = 50)
```

:::: {.columns .xsmall}
::: {.column width='50%'}
```{r}
lasso_grid |>
  collect_metrics()
```
:::

::: {.column width='50%' .fragment}
```{r}
lasso_grid |> 
  autoplot()
```
:::
::::


```{r include=FALSE}
options(width=60)
options(tibble.width = 60)
```

## "Best" models

::: {.small}
```{r}
lasso_grid |>
  show_best(metric = "roc_auc", n=5)
```
:::

## "Best" model

::: {.small}
```{r}
( lasso_best = lasso_grid |>
    collect_metrics() |>
    mutate(mean = round(mean, 2)) |>
    arrange(desc(mean), desc(penalty)) |>
    slice(1) )
```
:::

## Extracting predictions

Since we used `control_grid(save_pred = TRUE)` with `tune_grid()` we can recover the predictions for the out-of-sample values for each fold:

::: {.xsmall}
```{r}
( lasso_train_perf = lasso_grid |>
    collect_predictions(parameters = lasso_best) )
```
:::

##

::: {.xsmall}
```{r}
#| out-width: 40%
lasso_train_perf |>
  roc_curve(children, .pred_children) |>
  autoplot()

lasso_train_perf |>
  roc_auc(children, .pred_children)
```
:::


## Re-fitting

Typically with a tuned model we update the workflow (or model) with the optimal parameter values and then refit using the complete training data,

::: {.xsmall}
```{r}
lasso_work_tuned = finalize_workflow(
  lasso_work,
  lasso_best
)

( lasso_fit = lasso_work_tuned |>
    fit(data=hotel_train) )
```
:::


## Test Performance (out-of-sample)

:::: {.columns .xsmall}
::: {.column width='50%'}
```{r}
lasso_test_perf = lasso_fit |>
  augment(new_data = hotel_test) |>
  select(children, starts_with(".pred"))
```

```{r}
conf_mat(lasso_test_perf, children, .pred_class)

accuracy(lasso_test_perf, children, .pred_class)

precision(lasso_test_perf, children, .pred_class)
```
:::

::: {.column width='50%'}
```{r}
yardstick::roc_curve(lasso_test_perf, children, .pred_children) |>
  autoplot()

roc_auc(lasso_test_perf, children, .pred_children)
```
:::
::::

## Comparing models

::: {.xsmall}
```{r}
#| out-width: 75%
bind_rows(
  lr_test_perf |> mutate(name = "logistic - test"),
  lasso_test_perf |> mutate(name = "lasso - test")
) |>
  group_by(name) |>
  yardstick::roc_curve(children, .pred_children) |>
  autoplot()
```
:::


# Decision tree

## Decision tree models

::: {.small}
```{r}
show_engines("decision_tree")
```
:::

. . .

::: {.small}
```{r}
dt_model = decision_tree(
  tree_depth = tune(), 
  min_n = tune(),
  cost_complexity = tune()
) |> 
  set_engine("rpart") |> 
  set_mode("classification")
```
:::

## Recipe & workflow

We skip dummy coding in the recipe as it is not needed by rpart,

::: {.small}
```{r}
dt_recipe = recipe(children ~ ., data = hotel_train) |> 
  step_date(arrival_date) |> 
  step_holiday(arrival_date, holidays = holidays) |> 
  step_rm(arrival_date) |>
  step_rm(country)
```
:::

. . .

::: {.small}
```{r}
dt_work = workflow() |> 
  add_model(dt_model) |> 
  add_recipe(dt_recipe)
```
:::

## Tuning

:::: {.columns .xsmall}
::: {.column width='50%'}
```{r}
( dt_grid = grid_regular(
    cost_complexity(), 
    tree_depth(), 
    min_n(), 
    levels = 3
) )
```
:::

::: {.column width='50%' .fragment}
```{r}
doFuture::registerDoFuture()
future::plan(future::multisession, workers = 8)
```

```{r}
#| cache: true
dt_tune = dt_work |> 
  tune_grid(
    hotel_vf,
    grid = dt_grid,
    control = control_grid(save_pred = TRUE),
    metrics = metric_set(roc_auc)
  )
```
:::
::::

. . .

How many decision tree models were fit?


## Tuning results

::: {.small}
```{r}
dt_tune |>
  collect_metrics() |>
  arrange(desc(mean))
```
:::


## "Best" parameters

```{r include=FALSE}
options(width=50)
options(tibble.width = 50)
```

:::: {.columns .xsmall}
::: {.column width='50%'}
```{r}
dt_tune |> 
  show_best(metric = "roc_auc")
```
:::

::: {.column width='50%'}
```{r}
autoplot(dt_tune)
```
:::
::::


```{r include=FALSE}
options(width=60)
options(tibble.width = 60)
```

## Re-fitting

::: {.xsmall}
```{r}
(dt_best = dt_tune |>
  select_best(metric = "roc_auc"))
```
:::

. . .

::: {.xsmall}
```{r}
dt_work_tuned = finalize_workflow(
  dt_work,
  dt_best
)

( dt_fit = dt_work_tuned |>
    fit(data=hotel_train))
```
:::

## Model extraction

::: {.small}
```{r}
dt_fit |> 
  hardhat::extract_fit_engine() |> 
  plot()
```
:::

## Test Performance (out-of-sample)

:::: {.columns .xsmall}
::: {.column width='50%'}
```{r}
dt_test_perf = dt_fit |>
  augment(new_data = hotel_test) |>
  select(children, starts_with(".pred"))
```

```{r}
conf_mat(dt_test_perf, children, .pred_class)

accuracy(dt_test_perf, children, .pred_class)

precision(dt_test_perf, children, .pred_class)
```
:::

::: {.column width='50%'}
```{r}
yardstick::roc_curve(dt_test_perf, children, .pred_children) |>
  autoplot()

roc_auc(dt_test_perf, children, .pred_children)
```
:::
::::


## Comparing models

::: {.xsmall}
```{r}
#| out-width: 75%
bind_rows(
  lr_test_perf |> mutate(name = "logistic - test"),
  lasso_test_perf |> mutate(name = "lasso - test"),
  dt_test_perf |> mutate(name = "decision tree - test")
) |>
  group_by(name) |>
  yardstick::roc_curve(children, .pred_children) |>
  autoplot()
```
:::


# Random Forest

## Random forest models

::: {.small}
```{r}
show_engines("rand_forest")
```
:::

. . .

::: {.small}
```{r}
rf_model = rand_forest(mtry = tune(), min_n = tune(), trees = 100) |> 
  set_engine("ranger", num.threads = 8) |> 
  set_mode("classification")
```
:::

## Recipe & workflow

We skip dummy coding in the recipe as it is not needed by ranger,

::: {.small}
```{r}
rf_recipe = recipe(children ~ ., data = hotel_train) |> 
  step_date(arrival_date) |> 
  step_holiday(arrival_date, holidays = holidays) |> 
  step_rm(arrival_date) |>
  step_rm(country)
```
:::

. . .

::: {.small}
```{r}
rf_work = workflow() |> 
  add_model(rf_model) |> 
  add_recipe(rf_recipe)
```
:::


## Tuning - automatic grid search

:::: {.columns .xsmall}
::: {.column width='40%'}
```{r}
#| include: false
future::plan(future::multisession, workers = 1)
```

```{r}
#| cache: true
rf_tune = rf_work |> 
  tune_grid(
    hotel_vf,
    grid = 10,
    control = control_grid(save_pred = TRUE),
    metrics = metric_set(roc_auc)
  )
```
:::

::: {.column width='60%' .fragment}
```{r}
rf_tune |> 
  collect_metrics() |>
  arrange(desc(mean))
```
:::
::::


## "Best" parameters

```{r include=FALSE}
options(width=50)
options(tibble.width = 50)
```

:::: {.columns .xsmall}
::: {.column width='50%'}
```{r}
rf_tune |> 
  show_best(metric = "roc_auc")
```
:::

::: {.column width='50%'}
```{r}
autoplot(rf_tune)
```
:::
::::


```{r include=FALSE}
options(width=60)
options(tibble.width = 60)
```

## Re-fitting

::: {.xsmall}
```{r}
rf_best = rf_tune |>
  select_best(metric = "roc_auc")
```

```{r}
rf_work_tuned = finalize_workflow(
  rf_work, 
  rf_best
)

( rf_fit = rf_work_tuned |>
    fit(data=hotel_train) )
```
:::

## Test Performance (out-of-sample)

:::: {.columns .xsmall}
::: {.column width='50%'}
```{r}
rf_test_perf = rf_fit |>
  augment(new_data = hotel_test) |>
  select(children, starts_with(".pred"))
```

```{r}
conf_mat(rf_test_perf, children, .pred_class)

accuracy(rf_test_perf, children, .pred_class)

precision(rf_test_perf, children, .pred_class)
```
:::

::: {.column width='50%'}
```{r}
yardstick::roc_curve(rf_test_perf, children, .pred_children) |>
  autoplot()

roc_auc(rf_test_perf, children, .pred_children)
```
:::
::::


## Comparing models

::: {.xsmall}
```{r}
#| out-width: 75%
bind_rows(
  lr_test_perf |> mutate(name = "logistic - test"),
  lasso_test_perf |> mutate(name = "lasso - test"),
  dt_test_perf |> mutate(name = "decision tree - test"),
  rf_test_perf |> mutate(name = "random forest - test")
) |>
  group_by(name) |>
  yardstick::roc_curve(children, .pred_children) |>
  autoplot()
```
:::