本教程為腦機學習者Rose發表於公眾號：腦機接口社區(微信號：Brain_Computer).QQ交流群：903290195
案例介紹
本案例通過對多導睡眠圖（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()，可以方便地從Sleep 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進行多分類
下面展示了解決如何從愛麗絲的數據中預測鮑勃的睡眠階段並盡可能避免重復樣板代碼的問題。
這里將利用sckit-learn的Pipeline和FunctionTransformer。
擴展：[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))
預測的准確精度為84.7%
查看分類報告做進一步分析
print(classification_report(y_test, y_pred, target_names=event_id.keys()))
從分類報告中可以看出，
Bob的每個階段訓練測試樣本，以及對應的睡眠階段的精度。
比如W階段的精度為86%,測試樣本為1856。測試總樣本為2802。也可以看到其他一些指標比如召回率和F1值。這些指標的含義以后介紹。
腦機接口 QQ交流群：903290195
更多分享，請關注公眾號