---
title: "🗓️ Week 5 - Other Resampling Methods"
subtitle: "Theme: Supervised Learning"
author: "Dr. Jon Cardoso-Silva"
format: 
  html:
    from: markdown+emoji
    page-layout: full
    toc: true
    toc-depth: 2
self-contained: true
---

Welcome to Week 05!

Note that this is just a piece of the lecture's material. The first half was a live recap explanation on the whiteboard about logistic regression and classification metrics. To study those principles deeply, take a look at:

-   James, Gareth, Daniela Witten, Trevor Hastie, and Robert Tibshirani. 2021. [An Introduction to Statistical Learning: With Applications in R](https://www.statlearning.com/). 2nd edition. Springer Texts in Statistics. New York NY: Springer.
-   `parsnip::logistic_reg()` [documentation](https://parsnip.tidymodels.org/reference/logistic_reg.html)

# ⚙️ Setup

```{r, echo=FALSE, message=FALSE, warning=FALSE}
# Tidyverse packages
library(tidyr)
library(dplyr)
library(lubridate)
library(ggplot2)

# Tidymodels packages, spelled out individually for clarity
library(parsnip)
library(yardstick)
library(recipes)
library(rsample)
library(tune)
library(workflows)
```

# Part 1: Meet the dataset

**Let's use a different dataset**

To avoid the confusion that typically surrounds the lagged variables needed for forecasting, we will use a dataset that does not contain a temporal component. This will allow us to focus on the resampling methods themselves.

Here, I will use a fictional dataset created by IBM with "employee attrition" as the target variable, which we will interpret as indicating whether an employee intends to leave the company.

```{r}
library(modeldata)

# Run this to load the dataset into memory
# The dataset will be saved in a variable called `attrition`
data(attrition)

# I'd rather be safe than sorry
attrition$Attrition <- factor(attrition$Attrition, levels=c("No", "Yes"))
```

The dataset contains 1470 rows and 35 columns and is already tidy. Each row represents a single employee, and the columns are information about the employee. The target variable is `Attrition`, and all the rest are potential predictors.

🎯 **ACTION POINTS**

1.  **Explore the dataset**. What are the columns? What are the data types? Are there any missing values?

```{r}

```

2.  **Create a `recipe`**. Assume we want to use **all features** to predict `Attrition`. What are the steps you would take to create a recipe for this?

```{r}
rec <- 
    recipe(Attrition ~ ., data = attrition) %>% 
    step_dummy(all_nominal(), -all_outcomes()) %>% 
    step_zv(all_predictors()) %>% 
    step_nzv(all_predictors()) %>% 
    step_corr(all_predictors(), threshold = 0.9) %>%
    prep()

summary(rec)
```

3.  **Create a `workflow`**. Just run the code below to create a workflow that ties the recipe and the model together.

```{r}
log_reg_spec <- 
    logistic_reg() %>% 
    set_engine("glm") %>% 
    set_mode("classification")

log_rec_wf <- 
    workflow() %>% 
    add_recipe(rec) %>% 
    add_model(log_reg_spec)
```

# Part 2: A deeper look at training vs test splits

## Firstly, a train vs test split

Let's practice training vs testing splits again. We start like we've done before by splitting the data into training and test sets.

```{r}
dataset_split <- attrition %>% initial_split(prop=0.8)
```

Let me show you what's inside the resulting object:

```{r}
df_train <- training(dataset_split)
```

```{r}
df_test <- testing(dataset_split)
```

::: {style="width:50%;border: 1px solid #aaa; border-radius:1em; padding: 1em; margin: 1em 0;"}
📢 **THE CODE ABOVE IS THE SAME WE'VE BEEN DOING, BUT WE WILL THINK ABOUT THE SPLITS A BIT DIFFERENTLY HERE!**

This is because, although we've already done one split, we will still **further divide** the `df_train`.

Keep your eyes peeled to understand why!
:::

## Secondly, we *resample* the training data!

We will take the training data and further divide it into multiple training vs test sets, in a practice called **cross-validation**.

We need to specify the number of `v` **folds** (the number of splits), producing `v` disjoint test sets. We do like we were doing before: for each fold, we will reserve the test set and train the model on the remaining samples of the data. Then, at the end, we take a look at the average error across all of the folds.

💡 Sometimes, these test sets that are *internal* to cross-validation are also called 'validation sets'.

```{r}
# We further split df_train into v=10 folds
folds <- vfold_cv(df_train, v=10)

# Take a look at what this object looks like
folds
```

## Thirdly, we evaluate the model's performance under cross-validation

```{r}
# I would use just the ROC AUC here as an evaluating metric
# AUC is the metric that best indicates the OVERALL performance of the model (regardless of the threshold)
cv_results <- 
    log_rec_wf %>%
    fit_resamples(folds, metrics=metric_set(roc_auc))
```

Collect results:

```{r}
# Not bad! Finally, I found a good model in this course!! 🎉
# Well, of course, it is good. It's a fictional dataset!
cv_results %>% collect_metrics()
```

## Fourthly, we evaluate the model's performance on the test set

```{r}
# Fit the model to the training data
log_rec_fit <- 
    log_rec_wf %>%
    fit(data = df_train)

# Evaluate the model on the test set to see if it is very close to the cross-validation results
log_rec_fit %>% 
    augment(new_data = df_test) %>% 
    roc_auc(truth = Attrition, .pred_Yes, event_level="second")
```

⚠️ Note that I had to specify the event level to 'roc_auc' there. Otherwise, the metric will default to the first level of the factor! Inconsistencies... inconsistencies everywhere!

## Visually comparing the results

The first best thing to look at is the ROC AUC curve on the training set. This is a good way to see how well the model is doing across different thresholds.

```{r}
g <- 
    log_rec_fit %>% augment(new_data = df_train) %>%
    roc_curve(truth = Attrition, .pred_Yes, event_level="second") %>% 
    autoplot()

g <- g + 
  ## Theming. 
  ## Increase/decrease the font sizes according to your screen
  theme_bw(base_size=15) + 

  labs(title="A good ROC Curve on the Training Set",
       subtitle="This curve shows how the\n model performs under different thresholds",
       x="False Positive Rate\n(1 - Specificity)", y="True Positive Rate\n(Sensitivity)")

g
```

What about the test set?

```{r}
g <- 
    log_rec_fit %>% augment(new_data = df_test) %>%
    roc_curve(truth = Attrition, .pred_Yes, event_level="second") %>% 
    autoplot()

g <- g + 
  ## Theming. 
  ## Increase/decrease the font sizes according to your screen
  theme_bw(base_size=15) + 

  labs(title="Also good on the TESTING Set",
       subtitle="This curve shows how the\n model performs under different thresholds",
       x="False Positive Rate\n(1 - Specificity)", y="True Positive Rate\n(Sensitivity)")

g
```

::: {style="width:50%;border: 1px solid #aaa; border-radius:1em; padding: 1em; margin: 1em 0;"}
📝 **What can we learn from this?**

-   **Resampling by means of cross-validation** allows us to estimate the average performance of a model on unseen data.
-   We can confirm that the estimates were good by comparing the results of the test set.
    -   The Area Under the test set's curve was very close to the average AUC of the cross-validation.
-   
:::

# Part 3: 📉 Making decisions

The talk above is about model building, but what about **deploying** this model to make decisions?

Having established that this is a good model, we can ask:

## Is there an optimal threshold based on the ROC curve?

NOTE: only calculate optimal thresholds on a training set. The test set is for final evaluation only.

We can use a few helper functions from the `probably` package to help us with this. Here, I will frame it around a metric called [j_index](https://yardstick.tidymodels.org/reference/j_index.html), which measures the distance between the ROC curve and the diagonal 'dotted' line (see [figure on this page](https://www.wikiwand.com/en/Youden%27s_J_statistic)).

```{r}
library(probably) # also part of tidymodels nowadays

# 0.5 wins! Perfect!
log_rec_fit %>% 
    augment(new_data = df_train) %>% 
    threshold_perf(truth = Attrition, estimate = .pred_Yes, event_level="second") %>%
    pivot_wider(names_from=.metric, values_from=.estimate) %>% 
    arrange(desc(j_index))
```

## What does that mean in terms of the most common metrics?

See the full list of classification metrics on the [`yardstick` reference documentation page](https://yardstick.tidymodels.org/reference/index.html#classification-metrics). Here, I will show you a few of them.

| **Metric**  | **Description**                                                                                       | **`yardstick` function**                                                   |
|------------------|----------------------|--------------------------------|
| Accuracy    | Proportion of correct predictions                                                                     | [`accuracy()`](https://yardstick.tidymodels.org/reference/accuracy.html)   |
| Sensitivity | The proportion of positive results out of the number of samples that were actually positive.          | [`sens()`](https://yardstick.tidymodels.org/reference/sens.html)           |
| Specificity | The proportion of negatives that are correctly identified as negatives.                               | [`spec()`](https://yardstick.tidymodels.org/reference/spec.html)           |
| Precision   | The proportion of true positive results out of the number of samples that were predicted as positive. | [`precision()`](https://yardstick.tidymodels.org/reference/precision.html) |
| Recall      | The proportion of true positive results out of the number of samples that were actually positive.     | [`recall()`](https://yardstick.tidymodels.org/reference/recall.html)       |
| F-measure   | The harmonic mean of precision and recall.                                                            | [`f_meas()`](https://yardstick.tidymodels.org/reference/f_meas.html)       |

**Training set performance**

```{r}
# Identify the multiple metrics I care about
multi_metric <- metric_set(accuracy, sens, spec, precision, recall, f_meas)

log_rec_fit %>% 
    augment(new_data = df_train) %>%    
    multi_metric(truth = Attrition, estimate=.pred_class)
```

**Test set performance**

```{r}
log_rec_fit %>% 
    augment(new_data = df_test) %>%    
    multi_metric(truth = Attrition, estimate=.pred_class)
```

# Part 4: ⌛ Going back in time

Continuing from where we left off last week, we will now look at a different resampling method: **rolling window validation**.

When we are using ML methods to do forecasting, like the UK House Price Index, we need to be careful about how we split our data. One way to check this is by training **multiple models**, each for a separate period, then evaluating its performance on the following year/period (test).

In these situations, we can use the `rsample::rolling-window()` function to create multiple training vs test splits.

Back to the UK House Price Index dataset:

```{r}
uk_hpi <- 
    readr::read_csv("data/UK-HPI-full-file-2023-06.csv") %>%
    mutate(Date = dmy(Date)) %>%
    janitor::clean_names()

UK_countries <- c("England", "Wales", "Scotland", "Northern Ireland")

dataset <- 
    uk_hpi %>% 
    filter(region_name %in% UK_countries) %>%
    group_by(region_name) %>%
    arrange(date) %>%
    mutate(
      lag_1_average_price = lag(average_price, 1),
      lag_13_average_price = lag(average_price, 13),
      change_last_year = (lag_1_average_price/lag_13_average_price)
    ) %>%
    drop_na(lag_13_average_price) %>%
    mutate(
    sum_past_rate_change = cumsum(change_last_year) - change_last_year,
    n_past_rows = row_number() - 1,
    historical_avg_change = sum_past_rate_change/n_past_rows,
    price_up=factor(average_price > lag_1_average_price, 
                    labels=c("No", "Yes"),
                    levels=c(FALSE, TRUE),
                    ordered=TRUE)) %>%
    drop_na(historical_avg_change) %>%
    arrange(date, region_name) %>%
    select(date, region_name, 
           average_price, lag_1_average_price, lag_13_average_price, 
           change_last_year, 
          #  sum_past_rate_change, n_past_rows, # Uncomment to follow the calculations
           historical_avg_change, price_up)
```

## 4.1 Rolling window validation

Following a similar design as [this official tidymodels tutorial](https://www.tidymodels.org/learn/models/time-series/):

```{r}
# Training: 50 years of data (50*12 months)
# Testing:  the 12-month period after training 
roll_rs <- rolling_origin(dataset %>% group_by(region_name) %>% arrange(date),   
               initial = 4 * 12 * 50, # x4 countries x 12x50 months 
               assess = 4 * 12, # x4 countries x 12 months
               cumulative = FALSE)
```

What does this object look like?

There are a LOT of splits:

```{r}
# Also try roll_rs$splits to see each individual split
roll_rs %>% nrow()
```

The data was split into 149 different training vs test splits!!

**What did we do?**

The above was created **multiple samples** (a resampling of the original data) such that each sample (split) has:

-   an Analysis subset: the training data (50 years worth of data)
-   an Assess subset: the testing data (12 months worth of data)

The **training and test data are different for each split!** This is the key here.

## 4.2 Inspecting the splits

**How do I know what is in each split?**

```{r}
# Let me get the first split
roll_rs$splits[[1]]
```

This tells me there are 2400 rows in the analysis (training) set and 48 rows in the assessment (test) set.

I can look at what is in each of them by running `training(roll_rs$splits[[1]])` or `testing(roll_rs$splits[[1]])`.

Let me do a summary to understand the dates:

```{r}
# This is the training data
training(roll_rs$splits[[1]]) %>% summarise(min_date=min(date), max_date=max(date))
```

```{r}
# This is the test data
testing(roll_rs$splits[[1]]) %>% summarise(min_date=min(date), max_date=max(date))
```

Similarly, the second split will have training data from '06/1969 - 05/2019' and testing data from '06/2019 - 05/2020'. That is, we slide one month to the future.

## 4.3 Running the model on all splits

How do we test our baseline model for all data splits above? Instead of `fit()`, we would use `tune::fit_resamples()`.

```{r, message=FALSE, warning=FALSE}
# Define the recipe
rec <- 
    recipe(price_up ~ change_last_year, data = roll_rs$splits[[1]]) %>% 
    prep()

# Use a workflow to tie the recipe and the model together
log_rec_wf <- 
    workflow() %>% 
    add_recipe(rec) %>% 
    add_model(log_reg_spec)


# Use a workflow to tie the recipe and the model together
log_rec_wf <- 
    workflow() %>% 
    add_recipe(rec) %>% 
    add_model(log_reg_spec)

multi_metric <- metric_set(roc_auc, accuracy, sens, spec, precision, recall, f_meas)

# This takes a while to run
results <- tune::fit_resamples(log_rec_wf, roll_rs, 
                               metrics=multi_metric,
                               control=control_resamples(save_pred = TRUE, event_level="second"))
```

::: {style="width:50%;border: 1px solid #aaa; border-radius:1em; padding: 1em; margin: 1em 0;"}
📝 **Do you see a lot of errors?**

That's because the model is not good at all. We will see that in a moment.

IF WE HAVE TIME: fix the model above based on what we learned from the submissions to the W05 competition.
:::

You end up with a huge data frame with a lot of results per split:

```{r}
head(results)
```

Wanna see the result for a split, say (Slice100)?

```{r}
# See the AUC for a selected split
results$.metrics[[100]]
```

Luckily, tidymodels has a useful function to collect the metrics in a tidy way:

```{r}
collect_metrics(results)
```

# How to practice this

-   Tweak the UK HPI model, add more variables, and see if you can improve the results.
-   **Plot** the variables, colouring the dots by the `price_up` variable. This will give you a sense of which variables are good for predicting the target.