Time and Frequency Permutation Cluster Statistics¶

Below is a code sample for performing a permutation test with cluster correction

from ieeg.navigate import trial_ieeg, outliers_to_nan
from ieeg.calc.stats import time_perm_cluster
from ieeg.timefreq.utils import wavelet_scaleogram, crop_pad
import matplotlib.pyplot as plt
import mne

Load and Rereference Data¶

misc_path = mne.datasets.misc.data_path()
raw = mne.io.read_raw(misc_path / 'seeg' / 'sample_seeg_ieeg.fif')

# Exclude bad channels
raw.load_data()

# CAR (common average reference)
raw.set_eeg_reference(ref_channels="average", ch_type='seeg')

Opening raw data file /home/docs/mne_data/MNE-misc-data/seeg/sample_seeg_ieeg.fif...
    Range : 1310640 ... 1370605 =   1311.411 ...  1371.411 secs
Ready.
Reading 0 ... 59965  =      0.000 ...    60.000 secs...
Applying average reference.
Applying a custom ('sEEG',) reference.

	General
	Filename(s)	sample_seeg_ieeg.fif
	MNE object type	Raw
	Measurement date	2019-10-18 at 11:09:44 UTC
	Participant	sub-1
	Experimenter	Unknown
	Acquisition
	Duration	00:01:01 (HH:MM:SS)
	Sampling frequency	999.41 Hz
	Time points	59,966
	Channels
	sEEG
	Head & sensor digitization	122 points
	Filters
	Highpass	0.00 Hz
	Lowpass	499.71 Hz

Calculate Spectra¶

out = []
for epoch, t in zip(('Fixation', 'Response'),  # epochs to extract
                    ((-0.3, 0), (-0.1, 0.2))):  # times to extract
    times = [None, None]
    times[0] = t[0] - 0.5  # add 0.5s to the beginning
    times[1] = t[1] + 0.5  # add 0.5s to the end
    trials = trial_ieeg(raw, epoch, times, preload=True)
    # values greater than 10 standard deviations from the mean are set to NaN
    outliers_to_nan(trials, 10)
    spec = wavelet_scaleogram(trials,
                              n_jobs=1,
                              decim=20)
    # trim 0.5 seconds on the beginning and end of the data (edge artifacts)
    crop_pad(spec, "0.5s")
    out.append(spec)
resp = out[1]
base = out[0]

Used Annotations descriptions: [np.str_('Fixation'), np.str_('Go Cue'), np.str_('ISI Onset'), np.str_('Response')]
Not setting metadata
8 matching events found
No baseline correction applied
0 projection items activated
Using data from preloaded Raw for 8 events and 1301 original time points ...
0 bad epochs dropped
Data is self data: True
[Parallel(n_jobs=1)]: Done   1 tasks      | elapsed:    0.0s
[Parallel(n_jobs=1)]: Done   4 tasks      | elapsed:    0.1s
[Parallel(n_jobs=1)]: Done   7 tasks      | elapsed:    0.1s
[Parallel(n_jobs=1)]: Done  12 tasks      | elapsed:    0.1s
[Parallel(n_jobs=1)]: Done  17 tasks      | elapsed:    0.2s
[Parallel(n_jobs=1)]: Done  24 tasks      | elapsed:    0.3s
[Parallel(n_jobs=1)]: Done  31 tasks      | elapsed:    0.3s
[Parallel(n_jobs=1)]: Done  40 tasks      | elapsed:    0.4s
[Parallel(n_jobs=1)]: Done  49 tasks      | elapsed:    0.5s
[Parallel(n_jobs=1)]: Done  60 tasks      | elapsed:    0.6s
[Parallel(n_jobs=1)]: Done  71 tasks      | elapsed:    0.7s
[Parallel(n_jobs=1)]: Done  84 tasks      | elapsed:    0.9s
[Parallel(n_jobs=1)]: Done  97 tasks      | elapsed:    1.0s
[Parallel(n_jobs=1)]: Done 112 tasks      | elapsed:    1.1s
[Parallel(n_jobs=1)]: Done 119 out of 119 | elapsed:    1.2s finished
Used Annotations descriptions: [np.str_('Fixation'), np.str_('Go Cue'), np.str_('ISI Onset'), np.str_('Response')]
Not setting metadata
8 matching events found
No baseline correction applied
0 projection items activated
Using data from preloaded Raw for 8 events and 1301 original time points ...
1 bad epochs dropped
Data is self data: True
[Parallel(n_jobs=1)]: Done   1 tasks      | elapsed:    0.0s
[Parallel(n_jobs=1)]: Done   4 tasks      | elapsed:    0.0s
[Parallel(n_jobs=1)]: Done   7 tasks      | elapsed:    0.1s
[Parallel(n_jobs=1)]: Done  12 tasks      | elapsed:    0.1s
[Parallel(n_jobs=1)]: Done  17 tasks      | elapsed:    0.2s
[Parallel(n_jobs=1)]: Done  24 tasks      | elapsed:    0.2s
[Parallel(n_jobs=1)]: Done  31 tasks      | elapsed:    0.3s
[Parallel(n_jobs=1)]: Done  40 tasks      | elapsed:    0.4s
[Parallel(n_jobs=1)]: Done  49 tasks      | elapsed:    0.4s
[Parallel(n_jobs=1)]: Done  60 tasks      | elapsed:    0.5s
[Parallel(n_jobs=1)]: Done  71 tasks      | elapsed:    0.6s
[Parallel(n_jobs=1)]: Done  84 tasks      | elapsed:    0.7s
[Parallel(n_jobs=1)]: Done  97 tasks      | elapsed:    0.9s
[Parallel(n_jobs=1)]: Done 112 tasks      | elapsed:    1.0s
[Parallel(n_jobs=1)]: Done 119 out of 119 | elapsed:    1.0s finished

Time Cluster Statistics¶

mask, pvals = time_perm_cluster(resp._data, base._data,
                               p_thresh=0.1,
                               ignore_adjacency=1,  # ignore channel adjacency
                               n_perm=2000, n_jobs=1)

[Parallel(n_jobs=1)]: Done   1 tasks      | elapsed:    0.3s
[Parallel(n_jobs=1)]: Done   2 tasks      | elapsed:    0.6s
[Parallel(n_jobs=1)]: Done   3 tasks      | elapsed:    0.8s
[Parallel(n_jobs=1)]: Done   4 tasks      | elapsed:    1.1s
[Parallel(n_jobs=1)]: Done   5 tasks      | elapsed:    1.3s
[Parallel(n_jobs=1)]: Done   6 tasks      | elapsed:    1.5s
[Parallel(n_jobs=1)]: Done   7 tasks      | elapsed:    1.8s
[Parallel(n_jobs=1)]: Done   8 tasks      | elapsed:    2.0s
[Parallel(n_jobs=1)]: Done   9 tasks      | elapsed:    2.3s
[Parallel(n_jobs=1)]: Done  10 tasks      | elapsed:    2.5s
[Parallel(n_jobs=1)]: Done  11 tasks      | elapsed:    2.8s
[Parallel(n_jobs=1)]: Done  12 tasks      | elapsed:    3.0s
[Parallel(n_jobs=1)]: Done  13 tasks      | elapsed:    3.2s
[Parallel(n_jobs=1)]: Done  14 tasks      | elapsed:    3.5s
[Parallel(n_jobs=1)]: Done  15 tasks      | elapsed:    3.7s
[Parallel(n_jobs=1)]: Done  16 tasks      | elapsed:    4.0s
[Parallel(n_jobs=1)]: Done  17 tasks      | elapsed:    4.2s
[Parallel(n_jobs=1)]: Done  18 tasks      | elapsed:    4.5s
[Parallel(n_jobs=1)]: Done  19 tasks      | elapsed:    4.7s
[Parallel(n_jobs=1)]: Done  20 tasks      | elapsed:    5.0s
[Parallel(n_jobs=1)]: Done  21 tasks      | elapsed:    5.2s
[Parallel(n_jobs=1)]: Done  22 tasks      | elapsed:    5.4s
[Parallel(n_jobs=1)]: Done  23 tasks      | elapsed:    5.7s
[Parallel(n_jobs=1)]: Done  24 tasks      | elapsed:    5.9s
[Parallel(n_jobs=1)]: Done  25 tasks      | elapsed:    6.2s
[Parallel(n_jobs=1)]: Done  26 tasks      | elapsed:    6.4s
[Parallel(n_jobs=1)]: Done  27 tasks      | elapsed:    6.6s
[Parallel(n_jobs=1)]: Done  28 tasks      | elapsed:    6.9s
[Parallel(n_jobs=1)]: Done  29 tasks      | elapsed:    7.1s
[Parallel(n_jobs=1)]: Done  30 tasks      | elapsed:    7.3s
[Parallel(n_jobs=1)]: Done  31 tasks      | elapsed:    7.6s
[Parallel(n_jobs=1)]: Done  32 tasks      | elapsed:    7.8s
[Parallel(n_jobs=1)]: Done  33 tasks      | elapsed:    8.1s
[Parallel(n_jobs=1)]: Done  34 tasks      | elapsed:    8.3s
[Parallel(n_jobs=1)]: Done  35 tasks      | elapsed:    8.6s
[Parallel(n_jobs=1)]: Done  36 tasks      | elapsed:    8.8s
[Parallel(n_jobs=1)]: Done  37 tasks      | elapsed:    9.0s
[Parallel(n_jobs=1)]: Done  38 tasks      | elapsed:    9.3s
[Parallel(n_jobs=1)]: Done  39 tasks      | elapsed:    9.5s
[Parallel(n_jobs=1)]: Done  40 tasks      | elapsed:    9.7s
[Parallel(n_jobs=1)]: Done  41 tasks      | elapsed:   10.0s
[Parallel(n_jobs=1)]: Done  42 tasks      | elapsed:   10.2s
[Parallel(n_jobs=1)]: Done  43 tasks      | elapsed:   10.5s
[Parallel(n_jobs=1)]: Done  44 tasks      | elapsed:   10.7s
[Parallel(n_jobs=1)]: Done  45 tasks      | elapsed:   10.9s
[Parallel(n_jobs=1)]: Done  46 tasks      | elapsed:   11.2s
[Parallel(n_jobs=1)]: Done  47 tasks      | elapsed:   11.4s
[Parallel(n_jobs=1)]: Done  48 tasks      | elapsed:   11.7s
[Parallel(n_jobs=1)]: Done  49 tasks      | elapsed:   11.9s
[Parallel(n_jobs=1)]: Done  50 tasks      | elapsed:   12.1s
[Parallel(n_jobs=1)]: Done  51 tasks      | elapsed:   12.4s
[Parallel(n_jobs=1)]: Done  52 tasks      | elapsed:   12.6s
[Parallel(n_jobs=1)]: Done  53 tasks      | elapsed:   12.9s
[Parallel(n_jobs=1)]: Done  54 tasks      | elapsed:   13.1s
[Parallel(n_jobs=1)]: Done  55 tasks      | elapsed:   13.3s
[Parallel(n_jobs=1)]: Done  56 tasks      | elapsed:   13.6s
[Parallel(n_jobs=1)]: Done  57 tasks      | elapsed:   13.8s
[Parallel(n_jobs=1)]: Done  58 tasks      | elapsed:   14.1s
[Parallel(n_jobs=1)]: Done  59 tasks      | elapsed:   14.3s
[Parallel(n_jobs=1)]: Done  60 tasks      | elapsed:   14.5s
[Parallel(n_jobs=1)]: Done  61 tasks      | elapsed:   14.8s
[Parallel(n_jobs=1)]: Done  62 tasks      | elapsed:   15.0s
[Parallel(n_jobs=1)]: Done  63 tasks      | elapsed:   15.3s
[Parallel(n_jobs=1)]: Done  64 tasks      | elapsed:   15.5s
[Parallel(n_jobs=1)]: Done  65 tasks      | elapsed:   15.7s
[Parallel(n_jobs=1)]: Done  66 tasks      | elapsed:   16.0s
[Parallel(n_jobs=1)]: Done  67 tasks      | elapsed:   16.2s
[Parallel(n_jobs=1)]: Done  68 tasks      | elapsed:   16.4s
[Parallel(n_jobs=1)]: Done  69 tasks      | elapsed:   16.7s
[Parallel(n_jobs=1)]: Done  70 tasks      | elapsed:   16.9s
[Parallel(n_jobs=1)]: Done  71 tasks      | elapsed:   17.2s
[Parallel(n_jobs=1)]: Done  72 tasks      | elapsed:   17.4s
[Parallel(n_jobs=1)]: Done  73 tasks      | elapsed:   17.6s
[Parallel(n_jobs=1)]: Done  74 tasks      | elapsed:   17.9s
[Parallel(n_jobs=1)]: Done  75 tasks      | elapsed:   18.1s
[Parallel(n_jobs=1)]: Done  76 tasks      | elapsed:   18.4s
[Parallel(n_jobs=1)]: Done  77 tasks      | elapsed:   18.6s
[Parallel(n_jobs=1)]: Done  78 tasks      | elapsed:   18.9s
[Parallel(n_jobs=1)]: Done  79 tasks      | elapsed:   19.1s
[Parallel(n_jobs=1)]: Done  80 tasks      | elapsed:   19.3s
[Parallel(n_jobs=1)]: Done  81 tasks      | elapsed:   19.6s
[Parallel(n_jobs=1)]: Done  82 tasks      | elapsed:   19.8s
[Parallel(n_jobs=1)]: Done  83 tasks      | elapsed:   20.1s
[Parallel(n_jobs=1)]: Done  84 tasks      | elapsed:   20.3s
[Parallel(n_jobs=1)]: Done  85 tasks      | elapsed:   20.5s
[Parallel(n_jobs=1)]: Done  86 tasks      | elapsed:   20.8s
[Parallel(n_jobs=1)]: Done  87 tasks      | elapsed:   21.0s
[Parallel(n_jobs=1)]: Done  88 tasks      | elapsed:   21.3s
[Parallel(n_jobs=1)]: Done  89 tasks      | elapsed:   21.5s
[Parallel(n_jobs=1)]: Done  90 tasks      | elapsed:   21.8s
[Parallel(n_jobs=1)]: Done  91 tasks      | elapsed:   22.0s
[Parallel(n_jobs=1)]: Done  92 tasks      | elapsed:   22.2s
[Parallel(n_jobs=1)]: Done  93 tasks      | elapsed:   22.5s
[Parallel(n_jobs=1)]: Done  94 tasks      | elapsed:   22.7s
[Parallel(n_jobs=1)]: Done  95 tasks      | elapsed:   22.9s
[Parallel(n_jobs=1)]: Done  96 tasks      | elapsed:   23.2s
[Parallel(n_jobs=1)]: Done  97 tasks      | elapsed:   23.4s
[Parallel(n_jobs=1)]: Done  98 tasks      | elapsed:   23.7s
[Parallel(n_jobs=1)]: Done  99 tasks      | elapsed:   23.9s
[Parallel(n_jobs=1)]: Done 100 tasks      | elapsed:   24.1s
[Parallel(n_jobs=1)]: Done 101 tasks      | elapsed:   24.4s
[Parallel(n_jobs=1)]: Done 102 tasks      | elapsed:   24.6s
[Parallel(n_jobs=1)]: Done 103 tasks      | elapsed:   24.9s
[Parallel(n_jobs=1)]: Done 104 tasks      | elapsed:   25.1s
[Parallel(n_jobs=1)]: Done 105 tasks      | elapsed:   25.3s
[Parallel(n_jobs=1)]: Done 106 tasks      | elapsed:   25.6s
[Parallel(n_jobs=1)]: Done 107 tasks      | elapsed:   25.8s
[Parallel(n_jobs=1)]: Done 108 tasks      | elapsed:   26.1s
[Parallel(n_jobs=1)]: Done 109 tasks      | elapsed:   26.3s
[Parallel(n_jobs=1)]: Done 110 tasks      | elapsed:   26.6s
[Parallel(n_jobs=1)]: Done 111 tasks      | elapsed:   26.8s
[Parallel(n_jobs=1)]: Done 112 tasks      | elapsed:   27.0s
[Parallel(n_jobs=1)]: Done 113 tasks      | elapsed:   27.3s
[Parallel(n_jobs=1)]: Done 114 tasks      | elapsed:   27.5s
[Parallel(n_jobs=1)]: Done 115 tasks      | elapsed:   27.8s
[Parallel(n_jobs=1)]: Done 116 tasks      | elapsed:   28.0s
[Parallel(n_jobs=1)]: Done 117 tasks      | elapsed:   28.2s
[Parallel(n_jobs=1)]: Done 118 tasks      | elapsed:   28.5s
[Parallel(n_jobs=1)]: Done 119 tasks      | elapsed:   28.7s
[Parallel(n_jobs=1)]: Done 119 out of 119 | elapsed:   28.7s finished

Plot the Time-Frequency Clusters¶

fig, axs = plt.subplots(5, 24, figsize=(20, 20))
for i, ax in enumerate(axs.flat):
    if i >= mask.shape[0]:
        ax.axis('off')
        continue
    ax.imshow(mask[i])
    ax.set_title(resp.info['ch_names'][i])