Find global rejection threshold

This example demonstrates how to use autoreject to find global rejection thresholds.

# Author: Mainak Jas <mainak.jas@telecom-paristech.fr>
# License: BSD-3-Clause

Let us import the data using MNE-Python and epoch it.

import mne
from mne import io
from mne.datasets import sample

event_id = {'Visual/Left': 3}
tmin, tmax = -0.2, 0.5

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'

raw = io.read_raw_fif(raw_fname, preload=True)
events = mne.read_events(event_fname)

include = []
picks = mne.pick_types(raw.info, meg=True, eeg=True, stim=False,
                       eog=True, include=include, exclude='bads')
epochs = mne.Epochs(raw, events, event_id, tmin, tmax,
                    picks=picks, baseline=(None, 0), preload=True,
                    reject=None, verbose=False, detrend=1)

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...

Now we get the rejection dictionary

from autoreject import get_rejection_threshold  # noqa

# We can use the `decim` parameter to only take every nth time slice.
# This speeds up the computation time. Note however that for low sampling
# rates and high decimation parameters, you might not detect "peaky artifacts"
# (with a fast timecourse) in your data. A low amount of decimation however is
# almost always beneficial at no decrease of accuracy.
reject = get_rejection_threshold(epochs, decim=2)

/home/circleci/project/autoreject/autoreject.py:242: RuntimeWarning: The measurement information indicates a low-pass frequency of 40.0 Hz. The decim=2 parameter will result in a sampling frequency of 75.07687377929688 Hz, which can cause aliasing artifacts.
  epochs.decimate(decim=decim)
Estimating rejection dictionary for mag
Estimating rejection dictionary for grad
Estimating rejection dictionary for eeg
Estimating rejection dictionary for eog

and print it

print('The rejection dictionary is %s' % reject)

The rejection dictionary is {'mag': np.float64(3.284586001851895e-12), 'grad': np.float64(1.0357642173052411e-10), 'eeg': np.float64(5.844322373549318e-05), 'eog': np.float64(6.559782003227096e-05)}

Finally, the cleaned epochs

epochs.drop_bad(reject=reject)
epochs.average().plot()

    Rejecting  epoch based on EEG : ['EEG 001', 'EEG 002', 'EEG 003', 'EEG 006', 'EEG 007']
    Rejecting  epoch based on EEG : ['EEG 001', 'EEG 002', 'EEG 003', 'EEG 006', 'EEG 007', 'EEG 014', 'EEG 015']
    Rejecting  epoch based on EEG : ['EEG 007']
    Rejecting  epoch based on EOG : ['EOG 061']
    Rejecting  epoch based on EEG : ['EEG 002', 'EEG 003', 'EEG 007']
    Rejecting  epoch based on EEG : ['EEG 001', 'EEG 003', 'EEG 007']
    Rejecting  epoch based on EEG : ['EEG 025']
    Rejecting  epoch based on EEG : ['EEG 001', 'EEG 002', 'EEG 003', 'EEG 006', 'EEG 007', 'EEG 008']
    Rejecting  epoch based on EEG : ['EEG 001', 'EEG 002', 'EEG 003', 'EEG 007']
    Rejecting  epoch based on EEG : ['EEG 001', 'EEG 004', 'EEG 007', 'EEG 015']
    Rejecting  epoch based on EEG : ['EEG 001', 'EEG 002', 'EEG 003', 'EEG 007']
    Rejecting  epoch based on EEG : ['EEG 007']
    Rejecting  epoch based on EEG : ['EEG 001', 'EEG 002', 'EEG 003', 'EEG 006', 'EEG 007']
    Rejecting  epoch based on MAG : ['MEG 1331']
    Rejecting  epoch based on EEG : ['EEG 004']
    Rejecting  epoch based on EEG : ['EEG 001', 'EEG 002', 'EEG 003', 'EEG 007']
    Rejecting  epoch based on EEG : ['EEG 007', 'EEG 016', 'EEG 024']
17 bad epochs dropped

<Figure size 640x500 with 6 Axes>