Single-condition spatiotemporal cluster test on source epochs (dSPM): Test against 0 vs. late-window subtraction variants?

Support & Discussions
1

Hi everyone,

I am running a single-participant, single-condition source analysis in MNE-Python, and I want to use a non-parametric spatiotemporal cluster permutation test to find significant cortical activation within a post-stimulus ROI (0.1 to 0.5 seconds), I want to do this test between the source estimate epochs (see the code below) and that to apply the within subject anaylsis at the source level not the channel level.

Because I only have a single experimental condition, I am confused about the most mathematically sound way to test my active ROI against a baseline or a state where there is no stimulus activity. I want to avoid direct comparison with the pre-stimulus baseline (-0.2 to 0 s) because the noise covariance matrix was optimized on that exact window, which may introduces statistical circularity. Instead, I am looking at a late post-stimulus window (e.g., 2.2 to 2.7 seconds) where I expect no stimulus-driven activity to remain.

note: the stimulus is every 5 seconds approximately.

I am considering a few different approaches and want to know which one is statistically and scientifically valid:

  • Approach A (Subtraction): Subtract the late window from the active window at the channel level (epochs(active) - epochs(late)) and feed its source estimate difference array “X” into mne.stats.spatio_temporal_cluster_1samp_test .

  • Approach B (Mean Subtraction): Take the average activity over time for the late window subtract it from the active window feed its source estimate difference array “X” into mne.stats.spatio_temporal_cluster_1samp_test .

  • Approach C (Test Against Zero, shown in the code below): Ignore the late window entirely, take the estimated source epochs, and feed them directly into mne.stats.spatio_temporal_cluster_1samp_test to evaluate whether the active ROI significantly deviates from 0.

I would pick Approach C but I am wondering from the scientific statistical perspective whihc one is best or maybe some other approach can be suggested.

fwd_general = make_forward_solution(
    epochs.info,
    trans=str(trans),
    src=str(SRC),
    bem=str(BEM),
    eeg=True,
    mindist=5.0,
    n_jobs=None,
)

epochs_bl = epochs.copy().apply_baseline((-0.2, 0))

cov_matrix = compute_covariance(
    epochs_bl,
    tmin = -0.2,
    tmax = 0,
    method=["empirical", "shrunk", "diagonal_fixed"],
)
inv = make_inverse_operator(
    epochs_bl.average().info,
    fwd_general,
    cov_matrix,
    loose="auto"
)


stc_epochs_dspm = apply_inverse_epochs(epochs=epochs_bl,
                                            inverse_operator=inv,
                                            lambda2=LAMBDA2,
                                            method="dSPM",
                                            pick_ori="normal")

src_ = read_source_spaces(SRC)
adjacency = spatial_src_adjacency(src_)

# time window for the "active" task (expected time of activity)
tmin_active = 0.1
tmax_active = 0.5

n_epochs = len(stc_epochs_dspm)

diff_data = []

for stc in stc_epochs_dspm:
    
    stc_act = stc.copy().crop(tmin=tmin_active, tmax=tmax_active)
    act = stc_act.data
    
    diff = act # - base_mean_over_time
    
    diff_transposed = diff.T
    
    diff_data.append(diff_transposed)
    
X = np.array(diff_data)
    
print("Running cluster permutation test ...")
T_obs, clusters, cluster_p_values, H0 = spatio_temporal_cluster_1samp_test(
    X,
    adjacency=adjacency,
    n_permutations=1000, # max_perms = 2 ** (n_samples - (tail == 0)) - 1
    threshold=None,
    tail=0, 
    n_jobs=-1
)
2

Option C is definitely wrong. In MNE source estimation, there is non-zero activity at each source point, at all times. A null-hypothesis of there being 0 activity can never be true.

I’m afraid you will have to compare the amount of activity during your time window of interest, with the amount of activity during some other time. The baseline window seems like the most obvious choice to me, but your concern that the noise covariance matrix is also estimated on the baseline window may be valid (I’m not entirely sure). Could the noise covariance be computed over some period of rest, somewhere else in the recording?

3

Thanks for the reply!

Do you mean appraoch A or B? if you mean A, which one is the best way to do it:

  1. diff = channel_data(0.1->0.5) - channel_data(2.2->2.7) → baseline(diff) → Estimate sources → apply the test

or

  1. diff = baseline(channel_data(0.1->0.5)) - baseline(channel_data(2.2->2.7)) → Estimate sources → apply the test

In the recorded data, there is an auditory stimulus every 5 seconds (~300ms stimulus duration)
Maybe I can use other time window(e.g: 3s-3.5s post stimulus) but unfortuently I don’t have empty room recording or other period of rest.

4

option 2 makes more sense to me than option 1.