Python协方差矩阵处理脑电数据

在本教程中，我们将介绍传感器协方差计算的基础知识，并构建一个噪声协方差矩阵，该矩阵可用于计算最小范数逆解.

诸如MNE的源估计方法需要从记录中进行噪声估计。  

在本教程中，我们介绍了噪声协方差的基础知识，并构造了一个噪声协方差矩阵，该矩阵可在计算逆解时使用。

下面我们将结合代码来进行分析。

# 导入工具库
import os.path as op

import mne
from mne.datasets import sample
import matplotlib.pyplot as plt

读取数据

# 构建数据地址
data_path = sample.data_path()
raw_empty_room_fname = op.join(
    data_path, 'MEG', 'sample', 'ernoise_raw.fif')
raw_empty_room = mne.io.read_raw_fif(raw_empty_room_fname)
raw_fname = op.join(data_path, 'MEG', 'sample', 'sample_audvis_raw.fif')
# 读取原始数据
raw = mne.io.read_raw_fif(raw_fname)
#设置平均参考
raw.set_eeg_reference('average', projection=True)
#增加一个不良通道
raw.info['bads'] += ['EEG 053']  # bads + 1 more

噪声的定义取决于范式。在MEG中，通常使用空房间测量来估计传感器噪声。但是，如果要处理诱发的反应，可能还需要考虑将静息状态的大脑活动视为噪音。首先，我们使用空房间记录来计算噪声。请注意，您还可以仅将录制的一部分与tmin和tmax参数一起使用。如果您将静息状态用作噪声基线，这这将非常有用。在这里，我们使用整个空房间记录来计算噪声协方差。

请记住，在操作时要在处理方面将空房间数据集与实际的MEG数据进行匹配。确保过滤器都相同，并且如果使用ICA，则将其等效地应用于空房间和主题数据。在这种情况下，我们没有过滤数据，也没有使用ICA。但是，我们在MEG数据中确实存在错误的通道和投影，因此，我们要确保将它们存储在协方差对象中。

raw_empty_room.info['bads'] = [
    bb for bb in raw.info['bads'] if 'EEG' not in bb]
raw_empty_room.add_proj(
    [pp.copy() for pp in raw.info['projs'] if 'EEG' not in pp['desc']])
# 进行协方差均在计算
noise_cov = mne.compute_raw_covariance(
    raw_empty_room, tmin=0, tmax=None)

3 projection items deactivated
Using up to 550 segments
Number of samples used : 66000
[done]

现在，已经在MNE-Python对象中有了协方差矩阵，可以使用mne.write_cov()将其保存到文件中。  

稍后，可以使用mne.read_cov()将其读回。

还可以使用刺激前的基线来估计噪声协方差。  

首先，我们必须构建epoch。  

计算协方差时，应该在构建epochs时使用基线校正。否则协方差矩阵将不准确。 

在MNE中，默认情况下会完成此操作，但为了确定，我们在此处手动定义它。

# 在原始数据中找到events
events = mne.find_events(raw)
# 构建epoch
epochs = mne.Epochs(raw, events, event_id=1, tmin=-0.2, tmax=0.5,
                    baseline=(-0.2, 0.0), decim=3,  # we'll decimate for speed
                    verbose='error')  # and ignore the warning about aliasing

320 events found
Event IDs: [ 1 2 3 4 5 32]

"""
对epochs数据进行协方差估计
"""
noise_cov_baseline = mne.compute_covariance(epochs, tmax=0)

绘制协方差矩阵

尝试将proj设置为False以查看效果。
请注意，epochs中的投影机已经应用，因此proj参数无效。

noise_cov.plot(raw_empty_room.info, proj=True)

plt.show()

noise_cov_baseline.plot(epochs.info, proj=True)
plt.show()

应该如何规范协方差矩阵？

估计的协方差可能在数值上不稳定，并且倾向于在估计的源振幅和可用样本数之间引起相关性。 

因此，MNE手册建议对噪声协方差矩阵进行正则化（请参阅对噪声协方差矩阵进行正则化），尤其是在只有少量样本可用的情况下。 

然而，要说出样本的有效数量并不容易，因此要选择适当的正则化。  

在MNE-Python中，使用[1]中所述的高级正则化方法来完成正则化。为此，可以使用'auto'选项。使用此选项，交叉验证将用于学习最佳正则化：

noise_cov_reg = mne.compute_covariance(epochs, tmax=0., method='auto',
                                       rank=None)

此过程使用看不见的数据的负对数似然来量化噪声协方差，从而量化噪声协方差。最终结果也可以用肉眼检查。在假设基线不包含系统信号(对感兴趣的事件进行时间锁定)的情况下，白化基线信号应遵循多元高斯分布，即，在给定的条件下，白化基线信号应在-1.96和1.96之间 时间样本。基于相同的推理，全局场功率(GFP)的期望值为1(计算GFP时应考虑真实的自由度，如ddof=3，其中有2个活动的SSP向量):

evoked = epochs.average()
evoked.plot_white(noise_cov_reg, time_unit='s')
plt.show()

该图同时显示了每个通道的白化诱发信号和白化的GFP。GFP面板中的数字代表数据的估计等级，相当于计算白色GFP时跨传感器的平方和除以的有效自由度。变白的GFP还有助于检测虚假的晚期诱发成分，这可能是由于过度或欠正则化所致。

请注意，如果使用信号空间分离(SSS) 2处理数据，则会同时显示梯度仪和磁力仪，因为两者都是由相同的SSS基向量以相同的数值秩重建的。这也意味着这两种传感器类型不再是统计上独立的。这些评估方法可用于评估模型违规。

对于专家用例或调试，还可以将替代估计量进行比较，并演示白化对源估计的影响：

noise_covs = mne.compute_covariance(
    epochs, tmax=0., method=('empirical', 'shrunk'), return_estimators=True,
    rank=None)
evoked.plot_white(noise_covs, time_unit='s')
plt.show()

这将绘制出为最佳估计量而诱发的白化并在相关面板中将所有估计器的GFP显示为单独的行。

最后，让我们看一下空房间和与事件相关的协方差之间的区别，使用"method"选项，使它们的类型显示在图例中。

evoked_meg = evoked.copy().pick('meg')
noise_cov['method'] = 'empty_room'
noise_cov_baseline['method'] = 'baseline'
evoked_meg.plot_white([noise_cov_baseline, noise_cov], time_unit='s')
plt.show()

基于对数可能性的负值，基线协方差似乎更合适。

参考

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.

2.Taulu, S., Simola, J., Kajola, M., 2005. Applications of the signal space separation method. IEEE Trans. Signal Proc. 53, 3359-3372.

---

仅用于学术交流，不用于商业行为，若有侵权及疑问，请后台留言，管理员即时删侵！

更多阅读

---

---

点个在看祝你开心一整天!