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, method='iir', 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
add_eeg_ref defaults to True in 0.13, will default to False in 0.14, and will be removed in 0.15. We recommend to use add_eeg_ref=False and set_eeg_reference() instead.
Reading 0 ... 41699  =      0.000 ...   277.709 secs...
Band-pass filtering from 1 - 40 Hz
73 matching events found
No baseline correction applied
add_eeg_ref defaults to True in 0.13, will default to False in 0.14, and will be removed in 0.15. We recommend to use add_eeg_ref=False and set_eeg_reference() instead.
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 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 006', u'EEG 007']
    Rejecting  epoch based on EEG : [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 008', u'EEG 009']
    Rejecting  epoch based on EEG : [u'EEG 001', 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']
11 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: -1801.587
... rank: 10 - loglik: -1749.706
... rank: 15 - loglik: -1703.215
... rank: 20 - loglik: -1675.987
... rank: 25 - loglik: -1656.535
... rank: 30 - loglik: -1642.556
... rank: 35 - loglik: -1634.068
... rank: 40 - loglik: -1628.518
... rank: 45 - loglik: -1624.140
... rank: 50 - loglik: -1622.486
... rank: 55 - loglik: -1621.697
... rank: 60 - loglik: -1621.290
... rank: 65 - loglik: -1621.628
... rank: 70 - loglik: -1622.674
... rank: 75 - loglik: -1623.352
early stopping parameter search.
... best model at rank = 60
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 : 1922
[done]
Number of samples used : 1922
[done]
Number of samples used : 1922
[done]
Number of samples used : 1922
[done]
log-likelihood on unseen data (descending order):
   shrunk: -1604.635
   factor_analysis: -1621.290
   diagonal_fixed: -1690.067
   empirical: -1810.229
Covariance estimates sorted from best to worst
shrunk : -1604.63546573
factor_analysis : -1621.28967547
diagonal_fixed : -1690.06725296
empirical : -1810.22877908

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
Setting small MEG eigenvalues to zero.
Not doing PCA for MEG.
Setting small EEG eigenvalues to zero.
Not doing PCA for EEG.
No average EEG reference present in info["projs"], covariance may be adversely affected. Consider recomputing covariance using a raw file with an average eeg reference projector added.
Setting small MEG eigenvalues to zero.
Not doing PCA for MEG.
Setting small EEG eigenvalues to zero.
Not doing PCA for EEG.
No average EEG reference present in info["projs"], covariance may be adversely affected. Consider recomputing covariance using a raw file with an average eeg reference projector added.
Setting small MEG eigenvalues to zero.
Not doing PCA for MEG.
Setting small EEG eigenvalues to zero.
Not doing PCA for EEG.
No average EEG reference present in info["projs"], covariance may be adversely affected. Consider recomputing covariance using a raw file with an average eeg reference projector added.
Setting small MEG eigenvalues to zero.
Not doing PCA for MEG.
Setting small EEG eigenvalues to zero.
Not doing PCA for EEG.
No average EEG reference present in info["projs"], covariance may be adversely affected. Consider recomputing covariance using a raw file with an average eeg reference projector added.

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

Generated by Sphinx-Gallery