Repairing artifacts with ICA

This tutorial covers the basics of independent components analysis (ICA) and shows how ICA can be used for artifact repair; an extended example illustrates repair of ocular and heartbeat artifacts. :depth: 2

We begin as always by importing the necessary Python modules and loading some example data <sample-dataset>. Because ICA can be computationally intense, we'll also crop the data to 60 seconds; and to save ourselves from repeatedly typing mne.preprocessing we'll directly import a few functions and classes from that submodule:

What is ICA? ^^^^^^^^^^^^

Independent components analysis (ICA) is a technique for estimating independent source signals from a set of recordings in which the source signals were mixed together in unknown ratios. A common example of this is the problem of blind source separation_: with 3 musical instruments playing in the same room, and 3 microphones recording the performance (each picking up all 3 instruments, but at varying levels), can you somehow "unmix" the signals recorded by the 3 microphones so that you end up with a separate "recording" isolating the sound of each instrument?

It is not hard to see how this analogy applies to EEG/MEG analysis: there are many "microphones" (sensor channels) simultaneously recording many "instruments" (blinks, heartbeats, activity in different areas of the brain, muscular activity from jaw clenching or swallowing, etc). As long as these various source signals are statistically independent_ and non-gaussian, it is usually possible to separate the sources using ICA, and then re-construct the sensor signals after excluding the sources that are unwanted.

ICA in MNE-Python

MNE-Python implements three different ICA algorithms: ``fastica`` (the
default), ``picard``, and ``infomax``. FastICA and Infomax are both in fairly
widespread use; Picard is a newer (2017) algorithm that is expected to
converge faster than FastICA and Infomax, and is more robust than other
algorithms in cases where the sources are not completely independent, which
typically happens with real EEG/MEG data. See [1]_ for more information.

The ICA interface in MNE-Python is similar to the interface in
`scikit-learn`_: some general parameters are specified when creating an
:class:`~mne.preprocessing.ICA` object, then the
:class:`~mne.preprocessing.ICA` object is fit to the data using its
:meth:`~mne.preprocessing.ICA.fit` method. The results of the fitting are
added to the :class:`~mne.preprocessing.ICA` object as attributes that end in
an underscore (``_``), such as ``ica.mixing_matrix_`` and
``ica.unmixing_matrix_``. After fitting, the ICA component(s) that you want
to remove must be chosen, and the ICA fit must then be applied to the
:class:`~mne.io.Raw` or :class:`~mne.Epochs` object using the
:class:`~mne.preprocessing.ICA` object's :meth:`~mne.preprocessing.ICA.apply`
method.

.. sidebar:: ICA and dimensionality reduction

    If you want to perform ICA with no dimensionality reduction (other than
    the number of Independent Components (ICs) given in ``n_components``, and
    any subsequent exclusion of ICs you specify in ``ICA.exclude``), pass
    ``max_pca_components=None`` and ``n_pca_components=None`` (these are the
    default values). If you want to reduce dimensionality, consider this
    example: if you have 300 sensor channels and you set
    ``max_pca_components=200``, ``n_components=50`` and
    ``n_pca_components=None``, then the PCA step yields 200 PCs, the first 50
    PCs are sent to the ICA algorithm (yielding 50 ICs), and during
    reconstruction :meth:`~mne.preprocessing.ICA.apply` will use the 50 ICs
    plus PCs number 51-200 (the full PCA residual). If instead you specify
    ``n_pca_components=120`` then :meth:`~mne.preprocessing.ICA.apply` will
    reconstruct using the 50 ICs plus the first 70 PCs in the PCA residual
    (numbers 51-120).

As is typically done with ICA, the data are first scaled to unit variance and
whitened using principal components analysis (PCA) before performing the ICA
decomposition. You can impose an optional dimensionality reduction at this
step by specifying ``max_pca_components``. From the retained Principal
Components (PCs), the first ``n_components`` are then passed to the ICA
algorithm (``n_components`` may be an integer number of components to use, or
a fraction of explained variance that used components should capture).

After visualizing the Independent Components (ICs) and excluding any that
capture artifacts you want to repair, the sensor signal can be reconstructed
using the :class:`~mne.preprocessing.ICA` object's
:meth:`~mne.preprocessing.ICA.apply` method. By default, signal
reconstruction uses all of the ICs (less any ICs listed in ``ICA.exclude``)
plus all of the PCs that were not included in the ICA decomposition (i.e.,
the "PCA residual"). If you want to reduce the number of components used at
the reconstruction stage, it is controlled by the ``n_pca_components``
parameter (which will in turn reduce the rank of your data; by default
``n_pca_components = max_pca_components`` resulting in no additional
dimensionality reduction). The fitting and reconstruction procedures and the
parameters that control dimensionality at various stages are summarized in
the diagram below:

.. graphviz:: ../../_static/diagrams/ica.dot
   :alt: Diagram of ICA procedure in MNE-Python
   :align: left

See the Notes section of the :class:`~mne.preprocessing.ICA` documentation
for further details. Next we'll walk through an extended example that
illustrates each of these steps in greater detail.


Example: EOG and ECG artifact repair
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

Visualizing the artifacts
~~~~~~~~

Let's begin by visualizing the artifacts that we want to repair. In this dataset they are big enough to see easily in the raw data:

We can get a summary of how the ocular artifact manifests across each channel type using :func:~mne.preprocessing.create_eog_epochs like we did in the tut-artifact-overview tutorial:

Now we'll do the same for the heartbeat artifacts, using :func:~mne.preprocessing.create_ecg_epochs:

Filtering to remove slow drifts ~~~~~~~

Before we run the ICA, an important step is filtering the data to remove low-frequency drifts, which can negatively affect the quality of the ICA fit. The slow drifts are problematic because they reduce the independence of the assumed-to-be-independent sources (e.g., during a slow upward drift, the neural, heartbeat, blink, and other muscular sources will all tend to have higher values), making it harder for the algorithm to find an accurate solution. A high-pass filter with 1 Hz cutoff frequency is recommended. However, because filtering is a linear operation, the ICA solution found from the filtered signal can be applied to the unfiltered signal (see [2]_ for more information), so we'll keep a copy of the unfiltered :class:~mne.io.Raw object around so we can apply the ICA solution to it later.

Fitting and plotting the ICA solution ~~~~~~~~~

.. sidebar:: Ignoring the time domain

The ICA algorithms implemented in MNE-Python find patterns across
channels, but ignore the time domain. This means you can compute ICA on
discontinuous :class:`~mne.Epochs` or :class:`~mne.Evoked` objects (not
just continuous :class:`~mne.io.Raw` objects), or only use every Nth
sample by passing the ``decim`` parameter to ``ICA.fit()``.

Now we're ready to set up and fit the ICA. Since we know (from observing our raw data) that the EOG and ECG artifacts are fairly strong, we would expect those artifacts to be captured in the first few dimensions of the PCA decomposition that happens before the ICA. Therefore, we probably don't need a huge number of components to do a good job of isolating our artifacts (though it is usually preferable to include more components for a more accurate solution). As a first guess, we'll run ICA with n_components=15 (use only the first 15 PCA components to compute the ICA decomposition) — a very small number given that our data has over 300 channels, but with the advantage that it will run quickly and we will able to tell easily whether it worked or not (because we already know what the EOG / ECG artifacts should look like).

ICA fitting is not deterministic (e.g., the components may get a sign flip on different runs, or may not always be returned in the same order), so we'll also specify a random seed_ so that we get identical results each time this tutorial is built by our web servers.

Some optional parameters that we could have passed to the :meth:~mne.preprocessing.ICA.fit method include decim (to use only every Nth sample in computing the ICs, which can yield a considerable speed-up) and reject (for providing a rejection dictionary for maximum acceptable peak-to-peak amplitudes for each channel type, just like we used when creating epoched data in the tut-overview tutorial).

Now we can examine the ICs to see what they captured. :meth:~mne.preprocessing.ICA.plot_sources will show the time series of the ICs. Note that in our call to :meth:~mne.preprocessing.ICA.plot_sources we can use the original, unfiltered :class:~mne.io.Raw object:

Here we can pretty clearly see that the first component (ICA000) captures the EOG signal quite well, and the second component (ICA001) looks a lot like a heartbeat <qrs_> (for more info on visually identifying Independent Components, this EEGLAB tutorial is a good resource). We can also visualize the scalp field distribution of each component using :meth:~mne.preprocessing.ICA.plot_components. These are interpolated based on the values in the ICA unmixing matrix:

In the plots above it's fairly obvious which ICs are capturing our EOG and ECG artifacts, but there are additional ways visualize them anyway just to be sure. First, we can plot an overlay of the original signal against the reconstructed signal with the artifactual ICs excluded, using :meth:~mne.preprocessing.ICA.plot_overlay:

We can also plot some diagnostics of each IC using :meth:~mne.preprocessing.ICA.plot_properties:

In the remaining sections, we'll look at different ways of choosing which ICs to exclude prior to reconstructing the sensor signals.

Selecting ICA components manually ~~~~~~~~~

Once we're certain which components we want to exclude, we can specify that manually by setting the ica.exclude attribute. Similar to marking bad channels, merely setting ica.exclude doesn't do anything immediately (it just adds the excluded ICs to a list that will get used later when it's needed). Once the exclusions have been set, ICA methods like :meth:~mne.preprocessing.ICA.plot_overlay will exclude those component(s) even if no exclude parameter is passed, and the list of excluded components will be preserved when using :meth:mne.preprocessing.ICA.save and :func:mne.preprocessing.read_ica.

Now that the exclusions have been set, we can reconstruct the sensor signals with artifacts removed using the :meth:~mne.preprocessing.ICA.apply method (remember, we're applying the ICA solution from the filtered data to the original unfiltered signal). Plotting the original raw data alongside the reconstructed data shows that the heartbeat and blink artifacts are repaired.

Using an EOG channel to select ICA components ~~~~~~~~~

It may have seemed easy to review the plots and manually select which ICs to exclude, but when processing dozens or hundreds of subjects this can become a tedious, rate-limiting step in the analysis pipeline. One alternative is to use dedicated EOG or ECG sensors as a "pattern" to check the ICs against, and automatically mark for exclusion any ICs that match the EOG/ECG pattern. Here we'll use :meth:~mne.preprocessing.ICA.find_bads_eog to automatically find the ICs that best match the EOG signal, then use :meth:~mne.preprocessing.ICA.plot_scores along with our other plotting functions to see which ICs it picked. We'll start by resetting ica.exclude back to an empty list:

Note that above we used :meth:~mne.preprocessing.ICA.plot_sources on both the original :class:~mne.io.Raw instance and also on an :class:~mne.Evoked instance of the extracted EOG artifacts. This can be another way to confirm that :meth:~mne.preprocessing.ICA.find_bads_eog has identified the correct components.

Using a simulated channel to select ICA components ~~~~~~~~~~

If you don't have an EOG channel, :meth:~mne.preprocessing.ICA.find_bads_eog has a ch_name parameter that you can use as a proxy for EOG. You can use a single channel, or create a bipolar reference from frontal EEG sensors and use that as virtual EOG channel. This carries a risk however: you must hope that the frontal EEG channels only reflect EOG and not brain dynamics in the prefrontal cortex (or you must not care about those prefrontal signals).

For ECG, it is easier: :meth:~mne.preprocessing.ICA.find_bads_ecg can use cross-channel averaging of magnetometer or gradiometer channels to construct a virtual ECG channel, so if you have MEG channels it is usually not necessary to pass a specific channel name. :meth:~mne.preprocessing.ICA.find_bads_ecg also has two options for its method parameter: 'ctps' (cross-trial phase statistics [3]_) and 'correlation' (Pearson correlation between data and ECG channel).

The last of these plots is especially useful: it shows us that the heartbeat artifact is coming through on two ICs, and we've only caught one of them. In fact, if we look closely at the output of :meth:~mne.preprocessing.ICA.plot_sources (online, you can right-click → "view image" to zoom in), it looks like ICA014 has a weak periodic component that is in-phase with ICA001. It might be worthwhile to re-run the ICA with more components to see if that second heartbeat artifact resolves out a little better:

Much better! Now we've captured both ICs that are reflecting the heartbeat artifact (and as a result, we got two diagnostic plots: one for each IC that reflects the heartbeat). This demonstrates the value of checking the results of automated approaches like :meth:~mne.preprocessing.ICA.find_bads_ecg before accepting them.

Selecting ICA components using template matching ~~~~~~~~~~~~

When dealing with multiple subjects, it is also possible to manually select an IC for exclusion on one subject, and then use that component as a template for selecting which ICs to exclude from other subjects' data, using :func:mne.preprocessing.corrmap [4]_. The idea behind :func:~mne.preprocessing.corrmap is that the artifact patterns are similar enough across subjects that corresponding ICs can be identified by correlating the ICs from each ICA solution with a common template, and picking the ICs with the highest correlation strength. :func:~mne.preprocessing.corrmap takes a list of ICA solutions, and a template parameter that specifies which ICA object and which component within it to use as a template.

Since our sample dataset only contains data from one subject, we'll use a different dataset with multiple subjects: the EEGBCI dataset [5] [6]. The dataset has 109 subjects, we'll just download one run (a left/right hand movement task) from each of the first 4 subjects:

Now let's run :func:~mne.preprocessing.corrmap:

The first figure shows the template map, while the second figure shows all the maps that were considered a "match" for the template (including the template itself). There were only three matches from the four subjects; notice the output message No maps selected for subject(s) 1, consider a more liberal threshold. By default the threshold is set automatically by trying several values; here it may have chosen a threshold that is too high. Let's take a look at the ICA sources for each subject:

Notice that subject 1 does seem to have an IC that looks like it reflects blink artifacts (component ICA000). Notice also that subject 3 appears to have two components that are reflecting ocular artifacts (ICA000 and ICA002), but only one was caught by :func:~mne.preprocessing.corrmap. Let's try setting the threshold manually:

Now we get the message At least 1 IC detected for each subject (which is good). At this point we'll re-run :func:~mne.preprocessing.corrmap with parameters label=blink, show=False to label the ICs from each subject that capture the blink artifacts (without plotting them again).

Notice that the first subject has 3 different labels for the IC at index 0: "eog/0/Fpz", "eog", and "blink". The first two were added by :meth:~mne.preprocessing.ICA.find_bads_eog; the "blink" label was added by the last call to :func:~mne.preprocessing.corrmap. Notice also that each subject has at least one IC index labelled "blink", and subject 3 has two components (0 and 2) labelled "blink" (consistent with the plot of IC sources above). The labels_ attribute of :class:~mne.preprocessing.ICA objects can also be manually edited to annotate the ICs with custom labels. They also come in handy when plotting:

As a final note, it is possible to extract ICs numerically using the :meth:~mne.preprocessing.ICA.get_components method of :class:~mne.preprocessing.ICA objects. This will return a :class:NumPy array <numpy.ndarray> that can be passed to :func:~mne.preprocessing.corrmap instead of the :class:tuple of (subject_index, component_index) we passed before, and will yield the same result:

An advantage of using this numerical representation of an IC to capture a particular artifact pattern is that it can be saved and used as a template for future template-matching tasks using :func:~mne.preprocessing.corrmap without having to load or recompute the ICA solution that yielded the template originally. Put another way, when the template is a NumPy array, the :class:~mne.preprocessing.ICA object containing the template does not need to be in the list of ICAs provided to :func:~mne.preprocessing.corrmap.

References ^^^^^^^^^^

.. [1] Ablin P, Cardoso J, Gramfort A (2018). Faster Independent Component Analysis by Preconditioning With Hessian Approximations. IEEE Transactions on Signal Processing 66:4040–4049. https://doi.org/10.1109/TSP.2018.2844203

.. [2] Winkler I, Debener S, Müller K-R, Tangermann M (2015). On the influence of high-pass filtering on ICA-based artifact reduction in EEG-ERP. Proceedings of EMBC-2015, 4101–4105. https://doi.org/10.1109/EMBC.2015.7319296

.. [3] Dammers J, Schiek M, Boers F, Silex C, Zvyagintsev M, Pietrzyk U, Mathiak K (2008). Integration of amplitude and phase statistics for complete artifact removal in independent components of neuromagnetic recordings. IEEE Transactions on Biomedical Engineering 55(10):2353–2362. https://doi.org/10.1109/TBME.2008.926677

.. [4] Viola FC, Thorne J, Edmonds B, Schneider T, Eichele T, Debener S (2009). Semi-automatic identification of independent components representing EEG artifact. Clinical Neurophysiology 120(5):868–877. https://doi.org/10.1016/j.clinph.2009.01.015

.. [5] Schalk G, McFarland DJ, Hinterberger T, Birbaumer N, Wolpaw JR (2004). BCI2000: A General-Purpose Brain-Computer Interface (BCI) System. IEEE Transactions on Biomedical Engineering 51(6):1034-1043. https://doi.org/10.1109/TBME.2004.827072

.. [6] Goldberger AL, Amaral LAN, Glass L, Hausdorff JM, Ivanov PCh, Mark RG, Mietus JE, Moody GB, Peng C-K, Stanley HE (2000). PhysioBank, PhysioToolkit, and PhysioNet: Components of a New Research Resource for Complex Physiologic Signals. Circulation 101(23):e215-e220. https://doi.org/10.1161/01.CIR.101.23.e215

.. LINKS

https://en.wikipedia.org/wiki/Signal_separation https://en.wikipedia.org/wiki/Independence_(probability_theory)