Time Permutation Cluster Statistics¶

Below is a code sample for performing a permutation test with cluster correction, as well as using the window averaged test.

from ieeg.navigate import channel_outlier_marker, trial_ieeg
from ieeg.timefreq.utils import crop_pad
from ieeg.timefreq import gamma
from ieeg.calc import stats
from ieeg.calc.fast import mean_diff
import matplotlib.pyplot as plt
import mne

Load Data¶

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

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.

Remove bad channels¶

mark channel outliers as bad

raw.info['bads'] = channel_outlier_marker(raw, 3)

# Exclude bad channels
raw.drop_channels(raw.info['bads'])
raw.load_data()

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

outlier round 1 channels: ['LAMY 7']
outlier round 1 channels: ['LAMY 7', 'LPHG 6']
outlier round 1 channels: ['LAMY 7', 'LPHG 6', 'LBRI 3']
outlier round 1 channels: ['LAMY 7', 'LPHG 6', 'LBRI 3', 'RAHP 3']
outlier round 2 channels: ['LAMY 7', 'LPHG 6', 'LBRI 3', 'RAHP 3', 'LENT 3']
outlier round 2 channels: ['LAMY 7', 'LPHG 6', 'LBRI 3', 'RAHP 3', 'LENT 3', 'LPIT 5']
outlier round 3 channels: ['LAMY 7', 'LPHG 6', 'LBRI 3', 'RAHP 3', 'LENT 3', 'LPIT 5', 'LENT 2']
outlier round 3 channels: ['LAMY 7', 'LPHG 6', 'LBRI 3', 'RAHP 3', 'LENT 3', 'LPIT 5', 'LENT 2', 'LENT 7']
outlier round 4 channels: ['LAMY 7', 'LPHG 6', 'LBRI 3', 'RAHP 3', 'LENT 3', 'LPIT 5', 'LENT 2', 'LENT 7', 'LPLI 8']
outlier round 5 channels: ['LAMY 7', 'LPHG 6', 'LBRI 3', 'RAHP 3', 'LENT 3', 'LPIT 5', 'LENT 2', 'LENT 7', 'LPLI 8', 'RAHP 2']
outlier round 6 channels: ['LAMY 7', 'LPHG 6', 'LBRI 3', 'RAHP 3', 'LENT 3', 'LPIT 5', 'LENT 2', 'LENT 7', 'LPLI 8', 'RAHP 2', 'LPIT 6']
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

Gamma Filtering¶

out = []
for epoch, t in zip(('Fixation', 'Response'),
                    ((-0.4, 0), (-0.1, 0.3))):
    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)
    gamma.extract(trials, copy=False, n_jobs=1)
    # trim 0.5 seconds on the beginning and end of the data (edge artifacts)
    crop_pad(trials, "0.5s")
    out.append(trials)
resp = out[1]
resp.decimate(2)
base = out[0]
base.decimate(2)

del raw

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 1400 original time points ...
0 bad epochs dropped

  0%|          | 0/8 [00:00<?, ?it/s]
 75%|███████▌  | 6/8 [00:00<00:00, 51.34it/s]
100%|██████████| 8/8 [00:00<00:00, 51.99it/s]
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 1401 original time points ...
1 bad epochs dropped

  0%|          | 0/7 [00:00<?, ?it/s]
 57%|█████▋    | 4/7 [00:00<00:00, 36.99it/s]
100%|██████████| 7/7 [00:00<00:00, 37.00it/s]
/home/docs/checkouts/readthedocs.org/user_builds/ieeg-pipelines/checkouts/latest/examples/plot_stats.py:51: RuntimeWarning: The measurement information indicates a low-pass frequency of 499.7060546875 Hz. The decim=2 parameter will result in a sampling frequency of 499.7060546875 Hz, which can cause aliasing artifacts.
  resp.decimate(2)
/home/docs/checkouts/readthedocs.org/user_builds/ieeg-pipelines/checkouts/latest/examples/plot_stats.py:53: RuntimeWarning: The measurement information indicates a low-pass frequency of 499.7060546875 Hz. The decim=2 parameter will result in a sampling frequency of 499.7060546875 Hz, which can cause aliasing artifacts.
  base.decimate(2)

Time Cluster Stats¶

mask1, _ = stats.time_perm_cluster(resp._data, base._data,
                                p_thresh=0.05,
                                axis=0,
                                n_perm=1000,
                                n_jobs=1,
                                ignore_adjacency=1)
fig1 = plt.imshow(mask1)

[Parallel(n_jobs=1)]: Done   1 tasks      | elapsed:    0.1s
[Parallel(n_jobs=1)]: Done   2 tasks      | elapsed:    0.1s
[Parallel(n_jobs=1)]: Done   3 tasks      | elapsed:    0.1s
[Parallel(n_jobs=1)]: Done   4 tasks      | elapsed:    0.1s
[Parallel(n_jobs=1)]: Done   5 tasks      | elapsed:    0.2s
[Parallel(n_jobs=1)]: Done   6 tasks      | elapsed:    0.2s
[Parallel(n_jobs=1)]: Done   7 tasks      | elapsed:    0.2s
[Parallel(n_jobs=1)]: Done   8 tasks      | elapsed:    0.2s
[Parallel(n_jobs=1)]: Done   9 tasks      | elapsed:    0.3s
[Parallel(n_jobs=1)]: Done  10 tasks      | elapsed:    0.3s
[Parallel(n_jobs=1)]: Done  11 tasks      | elapsed:    0.3s
[Parallel(n_jobs=1)]: Done  12 tasks      | elapsed:    0.3s
[Parallel(n_jobs=1)]: Done  13 tasks      | elapsed:    0.4s
[Parallel(n_jobs=1)]: Done  14 tasks      | elapsed:    0.4s
[Parallel(n_jobs=1)]: Done  15 tasks      | elapsed:    0.4s
[Parallel(n_jobs=1)]: Done  16 tasks      | elapsed:    0.4s
[Parallel(n_jobs=1)]: Done  17 tasks      | elapsed:    0.5s
[Parallel(n_jobs=1)]: Done  18 tasks      | elapsed:    0.5s
[Parallel(n_jobs=1)]: Done  19 tasks      | elapsed:    0.5s
[Parallel(n_jobs=1)]: Done  20 tasks      | elapsed:    0.5s
[Parallel(n_jobs=1)]: Done  21 tasks      | elapsed:    0.6s
[Parallel(n_jobs=1)]: Done  22 tasks      | elapsed:    0.6s
[Parallel(n_jobs=1)]: Done  23 tasks      | elapsed:    0.6s
[Parallel(n_jobs=1)]: Done  24 tasks      | elapsed:    0.7s
[Parallel(n_jobs=1)]: Done  25 tasks      | elapsed:    0.7s
[Parallel(n_jobs=1)]: Done  26 tasks      | elapsed:    0.7s
[Parallel(n_jobs=1)]: Done  27 tasks      | elapsed:    0.7s
[Parallel(n_jobs=1)]: Done  28 tasks      | elapsed:    0.8s
[Parallel(n_jobs=1)]: Done  29 tasks      | elapsed:    0.8s
[Parallel(n_jobs=1)]: Done  30 tasks      | elapsed:    0.8s
[Parallel(n_jobs=1)]: Done  31 tasks      | elapsed:    0.8s
[Parallel(n_jobs=1)]: Done  32 tasks      | elapsed:    0.9s
[Parallel(n_jobs=1)]: Done  33 tasks      | elapsed:    0.9s
[Parallel(n_jobs=1)]: Done  34 tasks      | elapsed:    0.9s
[Parallel(n_jobs=1)]: Done  35 tasks      | elapsed:    0.9s
[Parallel(n_jobs=1)]: Done  36 tasks      | elapsed:    1.0s
[Parallel(n_jobs=1)]: Done  37 tasks      | elapsed:    1.0s
[Parallel(n_jobs=1)]: Done  38 tasks      | elapsed:    1.0s
[Parallel(n_jobs=1)]: Done  39 tasks      | elapsed:    1.0s
[Parallel(n_jobs=1)]: Done  40 tasks      | elapsed:    1.1s
[Parallel(n_jobs=1)]: Done  41 tasks      | elapsed:    1.1s
[Parallel(n_jobs=1)]: Done  42 tasks      | elapsed:    1.1s
[Parallel(n_jobs=1)]: Done  43 tasks      | elapsed:    1.1s
[Parallel(n_jobs=1)]: Done  44 tasks      | elapsed:    1.2s
[Parallel(n_jobs=1)]: Done  45 tasks      | elapsed:    1.2s
[Parallel(n_jobs=1)]: Done  46 tasks      | elapsed:    1.2s
[Parallel(n_jobs=1)]: Done  47 tasks      | elapsed:    1.2s
[Parallel(n_jobs=1)]: Done  48 tasks      | elapsed:    1.3s
[Parallel(n_jobs=1)]: Done  49 tasks      | elapsed:    1.3s
[Parallel(n_jobs=1)]: Done  50 tasks      | elapsed:    1.3s
[Parallel(n_jobs=1)]: Done  51 tasks      | elapsed:    1.3s
[Parallel(n_jobs=1)]: Done  52 tasks      | elapsed:    1.4s
[Parallel(n_jobs=1)]: Done  53 tasks      | elapsed:    1.4s
[Parallel(n_jobs=1)]: Done  54 tasks      | elapsed:    1.4s
[Parallel(n_jobs=1)]: Done  55 tasks      | elapsed:    1.5s
[Parallel(n_jobs=1)]: Done  56 tasks      | elapsed:    1.5s
[Parallel(n_jobs=1)]: Done  57 tasks      | elapsed:    1.5s
[Parallel(n_jobs=1)]: Done  58 tasks      | elapsed:    1.5s
[Parallel(n_jobs=1)]: Done  59 tasks      | elapsed:    1.6s
[Parallel(n_jobs=1)]: Done  60 tasks      | elapsed:    1.6s
[Parallel(n_jobs=1)]: Done  61 tasks      | elapsed:    1.6s
[Parallel(n_jobs=1)]: Done  62 tasks      | elapsed:    1.6s
[Parallel(n_jobs=1)]: Done  63 tasks      | elapsed:    1.7s
[Parallel(n_jobs=1)]: Done  64 tasks      | elapsed:    1.7s
[Parallel(n_jobs=1)]: Done  65 tasks      | elapsed:    1.7s
[Parallel(n_jobs=1)]: Done  66 tasks      | elapsed:    1.7s
[Parallel(n_jobs=1)]: Done  67 tasks      | elapsed:    1.8s
[Parallel(n_jobs=1)]: Done  68 tasks      | elapsed:    1.8s
[Parallel(n_jobs=1)]: Done  69 tasks      | elapsed:    1.8s
[Parallel(n_jobs=1)]: Done  70 tasks      | elapsed:    1.8s
[Parallel(n_jobs=1)]: Done  71 tasks      | elapsed:    1.9s
[Parallel(n_jobs=1)]: Done  72 tasks      | elapsed:    1.9s
[Parallel(n_jobs=1)]: Done  73 tasks      | elapsed:    1.9s
[Parallel(n_jobs=1)]: Done  74 tasks      | elapsed:    1.9s
[Parallel(n_jobs=1)]: Done  75 tasks      | elapsed:    2.0s
[Parallel(n_jobs=1)]: Done  76 tasks      | elapsed:    2.0s
[Parallel(n_jobs=1)]: Done  77 tasks      | elapsed:    2.0s
[Parallel(n_jobs=1)]: Done  78 tasks      | elapsed:    2.0s
[Parallel(n_jobs=1)]: Done  79 tasks      | elapsed:    2.1s
[Parallel(n_jobs=1)]: Done  80 tasks      | elapsed:    2.1s
[Parallel(n_jobs=1)]: Done  81 tasks      | elapsed:    2.1s
[Parallel(n_jobs=1)]: Done  82 tasks      | elapsed:    2.1s
[Parallel(n_jobs=1)]: Done  83 tasks      | elapsed:    2.2s
[Parallel(n_jobs=1)]: Done  84 tasks      | elapsed:    2.2s
[Parallel(n_jobs=1)]: Done  85 tasks      | elapsed:    2.2s
[Parallel(n_jobs=1)]: Done  86 tasks      | elapsed:    2.3s
[Parallel(n_jobs=1)]: Done  87 tasks      | elapsed:    2.3s
[Parallel(n_jobs=1)]: Done  88 tasks      | elapsed:    2.3s
[Parallel(n_jobs=1)]: Done  89 tasks      | elapsed:    2.3s
[Parallel(n_jobs=1)]: Done  90 tasks      | elapsed:    2.4s
[Parallel(n_jobs=1)]: Done  91 tasks      | elapsed:    2.4s
[Parallel(n_jobs=1)]: Done  92 tasks      | elapsed:    2.4s
[Parallel(n_jobs=1)]: Done  93 tasks      | elapsed:    2.4s
[Parallel(n_jobs=1)]: Done  94 tasks      | elapsed:    2.5s
[Parallel(n_jobs=1)]: Done  95 tasks      | elapsed:    2.5s
[Parallel(n_jobs=1)]: Done  96 tasks      | elapsed:    2.5s
[Parallel(n_jobs=1)]: Done  97 tasks      | elapsed:    2.5s
[Parallel(n_jobs=1)]: Done  98 tasks      | elapsed:    2.6s
[Parallel(n_jobs=1)]: Done  99 tasks      | elapsed:    2.6s
[Parallel(n_jobs=1)]: Done 100 tasks      | elapsed:    2.6s
[Parallel(n_jobs=1)]: Done 101 tasks      | elapsed:    2.6s
[Parallel(n_jobs=1)]: Done 102 tasks      | elapsed:    2.7s
[Parallel(n_jobs=1)]: Done 103 tasks      | elapsed:    2.7s
[Parallel(n_jobs=1)]: Done 104 tasks      | elapsed:    2.7s
[Parallel(n_jobs=1)]: Done 105 tasks      | elapsed:    2.7s
[Parallel(n_jobs=1)]: Done 106 tasks      | elapsed:    2.8s
[Parallel(n_jobs=1)]: Done 107 tasks      | elapsed:    2.8s
[Parallel(n_jobs=1)]: Done 108 tasks      | elapsed:    2.8s
[Parallel(n_jobs=1)]: Done 108 out of 108 | elapsed:    2.8s finished

Windowed Average Test (alternative test)¶

# response data points are length 201, while baseline is 200. We need to
# trim the response data to match the baseline with [..., :-1]
mask2 = stats.window_averaged_shuffle(resp._data[..., :-1], base._data,
                                      n_perm=1000,
                                      stat_func=mean_diff)
fig2 = plt.imshow(0.05 >= mask2[:, None])