How today works

Today you run a complete GAMM analysis yourselves, in pairs, on a new data set.

• Each task slide says what to do, with code to complete where ___ appears
• A checkpoint slide after some tasks shows what you should get
• The solutions are in the .Rmd file of these slides, on GitHub. Try first; look only when you are stuck

The study

A simulated pupillometry experiment on listening effort. 24 participants heard 12 sentences in quiet and 12 in background noise. Pupil size was recorded from sentence onset, baseline-corrected and averaged into 100 ms bins from 0 to 3,000 ms.

Single trials

Four trials per condition from participant P01. Single trials are noisy and differ in height: the random smooths and the AR1 model will have to deal with both.

Task 1: look at the data

1. Run the code on the previous slide and look at pdata. How many rows are there, and how many trajectories?
2. Compute the mean pupil size per participant, condition and time bin
3. Plot one line per participant, with a panel for each condition

Checkpoint 1

17856 rows and 576 trajectories. In noise, the response tends to be larger and later, with a lot of variation between participants.

Task 2: a first GAM

pdata$condition_o <- as.ordered(pdata$condition)
contrasts(pdata$condition_o) <- "contr.treatment"
m_gam <- bam(pupil ~ ___ + s(time, k = 20) + s(time, by = ___, k = 20),
             data = pdata, method = "ML")
summary(m_gam)

• Complete the formula with the ordered factor, as in lecture 3. Which level is the reference?
• What does each smooth in the summary describe?
• Is k = 20 enough? Check with k.check()

Task 3: is anything left in the residuals?

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()
}

The same function as in lecture 4. Compute the autocorrelation of the residuals of m_gam within trajectories, and plot it against the lag.

Task 4: random smooths

Every participant heard both conditions, so the comparison is within participants. Participants can differ in their curves and also in how noise changes their curves.

pdata$part_cond <- interaction(pdata$participant, pdata$condition)
m_rs <- bam(pupil ~ condition_o + s(time, k = 20) + s(time, by = condition_o, k = 20) +
              s(time, ___, bs = "fs", m = 1, k = 5),
            data = pdata, method = "ML")

• Why part_cond and not participant? Look back at “The random slope equivalent” in lecture 4
• How many levels does part_cond have, and how many coefficients do the random smooths add?

Task 5: an AR1 error model

rho <- itsadug::start_value_rho(m_rs)
m_ar <- bam(pupil ~ condition_o + s(time, k = 20) + s(time, by = condition_o, k = 20) +
              s(time, part_cond, bs = "fs", m = 1, k = 5),
            data = pdata, method = "ML", AR.start = ___, rho = rho)

1. Mark the first bin of every trajectory in AR.start. Is the data frame sorted the way the AR1 model needs?
2. Compare the autocorrelation of residuals(m_gam), residuals(m_rs) and m_ar$std.rsd

Checkpoint 2

Your rho should be close to 0.92. With the AR1 model, little autocorrelation should remain.

Task 6: check the model

• Run k.check(m_ar). Does any smooth need a higher k?
• Plot the standardised residuals (m_ar$std.rsd) against the fitted values, and make a QQ plot. Lecture 4 has the code
• Would you trust the confidence intervals of this model more than those of m_gam? Why?

Task 7: test the effect of noise

Fit a null model without any term for condition, and compare it with m_ar using itsadug::compareML().

Task 8: plot the results

With {tidygam}, plot:

1. the predicted curves for the two conditions, with 95% confidence intervals, for an average participant
2. the difference smooth, noise minus quiet, with a line at zero

Checkpoint 3

Task 9: interpret and report

• From roughly what time does the confidence interval of the difference exclude zero? Is that the same as “noise has an effect from that time on”?
• Write a reporting paragraph for a methods and results section, like the one at the end of lecture 4: the model, the random effects and error model with their settings, the test and its result, and a figure caption

If you finish early

The simulation has one more effect built in: the pupil response gets smaller over the 12 trials, as participants tire or get used to the task.

• Add s(trial, k = 5) to m_ar. Does the model improve?
• Then add ti(time, trial). Does fatigue only shrink the curve, or does it change its shape too?

Summary

Slides and code

These slides

mcanzi.de/ss27_advanced_statistics/05_gamm_practical.html

Code and all lectures

github.com/mcanzi/ss27_advanced_statistics