Note
Click here to download the full example code
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].
| [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
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')
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: ( 0 minutes 39.251 seconds)