Linear models

Intercepts, slopes and residuals

Dr Massimiliano Canzi

From one variable to two

Did you know: an Alpine ibex hoof has a hard, thin rim around a soft, spongy sole, and the ibex can move its two toes independently, placing each one on a different edge of the rock.

Source: ALPECOLE, University of Zurich

Today

  • From the mean to the conditional mean
  • Intercepts, slopes and residuals
  • Fitting a line with least squares and lm()
  • Reading a model and making predictions
  • Measuring how well the line fits with R²
  • Checking the model, and a warning about groups

Last week: one variable

The mean (201 mm) is the simplest model of a variable: it predicts the same flipper length for every penguin. Most penguins are some distance away from it.

Today: two variables

If we know a penguin’s body mass, we can do better: the conditional mean is the mean flipper length of penguins with a given body mass (orange). The conditional means lie close to a straight line.

Outcome and predictor

Simple linear regression models one continuous outcome with one predictor (Winter, 2019).

  • The outcome (dependent or response variable, y) is what you measure
  • The predictor (independent variable, x) is what you manipulate, or measure alongside, to see whether the outcome changes with it
Outcome Predictor
Reaction time Word frequency
Vowel duration Speaking rate
Pitch (f0) Speaker age
Reading speed Font size

Intercepts & slopes

Did you know: Alpine ibex climb near-vertical slopes. In the Italian Alps, females scale the concrete walls of dams such as the Cingino dam to lick the salts that form on them.

Source: Water Power Magazine

The equation of a line

y = b₀ + b₁ · x + e
  • The intercept b₀ is where the line crosses the vertical axis, at x = 0
  • The slope b₁ is how much y changes when x goes up by 1
  • The error e is the part of each observation the line misses

Moving and tilting the line

The intercept moves the line up or down. The slope tilts it: positive slopes rise, negative slopes fall.

Which line?

All three lines go through the middle of the data. Which one describes it best, and how would you decide?

Residuals

residual = observed − fitted

Each pink segment is a residual: how far a penguin’s flipper is from the length the line predicts. Points above the line have positive residuals, points below negative ones.

Least squares

The least squares method picks the line with the smallest sum of squared residuals.

sum(residuals(m)^2)
[1] 16250.3
guess <- 140 + 0.014 * pengs$body_mass_g
sum((pengs$flipper_length_mm - guess)^2)
[1] 18102.13

Squaring makes every residual positive, so misses above and below the line cannot cancel out. It also counts one large miss as worse than several small ones.

Searching for the best line

Each point on the curve is one candidate line. The sum of squared residuals is smallest at a slope of 0.0153. That is line B from the previous slide, and it is the line lm() finds.

Linear models in R

Did you know: a male Alpine ibex weighs about 90 kg, and his horns grow to about 92 cm on average. A female weighs 35 to 49 kg, with horns of about 34 cm. A model of horn length would need very different intercepts for the two sexes.

Source: Gran Paradiso National Park

lm() in R

lm() fits a linear model. The formula outcome ~ predictor reads “flipper length as a function of body mass”.

m <- lm(flipper_length_mm ~ body_mass_g, data = pengs)
m

Call:
lm(formula = flipper_length_mm ~ body_mass_g, data = pengs)

Coefficients:
(Intercept)  body_mass_g  
  136.72956      0.01528

pengs is penguins without the rows that lack either measurement, made with drop_na() from last week.

summary()

summary(m)

Call:
lm(formula = flipper_length_mm ~ body_mass_g, data = pengs)

Residuals:
     Min       1Q   Median       3Q      Max 
-23.7626  -4.9138   0.9891   5.1166  16.6392 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 1.367e+02  1.997e+00   68.47   <2e-16 ***
body_mass_g 1.528e-02  4.668e-04   32.72   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 6.913 on 340 degrees of freedom
Multiple R-squared:  0.759, Adjusted R-squared:  0.7583 
F-statistic:  1071 on 1 and 340 DF,  p-value: < 2.2e-16

Reading summary()

Part of the output What it tells you
Residuals How far the points are from the line: the median should be close to 0
Coefficients, Estimate The intercept b₀ and the slope b₁
Std. Error, t value, Pr(>\|t\|) How certain the estimates are: lecture 9
Residual standard error The typical size of a residual, in mm
Multiple R-squared How much of the variation the model accounts for: later today

Reading the coefficients

coefficients(m)
 (Intercept)  body_mass_g 
136.72955927   0.01527592
  • Each extra kilogram of body mass goes with flippers 15.3 mm longer (0.0153 mm per gram)
  • A penguin weighing 0 g would have flippers 136.7 mm long. No penguin weighs 0 g, so next week centering will give the intercept a meaning

Residuals average zero

The points above the line balance the points below it, so the residuals average zero:

mean(residuals(m))
[1] 5.02522e-16

e-16 means “times 10 to the power of −16”, so this mean is zero up to rounding error.

{broom}: models as tibbles

tidy() gives one row per coefficient and glance() one row per model:

tidy(m)
# A tibble: 2 × 5
  term        estimate std.error statistic   p.value
  <chr>          <dbl>     <dbl>     <dbl>     <dbl>
1 (Intercept) 137.      2.00          68.5 5.71e-201
2 body_mass_g   0.0153  0.000467      32.7 4.37e-107
glance(m)[, 1:4]
# A tibble: 1 × 4
  r.squared adj.r.squared sigma statistic
      <dbl>         <dbl> <dbl>     <dbl>
1     0.759         0.758  6.91     1071.

augment(): one row per penguin

augment(m)[, 1:4]
# A tibble: 342 × 4
  flipper_length_mm body_mass_g .fitted .resid
              <int>       <int>   <dbl>  <dbl>
1               181        3750    194. -13.0 
2               186        3800    195.  -8.78
3               195        3250    186.   8.62
4               193        3450    189.   3.57
# ℹ 338 more rows

augment() adds the fitted value and the residual of every penguin to the data. All three {broom} functions return tibbles, ready for {dplyr} and {ggplot2}.

Plotting the model

ggplot(pengs, aes(x = body_mass_g, y = flipper_length_mm)) +
  geom_point() +
  geom_smooth(method = "lm")

geom_smooth(method = "lm") fits the same model as lm(). The band shows how uncertain the line is: lecture 9.

Predictions

new_penguins <- tibble(body_mass_g = c(3000, 4000, 5000))
predict(m, newdata = new_penguins)
       1        2        3 
182.5573 197.8332 213.1091

Each prediction is b₀ + b₁ · x. For 4000 g: 136.7 + 0.0153 × 4000 = 197.8 mm.

Stay within the data

The data cover 2700 to 6300 g (shaded). Outside that range the dashed line is a guess: it gives a weightless penguin 137 mm flippers.

Model fit

Did you know: an Alpine ibex horn adds one segment a year, so you can count a male’s age on it. In males aged 1 to 6, age accounted for 50% of the variance in yearly horn growth. Warm springs with early plant growth gave longer segments.

Sources: Giacometti et al. (2002), Journal of Mammalogy; Büntgen et al. (2014), Ecology Letters

The 1 in a formula

The 1 in a formula stands for the intercept, and R adds it for you.

coef(lm(flipper_length_mm ~ 1 + body_mass_g, data = pengs))
 (Intercept)  body_mass_g 
136.72955927   0.01527592 
m0 <- lm(flipper_length_mm ~ 1, data = pengs)
coef(m0)
(Intercept) 
   200.9152

So y ~ x and y ~ 1 + x fit the same model. y ~ 1 fits an intercept-only model: its one estimate is the mean flipper length, the model from last week.

The mean as a baseline

The left panel shows each penguin’s distance from the mean, the residuals of y ~ 1. With our line (right, same 30 penguins) the residuals are much shorter.

R² step by step

R² = 1 − Σe²line / Σe²mean 
sq_mean <- sum(residuals(m0)^2)
sq_line <- sum(residuals(m)^2)
1 - sq_line / sq_mean
[1] 0.7589925

Σe² is a sum of squared residuals, as in least squares. R² is the proportion of the squared misses of the mean that the line removes, here 76%.

Reading R²

glance(m)$r.squared
[1] 0.7589925

This is the Multiple R-squared in summary(): 76% of the variation in flipper length goes with body mass.

It runs from 0, when the line does no better than the mean, to 1, when every point lies on the line. In linguistics much lower values are normal, because speakers and listeners vary for many reasons besides any one predictor.

Discussion

Body mass accounts for 76% of the variation in flipper length. Suppose a study finds that word frequency accounts for 5% of the variation in reaction times.

  • Is that model useless?
  • What might account for the other 95%?

Checking the model

Did you know: by the 19th century, Alpine ibex had been hunted out of the Alps except in the Italian king’s hunting reserve, today Gran Paradiso National Park. In the early 1900s, young ibex were smuggled from there into Switzerland, where the first reintroduction followed in 1911.

Sources: Gran Paradiso National Park; Springer Nature Research Communities

Are the residuals normal?

Linear models assume that the residuals are roughly normal, symmetric and centred on zero, and these residuals are close enough.

Is the spread constant?

The residuals should also have about the same spread all along the line. On the right, the misses grow with the fitted values. Chapter 4 of Winter (2019) covers both checks.

Grouped data

Did you know: outside the mating season, adult male and female Alpine ibex live in separate groups, and males join groups of males of a similar age.

Source: Villaret & Bon (1995), Ethology

A warning: three species

Bill length by bill depth, all penguins together:

The slope is negative: penguins with deeper bills seem to have shorter ones.

A warning: three species

Within every species the slope is positive. Leaving species out reversed the relationship: this is Simpson’s paradox. Lectures 6 and 7 show how to add species to the model.

Your turn

Exercises

With penguins, after drop_na() on the columns you use:

  1. Fit a model of bill_length_mm as a function of flipper_length_mm
  2. Report the slope in a sentence. What does the intercept mean here?
  3. What is R²? Compute it by hand with an intercept-only model
  4. Predict the bill length of penguins with 190 and 220 mm flippers
  5. Bonus: plot the residuals against the fitted values with augment(), coloured by species. What do you see?

Solutions

bill_data <- drop_na(penguins, bill_length_mm, flipper_length_mm)

m_bill <- lm(bill_length_mm ~ flipper_length_mm, data = bill_data)
coef(m_bill)
      (Intercept) flipper_length_mm 
       -7.2648678         0.2547682 
m_bill0 <- lm(bill_length_mm ~ 1, data = bill_data)
1 - sum(residuals(m_bill)^2) / sum(residuals(m_bill0)^2)
[1] 0.430574
glance(m_bill)$r.squared
[1] 0.430574

Solutions

predict(m_bill, newdata = tibble(flipper_length_mm = c(190, 220)))
       1        2 
41.14108 48.78413 
augment(m_bill, data = bill_data) |>
  ggplot(aes(x = .fitted, y = .resid, colour = species)) +
  geom_point()

Summary

  • y = b₀ + b₁ · x + e predicts the conditional mean of y from x
  • The intercept is the prediction at x = 0; the slope is the change in y when x goes up by 1
  • Least squares picks the line with the smallest sum of squared residuals, and lm(y ~ x) finds it
  • R² is the proportion of variation the line accounts for, compared with the intercept-only model y ~ 1
  • Check that the residuals are roughly normal, with the same spread along the line
  • Predict only within the range of the data, and look for groups that a single line hides

Slides and code

QR code linking to these slides

These slides

mcanzi.de/ws26_statistics/04_linear_models.html

Code, data and all lectures

github.com/mcanzi/ws26_statistics