Temporal Emotional Gaze in MDD: 16-GB Quasi-Binomial AR(1) GAMM
Tak Kwan Lam
2026-08-13
1 Purpose
This version is designed for a computer with about 16 GB RAM. It
preserves trial identity, participant-specific nonlinear time courses,
condition-specific group-by-time effects, and multiplicity correction,
but it deliberately avoids an observation-level random effect with one
level per row. Extra-binomial variation is handled with a quasi-binomial
variance model, and residual serial dependence within each trial is
handled by the AR(1) working-residual correction available in
mgcv::bam(..., discrete = TRUE).
The inferential plan is:
- A global full-versus-reduced group-by-time comparison across all conditions using an approximate F test under the quasi-binomial model.
- One condition-specific GAMM per expression-by-start-code condition, fitted sequentially to keep memory bounded.
- Benjamini-Hochberg correction across the complete planned condition family.
- Simultaneous MDD-minus-HC difference bands only for conditions that survive the multiplicity correction.
2 Settings
data_folder <- "C:/Users/psyuser/Desktop/fmri eye/data"
input_file <- file.path(data_folder, "all_subjects_full.csv")
output_dir <- file.path(data_folder, "2_eye_temporal_gamm_outputs")
# Keep raw AOI and start codes until their semantic mapping is fully verified.
neutral_aoi_code <- 1L
emotional_aoi_code <- 2L
expression_map <- c(F = "Fear", H = "Happy", S = "Sad")
sample_rate_hz <- 500
# The manuscript audit suggests `count` may be a sample index rather than ms.
# "auto" converts a median step near 1 into 2 ms/sample at 500 Hz, and treats a
# median step near 2 as already being milliseconds. For the final paper, verify
# this against the task program/codebook and then set explicitly to
# "sample_index" or "milliseconds".
count_time_unit <- "auto" # "auto", "sample_index", or "milliseconds"
analysis_start_ms <- 0
analysis_end_ms <- 2000
bin_width_ms <- 20 # 40 bins across 2 s; use 100 only if still needed
minimum_valid_samples_per_bin <- 3L
# Conservative basis sizes for a 16-GB machine. Increase only if k-check says
# the basis is too restrictive.
k_time <- 10L
k_subject_time <- 4L
# Fewer threads usually reduce peak RAM on a laptop/desktop.
n_threads <- min(2L, max(1L, parallel::detectCores(logical = TRUE) - 1L))
bam_chunk_size <- 5000L
bam_gc_level <- 2L
bam_samfrac <- 0.10
# AR(1) estimate is clipped to avoid unstable extreme values.
rho_min <- -0.80
rho_max <- 0.80
# Lower Monte-Carlo precision for confidence bands saves time and some RAM.
# Increase to 2000-5000 for the final manuscript if desired.
confidence_band_simulations <- 1000L
confidence_band_grid <- 120L
# Do not silently exclude a condition. Failed fits are reported and receive no
# inferential claim; BH correction still uses the full planned family size.
excluded_conditions <- character(0)
verbose_progress <- TRUE
dir.create(output_dir, recursive = TRUE, showWarnings = FALSE)3 Packages and helpers
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 required package(s): ", paste(missing_packages, collapse = ", "))
}
library(data.table)## Warning: package 'data.table' was built under R version 4.3.3library(mgcv)## Loading required package: nlme## This is mgcv 1.9-0. For overview type 'help("mgcv-package")'.library(ggplot2)
log_progress <- function(...) {
if (!isTRUE(verbose_progress)) return(invisible(NULL))
cat(sprintf("[%s] %s\n", format(Sys.time(), "%Y-%m-%d %H:%M:%S"),
paste0(..., collapse = "")))
}
time_step <- function(label, expr) {
expr_sub <- substitute(expr)
log_progress("START ", label)
t0 <- proc.time()[["elapsed"]]
out <- eval.parent(expr_sub)
elapsed <- proc.time()[["elapsed"]] - t0
log_progress("END ", label, " (", sprintf("%.2f", elapsed), " sec)")
out
}
safe_name <- function(x) gsub("[^A-Za-z0-9._-]", "_", x)
estimate_rho <- function(resid_vec, series_id) {
tmp <- data.table(series_id = series_id, r = as.numeric(resid_vec))
tmp[, lag_r := shift(r), by = series_id]
val <- suppressWarnings(cor(tmp$r, tmp$lag_r, use = "complete.obs"))
if (!is.finite(val)) val <- 0
max(rho_min, min(rho_max, val))
}
extract_difference_smooth_p <- function(fit) {
st <- summary(fit)$s.table
if (is.null(st) || nrow(st) == 0L) return(NA_real_)
idx <- grep("group_ordered", rownames(st), fixed = TRUE)
if (length(idx) == 0L) return(NA_real_)
pcol <- grep("p", colnames(st), ignore.case = TRUE)
if (length(pcol) == 0L) return(NA_real_)
as.numeric(st[idx[1L], pcol[length(pcol)]])
}4 Load, audit, and reduce raw samples without duplicating the full data
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 = ", "))
}
# Rename instead of creating duplicate copies of the original columns.
setnames(
raw,
required_columns,
c("subject", "group_raw", "start_code", "expression_raw",
"trial", "count_raw", "AOI_code")
)
raw[, `:=`(
subject = trimws(as.character(subject)),
group_raw = toupper(trimws(as.character(group_raw))),
start_code = toupper(trimws(as.character(start_code))),
expression_raw = toupper(trimws(as.character(expression_raw))),
trial = as.integer(trial),
count_raw = as.numeric(count_raw),
AOI_code = suppressWarnings(as.integer(AOI_code))
)]
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])]
if (anyNA(raw$group)) stop("At least one sub_type value could not be mapped.")
if (anyNA(raw$expression_label)) {
stop("Unmapped expression code(s): ",
paste(unique(raw[is.na(expression_label), expression_raw]), collapse = ", "))
}
if (!all(raw$start_code %in% c("C", "E", "N"))) {
stop("Unexpected trial_type code(s): ",
paste(setdiff(unique(raw$start_code), c("C", "E", "N")), collapse = ", "))
}
observed_aoi_codes <- sort(unique(na.omit(raw$AOI_code)))
if (!all(c(neutral_aoi_code, emotional_aoi_code) %in% observed_aoi_codes)) {
stop("Configured AOI codes were not both observed. Observed: ",
paste(observed_aoi_codes, collapse = ", "))
}
# Time-unit audit before aggregation.
time_step_audit <- raw[
is.finite(count_raw),
.(median_step = as.numeric(median(diff(sort(unique(count_raw))), na.rm = TRUE))),
by = .(subject, trial)
]
overall_median_step <- median(time_step_audit$median_step, na.rm = TRUE)
fwrite(time_step_audit, file.path(output_dir, "time_step_audit_by_trial.csv"))
if (count_time_unit == "sample_index") {
time_multiplier <- 1000 / sample_rate_hz
} else if (count_time_unit == "milliseconds") {
time_multiplier <- 1
} else if (count_time_unit == "auto") {
if (is.finite(overall_median_step) && abs(overall_median_step - 1) <= 0.25) {
time_multiplier <- 1000 / sample_rate_hz
warning("count behaves like a sample index; converting 1 unit to ",
time_multiplier, " ms. Verify against the task program before final reporting.")
} else if (is.finite(overall_median_step) &&
abs(overall_median_step - (1000 / sample_rate_hz)) <= 0.5) {
time_multiplier <- 1
} else {
stop("Could not safely infer time units from median step = ",
round(overall_median_step, 3),
". Set count_time_unit explicitly after checking the codebook.")
}
} else {
stop("count_time_unit must be 'auto', 'sample_index', or 'milliseconds'.")
}## Warning: count behaves like a sample index; converting 1 unit to 2 ms. Verify
## against the task program before final reporting.raw[, time_ms := count_raw * time_multiplier]
coding_audit <- data.table(
field = c("AOI code 1", "AOI code 2", "start codes", "time multiplier"),
configured_value = c(
neutral_aoi_code, emotional_aoi_code, "C/E/N kept as raw codes", time_multiplier
)
)
fwrite(coding_audit, file.path(output_dir, "coding_audit.csv"))
print(coding_audit)## field configured_value
## <char> <char>
## 1: AOI code 1 1
## 2: AOI code 2 2
## 3: start codes C/E/N kept as raw codes
## 4: time multiplier 25 Trial-level non-overlapping bins
# Aggregate directly from `raw`; do not create a second full-size `valid` copy.
trial_bins <- raw[
is.finite(time_ms) &
time_ms >= analysis_start_ms & time_ms < analysis_end_ms &
AOI_code %in% c(neutral_aoi_code, emotional_aoi_code),
.(
n_code2 = sum(AOI_code == emotional_aoi_code),
n_code1 = sum(AOI_code == neutral_aoi_code),
n_valid = .N
),
by = .(
subject, group, expression_label, start_code, trial,
time_bin = floor((time_ms - analysis_start_ms) / bin_width_ms)
)
]
trial_bins <- trial_bins[n_valid >= minimum_valid_samples_per_bin]
trial_bins[, time_center_ms :=
analysis_start_ms + (time_bin + 0.5) * bin_width_ms]
trial_bins[, probability_code2 := n_code2 / (n_code2 + n_code1)]
# The raw sample table is no longer needed. Releasing it here is a major RAM
# saving compared with retaining raw + valid + trial_bins simultaneously.
rm(raw)
gc(verbose = FALSE)## used (Mb) gc trigger (Mb) max used (Mb)
## Ncells 2032297 108.6 4025007 215.0 4025007 215.0
## Vcells 14730713 112.4 305963163 2334.4 316627955 2415.7trial_bins[, `:=`(
subject = factor(subject),
group = factor(group, levels = c("HC", "MDD")),
expression_label = factor(expression_label, levels = c("Fear", "Happy", "Sad")),
start_code = factor(start_code, levels = c("C", "E", "N"))
)]
trial_bins[, condition := interaction(
expression_label, start_code, drop = TRUE, sep = "__"
)]
if (length(excluded_conditions) > 0L) {
trial_bins <- trial_bins[!as.character(condition) %in% excluded_conditions]
trial_bins[, condition := droplevels(condition)]
}
# A unique trial series is needed for the AR(1) reset flag.
trial_bins[, trial_id := interaction(
subject, expression_label, start_code, trial, drop = TRUE, sep = "__"
)]
setorder(trial_bins, trial_id, time_bin)
trial_bins[, AR_start := !duplicated(trial_id)]
fwrite(
trial_bins[, .(subject, group, expression_label, start_code, trial,
time_bin, time_center_ms, n_code2, n_code1, n_valid)],
file.path(output_dir, "trial_time_bins_for_gamm.csv")
)
sample_audit <- trial_bins[, .(
participants = uniqueN(subject),
trials = uniqueN(trial_id),
trial_bins = .N,
median_valid_samples = as.numeric(median(n_valid))
), by = .(group, expression_label, start_code)]
fwrite(sample_audit, file.path(output_dir, "temporal_sample_audit.csv"))
print(sample_audit)## group expression_label start_code participants trials trial_bins
## <fctr> <fctr> <fctr> <int> <int> <int>
## 1: HC Fear C 47 729 54539
## 2: MDD Fear C 63 945 64068
## 3: HC Happy C 47 728 52125
## 4: MDD Happy C 63 942 62575
## 5: HC Sad C 47 717 51918
## 6: MDD Sad C 63 953 62901
## 7: HC Fear E 47 752 66830
## 8: MDD Fear E 63 1008 82217
## 9: HC Happy E 47 752 66416
## 10: MDD Happy E 63 1007 82312
## 11: HC Sad E 47 752 65333
## 12: MDD Sad E 63 1008 80942
## 13: HC Fear N 47 752 66441
## 14: MDD Fear N 63 1008 82871
## 15: HC Happy N 47 752 66823
## 16: MDD Happy N 63 1008 81702
## 17: HC Sad N 47 752 65811
## 18: MDD Sad N 63 1007 81534
## median_valid_samples
## <num>
## 1: 10
## 2: 10
## 3: 10
## 4: 10
## 5: 10
## 6: 10
## 7: 10
## 8: 10
## 9: 10
## 10: 10
## 11: 10
## 12: 10
## 13: 10
## 14: 10
## 15: 10
## 16: 10
## 17: 10
## 18: 106 Descriptive trajectories
# Average within participant first so participants, not trials, define the
# descriptive group mean.
subject_plot <- trial_bins[, .(
n_code2 = sum(n_code2),
n_code1 = sum(n_code1)
), by = .(subject, group, expression_label, start_code, time_center_ms)]
subject_plot[, probability_code2 := n_code2 / (n_code2 + n_code1)]
plot_data <- subject_plot[, .(
mean_probability = mean(probability_code2, na.rm = TRUE),
se_probability = sd(probability_code2, na.rm = TRUE) / sqrt(.N)
), by = .(group, expression_label, start_code, 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 ~ start_code) +
geom_hline(yintercept = 0.5, linetype = 2) +
theme_classic() +
labs(
x = "Time after stimulus onset (ms)",
y = "Probability assigned to AOI code 2",
colour = "Group", fill = "Group"
)
print(p)
ggsave(file.path(output_dir, "descriptive_temporal_trajectories.png"),
p, width = 12, height = 9, dpi = 300)
rm(subject_plot, plot_data)
gc(verbose = FALSE)## used (Mb) gc trigger (Mb) max used (Mb)
## Ncells 2263324 120.9 4025007 215.0 4025007 215.0
## Vcells 14717094 112.3 244770531 1867.5 316627955 2415.77 Global quasi-binomial AR(1) GAMM
condition_levels <- levels(trial_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]
trial_bins[[nm]] <- as.numeric(
trial_bins$group == "MDD" & trial_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_code2, n_code1) ~", paste(base_terms, collapse = " + ")
))
full_formula <- as.formula(paste(
"cbind(n_code2, n_code1) ~",
paste(c(base_terms, difference_terms), collapse = " + ")
))
# Pilot fit only estimates a reasonable lag-1 working-residual correlation.
global_pilot <- time_step(
"global pilot fit for rho",
bam(
reduced_formula,
family = quasibinomial(link = "logit"),
data = trial_bins,
method = "fREML",
discrete = TRUE,
rho = 0,
AR.start = trial_bins$AR_start,
nthreads = n_threads,
chunk.size = bam_chunk_size,
gc.level = bam_gc_level,
samfrac = bam_samfrac
)
)## [2026-08-14 11:29:25] START global pilot fit for rho
## [2026-08-14 12:00:54] END global pilot fit for rho (1888.84 sec)global_rho <- estimate_rho(
residuals(global_pilot, type = "pearson"),
trial_bins$trial_id
)
log_progress("Estimated AR(1) rho = ", sprintf("%.3f", global_rho))## [2026-08-14 12:00:54] Estimated AR(1) rho = 0.800rm(global_pilot)
gc(verbose = FALSE)## used (Mb) gc trigger (Mb) max used (Mb)
## Ncells 2326263 124.3 4025007 215 4025007 215.0
## Vcells 25942974 198.0 80206408 612 316627955 2415.7# Reduced and full models use the same data, smoothing-selection method, and rho.
global_reduced <- time_step(
"global reduced quasi-AR1 fit",
bam(
reduced_formula,
family = quasibinomial(link = "logit"),
data = trial_bins,
method = "fREML",
discrete = TRUE,
rho = global_rho,
AR.start = trial_bins$AR_start,
nthreads = n_threads,
chunk.size = bam_chunk_size,
gc.level = bam_gc_level,
samfrac = bam_samfrac
)
)## [2026-08-14 12:00:55] START global reduced quasi-AR1 fit
## [2026-08-14 12:50:56] END global reduced quasi-AR1 fit (3000.91 sec)global_full <- time_step(
"global full quasi-AR1 fit",
bam(
full_formula,
family = quasibinomial(link = "logit"),
data = trial_bins,
method = "fREML",
discrete = TRUE,
rho = global_rho,
AR.start = trial_bins$AR_start,
nthreads = n_threads,
chunk.size = bam_chunk_size,
gc.level = bam_gc_level,
samfrac = bam_samfrac
)
)## [2026-08-14 12:50:56] START global full quasi-AR1 fit
## [2026-08-14 13:43:17] END global full quasi-AR1 fit (3141.27 sec)# With a quasi family there is no ordinary likelihood-ratio test. This F test
# is approximate and should be reported as such.
global_comparison <- as.data.frame(
anova(global_reduced, global_full, test = "F")
)
print(global_comparison)## Resid. Df Resid. Dev Df Deviance F Pr(>F)
## 1 1236833 14556026 NA NA NA NA
## 2 1236778 14537315 54.98438 18710.35 105.0139 0write.csv(global_comparison,
file.path(output_dir, "global_group_by_time_approximate_F_test.csv"),
row.names = FALSE)
saveRDS(global_full, file.path(output_dir, "global_temporal_quasi_AR1_gamm.rds"))
writeLines(capture.output(summary(global_full)),
file.path(output_dir, "global_temporal_quasi_AR1_summary.txt"))
writeLines(capture.output(k.check(global_full)),
file.path(output_dir, "global_temporal_quasi_AR1_basis_check.txt"))
quasi_scale <- summary(global_full)$scale
fwrite(
data.table(
metric = c("estimated quasi dispersion", "AR1 rho", "rows", "conditions"),
value = c(quasi_scale, global_rho, nrow(trial_bins), length(condition_levels))
),
file.path(output_dir, "global_model_diagnostics.csv")
)
rm(global_reduced)
gc(verbose = FALSE)## used (Mb) gc trigger (Mb) max used (Mb)
## Ncells 2349356 125.5 4025007 215.0 4025007 215.0
## Vcells 54451233 415.5 160593235 1225.3 316627955 2415.78 Condition-specific tests, one model at a time
condition_result_list <- vector("list", length(condition_levels))
model_paths <- character(length(condition_levels))
for (i in seq_along(condition_levels)) {
condition_value <- condition_levels[i]
log_progress("Condition ", i, "/", length(condition_levels), ": ", condition_value)
dat <- copy(trial_bins[condition == condition_value])
dat[, `:=`(
subject = droplevels(subject),
condition = droplevels(condition),
trial_id = droplevels(trial_id),
group_ordered = ordered(group, levels = c("HC", "MDD"))
)]
setorder(dat, trial_id, time_bin)
dat[, AR_start := !duplicated(trial_id)]
formula_condition <- cbind(n_code2, n_code1) ~ 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)
fit_obj <- tryCatch(
time_step(
paste0("condition fit: ", condition_value),
bam(
formula_condition,
family = quasibinomial(link = "logit"),
data = dat,
method = "fREML",
discrete = TRUE,
rho = global_rho,
AR.start = dat$AR_start,
nthreads = n_threads,
chunk.size = bam_chunk_size,
gc.level = bam_gc_level,
samfrac = bam_samfrac
)
),
error = function(e) e
)
parts <- tstrsplit(condition_value, "__", fixed = TRUE)
if (inherits(fit_obj, "error")) {
condition_result_list[[i]] <- data.table(
condition = condition_value,
expression = parts[[1L]],
start_code = parts[[2L]],
p_group_by_time = NA_real_,
convergence_status = paste("FAILED:", conditionMessage(fit_obj))
)
model_paths[i] <- NA_character_
} else {
p_value <- extract_difference_smooth_p(fit_obj)
condition_result_list[[i]] <- data.table(
condition = condition_value,
expression = parts[[1L]],
start_code = parts[[2L]],
p_group_by_time = p_value,
convergence_status = "fit completed"
)
model_paths[i] <- file.path(
output_dir,
paste0("condition_", safe_name(condition_value), "_quasi_AR1.rds")
)
saveRDS(fit_obj, model_paths[i])
writeLines(
capture.output(summary(fit_obj)),
file.path(output_dir,
paste0("condition_", safe_name(condition_value), "_summary.txt"))
)
}
rm(dat, fit_obj)
gc(verbose = FALSE)
}## [2026-08-14 13:43:41] Condition 1/9: Fear__C
## [2026-08-14 13:43:41] START condition fit: Fear__C## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has repeated
## 1-d smooths of same variable.## [2026-08-14 13:45:33] END condition fit: Fear__C (111.98 sec)
## [2026-08-14 13:45:36] Condition 2/9: Happy__C
## [2026-08-14 13:45:36] START condition fit: Happy__C## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has repeated
## 1-d smooths of same variable.## [2026-08-14 13:47:24] END condition fit: Happy__C (107.35 sec)
## [2026-08-14 13:47:27] Condition 3/9: Sad__C
## [2026-08-14 13:47:27] START condition fit: Sad__C## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has repeated
## 1-d smooths of same variable.## [2026-08-14 13:49:13] END condition fit: Sad__C (106.91 sec)
## [2026-08-14 13:49:17] Condition 4/9: Fear__E
## [2026-08-14 13:49:17] START condition fit: Fear__E## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has repeated
## 1-d smooths of same variable.## [2026-08-14 13:53:00] END condition fit: Fear__E (223.64 sec)
## [2026-08-14 13:53:04] Condition 5/9: Happy__E
## [2026-08-14 13:53:04] START condition fit: Happy__E## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has repeated
## 1-d smooths of same variable.## [2026-08-14 13:56:09] END condition fit: Happy__E (184.93 sec)
## [2026-08-14 13:56:12] Condition 6/9: Sad__E
## [2026-08-14 13:56:12] START condition fit: Sad__E## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has repeated
## 1-d smooths of same variable.## [2026-08-14 13:58:42] END condition fit: Sad__E (150.70 sec)
## [2026-08-14 13:58:46] Condition 7/9: Fear__N
## [2026-08-14 13:58:46] START condition fit: Fear__N## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has repeated
## 1-d smooths of same variable.## [2026-08-14 14:01:21] END condition fit: Fear__N (155.53 sec)
## [2026-08-14 14:01:24] Condition 8/9: Happy__N
## [2026-08-14 14:01:25] START condition fit: Happy__N## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has repeated
## 1-d smooths of same variable.## [2026-08-14 14:04:11] END condition fit: Happy__N (166.33 sec)
## [2026-08-14 14:04:14] Condition 9/9: Sad__N
## [2026-08-14 14:04:14] START condition fit: Sad__N## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has repeated
## 1-d smooths of same variable.## [2026-08-14 14:10:45] END condition fit: Sad__N (390.98 sec)condition_tests <- rbindlist(condition_result_list, use.names = TRUE, fill = TRUE)
# Conservatively keep the planned family size even if a condition fit fails:
# an unestimable p value is treated as 1 for the adjustment, then displayed NA.
p_for_adjust <- condition_tests$p_group_by_time
missing_p <- !is.finite(p_for_adjust)
p_for_adjust[missing_p] <- 1
condition_tests[, q_bh_all_planned_conditions := p.adjust(p_for_adjust, method = "BH")]
condition_tests[missing_p, q_bh_all_planned_conditions := NA_real_]
fwrite(condition_tests,
file.path(output_dir, "condition_specific_group_by_time_tests.csv"))
print(condition_tests)## condition expression start_code p_group_by_time convergence_status
## <char> <char> <char> <num> <char>
## 1: Fear__C Fear C 9.085307e-05 fit completed
## 2: Happy__C Happy C 7.333480e-01 fit completed
## 3: Sad__C Sad C 0.000000e+00 fit completed
## 4: Fear__E Fear E 3.533447e-01 fit completed
## 5: Happy__E Happy E 6.247985e-01 fit completed
## 6: Sad__E Sad E 2.569658e-05 fit completed
## 7: Fear__N Fear N 5.395160e-01 fit completed
## 8: Happy__N Happy N 8.537635e-01 fit completed
## 9: Sad__N Sad N 9.999477e-01 fit completed
## q_bh_all_planned_conditions
## <num>
## 1: 0.0002725592
## 2: 0.9428760415
## 3: 0.0000000000
## 4: 0.7950256243
## 5: 0.9371977450
## 6: 0.0001156346
## 7: 0.9371977450
## 8: 0.9604839840
## 9: 0.99994773949 Simultaneous difference bands
make_difference_band <- function(fit, dat,
n_grid = confidence_band_grid,
n_sim = confidence_band_simulations,
seed = 20260813L) {
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")
# Exclude the participant-specific factor smooth so the contrast is a
# population-level MDD-minus-HC trajectory rather than a reference-subject
# trajectory.
smooth_labels <- vapply(fit$smooth, function(x) x$label, character(1L))
subject_smooth_label <- grep("subject", smooth_labels, value = TRUE, fixed = TRUE)
X_hc <- predict(fit, newdata = hc_new, type = "lpmatrix",
exclude = subject_smooth_label)
X_mdd <- predict(fit, newdata = mdd_new, type = "lpmatrix",
exclude = subject_smooth_label)
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(condition_levels)) {
condition_value <- condition_levels[i]
q_value <- condition_tests[
condition == condition_value, q_bh_all_planned_conditions
]
if (length(q_value) == 1L && is.finite(q_value) && q_value < .05 &&
is.character(model_paths[i]) && !is.na(model_paths[i]) &&
file.exists(model_paths[i])) {
dat <- copy(trial_bins[condition == condition_value])
dat[, `:=`(
subject = droplevels(subject),
group_ordered = ordered(group, levels = c("HC", "MDD"))
)]
fit_obj <- readRDS(model_paths[i])
band <- make_difference_band(fit_obj, dat)
band[, condition := condition_value]
band_tables[[condition_value]] <- band
rm(dat, fit_obj, band)
gc(verbose = FALSE)
}
}
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 Runtime / reproducibility metadata
fwrite(
data.table(
setting = c(
"bin_width_ms", "k_time", "k_subject_time", "n_threads",
"bam_chunk_size", "bam_gc_level", "bam_samfrac",
"confidence_band_simulations", "count_time_multiplier", "AR1_rho"
),
value = as.character(c(
bin_width_ms, k_time, k_subject_time, n_threads,
bam_chunk_size, bam_gc_level, bam_samfrac,
confidence_band_simulations, time_multiplier, global_rho
))
),
file.path(output_dir, "runtime_and_model_settings.csv")
)
sessionInfo()## R version 4.3.2 (2023-10-31 ucrt)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 11 x64 (build 26200)
##
## Matrix products: default
##
##
## locale:
## [1] LC_COLLATE=English_Hong Kong SAR.utf8
## [2] LC_CTYPE=English_Hong Kong SAR.utf8
## [3] LC_MONETARY=English_Hong Kong SAR.utf8
## [4] LC_NUMERIC=C
## [5] LC_TIME=English_Hong Kong SAR.utf8
##
## time zone: Asia/Hong_Kong
## tzcode source: internal
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] ggplot2_3.5.1 mgcv_1.9-0 nlme_3.1-163 data.table_1.15.4
##
## loaded via a namespace (and not attached):
## [1] Matrix_1.6-1.1 gtable_0.3.5 jsonlite_1.8.8 dplyr_1.1.4
## [5] compiler_4.3.2 tidyselect_1.2.1 parallel_4.3.2 jquerylib_0.1.4
## [9] textshaping_0.4.0 systemfonts_1.1.0 splines_4.3.2 scales_1.3.0
## [13] yaml_2.3.8 fastmap_1.2.0 lattice_0.21-9 R6_2.5.1
## [17] labeling_0.4.3 generics_0.1.3 knitr_1.49 MASS_7.3-60
## [21] tibble_3.2.1 munsell_0.5.1 bslib_0.7.0 pillar_1.9.0
## [25] rlang_1.1.4 utf8_1.2.4 cachem_1.1.0 xfun_0.49
## [29] sass_0.4.9 cli_3.6.2 withr_3.0.0 magrittr_2.0.3
## [33] digest_0.6.35 grid_4.3.2 rstudioapi_0.16.0 lifecycle_1.0.4
## [37] vctrs_0.6.5 evaluate_0.24.0 glue_1.7.0 farver_2.1.2
## [41] ragg_1.3.2 fansi_1.0.6 colorspace_2.1-0 rmarkdown_2.29
## [45] tools_4.3.2 pkgconfig_2.0.3 htmltools_0.5.8.111 Interpretation boundary
This lower-memory model is intended to answer the same temporal questions while avoiding the one-random-effect-per-row construction. The global quasi-model F test is approximate rather than a likelihood-ratio test. Condition-specific inference is based on the MDD-minus-HC difference smooth in each full model and BH correction across the planned condition family. Before directional claims about emotional versus neutral gaze are made, verify the AOI and C/E/N code mappings and the source time unit against the original task documentation.