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].
| [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)
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 30.376 seconds)