说明
转到末尾下载完整示例代码。
基于回归的基线校正#
本教程比较传统基线校正（对 epoch 中的每个时间点加上或减去同一个标量值）与基于回归的基线校正方法（允许基线时段的效应随时间点变化）。具体而言，本教程遵循 Alday [ 1 ] 提出的方法。
至少有两点理由值得考虑使用基于回归的基线校正：
与传统基线校正不同，基于回归的方法不假定基线的影响在不同实验条件之间是等同的。因此它在避免引入偏差方面更为安全。
假设试次前基线信号水平主要由数据中的慢漂移决定，那么在时间上离基线时段越远，信号水平的振幅与基线振幅相近的可能性就越小。因此，使用随时间变化的基线校正，在长时程 epoch 的较后时段引入信号失真或虚假效应的可能性更低。
无论是传统基线校正还是基于回归的基线校正，都存在一个共同问题：应选择哪个时间窗作为基线窗。
# Authors: Carina Forster
# Email: carinaforster0611@gmail.com
# License: BSD-3-Clause
# Copyright the MNE-Python contributors.
import numpy as np
import mne
加载数据#
我们首先加载 MNE-Python 示例数据集，并提取实验事件以获得试次位置和试次类型。由于本教程只关注 EEG 通道，我们可以丢弃其他通道类型以加快速度：
data_path = mne.datasets.sample.data_path()
raw_fname = data_path / "MEG" / "sample" / "sample_audvis_filt-0-40_raw.fif"
raw = mne.io.read_raw_fif(raw_fname, preload=True)
events = mne.find_events(raw)
raw.pick(picks="eeg")
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...
Finding events on: STI 014
319 events found on stim channel STI 014
Event IDs: [ 1 2 3 4 5 32]
这里我们把来自左右半球的视觉和听觉事件合并，并构建用于分段的 event_id 字典。
events = mne.merge_events(events, [1, 2], 1) # auditory events will be "1"
events = mne.merge_events(events, [3, 4], 2) # visual events will be "2"
event_id = {"auditory": 1, "visual": 2}
预处理#
接下来定义用于分段和预处理数据的一些变量。由于要合并左右两侧的刺激，因此我们将观察单个中央电极，以可视化听觉与视觉试次之间的差异。
我们进行一些标准预处理（带通滤波），然后对数据分段。注意我们没有对 epoch 做基线校正（指定 baseline=None）；我们只是通过剔除振幅过高或过低的通道来做最基本的清理。另请注意，我们在数据的副本上操作，以便之后与传统的基线校正方法进行比较。
raw_filtered = raw.copy().filter(highpass, lowpass)
epochs = mne.Epochs(
raw_filtered,
events,
event_id,
tmin=tmin,
tmax=tmax,
reject=dict(eeg=150e-6),
flat=dict(eeg=5e-6),
baseline=None,
preload=True,
)
Filtering raw data in 1 contiguous segment
Setting up band-pass filter from 0.1 - 40 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.10
- Lower transition bandwidth: 0.10 Hz (-6 dB cutoff frequency: 0.05 Hz)
- Upper passband edge: 40.00 Hz
- Upper transition bandwidth: 10.00 Hz (-6 dB cutoff frequency: 45.00 Hz)
- Filter length: 4957 samples (33.013 s)
Not setting metadata
288 matching events found
No baseline correction applied
Created an SSP operator (subspace dimension = 1)
1 projection items activated
Using data from preloaded Raw for 288 events and 106 original time points ...
Rejecting epoch based on EEG : ['EEG 001', 'EEG 002', 'EEG 003', 'EEG 007']
1 bad epochs dropped
传统基线校正#
首先用传统方式对数据做基线校正。我们对每个条件内的 epoch 求平均，并分别对听觉和视觉试次减去各自条件特异的基线。
Applying baseline correction (mode: mean)
Applying baseline correction (mode: mean)
基于回归的基线校正#
现在尝试基于回归的基线校正方法。我们将使用 mne.stats.linear_regression()，它需要一个设计矩阵来表示回归预测变量。我们将使用四个预测变量：每个实验条件各一个、一个表示基线效应、另一个表示基线与其中一个条件之间的交互作用（以考虑基线效应在两个条件之间可能存在的异质性）。下面是前两个：
aud_predictor = epochs.events[:, 2] == epochs.event_id["auditory"]
vis_predictor = epochs.events[:, 2] == epochs.event_id["visual"]
基线预测变量的计算稍复杂：我们将在基线期内分别对每个 epoch 求平均值，并将该值作为（试次级）预测变量。这里由于我们只关注一个特定通道，将使用该通道的基线值作为预测变量；但依据你的研究问题，你可能希望对每个通道分别进行此操作，或跨通道整合信息。
baseline_predictor = (
epochs.copy()
.crop(*baseline)
.pick([ch])
.get_data(copy=False) # convert to NumPy array
.mean(axis=-1) # average across timepoints
.squeeze() # only 1 channel, so remove singleton dimension
)
baseline_predictor *= 1e6 # convert V → μV
注意我们只把预测变量（而非 epochs 数据）从 Volts 转换为 microVolts。这样做是为了回归模型拟合（极小的数值可能使模型拟合不稳定）。
现在可以构建设计矩阵：把一维预测变量按行堆叠，再用 .T 转置使其成为列。将它们组合进一个 array() 中也会自动把布尔型 aud_predictor 和 vis_predictor 转换为 1 和 0：
最后，我们拟合回归模型：
reg_model = mne.stats.linear_regression(
epochs, design_matrix, names=["auditory", "visual", "baseline", "baseline:visual"]
)
Fitting linear model to epochs, (6360 targets, 4 regressors)
Done
该函数返回一个由 mne.stats.regression.lm 对象构成的字典，每个对象都是一个 namedtuple()，其中存储的各种估计值可以像 Evoked 对象一样访问。让我们检查一下：
dict_keys(['auditory', 'visual', 'baseline', 'baseline:visual'])
model attributes: ('beta', 'stderr', 't_val', 'p_val', 'mlog10_p_val')
values are stored in Evoked objects:
<Evoked | '' (average, N=1), -0.1998 – 0.49949 s, baseline off, 60 ch, ~3.0 MiB>
绘制基线回归量#
首先看基线时段的估计效应。我们关心的是 beta 值，它告诉我们基线值在每个时间点上的预测强度。模型会为每个通道估计其有效性，但由于我们只用一个通道构建了基线预测变量，我们只看该通道的结果。我们将在 β=1 处添加一条水平线，代表传统基线校正——它假定该效应在各时间点上恒定：
effect_of_baseline = reg_model["baseline"].beta
effect_of_baseline.plot(
picks=ch,
hline=[1.0],
units=dict(eeg=r"$\beta$ value"),
titles=dict(eeg=ch),
selectable=False,
)
不出所料，趋势是：离基线时段越远，基线振幅的预测价值越弱。换句话说，（在此数据中）较早的时间点应比之后的时间点接受更强的基线校正。
绘制 ERP#
现在看两个条件（auditory 和 visual）的 beta 值：这些系数表示在考虑了基线（随时间变化的！）效应之后，实验刺激对信号的“纯”影响。我们把它们绘制在一起，并与传统基线方法并排比较：
reg_aud = reg_model["auditory"].beta
reg_vis = reg_model["visual"].beta
kwargs = dict(picks=ch, show_sensors=False, truncate_yaxis=False)
mne.viz.plot_compare_evokeds(
dict(auditory=trad_aud, visual=trad_vis), title="Traditional", **kwargs
)
mne.viz.plot_compare_evokeds(
dict(auditory=reg_aud, visual=reg_vis), title="Regression-based", **kwargs
)
它们看起来相当相似，但两个条件之间的差距在一些细节上存在差异（例如在 400-500 ms 附近）。
绘制头皮地形图与差异波#
现在比较传统方法与基于回归方法的头皮地形图。我们通过计算条件之间的差异来实现：
diff_traditional = mne.combine_evoked([trad_aud, trad_vis], weights=[1, -1])
diff_regression = mne.combine_evoked([reg_aud, reg_vis], weights=[1, -1])
绘图之前，先确保两幅图使用相同的颜色刻度：
vmin = min(diff_traditional.get_data().min(), diff_regression.get_data().min()) * 1e6
vmax = max(diff_traditional.get_data().max(), diff_regression.get_data().max()) * 1e6
topo_kwargs = dict(vlim=(vmin, vmax), ch_type="eeg", times=np.linspace(0.05, 0.45, 9))
fig = diff_traditional.plot_topomap(**topo_kwargs)
fig.suptitle("Traditional")
fig = diff_regression.plot_topomap(**topo_kwargs)
fig.suptitle("Regression-based")
可以看到，基于回归的方法在早期（约 100-150 ms）显示出条件间更强的差异，而在后期（约 250-350 ms，以及 450 ms 附近）差异更弱。这一点也反映在差异波本身：注意基于回归的差异波在 150 ms 附近离零更远，而在 250-350 ms 附近更接近零。
title = "Difference in evoked potential (auditory minus visual)"
fig = mne.viz.plot_compare_evokeds(
dict(Traditional=diff_traditional, Regression=diff_regression),
title=title,
**kwargs,
)
考察交互项#
最后看回归模型中的交互项。它告诉我们基线时段的效应在视觉试次中是否不同于其在听觉试次中的效应。这里我们在零处添加一条水平线，表示“应当没有差异”的假设（即基线在某一类试次中不应系统性地更高，且基线效应在两类试次中随时间持续的时间也不应有差异）。
interaction_effect = reg_model["baseline:visual"].beta
interaction_effect.plot(
picks=ch,
hline=[0.0],
units=dict(eeg=r"$\beta$ value"),
titles=dict(eeg=ch),
selectable=False,
)
确实，交互项的 beta 权重相当小，且似乎在零附近随机波动，提示基线效应在我们的两类试次之间没有系统性差异。
参考文献#
脚本总运行时间：（0 分钟 3.038 秒）