Multiple regression

More than one continuous predictor

Dr Massimiliano Canzi

More than one predictor

Did you know: the word pangolin comes from the Malay pengguling: “one who rolls up”. Threatened pangolins curl into a ball of scales.

Sources: Born Free Foundation; National Geographic

Today

  • Why one predictor is often not enough
  • Fitting and reading a model with two continuous predictors
  • What “holding the other predictors constant” means
  • Comparing predictors with standardised coefficients
  • Model fit: R² and adjusted R²
  • Collinearity: when predictors overlap

So far

In lectures 4 and 5, every model had one predictor (y ~ x). We also saw a warning: with the penguins, leaving out species reversed a slope.

Outcomes in linguistics rarely depend on one thing. Reaction times depend on word frequency, word length, the speaker, the task…

Today: continuous predictors only. Categorical predictors (like species) come next week.

The data

A simulated lexical decision experiment: one row per word.

set.seed(6)
words <- tibble(
  log_freq  = runif(400, 1, 6),                       # log10 frequency
  length    = pmin(pmax(round(11 - 1.1 * log_freq +   # in letters:
                              rnorm(400, 0, 1.6)), 2), 14),  # frequent words are shorter
  syllables = pmax(round(length / 3 + rnorm(400, 0, 0.35)), 1),
  rt        = 760 - 30 * log_freq + 14 * length + rnorm(400, 0, 45)
)
words
# A tibble: 400 × 4
  log_freq length syllables    rt
     <dbl>  <dbl>     <dbl> <dbl>
1     4.03      9         3  778.
2     5.69      4         1  598.
3     2.32      8         2  754.
4     2.90      7         2  766.
# ℹ 396 more rows

We know the truth: frequent words are faster (−30 ms per log unit), longer words are slower (+14 ms per letter). Syllables have no effect of their own.

Two predictors, one outcome

Frequent words are faster and long words are slower.

The predictors are correlated

cor(words$log_freq, words$length)
[1] -0.7009648

Frequent words tend to be shorter. Each predictor carries part of the other’s information.

Two simple regressions

coef(lm(rt ~ log_freq, data = words))
(Intercept)    log_freq 
  909.92085   -44.73495 
coef(lm(rt ~ length,   data = words))
(Intercept)      length 
  549.30028    28.60034

Both slopes are too big compared with the truth (−30 and +14).

On its own, length also picks up the effect of frequency: long words are slow partly because they are rare.

One model, two predictors

Did you know: pangolin scales are made of keratin, like our fingernails. Scales and skin can make up about 25% of a pangolin’s body weight.

Sources: National Geographic; Born Free Foundation

The equation

y = b₀ + b₁ · x₁ + b₂ · x₂ + e
  • The intercept b₀ is the predicted y when both predictors are 0
  • Each slope is the change in y when its predictor goes up by 1 while the other stays the same

In R, add predictors with +:

lm(outcome ~ predictor1 + predictor2, data = ...)

lm() with two predictors

m <- lm(rt ~ log_freq + length, data = words)
coef(m)
(Intercept)    log_freq      length 
  749.78281   -29.13551    14.91379
  • Each extra log unit of frequency: -29.1 ms, for words of the same length
  • Each extra letter: 14.9 ms, for words of the same frequency

Both are now close to the true values. The multiple regression separated what the two predictors share.

Simple vs multiple

predictor true effect simple regression multiple regression
log_freq -30 -44.7 -29.1
length 14 28.6 14.9

With correlated predictors, a simple regression answers “what goes with x?”. A multiple regression answers “what goes with x, over and above the other predictors?”.

Holding the other predictor constant

Predicted RTs from the model, for words of three fixed frequencies. Each line has the same slope: 14.9 ms per letter.

Parallel lines

The lines are parallel because our model says the effect of length is the same at every frequency.

If the effect of length were different for frequent and rare words, the lines would not be parallel. That is an interaction: coming up in two weeks.

A regression plane

With two predictors, the model is a plane, seen here from above. The contour lines are straight and evenly spaced: no curves, no interactions.

summary()

summary(m)

Call:
lm(formula = rt ~ log_freq + length, data = words)

Residuals:
     Min       1Q   Median       3Q      Max 
-135.162  -28.989    0.268   29.255  141.362 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  749.783     16.551   45.30   <2e-16 ***
log_freq     -29.136      2.146  -13.58   <2e-16 ***
length        14.914      1.438   10.37   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 44.81 on 397 degrees of freedom
Multiple R-squared:  0.7079,    Adjusted R-squared:  0.7064 
F-statistic: 481.1 on 2 and 397 DF,  p-value: < 2.2e-16

{broom} and predictions

tidy(m)[, 1:3]
# A tibble: 3 × 3
  term        estimate std.error
  <chr>          <dbl>     <dbl>
1 (Intercept)    750.      16.6 
2 log_freq       -29.1      2.15
3 length          14.9      1.44

Predict for words we choose, giving every predictor a value:

new_words <- tibble(log_freq = c(2, 5), length = c(9, 4))
predict(m, newdata = new_words)
       1        2 
825.7359 663.7604

A rare, long word versus a frequent, short one: a difference of about 162 ms.

Standardised coefficients

Did you know: a pangolin’s tongue can be up to 70 cm long. Pangolins have no teeth, and can eat up to 70 million insects a year, mostly ants and termites.

Source: Born Free Foundation

Which predictor matters more?

coef(m)[-1]
 log_freq    length 
-29.13551  14.91379

-29 ms per log unit of frequency, or 15 ms per letter? The units differ, so the numbers cannot be compared directly.

From lecture 5: standardise the predictors, and every slope is “per standard deviation”.

Standardising the predictors

words <- words |>
  mutate(across(c(log_freq, length), \(x) as.numeric(scale(x)), .names = "{.col}_z"))

m_z <- lm(rt ~ log_freq_z + length_z, data = words)
coef(m_z)
(Intercept)  log_freq_z    length_z 
  748.35868   -42.71253    32.62463
  • One SD more frequency: -43 ms
  • One SD more length: 33 ms

The intercept is now the predicted RT for a word of average frequency and average length.

Same model, new numbers

c(raw = glance(m)$r.squared, standardised = glance(m_z)$r.squared)
         raw standardised 
   0.7079195    0.7079195

As in lecture 5: standardising changes what the coefficients mean, and the fit stays the same.

Standardised coefficients are useful to compare predictors within a model. Raw coefficients are easier to interpret (“ms per letter”): report the ones your reader needs.

Model fit

Did you know: Sunda pangolins climb trees to reach ant nests, holding on to branches with their prehensile tail. During the day they sleep in tree hollows or burrows.

Source: National Geographic

R² with more predictors

glance(lm(rt ~ log_freq,          data = words))$r.squared
[1] 0.6287711
glance(lm(rt ~ length,            data = words))$r.squared
[1] 0.5722565
glance(lm(rt ~ log_freq + length, data = words))$r.squared
[1] 0.7079195

Together, the predictors explain more than either alone, but less than the sum of the two: they overlap.

R² never goes down

Add a predictor of pure noise:

set.seed(1)
words <- words |> mutate(noise = rnorm(n()))

r2 <- c(glance(lm(rt ~ log_freq + length,         data = words))$r.squared,
        glance(lm(rt ~ log_freq + length + noise, data = words))$r.squared)
print(r2, digits = 6)
[1] 0.707919 0.707920

R² went up, a tiny bit. Every new predictor fits some of the noise in the sample.

Adjusted R²

Adjusted R² penalises each extra predictor (as.data.frame() prints more digits):

glance(lm(rt ~ log_freq + length,         data = words))[, 1:2] |> as.data.frame()
  r.squared adj.r.squared
1 0.7079195      0.706448
glance(lm(rt ~ log_freq + length + noise, data = words))[, 1:2] |> as.data.frame()
  r.squared adj.r.squared
1 0.7079199     0.7057072

Use adjusted R² to compare models with different numbers of predictors. A predictor that does not help lowers it.

Checking assumptions

The same checks as lecture 4, on one set of residuals for the whole model: roughly normal, constant spread.

Collinearity

Did you know: the Sunda pangolin is critically endangered. Pangolins are among the most trafficked mammals in the world; international commercial trade in them was banned in 2016.

Source: National Geographic

Overlapping predictors

Length in letters and length in syllables measure almost the same thing:

cor(words$length, words$syllables)
[1] 0.8443081

When two predictors are highly correlated, the model struggles to tell which one does the work. This is collinearity.

What collinearity does

tidy(lm(rt ~ log_freq + length,             data = words))[, 1:3]
# A tibble: 3 × 3
  term        estimate std.error
  <chr>          <dbl>     <dbl>
1 (Intercept)    750.      16.6 
2 log_freq       -29.1      2.15
3 length          14.9      1.44
tidy(lm(rt ~ log_freq + length + syllables, data = words))[, 1:3]
# A tibble: 4 × 3
  term        estimate std.error
  <chr>          <dbl>     <dbl>
1 (Intercept)   750.       16.6 
2 log_freq      -29.2       2.15
3 length         15.7       2.14
4 syllables      -2.40      4.93

Adding syllables makes the estimate for length less precise: its standard error grows. Syllables also gets a coefficient of its own, although in our simulation it has no effect, while the overall fit barely changes.

Collinearity, visually

The dotted line is the true effect of length (14 ms per letter). With syllables in the model, both intervals are much wider: the model cannot tell the two overlapping predictors apart.

Measuring collinearity: VIF

The variance inflation factor: how well can a predictor be predicted from the other predictors?

VIF = 1 / (1 − R²ⱼ) 
r2_len <- glance(lm(length ~ log_freq + syllables, data = words))$r.squared
1 / (1 - r2_len)          # length, with syllables in the model
[1] 4.346639
r2_len2 <- glance(lm(length ~ log_freq, data = words))$r.squared
1 / (1 - r2_len2)         # length, without syllables
[1] 1.965995

1 means no overlap. The higher the VIF, the more the standard error is inflated. Common rules of thumb flag values above 5 or 10.

What to do about it

  • Collinearity comes from your predictors themselves: R is working correctly
  • Think about theory: are the two predictors really different things?
  • Keep the predictor that best answers your research question, or combine them
  • Never add predictors just because the data happen to contain them

Discussion

Frequent words tend to be short, so word frequency and word length overlap. You want to know whether frequency affects reading times.

  • Would you put length in the model too? Why, or why not?
  • What does “the effect of frequency, holding length constant” mean for a real word?

Your turn

Exercises

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

  1. Predict body_mass_g from flipper_length_mm alone
  2. Add bill_length_mm and bill_depth_mm as predictors
  3. How did the flipper slope change? What is the adjusted R² of each model?
  4. Standardise the three predictors. Which one goes with the largest change in body mass?
  5. Bonus: compute the VIF of flipper_length_mm in the model from 2

Solutions

pengs <- drop_na(penguins, body_mass_g, flipper_length_mm, bill_length_mm, bill_depth_mm)

m1 <- lm(body_mass_g ~ flipper_length_mm, data = pengs)
m2 <- lm(body_mass_g ~ flipper_length_mm + bill_length_mm + bill_depth_mm, data = pengs)

coef(m1)[2]; coef(m2)[2]
flipper_length_mm 
         49.68557 
flipper_length_mm 
         50.26922 
c(m1 = glance(m1)$adj.r.squared, m2 = glance(m2)$adj.r.squared)
       m1        m2 
0.7582837 0.7593534

Solutions

pengs_z <- pengs |>
  mutate(across(c(flipper_length_mm, bill_length_mm, bill_depth_mm),
                \(x) as.numeric(scale(x))))

coef(lm(body_mass_g ~ flipper_length_mm + bill_length_mm + bill_depth_mm,
        data = pengs_z))[-1]
flipper_length_mm    bill_length_mm     bill_depth_mm 
        706.87140          22.72181          39.59368 
r2_flip <- glance(lm(flipper_length_mm ~ bill_length_mm + bill_depth_mm, data = pengs))$r.squared
1 / (1 - r2_flip)
[1] 2.673338

Summary

  • y ~ x1 + x2: each slope is the effect of one predictor, holding the others constant
  • With correlated predictors, simple and multiple slopes differ
  • No interactions: the effect of one predictor is the same at every value of the other (parallel lines)
  • Standardised coefficients compare predictors measured in different units
  • R² always rises with more predictors; adjusted R² does not
  • Collinearity inflates standard errors: check correlations and VIFs

Slides and code

QR code linking to these slides

These slides

mcanzi.de/ws26_statistics/06_multiple_regression.html

Code, data and all lectures

github.com/mcanzi/ws26_statistics