听觉P300实验与视觉P300相似,但使用听觉刺激来产生oddball
刺激时间为200ms,时间间隔400ms,随机抖动±100ms, 任务是计算玩奇数球刺激的次数,记录单个参与者进行的6次2分钟的实验。
非目标刺激的音调约为523Hz(C5),目标刺激的音调约为1174 Hz(D6)
下面利用muse库来进行数据分析
第一步:导入相应的包
代码语言:javascript复制import os
import sys
from collections import OrderedDict
import pandas as pd
import numpy as np
import seaborn as sns
from matplotlib import pyplot as plt
%matplotlib inline
from eeg_utils.utils import utils as utils
import warnings
warnings.filterwarnings("ignore", category=Warning)
第二步:加载数据
代码语言:javascript复制subject = 1
session = 1
raw = utils.load_data('auditory/P300', sfreq=256.,
subject_nb=subject, session_nb=session)
查看功率谱(Power Spectrum)
代码语言:javascript复制raw.plot_psd(tmax=np.inf);
第三步:过滤
过滤掉1到30Hz之间的数据
代码语言:javascript复制raw.filter(1, 30, method='iir')
Epoching
我们在刺激之后将数据从-100ms缩短到800ms,无需进行基线校正(信号经过带通滤波),并且在信号超过75uV时,我们拒绝每个epochs,主要是拒绝眨眼。
代码语言:javascript复制from mne import Epochs, find_events
events = find_events(raw)
event_id = {'Non-Target': 1, 'Target': 2}
epochs = Epochs(raw, events=events, event_id=event_id,
tmin=-0.1, tmax=0.8, baseline=None,
reject={'eeg': 75e-6}, preload=True,
verbose=False, picks=[0,1,2,3]
Epoch average
我们绘制两种情况下的平均ERP
代码语言:javascript复制conditions = OrderedDict()
conditions['Non-target'] = [1]
conditions['Target'] = [2]
fig, ax = utils.plot_conditions(epochs, conditions=conditions,
ci=97.5, n_boot=1000, title='',
diff_waveform=(1, 2))
在两种情况下都可以看到清晰的P300
解码
通过上述的平均epochs,可以很清楚识别ERP。 但如何了解有关P300的SNR的任何信息,可以通过分类管道(classification pipline)了解P300响应的强度。
下面我们将使用4个不同的管道。
Vect RL: 实验矢量化 Logistic回归。可以将其视为MEG/EEG中的标准解码管道。
Vect RegLDA:实验矢量化 正则化LDA.P300 BCI中经常利用这一点,但维数过多,则可能无法使用。
ERPCov MDM: ERP协方差 MDM.黎曼几何分类器,一种简单有效的方法(用于较少的通道数)
ERPCov TS: ERP协方差 切线空间映射。这是一种基于黎曼几何的管道之一。
以AUC作为度量,以交叉验证的方式进行评估(AUC可能是针对二进制和非平衡分类问题的最佳度量标准)
代码语言:javascript复制from sklearn.pipeline import make_pipeline
from mne.decoding import Vectorizer
from sklearn.linear_model import LogisticRegression
from sklearn.preprocessing import StandardScaler
from sklearn.discriminant_analysis import LinearDiscriminantAnalysis as LDA
from sklearn.model_selection import cross_val_score, StratifiedShuffleSplit
from pyriemann.estimation import ERPCovariances, XdawnCovariances
from pyriemann.tangentspace import TangentSpace
from pyriemann.classification import MDM
clfs = OrderedDict()
clfs['Vect LR'] = make_pipeline(Vectorizer(), StandardScaler(), LogisticRegression())
clfs['Vect RegLDA'] = make_pipeline(Vectorizer(), LDA(shrinkage='auto', solver='eigen'))
clfs['ERPCov TS'] = make_pipeline(ERPCovariances(estimator='oas'), TangentSpace(), LogisticRegression())
clfs['ERPCov MDM'] = make_pipeline(ERPCovariances(estimator='oas'), MDM())
clfs['XdawnCov TS'] = make_pipeline(XdawnCovariances(estimator='oas'), TangentSpace(), LogisticRegression())
clfs['XdawnCov MDM'] = make_pipeline(XdawnCovariances(estimator='oas'), MDM())
# format data
epochs.pick_types(eeg=True)
X = epochs.get_data() * 1e6
times = epochs.times
y = epochs.events[:, -1]
# 定义交叉验证
cv = StratifiedShuffleSplit(n_splits=20, test_size=0.25,
random_state=42)
# 对每个管道(pipeline)运行交叉验证
auc = []
methods = []
for m in clfs:
try:
res = cross_val_score(clfs[m], X, y==2, scoring='roc_auc',
cv=cv, n_jobs=-1)
auc.extend(res)
methods.extend([m]*len(res))
except:
pass
results = pd.DataFrame(data=auc, columns=['AUC'])
results['Method'] = methods
代码语言:javascript复制fig = plt.figure(figsize=[8,4])
sns.barplot(data=results, x='AUC', y='Method')
plt.xlim(0.4, 0.75)
sns.despine()