---
title: "🗓️ Week 5 - 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 05!

```{r, echo=FALSE, message=FALSE, warning=FALSE}

library(dplyr)
library(lubridate)
library(ggplot2)
library(readr)
library(tidymodels)

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

# Part 1: 🥇 Week 05 - Competition

Let's go over the solution that won the mini-competition, plus an interesting approach that will conveniently help us understand more complex aspects of feature selection and something called grid search!  

## 🥇 The winning solution: Bosco Hung

Let's take a look at the winning solution first. Bosco simply used the same baseline model but selected a threshold from a list of 100 different values. 

<details><summary> Click here to see the pre-processing</summary>

```{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_12_average_price = lag(average_price, 12),
      yearly_rate_increase = (average_price/lag_12_average_price) - 1
    ) %>%
    drop_na(lag_12_average_price) %>%
    mutate(
    sum_past_rate_increase = cumsum(yearly_rate_increase) - yearly_rate_increase,
    n_past_rows = row_number() - 1,
    historical_avg_increase = sum_past_rate_increase/n_past_rows,
    price_up=factor(yearly_rate_increase > historical_avg_increase, 
                    labels=c("No","Yes"),
                    levels=c(FALSE, TRUE),
                    ordered=TRUE)) %>%
    drop_na(historical_avg_increase) %>%
    arrange(desc(date), region_name) %>%
    select(-c(sum_past_rate_increase, n_past_rows))

dataset_train <- dataset %>% filter(year(date) <= 2018)
dataset_test <- dataset %>% filter(year(date) > 2018)
```

</details>

<details><summary>Click here to see the model</summary>

```{r}

log_rec_specification <- 
    logistic_reg() %>% 
    set_engine("glm") %>% 
    set_mode("classification") 

rec <- 
    recipe(price_up ~ ., data = dataset_train) %>% 
    update_role(-c(price_up), new_role = "ID") %>%
    step_lag(c(yearly_rate_increase, historical_avg_increase), lag = 1) %>%
    step_naomit(all_predictors(), skip=FALSE) %>%
    step_normalize() %>% # Bosco also added this step here, but normalisation doesn't change results in logistic regression.
    prep()

winning_model <-
    workflow() %>%
    add_recipe(rec) %>%
    add_model(log_rec_specification) %>% 
    fit(data = dataset_train) %>% 
    extract_fit_parsnip()

```

</details>

**What were the variables used?**

```{r}
summary(winning_model$fit)
```
Out of curiosity, let's plot them:

```{r echo=FALSE, fig.width=8, fig.height=5}
# Remember to bake the dataset!
plot_df <- 
    rec %>% bake(dataset_train) %>% 
    select(lag_1_historical_avg_increase, lag_1_yearly_rate_increase, price_up)

ggplot(plot_df, aes(x=lag_1_yearly_rate_increase, y=lag_1_historical_avg_increase)) +
  geom_bin_2d() +
  scale_fill_viridis_c(name="Frequency", option="E") +
  scale_x_continuous(name="% Yearly rate increase \n(lag: 1 month)", labels=scales::percent) +
  scale_y_continuous(name="Average Historical price increase \n(lag: 1 month)", labels=scales::percent) +
  theme_bw() +
  labs(title="How do the two variables relate to each other?",
       subtitle="Here I am ignoring the outcome variable. `price_up` is not represented here.") +
  theme(axis.text=element_text(size=rel(1.2)),
        axis.title = element_text(size=rel(1.3)),
        plot.title = element_text(size=rel(1.5)))

```

```{r echo=FALSE, fig.width=8, fig.height=6}
# Remember to bake the dataset!
plot_df <- 
    rec %>% bake(dataset_train) %>% 
    select(lag_1_historical_avg_increase, lag_1_yearly_rate_increase, price_up)

ggplot(plot_df, aes(x=lag_1_yearly_rate_increase, y=lag_1_historical_avg_increase)) +
  geom_bin_2d() +
  scale_fill_viridis_c(name="Frequency", option="E") +
  scale_x_continuous(name="% Yearly rate increase\n (as measured in previous month)", labels=scales::percent) +
  scale_y_continuous(name="Average Historical price increase \n(as measured in previous month)", labels=scales::percent) +
  theme_bw() +
  labs(title="Does the joint distribution change when priceUp = Yes vs No?",
       subtitle="Values in X-axis are higher when priceUp=Yes") +
  theme(axis.text=element_text(size=rel(1.2)),
        axis.title = element_text(size=rel(1.3)),
        plot.title = element_text(size=rel(1.5))) +
  facet_grid(price_up ~ .)

```

**NOTE:** What do we learn from plotting the two distributions together? Because the distributions are different when `priceUp = Yes` versus when `priceUp = No`, this reinforces our suspicion that these variables are somewhat useful in a model. **Plotting is a useful way to discover variables to include in a model**.

---

In equation form:

$$
\begin{align}
  Pr(\text{PriceUP}_t = \text{Yes}~|~\text{X1}_{t-1}, \text{X2}_{t-1}) &= \frac{e^{\beta_0 + \beta_1{\text{X1}}_{t-1} + \beta_2{\text{X2}}_{t-1}}}{1 + e^{\beta_0 + \beta_1\text{X1} + \beta_2\text{X2}}} \\
  \\
  &= \frac{e^{-0.7977 + 14.1426\text{X1}_{t-1} -12.9114\text{X2}_{t-1}}}{1 + e^{-0.7977 + 14.1426\text{X1}_{t-1} -12.9114\text{X2}_{t-1}}}

\end{align}

$$

where X1 represents the `yearly_rate_increase` and X2, `historical_avg_increase`.

**What does that formula represent?**

Many things:

- In the imaginary scenario in which both predictors are zero, the probability of the price going up in a year is $\frac{e^{-0.7977}}{1 + e^{-0.7977}}=\approx 31.05\%$. 
- Because $\hat{\beta_1}$ is positive, we know that any increase in $X1_{t-1}$ (yearly rate increase with a lag of 1) will lead to an increase in the probability of `priceUp`.
- By how much? In the imaginary scenario above, if $X1_{t-1}$ changed from 0 to $X1_{t-1} = 1$ and, **importantly, $X2_{t-1}$ remains zero**, we could expect the probability to climb to $\frac{e^{-0.7977  + 14.1426}}{1 + e^{-0.7977 + 14.1426}} = \approx 99.99%$! 
- Conversely, note that because $\hat{\beta_2}$ is negative, we know that the probability of `priceUp` will decrease when $X2_{t-1}$ (avg historical price increase with a lag of 1) increase.
- Note that, because both variables are positively linearly correlated and the regression coefficients attributed to them are opposite (one is positive and the other, negative), we are less likely to see huge jumps in the probability values.
  - Test this correlation by running:
  
    ```r
    rec %>% 
      bake(dataset_train) %>% 
      select(lag_1_historical_avg_increase, lag_1_yearly_rate_increase) %>% 
      cor()
    ```


### 📈 Another way to visualise the goodness-of-fit of this solution

```{r echo=FALSE, message=FALSE}
winning_model %>% 
    augment(rec %>% bake(dataset_train)) %>%
    roc_curve(truth=price_up, .pred_Yes, event_level="second") %>%
    autoplot() + 
    geom_point(aes(x=1-specificity, y=sensitivity, color=.threshold)) + 
    scale_color_gradient(name="Threshold", low = "#c6733c", high="#3cc6b8", limits=c(0, 1)) + 
    labs(title="Meet the ROC curve",
         subtitle="This model performs better than a random guess!", 
         x="(1 - specificity) = 1 - TN/N",
         y="(sensitivity) = TP/P")
```
There is a way to summarise the plot above in a single number, by calculating the Area Under the Curve:

```{r}
winning_model %>% 
    augment(rec %>% bake(dataset_train)) %>%
    roc_auc(truth=price_up, .pred_Yes, event_level="second") 
```


The ROC curve shows you what would happen to your model if you chose different thresholds. **The decision about which threshold is best is a 'business decision'**, it depends on what you want to do with this model afterwards. As I have been emphasizing in this course: _what do we want to do with this model? what metric is more relevant for my business model/research design?_

```{r message=FALSE}
# After looking at the ROC curve, 
#   this is roughly the threshold that represents the 'elbow' point in the curve
my_threshold <- 0.225

winning_model %>% 
    augment(rec %>% bake(dataset_train)) %>%
    mutate(.pred_class = .pred_Yes > my_threshold,
       .pred_class = factor(.pred_class, 
                            labels=c("No","Yes"), 
                            levels=c(FALSE, TRUE), 
                            ordered=TRUE)) %>%
    conf_mat(truth=price_up, estimate=.pred_class) %>%
    autoplot(type="heatmap") +
    scale_fill_gradient(low="#D6EAF8",high = "#2E86C1")
```

**Metrics-wise, what does the confusion matrix above mean?**

```{r, echo=FALSE, message=FALSE}
winning_model %>% 
    augment(rec %>% bake(dataset_train)) %>%
    mutate(.pred_class = .pred_Yes > my_threshold,
       .pred_class = factor(.pred_class, 
                            labels=c("No","Yes"), 
                            levels=c(FALSE, TRUE), 
                            ordered=TRUE)) %>%
    conf_mat(truth=price_up, estimate=.pred_class) %>%
    summary(estimator="binary", event_level="second") %>%
    knitr::kable()
```

💡 If you want to understand metrics more, check out Chapter 4 of the [ISLRv2 textbook](https://www.statlearning.com/), which is freely available online. But, honestly, even [Wikipedia](https://www.wikiwand.com/en/Confusion_matrix) has a good summary of classification metrics that stem out of a confusion matrix.

```{r}
winning_model %>% 
    augment(rec %>% bake(dataset_test)) %>%
    roc_curve(truth=price_up, .pred_Yes, event_level="second") %>%
    autoplot() +
    geom_point(aes(x=1-specificity, y=sensitivity, color=.threshold)) + 
    scale_color_gradient(low = "#c6733c", high="#3cc6b8", limits=c(0, 1)) + 
    labs(title="Ops...! Not the exact same story on the test set", 
         x="(1 - specificity) = 1 - TN/N",
         y="(sensitivity) = TP/P")
```
Here we see the same story as with when we were playing with training and testing sets for linear regression. Any model is more likely to fit the training data better than the test set. No surprises here.


## 🥈 Interesting solution: Louis

One of the submissions was very clever in that the student used many many variables in the hopes that putting them all together would lead to a better model. Let's see what happened!

<details><summary>Click here to see the pre-processing steps</summary>

```{r}
dataset <- 
    uk_hpi %>% 
    filter(region_name %in% UK_countries) %>%
    group_by(region_name) %>%
    arrange(date) %>%
    mutate(
      lag_12_average_price = lag(average_price, 12),
      yearly_rate_increase = (average_price/lag_12_average_price) - 1, 
      lag_12_detached_price = lag(detached_price, 12), 
      yearly_rate_detached_increase = (detached_price/lag_12_detached_price) - 1, 
      lag_12_semi_detached_price = lag(semi_detached_price, 12), 
      yearly_rate_semi_detached_increase = (semi_detached_price/lag_12_semi_detached_price) - 1,
      lag_12_terraced_price = lag(terraced_price, 12),
      yearly_rate_terraced_increase = (terraced_price/lag_12_terraced_price) - 1,
      lag_12_flat_price = lag(flat_price, 12), 
      yearly_rate_flat_increase = (flat_price/lag_12_flat_price) - 1,
      lag_12_old_price = lag(old_price, 12), 
      yearly_rate_old_increase = (old_price/lag_12_old_price) - 1
    ) %>%
    drop_na(lag_12_average_price, 
            lag_12_detached_price, 
            lag_12_semi_detached_price,  
            lag_12_terraced_price, 
            lag_12_flat_price, 
            lag_12_old_price) %>%
    mutate(
      sum_past_rate_increase = cumsum(yearly_rate_increase) - yearly_rate_increase,  # Cumsum = cumulative sum
      sum_past_rate_detached_increase = cumsum(yearly_rate_detached_increase) - yearly_rate_detached_increase, 
      sum_past_rate_semi_detached_increase = cumsum(yearly_rate_semi_detached_increase) - yearly_rate_semi_detached_increase, 
      sum_past_rate_terraced_increase = cumsum(yearly_rate_terraced_increase) - yearly_rate_terraced_increase,
      sum_past_rate_flat_increase = cumsum(yearly_rate_flat_increase) - yearly_rate_flat_increase, 
      sum_past_rate_old_increase = cumsum(yearly_rate_old_increase) - yearly_rate_old_increase,
      n_past_rows = row_number() - 1,  # ROW_NUMBER function returns the row number over a named or unnamed window specification.
      historical_avg_increase = sum_past_rate_increase/n_past_rows,
      historical_avg_detached_increase = sum_past_rate_detached_increase/n_past_rows,
      historical_avg_semi_detached_increase = sum_past_rate_semi_detached_increase/n_past_rows,
      historical_avg_terraced_increase = sum_past_rate_terraced_increase/n_past_rows,
      historical_avg_flat_increase = sum_past_rate_flat_increase/n_past_rows,
      historical_avg_old_increase = sum_past_rate_old_increase/n_past_rows,
      price_up=factor(yearly_rate_increase > historical_avg_increase, 
                      labels=c("No","Yes"),
                      levels=c(FALSE, TRUE),
                      ordered=FALSE)) %>%
    drop_na(historical_avg_increase, historical_avg_detached_increase, historical_avg_semi_detached_increase, historical_avg_terraced_increase, historical_avg_flat_increase, historical_avg_old_increase) %>%
    arrange(desc(date), region_name) %>%
    select(-c(sum_past_rate_increase, n_past_rows))

dataset_train <- dataset %>% filter(year(date) <= 2018)
dataset_test <- dataset %>% filter(year(date) > 2018)

```

</details>

<details><summary>Click here to see the model</summary>

```{r message=FALSE}
log_rec_specification <- 
    logistic_reg() %>% 
    set_engine("glm") %>% 
    set_mode("classification") 

rec_louis <- 
    recipe(price_up ~ ., data = dataset_train) %>% 
    update_role(-c(price_up), new_role = "ID") %>%
    step_lag(c(yearly_rate_increase, 
               historical_avg_increase, 
               yearly_rate_detached_increase, 
               historical_avg_detached_increase, 
               yearly_rate_semi_detached_increase, 
               historical_avg_semi_detached_increase,
               yearly_rate_terraced_increase, 
               historical_avg_terraced_increase, 
               yearly_rate_flat_increase, 
               historical_avg_flat_increase, 
               yearly_rate_old_increase, 
               historical_avg_old_increase), lag = 1:12) %>%
    step_normalize() %>%
    step_naomit(all_predictors(), skip=FALSE) %>%
    prep()

louis_model <-
    workflow() %>%
    add_recipe(rec_louis) %>%
    add_model(log_rec_specification) %>% 
    fit(data = dataset_train) %>% 
    extract_fit_parsnip()

```

</details>

The warning you get above should make you feel suspicious about this model, as it implies that our model was so good it fitted the training data perfectly! 

**How does this model look like?**

```{r}
summary(louis_model$fit)
```

First of all, it's a very **complex model**. Interpreting it is very tricky and, if you are very religious about p-values, all predictors are deemed statistically significant. 


```{r}
louis_model %>% 
    augment(rec_louis %>% bake(dataset_train)) %>%
    roc_curve(truth=price_up, .pred_Yes, event_level="second") %>%
    autoplot() + 
    geom_point(aes(x=1-specificity, y=sensitivity, color=.threshold)) + 
    scale_color_gradient(name="Threshold", low = "#c6733c", high="#3cc6b8", limits=c(0, 1)) + 
    labs(title="A near-perfect model on the training set!", 
         x="(1 - specificity) = 1 - TN/N",
         y="(sensitivity) = TP/P")
```

```{r}
louis_model %>% 
    augment(rec_louis %>% bake(dataset_train)) %>%
    roc_auc(truth=price_up, .pred_Yes, event_level="second") 
```
This model is nearly perfect!


### How good is this model on the test set?

```{r message=FALSE}
louis_model %>% 
    augment(rec_louis %>% bake(dataset_test)) %>%
    conf_mat(truth=price_up, estimate=.pred_class) %>%
    autoplot(type="heatmap") +
    scale_fill_gradient(low="#D6EAF8",high = "#2E86C1") +
    labs(title="It doesn't look all that perfect on the test set though")
```


```{r}
louis_model %>% 
    augment(rec_louis %>% bake(dataset_test)) %>%
    roc_curve(truth=price_up, .pred_Yes, event_level="second") %>%
    autoplot() + 
    geom_point(aes(x=1-specificity, y=sensitivity, color=.threshold)) + 
    scale_color_gradient(name="Threshold", low = "#c6733c", high="#3cc6b8", limits=c(0, 1)) + 
    labs(title="It's also awkward on the test set. What's going on?", 
         x="(1 - specificity) = 1 - TN/N",
         y="(sensitivity) = TP/P")
```

```{r}
louis_model %>% 
    augment(rec_louis %>% bake(dataset_test)) %>%
    roc_auc(truth=price_up, .pred_Yes, event_level="second") 
```
But it's not bad either! It's comparable to the performance of model1. 

**These examples show us two different approaches to modelling.** One could start from domain expertise, exploratory analysis, plotting and build a model by introducing more variables that seem important in isolation (read more about Feature Selection on Chapter 6 of [ISLRv2 textbook](https://www.statlearning.com/)). 

Alternatively, one could throw everything at the model and let it figure out what's most important! If this doesn't sound all that robust to you, that's because it isn't. But there are ways to make use of this strategy in a more principled way.

We will come back to this model in Part 3 of this tutorial! But first, I need you to understand the idea of **resampling**.

# Part 2: Resampling

## Performance of the model per year (rolling-window resampling)

So far, we have fixed the training and test sets, trained the algorithm once and evaluated it once. But how confident, really, can we be with our models? What if they're good just on that particular combination of train vs test years?

One way to solve this is by training **multiple models**, each for a separate period of time and then evaluate its performance on the following year/period (test).

<details><summary> Click here to see the pre-processing again</summary>

```{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_12_average_price = lag(average_price, 12),
      yearly_rate_increase = (average_price/lag_12_average_price) - 1
    ) %>%
    drop_na(lag_12_average_price) %>%
    mutate(
    sum_past_rate_increase = cumsum(yearly_rate_increase) - yearly_rate_increase,
    n_past_rows = row_number() - 1,
    historical_avg_increase = sum_past_rate_increase/n_past_rows,
    price_up=factor(yearly_rate_increase > historical_avg_increase, 
                    labels=c("No","Yes"),
                    levels=c(FALSE, TRUE),
                    ordered=TRUE)) %>%
    drop_na(historical_avg_increase) %>%
    arrange(desc(date), region_name) %>%
    select(-c(sum_past_rate_increase, n_past_rows))

```

</details>

This time, we don't create a `dataset_train` or `dataset_test`. But we will need to bake our dataset first, as our recipe wouldn't work for resampling:

```{r}
dataset_baked <- rec %>% bake(dataset)
```

With time-series data, the one we have, to create multiple training vs test splits, we have to use `rsample::rolling-window()` function. Following a similar design as [this official tidymodels tutorial](https://www.tidymodels.org/learn/models/time-series/):

```{r}

# A rolling window resampling
# Training: 50 years of data (50*12 months)
# Testing:  the 12-month period after training 
roll_rs <- rolling_origin(dataset_baked %>% 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)
```

How 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 152 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.

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

```{r}
# Let me get the first split
roll_rs$splits[[1]]
```
This tells me there 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 `analysis(roll_rs$splits[[1]])` or `assessment(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
assessment(roll_rs$splits[[1]]) %>% summarise(min_date=min(date), max_date=max(date))
```
Similarly, the second split will have a training data from '06/1969 - 05/2019' and a testing data from '06/2019 - 05/2020'. That is, we slide one month to the future.

## Running the same model multiple times (resampling)

How do we test our baseline model for all the different 76 data splits above? Instead of `fit()` we would use `fit_resamples()`. 


We specify the model as before:

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


Then we use the `fit_resamples()` function to run the model for all the different splits of our data:

```{r, message=FALSE, warning=FALSE}
# This takes a while to run
results <- fit_resamples(log_rec_specification, price_up ~ ., roll_rs)
```


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? That on average, the Area Under the Curve of this model is on average 0.66 (+- 0.019) on the test set. This is not too far, albeit a bit less than, the AUC we achieved on the single test set we had in the W05 mini competition and at the start of this notebook.

This is definitely a more robust way of assessing your model: running it multiple times on different samples of our data. Here, because time is an important factor of 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 3: 'Penalised' Logistic Regression

Now let's go back to `louis_model`, the complex model we had with 144 predictors.

Our question now is: What if we policed the model and got rid of the features (aka. predictors) that contribute less to the model?

The maths of how penalties work won't fit in this class, I might write a blog post about it later, but suffice to say that when you introduce a penalty, you **force** the model to reduce the values set to regression coefficients of the predictors ($\beta_1, \beta_2, \ldots$) to small numbers. Some will be forced entirely to zero. This is also often called _shrinkage_.

## Set up

Warning: you will need to install several other software packages to get this part to work, as we will not be using the default, `glm` engine that already comes pre-installed with tidymodels. There are many different ways (separate engines) to run logistic regresion models, you can read about each of them by going to the [`parsnip::logistic_reg()` documentation page](https://parsnip.tidymodels.org/reference/logistic_reg.html) and clicking on a desired engine.  

Initially I wanted to show you the example with the `LiblineaR` engine because it is a package that is easy to install. However, for some reason that is still unclear to me, RStudio breaks -- in any computer -- when I run a logistic regression on this data using this engine.

Therefore, for this tutorial, I will use the one in the `h2o` package. It allow us to tweak the parameters of the logistic regression, namely `penalty` and `mixture`, to control for the shrinkage we want to apply to this model. Penalty is a number that goes from 0 to infinity. The larger the penalty, the more we force the model to shrink and use fewer predictors. Mixture can only be set to 0 or to 1, and it changes how shrinkage penalty is applied (if you want to understand more, read the documentation `?logistic_reg` and research Lasso regression and Ridge Regression, or read Chapter 6 of the [ISLRv2 textbook](https://www.statlearning.com/)).

**(Optional) How to install `h2o`**

If you don't feel like going through the complication of installing extra software, you can just read the output of our models and try to understand what we are trying to achieve.

But if you want to go for the full experience:

1. Install [Java SDK 64-bit version](https://www.oracle.com/java/technologies/downloads/).
2. Install the agua package with `install.packages("agua")`
3. Initialise the H2o package:


```{r}
library(agua)

h2o::h2o.init()
```

4. This package doesn't like ordered factor (go figure!). So, if you used `factor(..., ordered=TRUE)` like I did in the baseline model, you would have to recreate it! When I recreated Louis' model, I removed the ordered parameter, so we'll be fine.

## Back to the data

We have to recreate the `dataset` object once again because it is different to the baseline model.

<details><summary>Click here to see the pre-processing steps for this model again</summary>

```{r}
dataset <- 
    uk_hpi %>% 
    filter(region_name %in% UK_countries) %>%
    group_by(region_name) %>%
    arrange(date) %>%
    mutate(
      lag_12_average_price = lag(average_price, 12),
      yearly_rate_increase = (average_price/lag_12_average_price) - 1, 
      lag_12_detached_price = lag(detached_price, 12), 
      yearly_rate_detached_increase = (detached_price/lag_12_detached_price) - 1, 
      lag_12_semi_detached_price = lag(semi_detached_price, 12), 
      yearly_rate_semi_detached_increase = (semi_detached_price/lag_12_semi_detached_price) - 1,
      lag_12_terraced_price = lag(terraced_price, 12),
      yearly_rate_terraced_increase = (terraced_price/lag_12_terraced_price) - 1,
      lag_12_flat_price = lag(flat_price, 12), 
      yearly_rate_flat_increase = (flat_price/lag_12_flat_price) - 1,
      lag_12_old_price = lag(old_price, 12), 
      yearly_rate_old_increase = (old_price/lag_12_old_price) - 1
    ) %>%
    drop_na(lag_12_average_price, 
            lag_12_detached_price, 
            lag_12_semi_detached_price,  
            lag_12_terraced_price, 
            lag_12_flat_price, 
            lag_12_old_price) %>%
    mutate(
      sum_past_rate_increase = cumsum(yearly_rate_increase) - yearly_rate_increase,  # Cumsum = cumulative sum
      sum_past_rate_detached_increase = cumsum(yearly_rate_detached_increase) - yearly_rate_detached_increase, 
      sum_past_rate_semi_detached_increase = cumsum(yearly_rate_semi_detached_increase) - yearly_rate_semi_detached_increase, 
      sum_past_rate_terraced_increase = cumsum(yearly_rate_terraced_increase) - yearly_rate_terraced_increase,
      sum_past_rate_flat_increase = cumsum(yearly_rate_flat_increase) - yearly_rate_flat_increase, 
      sum_past_rate_old_increase = cumsum(yearly_rate_old_increase) - yearly_rate_old_increase,
      n_past_rows = row_number() - 1,  # ROW_NUMBER function returns the row number over a named or unnamed window specification.
      historical_avg_increase = sum_past_rate_increase/n_past_rows,
      historical_avg_detached_increase = sum_past_rate_detached_increase/n_past_rows,
      historical_avg_semi_detached_increase = sum_past_rate_semi_detached_increase/n_past_rows,
      historical_avg_terraced_increase = sum_past_rate_terraced_increase/n_past_rows,
      historical_avg_flat_increase = sum_past_rate_flat_increase/n_past_rows,
      historical_avg_old_increase = sum_past_rate_old_increase/n_past_rows,
      price_up=factor(yearly_rate_increase > historical_avg_increase, 
                      labels=c("No","Yes"),
                      levels=c(FALSE, TRUE),
                      ordered=FALSE)) %>%
    drop_na(historical_avg_increase, historical_avg_detached_increase, historical_avg_semi_detached_increase, historical_avg_terraced_increase, historical_avg_flat_increase, historical_avg_old_increase) %>%
    arrange(desc(date), region_name) %>%
    select(-c(sum_past_rate_increase, n_past_rows))

```

</details>

Let's bake it again:

```{r}
dataset_baked <- rec_louis %>% bake(dataset)

# Let's get just the list of predictors or we will have issues with metadata columns later
list_predictors <- rec_louis$term_info %>% filter(role == "predictor") %>% pull(variable)

# I will get rid of metadata (keep just the date) to avoid issues with the recipes package
dataset_baked <- 
  dataset_baked %>% 
  group_by(region_name) %>% 
  arrange(date) %>% 
  ungroup() %>%
  select(date, all_of(list_predictors), price_up)

```

Because we used many many lags in this version of the data, we won't have 50-years worth of past data to use to train our model, our training sets in the resampling will have to do with less (say 20 years):


```{r}
# A rolling window resampling
# Training: 20 years of data (20*12 months)
# Testing:  the 12-month period after training 
roll_rs <- rolling_origin(dataset_baked,   
               initial = 4 * 12 * 20, # x4 countries x 12x50 months 
               assess = 4 * 12, # x4 countries x 12 months
               cumulative = FALSE)

roll_rs %>% nrow()
```
There are still a good number of resamples!


## Specifying the model

The way we specify this penalised model is almost the same as before, but now I will tweak the parameter `penalty` to show you what happens:

```{r}
log_reg_h2o_spec <- 
    logistic_reg(penalty = 0.1) %>% 
    set_engine("h2o") %>% 
    set_mode("classification")
```

Say I want to run this only on a single training set:

```{r}
# This gets the training data in the first split
training_data <- analysis(roll_rs$splits[[1]]) %>% select(-date)

penalised_model <- 
  log_reg_h2o_spec %>% 
  fit(price_up ~ ., data = training_data)
```

Because I set penalty to a nonzero value, I applied a penalty to my logistic regression and, therefore, the regression coefficients were shrinked. Some of them might have been set to zero. How do we know?

Let's look at the summary of the model:

```{r}
penalised_model
```

The second tab has information about how many predictors were actually active in the final model: 37, as opposed to 144 predictors that _could_ have been used. There is a lot more info on the other tabs: coefficients, confusion matrix, metrics, thresholds, etc.

**How to get the coefficients?**

This info is hidden inside the model attribute of the fitted model:

```{r}
penalised_model$fit@model$coefficients
```

**How do I get the nonzero coefficients?**

```{r}
all_coefs <- penalised_model$fit@model$coefficients
nonzero_coeffs <- all_coefs != 0
all_coefs[nonzero_coeffs]
```

Don't get me wrong, this is still a complex model, but less than the original one. Note, for example, how the regression coefficients are smaller in absolute number -- they range from around -2 to around + 2 -- than the ones in the original more complex model, which had extremely unreasonably large regression coefficients such as 5.476e+16.


## Further tweak penalty

Say, you want to **shrink** the model further. Just increase penalty! Instead of `penalty=0.1`, let me try `penalty=0.7`.

```{r}
log_reg_h2o_spec <- 
    logistic_reg(penalty = 0.7) %>% 
    set_engine("h2o") %>% 
    set_mode("classification")

penalised_model <- 
  log_reg_h2o_spec %>% 
  fit(price_up ~ ., data = training_data)

# What are the coefficients?
all_coefs <- penalised_model$fit@model$coefficients
nonzero_coeffs <- all_coefs != 0
all_coefs[nonzero_coeffs]
```

The **only** coefficient is $\hat{\beta_0} = -0.5085759$! What does that mean? This means that we forced the shrinkage so much that we forced ALL regression coefficients to zero. Essentially, our model has been reduced to a very simple equation:

$$
\text{price\_up} = \frac{e^{-0.484246}}{1 + e^{-0.484246}} = \approx 38.125\%
$$
In other words: all predictions will be the same! If your threshold is above the percentage above, all predictions will be `price_up = "Yes"`, and it will `price_up = "No"` otherwise -- in the case where your threshold is less than the only percentage you can get out of this model.

**This is not a useful model!**

How do I choose the right penalty that will produce a 'good' model, using all possible 144 predictors? It all depends on what metrics matter the most to you, of course. Let's assume that Area Under the Curve (AUC) is what we care about the most. We can follow what Bosco did and write a for loop with many different values of penalty, assess the metric on the training and test set, and pick the one that lead to the best result.

Our code would look somewhat like the following:

```r
# Vary from 0.10, 0.15, 0.20, 0.25, ..., 0.50
all_penalties <- seq(0.10, 0.50, 0.05)

best_auc <- NULL

for(penalty in all_penalties){

  ...fit the model...
  
  ...calculate the current auc...
  
  if(is.null(best_auc) or current_auc > best_auc){
    best_auc <- current_auc
  }

}
```

The above is valid and it would work. But we could make use of the [`tune` package](https://tune.tidymodels.org/) of tidymodels to 'tune' this parameter for us using fewer lines of code. More on that in the future!


