- MNE-Python version: 0.21.0
- operating system: Windows10
Hi,
my script below aims in sequence toperform the following actions:
import raw data > filtering > run the ICA and reconstruct the signal > apply the CAR > create epochs from a single event > find the epochs to reject because of artifacts. In my example, each epoch is equal to each second of the recording. Regarding the ICA, I used the ICA version suggested by @larsoner that generates the same decomposition each time the script is run.
MY QUESTIONS:
- Considering that Iâve only frontal channels (if you import the mne raw object through my script youâll see all info about channels and more), do you think that my approach for the ICA is correct? Please see the â# ICA SECTIONâ in the script.
- The final result of the script is represented by the percentage of epochs to rejects. In general I would expect that if I use the ICA (versus not using the ICA) the percentage of rejected epochs would be reduced, but it is not always like that. I did several tests running the analysis with all the following combinations, where for each combination I report the percentage of epochs rejected. From my tests, the better combination in order to have less epoch rejected is the âDâ, where I run the ICA with the reject parameter enabled in the ica.fit, along with the CAR. From your experience, do you think that such results are reasonable? Do you have some tips to improve my pipeline? In general, whatever suggestion will be really appreciated.
For the following 4 tests, if ICA was run, the ârejectâ ICA parameter WAS used:
- A) without ICA, without CAR = 80.1% epochs rejected
- B) ICA, without CAR = 40.3 % epochs rejected
- C) without ICA, CAR = 80.9% epochs rejected
- D) ICA, CAR = 2.2% epochs rejected
For the following 2 tests ICA was run excluding the ârejectâ ICA parameter from ica.fit :
- E) ICA, without CAR = 40.3 % epochs rejected
- F) ICA, CAR = 61% epochs rejected
N.B tests âBâ and âEâ led to same results, while âDâ and âFâ led to different results.
DATA AND SCRIPT
- Here mne raw object to download (the raw object of one subject, 4 MB)
- The entire script is shown below (the only change required is the path to the folder where you donwloaded the file)
import os
import mne
# LOAD THE EEG RECORDING
sample_data_folder = r"F:\folderContaningData" # change the path to the folder where you downloaded the file
sample_data_raw_file = os.path.join(sample_data_folder, 'exp_c3r1c1m4a1o_conc.fif')
raw = mne.io.read_raw_fif(sample_data_raw_file, preload=True)
# FILTERING
lowPass = 30 # filtering
highPass = 1 # filtering
raw_filtered = raw.copy().filter(l_freq=highPass, h_freq=lowPass)
# --------------------- START ICA SECTION
from mne.preprocessing import (ICA)
n_components = None # if float, select n_components by explained variance of PCA
method = 'fastica' # for comparison with EEGLAB try "extended-infomax" here
decim = 3 # we need sufficient statistics, not all time points -> saves time
# we will also set state of the random number generator - ICA is a
# non-deterministic algorithm, but we want to have the same decomposition
# and the same order of components each time this tutorial is run
random_state = 23
ica = ICA(n_components=n_components, method=method, random_state=random_state)
reject = dict(mag=5e-12, grad=4000e-13) # ???
ica.fit(raw_filtered, decim=decim, reject=reject)
# ica.fit(raw_filtered, decim=decim) # enable this line and disable the above line, if you want to run ICA without reject parameter
print(ica)
# plot components
raw_filtered.plot(title='EEG before ICA')
ica.plot_components()
ica.plot_sources(raw_filtered)
# select and exclude component, reconstruct and plot
ICAcomponents2Exlude = [0]
ica.plot_overlay(raw_filtered, exclude=ICAcomponents2Exlude, picks='eeg')
ica.exclude = ICAcomponents2Exlude
reconst_raw = raw_filtered.copy()
ica.apply(reconst_raw) # Apply ICA for remove the components
reconst_raw.plot(title='EEG after ICA')
raw_filtered.plot(title='EEG before ICA')
# --------------------- END ICA SECTION
# Enable the following line in case you don't run the ICA section
# reconst_raw = raw_filtered.copy()
# --------------------- START CAR (common average reference) SECTION
reconst_raw.set_eeg_reference('average', projection=True)
reconst_raw.apply_proj()
reconst_raw.plot(title='EEG after CAR')
print(reconst_raw.info)
# --------------------- END CAR SECTION
# Enable the following line in case you don't run BOTH the ICA section and the CAR section
# reconst_raw = raw_filtered.copy()
# --------------------- START EVENT-EPOCHS (ONLY ONE EVENT, see 'My Event': 1)
# each second of the recording will be equal to one epoch
import numpy as np
from autoreject import get_rejection_threshold # noqa
sfreq = 250 # eeg data sample frequency
event_id = 1 # This is used to identify the events.
events_raw = []
total_sec = len(reconst_raw) / sfreq
for i in range(0, int(total_sec)):
event_eachSecond = [i*sfreq, 0, event_id]
events_raw.append(event_eachSecond)
events = np.array(events_raw)
mne.viz.plot_events(events, sfreq=reconst_raw.info['sfreq']) # plot events
# from events to epochs
event_id = {'My Event': 1}
epochs = mne.Epochs(reconst_raw, events, event_id, tmin=0, tmax=1,baseline=(0, 0), preload=True)
reject = get_rejection_threshold(epochs, decim=2)
print('The rejection dictionary is %s' % reject)
epochs.drop_bad(reject=reject)
epochs.average().plot()
mne.viz.plot_epochs(epochs)
# FINAL RESULTS (PERCENTAGE OF REJECTED EPOCHS)
epochs.plot_drop_log() # A summary of the causes of rejected epochs can be generated with the plot_drop_log() method. If it gives the name of channel such as âF8â it means the epoch was rejected because âF8â exceeded the peak-to-peak rejection limit.
Thank you in advance