library(ISLR)
library(MASS)
library(tidyverse)Supervised learning: regression in R
Introduction
In this practical, you will learn how to perform regression analysis in R, how to create predictions, how to plot with confidence and prediction intervals, how to calculate MSE, perform train-test splits, and write a function for cross validation.
Just like in the practical at the end of chapter 3 of the ISLR book, we will use the Boston dataset, which is in the MASS package that comes with R.
Make sure to load MASS before tidyverse otherwise the function MASS::select() will overwrite dplyr::select()
Regression in R
Regression is performed through the lm() function. It requires two arguments: a formula and data. A formula is a specific type of object that can be constructed like so:
some_formula <- outcome ~ predictor_1 + predictor_2 You can read it as “the outcome variable is a function of predictors 1 and 2”. As with other objects, you can check its class and even convert it to other classes, such as a character vector:
class(some_formula)
as.character(some_formula)You can estimate a linear model using lm() by specifying the outcome variable and the predictors in a formula and by inputting the dataset these variables should be taken from.
You have now trained a regression model with medv (housing value) as the outcome/dependent variable and lstat (socio-economic status) as the predictor / independent variable.
Remember that a regression estimates \(\beta_0\) (the intercept) and \(\beta_1\) (the slope) in the following equation:
\[\boldsymbol{y} = \beta_0 + \beta_1\cdot \boldsymbol{x}_1 + \boldsymbol{\epsilon}\]
We now have a model object lm_ses that represents the formula
\[\text{medv}_i = 34.55 - 0.95 * \text{lstat}_i + \epsilon_i\]
With this object, we can predict a new medv value by inputting its lstat value. The predict() method enables us to do this for the lstat values in the original dataset.
5. Create a scatter plot with
y_pred mapped to the x position and the true y value (Boston$medv) mapped to the y value. What do you see? What would this plot look like if the fit were perfect?
tibble(pred = y_pred,
obs = Boston$medv) %>%
ggplot(aes(x = pred, y = obs)) +
geom_point() +
theme_minimal() +
geom_abline(slope = 1)
# I've added an ideal line where all the points would lie on if the
# fit were perfect.We can also generate predictions from new data using the newdata argument in the predict() method. For that, we need to prepare a data frame with new values for the original predictors.
Plotting lm() in ggplot
A good way of understanding your model is by visualizing it. We are going to walk through the construction of a plot with a fit line and prediction / confidence intervals from an lm object.
8. Create a scatter plot from the
Boston dataset with lstat mapped to the x position and medv mapped to the y position. Store the plot in an object called p_scatter.
p_scatter <-
Boston %>%
ggplot(aes(x = lstat, y = medv)) +
geom_point() +
theme_minimal()
p_scatter
Now we’re going to add a prediction line to this plot.
10. Add a geom_line() to
p_scatter, with pred_dat as the data argument. What does this line represent?
p_scatter + geom_line(data = pred_dat)
# This line represents predicted values of medv for the values of lstat
13. Add a
geom_ribbon() to the plot with the data frame you just made. The ribbon geom requires three aesthetics: x (lstat, already mapped), ymin (lower), and ymax (upper).
Add the ribbon first, then the geom_line() and the geom_points() from before to make sure they remain visible. Give it a nice colour and clean up the plot, too!
# Create the plot
Boston %>%
ggplot(aes(x = lstat, y = medv)) +
geom_ribbon(aes(ymin = lower, ymax = upper), data = gg_pred, fill = "#00008b44") +
geom_point(colour = "#883321") +
geom_line(data = pred_dat, colour = "#00008b", size = 1) +
theme_minimal() +
labs(x = "Proportion of low SES households",
y = "Median house value",
title = "Boston house prices")
15. Do the same thing, but now with the prediction interval instead of the confidence interval.
# pred with pred interval
y_pred_95 <- predict(lm_ses, newdata = pred_dat, interval = "prediction")
# create the df
gg_pred <- tibble(
lstat = pred_dat$lstat,
medv = y_pred_95[, 1],
l95 = y_pred_95[, 2],
u95 = y_pred_95[, 3]
)
# Create the plot
Boston %>%
ggplot(aes(x = lstat, y = medv)) +
geom_ribbon(aes(ymin = l95, ymax = u95), data = gg_pred, fill = "#00008b44") +
geom_point(colour = "#883321") +
geom_line(data = pred_dat, colour = "#00008b", size = 1) +
theme_minimal() +
labs(x = "Proportion of low SES households",
y = "Median house value",
title = "Boston house prices")
# You can look at ISLR p.82-83 for a discussion of prediction intervalsMean squared error
You have now calculated the mean squared length of the dashed lines below.
You have calculated the mean squared length of the dashed lines in the plot below.
Train-validation-test split
Now we will use the sample() function to randomly select observations from the Boston dataset to go into a training, test, and validation set. The training set will be used to fit our model, the validation set will be used to calculate the out-of sample prediction error during model building (this can also be called hyperparameter tuning), and the test set will be used to estimate the true out-of-sample MSE.
We will set aside the boston_test dataset for now.
This is the estimated out-of-sample mean squared error.
OPTIONAL Programming exercise: cross-validation
This is an advanced exercise. Some components we have seen before in this and previous practicals, but some things will be completely new. Try to complete it by yourself, but don’t worry if you get stuck. If you don’t know about for loops in R, read up on those before you start the exercise.
Use help in this order:
- R help files
- Internet search & stack exchange
- Your peers
- The answer, which shows one solution
You may also just read the answer and try to understand what happens in each step.