Lab 3: Linear regression as a machine learning algorithm

Author

The DS202 Team

Welcome to the third DS202A lab! This week is all about regression, in particular linear regression. So let’s get started 😉

🥅 Learning Objectives

Get exposed to the concepts of training and test sets
Learning how to train and test a linear regression model the machine learning way
Learning how to train and test a LASSO model

📚 Preparation

Click here to read about how to prepare for this lab.

Then use the link below to download the lab materials:

as well as the dataset we will use for this lab (if you haven’t already done so when doing the lab prep!):

Put this dataset file in the data folder within your DS202A project folder.

📋 Lab Tasks

No need to wait! Start reading the tasks and tackling the action points below when you come to the classroom.

Loading libraries and functions


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


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.


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


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


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


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.


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.


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 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.


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.

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

# 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.

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

# 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.


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.

# 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

# Code here

Residual plots


# 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.

# 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


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

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


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


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


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


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.


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.

# 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

# 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.

# 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?

# 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.