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, fir_design='firwin')
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:

Successfully extracted to: [u'/home/ubuntu/mne_data/MNE-sample-data']
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 496 samples (3.303 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: -1824.075
... rank: 10 - loglik: -1773.208
... rank: 15 - loglik: -1728.057
... rank: 20 - loglik: -1703.455
... rank: 25 - loglik: -1685.455
... rank: 30 - loglik: -1672.671
... rank: 35 - loglik: -1665.096
... rank: 40 - loglik: -1661.920
... rank: 45 - loglik: -1658.426
... rank: 50 - loglik: -1658.258
... rank: 55 - loglik: -1658.507
... rank: 60 - loglik: -1658.373
... rank: 65 - loglik: -1659.500
... rank: 70 - loglik: -1660.446
... rank: 75 - loglik: -1662.101
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: -1639.020
   factor_analysis: -1658.258
   diagonal_fixed: -1735.111
   empirical: -1859.569
Covariance estimates sorted from best to worst
shrunk : -1639.02038102
factor_analysis : -1658.25785868
diagonal_fixed : -1735.11111207
empirical : -1859.56900735

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 (eeg): 59
estimated rank (grad): 203
estimated rank (mag): 102
estimated rank (eeg): 58
estimated rank (grad): 203
estimated rank (mag): 102
estimated rank (eeg): 58
estimated rank (grad): 203
estimated rank (mag): 99
    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 rank is 362
    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 rank is 362
    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 rank is 362
    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 rank is 360

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

Gallery generated by Sphinx-Gallery