Mixed models 2

Fitting, testing and troubleshooting

Dr Massimiliano Canzi

Mixed models in practice

Did you know: octopus has three English plurals: octopuses, octopi (the oldest, from the idea that a Latin-looking word needs a Latin ending) and the rare octopodes (from the idea that it is Greek).

Source: Merriam-Webster

Today

Simulating a data set where we know the truth
Fitting and reading an lmer() model
Extracting estimates, random effects, confidence intervals and predictions
Significance testing: likelihood ratio tests
Common mistakes, made on purpose
R², convergence problems and singular fits

Last week

Repeated measures are not independent: data from the same speaker or word are alike
Mixed models add varying intercepts and varying slopes for each grouping factor
(1 + x | speaker) for predictors that vary within speakers; (1 | word) when x does not vary within words
Without varying slopes, standard errors are too small

Today: doing it properly, from simulation to write-up.

Simulating vowel durations

Why simulate?

With real data, we never know the true effects. With simulated data, we do, so we can check whether the model recovers them.

Like Winter (2019, section 15.2), we simulate vowel durations. Our design:

16 speakers each read the same 24 words. Half the words end in a voiced consonant, half in a voiceless one
True voicing effect: vowels 20 ms longer before voiced consonants

The simulation

set.seed(1)
n_speakers <- 16
n_words    <- 24

speakers <- tibble(speaker = sprintf("S%02d", 1:n_speakers),
                   s0 = rnorm(n_speakers, 0, 25),     # speaker's own speed
                   s1 = rnorm(n_speakers, 0, 12))     # speaker's own voicing effect
words <- tibble(word    = sprintf("w%02d", 1:n_words),
                voicing = rep(c("voiceless", "voiced"), each = n_words / 2),
                w0      = rnorm(n_words, 0, 15))     # word's own duration

vowels <- crossing(speakers, words) |>
  mutate(duration = round(150 + s0 + w0 +
                          (20 + s1) * (voicing == "voiced") +
                          rnorm(n(), 0, 20)),
         voicing = factor(voicing, levels = c("voiceless", "voiced"))) |>
  select(speaker, word, voicing, duration)

(voicing == "voiced") is TRUE or FALSE, which R treats as 1 or 0: a dummy variable, as in lecture 7.

The data

One line per speaker (blue), and the overall means (orange). Speakers differ in speed and in the size of their voicing effect.

Fitting the model

Did you know: octopuses have three hearts and blue blood: their blood carries oxygen with copper-based hemocyanin, not iron-based haemoglobin.

Source: Smithsonian Magazine

Which random effects?

Every speaker produces both voiced and voiceless words, so voicing varies within speakers and can get a slope: (1 + voicing | speaker). Each word is either voiced or voiceless, so words get only an intercept: (1 | word).

m <- lmer(duration ~ voicing + (1 + voicing | speaker) + (1 | word),
          data = vowels)

The maximal random effects structure justified by the design (Barr et al., 2013, lecture 10).

`summary()`

summary(m)

Linear mixed model fit by REML ['lmerMod']
Formula: duration ~ voicing + (1 + voicing | speaker) + (1 | word)
   Data: vowels

REML criterion at convergence: 3492.4

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.68928 -0.58558 -0.03095  0.65537  2.75815 

Random effects:
 Groups   Name          Variance Std.Dev. Corr  
 word     (Intercept)   158.3    12.583         
 speaker  (Intercept)   598.5    24.464         
          voicingvoiced  57.9     7.609   -0.06 
 Residual               402.0    20.049         
Number of obs: 384, groups:  word, 24; speaker, 16

Fixed effects:
              Estimate Std. Error t value
(Intercept)    155.516      7.259  21.423
voicingvoiced   20.516      5.848   3.508

Correlation of Fixed Effects:
            (Intr)
voicingvocd -0.376

Reading the output

Part of the output	What it tells you
`REML criterion`	a fit statistic for REML estimation (more below)
`Random effects`	how much speakers and words vary: variances, SDs and the correlation
`Fixed effects`	estimates, SEs and t values, read exactly as in `lm()`
`Correlation of Fixed Effects`	how the estimates covary: rarely of interest

The voicing effect is estimated at 20.5 ms (truth: 20). The intercept, 155.5 ms, is the mean duration before voiceless consonants (truth: 150).

Where are the p-values?

lme4 deliberately does not print p-values for the fixed effects.

For a t-test, we need the degrees of freedom. In a mixed model, it is not clear what they are: 384 observations? 16 speakers? 24 words? Something in between?

Several solutions exist. We will look at two: likelihood ratio tests, and approximate degrees of freedom with {lmerTest}.

Extracting information

Did you know: about two-thirds of an octopus’s neurons are in its arms, not its head, so each arm can act semi-independently.

Source: Smithsonian Magazine

Fixed effects and their SEs

fixef(m)

  (Intercept) voicingvoiced 
    155.51562      20.51563

coef(summary(m))

               Estimate Std. Error   t value
(Intercept)   155.51562   7.259179 21.423307
voicingvoiced  20.51563   5.847762  3.508287

Random effects

VarCorr(m)

 Groups   Name          Std.Dev. Corr   
 word     (Intercept)   12.5835         
 speaker  (Intercept)   24.4644         
          voicingvoiced  7.6089  -0.059 
 Residual               20.0491

sigma(m)              # the residual SD

[1] 20.04914

VarCorr() gives the SDs of the random effects. Speakers’ voicing effects vary with an SD of about 7.6 ms around the average of 20.5 ms.

Each speaker’s estimates

ranef(): each speaker’s adjustment to the fixed effects. coef(): fixed effects plus adjustments, that is each speaker’s own estimates.

head(ranef(m)$speaker, 3)

    (Intercept) voicingvoiced
S01  -17.635720    -0.4467663
S02    5.397587     4.2260787
S03  -16.590655     2.7863585

head(coef(m)$speaker, 3)

    (Intercept) voicingvoiced
S01    137.8799      20.06886
S02    160.9132      24.74170
S03    138.9250      23.30198

Speaker S01’s voicing effect: 20.52 + -0.45 = 20.07 ms.

Confidence intervals

confint(m, parm = "beta_", method = "profile")

                   2.5 %    97.5 %
(Intercept)   141.114192 169.91706
voicingvoiced   8.894863  32.13639

parm = "beta_" asks for the fixed effects only. Profile intervals can be slow for large models; method = "Wald" is a quick approximation:

confint(m, parm = "beta_", method = "Wald")

                   2.5 %    97.5 %
(Intercept)   141.287895 169.74336
voicingvoiced   9.054222  31.97703

Predictions

re.form = NA: the population-level prediction, for an average speaker and an average word:

new <- tibble(voicing = factor(c("voiceless", "voiced"), levels = levels(vowels$voicing)))
predict(m, newdata = new, re.form = NA)

       1        2 
155.5156 176.0312

Without re.form = NA, the prediction includes the random effects: here, speaker S01 reading word w01:

predict(m, newdata = tibble(speaker = "S01", word = "w01", voicing = new$voicing[1]))

       1 
137.7585

For a new speaker, not in the data, only the population-level prediction makes sense.

Significance testing

Did you know: octopus ink does more than hide: it contains tyrosinase, which irritates a predator’s eyes and garbles its senses of smell and taste.

Source: Smithsonian Magazine

Likelihood ratio tests

Compare two nested models: one with the predictor, one without it, and the same random effects.

m_null <- lmer(duration ~ 1 + (1 + voicing | speaker) + (1 | word),
               data = vowels, REML = FALSE)
m_full <- lmer(duration ~ voicing + (1 + voicing | speaker) + (1 | word),
               data = vowels, REML = FALSE)

Does adding voicing improve the fit by more than we would expect by chance?

`anova()`

anova(m_null, m_full)

Data: vowels
Models:
m_null: duration ~ 1 + (1 + voicing | speaker) + (1 | word)
m_full: duration ~ voicing + (1 + voicing | speaker) + (1 | word)
       npar    AIC    BIC  logLik -2*log(L)  Chisq Df Pr(>Chisq)   
m_null    6 3525.8 3549.5 -1756.9    3513.8                        
m_full    7 3517.3 3545.0 -1751.7    3503.3 10.506  1    0.00119 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

The full model fits much better: χ²(1) = 10.51, p = 0.0012.

What the test does

Each model has a log-likelihood, which says how probable the data are under that model’s estimates
The model with more parameters always fits at least as well, so the question is whether it fits better by more than chance would give
Twice the difference in log-likelihood is compared with a χ² distribution, whose degrees of freedom are the number of extra parameters

Here the models differ by one parameter, the voicing effect, so the test has 1 df.

REML or ML?

lmer() uses REML (restricted maximum likelihood) by default: better estimates of the variances
To compare models that differ in their fixed effects, the likelihoods must come from ML: REML = FALSE

If you forget, anova() refits the models with ML for you, and says so:

m_reml <- lmer(duration ~ voicing + (1 + voicing | speaker) + (1 | word), data = vowels)
m_reml_null <- update(m_reml, . ~ . - voicing)
anova(m_reml_null, m_reml)$"Pr(>Chisq)"

refitting model(s) with ML (instead of REML)

[1]          NA 0.001189648

p-values with `{lmerTest}`

{lmerTest} (Kuznetsova, Brockhoff & Christensen, 2017) adds approximate degrees of freedom and p-values to summary():

library(lmerTest)
m_lt <- lmer(duration ~ voicing + (1 + voicing | speaker) + (1 | word), data = vowels)
coef(summary(m_lt))

               Estimate Std. Error       df   t value     Pr(>|t|)
(Intercept)   155.51562   7.259179 24.21771 21.423307 2.982058e-17
voicingvoiced  20.51563   5.847762 25.08672  3.508287 1.722969e-03

The df are estimated with the Satterthwaite method: here about 25, much closer to the number of speakers and words than to 384 observations.

Reporting

Vowels were 20.5 ms longer before voiced than before voiceless consonants (95% CI [8.9, 32.1]). A likelihood ratio test against a model without voicing showed that this effect was reliable, χ²(1) = 10.51, p = 0.0012.

Also report the random effects structure: “The model included by-speaker random intercepts and random slopes for voicing, and by-word random intercepts.”

Common mistakes

A slope that cannot exist

What if we add a random slope for voicing by word?

m_bad <- lmer(duration ~ voicing + (1 + voicing | speaker) + (1 + voicing | word),
              data = vowels)

Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl =
control$checkConv, : unable to evaluate scaled gradient

Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model failed to converge: degenerate  Hessian with 1 negative eigenvalues
  See ?lme4::convergence and ?lme4::troubleshooting.

Each word is either voiced or voiceless, so a by-word voicing slope is impossible to estimate. The model cannot converge to a sensible answer: always check which slopes your design supports.

Forgetting the random slopes

m_int_only <- lmer(duration ~ voicing + (1 | speaker) + (1 | word), data = vowels)

c(with_slopes    = coef(summary(m))["voicingvoiced", "Std. Error"],
  without_slopes = coef(summary(m_int_only))["voicingvoiced", "Std. Error"])

   with_slopes without_slopes 
      5.847762       5.529815

Here the difference is small: speakers’ voicing effects vary little (SD ≈ 7.6 ms). The more slopes vary across speakers, the bigger the gap: in lecture 10, adding the slopes made the SE about three times larger.

Discussion

A paper reports a significant effect of a condition that every speaker took part in. The model has random intercepts for speakers, but no random slopes.

As a reviewer, what would you ask the authors, and why?
Would your trust in the result change if they added the slopes and the effect stayed significant?

Remaining issues

Did you know: in 2016, Inky, a common octopus at the National Aquarium of New Zealand in Napier, escaped: he squeezed through a gap in his tank lid, crossed the floor and slid down a drainpipe to the sea.

Source: RNZ

R² for mixed models

There is no single R². Nakagawa and Schielzeth (2013) distinguish:

Marginal R² is the variance explained by the fixed effects alone. Conditional R² also counts what the random effects explain.

For a model with random intercepts, both are simple ratios of variance components:

m_ri <- lmer(duration ~ voicing + (1 | speaker) + (1 | word), data = vowels)
var_fixed  <- var(predict(m_ri, re.form = NA))
var_random <- sum(as.data.frame(VarCorr(m_ri))$vcov[1:2])
var_resid  <- sigma(m_ri)^2

c(marginal    = var_fixed / (var_fixed + var_random + var_resid),
  conditional = (var_fixed + var_random) / (var_fixed + var_random + var_resid))

   marginal conditional 
 0.08233999  0.67447930

Reading R² for mixed models

Voicing alone accounts for only about 8% of the variation
With speakers and words, about 67%: most variation is between speakers and words

For models with random slopes, use a package, e.g. MuMIn::r.squaredGLMM().

Convergence warnings

Mixed models are fitted by an optimiser, which searches for the best parameter values. Sometimes it does not finish cleanly:

Model failed to converge with max|grad| = 0.0066 (tol = 0.002)

This is a warning rather than an error, so R still returns a model, but check it before you report anything.

What to do

Check the data and the design: is every random slope possible? (see above)
Centre or scale continuous predictors (lecture 5)
Try a different optimiser: lmer(..., control = lmerControl(optimizer = "bobyqa"))
Compare all optimisers with allFit(model): if they agree, the warning is usually harmless
Simplify the random effects: first remove the correlations with ||, then slopes with the least variance. || only works for numeric predictors: code a factor like voicing as 0/1 first

Barr et al. (2013) recommend starting maximal; Matuschek et al. (2017) show that maximal models can lose power, and recommend keeping only the random effects the data support.

Singular fits

boundary (singular) fit: see help('isSingular')

The model estimated a random-effect variance of zero, or a correlation of exactly ±1: the data cannot support that part of the random effects structure.

isSingular(m)

[1] FALSE

If a model is singular, consider simplifying the random effects (point 5 above), and report what you did.

Shrinkage and individual differences

Pink: each speaker’s raw voicing effect. Blue: the model’s estimate, shrunk towards the average. When interpreting individual speakers, remember that their estimates borrow strength from everyone else.

Your turn

Real data: `lexdec`

The lexdec data from {languageR} (Baayen, 2008): 21 participants made lexical decisions on 79 English words. RT is the log reaction time; 12 participants are native speakers of English, 9 are not.

library(languageR)
lexdec |> select(Subject, Word, RT, Frequency, NativeLanguage) |> head(3)

  Subject   Word       RT Frequency NativeLanguage
1      A1    owl 6.340359  4.859812        English
2      A1   mole 6.308098  4.605170        English
3      A1 cherry 6.349139  4.997212        English

Exercises

Which random slopes does the design allow for Frequency (a property of words) and NativeLanguage (a property of participants)?
Fit RT ~ Frequency + NativeLanguage with the maximal random effects structure justified by the design. What happens?
Refit with the bobyqa optimiser. Did the fixed effects change?
Test NativeLanguage with a likelihood ratio test.

Solutions

Frequency varies within participants, not within words: (1 + Frequency | Subject). NativeLanguage varies within words, not within participants: (1 + NativeLanguage | Word)

lex_max <- lmer(RT ~ Frequency + NativeLanguage +
                  (1 + Frequency | Subject) + (1 + NativeLanguage | Word),
                data = lexdec)

Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model failed to converge with max|grad| = 0.00767209 (tol = 0.002, component 1)
  See ?lme4::convergence and ?lme4::troubleshooting.

Solutions

lex_bob <- update(lex_max, control = lmerControl(optimizer = "bobyqa"))
cbind(default = fixef(lex_max), bobyqa = fixef(lex_bob))

                        default      bobyqa
(Intercept)          6.54085899  6.54087587
Frequency           -0.03987850 -0.03987869
NativeLanguageOther  0.07862839  0.07859110

The warning disappears, and the estimates barely change.

lex_full <- update(lex_bob, REML = FALSE)
lex_null <- update(lex_full, . ~ . - NativeLanguage)
anova(lex_null, lex_full)

Data: lexdec
Models:
lex_null: RT ~ Frequency + (1 + Frequency | Subject) + (1 + NativeLanguage | Word)
lex_full: RT ~ Frequency + NativeLanguage + (1 + Frequency | Subject) + (1 + NativeLanguage | Word)
         npar     AIC     BIC logLik -2*log(L)  Chisq Df Pr(>Chisq)
lex_null    9 -929.64 -880.92 473.82   -947.64                     
lex_full   10 -930.10 -875.96 475.05   -950.10 2.4608  1     0.1167

With random intercepts only, the same test gives p = 0.012. The by-word slope for NativeLanguage makes the test more cautious.

Summary

Simulate first: you can only check a model against a known truth
summary(), fixef(), ranef(), coef(), VarCorr(), confint(), predict(re.form = NA)
lme4 gives no p-values: use likelihood ratio tests (with ML), or {lmerTest}
Random slopes only where the predictor varies within the group
Convergence warnings and singular fits: check, change optimiser, then simplify
Marginal vs conditional R²

Slides and code

These slides

mcanzi.de/ws26_statistics/11_mixed_models_2.html

Code, data and all lectures

github.com/mcanzi/ws26_statistics

Mixed models 2

Fitting, testing and troubleshooting

Dr Massimiliano Canzi

Mixed models in practice

Today

Last week

Simulating vowel durations

Why simulate?

The simulation

The data

Fitting the model

Which random effects?

summary()

Reading the output

Where are the p-values?

Extracting information

Fixed effects and their SEs

Random effects

Each speaker’s estimates

Confidence intervals

Predictions

Significance testing

Likelihood ratio tests

anova()

What the test does

REML or ML?

p-values with {lmerTest}

Reporting

Common mistakes

A slope that cannot exist

Forgetting the random slopes

Discussion

Remaining issues

R² for mixed models

Reading R² for mixed models

Convergence warnings

What to do

Singular fits

Shrinkage and individual differences

Your turn

Real data: lexdec

Exercises

Solutions

Solutions

Summary

Slides and code

`summary()`

`anova()`

p-values with `{lmerTest}`

Real data: `lexdec`