全球每年約有1700萬人死於心血管疾病，當中主要表現為心肌梗死和心力衰竭。當心臟不能泵出足夠的血液來滿足人體的需要時，就會發生心力衰竭，通常由糖尿病、高血壓或其他心臟疾病引起。
在檢測心血管疾病的早期症狀時，機器學習就能派上用場了。通過患者的電子病歷，可以記錄患者的症狀、身體特徵、臨床實驗室測試值，從而進行生物統計分析，這能夠發現那些醫生無法檢測到的模式和相關性。
尤其通過機器學習，根據數據就能預測患者的存活率，今天我們就教大家如何用Python寫一個心血管疾病的預測模型。
我們用到的數據集來自Davide Chicco和Giuseppe Jurman發表的論文：《機器學習可以僅通過血肌酐和射血分數來預測心力衰竭患者的生存率》。
他們收集整理了299名心力衰竭患者的醫療記錄，這些患者數據來自2015年4月至12月間巴基斯坦費薩拉巴德心臟病研究所和費薩拉巴德聯合醫院。這些患者由105名女性和194名男性組成，年齡在40至95歲之間。所有299例患者均患有左心室收縮功能不全，並曾出現過心力衰竭。
Davide和Giuseppe應用了多個機器學習分類器來預測患者的生存率，並根據最重要的危險因素對特徵進行排序。同時還利用傳統的生物統計學測試進行了另一種特徵排序分析，並將這些結果與機器學習算法提供的結果進行比較。
他們分析對比了心力衰竭患者的一系列數據，最終發現根據血肌酐和射血分數這兩項數據能夠很好的預測心力衰竭患者的存活率。
今天我們就教教大家，如果根據這共13個字段的299 條病人診斷記錄，用Python寫出預測心力衰竭患者存活率的預測模型。
下面是具體的步驟和關鍵代碼。
數據理解
數據取自於kaggle平台分享的心血管疾病數據集，共有13個字段299 條病人診斷記錄。具體的字段概要如下：
數據讀入和初步處理
首先導入所需包。
#數據整理importnumpyasnpimportpandasaspd#可視化importmatplotlib.pyplotaspltimportseabornassnsimportplotlyaspyimportplotly.graph_objsasgoimportplotly.expressaspximportplotly.figure_factoryasff#模型建立fromsklearn.linear_modelimportLogisticRegressionfromsklearn.treeimportDecisionTreeClassifierfromsklearn.ensembleimportGradientBoostingClassifier,RandomForestClassifierimportlightgbm#前處理fromsklearn.preprocessingimportStandardScaler#模型評估fromsklearn.model_selectionimporttrain_test_split,GridSearchCVfromsklearn.metricsimportplot_confusion_matrix,confusion_matrix,f1_score加載並預覽數據集：
#讀入數據df=pd.read_csv('./data/heart_failure.csv')df.head()
探索性分析
1. 描述性分析
從上述描述性分析結果簡單總結如下：
2. 目標變量
#產生數據death_num=df['DEATH_EVENT'].value_counts()death_num=death_num.reset_index()#餅圖fig=px.pie(death_num,names='index',values='DEATH_EVENT')fig.update_layout(title_text='目標變量DEATH_EVENT的分布')py.offline.plot(fig,filename='./html/目標變量DEATH_EVENT的分布.html')總共有299人，其中隨訪期未存活人數96人，占總人數的32.1%
3. 貧血
從圖中可以看出，有貧血症狀的患者死亡概率較高，為35.66%。
bar1=draw_categorical_graph(df['anaemia'],df['DEATH_EVENT'],title='紅細胞、血紅蛋白減少和是否存活')bar1.render('./html/紅細胞血紅蛋白減少和是否存活.html')4. 年齡
從直方圖可以看出，在患心血管疾病的病人中年齡分布差異較大，表現趨勢為年齡越大，生存比例越低、死亡的比例越高。
#產生數據surv=df[df['DEATH_EVENT']==0]['age']not_surv=df[df['DEATH_EVENT']==1]['age']hist_data=[surv,not_surv]group_labels=['Survived','NotSurvived']#直方圖fig=ff.create_distplot(hist_data,group_labels,bin_size=0.5)fig.update_layout(title_text='年齡和生存狀態關係')py.offline.plot(fig,filename='./html/年齡和生存狀態關係.html')5. 年齡/性別
從分組統計和圖形可以看出，不同性別之間生存狀態沒有顯著性差異。在死亡的病例中，男性的平均年齡相對較高。
6. 年齡/抽煙
數據顯示，整體來看，是否抽煙與生存與否沒有顯著相關性。但是當我們關注抽煙的人群中，年齡在50歲以下生存概率較高。
7. 磷酸肌酸激酶(CPK)
從直方圖可以看出，血液中CPK酶的水平較高的人群死亡的概率較高。
8. 射血分數
射血分數代表了心臟的泵血功能，過高和過低水平下，生存的概率較低。
9. 血小板
血液中血小板（100~300）×10^9個/L，較高或較低的水平則代表不正常，存活的概率較低。
10. 血肌酐水平
血肌酐是檢測腎功能的最常用指標，較高的指數代表腎功能不全、腎衰竭，有較高的概率死亡。
11. 血清鈉水平
圖形顯示，血清鈉較高或較低往往伴隨着風險。
12. 相關性分析
從數值型屬性的相關性圖可以看出，變量之間沒有顯著的共線性關係。
num_df=df[['age','creatinine_phosphokinase','ejection_fraction','platelets','serum_creatinine','serum_sodium']]plt.figure(figsize=(12,12))sns.heatmap(num_df.corr(),vmin=-1,cmap='coolwarm',linewidths=0.1,annot=True)plt.title('Pearsoncorrelationcoefficientbetweennumericvariables',fontdict={'fontsize':15})plt.show()特徵篩選
我們使用統計方法進行特徵篩選，目標變量DEATH_EVENT是分類變量時，當自變量是分類變量，使用卡方鑑定，自變量是數值型變量，使用方差分析。
#劃分X和yX=df.drop('DEATH_EVENT',axis=1)y=df['DEATH_EVENT']fromfeature_selectionimportFeature_selectfs=Feature_select(num_method='anova',cate_method='kf')X_selected=fs.fit_transform(X,y)X_selected.head()202017:19:49INFOattrselectsuccess!Afterselectattr:['serum_creatinine','serum_sodium','ejection_fraction','age','time']
數據建模
首先劃分訓練集和測試集。
#劃分訓練集和測試集Features=X_selected.columnsX=df[Features]y=df["DEATH_EVENT"]X_train,X_test,y_train,y_test=train_test_split(X,y,test_size=0.2,stratify=y,random_state=2020)#標準化scaler=StandardScaler()scaler_Xtrain=scaler.fit_transform(X_train)scaler_Xtest=scaler.fit_transform(X_test)lr=LogisticRegression()lr.fit(scaler_Xtrain,y_train)test_pred=lr.predict(scaler_Xtest)#F1-scoreprint("F1_scoreofLogisticRegressionis:",round(f1_score(y_true=y_test,y_pred=test_pred),2))我們使用決策樹進行建模，設置特徵選擇標準為gini，樹的深度為5。輸出混淆矩陣圖：在這個案例中，1類是我們關注的對象。
#DecisionTreeClassifierclf=DecisionTreeClassifier(criterion='gini',max_depth=5,random_state=1)clf.fit(X_train,y_train)test_pred=clf.predict(X_test)#F1-scoreprint("F1_scoreofDecisionTreeClassifieris:",round(f1_score(y_true=y_test,y_pred=test_pred),2))#繪圖plt.figure(figsize=(10,7))plot_confusion_matrix(clf,X_test,y_test,cmap='Blues')plt.title("DecisionTreeClassifier-ConfusionMatrix",fontsize=15)plt.xticks(range(2),["HeartNotFailed","HeartFail"],fontsize=12)plt.yticks(range(2),["HeartNotFailed","HeartFail"],fontsize=12)plt.show()F1_scoreofDecisionTreeClassifieris:0.61<Figuresize720x504with0Axes>
使用網格搜索進行參數調優，優化標準為f1。
parameters={'splitter':('best','random'),'criterion':("gini","entropy"),"max_depth":[*range(1,20)],}clf=DecisionTreeClassifier(random_state=1)GS=GridSearchCV(clf,param_grid=parameters,cv=10,scoring='f1',n_jobs=-1)GS.fit(X_train,y_train)print(GS.best_params_)print(GS.best_score_){'criterion':'entropy','max_depth':3,'splitter':'best'}0.7638956305132776使用最優的模型重新評估測試集效果：
test_pred=GS.best_estimator_.predict(X_test)#F1-scoreprint("F1_scoreofDecisionTreeClassifieris:",round(f1_score(y_true=y_test,y_pred=test_pred),2))#繪圖plt.figure(figsize=(10,7))plot_confusion_matrix(GS,X_test,y_test,cmap='Blues')plt.title("DecisionTreeClassifier-ConfusionMatrix",fontsize=15)plt.xticks(range(2),["HeartNotFailed","HeartFail"],fontsize=12)plt.yticks(range(2),["HeartNotFailed","HeartFail"],fontsize=12)plt.show()
使用隨機森林
#RandomForestClassifierrfc=RandomForestClassifier(n_estimators=1000,random_state=1)parameters={'max_depth':np.arange(2,20,1)}GS=GridSearchCV(rfc,param_grid=parameters,cv=10,scoring='f1',n_jobs=-1)GS.fit(X_train,y_train)print(GS.best_params_)print(GS.best_score_)test_pred=GS.best_estimator_.predict(X_test)#F1-scoreprint("F1_scoreofRandomForestClassifieris:",round(f1_score(y_true=y_test,y_pred=test_pred),2)){'max_depth':3}0.791157747481277F1_scoreofRandomForestClassifieris:0.53使用Boosting
gbl=GradientBoostingClassifier(n_estimators=1000,random_state=1)parameters={'max_depth':np.arange(2,20,1)}GS=GridSearchCV(gbl,param_grid=parameters,cv=10,scoring='f1',n_jobs=-1)GS.fit(X_train,y_train)print(GS.best_params_)print(GS.best_score_)#測試集test_pred=GS.best_estimator_.predict(X_test)#F1-scoreprint("F1_scoreofGradientBoostingClassifieris:",round(f1_score(y_true=y_test,y_pred=test_pred),2)){'max_depth':3}0.7288420428900305F1_scoreofGradientBoostingClassifieris:0.65使用LGBMClassifier
lgb_clf=lightgbm.LGBMClassifier(boosting_type='gbdt',random_state=1)parameters={'max_depth':np.arange(2,20,1)}GS=GridSearchCV(lgb_clf,param_grid=parameters,cv=10,scoring='f1',n_jobs=-1)GS.fit(X_train,y_train)print(GS.best_params_)print(GS.best_score_)#測試集test_pred=GS.best_estimator_.predict(X_test)#F1-scoreprint("F1_scoreofLGBMClassifieris:",round(f1_score(y_true=y_test,y_pred=test_pred),2)){'max_depth':2}0.780378102289867F1_scoreofLGBMClassifieris:0.74以下為各模型在測試集上的表現效果對比：
LogisticRegression：0.63
DecisionTree Classifier：0.73
Random Forest Classifier: 0.53
GradientBoosting Classifier: 0.65
LGBM Classifier: 0.74
參考鏈接：
Machine learning can predict survival of patients with heart failure from serum creatinine and ejection fraction alone
https://bmcmedinformdecismak.biomedcentral.com/articles/10.1186/s12911-020-1023-5#Abs1
留言列表