Cross-Modal MVPA: Temporal Gaze Prediction of Resting-State Connectivity in MDD

1 Purpose and inferential boundary

This report reverses the original cross-modal prediction direction. Instead of using resting-state functional connectivity to predict temporal gaze, it uses a compact participant-level representation of the eye-tracking time courses to predict the 22 targeted resting-state connectivity edges.

The workflow deliberately consumes validated upstream outputs rather than sourcing analysis scripts with executable top-level code:

1. 1_resting_state_fc_prepare_validate_updated(5).R
2. 3_resting_state_fc_two_axis_targeted_updated(5).R
3. A_new_2_eye_temporal_two_axis_fixed(2).Rmd

The predictive model is a nested-cross-validated multi-response elastic net:

• Predictors (X): temporal eye-tracking B-spline coefficients.
• Outcomes (Y): the 22 targeted resting-state FC edges.
• Primary interpretation: incremental_beyond_group, which asks whether individual temporal gaze variation predicts individual FC variation beyond diagnosis and covariates.
• Secondary interpretation: total_without_group, which allows the eye predictors to exploit diagnostic-group structure shared by both modalities.

The reverse direction is a predictive association analysis. It does not imply that gaze causes resting-state connectivity, that connectivity causes gaze, or that either modality mediates the other.

2 Settings

data_dir <- params$data_dir

prepared_rds <- file.path(
  data_dir, "rsfc_prepared", "rsfc_prepared_data.rds"
)
manifest_file <- file.path(
  data_dir, "rsfc_two_axis", "two_axis_edge_manifest_verified.csv"
)

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_eye_to_fmri_mvpa")
dir.create(output_dir, recursive = TRUE, showWarnings = FALSE)

# Keep the same eye-trajectory definition as the original direction so the two
# analyses remain comparable. Because the three surviving conditions were
# identified in the same eye sample, the all-nine rerun remains an important
# 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 empirical log odds of emotional-face versus neutral-face gaze.
eye_basis_df <- 8L
eye_basis_ridge <- 1e-4
minimum_eye_bins_per_condition <- 60L

analysis_modes_to_run <- c("total_without_group", "incremental_beyond_group")

# Keep the original nested-CV and regularization settings. Do not make the
# reverse analysis more permissive merely because the original direction was null.
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

minimum_motion_fraction <- 0.80
require_motion <- FALSE
include_icv_sensitivity <- FALSE

# This is a fixed-prediction permutation of held-out predictions, not a full
# re-fitting permutation. It is substantially cheaper but should be described
# as such in the manuscript.
prediction_permutations <- 5000L
prediction_bootstraps <- 2000L

# Keep this FALSE unless the 22-edge manifest was fixed independently of the
# current-sample group effects.
targets_independently_prespecified <- FALSE

random_seed <- 20260816L
verbose_progress <- TRUE

3 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("Required input file(s) not found:\n", paste(missing, collapse = "\n"))
  }
}

stop_if_missing(c(prepared_rds, manifest_file, eye_bins_file))

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)
  )
}

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 BOTH modalities using nuisance coefficients estimated only in the
# training fold. For the reverse analysis, X is eye tracking and Y is fMRI.
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_eye_to_fmri <- function(
    x_train,
    y_train,
    meta_train,
    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 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.")

  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)

    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 = x_train[idx_training, , drop = FALSE],
      X_test = x_train[idx_validation, , drop = FALSE],
      Y_train = y_train[idx_training, , drop = FALSE],
      Y_test = y_train[idx_validation, , drop = FALSE],
      C_train = nuisance$train,
      C_test = nuisance$test
    )

    response_scale <- apply(res$Y_train, 2L, sd)
    response_scale[!is.finite(response_scale) | response_scale <= 0] <- 1

    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 = response_scale
        )
      }
    }
  }

  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_eye_to_fmri <- function(
    analysis_data,
    predictor_cols,
    edge_cols,
    mode,
    outer_fold_table,
    use_motion,
    include_icv,
    seed
) {
  predictor_matrix <- as.matrix(analysis_data[, ..predictor_cols])
  edge_matrix <- as.matrix(analysis_data[, ..edge_cols])
  storage.mode(predictor_matrix) <- "double"
  storage.mode(edge_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(
      "eye_temporal_basis -> targeted_edges22 | ", 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(
        "Skipping outer repeat ", current_repeat_id,
        ", fold ", current_outer_fold,
        " because training groups are too small."
      )
      next
    }

    tuned <- inner_tune_eye_to_fmri(
      x_train = predictor_matrix[train_idx, , drop = FALSE],
      y_train = edge_matrix[train_idx, , drop = FALSE],
      meta_train = analysis_data[train_idx],
      mode = mode,
      use_motion = use_motion,
      include_icv = include_icv,
      seed = seed + current_repeat_id * 10007L + current_outer_fold * 101L
    )

    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 = predictor_matrix[train_idx, , drop = FALSE],
      X_test = predictor_matrix[test_idx, , drop = FALSE],
      Y_train = edge_matrix[train_idx, , drop = FALSE],
      Y_test = edge_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(edge_matrix)
    )
    predicted_full <- res$Y_test_baseline + predicted_residual

    prediction_list[[pred_counter]] <- data.table(
      analysis_mode = mode,
      repeat_id = current_repeat_id,
      outer_fold = current_outer_fold,
      row_id = rep(test_idx, times = length(edge_cols)),
      subject = rep(analysis_data$subject[test_idx], times = length(edge_cols)),
      Group = rep(as.character(analysis_data$Group[test_idx]), times = length(edge_cols)),
      edge = rep(edge_cols, each = length(test_idx)),
      actual = as.vector(edge_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(
      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 = predictor_cols,
      outcome_names = edge_cols
    )
    fold_weights[, `:=`(
      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, nuisance, res, tuned)
    gc(verbose = FALSE)
  }

  list(
    predictions = rbindlist(prediction_list, fill = TRUE),
    hyperparameters = rbindlist(hyperparameter_list, fill = TRUE),
    weights = rbindlist(weight_list, fill = TRUE)
  )
}

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(
      analysis_mode = character(0L),
      row_id = integer(0L),
      subject = character(0L),
      Group = character(0L),
      edge = 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 = .(analysis_mode, row_id, subject, Group, edge)]
}

edge_metric_tables <- function(prediction_average) {
  source_long <- melt(
    prediction_average,
    id.vars = c("analysis_mode", "row_id", "subject", "Group", "edge", "actual"),
    measure.vars = c("baseline_prediction", "full_prediction"),
    variable.name = "prediction_source",
    value.name = "prediction"
  )

  by_edge <- source_long[, calculate_basic_metrics(actual, prediction),
                         by = .(analysis_mode, edge, prediction_source)]

  # Overall multivariate metric after equalizing edge scales.
  source_long[, edge_sd := sd(actual), by = .(analysis_mode, edge, prediction_source)]
  source_long[!is.finite(edge_sd) | edge_sd <= 0, edge_sd := 1]
  source_long[, `:=`(
    actual_scaled = actual / edge_sd,
    prediction_scaled = prediction / edge_sd
  )]
  overall <- source_long[, calculate_basic_metrics(actual_scaled, prediction_scaled),
                         by = .(analysis_mode, prediction_source)]

  list(by_edge = by_edge, overall = overall)
}

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
}

prediction_alignment_test_edges <- function(
    prediction_average,
    n_permutations,
    n_bootstrap,
    seed
) {
  modes <- unique(prediction_average$analysis_mode)
  output <- list()
  counter <- 1L

  for (md in modes) {
    dat_mode <- prediction_average[analysis_mode == md]
    targets <- c("ALL_EDGES", unique(dat_mode$edge))

    for (target_edge in targets) {
      use <- if (target_edge == "ALL_EDGES") dat_mode else dat_mode[edge == target_edge]

      actual_wide <- dcast(use, row_id + Group ~ edge, value.var = "actual")
      baseline_wide <- dcast(
        use, row_id + Group ~ edge, value.var = "baseline_prediction"
      )
      full_wide <- dcast(use, row_id + Group ~ edge, value.var = "full_prediction")
      setorder(actual_wide, row_id)
      setorder(baseline_wide, row_id)
      setorder(full_wide, row_id)

      edge_names <- setdiff(names(actual_wide), c("row_id", "Group"))
      A <- as.matrix(actual_wide[, ..edge_names])
      B <- as.matrix(baseline_wide[, ..edge_names])
      F <- as.matrix(full_wide[, ..edge_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 + counter * 1009L)
      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(
        analysis_mode = md,
        edge = target_edge,
        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_edges = ncol(Y_residual),
        n_permutations = n_permutations,
        n_bootstrap = n_bootstrap
      )
      counter <- counter + 1L
    }
  }

  rbindlist(output)
}

5 Load the 22 fMRI targets

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)
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
)
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 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     0

fwrite(fmri_audit, file.path(output_dir, "fmri_target_audit.csv"))

6 Create subject-level temporal eye predictors

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 = ", "))
}

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.0001156346

eye_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)
  )
}

eye_coeff_long <- rbindlist(coefficient_rows, use.names = TRUE, fill = FALSE)
eye_basis_audit <- rbindlist(audit_rows, use.names = TRUE, fill = FALSE)

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[, predictor := paste0(
  condition, "__b", sprintf("%02d", basis_index)
)]

eye_coeff_wide <- dcast(
  eye_coeff_long,
  subject + group ~ predictor,
  value.var = "coefficient"
)
eye_coeff_wide[, id_key := normalize_subject_id(subject)]
check_unique_key(eye_coeff_wide, "id_key", "Eye-trajectory table")

predictor_cols <- setdiff(names(eye_coeff_wide), c("subject", "group", "id_key"))
expected_predictors <- length(primary_eye_conditions) * eye_basis_df
if (length(predictor_cols) != expected_predictors) {
  stop("Expected ", expected_predictors,
       " eye predictors but found ", length(predictor_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_predictors.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.4839141

7 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, predictor_cols
)
complete_core <- complete.cases(matched[, ..core_columns])
analysis_data <- matched[complete_core]
analysis_data[, row_id := .I]

predictor_sd_check <- vapply(
  analysis_data[, ..predictor_cols], sd, numeric(1L), na.rm = TRUE
)
invalid_predictors <- names(predictor_sd_check)[
  !is.finite(predictor_sd_check) | predictor_sd_check <= 0
]
if (length(invalid_predictors) > 0L) {
  stop("Zero-variance or invalid eye predictor(s): ",
       paste(invalid_predictors, collapse = ", "))
}

edge_sd_check <- vapply(
  analysis_data[, ..target_edge_cols], sd, numeric(1L), na.rm = TRUE
)
invalid_edges <- names(edge_sd_check)[!is.finite(edge_sd_check) | edge_sd_check <= 0]
if (length(invalid_edges) > 0L) {
  stop("Zero-variance or invalid targeted edge(s): ",
       paste(invalid_edges, 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",
    "eye_predictors",
    "fmri_outcome_edges",
    "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"),
    length(predictor_cols),
    length(target_edge_cols),
    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:                             eye_predictors     24
##  7:                         fmri_outcome_edges     22
##  8:                  motion_available_fraction      0
##  9:                    motion_used_as_nuisance      0
## 10: target_manifest_independently_prespecified      0

fwrite(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 Repeated nested cross-validated prediction: eye tracking -> 22 FC edges

The model treats the 22 targeted FC edges as correlated responses in a multi-response Gaussian elastic net. Nuisance coefficients, eye-predictor residualization, fMRI-outcome residualization, and elastic-net tuning are all estimated strictly within training folds.

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()
for (i in seq_along(analysis_modes_to_run)) {
  mode <- analysis_modes_to_run[i]
  mvpa_runs[[i]] <- run_nested_eye_to_fmri(
    analysis_data = analysis_data,
    predictor_cols = predictor_cols,
    edge_cols = target_edge_cols,
    mode = mode,
    outer_fold_table = outer_fold_table,
    use_motion = use_motion,
    include_icv = include_icv_sensitivity,
    seed = random_seed + i * 1000003L
  )
}

## [2026-08-17 09:03:01] eye_temporal_basis -> targeted_edges22 | total_without_group | outer repeat 1, fold 1
## [2026-08-17 09:03:03] eye_temporal_basis -> targeted_edges22 | total_without_group | outer repeat 1, fold 2
## [2026-08-17 09:03:04] eye_temporal_basis -> targeted_edges22 | total_without_group | outer repeat 1, fold 3
## [2026-08-17 09:03:05] eye_temporal_basis -> targeted_edges22 | total_without_group | outer repeat 1, fold 4
## [2026-08-17 09:03:06] eye_temporal_basis -> targeted_edges22 | total_without_group | outer repeat 1, fold 5
## [2026-08-17 09:03:07] eye_temporal_basis -> targeted_edges22 | total_without_group | outer repeat 2, fold 1
## [2026-08-17 09:03:08] eye_temporal_basis -> targeted_edges22 | total_without_group | outer repeat 2, fold 2
## [2026-08-17 09:03:09] eye_temporal_basis -> targeted_edges22 | total_without_group | outer repeat 2, fold 3
## [2026-08-17 09:03:11] eye_temporal_basis -> targeted_edges22 | total_without_group | outer repeat 2, fold 4
## [2026-08-17 09:03:12] eye_temporal_basis -> targeted_edges22 | total_without_group | outer repeat 2, fold 5
## [2026-08-17 09:03:13] eye_temporal_basis -> targeted_edges22 | total_without_group | outer repeat 3, fold 1
## [2026-08-17 09:03:14] eye_temporal_basis -> targeted_edges22 | total_without_group | outer repeat 3, fold 2
## [2026-08-17 09:03:15] eye_temporal_basis -> targeted_edges22 | total_without_group | outer repeat 3, fold 3
## [2026-08-17 09:03:16] eye_temporal_basis -> targeted_edges22 | total_without_group | outer repeat 3, fold 4
## [2026-08-17 09:03:18] eye_temporal_basis -> targeted_edges22 | total_without_group | outer repeat 3, fold 5
## [2026-08-17 09:03:19] eye_temporal_basis -> targeted_edges22 | total_without_group | outer repeat 4, fold 1
## [2026-08-17 09:03:20] eye_temporal_basis -> targeted_edges22 | total_without_group | outer repeat 4, fold 2
## [2026-08-17 09:03:21] eye_temporal_basis -> targeted_edges22 | total_without_group | outer repeat 4, fold 3
## [2026-08-17 09:03:22] eye_temporal_basis -> targeted_edges22 | total_without_group | outer repeat 4, fold 4
## [2026-08-17 09:03:23] eye_temporal_basis -> targeted_edges22 | total_without_group | outer repeat 4, fold 5
## [2026-08-17 09:03:24] eye_temporal_basis -> targeted_edges22 | total_without_group | outer repeat 5, fold 1
## [2026-08-17 09:03:26] eye_temporal_basis -> targeted_edges22 | total_without_group | outer repeat 5, fold 2
## [2026-08-17 09:03:27] eye_temporal_basis -> targeted_edges22 | total_without_group | outer repeat 5, fold 3
## [2026-08-17 09:03:28] eye_temporal_basis -> targeted_edges22 | total_without_group | outer repeat 5, fold 4
## [2026-08-17 09:03:29] eye_temporal_basis -> targeted_edges22 | total_without_group | outer repeat 5, fold 5
## [2026-08-17 09:03:30] eye_temporal_basis -> targeted_edges22 | total_without_group | outer repeat 6, fold 1
## [2026-08-17 09:03:31] eye_temporal_basis -> targeted_edges22 | total_without_group | outer repeat 6, fold 2
## [2026-08-17 09:03:32] eye_temporal_basis -> targeted_edges22 | total_without_group | outer repeat 6, fold 3
## [2026-08-17 09:03:33] eye_temporal_basis -> targeted_edges22 | total_without_group | outer repeat 6, fold 4
## [2026-08-17 09:03:35] eye_temporal_basis -> targeted_edges22 | total_without_group | outer repeat 6, fold 5
## [2026-08-17 09:03:36] eye_temporal_basis -> targeted_edges22 | total_without_group | outer repeat 7, fold 1
## [2026-08-17 09:03:37] eye_temporal_basis -> targeted_edges22 | total_without_group | outer repeat 7, fold 2
## [2026-08-17 09:03:38] eye_temporal_basis -> targeted_edges22 | total_without_group | outer repeat 7, fold 3
## [2026-08-17 09:03:39] eye_temporal_basis -> targeted_edges22 | total_without_group | outer repeat 7, fold 4
## [2026-08-17 09:03:40] eye_temporal_basis -> targeted_edges22 | total_without_group | outer repeat 7, fold 5
## [2026-08-17 09:03:41] eye_temporal_basis -> targeted_edges22 | total_without_group | outer repeat 8, fold 1
## [2026-08-17 09:03:42] eye_temporal_basis -> targeted_edges22 | total_without_group | outer repeat 8, fold 2
## [2026-08-17 09:03:43] eye_temporal_basis -> targeted_edges22 | total_without_group | outer repeat 8, fold 3
## [2026-08-17 09:03:44] eye_temporal_basis -> targeted_edges22 | total_without_group | outer repeat 8, fold 4
## [2026-08-17 09:03:46] eye_temporal_basis -> targeted_edges22 | total_without_group | outer repeat 8, fold 5
## [2026-08-17 09:03:47] eye_temporal_basis -> targeted_edges22 | total_without_group | outer repeat 9, fold 1
## [2026-08-17 09:03:48] eye_temporal_basis -> targeted_edges22 | total_without_group | outer repeat 9, fold 2
## [2026-08-17 09:03:49] eye_temporal_basis -> targeted_edges22 | total_without_group | outer repeat 9, fold 3
## [2026-08-17 09:03:50] eye_temporal_basis -> targeted_edges22 | total_without_group | outer repeat 9, fold 4
## [2026-08-17 09:03:51] eye_temporal_basis -> targeted_edges22 | total_without_group | outer repeat 9, fold 5
## [2026-08-17 09:03:52] eye_temporal_basis -> targeted_edges22 | total_without_group | outer repeat 10, fold 1
## [2026-08-17 09:03:53] eye_temporal_basis -> targeted_edges22 | total_without_group | outer repeat 10, fold 2
## [2026-08-17 09:03:54] eye_temporal_basis -> targeted_edges22 | total_without_group | outer repeat 10, fold 3
## [2026-08-17 09:03:56] eye_temporal_basis -> targeted_edges22 | total_without_group | outer repeat 10, fold 4
## [2026-08-17 09:03:57] eye_temporal_basis -> targeted_edges22 | total_without_group | outer repeat 10, fold 5
## [2026-08-17 09:03:58] eye_temporal_basis -> targeted_edges22 | incremental_beyond_group | outer repeat 1, fold 1
## [2026-08-17 09:03:59] eye_temporal_basis -> targeted_edges22 | incremental_beyond_group | outer repeat 1, fold 2
## [2026-08-17 09:04:00] eye_temporal_basis -> targeted_edges22 | incremental_beyond_group | outer repeat 1, fold 3
## [2026-08-17 09:04:01] eye_temporal_basis -> targeted_edges22 | incremental_beyond_group | outer repeat 1, fold 4
## [2026-08-17 09:04:02] eye_temporal_basis -> targeted_edges22 | incremental_beyond_group | outer repeat 1, fold 5
## [2026-08-17 09:04:03] eye_temporal_basis -> targeted_edges22 | incremental_beyond_group | outer repeat 2, fold 1
## [2026-08-17 09:04:04] eye_temporal_basis -> targeted_edges22 | incremental_beyond_group | outer repeat 2, fold 2
## [2026-08-17 09:04:06] eye_temporal_basis -> targeted_edges22 | incremental_beyond_group | outer repeat 2, fold 3
## [2026-08-17 09:04:07] eye_temporal_basis -> targeted_edges22 | incremental_beyond_group | outer repeat 2, fold 4
## [2026-08-17 09:04:08] eye_temporal_basis -> targeted_edges22 | incremental_beyond_group | outer repeat 2, fold 5
## [2026-08-17 09:04:09] eye_temporal_basis -> targeted_edges22 | incremental_beyond_group | outer repeat 3, fold 1
## [2026-08-17 09:04:10] eye_temporal_basis -> targeted_edges22 | incremental_beyond_group | outer repeat 3, fold 2
## [2026-08-17 09:04:11] eye_temporal_basis -> targeted_edges22 | incremental_beyond_group | outer repeat 3, fold 3
## [2026-08-17 09:04:12] eye_temporal_basis -> targeted_edges22 | incremental_beyond_group | outer repeat 3, fold 4
## [2026-08-17 09:04:13] eye_temporal_basis -> targeted_edges22 | incremental_beyond_group | outer repeat 3, fold 5
## [2026-08-17 09:04:14] eye_temporal_basis -> targeted_edges22 | incremental_beyond_group | outer repeat 4, fold 1
## [2026-08-17 09:04:16] eye_temporal_basis -> targeted_edges22 | incremental_beyond_group | outer repeat 4, fold 2
## [2026-08-17 09:04:17] eye_temporal_basis -> targeted_edges22 | incremental_beyond_group | outer repeat 4, fold 3
## [2026-08-17 09:04:18] eye_temporal_basis -> targeted_edges22 | incremental_beyond_group | outer repeat 4, fold 4
## [2026-08-17 09:04:19] eye_temporal_basis -> targeted_edges22 | incremental_beyond_group | outer repeat 4, fold 5
## [2026-08-17 09:04:20] eye_temporal_basis -> targeted_edges22 | incremental_beyond_group | outer repeat 5, fold 1
## [2026-08-17 09:04:21] eye_temporal_basis -> targeted_edges22 | incremental_beyond_group | outer repeat 5, fold 2
## [2026-08-17 09:04:22] eye_temporal_basis -> targeted_edges22 | incremental_beyond_group | outer repeat 5, fold 3
## [2026-08-17 09:04:23] eye_temporal_basis -> targeted_edges22 | incremental_beyond_group | outer repeat 5, fold 4
## [2026-08-17 09:04:24] eye_temporal_basis -> targeted_edges22 | incremental_beyond_group | outer repeat 5, fold 5
## [2026-08-17 09:04:26] eye_temporal_basis -> targeted_edges22 | incremental_beyond_group | outer repeat 6, fold 1
## [2026-08-17 09:04:27] eye_temporal_basis -> targeted_edges22 | incremental_beyond_group | outer repeat 6, fold 2
## [2026-08-17 09:04:28] eye_temporal_basis -> targeted_edges22 | incremental_beyond_group | outer repeat 6, fold 3
## [2026-08-17 09:04:29] eye_temporal_basis -> targeted_edges22 | incremental_beyond_group | outer repeat 6, fold 4
## [2026-08-17 09:04:30] eye_temporal_basis -> targeted_edges22 | incremental_beyond_group | outer repeat 6, fold 5
## [2026-08-17 09:04:31] eye_temporal_basis -> targeted_edges22 | incremental_beyond_group | outer repeat 7, fold 1
## [2026-08-17 09:04:32] eye_temporal_basis -> targeted_edges22 | incremental_beyond_group | outer repeat 7, fold 2
## [2026-08-17 09:04:33] eye_temporal_basis -> targeted_edges22 | incremental_beyond_group | outer repeat 7, fold 3
## [2026-08-17 09:04:34] eye_temporal_basis -> targeted_edges22 | incremental_beyond_group | outer repeat 7, fold 4
## [2026-08-17 09:04:36] eye_temporal_basis -> targeted_edges22 | incremental_beyond_group | outer repeat 7, fold 5
## [2026-08-17 09:04:37] eye_temporal_basis -> targeted_edges22 | incremental_beyond_group | outer repeat 8, fold 1
## [2026-08-17 09:04:38] eye_temporal_basis -> targeted_edges22 | incremental_beyond_group | outer repeat 8, fold 2
## [2026-08-17 09:04:39] eye_temporal_basis -> targeted_edges22 | incremental_beyond_group | outer repeat 8, fold 3
## [2026-08-17 09:04:40] eye_temporal_basis -> targeted_edges22 | incremental_beyond_group | outer repeat 8, fold 4
## [2026-08-17 09:04:41] eye_temporal_basis -> targeted_edges22 | incremental_beyond_group | outer repeat 8, fold 5
## [2026-08-17 09:04:42] eye_temporal_basis -> targeted_edges22 | incremental_beyond_group | outer repeat 9, fold 1
## [2026-08-17 09:04:43] eye_temporal_basis -> targeted_edges22 | incremental_beyond_group | outer repeat 9, fold 2
## [2026-08-17 09:04:44] eye_temporal_basis -> targeted_edges22 | incremental_beyond_group | outer repeat 9, fold 3
## [2026-08-17 09:04:45] eye_temporal_basis -> targeted_edges22 | incremental_beyond_group | outer repeat 9, fold 4
## [2026-08-17 09:04:47] eye_temporal_basis -> targeted_edges22 | incremental_beyond_group | outer repeat 9, fold 5
## [2026-08-17 09:04:48] eye_temporal_basis -> targeted_edges22 | incremental_beyond_group | outer repeat 10, fold 1
## [2026-08-17 09:04:49] eye_temporal_basis -> targeted_edges22 | incremental_beyond_group | outer repeat 10, fold 2
## [2026-08-17 09:04:50] eye_temporal_basis -> targeted_edges22 | incremental_beyond_group | outer repeat 10, fold 3
## [2026-08-17 09:04:51] eye_temporal_basis -> targeted_edges22 | incremental_beyond_group | outer repeat 10, fold 4
## [2026-08-17 09:04:52] eye_temporal_basis -> targeted_edges22 | incremental_beyond_group | outer repeat 10, fold 5

all_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.")
}

fwrite(
  all_predictions,
  file.path(output_dir, "eye_to_fmri_nested_cv_all_outer_predictions_long.csv")
)
fwrite(
  all_hyperparameters,
  file.path(output_dir, "eye_to_fmri_nested_cv_selected_hyperparameters.csv")
)
fwrite(
  all_weights,
  file.path(output_dir, "eye_to_fmri_nested_cv_selected_model_weights_long.csv")
)

9 Out-of-fold FC prediction performance

prediction_average <- average_oof_predictions(all_predictions)
fwrite(
  prediction_average,
  file.path(output_dir, "eye_to_fmri_subject_averaged_predictions_long.csv")
)

edge_metrics <- edge_metric_tables(prediction_average)

edge_labels <- unique(manifest[, .(
  edge = FC,
  edge_label = if (all(c("node1_label", "node2_label") %in% names(manifest))) {
    paste(node1_label, node2_label, sep = " -- ")
  } else {
    FC
  },
  axis,
  subfamily
)])
edge_metrics$by_edge <- merge(
  edge_metrics$by_edge,
  edge_labels,
  by = "edge",
  all.x = TRUE,
  sort = FALSE
)
edge_metrics$by_edge[is.na(edge_label), edge_label := edge]

fwrite(
  edge_metrics$by_edge,
  file.path(output_dir, "eye_to_fmri_metrics_by_edge.csv")
)
fwrite(
  edge_metrics$overall,
  file.path(output_dir, "eye_to_fmri_metrics_overall.csv")
)

print(edge_metrics$overall)

##               analysis_mode   prediction_source     n     rmse       mae
##                      <char>              <fctr> <int>    <num>     <num>
## 1:      total_without_group baseline_prediction  2266 1.011808 0.7930222
## 2: incremental_beyond_group baseline_prediction  2266 1.006155 0.7874902
## 3:      total_without_group     full_prediction  2266 1.011808 0.7930222
## 4: incremental_beyond_group     full_prediction  2266 1.006155 0.7874902
##           q2 correlation
##        <num>       <num>
## 1: 0.3257908   0.5717652
## 2: 0.3333031   0.5789558
## 3: 0.3257908   0.5717652
## 4: 0.3333031   0.5789558

print(edge_metrics$by_edge)

##          edge            analysis_mode   prediction_source     n      rmse
##        <char>                   <char>              <fctr> <int>     <num>
##  1:  FC51_181      total_without_group baseline_prediction   103 0.2227894
##  2: FC125_181      total_without_group baseline_prediction   103 0.1824466
##  3:  FC48_155      total_without_group baseline_prediction   103 0.2595477
##  4: FC122_155      total_without_group baseline_prediction   103 0.2766308
##  5:  FC48_156      total_without_group baseline_prediction   103 0.2550055
##  6: FC122_156      total_without_group baseline_prediction   103 0.2516392
##  7:   FC48_80      total_without_group baseline_prediction   103 0.2531951
##  8:  FC80_122      total_without_group baseline_prediction   103 0.2579233
##  9:   FC80_89      total_without_group baseline_prediction   103 0.1873863
## 10:   FC80_90      total_without_group baseline_prediction   103 0.2514881
## 11:   FC89_90      total_without_group baseline_prediction   103 0.2235822
## 12: FC155_163      total_without_group baseline_prediction   103 0.2042408
## 13: FC155_164      total_without_group baseline_prediction   103 0.2957900
## 14:   FC16_43      total_without_group baseline_prediction   103 0.1391318
## 15:  FC16_117      total_without_group baseline_prediction   103 0.1554037
## 16:   FC16_60      total_without_group baseline_prediction   103 0.1780170
## 17:  FC16_134      total_without_group baseline_prediction   103 0.1623931
## 18:  FC12_160      total_without_group baseline_prediction   103 0.2260080
## 19:  FC12_143      total_without_group baseline_prediction   103 0.2101566
## 20:   FC12_69      total_without_group baseline_prediction   103 0.2024901
## 21:  FC12_118      total_without_group baseline_prediction   103 0.2183862
## 22:   FC12_44      total_without_group baseline_prediction   103 0.2133637
## 23:  FC51_181 incremental_beyond_group baseline_prediction   103 0.2119069
## 24: FC125_181 incremental_beyond_group baseline_prediction   103 0.1794127
## 25:  FC48_155 incremental_beyond_group baseline_prediction   103 0.2541106
## 26: FC122_155 incremental_beyond_group baseline_prediction   103 0.2722234
## 27:  FC48_156 incremental_beyond_group baseline_prediction   103 0.2564579
## 28: FC122_156 incremental_beyond_group baseline_prediction   103 0.2539940
## 29:   FC48_80 incremental_beyond_group baseline_prediction   103 0.2528686
## 30:  FC80_122 incremental_beyond_group baseline_prediction   103 0.2603784
## 31:   FC80_89 incremental_beyond_group baseline_prediction   103 0.1867590
## 32:   FC80_90 incremental_beyond_group baseline_prediction   103 0.2541986
## 33:   FC89_90 incremental_beyond_group baseline_prediction   103 0.2252028
## 34: FC155_163 incremental_beyond_group baseline_prediction   103 0.2042559
## 35: FC155_164 incremental_beyond_group baseline_prediction   103 0.2972534
## 36:   FC16_43 incremental_beyond_group baseline_prediction   103 0.1353064
## 37:  FC16_117 incremental_beyond_group baseline_prediction   103 0.1539635
## 38:   FC16_60 incremental_beyond_group baseline_prediction   103 0.1788390
## 39:  FC16_134 incremental_beyond_group baseline_prediction   103 0.1633798
## 40:  FC12_160 incremental_beyond_group baseline_prediction   103 0.2256895
## 41:  FC12_143 incremental_beyond_group baseline_prediction   103 0.2066155
## 42:   FC12_69 incremental_beyond_group baseline_prediction   103 0.1976441
## 43:  FC12_118 incremental_beyond_group baseline_prediction   103 0.2173570
## 44:   FC12_44 incremental_beyond_group baseline_prediction   103 0.2140016
## 45:  FC51_181      total_without_group     full_prediction   103 0.2227894
## 46: FC125_181      total_without_group     full_prediction   103 0.1824466
## 47:  FC48_155      total_without_group     full_prediction   103 0.2595477
## 48: FC122_155      total_without_group     full_prediction   103 0.2766308
## 49:  FC48_156      total_without_group     full_prediction   103 0.2550055
## 50: FC122_156      total_without_group     full_prediction   103 0.2516392
## 51:   FC48_80      total_without_group     full_prediction   103 0.2531951
## 52:  FC80_122      total_without_group     full_prediction   103 0.2579233
## 53:   FC80_89      total_without_group     full_prediction   103 0.1873863
## 54:   FC80_90      total_without_group     full_prediction   103 0.2514881
## 55:   FC89_90      total_without_group     full_prediction   103 0.2235822
## 56: FC155_163      total_without_group     full_prediction   103 0.2042408
## 57: FC155_164      total_without_group     full_prediction   103 0.2957900
## 58:   FC16_43      total_without_group     full_prediction   103 0.1391318
## 59:  FC16_117      total_without_group     full_prediction   103 0.1554037
## 60:   FC16_60      total_without_group     full_prediction   103 0.1780170
## 61:  FC16_134      total_without_group     full_prediction   103 0.1623931
## 62:  FC12_160      total_without_group     full_prediction   103 0.2260080
## 63:  FC12_143      total_without_group     full_prediction   103 0.2101566
## 64:   FC12_69      total_without_group     full_prediction   103 0.2024901
## 65:  FC12_118      total_without_group     full_prediction   103 0.2183862
## 66:   FC12_44      total_without_group     full_prediction   103 0.2133637
## 67:  FC51_181 incremental_beyond_group     full_prediction   103 0.2119069
## 68: FC125_181 incremental_beyond_group     full_prediction   103 0.1794127
## 69:  FC48_155 incremental_beyond_group     full_prediction   103 0.2541106
## 70: FC122_155 incremental_beyond_group     full_prediction   103 0.2722234
## 71:  FC48_156 incremental_beyond_group     full_prediction   103 0.2564579
## 72: FC122_156 incremental_beyond_group     full_prediction   103 0.2539940
## 73:   FC48_80 incremental_beyond_group     full_prediction   103 0.2528686
## 74:  FC80_122 incremental_beyond_group     full_prediction   103 0.2603784
## 75:   FC80_89 incremental_beyond_group     full_prediction   103 0.1867590
## 76:   FC80_90 incremental_beyond_group     full_prediction   103 0.2541986
## 77:   FC89_90 incremental_beyond_group     full_prediction   103 0.2252028
## 78: FC155_163 incremental_beyond_group     full_prediction   103 0.2042559
## 79: FC155_164 incremental_beyond_group     full_prediction   103 0.2972534
## 80:   FC16_43 incremental_beyond_group     full_prediction   103 0.1353064
## 81:  FC16_117 incremental_beyond_group     full_prediction   103 0.1539635
## 82:   FC16_60 incremental_beyond_group     full_prediction   103 0.1788390
## 83:  FC16_134 incremental_beyond_group     full_prediction   103 0.1633798
## 84:  FC12_160 incremental_beyond_group     full_prediction   103 0.2256895
## 85:  FC12_143 incremental_beyond_group     full_prediction   103 0.2066155
## 86:   FC12_69 incremental_beyond_group     full_prediction   103 0.1976441
## 87:  FC12_118 incremental_beyond_group     full_prediction   103 0.2173570
## 88:   FC12_44 incremental_beyond_group     full_prediction   103 0.2140016
##          edge            analysis_mode   prediction_source     n      rmse
##           mae           q2   correlation
##         <num>        <num>         <num>
##  1: 0.1728940  0.001970377  0.1212565926
##  2: 0.1458961  0.125325153  0.3565086731
##  3: 0.2023050 -0.020394828  0.0575956593
##  4: 0.2069974 -0.043462773 -0.1469635879
##  5: 0.2136048 -0.052485038 -0.1482850684
##  6: 0.1965177 -0.027710823  0.0328839490
##  7: 0.2023712 -0.020138464  0.0846140083
##  8: 0.2013727 -0.060739807 -0.1300439808
##  9: 0.1413366 -0.048369982 -0.1031858956
## 10: 0.1899375 -0.040859866 -0.0268098517
## 11: 0.1733683 -0.065817629 -0.3330384552
## 12: 0.1648401 -0.055417833 -0.3581890693
## 13: 0.2299909 -0.045670344 -0.2955569192
## 14: 0.1117729 -0.037834880 -0.0933464605
## 15: 0.1239696 -0.075976121 -0.1880866379
## 16: 0.1403859 -0.016136500  0.0423161026
## 17: 0.1273142 -0.020680263  0.0019983976
## 18: 0.1795022 -0.053028360 -0.2951754757
## 19: 0.1714989 -0.058285822 -0.4533313972
## 20: 0.1536814 -0.029009622  0.0017531714
## 21: 0.1681796 -0.041968814 -0.0871620991
## 22: 0.1626239 -0.056725381 -0.5191852699
## 23: 0.1611102  0.097089687  0.3231040475
## 24: 0.1434087  0.154173059  0.3966883610
## 25: 0.2037599  0.021909049  0.1956804243
## 26: 0.2076030 -0.010477355  0.0971296026
## 27: 0.2149560 -0.064508422 -0.1060010052
## 28: 0.1998399 -0.047034631  0.0009304731
## 29: 0.2034826 -0.017509425  0.1103989980
## 30: 0.2028417 -0.081029356 -0.1797599531
## 31: 0.1400247 -0.041362989  0.0186683792
## 32: 0.1921608 -0.063416883 -0.0790316092
## 33: 0.1733265 -0.081325082 -0.2678243047
## 34: 0.1634089 -0.055573912 -0.1323179918
## 35: 0.2307481 -0.056042341 -0.2238332957
## 36: 0.1074781  0.018451581  0.1841435078
## 37: 0.1215696 -0.056125307  0.0105038403
## 38: 0.1397571 -0.025541324  0.0296954586
## 39: 0.1296532 -0.033121064 -0.0017700302
## 40: 0.1773874 -0.050062521 -0.0233647751
## 41: 0.1663488 -0.022922494  0.0812724357
## 42: 0.1495426  0.019653954  0.1842185457
## 43: 0.1688101 -0.032170312  0.0383955606
## 44: 0.1622654 -0.063053192 -0.1672984225
## 45: 0.1728940  0.001970377  0.1212565926
## 46: 0.1458961  0.125325153  0.3565086731
## 47: 0.2023050 -0.020394828  0.0575956593
## 48: 0.2069974 -0.043462773 -0.1469635879
## 49: 0.2136048 -0.052485038 -0.1482850684
## 50: 0.1965177 -0.027710823  0.0328839490
## 51: 0.2023712 -0.020138464  0.0846140083
## 52: 0.2013727 -0.060739807 -0.1300439808
## 53: 0.1413366 -0.048369982 -0.1031858956
## 54: 0.1899375 -0.040859866 -0.0268098517
## 55: 0.1733683 -0.065817629 -0.3330384552
## 56: 0.1648401 -0.055417833 -0.3581890693
## 57: 0.2299909 -0.045670344 -0.2955569192
## 58: 0.1117729 -0.037834880 -0.0933464605
## 59: 0.1239696 -0.075976121 -0.1880866379
## 60: 0.1403859 -0.016136500  0.0423161026
## 61: 0.1273142 -0.020680263  0.0019983976
## 62: 0.1795022 -0.053028360 -0.2951754757
## 63: 0.1714989 -0.058285822 -0.4533313972
## 64: 0.1536814 -0.029009622  0.0017531714
## 65: 0.1681796 -0.041968814 -0.0871620991
## 66: 0.1626239 -0.056725381 -0.5191852699
## 67: 0.1611102  0.097089687  0.3231040475
## 68: 0.1434087  0.154173059  0.3966883610
## 69: 0.2037599  0.021909049  0.1956804243
## 70: 0.2076030 -0.010477355  0.0971296026
## 71: 0.2149560 -0.064508422 -0.1060010052
## 72: 0.1998399 -0.047034631  0.0009304731
## 73: 0.2034826 -0.017509425  0.1103989980
## 74: 0.2028417 -0.081029356 -0.1797599531
## 75: 0.1400247 -0.041362989  0.0186683792
## 76: 0.1921608 -0.063416883 -0.0790316092
## 77: 0.1733265 -0.081325082 -0.2678243047
## 78: 0.1634089 -0.055573912 -0.1323179918
## 79: 0.2307481 -0.056042341 -0.2238332957
## 80: 0.1074781  0.018451581  0.1841435078
## 81: 0.1215696 -0.056125307  0.0105038403
## 82: 0.1397571 -0.025541324  0.0296954586
## 83: 0.1296532 -0.033121064 -0.0017700302
## 84: 0.1773874 -0.050062521 -0.0233647751
## 85: 0.1663488 -0.022922494  0.0812724357
## 86: 0.1495426  0.019653954  0.1842185457
## 87: 0.1688101 -0.032170312  0.0383955606
## 88: 0.1622654 -0.063053192 -0.1672984225
##           mae           q2   correlation
##                                                       edge_label
##                                                           <char>
##  1:        ctx_lh_G_cingul-Post-ventral -- ctx_rh_S_pericallosal
##  2:        ctx_rh_G_cingul-Post-ventral -- ctx_rh_S_pericallosal
##  3: ctx_lh_G_and_S_cingul-Mid-Ant -- ctx_rh_Lat_Fis-ant-Vertical
##  4: ctx_rh_G_and_S_cingul-Mid-Ant -- ctx_rh_Lat_Fis-ant-Vertical
##  5:         ctx_lh_G_and_S_cingul-Mid-Ant -- ctx_rh_Lat_Fis-post
##  6:         ctx_rh_G_and_S_cingul-Mid-Ant -- ctx_rh_Lat_Fis-post
##  7: ctx_lh_G_and_S_cingul-Mid-Ant -- ctx_lh_Lat_Fis-ant-Horizont
##  8: ctx_rh_G_and_S_cingul-Mid-Ant -- ctx_lh_Lat_Fis-ant-Horizont
##  9:  ctx_lh_Lat_Fis-ant-Horizont -- ctx_lh_S_circular_insula_inf
## 10:  ctx_lh_Lat_Fis-ant-Horizont -- ctx_lh_S_circular_insula_sup
## 11: ctx_lh_S_circular_insula_inf -- ctx_lh_S_circular_insula_sup
## 12:  ctx_rh_Lat_Fis-ant-Vertical -- ctx_rh_S_circular_insula_inf
## 13:  ctx_rh_Lat_Fis-ant-Vertical -- ctx_rh_S_circular_insula_sup
## 14:          Left-Accumbens-area -- ctx_lh_G_and_S_occipital_inf
## 15:          Left-Accumbens-area -- ctx_rh_G_and_S_occipital_inf
## 16:             Left-Accumbens-area -- ctx_lh_G_occipital_middle
## 17:             Left-Accumbens-area -- ctx_rh_G_occipital_middle
## 18:                               Brain-Stem -- ctx_rh_S_central
## 19:                           Brain-Stem -- ctx_rh_G_postcentral
## 20:                           Brain-Stem -- ctx_lh_G_postcentral
## 21:                     Brain-Stem -- ctx_rh_G_and_S_paracentral
## 22:                     Brain-Stem -- ctx_lh_G_and_S_paracentral
## 23:        ctx_lh_G_cingul-Post-ventral -- ctx_rh_S_pericallosal
## 24:        ctx_rh_G_cingul-Post-ventral -- ctx_rh_S_pericallosal
## 25: ctx_lh_G_and_S_cingul-Mid-Ant -- ctx_rh_Lat_Fis-ant-Vertical
## 26: ctx_rh_G_and_S_cingul-Mid-Ant -- ctx_rh_Lat_Fis-ant-Vertical
## 27:         ctx_lh_G_and_S_cingul-Mid-Ant -- ctx_rh_Lat_Fis-post
## 28:         ctx_rh_G_and_S_cingul-Mid-Ant -- ctx_rh_Lat_Fis-post
## 29: ctx_lh_G_and_S_cingul-Mid-Ant -- ctx_lh_Lat_Fis-ant-Horizont
## 30: ctx_rh_G_and_S_cingul-Mid-Ant -- ctx_lh_Lat_Fis-ant-Horizont
## 31:  ctx_lh_Lat_Fis-ant-Horizont -- ctx_lh_S_circular_insula_inf
## 32:  ctx_lh_Lat_Fis-ant-Horizont -- ctx_lh_S_circular_insula_sup
## 33: ctx_lh_S_circular_insula_inf -- ctx_lh_S_circular_insula_sup
## 34:  ctx_rh_Lat_Fis-ant-Vertical -- ctx_rh_S_circular_insula_inf
## 35:  ctx_rh_Lat_Fis-ant-Vertical -- ctx_rh_S_circular_insula_sup
## 36:          Left-Accumbens-area -- ctx_lh_G_and_S_occipital_inf
## 37:          Left-Accumbens-area -- ctx_rh_G_and_S_occipital_inf
## 38:             Left-Accumbens-area -- ctx_lh_G_occipital_middle
## 39:             Left-Accumbens-area -- ctx_rh_G_occipital_middle
## 40:                               Brain-Stem -- ctx_rh_S_central
## 41:                           Brain-Stem -- ctx_rh_G_postcentral
## 42:                           Brain-Stem -- ctx_lh_G_postcentral
## 43:                     Brain-Stem -- ctx_rh_G_and_S_paracentral
## 44:                     Brain-Stem -- ctx_lh_G_and_S_paracentral
## 45:        ctx_lh_G_cingul-Post-ventral -- ctx_rh_S_pericallosal
## 46:        ctx_rh_G_cingul-Post-ventral -- ctx_rh_S_pericallosal
## 47: ctx_lh_G_and_S_cingul-Mid-Ant -- ctx_rh_Lat_Fis-ant-Vertical
## 48: ctx_rh_G_and_S_cingul-Mid-Ant -- ctx_rh_Lat_Fis-ant-Vertical
## 49:         ctx_lh_G_and_S_cingul-Mid-Ant -- ctx_rh_Lat_Fis-post
## 50:         ctx_rh_G_and_S_cingul-Mid-Ant -- ctx_rh_Lat_Fis-post
## 51: ctx_lh_G_and_S_cingul-Mid-Ant -- ctx_lh_Lat_Fis-ant-Horizont
## 52: ctx_rh_G_and_S_cingul-Mid-Ant -- ctx_lh_Lat_Fis-ant-Horizont
## 53:  ctx_lh_Lat_Fis-ant-Horizont -- ctx_lh_S_circular_insula_inf
## 54:  ctx_lh_Lat_Fis-ant-Horizont -- ctx_lh_S_circular_insula_sup
## 55: ctx_lh_S_circular_insula_inf -- ctx_lh_S_circular_insula_sup
## 56:  ctx_rh_Lat_Fis-ant-Vertical -- ctx_rh_S_circular_insula_inf
## 57:  ctx_rh_Lat_Fis-ant-Vertical -- ctx_rh_S_circular_insula_sup
## 58:          Left-Accumbens-area -- ctx_lh_G_and_S_occipital_inf
## 59:          Left-Accumbens-area -- ctx_rh_G_and_S_occipital_inf
## 60:             Left-Accumbens-area -- ctx_lh_G_occipital_middle
## 61:             Left-Accumbens-area -- ctx_rh_G_occipital_middle
## 62:                               Brain-Stem -- ctx_rh_S_central
## 63:                           Brain-Stem -- ctx_rh_G_postcentral
## 64:                           Brain-Stem -- ctx_lh_G_postcentral
## 65:                     Brain-Stem -- ctx_rh_G_and_S_paracentral
## 66:                     Brain-Stem -- ctx_lh_G_and_S_paracentral
## 67:        ctx_lh_G_cingul-Post-ventral -- ctx_rh_S_pericallosal
## 68:        ctx_rh_G_cingul-Post-ventral -- ctx_rh_S_pericallosal
## 69: ctx_lh_G_and_S_cingul-Mid-Ant -- ctx_rh_Lat_Fis-ant-Vertical
## 70: ctx_rh_G_and_S_cingul-Mid-Ant -- ctx_rh_Lat_Fis-ant-Vertical
## 71:         ctx_lh_G_and_S_cingul-Mid-Ant -- ctx_rh_Lat_Fis-post
## 72:         ctx_rh_G_and_S_cingul-Mid-Ant -- ctx_rh_Lat_Fis-post
## 73: ctx_lh_G_and_S_cingul-Mid-Ant -- ctx_lh_Lat_Fis-ant-Horizont
## 74: ctx_rh_G_and_S_cingul-Mid-Ant -- ctx_lh_Lat_Fis-ant-Horizont
## 75:  ctx_lh_Lat_Fis-ant-Horizont -- ctx_lh_S_circular_insula_inf
## 76:  ctx_lh_Lat_Fis-ant-Horizont -- ctx_lh_S_circular_insula_sup
## 77: ctx_lh_S_circular_insula_inf -- ctx_lh_S_circular_insula_sup
## 78:  ctx_rh_Lat_Fis-ant-Vertical -- ctx_rh_S_circular_insula_inf
## 79:  ctx_rh_Lat_Fis-ant-Vertical -- ctx_rh_S_circular_insula_sup
## 80:          Left-Accumbens-area -- ctx_lh_G_and_S_occipital_inf
## 81:          Left-Accumbens-area -- ctx_rh_G_and_S_occipital_inf
## 82:             Left-Accumbens-area -- ctx_lh_G_occipital_middle
## 83:             Left-Accumbens-area -- ctx_rh_G_occipital_middle
## 84:                               Brain-Stem -- ctx_rh_S_central
## 85:                           Brain-Stem -- ctx_rh_G_postcentral
## 86:                           Brain-Stem -- ctx_lh_G_postcentral
## 87:                     Brain-Stem -- ctx_rh_G_and_S_paracentral
## 88:                     Brain-Stem -- ctx_lh_G_and_S_paracentral
##                                                       edge_label
##                                axis              subfamily
##                              <char>                 <char>
##  1:      Axis1_Cortical_Integration       posterior_medial
##  2:      Axis1_Cortical_Integration       posterior_medial
##  3:      Axis1_Cortical_Integration    cingulo_perisylvian
##  4:      Axis1_Cortical_Integration    cingulo_perisylvian
##  5:      Axis1_Cortical_Integration    cingulo_perisylvian
##  6:      Axis1_Cortical_Integration    cingulo_perisylvian
##  7:      Axis1_Cortical_Integration    cingulo_perisylvian
##  8:      Axis1_Cortical_Integration    cingulo_perisylvian
##  9:      Axis1_Cortical_Integration      insular_opercular
## 10:      Axis1_Cortical_Integration      insular_opercular
## 11:      Axis1_Cortical_Integration      insular_opercular
## 12:      Axis1_Cortical_Integration      insular_opercular
## 13:      Axis1_Cortical_Integration      insular_opercular
## 14: Axis2_Perceptual_Motor_Coupling          reward_visual
## 15: Axis2_Perceptual_Motor_Coupling          reward_visual
## 16: Axis2_Perceptual_Motor_Coupling          reward_visual
## 17: Axis2_Perceptual_Motor_Coupling          reward_visual
## 18: Axis2_Perceptual_Motor_Coupling brainstem_sensorimotor
## 19: Axis2_Perceptual_Motor_Coupling brainstem_sensorimotor
## 20: Axis2_Perceptual_Motor_Coupling brainstem_sensorimotor
## 21: Axis2_Perceptual_Motor_Coupling brainstem_sensorimotor
## 22: Axis2_Perceptual_Motor_Coupling brainstem_sensorimotor
## 23:      Axis1_Cortical_Integration       posterior_medial
## 24:      Axis1_Cortical_Integration       posterior_medial
## 25:      Axis1_Cortical_Integration    cingulo_perisylvian
## 26:      Axis1_Cortical_Integration    cingulo_perisylvian
## 27:      Axis1_Cortical_Integration    cingulo_perisylvian
## 28:      Axis1_Cortical_Integration    cingulo_perisylvian
## 29:      Axis1_Cortical_Integration    cingulo_perisylvian
## 30:      Axis1_Cortical_Integration    cingulo_perisylvian
## 31:      Axis1_Cortical_Integration      insular_opercular
## 32:      Axis1_Cortical_Integration      insular_opercular
## 33:      Axis1_Cortical_Integration      insular_opercular
## 34:      Axis1_Cortical_Integration      insular_opercular
## 35:      Axis1_Cortical_Integration      insular_opercular
## 36: Axis2_Perceptual_Motor_Coupling          reward_visual
## 37: Axis2_Perceptual_Motor_Coupling          reward_visual
## 38: Axis2_Perceptual_Motor_Coupling          reward_visual
## 39: Axis2_Perceptual_Motor_Coupling          reward_visual
## 40: Axis2_Perceptual_Motor_Coupling brainstem_sensorimotor
## 41: Axis2_Perceptual_Motor_Coupling brainstem_sensorimotor
## 42: Axis2_Perceptual_Motor_Coupling brainstem_sensorimotor
## 43: Axis2_Perceptual_Motor_Coupling brainstem_sensorimotor
## 44: Axis2_Perceptual_Motor_Coupling brainstem_sensorimotor
## 45:      Axis1_Cortical_Integration       posterior_medial
## 46:      Axis1_Cortical_Integration       posterior_medial
## 47:      Axis1_Cortical_Integration    cingulo_perisylvian
## 48:      Axis1_Cortical_Integration    cingulo_perisylvian
## 49:      Axis1_Cortical_Integration    cingulo_perisylvian
## 50:      Axis1_Cortical_Integration    cingulo_perisylvian
## 51:      Axis1_Cortical_Integration    cingulo_perisylvian
## 52:      Axis1_Cortical_Integration    cingulo_perisylvian
## 53:      Axis1_Cortical_Integration      insular_opercular
## 54:      Axis1_Cortical_Integration      insular_opercular
## 55:      Axis1_Cortical_Integration      insular_opercular
## 56:      Axis1_Cortical_Integration      insular_opercular
## 57:      Axis1_Cortical_Integration      insular_opercular
## 58: Axis2_Perceptual_Motor_Coupling          reward_visual
## 59: Axis2_Perceptual_Motor_Coupling          reward_visual
## 60: Axis2_Perceptual_Motor_Coupling          reward_visual
## 61: Axis2_Perceptual_Motor_Coupling          reward_visual
## 62: Axis2_Perceptual_Motor_Coupling brainstem_sensorimotor
## 63: Axis2_Perceptual_Motor_Coupling brainstem_sensorimotor
## 64: Axis2_Perceptual_Motor_Coupling brainstem_sensorimotor
## 65: Axis2_Perceptual_Motor_Coupling brainstem_sensorimotor
## 66: Axis2_Perceptual_Motor_Coupling brainstem_sensorimotor
## 67:      Axis1_Cortical_Integration       posterior_medial
## 68:      Axis1_Cortical_Integration       posterior_medial
## 69:      Axis1_Cortical_Integration    cingulo_perisylvian
## 70:      Axis1_Cortical_Integration    cingulo_perisylvian
## 71:      Axis1_Cortical_Integration    cingulo_perisylvian
## 72:      Axis1_Cortical_Integration    cingulo_perisylvian
## 73:      Axis1_Cortical_Integration    cingulo_perisylvian
## 74:      Axis1_Cortical_Integration    cingulo_perisylvian
## 75:      Axis1_Cortical_Integration      insular_opercular
## 76:      Axis1_Cortical_Integration      insular_opercular
## 77:      Axis1_Cortical_Integration      insular_opercular
## 78:      Axis1_Cortical_Integration      insular_opercular
## 79:      Axis1_Cortical_Integration      insular_opercular
## 80: Axis2_Perceptual_Motor_Coupling          reward_visual
## 81: Axis2_Perceptual_Motor_Coupling          reward_visual
## 82: Axis2_Perceptual_Motor_Coupling          reward_visual
## 83: Axis2_Perceptual_Motor_Coupling          reward_visual
## 84: Axis2_Perceptual_Motor_Coupling brainstem_sensorimotor
## 85: Axis2_Perceptual_Motor_Coupling brainstem_sensorimotor
## 86: Axis2_Perceptual_Motor_Coupling brainstem_sensorimotor
## 87: Axis2_Perceptual_Motor_Coupling brainstem_sensorimotor
## 88: Axis2_Perceptual_Motor_Coupling brainstem_sensorimotor
##                                axis              subfamily

edge_metric_wide <- dcast(
  edge_metrics$by_edge,
  analysis_mode + edge + edge_label + axis + subfamily ~ prediction_source,
  value.var = c("rmse", "mae", "q2", "correlation")
)
edge_metric_wide[, `:=`(
  delta_q2 = q2_full_prediction - q2_baseline_prediction,
  delta_rmse = rmse_baseline_prediction - rmse_full_prediction
)]
fwrite(
  edge_metric_wide,
  file.path(output_dir, "eye_to_fmri_edge_prediction_improvement.csv")
)

plot_data <- edge_metric_wide[analysis_mode == "incremental_beyond_group"]
setorder(plot_data, delta_q2)
plot_data[, edge_label := factor(edge_label, levels = unique(edge_label))]

p_edge_q2 <- ggplot(plot_data, aes(x = delta_q2, y = edge_label)) +
  geom_vline(xintercept = 0, linetype = 3) +
  geom_point() +
  theme_classic() +
  labs(
    title = "Eye-tracking increment in held-out FC prediction beyond diagnosis",
    x = "Delta Q2: full eye model minus nuisance baseline",
    y = NULL
  )
print(p_edge_q2)

ggsave(
  file.path(output_dir, "eye_to_fmri_delta_q2_by_edge.png"),
  p_edge_q2,
  width = 11,
  height = 8,
  dpi = 300
)

10 Prediction-alignment permutation test and uncertainty

The statistic is the reduction in standardized held-out mean squared error when the eye-tracking increment is added to the nuisance-only prediction. Positive values favor the eye-to-fMRI model. For incremental_beyond_group, permutation is restricted within diagnosis. The held-out predictions are kept fixed, so this is not a full pipeline re-fitting permutation.

prediction_tests <- prediction_alignment_test_edges(
  prediction_average = prediction_average,
  n_permutations = prediction_permutations,
  n_bootstrap = prediction_bootstraps,
  seed = random_seed
)

# Multiplicity correction is applied across the 22 individual edges within each
# analysis mode. The ALL_EDGES omnibus row is reported separately.
prediction_tests[, p_bh_22_edges := NA_real_]
prediction_tests[edge != "ALL_EDGES", p_bh_22_edges := p.adjust(
  p_fixed_prediction_permutation,
  method = "BH"
), by = analysis_mode]

prediction_tests <- merge(
  prediction_tests,
  edge_labels,
  by = "edge",
  all.x = TRUE,
  sort = FALSE
)
prediction_tests[edge == "ALL_EDGES", edge_label := "ALL_EDGES"]

print(prediction_tests)

## Index: <edge>
##          edge            analysis_mode     delta_mse proportional_mse_reduction
##        <char>                   <char>         <num>                      <num>
##  1: ALL_EDGES      total_without_group  0.000000e+00               0.000000e+00
##  2:  FC51_181      total_without_group  0.000000e+00               0.000000e+00
##  3: FC125_181      total_without_group  0.000000e+00               0.000000e+00
##  4:  FC48_155      total_without_group  0.000000e+00               0.000000e+00
##  5: FC122_155      total_without_group  0.000000e+00               0.000000e+00
##  6:  FC48_156      total_without_group  0.000000e+00               0.000000e+00
##  7: FC122_156      total_without_group  2.220446e-16               2.181757e-16
##  8:   FC48_80      total_without_group  0.000000e+00               0.000000e+00
##  9:  FC80_122      total_without_group  0.000000e+00               0.000000e+00
## 10:   FC80_89      total_without_group  0.000000e+00               0.000000e+00
## 11:   FC80_90      total_without_group  0.000000e+00               0.000000e+00
## 12:   FC89_90      total_without_group  0.000000e+00               0.000000e+00
## 13: FC155_163      total_without_group  0.000000e+00               0.000000e+00
## 14: FC155_164      total_without_group  0.000000e+00               0.000000e+00
## 15:   FC16_43      total_without_group  0.000000e+00               0.000000e+00
## 16:  FC16_117      total_without_group  0.000000e+00               0.000000e+00
## 17:   FC16_60      total_without_group  0.000000e+00               0.000000e+00
## 18:  FC16_134      total_without_group  0.000000e+00               0.000000e+00
## 19:  FC12_160      total_without_group  0.000000e+00               0.000000e+00
## 20:  FC12_143      total_without_group  0.000000e+00               0.000000e+00
## 21:   FC12_69      total_without_group  0.000000e+00               0.000000e+00
## 22:  FC12_118      total_without_group  0.000000e+00               0.000000e+00
## 23:   FC12_44      total_without_group  0.000000e+00               0.000000e+00
## 24: ALL_EDGES incremental_beyond_group -2.220446e-16              -2.193363e-16
## 25:  FC51_181 incremental_beyond_group  0.000000e+00               0.000000e+00
## 26: FC125_181 incremental_beyond_group  0.000000e+00               0.000000e+00
## 27:  FC48_155 incremental_beyond_group  0.000000e+00               0.000000e+00
## 28: FC122_155 incremental_beyond_group  0.000000e+00               0.000000e+00
## 29:  FC48_156 incremental_beyond_group  0.000000e+00               0.000000e+00
## 30: FC122_156 incremental_beyond_group  0.000000e+00               0.000000e+00
## 31:   FC48_80 incremental_beyond_group  0.000000e+00               0.000000e+00
## 32:  FC80_122 incremental_beyond_group  0.000000e+00               0.000000e+00
## 33:   FC80_89 incremental_beyond_group  0.000000e+00               0.000000e+00
## 34:   FC80_90 incremental_beyond_group  0.000000e+00               0.000000e+00
## 35:   FC89_90 incremental_beyond_group  0.000000e+00               0.000000e+00
## 36: FC155_163 incremental_beyond_group  0.000000e+00               0.000000e+00
## 37: FC155_164 incremental_beyond_group  0.000000e+00               0.000000e+00
## 38:   FC16_43 incremental_beyond_group  0.000000e+00               0.000000e+00
## 39:  FC16_117 incremental_beyond_group  0.000000e+00               0.000000e+00
## 40:   FC16_60 incremental_beyond_group  0.000000e+00               0.000000e+00
## 41:  FC16_134 incremental_beyond_group  0.000000e+00               0.000000e+00
## 42:  FC12_160 incremental_beyond_group  0.000000e+00               0.000000e+00
## 43:  FC12_143 incremental_beyond_group  0.000000e+00               0.000000e+00
## 44:   FC12_69 incremental_beyond_group  0.000000e+00               0.000000e+00
## 45:  FC12_118 incremental_beyond_group  0.000000e+00               0.000000e+00
## 46:   FC12_44 incremental_beyond_group  0.000000e+00               0.000000e+00
##          edge            analysis_mode     delta_mse proportional_mse_reduction
##     bootstrap_ci_low bootstrap_ci_high p_fixed_prediction_permutation
##                <num>             <num>                          <num>
##  1:    -1.137979e-16      1.110223e-16                      1.0000000
##  2:    -2.220446e-16      2.220446e-16                      0.9696061
##  3:    -2.220446e-16      1.110223e-16                      0.9536093
##  4:    -2.220446e-16      0.000000e+00                      0.9998000
##  5:    -2.220446e-16      0.000000e+00                      1.0000000
##  6:    -2.220446e-16      2.220446e-16                      0.9988002
##  7:    -2.220446e-16      3.330669e-16                      0.4341132
##  8:    -2.220446e-16      2.220446e-16                      1.0000000
##  9:    -2.220446e-16      2.220446e-16                      0.9936013
## 10:    -2.220446e-16      2.220446e-16                      0.9250150
## 11:    -2.220446e-16      2.220446e-16                      1.0000000
## 12:    -2.220446e-16      2.220446e-16                      0.9930014
## 13:    -2.220446e-16      0.000000e+00                      1.0000000
## 14:     0.000000e+00      2.220446e-16                      1.0000000
## 15:    -2.220446e-16      0.000000e+00                      0.9946011
## 16:    -2.220446e-16      2.220446e-16                      0.9822036
## 17:     0.000000e+00      2.220446e-16                      1.0000000
## 18:    -2.220446e-16      2.220446e-16                      1.0000000
## 19:    -2.220446e-16      0.000000e+00                      0.9996001
## 20:     0.000000e+00      2.220446e-16                      1.0000000
## 21:    -2.220446e-16      2.220446e-16                      0.9760048
## 22:    -2.220446e-16      2.220446e-16                      0.9996001
## 23:    -1.110223e-16      1.110223e-16                      1.0000000
## 24:    -1.110223e-16      1.110223e-16                      1.0000000
## 25:    -2.220446e-16      2.220446e-16                      0.9908018
## 26:    -1.110223e-16      1.110223e-16                      0.9214157
## 27:    -1.110223e-16      1.110223e-16                      0.9724055
## 28:    -2.220446e-16      0.000000e+00                      1.0000000
## 29:    -2.220446e-16      2.220446e-16                      0.9988002
## 30:    -2.220446e-16      2.220446e-16                      0.9942012
## 31:    -3.330669e-16      2.220446e-16                      1.0000000
## 32:    -2.220446e-16      2.220446e-16                      0.9978004
## 33:    -2.220446e-16      2.220446e-16                      0.8808238
## 34:    -2.220446e-16      2.220446e-16                      1.0000000
## 35:    -2.220446e-16      0.000000e+00                      0.9998000
## 36:    -2.220446e-16      0.000000e+00                      1.0000000
## 37:    -1.110223e-16      2.220446e-16                      0.8024395
## 38:    -2.220446e-16      1.110223e-16                      0.9986003
## 39:    -2.220446e-16      1.110223e-16                      1.0000000
## 40:     0.000000e+00      2.220446e-16                      0.9188162
## 41:     0.000000e+00      2.220446e-16                      0.9890022
## 42:    -2.220446e-16      0.000000e+00                      1.0000000
## 43:     0.000000e+00      2.220446e-16                      1.0000000
## 44:    -2.220446e-16      2.220446e-16                      0.9996001
## 45:    -2.220446e-16      2.220446e-16                      1.0000000
## 46:    -2.220446e-16      0.000000e+00                      0.9706059
##     bootstrap_ci_low bootstrap_ci_high p_fixed_prediction_permutation
##     n_subjects n_edges n_permutations n_bootstrap p_bh_22_edges
##          <int>   <int>          <int>       <int>         <num>
##  1:        103      22           5000        2000            NA
##  2:        103       1           5000        2000             1
##  3:        103       1           5000        2000             1
##  4:        103       1           5000        2000             1
##  5:        103       1           5000        2000             1
##  6:        103       1           5000        2000             1
##  7:        103       1           5000        2000             1
##  8:        103       1           5000        2000             1
##  9:        103       1           5000        2000             1
## 10:        103       1           5000        2000             1
## 11:        103       1           5000        2000             1
## 12:        103       1           5000        2000             1
## 13:        103       1           5000        2000             1
## 14:        103       1           5000        2000             1
## 15:        103       1           5000        2000             1
## 16:        103       1           5000        2000             1
## 17:        103       1           5000        2000             1
## 18:        103       1           5000        2000             1
## 19:        103       1           5000        2000             1
## 20:        103       1           5000        2000             1
## 21:        103       1           5000        2000             1
## 22:        103       1           5000        2000             1
## 23:        103       1           5000        2000             1
## 24:        103      22           5000        2000            NA
## 25:        103       1           5000        2000             1
## 26:        103       1           5000        2000             1
## 27:        103       1           5000        2000             1
## 28:        103       1           5000        2000             1
## 29:        103       1           5000        2000             1
## 30:        103       1           5000        2000             1
## 31:        103       1           5000        2000             1
## 32:        103       1           5000        2000             1
## 33:        103       1           5000        2000             1
## 34:        103       1           5000        2000             1
## 35:        103       1           5000        2000             1
## 36:        103       1           5000        2000             1
## 37:        103       1           5000        2000             1
## 38:        103       1           5000        2000             1
## 39:        103       1           5000        2000             1
## 40:        103       1           5000        2000             1
## 41:        103       1           5000        2000             1
## 42:        103       1           5000        2000             1
## 43:        103       1           5000        2000             1
## 44:        103       1           5000        2000             1
## 45:        103       1           5000        2000             1
## 46:        103       1           5000        2000             1
##     n_subjects n_edges n_permutations n_bootstrap p_bh_22_edges
##                                                       edge_label
##                                                           <char>
##  1:                                                    ALL_EDGES
##  2:        ctx_lh_G_cingul-Post-ventral -- ctx_rh_S_pericallosal
##  3:        ctx_rh_G_cingul-Post-ventral -- ctx_rh_S_pericallosal
##  4: ctx_lh_G_and_S_cingul-Mid-Ant -- ctx_rh_Lat_Fis-ant-Vertical
##  5: ctx_rh_G_and_S_cingul-Mid-Ant -- ctx_rh_Lat_Fis-ant-Vertical
##  6:         ctx_lh_G_and_S_cingul-Mid-Ant -- ctx_rh_Lat_Fis-post
##  7:         ctx_rh_G_and_S_cingul-Mid-Ant -- ctx_rh_Lat_Fis-post
##  8: ctx_lh_G_and_S_cingul-Mid-Ant -- ctx_lh_Lat_Fis-ant-Horizont
##  9: ctx_rh_G_and_S_cingul-Mid-Ant -- ctx_lh_Lat_Fis-ant-Horizont
## 10:  ctx_lh_Lat_Fis-ant-Horizont -- ctx_lh_S_circular_insula_inf
## 11:  ctx_lh_Lat_Fis-ant-Horizont -- ctx_lh_S_circular_insula_sup
## 12: ctx_lh_S_circular_insula_inf -- ctx_lh_S_circular_insula_sup
## 13:  ctx_rh_Lat_Fis-ant-Vertical -- ctx_rh_S_circular_insula_inf
## 14:  ctx_rh_Lat_Fis-ant-Vertical -- ctx_rh_S_circular_insula_sup
## 15:          Left-Accumbens-area -- ctx_lh_G_and_S_occipital_inf
## 16:          Left-Accumbens-area -- ctx_rh_G_and_S_occipital_inf
## 17:             Left-Accumbens-area -- ctx_lh_G_occipital_middle
## 18:             Left-Accumbens-area -- ctx_rh_G_occipital_middle
## 19:                               Brain-Stem -- ctx_rh_S_central
## 20:                           Brain-Stem -- ctx_rh_G_postcentral
## 21:                           Brain-Stem -- ctx_lh_G_postcentral
## 22:                     Brain-Stem -- ctx_rh_G_and_S_paracentral
## 23:                     Brain-Stem -- ctx_lh_G_and_S_paracentral
## 24:                                                    ALL_EDGES
## 25:        ctx_lh_G_cingul-Post-ventral -- ctx_rh_S_pericallosal
## 26:        ctx_rh_G_cingul-Post-ventral -- ctx_rh_S_pericallosal
## 27: ctx_lh_G_and_S_cingul-Mid-Ant -- ctx_rh_Lat_Fis-ant-Vertical
## 28: ctx_rh_G_and_S_cingul-Mid-Ant -- ctx_rh_Lat_Fis-ant-Vertical
## 29:         ctx_lh_G_and_S_cingul-Mid-Ant -- ctx_rh_Lat_Fis-post
## 30:         ctx_rh_G_and_S_cingul-Mid-Ant -- ctx_rh_Lat_Fis-post
## 31: ctx_lh_G_and_S_cingul-Mid-Ant -- ctx_lh_Lat_Fis-ant-Horizont
## 32: ctx_rh_G_and_S_cingul-Mid-Ant -- ctx_lh_Lat_Fis-ant-Horizont
## 33:  ctx_lh_Lat_Fis-ant-Horizont -- ctx_lh_S_circular_insula_inf
## 34:  ctx_lh_Lat_Fis-ant-Horizont -- ctx_lh_S_circular_insula_sup
## 35: ctx_lh_S_circular_insula_inf -- ctx_lh_S_circular_insula_sup
## 36:  ctx_rh_Lat_Fis-ant-Vertical -- ctx_rh_S_circular_insula_inf
## 37:  ctx_rh_Lat_Fis-ant-Vertical -- ctx_rh_S_circular_insula_sup
## 38:          Left-Accumbens-area -- ctx_lh_G_and_S_occipital_inf
## 39:          Left-Accumbens-area -- ctx_rh_G_and_S_occipital_inf
## 40:             Left-Accumbens-area -- ctx_lh_G_occipital_middle
## 41:             Left-Accumbens-area -- ctx_rh_G_occipital_middle
## 42:                               Brain-Stem -- ctx_rh_S_central
## 43:                           Brain-Stem -- ctx_rh_G_postcentral
## 44:                           Brain-Stem -- ctx_lh_G_postcentral
## 45:                     Brain-Stem -- ctx_rh_G_and_S_paracentral
## 46:                     Brain-Stem -- ctx_lh_G_and_S_paracentral
##                                                       edge_label
##                                axis              subfamily
##                              <char>                 <char>
##  1:                            <NA>                   <NA>
##  2:      Axis1_Cortical_Integration       posterior_medial
##  3:      Axis1_Cortical_Integration       posterior_medial
##  4:      Axis1_Cortical_Integration    cingulo_perisylvian
##  5:      Axis1_Cortical_Integration    cingulo_perisylvian
##  6:      Axis1_Cortical_Integration    cingulo_perisylvian
##  7:      Axis1_Cortical_Integration    cingulo_perisylvian
##  8:      Axis1_Cortical_Integration    cingulo_perisylvian
##  9:      Axis1_Cortical_Integration    cingulo_perisylvian
## 10:      Axis1_Cortical_Integration      insular_opercular
## 11:      Axis1_Cortical_Integration      insular_opercular
## 12:      Axis1_Cortical_Integration      insular_opercular
## 13:      Axis1_Cortical_Integration      insular_opercular
## 14:      Axis1_Cortical_Integration      insular_opercular
## 15: Axis2_Perceptual_Motor_Coupling          reward_visual
## 16: Axis2_Perceptual_Motor_Coupling          reward_visual
## 17: Axis2_Perceptual_Motor_Coupling          reward_visual
## 18: Axis2_Perceptual_Motor_Coupling          reward_visual
## 19: Axis2_Perceptual_Motor_Coupling brainstem_sensorimotor
## 20: Axis2_Perceptual_Motor_Coupling brainstem_sensorimotor
## 21: Axis2_Perceptual_Motor_Coupling brainstem_sensorimotor
## 22: Axis2_Perceptual_Motor_Coupling brainstem_sensorimotor
## 23: Axis2_Perceptual_Motor_Coupling brainstem_sensorimotor
## 24:                            <NA>                   <NA>
## 25:      Axis1_Cortical_Integration       posterior_medial
## 26:      Axis1_Cortical_Integration       posterior_medial
## 27:      Axis1_Cortical_Integration    cingulo_perisylvian
## 28:      Axis1_Cortical_Integration    cingulo_perisylvian
## 29:      Axis1_Cortical_Integration    cingulo_perisylvian
## 30:      Axis1_Cortical_Integration    cingulo_perisylvian
## 31:      Axis1_Cortical_Integration    cingulo_perisylvian
## 32:      Axis1_Cortical_Integration    cingulo_perisylvian
## 33:      Axis1_Cortical_Integration      insular_opercular
## 34:      Axis1_Cortical_Integration      insular_opercular
## 35:      Axis1_Cortical_Integration      insular_opercular
## 36:      Axis1_Cortical_Integration      insular_opercular
## 37:      Axis1_Cortical_Integration      insular_opercular
## 38: Axis2_Perceptual_Motor_Coupling          reward_visual
## 39: Axis2_Perceptual_Motor_Coupling          reward_visual
## 40: Axis2_Perceptual_Motor_Coupling          reward_visual
## 41: Axis2_Perceptual_Motor_Coupling          reward_visual
## 42: Axis2_Perceptual_Motor_Coupling brainstem_sensorimotor
## 43: Axis2_Perceptual_Motor_Coupling brainstem_sensorimotor
## 44: Axis2_Perceptual_Motor_Coupling brainstem_sensorimotor
## 45: Axis2_Perceptual_Motor_Coupling brainstem_sensorimotor
## 46: Axis2_Perceptual_Motor_Coupling brainstem_sensorimotor
##                                axis              subfamily

fwrite(
  prediction_tests,
  file.path(output_dir, "eye_to_fmri_prediction_permutation_and_bootstrap.csv")
)

11 Eye-temporal predictor stability

The coefficients below are conditional multivariate predictive weights. They should not be interpreted as isolated causal effects of a time-course component on an FC edge.

model_predictor_stability <- all_weights[, .(
  any_edge_selected = any(abs(weight) > 1e-10),
  selected_edge_fraction = mean(abs(weight) > 1e-10),
  feature_l2_norm = feature_l2_norm[1L]
), by = .(analysis_mode, repeat_id, outer_fold, feature)]

predictor_stability <- model_predictor_stability[, .(
  outer_model_selection_frequency = mean(any_edge_selected),
  mean_selected_edge_fraction = mean(selected_edge_fraction),
  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 = .(analysis_mode, feature)]

predictor_stability[, condition := sub("__b[0-9]+$", "", feature)]
predictor_stability[, basis_index := as.integer(sub("^.*__b", "", feature))]

edge_specific_weights <- all_weights[, .(
  selection_frequency = mean(abs(weight) > 1e-10),
  mean_weight = mean(weight),
  mean_abs_weight = mean(abs(weight)),
  median_weight = median(weight)
), by = .(analysis_mode, outcome, feature)]
setnames(edge_specific_weights, "outcome", "edge")
edge_specific_weights <- merge(
  edge_specific_weights,
  edge_labels,
  by = "edge",
  all.x = TRUE,
  sort = FALSE
)

fwrite(
  predictor_stability,
  file.path(output_dir, "eye_temporal_predictor_stability.csv")
)
fwrite(
  edge_specific_weights,
  file.path(output_dir, "eye_to_fmri_edge_specific_weights.csv")
)

plot_predictors <- predictor_stability[
  analysis_mode == "incremental_beyond_group"
]
setorder(plot_predictors, -mean_l2_norm)
plot_predictors <- head(plot_predictors, 20L)
plot_predictors[, feature := factor(
  feature,
  levels = unique(rev(as.character(feature)))
)]

if (nrow(plot_predictors) > 0L) {
  p_weights <- ggplot(plot_predictors, aes(x = mean_l2_norm, y = feature)) +
    geom_point() +
    geom_segment(aes(
      x = l2_norm_q25,
      xend = l2_norm_q75,
      y = feature,
      yend = feature
    )) +
    theme_classic() +
    labs(
      title = "Eye-temporal predictor stability for FC prediction beyond diagnosis",
      x = "Mean L2 norm across 22 edge outcomes and outer models",
      y = NULL
    )
  print(p_weights)
  ggsave(
    file.path(output_dir, "eye_temporal_predictor_stability.png"),
    p_weights,
    width = 11,
    height = 8,
    dpi = 300
  )
}

12 Null-model diagnostic

Because the original fMRI-to-eye model selected no incremental predictive signal, this diagnostic explicitly reports how often the reverse model also collapses to an all-zero eye increment. This is useful for distinguishing “poor but nonzero” prediction from a regularized null solution.

outer_model_nonzero <- all_weights[, .(
  n_nonzero_coefficients = sum(abs(weight) > 1e-10),
  any_nonzero = any(abs(weight) > 1e-10)
), by = .(analysis_mode, repeat_id, outer_fold)]

null_model_summary <- outer_model_nonzero[, .(
  n_outer_models = .N,
  n_all_zero_models = sum(!any_nonzero),
  proportion_all_zero_models = mean(!any_nonzero),
  median_nonzero_coefficients = median(n_nonzero_coefficients)
), by = analysis_mode]

print(null_model_summary)

##               analysis_mode n_outer_models n_all_zero_models
##                      <char>          <int>             <int>
## 1:      total_without_group             50                50
## 2: incremental_beyond_group             50                50
##    proportion_all_zero_models median_nonzero_coefficients
##                         <num>                       <num>
## 1:                          1                           0
## 2:                          1                           0

fwrite(
  null_model_summary,
  file.path(output_dir, "eye_to_fmri_null_model_diagnostic.csv")
)

13 Interpretation guide

1. Matched sample first. Report the exact number of participants contributing both eligible eye trajectories and all 22 targeted FC edges.
2. Primary result. Use the incremental_beyond_group / ALL_EDGES permutation row as the main test of whether temporal gaze adds held-out multivariate FC prediction beyond diagnosis and covariates.
3. Edge-level results. Treat individual-edge tests as secondary and use the BH-adjusted p values across the 22 edges.
4. Group-structure result. If total_without_group predicts but incremental_beyond_group does not, the predictive association is largely attributable to shared MDD-HC structure rather than participant-level cross-modal coupling beyond diagnosis.
5. Null shrinkage matters. If many outer models are all-zero under the one-SE rule, report that the nested CV preferred the nuisance-only model. Do not switch to weaker regularization solely to obtain nonzero coefficients.
6. Direction is not causality. Eye-to-fMRI prediction does not establish that gaze changes resting-state connectivity.
7. External validity. Repeated nested CV reduces leakage but does not replace independent replication.

14 Reproducibility metadata

settings_table <- data.table(
  setting = c(
    "data_dir",
    "prediction_direction",
    "eye_condition_set",
    "primary_eye_conditions",
    "eye_basis_df",
    "n_eye_predictors",
    "n_fmri_outcome_edges",
    "minimum_eye_bins_per_condition",
    "outer_folds",
    "outer_repeats",
    "inner_folds",
    "alpha_grid",
    "lambda_selection_rule",
    "motion_available_fraction",
    "motion_used",
    "include_icv_sensitivity",
    "prediction_permutations",
    "prediction_bootstraps",
    "targets_independently_prespecified",
    "random_seed"
  ),
  value = c(
    data_dir,
    "eye_temporal_to_fmri_targeted_edges22",
    eye_condition_set,
    paste(primary_eye_conditions, collapse = ";"),
    eye_basis_df,
    length(predictor_cols),
    length(target_edge_cols),
    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,
    prediction_permutations,
    prediction_bootstraps,
    targets_independently_prespecified,
    random_seed
  )
)
fwrite(settings_table, file.path(output_dir, "eye_to_fmri_mvpa_settings.csv"))

writeLines(
  c(
    paste0("Created: ", Sys.time()),
    "Prediction direction: eye temporal basis -> 22 targeted fMRI edges",
    paste0("Matched participants: ", nrow(analysis_data)),
    paste0("HC: ", sum(analysis_data$Group == "HC")),
    paste0("MDD: ", sum(analysis_data$Group == "MDD")),
    paste0("Eye predictors: ", length(predictor_cols)),
    paste0("fMRI outcomes: ", length(target_edge_cols)),
    paste0("Motion used: ", use_motion),
    paste0("Output directory: ", output_dir),
    "",
    capture.output(sessionInfo())
  ),
  file.path(output_dir, "eye_to_fmri_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