Whitening evoked data with a noise covariance

Evoked data are loaded and then whitened using a given noise covariance matrix. It’s an excellent quality check to see if baseline signals match the assumption of Gaussian white noise from which we expect values around 0 with less than 2 standard deviations. Covariance estimation and diagnostic plots are based on [1].

References

[1]Engemann D. and Gramfort A. (2015) Automated model selection in covariance estimation and spatial whitening of MEG and EEG signals, vol. 108, 328-342, NeuroImage.
# Authors: Alexandre Gramfort <alexandre.gramfort@telecom-paristech.fr>
#          Denis A. Engemann <denis.engemann@gmail.com>
#
# License: BSD (3-clause)

import mne

from mne import io
from mne.datasets import sample
from mne.cov import compute_covariance

print(__doc__)

Set parameters

data_path = sample.data_path()
raw_fname = data_path + '/MEG/sample/sample_audvis_filt-0-40_raw.fif'
event_fname = data_path + '/MEG/sample/sample_audvis_filt-0-40_raw-eve.fif'

raw = io.read_raw_fif(raw_fname, preload=True)
raw.filter(1, 40, n_jobs=1)
raw.info['bads'] += ['MEG 2443']  # bads + 1 more
events = mne.read_events(event_fname)

# let's look at rare events, button presses
event_id, tmin, tmax = 2, -0.2, 0.5
picks = mne.pick_types(raw.info, meg=True, eeg=True, eog=True, exclude='bads')
reject = dict(mag=4e-12, grad=4000e-13, eeg=80e-6)

epochs = mne.Epochs(raw, events, event_id, tmin, tmax, picks=picks,
                    baseline=None, reject=reject, preload=True)

# Uncomment next line to use fewer samples and study regularization effects
# epochs = epochs[:20]  # For your data, use as many samples as you can!

Out:

Opening raw data file /home/ubuntu/mne_data/MNE-sample-data/MEG/sample/sample_audvis_filt-0-40_raw.fif...
    Read a total of 4 projection items:
        PCA-v1 (1 x 102)  idle
        PCA-v2 (1 x 102)  idle
        PCA-v3 (1 x 102)  idle
        Average EEG reference (1 x 60)  idle
    Range : 6450 ... 48149 =     42.956 ...   320.665 secs
Ready.
Current compensation grade : 0
Reading 0 ... 41699  =      0.000 ...   277.709 secs...
Setting up band-pass filter from 1 - 40 Hz
l_trans_bandwidth chosen to be 1.0 Hz
h_trans_bandwidth chosen to be 10.0 Hz
Filter length of 991 samples (6.600 sec) selected
73 matching events found
Created an SSP operator (subspace dimension = 4)
4 projection items activated
Loading data for 73 events and 106 original time points ...
    Rejecting  epoch based on EEG : [u'EEG 001', u'EEG 002', u'EEG 003', u'EEG 007']
    Rejecting  epoch based on EEG : [u'EEG 001', u'EEG 002', u'EEG 003', u'EEG 004', u'EEG 005', u'EEG 006', u'EEG 007']
    Rejecting  epoch based on EEG : [u'EEG 003', u'EEG 007']
    Rejecting  epoch based on EEG : [u'EEG 001', u'EEG 002', u'EEG 003', u'EEG 006', u'EEG 007']
    Rejecting  epoch based on MAG : [u'MEG 1711']
    Rejecting  epoch based on EEG : [u'EEG 007']
    Rejecting  epoch based on EEG : [u'EEG 008', u'EEG 009']
    Rejecting  epoch based on EEG : [u'EEG 001', u'EEG 002', u'EEG 003', u'EEG 007']
    Rejecting  epoch based on EEG : [u'EEG 001', u'EEG 002', u'EEG 003', u'EEG 006', u'EEG 007']
    Rejecting  epoch based on EEG : [u'EEG 001', u'EEG 002', u'EEG 003', u'EEG 007']
    Rejecting  epoch based on EEG : [u'EEG 001', u'EEG 002', u'EEG 003', u'EEG 007']
    Rejecting  epoch based on EEG : [u'EEG 001', u'EEG 007']
12 bad epochs dropped

Compute covariance using automated regularization

noise_covs = compute_covariance(epochs, tmin=None, tmax=0, method='auto',
                                return_estimators=True, verbose=True, n_jobs=1,
                                projs=None)

# With "return_estimator=True" all estimated covariances sorted
# by log-likelihood are returned.

print('Covariance estimates sorted from best to worst')
for c in noise_covs:
    print("%s : %s" % (c['method'], c['loglik']))

Out:

Estimating covariance using SHRUNK
Done.
Estimating covariance using DIAGONAL_FIXED
    EEG regularization : None
    MAG regularization : 0.01
    GRAD regularization : 0.01
Done.
Estimating covariance using EMPIRICAL
Done.
Estimating covariance using FACTOR_ANALYSIS
... rank: 5 - loglik: -1826.779
... rank: 10 - loglik: -1776.750
... rank: 15 - loglik: -1733.126
... rank: 20 - loglik: -1709.022
... rank: 25 - loglik: -1691.479
... rank: 30 - loglik: -1679.398
... rank: 35 - loglik: -1672.391
... rank: 40 - loglik: -1669.625
... rank: 45 - loglik: -1666.613
... rank: 50 - loglik: -1666.548
... rank: 55 - loglik: -1667.393
... rank: 60 - loglik: -1667.699
... rank: 65 - loglik: -1668.827
early stopping parameter search.
... best model at rank = 50
Done.
Using cross-validation to select the best estimator.
    EEG regularization : None
    MAG regularization : 0.01
    GRAD regularization : 0.01
    EEG regularization : None
    MAG regularization : 0.01
    GRAD regularization : 0.01
    EEG regularization : None
    MAG regularization : 0.01
    GRAD regularization : 0.01
Number of samples used : 1891
[done]
Number of samples used : 1891
[done]
Number of samples used : 1891
[done]
Number of samples used : 1891
[done]
log-likelihood on unseen data (descending order):
   shrunk: -1644.818
   factor_analysis: -1666.548
   diagonal_fixed: -1743.812
   empirical: -1871.579
Covariance estimates sorted from best to worst
shrunk : -1644.81829922
factor_analysis : -1666.5484432
diagonal_fixed : -1743.81197537
empirical : -1871.57922919

Show whitening

evoked = epochs.average()

evoked.plot()  # plot evoked response

# plot the whitened evoked data for to see if baseline signals match the
# assumption of Gaussian white noise from which we expect values around
# 0 with less than 2 standard deviations. For the Global field power we expect
# a value of 1.

evoked.plot_white(noise_covs)
  • ../../_images/sphx_glr_plot_evoked_whitening_001.png
  • ../../_images/sphx_glr_plot_evoked_whitening_002.png

Out:

estimated rank (eeg): 59
estimated rank (grad): 203
estimated rank (mag): 102
estimated rank (mag + grad): 305
estimated rank (eeg): 59
estimated rank (eeg): 59
estimated rank (grad): 203
estimated rank (mag): 102
estimated rank (mag + grad): 305
estimated rank (eeg): 59
estimated rank (eeg): 58
estimated rank (grad): 203
estimated rank (mag): 102
estimated rank (mag + grad): 305
estimated rank (eeg): 58
estimated rank (eeg): 58
estimated rank (grad): 203
estimated rank (mag): 99
estimated rank (mag + grad): 302
estimated rank (eeg): 58
    Created an SSP operator (subspace dimension = 4)
Setting small MEG eigenvalues to zero.
Not doing PCA for MEG.
Setting small EEG eigenvalues to zero.
Not doing PCA for EEG.
    Created an SSP operator (subspace dimension = 4)
Setting small MEG eigenvalues to zero.
Not doing PCA for MEG.
Setting small EEG eigenvalues to zero.
Not doing PCA for EEG.
    Created an SSP operator (subspace dimension = 4)
Setting small MEG eigenvalues to zero.
Not doing PCA for MEG.
Setting small EEG eigenvalues to zero.
Not doing PCA for EEG.
    Created an SSP operator (subspace dimension = 4)
Setting small MEG eigenvalues to zero.
Not doing PCA for MEG.
Setting small EEG eigenvalues to zero.
Not doing PCA for EEG.

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

Generated by Sphinx-Gallery