Add regression-based approach to removing EOG artifacts (#11046)

Co-authored-by: Eric Larson <larson.eric.d@gmail.com>
Co-authored-by: Daniel McCloy <dan@mccloy.info>
Co-authored-by: Alexandre Gramfort <alexandre.gramfort@m4x.org>

doc/changes/latest.inc

@@ -85,3 +85,5 @@ API changes
 - New classes :class:`~mne.time_frequency.Spectrum` and :class:`~mne.time_frequency.EpochsSpectrum`, created via new methods :meth:`Raw.compute_psd()<mne.io.Raw.compute_psd>`, :meth:`Epochs.compute_psd()<mne.Epochs.compute_psd>`, and :meth:`Evoked.compute_psd()<mne.Evoked.compute_psd>` (:gh:`10184` by `Daniel McCloy`_)
 - The ``mne.epochs.add_channels_epochs`` function has been deprecated in favor of :meth:`epochs.add_channels <mne.Epochs.add_channels>` (:gh:`11180` by `Eric Larson`_)
 - The PSD functions that operate on Raw/Epochs/Evoked instances (``mne.time_frequency.psd_welch`` and ``mne.time_frequency.psd_multitaper``) are deprecated; for equivalent functionality create :class:`~mne.time_frequency.Spectrum` or :class:`~mne.time_frequency.EpochsSpectrum` objects instead and then run ``spectrum.get_data(return_freqs=True)`` (:gh:`10184` by `Daniel McCloy`_)
+- Added new class  :class:`mne.preprocessing.EOGRegression` to allow more flexibility when using regression to reduce EOG artifacts (:gh:`11046` by `Marijn van Vliet`_)
+- New parameter ``exclude`` added to :func:`mne.preprocessing.regress_artifact` to prevent regression from being applied to certain channels (:gh:`11046` by `Marijn van Vliet`_)

doc/conf.py

@@ -237,6 +237,7 @@
     'Transform': 'mne.transforms.Transform',
     'Coregistration': 'mne.coreg.Coregistration',
     'Figure3D': 'mne.viz.Figure3D',
+    'EOGRegression': 'mne.preprocessing.EOGRegression',
     'Spectrum': 'mne.time_frequency.Spectrum',
     'EpochsSpectrum': 'mne.time_frequency.EpochsSpectrum',
     # dipy

doc/preprocessing.rst

@@ -70,6 +70,7 @@ Projections:
 
    ICA
    Xdawn
+   EOGRegression
    annotate_amplitude
    annotate_break
    annotate_movement
@@ -99,6 +100,7 @@ Projections:
    oversampled_temporal_projection
    peak_finder
    read_ica
+   read_eog_regression
    realign_raw
    regress_artifact
    corrmap

doc/references.bib

@@ -2354,3 +2354,14 @@ @article{StenroosHauk2013
 	pages = {265--272},
 	file = {Full Text:/Users/larsoner/Zotero/storage/VJAIPDFL/Stenroos and Hauk - 2013 - Minimum-norm cortical source estimation in layered.pdf:application/pdf;ScienceDirect Snapshot:/Users/larsoner/Zotero/storage/FR9YZCJZ/S1053811913004333.html:text/html},
 }
+
+@article{CroftBarry2000,
+  title = {Removal of ocular artifact from the {EEG}: a review},
+  author = {Croft, R. J. and Barry, R. J.},
+  year = {2000},
+  journal = {Clinical Neurophysiology},
+  volume = {30},
+  number = {1},
+  pages = {5--19},
+  doi = {10.1016/S0987-7053(00)00055-1}
+}

doc/visualization.rst

@@ -59,6 +59,7 @@ Visualization
    plot_projs_joint
    plot_raw
    plot_raw_psd
+   plot_regression_weights
    plot_sensors
    plot_snr_estimate
    plot_source_estimates

examples/preprocessing/eog_regression.py

@@ -0,0 +1,78 @@
+# -*- coding: utf-8 -*-
+"""
+=======================================
+Reduce EOG artifacts through regression
+=======================================
+
+Reduce artifacts by regressing the EOG channels onto the rest of the channels
+and then subtracting the EOG signal.
+
+This is a quick example to show the most basic application of the technique.
+See the :ref:`tutorial <tut-artifact-regression>` for a more thorough
+explanation that demonstrates more advanced approaches.
+"""
+
+# Author: Marijn van Vliet <w.m.vanvliet@gmail.com>
+#
+# License: BSD (3-clause)
+
+# %%
+# Import packages and load data
+# -----------------------------
+#
+# We begin as always by importing the necessary Python modules and loading some
+# data, in this case the :ref:`MNE sample dataset <sample-dataset>`.
+
+import mne
+from mne.datasets import sample
+from mne.preprocessing import EOGRegression
+from matplotlib import pyplot as plt
+
+print(__doc__)
+
+data_path = sample.data_path()
+raw_fname = data_path / 'MEG' / 'sample' / 'sample_audvis_filt-0-40_raw.fif'
+
+# Read raw data
+raw = mne.io.read_raw_fif(raw_fname, preload=True)
+events = mne.find_events(raw, 'STI 014')
+
+# Highpass filter to eliminate slow drifts
+raw.filter(0.3, None, picks='all')
+
+# %%
+# Perform regression and remove EOG
+# ---------------------------------
+
+# Fit the regression model
+weights = EOGRegression().fit(raw)
+raw_clean = weights.apply(raw, copy=True)
+
+# Show the filter weights in a topomap
+weights.plot()
+
+# %%
+# Before/after comparison
+# -----------------------
+# Let's compare the signal before and after cleaning with EOG regression. This
+# is best visualized by extracting epochs and plotting the evoked potential.
+
+tmin, tmax = -0.1, 0.5
+event_id = {'visual/left': 3, 'visual/right': 4}
+evoked_before = mne.Epochs(raw, events, event_id, tmin, tmax,
+                           baseline=(tmin, 0)).average()
+evoked_after = mne.Epochs(raw_clean, events, event_id, tmin, tmax,
+                          baseline=(tmin, 0)).average()
+
+# Create epochs after EOG correction
+epochs_after = mne.Epochs(raw_clean, events, event_id, tmin, tmax,
+                          baseline=(tmin, 0))
+evoked_after = epochs_after.average()
+
+fig, ax = plt.subplots(nrows=3, ncols=2, figsize=(10, 7),
+                       sharex=True, sharey='row')
+evoked_before.plot(axes=ax[:, 0], spatial_colors=True)
+evoked_after.plot(axes=ax[:, 1], spatial_colors=True)
+fig.subplots_adjust(top=0.905, bottom=0.09, left=0.08, right=0.975,
+                    hspace=0.325, wspace=0.145)
+fig.suptitle('Before --> After')

mne/preprocessing/__init__.py

@@ -26,7 +26,7 @@
 from . import nirs
 from .artifact_detection import (annotate_movement, compute_average_dev_head_t,
                                  annotate_muscle_zscore, annotate_break)
-from ._regress import regress_artifact
+from ._regress import regress_artifact, EOGRegression, read_eog_regression
 from ._fine_cal import (compute_fine_calibration, read_fine_calibration,
                         write_fine_calibration)
 from .annotate_nan import annotate_nan