Spatiotemporal permutation F-test - No sign. Clusters found

💬 Support & Discussions
1

- MNE version: 1.9
- operating system: Windows 11
- 64 electrodes used (data already preprocessed: ICA, Bandpass filtered, Re-referenced.)

Dear MNE family,

I have been exploring the spatio-temporal permutation F-test in MNE-Python, following the well-written tutorial. I am using EEG data from a previous experiment with two conditions: standing and walking. In the experiment, participants stood for 1 minute to establish a baseline, then walked. I created 30 events for each condition, spaced 2 seconds apart, with each event lasting 1.5 seconds. I then extracted epochs for these conditions and applied the spatio-temporal permutation F-test as outlined in the tutorial.

montage = mne.channels.make_standard_montage("easycap-M1")
raw.set_montage(montage)
sfreq = raw.info['sfreq'] 
events, event_id = mne.events_from_annotations(raw)

start_standing = event_id['Photosensor 255']
start_walking = event_id['Photosensor 2']

event_samples = events[:, 0]
event_types = events[:, 2] 

marker_sample_start_standing = event_samples[event_types == start_standing][0]
marker_sample_start_walking = event_samples[event_types == start_walking][0]

marker_time_standing = marker_sample_start_standing / sfreq
marker_time_walking = marker_sample_start_walking / sfreq

n_annotations = 30
spacing = 2 

def create_annotations(timepoint, label):
    onsets = np.array([timepoint + i * spacing for i in range(n_annotations)])
    durations = np.zeros(n_annotations) 
    descriptions = [label] * n_annotations
    return mne.Annotations(onset=onsets, duration=durations, description=descriptions)

annotations_standing = create_annotations(marker_time_standing, 'Standing')
annotations_walking = create_annotations(marker_time_walking, 'Walking')

combined_annotations = annotations_standing + annotations_walking

raw.set_annotations(combined_annotations)

events, event_id = mne.events_from_annotations(raw)

tmin = 0.0
tmax = 1.5 

epochs = mne.Epochs(raw, 
                        events,
                        event_id,
                        tmin=tmin,
                        tmax=tmax,
                        baseline=None,
                        preload=True,
                        picks='eeg')

Until here everything works fine.


X = [epochs[event_name].get_data(copy=False) for event_name in event_id]
X = [np.transpose(x, (0, 2, 1)) for x in X]


adjacency, ch_names = find_ch_adjacency(epochs.info, ch_type="eeg")
mne.viz.plot_ch_adjacency(epochs.info, adjacency, ch_names, kind="3d")

alpha_cluster_forming = 0.05
tail = 1


n_conditions = len(event_id)
n_observations = len(X[0])
dfn = n_conditions - 1
dfd = n_observations - n_conditions


f_thresh = scipy.stats.f.ppf(1 - alpha_cluster_forming, dfn=dfn, dfd=dfd)

cluster_stats = spatio_temporal_cluster_test(
    X,
    n_permutations=1000,
    threshold=f_thresh,
    tail=tail,
    n_jobs=None,
    buffer_size=None,
    adjacency=adjacency,
)
F_obs, clusters, p_values, _ = cluster_stats

X is a list, containing two arrays with a size of 30, 751, 65. Looking fine and being in line with the documentation. Adjacency also looks regular.

The permutations are done, and everything seems to be working. I am finding 377 clusters. But I’m puzzled by the p-values. They are consistently either 1 or very close to 1 (e.g., 0.999). I’m new to this field of data analysis, but I would expect some differences between the standing and walking conditions. I have experimented with various values regarding thresholds, number of permutations, one-sided vs. two-sided tails, different tmin and tmax settings for the epochs and varying the number of epochs, but still no significant p-values.

Where could I have gone wrong, or what might I be missing?

Thank you! :slight_smile:

Cheers,
Tim

2

The way to proceed is to plot the raw f-values. If you can somehow get a matrix of size (channels x times) filled with the f-values across your conditions, you can load it into an EvokedArray object and plot it:

fvalues = scipy.stats.f_oneway(X)  # check that is (channels x times)
f_evoked = mne.EvokedArray(fvalues, info=epochs.info, tmin=epochs.times[0])
f_evoked.plot_topo()

…and see if these f-values makes sense to you: high values where you indeed expect differences and low value here you don’t expect any differences.

3 Likes
3

Dear Marijn,

thank you for your reply and sorry for my late response.
I did check the raw f-values as you mentioned.

fvalues = scipy.stats.f_oneway(*X)
X_reshaped = [x.transpose(0, 2, 1) for x in X]  # Now (epochs, channels, timepoints)
fvalues = np.array([scipy.stats.f_oneway(*[X_reshaped[c][:, ch, :] for c in range(len(X))])[0]
                    for ch in range(X_reshaped[0].shape[1])])

f_evoked = mne.EvokedArray(fvalues, info=epochs.info, tmin=epochs.times[0])
f_evoked.plot_topo()

The f-values are looking quiet similiar and I would expect higher values in central regions.
Nevertheless, if I take a look at critical F-values:

f_threshold = scipy.stats.f.ppf(0.95, df_between, df_within)
print(f"F-value threshold for p < 0.05: {f_threshold}")

I do find several. However, not in the regions I would expect them to be.

I now tried to align the ‘Walking’ annotation to the heel strike. This at least led to significant clusters being found that would make sense. Still have to figure out, if this is only driven by artefacts. But atleast it worked somehow :-).

1 Like