lme4gammgcv

How to Plot Difference Smooth with confidence Intervals using gamm4 in R

I am currently working on fitting splines with gamm4 and I have hierarchical data. In my dataset, Trial_no is nested within VP, and SoundType (levels: STA and DEV) is a level-1 predictor. My model should include a random intercept and a random slope for SoundType, and the spline curve should be allowed to differ between SoundType levels. Here is an overview of my data:

head(data)
#     VP SoundType       RT Trial_no
# 2 1001       STA 1008.610        5
# 3 1001       DEV  847.140        6
# 5 1001       STA  982.846        8
# 6 1001       STA  749.330        9
# 7 1001       STA  965.735       10
# 8 1001       STA 1198.200       11

I have been reading the passages in the paper by Sóskuthy, M. (2017) (https://doi.org/10.48550/arXiv.1703.05339 ) about different inference methods (pages 18/19) and would like to plot the difference smooth with a confidence interval. The implementation for a specific dataset is explained on page 23 and pages 28/29. Based on my understanding, here is how I would implement it for my data using bam:

data$SoundType.ordered <- as.ordered(data$SoundType) 
contrasts(data$SoundType.ordered) <- "contr.treatment"
ME_spline.diff.bam <- bam(RT ~ SoundType.ordered +
                          s(Trial_no, fx = FALSE, bs = "bs) +
                          s(Trial_no, by = SoundType.ordered, fx = FALSE, bs = "bs") +
                          s(VP, bs = "re") +
                          s(VP, SoundType.ordered, bs = "re"),
                          method = "REML",
                          data = data)

library(itsadug)
plot_diff(ME_spline.diff.bam, view = "Trial_no",
          comp = list(SoundType.ordered = c("DEV", "STA")))

Now, I want to do this using gamm4 instead of bam.

  1. Is there a way to plot the difference smooth with a confidence interval using gamm4?

  2. I have already tried fitting the model with the difference smooth. Is this correct?

ME_spline.diff <- gamm4(RT ~ SoundType.ordered +
                        s(Trial_no, fx = FALSE, bs = "bs") +
                        s(Trial_no, by = SoundType.ordered, fx = FALSE, bs = "bs"),
                        random = ~(1 + SoundType.ordered | VP),
                        data = data)

This is my first time asking a question here. If I forgot to include any necessary information, please let me know.

Thank you!

Solution

  • One option (although I'd be surprised it itsadug::plot_diff() didn't work) is the difference_smooths() function in my gratia package:

    library("mgcv")
library("gamm4")
library("gratia")

df <- data_sim("eg4", seed = 42)
m <- gamm4(y ~ fac + s(x2, by = fac) + s(x0), data = df, REML = TRUE)

sm_diff <- difference_smooths(m, select = "s(x2)")

sm_diff

sm_diff |>
  draw()

    enter image description here

    If you don't like the way I chose to plot this, you can always cook up your own plot using the output from difference_smooths() as it is designed to work with ggplot2

    # A tibble: 300 × 9
   .smooth .by   .level_1 .level_2 .diff   .se .lower_ci .upper_ci      x2
   <chr>   <chr> <chr>    <chr>    <dbl> <dbl>     <dbl>     <dbl>   <dbl>
 1 s(x2)   fac   1        2        0.386 0.618  -0.824        1.60 0.00359
 2 s(x2)   fac   1        2        0.479 0.574  -0.646        1.60 0.0136
 3 s(x2)   fac   1        2        0.572 0.534  -0.474        1.62 0.0237
 4 s(x2)   fac   1        2        0.665 0.497  -0.308        1.64 0.0338
 5 s(x2)   fac   1        2        0.758 0.464  -0.151        1.67 0.0438
 6 s(x2)   fac   1        2        0.850 0.435  -0.00342      1.70 0.0539
 7 s(x2)   fac   1        2        0.941 0.412   0.134        1.75 0.0639
 8 s(x2)   fac   1        2        1.03  0.393   0.262        1.80 0.0740
 9 s(x2)   fac   1        2        1.12  0.378   0.380        1.86 0.0841
10 s(x2)   fac   1        2        1.21  0.367   0.489        1.93 0.0941
# ℹ 290 more rows
# ℹ Use `print(n = ...)` to see more rows

    For example, to put all the difference smooths on a single plot we might use:

    library("dplyr")
library("ggplot2")

sm_diff |>
  mutate(.comp = factor(paste(.level_1, "vs", .level_2))) |>
  ggplot(aes(x = x2, y = .diff, .group = .comp)) +
  geom_ribbon(aes(ymin = .lower_ci, ymax = .upper_ci, fill = .comp),
    alpha = 0.2) +
  geom_line(aes(colour = .comp)) +
  labs(y = "Difference", colour = "Comparison", fill = "Comparison")

    enter image description here