Cross-Modal MVPA: Resting-State Connectivity Prediction of Temporal Gaze in MDD
Tak Kwan Lam
2026-08-16
1 Purpose and inferential boundary
This report performs a participant-level cross-modal analysis. It asks whether resting-state functional-connectivity patterns predict individual temporal gaze trajectories in the participants who have both modalities.
The workflow deliberately consumes the validated outputs of these upstream files rather than sourcing them, because the upstream files contain executable top-level analysis code:
1_resting_state_fc_prepare_validate_updated(5).R3_resting_state_fc_two_axis_targeted_updated(5).RA_new_2_eye_temporal_two_axis_fixed(2).Rmd
The analysis has two complementary parts:
- A theory-guided multivariate association test asks whether the two manuscript-defined connectivity axes account for additional variation in a compact representation of the gaze trajectories after diagnosis and covariates are included.
- A nested-cross-validated multi-response elastic-net model asks whether either the two-axis pattern or the complete 22-edge targeted pattern predicts gaze trajectories in held-out participants.
The primary predictive interpretation should be attached to the
incremental_beyond_group analysis. The
total_without_group analysis is useful but can succeed
simply because both fMRI and eye tracking differ between MDD and HC.
Neither analysis establishes causality, mediation, or a mechanism by
which resting-state connectivity causes gaze behavior.
2 Settings
data_dir <- params$data_dir
# Outputs created by the two resting-state scripts.
prepared_rds <- file.path(
data_dir, "rsfc_prepared", "rsfc_prepared_data.rds"
)
axis_scores_file <- file.path(
data_dir, "rsfc_two_axis", "two_axis_subject_scores_for_temporal_integration.csv"
)
manifest_file <- file.path(
data_dir, "rsfc_two_axis", "two_axis_edge_manifest_verified.csv"
)
# Output created by the temporal eye-tracking R Markdown file.
eye_bins_file <- file.path(
data_dir, "2_eye_temporal_gamm_outputs", "trial_time_bins_for_gamm.csv"
)
eye_condition_tests_file <- file.path(
data_dir, "2_eye_temporal_gamm_outputs",
"condition_specific_group_by_time_tests.csv"
)
output_dir <- file.path(data_dir, "cross_modal_fmri_eye_mvpa")
dir.create(output_dir, recursive = TRUE, showWarnings = FALSE)
# Choose the eye-trajectory family. The three surviving conditions are the
# manuscript-focused analysis, but they were identified in the same eye sample
# and therefore remain secondary/exploratory. Re-run with "all_nine" as a
# selection-bias sensitivity analysis.
eye_condition_set <- "surviving_three" # "surviving_three" or "all_nine"
surviving_eye_conditions <- c("Fear__C", "Sad__C", "Sad__E")
all_nine_eye_conditions <- as.vector(outer(
c("Fear", "Happy", "Sad"),
c("C", "E", "N"),
paste,
sep = "__"
))
primary_eye_conditions <- switch(
eye_condition_set,
surviving_three = surviving_eye_conditions,
all_nine = all_nine_eye_conditions,
stop("eye_condition_set must be 'surviving_three' or 'all_nine'.")
)
# Each subject-condition trajectory is represented by a fixed cubic B-spline
# basis fitted to the empirical log odds of emotional-face versus neutral-face
# gaze. This is a compact phenotype, not a replacement for the trial-level GAMM.
eye_basis_df <- 8L
eye_basis_ridge <- 1e-4
minimum_eye_bins_per_condition <- 60L
# Predictive feature sets. The axes are the theory-driven primary pattern; the
# 22-edge model is an exploratory targeted MVPA.
feature_sets_to_run <- c("axes2", "targeted_edges22")
analysis_modes_to_run <- c("total_without_group", "incremental_beyond_group")
# Repeated nested cross-validation. All preprocessing, nuisance regression,
# axis construction, and model tuning are refitted within the appropriate fold.
outer_folds <- 5L
outer_repeats <- 10L
inner_folds <- 5L
alpha_grid <- c(0, 0.25, 0.50, 0.75, 1)
lambda_ratio_grid <- exp(seq(log(1), log(0.01), length.out = 40L))
lambda_selection_rule <- "one_se" # "one_se" or "minimum"
glmnet_nlambda <- 60L
glmnet_lambda_min_ratio <- 0.01
# Motion is strongly preferred. If it is unavailable for enough participants,
# the analysis continues but is explicitly flagged as provisional.
minimum_motion_fraction <- 0.80
require_motion <- FALSE
include_icv_sensitivity <- FALSE
# Permutation and bootstrap counts. The prediction permutation holds the
# out-of-fold predictions fixed; it is not a full retraining permutation.
axis_permutations <- 2000L
prediction_permutations <- 5000L
prediction_bootstraps <- 2000L
# Set TRUE only after documenting that the 22-edge manifest was specified
# independently of the group effects in the current sample.
targets_independently_prespecified <- FALSE
random_seed <- 20260816L
verbose_progress <- TRUE3 Packages
required_packages <- c("data.table", "ggplot2", "glmnet", "knitr")
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)
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 = "")
))
}
stop_if_missing <- function(paths) {
missing <- paths[!file.exists(paths)]
if (length(missing) > 0L) {
stop("Missing required file(s): ", paste(missing, collapse = ", "))
}
}
required_files <- c(
prepared_rds,
axis_scores_file,
manifest_file,
eye_bins_file
)
stop_if_missing(required_files)
if (!isTRUE(targets_independently_prespecified)) {
warning(
"The 22-edge model must be described as exploratory until independent ",
"prespecification of the target manifest is documented."
)
}## Warning: The 22-edge model must be described as exploratory until independent
## prespecification of the target manifest is documented.4 Helper functions
normalize_subject_id <- function(x) {
missing <- is.na(x) | trimws(as.character(x)) == ""
out <- toupper(trimws(as.character(x)))
out <- sub("^SUB[-_ ]*", "", out)
out <- gsub("[^A-Z0-9]", "", out)
out[missing | out == ""] <- NA_character_
out
}
check_unique_key <- function(data, key_col, label) {
key <- data[[key_col]]
if (anyNA(key) || any(key == "")) {
stop(label, " has a missing or empty normalized subject ID.")
}
duplicate_rows <- duplicated(key) | duplicated(key, fromLast = TRUE)
dup <- data[duplicate_rows]
if (nrow(dup) > 0L) {
stop(
label, " has duplicated normalized IDs. First duplicated key: ",
dup[[key_col]][1L]
)
}
}
make_eye_basis_function <- function(time_min, time_max, basis_df) {
force(time_min)
force(time_max)
force(basis_df)
function(time_ms) {
if (!all(is.finite(time_ms))) stop("Non-finite eye-tracking time value.")
time_01 <- (time_ms - time_min) / (time_max - time_min)
splines::bs(
time_01,
df = basis_df,
degree = 3,
intercept = TRUE,
Boundary.knots = c(0, 1)
)
}
}
fit_subject_eye_basis <- function(dat, basis_fun, ridge) {
setorder(dat, time_center_ms)
B <- basis_fun(dat$time_center_ms)
y <- log((dat$n_code2 + 0.5) / (dat$n_code1 + 0.5))
w <- dat$n_code2 + dat$n_code1
gram <- crossprod(B, B * w)
rhs <- crossprod(B, y * w)
beta <- as.numeric(solve(
gram + diag(ridge, ncol(B)),
rhs
))
fitted_logit <- as.numeric(B %*% beta)
weighted_rmse <- sqrt(weighted.mean((y - fitted_logit)^2, w = w))
list(
beta = beta,
n_bins = nrow(dat),
total_valid_samples = sum(w),
weighted_logit_rmse = weighted_rmse,
first_time_ms = min(dat$time_center_ms),
last_time_ms = max(dat$time_center_ms)
)
}axis1_subfamilies <- c(
"posterior_medial",
"cingulo_perisylvian",
"insular_opercular"
)
axis2_subfamilies <- c(
"reward_visual",
"brainstem_sensorimotor"
)
all_subfamilies <- c(axis1_subfamilies, axis2_subfamilies)
standardize_from_training <- function(train_matrix, test_matrix) {
mu <- colMeans(train_matrix)
sig <- apply(train_matrix, 2L, sd)
if (any(!is.finite(sig) | sig <= 0)) {
bad <- colnames(train_matrix)[!is.finite(sig) | sig <= 0]
stop("Zero-variance or invalid training feature(s): ", paste(bad, collapse = ", "))
}
list(
train = sweep(sweep(train_matrix, 2L, mu, "-"), 2L, sig, "/"),
test = sweep(sweep(test_matrix, 2L, mu, "-"), 2L, sig, "/"),
mean = mu,
sd = sig
)
}
build_fold_features <- function(
edge_train,
edge_test,
feature_spec,
manifest
) {
edge_train <- as.matrix(edge_train)
edge_test <- as.matrix(edge_test)
storage.mode(edge_train) <- "double"
storage.mode(edge_test) <- "double"
if (feature_spec == "targeted_edges22") {
return(list(
train = edge_train,
test = edge_test,
feature_names = colnames(edge_train)
))
}
edge_z <- standardize_from_training(edge_train, edge_test)
make_subfamily_matrix <- function(zmat) {
out <- vapply(all_subfamilies, function(sf) {
sf_edges <- unique(manifest[subfamily == sf, FC])
missing <- setdiff(sf_edges, colnames(zmat))
if (length(missing) > 0L) {
stop("Missing edge(s) for subfamily ", sf, ": ", paste(missing, collapse = ", "))
}
rowMeans(zmat[, sf_edges, drop = FALSE])
}, numeric(nrow(zmat)))
colnames(out) <- paste0("subfamily_", all_subfamilies)
out
}
sf_train <- make_subfamily_matrix(edge_z$train)
sf_test <- make_subfamily_matrix(edge_z$test)
if (feature_spec == "subfamilies5") {
return(list(
train = sf_train,
test = sf_test,
feature_names = colnames(sf_train)
))
}
if (feature_spec != "axes2") {
stop("Unknown feature specification: ", feature_spec)
}
axis1_train_raw <- rowMeans(
sf_train[, paste0("subfamily_", axis1_subfamilies), drop = FALSE]
)
axis1_test_raw <- rowMeans(
sf_test[, paste0("subfamily_", axis1_subfamilies), drop = FALSE]
)
axis2_train_raw <- rowMeans(
sf_train[, paste0("subfamily_", axis2_subfamilies), drop = FALSE]
)
axis2_test_raw <- rowMeans(
sf_test[, paste0("subfamily_", axis2_subfamilies), drop = FALSE]
)
axis_raw_train <- cbind(axis1_train_raw, axis2_train_raw)
axis_raw_test <- cbind(axis1_test_raw, axis2_test_raw)
colnames(axis_raw_train) <- c(
"axis1_cortical_integration",
"axis2_perceptual_motor_coupling"
)
colnames(axis_raw_test) <- colnames(axis_raw_train)
axis_z <- standardize_from_training(axis_raw_train, axis_raw_test)
axis_z$train[, 1L] <- -axis_z$train[, 1L]
axis_z$test[, 1L] <- -axis_z$test[, 1L]
colnames(axis_z$train) <- c("axis1_disconnection", "axis2_hypercoupling")
colnames(axis_z$test) <- colnames(axis_z$train)
list(
train = axis_z$train,
test = axis_z$test,
feature_names = colnames(axis_z$train)
)
}
make_nuisance_design <- function(
train_meta,
test_meta,
mode,
use_motion,
include_icv = FALSE
) {
train_meta <- as.data.frame(copy(train_meta))
test_meta <- as.data.frame(copy(test_meta))
train_meta$Group <- factor(train_meta$Group, levels = c("HC", "MDD"))
test_meta$Group <- factor(test_meta$Group, levels = c("HC", "MDD"))
train_meta$Sex <- factor(train_meta$Sex, levels = c("Male", "Female"))
test_meta$Sex <- factor(test_meta$Sex, levels = c("Male", "Female"))
terms <- c("age", "Sex")
if (isTRUE(use_motion)) {
motion_median <- median(train_meta$Motion[is.finite(train_meta$Motion)], na.rm = TRUE)
if (!is.finite(motion_median)) motion_median <- 0
train_meta$Motion_missing <- as.numeric(!is.finite(train_meta$Motion))
test_meta$Motion_missing <- as.numeric(!is.finite(test_meta$Motion))
train_meta$Motion[!is.finite(train_meta$Motion)] <- motion_median
test_meta$Motion[!is.finite(test_meta$Motion)] <- motion_median
terms <- c(terms, "Motion", "Motion_missing")
}
if (isTRUE(include_icv)) {
icv_median <- median(train_meta$ICV[is.finite(train_meta$ICV)], na.rm = TRUE)
if (!is.finite(icv_median)) icv_median <- 0
train_meta$ICV_missing <- as.numeric(!is.finite(train_meta$ICV))
test_meta$ICV_missing <- as.numeric(!is.finite(test_meta$ICV))
train_meta$ICV[!is.finite(train_meta$ICV)] <- icv_median
test_meta$ICV[!is.finite(test_meta$ICV)] <- icv_median
terms <- c(terms, "ICV", "ICV_missing")
}
if (mode == "incremental_beyond_group") {
terms <- c("Group", terms)
} else if (mode != "total_without_group") {
stop("Unknown analysis mode: ", mode)
}
f <- as.formula(paste("~", paste(terms, collapse = " + ")))
C_train <- model.matrix(f, data = train_meta)
C_test <- model.matrix(f, data = test_meta)
missing_test_cols <- setdiff(colnames(C_train), colnames(C_test))
if (length(missing_test_cols) > 0L) {
for (nm in missing_test_cols) {
new_col <- matrix(
0,
nrow = nrow(C_test),
ncol = 1L,
dimnames = list(NULL, nm)
)
C_test <- cbind(C_test, new_col)
}
}
extra_test_cols <- setdiff(colnames(C_test), colnames(C_train))
if (length(extra_test_cols) > 0L) {
C_test <- C_test[, setdiff(colnames(C_test), extra_test_cols), drop = FALSE]
}
C_test <- C_test[, colnames(C_train), drop = FALSE]
list(
train = C_train,
test = C_test,
formula = deparse(f)
)
}
fit_matrix_ols <- function(C, Z) {
fit <- lm.fit(x = C, y = Z)
beta <- fit$coefficients
if (is.null(dim(beta))) beta <- matrix(beta, ncol = 1L)
beta[!is.finite(beta)] <- 0
beta
}
residualize_fold <- function(
X_train,
X_test,
Y_train,
Y_test,
C_train,
C_test
) {
beta_x <- fit_matrix_ols(C_train, X_train)
beta_y <- fit_matrix_ols(C_train, Y_train)
X_train_resid <- X_train - C_train %*% beta_x
X_test_resid <- X_test - C_test %*% beta_x
Y_train_resid <- Y_train - C_train %*% beta_y
Y_test_baseline <- C_test %*% beta_y
list(
X_train = X_train_resid,
X_test = X_test_resid,
Y_train = Y_train_resid,
Y_test_actual = Y_test,
Y_test_baseline = Y_test_baseline,
beta_x = beta_x,
beta_y = beta_y
)
}make_stratified_foldid <- function(strata, v, seed) {
strata <- as.character(strata)
if (v < 2L) stop("At least two folds are required.")
set.seed(seed)
foldid <- integer(length(strata))
for (lev in sort(unique(strata))) {
idx <- which(strata == lev)
idx <- sample(idx, length(idx), replace = FALSE)
foldid[idx] <- rep(seq_len(v), length.out = length(idx))
}
foldid
}
make_outer_fold_table <- function(group, v, repeats, seed) {
rbindlist(lapply(seq_len(repeats), function(r) {
data.table(
row_id = seq_along(group),
repeat_id = r,
outer_fold = make_stratified_foldid(
group,
v = v,
seed = seed + r * 1009L
)
)
}))
}
fit_glmnet_path <- function(x, y, alpha_value) {
glmnet::glmnet(
x = x,
y = y,
family = "mgaussian",
alpha = alpha_value,
nlambda = glmnet_nlambda,
lambda.min.ratio = glmnet_lambda_min_ratio,
standardize = TRUE,
standardize.response = TRUE,
intercept = TRUE,
thresh = 1e-7
)
}
coerce_mgaussian_prediction <- function(pred, n_obs, n_resp) {
d <- dim(pred)
if (is.null(d)) {
stop("glmnet returned a prediction without dimensions.")
}
if (length(d) == 2L) {
if (d[1L] == n_obs && d[2L] == n_resp) {
return(array(pred, dim = c(n_obs, n_resp, 1L)))
}
if (d[1L] == n_resp && d[2L] == n_obs) {
return(array(t(pred), dim = c(n_obs, n_resp, 1L)))
}
}
if (length(d) == 3L) {
if (d[1L] == n_obs && d[2L] == n_resp) return(pred)
if (d[1L] == n_resp && d[2L] == n_obs) return(aperm(pred, c(2L, 1L, 3L)))
if (d[1L] == n_obs && d[3L] == n_resp) return(aperm(pred, c(1L, 3L, 2L)))
}
stop(
"Unexpected multi-response prediction dimensions: ",
paste(d, collapse = " x ")
)
}
path_mse_on_ratio_grid <- function(
fit,
x_validation,
y_validation,
ratio_grid,
response_scale
) {
pred_raw <- predict(
fit,
newx = x_validation,
s = fit$lambda,
type = "response"
)
pred <- coerce_mgaussian_prediction(
pred_raw,
n_obs = nrow(x_validation),
n_resp = ncol(y_validation)
)
response_scale <- as.numeric(response_scale)
response_scale[!is.finite(response_scale) | response_scale <= 0] <- 1
n_lambda <- dim(pred)[3L]
mse <- vapply(seq_len(n_lambda), function(j) {
error <- y_validation - pred[, , j, drop = TRUE]
error <- sweep(error, 2L, response_scale, "/")
mean(error^2)
}, numeric(1L))
lambda_ratio <- fit$lambda / max(fit$lambda)
ord <- order(lambda_ratio)
approx(
x = log(lambda_ratio[ord]),
y = mse[ord],
xout = log(ratio_grid),
rule = 2,
ties = "ordered"
)$y
}
predict_one_lambda <- function(fit, x_new, lambda_value, n_resp) {
raw <- predict(
fit,
newx = x_new,
s = lambda_value,
type = "response"
)
arr <- coerce_mgaussian_prediction(raw, nrow(x_new), n_resp)
arr[, , 1L, drop = TRUE]
}
select_tuning_candidate <- function(tuning_table, rule) {
tuning_table <- copy(tuning_table[is.finite(mean_mse)])
if (nrow(tuning_table) == 0L) stop("All inner-loop models failed.")
min_row <- tuning_table[which.min(mean_mse)]
if (rule == "minimum") return(min_row)
if (rule != "one_se") stop("Unknown lambda selection rule: ", rule)
threshold <- min_row$mean_mse + min_row$se_mse
candidates <- tuning_table[mean_mse <= threshold]
# Prefer stronger regularization; if tied, prefer the more ridge-like model.
setorder(candidates, -lambda_ratio, alpha, mean_mse)
candidates[1L]
}
inner_tune <- function(
edge_train,
y_train,
meta_train,
feature_spec,
manifest,
mode,
use_motion,
include_icv,
seed
) {
min_group_n <- min(table(meta_train$Group))
if (min_group_n < 2L) {
stop("Too few participants per group for any valid inner cross-validation.")
}
v <- min(inner_folds, as.integer(min_group_n))
if (v < 2L) {
stop("Too few participants per group for inner cross-validation.")
}
if (v < 3L) {
warning(
"Only 2 participants per group remain in the training split; using 2-fold inner validation instead of the default 3+ folds.",
call. = FALSE
)
}
foldid <- make_stratified_foldid(meta_train$Group, v = v, seed = seed)
error_array <- array(
NA_real_,
dim = c(length(alpha_grid), v, length(lambda_ratio_grid)),
dimnames = list(
alpha = as.character(alpha_grid),
fold = as.character(seq_len(v)),
lambda_ratio = sprintf("%.8f", lambda_ratio_grid)
)
)
for (fold in seq_len(v)) {
idx_validation <- which(foldid == fold)
idx_training <- which(foldid != fold)
feat <- build_fold_features(
edge_train[idx_training, , drop = FALSE],
edge_train[idx_validation, , drop = FALSE],
feature_spec = feature_spec,
manifest = manifest
)
nuisance <- make_nuisance_design(
meta_train[idx_training],
meta_train[idx_validation],
mode = mode,
use_motion = use_motion,
include_icv = include_icv
)
res <- residualize_fold(
X_train = feat$train,
X_test = feat$test,
Y_train = y_train[idx_training, , drop = FALSE],
Y_test = y_train[idx_validation, , drop = FALSE],
C_train = nuisance$train,
C_test = nuisance$test
)
for (a in seq_along(alpha_grid)) {
fit <- tryCatch(
fit_glmnet_path(res$X_train, res$Y_train, alpha_grid[a]),
error = function(e) e
)
if (inherits(fit, "error")) {
error_array[a, fold, ] <- Inf
} else {
error_array[a, fold, ] <- path_mse_on_ratio_grid(
fit,
x_validation = res$X_test,
y_validation = y_train[idx_validation, , drop = FALSE] -
res$Y_test_baseline,
ratio_grid = lambda_ratio_grid,
response_scale = apply(res$Y_train, 2L, sd)
)
}
}
}
tuning_rows <- list()
counter <- 1L
for (a in seq_along(alpha_grid)) {
fold_errors <- error_array[a, , , drop = TRUE]
if (is.null(dim(fold_errors))) {
fold_errors <- matrix(fold_errors, nrow = v)
}
for (j in seq_along(lambda_ratio_grid)) {
values <- fold_errors[, j]
values <- values[is.finite(values)]
tuning_rows[[counter]] <- data.table(
alpha = alpha_grid[a],
lambda_ratio = lambda_ratio_grid[j],
mean_mse = if (length(values) == 0L) Inf else mean(values),
se_mse = if (length(values) <= 1L) 0 else sd(values) / sqrt(length(values)),
successful_inner_folds = length(values)
)
counter <- counter + 1L
}
}
tuning_table <- rbindlist(tuning_rows)
selected <- select_tuning_candidate(tuning_table, lambda_selection_rule)
list(selected = selected, tuning_table = tuning_table)
}
extract_selected_weights <- function(
fit,
lambda_index,
feature_names,
outcome_names
) {
if (!is.list(fit$beta)) {
stop("Expected a list of coefficient matrices for family='mgaussian'.")
}
coefficient_matrix <- do.call(cbind, lapply(fit$beta, function(b) {
as.numeric(b[, lambda_index])
}))
rownames(coefficient_matrix) <- feature_names
colnames(coefficient_matrix) <- outcome_names
weight_long <- as.data.table(as.table(coefficient_matrix))
setnames(weight_long, c("feature", "outcome", "weight"))
weight_long[, feature_l2_norm := sqrt(sum(weight^2)), by = feature]
weight_long
}
run_nested_mvpa <- function(
analysis_data,
edge_cols,
outcome_cols,
feature_spec,
mode,
manifest,
outer_fold_table,
use_motion,
include_icv,
seed
) {
edge_matrix <- as.matrix(analysis_data[, ..edge_cols])
outcome_matrix <- as.matrix(analysis_data[, ..outcome_cols])
storage.mode(edge_matrix) <- "double"
storage.mode(outcome_matrix) <- "double"
prediction_list <- list()
hyperparameter_list <- list()
weight_list <- list()
pred_counter <- 1L
hyper_counter <- 1L
weight_counter <- 1L
fold_keys <- unique(outer_fold_table[, .(repeat_id, outer_fold)])
setorder(fold_keys, repeat_id, outer_fold)
for (k in seq_len(nrow(fold_keys))) {
current_repeat_id <- fold_keys$repeat_id[k]
current_outer_fold <- fold_keys$outer_fold[k]
test_idx <- outer_fold_table[
repeat_id == current_repeat_id &
outer_fold == current_outer_fold,
row_id
]
train_idx <- setdiff(seq_len(nrow(analysis_data)), test_idx)
log_progress(
feature_spec, " | ", mode,
" | outer repeat ", current_repeat_id,
", fold ", current_outer_fold
)
train_group_counts <- table(analysis_data$Group[train_idx])
if (length(train_group_counts) < 2L || min(train_group_counts) < 2L) {
log_progress(
feature_spec, " | ", mode,
" | skipping outer repeat ", current_repeat_id,
", fold ", current_outer_fold,
" because training groups are too small: ",
paste(names(train_group_counts), train_group_counts, sep = "=", collapse = "; ")
)
next
}
tuned <- inner_tune(
edge_train = edge_matrix[train_idx, , drop = FALSE],
y_train = outcome_matrix[train_idx, , drop = FALSE],
meta_train = analysis_data[train_idx],
feature_spec = feature_spec,
manifest = manifest,
mode = mode,
use_motion = use_motion,
include_icv = include_icv,
seed = seed + current_repeat_id * 10007L + current_outer_fold * 101L
)
feat <- build_fold_features(
edge_matrix[train_idx, , drop = FALSE],
edge_matrix[test_idx, , drop = FALSE],
feature_spec = feature_spec,
manifest = manifest
)
nuisance <- make_nuisance_design(
analysis_data[train_idx],
analysis_data[test_idx],
mode = mode,
use_motion = use_motion,
include_icv = include_icv
)
res <- residualize_fold(
X_train = feat$train,
X_test = feat$test,
Y_train = outcome_matrix[train_idx, , drop = FALSE],
Y_test = outcome_matrix[test_idx, , drop = FALSE],
C_train = nuisance$train,
C_test = nuisance$test
)
final_fit <- fit_glmnet_path(
res$X_train,
res$Y_train,
alpha_value = tuned$selected$alpha
)
final_ratio <- final_fit$lambda / max(final_fit$lambda)
lambda_index <- which.min(abs(
log(final_ratio) - log(tuned$selected$lambda_ratio)
))
selected_lambda <- final_fit$lambda[lambda_index]
predicted_residual <- predict_one_lambda(
final_fit,
x_new = res$X_test,
lambda_value = selected_lambda,
n_resp = ncol(outcome_matrix)
)
predicted_full <- res$Y_test_baseline + predicted_residual
prediction_list[[pred_counter]] <- data.table(
feature_set = feature_spec,
analysis_mode = mode,
repeat_id = current_repeat_id,
outer_fold = current_outer_fold,
row_id = rep(test_idx, times = length(outcome_cols)),
subject = rep(analysis_data$subject[test_idx], times = length(outcome_cols)),
Group = rep(as.character(analysis_data$Group[test_idx]), times = length(outcome_cols)),
outcome = rep(outcome_cols, each = length(test_idx)),
actual = as.vector(outcome_matrix[test_idx, , drop = FALSE]),
baseline_prediction = as.vector(res$Y_test_baseline),
full_prediction = as.vector(predicted_full)
)
pred_counter <- pred_counter + 1L
hyperparameter_list[[hyper_counter]] <- data.table(
feature_set = feature_spec,
analysis_mode = mode,
repeat_id = current_repeat_id,
outer_fold = current_outer_fold,
alpha = tuned$selected$alpha,
selected_lambda_ratio = tuned$selected$lambda_ratio,
selected_lambda = selected_lambda,
inner_mean_mse = tuned$selected$mean_mse,
inner_se_mse = tuned$selected$se_mse,
nuisance_formula = nuisance$formula,
n_train = length(train_idx),
n_test = length(test_idx)
)
hyper_counter <- hyper_counter + 1L
fold_weights <- extract_selected_weights(
final_fit,
lambda_index = lambda_index,
feature_names = feat$feature_names,
outcome_names = outcome_cols
)
fold_weights[, `:=`(
feature_set = feature_spec,
analysis_mode = mode,
repeat_id = current_repeat_id,
outer_fold = current_outer_fold,
alpha = tuned$selected$alpha,
selected_lambda = selected_lambda
)]
weight_list[[weight_counter]] <- fold_weights
weight_counter <- weight_counter + 1L
rm(final_fit, feat, nuisance, res, tuned)
gc(verbose = FALSE)
}
list(
predictions = rbindlist(prediction_list),
hyperparameters = rbindlist(hyperparameter_list),
weights = rbindlist(weight_list)
)
}multivariate_sse <- function(design, Y) {
fit <- lm.fit(design, Y)
sum(fit$residuals^2)
}
permute_indices_within_strata <- function(strata) {
strata <- as.character(strata)
out <- seq_along(strata)
for (lev in unique(strata)) {
idx <- which(strata == lev)
out[idx] <- sample(idx, length(idx), replace = FALSE)
}
out
}
freedman_lane_partial_r2 <- function(
Y,
data,
target,
reduced_terms,
strata,
n_permutations,
seed
) {
reduced_formula <- as.formula(paste("~", paste(reduced_terms, collapse = " + ")))
full_formula <- as.formula(paste(
"~", paste(c(reduced_terms, target), collapse = " + ")
))
C_reduced <- model.matrix(reduced_formula, data = data)
C_full <- model.matrix(full_formula, data = data)
reduced_fit <- lm.fit(C_reduced, Y)
full_fit <- lm.fit(C_full, Y)
sse_reduced <- sum(reduced_fit$residuals^2)
sse_full <- sum(full_fit$residuals^2)
observed <- (sse_reduced - sse_full) / sse_reduced
set.seed(seed)
permutation_stats <- numeric(n_permutations)
for (b in seq_len(n_permutations)) {
perm_idx <- permute_indices_within_strata(strata)
Y_perm <- reduced_fit$fitted.values + reduced_fit$residuals[perm_idx, , drop = FALSE]
sse_r <- multivariate_sse(C_reduced, Y_perm)
sse_f <- multivariate_sse(C_full, Y_perm)
permutation_stats[b] <- (sse_r - sse_f) / sse_r
}
data.table(
target = target,
partial_r2 = observed,
p_permutation = (1 + sum(permutation_stats >= observed)) /
(n_permutations + 1),
n_permutations = n_permutations,
reduced_formula = deparse(reduced_formula),
full_formula = deparse(full_formula)
)
}
calculate_basic_metrics <- function(observed, predicted) {
keep <- is.finite(observed) & is.finite(predicted)
observed <- observed[keep]
predicted <- predicted[keep]
if (length(observed) < 3L) {
return(data.table(n = length(observed), rmse = NA_real_, mae = NA_real_,
q2 = NA_real_, correlation = NA_real_))
}
sst <- sum((observed - mean(observed))^2)
data.table(
n = length(observed),
rmse = sqrt(mean((observed - predicted)^2)),
mae = mean(abs(observed - predicted)),
q2 = if (sst <= 0) NA_real_ else 1 - sum((observed - predicted)^2) / sst,
correlation = suppressWarnings(cor(observed, predicted))
)
}
average_oof_predictions <- function(predictions) {
if (nrow(predictions) == 0L) {
return(data.table(
feature_set = character(0L),
analysis_mode = character(0L),
row_id = integer(0L),
subject = character(0L),
Group = character(0L),
outcome = character(0L),
actual = numeric(0L),
baseline_prediction = numeric(0L),
full_prediction = numeric(0L),
n_outer_predictions = integer(0L)
))
}
predictions[, .(
actual = actual[1L],
baseline_prediction = mean(baseline_prediction),
full_prediction = mean(full_prediction),
n_outer_predictions = .N
), by = .(
feature_set, analysis_mode, row_id, subject, Group, outcome
)]
}
coefficient_metric_tables <- function(prediction_average) {
source_long <- melt(
prediction_average,
id.vars = c("feature_set", "analysis_mode", "row_id", "subject", "Group", "outcome", "actual"),
measure.vars = c("baseline_prediction", "full_prediction"),
variable.name = "prediction_source",
value.name = "prediction"
)
by_outcome <- source_long[, calculate_basic_metrics(actual, prediction),
by = .(feature_set, analysis_mode, outcome, prediction_source)]
by_outcome[, condition := sub("__b[0-9]+$", "", outcome)]
# Equalize response scales for an overall multivariate score.
source_long[, outcome_sd := sd(actual),
by = .(outcome, prediction_source)]
source_long[!is.finite(outcome_sd) | outcome_sd <= 0, outcome_sd := 1]
source_long[, `:=`(
actual_scaled = actual / outcome_sd,
prediction_scaled = prediction / outcome_sd
)]
overall <- source_long[, calculate_basic_metrics(actual_scaled, prediction_scaled),
by = .(feature_set, analysis_mode, prediction_source)]
list(by_outcome = by_outcome, overall = overall)
}
reconstruct_curve_predictions <- function(
prediction_average,
basis_grid,
conditions,
basis_df
) {
B_grid <- as.matrix(basis_grid[, paste0("B", seq_len(basis_df)), with = FALSE])
output <- list()
counter <- 1L
model_keys <- unique(prediction_average[, .(feature_set, analysis_mode)])
for (m in seq_len(nrow(model_keys))) {
fs <- model_keys$feature_set[m]
md <- model_keys$analysis_mode[m]
model_data <- prediction_average[
feature_set == fs & analysis_mode == md
]
for (condition_value in conditions) {
coefficient_names <- paste0(
condition_value, "__b", sprintf("%02d", seq_len(basis_df))
)
dat <- model_data[outcome %in% coefficient_names]
if (nrow(dat) == 0L) next
make_wide <- function(value_col) {
dcast(
dat,
row_id + subject + Group ~ outcome,
value.var = value_col
)
}
observed_wide <- make_wide("actual")
baseline_wide <- make_wide("baseline_prediction")
full_wide <- make_wide("full_prediction")
setorder(observed_wide, row_id)
setorder(baseline_wide, row_id)
setorder(full_wide, row_id)
observed_coef <- as.matrix(observed_wide[, ..coefficient_names])
baseline_coef <- as.matrix(baseline_wide[, ..coefficient_names])
full_coef <- as.matrix(full_wide[, ..coefficient_names])
observed_eta <- observed_coef %*% t(B_grid)
baseline_eta <- baseline_coef %*% t(B_grid)
full_eta <- full_coef %*% t(B_grid)
n_subjects <- nrow(observed_wide)
n_times <- nrow(basis_grid)
output[[counter]] <- data.table(
feature_set = fs,
analysis_mode = md,
condition = condition_value,
row_id = rep(observed_wide$row_id, each = n_times),
subject = rep(observed_wide$subject, each = n_times),
Group = rep(observed_wide$Group, each = n_times),
time_ms = rep(basis_grid$time_ms, times = n_subjects),
observed_probability = as.vector(t(plogis(observed_eta))),
baseline_probability = as.vector(t(plogis(baseline_eta))),
full_probability = as.vector(t(plogis(full_eta)))
)
counter <- counter + 1L
}
}
rbindlist(output)
}
curve_metric_tables <- function(curve_predictions) {
long <- melt(
curve_predictions,
id.vars = c(
"feature_set", "analysis_mode", "condition", "row_id",
"subject", "Group", "time_ms", "observed_probability"
),
measure.vars = c("baseline_probability", "full_probability"),
variable.name = "prediction_source",
value.name = "predicted_probability"
)
by_condition <- long[, calculate_basic_metrics(
observed_probability,
predicted_probability
), by = .(feature_set, analysis_mode, condition, prediction_source)]
overall <- long[, calculate_basic_metrics(
observed_probability,
predicted_probability
), by = .(feature_set, analysis_mode, prediction_source)]
list(by_condition = by_condition, overall = overall)
}
make_group_difference_table <- function(curve_predictions) {
group_means <- curve_predictions[, .(
observed_probability = mean(observed_probability),
baseline_probability = mean(baseline_probability),
full_probability = mean(full_probability)
), by = .(feature_set, analysis_mode, condition, time_ms, Group)]
long <- melt(
group_means,
id.vars = c("feature_set", "analysis_mode", "condition", "time_ms", "Group"),
measure.vars = c(
"observed_probability", "baseline_probability", "full_probability"
),
variable.name = "series",
value.name = "probability"
)
wide <- dcast(
long,
feature_set + analysis_mode + condition + time_ms + series ~ Group,
value.var = "probability"
)
if (!all(c("HC", "MDD") %in% names(wide))) {
stop("Both HC and MDD are required for group-difference curves.")
}
wide[, difference_mdd_minus_hc := MDD - HC]
result <- dcast(
wide,
feature_set + analysis_mode + condition + time_ms ~ series,
value.var = "difference_mdd_minus_hc"
)
setnames(
result,
old = c("observed_probability", "baseline_probability", "full_probability"),
new = c("observed_difference", "baseline_difference", "full_difference")
)
result
}
prediction_alignment_test <- function(
prediction_average,
n_permutations,
n_bootstrap,
seed
) {
model_keys <- unique(prediction_average[, .(feature_set, analysis_mode)])
output <- list()
counter <- 1L
for (m in seq_len(nrow(model_keys))) {
fs <- model_keys$feature_set[m]
md <- model_keys$analysis_mode[m]
dat <- prediction_average[feature_set == fs & analysis_mode == md]
conditions <- unique(sub("__b[0-9]+$", "", dat$outcome))
target_sets <- c("ALL_CONDITIONS", conditions)
for (target_condition in target_sets) {
use <- if (target_condition == "ALL_CONDITIONS") {
dat
} else {
dat[sub("__b[0-9]+$", "", outcome) == target_condition]
}
actual_wide <- dcast(use, row_id + Group ~ outcome, value.var = "actual")
baseline_wide <- dcast(
use, row_id + Group ~ outcome, value.var = "baseline_prediction"
)
full_wide <- dcast(use, row_id + Group ~ outcome, value.var = "full_prediction")
setorder(actual_wide, row_id)
setorder(baseline_wide, row_id)
setorder(full_wide, row_id)
outcome_names <- setdiff(names(actual_wide), c("row_id", "Group"))
A <- as.matrix(actual_wide[, ..outcome_names])
B <- as.matrix(baseline_wide[, ..outcome_names])
F <- as.matrix(full_wide[, ..outcome_names])
scale_sd <- apply(A, 2L, sd)
scale_sd[!is.finite(scale_sd) | scale_sd <= 0] <- 1
Y_residual <- sweep(A - B, 2L, scale_sd, "/")
P_increment <- sweep(F - B, 2L, scale_sd, "/")
stat_fun <- function(Ymat, Pmat) {
mean(Ymat^2) - mean((Ymat - Pmat)^2)
}
observed_stat <- stat_fun(Y_residual, P_increment)
null_mse <- mean(Y_residual^2)
percent_reduction <- if (null_mse <= 0) NA_real_ else observed_stat / null_mse
set.seed(seed + m * 1009L + counter)
permutation_stats <- numeric(n_permutations)
for (b in seq_len(n_permutations)) {
perm_idx <- if (md == "incremental_beyond_group") {
permute_indices_within_strata(actual_wide$Group)
} else {
sample(seq_len(nrow(actual_wide)))
}
permutation_stats[b] <- stat_fun(
Y_residual[perm_idx, , drop = FALSE],
P_increment
)
}
bootstrap_stats <- numeric(n_bootstrap)
for (b in seq_len(n_bootstrap)) {
idx <- sample(seq_len(nrow(actual_wide)), replace = TRUE)
bootstrap_stats[b] <- stat_fun(
Y_residual[idx, , drop = FALSE],
P_increment[idx, , drop = FALSE]
)
}
output[[counter]] <- data.table(
feature_set = fs,
analysis_mode = md,
condition = target_condition,
delta_mse = observed_stat,
proportional_mse_reduction = percent_reduction,
bootstrap_ci_low = unname(quantile(bootstrap_stats, 0.025, na.rm = TRUE)),
bootstrap_ci_high = unname(quantile(bootstrap_stats, 0.975, na.rm = TRUE)),
p_fixed_prediction_permutation =
(1 + sum(permutation_stats >= observed_stat)) /
(n_permutations + 1),
n_subjects = nrow(actual_wide),
n_outcomes = ncol(Y_residual),
n_permutations = n_permutations,
n_bootstrap = n_bootstrap
)
counter <- counter + 1L
}
}
rbindlist(output)
}5 Load and assemble fMRI features
prepared <- readRDS(prepared_rds)
if (!all(c("analysis_long", "subjects", "edge_map", "roi") %in% names(prepared))) {
stop("The prepared RDS does not have the expected Step-1 structure.")
}
manifest <- fread(manifest_file)
axis_scores <- fread(axis_scores_file)
analysis_long <- as.data.table(prepared$analysis_long)
subjects_fmri <- as.data.table(prepared$subjects)
required_manifest_columns <- c("FC", "axis", "subfamily")
if (!all(required_manifest_columns %in% names(manifest))) {
stop("Manifest is missing: ", paste(
setdiff(required_manifest_columns, names(manifest)), collapse = ", "
))
}
target_edge_cols <- unique(manifest$FC)
if (length(target_edge_cols) != 22L) {
warning("Expected 22 unique targeted edges but found ", length(target_edge_cols), ".")
}
edge_long <- analysis_long[FC %in% target_edge_cols]
edge_wide <- dcast(
edge_long,
eye_subject_id + fmri_subject_id ~ FC,
value.var = "fc_value"
)
missing_edge_columns <- setdiff(target_edge_cols, names(edge_wide))
if (length(missing_edge_columns) > 0L) {
stop(
"Targeted FC column(s) absent after reshaping: ",
paste(missing_edge_columns, collapse = ", ")
)
}
subject_columns <- c(
"eye_subject_id", "fmri_subject_id", "Group", "Sex",
"age", "ICV", "Motion"
)
if (!all(subject_columns %in% names(subjects_fmri))) {
stop("Prepared subject table is missing: ", paste(
setdiff(subject_columns, names(subjects_fmri)), collapse = ", "
))
}
fmri_table <- merge(
unique(subjects_fmri[, ..subject_columns]),
edge_wide,
by = c("eye_subject_id", "fmri_subject_id"),
all = FALSE,
sort = FALSE
)
axis_score_columns <- c(
"axis1_cortical_integration",
"axis2_perceptual_motor_coupling",
"axis1_disconnection",
"axis2_hypercoupling",
"cross_axis_imbalance",
paste0("subfamily_", all_subfamilies)
)
missing_axis_cols <- setdiff(axis_score_columns, names(axis_scores))
if (length(missing_axis_cols) > 0L) {
stop("Axis-score export is missing: ", paste(missing_axis_cols, collapse = ", "))
}
fmri_table <- merge(
fmri_table,
axis_scores[, c("eye_subject_id", "fmri_subject_id", axis_score_columns), with = FALSE],
by = c("eye_subject_id", "fmri_subject_id"),
all.x = TRUE,
sort = FALSE
)
fmri_table[, id_key := normalize_subject_id(eye_subject_id)]
check_unique_key(fmri_table, "id_key", "fMRI table")
if (anyNA(fmri_table[, ..target_edge_cols])) {
stop("At least one matched fMRI participant is missing a targeted FC edge.")
}
fmri_audit <- data.table(
metric = c(
"prepared_fmri_participants",
"targeted_edges",
"participants_with_all_targeted_edges",
"motion_available"
),
value = c(
nrow(subjects_fmri),
length(target_edge_cols),
nrow(fmri_table),
sum(is.finite(fmri_table$Motion))
)
)
print(fmri_audit)## metric value
## <char> <int>
## 1: prepared_fmri_participants 203
## 2: targeted_edges 22
## 3: participants_with_all_targeted_edges 203
## 4: motion_available 0fwrite(fmri_audit, file.path(output_dir, "fmri_input_audit.csv"))6 Create subject-level temporal eye phenotypes
eye_bins <- fread(eye_bins_file)
required_eye_columns <- c(
"subject", "group", "expression_label", "start_code",
"trial", "time_center_ms", "n_code2", "n_code1", "n_valid"
)
if (!all(required_eye_columns %in% names(eye_bins))) {
stop("Eye-bin file is missing: ", paste(
setdiff(required_eye_columns, names(eye_bins)), collapse = ", "
))
}
eye_bins[, `:=`(
subject = trimws(as.character(subject)),
group = toupper(trimws(as.character(group))),
expression_label = trimws(as.character(expression_label)),
start_code = toupper(trimws(as.character(start_code))),
time_center_ms = as.numeric(time_center_ms),
n_code2 = as.numeric(n_code2),
n_code1 = as.numeric(n_code1)
)]
eye_bins[, condition := paste(expression_label, start_code, sep = "__")]
missing_conditions <- setdiff(primary_eye_conditions, unique(eye_bins$condition))
if (length(missing_conditions) > 0L) {
stop("Primary eye condition(s) absent from trial-bin file: ",
paste(missing_conditions, collapse = ", "))
}
# Audit that the frozen conditions still correspond to the prior BH results.
if (file.exists(eye_condition_tests_file)) {
condition_tests <- fread(eye_condition_tests_file)
if (all(c("condition", "q_bh_all_planned_conditions") %in% names(condition_tests))) {
print(condition_tests[condition %in% primary_eye_conditions])
}
}## 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: Sad__C Sad C 0.000000e+00 fit completed
## 3: Sad__E Sad E 2.569658e-05 fit completed
## q_bh_all_planned_conditions
## <num>
## 1: 0.0002725592
## 2: 0.0000000000
## 3: 0.0001156346eye_aggregated <- eye_bins[
condition %in% primary_eye_conditions,
.(
n_code2 = sum(n_code2),
n_code1 = sum(n_code1),
n_valid = sum(n_code2 + n_code1)
),
by = .(subject, group, condition, time_center_ms)
]
time_min <- min(eye_aggregated$time_center_ms)
time_max <- max(eye_aggregated$time_center_ms)
eye_basis_fun <- make_eye_basis_function(time_min, time_max, eye_basis_df)
basis_time_grid <- sort(unique(eye_aggregated$time_center_ms))
basis_matrix_grid <- eye_basis_fun(basis_time_grid)
basis_grid <- data.table(time_ms = basis_time_grid)
for (j in seq_len(ncol(basis_matrix_grid))) {
basis_grid[[paste0("B", j)]] <- basis_matrix_grid[, j]
}
fwrite(basis_grid, file.path(output_dir, "eye_temporal_basis_grid.csv"))
split_eye <- split(
eye_aggregated,
by = c("subject", "group", "condition"),
keep.by = TRUE,
drop = TRUE
)
coefficient_rows <- list()
audit_rows <- list()
for (i in seq_along(split_eye)) {
dat <- as.data.table(split_eye[[i]])
fitted <- fit_subject_eye_basis(dat, eye_basis_fun, eye_basis_ridge)
coefficient_rows[[i]] <- data.table(
subject = as.character(dat$subject[1L]),
group = as.character(dat$group[1L]),
condition = as.character(dat$condition[1L]),
basis_index = as.integer(seq_along(fitted$beta)),
coefficient = as.numeric(fitted$beta)
)
audit_rows[[i]] <- data.table(
subject = as.character(dat$subject[1L]),
group = as.character(dat$group[1L]),
condition = as.character(dat$condition[1L]),
n_bins = as.integer(fitted$n_bins),
total_valid_samples = as.integer(fitted$total_valid_samples),
weighted_logit_rmse = as.numeric(fitted$weighted_logit_rmse),
first_time_ms = as.numeric(fitted$first_time_ms),
last_time_ms = as.numeric(fitted$last_time_ms)
)
}
# Force consistent column types before stacking across subjects/conditions.
eye_coeff_long <- rbindlist(coefficient_rows, use.names = TRUE, fill = FALSE)
eye_coeff_long[, `:=`(
subject = as.character(subject),
group = as.character(group),
condition = as.character(condition),
basis_index = as.integer(basis_index),
coefficient = as.numeric(coefficient)
)]
eye_basis_audit <- rbindlist(audit_rows, use.names = TRUE, fill = FALSE)
eye_basis_audit[, `:=`(
subject = as.character(subject),
group = as.character(group),
condition = as.character(condition),
n_bins = as.integer(n_bins),
total_valid_samples = as.integer(total_valid_samples),
weighted_logit_rmse = as.numeric(weighted_logit_rmse),
first_time_ms = as.numeric(first_time_ms),
last_time_ms = as.numeric(last_time_ms)
)]
eligible_subjects <- eye_basis_audit[
n_bins >= minimum_eye_bins_per_condition,
.(n_conditions = uniqueN(condition)),
by = subject
][n_conditions == length(primary_eye_conditions), subject]
eye_coeff_long <- eye_coeff_long[subject %in% eligible_subjects]
eye_basis_audit <- eye_basis_audit[subject %in% eligible_subjects]
eye_coeff_long[, outcome := paste0(
condition, "__b", sprintf("%02d", basis_index)
)]
eye_coeff_wide <- dcast(
eye_coeff_long,
subject + group ~ outcome,
value.var = "coefficient"
)
eye_coeff_wide[, id_key := normalize_subject_id(subject)]
check_unique_key(eye_coeff_wide, "id_key", "Eye-trajectory table")
outcome_cols <- setdiff(
names(eye_coeff_wide),
c("subject", "group", "id_key")
)
expected_outcomes <- length(primary_eye_conditions) * eye_basis_df
if (length(outcome_cols) != expected_outcomes) {
stop("Expected ", expected_outcomes, " eye outcomes but found ", length(outcome_cols), ".")
}
fwrite(eye_basis_audit, file.path(output_dir, "eye_basis_fit_audit.csv"))
fwrite(eye_coeff_wide, file.path(output_dir, "eye_subject_temporal_basis_coefficients.csv"))
print(eye_basis_audit[, .(
participants = as.integer(uniqueN(subject)),
median_bins = as.numeric(median(n_bins)),
min_bins = as.numeric(min(n_bins)),
median_weighted_logit_rmse = as.numeric(median(weighted_logit_rmse))
), by = .(group, condition)])## group condition participants median_bins min_bins
## <char> <char> <int> <num> <num>
## 1: HC Fear__C 47 91 84
## 2: MDD Fear__C 62 90 78
## 3: HC Sad__C 47 90 80
## 4: MDD Sad__C 62 90 76
## 5: HC Sad__E 47 100 100
## 6: MDD Sad__E 62 100 100
## median_weighted_logit_rmse
## <num>
## 1: 0.2422215
## 2: 0.2628336
## 3: 0.2370244
## 4: 0.2414341
## 5: 0.4362074
## 6: 0.48391417 Match modalities and freeze the analytic sample
matched <- merge(
fmri_table,
eye_coeff_wide,
by = "id_key",
all = FALSE,
suffixes = c("_fmri", "_eye"),
sort = FALSE
)
matched[, `:=`(
Group = factor(as.character(Group), levels = c("HC", "MDD")),
Sex = factor(as.character(Sex), levels = c("Male", "Female")),
eye_group = factor(group, levels = c("HC", "MDD")),
subject = subject
)]
if (any(as.character(matched$Group) != as.character(matched$eye_group))) {
bad <- matched[as.character(Group) != as.character(eye_group)][1L]
stop(
"Group mismatch after modality merge for ID ", bad$id_key,
": fMRI=", bad$Group, ", eye=", bad$eye_group
)
}
core_columns <- c(
"subject", "Group", "Sex", "age",
target_edge_cols, axis_score_columns, outcome_cols
)
complete_core <- complete.cases(matched[, ..core_columns])
analysis_data <- matched[complete_core]
analysis_data[, row_id := .I]
outcome_sd_check <- vapply(
analysis_data[, ..outcome_cols],
sd,
numeric(1L),
na.rm = TRUE
)
invalid_outcomes <- names(outcome_sd_check)[
!is.finite(outcome_sd_check) | outcome_sd_check <= 0
]
if (length(invalid_outcomes) > 0L) {
stop(
"Zero-variance or invalid eye-trajectory outcome(s): ",
paste(invalid_outcomes, collapse = ", ")
)
}
motion_fraction <- mean(is.finite(analysis_data$Motion))
use_motion <- motion_fraction >= minimum_motion_fraction
if (isTRUE(require_motion) && !use_motion) {
stop(
"Motion is available for only ", round(100 * motion_fraction, 1),
"% of the matched sample; require_motion=TRUE."
)
}
if (!use_motion) {
warning(
"Motion is available for only ", round(100 * motion_fraction, 1),
"% of the matched sample. Predictive results are provisional."
)
}## Warning: Motion is available for only 0% of the matched sample. Predictive
## results are provisional.match_audit <- data.table(
metric = c(
"eye_participants_with_eligible_curves",
"fmri_participants_with_targeted_edges",
"cross_modal_matched_participants",
"matched_hc",
"matched_mdd",
"motion_available_fraction",
"motion_used_as_nuisance",
"target_manifest_independently_prespecified"
),
value = as.character(c(
nrow(eye_coeff_wide),
nrow(fmri_table),
nrow(analysis_data),
sum(analysis_data$Group == "HC"),
sum(analysis_data$Group == "MDD"),
motion_fraction,
use_motion,
targets_independently_prespecified
))
)
print(match_audit)## metric value
## <char> <char>
## 1: eye_participants_with_eligible_curves 109
## 2: fmri_participants_with_targeted_edges 203
## 3: cross_modal_matched_participants 103
## 4: matched_hc 44
## 5: matched_mdd 59
## 6: motion_available_fraction 0
## 7: motion_used_as_nuisance 0
## 8: target_manifest_independently_prespecified 0fwrite(match_audit, file.path(output_dir, "cross_modal_match_audit.csv"))
if (nrow(analysis_data) < 60L || min(table(analysis_data$Group)) < 20L) {
stop(
"The matched sample is too small for the prespecified 5-fold nested CV. ",
"Reduce model complexity only with a documented revised analysis plan."
)
}8 Theory-guided multivariate association of the two axes
This test treats the B-spline coefficients jointly. It asks whether
each axis adds multivariate explanatory variance after diagnosis, age,
sex, motion when available, and the other axis. The permutation is a
Freedman-Lane residual permutation within diagnosis. When
eye_condition_set = "surviving_three", the conditions were
selected in a prior analysis of the same eye data, so these p values
remain secondary and exploratory. The all_nine rerun is the
recommended selection-bias sensitivity analysis.
association_data <- copy(analysis_data)
association_data[, Group := factor(Group, levels = c("HC", "MDD"))]
association_data[, Sex := factor(Sex, levels = c("Male", "Female"))]
nuisance_terms <- c("Group", "age", "Sex")
if (use_motion) {
motion_median <- median(
association_data$Motion[is.finite(association_data$Motion)],
na.rm = TRUE
)
association_data[, Motion_missing := as.numeric(!is.finite(Motion))]
association_data[!is.finite(Motion), Motion := motion_median]
nuisance_terms <- c(nuisance_terms, "Motion", "Motion_missing")
}
if (include_icv_sensitivity) {
icv_median <- median(
association_data$ICV[is.finite(association_data$ICV)],
na.rm = TRUE
)
association_data[, ICV_missing := as.numeric(!is.finite(ICV))]
association_data[!is.finite(ICV), ICV := icv_median]
nuisance_terms <- c(nuisance_terms, "ICV", "ICV_missing")
}
Y_association <- scale(as.matrix(association_data[, ..outcome_cols]))
axis_targets <- c("axis1_disconnection", "axis2_hypercoupling")
axis_association_results <- rbindlist(lapply(seq_along(axis_targets), function(i) {
target <- axis_targets[i]
other_axis <- setdiff(axis_targets, target)
freedman_lane_partial_r2(
Y = Y_association,
data = association_data,
target = target,
reduced_terms = c(nuisance_terms, other_axis),
strata = association_data$Group,
n_permutations = axis_permutations,
seed = random_seed + i * 100003L
)
}))
axis_association_results[, p_holm_two_axes := p.adjust(
p_permutation,
method = "holm"
)]
print(axis_association_results)## target partial_r2 p_permutation n_permutations
## <char> <num> <num> <int>
## 1: axis1_disconnection 0.007614068 0.7106447 2000
## 2: axis2_hypercoupling 0.010170836 0.4482759 2000
## reduced_formula
## <char>
## 1: ~Group + age + Sex + axis2_hypercoupling
## 2: ~Group + age + Sex + axis1_disconnection
## full_formula
## <char>
## 1: ~Group + age + Sex + axis2_hypercoupling + axis1_disconnection
## 2: ~Group + age + Sex + axis1_disconnection + axis2_hypercoupling
## p_holm_two_axes
## <num>
## 1: 0.8965517
## 2: 0.8965517fwrite(
axis_association_results,
file.path(output_dir, "two_axis_multivariate_temporal_association.csv")
)9 Repeated nested cross-validated multi-response prediction
The outcome is the complete set of B-spline coefficients for the
condition family chosen by eye_condition_set (default:
fearful center, sad center, and sad emotional-start).
glmnet uses a multi-response Gaussian elastic net, so each
fMRI feature is selected or shrunk jointly across the correlated
temporal outcomes. Axis construction and every nuisance-regression step
are repeated inside the training fold.
outer_fold_table <- make_outer_fold_table(
group = analysis_data$Group,
v = outer_folds,
repeats = outer_repeats,
seed = random_seed
)
fwrite(outer_fold_table, file.path(output_dir, "outer_cross_validation_folds.csv"))
mvpa_runs <- list()
run_counter <- 1L
for (feature_spec in feature_sets_to_run) {
for (mode in analysis_modes_to_run) {
mvpa_runs[[run_counter]] <- run_nested_mvpa(
analysis_data = analysis_data,
edge_cols = target_edge_cols,
outcome_cols = outcome_cols,
feature_spec = feature_spec,
mode = mode,
manifest = manifest,
outer_fold_table = outer_fold_table,
use_motion = use_motion,
include_icv = include_icv_sensitivity,
seed = random_seed + run_counter * 1000003L
)
run_counter <- run_counter + 1L
}
}## [2026-08-16 13:04:17] axes2 | total_without_group | outer repeat 1, fold 1
## [2026-08-16 13:04:19] axes2 | total_without_group | outer repeat 1, fold 2
## [2026-08-16 13:04:20] axes2 | total_without_group | outer repeat 1, fold 3
## [2026-08-16 13:04:21] axes2 | total_without_group | outer repeat 1, fold 4
## [2026-08-16 13:04:22] axes2 | total_without_group | outer repeat 1, fold 5
## [2026-08-16 13:04:23] axes2 | total_without_group | outer repeat 2, fold 1
## [2026-08-16 13:04:24] axes2 | total_without_group | outer repeat 2, fold 2
## [2026-08-16 13:04:25] axes2 | total_without_group | outer repeat 2, fold 3
## [2026-08-16 13:04:26] axes2 | total_without_group | outer repeat 2, fold 4
## [2026-08-16 13:04:27] axes2 | total_without_group | outer repeat 2, fold 5
## [2026-08-16 13:04:28] axes2 | total_without_group | outer repeat 3, fold 1
## [2026-08-16 13:04:29] axes2 | total_without_group | outer repeat 3, fold 2
## [2026-08-16 13:04:30] axes2 | total_without_group | outer repeat 3, fold 3
## [2026-08-16 13:04:31] axes2 | total_without_group | outer repeat 3, fold 4
## [2026-08-16 13:04:32] axes2 | total_without_group | outer repeat 3, fold 5
## [2026-08-16 13:04:33] axes2 | total_without_group | outer repeat 4, fold 1
## [2026-08-16 13:04:34] axes2 | total_without_group | outer repeat 4, fold 2
## [2026-08-16 13:04:35] axes2 | total_without_group | outer repeat 4, fold 3
## [2026-08-16 13:04:36] axes2 | total_without_group | outer repeat 4, fold 4
## [2026-08-16 13:04:37] axes2 | total_without_group | outer repeat 4, fold 5
## [2026-08-16 13:04:39] axes2 | total_without_group | outer repeat 5, fold 1
## [2026-08-16 13:04:40] axes2 | total_without_group | outer repeat 5, fold 2
## [2026-08-16 13:04:41] axes2 | total_without_group | outer repeat 5, fold 3
## [2026-08-16 13:04:42] axes2 | total_without_group | outer repeat 5, fold 4
## [2026-08-16 13:04:43] axes2 | total_without_group | outer repeat 5, fold 5
## [2026-08-16 13:04:44] axes2 | total_without_group | outer repeat 6, fold 1
## [2026-08-16 13:04:45] axes2 | total_without_group | outer repeat 6, fold 2
## [2026-08-16 13:04:46] axes2 | total_without_group | outer repeat 6, fold 3
## [2026-08-16 13:04:47] axes2 | total_without_group | outer repeat 6, fold 4
## [2026-08-16 13:04:48] axes2 | total_without_group | outer repeat 6, fold 5
## [2026-08-16 13:04:49] axes2 | total_without_group | outer repeat 7, fold 1
## [2026-08-16 13:04:50] axes2 | total_without_group | outer repeat 7, fold 2
## [2026-08-16 13:04:51] axes2 | total_without_group | outer repeat 7, fold 3
## [2026-08-16 13:04:52] axes2 | total_without_group | outer repeat 7, fold 4
## [2026-08-16 13:04:53] axes2 | total_without_group | outer repeat 7, fold 5
## [2026-08-16 13:04:54] axes2 | total_without_group | outer repeat 8, fold 1
## [2026-08-16 13:04:55] axes2 | total_without_group | outer repeat 8, fold 2
## [2026-08-16 13:04:56] axes2 | total_without_group | outer repeat 8, fold 3
## [2026-08-16 13:04:57] axes2 | total_without_group | outer repeat 8, fold 4
## [2026-08-16 13:04:59] axes2 | total_without_group | outer repeat 8, fold 5
## [2026-08-16 13:05:00] axes2 | total_without_group | outer repeat 9, fold 1
## [2026-08-16 13:05:01] axes2 | total_without_group | outer repeat 9, fold 2
## [2026-08-16 13:05:02] axes2 | total_without_group | outer repeat 9, fold 3
## [2026-08-16 13:05:03] axes2 | total_without_group | outer repeat 9, fold 4
## [2026-08-16 13:05:04] axes2 | total_without_group | outer repeat 9, fold 5
## [2026-08-16 13:05:05] axes2 | total_without_group | outer repeat 10, fold 1
## [2026-08-16 13:05:06] axes2 | total_without_group | outer repeat 10, fold 2
## [2026-08-16 13:05:07] axes2 | total_without_group | outer repeat 10, fold 3
## [2026-08-16 13:05:08] axes2 | total_without_group | outer repeat 10, fold 4
## [2026-08-16 13:05:09] axes2 | total_without_group | outer repeat 10, fold 5
## [2026-08-16 13:05:10] axes2 | incremental_beyond_group | outer repeat 1, fold 1
## [2026-08-16 13:05:11] axes2 | incremental_beyond_group | outer repeat 1, fold 2
## [2026-08-16 13:05:12] axes2 | incremental_beyond_group | outer repeat 1, fold 3
## [2026-08-16 13:05:13] axes2 | incremental_beyond_group | outer repeat 1, fold 4
## [2026-08-16 13:05:14] axes2 | incremental_beyond_group | outer repeat 1, fold 5
## [2026-08-16 13:05:15] axes2 | incremental_beyond_group | outer repeat 2, fold 1
## [2026-08-16 13:05:17] axes2 | incremental_beyond_group | outer repeat 2, fold 2
## [2026-08-16 13:05:18] axes2 | incremental_beyond_group | outer repeat 2, fold 3
## [2026-08-16 13:05:19] axes2 | incremental_beyond_group | outer repeat 2, fold 4
## [2026-08-16 13:05:20] axes2 | incremental_beyond_group | outer repeat 2, fold 5
## [2026-08-16 13:05:21] axes2 | incremental_beyond_group | outer repeat 3, fold 1
## [2026-08-16 13:05:22] axes2 | incremental_beyond_group | outer repeat 3, fold 2
## [2026-08-16 13:05:23] axes2 | incremental_beyond_group | outer repeat 3, fold 3
## [2026-08-16 13:05:24] axes2 | incremental_beyond_group | outer repeat 3, fold 4
## [2026-08-16 13:05:25] axes2 | incremental_beyond_group | outer repeat 3, fold 5
## [2026-08-16 13:05:26] axes2 | incremental_beyond_group | outer repeat 4, fold 1
## [2026-08-16 13:05:27] axes2 | incremental_beyond_group | outer repeat 4, fold 2
## [2026-08-16 13:05:28] axes2 | incremental_beyond_group | outer repeat 4, fold 3
## [2026-08-16 13:05:29] axes2 | incremental_beyond_group | outer repeat 4, fold 4
## [2026-08-16 13:05:30] axes2 | incremental_beyond_group | outer repeat 4, fold 5
## [2026-08-16 13:05:31] axes2 | incremental_beyond_group | outer repeat 5, fold 1
## [2026-08-16 13:05:32] axes2 | incremental_beyond_group | outer repeat 5, fold 2
## [2026-08-16 13:05:33] axes2 | incremental_beyond_group | outer repeat 5, fold 3
## [2026-08-16 13:05:34] axes2 | incremental_beyond_group | outer repeat 5, fold 4
## [2026-08-16 13:05:36] axes2 | incremental_beyond_group | outer repeat 5, fold 5
## [2026-08-16 13:05:37] axes2 | incremental_beyond_group | outer repeat 6, fold 1
## [2026-08-16 13:05:38] axes2 | incremental_beyond_group | outer repeat 6, fold 2
## [2026-08-16 13:05:39] axes2 | incremental_beyond_group | outer repeat 6, fold 3
## [2026-08-16 13:05:40] axes2 | incremental_beyond_group | outer repeat 6, fold 4
## [2026-08-16 13:05:41] axes2 | incremental_beyond_group | outer repeat 6, fold 5
## [2026-08-16 13:05:42] axes2 | incremental_beyond_group | outer repeat 7, fold 1
## [2026-08-16 13:05:43] axes2 | incremental_beyond_group | outer repeat 7, fold 2
## [2026-08-16 13:05:44] axes2 | incremental_beyond_group | outer repeat 7, fold 3
## [2026-08-16 13:05:45] axes2 | incremental_beyond_group | outer repeat 7, fold 4
## [2026-08-16 13:05:46] axes2 | incremental_beyond_group | outer repeat 7, fold 5
## [2026-08-16 13:05:47] axes2 | incremental_beyond_group | outer repeat 8, fold 1
## [2026-08-16 13:05:48] axes2 | incremental_beyond_group | outer repeat 8, fold 2
## [2026-08-16 13:05:49] axes2 | incremental_beyond_group | outer repeat 8, fold 3
## [2026-08-16 13:05:50] axes2 | incremental_beyond_group | outer repeat 8, fold 4
## [2026-08-16 13:05:51] axes2 | incremental_beyond_group | outer repeat 8, fold 5
## [2026-08-16 13:05:52] axes2 | incremental_beyond_group | outer repeat 9, fold 1
## [2026-08-16 13:05:53] axes2 | incremental_beyond_group | outer repeat 9, fold 2
## [2026-08-16 13:05:55] axes2 | incremental_beyond_group | outer repeat 9, fold 3
## [2026-08-16 13:05:56] axes2 | incremental_beyond_group | outer repeat 9, fold 4
## [2026-08-16 13:05:57] axes2 | incremental_beyond_group | outer repeat 9, fold 5
## [2026-08-16 13:05:58] axes2 | incremental_beyond_group | outer repeat 10, fold 1
## [2026-08-16 13:05:59] axes2 | incremental_beyond_group | outer repeat 10, fold 2
## [2026-08-16 13:06:00] axes2 | incremental_beyond_group | outer repeat 10, fold 3
## [2026-08-16 13:06:01] axes2 | incremental_beyond_group | outer repeat 10, fold 4
## [2026-08-16 13:06:02] axes2 | incremental_beyond_group | outer repeat 10, fold 5
## [2026-08-16 13:06:03] targeted_edges22 | total_without_group | outer repeat 1, fold 1
## [2026-08-16 13:06:04] targeted_edges22 | total_without_group | outer repeat 1, fold 2
## [2026-08-16 13:06:05] targeted_edges22 | total_without_group | outer repeat 1, fold 3
## [2026-08-16 13:06:07] targeted_edges22 | total_without_group | outer repeat 1, fold 4
## [2026-08-16 13:06:08] targeted_edges22 | total_without_group | outer repeat 1, fold 5
## [2026-08-16 13:06:09] targeted_edges22 | total_without_group | outer repeat 2, fold 1
## [2026-08-16 13:06:11] targeted_edges22 | total_without_group | outer repeat 2, fold 2
## [2026-08-16 13:06:12] targeted_edges22 | total_without_group | outer repeat 2, fold 3
## [2026-08-16 13:06:13] targeted_edges22 | total_without_group | outer repeat 2, fold 4
## [2026-08-16 13:06:14] targeted_edges22 | total_without_group | outer repeat 2, fold 5
## [2026-08-16 13:06:16] targeted_edges22 | total_without_group | outer repeat 3, fold 1
## [2026-08-16 13:06:17] targeted_edges22 | total_without_group | outer repeat 3, fold 2
## [2026-08-16 13:06:18] targeted_edges22 | total_without_group | outer repeat 3, fold 3
## [2026-08-16 13:06:20] targeted_edges22 | total_without_group | outer repeat 3, fold 4
## [2026-08-16 13:06:21] targeted_edges22 | total_without_group | outer repeat 3, fold 5
## [2026-08-16 13:06:22] targeted_edges22 | total_without_group | outer repeat 4, fold 1
## [2026-08-16 13:06:24] targeted_edges22 | total_without_group | outer repeat 4, fold 2
## [2026-08-16 13:06:25] targeted_edges22 | total_without_group | outer repeat 4, fold 3
## [2026-08-16 13:06:26] targeted_edges22 | total_without_group | outer repeat 4, fold 4
## [2026-08-16 13:06:27] targeted_edges22 | total_without_group | outer repeat 4, fold 5
## [2026-08-16 13:06:29] targeted_edges22 | total_without_group | outer repeat 5, fold 1
## [2026-08-16 13:06:30] targeted_edges22 | total_without_group | outer repeat 5, fold 2
## [2026-08-16 13:06:31] targeted_edges22 | total_without_group | outer repeat 5, fold 3
## [2026-08-16 13:06:32] targeted_edges22 | total_without_group | outer repeat 5, fold 4
## [2026-08-16 13:06:34] targeted_edges22 | total_without_group | outer repeat 5, fold 5
## [2026-08-16 13:06:35] targeted_edges22 | total_without_group | outer repeat 6, fold 1
## [2026-08-16 13:06:36] targeted_edges22 | total_without_group | outer repeat 6, fold 2
## [2026-08-16 13:06:38] targeted_edges22 | total_without_group | outer repeat 6, fold 3
## [2026-08-16 13:06:39] targeted_edges22 | total_without_group | outer repeat 6, fold 4
## [2026-08-16 13:06:40] targeted_edges22 | total_without_group | outer repeat 6, fold 5
## [2026-08-16 13:06:42] targeted_edges22 | total_without_group | outer repeat 7, fold 1
## [2026-08-16 13:06:43] targeted_edges22 | total_without_group | outer repeat 7, fold 2
## [2026-08-16 13:06:44] targeted_edges22 | total_without_group | outer repeat 7, fold 3
## [2026-08-16 13:06:45] targeted_edges22 | total_without_group | outer repeat 7, fold 4
## [2026-08-16 13:06:47] targeted_edges22 | total_without_group | outer repeat 7, fold 5
## [2026-08-16 13:06:48] targeted_edges22 | total_without_group | outer repeat 8, fold 1
## [2026-08-16 13:06:49] targeted_edges22 | total_without_group | outer repeat 8, fold 2
## [2026-08-16 13:06:51] targeted_edges22 | total_without_group | outer repeat 8, fold 3
## [2026-08-16 13:06:52] targeted_edges22 | total_without_group | outer repeat 8, fold 4
## [2026-08-16 13:06:53] targeted_edges22 | total_without_group | outer repeat 8, fold 5
## [2026-08-16 13:06:54] targeted_edges22 | total_without_group | outer repeat 9, fold 1
## [2026-08-16 13:06:56] targeted_edges22 | total_without_group | outer repeat 9, fold 2
## [2026-08-16 13:06:57] targeted_edges22 | total_without_group | outer repeat 9, fold 3
## [2026-08-16 13:06:58] targeted_edges22 | total_without_group | outer repeat 9, fold 4
## [2026-08-16 13:07:00] targeted_edges22 | total_without_group | outer repeat 9, fold 5
## [2026-08-16 13:07:01] targeted_edges22 | total_without_group | outer repeat 10, fold 1
## [2026-08-16 13:07:02] targeted_edges22 | total_without_group | outer repeat 10, fold 2
## [2026-08-16 13:07:04] targeted_edges22 | total_without_group | outer repeat 10, fold 3
## [2026-08-16 13:07:05] targeted_edges22 | total_without_group | outer repeat 10, fold 4
## [2026-08-16 13:07:06] targeted_edges22 | total_without_group | outer repeat 10, fold 5
## [2026-08-16 13:07:07] targeted_edges22 | incremental_beyond_group | outer repeat 1, fold 1
## [2026-08-16 13:07:09] targeted_edges22 | incremental_beyond_group | outer repeat 1, fold 2
## [2026-08-16 13:07:10] targeted_edges22 | incremental_beyond_group | outer repeat 1, fold 3
## [2026-08-16 13:07:11] targeted_edges22 | incremental_beyond_group | outer repeat 1, fold 4
## [2026-08-16 13:07:13] targeted_edges22 | incremental_beyond_group | outer repeat 1, fold 5
## [2026-08-16 13:07:14] targeted_edges22 | incremental_beyond_group | outer repeat 2, fold 1
## [2026-08-16 13:07:15] targeted_edges22 | incremental_beyond_group | outer repeat 2, fold 2
## [2026-08-16 13:07:16] targeted_edges22 | incremental_beyond_group | outer repeat 2, fold 3
## [2026-08-16 13:07:18] targeted_edges22 | incremental_beyond_group | outer repeat 2, fold 4
## [2026-08-16 13:07:19] targeted_edges22 | incremental_beyond_group | outer repeat 2, fold 5
## [2026-08-16 13:07:20] targeted_edges22 | incremental_beyond_group | outer repeat 3, fold 1
## [2026-08-16 13:07:21] targeted_edges22 | incremental_beyond_group | outer repeat 3, fold 2
## [2026-08-16 13:07:23] targeted_edges22 | incremental_beyond_group | outer repeat 3, fold 3
## [2026-08-16 13:07:24] targeted_edges22 | incremental_beyond_group | outer repeat 3, fold 4
## [2026-08-16 13:07:25] targeted_edges22 | incremental_beyond_group | outer repeat 3, fold 5
## [2026-08-16 13:07:27] targeted_edges22 | incremental_beyond_group | outer repeat 4, fold 1
## [2026-08-16 13:07:28] targeted_edges22 | incremental_beyond_group | outer repeat 4, fold 2
## [2026-08-16 13:07:29] targeted_edges22 | incremental_beyond_group | outer repeat 4, fold 3
## [2026-08-16 13:07:31] targeted_edges22 | incremental_beyond_group | outer repeat 4, fold 4
## [2026-08-16 13:07:32] targeted_edges22 | incremental_beyond_group | outer repeat 4, fold 5
## [2026-08-16 13:07:33] targeted_edges22 | incremental_beyond_group | outer repeat 5, fold 1
## [2026-08-16 13:07:34] targeted_edges22 | incremental_beyond_group | outer repeat 5, fold 2
## [2026-08-16 13:07:36] targeted_edges22 | incremental_beyond_group | outer repeat 5, fold 3
## [2026-08-16 13:07:37] targeted_edges22 | incremental_beyond_group | outer repeat 5, fold 4
## [2026-08-16 13:07:38] targeted_edges22 | incremental_beyond_group | outer repeat 5, fold 5
## [2026-08-16 13:07:39] targeted_edges22 | incremental_beyond_group | outer repeat 6, fold 1
## [2026-08-16 13:07:41] targeted_edges22 | incremental_beyond_group | outer repeat 6, fold 2
## [2026-08-16 13:07:42] targeted_edges22 | incremental_beyond_group | outer repeat 6, fold 3
## [2026-08-16 13:07:43] targeted_edges22 | incremental_beyond_group | outer repeat 6, fold 4
## [2026-08-16 13:07:45] targeted_edges22 | incremental_beyond_group | outer repeat 6, fold 5
## [2026-08-16 13:07:46] targeted_edges22 | incremental_beyond_group | outer repeat 7, fold 1
## [2026-08-16 13:07:47] targeted_edges22 | incremental_beyond_group | outer repeat 7, fold 2
## [2026-08-16 13:07:49] targeted_edges22 | incremental_beyond_group | outer repeat 7, fold 3
## [2026-08-16 13:07:50] targeted_edges22 | incremental_beyond_group | outer repeat 7, fold 4
## [2026-08-16 13:07:51] targeted_edges22 | incremental_beyond_group | outer repeat 7, fold 5
## [2026-08-16 13:07:52] targeted_edges22 | incremental_beyond_group | outer repeat 8, fold 1
## [2026-08-16 13:07:54] targeted_edges22 | incremental_beyond_group | outer repeat 8, fold 2
## [2026-08-16 13:07:55] targeted_edges22 | incremental_beyond_group | outer repeat 8, fold 3
## [2026-08-16 13:07:56] targeted_edges22 | incremental_beyond_group | outer repeat 8, fold 4
## [2026-08-16 13:07:57] targeted_edges22 | incremental_beyond_group | outer repeat 8, fold 5
## [2026-08-16 13:07:59] targeted_edges22 | incremental_beyond_group | outer repeat 9, fold 1
## [2026-08-16 13:08:00] targeted_edges22 | incremental_beyond_group | outer repeat 9, fold 2
## [2026-08-16 13:08:01] targeted_edges22 | incremental_beyond_group | outer repeat 9, fold 3
## [2026-08-16 13:08:03] targeted_edges22 | incremental_beyond_group | outer repeat 9, fold 4
## [2026-08-16 13:08:04] targeted_edges22 | incremental_beyond_group | outer repeat 9, fold 5
## [2026-08-16 13:08:05] targeted_edges22 | incremental_beyond_group | outer repeat 10, fold 1
## [2026-08-16 13:08:06] targeted_edges22 | incremental_beyond_group | outer repeat 10, fold 2
## [2026-08-16 13:08:08] targeted_edges22 | incremental_beyond_group | outer repeat 10, fold 3
## [2026-08-16 13:08:09] targeted_edges22 | incremental_beyond_group | outer repeat 10, fold 4
## [2026-08-16 13:08:10] targeted_edges22 | incremental_beyond_group | outer repeat 10, fold 5all_predictions <- rbindlist(lapply(mvpa_runs, `[[`, "predictions"), fill = TRUE)
all_hyperparameters <- rbindlist(lapply(mvpa_runs, `[[`, "hyperparameters"), fill = TRUE)
all_weights <- rbindlist(lapply(mvpa_runs, `[[`, "weights"), fill = TRUE)
if (nrow(all_predictions) == 0L) {
stop(
"No valid nested-CV predictions were produced. This usually means the matched sample or group balance is too small for the planned outer/inner folds."
)
}
fwrite(
all_predictions,
file.path(output_dir, "nested_cv_all_outer_predictions_long.csv")
)
fwrite(
all_hyperparameters,
file.path(output_dir, "nested_cv_selected_hyperparameters.csv")
)
fwrite(
all_weights,
file.path(output_dir, "nested_cv_selected_model_weights_long.csv")
)10 Out-of-fold performance
prediction_average <- average_oof_predictions(all_predictions)
fwrite(
prediction_average,
file.path(output_dir, "nested_cv_subject_averaged_predictions_long.csv")
)
coefficient_metrics <- coefficient_metric_tables(prediction_average)
fwrite(
coefficient_metrics$by_outcome,
file.path(output_dir, "nested_cv_coefficient_metrics_by_outcome.csv")
)
fwrite(
coefficient_metrics$overall,
file.path(output_dir, "nested_cv_coefficient_metrics_overall.csv")
)
print(coefficient_metrics$overall)## feature_set analysis_mode prediction_source n rmse
## <char> <char> <fctr> <int> <num>
## 1: axes2 total_without_group baseline_prediction 2472 1.020057
## 2: axes2 incremental_beyond_group baseline_prediction 2472 1.023856
## 3: targeted_edges22 total_without_group baseline_prediction 2472 1.020057
## 4: targeted_edges22 incremental_beyond_group baseline_prediction 2472 1.023856
## 5: axes2 total_without_group full_prediction 2472 1.020057
## 6: axes2 incremental_beyond_group full_prediction 2472 1.023856
## 7: targeted_edges22 total_without_group full_prediction 2472 1.020057
## 8: targeted_edges22 incremental_beyond_group full_prediction 2472 1.023856
## mae q2 correlation
## <num> <num> <num>
## 1: 0.7454293 0.6887832 0.8300115
## 2: 0.7496451 0.6864611 0.8286699
## 3: 0.7454293 0.6887832 0.8300115
## 4: 0.7496451 0.6864611 0.8286699
## 5: 0.7454293 0.6887832 0.8300115
## 6: 0.7496451 0.6864611 0.8286699
## 7: 0.7454293 0.6887832 0.8300115
## 8: 0.7496451 0.6864611 0.8286699curve_predictions <- reconstruct_curve_predictions(
prediction_average = prediction_average,
basis_grid = basis_grid,
conditions = primary_eye_conditions,
basis_df = eye_basis_df
)
fwrite(
curve_predictions,
file.path(output_dir, "nested_cv_reconstructed_curve_predictions.csv")
)
curve_metrics <- curve_metric_tables(curve_predictions)
fwrite(
curve_metrics$by_condition,
file.path(output_dir, "nested_cv_curve_metrics_by_condition.csv")
)
fwrite(
curve_metrics$overall,
file.path(output_dir, "nested_cv_curve_metrics_overall.csv")
)
print(curve_metrics$by_condition)## feature_set analysis_mode condition prediction_source
## <char> <char> <char> <fctr>
## 1: axes2 total_without_group Fear__C baseline_probability
## 2: axes2 total_without_group Sad__C baseline_probability
## 3: axes2 total_without_group Sad__E baseline_probability
## 4: axes2 incremental_beyond_group Fear__C baseline_probability
## 5: axes2 incremental_beyond_group Sad__C baseline_probability
## 6: axes2 incremental_beyond_group Sad__E baseline_probability
## 7: targeted_edges22 total_without_group Fear__C baseline_probability
## 8: targeted_edges22 total_without_group Sad__C baseline_probability
## 9: targeted_edges22 total_without_group Sad__E baseline_probability
## 10: targeted_edges22 incremental_beyond_group Fear__C baseline_probability
## 11: targeted_edges22 incremental_beyond_group Sad__C baseline_probability
## 12: targeted_edges22 incremental_beyond_group Sad__E baseline_probability
## 13: axes2 total_without_group Fear__C full_probability
## 14: axes2 total_without_group Sad__C full_probability
## 15: axes2 total_without_group Sad__E full_probability
## 16: axes2 incremental_beyond_group Fear__C full_probability
## 17: axes2 incremental_beyond_group Sad__C full_probability
## 18: axes2 incremental_beyond_group Sad__E full_probability
## 19: targeted_edges22 total_without_group Fear__C full_probability
## 20: targeted_edges22 total_without_group Sad__C full_probability
## 21: targeted_edges22 total_without_group Sad__E full_probability
## 22: targeted_edges22 incremental_beyond_group Fear__C full_probability
## 23: targeted_edges22 incremental_beyond_group Sad__C full_probability
## 24: targeted_edges22 incremental_beyond_group Sad__E full_probability
## feature_set analysis_mode condition prediction_source
## n rmse mae q2 correlation
## <int> <num> <num> <num> <num>
## 1: 10300 0.2150084 0.1612299 -0.06198417 0.21112062
## 2: 10300 0.2045974 0.1472725 -0.13590913 0.09881540
## 3: 10300 0.1820177 0.1367654 0.57433762 0.76522265
## 4: 10300 0.2151015 0.1614665 -0.06290443 0.21417921
## 5: 10300 0.2069233 0.1486674 -0.16188263 0.07831281
## 6: 10300 0.1818246 0.1365719 0.57524030 0.76611387
## 7: 10300 0.2150084 0.1612299 -0.06198417 0.21112062
## 8: 10300 0.2045974 0.1472725 -0.13590913 0.09881540
## 9: 10300 0.1820177 0.1367654 0.57433762 0.76522265
## 10: 10300 0.2151015 0.1614665 -0.06290443 0.21417921
## 11: 10300 0.2069233 0.1486674 -0.16188263 0.07831281
## 12: 10300 0.1818246 0.1365719 0.57524030 0.76611387
## 13: 10300 0.2150084 0.1612299 -0.06198417 0.21112062
## 14: 10300 0.2045974 0.1472725 -0.13590913 0.09881540
## 15: 10300 0.1820177 0.1367654 0.57433762 0.76522265
## 16: 10300 0.2151015 0.1614665 -0.06290443 0.21417921
## 17: 10300 0.2069233 0.1486674 -0.16188263 0.07831281
## 18: 10300 0.1818246 0.1365719 0.57524030 0.76611387
## 19: 10300 0.2150084 0.1612299 -0.06198417 0.21112062
## 20: 10300 0.2045974 0.1472725 -0.13590913 0.09881540
## 21: 10300 0.1820177 0.1367654 0.57433762 0.76522265
## 22: 10300 0.2151015 0.1614665 -0.06290443 0.21417921
## 23: 10300 0.2069233 0.1486674 -0.16188263 0.07831281
## 24: 10300 0.1818246 0.1365719 0.57524030 0.76611387
## n rmse mae q2 correlation11 Can fMRI reproduce the held-out MDD-HC temporal difference?
These are matched-sample reconstructions from the subject-level basis phenotypes. They do not replace the original trial-level GAMM or its simultaneous difference bands. Their purpose is to show whether predictions made for held-out participants recover the shape of the MDD-HC contrast.
group_difference <- make_group_difference_table(curve_predictions)
fwrite(
group_difference,
file.path(output_dir, "observed_and_predicted_mdd_hc_difference_curves.csv")
)
group_difference_long <- melt(
group_difference,
id.vars = c(
"feature_set", "analysis_mode", "condition", "time_ms",
"observed_difference"
),
measure.vars = c("baseline_difference", "full_difference"),
variable.name = "prediction_source",
value.name = "predicted_difference"
)
group_difference_metrics <- group_difference_long[,
calculate_basic_metrics(observed_difference, predicted_difference),
by = .(feature_set, analysis_mode, condition, prediction_source)
]
fwrite(
group_difference_metrics,
file.path(output_dir, "mdd_hc_difference_curve_metrics.csv")
)
print(group_difference_metrics)## feature_set analysis_mode condition prediction_source
## <char> <char> <char> <fctr>
## 1: axes2 incremental_beyond_group Fear__C baseline_difference
## 2: axes2 incremental_beyond_group Sad__C baseline_difference
## 3: axes2 incremental_beyond_group Sad__E baseline_difference
## 4: axes2 total_without_group Fear__C baseline_difference
## 5: axes2 total_without_group Sad__C baseline_difference
## 6: axes2 total_without_group Sad__E baseline_difference
## 7: targeted_edges22 incremental_beyond_group Fear__C baseline_difference
## 8: targeted_edges22 incremental_beyond_group Sad__C baseline_difference
## 9: targeted_edges22 incremental_beyond_group Sad__E baseline_difference
## 10: targeted_edges22 total_without_group Fear__C baseline_difference
## 11: targeted_edges22 total_without_group Sad__C baseline_difference
## 12: targeted_edges22 total_without_group Sad__E baseline_difference
## 13: axes2 incremental_beyond_group Fear__C full_difference
## 14: axes2 incremental_beyond_group Sad__C full_difference
## 15: axes2 incremental_beyond_group Sad__E full_difference
## 16: axes2 total_without_group Fear__C full_difference
## 17: axes2 total_without_group Sad__C full_difference
## 18: axes2 total_without_group Sad__E full_difference
## 19: targeted_edges22 incremental_beyond_group Fear__C full_difference
## 20: targeted_edges22 incremental_beyond_group Sad__C full_difference
## 21: targeted_edges22 incremental_beyond_group Sad__E full_difference
## 22: targeted_edges22 total_without_group Fear__C full_difference
## 23: targeted_edges22 total_without_group Sad__C full_difference
## 24: targeted_edges22 total_without_group Sad__E full_difference
## feature_set analysis_mode condition prediction_source
## n rmse mae q2 correlation
## <int> <num> <num> <num> <num>
## 1: 100 0.03226227 0.01444331 0.48768999 0.9199588
## 2: 100 0.01804955 0.01448516 0.67421660 0.9181206
## 3: 100 0.01350200 0.01094395 0.91305136 0.9863743
## 4: 100 0.04372680 0.03905330 0.05889427 0.2500319
## 5: 100 0.02410158 0.01887337 0.41911916 0.6957355
## 6: 100 0.05147204 0.03876492 -0.26359970 -0.7359776
## 7: 100 0.03226227 0.01444331 0.48768999 0.9199588
## 8: 100 0.01804955 0.01448516 0.67421660 0.9181206
## 9: 100 0.01350200 0.01094395 0.91305136 0.9863743
## 10: 100 0.04372680 0.03905330 0.05889427 0.2500319
## 11: 100 0.02410158 0.01887337 0.41911916 0.6957355
## 12: 100 0.05147204 0.03876492 -0.26359970 -0.7359776
## 13: 100 0.03226227 0.01444331 0.48768999 0.9199588
## 14: 100 0.01804955 0.01448516 0.67421660 0.9181206
## 15: 100 0.01350200 0.01094395 0.91305136 0.9863743
## 16: 100 0.04372680 0.03905330 0.05889427 0.2500319
## 17: 100 0.02410158 0.01887337 0.41911916 0.6957355
## 18: 100 0.05147204 0.03876492 -0.26359970 -0.7359776
## 19: 100 0.03226227 0.01444331 0.48768999 0.9199588
## 20: 100 0.01804955 0.01448516 0.67421660 0.9181206
## 21: 100 0.01350200 0.01094395 0.91305136 0.9863743
## 22: 100 0.04372680 0.03905330 0.05889427 0.2500319
## 23: 100 0.02410158 0.01887337 0.41911916 0.6957355
## 24: 100 0.05147204 0.03876492 -0.26359970 -0.7359776
## n rmse mae q2 correlationfor (mode_value in analysis_modes_to_run) {
plot_data <- group_difference[analysis_mode == mode_value]
plot_long <- melt(
plot_data,
id.vars = c("feature_set", "analysis_mode", "condition", "time_ms"),
measure.vars = c(
"observed_difference", "baseline_difference", "full_difference"
),
variable.name = "series",
value.name = "mdd_minus_hc_probability"
)
p_difference <- ggplot(
plot_long,
aes(
x = time_ms,
y = mdd_minus_hc_probability,
colour = series,
linetype = series
)
) +
geom_hline(yintercept = 0, linetype = 3) +
geom_line(linewidth = 0.8) +
facet_grid(condition ~ feature_set) +
theme_classic() +
labs(
title = paste("Held-out MDD-HC temporal difference:", mode_value),
x = "Time after stimulus onset (ms)",
y = "MDD minus HC emotional-face probability",
colour = "Curve",
linetype = "Curve"
)
print(p_difference)
ggsave(
file.path(
output_dir,
paste0("mdd_hc_difference_", mode_value, ".png")
),
p_difference,
width = 13,
height = 9,
dpi = 300
)
}
12 Incremental prediction test and uncertainty
The statistic below is the reduction in standardized out-of-fold mean squared error when the fMRI increment is added to the covariate-only prediction. A positive value favors the fMRI model. The permutation keeps the already-created out-of-fold predictions fixed and shuffles the multivariate residual phenotype at the participant level. This is a useful prediction-alignment test, but a full-pipeline permutation that repeats every nested-CV fit would be stronger and much more computationally expensive.
prediction_tests <- prediction_alignment_test(
prediction_average = prediction_average,
n_permutations = prediction_permutations,
n_bootstrap = prediction_bootstraps,
seed = random_seed
)
prediction_tests[, p_bh_within_mode := p.adjust(
p_fixed_prediction_permutation,
method = "BH"
), by = .(analysis_mode)]
print(prediction_tests)## feature_set analysis_mode condition delta_mse
## <char> <char> <char> <num>
## 1: axes2 total_without_group ALL_CONDITIONS 0
## 2: axes2 total_without_group Fear__C 0
## 3: axes2 total_without_group Sad__C 0
## 4: axes2 total_without_group Sad__E 0
## 5: axes2 incremental_beyond_group ALL_CONDITIONS 0
## 6: axes2 incremental_beyond_group Fear__C 0
## 7: axes2 incremental_beyond_group Sad__C 0
## 8: axes2 incremental_beyond_group Sad__E 0
## 9: targeted_edges22 total_without_group ALL_CONDITIONS 0
## 10: targeted_edges22 total_without_group Fear__C 0
## 11: targeted_edges22 total_without_group Sad__C 0
## 12: targeted_edges22 total_without_group Sad__E 0
## 13: targeted_edges22 incremental_beyond_group ALL_CONDITIONS 0
## 14: targeted_edges22 incremental_beyond_group Fear__C 0
## 15: targeted_edges22 incremental_beyond_group Sad__C 0
## 16: targeted_edges22 incremental_beyond_group Sad__E 0
## proportional_mse_reduction bootstrap_ci_low bootstrap_ci_high
## <num> <num> <num>
## 1: 0 -1.110223e-16 1.110223e-16
## 2: 0 -2.220446e-16 0.000000e+00
## 3: 0 -1.110223e-16 2.220446e-16
## 4: 0 -1.110223e-16 2.220446e-16
## 5: 0 -2.220446e-16 0.000000e+00
## 6: 0 -2.220446e-16 0.000000e+00
## 7: 0 -1.110223e-16 2.220446e-16
## 8: 0 -2.220446e-16 0.000000e+00
## 9: 0 -1.110223e-16 2.220446e-16
## 10: 0 -2.220446e-16 0.000000e+00
## 11: 0 -1.110223e-16 2.220446e-16
## 12: 0 -1.110223e-16 2.220446e-16
## 13: 0 -2.220446e-16 0.000000e+00
## 14: 0 -2.220446e-16 0.000000e+00
## 15: 0 -1.110223e-16 2.220446e-16
## 16: 0 -2.220446e-16 0.000000e+00
## p_fixed_prediction_permutation n_subjects n_outcomes n_permutations
## <num> <int> <int> <int>
## 1: 0.5234953 103 24 5000
## 2: 1.0000000 103 8 5000
## 3: 0.8614277 103 8 5000
## 4: 1.0000000 103 8 5000
## 5: 0.9958008 103 24 5000
## 6: 1.0000000 103 8 5000
## 7: 1.0000000 103 8 5000
## 8: 1.0000000 103 8 5000
## 9: 0.5254949 103 24 5000
## 10: 1.0000000 103 8 5000
## 11: 0.8684263 103 8 5000
## 12: 1.0000000 103 8 5000
## 13: 0.9980004 103 24 5000
## 14: 1.0000000 103 8 5000
## 15: 1.0000000 103 8 5000
## 16: 1.0000000 103 8 5000
## n_bootstrap p_bh_within_mode
## <int> <num>
## 1: 2000 1
## 2: 2000 1
## 3: 2000 1
## 4: 2000 1
## 5: 2000 1
## 6: 2000 1
## 7: 2000 1
## 8: 2000 1
## 9: 2000 1
## 10: 2000 1
## 11: 2000 1
## 12: 2000 1
## 13: 2000 1
## 14: 2000 1
## 15: 2000 1
## 16: 2000 1fwrite(
prediction_tests,
file.path(output_dir, "incremental_prediction_permutation_and_bootstrap.csv")
)13 Feature stability
Predictive weights are conditional multivariate weights. They are not isolated edge effects and should not be interpreted as causal or independently significant connections.
feature_stability <- all_weights[, .(
selection_frequency = mean(abs(weight) > 1e-10),
mean_l2_norm = mean(feature_l2_norm),
median_l2_norm = median(feature_l2_norm),
l2_norm_q25 = quantile(feature_l2_norm, 0.25),
l2_norm_q75 = quantile(feature_l2_norm, 0.75)
), by = .(feature_set, analysis_mode, feature)]
edge_labels <- unique(manifest[, .(
feature = FC,
edge_label = if (all(c("node1_label", "node2_label") %in% names(manifest))) {
paste(node1_label, node2_label, sep = " -- ")
} else {
FC
},
axis,
subfamily
)])
feature_stability <- merge(
feature_stability,
edge_labels,
by = "feature",
all.x = TRUE,
sort = FALSE
)
feature_stability[is.na(edge_label), edge_label := feature]
fwrite(
feature_stability,
file.path(output_dir, "nested_cv_feature_stability.csv")
)
plot_edges <- feature_stability[
feature_set == "targeted_edges22" &
analysis_mode == "incremental_beyond_group"
]
setorder(plot_edges, -mean_l2_norm)
plot_edges <- head(plot_edges, 15L)
plot_edges[, edge_label := factor(
edge_label,
levels = unique(rev(as.character(edge_label)))
)]
if (nrow(plot_edges) > 0L) {
p_weights <- ggplot(
plot_edges,
aes(x = mean_l2_norm, y = edge_label)
) +
geom_point() +
geom_segment(
aes(
x = l2_norm_q25,
xend = l2_norm_q75,
y = edge_label,
yend = edge_label
)
) +
theme_classic() +
labs(
title = "Targeted-edge weight stability beyond diagnosis",
x = "Mean L2 norm across outer models",
y = NULL
)
print(p_weights)
ggsave(
file.path(output_dir, "targeted_edge_weight_stability.png"),
p_weights,
width = 11,
height = 7,
dpi = 300
)
}
14 Interpretation guide
Use the following hierarchy when writing the manuscript:
- Matched-sample requirement. Report the exact number of participants with both valid fMRI and eligible eye trajectories. Do not combine the larger modality-specific sample sizes as though all participants contributed to the cross-modal analysis.
- Primary cross-modal evidence. Emphasize the
two-axis multivariate test and the
axes2 / incremental_beyond_groupheld-out performance. This asks whether individual variation in the MDD-oriented connectivity axes predicts temporal gaze variation beyond diagnosis and covariates. - Exploratory pattern evidence. Treat the 22-edge elastic-net model as exploratory unless the edge manifest was fixed independently of the current group results. Report prediction, not edge-wise significance.
- Diagnostic-group ambiguity. If
total_without_grouppredicts well butincremental_beyond_groupdoes not, the data show shared diagnostic-group structure across modalities, not a participant-level connectivity-gaze link. - External validity. Repeated nested cross-validation reduces leakage but does not substitute for an independent replication sample. With a matched sample near 100, uncertainty can remain large.
- Causal language. Use “predicts,” “is associated with,” or “accounts for out-of-sample variation.” Do not use “explains the mechanism,” “mediates,” or “causes” without a separate causal design.
15 Reproducibility metadata
settings_table <- data.table(
setting = c(
"data_dir",
"eye_condition_set",
"primary_eye_conditions",
"eye_basis_df",
"minimum_eye_bins_per_condition",
"outer_folds",
"outer_repeats",
"inner_folds",
"alpha_grid",
"lambda_selection_rule",
"motion_available_fraction",
"motion_used",
"include_icv_sensitivity",
"axis_permutations",
"prediction_permutations",
"prediction_bootstraps",
"targets_independently_prespecified",
"random_seed"
),
value = c(
data_dir,
eye_condition_set,
paste(primary_eye_conditions, collapse = ";"),
eye_basis_df,
minimum_eye_bins_per_condition,
outer_folds,
outer_repeats,
inner_folds,
paste(alpha_grid, collapse = ";"),
lambda_selection_rule,
motion_fraction,
use_motion,
include_icv_sensitivity,
axis_permutations,
prediction_permutations,
prediction_bootstraps,
targets_independently_prespecified,
random_seed
)
)
fwrite(settings_table, file.path(output_dir, "cross_modal_mvpa_settings.csv"))
writeLines(
c(
paste0("Created: ", Sys.time()),
paste0("Matched participants: ", nrow(analysis_data)),
paste0("HC: ", sum(analysis_data$Group == "HC")),
paste0("MDD: ", sum(analysis_data$Group == "MDD")),
paste0("Motion used: ", use_motion),
paste0("Output directory: ", output_dir),
"",
capture.output(sessionInfo())
),
file.path(output_dir, "cross_modal_mvpa_session_log.txt")
)
sessionInfo()## R version 4.4.2 (2024-10-31 ucrt)
## Platform: x86_64-w64-mingw32/x64
## 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 data.table_1.16.2
##
## loaded via a namespace (and not attached):
## [1] Matrix_1.7-1 glmnet_4.1-10 gtable_0.3.6 jsonlite_1.8.9
## [5] dplyr_1.1.4 compiler_4.4.2 tidyselect_1.2.1 Rcpp_1.0.13-1
## [9] jquerylib_0.1.4 textshaping_0.4.0 systemfonts_1.3.1 splines_4.4.2
## [13] scales_1.3.0 yaml_2.3.10 fastmap_1.2.0 lattice_0.22-6
## [17] R6_2.5.1 labeling_0.4.3 generics_0.1.3 shape_1.4.6.1
## [21] knitr_1.49 iterators_1.0.14 tibble_3.2.1 munsell_0.5.1
## [25] bslib_0.8.0 pillar_1.9.0 rlang_1.2.0 utf8_1.2.4
## [29] cachem_1.1.0 xfun_0.49 sass_0.4.9 cli_3.6.3
## [33] withr_3.0.2 magrittr_2.0.3 digest_0.6.37 foreach_1.5.2
## [37] grid_4.4.2 rstudioapi_0.19.0 lifecycle_1.0.4 vctrs_0.6.5
## [41] evaluate_1.0.1 glue_1.8.0 farver_2.1.2 ragg_1.3.3
## [45] codetools_0.2-20 survival_3.7-0 fansi_1.0.6 colorspace_2.1-1
## [49] rmarkdown_2.31 tools_4.4.2 pkgconfig_2.0.3 htmltools_0.5.8.1