Generalized Addditive Mixed Models (GAMMs)
GAMMSs can used to characterize non-linear changes in pupil diameter across trial’s to identify if changes in pupil diameter are unique to the type of stimuli presented. In this study we were interested in testing if pupil diameter and essentially if cognitive effort is different between trials with images with faces and trials with images without faces. Next a multiple comparison with false discovery rate correction is used to better interpret the GAMM and how changes in pupil diameters differed between trial types. This approach identifies specific time segments where pupil diameters were significantly different between both image types by performing pairwise comparisons at each time point, while also controlling for Type 1 errors that arise when computing multiple comparisons simultaneously (Benjamini & Hochberg, 1995; Muncy et al., 2022).
The final GAMM model was fitted to predict mean pupil diameter as a function of image type (non-social versus social), child age, sex, accuracy and smooths of time by image type, and participants included as a random effect. A scaled t-distribution (scat) was identified as the appropriate family type for the heavy-tailed distribution seen in the data. The fast residual estimates of maximum likelihood (fREML) method, which completes computations faster with large datasets, was used to assess overall model fit. Lastly, an AR1 model when rho is set to an optimal value of .1 was used to address autocorrelation. For a complete overview of GAMMs refer to Wood’s (2017) book on generalized additive models.
#Load libraries
library("mgcv")
library("dplyr")
library("ggplot2")
library("scales")
library("gridExtra")
library("fitdistrplus")Identifying appropriate family type:
The function descdist() from fitdistrplus package, generates a Cullen and Frey skewness-kurtosis plot to help guide the choice of which parametric distributions might be suitable to fit data. The data used in the example below does not fall within the listed theroetical parametric distributions. Therefore, we needed to run further tests to determeing if the data is heavey tailed.
#Cullen and Frey skewness-kurtosis plot
descdist(data$mean_pupil, discrete=F)
summary statistics
------
min: -13.88479 max: 9.665689
median: 0.4280267
mean: 0.5260173
estimated sd: 0.9952966
estimated skewness: -0.9272713
estimated kurtosis: 29.195
> 
The function qqnorm(), built-in part of R’s base stats package, can be used to generates a normal quantile-quantile (Q-Q) plot to visually assess whether the dataset follows a theoretical normal distribution.
If the points fall approximately along a straight line, it suggests the data is likely normally distributed.”S” shape suggests the distribution has longer tails (heavy tails) than a normal distribution.”U” shape implies one distribution is skewed relative to the other. Points bent down on the left and up on the right specifically mean the data has heavier tails than a Gaussian distribution (bell shaped).
The data used in this study is heavy tailed and should be addressed by using scaled t family parameters.
#Is the data heavy tailed?
qqnorm(data_1$mean_pupil)
Addressing autocorrelation:
1: Checking for autocorrelation
Autocorrelation is what happens when current values depend on the previous values, which is common in time sereis data. One way to check for autocorrelation in the data is by running a Box-Ljung test. If p.vlaue is >.05 there is no evidence of autocorrelation If p.value is <.05 there is evidice of significant autocorrelation. X-square is that test statistic and df is number of lags tested. Therefore in this case, there is signigicant autocorrelation and needs to be addressed. THere are varipsu ways to address autocorrelaton in GAMMs, but for this one an AR.Start time vairable was created.
# Box - Ljung test
Box.test(data_1$mean_pupil, lag = 10, type = "Ljung-Box")
#Output
>
Box-Ljung test
data: data_1$mean_pupil
X-squared = 14017, df = 10, p-value < 2.2e-16
> 2: Creating an AR.Start time variable
An AR.start time variable is a logical indicator variable that informs the model which row starts a new independent time series. Without this the model would assume the last time point in trial 1 is correlated with the first time point in trial 2. With it each new trial starts a new time series and autocorrelation should reset. Next we need to determine how much the previpus time point should influcne the next one within a trial by calcuating the optimal rho value.
#Create AR. Start time variable
data_1 <- data_1 %>%
group_by(id, trial) %>%
arrange(time) %>%
mutate(AR.start = row_number() == 1) %>%
ungroup()3: calcuating optimal rho value
The code below is a loop that tests all values between 0 to .9 by .1 to determine which rho value is optimal based on fREML scores. Lower fREML scores indicate better model fits. In this case we set Rho to .1.
> # Get Rho value for model
> rhos <- seq(0, .9, by = 0.1)
>
> models <- lapply(rhos, function(r) {
+ bam(
+ mean_pupil ~ social + ageinmonth + sex + correct +
+ s(time, by = social) + #separate smooth function of time for each condition
+ s(id, bs = "re"),
+ data = data_1,
+ method = "fREML",
+ rho = r,
+ AR.start = data$AR.start
+ )
+ })
>
> fREML_scores <- sapply(models, function(m) m$gcv.ubre)
>
> rho_comparison <- data.frame(
+ rho = rhos,
+ fREML = fREML_scores
+ )
>
> print(rho_comparison)
rho fREML
1 0.0 236775.8
2 0.1 238546.5
3 0.2 241916.1
4 0.3 246684.4
5 0.4 252602.0
6 0.5 259405.2
7 0.6 266847.8
8 0.7 274725.9
9 0.8 282904.0
10 0.9 291387.4
> Final GAMM Model
The model’s overall fit was evaluated with the gam.check(aim2_model, rep =100) function (see Figure 3). The model’s rank was adequate (rank = 180/180) and successfully converged. Basis dimension (k) checks identified that the specified basis dimensions (k = 8) are sufficient for all smooth terms, with no evidence of under or over smoothing. Specifically, the smooths for time by image type have large enough p-values to confirm that residuals are randomly distributed (non-social: edf = 6.97, k′ = 7, p = 0.40; social: edf = 6.96, k′ = 7, p = 0.42). The random smooth for participants (edf = 156.91, k′ = 161) is also sufficiently flexibility. Overall, the model was well-specified, with no evidence that the basis dimensions were too low or too high. Next, we check the model’s output statistics using the summary() function.
# Recommend setting plotting window to 2x2 to be able to see 4 graphs at a time
par(mfrow = c(2, 2))
#set seed to make make random processes reproducible
> set.seed(456)
> mod <- bam(mean_pupil ~ social + ageinmonth + sex + correct +
+ s(time, by = social, k = 8) + #separate smooth function of time for each condition
+ s(id, bs = "re"), #participant as a random intercept allows pupil size to vary across individuals
+ data = data_1,
+ discrete = TRUE,
+ method = "fREML",
+ family = scat(), #scaled t-distributions since data is heavy tailed
+ rho = .1, # autocorrelation parameter for residuals, expected correlation between consecutive time points
+ AR.start = data_1$AR.start) #start of each trial
>
> set.seed (456)
> gam.check(mod, rep =100)
Method: fREML Optimizer: perf chol
$grad
[1] 7.322859e-06 -1.754171e-05 2.343872e-03
$hess
[,1] [,2] [,3]
[1,] 2.9235275292 -9.651005e-03 3.619741e-04
[2,] -0.0096510054 2.894344e+00 -2.719821e-05
[3,] 0.0003619741 -2.719821e-05 7.762774e+01
Model rank = 180 / 180
Basis dimension (k) checking results. Low p-value (k-index<1) may
indicate that k is too low, especially if edf is close to k'.
k' edf k-index p-value
s(time):socialnonsocial 7.00 6.97 1 0.41
s(time):socialsocial 7.00 6.95 1 0.43
s(id) 161.00 156.71 NA NA
The model indicates significant effects of trial type, correctness, and the smooth terms for time. Compared to non-social images (intercept), social images had a significantly larger change in pupil diameter across time (b = 0.022, SE = 0.003, t = 7.379, p < .001). Accuracy was a significant predictor (b = -0.034, SE = 0.004, t = -9.060, p < .001), indicating that mean change in pupil diameter was smaller for correct trials across time. Age (b = 0.001, SE = 0.002, t = -0.360, p = .719) and sex (b = -0.062, SE = 0.039, t = -1.599, p = .110) were not significant predictors. Additionally, regardless of image type, time has a significant non-linear effect on pupil diameter (non-social: edf = 6.97, F = 1856.5, p < .001; social: edf = 6.96, F = 1322.3, p < .001). The model also found strong evidence for individual differences (edf = 156.906, F = 165.5, p < .001). Overall, the model accounted for 14.1% of the deviance explained in mean change in pupil diameter (adjusted R² = .126).
> summary(mod)
Family: Scaled t(3,0.509)
Link function: identity
Formula:
mean_pupil ~ social + ageinmonth + sex + correct + s(time, by = social,
k = 8) + s(id, bs = "re")
Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.5903756 0.1047504 5.636 1.74e-08 ***
socialsocial 0.0252250 0.0030789 8.193 2.57e-16 ***
ageinmonth -0.0005881 0.0015437 -0.381 0.703
sex1 -0.0623442 0.0382230 -1.631 0.103
correct1 -0.0310008 0.0036887 -8.404 < 2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Approximate significance of smooth terms:
edf Ref.df F p-value
s(time):socialnonsocial 6.971 7.000 1671.1 <2e-16 ***
s(time):socialsocial 6.954 6.999 1223.2 <2e-16 ***
s(id) 156.709 159.000 717.9 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
R-sq.(adj) = 0.125 Deviance explained = 14.1%
fREML = 3.0505e+05 Scale est. = 1 n = 177874Plotting GAMM
The function plot.gam() can be used to quickly visulize the model. The code below was created to create a plot that can be used for presentation and publications. # Generate a prediction grid
newdata <- expand.grid(
time = seq(0, 6500, length.out = 200),
social = c("social", "nonsocial"),
ageinmonth = mean(data_1$ageinmonth, na.rm = TRUE),
sex = names(sort(table(data_1$sex), decreasing = TRUE))[1],
correct = names(sort(table(data_1$correct), decreasing = TRUE))[1],
id = data_1$id[1] # dummy id for prediction
)
# Predict values
preds <- predict(mod, newdata = newdata, se.fit = TRUE, exclude = "s(id)")
# Add predictions to the new data
newdata$fit <- as.numeric(preds$fit)
newdata$se <- as.numeric(preds$se.fit)
newdata$lower <- newdata$fit - 1.96 * newdata$se
newdata$upper <- newdata$fit + 1.96 * newdata$se
newdata$fitR <- round(as.numeric(newdata$fit), 3)
newdata$lowerR <- round(as.numeric(newdata$lower), 3)
newdata$upperR <- round(as.numeric(newdata$upper), 3)
facet_labels <- c(
"Social" = "Social\u00B0",
"NonSocial" = "NonSocial\u00B0"
)
predicted_plot_facets <- ggplot(newdata, aes(x = time, y = fit, color = social, fill = social)) +
geom_line(linewidth = 1.2) +
geom_ribbon(aes(ymin = lower, ymax = upper), alpha = 0.2, color = NA) +
labs(
x = "Time(ms)",
y = ("Δ Pupil Size from Baseline (mm)")) +
scale_color_manual(
values = c("social" = "#C3593F", "nonsocial" = "#F9E394")) +
scale_fill_manual(
values = c("social" = "#C3593F", "nonsocial" = "#F9E394")) +
scale_y_continuous(
breaks = seq(-0.2,0.8, by = 0.2),
labels = number_format(accuracy = 0.2)
) +
theme_minimal() +
facet_wrap(~social, ncol =1, labeller = labeller(social = facet_labels)) +
theme(panel.grid.minor = element_blank(),
axis.title.x = element_text(size = 12), # no x-axis title on first plot
axis.title.y = element_text(size =12),
legend.position = "none",
strip.text = element_text(size = 12, hjust = 0.5),
panel.spacing = unit(0.1, "lines"),
plot.margin = unit(c(0.2, 0.1, 0.2, 0.1), "cm"))
predicted_plot_facets
ggsave("predicted_plot_facets10.23.png", plot = predicted_plot_facets,
width = 4, height = 8, dpi = 600, units = "in")
# save as SVG file
ggplot2:: ggsave("predicted_plot_facets11.05.svg",
plot = predicted_plot_facets,
device = "svg",
width = 4, height = 8, units = "in")
predicted_plot_combined <- ggplot(newdata, aes(x = time, y = fit, color = social, fill = social)) +
geom_line(linewidth = 1.2) +
geom_ribbon(aes(ymin = lower, ymax = upper), alpha = 0.2, color = NA) +
labs(
x = "Time(ms)",
y = "NA",
subtitle = "All Trial Types" ) +
scale_color_manual(
name = "Trial Type",
values = c("social" = "#C3593F", "nonsocial" = "#F9E394")) +
scale_fill_manual(
name = "Trial Type",
values = c("social" = "#C3593F", "nonsocial" = "#F9E394")) +
scale_y_continuous(
breaks = seq(-0.2,0.8, by = 0.2),
labels = number_format(accuracy = 0.2)) +
theme_minimal() +
theme(panel.grid.minor = element_blank(),
axis.title.x = element_text(size = 12),
axis.title.y = element_blank(),
plot.subtitle = element_text(hjust = 0.5, vjust = -1, size = 12),
legend.position = "none",
plot.margin = unit(c(0.2, 0.2, 0.2, 0.1), "cm"))
predicted_plot_combined
ggsave("predicted_plot_combined.10.23.png", plot = predicted_plot_combined,
width = 6, height = 8, dpi = 600, units = "in")
# save as SVG file
ggplot2:: ggsave("predicted_plot_combined11.04.svg",
plot = predicted_plot_combined,
device = "svg",
width = 6, height = 8, units = "in")
Multiple Comparison / FDR correction
A multiple comparison with false discovery rate correction was used to better interpret the model and how changes in pupil diameters differed between image types. This approach identifies specific time segments where pupil diameters were significantly different between both image types by performing pairwise comparisons at each time point, while also controlling for Type 1 errors that arise when computing multiple comparisons simultaneously (Benjamini & Hochberg, 1995; Curtis et al., n.d.; Muncy et al., 2022). To minimize the likelihood of falsely identifying significant time periods, time points after 4500 ms were excluded because of the overfitted predictions observed in Figure 3 and theoretical basis that changes greater than .5 mm are due to light. Results shown in Figure 4, indicate when comparing social and non-social images there is significant time cluster spanning from 1500 ms to 1980 ms (mean estimate = -0.116, FDR-corrected p < .05).
> # Filter time and ensure social is a factor
> data_filtered <- data_1 %>%
+ filter(time >= 0 & time <= 4500) %>%
+ mutate(social = as.factor(social))
>
> # Calculate sample size and degrees of freedom
> n_participants <- length(unique(data_filtered$id))
> n_social_levels <- length(unique(data_filtered$social))
> df_between <- n_participants - n_social_levels
>
> print(paste("Number of unique participants:", n_participants))
[1] "Number of unique participants: 161"
> print(paste("Degrees of freedom for between-subjects comparisons:", df_between))
[1] "Degrees of freedom for between-subjects comparisons: 159"
>
> # Create time points for comparison
> time_points <- seq(0, 4500, by = 45)
>
> # Reference values for covariates
> ref_ageinmonth <- mean(data_filtered$ageinmonth, na.rm = TRUE)
> ref_sex <- names(sort(table(data_filtered$sex), decreasing = TRUE))[1]
> ref_correct <- names(sort(table(data_filtered$correct), decreasing = TRUE))[1]
> dummy_id <- data_filtered$id[1]
>
> # Ensure social levels match data
> social_levels <- levels(factor(data_filtered$social))
>
> # Prediction grid
> pred_grid <- expand.grid(
+ social = factor(social_levels, levels = social_levels),
+ time = time_points,
+ ageinmonth = ref_ageinmonth,
+ sex = factor(ref_sex, levels = levels(factor(data_filtered$sex))),
+ correct = factor(ref_correct, levels = levels(factor(data_filtered$correct))),
+ id = factor(dummy_id, levels = levels(factor(data_filtered$id)))
+ )
>
> # Predictions from model (exclude random effect)
> predictions <- predict(mod, newdata = pred_grid, exclude = "s(id)", se.fit = TRUE)
>
> # Add predicted values and SE
> pred_grid$fit <- predictions$fit
> pred_grid$se <- predictions$se.fit
>
> # Social levels
> social_conditions <- sort(unique(data_filtered$social))
> print(paste("Social levels found:", paste(social_conditions, collapse = ", ")))
[1] "Social levels found: nonsocial, social"
>
> # Function for pairwise difference
> calc_pairwise_diff <- function(fit1, se1, fit2, se2) {
+ diff <- fit1 - fit2
+ se_diff <- sqrt(se1^2 + se2^2)
+ return(list(diff = diff, se = se_diff))
+ }
>
> # Pairwise comparisons at each time point
> comparison_results <- list()
>
> for(i in seq_along(time_points)) {
+ current_time <- time_points[i]
+ time_data <- pred_grid[pred_grid$time == current_time, ]
+
+ comparisons_this_time <- list()
+
+ # Only two levels: social vs nonsocial
+ fit1 <- time_data$fit[time_data$social == social_conditions[1]]
+ se1 <- time_data$se[time_data$social == social_conditions[1]]
+ fit2 <- time_data$fit[time_data$social == social_conditions[2]]
+ se2 <- time_data$se[time_data$social == social_conditions[2]]
+
+ diff_result <- calc_pairwise_diff(fit1, se1, fit2, se2)
+ t_stat <- diff_result$diff / diff_result$se
+ p_value <- 2 * pt(abs(t_stat), df = df_between, lower.tail = FALSE)
+
+ comparisons_this_time[[1]] <- data.frame(
+ time = current_time,
+ comparison = paste(social_conditions[1], "vs", social_conditions[2]),
+ estimate = diff_result$diff,
+ se = diff_result$se,
+ t_value = t_stat,
+ df = df_between,
+ p_value = p_value
+ )
+
+ comparison_results[[i]] <- do.call(rbind, comparisons_this_time)
+ }
>
> # Combine all results
> all_comparisons <- do.call(rbind, comparison_results)
>
> # Apply FDR correction
> all_comparisons$p_fdr <- p.adjust(all_comparisons$p_value, method = "fdr")
> all_comparisons$significant_fdr <- all_comparisons$p_fdr < 0.05
>
> # Summary of significant time windows
> significant_periods <- all_comparisons %>%
+ filter(significant_fdr == TRUE) %>%
+ arrange(time)
>
> print("=== SIGNIFICANT DIFFERENCES AFTER FDR CORRECTION ===")
[1] "=== SIGNIFICANT DIFFERENCES AFTER FDR CORRECTION ==="
> print(significant_periods)
time comparison estimate se t_value df p_value p_fdr significant_fdr
1 1170 nonsocial vs social -0.1038711 0.03811660 -2.725088 159 0.007148771 0.04011255 TRUE
2 1215 nonsocial vs social -0.1080271 0.03808899 -2.836176 159 0.005159981 0.03257238 TRUE
3 1260 nonsocial vs social -0.1118272 0.03806812 -2.937557 159 0.003800094 0.02741496 TRUE
4 1305 nonsocial vs social -0.1152292 0.03805536 -3.027936 159 0.002873743 0.02418734 TRUE
5 1350 nonsocial vs social -0.1181928 0.03805134 -3.106140 159 0.002245181 0.02205680 TRUE
6 1395 nonsocial vs social -0.1206811 0.03805594 -3.171149 159 0.001822216 0.02156808 TRUE
7 1440 nonsocial vs social -0.1226608 0.03806831 -3.222123 159 0.001543641 0.02156808 TRUE
8 1485 nonsocial vs social -0.1241032 0.03808701 -3.258414 159 0.001370029 0.02156808 TRUE
9 1530 nonsocial vs social -0.1249847 0.03811012 -3.279566 159 0.001277405 0.02156808 TRUE
10 1575 nonsocial vs social -0.1252869 0.03813544 -3.285315 159 0.001253260 0.02156808 TRUE
11 1620 nonsocial vs social -0.1249979 0.03816068 -3.275567 159 0.001294457 0.02156808 TRUE
12 1665 nonsocial vs social -0.1241117 0.03818367 -3.250386 159 0.001406789 0.02156808 TRUE
13 1710 nonsocial vs social -0.1226290 0.03820249 -3.209973 159 0.001606194 0.02156808 TRUE
14 1755 nonsocial vs social -0.1205573 0.03821568 -3.154655 159 0.001921908 0.02156808 TRUE
15 1800 nonsocial vs social -0.1179106 0.03822228 -3.084866 159 0.002402226 0.02205680 TRUE
16 1845 nonsocial vs social -0.1147097 0.03822196 -3.001148 159 0.003123917 0.02427043 TRUE
17 1890 nonsocial vs social -0.1109818 0.03821495 -2.904146 159 0.004206870 0.02832626 TRUE
18 1935 nonsocial vs social -0.1067600 0.03820206 -2.794614 159 0.005835926 0.03467227 TRUE
19 1980 nonsocial vs social -0.1020834 0.03818460 -2.673417 159 0.008291943 0.04407822 TRUE
>Plotting Multiple Comparison with FDR Correction Results
Blue shaded area indicates significant timepoints that pupil diameters were different after FDR corrections between social and non-social images. Black line represents the observed comparison results across time, and the dashed red line demonstrates the significance threshold at −log10(0.05).
# Function to identify significant time sections
get_sig_sections <- function(df) {
df <- df %>% arrange(time)
sig_rle <- rle(df$significant_fdr)
ends <- cumsum(sig_rle$lengths)
starts <- c(1, head(ends, -1) + 1)
sections <- tibble(
start = df$time[starts],
end = df$time[ends],
sig = sig_rle$values
) %>%
mutate(comparison = unique(df$comparison))
return(sections)
}
# Apply function per comparison
sig_sections <- all_comparisons %>%
group_by(comparison) %>%
group_split() %>%
purrr::map_dfr(get_sig_sections)
# Split data for plotting
comparison_dfs <- split(all_comparisons, all_comparisons$comparison)
sig_section_split <- split(sig_sections, sig_sections$comparison)
# Define labels and colors for social comparison
comp_labels <- c("nonsocial vs social" = "Nonsocial vs Social")
comp_colors <- c("nonsocial vs social" = "#6BAED6") # blue tone for social comparison
# Generate plots
comparison_plots <- Map(function(df, rects) {
comp <- unique(df$comparison)
comp_title <- comp_labels[comp]
fill_color <- comp_colors[comp]
p <- ggplot(df, aes(x = time, y = -log10(p_fdr))) +
geom_line(color = "black", linewidth = 1) +
geom_hline(yintercept = -log10(0.05), linetype = "dashed", color = "red", linewidth = .9) +
labs(
x = "Time (ms)",
y = "FDR-corrected p-values (-log10 scale)",
title = comp_title
) +
theme_classic() +
theme(
axis.title.x = element_text(size = 10),
axis.title.y = element_text(size = 10),
plot.title = element_text(size = 12, hjust = 0.5)
)
# Add shaded regions for significant time windows
if (nrow(rects) > 0 && any(rects$sig == TRUE, na.rm = TRUE)) {
p <- p + geom_rect(
data = dplyr::filter(rects, sig == TRUE),
aes(xmin = start, xmax = end, ymin = -Inf, ymax = Inf),
fill = fill_color, alpha = 0.4, inherit.aes = FALSE
)
}
return(p)
}, comparison_dfs, sig_section_split)
# Combine plot (only one comparison here)
final_plot <- comparison_plots[[1]]
print(final_plot)
ggsave("MC_FDR_Ch3.12.12.png", plot = final_plot,
width = 10, height = 8, dpi = 600)
ggplot2:: ggsave("MC_FDR_Ch3.12.12.svg",
plot = final_plot,
device = "svg",
width = 8, height = 10, units = "in")References
Wood, S. N. (2017). Generalized Additive Models An Introduction with R (J. K. Blitzstein, J. J. Faraway, M. Tanner, & J. Zidek, Eds.; Second). Taylor & Francis Group.