Temporal Emotional Gaze in MDD: GAMM and Two-Axis rsFC Integration
Tak Kwan Lam
2026-07-28
1 Analysis purpose
This analysis replaces the ten-bin ANOVA and the almost completely overlapping 200-ms moving windows with a hierarchical generalized additive modelling strategy. The primary outcome is the number of valid gaze samples directed to the emotional versus neutral face in non-overlapping time bins. A binomial GAMM therefore models the nonlinear probability of emotional-face gaze over time.
The script first tests the nine expression-by-starting-position
trajectories, with multiplicity control across the nine planned
group-by-time comparisons. An optional second stage merges the
subject-level resting-state scores produced by
3_resting_state_fc_two_axis_targeted_updated.R.
2 Settings
data_folder <- "C:/Users/quent/Desktop/fmri eye/data"
input_file <- file.path(data_folder, "all_subjects_full.csv")
output_dir <- file.path(data_folder, "eye_temporal_gamm_outputs")
# The original codebook states:
# AOI_tem_emo == 1 : neutral face
# AOI_tem_emo == 2 : emotional face
# The old executable code accidentally reversed these assignments.
neutral_aoi_code <- 1L
emotional_aoi_code <- 2L
sample_rate_hz <- 500
bin_width_ms <- 20
analysis_start_ms <- 0
analysis_end_ms <- 2000
minimum_valid_samples_per_bin <- 3L
k_time <- 15L
k_subject_time <- 5L
n_threads <- max(1L, parallel::detectCores(logical = TRUE) - 1L)
# Verify these labels against the task program before manuscript submission.
expression_map <- c(F = "Fear", H = "Happy", S = "Sad")
starting_position_map <- c(
C = "Center",
E = "Emotional_start",
N = "Neutral_start"
)
# Optional multimodal stage. Keep FALSE until participant overlap, AOI coding,
# starting-position coding, and the rsFC analysis have all been verified.
run_multimodal_models <- FALSE
fc_score_file <- file.path(
data_folder, "rsfc_two_axis",
"two_axis_subject_scores_for_temporal_integration.csv"
)
dir.create(output_dir, recursive = TRUE, showWarnings = FALSE)3 Packages
required_packages <- c("data.table", "mgcv", "ggplot2", "MASS")
missing_packages <- required_packages[
!vapply(required_packages, requireNamespace, logical(1L), quietly = TRUE)
]
if (length(missing_packages) > 0L) {
stop(
"Install the required package(s) before knitting: ",
paste(missing_packages, collapse = ", ")
)
}
library(data.table)
library(mgcv)## Loading required package: nlme## This is mgcv 1.9-1. For overview type 'help("mgcv-package")'.library(ggplot2)4 Load and audit the raw gaze samples
if (!file.exists(input_file)) stop("Missing eye-tracking file: ", input_file)
required_columns <- c(
"subject_id", "sub_type", "trial_type", "expression",
"Ntrial", "count", "AOI_tem_emo"
)
raw <- fread(
input_file,
select = required_columns,
showProgress = TRUE
)
if (!all(required_columns %in% names(raw))) {
stop("Input file is missing: ",
paste(setdiff(required_columns, names(raw)), collapse = ", "))
}
raw[, `:=`(
subject = trimws(as.character(subject_id)),
group_raw = toupper(trimws(as.character(sub_type))),
expression_raw = toupper(trimws(as.character(expression))),
start_raw = toupper(trimws(as.character(trial_type))),
trial = as.integer(Ntrial),
time_ms = as.numeric(count),
AOI_code = suppressWarnings(as.integer(AOI_tem_emo))
)]
raw[, group := fcase(
group_raw %in% c("HC", "1", "CONTROL"), "HC",
group_raw %in% c("MDD", "MD", "2", "PATIENT"), "MDD",
default = NA_character_
)]
raw[, expression_label := unname(expression_map[expression_raw])]
raw[, starting_position := unname(starting_position_map[start_raw])]
unmapped_expression <- unique(raw[is.na(expression_label), expression_raw])
unmapped_start <- unique(raw[is.na(starting_position), start_raw])
if (length(unmapped_expression) > 0L) {
stop("Unmapped expression code(s): ", paste(unmapped_expression, collapse = ", "))
}
if (length(unmapped_start) > 0L) {
stop(
"Unmapped trial_type code(s): ", paste(unmapped_start, collapse = ", "),
". Edit starting_position_map only after checking the task documentation."
)
}
if (anyNA(raw$group)) stop("At least one sub_type value could not be mapped to HC or MDD.")
observed_aoi_codes <- sort(unique(na.omit(raw$AOI_code)))
if (!all(c(neutral_aoi_code, emotional_aoi_code) %in% observed_aoi_codes)) {
stop(
"The configured neutral/emotional AOI codes were not both observed. ",
"Observed codes: ", paste(observed_aoi_codes, collapse = ", ")
)
}
coding_audit <- data.table(
field = c("neutral AOI", "emotional AOI", "expression codes", "trial_type codes"),
configured_value = c(
neutral_aoi_code,
emotional_aoi_code,
paste(names(expression_map), expression_map, sep = "=", collapse = ";"),
paste(names(starting_position_map), starting_position_map, sep = "=", collapse = ";")
)
)
fwrite(coding_audit, file.path(output_dir, "coding_audit.csv"))
print(coding_audit)## field configured_value
## <char> <char>
## 1: neutral AOI 1
## 2: emotional AOI 2
## 3: expression codes F=Fear;H=Happy;S=Sad
## 4: trial_type codes C=Center;E=Emotional_start;N=Neutral_start# Confirm that count behaves like milliseconds. At 500 Hz, consecutive valid
# samples should usually be about 2 ms apart within a trial.
time_step_audit <- raw[
is.finite(time_ms),
.(median_step_ms = as.numeric(median(diff(sort(unique(time_ms))), na.rm = TRUE))),
by = .(subject, trial)
]
overall_median_step_ms <- median(time_step_audit$median_step_ms, na.rm = TRUE)
fwrite(time_step_audit, file.path(output_dir, "time_step_audit_by_trial.csv"))
if (!is.finite(overall_median_step_ms) ||
abs(overall_median_step_ms - 1000 / sample_rate_hz) > 0.5) {
warning(
"The median within-trial increment of count is ",
round(overall_median_step_ms, 3),
" ms, not approximately ", 1000 / sample_rate_hz,
" ms. Confirm whether count is milliseconds or a sample index."
)
}## Warning: The median within-trial increment of count is 1 ms, not approximately
## 2 ms. Confirm whether count is milliseconds or a sample index.5 Create non-overlapping time-bin counts
A GAM already estimates a smooth trajectory, so the old 200-ms window advanced in 2-ms steps is unnecessary and creates severe redundancy. Here, each raw gaze sample contributes to one 20-ms bin only. The response retains the numbers of emotional and neutral samples rather than treating each binned mean as a Gaussian observation.
valid <- raw[
time_ms >= analysis_start_ms & time_ms < analysis_end_ms &
AOI_code %in% c(neutral_aoi_code, emotional_aoi_code)
]
valid[, time_bin := floor((time_ms - analysis_start_ms) / bin_width_ms)]
valid[, time_center_ms := analysis_start_ms + (time_bin + 0.5) * bin_width_ms]
# Aggregate first within trial, then across trials for each participant. This
# preserves a valid-sample denominator and avoids overweighting a trial merely
# because it has more rows after merging.
trial_bins <- valid[, .(
n_emotional = sum(AOI_code == emotional_aoi_code),
n_neutral = sum(AOI_code == neutral_aoi_code),
n_valid = .N
), by = .(
subject, group, expression_label, starting_position,
trial, time_bin, time_center_ms
)]
trial_bins <- trial_bins[n_valid >= minimum_valid_samples_per_bin]
subject_bins <- trial_bins[, .(
n_emotional = sum(n_emotional),
n_neutral = sum(n_neutral),
n_valid = sum(n_valid),
n_trials = uniqueN(trial)
), by = .(
subject, group, expression_label, starting_position,
time_bin, time_center_ms
)]
subject_bins[, emotional_probability := n_emotional / (n_emotional + n_neutral)]
subject_bins[, bias_proportion := (n_emotional - n_neutral) /
(n_emotional + n_neutral)]
subject_bins[, `:=`(
subject = factor(subject),
group = factor(group, levels = c("HC", "MDD")),
expression_label = factor(expression_label,
levels = c("Fear", "Happy", "Sad")),
starting_position = factor(starting_position,
levels = c("Center", "Emotional_start", "Neutral_start"))
)]
subject_bins[, condition := interaction(
expression_label, starting_position, drop = TRUE, sep = "__"
)]
fwrite(subject_bins, file.path(output_dir, "subject_time_bins_for_gamm.csv"))
sample_audit <- subject_bins[, .(
participants = uniqueN(subject),
median_trials_per_time_bin = as.numeric(median(n_trials)),
minimum_trials_per_time_bin = as.numeric(min(n_trials)),
median_valid_samples = as.numeric(median(n_valid))
), by = .(group, expression_label, starting_position)]
fwrite(sample_audit, file.path(output_dir, "temporal_sample_audit.csv"))
print(sample_audit)## group expression_label starting_position participants
## <fctr> <fctr> <fctr> <int>
## 1: HC Happy Neutral_start 47
## 2: HC Happy Center 47
## 3: HC Sad Emotional_start 47
## 4: HC Fear Neutral_start 47
## 5: HC Sad Center 47
## 6: HC Happy Emotional_start 47
## 7: HC Fear Center 47
## 8: HC Sad Neutral_start 47
## 9: HC Fear Emotional_start 47
## 10: MDD Fear Emotional_start 63
## 11: MDD Sad Neutral_start 63
## 12: MDD Fear Center 63
## 13: MDD Sad Emotional_start 63
## 14: MDD Happy Neutral_start 63
## 15: MDD Happy Center 63
## 16: MDD Sad Center 63
## 17: MDD Fear Neutral_start 63
## 18: MDD Happy Emotional_start 63
## median_trials_per_time_bin minimum_trials_per_time_bin median_valid_samples
## <num> <num> <num>
## 1: 15 4 288
## 2: 13 1 253
## 3: 15 4 280
## 4: 15 3 288
## 5: 14 1 260
## 6: 15 5 284
## 7: 14 1 266
## 8: 15 3 280
## 9: 15 5 286
## 10: 14 1 273
## 11: 14 1 269
## 12: 13 1 240
## 13: 14 1 262
## 14: 14 1 266
## 15: 12 1 234
## 16: 12 1 238
## 17: 14 1 274
## 18: 14 1 2736 Descriptive trajectories
plot_data <- subject_bins[, .(
mean_probability = mean(emotional_probability, na.rm = TRUE),
se_probability = sd(emotional_probability, na.rm = TRUE) / sqrt(.N)
), by = .(group, expression_label, starting_position, time_center_ms)]
p <- ggplot(
plot_data,
aes(time_center_ms, mean_probability, colour = group, fill = group)
) +
geom_ribbon(
aes(ymin = pmax(0, mean_probability - 1.96 * se_probability),
ymax = pmin(1, mean_probability + 1.96 * se_probability)),
alpha = 0.15, colour = NA
) +
geom_line(linewidth = 0.8) +
facet_grid(expression_label ~ starting_position) +
geom_hline(yintercept = 0.5, linetype = 2) +
theme_classic() +
labs(
x = "Time after stimulus onset (ms)",
y = "Probability of gaze to the emotional face",
colour = "Group", fill = "Group"
)
print(p)
ggsave(file.path(output_dir, "descriptive_temporal_trajectories.png"),
p, width = 12, height = 9, dpi = 300)7 Hierarchical GAMM analysis
The omnibus model allows each expression-by-starting-position condition to have its own baseline smooth. Nine numeric indicator variables then represent the additional MDD-specific smooth for each condition. The full-versus-reduced model comparison is the global group-by-time test. Condition-specific tests are examined only after the global test and are adjusted across the nine conditions.
condition_levels <- levels(subject_bins$condition)
indicator_names <- paste0("mdd_diff_", make.names(condition_levels))
for (i in seq_along(condition_levels)) {
nm <- indicator_names[i]
lev <- condition_levels[i]
subject_bins[[nm]] <- as.numeric(
subject_bins$group == "MDD" & subject_bins$condition == lev
)
}
base_terms <- c(
"group * condition",
sprintf("s(time_center_ms, by = condition, bs = 'cr', k = %d)", k_time),
sprintf("s(time_center_ms, subject, bs = 'fs', m = 1, k = %d)", k_subject_time)
)
difference_terms <- sprintf(
"s(time_center_ms, by = %s, bs = 'cr', k = %d)",
indicator_names, k_time
)
reduced_formula <- as.formula(paste(
"cbind(n_emotional, n_neutral) ~",
paste(base_terms, collapse = " + ")
))
full_formula <- as.formula(paste(
"cbind(n_emotional, n_neutral) ~",
paste(c(base_terms, difference_terms), collapse = " + ")
))
# ML is used for the nested model comparison; fREML is used for the final model.
global_reduced_ml <- bam(
reduced_formula,
family = binomial(link = "logit"),
data = subject_bins,
method = "ML",
discrete = FALSE,
nthreads = n_threads
)
global_full_ml <- bam(
full_formula,
family = binomial(link = "logit"),
data = subject_bins,
method = "ML",
discrete = FALSE,
nthreads = n_threads
)
global_comparison <- anova(global_reduced_ml, global_full_ml, test = "Chisq")
print(global_comparison)## Analysis of Deviance Table
##
## Model 1: cbind(n_emotional, n_neutral) ~ group * condition + s(time_center_ms,
## by = condition, bs = "cr", k = 15) + s(time_center_ms, subject,
## bs = "fs", m = 1, k = 5)
## Model 2: cbind(n_emotional, n_neutral) ~ group * condition + s(time_center_ms,
## by = condition, bs = "cr", k = 15) + s(time_center_ms, subject,
## bs = "fs", m = 1, k = 5) + s(time_center_ms, by = mdd_diff_Fear__Center,
## bs = "cr", k = 15) + s(time_center_ms, by = mdd_diff_Happy__Center,
## bs = "cr", k = 15) + s(time_center_ms, by = mdd_diff_Sad__Center,
## bs = "cr", k = 15) + s(time_center_ms, by = mdd_diff_Fear__Emotional_start,
## bs = "cr", k = 15) + s(time_center_ms, by = mdd_diff_Happy__Emotional_start,
## bs = "cr", k = 15) + s(time_center_ms, by = mdd_diff_Sad__Emotional_start,
## bs = "cr", k = 15) + s(time_center_ms, by = mdd_diff_Fear__Neutral_start,
## bs = "cr", k = 15) + s(time_center_ms, by = mdd_diff_Happy__Neutral_start,
## bs = "cr", k = 15) + s(time_center_ms, by = mdd_diff_Sad__Neutral_start,
## bs = "cr", k = 15)
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 60125 1861267
## 2 60003 1838074 122.17 23193 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1write.csv(
as.data.frame(global_comparison),
file.path(output_dir, "global_group_by_time_model_comparison.csv")
)
global_final <- bam(
full_formula,
family = binomial(link = "logit"),
data = subject_bins,
method = "fREML",
discrete = TRUE,
nthreads = n_threads
)
saveRDS(global_final, file.path(output_dir, "global_temporal_gamm.rds"))
writeLines(capture.output(summary(global_final)),
file.path(output_dir, "global_temporal_gamm_summary.txt"))
pearson_residuals <- residuals(global_final, type = "pearson")
dispersion_ratio <- sum(pearson_residuals^2, na.rm = TRUE) /
max(1, df.residual(global_final))
fwrite(
data.table(metric = "Pearson dispersion ratio", value = dispersion_ratio),
file.path(output_dir, "global_gamm_dispersion_check.csv")
)
if (is.finite(dispersion_ratio) && dispersion_ratio > 1.5) {
warning(
"The binomial GAMM shows substantial overdispersion (ratio = ",
round(dispersion_ratio, 2),
"). Treat p values cautiously and fit a trial-level or overdispersion-robust sensitivity model."
)
}## Warning: The binomial GAMM shows substantial overdispersion (ratio = 29.53).
## Treat p values cautiously and fit a trial-level or overdispersion-robust
## sensitivity model.writeLines(capture.output(k.check(global_final)),
file.path(output_dir, "global_temporal_gamm_basis_check.txt"))8 Condition-specific group-by-time tests
fit_condition_models <- function(condition_value) {
dat <- subject_bins[condition == condition_value]
dat[, group_ordered := ordered(group, levels = c("HC", "MDD"))]
reduced <- bam(
cbind(n_emotional, n_neutral) ~ group +
s(time_center_ms, bs = "cr", k = k_time) +
s(time_center_ms, subject, bs = "fs", m = 1, k = k_subject_time),
family = binomial(link = "logit"),
data = dat,
method = "ML",
discrete = FALSE,
nthreads = n_threads
)
# For an ordered factor, this smooth estimates the MDD-minus-HC difference
# over time while the common smooth represents the HC trajectory.
full <- bam(
cbind(n_emotional, n_neutral) ~ group +
s(time_center_ms, bs = "cr", k = k_time) +
s(time_center_ms, by = group_ordered, bs = "cr", k = k_time) +
s(time_center_ms, subject, bs = "fs", m = 1, k = k_subject_time),
family = binomial(link = "logit"),
data = dat,
method = "ML",
discrete = FALSE,
nthreads = n_threads
)
cmp <- anova(reduced, full, test = "Chisq")
cmp_df <- as.data.frame(cmp)
p_col <- grep("Pr", names(cmp_df), value = TRUE)
p_value <- if (length(p_col) > 0L) tail(cmp_df[[p_col[1L]]], 1L) else NA_real_
final <- bam(
formula(full),
family = binomial(link = "logit"),
data = dat,
method = "fREML",
discrete = TRUE,
nthreads = n_threads
)
list(
condition = condition_value,
p_value = p_value,
reduced = reduced,
full_ml = full,
final = final,
comparison = cmp_df
)
}
cell_fits <- lapply(condition_levels, fit_condition_models)## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has repeated
## 1-d smooths of same variable.
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has repeated
## 1-d smooths of same variable.
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has repeated
## 1-d smooths of same variable.
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has repeated
## 1-d smooths of same variable.
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has repeated
## 1-d smooths of same variable.
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has repeated
## 1-d smooths of same variable.
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has repeated
## 1-d smooths of same variable.
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has repeated
## 1-d smooths of same variable.
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has repeated
## 1-d smooths of same variable.
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has repeated
## 1-d smooths of same variable.
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has repeated
## 1-d smooths of same variable.## Warning in bgam.fit(G, mf, chunk.size, gp, scale, gamma, method = method, :
## algorithm did not converge## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has repeated
## 1-d smooths of same variable.
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has repeated
## 1-d smooths of same variable.
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has repeated
## 1-d smooths of same variable.
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has repeated
## 1-d smooths of same variable.
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has repeated
## 1-d smooths of same variable.
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has repeated
## 1-d smooths of same variable.
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has repeated
## 1-d smooths of same variable.
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has repeated
## 1-d smooths of same variable.
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has repeated
## 1-d smooths of same variable.
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has repeated
## 1-d smooths of same variable.
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has repeated
## 1-d smooths of same variable.
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has repeated
## 1-d smooths of same variable.
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has repeated
## 1-d smooths of same variable.
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has repeated
## 1-d smooths of same variable.
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has repeated
## 1-d smooths of same variable.
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has repeated
## 1-d smooths of same variable.cell_tests <- rbindlist(lapply(cell_fits, function(x) {
parts <- tstrsplit(x$condition, "__", fixed = TRUE)
data.table(
condition = x$condition,
expression = parts[[1L]],
starting_position = parts[[2L]],
p_group_by_time = x$p_value
)
}))
cell_tests[, q_bh_nine_conditions := p.adjust(p_group_by_time, method = "BH")]
cell_tests[, global_gate_note :=
"Interpret only after the global full-vs-reduced group-by-time test."]
fwrite(cell_tests, file.path(output_dir, "condition_specific_group_by_time_tests.csv"))
print(cell_tests)## condition expression starting_position p_group_by_time
## <char> <char> <char> <num>
## 1: Fear__Center Fear Center 3.332404e-78
## 2: Happy__Center Happy Center 1.212524e-15
## 3: Sad__Center Sad Center 1.769290e-35
## 4: Fear__Emotional_start Fear Emotional_start 2.000604e-01
## 5: Happy__Emotional_start Happy Emotional_start 1.857202e-22
## 6: Sad__Emotional_start Sad Emotional_start 5.009085e-85
## 7: Fear__Neutral_start Fear Neutral_start 2.241873e-60
## 8: Happy__Neutral_start Happy Neutral_start 3.937832e-36
## 9: Sad__Neutral_start Sad Neutral_start 3.413310e-01
## q_bh_nine_conditions
## <num>
## 1: 1.499582e-77
## 2: 1.558959e-15
## 3: 3.184722e-35
## 4: 2.250680e-01
## 5: 2.785803e-22
## 6: 4.508176e-84
## 7: 6.725620e-60
## 8: 8.860122e-36
## 9: 3.413310e-01
## global_gate_note
## <char>
## 1: Interpret only after the global full-vs-reduced group-by-time test.
## 2: Interpret only after the global full-vs-reduced group-by-time test.
## 3: Interpret only after the global full-vs-reduced group-by-time test.
## 4: Interpret only after the global full-vs-reduced group-by-time test.
## 5: Interpret only after the global full-vs-reduced group-by-time test.
## 6: Interpret only after the global full-vs-reduced group-by-time test.
## 7: Interpret only after the global full-vs-reduced group-by-time test.
## 8: Interpret only after the global full-vs-reduced group-by-time test.
## 9: Interpret only after the global full-vs-reduced group-by-time test.9 Simultaneous confidence bands for MDD-minus-HC differences
The function below simulates from the fitted coefficient covariance matrix and constructs a simultaneous band for the difference in predicted emotional-gaze probability. It should be used only for conditions that pass the hierarchical inferential plan.
make_difference_band <- function(fit, dat, n_grid = 200L, n_sim = 2000L,
seed = 20260728L) {
stopifnot(length(unique(dat$condition)) == 1L)
time_grid <- seq(min(dat$time_center_ms), max(dat$time_center_ms), length.out = n_grid)
reference_subject <- levels(dat$subject)[1L]
make_new <- function(group_value) {
data.frame(
group = factor(group_value, levels = c("HC", "MDD")),
group_ordered = ordered(group_value, levels = c("HC", "MDD")),
time_center_ms = time_grid,
subject = factor(reference_subject, levels = levels(dat$subject))
)
}
hc_new <- make_new("HC")
mdd_new <- make_new("MDD")
X_hc <- predict(fit, newdata = hc_new, type = "lpmatrix")
X_mdd <- predict(fit, newdata = mdd_new, type = "lpmatrix")
beta <- coef(fit)
V <- vcov(fit, unconditional = TRUE)
set.seed(seed)
beta_draws <- MASS::mvrnorm(n_sim, mu = beta, Sigma = V)
eta_hc <- as.numeric(X_hc %*% beta)
eta_mdd <- as.numeric(X_mdd %*% beta)
point_difference <- plogis(eta_mdd) - plogis(eta_hc)
draw_hc <- plogis(X_hc %*% t(beta_draws))
draw_mdd <- plogis(X_mdd %*% t(beta_draws))
draw_diff <- draw_mdd - draw_hc
point_sd <- apply(draw_diff, 1L, sd)
standardized <- sweep(draw_diff, 1L, point_difference, "-")
standardized <- sweep(standardized, 1L, pmax(point_sd, 1e-10), "/")
critical <- unname(quantile(apply(abs(standardized), 2L, max), .95))
data.table(
time_ms = time_grid,
difference_mdd_minus_hc = point_difference,
simultaneous_low = point_difference - critical * point_sd,
simultaneous_high = point_difference + critical * point_sd
)
}
band_tables <- list()
for (i in seq_along(cell_fits)) {
condition_name <- cell_fits[[i]]$condition
if (is.finite(cell_tests[condition == condition_name, q_bh_nine_conditions]) &&
cell_tests[condition == condition_name, q_bh_nine_conditions] < .05) {
dat <- subject_bins[condition == condition_name]
band <- make_difference_band(cell_fits[[i]]$final, dat)
band[, condition := condition_name]
band_tables[[condition_name]] <- band
}
}
if (length(band_tables) > 0L) {
bands <- rbindlist(band_tables, use.names = TRUE, fill = TRUE)
fwrite(bands, file.path(output_dir, "simultaneous_group_difference_bands.csv"))
}10 Optional integration with the two rsFC axes
The multimodal models are intentionally disabled by default. Axis 1 is first related to fearful emotional-start and neutral-start trajectories, where the main theoretical question is start-dependent switching. Axis 2 is first related to the happy emotional-start trajectory. A significant axis-level result can be decomposed into its prespecified subfamilies using the exported score columns.
if (!file.exists(fc_score_file)) stop("Missing rsFC score file: ", fc_score_file)
fc_scores <- fread(fc_score_file)
fc_scores[, eye_subject_id := trimws(as.character(eye_subject_id))]
multimodal <- merge(
subject_bins,
fc_scores,
by.x = "subject",
by.y = "eye_subject_id",
all = FALSE
)
if (uniqueN(multimodal$subject) < 50L) {
warning("The overlapping multimodal sample contains fewer than 50 participants.")
}
multimodal[, axis1_z := as.numeric(scale(axis1_cortical_integration))]
multimodal[, axis2_z := as.numeric(scale(axis2_perceptual_motor_coupling))]
fit_axis_time_model <- function(dat, axis_variable, include_start = FALSE) {
dat <- copy(dat)
dat[, group_start := interaction(group, starting_position, drop = TRUE)]
if (include_start) {
reduced_text <- paste0(
"cbind(n_emotional, n_neutral) ~ group * starting_position * ", axis_variable,
" + s(time_center_ms, by = group_start, bs = 'cr', k = ", k_time, ")",
" + s(time_center_ms, subject, bs = 'fs', m = 1, k = ", k_subject_time, ")"
)
full_text <- paste0(
reduced_text,
" + ti(time_center_ms, ", axis_variable,
", by = group_start, k = c(", k_time, ", 5))"
)
} else {
reduced_text <- paste0(
"cbind(n_emotional, n_neutral) ~ group * ", axis_variable,
" + s(time_center_ms, by = group, bs = 'cr', k = ", k_time, ")",
" + s(time_center_ms, subject, bs = 'fs', m = 1, k = ", k_subject_time, ")"
)
full_text <- paste0(
reduced_text,
" + ti(time_center_ms, ", axis_variable,
", by = group, k = c(", k_time, ", 5))"
)
}
reduced_ml <- bam(
as.formula(reduced_text),
family = binomial(link = "logit"),
data = dat,
method = "ML",
discrete = FALSE,
nthreads = n_threads
)
full_ml <- bam(
as.formula(full_text),
family = binomial(link = "logit"),
data = dat,
method = "ML",
discrete = FALSE,
nthreads = n_threads
)
comparison <- as.data.frame(anova(reduced_ml, full_ml, test = "Chisq"))
p_col <- grep("Pr", names(comparison), value = TRUE)
p_value <- if (length(p_col) > 0L) tail(comparison[[p_col[1L]]], 1L) else NA_real_
final <- bam(
as.formula(full_text),
family = binomial(link = "logit"),
data = dat,
method = "fREML",
discrete = TRUE,
nthreads = n_threads
)
list(final = final, comparison = comparison, p_value = p_value,
reduced_formula = reduced_text, full_formula = full_text)
}
# Axis 1 primary multimodal condition: fearful E/N starts.
fear_switch <- multimodal[
expression_label == "Fear" &
starting_position %in% c("Emotional_start", "Neutral_start")
]
axis1_fear <- fit_axis_time_model(fear_switch, "axis1_z", include_start = TRUE)
saveRDS(axis1_fear$final, file.path(output_dir, "axis1_fear_switching_gamm.rds"))
writeLines(capture.output(summary(axis1_fear$final)),
file.path(output_dir, "axis1_fear_switching_gamm_summary.txt"))
# Axis 2 primary multimodal condition: happy emotional-location start.
happy_maintenance <- multimodal[
expression_label == "Happy" & starting_position == "Emotional_start"
]
axis2_happy <- fit_axis_time_model(happy_maintenance, "axis2_z", include_start = FALSE)
saveRDS(axis2_happy$final, file.path(output_dir, "axis2_happy_maintenance_gamm.rds"))
writeLines(capture.output(summary(axis2_happy$final)),
file.path(output_dir, "axis2_happy_maintenance_gamm_summary.txt"))
multimodal_tests <- data.table(
hypothesis = c(
"Axis1 x fearful start-dependent temporal trajectory",
"Axis2 x happy emotional-start temporal trajectory"
),
p_likelihood_ratio = c(axis1_fear$p_value, axis2_happy$p_value)
)
multimodal_tests[, p_holm_two_primary_tests := p.adjust(
p_likelihood_ratio, method = "holm"
)]
fwrite(multimodal_tests, file.path(output_dir, "multimodal_primary_tests.csv"))
write.csv(axis1_fear$comparison,
file.path(output_dir, "axis1_fear_model_comparison.csv"))
write.csv(axis2_happy$comparison,
file.path(output_dir, "axis2_happy_model_comparison.csv"))
fwrite(unique(multimodal[, .(subject, group)]),
file.path(output_dir, "multimodal_overlap_subjects.csv"))11 Reproducibility record
writeLines(c(
paste0("Created: ", Sys.time()),
paste0("Input file: ", input_file),
paste0("AOI coding: neutral = ", neutral_aoi_code,
"; emotional = ", emotional_aoi_code),
paste0("Bin width: ", bin_width_ms, " ms"),
paste0("Time range: ", analysis_start_ms, " to ", analysis_end_ms, " ms"),
paste0("Participants: ", uniqueN(subject_bins$subject)),
paste0("Multimodal models run: ", run_multimodal_models),
"",
capture.output(sessionInfo())
), file.path(output_dir, "temporal_gamm_analysis_log.txt"))