案例介绍 本案例对多导睡眠图(Polysomnography,PSG)数据进行睡眠阶段的分类来判断睡眠类型 训练:对Alice的睡眠数据进行训练 测试:利用训练结果对Bob的睡眠数据进行测试,判断其睡眠类型 在分析之前,先简单介绍一下多导睡眠图,多导睡眠图(Polysomnography,PSG)又称睡眠脑电图。主要用于睡眠和梦境研究以及抑郁症和呼吸暂停综合征的诊断 多导睡眠图是通过不同部位的生物电或通过不同传感获得生物讯号,经前置放大,输出为不同的电讯号,记录出不同的图形以供分析。
数据集介绍: 本案例用的数据是来自于PhysioNet上关于健康受试者的年龄对睡眠影响研究的公开数据集的一个子集。 mne.datasets.sleep_physionet.age.fetch_data可以下载PhysioNet数据集的子数据集 该子数据集中包含20位受试者的实验数据,记录当时年龄为25-34岁的10位男性和10位女性的实验数据。 由于受试者13的第二个记录遗失了,所以除了受试者13以外,每个受试者都有两次夜间记录
Sleep Physionet数据集使用8个标签进行标注,代表8个阶段: Wake(W), Stage 1, Stage 2, Stage 3, Stage 4, REM(R), Movement time(M), Stage(?).Stage(?)-(not scored) 唤醒-Wake(W)、第1阶段、第2阶段、第3阶段、第4阶段、对应于轻度睡眠到深度睡眠的范围;REM睡眠(R),其中REM是Rapid Eye的缩写,表示快速眼运动睡眠,运动(M)和阶段(?)的任何未得分部分
#导入工具库 import numpy as np import matplotlib.pyplot as plt import mne from mne.datasets.sleep_physionet.age import fetch_data from mne.time_frequency import psd_welch from sklearn.ensemble import RandomForestClassifier from sklearn.metrics import accuracy_score from sklearn.metrics import confusion_matrix from sklearn.metrics import classification_report from sklearn.pipeline import make_pipeline from sklearn.preprocessing import FunctionTransformer第二步:加载数据 在这里,我们从两个主题下载数据,最终目标是获得时间片段(epochs)及其相关的地面真理
MNE-Python为我们提供了mne.datasets.sleep_physionet.age.fetch_data(),可以方便地从Slepp Physionet数据集下载数据,给定主题和记录的列表,提取程序将下载数据并为每个主题提供数据,一对文件: -PSG.edf 包含多导睡眠图。来自EEG头盔的原始数据 -Hypnogram.edf 包含专家记录的注释
然后,将这两个对象合并到mne.io.Raw对象中,就可以根据注释的描述提取事件以获得时间片段(epochs)
可以通过 mne.datasets.sleep_physionet.age.fetch_data(subjects,recording,path) 来获取PhysioNet多导睡眠图数据集文件。
subjects:表示想要使用哪些受试者对象,可供选择的受试者对象范围0-19 recording:表示夜间记录的编号(索引),有效值为[1]、[2]或[1、2]
path:PhysioNet数据的存放地址,如果没有给定,则加载默认存放数据的地址: 如果默认存放数据集的地址不存在数据,则从网络中下载相关数据。
#选择两个受试者实验对象ALICE,BOB(改名字并非实验中的真实名,这里是为了方便临时取的名字) ALICE,BOB=0,1 #加载ALICE,BOB的实验数据文件 [alice_files,bob_files]=fetch_data(subjects=[ALICE,BOB],recording=[1]) #通道名称映射 mapping={'EOG horizontal':'eog', 'Resp oro-nasal':'misc', 'EMG submental':'misc', 'Temp rectal':'misc', 'Event marker':'misc'} #读取ALICE的edf文件,和其对应的注释文件 raw_train=mne.io.read_raw_edf(alice_files[0]) annot_train=mne.read_annotations(alice_files[1]) raw_train.set_annotations(annot_train,emit_warning=False) raw_train.set_channel_types(mapping) #绘制空0s开始,时间窗口长度为40s的连续通道数据波形图 raw_train.plot(duration=40,scalings='auto') plt.show()输出结果: 不知道为什么,我这个图出现了这样子的问题,希望有大佬可以指点一下 这里仅使用5个阶段:唤醒(W),阶段1,阶段2,阶段3/4和REM睡眠(R)。 为此,这里使用mne.events_from_annotations()中的event_id参数来选择我们感兴趣的事件,并将事件标识符与每个事件相关联。
''' 睡眠表示与事件映射 ''' annotation_desc_2_event_id={'Sleep stage W':1, 'Sleep stage 1':2, 'Sleep stage 2':3, 'Sleep stage 3':4, 'Sleep stage 4':4, 'Sleep stage R':5} events_train,_=mne.events_from_annotations(raw_train,event_id=annotation_desc_2_event_id, chunk_duration=30) #创建一个新的event_id以统一阶段3和4 event_id={'Sleep stage W':1, 'Sleep stage 1':2, 'Sleep stage 2':3, 'Sleep stage 3/4':4, 'Sleep stage R':5} #绘制事件数据 mne.viz.plot_events(events_train,event_id=event_id, sfreq=raw_train.info['sfreq']) #保留颜色代码以便进一步绘制 stage_colors=plt.rcParams['axes.prop_cycle'].by_key()['color'] #根据注释中的事件从数据创建epochs(时间片段) tmax=30.-1./raw_train.info['sfreq']#tmax in included ''' 所创建的是时间从tmin=0开始,到tmax为止的epochs ''' epochs_train=mne.Epochs(raw=raw_train,events=events_train, event_id=event_id,tmin=0,tmax=tmax,baseline=None) print(epochs_train) ###第三步:加载Bob的数据作为测试数据 #按照上述相同的步骤来获取Bob的测试数据 raw_test=mne.io.read_raw_edf(bob_files[0]) annot_test=mne.read_annotations(bob_files[1]) raw_test.set_annotations(annot_test,emit_warning=False) raw_test.set_channel_types(mapping) events_test,_=mne.events_from_annotations(raw_test,event_id=annotation_desc_2_event_id, chunk_duration=30.) epochs_test=mne.Epochs(raw=raw_test,events=events_test,event_id=event_id, tmin=0.,tmax=tmax,baseline=None) print(epochs_test)特征工程 观察不同睡眠阶段的功率谱密度(PSD)图,可以看到不同睡眠阶段具有不同的特征。这些签名在Alice和Bob的数据中保持相似
在本节的其余部分中,将基于特定频带的相对功率来创建EEG特征,以捕获数据中睡眠阶段之间的差异
fig,(ax1,ax2)=plt.subplots(ncols=2) #iterate over the subjects stages=sorted(event_id.keys()) for ax,title,epochs in zip([ax1,ax2], ['Alice','Bob'], [epochs_train,epochs_test]): for stage,color in zip(stages,stage_colors): epochs[stage].plot_psd(area_mode=None,color=color,ax=ax, fmin=0.1,fmax=20.,show=False, average=True,spatial_colors=False) ax.set(title=title,xlabel='Frequency(Hz)') ax2.set(ylabel='uv^2/hz(dB)') ax2.legend(ax2.lines[2::3],stages) plt.tight_layout() plt.show()第四步:设计scikit-learn转换器 创建一个函数,根据特定频带中的相对功率提取脑电图特征,从而能够根据脑电图信号预测睡眠阶段
def eeg_power_band(epochs): ''' 脑电相对功率带特征提取 该函数接受一个‘mne.Epochs’对象, 并基于与scikit-learn兼容的特定频带中的相对功率创建EEG特征 Parameters ------------ epochs:Epochs The data. Returns ------------ X:numpy array of shape[n_samples,5] Transformed data. ''' #特定频带 FREQ_BANDS={'delta':[0.5,4.5], 'theta':[4.5,8.5], 'alpha':[8.5,11.5], 'sigma':[11.5,15.5], 'beta':[15.5,30]} psds,freqs=psd_welch(epochs,picks='eeg',fmin=0.5,fmax=30.) #归一化 PSDs psds/=np.sum(psds,axis=-1,keepdims=True) X=[] for fmin,fmax in FREQ_BANDS.values(): psds_band=psds[:,:,(freqs>=fmin)&(freqs<fmax)].mean(axis=-1) X.append(psds_band.reshape(len(psds),-1)) return np.concatenate(X,axis=1)第五步:根据Alice的数据来预测Bob的睡眠阶段 使用scikit-learn进行多分类 下面展示了解决如何从爱丽丝的数据中预测Bob的睡眠阶段并尽可能避免重复样板代码的问题,这里将利用scikit-learn的Pipeline和Function Transformer 扩展:[Pipeline可以将许多算法模型串联起来,可以用于把多个estamitors级联成一个estamitor,比如将特征提取、归一化、分类组织在一起形成一个典型的机器学习问题工作流。 FunctionTransformer将python函数转换为与estamitor兼容的对象。]
pipe=make_pipeline(FunctionTransformer(eeg_power_band,validate=False), RandomForestClassifier(n_estimators=100,random_state=42)) #训练 y_train=epochs_train.events[:,2] pipe.fit(epochs_train,y_train) #预测 y_pred=pipe.predict(epochs_test) #评估准确率 y_test=epochs_test.events[:,2] acc=accuracy_score(y_test,y_pred) print('Accuracy score:{}'.format(acc))准确率达到了85.3%
####查看分类报告做进一步分析 print(classification_report(y_test,y_pred,target_names=event_id.keys()))学起来还是有一些懵懵懂懂,继续加油!
