计算效应匹配空间滤波（EMS）#

本示例计算 EMS，以重建 [1] 中描述的实验效应时间过程。

该技术基于两个条件之间的差异构建空间滤波器。把试次投影到相应的空间滤波器上，可得到替代性的单试次数据：多传感器活动被压缩为一条时间序列，若存在实验效应便会显现出来。

首先绘制单试次的“试次 × 时间”图像，并按条件排列试次；第二幅图显示各条件的平均时间序列；最后绘制地形图，展示空间滤波器随时间的演变。

# Author: Denis Engemann <denis.engemann@gmail.com>
#         Jean-Rémi King <jeanremi.king@gmail.com>
#
# License: BSD-3-Clause
# Copyright the MNE-Python contributors.

import matplotlib.pyplot as plt
import numpy as np
from sklearn.model_selection import StratifiedKFold

import mne
from mne import EvokedArray, io
from mne.datasets import sample
from mne.decoding import EMS, compute_ems

print(__doc__)

data_path = sample.data_path()

# Preprocess the data
meg_path = data_path / "MEG" / "sample"
raw_fname = meg_path / "sample_audvis_filt-0-40_raw.fif"
event_fname = meg_path / "sample_audvis_filt-0-40_raw-eve.fif"
event_ids = {"AudL": 1, "VisL": 3}

# Read data and create epochs
raw = io.read_raw_fif(raw_fname, preload=True)
raw.filter(0.5, 45, fir_design="firwin")
events = mne.read_events(event_fname)

raw.pick(["grad", "eog"], exclude="bads")

epochs = mne.Epochs(
    raw,
    events,
    event_ids,
    tmin=-0.2,
    tmax=0.5,
    baseline=None,
    reject=dict(grad=4000e-13, eog=150e-6),
    preload=True,
)
epochs.pick("grad")

# Setup the data to use it a scikit-learn way:
X = epochs.get_data(copy=False)  # The MEG data
y = epochs.events[:, 2]  # The conditions indices
n_epochs, n_channels, n_times = X.shape

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.
Reading 0 ... 41699  =      0.000 ...   277.709 secs...
Filtering raw data in 1 contiguous segment
Setting up band-pass filter from 0.5 - 45 Hz

FIR filter parameters
---------------------
Designing a one-pass, zero-phase, non-causal bandpass filter:
- Windowed time-domain design (firwin) method
- Hamming window with 0.0194 passband ripple and 53 dB stopband attenuation
- Lower passband edge: 0.50
- Lower transition bandwidth: 0.50 Hz (-6 dB cutoff frequency: 0.25 Hz)
- Upper passband edge: 45.00 Hz
- Upper transition bandwidth: 11.25 Hz (-6 dB cutoff frequency: 50.62 Hz)
- Filter length: 993 samples (6.613 s)

Not setting metadata
145 matching events found
No baseline correction applied
0 projection items activated
Using data from preloaded Raw for 145 events and 106 original time points ...
    Rejecting  epoch based on EOG : ['EOG 061']
    Rejecting  epoch based on EOG : ['EOG 061']
    Rejecting  epoch based on EOG : ['EOG 061']
    Rejecting  epoch based on EOG : ['EOG 061']
    Rejecting  epoch based on EOG : ['EOG 061']
    Rejecting  epoch based on EOG : ['EOG 061']
    Rejecting  epoch based on EOG : ['EOG 061']
    Rejecting  epoch based on EOG : ['EOG 061']
    Rejecting  epoch based on EOG : ['EOG 061']
    Rejecting  epoch based on EOG : ['EOG 061']
    Rejecting  epoch based on EOG : ['EOG 061']
    Rejecting  epoch based on EOG : ['EOG 061']
    Rejecting  epoch based on EOG : ['EOG 061']
    Rejecting  epoch based on EOG : ['EOG 061']
    Rejecting  epoch based on EOG : ['EOG 061']
    Rejecting  epoch based on EOG : ['EOG 061']
    Rejecting  epoch based on EOG : ['EOG 061']
    Rejecting  epoch based on EOG : ['EOG 061']
    Rejecting  epoch based on EOG : ['EOG 061']
    Rejecting  epoch based on EOG : ['EOG 061']
    Rejecting  epoch based on EOG : ['EOG 061']
    Rejecting  epoch based on EOG : ['EOG 061']
22 bad epochs dropped

# Initialize EMS transformer
ems = EMS()

# Initialize the variables of interest
X_transform = np.zeros((n_epochs, n_times))  # Data after EMS transformation
filters = list()  # Spatial filters at each time point

# In the original paper, the cross-validation is a leave-one-out. However,
# we recommend using a Stratified KFold, because leave-one-out tends
# to overfit and cannot be used to estimate the variance of the
# prediction within a given fold.

for train, test in StratifiedKFold(n_splits=5).split(X, y):
    # In the original paper, the z-scoring is applied outside the CV.
    # However, we recommend to apply this preprocessing inside the CV.
    # Note that such scaling should be done separately for each channels if the
    # data contains multiple channel types.
    X_scaled = X / np.std(X[train])

    # Fit and store the spatial filters
    ems.fit(X_scaled[train], y[train])

    # Store filters for future plotting
    filters.append(ems.filters_)

    # Generate the transformed data
    X_transform[test] = ems.transform(X_scaled[test])

# Average the spatial filters across folds
filters = np.mean(filters, axis=0)

# Plot individual trials
plt.figure()
plt.title("single trial surrogates")
plt.imshow(
    X_transform[y.argsort()],
    origin="lower",
    aspect="auto",
    extent=[epochs.times[0], epochs.times[-1], 1, len(X_transform)],
    cmap="RdBu_r",
)
plt.xlabel("Time (ms)")
plt.ylabel("Trials (reordered by condition)")

# Plot average response
plt.figure()
plt.title("Average EMS signal")
mappings = [(key, value) for key, value in event_ids.items()]
for key, value in mappings:
    ems_ave = X_transform[y == value]
    plt.plot(epochs.times, ems_ave.mean(0), label=key)
plt.xlabel("Time (ms)")
plt.ylabel("a.u.")
plt.legend(loc="best")
plt.show()

# Visualize spatial filters across time
evoked = EvokedArray(filters, epochs.info, tmin=epochs.tmin)
evoked.plot_topomap(scalings=1)

注意，compute_ems 也能施加类似的变换。但该函数复现的是 Schurger 等人原论文的做法，因此把归一化放在留一交叉验证之外，我们不建议这样使用。

epochs.equalize_event_counts(event_ids)
X_transform, filters, classes = compute_ems(epochs)

Dropped 11 epochs: 26, 47, 66, 75, 76, 79, 80, 81, 100, 103, 108
...computing surrogate time series. This can take some time

参考文献#

脚本总运行时间：（0 分 1.296 秒）

图库由 Sphinx-Gallery 生成