Linear regression: Scottish hill races
In fdm2id: Data Mining and R Programming for Beginners

knitr::opts_chunk$set (collapse = TRUE, comment = "#>", fig.width = 6, fig.height = 4.2,
                       fig.align = "center")
optional = c ("MASS")
available = all (sapply (optional, requireNamespace, quietly = TRUE))
knitr::opts_chunk$set (eval = available)

cat ("**Note.** This vignette needs the following packages, some of which are missing:",
     paste (optional, collapse = ", "), "-- the code is shown but not run.\n")

One of the case studies of the Analyse de données (L3 Informatique) course, for which fdm2id was written. Choosing predictors by cross-validation, reading residuals, and what "linear" actually constrains in a linear model.

The other case studies are listed by vignette (package = "fdm2id"); they use the same handful of functions on other data, and can be read in any order.

library (fdm2id)

The data

Records of hill races in Scotland. Each race is described by three variables: its distance (in miles), its total climb (in feet), and the record time set in 1984 (in minutes). The question is whether the time can be predicted from the distance and the climb.

data (hills, package = "MASS")
summary (hills)

plotdata (hills)

Question 1. Which single variable predicts best? Was it foreseeable?

# Reproducible without a seed: leave-one-out builds n folds of one observation each, so there
# is nothing to draw. It is the one protocol of the package that needs no 'seed'.
performance (LINREG, hills [, 1], hills [, 3], protocol = "loocv", eval = c ("adjr2", "msep"))
performance (LINREG, hills [, 2], hills [, 3], protocol = "loocv", eval = c ("adjr2", "msep"))

Answer. Under leave-one-out cross-validation, dist predicts more precisely than climb. It was foreseeable: its linear correlation with time is the higher of the two.

round (cor (hills), 3)

Question 2. One variable or two?

performance (LINREG, hills [, -3], hills [, 3], protocol = "loocv", eval = c ("adjr2", "msep"))

Answer. Two. The mean squared error of prediction drops by nearly half.

Question 3. What do the residuals say?

model = LINREG (hills [, -3], hills [, 3])
resplot (model)
resplot (model, index = 0)
resplot (model, index = 1)
resplot (model, index = 2)

Answer. Two races stand out: one of the flattest, and the steepest of them all.

head (sort (abs (residuals (model$model)), decreasing = TRUE), 3)
hills [c ("Knock Hill", "Bens of Jura"), ]

Knock Hill is three miles with 350 feet of climb and a record of 78 minutes, which is not a record but a recording error -- the accepted reading is 18 minutes. Bens of Jura is the longest and steepest race in the table, and the model has nothing else like it to learn from.

Question 4. Can the predictions be improved?

The residuals plotted against climb are vaguely parabolic, which suggests adding a climb² variable:

hills2 = cbind (hills, hills$climb^2)
colnames (hills2) = c (colnames (hills), "climb2")
performance (LINREG, hills2 [, -3], hills2 [, 3], protocol = "loocv",
             eval = c ("adjr2", "msep"))

Answer. Under leave-one-out cross-validation this improves the predictions markedly -- the error falls again by more than a third. A linear model is linear in its coefficients, not in the variables it is given.