Note
Go to the end to download the full example code.
Spatiotemporal permutation F-test on full sensor data#
Tests for differential evoked responses in at least one condition using a permutation clustering test. The FieldTrip neighbor templates will be used to determine the adjacency between sensors. This serves as a spatial prior to the clustering. Spatiotemporal clusters will then be visualized using custom matplotlib code.
Here, the unit of observation is epochs from a specific study subject. However, the same logic applies when the unit observation is a number of study subject each of whom contribute their own averaged data (i.e., an average of their epochs). This would then be considered an analysis at the “2nd level”.
See the FieldTrip tutorial for a caveat regarding the possible interpretation of “significant” clusters.
For more information on cluster-based permutation testing in MNE-Python, see also: Non-parametric 1 sample cluster statistic on single trial power.
# Authors: Denis Engemann <denis.engemann@gmail.com>
# Jona Sassenhagen <jona.sassenhagen@gmail.com>
# Alex Rockhill <aprockhill@mailbox.org>
# Stefan Appelhoff <stefan.appelhoff@mailbox.org>
#
# License: BSD-3-Clause
# Copyright the MNE-Python contributors.
import matplotlib.pyplot as plt
import numpy as np
import scipy.stats
from mpl_toolkits.axes_grid1 import make_axes_locatable
import mne
from mne.channels import find_ch_adjacency
from mne.datasets import sample
from mne.stats import combine_adjacency, spatio_temporal_cluster_test
from mne.viz import plot_compare_evokeds
Set parameters#
data_path = sample.data_path()
meg_path = data_path / "MEG" / "sample"
raw_fname = meg_path / "sample_audvis_filt-0-40_raw.fif"
event_fname = meg_path / "sample_audvis_filt-0-40_raw-eve.fif"
event_id = {"Aud/L": 1, "Aud/R": 2, "Vis/L": 3, "Vis/R": 4}
tmin = -0.2
tmax = 0.5
# Setup for reading the raw data
raw = mne.io.read_raw_fif(raw_fname, preload=True)
raw.filter(1, 25)
events = mne.read_events(event_fname)
Opening raw data file /home/circleci/mne_data/MNE-sample-data/MEG/sample/sample_audvis_filt-0-40_raw.fif...
Read a total of 4 projection items:
PCA-v1 (1 x 102) idle
PCA-v2 (1 x 102) idle
PCA-v3 (1 x 102) idle
Average EEG reference (1 x 60) idle
Range : 6450 ... 48149 = 42.956 ... 320.665 secs
Ready.
Reading 0 ... 41699 = 0.000 ... 277.709 secs...
Filtering raw data in 1 contiguous segment
Setting up band-pass filter from 1 - 25 Hz
FIR filter parameters
---------------------
Designing a one-pass, zero-phase, non-causal bandpass filter:
- Windowed time-domain design (firwin) method
- Hamming window with 0.0194 passband ripple and 53 dB stopband attenuation
- Lower passband edge: 1.00
- Lower transition bandwidth: 1.00 Hz (-6 dB cutoff frequency: 0.50 Hz)
- Upper passband edge: 25.00 Hz
- Upper transition bandwidth: 6.25 Hz (-6 dB cutoff frequency: 28.12 Hz)
- Filter length: 497 samples (3.310 s)
[Parallel(n_jobs=1)]: Done 17 tasks | elapsed: 0.0s
[Parallel(n_jobs=1)]: Done 71 tasks | elapsed: 0.1s
[Parallel(n_jobs=1)]: Done 161 tasks | elapsed: 0.2s
[Parallel(n_jobs=1)]: Done 287 tasks | elapsed: 0.4s
Read epochs for the channel of interest#
picks = mne.pick_types(raw.info, meg="mag", eog=True)
reject = dict(mag=4e-12, eog=150e-6)
epochs = mne.Epochs(
raw,
events,
event_id,
tmin,
tmax,
picks=picks,
decim=2, # just for speed!
baseline=None,
reject=reject,
preload=True,
)
epochs.drop_channels(["EOG 061"])
epochs.equalize_event_counts(event_id)
# Obtain the data as a 3D matrix and transpose it such that
# the dimensions are as expected for the cluster permutation test:
# n_epochs × n_times × n_channels
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]
Not setting metadata
288 matching events found
No baseline correction applied
Created an SSP operator (subspace dimension = 3)
3 projection items activated
Using data from preloaded Raw for 288 events and 106 original time points (prior to decimation) ...
Rejecting epoch based on EOG : ['EOG 061']
Rejecting epoch based on EOG : ['EOG 061']
Rejecting epoch based on EOG : ['EOG 061']
Rejecting epoch based on EOG : ['EOG 061']
Rejecting epoch based on EOG : ['EOG 061']
Rejecting epoch based on EOG : ['EOG 061']
Rejecting epoch based on EOG : ['EOG 061']
Rejecting epoch based on EOG : ['EOG 061']
Rejecting epoch based on EOG : ['EOG 061']
Rejecting epoch based on EOG : ['EOG 061']
Rejecting epoch based on EOG : ['EOG 061']
Rejecting epoch based on EOG : ['EOG 061']
Rejecting epoch based on EOG : ['EOG 061']
Rejecting epoch based on EOG : ['EOG 061']
Rejecting epoch based on EOG : ['EOG 061']
Rejecting epoch based on MAG : ['MEG 1711']
Rejecting epoch based on EOG : ['EOG 061']
Rejecting epoch based on EOG : ['EOG 061']
Rejecting epoch based on EOG : ['EOG 061']
Rejecting epoch based on EOG : ['EOG 061']
Rejecting epoch based on EOG : ['EOG 061']
Rejecting epoch based on MAG : ['MEG 1711']
Rejecting epoch based on EOG : ['EOG 061']
Rejecting epoch based on EOG : ['EOG 061']
Rejecting epoch based on EOG : ['EOG 061']
Rejecting epoch based on EOG : ['EOG 061']
Rejecting epoch based on EOG : ['EOG 061']
Rejecting epoch based on EOG : ['EOG 061']
Rejecting epoch based on EOG : ['EOG 061']
Rejecting epoch based on EOG : ['EOG 061']
Rejecting epoch based on EOG : ['EOG 061']
Rejecting epoch based on EOG : ['EOG 061']
Rejecting epoch based on EOG : ['EOG 061']
Rejecting epoch based on EOG : ['EOG 061']
Rejecting epoch based on EOG : ['EOG 061']
Rejecting epoch based on EOG : ['EOG 061']
Rejecting epoch based on EOG : ['EOG 061']
Rejecting epoch based on EOG : ['EOG 061']
Rejecting epoch based on EOG : ['EOG 061']
Rejecting epoch based on EOG : ['EOG 061']
Rejecting epoch based on EOG : ['EOG 061']
Rejecting epoch based on EOG : ['EOG 061']
Rejecting epoch based on EOG : ['EOG 061']
Rejecting epoch based on EOG : ['EOG 061']
Rejecting epoch based on EOG : ['EOG 061']
Rejecting epoch based on EOG : ['EOG 061']
Rejecting epoch based on EOG : ['EOG 061']
Rejecting epoch based on EOG : ['EOG 061']
48 bad epochs dropped
Dropped 16 epochs: 50, 51, 84, 93, 95, 96, 129, 146, 149, 154, 158, 195, 201, 203, 211, 212
Find the FieldTrip neighbor definition to setup sensor adjacency#
adjacency, ch_names = find_ch_adjacency(epochs.info, ch_type="mag")
print(type(adjacency)) # it's a sparse matrix!
mne.viz.plot_ch_adjacency(epochs.info, adjacency, ch_names)
Reading adjacency matrix for neuromag306mag.
<class 'scipy.sparse._csr.csr_array'>
Compute permutation statistic#
How does it work? We use clustering to “bind” together features which are similar. Our features are the magnetic fields measured over our sensor array at different times. This reduces the multiple comparison problem. To compute the actual test-statistic, we first sum all F-values in all clusters. We end up with one statistic for each cluster. Then we generate a distribution from the data by shuffling our conditions between our samples and recomputing our clusters and the test statistics. We test for the significance of a given cluster by computing the probability of observing a cluster of that size [1][2].
# We are running an F test, so we look at the upper tail
# see also: https://stats.stackexchange.com/a/73993
tail = 1
# We want to set a critical test statistic (here: F), to determine when
# clusters are being formed. Using Scipy's percent point function of the F
# distribution, we can conveniently select a threshold that corresponds to
# some alpha level that we arbitrarily pick.
alpha_cluster_forming = 0.001
# For an F test we need the degrees of freedom for the numerator
# (number of conditions - 1) and the denominator (number of observations
# - number of conditions):
n_conditions = len(event_id)
n_observations = len(X[0])
dfn = n_conditions - 1
dfd = n_observations - n_conditions
# Note: we calculate 1 - alpha_cluster_forming to get the critical value
# on the right tail
f_thresh = scipy.stats.f.ppf(1 - alpha_cluster_forming, dfn=dfn, dfd=dfd)
# run the cluster based permutation analysis
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
stat_fun(H1): min=0.0033787374718191373 max=207.77535362731717
Running initial clustering …
Found 22 clusters
0%| | Permuting : 0/999 [00:00<?, ?it/s]
1%| | Permuting : 5/999 [00:00<00:06, 147.41it/s]
1%| | Permuting : 11/999 [00:00<00:06, 162.87it/s]
2%|▏ | Permuting : 18/999 [00:00<00:05, 178.53it/s]
3%|▎ | Permuting : 25/999 [00:00<00:05, 186.50it/s]
3%|▎ | Permuting : 32/999 [00:00<00:05, 191.25it/s]
4%|▍ | Permuting : 38/999 [00:00<00:05, 188.82it/s]
4%|▍ | Permuting : 43/999 [00:00<00:05, 182.10it/s]
5%|▍ | Permuting : 49/999 [00:00<00:05, 181.44it/s]
6%|▌ | Permuting : 57/999 [00:00<00:04, 189.01it/s]
6%|▋ | Permuting : 64/999 [00:00<00:04, 191.36it/s]
7%|▋ | Permuting : 69/999 [00:00<00:04, 186.39it/s]
8%|▊ | Permuting : 77/999 [00:00<00:04, 191.95it/s]
9%|▊ | Permuting : 86/999 [00:00<00:04, 199.62it/s]
9%|▉ | Permuting : 94/999 [00:00<00:04, 203.27it/s]
10%|█ | Permuting : 104/999 [00:00<00:04, 211.98it/s]
11%|█▏ | Permuting : 114/999 [00:00<00:04, 219.56it/s]
12%|█▏ | Permuting : 124/999 [00:00<00:03, 226.20it/s]
13%|█▎ | Permuting : 133/999 [00:00<00:03, 229.59it/s]
14%|█▍ | Permuting : 142/999 [00:00<00:03, 232.59it/s]
15%|█▌ | Permuting : 151/999 [00:00<00:03, 235.29it/s]
16%|█▌ | Permuting : 161/999 [00:00<00:03, 239.92it/s]
17%|█▋ | Permuting : 170/999 [00:00<00:03, 241.93it/s]
18%|█▊ | Permuting : 179/999 [00:00<00:03, 243.76it/s]
19%|█▉ | Permuting : 188/999 [00:00<00:03, 245.41it/s]
20%|█▉ | Permuting : 199/999 [00:00<00:03, 251.03it/s]
21%|██ | Permuting : 208/999 [00:00<00:03, 252.12it/s]
22%|██▏ | Permuting : 218/999 [00:00<00:03, 255.12it/s]
23%|██▎ | Permuting : 227/999 [00:00<00:03, 255.91it/s]
24%|██▎ | Permuting : 236/999 [00:00<00:02, 256.62it/s]
25%|██▍ | Permuting : 245/999 [00:01<00:02, 257.15it/s]
25%|██▌ | Permuting : 254/999 [00:01<00:02, 257.80it/s]
26%|██▋ | Permuting : 263/999 [00:01<00:02, 258.38it/s]
27%|██▋ | Permuting : 273/999 [00:01<00:02, 260.73it/s]
28%|██▊ | Permuting : 283/999 [00:01<00:02, 262.93it/s]
29%|██▉ | Permuting : 292/999 [00:01<00:02, 263.17it/s]
30%|███ | Permuting : 302/999 [00:01<00:02, 265.18it/s]
31%|███ | Permuting : 311/999 [00:01<00:02, 265.31it/s]
32%|███▏ | Permuting : 321/999 [00:01<00:02, 267.17it/s]
33%|███▎ | Permuting : 331/999 [00:01<00:02, 268.87it/s]
34%|███▍ | Permuting : 341/999 [00:01<00:02, 270.50it/s]
35%|███▌ | Permuting : 351/999 [00:01<00:02, 272.01it/s]
36%|███▌ | Permuting : 361/999 [00:01<00:02, 273.43it/s]
37%|███▋ | Permuting : 370/999 [00:01<00:02, 273.09it/s]
38%|███▊ | Permuting : 378/999 [00:01<00:02, 271.10it/s]
39%|███▉ | Permuting : 388/999 [00:01<00:02, 272.53it/s]
40%|███▉ | Permuting : 399/999 [00:01<00:02, 275.51it/s]
41%|████ | Permuting : 408/999 [00:01<00:02, 275.05it/s]
42%|████▏ | Permuting : 418/999 [00:01<00:02, 276.24it/s]
43%|████▎ | Permuting : 427/999 [00:01<00:02, 275.75it/s]
44%|████▍ | Permuting : 438/999 [00:01<00:02, 278.34it/s]
45%|████▍ | Permuting : 448/999 [00:01<00:01, 279.35it/s]
46%|████▌ | Permuting : 459/999 [00:01<00:01, 281.90it/s]
47%|████▋ | Permuting : 468/999 [00:01<00:01, 281.13it/s]
48%|████▊ | Permuting : 477/999 [00:01<00:01, 280.38it/s]
49%|████▊ | Permuting : 487/999 [00:01<00:01, 281.27it/s]
50%|████▉ | Permuting : 497/999 [00:01<00:01, 282.10it/s]
51%|█████ | Permuting : 507/999 [00:01<00:01, 282.88it/s]
52%|█████▏ | Permuting : 516/999 [00:01<00:01, 282.06it/s]
53%|█████▎ | Permuting : 525/999 [00:01<00:01, 281.32it/s]
54%|█████▎ | Permuting : 536/999 [00:02<00:01, 283.66it/s]
55%|█████▍ | Permuting : 547/999 [00:02<00:01, 285.91it/s]
56%|█████▌ | Permuting : 557/999 [00:02<00:01, 286.54it/s]
57%|█████▋ | Permuting : 567/999 [00:02<00:01, 287.16it/s]
58%|█████▊ | Permuting : 577/999 [00:02<00:01, 287.75it/s]
59%|█████▉ | Permuting : 587/999 [00:02<00:01, 288.30it/s]
60%|█████▉ | Permuting : 597/999 [00:02<00:01, 288.82it/s]
61%|██████ | Permuting : 608/999 [00:02<00:01, 290.86it/s]
62%|██████▏ | Permuting : 617/999 [00:02<00:01, 289.72it/s]
63%|██████▎ | Permuting : 626/999 [00:02<00:01, 288.64it/s]
64%|██████▎ | Permuting : 636/999 [00:02<00:01, 289.15it/s]
65%|██████▍ | Permuting : 645/999 [00:02<00:01, 288.09it/s]
66%|██████▌ | Permuting : 655/999 [00:02<00:01, 288.62it/s]
67%|██████▋ | Permuting : 665/999 [00:02<00:01, 289.13it/s]
68%|██████▊ | Permuting : 675/999 [00:02<00:01, 289.61it/s]
69%|██████▊ | Permuting : 685/999 [00:02<00:01, 289.95it/s]
70%|██████▉ | Permuting : 695/999 [00:02<00:01, 290.39it/s]
71%|███████ | Permuting : 707/999 [00:02<00:00, 293.84it/s]
72%|███████▏ | Permuting : 717/999 [00:02<00:00, 294.09it/s]
73%|███████▎ | Permuting : 727/999 [00:02<00:00, 294.31it/s]
74%|███████▎ | Permuting : 736/999 [00:02<00:00, 293.00it/s]
75%|███████▍ | Permuting : 746/999 [00:02<00:00, 293.27it/s]
76%|███████▌ | Permuting : 755/999 [00:02<00:00, 292.03it/s]
77%|███████▋ | Permuting : 765/999 [00:02<00:00, 292.36it/s]
78%|███████▊ | Permuting : 775/999 [00:02<00:00, 292.67it/s]
79%|███████▊ | Permuting : 785/999 [00:02<00:00, 292.96it/s]
79%|███████▉ | Permuting : 794/999 [00:02<00:00, 291.73it/s]
80%|████████ | Permuting : 804/999 [00:02<00:00, 292.07it/s]
81%|████████▏ | Permuting : 813/999 [00:02<00:00, 290.89it/s]
82%|████████▏ | Permuting : 823/999 [00:02<00:00, 291.26it/s]
83%|████████▎ | Permuting : 834/999 [00:03<00:00, 293.14it/s]
84%|████████▍ | Permuting : 843/999 [00:03<00:00, 291.90it/s]
85%|████████▌ | Permuting : 853/999 [00:03<00:00, 292.23it/s]
86%|████████▋ | Permuting : 862/999 [00:03<00:00, 291.04it/s]
87%|████████▋ | Permuting : 872/999 [00:03<00:00, 291.41it/s]
88%|████████▊ | Permuting : 882/999 [00:03<00:00, 291.76it/s]
89%|████████▉ | Permuting : 892/999 [00:03<00:00, 292.10it/s]
90%|█████████ | Permuting : 902/999 [00:03<00:00, 292.42it/s]
91%|█████████▏| Permuting : 912/999 [00:03<00:00, 292.71it/s]
92%|█████████▏| Permuting : 921/999 [00:03<00:00, 291.50it/s]
93%|█████████▎| Permuting : 930/999 [00:03<00:00, 290.35it/s]
94%|█████████▍| Permuting : 940/999 [00:03<00:00, 290.76it/s]
95%|█████████▌| Permuting : 950/999 [00:03<00:00, 291.15it/s]
96%|█████████▌| Permuting : 959/999 [00:03<00:00, 290.02it/s]
97%|█████████▋| Permuting : 969/999 [00:03<00:00, 290.44it/s]
98%|█████████▊| Permuting : 979/999 [00:03<00:00, 290.85it/s]
99%|█████████▉| Permuting : 989/999 [00:03<00:00, 291.23it/s]
100%|█████████▉| Permuting : 998/999 [00:03<00:00, 290.09it/s]
100%|██████████| Permuting : 999/999 [00:03<00:00, 277.61it/s]
Note
Note how we only specified an adjacency for sensors! However,
because we used mne.stats.spatio_temporal_cluster_test(),
an adjacency for time points was automatically taken into
account. That is, at time point N, the time points N - 1 and
N + 1 were considered as adjacent (this is also called “lattice
adjacency”). This is only possible because we ran the analysis on
2D data (times × channels) per observation … for 3D data per
observation (e.g., times × frequencies × channels), we will need
to use mne.stats.combine_adjacency(), as shown further
below.
Note also that the same functions work with source estimates. The only differences are the origin of the data, the size, and the adjacency definition. It can be used for single trials or for groups of subjects.
Visualize clusters#
# We subselect clusters that we consider significant at an arbitrarily
# picked alpha level: "p_accept".
# NOTE: remember the caveats with respect to "significant" clusters that
# we mentioned in the introduction of this tutorial!
p_accept = 0.01
good_cluster_inds = np.where(p_values < p_accept)[0]
# configure variables for visualization
colors = {"Aud": "crimson", "Vis": "steelblue"}
linestyles = {"L": "-", "R": "--"}
# organize data for plotting
evokeds = {cond: epochs[cond].average() for cond in event_id}
# loop over clusters
for i_clu, clu_idx in enumerate(good_cluster_inds):
# unpack cluster information, get unique indices
time_inds, space_inds = np.squeeze(clusters[clu_idx])
ch_inds = np.unique(space_inds)
time_inds = np.unique(time_inds)
# get topography for F stat
f_map = F_obs[time_inds, ...].mean(axis=0)
# get signals at the sensors contributing to the cluster
sig_times = epochs.times[time_inds]
# create spatial mask
mask = np.zeros((f_map.shape[0], 1), dtype=bool)
mask[ch_inds, :] = True
# initialize figure
fig, ax_topo = plt.subplots(1, 1, figsize=(10, 3), layout="constrained")
# plot average test statistic and mark significant sensors
f_evoked = mne.EvokedArray(f_map[:, np.newaxis], epochs.info, tmin=0)
f_evoked.plot_topomap(
times=0,
mask=mask,
axes=ax_topo,
cmap="Reds",
vlim=(np.min, np.max),
show=False,
colorbar=False,
mask_params=dict(markersize=10),
)
image = ax_topo.images[0]
# remove the title that would otherwise say "0.000 s"
ax_topo.set_title("")
# create additional axes (for ERF and colorbar)
divider = make_axes_locatable(ax_topo)
# add axes for colorbar
ax_colorbar = divider.append_axes("right", size="5%", pad=0.05)
plt.colorbar(image, cax=ax_colorbar)
ax_topo.set_xlabel(
"Averaged F-map ({:0.3f} - {:0.3f} s)".format(*sig_times[[0, -1]])
)
# add new axis for time courses and plot time courses
ax_signals = divider.append_axes("right", size="300%", pad=1.2)
title = f"Cluster #{i_clu + 1}, {len(ch_inds)} sensor"
if len(ch_inds) > 1:
title += "s (mean)"
plot_compare_evokeds(
evokeds,
title=title,
picks=ch_inds,
axes=ax_signals,
colors=colors,
linestyles=linestyles,
show=False,
split_legend=True,
truncate_yaxis="auto",
)
# plot temporal cluster extent
ymin, ymax = ax_signals.get_ylim()
ax_signals.fill_betweenx(
(ymin, ymax), sig_times[0], sig_times[-1], color="orange", alpha=0.3
)
plt.show()
combining channels using RMS (mag channels)
combining channels using RMS (mag channels)
combining channels using RMS (mag channels)
combining channels using RMS (mag channels)
combining channels using RMS (mag channels)
combining channels using RMS (mag channels)
combining channels using RMS (mag channels)
combining channels using RMS (mag channels)
combining channels using RMS (mag channels)
combining channels using RMS (mag channels)
combining channels using RMS (mag channels)
combining channels using RMS (mag channels)
combining channels using RMS (mag channels)
combining channels using RMS (mag channels)
combining channels using RMS (mag channels)
combining channels using RMS (mag channels)
combining channels using RMS (mag channels)
combining channels using RMS (mag channels)
combining channels using RMS (mag channels)
combining channels using RMS (mag channels)
combining channels using RMS (mag channels)
combining channels using RMS (mag channels)
combining channels using RMS (mag channels)
combining channels using RMS (mag channels)
combining channels using RMS (mag channels)
combining channels using RMS (mag channels)
combining channels using RMS (mag channels)
combining channels using RMS (mag channels)
Permutation statistic for time-frequencies#
Let’s do the same thing with the time-frequency decomposition of the data (see Frequency and time-frequency sensor analysis for a tutorial and Time-frequency on simulated data (Multitaper vs. Morlet vs. Stockwell vs. Hilbert) for a comparison of time-frequency methods) to show how cluster permutations can be done on higher-dimensional data.
decim = 4
freqs = np.arange(7, 30, 3) # define frequencies of interest
n_cycles = freqs / freqs[0]
epochs_power = list()
for condition in [epochs[k] for k in ("Aud/L", "Vis/L")]:
this_tfr = condition.compute_tfr(
method="morlet",
freqs=freqs,
n_cycles=n_cycles,
decim=decim,
average=False,
return_itc=False,
)
this_tfr.apply_baseline(mode="ratio", baseline=(None, 0))
epochs_power.append(this_tfr.data)
# transpose again to (epochs, frequencies, times, channels)
X = [np.transpose(x, (0, 2, 3, 1)) for x in epochs_power]
[Parallel(n_jobs=1)]: Done 17 tasks | elapsed: 0.1s
[Parallel(n_jobs=1)]: Done 71 tasks | elapsed: 0.4s
Applying baseline correction (mode: ratio)
[Parallel(n_jobs=1)]: Done 17 tasks | elapsed: 0.1s
[Parallel(n_jobs=1)]: Done 71 tasks | elapsed: 0.4s
Applying baseline correction (mode: ratio)
Remember the note on the adjacency matrix from above: For 3D data, as here,
we must use mne.stats.combine_adjacency() to extend the
sensor-based adjacency to incorporate the time-frequency plane as well.
Here, the integer inputs are converted into a lattice and combined with the sensor adjacency matrix so that data at similar times and with similar frequencies and at close sensor locations are clustered together.
# our data at each observation is of shape frequencies × times × channels
tfr_adjacency = combine_adjacency(len(freqs), len(this_tfr.times), adjacency)
Now we can run the cluster permutation test, but first we have to set a threshold. This example decimates in time and uses few frequencies so we need to increase the threshold from the default value in order to have differentiated clusters (i.e., so that our algorithm doesn’t just find one large cluster). For a more principled method of setting this parameter, threshold-free cluster enhancement may be used. See Statistical inference for a discussion.
# This time we don't calculate a threshold based on the F distribution.
# We might as well select an arbitrary threshold for cluster forming
tfr_threshold = 15.0
# run cluster based permutation analysis
cluster_stats = spatio_temporal_cluster_test(
X,
n_permutations=1000,
threshold=tfr_threshold,
tail=1,
n_jobs=None,
buffer_size=None,
adjacency=tfr_adjacency,
)
stat_fun(H1): min=6.752417666025854e-07 max=38.408365559505015
Running initial clustering …
Found 3 clusters
0%| | Permuting : 0/999 [00:00<?, ?it/s]
0%| | Permuting : 3/999 [00:00<01:17, 12.89it/s]
1%| | Permuting : 8/999 [00:00<00:32, 30.19it/s]
1%| | Permuting : 12/999 [00:00<00:24, 40.05it/s]
2%|▏ | Permuting : 19/999 [00:00<00:16, 57.82it/s]
2%|▏ | Permuting : 23/999 [00:00<00:15, 64.10it/s]
3%|▎ | Permuting : 26/999 [00:00<00:14, 65.60it/s]
3%|▎ | Permuting : 29/999 [00:00<00:14, 67.77it/s]
4%|▎ | Permuting : 36/999 [00:00<00:12, 79.32it/s]
4%|▍ | Permuting : 40/999 [00:00<00:11, 82.67it/s]
4%|▍ | Permuting : 42/999 [00:00<00:11, 80.56it/s]
4%|▍ | Permuting : 44/999 [00:00<00:12, 77.99it/s]
5%|▍ | Permuting : 48/999 [00:00<00:11, 81.10it/s]
5%|▌ | Permuting : 53/999 [00:00<00:10, 86.11it/s]
6%|▌ | Permuting : 56/999 [00:00<00:14, 65.17it/s]
6%|▌ | Permuting : 60/999 [00:00<00:13, 68.09it/s]
6%|▋ | Permuting : 63/999 [00:00<00:13, 68.68it/s]
7%|▋ | Permuting : 69/999 [00:00<00:12, 73.96it/s]
7%|▋ | Permuting : 71/999 [00:01<00:12, 72.62it/s]
8%|▊ | Permuting : 75/999 [00:01<00:12, 75.06it/s]
8%|▊ | Permuting : 77/999 [00:01<00:12, 73.87it/s]
8%|▊ | Permuting : 80/999 [00:01<00:12, 74.07it/s]
8%|▊ | Permuting : 83/999 [00:01<00:12, 74.86it/s]
9%|▊ | Permuting : 85/999 [00:01<00:12, 73.77it/s]
9%|▊ | Permuting : 87/999 [00:01<00:12, 72.50it/s]
10%|▉ | Permuting : 96/999 [00:01<00:11, 81.66it/s]
10%|▉ | Permuting : 99/999 [00:01<00:13, 68.10it/s]
11%|█ | Permuting : 105/999 [00:01<00:12, 72.74it/s]
11%|█ | Permuting : 109/999 [00:01<00:11, 74.23it/s]
11%|█▏ | Permuting : 114/999 [00:01<00:11, 76.90it/s]
12%|█▏ | Permuting : 118/999 [00:01<00:11, 78.69it/s]
12%|█▏ | Permuting : 123/999 [00:01<00:10, 81.69it/s]
13%|█▎ | Permuting : 128/999 [00:01<00:10, 84.59it/s]
13%|█▎ | Permuting : 133/999 [00:01<00:09, 87.39it/s]
14%|█▍ | Permuting : 139/999 [00:01<00:09, 91.37it/s]
14%|█▍ | Permuting : 144/999 [00:01<00:09, 93.90it/s]
15%|█▍ | Permuting : 149/999 [00:01<00:08, 96.34it/s]
15%|█▌ | Permuting : 154/999 [00:01<00:08, 98.68it/s]
16%|█▌ | Permuting : 158/999 [00:01<00:08, 98.94it/s]
16%|█▌ | Permuting : 161/999 [00:01<00:08, 98.52it/s]
17%|█▋ | Permuting : 165/999 [00:01<00:08, 98.64it/s]
17%|█▋ | Permuting : 171/999 [00:01<00:08, 101.54it/s]
17%|█▋ | Permuting : 172/999 [00:02<00:09, 84.20it/s]
18%|█▊ | Permuting : 178/999 [00:02<00:09, 87.86it/s]
18%|█▊ | Permuting : 180/999 [00:02<00:09, 86.23it/s]
19%|█▊ | Permuting : 186/999 [00:02<00:09, 89.88it/s]
19%|█▉ | Permuting : 190/999 [00:02<00:08, 91.05it/s]
20%|██ | Permuting : 202/999 [00:02<00:07, 101.77it/s]
20%|██ | Permuting : 204/999 [00:02<00:07, 99.45it/s]
21%|██ | Permuting : 208/999 [00:02<00:07, 100.26it/s]
21%|██▏ | Permuting : 213/999 [00:02<00:07, 102.27it/s]
22%|██▏ | Permuting : 217/999 [00:02<00:07, 102.26it/s]
22%|██▏ | Permuting : 220/999 [00:02<00:07, 101.73it/s]
22%|██▏ | Permuting : 224/999 [00:02<00:07, 102.46it/s]
23%|██▎ | Permuting : 226/999 [00:02<00:07, 100.34it/s]
23%|██▎ | Permuting : 227/999 [00:02<00:07, 96.72it/s]
23%|██▎ | Permuting : 231/999 [00:02<00:09, 83.70it/s]
23%|██▎ | Permuting : 234/999 [00:02<00:09, 83.47it/s]
24%|██▍ | Permuting : 238/999 [00:02<00:08, 84.80it/s]
24%|██▍ | Permuting : 243/999 [00:02<00:08, 87.23it/s]
25%|██▍ | Permuting : 245/999 [00:02<00:08, 85.96it/s]
25%|██▍ | Permuting : 249/999 [00:02<00:08, 87.24it/s]
26%|██▌ | Permuting : 255/999 [00:02<00:08, 90.81it/s]
26%|██▌ | Permuting : 260/999 [00:02<00:07, 93.11it/s]
27%|██▋ | Permuting : 268/999 [00:03<00:07, 98.32it/s]
28%|██▊ | Permuting : 275/999 [00:03<00:07, 102.11it/s]
28%|██▊ | Permuting : 280/999 [00:03<00:06, 104.00it/s]
28%|██▊ | Permuting : 284/999 [00:03<00:06, 104.62it/s]
29%|██▉ | Permuting : 292/999 [00:03<00:06, 109.43it/s]
30%|██▉ | Permuting : 296/999 [00:03<00:06, 109.84it/s]
30%|██▉ | Permuting : 297/999 [00:03<00:07, 91.20it/s]
30%|███ | Permuting : 300/999 [00:03<00:07, 90.63it/s]
30%|███ | Permuting : 304/999 [00:03<00:07, 91.66it/s]
31%|███ | Permuting : 307/999 [00:03<00:07, 91.04it/s]
31%|███ | Permuting : 310/999 [00:03<00:07, 90.97it/s]
31%|███ | Permuting : 312/999 [00:03<00:07, 89.27it/s]
31%|███▏ | Permuting : 314/999 [00:03<00:07, 87.62it/s]
32%|███▏ | Permuting : 317/999 [00:03<00:07, 87.68it/s]
32%|███▏ | Permuting : 321/999 [00:03<00:07, 88.86it/s]
32%|███▏ | Permuting : 324/999 [00:03<00:07, 88.35it/s]
33%|███▎ | Permuting : 328/999 [00:03<00:07, 89.22it/s]
33%|███▎ | Permuting : 332/999 [00:03<00:07, 90.38it/s]
34%|███▎ | Permuting : 337/999 [00:03<00:07, 92.71it/s]
34%|███▍ | Permuting : 342/999 [00:03<00:06, 94.97it/s]
35%|███▌ | Permuting : 353/999 [00:03<00:06, 104.42it/s]
36%|███▌ | Permuting : 355/999 [00:04<00:07, 87.44it/s]
36%|███▌ | Permuting : 358/999 [00:04<00:07, 87.04it/s]
36%|███▌ | Permuting : 360/999 [00:04<00:07, 85.58it/s]
37%|███▋ | Permuting : 369/999 [00:04<00:06, 91.68it/s]
37%|███▋ | Permuting : 371/999 [00:04<00:06, 90.50it/s]
38%|███▊ | Permuting : 376/999 [00:04<00:06, 92.66it/s]
38%|███▊ | Permuting : 378/999 [00:04<00:06, 90.84it/s]
38%|███▊ | Permuting : 382/999 [00:04<00:06, 91.37it/s]
39%|███▉ | Permuting : 390/999 [00:04<00:06, 96.95it/s]
40%|███▉ | Permuting : 397/999 [00:04<00:05, 101.01it/s]
40%|████ | Permuting : 401/999 [00:04<00:05, 101.69it/s]
41%|████ | Permuting : 407/999 [00:04<00:05, 104.70it/s]
41%|████ | Permuting : 412/999 [00:04<00:05, 106.43it/s]
41%|████▏ | Permuting : 414/999 [00:04<00:05, 103.91it/s]
42%|████▏ | Permuting : 421/999 [00:04<00:05, 107.46it/s]
43%|████▎ | Permuting : 426/999 [00:04<00:06, 91.39it/s]
43%|████▎ | Permuting : 432/999 [00:04<00:06, 93.90it/s]
44%|████▍ | Permuting : 441/999 [00:04<00:05, 99.95it/s]
45%|████▍ | Permuting : 447/999 [00:04<00:05, 102.71it/s]
45%|████▍ | Permuting : 449/999 [00:04<00:05, 100.60it/s]
46%|████▌ | Permuting : 455/999 [00:04<00:05, 102.85it/s]
46%|████▌ | Permuting : 459/999 [00:04<00:05, 103.44it/s]
46%|████▌ | Permuting : 462/999 [00:05<00:05, 102.34it/s]
47%|████▋ | Permuting : 466/999 [00:05<00:05, 102.95it/s]
47%|████▋ | Permuting : 471/999 [00:05<00:05, 104.65it/s]
47%|████▋ | Permuting : 474/999 [00:05<00:05, 103.46it/s]
48%|████▊ | Permuting : 478/999 [00:05<00:05, 104.05it/s]
48%|████▊ | Permuting : 483/999 [00:05<00:04, 105.77it/s]
49%|████▊ | Permuting : 486/999 [00:05<00:05, 89.72it/s]
49%|████▉ | Permuting : 488/999 [00:05<00:05, 88.60it/s]
49%|████▉ | Permuting : 491/999 [00:05<00:05, 88.13it/s]
50%|████▉ | Permuting : 495/999 [00:05<00:05, 89.03it/s]
50%|█████ | Permuting : 501/999 [00:05<00:05, 91.70it/s]
51%|█████ | Permuting : 506/999 [00:05<00:05, 93.74it/s]
51%|█████ | Permuting : 511/999 [00:05<00:05, 95.72it/s]
52%|█████▏ | Permuting : 516/999 [00:05<00:04, 97.67it/s]
52%|█████▏ | Permuting : 520/999 [00:05<00:04, 97.91it/s]
52%|█████▏ | Permuting : 523/999 [00:05<00:04, 97.58it/s]
53%|█████▎ | Permuting : 527/999 [00:05<00:04, 98.37it/s]
53%|█████▎ | Permuting : 533/999 [00:05<00:04, 101.45it/s]
54%|█████▍ | Permuting : 537/999 [00:05<00:04, 102.10it/s]
54%|█████▍ | Permuting : 538/999 [00:05<00:05, 85.78it/s]
54%|█████▍ | Permuting : 544/999 [00:05<00:05, 88.99it/s]
55%|█████▍ | Permuting : 546/999 [00:06<00:05, 87.48it/s]
55%|█████▌ | Permuting : 550/999 [00:06<00:05, 88.60it/s]
56%|█████▌ | Permuting : 555/999 [00:06<00:04, 90.75it/s]
56%|█████▌ | Permuting : 561/999 [00:06<00:04, 93.95it/s]
57%|█████▋ | Permuting : 566/999 [00:06<00:04, 95.97it/s]
57%|█████▋ | Permuting : 568/999 [00:06<00:04, 94.07it/s]
58%|█████▊ | Permuting : 575/999 [00:06<00:04, 97.82it/s]
58%|█████▊ | Permuting : 580/999 [00:06<00:04, 99.75it/s]
59%|█████▊ | Permuting : 585/999 [00:06<00:04, 101.61it/s]
59%|█████▉ | Permuting : 587/999 [00:06<00:04, 99.38it/s]
60%|█████▉ | Permuting : 599/999 [00:06<00:03, 109.45it/s]
60%|██████ | Permuting : 601/999 [00:06<00:03, 106.83it/s]
61%|██████ | Permuting : 605/999 [00:06<00:03, 107.31it/s]
61%|██████ | Permuting : 610/999 [00:06<00:03, 108.97it/s]
62%|██████▏ | Permuting : 617/999 [00:06<00:03, 112.87it/s]
62%|██████▏ | Permuting : 623/999 [00:06<00:03, 97.35it/s]
63%|██████▎ | Permuting : 625/999 [00:06<00:03, 95.52it/s]
63%|██████▎ | Permuting : 630/999 [00:06<00:03, 97.39it/s]
63%|██████▎ | Permuting : 632/999 [00:06<00:03, 95.49it/s]
64%|██████▎ | Permuting : 635/999 [00:06<00:03, 94.75it/s]
64%|██████▍ | Permuting : 639/999 [00:06<00:03, 95.39it/s]
64%|██████▍ | Permuting : 642/999 [00:06<00:03, 94.62it/s]
65%|██████▍ | Permuting : 646/999 [00:06<00:03, 94.98it/s]
66%|██████▌ | Permuting : 657/999 [00:06<00:03, 103.68it/s]
66%|██████▌ | Permuting : 659/999 [00:07<00:03, 101.54it/s]
66%|██████▋ | Permuting : 664/999 [00:07<00:03, 103.32it/s]
67%|██████▋ | Permuting : 668/999 [00:07<00:03, 103.90it/s]
68%|██████▊ | Permuting : 678/999 [00:07<00:02, 110.63it/s]
68%|██████▊ | Permuting : 683/999 [00:07<00:02, 111.46it/s]
69%|██████▊ | Permuting : 686/999 [00:07<00:03, 95.62it/s]
69%|██████▉ | Permuting : 692/999 [00:07<00:03, 98.49it/s]
70%|██████▉ | Permuting : 695/999 [00:07<00:03, 97.62it/s]
70%|███████ | Permuting : 703/999 [00:07<00:02, 102.61it/s]
71%|███████ | Permuting : 707/999 [00:07<00:02, 102.53it/s]
71%|███████ | Permuting : 709/999 [00:07<00:02, 100.96it/s]
71%|███████▏ | Permuting : 712/999 [00:07<00:02, 99.97it/s]
72%|███████▏ | Permuting : 716/999 [00:07<00:02, 100.68it/s]
72%|███████▏ | Permuting : 721/999 [00:07<00:02, 102.47it/s]
73%|███████▎ | Permuting : 726/999 [00:07<00:02, 104.22it/s]
73%|███████▎ | Permuting : 731/999 [00:07<00:02, 105.92it/s]
73%|███████▎ | Permuting : 734/999 [00:07<00:02, 104.65it/s]
74%|███████▍ | Permuting : 738/999 [00:07<00:02, 105.21it/s]
74%|███████▍ | Permuting : 743/999 [00:07<00:02, 106.91it/s]
75%|███████▍ | Permuting : 745/999 [00:07<00:02, 89.77it/s]
75%|███████▍ | Permuting : 749/999 [00:08<00:02, 90.31it/s]
75%|███████▌ | Permuting : 752/999 [00:08<00:02, 90.27it/s]
75%|███████▌ | Permuting : 754/999 [00:08<00:02, 88.92it/s]
76%|███████▌ | Permuting : 758/999 [00:08<00:02, 89.98it/s]
76%|███████▋ | Permuting : 763/999 [00:08<00:02, 92.12it/s]
77%|███████▋ | Permuting : 770/999 [00:08<00:02, 96.29it/s]
77%|███████▋ | Permuting : 774/999 [00:08<00:02, 97.11it/s]
78%|███████▊ | Permuting : 782/999 [00:08<00:02, 101.86it/s]
79%|███████▉ | Permuting : 790/999 [00:08<00:01, 107.06it/s]
79%|███████▉ | Permuting : 792/999 [00:08<00:01, 104.79it/s]
80%|███████▉ | Permuting : 795/999 [00:08<00:01, 103.55it/s]
80%|███████▉ | Permuting : 799/999 [00:08<00:01, 103.85it/s]
80%|████████ | Permuting : 800/999 [00:08<00:02, 86.67it/s]
81%|████████ | Permuting : 806/999 [00:08<00:02, 89.83it/s]
81%|████████ | Permuting : 811/999 [00:08<00:02, 91.89it/s]
82%|████████▏ | Permuting : 816/999 [00:08<00:01, 93.89it/s]
82%|████████▏ | Permuting : 821/999 [00:08<00:01, 95.86it/s]
83%|████████▎ | Permuting : 825/999 [00:08<00:01, 96.06it/s]
83%|████████▎ | Permuting : 830/999 [00:08<00:01, 97.50it/s]
84%|████████▎ | Permuting : 835/999 [00:08<00:01, 99.39it/s]
84%|████████▍ | Permuting : 840/999 [00:08<00:01, 101.21it/s]
84%|████████▍ | Permuting : 843/999 [00:08<00:01, 100.18it/s]
85%|████████▍ | Permuting : 846/999 [00:08<00:01, 99.13it/s]
85%|████████▌ | Permuting : 850/999 [00:08<00:01, 99.69it/s]
85%|████████▌ | Permuting : 851/999 [00:09<00:01, 96.39it/s]
86%|████████▌ | Permuting : 856/999 [00:09<00:01, 98.43it/s]
86%|████████▌ | Permuting : 861/999 [00:09<00:01, 100.39it/s]
87%|████████▋ | Permuting : 867/999 [00:09<00:01, 103.50it/s]
87%|████████▋ | Permuting : 869/999 [00:09<00:01, 88.51it/s]
87%|████████▋ | Permuting : 872/999 [00:09<00:01, 88.08it/s]
88%|████████▊ | Permuting : 875/999 [00:09<00:01, 87.63it/s]
88%|████████▊ | Permuting : 884/999 [00:09<00:01, 94.16it/s]
89%|████████▉ | Permuting : 892/999 [00:09<00:01, 98.98it/s]
90%|████████▉ | Permuting : 897/999 [00:09<00:01, 100.81it/s]
91%|█████████ | Permuting : 907/999 [00:09<00:00, 108.18it/s]
91%|█████████▏| Permuting : 913/999 [00:09<00:00, 110.84it/s]
92%|█████████▏| Permuting : 918/999 [00:09<00:00, 112.28it/s]
93%|█████████▎| Permuting : 929/999 [00:09<00:00, 120.64it/s]
94%|█████████▍| Permuting : 937/999 [00:09<00:00, 124.50it/s]
95%|█████████▍| Permuting : 945/999 [00:09<00:00, 128.23it/s]
95%|█████████▍| Permuting : 949/999 [00:09<00:00, 127.53it/s]
95%|█████████▌| Permuting : 953/999 [00:09<00:00, 127.17it/s]
96%|█████████▌| Permuting : 959/999 [00:09<00:00, 129.23it/s]
96%|█████████▌| Permuting : 961/999 [00:09<00:00, 125.57it/s]
97%|█████████▋| Permuting : 966/999 [00:09<00:00, 126.49it/s]
97%|█████████▋| Permuting : 968/999 [00:09<00:00, 122.91it/s]
98%|█████████▊| Permuting : 977/999 [00:09<00:00, 128.91it/s]
98%|█████████▊| Permuting : 978/999 [00:10<00:00, 108.09it/s]
98%|█████████▊| Permuting : 980/999 [00:10<00:00, 106.11it/s]
99%|█████████▊| Permuting : 985/999 [00:10<00:00, 107.64it/s]
99%|█████████▉| Permuting : 989/999 [00:10<00:00, 108.05it/s]
99%|█████████▉| Permuting : 994/999 [00:10<00:00, 109.59it/s]
100%|██████████| Permuting : 999/999 [00:10<00:00, 111.26it/s]
100%|██████████| Permuting : 999/999 [00:10<00:00, 97.44it/s]
Finally, we can plot our results. It is difficult to visualize clusters in time-frequency-sensor space; plotting time-frequency spectrograms and plotting topomaps display time-frequency and sensor space respectively but they are difficult to combine. We will plot topomaps with the clustered sensors colored in white adjacent to spectrograms in order to provide a visualization of the results. This is a dimensionally limited view, however. Each sensor has its own significant time-frequencies, but, in order to display a single spectrogram, all the time-frequencies that are significant for any sensor in the cluster are plotted as significant. This is a difficulty inherent to visualizing high-dimensional data and should be taken into consideration when interpreting results.
F_obs, clusters, p_values, _ = cluster_stats
good_cluster_inds = np.where(p_values < p_accept)[0]
for i_clu, clu_idx in enumerate(good_cluster_inds):
# unpack cluster information, get unique indices
freq_inds, time_inds, space_inds = clusters[clu_idx]
ch_inds = np.unique(space_inds)
time_inds = np.unique(time_inds)
freq_inds = np.unique(freq_inds)
# get topography for F stat
f_map = F_obs[freq_inds].mean(axis=0)
f_map = f_map[time_inds].mean(axis=0)
# get signals at the sensors contributing to the cluster
sig_times = epochs.times[time_inds]
# initialize figure
fig, ax_topo = plt.subplots(1, 1, figsize=(10, 3), layout="constrained")
# create spatial mask
mask = np.zeros((f_map.shape[0], 1), dtype=bool)
mask[ch_inds, :] = True
# plot average test statistic and mark significant sensors
f_evoked = mne.EvokedArray(f_map[:, np.newaxis], epochs.info, tmin=0)
f_evoked.plot_topomap(
times=0,
mask=mask,
axes=ax_topo,
cmap="Reds",
vlim=(np.min, np.max),
show=False,
colorbar=False,
mask_params=dict(markersize=10),
)
image = ax_topo.images[0]
# create additional axes (for ERF and colorbar)
divider = make_axes_locatable(ax_topo)
# add axes for colorbar
ax_colorbar = divider.append_axes("right", size="5%", pad=0.05)
plt.colorbar(image, cax=ax_colorbar)
ax_topo.set_xlabel(
"Averaged F-map ({:0.3f} - {:0.3f} s)".format(*sig_times[[0, -1]])
)
# remove the title that would otherwise say "0.000 s"
ax_topo.set_title("")
# add new axis for spectrogram
ax_spec = divider.append_axes("right", size="300%", pad=1.2)
title = f"Cluster #{i_clu + 1}, {len(ch_inds)} spectrogram"
if len(ch_inds) > 1:
title += " (max over channels)"
F_obs_plot = F_obs[..., ch_inds].max(axis=-1)
F_obs_plot_sig = np.zeros(F_obs_plot.shape) * np.nan
F_obs_plot_sig[tuple(np.meshgrid(freq_inds, time_inds))] = F_obs_plot[
tuple(np.meshgrid(freq_inds, time_inds))
]
for f_image, cmap in zip([F_obs_plot, F_obs_plot_sig], ["gray", "autumn"]):
c = ax_spec.imshow(
f_image,
cmap=cmap,
aspect="auto",
origin="lower",
extent=[epochs.times[0], epochs.times[-1], freqs[0], freqs[-1]],
)
ax_spec.set_xlabel("Time (ms)")
ax_spec.set_ylabel("Frequency (Hz)")
ax_spec.set_title(title)
# add another colorbar
ax_colorbar2 = divider.append_axes("right", size="5%", pad=0.05)
plt.colorbar(c, cax=ax_colorbar2)
ax_colorbar2.set_ylabel("F-stat")
# clean up viz
plt.show()
Exercises#
References#
Total running time of the script: (0 minutes 21.932 seconds)