---
title: "🗓️ Week 4 - Metrics & 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 04!

We continue our saga of finding good models for the data we have. This week, we introduce new metrics and tools to assess your algorithms, but note that **we will show you bad models**. There is a competition next week for whoever can come up with the best model using resampling methods.

::: callout-important
## 📝 Errata

Guess who saw the future? 🔮 I did.

I made a mistake in the notebook shared last week, which resulted in a sneaky data leakage.

The model `rate_change_1m ~ rate_change_2m` was not a correct example for predicting the future. It defined `rate_change_2m` using the current ('future') value of `average_price` (the target variable), which is not feasible in reality since we don't have the future value of the target variable when predicting it.

The intended definition was:

``` r
rate_change_2m =  (lag_1_sales_volume / lag_2_sales_volume - 1)
```

Apologies for any confusion caused! A revised version of the notebook will be issued soon.
:::

# ⚙️ 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)    # This one is new today
library(tune)       # This one is new today
library(workflows)

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

# Part 1: Building a temporally-aware dataset

If you are revising this later, remember what we did in the lecture: instead of running the entire code below, break it down into smaller pieces and run each separately to understand what is happening.

Here, we will work with the assumption that the **average yearly increase** in house prices is a good predictor of whether the price will go up in the following month.

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

```

# Part 2: Resampling

## 2.1 Training vs Test sets with `rsample` package

So far, we have fixed the date at which we split the data into training and test sets. Now I want to show you how to do this in a more automated way using the `rsample` package.

::: callout-important
## ☢️ Warning

The code snippet below illustrates a process that works only when your data lacks a temporal component. For instance, if you were predicting cats vs dogs, you could employ the `initial_split()` function to divide your data into training and test sets.
:::

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

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

```{r}
training(dataset_split)
```

```{r}
testing(dataset_split)
```

Do you see why the above wouldn't work for our dataset?

## 2.2 When your data has a time-series dimension

As mentioned above, we must be aware of time when splitting our datasets. Thankfully, there is a function for that: `rsample::initial_time_split()` .

```{r}
dataset_split <- dataset %>% initial_time_split(prop=4/5)
```

Were the proportions respected?

```{r}
# Print dimensions of the training and test sets
cat(paste0("Our Training set has: ", training(dataset_split) %>% nrow(), " rows\n"))
cat(paste0("Our Testing set has : ", testing(dataset_split) %>% nrow(), " rows"))
```

What is the range of dates in both sets?

```{r}
training(dataset_split) %>% summarise(min_date=min(date), max_date=max(date))
```

```{r}
testing(dataset_split) %>% summarise(min_date=min(date), max_date=max(date))
```

**AVOID DATA LEAKAGE! Always check the range of dates in your training and test sets!**

## 2.3 Running a logistic regression model

Let's run a model on this data. We will use a logistic regression model, which is a simple model that is easy to understand and interpret. We will use the `tidymodels` package to do this.

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

# Define the recipe
rec <- 
    recipe(price_up ~ change_last_year, data = training(dataset_split)) %>% 
    prep()

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

# Fit the model
log_rec_fit <- log_rec_wf %>% fit(data = training(dataset_split)) %>% extract_fit_parsnip()
```

How should we express this model in equation form?

```{r}
log_rec_fit
```

yields:

``` text
parsnip model object


Call:  stats::glm(formula = ..y ~ ., family = stats::binomial, data = data)

Coefficients:
     (Intercept)  change_last_year  
         -1.5580            0.8841  

Degrees of Freedom: 2075 Total (i.e. Null);  2074 Residual
Null Deviance:      2706 
Residual Deviance: 2701         AIC: 2705
```

How should we express this model in equation form?

$$
\text{price_up} = \frac{e^{-1.5580 + 0.8841 \times \text{change_last_year}}}{1 + e^{-1.5580 + 0.8841 \times \text{change_last_year}}}
$$

## 2.4 Interpreting the model

First, let's go a bit deeper on logistic regression models. What do the regression coefficients mean?

**Intercept:** the intercept changes the [log of the odds ratio](https://www.statisticshowto.com/log-odds/) of the outcome when the predictor is zero.

- The probability of `price_up = "Yes"` when `change_last_year = 0` is:

$$
\frac{e^{-1.5580}}{1 + e^{-1.5580}} = 17.4\%
$$

 This coefficient indicates how much the log of the odds ratio changes for a one-unit increase in the predictor variable. In your example:

    When change_last_year increases by 1, the log odds of price_up = "Yes" increases by 0.8841.

**The coefficient of `change_last_year`** indicates how much the log of the odds ratio changes for a one-unit increase in the predictor variable. In your example:

- When `change_last_year` increases by 1, the probability of `price_up = "Yes"` increases by log-odds of 0.8841.
- The probability of `price_up = "Yes"` when `change_last_year = 1` is:

$$
\frac{e^{-1.5580 + 0.8841}}{1 + e^{-1.5580 + 0.8841}} = 70.3\%
$$

🎯 **ACTION POINTS**

Test it out! Ideally, in groups of 2-3, try to:

1. Create a plot of the probability of `price_up = "Yes"` as a function of `change_last_year`, using the equation produced by the model.

```{r}

```

## 2.5 Assessing the model

How well did this model perform on the test set?

```{r}
# Predictions
predictions <- log_rec_fit %>% augment(new_data = testing(dataset_split))

# Confusion matrix
conf_mat <- predictions %>% conf_mat(truth = price_up, estimate = .pred_class)

g <- conf_mat %>% autoplot(type = "heatmap") 

g <- g + 
  ## Customise scales
  scale_fill_gradient(low="#D6EAF8",high = "#2E86C1") +

  ## Theming. 
  ## Increase/decrease the font sizes according to your screen
  theme_bw(base_size=30)
g
```

What does that mean in terms of metrics?

**Accuracy**

Accuracy is the proportion of correct predictions over the total number of predictions. It is a good metric when the classes are balanced.

```{r}
accuracy(predictions, truth = price_up, estimate = .pred_class)
```

How do you check if classes are balanced?

```{r}
table(dataset$price_up)
```

**Precision or Positive Predictive Value (PPV)**

Precision is the proportion of true positive predictions over the total number of positive predictions. It is a good metric when the cost of false positives is high.

Note: what is considered **positive** is defined by the factor levels of the `price_up` variable. Here, we set `price_up = "Yes"` as the positive class with the parameter `event_level = "second"`.

```{r}
# event level defines what we mean by 'positive'
precision(predictions, truth = price_up, estimate = .pred_class, event_level = "second")
```

$$
\text{Precision} = \frac{\text{True Positives}}{\text{True Positives} + \text{False Positives}}
$$

🗣️ **CLASSROOM DISCUSSION:** Why is this model SO bad??

**Recall or Sensitivity**

Recall is the proportion of true positive predictions over the total number of actual positive cases. It is a good metric when the cost of false negatives is high.

$$
\text{Recall} = \frac{\text{True Positives}}{\text{True Positives} + \text{False Negatives}}
$$

```{r}
# event level defines what we mean by 'positive'
recall(predictions, truth = price_up, estimate = .pred_class, event_level = "second")
```

**F-score**

The F-score is the harmonic mean of precision and recall. It is a good metric when you want to balance precision and recall.

$$
F_1 = 2 \times \frac{\text{precision} \times \text{recall}}{\text{precision} + \text{recall}}
$$

```{r}
f_meas(predictions, truth = price_up, estimate = .pred_class, event_level = "second")
```

**ROC Curve**

The ROC curve is a plot of the true positive rate (recall) against the false positive rate. It is a good metric when you want to visualise the trade-off between precision and recall.

```{r}
roc_curve(predictions, truth = price_up, estimate = .pred_Yes, event_level="second") %>% autoplot()
```

**Area Under the Curve (AUC)**

The AUC is a metric that measures the performance of a model across all classification thresholds. It is a good metric when you want to compare models across different thresholds.

```{r}
roc_auc(predictions, truth = price_up, estimate=.pred_Yes, event_level="second")
```

Our model is worse than a random guess (AUC = 0.5). This is a very bad model.

# Part 3: One split is not enough

Is that a fluke? Or is this model really bad, no matter how we split the 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).

## 3.1 Rolling window validation

With data with a strong temporal component, like ours, we can use the `rsample::rolling-window()` function to create multiple training vs test splits. 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.

## 3.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 assess (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
analysis(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.

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

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

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

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

What does the above represent? The Area Under the Curve of this model is, on average 0.447 (+- 0.00899) on the test set. This shows that the bad performance we saw when we ran a single model was not a fluke: it is consistently bad across all splits.

This is definitely a more robust way of assessing your model: running it multiple times on different samples of our data. Because time is an important factor in our data, we used a rolling window validation, but there are other ways to perform resampling. See [rsample Common patterns page](https://rsample.tidymodels.org/articles/Common_Patterns.html) for more.

# Part 4: 🏆 Competition

**We created some awfully bad models. Can you do better?**

Using the same target variable I used in this notebook, can you come up with a set of features that do not contain any data leakage and can beat the awful performance of the one I showed today? That is, a model specification with a non-null F1-score for the class 'Yes' and an F1-score and accuracy as high as possible?

- Use the same logistic regression as the algorithm specification
- Use the same target variable
- Use the same rolling window resampling technique
- Feel free to play with any combination of features you can think of, and to even bring in additional data if you want.


Winners will be put in our Hall of Fame and receive DSI tote bags! 🎉 