crossnobis_plan.md

In bbuchsbaum/rMVPA: Multivoxel Pattern Analysis in R

Implementation Plan – Crossnobis Support for MS-ReVE / Contrast-RSA

0 Paper Summary & Motivation

The cross-nobis (a.k.a. cross-validated Mahalanobis) distance provides an unbiased estimate of the true squared pattern distance by multiplying pattern differences that were estimated from independent data partitions. Ordinary (same-fold) distances are positively biased because noise is multiplied by itself.

Formal definition for condition pair (k=(i,j)) with (M) partitions (Eq. 5, Diedrichsen & Kriegeskorte 2017): [ \tilde d_k = \frac{1}{M(M-1)\,P}\sum_{m\neq n} \hat\delta_{k,m}^{\mathsf T}\,\hat\delta_{k,n}, \qquad \hat\delta_{k,m}=\hat b_{i,m}-\hat b_{j,m}, ] where (P) is the number of voxels/channels. Equivalent "single-fold-vs-rest" algebra can also be used.

Key properties Unbiased: (\mathbb E[\tilde d_k]=\delta_k^{\mathsf T}\delta_k/P). Can be negative when the true distance is ~0 – keep these values, do not truncate. Variance is (\tfrac{M}{M-1}) times larger than the biased estimator but decreases with more partitions. Whitening: pre-multiply patterns by (\mathbf W=\Sigma_{\text{noise}}^{-1/2}) to obtain Mahalanobis distances; crucial for fMRI.

(Sections 1–9 below incorporate this theory.)

1 High-Level Requirements

Produce unbiased squared distances … (unchanged from previous draft)

// ... existing content of sections 1–9 from crossnobis.md retained unchanged ...

9 Timeline (suggestion)

| Week | Deliverable | |------|-------------| | 1 | Helper function + unit tests pass | | 2 | Integration into train_model + searchlight plumbing | | 3 | Documentation, vignette, example rebuild | | 4 | Profiling & optional Rcpp optimisation |

Appendix A Detailed Implementation Table & Sanity-Check Code

A.1 Step-wise Implementation Table (adapted from user note)

| Step | What to compute | Suggested R objects / code changes | |------|-----------------|------------------------------------| | A | Return fold-wise condition means. | compute_crossvalidated_means_sl(return_folds = TRUE) returns a list with mean = U_hat_sl, folds = U_folds (K × V × M). | | B | Pre-whiten each fold. | Within the same helper: if W provided, U_folds[,,m] %*% W. | C | Cross-fold inner products only for m ≠ n. | New helper crossnobis_distance(U_folds, P) computes distance vector; pseudo-code shown below. | | D | Expose to RSA regression. | train_model.contrast_rsa_model() uses this distance vector when estimation_method == "crossnobis". include_vec still applied. | | E | Allow negative outputs. | No truncation; document in Roxygen. | | F | Unit tests. | Noise-only simulation: mean distance ≈ 0; biased Euclidean positive. | | G | Performance tips. | Vectorise across pairs; optionally migrate to Rcpp. |

crossnobis_distance <- function(U_folds) {
  K <- dim(U_folds)[1]; V <- dim(U_folds)[2]; M <- dim(U_folds)[3]
  if (M < 2) rlang::abort("Need at least two partitions for cross-nobis.")
  pair_mat <- utils::combn(K, 2)
  D <- ncol(pair_mat)
  out <- numeric(D)
  P <- V
  for (p in seq_len(D)) {
    i <- pair_mat[1, p]; j <- pair_mat[2, p]
    deltas <- U_folds[i, , ] - U_folds[j, , ]   # V × M
    ip <- tcrossprod(t(deltas))                 # M × M inner products
    diag(ip) <- 0
    out[p] <- sum(ip) / (P * M * (M - 1))
  }
  names(out) <- paste0(dimnames(U_folds)[[1]][pair_mat[1, ]], "_vs_", 
                       dimnames(U_folds)[[1]][pair_mat[2, ]])
  out
}

A.2 Quick Sanity-Check

set.seed(1)
K <- 4; V <- 30; M <- 8
true <- matrix(rnorm(K*V), K, V)          # ground-truth patterns
noise <- array(rnorm(K*V*M, sd = 3), c(K, V, M))
U_folds <- true + noise                   # simulated noisy means per fold
cross_d <- crossnobis_distance(U_folds)
mean(cross_d)   # ≈ 0 for pure noise case

10 Checkable Implementation Tickets

This section breaks down the implementation plan into granular, checkable tasks for the specified R files and the new helper function.

Always check off each task as you complete it.

New Helper Function (e.g., in R/crossnobis_helpers.R or similar)

• [ ] Implement compute_crossnobis_distances_sl function:
  • [ ] Define function signature: compute_crossnobis_distances_sl(U_folds, P_voxels).
    • U_folds: K x V x M array of condition means per fold (K=conditions, V=voxels, M=folds).
    • P_voxels: Number of voxels (V).
  • [ ] Handle edge case: If M < 2, abort with an informative message.
  • [ ] Generate all unique unordered condition pairs (e.g., using utils::combn(K, 2)).
  • [ ] For each condition pair k = (i,j):
    • [ ] Initialize sum_cross_prods_k = 0.
    • [ ] For each pair of distinct folds (m, n) where m != n:
      • [ ] Calculate delta_k_m = U_folds[i,,m] - U_folds[j,,m].
      • [ ] Calculate delta_k_n = U_folds[i,,n] - U_folds[j,,n].
      • [ ] Add crossprod(delta_k_m, delta_k_n) to sum_cross_prods_k.
    • [ ] Calculate d_crossnobis_k = sum_cross_prods_k / (P_voxels * M * (M - 1)).
  • [ ] Return a named numeric vector of these d_crossnobis_k values, with names like "condA_vs_condB", ordered to match lower-triangle vectorization.
  • [ ] Handle potential NAs in U_folds (e.g., if a condition is missing in a fold): ensure resulting distances involving such conditions are NA, or relevant pairs are skipped.
  • [ ] Add Roxygen documentation for this new function.

R/compute_cv_means.R (Changes to compute_crossvalidated_means_sl)

• [ ] Argument and Return Value Modification:
  • [ ] Add a new boolean argument return_folds = FALSE to compute_crossvalidated_means_sl.
  • [ ] If return_folds = TRUE, change the return value to a list: list(mean_estimate = U_hat_sl, fold_estimates = U_folds_array), where U_folds_array is (K x V x M).
  • [ ] If return_folds = FALSE (default), return U_hat_sl as before (or mean_estimate from the list).
• [ ] Per-Fold Mean Calculation (when return_folds = TRUE):
  • [ ] Inside the loop over folds, store the fold_means_mat_processed (K_fold x V matrix for conditions present in the current fold) before re-indexing to full_fold_means.
  • [ ] Collect these fold_means_mat_processed (or their full_fold_means equivalent containing NAs for absent conditions) into the U_folds_array (K x V x M).
• [ ] Whitening Logic (when return_folds = TRUE and whitening_matrix_W is provided):
  • [ ] Ensure that if estimation_method == "crossnobis" (as passed to compute_cv_means) and whitening_matrix_W is provided, the whitening (fold_means_mat %*% whitening_matrix_W) is applied to the patterns that will form U_folds_array. (Current logic already does this for fold_means_mat_processed).
• [ ] Documentation:
  • [ ] Update Roxygen for compute_crossvalidated_means_sl detailing return_folds and the new list return type.

R/contrast_rsa_model.R

• contrast_rsa_model Constructor:

  • [ ] Re-enable estimation_method = "crossnobis" in match.arg(). (Remove the rlang::abort for it).
  • [ ] Add argument whitening_matrix_W = NULL to the function signature.
  • [ ] Store whitening_matrix_W in the model specification object (e.g., x$whitening_matrix_W <- whitening_matrix_W).
  • [ ] Update Roxygen:
    • [ ] Reflect "crossnobis" as a valid estimation_method.
    • [ ] Document the new whitening_matrix_W argument, its purpose for crossnobis, and that it's passed to compute_crossvalidated_means_sl.

• train_model.contrast_rsa_model Method:

  • [ ] Initialization:
    • [ ] Declare U_hat_for_delta_calc (will hold the K x V matrix for Δ projections).
    • [ ] Declare dvec_sl (will hold the vector of (dis)similarities for regression).
  • [ ] Conditional Logic for estimation_method:
    • [ ] If obj$estimation_method == "crossnobis":
      • [ ] Call compute_crossvalidated_means_sl with:
        • estimation_method = "crossnobis" (this tells compute_cv_means to handle whitening internally if W is passed).
        • whitening_matrix_W = obj$whitening_matrix_W.
        • return_folds = TRUE.
        • Store result in cv_outputs.
      • [ ] Assign U_hat_for_delta_calc <- cv_outputs$mean_estimate.
      • [ ] Let U_folds_data <- cv_outputs$fold_estimates.
      • [ ] P_voxels <- ncol(sl_data).
      • [ ] dvec_sl <- compute_crossnobis_distances_sl(U_folds_data, P_voxels).
    • [ ] Else (for "average", "L2_norm"):
      • [ ] U_hat_for_delta_calc <- compute_crossvalidated_means_sl(sl_data, mvpa_des, cv_spec, obj$estimation_method, whitening_matrix_W = NULL) (Note: compute_cv_means correctly ignores W if its own estimation_method is not "crossnobis").
      • [ ] G_hat_sl <- U_hat_for_delta_calc %*% t(U_hat_for_delta_calc).
      • [ ] dvec_sl <- G_hat_sl[lower.tri(G_hat_sl)].
  • [ ] Common Post-Processing for dvec_sl (after if/else block):
    • [ ] If obj$estimation_method != "crossnobis", apply include_vec logic to dvec_sl based on G_hat_sl as currently done.
    • [ ] If obj$estimation_method == "crossnobis", apply include_vec logic directly to the dvec_sl output by compute_crossnobis_distances_sl (ensure mapping from KxK is_intra_block_pair to K(K-1)/2 dvec_sl elements is correct).
  • [ ] Voxel Projections (Δ_sl):
    • [ ] Calculate Delta_sl <- t(U_hat_for_delta_calc) %*% C_ord using the consistently derived U_hat_for_delta_calc.
  • [ ] (Rest of the RSA regression, metric calculation proceeds using dvec_sl and beta_sl, and Delta_sl as appropriate).
• run_searchlight.contrast_rsa_model Method:
  • [ ] No specific changes required here if whitening_matrix_W is stored in and accessed from model_spec by train_model.

R/msreve_design.R

• [ ] No changes anticipated for this file based on the current crossnobis_plan.md.

bbuchsbaum/rMVPA documentation built on June 10, 2025, 8:23 p.m.