---
title: "🗓️ Week 04 -Fundamentals of Classification ([DS202](https://lse-dsi.github.io/DS202))"
subtitle: "Theme: Supervised Learning"
author: "Dr. Jon Cardoso-Silva"
format: 
  html:
    page-layout: full
    toc: true
    toc-depth: 2
self-contained: true
---

Welcome to Week 04!

Instead of slides, we are going to use a Quarto markdown as the lecture material for this week. This is to encourage you to follow along on your own computers and help you practice taking notes of concepts, equations and code using markdown.

⚙️ **Set up**

```{r import, message=FALSE, warning=FALSE}
library(dplyr)
library(lubridate)
library(ggplot2)

library(parsnip)
library(yardstick)
library(recipes)
library(workflows)
```

```{r read_data, message=FALSE, warning=FALSE, echo=FALSE}
uk_hpi <- 
  readr::read_csv("../../../../data/UK-HPI-full-file-2023-06.csv") %>%
    mutate(Date = dmy(Date)) %>%
    filter(Date >= ymd(20050101)) %>%
    ungroup()
```

**💡 Tip: the janitor package**

Before we begin with the proper content, let me give you a tip: the [`janitor` package](https://sfirke.github.io/janitor/) has very useful functions!

Here are our old columns:

```{r}
uk_hpi %>% colnames()
```

I can use `janitor::clean_names()` to make them more standardised and avoid having to use the backticks:

```{r}
uk_hpi <- uk_hpi %>% janitor::clean_names()
uk_hpi %>% colnames()
```

(I will use these columns from now on)

# 🥅 Learning Objectives

<details>

<summary>Click here to read if you are seeing this document as a web page</summary>

What is the point of today?

The hope is that, by the end of today, you will learn the following things:

**1st Hour**

-   Understand the term 'goodness-of-fit'
-   Use (Adjusted) R-squared, Pearson correlation coefficient to compare linear models
-   Obtain a deeper insight into training and testing splits
-   Notice that transforming variables can be a useful tool when modelling
-   Discover the usefulness of dummy variables
-   Apply different metrics to assess the of regression algorithms

**2nd Hour**

-   Distinguish regression versus classification
-   Run a logistic regression (a classification method!) model
-   The very basics of classification metrics (confusion matrix)

</details>

# 1: Assumptions + Bad vs Good Models

Sometimes, we can achieve the same thing if we change the way we frame the problem. Along this section, I want to show you that I could predict `sales_volume` in the UK without using the `sales_volume` column at all!

## 1.1 Pre-processing: build the dataset

For this example, I will create the following predictors:

-   `lag_1_sales_volume`: simply the lagged sales volume. You already know this one well.

-   `rate_change_1m`: the rate of change in the sales volume from the current to the previous month.

-   `rate_change_2m`: similar to above but represents the rate of change relative to two months prior.

Let's build a dataset with UK countries that contains both:

```{r message=FALSE}
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_sales_volume = lag(sales_volume, 1),
      lag_2_sales_volume = lag(sales_volume, 2), 
      rate_change_1m = (sales_volume / lag_1_sales_volume - 1),
      rate_change_2m = (sales_volume / lag_2_sales_volume - 1)
    ) %>%
    tidyr::drop_na(sales_volume, lag_1_sales_volume, lag_2_sales_volume, 
                   rate_change_1m, rate_change_2m) %>%
    arrange(desc(date), region_name) %>%
    select(date, region_name, sales_volume, 
           lag_1_sales_volume, lag_2_sales_volume,
           rate_change_1m, rate_change_2m)

# This is what we have
dataset %>% head(8) 
```

**Another training vs testing split**...

To assess the capabilities of this model to generalise to unseen data, I will use a temporal train-test split. I will train the model on data from 2005 to 2018 and test it on data from 2019 onwards.

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

👉 In the past labs and exercises, I invited you to set end of 2019 as a cut-off point for the test set. This is because I wanted you to see if anyone would spot that 2020 was the year of the pandemic! The year were all Machine Learning algorithms trained on historical data failed to predict the future! Here I will stop one year earlier in case we have time to make comparisons between the test set performance in 2019 and 2020.

💡 This split is the same as pretending we are in December 2018 and we have all the data until now. We just don't know what the future holds.

## 1.2 Don't forget to visualise your assumptions

Even before I start investigating a model, let me actually visually inspect if it makes sense:

```{r message=FALSE}
# Remember: ggplot uses + instead of %>%
ggplot(dataset_train, 
       aes(x = lag_1_sales_volume, y = sales_volume, color=region_name)) +
  geom_point(size=2, alpha=0.3, stroke=1) +
  # geom_smooth fits a linear model just for the sake of the plot
  geom_smooth(method = "lm", se = FALSE, color="black") + 
  scale_color_brewer(name="Country", palette = "Set1") +
  scale_x_continuous(labels=scales::unit_format(unit = "K", scale = 1e-3)) +
  scale_y_continuous(labels=scales::unit_format(unit = "K", scale = 1e-3)) +
  labs(x = "Last month's sales volume", y = "This month's sales volume",
       title="It is fair to say there is a linear relationship",
       subtitle="However, the variance is much higher for England data.") + 
  theme_bw()
```

🤔 I mean, the linear assumption seems to hold but the England data will always have more 'leverage' in the model because it is such a clear outlier.

What about the rate of change?

```{r message=FALSE}
ggplot(dataset_train, 
       aes(x = rate_change_1m, y = sales_volume, color=region_name)) +
  geom_point(size=2, alpha=0.3, stroke=1) +
  # geom_smooth fits a linear model just for the sake of the plot
  geom_smooth(method = "lm", se = FALSE, color="black") + 
  scale_x_continuous(labels=scales::percent_format()) +
  scale_y_continuous(labels=scales::unit_format(unit = "K", scale = 1e-3)) +
  scale_color_brewer(name="Country", palette = "Set1") +
  labs(x = "Rate of change in sales volume", y = "This month's sales volume",
       title = "This here, ladies and gentlemen, theys and thems, is a bad linear model",
       subtitle = "We don't even need to look at the metrics, we can SEE that it is bad") + 
  theme_bw()
```

🤔 This combination does not make a lot of sense. Our linear regression would be way off.

💡 But what if we put the two variables in similar ranges? I could model the percentage change in sales volume instead of the absolute change. Let's try that!

Let's create a plot of `rate_change_1m` versus `rate_change_2m`:

```{r message=FALSE, warning=FALSE}
ggplot(dataset_train, 
         aes(x = rate_change_2m, y = rate_change_1m, color=region_name)) +
  geom_point(size=2, alpha=0.3, stroke=1) +
  # geom_smooth fits a linear model just for the sake of the plot
  geom_smooth(method = "lm", se = FALSE, color="black") + 
  scale_x_continuous(labels=scales::percent_format()) +
  scale_y_continuous(labels=scales::percent_format()) +
  scale_color_brewer(name="Country", palette = "Set1") +
  labs(x = "Rate of change in sales volume (2 months ago)", 
       y = "Rate of change in sales volume (1 month ago)",
       title = "Our linear assumptions seem to make more sense here!") + 
  theme_bw()
```

👉 Visually, this looks like a good assumption to make. It is definitely a better model than the previous one, **but how good is it compared to `sales_volume ~ lag_1_sales_volume`**?

How do we know if one linear model is better than another?

-   By comparing the **goodness-of-fit** of the models using standardised metrics
-   By comparing predictions of the models on a test set, **using the same unit of measurement**

## 1.3 Comparing models

Let's create one first model `model1` using `rate_change_1m` as a direct predictor of `sales_volume`:

```{r}
model1 <- 
  linear_reg() %>%
  set_engine("lm") %>%
  fit(sales_volume ~ lag_1_sales_volume, data = dataset_train)
```

Then, to align with our last plot, let's create a `model2` that uses `rate_change_2m` as a predictor of `rate_change_1m` (not `sales_volume`):

```{r}
model2 <- 
  linear_reg() %>%
  set_engine("lm") %>%
  fit(rate_change_1m ~ rate_change_2m, data = dataset_train)
```

If we use the MAE, we would be comparing apples with oranges because the units of measurement are different. But there's one metric that appears when you run `summary()` on a linear model that is agnostic to the unit of measurement: the $R^2$ metric.

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

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

## 1.4 $R^2$ as a measures of goodness-of-fit

When you were asked to calculate, for example, the MAE of a model when applied to the training set, you were calculating the goodness-of-fit.

> 📕 Goodness-of-fit is the concept of how well a model fits the data it was trained on.

The MAE has the same unit of measurement as the variable you are trying to predict, but there are some measures that are agnostic to the unit of measurement and it might be useful to know them.

The R-squared metric uses the concept of residuals, which you have explored in the previous weeks. The $R^2$ metric is defined as:

$$
R^2 = 1 - \frac{RSS}{TSS}
$$

where $RSS = \sum_{i=1}^n (y_i - \hat{y}_i)^2$ (the sum of squares of the residuals) and $TSS = \sum_{i=1}^n (y_i - \bar{y})^2$ is the total sum of squares of all the data points used to fit the model [^1]. The closer $R^2$ is to 1, the better the model *explains the variance* in the data. The closer it is to 0, the worse the model is.

[^1]: See [An Introduction to Statistical Learning](https://www.statlearning.com/) by Hastie et al. Page 70 of the 2nd edition (freely available online).

The $R^2$ metric is closely related to the Pearson correlation coefficient ($\rho$):

$$
\rho = \sqrt{R^2}
$$

where $\rho$ is the Pearson correlation coefficient.

The interpretation is a little different: $\rho$ ranges from -1 to 1. Values of $|\rho|$ closer to 1 indicate that the two variables are more strongly **linearly correlated**. Values of $|\rho|$ closer to 0 indicate that the two variables are less strongly **linearly correlated**.

```{r}
# Function cor returns the Pearson correlation coefficient
cor(dataset_train$sales_volume, dataset_train$lag_1_sales_volume)
```

vs

```{r}
cor(dataset_train$rate_change_1m, dataset_train$rate_change_2m)
```

⚠️ **Note:** these two metrics are only useful in the context of linear models. It is definitely possible for two variables to be strongly associated with each other but in a *non-linear* way. For example, the relationship between the number of hours you spend studying and your exam score is likely to be non-linear. The more you study, the better your exam score will be, but only up to a point. After that point, you will start to get diminishing returns.

What can we learn from this? Because I am still focused in linear associations, these metrics are useful to me and therefore I can be more confident that `model2` is a better model than `model1`.

But *'How do I know **for sure**'*, I hear you ask? Well, this is more of a deep epistemological question. I will adopt the convenient philosophy that all models are wrong, but some are useful.

But we can do more! We can compare the predictions of both models on a test set. Let's do that.

# 2. Don't stick just to goodness-of-fit though, go deeper

It is very easy to conclude by looking at measures of goodness-of-fit that `model1` is superior to `model2`. It is true that the first model is a better linear fit. However, what I want you to remember when comparing different models is the following:

> *All models are wrong, but just how useful these models are?*

## 2.1 Why use a test set?

The whole point of having a test set is to check if our model is as good as it sounds on paper. We want our models to be useful in the real world, not just when used on the training data.

Thinking this way, let's compare the predictions made by `model1` and `model2` and although I could still continue to compare them with unit-agnostic measures like $R^2$, I want to bring you back to the Mean Absolute Error - it will be easier to interpret and make it more 'real-world'.

Predictions made by `model1` are already made in unit of sales:

```{r}
model1 %>% predict(dataset_test) %>% head()
```

Whereas the predictions made by `model2` are in the unit of rate of change:

```{r}
model2 %>% predict(dataset_test) %>% head()
```

## 2.2 Residuals plot (on test set)

### 2.2.1 Easy to do for `model1`

```{r message=FALSE}
plot_df <- model1 %>% augment(dataset_test)

g <- (
  ggplot(plot_df, aes(x = sales_volume, y=.resid, color=region_name)) +
  geom_hline(yintercept=0, linetype="dashed", color="#212121") +
  geom_point(size=2, alpha=0.3, stroke=1) +
  scale_color_brewer(name="Country", palette = "Set1") +
  scale_x_continuous(labels=scales::unit_format(unit = "K", scale = 1e-3)) +
  scale_y_continuous(labels=scales::unit_format(unit = "K", scale = 1e-3),
                     limits=c(-1e5, 1e5)) +
  labs(x = "Sales volume (data)", 
       y = "Residual sales volume",
       title="Residuals plot for model1") + 
  theme_bw()
) 

g
```

### 2.2.2 We gotta do some work to get the same plot for `model2`

Let's convert the predictions of `model2` to units of sales. Here is how `rate_change_1m` was created:

```         
rate_change_1m = (sales_volume / lag(sales_volume, 1) - 1)
```

After I predict (or augment) from `model2`, I will do a little calculation to create a column I will call `pred_sales_volume`:

```{r}
model2 %>%
  augment(dataset_test) %>%
  mutate(
    pred_sales_volume = (1 + .pred) * lag_1_sales_volume
  ) %>%
  select(lag_1_sales_volume, lag_2_sales_volume,
         rate_change_1m, pred_sales_volume, sales_volume) %>%
  head()
```

Then, let's create recreate the residuals plot but using `sales_volume` as reference.

```{r message=FALSE}
plot_df <- 
  model2 %>%
  augment(dataset_test) %>%
  mutate(
    pred_sales_volume = (1 + .pred) * lag_1_sales_volume,
    resid_sales_volume = sales_volume - pred_sales_volume
  )

g <- (
  ggplot(plot_df, aes(x = sales_volume, y=resid_sales_volume, color=region_name)) +
  geom_hline(yintercept=0, linetype="dashed", color="#212121") +
  geom_point(size=2, alpha=0.3, stroke=1) +
  scale_color_brewer(name="Country", palette = "Set1") +
  scale_x_continuous(labels=scales::unit_format(unit = "K", scale = 1e-3)) +
  scale_y_continuous(labels=scales::unit_format(unit = "K", scale = 1e-3),
                     limits=c(-1e5, 1e5)) +
  labs(x = "Sales volume (data)", 
       y = "Residual sales volume",
       title="Residuals plot for model2") + 
  theme_bw()
) 

g
  
```

## 2.3 Other metrics (test set)

### 2.3.1 Correlation coefficient (test set)

Compare the correlation coefficient between `sales_volume` and `pred_sales_volume` in the two models:

**Model 1**

```{r}
model1 %>%
  augment(dataset_test) %>%
  summarise(pearson_correlation=cor(.pred, sales_volume))
```

**Model 2**

```{r}
model2 %>%
  augment(dataset_test) %>%
  mutate(
    pred_sales_volume = (1 + .pred) * lag_1_sales_volume
  ) %>%
  summarise(pearson_correlation=cor(sales_volume, pred_sales_volume))
```

### 2.3.2 MAE

You know this one already!

**Model 1**

``` r
model1 %>%
  augment(dataset_test) %>%
  group_by(region_name, year=year(date)) %>%
  mae(sales_volume, .pred) %>%
  rename(mae=.estimate) %>%
  select(region_name, year, mae) %>%
  tidyr::pivot_wider(names_from=region_name, values_from=mae) %>%
  knitr::kable()
```

| year |   England | Northern Ireland |  Scotland |     Wales |
|-----:|----------:|-----------------:|----------:|----------:|
| 2019 |  6427.576 |         392.8642 | 1148.5991 |  387.3800 |
| 2020 |  9645.571 |         574.6722 | 1466.5728 |  529.7189 |
| 2021 | 37582.779 |         429.3154 | 1879.7946 | 1016.1121 |
| 2022 |  4989.231 |         395.9884 |  775.4967 |  480.3977 |
| 2023 |  7305.470 |         555.8381 | 1682.9940 |  769.4262 |

: **\[Model 1\]** Mean Absolute Error per country per year (test set)

**Model 2**

``` r
model2 %>%
  augment(dataset_test) %>%
  mutate(
    pred_sales_volume = (1 + .pred) * lag_1_sales_volume
  ) %>%
  group_by(region_name, year=year(date)) %>%
  mae(sales_volume, pred_sales_volume) %>%
  rename(mae=.estimate) %>%
  select(region_name, year, mae) %>%
  tidyr::pivot_wider(names_from=region_name, values_from=mae) %>%
  knitr::kable()
```

| year |   England | Northern Ireland |  Scotland |    Wales |
|-----:|----------:|-----------------:|----------:|---------:|
| 2019 |  4310.027 |         88.78054 |  991.2044 | 268.4131 |
| 2020 |  6238.819 |        383.18635 |  783.3222 | 372.3814 |
| 2021 | 26500.811 |         96.10372 | 1366.5793 | 798.5002 |
| 2022 |  3731.289 |         47.68254 |  515.5830 | 270.3233 |
| 2023 |  4237.989 |        158.60761 | 1229.1318 | 319.4935 |

: **\[Model 2\]** Mean Absolute Error per country per year (test set)

Model 2 is pretty much always better...

### 2.3.2 MAPE

As I have reinforced elsewhere, it will help if you think of MAE in comparison to some other baseline. For example: Is MAE ok, in comparison to the usual variance of the sales volume?

Another way to look at this is by comparing the Mean Absolute Percentage Error (MAPE). There is a handy [`yardstick::mape()` function](https://yardstick.tidymodels.org/reference/mape.html) to calculate this for us, but the equation is as follows:

$$
\text{MAPE} = \frac{1}{n}\sum_{i=1}^n{\left|\frac{y_i - \hat{y_i}}{y_i}\right|}
$$

That is: on average, how much does the residual represent from the actual baseline value?

| region_name      | year |  model1 |  model2 |
|:-----------------|-----:|--------:|--------:|
| England          | 2019 |  9.82 % |  6.61 % |
| England          | 2020 | 24.03 % | 15.34 % |
| England          | 2021 | 54.28 % | 42.75 % |
| England          | 2022 |  8.42 % |  6.37 % |
| England          | 2023 | 20.44 % | 11.89 % |
| Northern Ireland | 2019 | 19.66 % |  4.60 % |
| Northern Ireland | 2020 | 48.11 % | 25.08 % |
| Northern Ireland | 2021 | 17.21 % |  3.92 % |
| Northern Ireland | 2022 | 19.05 % |  2.34 % |
| Northern Ireland | 2023 | 36.62 % | 10.50 % |
| Scotland         | 2019 | 14.27 % | 11.69 % |
| Scotland         | 2020 | 26.93 % | 14.78 % |
| Scotland         | 2021 | 20.91 % | 15.62 % |
| Scotland         | 2022 |  9.62 % |  6.06 % |
| Scotland         | 2023 | 25.45 % | 17.86 % |
| Wales            | 2019 | 10.90 % |  7.22 % |
| Wales            | 2020 | 24.27 % | 15.82 % |
| Wales            | 2021 | 25.34 % | 20.35 % |
| Wales            | 2022 | 14.37 % |  7.81 % |
| Wales            | 2023 | 37.37 % | 15.46 % |

: Comparison of `model1` and `model2` in terms of Mean Absolute Percentage Error

<details>

<summary>Click here to see the code that produced this summary table</summary>

``` r
model1_test_results <-
    model1 %>%
    augment(dataset_test) %>%
    group_by(region_name, year=year(date)) %>%
    mape(sales_volume, .pred) %>%
    rename(mape=.estimate) %>%
    select(region_name, year, mape) %>%
    mutate(model="model1")

model2_test_results <-
  model2 %>%
  augment(dataset_test) %>%
  mutate(
      pred_sales_volume = (1 + .pred) * lag_1_sales_volume
  ) %>%
  group_by(region_name, year=year(date)) %>%
  mape(sales_volume, pred_sales_volume) %>%
  rename(mape=.estimate) %>%
  select(region_name, year, mape) %>%
  mutate(model="model2")

bind_rows(model1_test_results, model2_test_results) %>% 
    mutate(mape=(function(x){sprintf("%.2f %%", x)})(mape)) %>%
    tidyr::pivot_wider(names_from=model, values_from=mape) %>%
    knitr::kable()
```

</details>

## 2.4 Go deeper: Compare yearly predictions (2019-now)

Let's create one big plot where we have a `sales_volume` line, the `pred_sales_volume` line from `model1` and the `pred_sales_volume` line from `model2`. Let's use the months as the x-axis and the sales volume as the y-axis. Let's also facet by `region_name`:

```{r}
# I will start by copying my dataset_test to plot_df
plot_df <- dataset_test

# Add a column with the predictions from model1
plot_df$model1 <- model1 %>% augment(dataset_test) %>% pull(.pred)

# Add a column with the predictions from model2
plot_df$model2 <- 
  model2 %>%
  augment(dataset_test) %>%
  mutate(pred_sales_volume = (1 + .pred) * lag_1_sales_volume) %>%
  pull(pred_sales_volume)

# Let's see how the data looks like
plot_df %>% head()
```

**I need to use some advanced tricks now**

```{r}
# Now I will use and advanced technique called 'pivot_longer' to make the data 'tidy'
plot_df <- 
  plot_df %>%
  tidyr::pivot_longer(
    cols = c(sales_volume, model1, model2),
    names_to = "model",
    values_to = "sales_volume"
  )

plot_df %>% head()
```

Now I can plot the data:

```{r fig.width=8, fig.height=7}
# You will likely have to tweak the width and height (or even size of geom_point)
# to adjust to your monitor

g <- (
  ggplot(plot_df,
         aes(x = date, y = sales_volume, color=model, linetype=model)) +
  geom_line(linewidth=0.75, alpha=1) +
  geom_point(alpha=0.6, size=2) + 
  scale_linetype_manual(name="Data or model", 
                        values=c("model1"="dashed", 
                                 "model2"="dashed", 
                                 "sales_volume"="solid")) +  
  scale_colour_manual(name="Data or model",
                      values=c("model1"="#FE6F5E", 
                               "model2"="#00939D", 
                               "sales_volume"="#595959")) + 
  scale_y_continuous(labels=scales::unit_format(unit = "K", scale = 1e-3),
                     expand=c(0.25, 0.25)) + 
  scale_x_date(breaks="6 months", 
               date_labels = "%b %Y",
               date_minor_breaks="1 month") + 
  facet_wrap(~region_name, scales = "free_y", ncol=1) +
  labs(x = "Date", y = "Sales volume and Predictions",
      title="Model 2 can handle even the fluctuations of the market",
      subtitle="Plotting is a very useful tool to compare models!") + 
  theme_bw()
)

#ggsave(g, filename="./figures/comparison_plot.png", width=8, height=4)

g
```

🎯 **ACTION POINTS**

Later at home, tweak the code above to see different versions of this plot:

-   Zoom in on a country

-   Zoom in on a particular year

💡 I am focusing on simple linear regression, with just one predictor, here because it is easier to communicate these points. It should be easy for you to translate the comparisons above for multiple linear regression.

⚠️ **WARNING:** R-squared [will always improve](https://lse-dsi.github.io/lse-ds202-course-notes/blog/posts/r-squared.html) if you add more predictors, even if the predictors are random. So, if you are comparing models with different predictors, use the adjusted R-squared instead.

# 3. Other tricks

## 3.1 Dummy variables

What if let my model tweak the intercept and slope according to the country?

Conceptually, I want to create a model that behaves sort of like:

$$
\text{rate\_change\_1m} = \beta_0 + \beta_1 \times \text{rate\_change\_2m} + \beta_2^*\times\text{region\_name}
$$ 

(Tip: you won't see these `\` characters in the rendered (final) version of this notebook.)

This is what recipes are for. Let's hope they lead us astray today! This is how dummy variables actually look like in a dataset:

```{r}

rec_dummy <-
  recipe(rate_change_1m ~ rate_change_2m + region_name, data=dataset_train) %>%
  step_dummy(region_name) %>%
  prep()

# I will run this in console (not suitable for qmd):
# rec_dummy %>% bake(dataset_train) %>% View()
```

Let's create a `model3` for that:

```{r}
# I will do all in one pipe chain
model3_workflow <-
  workflow() %>%
  add_recipe(rec_dummy) %>%
  add_model(linear_reg() %>% set_engine("lm")) %>%
  fit(dataset_train)

model3 <- model3_workflow %>% extract_fit_parsnip()

model3$fit
```

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

This is the MAPE comparison:

| region_name      | year |  model1 |  model2 |  model3 |
|:-----------------|-----:|--------:|--------:|--------:|
| England          | 2019 |  9.82 % |  6.61 % |  6.62 % |
| England          | 2020 | 24.03 % | 15.34 % | 15.34 % |
| England          | 2021 | 54.28 % | 42.75 % | 42.78 % |
| England          | 2022 |  8.42 % |  6.37 % |  6.39 % |
| England          | 2023 | 20.44 % | 11.89 % | 11.90 % |
| Northern Ireland | 2019 | 19.66 % |  4.60 % |  4.52 % |
| Northern Ireland | 2020 | 48.11 % | 25.08 % | 24.97 % |
| Northern Ireland | 2021 | 17.21 % |  3.92 % |  3.87 % |
| Northern Ireland | 2022 | 19.05 % |  2.34 % |  2.28 % |
| Northern Ireland | 2023 | 36.62 % | 10.50 % | 10.49 % |
| Scotland         | 2019 | 14.27 % | 11.69 % | 11.70 % |
| Scotland         | 2020 | 26.93 % | 14.78 % | 14.85 % |
| Scotland         | 2021 | 20.91 % | 15.62 % | 15.64 % |
| Scotland         | 2022 |  9.62 % |  6.06 % |  6.06 % |
| Scotland         | 2023 | 25.45 % | 17.86 % | 17.89 % |
| Wales            | 2019 | 10.90 % |  7.22 % |  7.23 % |
| Wales            | 2020 | 24.27 % | 15.82 % | 15.83 % |
| Wales            | 2021 | 25.34 % | 20.35 % | 20.35 % |
| Wales            | 2022 | 14.37 % |  7.81 % |  7.83 % |
| Wales            | 2023 | 37.37 % | 15.46 % | 15.47 % |

: TABLE 3. things got slightly better for the Northern Ireland prediction, but other than that, `model3` barely made any difference. At least now you know a way to create 'dummy' variables if the opportunity arrives.

## 3.2 Feature Selection

# 4. Classification: a different type of supervised learning model

## 4.1 Making binary predictions

Now say I have a brand new ambition: to predict whether average prices will go up (`Yes` or `No`) a year from now. It is a very simplistic model because it doesn't take into account any fluctuations or inflation adjustments, but let's use it as a toy example for now.

Allow me to start over with my pre-processing and to create a target variable `price_up`:

```{r message=FALSE, error=FALSE}
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),
      price_up = if_else(average_price > lag_12_average_price, "Yes", "No")
    ) %>%
    tidyr::drop_na(price_up) %>%
    arrange(desc(date), region_name) %>%
    select(date, region_name, average_price, lag_12_average_price, price_up)

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

# Run View(dataset) to browse the new data
```

```{r}
 dataset_train %>% group_by(price_up) %>% tally()
```

It's pretty much always up!

OK, but how do I write a model for the `price_up` variable?

## 4.2 Could I do linear regression?

In theory, yes, but only if you convert the `Yes` and `No` **labels** to numbers. Say: 1 and 0, respectively.

```{r}
model_yes_no <-
  linear_reg() %>%
  set_engine("lm") %>%
  fit(price_up ~ lag_12_average_price, 
      data=dataset_train %>% mutate(price_up=if_else(price_up == "Yes", 1, 0)))
model_yes_no
```

What does this mean? Well, since linear regression is a **regression** algorithm, i.e., it operates on numbers, I can do as before and plot the predictions versus the actual values:

```{r}
plot_df <- model_yes_no %>% augment(dataset_test %>% mutate(price_up=if_else(price_up == "Yes", 1, 0)))

min_x <- plot_df %>% filter(lag_12_average_price == min(lag_12_average_price)) %>% slice(1)
max_x <- plot_df %>% filter(lag_12_average_price == max(lag_12_average_price)) %>% slice(1)

g <- (
  ggplot(plot_df, aes(x = lag_12_average_price, y=.pred, color=region_name)) +
  geom_point(size=2, alpha=0.3, stroke=1) +
  annotate(geom="segment", 
           x=min_x$lag_12_average_price, 
           y=min_x$.pred,
           xend=max_x$lag_12_average_price,
           yend=max_x$.pred) +
  scale_color_brewer(name="Country", palette = "Set1") +
  scale_x_continuous(labels=scales::label_dollar(prefix = "£", suffix="K", scale=1e-3)) +
  scale_y_continuous(limits=c(0, NA)) +
  labs(x = "Average price 12 months ago", 
       y = "Prediction for variable 'price up' in the current month",
       title="Linear regression is not a good model for this problem",
       subtitle="The predictions are not bounded between 0 and 1",
       caption="I am not even sure what to make of this!") + 
  theme_bw()
)
g
```

NO. Linear Regression is not suitable for predicting labels. It would be more useful if we were to predict **probabilities** instead of single numbers.

## 4.3. Logistic Regression

Consider a **binary** response:

$$
Y = \begin{cases}
0 \\
1
\end{cases}
$$

We model the probability that $Y = 1$ using the **logistic function** (aka. sigmoid curve):

$$
Pr(Y = 1|X) = p(X) = \frac{e^{\beta_0 + \beta_1X}}{1 + e^{\beta_0 + \beta_1 X}}
$$

This function is mathematically guaranteed to be bounded by the interval $[0, 1]$.

### 4.3.1 Insight into the logistic function

Let's see how this mathematical function looks like:

```{r}
plot_logistic_function <- function(X_range, beta0, beta1, color){
  
  df <- expand.grid(X=X_range, beta0=beta0, beta1=beta1)
  
  df <- 
    df %>% group_by(beta0, beta1) %>% 
    transform(X=X, Y=exp(beta0+beta1*X)/(1+exp(beta0+beta1*X)))
  
  df <- df %>% 
    transform(beta0=factor(beta0, ordered = TRUE),
              beta1=factor(beta1, ordered = TRUE))
  
  g <- 
    if(missing(color)){
      ggplot(df, aes(x=X, y=Y))
    }else if(color=="beta0"){
      ggplot(df, aes(x=X, y=Y, color=beta0, linetype=beta0))
    }else{
      ggplot(df, aes(x=X, y=Y, color=beta1, linetype=beta1))
    }
  
  g <- 
    g + 
    geom_line(linewidth=1.4) +
  
    xlab("X") + ylab("Y") +
  
    # THEMING
    theme_bw() +
  	theme(
          axis.title.x = element_text(size=24,margin = margin(t = 10, r = 0, b = 0, l = 0)),
          axis.title.y = element_text(size=24, margin = margin(t = 0,  r = 10, b = 0, l = 0)),
          legend.text = element_text(size=18),
          legend.position="right",
          axis.text.x = element_text(size=18),
          axis.text.y = element_text(size=18),
          panel.grid.minor = element_blank(),
          aspect.ratio = 1/2)
  return(g)
}
```

```{r}
plot_logistic_function(X_range=(-400:400)/100, beta0=c(-2, 0, 1), beta1=c(1), color="beta0") + 
labs(title=expression(paste("If I keep ", beta[1]==1," fixed and vary ", beta[0] )))
```

```{r}
plot_logistic_function(X_range=seq(-4, 4, 0.1), beta0=c(0), beta1=c(-1, 0, 1, 2), color="beta1") + 
labs(title=expression(paste("If I keep ", beta[0]==0," fixed and vary ", beta[1] )))
```

### 4.3.2 Building a model

Looking up the [list of models available on `tidymodels`](https://www.tidymodels.org/find/parsnip/), I find that there are many implementations of logistic regression. Here I am going to use the one from the `glm` package.

That is, something like:

``` r
logistic_model1 <-
  logistic_reg() %>%
  set_engine("glm") %>%
  fit(price_up ~ lag_12_average_price, data=dataset_train)
```

However, if you try to run the code above, you will get an error saying:

> Error in `check_outcome()`: ! For a classification model, the outcome should be a `factor`, not a `character`.

So we need to convert our `price_up` from a character ("Yes" vs "No") variable to a factor, a categorical variable.

Let's go back to creating our dataset:

```{r message=FALSE, error=FALSE}
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),
      price_up = if_else(average_price > lag_12_average_price, "Yes", "No"),
      price_up = factor(price_up, levels=c("No", "Yes"))
    ) %>%
    tidyr::drop_na(price_up) %>%
    arrange(desc(date), region_name) %>%
    select(date, region_name, average_price, lag_12_average_price, price_up)

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

# Run View(dataset) to browse the new data
```

```{r}
logistic_model1 <-
  logistic_reg() %>%
  set_engine("glm") %>%
  fit(price_up ~ lag_12_average_price, data=dataset_train)

# Look at the coefficients
logistic_model1$fit
```

It worked! And the output is quite similar to that of a linear regression.

Let's see what happens when I augment the dataset with predictions:

```{r}
logistic_model1 %>% augment(dataset_train) %>% arrange() %>% head()
```

I have several columns that relate to predictions!

### 4.3.3 Plot the model

Let's plot the output of the model so we can better understand what it represents:

```{r}
plot_df <- logistic_model1 %>% augment(dataset_train)

min_x <- plot_df %>% filter(lag_12_average_price == min(lag_12_average_price)) %>% slice(1)
max_x <- plot_df %>% filter(lag_12_average_price == max(lag_12_average_price)) %>% slice(1)

g <- (
  ggplot(plot_df, aes(x = lag_12_average_price, y=.pred_Yes, color=region_name)) +
  geom_point(size=2, alpha=0.3, stroke=1) +
  scale_color_brewer(name="Country", palette = "Set1") +
  scale_x_continuous(labels=scales::label_dollar(prefix = "£", suffix="K", scale=1e-3)) +
  labs(x = "Average price 12 months ago", 
       y = "Probability of price going up",
       title="Logistic regression is a good model for this problem",
       subtitle="The predictions are bounded between 0 and 1",
       caption="I am not even sure what to make of this!") + 
  theme_bw()
)
g
```

### 4.3.4 How do we interpret the coefficients?

We know that:

$$
\beta_0 \approx 0.9491292 \\
\beta_1 \approx -8.877\times10^{-7}
$$

This means that:

1.  In the case where the average price 12 months ago was zero, the probability of the price going up in the current month is:

$$
\frac{e^{0.9491292}}{1 + e^{0.9491292}} = 0.721
$$

2.  For every £1 increase in the average price 12 months ago, the probability of the price going up increases by a factor of $-8.877\times10^{-7}$

🤔 "*increased by a factor of... a number that is much lower than 1".* Yes, I know, this can be very confusing.

Other ways to think about this:

- If the average price is £ 200k today, I can plug this into the equation above and conclude that there is $68.38583 \%$ chance that this average price will go up a year from now.

- To say that the probability will 'increase by a factor of $\beta_1$' means that: if instead of £200k, the average price today was £200,001, I don't need to plug the numbers into the equation again. I can reuse what I calculated for the case where price = £ 200k.

  - How? Multiply the probability of price = £ 200k by $\beta_1$ =
  
$$
0.6838583 + 0.6838583 \times -8.877\times10^{-7} \approx 68.38575 \%
$$

## 4.4 Making predictions

At what point do you say: Given these probabilities, I can confidently say the price will go up!

By default, the **threshold** is set at probability=0.5.

```{r}
logistic_model1 %>% 
  augment(dataset_train) %>%
  summarise(min_prob=min(.pred_Yes),
            q1=quantile(.pred_Yes, 0.25),
            median=quantile(.pred_Yes, 0.5),
            q3=quantile(.pred_Yes, 0.75),
            max=max(.pred_Yes))
```

So... what prediction will I make if my threshold is 0.5?

```{r}
logistic_model1 %>% 
  augment(dataset_train) %>%
  group_by(.pred_class, price_up) %>%
  tally()
```
What If I change the threshold?

```{r}
logistic_model1 %>% 
  augment(dataset_train) %>%
  # I made up a rule entirely:
  mutate(.pred_class=if_else(.pred_Yes < quantile(.pred_Yes, 0.25), "Yes", "No"),
         .pred_class=factor(.pred_class, levels=c("No", "Yes"))) %>%
  group_by(.pred_class, price_up) %>%
  tally()
```

If I rearrange, I get what is called a **confusion matrix**:

```{r}
logistic_model1 %>% 
  augment(dataset_train) %>%
  # I made up a rule entirely:
  mutate(.pred_class=if_else(.pred_Yes < quantile(.pred_Yes, 0.25), "Yes", "No"),
         .pred_class=factor(.pred_class, levels=c("No", "Yes"))) %>%
  group_by(.pred_class, price_up) %>%
  tally() %>%
  tidyr::pivot_wider(values_from=n, names_from="price_up")
```

## 4.4 Adding dummy variables

```{r}
rec_dummy <-
  recipe(price_up ~ lag_12_average_price + region_name, data=dataset_train) %>%
  step_dummy(region_name) %>%
  prep()
```

Let's create a `logistic_model2` for that:

```{r}
# I will do all in one pipe chain
logistic_model2_workflow <-
  workflow() %>%
  add_recipe(rec_dummy) %>%
  add_model(logistic_reg() %>% set_engine("glm")) %>%
  fit(dataset_train)

logistic_model2 <- logistic_model2_workflow %>% extract_fit_parsnip()

logistic_model2$fit
```
**Is this a better or a worse model?**

```{r}
logistic_model2_workflow %>% 
  augment(dataset_train) %>%
  # I made up a rule entirely:
  mutate(.pred_class=if_else(.pred_Yes < quantile(.pred_Yes, 0.25), "Yes", "No"),
         .pred_class=factor(.pred_class, levels=c("No", "Yes"))) %>%
  group_by(.pred_class, price_up) %>%
  tally() %>%
  tidyr::pivot_wider(values_from=n, names_from="price_up")
```

# What's Next

In the lab next week, you will be tasked with creating a logistic regression. To make things more exciting, we want to see who will be able to create the most accurate model, as shown by the confusion matrix!

