---
title: "LSE DS202A 2024: Week 03 - Lab"
subtitle: "Linear regression as a machine learning algorithm"
author: "MY NAME (MY CANDIDATE NUMBER)"
date: 14 Oct 2024
format: html
editor: visual
self-contained: true
---

## Loading libraries and functions

```{r, warning=FALSE,message=FALSE}
library("ggsci")       # For pretty colour schemes
library("MASS")        # For simulating data  
library("scales")      # For number formatting in ggplot2
library("tidymodels")  # For train / test splits
library("tidyverse")   # For data wrangling / visualisation

# Let's use some of the theme functions we created in lab 2

theme_histogram <- function() {
  
  theme_minimal() +
  theme(panel.grid.minor = element_blank(),
        panel.grid.major.x = element_blank())
  
}

theme_boxplot <- function() {
  
  theme_minimal() +
  theme(panel.grid.minor = element_blank(),
        panel.grid.major.x = element_blank(),
        legend.position = "none")
  
}

theme_bar <- function() {
  
  theme_minimal() +
  theme(panel.grid.minor = element_blank(),
        panel.grid.major.y = element_blank())
  
}

theme_scatter <- function() {
  
  theme_minimal() +
  theme(panel.grid.minor = element_blank(),
        legend.position = "bottom")
  
}

```

## Before we do anything more

Please create a data folder called `data` to store all the different data sets in this course.

## The World Values Survey

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

wvs <- read_csv("data/WVS_Wave_7.csv")

```

We start our machine learning journey with Wave 7 of the World Values Survey (`wvs`), which contains information on `r nrow(wvs)` from `r length(unique(wvs$iso3c))` countries. We have cleaned the data to only include non-missing values. The columns include:

-   `iso3c` 3-letter country iso code
-   `satisfaction` 1 to 10 rating of life satisfaction (**the outcome**)
-   `social_trust` TRUE / FALSE as to whether or not someone expresses social trust
-   `male` respondent is male (reference: female)
-   `age` age of respondent
-   `post_second_edu` respondent has post-secondary education
-   `rural` respondent lives in a rural area (reference: urban)
-   `employment` categorical employment variable (try `count(wvs, employment)`)
-   `financial_situ` 1 to 10 rating of respondent's financial situation
-   `married` the respondent is married (reference: other marital status)
-   `relig_import` 1 to 10 rating of how much importance the respondent attaches to religion.
-   `better_living` trichotomous rating of whether or not the respondent sees their lives as better off than their parents (try `count(wvs, better_living)`)
-   `no_food`: In the last 12 months, how often have your or your family gone without enough food to eat? (Sometimes / Often = TRUE, Rarely / Never = FALSE)
-   `no_safety` In the last 12 months, how often have your or your family felt unsafe from crime in your home? (Sometimes / Often = TRUE, Rarely / Never = FALSE)
-   `no_medical` In the last 12 months, how often have your or your family gone without medicine or medical treatment that you needed? (Sometimes / Often = TRUE, Rarely / Never = FALSE)
-   `no_cash` In the last 12 months, how often have your or your family gone without a cash income? (Sometimes / Often = TRUE, Rarely / Never = FALSE)
-   `no_shelter` In the last 12 months, how often have your or your family gone without a safe shelter over your head? (Sometimes / Often = TRUE, Rarely / Never = FALSE)

👉 **NOTE:** With variables such as `no_food`, we recoded the original question mainly for the sake of simplicity. This simplification can be useful as it reduces the number of parameters in our model. However, there may be distinct differences between each level that may produce different results.

If you find that your laptop is unable to handle the full data set without running slowly, try experimenting with the following code. 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, eval=FALSE}

# 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()

```

## Understanding life satisfaction: some exploratory data analysis (EDA) (5 minutes)

Here are a couple of graphs that tell a few stories about `wvs`.

### Median life satisfaction appears to be relatively high in the sample (7 out of 10)

```{r}

med_satisfaction <- median(wvs$satisfaction)

wvs %>% 
  ggplot(aes(satisfaction)) +
  geom_histogram(fill = "midnightblue", colour = "black", bins = 10, alpha = 0.5) +
  geom_vline(xintercept = med_satisfaction, linetype = "dashed", size = 2, colour = "red") +
  theme_histogram() +
  scale_x_continuous(breaks = c(1, med_satisfaction, 10)) +
  scale_y_continuous(labels = comma) +
  labs(x = "Life satisfaction", y = "Number of respondents",
       caption = "Note: dotted line represents median life satisfaction")

```

### Life satisfaction tracks positively with financial situation

```{r}

wvs %>% 
  count(financial_situ, satisfaction) %>% 
  ggplot(aes(financial_situ, satisfaction, size = n)) +
  geom_point() +
  scale_x_continuous(breaks = seq(2, 10, 2)) +
  scale_y_continuous(breaks = seq(2, 10, 2)) +
  scale_size_continuous(labels = comma) +
  theme(panel.grid = element_blank(),
        panel.background = element_rect(fill = "white"),
        legend.position = "bottom") +
  labs(x = "Financial situation", y = "Life satisfaction")

```

### The median individual with post-secondary education has a higher life satisfaction than the median individual without

```{r}

wvs %>% 
  ggplot(aes(post_second_edu, satisfaction, fill = post_second_edu)) +
  geom_boxplot() +
  theme_boxplot() +
  scale_fill_jco() +
  labs(x = "Post-secondary education?", y = "Life satisfaction")

```

## Understanding life satisfaction: the hypothesis-testing approach (5 minutes)

Why do some people have a higher life satisfaction than others? This is one question that a quantitative social scientist might answer by exploring the magnitude and precision of a series of variables. Suppose we hypothesised that individuals with post-secondary level education have greater life satisfaction. We can estimate a linear regression model by using `satisfaction` as the dependent variable and `post_second_edu` as the independent variable.

To run a linear regression, we can use the `lm` function, which requires two things:

-   A model formula (a.k.a. equation)
-   The data used to estimate the model

Let's do this now. We can call the summary function to get information on the coefficient estimate for `post_second_edu`.

```{r}

lm(satisfaction ~ post_second_edu, data = wvs) %>% 
  summary()

```

We see that individuals with post-secondary education have a positive and statistically significant (p \< 0.001) increase in life satisfaction of about 0.2 points.

👉 **NOTE:** The process of hypothesis testing is obviously more involved when using observational data than is portrayed by this simple example. Control variables will almost always be incorporated and, increasingly, identification strategies will be used to uncover causal effects. The end result, however, will involve as rigorous an attempt at falsifying a hypothesis as can be provided with the data.

For an example of how multivariate regression is used, we can run the following code.

```{r}

lm(satisfaction ~ . -iso3c, data = wvs) %>% 
  summary()

```

The `.` placeholder indicates that we want to use all other variables in the data set. `-iso3c` indicates that we are omitting country dummy variables. Interestingly, we can see that the coefficient estimate for `post_second_edu` is now negative and significant, albeit at a lower level than before (p \< 0.05).

👉 **NOTE:** p-values are useful to machine learning scientists as they indicate which variables may yield a significant increase in model performance. However, p-hacking where researchers manipulate data to find results that support their hypothesis make it hard to tell whether or not a relationship held up after honest attempts at falsification. This can range from using a specific modelling approach that produces statistically significant (while failing to report others that do not) findings to outright manipulation of the data. For a recent egregious case of the latter, we recommend the [**Data Falsificada**](https://datacolada.org/archives/category/fake-data) series.

## Predicting life satisfaction: the machine learning approach (30 minutes)

Machine learning scientists take a different approach. Our aim, in this context, is to build a model that can be used to accurately predict how happy a person is using a mixture of features and, for some models, hyperparameters (which we will address in Lab 5).

*Thus, rather than attempting to falsify the effects of causes, we are more concerned about the fit of the model in the aggregate when applied to unforeseen data.*

To achieve this, we do the following:

-   Split the data into training and test sets
-   Build a model using the training set
-   Evaluate the model on the test set

Let's look at each of these in turn.

### Split the data into training and test sets

It is worth considering what a training and test set is and why we might split the data this way.

A **training set** is data that we use to build (or "train") a model. In the case of multivariate linear regression, we are using the training data to estimate a series of coefficients. Here is a made-up multivariate linear model with three coefficients derived from (non-existent) data to illustrate things.

```{r}

sim_model_preds <- function(x1, x2, x3) {
  
  y <- 1.1*x1 + 2.2*x2 + 3.3*x3
  y
  
}

```

A **test set** is data that the model has not yet seen. We then apply the model to this data set and use an evaluation metric to find out how accurate our predictions are. For example, suppose we had a new observation where `x1 = 10`, `x2 = 20` and `x3 = 30` and `y = 150`. We can use the above model to develop a prediction.

```{r}
# Code here
```

We get a prediction of `r sim_model_preds(10, 20, 30)` points!

We can also calculate the amount of error we make by calculating residuals (actual value - predicted value).

```{r}
# Code here
```

We can see that our model is `r abs(150 - sim_model_preds(10, 20, 30))` points off the real answer!

Why do we evaluate our models using different data? Because, as stated earlier, machine learning scientists care about the applicability of a model to *unforeseen* data. If we were to evaluate the model using the training data, we obviously cannot do this to begin with. Furthermore, we cannot ascertain whether the model we have built can generalise to other data sets or if the model has simply learned the idiosyncrasies of the data it was used to train on. We will discuss the concept of **overfitting** throughout this course.

We can use the `rsample` package in the `tidymodels` ecosystem to split the data into training and test sets.

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

set.seed(123)

# 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)
```

👉 **NOTE:** Our data are purely cross-sectional, so we can use this approach. However, when working with more complex data structures (e.g. time series cross sectional), different approaches to splitting the data will need to be used.

### Build a model using the training set

This is remarkably simple. We will use almost exactly the same code we used to build a multivariate linear model, but with one exception. Instead of using the whole of the data, we will only use `wvs_train`. We will also only create a model object (`mv_model`, short for multivariate model).

```{r}
# Code here
```

### Evaluate the model using the test set

Now that we have trained a model, we can then evaluate its performance on the test set. We will look at two evaluation metrics:

-   **R-squared**: the proportion of variance in the outcome explained by the model.
-   **Root mean squared error (RMSE)**: the amount of error a typical observation parameterised as the units used in the initial measurement.

```{r}

reg_metrics <- metric_set(rsq, rmse)

mv_model %>% 
  augment(new_data = wvs_test) %>% 
  reg_metrics(truth = satisfaction, estimate = .fitted)

```

🗣️ **CLASSROOM DISCUSSION:**

How can we interpret these results?

## Graphically exploring where we make errors

We are going to build some residual scatter plots which look at the relationship between the values fitted by the model for each observation and the residuals (actual - predicted values). Before we do this for our data, let's take a look at an example where there is a near perfect relationship between two variables. As this very rarely exists in the social world, we will rely upon simulated data.

We adapted this code from [here](https://www.r-bloggers.com/2021/05/how-to-generate-correlated-data-in-r/).

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

set.seed(123)

# Create the variance covariance matrix

sigma <- rbind(c(1,0.99), c(0.99,1))

# Create the mean vector

mu <- c(10, 5) 

# Generate a multivariate normal distribution using 1,000 samples

sim_data <-
  mvrnorm(n=1000, mu=mu, Sigma=sigma) %>% 
  as.data.frame() %>% 
  as_tibble()

```

### Plot the correlation

```{r}
# Code here
```

### Residual plots

```{r}

# Build a linear model and plot the fitted versus residual values

lm(V2 ~ V1, data = sim_data) %>% 
  augment() %>% 
  ggplot(aes(.fitted, .resid)) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  geom_point() +
  theme_scatter() +
  labs(x = "Fitted values", y = "Residuals")

```

Now let's run this code for our model.

```{r}
# Code here
```

🎯 **ACTION POINTS** why does the graph of the simulated data illustrate a more well-fitting model when compared to our actual data?

## Introduction to using nested tibbles to aid feature selection (30 minutes)

👨🏻‍🏫 **TEACHING MOMENT:** Your tutor will take you through the code, so sit back, relax and enjoy!

Remember our univariate model earlier? We are going to do the same for all features so see which ones show the best improvements in predictive power.

We could build 15 different model objects, but this would be very inefficient. Instead, we are going to take advantage of a unique feature of tibbles, the **list column**.

So far, we have looked at columns that are of numeric, integer, Boolean, factor and character class. However, with list columns, we can nest and unnest any class of object within a single cell. This means we can do things like apply functions over list columns simply by adding another list column to our tibble.

We are going to leverage this by creating a series of formulas and build a linear model using a combination of each formula and the nested training data.

### Create a series of univariate regression formulas

```{r}

formulas <- paste("satisfaction ~", colnames(wvs)[3:17])
formulas

```

### Create a tibble that combines these formulas with the training and test sets

```{r}

wvs_tbl <-
  # Use crossing to find all key combinations between the
  # formulas and the training and test sets
  crossing(formula = formulas,
           # We use the nest function to create a list column
           # for both the training and test sets
           nest(wvs_train, .key = "train_set"),
           nest(wvs_test, .key = "test_set")) 

wvs_tbl

```

### Build a linear model using the training set and apply it to the test set

```{r}

models <-
  wvs_tbl %>% 
         # Add the features (this will help with plotting!)
  mutate(feature = sort(colnames(wvs)[3:17]),
         # Build a linear model using the training set
         model = map2(formula, train_set, ~ lm(.x, data = .y)),
         # Apply the predictions to the test set
         augmented = map2(model, test_set, ~ augment(.x, new_data = .y)))

models

```

### Unnest the data frame and calculate the r-squared value for each univariate regression

```{r}

preds <- 
  models %>% 
  # To unnest a list column, we use the unnest function
  unnest(augmented) %>% 
  # We use group_by to perform grouped calculations of the r-squared
  # by feature.
  group_by(feature) %>% 
  rsq(truth = satisfaction, estimate = .fitted) %>% 
  # We reorder the features so the "best" / "worst" is on the top / 
  # bottom
  mutate(feature = fct_reorder(feature, .estimate))

preds

```

### Plot the results

```{r}

preds %>% 
  # We the use a bar plot (see Lab 2)
  ggplot(aes(.estimate, feature)) +
  geom_col() +
  theme_bar() +
  labs(x = "Test set r-squared", y = NULL)

```

## Using penalised linear regression to perform feature selection (20 minutes)

We are now going to experiment with a lasso regression which, in this case, is a linear regression that uses a so-called **hyperparameter** - a "dial" built into a given model that can be experimented with to improve model performance. The hyperparameter in this case is a regularisation penalty which takes the value of a non-negative number. This penalty can shrink the magnitude of coefficients down to zero and the larger the penalty, the more shrinkage occurs.

### Step 1: Create a lasso model

Run the following code. This builds a lasso model with the penalty parameter set to `0.01`.

```{r}

lasso_model <-
  linear_reg(penalty = 0.01, mixture = 1) %>%
  set_engine("glmnet") %>% 
  fit(satisfaction ~ . -iso3c, data = wvs_train)

```

### Step 2: Extract lasso coefficients

Use the `tidy` function on the lasso model to get the coefficients.

```{r}
# Code here
```

🎯 **ACTION POINTS** What is the output? Which coefficients have been shrunk to zero? What is the most important feature?

### Step 3: Create a bar plot

```{r}
# Code here
```

### Step 4: Evaluate on the test set

Although a different model is used, the code for evaluating the model on the test set is exactly the same as earlier.

```{r}
# Code here
```

🗣️ **CLASSROOM DISCUSSION:**

Does this model represent an improvement on the linear model?

### Step 5: Experiment with different penalties

This is your chance to try out different penalties. Can you find a penalty that improves test set performance?

```{r}
# Code here
```

👉 **NOTE:** In labs 4 and 5, we are going to use a method called k-fold cross validation to systematically test different combinations of hyperparameters for models such as the lasso.
