DOC: EEG eye-tracking alignment tutorial

doc/changes/latest.inc

@@ -30,6 +30,7 @@ Enhancements
 - Added :class:`mne.preprocessing.eyetracking.Calibration` to store eye-tracking calibration info, and :func:`mne.preprocessing.eyetracking.read_eyelink_calibration` to read calibration data from EyeLink systems (:gh:`11719` by `Scott Huberty`_)
 - Ocular :class:`mne.Annotations` read in by :func:`mne.io.read_raw_eyelink` are now channel aware. This means if the left eye blinked, the associated annotation will store this in the ``'ch_names'`` key. (:gh:`11746` by `Scott Huberty`_)
 - Added :func:`mne.preprocessing.eyetracking.interpolate_blinks` to linear interpolate eyetrack signals during blink periods. (:gh:`11740` by `Scott Huberty`_)
+- Added a section for combining eye-tracking and EEG data to the preprocessing tutorial "working with eye tracker data in MNE-Python" (:gh:`11770` by `Scott Huberty`_)
 
 Bugs
 ~~~~

doc/overview/datasets_index.rst

@@ -475,14 +475,17 @@ standard.
 
     * :ref:`tut-ssvep`
 
+.. _eyelink-dataset:
+
 EYELINK
 =======
 :func:`mne.datasets.eyelink.data_path`
 
-A small example dataset in SR research's proprietary .asc format.
-1 participant fixated on the screen while short light flashes appeared.
-Monocular recording of gaze position and pupil size, 1000 Hz sampling 
-frequency.
+A small example dataset from a pupillary light reflex experiment. Both EEG (EGI) and
+eye-tracking (SR Research EyeLink; ASCII format) data were recorded and stored in
+separate files. 1 participant fixated on the screen while short light flashes appeared.
+Event onsets were recorded by a photodiode attached to the screen and were
+sent to both the EEG and eye-tracking systems.
 
 .. topic:: Examples
 

mne/datasets/config.py

@@ -345,9 +345,9 @@
 
 # eyelink dataset
 MNE_DATASETS["eyelink"] = dict(
-    archive_name="eyelink_example_data.zip",
-    hash="md5:081950c05f35267458d9c751e178f161",
-    url=("https://osf.io/r5ndq/download?version=1"),
+    archive_name="eeg-eyetrack_data.zip",
+    hash="md5:c4fc788fe01737e08e9086c90cab642d",
+    url=("https://osf.io/63fjm/download?version=1"),
     folder_name="eyelink-example-data",
     config_key="MNE_DATASETS_EYELINK_PATH",
 )

tutorials/preprocessing/90_eyetracking_data.py

@@ -6,67 +6,74 @@
 Working with eye tracker data in MNE-Python
 ===========================================
 
-In this tutorial we will load some eye tracker data and plot the average
-pupil response to light flashes (i.e. the pupillary light reflex).
+In this tutorial we will explore simultaneously recorded eye-tracking and EEG data from
+a pupillary light reflex task. We will combine the eye-tracking and EEG data, and plot
+the ERP and pupil response to the light flashes (i.e. the pupillary light reflex).
 
 """  # noqa: E501
-# Authors: Dominik Welke <[email protected]>
-#          Scott Huberty <[email protected]>
+# Authors: Scott Huberty <[email protected]>
+#          Dominik Welke <[email protected]>
+#
 #
 # License: BSD-3-Clause
 
 # %%
 # Data loading
 # ------------
 #
-# First we will load an eye tracker recording from SR research's proprietary
-# ``'.asc'`` file format.
+# As usual we start by importing the modules we need and loading some
+# :ref:`example data <eyelink-dataset>`: eye-tracking data recorded from SR research's
+# ``'.asc'`` file format, and EEG data recorded from EGI's ``'.mff'`` file format.
 #
-# The info structure tells us we loaded a monocular recording with 2
-# ``'eyegaze'``, channels (X/Y), 1 ``'pupil'`` channel, and 1 ``'stim'``
-# channel.
+# The info structure of the eye-tracking data tells us we loaded a monocular recording
+# with 2 ``'eyegaze'``, channels (X/Y), 1 ``'pupil'`` channel, 1 ``'stim'`` channel, and
+# 3 channels for the head distance and position (since this data was collected using
+# EyeLink's Remote mode).
 
 import mne
 from mne.datasets.eyelink import data_path
 from mne.preprocessing.eyetracking import read_eyelink_calibration
 
-eyelink_fname = data_path() / "mono_multi-block_multi-DINS.asc"
+et_fpath = data_path() / "sub-01_task-plr_eyetrack.asc"
+eeg_fpath = data_path() / "sub-01_task-plr_eeg.mff"
+
+raw_et = mne.io.read_raw_eyelink(et_fpath, create_annotations=["blinks"])
+raw_eeg = mne.io.read_raw_egi(eeg_fpath, preload=True, verbose="warning")
+raw_eeg.filter(1, 30, verbose="warning")
+print(f"RawEyelink info: {raw_et.info}")
 
-raw = mne.io.read_raw_eyelink(eyelink_fname, create_annotations=["blinks", "messages"])
-raw.crop(tmin=0, tmax=130)  # for this demonstration, let's take a subset of the data
+# %%
+# .. seealso:: :ref:`tut-importing-eyetracking-data`
 
 # %%
 # Ocular annotations
 # ------------------
-# By default, Eyelink files will output events for ocular events (blinks,
+# By default, EyeLink files will output events for ocular events (blinks,
 # saccades, fixations), and experiment messages. MNE will store these events
 # as `mne.Annotations`. Ocular annotations contain channel information, in the
-# ``'ch_names'``` key. This means that we can see which eye an ocular event occurred in:
+# ``'ch_names'`` key. This means that we can see which eye an ocular event occurred in,
+# which can be useful for binocular recordings:
 
-print(raw.annotations[0])  # a blink in the right eye
+print(raw_et.annotations[0])  # a blink in the right eye
 
 # %%
 # If we are only interested in certain event types from
-# the Eyelink file, we can select for these using the ``'create_annotations'``
-# argument of `mne.io.read_raw_eyelink`. above, we only created annotations
-# for blinks, and experiment messages.
-#
-# Note that ``'blink'`` annotations are read in as ``'BAD_blink'``, and MNE will treat
-# these as bad segments of data. This means that blink periods will be dropped during
-# epoching by default.
+# the EyeLink file, we can select for these using the ``'create_annotations'``
+# argument of `~mne.io.read_raw_eyelink`. above, we only created annotations
+# for blinks, which are read in as ``'BAD_blink'`` so that MNE will treat
+# these as bad segments of data.
 
 # %%
 # Checking the calibration
 # ------------------------
 #
-# We can also load the calibrations from the recording and visualize them.
-# Checking the quality of the calibration is a useful first step in assessing
-# the quality of the eye tracking data. Note that
+# We can also load and visualize eye-tracking calibrations, which is a useful first step
+# in assessing the quality of the eye-tracking data.
 # :func:`~mne.preprocessing.eyetracking.read_eyelink_calibration`
 # will return a list of :class:`~mne.preprocessing.eyetracking.Calibration` instances,
 # one for each calibration. We can index that list to access a specific calibration.
 
-cals = read_eyelink_calibration(eyelink_fname)
+cals = read_eyelink_calibration(et_fpath)
 print(f"number of calibrations: {len(cals)}")
 first_cal = cals[0]  # let's access the first (and only in this case) calibration
 print(first_cal)
@@ -101,123 +108,107 @@
 first_cal.plot(show_offsets=True)
 
 # %%
-# Get stimulus events from DIN channel
-# ------------------------------------
+# Extract common stimulus events from the data
+# --------------------------------------------
 #
-# Eyelink eye trackers have a DIN port that can be used to feed in stimulus
-# or response timings. :func:`mne.io.read_raw_eyelink` loads this data as a
-# ``'stim'`` channel. Alternatively, the onset of stimulus events could be sent
-# to the eyetracker as ``messages`` - these can be read in as
-# `mne.Annotations`.
+# In this experiment, a photodiode attached to the display screen was connected to both
+# the EEG and eye-tracking systems. The photodiode was triggered by the the light flash
+# stimuli, causing a signal to be sent to both systems simultaneously, signifying the
+# onset of the flash. The photodiode signal was recorded as a digital input channel in
+# the EEG and eye-tracking data. MNE loads these data as ``'stim'`` channels.
 #
-# In the example data, the DIN channel contains the onset of light flashes on
-# the screen. We now extract these events to visualize the pupil response. We will use
-# these later in this tutorial.
-
-events = mne.find_events(raw, shortest_event=1, min_duration=0.02, uint_cast=True)
-event_dict = {"flash": 3}
+# We'll extract the flash event onsets from both the EEG and eye-tracking data, as they
+# are necessary for aligning the EEG and eye-tracking data.
 
+et_events = mne.find_events(raw_et, min_duration=0.01, shortest_event=1, uint_cast=True)
+event_dict = {"Flash": 2}
+eeg_events = mne.find_events(raw_eeg, stim_channel="DIN3")
 
 # %%
-# Plot raw data
-# -------------
+# Plot the raw eye-tracking data
+# ------------------------------
 #
-# As the following plot shows, we now have a raw object with the eye tracker
-# data, eyeblink annotations and stimulus events (from the DIN channel).
-#
-# The plot also shows us that there is some noise in the data (not always
-# categorized as blinks). Also, notice that we have passed a custom `dict` into
-# the scalings argument of ``raw.plot``. This is necessary to make the eyegaze
-# channel traces legible when plotting, since the file contains pixel position
-# data (as opposed to eye angles, which are reported in radians). We also could
-# have simply passed ``scalings='auto'``.
-
-raw.plot(
-    events=events,
-    event_id={"Flash": 3},
+# Let's plot the raw eye-tracking data. We'll pass a custom `dict` into
+# the scalings argument to make the eyegaze channel traces legible when plotting,
+# since this file contains pixel position data (as opposed to eye angles,
+# which are reported in radians). We also could have simply passed ``scalings='auto'``.
+
+raw_et.plot(
+    events=et_events,
+    event_id=event_dict,
     event_color="g",
-    start=25,
-    duration=45,
+    duration=15,
     scalings=dict(eyegaze=1e3),
 )
 
 # %%
-# handling blink artifacts
+# Handling blink artifacts
 # ------------------------
-# We also notice that, naturally, there are blinks in our data, and these blink periods
+#
+# Naturally, there are blinks in our data, and these blink periods
 # occur within ``"BAD_blink"``  annotations. During blink periods, ``"eyegaze"``
 # coordinates are not reported, and ``"pupil"`` size data are ``0``. We don't want these
-# blink artifacts biasing our analysis, so we have two options: We can either drop the
-# blink periods from our data during epoching, or we can interpolate the missing data
-# during the blink periods. For this tutorial, let's interpolate the missing samples:
+# blink artifacts biasing our analysis, so we have two options: Drop the blink periods
+# from our data during epoching, or interpolate the missing data during the blink
+# periods. For this tutorial, let's interpolate the blink samples:
 
-mne.preprocessing.eyetracking.interpolate_blinks(raw, buffer=0.05)
-# Let's plot our data again to see the result of the interpolation:
-raw.pick(["pupil_right"])  # Let's pick just the pupil channel
-raw.plot(events=events, event_id={"Flash": 3}, event_color="g")
+mne.preprocessing.eyetracking.interpolate_blinks(raw_et, buffer=(0.05, 0.2))
 
 # %%
-# :func:`~mne.preprocessing.eyetracking.interpolate_blinks` performs a simple linear
-# interpolation of the pupil size data during blink periods. the ``buffer`` keyword
-# argument specifies the amount of time (in seconds) before and after the blinks to
-# include in the interpolation. This is helpful because the ``blink`` annotations
-# do not always capture the entire blink in the signal. We specified a value of ``.05``
-# seconds (50 ms), which is slightly more than the default value of ``.025``.
-
-# %%
-# Dealing with high frequency noise
-# ---------------------------------
-# From the plot above, we notice that there is some high frequency noise in the pupil
-# signal. We can remove this noise by low-pass filtering the data:
+# Aligning the eye-tracking data with EEG data
+# --------------------------------------------
+#
+# Eye-tracking and EEG data were recorded simultaneously, but on different
+# systems, so we'll need to align the data before we can analyze them together. We can
+# do this using the :func:`~mne.preprocessing.realign_raw` function, which will align
+# the data based on the timing of the shared events that are present in both
+# :class:`~mne.io.Raw` objects. We'll use the shared photodiode events we extracted
+# above, but first we need to convert the event onsets from samples to seconds. Once the
+# data have been aligned, we'll add the EEG channels to the eye-tracking raw object.
+
+# Convert event onsets from samples to seconds
+et_din_times = et_events[:, 0] / raw_et.info["sfreq"]
+eeg_din_times = eeg_events[:, 0] / raw_eeg.info["sfreq"]
+# Align the data
+mne.preprocessing.realign_raw(
+    raw_et, raw_eeg, et_din_times, eeg_din_times, verbose="error"
+)
+# Add EEG channels to the eye-tracking raw object
+raw_et.add_channels([raw_eeg], force_update_info=True)
+
+# Define a few channel groups of interest and plot the data
+frontal = ["E19", "E11", "E4", "E12", "E5"]
+occipital = ["E61", "E62", "E78", "E67", "E72", "E77"]
+pupil = ["pupil_right"]
+picks_idx = mne.pick_channels(
+    raw_et.ch_names, frontal + occipital + pupil, ordered=True
+)
+raw_et.plot(events=et_events, event_id=event_dict, event_color="g", order=picks_idx)
 
-# Apply a low pass filter to the pupil channel
-raw.filter(l_freq=None, h_freq=40, picks=["pupil_right"])
 
 # %%
-# Rejecting bad spans of data
-# ^^^^^^^^^^^^^^^^^^^^^^^^^^^
-# Even after filtering the pupil data and interpolating the blink periods, we still see
-# some artifacts in the data (the large spikes) that we don't want to include in our
-# analysis. Let's epoch our data and then reject any epochs that might contain these
-# artifacts. We'll use :class:`mne.Epochs` to epoch our data, and pass in the
-# ``events`` array and ``event_dict`` that we created earlier. We'll also pass in the
-# ``reject`` keyword argument to reject any epochs that contain data that exceeds a
-# peak-to-peak signal amplitude threshold of ``1500`` in the ``"pupil"`` channel.
-# Note that this threshold is arbitrary, and should be adjusted based on the data.
-# We chose 1500 because eyelink reports pupil size in arbitrary units (AU), which
-# typically ranges from 800 to 3000 units. Our epochs already contains large
-# signal fluctuations due to the pupil response, so a threshold of 1500 is conservative
-# enough to reject epochs only with large artifacts.
+# Epoching the data
+# -----------------
 
 epochs = mne.Epochs(
-    raw,
-    events,
-    tmin=-0.3,
-    tmax=5,
+    raw_et,
+    events=et_events,
     event_id=event_dict,
+    tmin=-0.3,
+    tmax=3,
     preload=True,
-    reject=dict(pupil=1500),
 )
-epochs.plot(events=events, event_id=event_dict)
+epochs[:8].plot(events=et_events, event_id=event_dict, order=picks_idx)
 
 # %%
 # We can clearly see the prominent decrease in pupil size following the
 # stimulation.
 
 # %%
-# Plot average pupil response
-# ---------------------------
+# Plot the evoked response
+# ------------------------
 #
-# Finally, let's plot the evoked response to the light flashes to get a sense of the
-# average pupillary light response.
+# Finally, let's plot the evoked responses to the light flashes to get a sense of the
+# average pupillary light response, and the ERP in the EEG data.
 
-epochs.average().plot()
-
-# %%
-# Again, it is important to note that pupil size data are reported by Eyelink (and
-# stored internally by MNE) as arbitrary units (AU). While it often can be
-# preferable to convert pupil size data to millimeters, this requires
-# information that is not present in the file. MNE does not currently
-# provide methods to convert pupil size data.
-# See :ref:`tut-importing-eyetracking-data` for more information on pupil size
-# data.
+epochs.average().plot(picks=occipital + pupil)