---
title: "Lec 25 - More Tidymodels"
subtitle: "<br/> Statistical Programming"
author: "Sta 323 | Spring 2022"
date: "<br/> Dr. Colin Rundel"
output:
  xaringan::moon_reader:
    css: ["slides.css"]
    lib_dir: libs
    nature:
      highlightStyle: solarized-light
      highlightLines: true
      countIncrementalSlides: false
      ratio: "16:9"
---
exclude: true

```{r, message=FALSE, warning=FALSE, include=FALSE}
options(
  htmltools.dir.version = FALSE, # for blogdown
  width = 80,
  tibble.width = 80
)

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

options(width=60)
options(tibble.width = 60)

library(tidyverse)
library(tidymodels)
```


---
background-image: url(imgs/hex_tidymodels.png)
background-position: center
background-size: contain

---

## Hotels Data

Original data from [Antonio, Almeida, and Nunes (2019)](https://doi.org/10.1016/j.dib.2018.11.126), [Data dictionary](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)
  )
```

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

---

## The data

```{r}
glimpse(hotels)
```

---

## 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))
```

---

## Clustering the test/train split

.pull-left[
```{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)
```
]

--

.pull-right[
```{r}
hotel_train %>%
  count(children) %>%
  mutate(prop = n/sum(n))

hotel_test %>%
  count(children) %>%
  mutate(prop = n/sum(n))
```
]

---

## Logistic Regression model

```{r}
show_engines("logistic_reg")
```

--

```{r}
lr_model = logistic_reg() %>%
  set_engine("glm")

lr_model %>%
  translate()
```

---

## Recipe

```{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())

lr_recipe
```

---

```{r}
lr_recipe %>%
  prep() %>%
  bake(new_data = hotel_train)
```

---

## Workflow

```{r}
lr_work = workflow() %>%
  add_model(lr_model) %>%
  add_recipe(lr_recipe)
```

---

## Fit

```{r}
lr_fit = lr_work %>%
  fit(data = hotel_train)

lr_fit
```

---

## Logistic regression predictions

```{r}
lr_perf = lr_fit %>%
  augment(new_data = hotel_train) %>%
  select(children, starts_with(".pred"))

lr_perf
```

---

## Performance metrics (within-sample)

.pull-left[
```{r}
lr_perf %>%
  conf_mat(children, .pred_class)

lr_perf %>%
  precision(children, .pred_class)

lr_perf %>%
  roc_auc(children, .pred_children)
```
]

.pull-right[
```{r}
lr_perf %>%
  yardstick::roc_curve(
    children,
    .pred_children
  ) %>%
  autoplot()
```
]

---

## Performance metrics (out-of-sample)

.pull-left[
```{r}
lr_test_perf = lr_fit %>%
  augment(new_data = hotel_test) %>%
  select(children, starts_with(".pred"))

lr_test_perf %>%
  conf_mat(children, .pred_class)

lr_test_perf %>%
  precision(children, .pred_class)

lr_test_perf %>%
  roc_auc(children, .pred_children)
```
]

.pull-right[
```{r}
lr_test_perf %>%
  yardstick::roc_curve(
    children,
    .pred_children
  ) %>%
  autoplot()
```
]

---

## Combining ROC curves

.pull-left[
```{r}
lr_test_perf %>%
  yardstick::roc_curve(
    children,
    .pred_children
  ) 
```

```{r}
lr_roc_train = lr_perf %>%
  yardstick::roc_curve(children, .pred_children) %>% 
  mutate(name="logistic train")

lr_roc_test = lr_test_perf %>%
  yardstick::roc_curve(children, .pred_children) %>% 
  mutate(name="logistic test")
```
]

--

.pull-right[
```{r}
bind_rows(
  lr_roc_train,
  lr_roc_test
) %>%
  ggplot(aes(x = 1 - specificity, y = sensitivity, col = name)) + 
    geom_path(lwd = 1.5, alpha = 0.8) +
    geom_abline(lty = 3) + 
    coord_equal()
```
]

---
class: center, middle

## Lasso 

---

## Lasso Model

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

The mixture argument determines the type of model fit with: `1` -> Lasso, `0` -> Ridge, other -> elastic net.

```{r}
lasso_model = logistic_reg(penalty = tune(), mixture = 1) %>%
  set_engine("glmnet")

lasso_model %>%
  translate()

lasso_model %>% 
  parameters()
```

---

## 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.

```{r}
lasso_recipe = lr_recipe %>%
  step_normalize(all_predictors())
```

```{r}
lasso_recipe %>%
  prep() %>%
  bake(new_data = hotel_train)
```

---

## Lasso workflow

```{r}
lasso_work = workflow() %>%
  add_model(lasso_model) %>%
  add_recipe(lasso_recipe)
```


---

## k-Folds for tuning

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

---

## grid search

```{r}
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)
  )

lasso_grid
```

---

## Results

.pull-left[
```{r}
lasso_grid %>%
  collect_metrics()
```
]

--

.pull-right[
```{r}
lasso_grid %>% 
  collect_metrics() %>% 
  ggplot(aes(x = penalty, y = mean)) + 
    geom_point() + 
    geom_line() + 
    ylab("Area under the ROC Curve") +
    scale_x_log10(labels = scales::label_number())
```
]


---

## "Best" models

```{r}
lasso_grid %>%
  show_best("roc_auc", n=10)
```

---

## "Best" model

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

lasso_best
```

---

## 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:

```{r}
lasso_train_perf = lasso_grid %>%
  collect_predictions(parameters = lasso_best)
lasso_train_perf
```

--

```{r}
lasso_train_perf %>%
  roc_auc(children, .pred_children)
```

---

## Re-fitting

Typically with a tuned model we will refit using the complete test data and the "best" parameter value(s),

```{r}
lasso_work_tuned = update_model(
  lasso_work, 
  logistic_reg(
    mixture = 1, 
    penalty = lasso_best$penalty
  ) %>%
    set_engine("glmnet")
)

lasso_fit = lasso_work_tuned %>%
  fit(data=hotel_train)
```


---

## Test Performance (out-of-sample)

.pull-left[
```{r}
lasso_test_perf = lasso_fit %>%
  augment(new_data = hotel_test) %>%
  select(children, starts_with(".pred"))

lasso_test_perf %>%
  conf_mat(children, .pred_class)

lasso_test_perf %>%
  precision(children, .pred_class)

lasso_test_perf %>%
  roc_auc(children, .pred_children)
```
]

.pull-right[
```{r}
lasso_roc = lasso_test_perf %>%
  yardstick::roc_curve(
    children,
    .pred_children
  ) %>%
  mutate(name = "lasso - test")
lasso_roc %>%
  autoplot()
```
]

---

## Comparing models

```{r out.width="55%", echo=FALSE }
bind_rows(
  lr_roc_test,
  lasso_roc
) %>%
  ggplot(aes(x = 1 - specificity, y = sensitivity, col = name)) + 
    geom_path(lwd = 1.5, alpha = 0.8) +
    geom_abline(lty = 3) + 
    coord_equal()
```

---
class: center, middle

## Random Forest

---

## Random forest models

```{r}
show_engines("rand_forest")
```

--

```{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,

```{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)
```

--

```{r}
rf_work = workflow() %>% 
  add_model(rf_model) %>% 
  add_recipe(rf_recipe)
```

---

## Tuning

.pull-left[
```{r}
rf_work %>%
  parameters()
```
]

.pull-right[
```{r}
rf_grid = rf_work %>% 
  tune_grid(
    hotel_vf,
    grid = 10,
    control = control_grid(save_pred = TRUE),
    metrics = metric_set(roc_auc)
  )
```
]


---

## "Best" parameters

.pull-left[
```{r}
rf_grid %>% 
  show_best(metric = "roc_auc")
```
]

.pull-right[
```{r}
autoplot(rf_grid)
```
]

---

## Refitting

```{r}
(rf_best = rf_grid %>%
  select_best(metric = "roc_auc"))
```

```{r}
rf_work_tuned = update_model(
  rf_work, 
  rand_forest(
    trees=100,
    mtry = rf_best$mtry, 
    min_n = rf_best$min_n
  ) %>%
    set_engine("ranger", num.threads = 8) %>%
    set_mode("classification")
)

rf_fit = rf_work_tuned %>%
  fit(data=hotel_train)
```


---

## Test Performance (out-of-sample)

.pull-left[
```{r}
rf_test_perf = rf_fit %>%
  augment(new_data = hotel_test) %>%
  select(children, starts_with(".pred"))

rf_test_perf %>%
  conf_mat(children, .pred_class)

rf_test_perf %>%
  precision(children, .pred_class)

rf_test_perf %>%
  roc_auc(children, .pred_children)
```
]

.pull-right[
```{r}
rf_roc = rf_test_perf %>%
  yardstick::roc_curve(
    children,
    .pred_children
  ) %>%
  mutate(name = "RF - test")
rf_roc %>%
  autoplot()
```
]

---

## Comparing models

```{r out.width="55%", echo=FALSE }
bind_rows(
  lr_roc_test,
  lasso_roc,
  rf_roc
) %>%
  ggplot(aes(x = 1 - specificity, y = sensitivity, col = name)) + 
    geom_path(lwd = 1.5, alpha = 0.8) +
    geom_abline(lty = 3) + 
    coord_equal()
```