Autorejection on resting state EEG data

💬 Support & Discussions
, ,
1

I am trying to perform autorejection on resting state EEG data (continuous) . I created epochs on my data and tried to identify the rejection threshold. I wanted to know how the rejection threshold is automatically created or what parameters it includes when I use – get_rejection_threshold(epochs). With my code I got eeg threshold of {‘eeg’: 0.0005535939698374599}. I thought that was too high ?

Also, is autorejection preferably performed after removing slow drifts , or before ICA or after ICA ?

Here is my code, is this the way it is done ? :

raw=mne.io.read_raw_eeglab('samplefile.set',preload=True)
raw_avg = raw.set_eeg_reference("average")
raw_drift=raw_avg.filter(0.5, None, fir_design='firwin') #removedrifts
epochs=mne.make_fixed_length_epochs(raw_drift, duration=1.0, preload=False, proj=True)
reject = get_rejection_threshold(epochs)
print('The rejection dictionary is %s' % reject)
epochs.drop_bad(reject=reject)
#reject_criteria = dict( eeg=150e-6)e
2

I was trying to mimic the autorejection from eeglab, where you make 1s epochs and then perform rejkurt (https://sccn.ucsd.edu/~arno/eeglab/auto/pop_rejkurt.html) at 5 std threshold, and jointprob (https://sccn.ucsd.edu/~arno/eeglab/auto/pop_jointprob.html) at 5std threshold.

3

I’m just tagging @mainakjas, who would probably be the most competent person around here when it comes to autoreject :ok_hand:

4

Hi Apoorva,

If you can share the file, I can take a quick look. 500 uV does appear too high. However, I normally suggest looking at the data to ensure what autoreject makes sense or not.

It’s better to apply autoreject after removing slow drifts as otherwise you might drop too many epochs. As for ICA, I would suggest after as you don’t want to drop epochs due to eyeblinks. However, you can estimate the thresholds before and supply to ICA if you want to help improve the ICA decomposition. See detailed FAQ here: http://autoreject.github.io/faq.html#should-i-apply-ica-first-or-autoreject-first

The EEGLAB links you shared seem to be for rejecting ICA components automatically. That is a different usecase than what autoreject targets, that is to remove sensor artifacts.

Mainak

5

Hi @mainakjas ,

Thank you so much. I can send a different one. This is a resting state eeg file. By running this code, it should download the file and run the autoreject. For this one I got a threshold of {‘eeg’: 0.00012221068272366844}. I just wanted more clarification on what you mean by “looking at the data to ensure what autoreject makes sense or not.” .

import os
import os.path as op
import openneuro
import mne
import numpy as np
from autoreject import AutoReject
from autoreject import get_rejection_threshold

dataset = 'ds002778'
subject = 'pd23'
openneuro.download(dataset=dataset, target_dir='.',
                   include=[f'sub-{subject}'])
raw = mne.io.read_raw_bdf('./sub-pd23/ses-off/eeg/sub-pd23_ses-off_task-rest_eeg.bdf', preload=True)
raw.drop_channels(['EXG1', 'EXG2', 'EXG3', 'EXG4','EXG5', 'EXG6', 'EXG7', 'EXG8',
                                  'Status']) #drop extra channels
raw.set_eeg_reference(ref_channels='average')
raw.filter(0.5, None, fir_design='firwin',phase='zero-double') #remove drifts

epochs=mne.make_fixed_length_epochs(raw, duration=1.0, preload=False, proj=True)
reject = get_rejection_threshold(epochs)
print('The rejection dictionary is %s' % reject)
epochs.drop_bad(reject=reject)

With regards to EEGLAB, the documentation here shows to use joinprob and rejkurt for abnormally distributed data as method for rejecting artifacts in continuous or epoched data. Legacy rejection - EEGLAB Wiki

6

I usually subtract blink and saccade ICA components before running autoreject so that epochs contaminated by those artifacts can be saved rather than being rejected.

I think this is a good example of how to check that autoreject is working Visualize bad sensors per trial — autoreject 0.2.1 documentation.

7

Hi Apoorva,

Apologies for my slow response. Thanks for the reproducible script, that was very helpful. Try the following:


epochs = mne.make_fixed_length_epochs(raw, duration=1.0, 
                                                                     preload=False, proj=True)
epochs_orig = epochs.copy()
reject = get_rejection_threshold(epochs)
print('The rejection dictionary is %s' % reject)
epochs.drop_bad(reject=reject)
bad_idx = [idx for idx, ch_log in 
                  enumerate(epochs.drop_log) if ch_log != ()]
epochs_orig[bad_idx].plot()

I just replaced the last few lines of your code. As you can see from the plot, most of the epochs do indeed look bad:

image
image1374×777 551 KB

What Alex suggested will also be useful if you use autoreject (local). Perhaps we can also add this small code snippet I shared as part of an example.

Best,
Mainak

2 Likes