Linear regression: Scottish hill races

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)
#>       dist            climb           time       
#>  Min.   : 2.000   Min.   : 300   Min.   : 15.95  
#>  1st Qu.: 4.500   1st Qu.: 725   1st Qu.: 28.00  
#>  Median : 6.000   Median :1000   Median : 39.75  
#>  Mean   : 7.529   Mean   :1815   Mean   : 57.88  
#>  3rd Qu.: 8.000   3rd Qu.:2200   3rd Qu.: 68.62  
#>  Max.   :28.000   Max.   :7500   Max.   :204.62
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"))
#>       adjr2        msep 
#>   0.7835907 510.9387160
performance (LINREG, hills [, 2], hills [, 3], protocol = "loocv", eval = c ("adjr2", "msep"))
#>       adjr2        msep 
#>   0.6035177 936.0880207

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)
#>        dist climb  time
#> dist  1.000 0.652 0.920
#> climb 0.652 1.000 0.805
#> time  0.920 0.805 1.000

Question 2. One variable or two?

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

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)
#>   Knock Hill Bens of Jura    Ben Nevis 
#>     65.12140     31.26242     16.21532
hills [c ("Knock Hill", "Bens of Jura"), ]
#>              dist climb    time
#> Knock Hill      3   350  78.650
#> Bens of Jura   16  7500 204.617

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"))
#>       adjr2        msep 
#>   0.9263591 163.3276285

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.