Python機(jī)器學(xué)習(xí)預(yù)測(cè)心血管疾病
出品:CDA數(shù)據(jù)分析師(cdacdacda)
大家好,歡迎來到 Crossin的編程教室 !
全球每年約有1700萬人死于心血管疾病,當(dāng)中主要表現(xiàn)為心肌梗死和心力衰竭。當(dāng)心臟不能泵出足夠的血液來滿足人體的需要時(shí),就會(huì)發(fā)生心力衰竭,通常由糖尿病、高血壓或其他心臟疾病引起。
在檢測(cè)心血管疾病的早期癥狀時(shí),機(jī)器學(xué)習(xí)就能派上用場(chǎng)了。通過患者的電子病歷,可以記錄患者的癥狀、身體特征、臨床實(shí)驗(yàn)室測(cè)試值,從而進(jìn)行生物統(tǒng)計(jì)分析,這能夠發(fā)現(xiàn)那些醫(yī)生無法檢測(cè)到的模式和相關(guān)性。
尤其通過機(jī)器學(xué)習(xí),根據(jù)數(shù)據(jù)就能預(yù)測(cè)患者的存活率,今天的文章我們分享就是如何用Python寫一個(gè)心血管疾病的預(yù)測(cè)模型。
我們用到的數(shù)據(jù)集來自Davide Chicco和Giuseppe Jurman發(fā)表的論文:《機(jī)器學(xué)習(xí)可以僅通過血肌酐和射血分?jǐn)?shù)來預(yù)測(cè)心力衰竭患者的生存率》。
他們收集整理了299名心力衰竭患者的醫(yī)療記錄,這些患者數(shù)據(jù)來自2015年4月至12月間巴基斯坦費(fèi)薩拉巴德心臟病研究所和費(fèi)薩拉巴德聯(lián)合醫(yī)院。這些患者由105名女性和194名男性組成,年齡在40至95歲之間。所有299例患者均患有左心室收縮功能不全,并曾出現(xiàn)過心力衰竭。
Davide和Giuseppe應(yīng)用了多個(gè)機(jī)器學(xué)習(xí)分類器來預(yù)測(cè)患者的生存率,并根據(jù)最重要的危險(xiǎn)因素對(duì)特征進(jìn)行排序。同時(shí)還利用傳統(tǒng)的生物統(tǒng)計(jì)學(xué)測(cè)試進(jìn)行了另一種特征排序分析,并將這些結(jié)果與機(jī)器學(xué)習(xí)算法提供的結(jié)果進(jìn)行比較。
他們分析對(duì)比了心力衰竭患者的一系列數(shù)據(jù),最終發(fā)現(xiàn)根據(jù)血肌酐和射血分?jǐn)?shù)這兩項(xiàng)數(shù)據(jù)能夠很好的預(yù)測(cè)心力衰竭患者的存活率。
今天我們就教教大家,如果根據(jù)這共13個(gè)字段的299 條病人診斷記錄,用Python寫出預(yù)測(cè)心力衰竭患者存活率的預(yù)測(cè)模型。
下面是具體的步驟和關(guān)鍵代碼。
01
數(shù)據(jù)理解
數(shù)據(jù)取自于kaggle平臺(tái)分享的心血管疾病數(shù)據(jù)集,共有13個(gè)字段299 條病人診斷記錄。具體的字段概要如下:

02
數(shù)據(jù)讀入和初步處理
首先導(dǎo)入所需包。
# 數(shù)據(jù)整理
import numpy as np
import pandas as pd
# 可視化
import matplotlib.pyplot as plt
import seaborn as sns
import plotly as py
import plotly.graph_objs as go
import plotly.express as px
import plotly.figure_factory as ff
# 模型建立
from sklearn.linear_model import LogisticRegression
from sklearn.tree import DecisionTreeClassifier
from sklearn.ensemble import GradientBoostingClassifier, RandomForestClassifier
import lightgbm
# 前處理
from sklearn.preprocessing import StandardScaler
# 模型評(píng)估
from sklearn.model_selection import train_test_split, GridSearchCV
from sklearn.metrics import plot_confusion_matrix, confusion_matrix, f1_score
加載并預(yù)覽數(shù)據(jù)集:
# 讀入數(shù)據(jù)
df = pd.read_csv('./data/heart_failure.csv')
df.head()

03
探索性分析
1. 描述性分析
df.describe().T

從上述描述性分析結(jié)果簡(jiǎn)單總結(jié)如下:
是否死亡:平均的死亡率為32%; 年齡分布:平均年齡60歲,最小40歲,最大95歲 是否有糖尿?。河?1.8%患有糖尿病 是否有高血壓:有35.1%患有高血壓 是否抽煙:有32.1%有抽煙
2. 目標(biāo)變量
# 產(chǎn)生數(shù)據(jù)
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='目標(biāo)變量DEATH_EVENT的分布')
py.offline.plot(fig, filename='./html/目標(biāo)變量DEATH_EVENT的分布.html')
總共有299人,其中隨訪期未存活人數(shù)96人,占總?cè)藬?shù)的32.1%
3. 貧血

從圖中可以看出,有貧血癥狀的患者死亡概率較高,為35.66%。
bar1 = draw_categorical_graph(df['anaemia'], df['DEATH_EVENT'], title='紅細(xì)胞、血紅蛋白減少和是否存活')
bar1.render('./html/紅細(xì)胞血紅蛋白減少和是否存活.html')
4. 年齡

從直方圖可以看出,在患心血管疾病的病人中年齡分布差異較大,表現(xiàn)趨勢(shì)為年齡越大,生存比例越低、死亡的比例越高。
# 產(chǎn)生數(shù)據(jù)
surv = df[df['DEATH_EVENT'] == 0]['age']
not_surv = df[df['DEATH_EVENT'] == 1]['age']
hist_data = [surv, not_surv]
group_labels = ['Survived', 'Not Survived']
# 直方圖
fig = ff.create_distplot(hist_data, group_labels, bin_size=0.5)
fig.update_layout(title_text='年齡和生存狀態(tài)關(guān)系')
py.offline.plot(fig, filename='./html/年齡和生存狀態(tài)關(guān)系.html')
5. 年齡/性別

從分組統(tǒng)計(jì)和圖形可以看出,不同性別之間生存狀態(tài)沒有顯著性差異。在死亡的病例中,男性的平均年齡相對(duì)較高。
6. 年齡/抽煙

數(shù)據(jù)顯示,整體來看,是否抽煙與生存與否沒有顯著相關(guān)性。但是當(dāng)我們關(guān)注抽煙的人群中,年齡在50歲以下生存概率較高。
7. 磷酸肌酸激酶(CPK)

從直方圖可以看出,血液中CPK酶的水平較高的人群死亡的概率較高。
8. 射血分?jǐn)?shù)

射血分?jǐn)?shù)代表了心臟的泵血功能,過高和過低水平下,生存的概率較低。
9. 血小板

血液中血小板(100~300)×10^9個(gè)/L,較高或較低的水平則代表不正常,存活的概率較低。
10. 血肌酐水平

血肌酐是檢測(cè)腎功能的最常用指標(biāo),較高的指數(shù)代表腎功能不全、腎衰竭,有較高的概率死亡。
11. 血清鈉水平

圖形顯示,血清鈉較高或較低往往伴隨著風(fēng)險(xiǎn)。
12. 相關(guān)性分析

從數(shù)值型屬性的相關(guān)性圖可以看出,變量之間沒有顯著的共線性關(guān)系。
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('Pearson correlation coefficient between numeric variables', fontdict={'fontsize': 15})
plt.show()
04
特征篩選
我們使用統(tǒng)計(jì)方法進(jìn)行特征篩選,目標(biāo)變量DEATH_EVENT是分類變量時(shí),當(dāng)自變量是分類變量,使用卡方鑒定,自變量是數(shù)值型變量,使用方差分析。
# 劃分X和y
X = df.drop('DEATH_EVENT', axis=1)
y = df['DEATH_EVENT']
from feature_selection import Feature_select
fs = Feature_select(num_method='anova', cate_method='kf')
X_selected = fs.fit_transform(X, y)
X_selected.head()
2020 17:19:49 INFO attr select success!
After select attr: ['serum_creatinine', 'serum_sodium', 'ejection_fraction', 'age', 'time']

05
數(shù)據(jù)建模
首先劃分訓(xùn)練集和測(cè)試集。
# 劃分訓(xùn)練集和測(cè)試集
Features = X_selected.columns
X = 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)
# 標(biāo)準(zhǔn)化
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-score
print("F1_score of LogisticRegression is : ", round(f1_score(y_true=y_test, y_pred=test_pred),2))
我們使用決策樹進(jìn)行建模,設(shè)置特征選擇標(biāo)準(zhǔn)為gini,樹的深度為5。輸出混淆矩陣圖:在這個(gè)案例中,1類是我們關(guān)注的對(duì)象。
# DecisionTreeClassifier
clf = DecisionTreeClassifier(criterion='gini', max_depth=5, random_state=1)
clf.fit(X_train, y_train)
test_pred = clf.predict(X_test)
# F1-score
print("F1_score of DecisionTreeClassifier is : ", 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 - Confusion Matrix", fontsize=15)
plt.xticks(range(2), ["Heart Not Failed","Heart Fail"], fontsize=12)
plt.yticks(range(2), ["Heart Not Failed","Heart Fail"], fontsize=12)
plt.show()
F1_score of DecisionTreeClassifier is : 0.61
<Figure size 720x504 with 0 Axes>

使用網(wǎng)格搜索進(jìn)行參數(shù)調(diào)優(yōu),優(yōu)化標(biāo)準(zhǔn)為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
使用最優(yōu)的模型重新評(píng)估測(cè)試集效果:
test_pred = GS.best_estimator_.predict(X_test)
# F1-score
print("F1_score of DecisionTreeClassifier is : ", 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 - Confusion Matrix", fontsize=15)
plt.xticks(range(2), ["Heart Not Failed","Heart Fail"], fontsize=12)
plt.yticks(range(2), ["Heart Not Failed","Heart Fail"], fontsize=12)
plt.show()

使用隨機(jī)森林
# RandomForestClassifier
rfc = 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-score
print("F1_score of RandomForestClassifier is : ", round(f1_score(y_true=y_test, y_pred=test_pred),2))
{'max_depth': 3}
0.791157747481277
F1_score of RandomForestClassifier is : 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_)
# 測(cè)試集
test_pred = GS.best_estimator_.predict(X_test)
# F1-score
print("F1_score of GradientBoostingClassifier is : ", round(f1_score(y_true=y_test, y_pred=test_pred),2))
{'max_depth': 3}
0.7288420428900305
F1_score of GradientBoostingClassifier is : 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_)
# 測(cè)試集
test_pred = GS.best_estimator_.predict(X_test)
# F1-score
print("F1_score of LGBMClassifier is : ", round(f1_score(y_true=y_test, y_pred=test_pred),2))
{'max_depth': 2}
0.780378102289867
F1_score of LGBMClassifier is : 0.74
以下為各模型在測(cè)試集上的表現(xiàn)效果對(duì)比:
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
獲取詳細(xì)數(shù)據(jù)及代碼,請(qǐng)?jiān)诠娞?hào)對(duì)話內(nèi)回復(fù)關(guān)鍵字 心血管
如果文章對(duì)你有幫助,歡迎轉(zhuǎn)發(fā)/點(diǎn)贊/收藏~
_往期文章推薦_
