Draft new dependence metrics

[tsalo]

Closes none, but adds a method that we discussed in our most recent devs call. The two functions that should be reviewed are compute_te_variance and compute_te_variance_permutation.

I'm not going to fix the unit tests unless we decide that these metrics are something we actually want to add to tedana.

Changes proposed in this pull request:

[tsalo]

Here's the Notes section from compute_te_variance_permutation's docstring.

Important: Inference vs. Description

• kappa_star and rho_star are descriptive variance fractions and are not used for statistical inference.
• Statistical inference is based solely on the permutation p-values (p_t2, p_s0), which are empirical tail probabilities from the null distribution.

Null hypothesis

The component's echo-wise parameter estimates show echo-wise structure, but this structure is not specifically aligned with voxel-local T2* or S0 sensitivity. Under the null, permuting which voxel's basis functions are used for fitting should not systematically change model fit.

Permutation scheme

This test uses independent permutations for T2* and S0 to isolate spatial specificity of each component:

• For p_t2: Permute only φ_T2* while keeping φ_S0 fixed (local). This tests whether T2* explains more unique variance with its local basis than with a randomly assigned basis.
• For p_s0: Permute only φ_S0 while keeping φ_T2* fixed (local). This tests whether S0 explains more unique variance with its local basis than with a randomly assigned basis.

Within groups of voxels sharing the same number of valid echoes (n_e), shuffle the assignment of the target basis function to voxels while keeping echo-wise PEs and the other basis fixed. This preserves:

• Total variance structure of each component
• Echo-wise correlation structure within voxels
• Marginal distribution of basis function shapes
• The contribution from the non-permuted basis (isolation of effects)

While breaking:

• Spatial alignment between the target basis and voxel-local sensitivity

Note that permutations are performed within echo-availability strata only. The null is therefore: "alignment no better than random assignment within echo-availability groups."

Test statistics

• ss_t2: Unique variance explained by T2* (Type III sum of squares). This is the variance that T2* explains beyond what S0 explains. Higher values indicate better spatial alignment with local T2*.
• ss_s0: Unique variance explained by S0 (Type III sum of squares). This is the variance that S0 explains beyond what T2* explains. Higher values indicate better spatial alignment with local S0.

P-values are the proportion of permutations achieving unique variance as high or higher than observed (one-tailed test, higher is better).

Interpretation

• Low p_t2: Echo-wise pattern is spatially aligned with local T2* sensitivity → consistent with BOLD-like signal
• Low p_s0: Echo-wise pattern is spatially aligned with local S0 sensitivity → consistent with S0-driven fluctuations
• Both low: Mixed signal with both T2* and S0 spatial specificity
• Both high: Non-specific echo-wise structure

Limitations and caveats

1. Model dependence: This test evaluates alignment conditional on the linearized monoexponential model δS ≈ (∂S/∂S0)δS0 + (∂S/∂T2*)δT2*. If the linearization is poor (large fluctuations, nonlinear effects), the interpretation weakens.

2. Map quality dependence: The test assumes the provided s0_hat and t2s_hat are meaningful voxelwise quantities. If T2*/S0 maps are noisy or biased, the null distribution becomes wider and power decreases. P-values are conditional on map quality.

3. Physiological interpretation: This test identifies components whose echo-wise structure is consistent with voxel-local T2* sensitivity. This is a necessary but not sufficient condition for neuronal BOLD signal. A component could show T2*-alignment due to non-neuronal effects (e.g., motion-correlated susceptibility artifacts).

Why this approach

Traditional kappa/rho metrics and parametric F-tests for component-level TE-dependence suffer from invalid aggregation (averaging ratios) or incorrect degrees of freedom. This permutation test avoids parametric assumptions entirely and directly tests the spatial specificity that distinguishes BOLD from other echo-wise structure.

[handwerkerd]

@tsalo & @eurunuela Just following up on a comment from the dev call & figure this is a reasonable place to share. I was playing a bit with estimating a T2* and R0 fit for a component across all voxels instead of doing the fit on each voxel separate and then calculating some type of weighted sum. I realized that sklearn.linear_mode.RANSACRegressor (manual) might be a useful too for this. It both allowed for a weighted linear regression so that voxels that contribute more to the component are weighted more and it also takes a nice approach to make sure outlier voxels don't overly skew the fit. I didn't love my first results & it's clear that how we weight matters (i.e. ICA weight maps have bright blobs. How much should those blobs be weights vs everything else?). Still, this seemed like a potentially promising direction.