---
title: "Linear regression: Scottish hill races"
output: rmarkdown::html_vignette
vignette: >
%\VignetteIndexEntry{Linear regression: Scottish hill races}
%\VignetteEngine{knitr::rmarkdown}
%\VignetteEncoding{UTF-8}
---
Choosing predictors by cross-validation, reading residuals, and what "linear" actually
constrains in a linear model.
The other vignettes are listed by `vignette (package = "fdm2id")`; they use the same
handful of functions on other data, and can be read in any order.
``` r
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.
``` r
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
```
``` r
plotdata (hills)
```
# Question 1. Which single variable predicts best? Was it foreseeable?
``` r
# 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.*
``` r
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?
``` r
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?
``` r
model = LINREG (hills [, -3], hills [, 3])
resplot (model)
```
``` r
resplot (model, index = 0)
```
``` r
resplot (model, index = 1)
```
``` r
resplot (model, index = 2)
```
**Answer.** *Two races stand out: one of the flattest, and the steepest of them all.*
``` r
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:
``` r
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.*
