Note
Go to the end to download the full example code.
Oddball Tones Example¶
Estimate NCRFs for standard and oddball tones.
For this tutorial, we use the auditory Brainstorm tutorial dataset [Tadel et al., 2011] that is available as a part of the Brainstorm software.
Note
Downloading the dataset requires answering an interactive prompt (see
mne.datasets.brainstorm.bst_auditory.data_path()).
# Authors: Proloy Das <proloy@umd.edu>
# Christian Brodbeck <brodbecc@mcmaster.ca>
#
# sphinx_gallery_thumbnail_number = 3
import numpy as np
import pandas as pd
import eelbrain
import mne
from ncrf import fit_ncrf
Preprocessing¶
Preprocess MEG Data: low pass filtering, power line attenuation, downsampling, etc. We broadly follow this mne-python tutorial.
data_path = mne.datasets.brainstorm.bst_auditory.data_path()
raw_fname = data_path / 'MEG' / 'bst_auditory' / 'S01_AEF_20131218_01.ds'
raw = mne.io.read_raw_ctf(raw_fname, preload=False)
n_times_run1 = raw.n_times
# We mark a set of bad channels that seem noisier than others.
raw.info['bads'] = ['MLO52-4408', 'MRT51-4408', 'MLO42-4408', 'MLO43-4408']
annotations_df = pd.DataFrame()
offset = n_times_run1
for idx in [1]:
csv_fname = data_path / 'MEG' / 'bst_auditory' / f'events_bad_0{idx}.csv'
df = pd.read_csv(csv_fname, header=None, names=['onset', 'duration', 'id', 'label'])
print('Events from run {0}:'.format(idx))
print(df)
df['onset'] += offset * (idx - 1)
annotations_df = pd.concat([annotations_df, df], axis=0)
# Conversion from samples to times:
onsets = annotations_df['onset'].values / raw.info['sfreq']
durations = annotations_df['duration'].values / raw.info['sfreq']
descriptions = annotations_df['label'].values
annotations = mne.Annotations(onsets, durations, descriptions)
raw.set_annotations(annotations)
del onsets, durations, descriptions
# events are the presentation times of the audio stimuli: UPPT001
event_fname = data_path / 'MEG' / 'bst_auditory' / 'S01_AEF_20131218_01-eve.fif'
events = mne.find_events(raw, stim_channel='UPPT001')
# The event timing is adjusted by comparing the trigger times on detected sound onsets on channel UADC001-4408.
sound_data = raw[raw.ch_names.index('UADC001-4408')][0][0]
onsets = np.where(np.abs(sound_data) > 2. * np.std(sound_data))[0]
min_diff = int(0.5 * raw.info['sfreq'])
diffs = np.concatenate([[min_diff + 1], np.diff(onsets)])
onsets = onsets[diffs > min_diff]
assert len(onsets) == len(events)
diffs = 1000. * (events[:, 0] - onsets) / raw.info['sfreq']
print('Trigger delay removed (μ ± σ): %0.1f ± %0.1f ms'
% (np.mean(diffs), np.std(diffs)))
# events times are rescaled according to new sampling freq, 100 Hz
events[:, 0] = np.int64(onsets * 100 / raw.info['sfreq'])
mne.write_events(event_fname, events, overwrite=True)
del sound_data, diffs
## set EOG channel
raw.set_eeg_reference('average', projection=True)
# raw_AEF.plot_psd(tmax=60., average=False)
raw.load_data()
raw.notch_filter(np.arange(60, 181, 60), fir_design='firwin')
# band pass filtering 1-8 Hz
raw.filter(1.0, 8.0, fir_design='firwin')
# resample to 100 Hz
raw.resample(100, npad="auto")
### LOAD RELEVANT VARIABLES AS eelbrain.NDVar
# load as epochs for plot only
ds = eelbrain.load.fiff.events(raw=raw, proj=True, stim_channel='UPPT001', events=event_fname)
epochs = eelbrain.load.fiff.epochs(ds, tmin=-0.1, tmax=0.5, baseline=(None, 0))
eelbrain.plot.Butterfly(epochs)
# pick MEG channels
picks = mne.pick_types(raw.info, meg=True, eeg=False, stim=False, eog=False,
ref_meg=False, exclude='bads')
# Read as a single chunk of data
y, t = raw.get_data(picks, return_times=True)
sensor_dim = eelbrain.load.fiff.sensor_dim(raw.info, picks=picks)
time = eelbrain.UTS.from_int(0, t.size - 1, raw.info['sfreq'])
meg = eelbrain.NDVar(y, dims=(sensor_dim, time))
print(meg)
License
-------
This tutorial dataset (EEG and MRI data) remains a property of the MEG Lab,
McConnell Brain Imaging Center, Montreal Neurological Institute,
McGill University, Canada. Its use and transfer outside the Brainstorm
tutorial, e.g. for research purposes, is prohibited without written consent
from the MEG Lab.
If you reference this dataset in your publications, please:
1) acknowledge its authors: Elizabeth Bock, Esther Florin, Francois Tadel
and Sylvain Baillet, and
2) cite Brainstorm as indicated on the website:
http://neuroimage.usc.edu/brainstorm
For questions, please contact Francois Tadel (francois.tadel@mcgill.ca).
Agree (y/[n])? Downloading file 'bst_auditory.tar.gz' from 'https://osf.io/5t9n8/download?version=1' to '/home/runner/mne_data'.
0%| | 0.00/1.64G [00:00<?, ?B/s]
0%| | 945k/1.64G [00:00<02:53, 9.44MB/s]
0%| | 2.64M/1.64G [00:00<02:01, 13.4MB/s]
0%| | 4.47M/1.64G [00:00<01:50, 14.7MB/s]
0%|▏ | 6.29M/1.64G [00:00<01:43, 15.7MB/s]
0%|▏ | 7.86M/1.64G [00:00<01:45, 15.4MB/s]
1%|▏ | 9.69M/1.64G [00:00<01:40, 16.1MB/s]
1%|▎ | 11.3M/1.64G [00:00<01:44, 15.6MB/s]
1%|▎ | 12.9M/1.64G [00:00<01:55, 14.1MB/s]
1%|▎ | 14.3M/1.64G [00:00<01:55, 14.0MB/s]
1%|▎ | 15.7M/1.64G [00:01<02:01, 13.3MB/s]
1%|▍ | 17.1M/1.64G [00:01<02:29, 10.9MB/s]
1%|▍ | 18.2M/1.64G [00:01<02:26, 11.0MB/s]
1%|▍ | 19.4M/1.64G [00:01<02:28, 10.9MB/s]
1%|▍ | 20.7M/1.64G [00:01<02:26, 11.1MB/s]
1%|▍ | 22.0M/1.64G [00:01<02:22, 11.3MB/s]
1%|▌ | 23.3M/1.64G [00:01<02:16, 11.8MB/s]
2%|▌ | 25.5M/1.64G [00:01<01:49, 14.7MB/s]
2%|▊ | 33.6M/1.64G [00:02<00:47, 33.4MB/s]
3%|▉ | 42.8M/1.64G [00:02<00:35, 45.0MB/s]
3%|█ | 47.1M/1.64G [00:02<00:38, 41.1MB/s]
3%|█▎ | 56.1M/1.64G [00:02<00:29, 53.7MB/s]
4%|█▍ | 65.1M/1.64G [00:02<00:24, 63.7MB/s]
5%|█▋ | 74.3M/1.64G [00:02<00:21, 71.6MB/s]
5%|█▉ | 83.4M/1.64G [00:02<00:20, 77.2MB/s]
6%|██ | 92.6M/1.64G [00:02<00:18, 81.3MB/s]
6%|██▎ | 102M/1.64G [00:02<00:18, 84.3MB/s]
7%|██▌ | 111M/1.64G [00:02<00:17, 87.0MB/s]
7%|██▊ | 121M/1.64G [00:03<00:16, 89.8MB/s]
8%|███ | 130M/1.64G [00:03<00:19, 75.5MB/s]
8%|███▏ | 138M/1.64G [00:04<01:05, 22.8MB/s]
9%|███▎ | 144M/1.64G [00:04<01:15, 19.7MB/s]
9%|███▍ | 148M/1.64G [00:05<01:25, 17.4MB/s]
9%|███▌ | 152M/1.64G [00:05<01:32, 16.1MB/s]
9%|███▌ | 154M/1.64G [00:05<01:37, 15.2MB/s]
10%|███▋ | 157M/1.64G [00:05<01:43, 14.3MB/s]
10%|███▋ | 159M/1.64G [00:06<02:00, 12.2MB/s]
10%|███▋ | 160M/1.64G [00:06<01:58, 12.5MB/s]
10%|███▊ | 162M/1.64G [00:06<02:01, 12.2MB/s]
10%|███▊ | 163M/1.64G [00:06<02:03, 11.9MB/s]
10%|███▊ | 164M/1.64G [00:06<02:35, 9.48MB/s]
10%|███▊ | 166M/1.64G [00:06<02:30, 9.75MB/s]
10%|███▊ | 167M/1.64G [00:06<02:27, 9.94MB/s]
10%|███▉ | 168M/1.64G [00:07<02:14, 10.9MB/s]
10%|███▉ | 170M/1.64G [00:07<02:04, 11.8MB/s]
10%|███▉ | 172M/1.64G [00:07<01:51, 13.1MB/s]
11%|████ | 173M/1.64G [00:07<01:42, 14.2MB/s]
11%|████ | 175M/1.64G [00:07<01:34, 15.5MB/s]
11%|████▎ | 183M/1.64G [00:07<00:44, 33.0MB/s]
12%|████▍ | 188M/1.64G [00:07<00:37, 38.3MB/s]
12%|████▌ | 195M/1.64G [00:07<00:33, 42.7MB/s]
12%|████▋ | 203M/1.64G [00:07<00:27, 51.8MB/s]
13%|████▉ | 212M/1.64G [00:07<00:22, 63.0MB/s]
13%|█████▏ | 221M/1.64G [00:08<00:20, 70.5MB/s]
14%|█████▎ | 230M/1.64G [00:08<00:18, 77.0MB/s]
15%|█████▌ | 238M/1.64G [00:08<00:19, 71.8MB/s]
15%|█████▋ | 245M/1.64G [00:08<00:39, 35.0MB/s]
15%|█████▊ | 251M/1.64G [00:09<00:44, 30.9MB/s]
16%|██████ | 259M/1.64G [00:09<00:35, 39.0MB/s]
16%|██████▏ | 268M/1.64G [00:09<00:28, 48.5MB/s]
17%|██████▍ | 278M/1.64G [00:09<00:23, 57.8MB/s]
17%|██████▋ | 285M/1.64G [00:09<00:21, 62.1MB/s]
18%|██████▊ | 293M/1.64G [00:09<00:22, 59.4MB/s]
18%|███████ | 302M/1.64G [00:09<00:19, 67.3MB/s]
19%|███████▏ | 311M/1.64G [00:09<00:17, 73.9MB/s]
20%|███████▍ | 319M/1.64G [00:09<00:19, 68.5MB/s]
20%|███████▋ | 328M/1.64G [00:10<00:17, 74.0MB/s]
21%|███████▊ | 338M/1.64G [00:10<00:16, 79.8MB/s]
21%|████████ | 348M/1.64G [00:10<00:15, 84.7MB/s]
22%|████████▎ | 357M/1.64G [00:10<00:18, 69.2MB/s]
22%|████████▍ | 365M/1.64G [00:10<00:17, 73.0MB/s]
23%|████████▋ | 374M/1.64G [00:10<00:16, 78.5MB/s]
23%|████████▉ | 383M/1.64G [00:10<00:17, 71.8MB/s]
24%|█████████ | 392M/1.64G [00:10<00:16, 76.1MB/s]
24%|█████████▎ | 401M/1.64G [00:10<00:15, 80.2MB/s]
25%|█████████▌ | 410M/1.64G [00:11<00:14, 83.6MB/s]
26%|█████████▋ | 419M/1.64G [00:11<00:14, 84.7MB/s]
26%|█████████▉ | 428M/1.64G [00:11<00:29, 40.3MB/s]
27%|██████████▏ | 436M/1.64G [00:11<00:24, 48.2MB/s]
27%|██████████▎ | 445M/1.64G [00:11<00:21, 56.1MB/s]
28%|██████████▌ | 455M/1.64G [00:11<00:18, 64.2MB/s]
28%|██████████▊ | 464M/1.64G [00:12<00:16, 70.0MB/s]
29%|██████████▉ | 473M/1.64G [00:12<00:15, 75.7MB/s]
29%|███████████▏ | 482M/1.64G [00:12<00:14, 80.5MB/s]
30%|███████████▍ | 492M/1.64G [00:12<00:13, 83.8MB/s]
31%|███████████▋ | 501M/1.64G [00:12<00:13, 85.0MB/s]
31%|███████████▊ | 510M/1.64G [00:12<00:24, 46.2MB/s]
32%|████████████ | 517M/1.64G [00:12<00:21, 51.5MB/s]
32%|████████████▏ | 527M/1.64G [00:13<00:18, 60.1MB/s]
33%|████████████▍ | 536M/1.64G [00:13<00:16, 68.0MB/s]
33%|████████████▋ | 546M/1.64G [00:13<00:14, 74.3MB/s]
34%|████████████▉ | 555M/1.64G [00:13<00:13, 80.2MB/s]
35%|█████████████▏ | 565M/1.64G [00:13<00:12, 85.1MB/s]
35%|█████████████▎ | 575M/1.64G [00:13<00:11, 88.8MB/s]
36%|█████████████▌ | 585M/1.64G [00:13<00:11, 91.9MB/s]
36%|█████████████▊ | 595M/1.64G [00:13<00:11, 93.6MB/s]
37%|██████████████ | 604M/1.64G [00:13<00:10, 94.5MB/s]
38%|██████████████▎ | 614M/1.64G [00:13<00:10, 94.5MB/s]
38%|██████████████▍ | 624M/1.64G [00:14<00:10, 94.6MB/s]
39%|██████████████▋ | 633M/1.64G [00:14<00:10, 94.3MB/s]
39%|██████████████▉ | 643M/1.64G [00:14<00:10, 94.4MB/s]
40%|███████████████▏ | 652M/1.64G [00:14<00:18, 53.2MB/s]
40%|███████████████▎ | 660M/1.64G [00:15<00:35, 27.8MB/s]
41%|███████████████▍ | 665M/1.64G [00:15<00:43, 22.2MB/s]
41%|███████████████▌ | 669M/1.64G [00:15<00:40, 23.8MB/s]
41%|███████████████▋ | 673M/1.64G [00:15<00:37, 26.0MB/s]
41%|███████████████▋ | 677M/1.64G [00:16<00:34, 28.2MB/s]
42%|███████████████▉ | 684M/1.64G [00:16<00:27, 34.5MB/s]
42%|████████████████ | 692M/1.64G [00:16<00:20, 45.1MB/s]
43%|████████████████▎ | 702M/1.64G [00:16<00:16, 56.1MB/s]
43%|████████████████▍ | 708M/1.64G [00:16<00:29, 31.7MB/s]
44%|████████████████▌ | 714M/1.64G [00:17<00:43, 21.2MB/s]
44%|████████████████▋ | 718M/1.64G [00:17<00:47, 19.2MB/s]
44%|████████████████▊ | 721M/1.64G [00:17<00:50, 18.1MB/s]
44%|████████████████▊ | 724M/1.64G [00:18<00:53, 16.9MB/s]
44%|████████████████▊ | 726M/1.64G [00:18<01:02, 14.5MB/s]
45%|████████████████▉ | 728M/1.64G [00:18<01:08, 13.3MB/s]
45%|████████████████▉ | 729M/1.64G [00:18<01:08, 13.2MB/s]
45%|█████████████████ | 732M/1.64G [00:18<00:59, 15.3MB/s]
45%|█████████████████ | 734M/1.64G [00:18<00:59, 15.2MB/s]
45%|█████████████████ | 735M/1.64G [00:18<00:59, 15.1MB/s]
45%|█████████████████▏ | 741M/1.64G [00:19<00:37, 23.9MB/s]
45%|█████████████████▎ | 744M/1.64G [00:19<00:36, 24.6MB/s]
46%|█████████████████▍ | 749M/1.64G [00:19<00:30, 28.9MB/s]
46%|█████████████████▍ | 752M/1.64G [00:19<00:34, 25.4MB/s]
46%|█████████████████▌ | 754M/1.64G [00:19<00:41, 21.2MB/s]
46%|█████████████████▌ | 759M/1.64G [00:19<00:34, 25.6MB/s]
47%|█████████████████▊ | 766M/1.64G [00:19<00:23, 36.9MB/s]
47%|█████████████████▉ | 774M/1.64G [00:19<00:17, 48.6MB/s]
48%|██████████████████▏ | 781M/1.64G [00:19<00:15, 53.7MB/s]
48%|██████████████████▎ | 788M/1.64G [00:20<00:14, 57.3MB/s]
49%|██████████████████▌ | 797M/1.64G [00:20<00:12, 66.8MB/s]
49%|██████████████████▋ | 806M/1.64G [00:20<00:11, 74.5MB/s]
50%|██████████████████▉ | 815M/1.64G [00:20<00:10, 80.2MB/s]
50%|███████████████████▏ | 825M/1.64G [00:20<00:09, 84.4MB/s]
51%|███████████████████▎ | 833M/1.64G [00:20<00:15, 53.0MB/s]
51%|███████████████████▌ | 840M/1.64G [00:21<00:36, 21.9MB/s]
52%|███████████████████▋ | 845M/1.64G [00:22<00:44, 17.7MB/s]
52%|███████████████████▋ | 849M/1.64G [00:22<00:47, 16.7MB/s]
52%|███████████████████▊ | 852M/1.64G [00:22<00:47, 16.3MB/s]
52%|███████████████████▊ | 855M/1.64G [00:22<00:47, 16.6MB/s]
52%|███████████████████▉ | 857M/1.64G [00:22<00:47, 16.3MB/s]
53%|███████████████████▉ | 859M/1.64G [00:23<00:47, 16.5MB/s]
53%|████████████████████ | 861M/1.64G [00:23<00:52, 14.8MB/s]
53%|████████████████████ | 863M/1.64G [00:23<00:51, 15.0MB/s]
53%|████████████████████ | 865M/1.64G [00:23<00:53, 14.5MB/s]
53%|████████████████████▏ | 866M/1.64G [00:23<00:50, 15.2MB/s]
53%|████████████████████▏ | 868M/1.64G [00:23<00:50, 15.2MB/s]
54%|████████████████████▎ | 875M/1.64G [00:23<00:26, 29.2MB/s]
54%|████████████████████▍ | 879M/1.64G [00:23<00:25, 30.1MB/s]
54%|████████████████████▍ | 882M/1.64G [00:24<00:27, 27.2MB/s]
54%|████████████████████▌ | 888M/1.64G [00:24<00:21, 35.2MB/s]
55%|████████████████████▋ | 891M/1.64G [00:24<00:26, 28.3MB/s]
55%|████████████████████▊ | 897M/1.64G [00:24<00:21, 34.6MB/s]
55%|█████████████████████ | 906M/1.64G [00:24<00:14, 48.8MB/s]
56%|█████████████████████▏ | 912M/1.64G [00:24<00:15, 45.5MB/s]
56%|█████████████████████▎ | 919M/1.64G [00:24<00:13, 53.1MB/s]
57%|█████████████████████▌ | 928M/1.64G [00:24<00:11, 62.9MB/s]
57%|█████████████████████▊ | 938M/1.64G [00:24<00:09, 71.4MB/s]
58%|██████████████████████ | 947M/1.64G [00:25<00:08, 77.1MB/s]
58%|██████████████████████▏ | 955M/1.64G [00:25<00:13, 50.0MB/s]
59%|██████████████████████▎ | 961M/1.64G [00:25<00:26, 25.7MB/s]
59%|██████████████████████▍ | 966M/1.64G [00:26<00:35, 18.6MB/s]
59%|██████████████████████▌ | 970M/1.64G [00:26<00:39, 16.7MB/s]
59%|██████████████████████▌ | 973M/1.64G [00:27<00:44, 14.9MB/s]
60%|██████████████████████▋ | 975M/1.64G [00:27<00:44, 14.9MB/s]
60%|██████████████████████▋ | 977M/1.64G [00:27<00:47, 13.8MB/s]
60%|██████████████████████▋ | 979M/1.64G [00:27<00:46, 14.1MB/s]
60%|██████████████████████▊ | 981M/1.64G [00:27<00:46, 14.2MB/s]
60%|██████████████████████▊ | 982M/1.64G [00:27<00:44, 14.7MB/s]
60%|██████████████████████▉ | 989M/1.64G [00:27<00:26, 24.5MB/s]
61%|███████████████████████ | 995M/1.64G [00:28<00:19, 33.2MB/s]
61%|██████████████████████▋ | 1.00G/1.64G [00:28<00:14, 42.8MB/s]
62%|██████████████████████▉ | 1.01G/1.64G [00:28<00:11, 55.4MB/s]
62%|███████████████████████ | 1.02G/1.64G [00:28<00:10, 60.9MB/s]
63%|███████████████████████▏ | 1.03G/1.64G [00:28<00:09, 63.1MB/s]
63%|███████████████████████▍ | 1.03G/1.64G [00:28<00:09, 64.0MB/s]
64%|███████████████████████▌ | 1.04G/1.64G [00:28<00:08, 71.2MB/s]
64%|███████████████████████▊ | 1.05G/1.64G [00:28<00:07, 77.5MB/s]
65%|████████████████████████ | 1.06G/1.64G [00:28<00:06, 82.1MB/s]
65%|████████████████████████▏ | 1.07G/1.64G [00:28<00:06, 85.6MB/s]
66%|████████████████████████▍ | 1.08G/1.64G [00:29<00:06, 86.6MB/s]
67%|████████████████████████▋ | 1.09G/1.64G [00:29<00:06, 88.6MB/s]
67%|████████████████████████▊ | 1.10G/1.64G [00:29<00:06, 89.0MB/s]
68%|█████████████████████████ | 1.11G/1.64G [00:29<00:05, 89.8MB/s]
68%|█████████████████████████▎ | 1.12G/1.64G [00:29<00:05, 90.7MB/s]
69%|█████████████████████████▍ | 1.13G/1.64G [00:29<00:05, 90.2MB/s]
69%|█████████████████████████▋ | 1.14G/1.64G [00:29<00:05, 91.0MB/s]
70%|█████████████████████████▉ | 1.14G/1.64G [00:29<00:05, 91.3MB/s]
71%|██████████████████████████ | 1.15G/1.64G [00:29<00:05, 91.6MB/s]
71%|██████████████████████████▎ | 1.16G/1.64G [00:29<00:05, 91.6MB/s]
72%|██████████████████████████▌ | 1.17G/1.64G [00:30<00:05, 91.2MB/s]
72%|██████████████████████████▋ | 1.18G/1.64G [00:30<00:05, 90.8MB/s]
73%|██████████████████████████▉ | 1.19G/1.64G [00:30<00:04, 90.5MB/s]
73%|███████████████████████████▏ | 1.20G/1.64G [00:30<00:04, 90.4MB/s]
74%|███████████████████████████▎ | 1.21G/1.64G [00:30<00:04, 90.3MB/s]
74%|███████████████████████████▌ | 1.22G/1.64G [00:30<00:04, 90.3MB/s]
75%|███████████████████████████▋ | 1.23G/1.64G [00:30<00:04, 90.2MB/s]
76%|███████████████████████████▉ | 1.24G/1.64G [00:30<00:07, 54.2MB/s]
76%|████████████████████████████ | 1.24G/1.64G [00:31<00:15, 25.6MB/s]
76%|████████████████████████████▏ | 1.25G/1.64G [00:32<00:19, 19.9MB/s]
77%|████████████████████████████▎ | 1.25G/1.64G [00:32<00:21, 17.8MB/s]
77%|████████████████████████████▍ | 1.25G/1.64G [00:32<00:22, 17.2MB/s]
77%|████████████████████████████▍ | 1.26G/1.64G [00:32<00:22, 16.8MB/s]
77%|████████████████████████████▌ | 1.26G/1.64G [00:33<00:22, 16.6MB/s]
77%|████████████████████████████▌ | 1.26G/1.64G [00:33<00:23, 16.0MB/s]
77%|████████████████████████████▌ | 1.26G/1.64G [00:33<00:26, 14.2MB/s]
77%|████████████████████████████▋ | 1.27G/1.64G [00:33<00:26, 14.1MB/s]
77%|████████████████████████████▋ | 1.27G/1.64G [00:33<00:25, 14.6MB/s]
78%|████████████████████████████▋ | 1.27G/1.64G [00:33<00:24, 14.9MB/s]
78%|████████████████████████████▋ | 1.27G/1.64G [00:33<00:23, 15.3MB/s]
78%|████████████████████████████▊ | 1.27G/1.64G [00:33<00:21, 17.2MB/s]
78%|████████████████████████████▉ | 1.28G/1.64G [00:34<00:15, 23.2MB/s]
78%|████████████████████████████▉ | 1.28G/1.64G [00:34<00:16, 21.9MB/s]
78%|████████████████████████████▉ | 1.28G/1.64G [00:34<00:18, 19.0MB/s]
78%|█████████████████████████████ | 1.28G/1.64G [00:34<00:22, 15.7MB/s]
79%|█████████████████████████████ | 1.29G/1.64G [00:34<00:17, 19.8MB/s]
79%|█████████████████████████████▎ | 1.29G/1.64G [00:34<00:10, 32.3MB/s]
79%|█████████████████████████████▎ | 1.30G/1.64G [00:34<00:10, 33.5MB/s]
80%|█████████████████████████████▍ | 1.30G/1.64G [00:34<00:08, 38.8MB/s]
80%|█████████████████████████████▌ | 1.31G/1.64G [00:35<00:07, 41.1MB/s]
80%|█████████████████████████████▋ | 1.31G/1.64G [00:35<00:06, 47.3MB/s]
81%|█████████████████████████████▉ | 1.32G/1.64G [00:35<00:05, 56.0MB/s]
81%|██████████████████████████████ | 1.33G/1.64G [00:35<00:04, 62.9MB/s]
82%|██████████████████████████████▎ | 1.34G/1.64G [00:35<00:04, 66.3MB/s]
82%|██████████████████████████████▍ | 1.34G/1.64G [00:35<00:09, 29.4MB/s]
83%|██████████████████████████████▌ | 1.35G/1.64G [00:36<00:19, 14.7MB/s]
83%|██████████████████████████████▌ | 1.35G/1.64G [00:37<00:18, 15.2MB/s]
83%|██████████████████████████████▋ | 1.36G/1.64G [00:37<00:18, 15.1MB/s]
83%|██████████████████████████████▋ | 1.36G/1.64G [00:37<00:18, 15.0MB/s]
83%|██████████████████████████████▊ | 1.36G/1.64G [00:37<00:18, 15.2MB/s]
83%|██████████████████████████████▊ | 1.36G/1.64G [00:37<00:17, 15.5MB/s]
83%|██████████████████████████████▉ | 1.37G/1.64G [00:37<00:19, 14.0MB/s]
84%|██████████████████████████████▉ | 1.37G/1.64G [00:38<00:19, 13.8MB/s]
84%|██████████████████████████████▉ | 1.37G/1.64G [00:38<00:18, 14.1MB/s]
84%|██████████████████████████████▉ | 1.37G/1.64G [00:38<00:18, 14.3MB/s]
84%|███████████████████████████████ | 1.37G/1.64G [00:38<00:18, 14.4MB/s]
84%|███████████████████████████████ | 1.37G/1.64G [00:38<00:17, 14.7MB/s]
84%|███████████████████████████████ | 1.37G/1.64G [00:38<00:20, 12.9MB/s]
84%|███████████████████████████████▏ | 1.38G/1.64G [00:38<00:19, 13.3MB/s]
84%|███████████████████████████████▏ | 1.38G/1.64G [00:38<00:19, 13.0MB/s]
84%|███████████████████████████████▏ | 1.38G/1.64G [00:38<00:19, 13.2MB/s]
84%|███████████████████████████████▏ | 1.38G/1.64G [00:39<00:19, 13.4MB/s]
85%|███████████████████████████████▎ | 1.38G/1.64G [00:39<00:19, 13.3MB/s]
85%|███████████████████████████████▎ | 1.38G/1.64G [00:39<00:20, 12.0MB/s]
85%|███████████████████████████████▍ | 1.39G/1.64G [00:39<00:08, 28.1MB/s]
85%|███████████████████████████████▌ | 1.39G/1.64G [00:39<00:08, 29.3MB/s]
85%|███████████████████████████████▌ | 1.40G/1.64G [00:39<00:08, 28.9MB/s]
86%|███████████████████████████████▋ | 1.40G/1.64G [00:39<00:07, 32.9MB/s]
86%|███████████████████████████████▉ | 1.41G/1.64G [00:39<00:04, 45.9MB/s]
87%|████████████████████████████████ | 1.42G/1.64G [00:39<00:03, 58.2MB/s]
87%|████████████████████████████████▎ | 1.43G/1.64G [00:40<00:03, 67.4MB/s]
88%|████████████████████████████████▍ | 1.44G/1.64G [00:40<00:02, 74.0MB/s]
88%|████████████████████████████████▋ | 1.45G/1.64G [00:40<00:02, 78.7MB/s]
89%|████████████████████████████████▉ | 1.45G/1.64G [00:40<00:02, 82.0MB/s]
89%|█████████████████████████████████ | 1.46G/1.64G [00:40<00:02, 84.5MB/s]
90%|█████████████████████████████████▎ | 1.47G/1.64G [00:40<00:01, 86.1MB/s]
91%|█████████████████████████████████▌ | 1.48G/1.64G [00:40<00:01, 86.0MB/s]
91%|█████████████████████████████████▋ | 1.49G/1.64G [00:40<00:01, 85.6MB/s]
92%|█████████████████████████████████▉ | 1.50G/1.64G [00:40<00:01, 87.3MB/s]
92%|██████████████████████████████████ | 1.51G/1.64G [00:40<00:01, 88.1MB/s]
93%|██████████████████████████████████▎ | 1.52G/1.64G [00:41<00:01, 88.5MB/s]
93%|██████████████████████████████████▌ | 1.53G/1.64G [00:41<00:01, 88.7MB/s]
94%|██████████████████████████████████▋ | 1.53G/1.64G [00:41<00:03, 32.1MB/s]
94%|██████████████████████████████████▊ | 1.54G/1.64G [00:42<00:04, 23.3MB/s]
95%|██████████████████████████████████▉ | 1.55G/1.64G [00:42<00:04, 21.3MB/s]
95%|███████████████████████████████████ | 1.55G/1.64G [00:42<00:04, 19.4MB/s]
95%|███████████████████████████████████▏ | 1.55G/1.64G [00:43<00:04, 18.5MB/s]
95%|███████████████████████████████████▏ | 1.56G/1.64G [00:43<00:04, 17.7MB/s]
95%|███████████████████████████████████▎ | 1.56G/1.64G [00:43<00:05, 14.2MB/s]
95%|███████████████████████████████████▎ | 1.56G/1.64G [00:43<00:05, 12.6MB/s]
95%|███████████████████████████████████▎ | 1.56G/1.64G [00:43<00:05, 12.5MB/s]
96%|███████████████████████████████████▎ | 1.56G/1.64G [00:44<00:05, 12.7MB/s]
96%|███████████████████████████████████▍ | 1.57G/1.64G [00:44<00:05, 13.0MB/s]
96%|███████████████████████████████████▍ | 1.57G/1.64G [00:44<00:05, 13.8MB/s]
96%|███████████████████████████████████▍ | 1.57G/1.64G [00:44<00:04, 14.5MB/s]
96%|███████████████████████████████████▌ | 1.57G/1.64G [00:44<00:04, 15.1MB/s]
96%|███████████████████████████████████▌ | 1.57G/1.64G [00:44<00:03, 18.4MB/s]
97%|███████████████████████████████████▋ | 1.58G/1.64G [00:44<00:02, 27.2MB/s]
97%|███████████████████████████████████▊ | 1.59G/1.64G [00:44<00:01, 38.4MB/s]
97%|████████████████████████████████████ | 1.59G/1.64G [00:44<00:00, 52.1MB/s]
98%|████████████████████████████████████▏| 1.60G/1.64G [00:45<00:00, 44.9MB/s]
98%|████████████████████████████████████▎| 1.61G/1.64G [00:45<00:00, 53.5MB/s]
99%|████████████████████████████████████▌| 1.62G/1.64G [00:45<00:00, 64.1MB/s]
99%|████████████████████████████████████▊| 1.63G/1.64G [00:45<00:00, 71.9MB/s]
100%|████████████████████████████████████▉| 1.63G/1.64G [00:46<00:00, 29.4MB/s]
0%| | 0.00/1.64G [00:00<?, ?B/s]
100%|█████████████████████████████████████| 1.64G/1.64G [00:00<00:00, 10.4TB/s]
Untarring contents of '/home/runner/mne_data/bst_auditory.tar.gz' to '/home/runner/mne_data/MNE-brainstorm-data'
Events from run 1:
onset duration id label
0 7625 2776 1 BAD
1 142459 892 1 BAD
2 216954 460 1 BAD
3 345135 5816 1 BAD
4 357687 1053 1 BAD
5 409101 3736 1 BAD
6 461110 179 1 BAD
7 479866 426 1 BAD
8 764914 11500 1 BAD
9 798174 6589 1 BAD
10 846880 5383 1 BAD
11 858863 5136 1 BAD
Trigger delay removed (μ ± σ): -13.9 ± 0.3 ms
/home/runner/mne_data/MNE-brainstorm-data/bst_auditory/MEG/bst_auditory/S01_AEF_20131218_01.ds/S01_AEF_20131218_01.meg4: MNE generated only 233 Epochs for 240 events. The raw file might end before the end of the last epoch.
<NDVar: 270 sensor, 36000 time>
Continuous stimulus variable construction¶
After loading and processing the raw data, we will construct the predictor variable for this particular experiment (by putting an impulse at every event time-point). Note that, the predictor variable and meg response should be of same length.
In case of repetitive trials (where you will have a eelbrain.Case dimension), supply one predictor variable for each trial. Different predictor variables for a single trial can be nested (see ncrf.fit_ncrf()).
In this example, we use two different predictor variables for a single trial
# For the common response, we put impulses at the presentation times of both the audio stimuli (i.e., all beeps).
stim1 = np.zeros(len(time))
stim1[events[:, 0]] = 1.
# To distinguish between standard and deviant beeps, we assign 1 and -1 impulses respectively.
stim2 = stim1.copy()
stim2[events[np.where(events[:, 2] == 2), 0]] = -1.
stim1 = eelbrain.NDVar(stim1, time)
stim2 = eelbrain.NDVar(stim2, time)
# Visualize the stimulus
# p = eelbrain.plot.LineStack(eelbrain.combine([stim1, stim2]), w=10, h=2.5, legend=False)
p = eelbrain.plot.UTS([stim1, stim2], color='black', stem=True, frame='none', w=10, h=2.5, legend=False)
Noise covariance estimation¶
Here we estimate the noise covariance from empty room data. Instead, you can also use pre-stimulus recordings to compute noise covariance.
noise_path = data_path / 'MEG' / 'bst_auditory' / 'S01_Noise_20131218_01.ds'
raw_empty_room = mne.io.read_raw_ctf(noise_path, preload=True)
# Apply the same pre-processing steps to empty room data
raw_empty_room.notch_filter(np.arange(60, 181, 60), fir_design='firwin')
raw_empty_room.filter(1.0, 8.0, fir_design='firwin')
raw_empty_room.resample(100, npad="auto")
# Compute the noise covariance matrix
noise_cov = mne.compute_raw_covariance(raw_empty_room, tmin=0, tmax=None, method='shrunk', rank=None)
Forward model (aka lead-field matrix)¶
Now is the time for forward modeling. ‘ico-4’ should be sufficient resolution if working with surface source space. You can choose to work with free or constrained lead fields. :func`ncrf.fit_ncrf` will choose the appropriate regularizer by looking at the provided lead-field matrix.
# The paths to FreeSurfer reconstructions
subjects_dir = data_path / 'subjects'
subject = 'bst_auditory'
# mne.viz.plot_bem(subject=subject, subjects_dir=subjects_dir,
# brain_surfaces='white', orientation='coronal')
# The transformation file obtained by coregistration
trans = data_path / 'MEG' / 'bst_auditory' / 'bst_auditory-trans.fif'
# Here we look at the head only.
# mne.viz.plot_alignment(raw.info, trans, subject=subject, dig=True,
# meg=['helmet', 'sensors'], subjects_dir=subjects_dir,
# surfaces='head')
srcfile = subjects_dir / 'bst_auditory' / 'bem' / 'bst_auditory-ico-4-src.fif'
if srcfile.is_file():
src = mne.read_source_spaces(srcfile)
else:
src = mne.setup_source_space(subject, spacing='ico4',
subjects_dir=subjects_dir, add_dist=False)
mne.add_source_space_distances(src)
mne.write_source_spaces(srcfile, src, overwrite=True) # needed for smoothing
src
<SourceSpaces: [<surface (lh), n_vertices=163080, n_used=2562>, <surface (rh), n_vertices=163816, n_used=2562>] MRI (surface RAS) coords, subject 'bst_auditory', ~34.5 MiB>
Compute the forward solution:
fwdfile = subjects_dir / 'bst_auditory' / 'bem' / 'bst_auditory-ico-4-fwd.fif'
if fwdfile.is_file():
fwd = mne.read_forward_solution(fwdfile)
else:
conductivity = (0.3,) # for single layer
# conductivity = (0.3, 0.006, 0.3) # for three layers
model = mne.make_bem_model(subject=subject, ico=4,
conductivity=conductivity,
subjects_dir=subjects_dir)
bem = mne.make_bem_solution(model)
fwd = mne.make_forward_solution(raw.info, trans=trans, src=src, bem=bem,
meg=True, eeg=False, mindist=5.0, n_jobs=2)
mne.write_forward_solution(fwdfile, fwd)
fwd
[Parallel(n_jobs=2)]: Using backend LokyBackend with 2 concurrent workers.
[Parallel(n_jobs=2)]: Done 2 out of 2 | elapsed: 1.5s finished
[Parallel(n_jobs=2)]: Using backend LokyBackend with 2 concurrent workers.
[Parallel(n_jobs=2)]: Done 2 out of 2 | elapsed: 1.0s finished
[Parallel(n_jobs=2)]: Using backend LokyBackend with 2 concurrent workers.
[Parallel(n_jobs=2)]: Done 2 out of 2 | elapsed: 1.0s finished
[Parallel(n_jobs=2)]: Using backend LokyBackend with 2 concurrent workers.
[Parallel(n_jobs=2)]: Done 2 out of 2 | elapsed: 0.9s finished
[Parallel(n_jobs=2)]: Using backend LokyBackend with 2 concurrent workers.
[Parallel(n_jobs=2)]: Done 2 out of 2 | elapsed: 0.9s finished
[Parallel(n_jobs=2)]: Using backend LokyBackend with 2 concurrent workers.
[Parallel(n_jobs=2)]: Done 2 out of 2 | elapsed: 0.9s finished
[Parallel(n_jobs=2)]: Using backend LokyBackend with 2 concurrent workers.
[Parallel(n_jobs=2)]: Done 2 out of 2 | elapsed: 0.6s finished
Extract the fixed orientation lead field matrix:
fwd_fixed = mne.convert_forward_solution(
fwd, surf_ori=True, force_fixed=True, use_cps=True)
# leadfield matrix
lf = eelbrain.load.fiff.forward_operator(fwd_fixed, src='ico-4', subjects_dir=subjects_dir)
NCRF estimation¶
Now that we have all the required data to estimate NCRFs.
Note
This example uses simplified settings to speed up estimation:
1) For this example, we use a fixed regularization parameter (mu).
For a real experiment, the optimal mu would be determined by
cross-validation (set mu='auto', which is the default).
The optimal mu will then be stored in model.mu
(this is how the mu used here was determined).
2) The example forces the estimation to stop after fewer iterations than
is recommended (n_iter). For stable models, we recommend to use the
default setting (n_iter=10).
# To speed up the example, we cache the NCRF:
ncrf_file = data_path / 'MEG' / 'bst_auditory' / 'oddball_ncrf.pickle'
if ncrf_file.exists():
model = eelbrain.load.unpickle(ncrf_file)
else:
model = fit_ncrf(
meg, [stim1, stim2], lf, noise_cov, tstart=0, tstop=0.5,
mu=0.0001756774187547859, n_iter=5,
)
eelbrain.save.pickle(model, ncrf_file)
The learned kernel/filter (the NCRF) can be accessed as an attribute of the
model.
NCRFs are stored as eelbrain.NDVar. Here, the two NCRFs correspond
to the two different predictor variables:
[<NDVar: 5107 source, 51 time>, <NDVar: 5107 source, 51 time>]
Visualization¶
A butterfly plot shows weights in all sources over time. This is good for forming a quick impression of important time lags, or peaks in the response:
Note
Since the estimates are sparse over cortical locations, smoothing the NCRFs over sources to make the visualization more intuitive.
The following code for plotting the anatomical localization is commented because the Mayavi based plots do not work reliably in the automatic documentation. Uncomment it to create anatomical plots.
A single time point can be visualized with the PySurfer (surfer)
based eelbrain.plot.brain.brain():
# brain = eelbrain.plot.brain.brain(h[0].sub(time=0.140), vmax=2e-11, surf='pial')
An eelbrain.plot.brain.SequencePlotter can be used to plot a
sequence of brain images, for example in a jupyter notebook:
# h_binned = h0.bin(0.1, 0.1, 0.4, 'extrema')
# sp = eelbrain.plot.brain.SequencePlotter()
# sp.set_brain_args(surf='inflated')
# sp.add_ndvar(h_binned)
# p = sp.plot_table(view='lateral')
In an interactive iPython session, we can also use interactive time-linked
plots with eelbrain.plot.brain.butterfly():
# brain, butterfly = eelbrain.plot.brain.butterfly(h0)
Total running time of the script: (11 minutes 16.436 seconds)