---
title: "🥇 Week 05 - Competition annotated"
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
---

# ⚙️ Setup

Important: make sure your `tidymodel` is up-to-date. If not, open the R Console and run the following command:

```r
install.packages("tidymodels")
```

You might also need to install:

```r    
install.packages("themis")
install.packages("probably")
```

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

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

# Submission #1: Spark

## Data Preparation

**The dataset and the splits**

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

    ## THESE ARE OUR FEATURES
    mutate(
        lag_4_yearly_rate_increase = lag(yearly_rate_increase, 4), 
        lag_4_historical_avg_increase = lag(historical_avg_increase, 4)
    ) %>%
    drop_na(historical_avg_increase, lag_4_yearly_rate_increase, lag_4_historical_avg_increase) %>%
    arrange(desc(date), region_name) %>%
    dplyr::select(-c(sum_past_rate_increase, n_past_rows))

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

## Specifying and fitting the model

**Model specification**

Spark went for a simple choice of features for his model: a 4-month lag of the yearly rate increase and the historical average increase:

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

rec_baseline <- 
    recipe(price_up ~ lag_4_yearly_rate_increase + lag_4_historical_avg_increase, data = dataset_train) %>% 
    step_naomit(all_predictors(), skip=FALSE) 

baseline_wf <-
    workflow() %>%
    add_recipe(rec_baseline) %>%
    add_model(log_rec_specification)
```

```{r}
summary(rec_baseline) %>% filter(role %in% c("predictor", "outcome"))
```

**Model fitting**

```{r}
baseline_fit <- baseline_wf %>% fit(data = dataset_train)
baseline_model <- baseline_fit %>% extract_fit_parsnip()
```

**Model evaluation**

```{r}
train_augmented <-    
    baseline_model %>% 
    augment(dataset_train)

g <- 
    train_augmented %>%
    roc_curve(truth = price_up, .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
```

## Decision-making & Metrics (Training Set)

**Is there an optimal threshold?**

With a good model, there is no need to fiddle with the threshold. The ROC curve is already very close to the top-left corner.

```{r}
# See the optimal threshold based on the J-index
train_augmented %>%
    threshold_perf(truth = price_up, estimate = .pred_Yes, event_level="second") %>%
    pivot_wider(names_from=.metric, values_from=.estimate) %>% 
    arrange(desc(j_index)) %>%
    head(2)
```

0.5 is already the optimal threshold.

**Metrics**

```{r}
cm <- 
    train_augmented %>%
    conf_mat(truth=price_up, estimate=.pred_class)

g <- cm %>% 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
```

```{r}
multi_metric <- metric_set(accuracy, sens, spec, precision, recall, f_meas)

train_augmented %>%
    multi_metric(truth = price_up, estimate=.pred_class)
```

## The real test: the test set

Does it perform as well in the test set?

```{r}
test_augmented <-    
    baseline_model %>% 
    augment(dataset_test)

test_augmented %>%
    multi_metric(truth = price_up, estimate=.pred_class)
```

Not bad!

**Confusion matrix on the test set**

```{r}
cm <- 
    test_augmented %>%
    conf_mat(truth=price_up, estimate=.pred_class)

g <- cm %>% 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
```

# Submission #2: Jacinta

## Data Preparation

**The dataset and the splits**

Jacinta used several lagged variables. I am replicating it as she did, with the data pre-processing steps such as lagging done in `dplyr`, not in the `recipes`.

```{r}
UK_countries <- c("England", "Wales", "Scotland", "Northern Ireland")

dataset <- 
    uk_hpi %>% 
    filter(region_name %in% UK_countries) %>%
    group_by(region_name) %>%
    arrange(date) %>%
    
    # define percentage change in price
    mutate(
        lag_12_average_price = lag(average_price, 12),
        lag_3_average_price = lag(average_price, 3),
        lag_1_average_price = lag(average_price, 1),
        monthly_rate_increase =  (average_price/lag_1_average_price) - 1,
        quarterly_rate_increase = (average_price/lag_3_average_price) - 1,
        yearly_rate_increase = (average_price/lag_12_average_price) - 1
    ) %>%
  
    drop_na(lag_12_average_price) %>%
    drop_na(lag_3_average_price) %>%
    drop_na(lag_1_average_price) %>%
  
    # lag one period for all regressors 
    mutate(
        lag_1_monthly_rate_increase = lag(monthly_rate_increase, 1),
        lag_1_quarterly_rate_increase = lag(quarterly_rate_increase, 1),
        lag_1_new12m_percent_change = lag(new12m_percent_change, 1),
     ) %>%

  
    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)
    ) %>%
   
    mutate(
        lag_1_yearly_rate_increase = lag(yearly_rate_increase, 1),
        lag_1_historical_avg_increase = lag(historical_avg_increase, 1))%>%
    
    drop_na(lag_1_monthly_rate_increase) %>%
    drop_na(lag_1_quarterly_rate_increase) %>%
    drop_na(lag_1_new12m_percent_change) %>%
    drop_na(lag_1_yearly_rate_increase) %>%
    drop_na(lag_1_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)
```

Jacinta framed her thought process around demand and supply (**domain knowledge**), and as far as I can tell, there is no data leakage.

She wrote:

> Relevant factors:
>
> demand and supply: price factors: (lag dependent variables) - monthly_rate_increase - quarterly_rate_increase - yearly_rate_increase - price of new and old house: new12m_percent_change, old12m_percent_change
>
> demand - finance type: mortgage12m_percent_change, mortgage1m_percent_change, cash1m_percent_change, cash12m_percent_change - two types of buyer: first-time buyer and former owner occupier index: ftb1m_percent_change, ftb12m_percent_change, foo12m_percent_change, foo1m_percent_change supply: monthly_sales_volume_increase

## Specifying and fitting the model

**Model specification**

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

rec_baseline <- 
    recipe(price_up ~ ., data = dataset_train) %>% 
    update_role(-c(price_up), new_role = "ID") %>%
    update_role(lag_1_monthly_rate_increase, 
                lag_1_quarterly_rate_increase, 
                lag_1_yearly_rate_increase, 
                lag_1_historical_avg_increase,
                lag_1_new12m_percent_change,
                new_role = "predictor") %>%
    step_naomit(all_predictors(), skip=FALSE)

baseline_wf <-
    workflow() %>%
    add_recipe(rec_baseline) %>%
    add_model(log_rec_specification)
```

In the end, her model had the following features:

-   lag_1_monthly_rate_increase
-   lag_1_quarterly_rate_increase
-   lag_1_new12m_percent_change
-   lag_1_yearly_rate_increase
-   lag_1_historical_avg_increase

```{r}
summary(rec_baseline) %>% filter(role %in% c("predictor", "outcome"))
```

However, the `dataset` object had several other lagged variables (including some lag_12_*). I suspect Jacinta played around with different combinations of the resulting features before settling on this one.

```{r}
dataset %>% select(starts_with("lag")) %>% colnames()
```

**Model fitting**

```{r}
baseline_fit <- baseline_wf %>% fit(data = dataset_train)
baseline_model <- baseline_fit %>% extract_fit_parsnip()
```

Here we get a warning message:

> Warning message:
> glm.fit: fitted probabilities numerically 0 or 1 occurred 

Typically, this means that the model has regression coefficients that are too large. This can be due to multicollinearity or overfitting. It is important to check the model's performance.


**Assessing the goodness-of-fit**

```{r}
train_augmented <-    
    baseline_model %>% 
    augment(dataset_train)

g <- 
    train_augmented %>%
    roc_curve(truth = price_up, .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 super ROC Curve on the Training Set",
       subtitle="This curve shows how the model performs under different thresholds",
       x="False Positive Rate\n(1 - Specificity)", y="True Positive Rate\n(Sensitivity)")

g
```

## Decision-making & Metrics (Training Set)

**Is there an optimal threshold?**

With a good model, there is no need to fiddle with the threshold. The ROC curve is already very close to the top-left corner.

```{r}
# See the optimal threshold based on the J-index
train_augmented %>%
    threshold_perf(truth = price_up, estimate = .pred_Yes, event_level="second") %>%
    pivot_wider(names_from=.metric, values_from=.estimate) %>% 
    arrange(desc(j_index)) %>%
    head(2)
```

(Technically, J-index suggests 0.55 is the optimal threshold, but the ROC curve is already very close to the top-left corner, so 0.5 would also be a good choice.)

**Update `train_augmented`**

```{r}
threshold <- 0.55

train_augmented <-    
    baseline_model %>% 
    augment(dataset_train) %>% 
    mutate(
        .pred_class = .pred_Yes > threshold,
        .pred_class = factor(.pred_class, 
                             labels=c("No", "Yes"), 
                             levels=c(FALSE, TRUE), 
                             ordered=TRUE))
```

**Metrics**

```{r}
cm <- 
    train_augmented %>%
    conf_mat(truth=price_up, estimate=.pred_class)

g <- cm %>% 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
```

```{r}
multi_metric <- metric_set(accuracy, sens, spec, precision, recall, f_meas)

train_augmented %>%
    multi_metric(truth = price_up, estimate=.pred_class)
```

## The real test: the test set

Does it perform as well in the test set?

```{r}
test_augmented <-    
    baseline_model %>% 
    augment(dataset_test) %>%
    mutate(
        .pred_class = .pred_Yes > threshold,
        .pred_class = factor(.pred_class, 
                             labels=c("No", "Yes"), 
                             levels=c(FALSE, TRUE), 
                             ordered=TRUE))

test_augmented %>%
    multi_metric(truth = price_up, estimate=.pred_class)
```

**Confusion matrix on the test set**

```{r}
cm <- 
    test_augmented %>%
    conf_mat(truth=price_up, estimate=.pred_class)

g <- cm %>% 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
```

# Submission #3: Terry

## Data Preparation

**The dataset and the splits**

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

    ## THESE ARE OUR FEATURES
    mutate(
        lag_1_average_price = lag(average_price, 1), 
        lag_1_index = lag(index, 1),
        lag_1_sales_volume = lag(sales_volume, 1),
        lag_1_x1m_percent_change = lag(x1m_percent_change, 1), 
        lag_1_x12m_percent_change = lag(x12m_percent_change, 1),
        lag_1_yearly_rate_increase = lag(yearly_rate_increase, 1),
        lag_1_historical_avg_increase = lag(historical_avg_increase, 1)
    ) %>%

    drop_na(historical_avg_increase, starts_with("lag")) %>%

    arrange(desc(date), region_name) %>%
    select(-c(sum_past_rate_increase, n_past_rows, lag_12_average_price))

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

## Specifying and fitting the model

**Model specification**

Terry chose to use the `recipe` package to specify the model. He also included some advanced and very clever data preparation steps, such as down-sampling the majority class and creating lagged variables.

```{r}
library(themis) # Required for step_downsample

rec_train <- 
  recipe(price_up ~ ., data = dataset_train) %>% 
  update_role(-c(price_up), new_role = "ID") %>%
  update_role(c(starts_with("lag"), contains("region_name")), new_role = "predictor") %>%
  step_dummy(region_name, skip=FALSE) %>%
  step_naomit(all_predictors(), skip=FALSE) %>%
  step_downsample(price_up, under_ratio = 1) # Downsample the majority class

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

train_wf <- workflow() %>% 
  add_recipe(rec_train) %>% 
  add_model(log_rec_spec)
```

```{r}
summary(rec_train %>% prep()) %>% filter(role %in% c("predictor", "outcome"))
```

We end up with a model that includes the following features:

- region_name_Northern.Ireland 
- region_name_Scotland 
- region_name_Wales 
- lag_1_average_price 
- lag_1_index 
- lag_1_sales_volume
- lag_1_x1m_percent_change
- lag_1_x12m_percent_change
- lag_1_yearly_rate_increase
- lag_1_historical_avg_increase

**Model fitting**

```{r}
baseline_fit <- train_wf %>% fit(data = dataset_train)
baseline_model <- baseline_fit %>% extract_fit_parsnip()
```

**Model evaluation**

```{r}
train_augmented <-    
    baseline_model %>% 
    augment(rec_train %>% prep() %>% bake(dataset_train))

g <- 
    train_augmented %>%
    roc_curve(truth = price_up, .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="ROC Curve on the Training Set",
       subtitle="This curve shows how the model performs under different thresholds",
       x="False Positive Rate\n(1 - Specificity)", y="True Positive Rate\n(Sensitivity)")

g
```

## Decision-making & Metrics (Training Set)

**Is there an optimal threshold?**

Yes, there is. But, with a good model, there is really no need to fiddle with the threshold. The ROC curve is already very close to the top-left corner. The gain is marginal.

```{r}
# See the optimal threshold based on the J-index
train_augmented %>%
    threshold_perf(truth = price_up, estimate = .pred_Yes, event_level="second") %>%
    pivot_wider(names_from=.metric, values_from=.estimate) %>% 
    arrange(desc(j_index)) %>%
    head(3)
```



**Metrics**

```{r}
cm <- 
    train_augmented %>%
    conf_mat(truth=price_up, estimate=.pred_class)

g <- cm %>% 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
```

```{r}
multi_metric <- metric_set(accuracy, sens, spec, precision, recall, f_meas)

train_augmented %>%
    multi_metric(truth = price_up, estimate=.pred_class)
```

## The real test: the test set

Does it perform as well in the test set?

```{r}
test_augmented <-    
    baseline_model %>% 
    augment(rec_train %>% prep() %>% bake(dataset_test))

test_augmented %>%
    multi_metric(truth = price_up, estimate=.pred_class)
```

Brilliantly!

**Confusion matrix on the test set**

```{r}
cm <- 
    test_augmented %>%
    conf_mat(truth=price_up, estimate=.pred_class)

g <- cm %>% 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
```