Epoching and averaging (ERP/ERF)

import os.path as op
import numpy as np

import mne

In MNE, epochs refers to a collection of single trials or short segments of time locked raw data. If you haven’t already, you might want to check out The Epochs data structure: epoched data. In this tutorial we take a deeper look into construction of epochs and averaging the epoch data to evoked instances. First let’s read in the raw sample data.

data_path = mne.datasets.sample.data_path()
fname = op.join(data_path, 'MEG', 'sample', 'sample_audvis_raw.fif')
raw = mne.io.read_raw_fif(fname)
raw.set_eeg_reference('average', projection=True)  # set EEG average reference

To create time locked epochs, we first need a set of events that contain the information about the times. In this tutorial we use the stimulus channel to define the events. Let’s look at the raw data.

order = np.arange(raw.info['nchan'])
order[9] = 312  # We exchange the plotting order of two channels
order[312] = 9  # to show the trigger channel as the 10th channel.
raw.plot(n_channels=10, order=order, block=True)
../_images/sphx_glr_plot_epoching_and_averaging_001.png

Notice channel STI 014 at the bottom. It is the trigger channel that was used for combining all the events to a single channel. We can see that it has several pulses of different amplitude throughout the recording. These pulses correspond to different stimuli presented to the subject during the acquisition. The pulses have values of 1, 2, 3, 4, 5 and 32. These are the events we are going to align the epochs to. To create an event list from raw data, we simply call a function dedicated just for that. Since the event list is simply a numpy array, you can also manually create one. If you create one from an outside source (like a separate file of events), pay special attention in aligning the events correctly with the raw data.

events = mne.find_events(raw)
print('Found %s events, first five:' % len(events))
print(events[:5])

# Plot the events to get an idea of the paradigm
# Specify colors and an event_id dictionary for the legend.
event_id = {'Auditory/Left': 1, 'Auditory/Right': 2,
            'Visual/Left': 3, 'Visual/Right': 4,
            'smiley': 5, 'button': 32}
color = {1: 'green', 2: 'yellow', 3: 'red', 4: 'c', 5: 'black', 32: 'blue'}

mne.viz.plot_events(events, raw.info['sfreq'], raw.first_samp, color=color,
                    event_id=event_id)
../_images/sphx_glr_plot_epoching_and_averaging_002.png

Out:

Found 320 events, first five:
[[27977     0     2]
 [28345     0     3]
 [28771     0     1]
 [29219     0     4]
 [29652     0     2]]

The event list contains three columns. The first column corresponds to sample number. To convert this to seconds, you should divide the sample number by the used sampling frequency. The second column is reserved for the old value of the trigger channel at the time of transition, but is currently not in use. The third column is the trigger id (amplitude of the pulse).

You might wonder why the samples don’t seem to align with the plotted data. For instance, the first event has a sample number of 27977 which should translate to roughly 46.6 seconds (27977 / 600). However looking at the pulses we see the first pulse at 3.6 seconds. This is because Neuromag recordings have an attribute first_samp which refers to the offset between the system start and the start of the recording. Our data has a first_samp equal to 25800. This means that the first sample you see with raw.plot is the sample number 25800. Generally you don’t need to worry about this offset as it is taken into account with MNE functions, but it is good to be aware of. Just to confirm, let’s plot the events together with the raw data. Notice how the vertical lines (events) align nicely with the pulses on STI 014.

raw.plot(events=events, n_channels=10, order=order)
../_images/sphx_glr_plot_epoching_and_averaging_003.png

In this tutorial we are only interested in triggers 1, 2, 3 and 4. These triggers correspond to auditory and visual stimuli. The event_id here can be an int, a list of ints or a dict. With dicts it is possible to assign these ids to distinct categories. When using ints or lists this information is lost. First we shall define some parameters to feed to the mne.Epochs constructor. The values tmin and tmax refer to offsets in relation to the events. Here we make epochs that collect the data from 200 ms before to 500 ms after the event.

tmin, tmax = -0.2, 0.5
event_id = {'Auditory/Left': 1, 'Auditory/Right': 2,
            'Visual/Left': 3, 'Visual/Right': 4}
# Only pick MEG and EOG channels.
picks = mne.pick_types(raw.info, meg=True, eeg=False, eog=True)

Now we have everything we need to construct the epochs. To get some meaningful results, we also want to baseline the epochs. Baselining computes the mean over the baseline period and adjusts the data accordingly. The epochs constructor uses a baseline period from tmin to 0.0 seconds by default, but it is wise to be explicit. That way you are less likely to end up with surprises along the way. None as the first element of the tuple refers to the start of the time window (-200 ms in this case). See mne.Epochs for more.

We also define rejection thresholds to get rid of noisy epochs. The rejection thresholds are defined as peak-to-peak values within the epoch time window. They are defined as T/m for gradiometers, T for magnetometers and V for EEG and EOG electrodes.

Note

In this tutorial, we don’t preprocess the data. This is not something you would normally do. See our Documentation on preprocessing for more.

baseline = (None, 0.0)
reject = {'mag': 4e-12, 'eog': 200e-6}
epochs = mne.Epochs(raw, events=events, event_id=event_id, tmin=tmin,
                    tmax=tmax, reject=reject, picks=picks)

Let’s plot the epochs to see the results. The number at the top refers to the id number. We can see that 128 good epochs out of total of 145 events got through the rejection process. Visual inspection also reveals that some epochs containing saccades or blinks got through. You can also reject epochs by hand by clicking on the epoch in the browser window. The selected epochs get rejected when you close the epochs browser. How you should reject the epochs and which thresholds to use is not a trivial question and this tutorial takes no stand on that matter.

To see all the interactive features of the epochs browser, click ‘Help’ in the lower left corner of the browser window.

epochs.plot(block=True)
../_images/sphx_glr_plot_epoching_and_averaging_004.png

To see why the epochs were rejected, we can plot the drop log.

epochs.plot_drop_log()
../_images/sphx_glr_plot_epoching_and_averaging_005.png

To get the evoked response you can simply do epochs.average(). It includes only the data channels by default. For the sake of example, we use picks to include the EOG channels as well. Notice that we cannot use the same picks as before as the indices are different. ‘Why are they different?’ you might ask. They’re different because picks is simply a list of channel indices and as the epochs were constructed, also a new info structure is created where the channel indices run from 0 to epochs.info['nchan']. See The Info data structure for more information.

picks = mne.pick_types(epochs.info, meg=True, eog=True)
evoked_left = epochs['Auditory/Left'].average(picks=picks)
evoked_right = epochs['Auditory/Right'].average(picks=picks)

Notice we have used forward slashes (‘/’) to separate the factors of the conditions of the experiment. We can use these ‘tags’ to select for example all left trials (both visual left and auditory right) …

epochs_left = epochs['Left']

# ... or to select a very specific subset. This is the same as above:
evoked_left = epochs['Left/Auditory'].average(picks=picks)

Finally, let’s plot the evoked responses.

evoked_left.plot()
evoked_right.plot()
  • ../_images/sphx_glr_plot_epoching_and_averaging_006.png
  • ../_images/sphx_glr_plot_epoching_and_averaging_007.png

Total running time of the script: ( 0 minutes 15.712 seconds)

Gallery generated by Sphinx-Gallery