---
title: "💻 LSE DS202A 2024: Week 03 - Homework"
author: "MY NAME (MY CANDIDATE NUMBER)"
date: 16 October 2024
format: html
self-contained: true
---

# Tasks
  
## Part I - Getting Ready 

### ⚙️ Setup

**Import required libraries:**

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

# Tidyverse packages we will use
library(dplyr)     
library(ggplot2)
library(tidyr)     
library(readr)
library(stringr)
library(rsample) #this is what allows you to split data in training and test set!

# Tidymodel packages we will use
library(yardstick) # this is where metrics like mae() are
library(parsnip)   # This is where models like linear_reg() are
library(recipes)   # You will learn about this today
library(workflows) # You will learn about this today

# And finally a package for cleaning variable names
library(janitor)
library(utils)

```

**Read the data set:**

It is the same data set you used in the previous lab:

It is the same data set you used in this week's lab:

```{r}
# Modify the filepath if needed
filepath <- "data/WVS_Wave_7.csv"
wvs <- read_csv(filepath)
```

As in this week's lab, if you find that your laptop is unable to handle the full data set without running slowly, try experimenting with the following code (reminder: this takes a random sample of the data, stratifying by country so the sampling algorithm doesn't take more data from one country and less from another).

```{r}
# Set a seed for reproducibility

set.seed(123)
  
wvs <-
  # Load the .csv
  read_csv("data/WVS_Wave_7.csv") %>%
  # Sample a proportion of the data set for each country.
  # We have used 25% but you can experiment depending on
  # the capability of your machine.
  group_by(iso3c) %>% 
  slice_sample(prop = 0.25) %>% 
  ungroup()

```


## Part II - Tidymodels recipes

Say our goal is to predict `satisfaction` using the variables available in the dataset with the exception of `-iso3c` (just as you did in this week's lab!). That is:

``` {r}
# Split the data with 75% being used to train the model

wvs_split <- initial_split(wvs, prop = 0.75)

# Create tibbles of the training and test set

wvs_train <- training(wvs_split)
wvs_test <- testing(wvs_split)
    

# Then train it with tidymodels:
model <- 
    linear_reg() %>%
    set_engine("lm") %>%
    fit(satisfaction ~ .-iso3c, data = wvs_train)
```

But today, we want to show you a different way of pre-processing your data before fitting a model!

**We start with a recipe**

The problem with the code above is that, if you need to use the fitted `model` to make predictions on new data, you will need to apply the same pre-processing steps to that new data. This is not a problem if you are only doing it once, but if you need to do it many times, it can become a bit of a pain. This is where the [recipes](https://recipes.tidymodels.org/) package (of tidymodels) comes in handy.

A recipe is a set of instructions of how the data should be *pre-processed* before it is used in a model. It is a bit like a cooking recipe, where you have a list of ingredients and a set of instructions on how to prepare them. In our case, the ingredient is the data frame, and the instructions are the pre-processing steps we want to apply to them.

Here's how we would construct a recipe to pre-process our data before fitting a linear regression model:

```{r}
# Create a recipe
rec <- 
  recipe(satisfaction ~ . , data = wvs) %>%
  step_rm(iso3c)


summary(rec)

```

In this particular case, the only pre-processing we do is remove the variable `iso3c` from consideration in the model: that's what's `step_rm` does. Other examples of pre-processing could include removing the observations/rows that contain missing values ((`step_naomit()`)[https://recipes.tidymodels.org//reference/step_naomit.html]) or imputing (i.e 'taking a best guess' at) them (e.g (`step_impute_median`)[https://recipes.tidymodels.org//reference/step_impute_median.html] or (`step_impute_linear`)[https://recipes.tidymodels.org//reference/step_impute_linear.html]).

:::callout-warning
The previous syntax is enough if all you're going to do is use your recipe to train a model. However, if you want to, say check the correctness of your pre-processing by printing out the actual dataset that you obtain after pre-processing, then you'll need a few extra steps to get there.

The previous code will look like this:

```{r}
rec <- 
  recipe(satisfaction ~ . , data = wvs) %>%
  step_rm(iso3c) %>%
  prep() %>%
  bake(wvs_train)
```

`prep()` estimates all the quantities and statistics needed for the recipe you've specified (i.e preprocessing steps) to be applied by the `bake` function (note that you specify which dataset the `bake` function, and ultimately the pre-processing steps, apply on - in this case the training set). You can think of it as the (ingredients) prepping step in a cooking recipe before you bake/cook everything and get the final meal! 

⚠️ If you don't go through the prep and bake steps, the pre-processing steps will not be applied. The recipe is just a recipe, not the final cooked meal. You need to bake it first.

:::

**How to use this in a model?**

So we have our recipe from before:

```{r}
#Create a recipe
rec <- 
  recipe(satisfaction ~ . , data = wvs) %>%
  step_rm(iso3c)
```

How do we now use this in a model?

For this we need a **workflow** to which we can attach the recipe and the model. Creating a workflow of a linear regression looks like this:

```{r}
# Create the specification of a model but don't fit it yet
lm_model <- 
  linear_reg() %>%
  set_engine("lm")

# Create a workflow to add the recipe and the model
wflow <- 
  workflow() %>% 
  add_recipe(rec) %>% 
  add_model(lm_model)

print(wflow)
```

Then, you can fit the workflow to the data:

```{r}
model <- 
  wflow %>% 
  fit(data = wvs_train)

model
```

Do you want to **confirm** what recipe was used to pre-process the data? You can `extract_recipe()`:

```{r}
model %>% 
  extract_recipe()
```

You can use your model to `augment()` the data with predictions:

```{r}
model %>% 
  augment(new_data = wvs_test) %>% 
  head()
```

Notice that this augments our data frame with `.pred` and `.resid` columns.

How well did our model perform?

```{r}
reg_metrics <- metric_set(rsq, rmse, mae, mape)
model %>%
    augment(new_data=wvs_test)%>%
    reg_metrics(truth = satisfaction, estimate = .pred)
```

You already know about the $Rˆ2$ (`rsq`) and $RMSE$ (`rmse`) metrics: can you make sense of the other two metrics (`mae` for MAE and `mape` for MAPE)?

> Based on $R^2$ (the only scale-agnostic metric in our set of metrics), our model doesn't seem like a good fit. For the model to be a good fit, you would ideally want $R^2$ to be closer to 1. For more on $R^2$, see [this page](https://statisticsbyjim.com/regression/interpret-r-squared-regression/), [this one](https://people.duke.edu/~rnau/rsquared.htm) and finally, [this one](https://www.investopedia.com/terms/r/r-squared.asp).

> RMSE represents the square root of the variance of the residuals. The lower the RMSE, the better a given model fits your dataset. For our given dataset, the RMSE turns out to be 1.82, which looks relatively low and decent. However, RMSE is a metric that is much more useful when used to compare different models on the same dataset: the interpretation of a single RMSE measurement is rather tricky, however comparing the RMSE of different models on the same dataset can tell you which one is a better fit. For more information on RMSE, see [this page](https://www.statology.org/how-to-interpret-rmse/) and [this page](https://statisticsbyjim.com/regression/root-mean-square-error-rmse/).

> MAPE (or mean average percentage error) for our model is 32.8. This means that, on average, the error between the life satisfaction predicted by the model and the actual satisfaction is 32.8%. Do you think this is a low or high number? As usual with these sorts of metrics, their value lies more in comparisons with other models: a model with a lower MAPE than another is likely to be a better fit. For more on MAPE, see [this page](https://www.statology.org/how-to-interpret-mape/) or [this page](https://medium.com/@pirthipalsingh138/mastering-mape-a-guide-to-understanding-and-using-mean-absolute-percentage-error-8fd88f347eaa) or [this one](https://www.aporia.com/learn/a-comprehensive-guide-to-mean-absolute-percentage-error-mape/)

> The last metric is MAE (mean absolute error). It simply tells you on average, how far off the model is from the actual data. The actual formula for this metric is 
$$\frac{1}{n}\sum_{i=1}^n|y_i-\hat{y_i}|$$

> This metric is expressed in the scale of the variable you are predicting (i.e if you are predicting a salary in GBP, MAE is in GBP). The lower this metric, the better. It treats all error values equally (i.e it doesn't emphasize low or high errors). However, you can't use it to compare models trained on different datasets since it reflects an absolute error. For more on MAE, see [this page](https://docs.kolena.com/metrics/mean-absolute-error/) or [this page](https://medium.com/@ooemma83/interpretation-of-evaluation-metrics-for-regression-analysis-mae-mse-rmse-mape-r-squared-and-5693b61a9833) or [this one](https://www.statology.org/mae-vs-rmse/).


(Bonus) See if you can plot the metrics for a few countries of your choice.

Let's get my metrics per countries (for the sake of example, I am choosing the US (country code: USA), the UK (country code: GBR), Morocco (country code: MAR) and the Netherlands (country code: NLD)):
```{r}
model %>%
    augment(new_data=wvs_test)%>%
      filter(iso3c %in% c('GBR','MAR','USA','NLD'))%>%
      group_by(iso3c)%>%
      reg_metrics(truth = satisfaction, estimate = .pred)
```

```{r}
metrics_per_country <- model %>%
                             augment(new_data=wvs_test)%>%
                             filter(iso3c %in% c('GBR','MAR','USA','NLD'))%>%
                             group_by(iso3c)%>%
                             reg_metrics(truth = satisfaction, estimate = .pred)
```

Our various metrics have different scales so let's plot them separately in bar charts. 

Let's start with $R^2$
```{r}
metrics_per_country %>%
          rename(country_code=iso3c,Value=.estimate,Metric=.metric)%>%
          filter(Metric=='rsq')%>%
          ggplot(aes(Metric,Value,fill=country_code))+
          geom_bar(stat="identity", position = "dodge")+
          scale_fill_brewer(palette='Pastel1')
```

Then RMSE:

```{r}
metrics_per_country %>%
          rename(country_code=iso3c,Value=.estimate,Metric=.metric)%>%
          filter(Metric=='rmse')%>%
          ggplot(aes(Metric,Value,fill=country_code))+
          geom_bar(stat="identity", position = "dodge")+
          scale_fill_brewer(palette='Pastel1')
```

Then MAE:

```{r}
metrics_per_country %>%
          rename(country_code=iso3c,Value=.estimate,Metric=.metric)%>%
          filter(Metric=='mae')%>%
          ggplot(aes(Metric,Value,fill=country_code))+
          geom_bar(stat="identity", position = "dodge")+
          scale_fill_brewer(palette='Pastel1')
```

And finally MAPE:

```{r}
metrics_per_country %>%
          rename(country_code=iso3c,Value=.estimate,Metric=.metric)%>%
          filter(Metric=='mape')%>%
          ggplot(aes(Metric,Value,fill=country_code))+
          geom_bar(stat="identity", position = "dodge")+
          scale_fill_brewer(palette='Pastel1')
```

If we'd stopped by looking at $R^2$, we would have concluded the model performs best for Morocco. However, the other metrics, particularly MAPE, paint a different picture: errors are, on average, smaller for the Netherlands. The picture for Morocco diverges so much in between $R^2$ and MAPE that it might be worth having a look at the residual plot here to ascertain the underlying picture. Did we stumble upon any of the limitations of these two metrics in this particular case?

Let's look at the residuals plot for our model!

```{r}
model %>% 
  augment(new_data = wvs_test ) %>%
  ggplot(aes(.pred, .resid)) +
  geom_point() +
  geom_hline(yintercept = 0, linetype = "dashed") +
  theme_bw() +
  labs(x = "Predicted life satisfaction", y = "Residuals")
```

It's the same residuals plot you saw in the W03 lab! No surprises there since we fit the same model 😉

## Part III - Practice with recipes and workflows

0.  Consider the following simulated data frame.

```{r}

set.seed(123)

df_sim <-
  tibble(outcome = rnorm(100, 0, 1),
         feature1 = rnorm(100, 0, 1),
         feature2 = rnorm(100, 0, 20))

```

Let's have a quick look at the first lines of `df_sim` before we proceed:

```{r}
df_sim  %>% head(10) %>% knitr::kable()
```

We have created a data set with one outcome (creatively labeled `outcome`) and two features `feature1` and `feature2` where one variable has been assigned a standard deviation 20 times the size of the other variable. In some cases like linear regression, it will be okay to use `outcome ~ feature1 + feature2` as a formula. However, with other models, variables with larger variation will play a disproportionate role. As a result, we can use standardization (a.k.a. normalization) to counter this issue. This involves subtracting the mean of the variable and dividing the result by the standard deviation. Create a recipe for this data frame using a step from the `tidymodels` manual, and make use of graphs to demonstrate that the standardization process "worked". 

```{r}

# Let's look at the distributions of each feature overlaid!

df_sim %>% 
  pivot_longer(-outcome, names_to = 'feature') %>% 
  ggplot(aes(value, fill = str_remove(feature, "feature"))) +
  geom_histogram(alpha = 0.75) +
  theme_bw() +
  labs(x = "Feature value", y = "Frequency", fill = "Feature")

```
```{r}
# Now let's normalize feature2 using step_normalize
df_normalized <-
  df_sim %>% 
  recipe(outcome ~ ., data = .) %>% 
  step_normalize(feature2) %>% 
  prep() %>% 
  bake(new_data = NULL)
```

```{r}
# Now let's plot the distributions of the features
df_normalized %>% 
  pivot_longer(-outcome, names_to = 'feature') %>% 
  ggplot(aes(value, fill = str_remove(feature, "feature"))) +
  geom_histogram(alpha = 0.75) +
  theme_bw() +
  labs(x = "Feature value", y = "Frequency", fill = "Feature")
```

1.  Some models will only require numeric input, meaning that variables like `employment` or `better_living` into several dummy variables. Using the `tidymodels` manual, can you find a step that can help us achieve this end?
   
```{r}
recipe2 <- 
  recipe(satisfaction ~ . , data = wvs) %>%
  step_rm(iso3c) %>%
  step_dummy(employment,better_living)  
```

I can stop here if all I am doing is using the recipe for modeling (which is what we'll do in Question 2).
To be able to visualise the result of my pre-processing and check I've performed all the correct steps, here is an extra piece of code:

```{r}
  recipe2 %>%
        prep() %>%
        bake(wvs_train)
```

If I wanted to view my results even better, I could do this (try this code in your console - don't render!):

```r
  baked_recipe <- recipe2 %>% 
                  prep() %>%
                  bake(wvs_train)
   View(baked_recipe)
```

2.  Can you tweak the model we wrote in Part I to include the dummy step from Question 1? Evaluate the performance of this new model (first on the training then the test set).

The model specification from Part I doesn't change: we're still using a linear regression!
```{r}
# Create the specification of a model but don't fit it yet
lm_model <- 
  linear_reg() %>%
  set_engine("lm")
```

We tweak our workflow slightly to use the recipe `recipe2` that we've just defined in Question 2 (and that uses the dummy steps!) 
```{r}
# Create a workflow to add the recipe and the model
wflow <- 
  workflow() %>% 
  add_recipe(recipe2) %>% 
  add_model(lm_model)

print(wflow)
```

Then, you can fit the workflow to the data:

```{r}
model <- 
  wflow %>% 
  fit(data = wvs_train)

model
```

After that, you can use your model to `augment()` the data with predictions and compute metrics just as you did before in Part I:

```{r}
model %>% 
  augment(new_data = wvs_test) %>% 
  head()
```


```{r}
reg_metrics <- metric_set(rsq, rmse, mae, mape)
model %>%
    augment(new_data=wvs_test)%>%
    reg_metrics(truth = satisfaction, estimate = .pred)
```
  
How well did our model perform?

The dummy steps don't seem to have made a difference in our model's performance.

3.  Write a function called `plot_residuals()` that takes a fitted workflow model and a data frame in its original form (that is, before any pre-processing) and plot the residuals against the fitted values. The function should return a ggplot object.


```{r}
  plot_residuals <- function(wflow_model, data) {
    plot_df <-
      wflow_model %>% # This will be the model fitted with
                      # the workflow, not the workflow 
                      # itself!
      augment(new_data = data)

    g <- 
      ggplot(plot_df, aes(x=.pred, y=.resid)) +
      geom_point(alpha = 0.75) +
      geom_hline(yintercept = 0, linetype = "dashed") +
      theme_bw() +
      labs(x = "Predicted life satisfaction", y = "Residuals")

    return(g)
    
}
```

Using the `plot_residuals` function you've just created, can you create:

 (a) one residuals against the fitted values plot for the training data used in Question 2 
 (b) and another, separate residuals against the fitted values plot for the test data used in Question 2?

```{r}
#(a)
plot_residuals(model, wvs_train)
```
This is the residual plot we plotted previously in Part II!
```{r}
#(b)
plot_residuals(model, wvs_test)
```

4. (Bonus) Can you train a LASSO model on the data and evaluate its performance using recipes and workflows? How does this model compare to linear regression?
  

Let's modify `recipe2` slightly to make sure that all variables within our dataset are numeric

```{r}
recipe3 <- 
  recipe(satisfaction ~ . , data = wvs) %>%
  step_rm(iso3c) %>%
  step_dummy(employment,better_living)%>%
  step_mutate(across(c(social_trust,male,post_second_edu,rural,married,no_food,no_safety,no_medical,no_cash,no_shelter), as.numeric))
```  

Now, we need to modify our model specification since we are no longer using a linear model but a LASSO model.
```{r}
lasso_model <-
  linear_reg(penalty = 0.01, mixture = 1) %>%
  set_engine("glmnet") 
```

We add recipe and model specification to a workflow just as before:

```{r}
wflow_lasso <- 
  workflow() %>% 
  add_recipe(recipe3) %>% 
  add_model(lasso_model)

print(wflow_lasso)
```

```{r}
model_lasso <- 
  wflow_lasso %>% 
  fit(data = wvs_train)

model_lasso
```

After that, you can use your model to `augment()` the data with predictions and compute metrics just as you did before with the linear model:

```{r}
model_lasso %>% 
  augment(new_data = wvs_test) %>% 
  head()
```


```{r}
reg_metrics <- metric_set(rsq, rmse, mae, mape)
model_lasso %>%
    augment(new_data=wvs_test)%>%
    reg_metrics(truth = satisfaction, estimate = .pred)
```
  
How well did our model perform?
 
Overall, our metrics (in particular MAPE) show LASSO to performing slightly worse than linear regression for our particular task.

