Statistics refresher

Linear models, interactions and mixed models

Dr Massimiliano Canzi

From lm() to lmer()

Did you know: king penguins breed in colonies of up to 300,000 birds and build no nests. When a parent comes back from the sea, it and its chick find each other by call alone.

Sources: Aubin & Jouventin (1998), Proceedings of the Royal Society B; Howlett (1998), Nature

Why a refresher

This course builds on Statistics. Before we add anything new, we go once through the models from that course, in the order we met them: one predictor, several predictors, interactions, and mixed models.

Nothing here is new. If a slide feels unfamiliar, the matching lecture of the first course has the details, and so does Winter (2019).

Today

  1. The linear model
  2. Categorical predictors and several predictors
  3. Interactions
  4. Mixed models

The data

A simulated lexical decision task: 30 participants (15 L1 and 15 L2 speakers of English) each respond to the same 40 words.

set.seed(2027)
participants <- tibble(participant = sprintf("P%02d", 1:30),
                       group = rep(c("L1", "L2"), each = 15),
                       p0 = rnorm(30, 0, 60),
                       p1 = rnorm(30, 0, 12))
words <- tibble(word = sprintf("w%02d", 1:40),
                log_freq = runif(40, 1, 6),
                length = pmin(pmax(round(11 - 1.1 * log_freq + rnorm(40, 0, 1.6)), 3), 13),
                w0 = rnorm(40, 0, 30))
lex <- crossing(participants, words) |>
  mutate(slope = -25 + p1 + if_else(group == "L2", -25, 0),
         rt = round(700 + p0 + w0 + if_else(group == "L2", 120, 0) +
                    slope * (log_freq - 3.5) + 12 * (length - 7) + rnorm(n(), 0, 70))) |>
  select(participant, group, word, log_freq, length, rt)

The data

lex
# A tibble: 1,200 × 6
  participant group word  log_freq length    rt
  <chr>       <chr> <chr>    <dbl>  <dbl> <dbl>
1 P01         L1    w01       4.38      3   552
2 P01         L1    w02       3.84      8   628
3 P01         L1    w03       2.64      8   629
4 P01         L1    w04       4.85      5   659
# ℹ 1,196 more rows

1200 rows: one per participant and word. Because we simulated the data, we know the true effects: frequent words are faster, long words are slower, L2 speakers are slower overall, and frequency matters more for them.

A first look

Grey: every response. Salmon: the mean RT for each word. Responses are noisy, but frequent words are faster on average.

The linear model

Did you know: in playback experiments, king penguin chicks picked out their parent’s call even when it was 6 dB quieter than the calls of other adults played at the same time.

Source: Aubin & Jouventin (1998), Proceedings of the Royal Society B

The equation

y = b₀ + b₁ · x + e
  • The intercept b₀ is the predicted y when x = 0
  • The slope b₁ is how much y changes when x increases by one unit
  • The residual e is the distance between each observation and the line

lm() finds the intercept and slope that minimise the sum of squared residuals: least squares.

Intercept, slope and residuals

Word means, the regression line, and each word’s residual (vertical lines). The intercept sits at log frequency 0, far from the data.

One predictor

m1 <- lm(rt ~ log_freq, data = lex)
coef(m1)
(Intercept)    log_freq 
   926.2022    -45.9173

Each extra log unit of frequency goes with 46 ms faster responses. The intercept is the predicted RT for a word with log frequency 0, which none of our words has.

Centering

lex <- lex |>
  mutate(freq_c = log_freq - mean(log_freq),
         length_c = length - mean(length))
coef(lm(rt ~ freq_c, data = lex))
(Intercept)      freq_c 
   769.6467    -45.9173

Centering moves 0 to the mean. The slope stays the same; the intercept becomes the predicted RT for a word of average frequency.

Centering, visually

The line and its slope stay where they were. Centering moves 0, and with it the intercept (gold).

tidy() and glance()

tidy(lm(rt ~ freq_c, data = lex))
# A tibble: 2 × 5
  term        estimate std.error statistic  p.value
  <chr>          <dbl>     <dbl>     <dbl>    <dbl>
1 (Intercept)    770.       3.69     209.  0       
2 freq_c         -45.9      2.50     -18.4 1.41e-66
glance(lm(rt ~ freq_c, data = lex))[, c("r.squared", "sigma", "nobs")]
# A tibble: 1 × 3
  r.squared sigma  nobs
      <dbl> <dbl> <int>
1     0.220  128.  1200

R² is the share of the variance in RT that the model accounts for. sigma is the residual standard deviation, in ms.

Checking the residuals

Linear models assume residuals that are roughly normal and equally spread along the fitted values. These residuals look fine on both counts.

More predictors

Did you know: king penguin chicks need very little of a call to recognise a parent. The first half of one syllable (0.23 s), with only its three lowest harmonics, was enough.

Source: Jouventin, Aubin & Lengagne (1999), Animal Behaviour

A categorical predictor

coef(lm(rt ~ group, data = lex))
(Intercept)     groupL2 
   692.6400    154.0133

With treatment coding, R’s default, group becomes a 0/1 dummy. The intercept is the mean of the reference level (L1, first in the alphabet), and groupL2 is the difference between L2 and L1.

Treatment and sum coding, visually

Grey: each participant’s mean. Salmon: the group means.

Sum coding

lex <- lex |> mutate(group = factor(group))
contrasts(lex$group) <- contr.sum(2)
coef(lm(rt ~ group, data = lex))
(Intercept)      group1 
  769.64667   -77.00667

With sum coding (−1 and 1 here, R’s contr.sum() gives L1 = 1), the intercept is the mean of the two group means, and group1 is half the difference between them. This matters once we add interactions.

Two continuous predictors

m2 <- lm(rt ~ freq_c + length_c, data = lex)
coef(m2)
(Intercept)      freq_c    length_c 
  769.64667   -34.09109    10.15483

Each slope is now the effect of one predictor holding the other constant. The frequency slope shrinks from -46 to -34 ms: in the simple model, frequency also stood in for length, because frequent words are shorter.

Simple and multiple slopes

Frequent words are also short (blue), so the frequency-only line (dashed) mixes the two effects. Holding length constant gives the shallower slope.

Correlated predictors

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

Frequent words tend to be short.

Collinearity

  • When predictors overlap, each coefficient rests on the part of the predictor the others do not share, so standard errors grow. This is collinearity
  • Check it with a correlation matrix or variance inflation factors (car::vif()), and think about which predictors the design really needs

Interactions

Did you know: king penguin chicks still recognised their parent’s call when its frequency was shifted 100 Hz down or 75 Hz up. Their recognition does not depend on exact pitch.

Source: Jouventin, Aubin & Lengagne (1999), Animal Behaviour

Two groups, two slopes

Frequency matters more for L2 speakers: the lines are not parallel.

An interaction in R

m3 <- lm(rt ~ freq_c * group, data = lex)
coef(m3)
  (Intercept)        freq_c        group1 freq_c:group1 
    769.64667     -45.91730     -77.00667      17.33990

freq_c * group is short for freq_c + group + freq_c:group. The interaction term freq_c:group1 says how much the frequency slope changes between the groups.

Reading the coefficients

With sum-coded group and centred freq_c:

  • (Intercept) is the predicted RT for an average-frequency word, averaged over the two groups
  • freq_c is the frequency slope, averaged over the groups
  • group1 is half the L1–L2 difference at average frequency
  • freq_c:group1 is half the difference between the two slopes

With treatment coding and uncentred predictors, the same model is correct but group would be the L1–L2 difference at log frequency 0. Each lower-order term in an interaction model is a simple effect: the effect when the other predictor is 0.

Reading interactions

  • Plot the data and the fitted lines first
  • Check how each predictor is coded and where its 0 is
  • Read the interaction as a difference in slopes, or a difference in differences
  • For three-way interactions, plot again before reading any number

Mixed models

Did you know: king and emperor penguins call with two voices at once, one from each side of the syrinx. When one voice was removed from an emperor penguin’s call in playbacks, neither its mate nor its chick replied.

Source: Aubin, Jouventin & Hildebrand (2000), Proceedings of the Royal Society B

Not independent

So far, lm() treated our 1200 rows as 1200 independent observations. They are not:

  • each participant gave 40 responses, and some participants are faster than others
  • each word was seen by 30 participants, and some words are easier than others

Ignoring this makes standard errors too small, so effects look more certain than they are.

Participants differ

One line per participant: they differ in overall speed (intercepts) and in how much frequency helps them (slopes).

Fixed and random effects

  • Fixed effects are the predictors we want to estimate: frequency, length, group, their interaction
  • Random effects model how grouping units vary around them: participants and words
  • A random intercept lets each participant or word have its own baseline
  • A random slope lets the effect of a predictor vary by participant (or word)

A random slope needs a predictor that varies within the grouping unit. Frequency and length vary within participants; group varies within words, because every word is seen by both groups.

lmer() syntax

Formula term Meaning
(1 | participant) random intercepts by participant
(1 + freq_c | participant) random intercepts and frequency slopes, with their correlation
(1 + freq_c || participant) the same, without the correlation
(1 | word) random intercepts by word

We fit rt ~ freq_c * group + length_c + (1 + freq_c | participant) + (1 | word). The maximal structure for this design would also include by-participant slopes for length and by-word slopes for group (Barr et al., 2013); we leave them out here to keep the example small.

Fitting the model

m4 <- lmer(rt ~ freq_c * group + length_c + (1 + freq_c | participant) + (1 | word),
           data = lex)
round(fixef(m4), 1)
  (Intercept)        freq_c        group1      length_c freq_c:group1 
        769.6         -34.1         -77.0          10.2          17.3

The fixed effects tell the same story as before: frequency speeds up responses, length slows them down, L2 speakers are slower, and the frequency slope differs between groups.

Random effects

VarCorr(m4)
 Groups      Name        Std.Dev. Corr   
 word        (Intercept) 26.428          
 participant (Intercept) 63.453          
             freq_c      10.575   -0.078 
 Residual                69.808
  • Participants’ baselines vary with a standard deviation of about 63 ms, words’ baselines with about 26 ms
  • Participants’ frequency slopes vary around the fixed slope with a standard deviation of about 11 ms
  • Residual is the variation left over within participants and words

Random effects, visually

Each participant’s and each word’s random intercept, with a 95% interval. Participants spread more widely than words, as the standard deviations on the previous slide say.

Standard errors

rbind(lm = coef(summary(lm(rt ~ freq_c * group + length_c, data = lex)))["freq_c", 1:2],
      lmer = coef(summary(m4))["freq_c", 1:2])
      Estimate Std. Error
lm   -34.09109   3.004794
lmer -34.09109   5.311654

The estimate barely changes, but its standard error grows once the model knows that responses from the same participant and word belong together.

Wider intervals

Estimate ± 1.96 SE for each fixed effect. The intervals widen most for group1, which varies between participants: only 30 participants stand behind it, not 1,200 responses.

Partial pooling

Each participant’s slope moves toward the overall slope (dashed). This shrinkage, also called partial pooling, lets extreme estimates borrow strength from the rest of the sample.

Testing an effect

A likelihood ratio test compares two nested models fitted with maximum likelihood:

m4_ml <- update(m4, REML = FALSE)
m4_noint <- update(m4_ml, . ~ . - freq_c:group)
anova(m4_noint, m4_ml)
Data: lex
Models:
m4_noint: rt ~ freq_c + group + length_c + (1 + freq_c | participant) + (1 | word)
m4_ml: rt ~ freq_c * group + length_c + (1 + freq_c | participant) + (1 | word)
         npar   AIC   BIC  logLik -2*log(L)  Chisq Df Pr(>Chisq)    
m4_noint    9 13843 13889 -6912.4     13825                         
m4_ml      10 13813 13864 -6896.6     13793 31.647  1  1.849e-08 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

p-values with {lmerTest}

library(lmerTest)
m4 <- lmer(rt ~ freq_c * group + length_c + (1 + freq_c | participant) + (1 | word),
           data = lex)
round(coef(summary(m4)), 3)
              Estimate Std. Error     df t value Pr(>|t|)
(Intercept)    769.647     12.479 34.879  61.674    0.000
freq_c         -34.091      5.312 45.943  -6.418    0.000
group1         -77.007     11.759 28.001  -6.549    0.000
length_c        10.155      3.279 37.000   3.097    0.004
freq_c:group1   17.340      2.366 27.999   7.330    0.000

{lmerTest} adds degrees of freedom (Satterthwaite’s approximation) and p-values to the fixed effects.

When things go wrong

  • A convergence warning means the optimiser may not have found the best estimates. Check the data and scaling, try another optimiser (lmerControl(optimizer = "bobyqa")), and simplify only if needed
  • A singular fit means some variance or correlation is estimated at its boundary, often a random-effect correlation of ±1 or a standard deviation of 0. The data cannot support that part of the structure

In both cases, simplify step by step and report what you removed: first correlations, then the random slopes that explain least.

Reporting

A sentence for our model:

We fitted a linear mixed model to RTs with centred log frequency, sum-coded group, their interaction and centred length as fixed effects, by-participant random intercepts and frequency slopes, and by-word random intercepts (lme4; Bates et al., 2015). The frequency effect was stronger for L2 than for L1 speakers (interaction, with L1 coded 1 and L2 −1: b = 17.3 ms, SE = 2.4, p < .001).

Report the full model structure, the coding of each predictor, and how the model was simplified, if it was.

Summary

  • A linear model estimates an intercept and slopes; centering makes the intercept meaningful
  • Treatment and sum coding answer different questions about categorical predictors
  • With several predictors, each slope holds the others constant; collinearity makes slopes uncertain
  • In an interaction model, each lower-order term is a simple effect, so coding and centering decide what it means
  • Mixed models add random intercepts and slopes for participants and items, which gives honest standard errors and shrinks individual estimates toward the group

References

  • Barr, D. J., Levy, R., Scheepers, C., & Tily, H. J. (2013). Random effects structure for confirmatory hypothesis testing: Keep it maximal. Journal of Memory and Language, 68(3), 255–278.
  • Bates, D., Mächler, M., Bolker, B., & Walker, S. (2015). Fitting linear mixed-effects models using lme4. Journal of Statistical Software, 67(1), 1–48.
  • Kuznetsova, A., Brockhoff, P. B., & Christensen, R. H. B. (2017). lmerTest package: Tests in linear mixed effects models. Journal of Statistical Software, 82(13), 1–26.
  • Winter, B. (2019). Statistics for Linguists: An Introduction Using R. Routledge. Chapters 4–8, 14 and 15.

Sources of the animal facts

Slides and code

QR code linking to these slides

These slides

mcanzi.de/ss27_advanced_statistics/01_refresher.html

Code and all lectures

github.com/mcanzi/ss27_advanced_statistics