Hello everyone,
I am analyzing MEG time-frequency data and I am trying to understand a mismatch between sensor-space and source-space group statistics.
In sensor space, I observe positive significant clusters over frontal sensors for some contrasts and frequency bands. For example, for the contrast Condition A - Condition B, I see positive frontal sensor-level clusters in the alpha and beta bands.
However, in source space, the significant clusters are consistently visualized as negative/blue clusters, and I do not see corresponding positive/red clusters in alpha and beta bands. As an example, I am attaching images of comparisons, questions to the first image.
I would like to understand whether this reflects a real difference between sensor-level and source-level statistics, or whether something in my source-space plotting/statistics pipeline may be wrong.
For sensor-space analysis, the pipeline is:
- Load preprocessed MEG data and keep gradiometers only.
- Create condition-specific epochs from event markers. The epochs used for the CSD analysis are correct-response epochs.
- Compute CSDs in sensor space using Morlet wavelets:
csd = mne.time_frequency.csd_morlet(
epochs,
freqs,
tmin=tmin,
tmax=tmax,
n_cycles=n_cycles,
decim=decim,
n_jobs=n_jobs,
)
- For each subject and contrast, compute condition-specific CSDs in the ROI time window and a paired baseline CSD time window.
For example, for CoVi - SiVi:
A = CoVi Storage2 CSD
B = SiVi Storage2 CSD
baseline = CoVi_SiVi Storage1 paired-baseline CSD
The individual normalized sensor contrast is:
normalized_csd = (A - B) / paired_baseline
- For group-level sensor statistics, the already-computed individual normalized CSD contrast files are loaded:
S1_CoVi_vs_SiVi_normalized_csd.npy
S2_CoVi_vs_SiVi_normalized_csd.npy
...
Each file has shape:
channels × frequencies
The subject-level contrast arrays are stacked and transposed before statistics:
p = list_norm_csd # subjects × channels × frequencies
obj = np.transpose(p, (0, 2, 1)) # subjects × frequencies × channels
Sensor adjacency is computed from a representative epochs object using gradiometers:
adj, ch_names = mne.channels.find_ch_adjacency(info, ch_type="grad")
Group-level cluster statistics are performed as a one-sample test against zero:
T_obs, clusters, cluster_p_values, H0 = mne.stats.spatio_temporal_cluster_1samp_test(
obj,
out_type=out_type,
adjacency=adj,
n_permutations=n_permutations,
threshold=t_threshold,
tail=tail,
)
So, at the sensor level, the group analysis does not subtract condition B from condition A again. It tests whether the individual normalized CSD contrasts differ from zero across subjects.
For the group-level source-space analysis, I use morphed source estimates that were already computed for each subject.
At the group level, I load three morphed STC files per subject:
stc_s # condition A
stc_t # condition B
stc_a # paired baseline
For example, for CoVi_vs_SiVi, the intended contrast is:
CoVi - SiVi
The normalized source contrast is computed as:
group_1 = np.array([stc.data for stc in stc_s])
group_2 = np.array([stc.data for stc in stc_t])
group_a = np.array([stc.data for stc in stc_a])
diff = (group_1 - group_2) / group_a
STAT = np.transpose(diff, [0, 2, 1]) # subjects × time/frequency × sources
Then I run a one-sample source-space cluster permutation test against zero:
adjacency = mne.spatial_src_adjacency(src_fs)
T_obs, clusters, cluster_p_values, H0 = spatio_temporal_cluster_1samp_test(
STAT,
adjacency=adjacency,
n_jobs=n_jobs,
threshold=t_threshold,
buffer_size=buffer_size,
verbose=True,
n_permutations=n_permutations,
out_type=out_type,
)
For the final figures shown here, I use a signed T-value visualization of all significant clusters, not the MNE cluster-extent summary.
The function keeps T-values only inside significant clusters:
T_obs_sig = np.zeros_like(T_obs)
for idx in good_cluster_inds:
cluster_mask = _cluster_to_mask(clusters[idx], T_obs.shape)
T_obs_sig[cluster_mask] = T_obs[cluster_mask]
If T_obs has more than one time/frequency point, I collapse it into one source map by selecting, for each vertex, the signed T-value with the maximum absolute magnitude:
best_time = np.argmax(np.abs(T_obs_sig), axis=0)
signed_map = T_obs_sig[best_time, np.arange(T_obs_sig.shape[1])]
data = signed_map[:, np.newaxis]
Then I create a SourceEstimate and plot it using:
stc_signed.plot(
subject="fsaverage",
subjects_dir=subjects_dir,
hemi=hemi,
surface=surface,
views=views,
colormap="RdBu_r",
clim=dict(kind="value", pos_lims=[absmax * 0.4, absmax * 0.7, absmax]),
smoothing_steps=10,
transparent=True,
time_label="signed T-values",
time_viewer=False,
show_traces=False,
background="white",
foreground="black",
)
Given this pipeline, my main questions are:
-
Is it expected that a positive frontal cluster in MEG sensor space does not necessarily appear as a positive frontal cluster in source space?
-
Does the source-space statistical pipeline look conceptually correct: loading morphed STCs for condition A, condition B, and paired baseline, computing
(A - B) / baseline, and then running a one-sample cluster permutation test against zero? -
Is my signed T-value visualization of source clusters appropriate, especially the step where I collapse multiple frequency/time points into one source map by selecting the signed T-value with the maximum absolute magnitude at each vertex?
-
Could my plotting settings, especially
RdBu_r,clim=dict(kind="value", pos_lims=[...]), smoothing, or transparency, cause positive/red clusters to be hidden or underrepresented?
In short, I would like to know whether the absence of positive/red frontal clusters in source space likely reflects the actual source-space statistics, or whether it may be caused by an issue in the source-space visualization/statistics code.
I am running the analysis in Jupyter Notebook with a Python [conda env:mne] kernel.
MNE version: 1.11.0
