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 during the baseline period.

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:

Opening raw data file /home/circleci/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 497 samples (3.310 sec) selected
73 matching events found
No baseline correction applied
Not setting metadata
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 : ['EEG 001', 'EEG 002', 'EEG 003', 'EEG 007']
    Rejecting  epoch based on EEG : ['EEG 001', 'EEG 002', 'EEG 003', 'EEG 004', 'EEG 005', 'EEG 006', 'EEG 007']
    Rejecting  epoch based on EEG : ['EEG 003', 'EEG 007']
    Rejecting  epoch based on EEG : ['EEG 001', 'EEG 002', 'EEG 003', 'EEG 006', 'EEG 007']
    Rejecting  epoch based on MAG : ['MEG 1711']
    Rejecting  epoch based on EEG : ['EEG 007']
    Rejecting  epoch based on EEG : ['EEG 008', 'EEG 009']
    Rejecting  epoch based on EEG : ['EEG 001', 'EEG 002', 'EEG 003', 'EEG 007']
    Rejecting  epoch based on EEG : ['EEG 001', 'EEG 002', 'EEG 003', 'EEG 006', 'EEG 007']
    Rejecting  epoch based on EEG : ['EEG 001', 'EEG 002', 'EEG 003', 'EEG 007']
    Rejecting  epoch based on EEG : ['EEG 001', 'EEG 002', 'EEG 003', 'EEG 007']
    Rejecting  epoch based on EEG : ['EEG 001', '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
    MAG regularization : 0.01
    GRAD regularization : 0.01
    EEG regularization : None
    SEEG regularization : None
    ECOG regularization : None
    HBO regularization : None
    HBR regularization : None
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.
    MAG regularization : 0.01
    GRAD regularization : 0.01
    EEG regularization : None
    SEEG regularization : None
    ECOG regularization : None
    HBO regularization : None
    HBR regularization : None
    MAG regularization : 0.01
    GRAD regularization : 0.01
    EEG regularization : None
    SEEG regularization : None
    ECOG regularization : None
    HBO regularization : None
    HBR regularization : None
    MAG regularization : 0.01
    GRAD regularization : 0.01
    EEG regularization : None
    SEEG regularization : None
    ECOG regularization : None
    HBO regularization : None
    HBR regularization : None
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.0203810171752
factor_analysis : -1658.2578586816837
diagonal_fixed : -1735.1111120706043
empirical : -1859.5690073462858

Show the evoked data:

evoked = epochs.average()

evoked.plot(time_unit='s')  # plot evoked response
../../_images/sphx_glr_plot_evoked_whitening_001.png

We can then show whitening for our various noise covariance estimates.

Here we should look to see if baseline signals match the assumption of Gaussian white noise. we expect values centered at 0 within 2 standard deviations for 95% of the time points.

For the Global field power we expect a value of 1.

evoked.plot_white(noise_covs, time_unit='s')
../../_images/sphx_glr_plot_evoked_whitening_002.png

Out:

Created an SSP operator (subspace dimension = 1)
8 projection items activated
estimated rank (eeg): 58
8 projection items activated
estimated rank (grad): 203
Created an SSP operator (subspace dimension = 3)
8 projection items activated
estimated rank (mag): 99
Created an SSP operator (subspace dimension = 1)
8 projection items activated
estimated rank (eeg): 58
8 projection items activated
estimated rank (grad): 203
Created an SSP operator (subspace dimension = 3)
8 projection items activated
estimated rank (mag): 99
Created an SSP operator (subspace dimension = 1)
8 projection items activated
estimated rank (eeg): 58
8 projection items activated
estimated rank (grad): 203
Created an SSP operator (subspace dimension = 3)
8 projection items activated
estimated rank (mag): 99
Created an SSP operator (subspace dimension = 1)
8 projection items activated
estimated rank (eeg): 58
8 projection items activated
estimated rank (grad): 203
Created an SSP operator (subspace dimension = 3)
8 projection items activated
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 360
    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
    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
    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: ( 1 minutes 8.333 seconds)

Gallery generated by Sphinx-Gallery