---
title: "📓 LSE DS202A 2025: Week 03 - Student Notebook"
author: "MY NAME (MY CANDIDATE NUMBER)"
date: 15 September 2025
format: html
self-contained: true
---

## 🎯 Learning Outcomes

By the end of this lab, you will be able to:

1.  **Distinguish between hypothesis testing and machine learning approaches** - Understand the fundamental difference between using statistical models to test causal hypotheses versus building predictive models for accurate forecasting on unseen data.

2.  **Implement proper machine learning workflow practices** - Split datasets into training and test sets using `tidymodels`, build models on training data, and evaluate performance using appropriate metrics (R-squared, RMSE) on held-out test data.

3.  **Apply and interpret penalized regression techniques** - Use lasso regression with regularization penalties to perform automatic feature selection, understand how penalty parameters affect coefficient shrinkage, and compare model performance across different penalty values.

4.  **Evaluate model performance and identify prediction errors** - Create and interpret residual plots to diagnose model fit, calculate prediction errors, and systematically compare multiple modeling approaches to determine the most effective predictive strategy.

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

## Before we do anything more

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

## Student Performance

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

students <- read_rds("data/students-data-cleaned.rds")

```

We start our machine learning journey with a student performance dataset (`students`), which contains information on students from two Portuguese schools. We have cleaned the data to only include statistically significant predictors. The columns include:

-   `final_grade` final grade from 0 to 20 (**the outcome**)
-   `school` student's school (binary: 'GP' - Gabriel Pereira or 'MS' - Mousinho da Silveira, reference: GP)
-   `sex` student's sex (binary: 'F' - female or 'M' - male, reference: female)
-   `age` student's age (numeric: from 15 to 22 years)
-   `studytime` weekly study time (categorical: '\<2hrs', '2-5hrs', '5-10hrs', '\>10hrs', reference: \<2hrs)
-   `failures` number of past class failures (categorical: '0', '1', '2', '3', '4+', reference: 0 failures)
-   `schoolsup` extra educational support (binary: yes or no, reference: no)
-   `higher` wants to take higher education (binary: yes or no, reference: no)
-   `goout` going out with friends frequency (categorical: 'VeryLow', 'Low', 'Average', 'High', 'VeryHigh', reference: VeryLow)
-   `dalc` workday alcohol consumption (categorical: 'VeryLow', 'Low', 'Average', 'High', 'VeryHigh', reference: VeryLow)
-   `health` current health status (categorical: 'VeryBad', 'Bad', 'Average', 'Good', 'VeryGood', reference: VeryBad)
-   `romantic` in a romantic relationship (binary: yes or no, reference: no)

## Understanding student performance: some exploratory data analysis (EDA) (10 minutes)

Now it's your turn to explore the `students` dataset! Use this time to create visualizations and discover patterns in the data that might help explain what drives student success.

### Some ideas to get you started:

-   How are final grades distributed? Are they normally distributed or skewed?
-   Do students who want higher education perform better than those who don't?
-   Is there a relationship between study time and final grades?
-   How does past academic failure affect current performance?
-   Are there differences in performance between the two schools?
-   Does going out with friends impact academic performance?
-   What's the relationship between health status and grades?
-   Do students receiving extra educational support perform differently?

### Challenge yourself:

-   Can you find any surprising relationships in the data?
-   What patterns emerge when you look at combinations of variables?
-   Are there any outliers or interesting edge cases?

**Share your most interesting findings on Slack!** We'd love to see what patterns you discover and which visualizations tell the most compelling stories about student performance.

```{r}
# Your EDA code here!
# Try different plot types: histograms, boxplots, scatter plots, bar charts
# Example starter code:

# students |>
#   ggplot(aes(x = final_grade)) +
#   geom_histogram()

```

## Understanding student performance: the hypothesis-testing approach (5 minutes)

Why do some students perform better 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 students who want to pursue higher education have better academic performance. We can estimate a linear regression model by using `final_grade` as the dependent variable and `higher` 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 `higher`.

```{r}

lm(final_grade ~ higher, data = students) |>
  summary()

```

We see that students who want to pursue higher education have a positive and statistically significant (p \< 0.001) increase in final grades of about 3.5 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(final_grade ~ ., data = students) |>
  summary()

```

The `.` placeholder indicates that we want to use all other variables in the data set. We can see that the coefficient estimate for `higher` approximately halves in size yet remains positive and highly significant (p \< 0.001), suggesting this relationship holds even when controlling for other factors like study time, past failures, and demographic characteristics.

👉 **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 student grades: 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}

sim_model_preds(10, 20, 30)

```

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

students_split <- initial_split(students, prop = 0.75)

# Create tibbles of the training and test set

students_train <- training(students_split)
students_test <- testing(students_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 `students_train`. We will also only create a model object (`mv_model`, short for multivariate model).

```{r}

mv_model <- lm(final_grade ~ ., data = students_train)

```

### 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 = students_test) |>
  reg_metrics(truth = final_grade, 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?

## Challenge: Running multiple univariate regressions efficiently (30 minutes)

🎯 **YOUR CHALLENGE:** Can you figure out how to run a univariate model for all features in the student dataset without creating 11 separate model objects?

Remember our univariate model earlier where we looked at `final_grade ~ higher`? Now we want to do the same thing for ALL features to see which individual variables are the strongest predictors of student performance.

The naive approach would be to copy-paste code 11 times and create 11 different model objects, but that's inefficient and error-prone. Your job is to find a more elegant solution!

### Your mission:

**Figure out how to systematically run univariate regressions for each predictor variable and compare their performance.**

Some things to think about: - How can you automate the process of fitting multiple models? - How will you store and compare the results? - What's the best way to visualize which features are most predictive?

### Work in pairs and brainstorm different approaches!

You might consider: - Loops (for loops, while loops) - Functional programming (map functions) - Advanced tidyverse techniques (nested tibbles) - Or something completely different!

### Questions to explore:

-   Which single variable is the strongest predictor of student performance?
-   Are there any surprisingly weak predictors?
-   How do these univariate results compare to what you found significant in the multivariate model?
-   Do the "best" univariate predictors match what you'd expect from theory?

### Get creative and experiment!

There's no single "right" way to solve this problem. Try different approaches, see what works, and don't be afraid to ask for help when you get stuck.

**Share your approach and findings on Slack!** We want to see the different methods you came up with and which features you found most predictive.

```{r}

# Code here

```

## 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(final_grade ~ ., data = students_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?

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

*Let's employ a nested tibbles approach to find better penalty values.*

```{r}

# Code here

```

*Which penalty value gives you the lowest RMSE? How do the coefficients change as you increase the penalty?*

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