Generalised additive mixed models

Random smooths, autocorrelation and non-linear interactions

Dr Massimiliano Canzi

Generalised additive mixed models

Did you know: Egyptian fruit bats live in colonies of hundreds to thousands of individuals, can live for at least 25 years in captivity, and are extremely vocal in their social interactions.

Source: Prat, Taub & Yovel (2015), Science Advances

Where we are

Last time: smooth terms, difference smooths and significance testing with GAMs, in models that treated every measurement as independent. Today: the “mixed” part, following the second half of Sóskuthy’s (2018) workshop.

Today

  1. Dependencies in curve data
  2. Random effects in GAMMs
  3. Autocorrelation
  4. Non-linear interactions
  5. Checking and reporting

The data

The PRICE trajectories from lecture 3, now with an ordered factor for L1.

set.seed(2029)
speakers <- tibble(speaker = sprintf("S%02d", 1:20),
                   l1 = rep(c("English", "Italian"), each = 10),
                   s0 = rnorm(20, 0, 60),
                   s_mid = rnorm(20, 0, 5))
tokens <- crossing(speakers, token = 1:10) |>
  mutate(traj = paste(speaker, token, sep = "."),
         duration = round(runif(n(), 120, 300)),
         t0 = rnorm(n(), 0, 40))
price <- crossing(tokens, time = seq(0, 100, by = 10)) |>
  arrange(traj, time) |>
  mutate(onset = if_else(l1 == "English", 1250, 1320),
         offset = if_else(l1 == "English", 2050, 1900) - 2 * (250 - duration),
         mid = if_else(l1 == "English", 62, 50) + s_mid,
         width = if_else(l1 == "English", 8, 14),
         f2 = onset + (offset - onset) * plogis((time - mid) / width) + s0 + t0) |>
  group_by(traj) |>
  mutate(f2 = round(f2 + as.numeric(arima.sim(list(ar = 0.6), n = n(), sd = 30)))) |>
  ungroup() |>
  mutate(across(c(speaker, l1, traj), factor), l1_o = as.ordered(l1)) |>
  select(speaker, l1, l1_o, traj, duration, time, f2)
contrasts(price$l1_o) <- "contr.treatment"

Speakers and tokens

Ten tokens from each of four speakers. Each speaker’s tokens sit at their own height, and within a token, neighbouring measurements move up and down together.

Dependencies

Did you know: Prat and colleagues recorded 22 Egyptian fruit bats for 75 days. Their calls carried information about the caller, the context and even the addressee: the bat the call was directed at.

Source: Prat, Taub & Yovel (2016), Scientific Reports

Four kinds of model

Fixed effects only With random effects
Parametric terms LM, lm() LMER, lmer()
Smooth terms GAM, bam() GAMM, bam()

In {mgcv}, random effects are written as smooth terms, so the same function fits GAMs and GAMMs.

Two levels of dependency

  • Neighbouring measurements of the same token are more similar than measurements far apart: dependency within a trajectory
  • The 110 measurements of a speaker share that speaker’s vocal tract and habits: dependency within a speaker, and in real data within a word too

Lecture 3’s model ignored both, so its confidence intervals were too narrow.

Residual autocorrelation

resid_acf <- function(res, traj, max_lag = 8) {
  d <- tibble(res, traj)
  map(0:max_lag, \(l) d |>
    mutate(res_lag = lag(res, l), .by = traj) |>
    summarise(acf = cor(res, res_lag, use = "complete.obs")) |>
    mutate(lag = l)) |>
    list_rbind()
}
m_gam <- bam(f2 ~ l1_o + s(time, k = 10) + s(time, by = l1_o, k = 10),
             data = price, method = "ML")
acf_gam <- resid_acf(residuals(m_gam), price$traj)

The autocorrelation at lag 1 is the correlation between each residual and the residual one measurement earlier in the same trajectory; lag 2 skips one measurement, and so on. itsadug::acf_resid() does the same with base R graphics.

Residual autocorrelation

ggplot(acf_gam, aes(x = lag, y = acf)) +
  geom_hline(yintercept = 0, colour = muted) +
  geom_segment(aes(xend = lag, yend = 0), colour = kingfisher, linewidth = 1.2) +
  geom_point(colour = kingfisher, size = 3) +
  labs(x = "Lag (measurements)", y = "Autocorrelation")

Neighbouring residuals correlate at about 0.94, and the correlation is still clearly positive several measurements apart.

Random effects in GAMMs

Did you know: fruit bat pups raised without hearing adult bats kept producing underdeveloped, pup-like calls far into adolescence.

Source: Prat, Taub & Yovel (2015), Science Advances

Three kinds of random effect

{mgcv} term Random effect
s(speaker, bs = "re") random intercepts by speaker
s(speaker, time, bs = "re") random slopes over time by speaker
s(time, speaker, bs = "fs", m = 1) random smooths over time by speaker

Watch the order: random slopes put the grouping factor first, random smooths put the continuous predictor first.

Fitting the three

m_int <- bam(f2 ~ l1_o + s(time, k = 10) + s(time, by = l1_o, k = 10) +
               s(speaker, bs = "re"),
             data = price, method = "ML")
m_slope <- bam(f2 ~ l1_o + s(time, k = 10) + s(time, by = l1_o, k = 10) +
                 s(speaker, bs = "re") + s(speaker, time, bs = "re"),
               data = price, method = "ML")
m_smooth <- bam(f2 ~ l1_o + s(time, k = 10) + s(time, by = l1_o, k = 10) +
                  s(time, speaker, bs = "fs", m = 1, k = 5),
                data = price, method = "ML")

What each one allows

One line per speaker. Intercepts shift each speaker’s curve up or down, slopes tilt it, and random smooths let each speaker’s curve take its own shape.

Random smooths

  • bs = "fs" fits one smooth per speaker, with a shared smoothing parameter, much like random slopes share one variance
  • m = 1 changes the penalty so that each speaker’s smooth is pulled toward a flat line, as a random effect should be (Sóskuthy, 2017)
  • k applies to each speaker: 20 speakers with k = 5 add 100 coefficients. Keep k low; too many coefficients for the data stops the model with an error
  • Many levels make random smooths slow and memory-hungry

The random slope equivalent

  • Our groups differ between speakers: each speaker is English or Italian, so the difference smooth cannot vary by speaker
  • When a comparison is within speakers, say two vowels from every speaker, the model also needs speakers to differ in the vowel difference, like a by-speaker random slope in lmer()

One way is a random smooth for every speaker and vowel combination, s(time, speaker_vowel, bs = "fs", m = 1). Sóskuthy (2021) compares this and other options in simulations.

Autocorrelation

Did you know: isolated fruit bat pups that heard recordings of adult bats developed a repertoire that copied the playbacks.

Source: Prat, Taub & Yovel (2015), Science Advances

Two solutions

Random smooths by speaker handle dependencies within speakers. For the dependencies within each trajectory, there are two options:

  • Random smooths by trajectory, s(time, traj, bs = "fs", m = 1), give each of the 200 tokens its own smooth
  • An AR1 error model predicts the error at each measurement partly from the error at the previous one, with the correlation rho

An AR1 model

m_noar <- bam(f2 ~ l1_o + s(time, k = 10) + s(time, by = l1_o, k = 10) +
                s(time, speaker, bs = "fs", m = 1, k = 5),
              data = price, method = "ML")
rho <- itsadug::start_value_rho(m_noar)
m_ar <- bam(f2 ~ l1_o + s(time, k = 10) + s(time, by = l1_o, k = 10) +
              s(time, speaker, bs = "fs", m = 1, k = 5),
            data = price, method = "ML", AR.start = price$time == 0, rho = rho)
  • start_value_rho() estimates rho from the lag-1 autocorrelation of a model without the AR1 component
  • AR.start marks the first measurement of each trajectory, so the error model does not link one token to the next. The data must be sorted by trajectory and time

Did it work?

For the AR1 model, the autocorrelation is computed on the standardised residuals, m_ar$std.rsd, which take the error model into account.

Random smooths or AR1?

  • Random smooths by trajectory become very slow with many trajectories: Sóskuthy (2018) puts the practical limit at around 500 to 800
  • AR1 adds almost no computing time, so for large data sets it is the only option
  • An AR1 model may not remove all autocorrelation for some kinds of curves: always check the residuals
  • rho is shared by all trajectories, and in {mgcv} the AR1 option is available in bam() only

Non-linear interactions

Did you know: fruit bat pups raised in three groups, each hearing a different playback, developed distinct dialects from about six months of age, each closer to its own playback.

Source: Prat, Azoulay, Dor & Yovel (2017), PLoS Biology

Duration and shape

ggplot(price, aes(x = time, y = f2, group = traj, colour = duration)) +
  geom_line(alpha = 0.6) +
  scale_colour_viridis_c(option = "magma", begin = 0.2) +
  facet_wrap(~ l1) +
  labs(x = "Time (%)", y = "F2 (Hz)", colour = "Duration (ms)") +
  theme(legend.position = "right")

Short vowels (dark) stop lower: the glide undershoots its target. Duration changes the shape of the curve as well as its height.

s(), te() and ti()

  • s(duration) adds a smooth effect of duration on overall F2: it can shift whole curves, but not change their shape
  • te(time, duration) is a tensor product smooth: one surface over time and duration, including both main effects
  • ti(time, duration) is the tensor product interaction alone. Combined with s(time) and s(duration), it separates the main effects from the interaction, so each gets its own row in the summary

The model

m_dur <- bam(f2 ~ l1_o + s(time, k = 10) + s(time, by = l1_o, k = 10) +
               s(duration) + ti(time, duration) +
               s(time, speaker, bs = "fs", m = 1, k = 5),
             data = price, method = "ML", AR.start = price$time == 0, rho = rho)
round(summary(m_dur)$s.table, 2)
                      edf Ref.df      F p-value
s(time)              8.88   8.97 190.31       0
s(time):l1_oItalian  8.67   8.91  36.40       0
s(duration)          1.00   1.00 304.45       0
ti(time,duration)    4.51   5.17  83.45       0
s(time,speaker)     78.30  98.00  10.41       0

Predictions at three durations

predict_gam(m_dur, series = "time", length_out = 50,
            values = list(duration = c(150, 220, 290)),
            exclude_terms = "s(time,speaker)") |>
  ggplot(aes(x = time, y = f2, colour = factor(duration), fill = factor(duration))) +
  geom_ribbon(aes(ymin = lower_ci, ymax = upper_ci), alpha = 0.2, colour = NA) +
  geom_line(linewidth = 1.2) +
  facet_wrap(~ l1_o) +
  labs(x = "Time (%)", y = "Predicted F2 (Hz)", colour = "Duration (ms)", fill = "Duration (ms)") +
  theme(legend.position = "right")

exclude_terms removes the random smooths, so the curves are for an average speaker.

The surface

predict_gam(m_dur, series = c("time", "duration"), length_out = 40,
            values = list(l1_o = "English"), exclude_terms = "s(time,speaker)") |>
  ggplot(aes(x = time, y = duration, fill = f2, z = f2)) +
  geom_raster() +
  geom_contour(colour = floor_col, alpha = 0.6, bins = 8) +
  scale_fill_viridis_c(option = "mako", begin = 0.15) +
  labs(x = "Time (%)", y = "Duration (ms)", fill = "Predicted\nF2 (Hz)") +
  theme(legend.position = "right")

A heat map shows the whole interaction for English speakers, but prediction curves at a few durations are easier to read.

Checking and reporting

Did you know: the pups in the dialect experiment heard their mothers too, but their calls moved toward the playback of the whole colony. Prat and colleagues called this crowd vocal learning.

Source: Prat, Azoulay, Dor & Yovel (2017), PLoS Biology

Checking the model

k.check(m_dur)
                     k'       edf   k-index p-value
s(time)               9  8.880414 1.0565009  0.9975
s(time):l1_oItalian   9  8.669315 1.0565009  0.9950
s(duration)           9  1.001933 0.2914938  0.0000
ti(time,duration)    16  4.514466 1.0488305  0.9975
s(time,speaker)     100 78.296626 1.0565009  0.9925
p_resid <- ggplot(tibble(fitted = fitted(m_dur), res = m_dur$std.rsd), aes(fitted, res)) +
  geom_point(colour = kingfisher, alpha = 0.2) + labs(x = "Fitted F2 (Hz)", y = "Std. residual")
p_qq <- ggplot(tibble(res = m_dur$std.rsd), aes(sample = res)) +
  geom_qq(colour = kingfisher, alpha = 0.4) + geom_qq_line(colour = mist) + labs(x = "Theoretical", y = "Sample")
p_resid + p_qq

Testing the group difference

m_dur_null <- bam(f2 ~ s(time, k = 10) + s(duration) + ti(time, duration) +
                    s(time, speaker, bs = "fs", m = 1, k = 5),
                  data = price, method = "ML", AR.start = price$time == 0, rho = rho)
itsadug::compareML(m_dur_null, m_dur)
m_dur_null: f2 ~ s(time, k = 10) + s(duration) + ti(time, duration) + s(time, 
    speaker, bs = "fs", m = 1, k = 5)

m_dur: f2 ~ l1_o + s(time, k = 10) + s(time, by = l1_o, k = 10) + s(duration) + 
    ti(time, duration) + s(time, speaker, bs = "fs", m = 1, k = 5)

Chi-square test of ML scores
-----
       Model    Score Edf Difference    Df  p.value Sig.
1 m_dur_null 11087.91  10                               
2      m_dur 11035.71  13     52.200 3.000  < 2e-16  ***

AIC difference: 18.90, model m_dur has lower AIC.

Both models have the same random effects, error model and rho, and are fitted with method = "ML".

GAM versus GAMM

compare <- list(l1_o = c("Italian", "English"))
bind_rows(`GAM (lecture 3)` = as_tibble(get_difference(m_gam, "time", compare, length_out = 50)),
          `GAMM with AR1` = as_tibble(get_difference(m_dur, "time", compare, length_out = 50,
                                                     exclude_terms = "s(time,speaker)")),
          .id = "model") |>
  ggplot(aes(x = time, y = diff, colour = model, fill = model)) +
  geom_hline(yintercept = 0, colour = muted, linetype = "dashed") +
  geom_ribbon(aes(ymin = lower_ci, ymax = upper_ci), alpha = 0.25, colour = NA) +
  geom_line(linewidth = 1.1) +
  labs(x = "Time (%)", y = "Italian − English F2 (Hz)", colour = NULL, fill = NULL) +
  theme(legend.position = "right")

The GAM treats the 2,200 measurements as independent; the GAMM knows they come from 200 tokens and 20 speakers.

Reporting

We fitted a GAMM to F2 (bam(), {mgcv}; Wood, 2017) with an ordered factor for L1 (reference: English), a reference smooth over time, a difference smooth for Italian speakers, a smooth for duration and a tensor product interaction between time and duration. Random smooths over time by speaker (m = 1, k = 5) and an AR1 error model (ρ = 0.8) accounted for dependencies within speakers and trajectories.

Then report the model comparison, and show the prediction curves and the difference smooth with their confidence intervals.

Summary

  • Curve data have dependencies within trajectories and within speakers (and words); ignoring them makes intervals too narrow
  • In {mgcv}, random intercepts and slopes use bs = "re", and random smooths use bs = "fs" with m = 1
  • An AR1 error model handles autocorrelation within trajectories cheaply; check the standardised residuals
  • ti() adds a non-linear interaction between two continuous predictors; show it with prediction curves at a few values
  • Test effects by comparing ML-fitted models with the same random effects and error model

References

  • Coretta, S. (2024). tidygam: Tidy prediction and plotting of generalised additive models (R package version 1.0.0). https://github.com/stefanocoretta/tidygam
  • Sóskuthy, M. (2017). Generalised additive mixed models for dynamic analysis in linguistics: A practical introduction. arXiv, 1703.05339. https://arxiv.org/abs/1703.05339
  • Sóskuthy, M. (2018). Introducing generalised additive models and Becoming a GAM(M) warrior [Workshop slides]. Analysing Curves, University of Manchester, 29 May 2018.
  • Sóskuthy, M. (2021). Evaluating generalised additive mixed modelling strategies for dynamic speech analysis. Journal of Phonetics, 84, 101017. https://doi.org/10.1016/j.wocn.2020.101017
  • Wieling, M. (2018). Analyzing dynamic phonetic data using generalized additive mixed modeling: A tutorial focusing on articulatory differences between L1 and L2 speakers of English. Journal of Phonetics, 70, 86–116.
  • Wood, S. N. (2017). Generalized additive models: An introduction with R (2nd ed.). CRC Press.

Sources of the animal facts

  • Prat, Y., Azoulay, L., Dor, R., & Yovel, Y. (2017). Crowd vocal learning induces vocal dialects in bats: Playback of conspecifics shapes fundamental frequency usage by pups. PLoS Biology, 15(10), e2002556. https://doi.org/10.1371/journal.pbio.2002556
  • Prat, Y., Taub, M., & Yovel, Y. (2015). Vocal learning in a social mammal: Demonstrated by isolation and playback experiments in bats. Science Advances, 1(2), e1500019. https://doi.org/10.1126/sciadv.1500019
  • Prat, Y., Taub, M., & Yovel, Y. (2016). Everyday bat vocalizations contain information about emitter, addressee, context, and behavior. Scientific Reports, 6, 39419. https://doi.org/10.1038/srep39419

Slides and code

QR code linking to these slides

These slides

mcanzi.de/ss27_advanced_statistics/04_gamms.html

Code and all lectures

github.com/mcanzi/ss27_advanced_statistics