---
title: "💻 LSE 202A 2023: Week 03 - Lab"
author: "MY NAME (MY CANDIDATE NUMBER)"
date: 29 January 2024
format: html
self-contained: true
---

# ⚙️ Setup

**Import required libraries:**

```{r setup, message=FALSE, warning=FALSE}
# Remember to install.packages() any missing packages
library(dplyr)
library(readr)
library(lubridate)
library(tidymodels)
```

**Read the data set:**

It is the same data set you used in the previous lab (as well as in the first formative assignment):

```{r}
# Modify the filepath if needed
filepath <- "data/UK-HPI-full-file-2023-06.csv"
uk_hpi <- read_csv(filepath)
```

**Pre-process the dataset:**

```{r}
# Replicate roughly the same steps you did in the previous lab:

# 1. Select the columns you need
# 2. Rename the columns
# 3. Filter the data set to include only the regions you need (let's focus on England for now)
# 4. Preprocess the date column

selected_cols    <- c("Date", "RegionName", "1m%Change", "SalesVolume")

df_england <-
    uk_hpi %>%
    select(all_of(selected_cols)) %>%
    rename(date = Date, region = RegionName, monthlyChangeHPI = `1m%Change`) %>%
    filter(region == "England") %>%
    drop_na(monthlyChangeHPI) %>% 
    mutate(date = lubridate::dmy(date)) %>% 
    arrange(desc(date))

df_england  %>% head(8) %>% knitr::kable()
```

# 📋 Lab Tasks

## Part I - the `lm` command (25 min)

👨🏻‍🏫 **TEACHING MOMENT:** Your class teacher will demonstrate how to run the `lm` function.

**Task 1:** Run a linear regression model with `SalesVolume` as the dependent variable and `monthlyChangeHPI` as the independent variable. Save the model as `model1`.

```{r}
model1 <- lm(SalesVolume ~ monthlyChangeHPI, data = df_england)
```

Print the model summary:

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

🗣️ **DISCUSSION:**

- What is the value of the intercept?

> $\approx 69427.6$

- What is the value of the slope?

> $\approx 12008.9$

- How should you interpret this model?

>[This model suggests that whenever the **house prices index** increases by 1 **percentage point** on a period of **one month** (this corresponds to our `monthlyChangeHPI` independent variable), then the **sales volume** (this corresponds to our `SalesVolume` dependent variable) increases by **a factor of 12008 approximately**.]{style="font-size:0.65em"}
> If the model is a good fit, we can also expect the slope to be correct within a 95% confidence interval of 12008.9 +/- 943.6 (given by the Std. Error of the slope coefficient).


**👨🏻‍🏫 TEACHING MOMENT:** Your class teacher will demonstrate how to plot the model.

**Task 2:** Plot the residuals against the fitted values.

```{r}
plot(model1, 1)
```

Alternatively, see all four basic lm diagnostic plots:

```{r}
par(mfrow = c(2, 2))
plot(model1)
```

🗣️ **DISCUSSION:**

Is this a good fit? Why? Why not?

> The residual plot suggests the relationship is not exactly linear. This model is not a good fit!

## Part II - the `tidymodels` way (20 min)

👨🏻‍🏫 **TEACHING MOMENT:** Your class teacher will demonstrate how to get the same thing done using the `tidymodels` framework.

**Task 3:** Run a linear regression model with `SalesVolume` as the dependent variable and `monthlyChangeHPI` as the independent variable. Save the model as `model2`.

```{r}
model2 <-
    linear_reg() %>%
    set_engine("lm") %>%
    fit(SalesVolume ~ monthlyChangeHPI, data = df_england)
```

Print the model coefficients in a tidy way:

```{r}
model2 %>% tidy()
```

See the full model summary:

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

Overall statistics:

```{r}
model2 %>% glance()
```

**Task 4:** Augument the model with the residuals and the fitted values and use the resulting data frame to plot the residuals against the fitted values.

```{r}
model2 %>% augment(df_england)
```

Then use `ggplot2` to plot the residuals against the fitted values:

```{r}
ggplot(model2 %>% augment(df_england), aes(.pred, .resid)) + 
    geom_point()
```

You can further customise the plot:

```{r}
plot_df <- model2 %>% augment(df_england)
g <- ggplot(plot_df, aes(x = .pred, y = .resid)) +
    geom_point(alpha=0.2, size=3, color="red", stroke=1) +
    geom_hline(yintercept = 0, linetype = "dashed") +
    labs(x = "Fitted values", y = "Residuals", title="Residuals vs Fitted") + 
    theme_bw() + 
    theme(axis.title.x = element_text(size = rel(1.2)), 
          axis.text.x = element_text(size = rel(1.2)),
          axis.title.y = element_text(size = rel(1.2)),
          axis.text.y = element_text(size = rel(1.2)),
          plot.title = element_text(size = rel(1.5), face = "bold"))
g
```

## Part III - evaluating the model

**STEP 01**

> 1. Separate the data into training and testing sets. The training set must contain data up until Dec 2020. The testing set must contain data from Jan 2021 onwards.

The reasoning here is as follows:

- I need to split my data into two separate data frames.
- The questions specifies what it calls 'training set' as all the data up until Dec 2020. 
- Since we are subsetting the data by selecting rows, we can use the `filter` function to select the rows we need.
- To make date comparisons, it is imperative that all date objects and columns are stored as proper date objects. 
- I have already properly modified the data type of the `date` column using [lubridate](https://lubridate.tidyverse.org/). Therefore, I don't need to modify anything else in that column.
- However, I need to convert the string of text `"2021-01-01"` to a proper date object. I can use lubridate's `ymd` function to convert a string to a date object.
- Great, now I can use the `filter` function to select the rows I need.

```{r task_03_step01}
df_england_train <- df_england %>% filter(date <   ymd("2021-01-01"))
df_england_test  <- df_england %>% filter(date >=  ymd("2021-01-01"))
```

Why do you need to convert date-related columns and objects to a date object using `lubridate`? The original data set has date in the format "DD/MM/YYYY" as a string. If you try to compare strings, you will get the wrong result because the comparison is done lexicographically (i.e. by alphabetical order).

For example, the following code will return `TRUE`, implying that 1 January 2021 happens before 1 December 2020:

```{r}
"01/01/2021" < "01/12/2020"
```

**STEP 02**

> 2. Fit a linear regression model with `tidymodels` using just the training set. Call it `model3`. How does this model compare to the one you fitted in Part II?

How would you reason about this?

- I need to fit a linear regression model using `tidymodels`. That is, copying the code from Part II (not with `lm()`).
- I need to modify just the source of the data. Instead of using the entire data set, I need to use just the training set I created in STEP 01.

```{r task_03_step02}
model3 <-
    linear_reg() %>%
    set_engine("lm") %>%
    fit(SalesVolume ~ monthlyChangeHPI, data = df_england_train)
```

In terms of comparison, if I tidy up the output of this new model, I will notice that it is very similar to the one we fitted in Part II. The coefficients are almost identical:

```{r}
model3 %>% tidy()
```

**STEP 03**

> 3. Calculate the Mean Absolute Error (MAE) for the **training set**. That is, on average, how far off is the model from the actual data?

How would you reason about this?

- I don't know what the MAE is. I need to look it up.
- I go to the [tidymodels](https://www.tidymodels.org/) website and search for "MAE" and "Mean Absolute Error"
- The results didn't exactly showed me how to calculate the MAE. I wonder if there is something in one of the packages that are part of tidymodels.
- Eventually, you should find the `yardstick` package. You can read the documentation of the `mae()` function [here](https://yardstick.tidymodels.org/reference/mae.html).
- From the documentation of the `mae()` function, I look at the examples and deduce that I need to compare the model predictions with the actual values used to fit it.
- I can use the `augment()` function to get the predictions of the model on the training set. Then I can use the `mae()` function to calculate the MAE.

```{r task_03_step03}
model3 %>% 
    augment(df_england_train) %>%
    mae(truth = SalesVolume, estimate = .pred)
```

I am just printing out the value, but - alternatively - I can save it in a variable if I need to use it later:

```{r task_03_step03_save}
mae_train <- 
    model3 %>% 
    augment(df_england_train) %>%
    mae(truth = SalesVolume, estimate = .pred)
mae_train
```

The MAE is approximately 14213. That means on average, the model is off by 14213.

**STEP 04**

> 4. Calculate the Mean Absolute Error (MAE) for the **testing set**. That is, on average, how far off is the model from the actual data?

Here is when things start to get interesting. How would you reason about this?

- I have already trained a model, `model3` on two variables: `SalesVolume` and `monthlyChangeHPI`.
- This means I can use the model to predict the `SalesVolume` for any given value of `monthlyChangeHPI`, **on any data set** -- it does not have to be the same data I used to train the model.
- Therefore, I can use the `augment()` function to get the predictions of the model on the testing set. Then I can use the `mae()` function to calculate the MAE, just the same way I did before.

```{r task_03_step04}
model3 %>% 
    augment(df_england_test) %>%
    mae(truth = SalesVolume, estimate = .pred)
```

OH WOW! This is a much higher error than the one we got for the training set. The model is wrong on average almost twice as much as it was on the training set. This can happen for many reasons: it could be that the model is so attuned to the training set that it is not able to **generalise** to the testing set, or it could mean that something changed in the dynamics that produced this data, and the model -- fitted on older data -- is not able to capture it. Which one do you think it is?

💡 A common source of confusion in this question was that many people assumed that we would have to train a new model. That is not the case, and this is the beauty of modelling. Once you have a model, you can use it to predict the outcome on any data set, as long as the data set has the same variables as the ones you used to train the model. (We will explore more of the implications and limitations of thinking this way in the future, in this course)

**STEP 05**

5. Now, add a second independent variable, call it `monthlyChangeHPI_lag1`, which is essentially a **lagged variable** of the `monthlyChangeHPI` variable you already have. Build a new model, call it `model4`, and compare it to `model3`. Which one is better?

**Add lagged variables:**

```{r}
# Add lagged variables for the dependent variable
df_england <- 
    df_england %>%
    arrange(date) %>% # Sort in ascending order so we can use lag later
    mutate(monthlyChangeHPI_lag1 = lag(monthlyChangeHPI, 1)) %>%
    drop_na() %>%
    arrange(desc(date)) # sort by date in descending order just to make it easier to read
```

**Split dataset in training and test sets:**

```{r}
df_england_train <- df_england %>% filter(date <   ymd("2021-01-01"))
df_england_test  <- df_england %>% filter(date >=  ymd("2021-01-01"))
```

**Building new regression model `model4`:**

```{r}
model4 <-
  linear_reg() %>%
  set_engine("lm") %>%
  fit(SalesVolume ~ monthlyChangeHPI + monthlyChangeHPI_lag1, data = df_england_train)
```

If we tidy the output of this new model, here is what we get:

```{r}
model4 %>% tidy()
```

**Evaluating the new model**

```{r}
model4 %>% 
  augment(df_england_train) %>%
  mae(truth = SalesVolume, estimate = .pred)
```

For the training set, the MAE for `model4` is slightly lower than for `model3`.
```{r}
model4 %>% 
  augment(df_england_test) %>%
  mae(truth = SalesVolume, estimate = .pred)
```

When it comes to test set too, the performance of `model4` is better than that of `model3`.

It looks as if `model4` fits the data better than `model3`!

**STEP 06**

> 6. Now create a `df_scotland` with the same variables as `df`, but only for Scotland (no need to separate into training and testing sets).

I can use the same code I used to create `df_england` in the beginning of this lab, but I need to change the filter condition to select only Scotland:

```{r task_03_step05}
df_scotland <-
    uk_hpi %>%
    select(all_of(selected_cols)) %>%
    rename(date = Date, region = RegionName, monthlyChangeHPI = `1m%Change`) %>%
    filter(region == "Scotland") %>%
    drop_na(monthlyChangeHPI) %>% 
    mutate(date = lubridate::dmy(date)) %>% 
    arrange(date) %>%
    mutate(monthlyChangeHPI_lag1 = lag(monthlyChangeHPI, 1)) %>%
    drop_na() %>%
    arrange(desc(date))
```

Checking that it is correct:

```{r task_03_step05_check}
df_scotland %>% head(8)
```

Further sanity check:

```{r task_03_step05_check2}
df_scotland$region %>% unique() # Retrieve all unique values in the `region` column
```

**STEP 07**

> 8. How well can `model3` and `model4` predict monthly changes in house prices in Scotland?

Now that you have realised that `model3` can be used in any dataset, you can use it to predict the monthly changes in house prices in Scotland the same way you did for the England test set. We can also calculate the MAE so we can compare it with the MAE we calculated for the England test set.


```{r task_03_step06}
model3 %>% 
    augment(df_scotland) %>%
    mae(truth = SalesVolume, estimate = .pred)
```

The error is also higher, on average, than the one we got for the England training set. This is not surprising, since the model was trained on England data, and Scotland is a different country with different dynamics.

(Extra) Out of curiosity, we could see if the MAE changes over time. I could group the data by year and calculate the MAE for each year:

```{r task_03_step06_extra}
model3 %>%
    augment(df_scotland) %>%
    group_by(year = lubridate::year(date)) %>%
    mae(yearly_change, .pred) # MAE is already a summary statistics, so I don't need to use summarise()
```

While I am at it, I can plot it:

```{r}
plot_df <- 
    model3 %>%
    augment(df_scotland) %>%
    group_by(year = lubridate::year(date)) %>%
    mae(yearly_change, .pred)

g <- 
    (
        ggplot(plot_df, aes(x=year, y=.estimate, color=.metric))
            + geom_line()
            + geom_point()
            + ylim(c(0, NA)) # Ensure the Y axis starts at zero. NA means "whatever the maximum value is"
            + labs(x = "Year", y = "MAE", title="MAE over time")
    )
g
```

nterestingly, the worst predictions were for 2004, not during the pandemic. This is indicative that the dynamics in Scotland changed around that time, and the model is not able to capture it as it was trained on data from a different country!

Let's do the same for `model4`. Is `model4` a better fit for Scotland than `model3` was?

```{r}
model4 %>% 
  augment(df_scotland) %>%
  mae(truth = SalesVolume, estimate = .pred)
```

The error is still much higher, on average, than for this model the one we got for the England training set. The same conclusions that applied for `model3` apply here as well!

Let's, as before, plot MAE over time:

```{r}
plot_df <- 
  model4 %>%
  augment(df_scotland) %>%
  group_by(year = lubridate::year(date)) %>%
  mae(SalesVolume, .pred)

g <- 
  (
    ggplot(plot_df, aes(x=year, y=.estimate, color=.metric))
    + geom_line()
    + geom_point()
    + ylim(c(0, NA)) # Ensure the Y axis starts at zero. NA means "whatever the maximum value is"
    + labs(x = "Year", y = "MAE", title="MAE over time")
  )
g
```

MAE is highest, as before, in 2004: this is the second model to suggest that the dynamics of the housing market in Scotland experienced some changes around that time, and the model is not able to capture it as it was trained on data from a different country! Some deeper investigations are needed here!

