基于随机森林模型的心脏病人预测分类

公众号：尤而小屋

作者：Peter

编辑：Peter

大家好，我是 Peter~

今天给大家分享一个新的 kaggle 案例：基于随机森林模型（RandomForest）的心脏病人预测分类。本文涉及到的知识点主要包含：

• 数据预处理和类型转化

• 随机森林模型建立与解释

• 决策树的可视化

• 部分依赖图 PDP 的绘制和解释

• AutoML 机器学习 SHAP 库的使用和解释（个人待提升）

导读

Of all the applications of machine-learning, diagnosing any serious disease using a black box is always going to be a hard sell. If the output from a model is the particular course of treatment (potentially with side-effects), or surgery, or the absence of treatment, people are going to want to know why.

在机器学习的所有应用中，使用黑匣子诊断任何严重疾病总是很难的。如果模型的输出是特定的治疗过程（可能有副作用）、手术或是否有疗效，人们会想知道为什么。

This dataset gives a number of variables along with a target condition of having or not having heart disease. Below, the data is first used in a simple random forest model, and then the model is investigated using ML explainability tools and techniques.

该数据集提供了许多变量以及患有或不患有心脏病的目标条件。下面，数据首先用于一个简单的随机森林模型，然后使用 ML 可解释性工具和技术对该模型进行研究。

感兴趣的可以参考原文学习，notebook 地址为：https://www.kaggle.com/tentotheminus9/what-causes-heart-disease-explaining-the-model

导入库

本案例中涉及到多个不同方向的库：

• 数据预处理

• 多种可视化绘图；尤其是 shap 的可视化，模型可解释性的使用（后面会专门写这个库）

• 随机森林模型

• 模型评价等

import numpy as npimport pandas as pdimport matplotlib.pyplot as pltimport seaborn as sns from sklearn.ensemble import RandomForestClassifier from sklearn.tree import DecisionTreeClassifierfrom sklearn.tree import export_graphviz from sklearn.metrics import roc_curve, auc from sklearn.metrics import classification_report from sklearn.metrics import confusion_matrix from sklearn.model_selection import train_test_split import eli5 from eli5.sklearn import PermutationImportanceimport shap from pdpbox import pdp, info_plots np.random.seed(123) 
pd.options.mode.chained_assignment = None  

数据探索 EDA

1、导入数据

2、缺失值情况

数据比较完美，没有任何缺失值！

字段含义

在这里重点介绍下各个字段的含义。Peter 近期导出的数据集中的额字段和原 notebook 中的字段名字写法稍有差异（时间原因导致），还好 Peter 已经为大家做了一一对应的关系，下面是具体的中文含义：

1. age：年龄

2. sex 性别 1=male 0=female

3. cp 胸痛类型；4 种取值情况

4. 1：典型心绞痛

5. 2：非典型心绞痛

6. 3：非心绞痛

7. 4：无症状

8. trestbps 静息血压

9. chol 血清胆固醇

10. fbs 空腹血糖 >120mg/dl ：1=true； 0=false

11. restecg 静息心电图(值 0,1,2)

12. thalach 达到的最大心率

13. exang 运动诱发的心绞痛(1=yes;0=no)

14. oldpeak 相对于休息的运动引起的 ST 值(ST 值与心电图上的位置有关)

15. slope 运动高峰 ST 段的坡度

16. 1：upsloping 向上倾斜

17. 2： flat 持平

18. 3：downsloping 向下倾斜

19. ca The number of major vessels(血管) (0-3)

20. thal A blood disorder called thalassemia ，一种叫做地中海贫血的血液疾病(3 = normal；6 = fixed defect;；7 = reversable defect)

21. target 生病没有(0=no；1=yes)

原 notebook 中的英文含义；

下面是 Peter 整理的对应关系。本文中以当前的版本为标准：

字段转化

转化编码

对部分字段进行一一的转化。以 sex 字段为例：将数据中的 0 变成 female，1 变成 male

# 1、sex
df["sex"][df["sex"] == 0] = "female"df["sex"][df["sex"] == 1] = "male"

字段类型转化

# 指定数据类型df["sex"] = df["sex"].astype("object")df["cp"] = df["cp"].astype("object")df["fbs"] = df["fbs"].astype("object")df["restecg"] = df["restecg"].astype("object")df["exang"] = df["exang"].astype("object")df["slope"] = df["slope"].astype("object")df["thal"] = df["thal"].astype("object")

生成哑变量

# 生成哑变量df = pd.get_dummies(df,drop_first=True)df

随机森林 RandomForest

切分数据

# 生成特征变量数据集和因变量数据集X = df.drop("target",1)y = df["target"]
# 切分比例为8：2X_train, X_test, y_train, y_test = train_test_split(X,y,test_size=0.2,random_state=10)
X_train

建模

rf = RandomForestClassifier(max_depth=5)rf.fit(X_train, y_train)

3 个重要属性

随机森林中 3 个重要的属性：

• 查看森林中树的状况：estimators_

• 袋外估计准确率得分：oob_score_，必须是oob_score参数选择 True 的时候才可用

• 变量的重要性：feature_importances_

决策树可视化

在这里我们选择的第二棵树的可视化过程：

# 查看第二棵树的状况estimator = rf.estimators_[1]
# 全部属性feature_names = [i for i in X_train.columns]#print(feature_names)

# 指定数据类型y_train_str = y_train.astype('str')# 0-no 1-diseasey_train_str[y_train_str == '0'] = 'no disease'y_train_str[y_train_str == '1'] = 'disease'# 训练数据的取值y_train_str = y_train_str.valuesy_train_str[:5]

绘图的具体代码为：

# 绘图显示
export_graphviz(    estimator,   # 传入第二颗树    out_file='tree.dot',   # 导出文件名    feature_names = feature_names,  # 属性名    class_names = y_train_str,  # 最终的分类数据    rounded = True,     proportion = True,     label='root',    precision = 2,     filled = True)
from subprocess import callcall(['dot', '-Tpng', 'tree.dot', '-o', 'tree.png', '-Gdpi=600'])
from IPython.display import ImageImage(filename = 'tree.png')

决策树的可视化过程能够让我们看到具体的分类过程，但是并不能解决哪些特征或者属性比较重要。后面会对部分属性的特征重要性进行探索

模型得分验证

关于混淆矩阵和使用特异性（specificity）以及灵敏度（sensitivity）这两个指标来描述分类器的性能：

# 模型预测y_predict = rf.predict(X_test)y_pred_quant = rf.predict_proba(X_test)[:,1]y_pred_bin = rf.predict(X_test)
# 混淆矩阵confusion_matrix = confusion_matrix(y_test,y_pred_bin)confusion_matrix
# 计算sensitivity and specificity total=sum(sum(confusion_matrix))sensitivity = confusion_matrix[0,0]/(confusion_matrix[0,0]+confusion_matrix[1,0])specificity = confusion_matrix[1,1]/(confusion_matrix[1,1]+confusion_matrix[0,1])

绘制 ROC 曲线

fpr, tpr, thresholds = roc_curve(y_test, y_pred_quant)
fig, ax = plt.subplots()ax.plot(fpr, tpr)
ax.plot([0,1],[0,1],        transform = ax.transAxes,        ls = "--",        c = ".3"       )
plt.xlim([0.0, 1.0])plt.ylim([0.0, 1.0])
plt.rcParams['font.size'] = 12
# 标题plt.title('ROC Curve')# 两个轴的名称plt.xlabel('False Positive Rate (1 - Specificity)')plt.ylabel('True Positive Rate (Sensitivity)')# 网格线plt.grid(True)

本案例中的 ROC 曲线值：

auc(fpr, tpr)# 结果0.9076923076923078

根据一般 ROC 曲线的评价标准，案例的表现结果还是不错的：

• 0.90 - 1.00 = excellent

• 0.80 - 0.90 = good

• 0.70 - 0.80 = fair

• 0.60 - 0.70 = poor

• 0.50 - 0.60 = fail

补充知识点：分类器的评价指标

考虑一个二分类的情况，类别为 1 和 0，我们将 1 和 0 分别作为正类（positive）和负类（negative），根据实际的结果和预测的结果，则最终的结果有 4 种，表格如下：

常见的评价指标：

1、ACC：classification accuracy，描述分类器的分类准确率

计算公式为：ACC=(TP+TN)/(TP+FP+FN+TN)

2、BER：balanced error rate 计算公式为：BER=1/2*(FPR+FN/(FN+TP))

3、TPR：true positive rate，描述识别出的所有正例占所有正例的比例计算公式为：TPR=TP/ (TP+ FN)

4、FPR：false positive rate，描述将负例识别为正例的情况占所有负例的比例计算公式为：FPR= FP / (FP + TN)

5、TNR：true negative rate，描述识别出的负例占所有负例的比例计算公式为：TNR= TN / (FP + TN)

6、PPV：Positive predictive value计算公式为：PPV=TP / (TP + FP)

7、NPV：Negative predictive value计算公式：NPV=TN / (FN + TN)

其中 TPR 即为敏感度（sensitivity），TNR 即为特异度（specificity）。

来自维基百科的经典图形：

可解释性

排列重要性-Permutation Importance

下面的内容是关于机器学习模型的结果可解释性。首先考察的是每个变量对模型的重要性。重点考量的排列重要性 Permutation Importance：

部分依赖图（ Partial dependence plots ，PDP）

一维 PDP

Partial Dependence 就是用来解释某个特征和目标值 y 的关系的，一般是通过画出 Partial Dependence Plot(PDP)来体现。 也就是说 PDP 在 X1 的值，就是把训练集中第一个变量换成 X1 之后，原模型预测出来的平均值。

重点：查看单个特征和目标值的关系

字段 ca

base_features = df.columns.values.tolist()base_features.remove("target")
feat_name = 'ca'  # ca-num_major_vessels 原文pdp_dist = pdp.pdp_isolate(    model=rf,  # 模型    dataset=X_test,  # 测试集    model_features=base_features,  # 特征变量；除去目标值     feature=feat_name  # 指定单个字段)
pdp.pdp_plot(pdp_dist, feat_name)  # 传入两个参数plt.show()

通过下面的图形我们观察到：当 ca 字段增加的时候，患病的几率在下降。ca 字段的含义是血管数量（num_major_vessels），也就是说当血管数量增加的时候，患病率随之降低

字段 age

feat_name = 'age'
pdp_dist = pdp.pdp_isolate(    model=rf,     dataset=X_test,     model_features=base_features,     feature=feat_name)
pdp.pdp_plot(pdp_dist, feat_name)plt.show()

关于年龄字段 age，原文的描述：

That's a bit odd. The higher the age, the lower the chance of heart disease? Althought the blue confidence regions show that this might not be true (the red baseline is within the blue zone).

翻译：这有点奇怪。年龄越大，患心脏病的几率越低？尽管蓝色置信区间表明这可能不是真的（红色基线在蓝色区域内）

字段 oldpeak

feat_name = 'oldpeak'
pdp_dist = pdp.pdp_isolate(    model=rf,     dataset=X_test,     model_features=base_features,     feature=feat_name)
pdp.pdp_plot(pdp_dist, feat_name)plt.show()

oldpeak 字段同样表明：取值越大，患病几率越低。

这个变量称之为“相对休息运动引起的 ST 压低值”。正常的状态下，该值越高，患病几率越高。但是上面的图像却显示了相反的结果。

作者推断：造成这个结果的原因除了 the depression amount，可能还和 slope type 有关系。原文摘录如下，于是作者绘制了 2D-PDP 图形

Perhaps it's not just the depression amount that's important, but the interaction with the slope type? Let's check with a 2D PDP

2D-PDP 图

查看的是 slope_upsloping 、slope_flat 和 oldpeak 的关系：

inter1  =  pdp.pdp_interact(    model=rf,  # 模型    dataset=X_test,  # 特征数据集    model_features=base_features,  # 特征    features=['slope_upsloping', 'oldpeak'])
pdp.pdp_interact_plot(    pdp_interact_out=inter1,     feature_names=['slope_upsloping', 'oldpeak'],     plot_type='contour')plt.show()
## ------------
inter1  =  pdp.pdp_interact(    model=rf,     dataset=X_test,    model_features=base_features,     features=['slope_flat', 'oldpeak'])
pdp.pdp_interact_plot(    pdp_interact_out=inter1,     feature_names=['slope_flat', 'oldpeak'],     plot_type='contour')plt.show()

从两张图形中我们可以观察到：在 oldpeak 取值较低的时候，患病几率都比较高（黄色），这是一个奇怪的现象。于是作者进行了如下的 SHAP 可视化探索：针对单个变量进行分析。

SHAP 可视化

关于 SHAP 的介绍可以参考文章：https://zhuanlan.zhihu.com/p/83412330 和 https://blog.csdn.net/sinat_26917383/article/details/115400327

SHAP 是 Python 开发的一个"模型解释"包，可以解释任何机器学习模型的输出。下面 SHAP 使用的部分功能：

Explainer

在 SHAP 中进行模型解释之前需要先创建一个 explainer，SHAP 支持很多类型的 explainer，例如 deep, gradient, kernel, linear, tree, sampling)。在这个案例我们以 tree 为例：

# 传入随机森林模型rfexplainer = shap.TreeExplainer(rf)  # 在explainer中传入特征值的数据，计算shap值shap_values = explainer.shap_values(X_test)  shap_values

Feature Importance

取每个特征 SHAP 值的绝对值的平均值作为该特征的重要性，得到一个标准的条形图(multi-class 则生成堆叠的条形图：

结论：能够直观地观察到 ca 字段的 SHAP 值最高

summary_plot

summary plot 为每个样本绘制其每个特征的 SHAP 值，这可以更好地理解整体模式，并允许发现预测异常值。

• 每一行代表一个特征，横坐标为 SHAP 值

• 一个点代表一个样本，颜色表示特征值的高低(红色高，蓝色低)

个体差异

查看单个病人的不同特征属性对其结果的影响，原文描述为：

Next, let's pick out individual patients and see how the different variables are affecting their outcomes

def heart_disease_risk_factors(model, patient):
    explainer = shap.TreeExplainer(model)  # 建立explainer    shap_values = explainer.shap_values(patient)  # 计算shape值    shap.initjs()    return shap.force_plot(      explainer.expected_value[1],      shap_values[1],       patient)

从两个病人的结果中显示：

• P1：预测准确率高达 29%（baseline 是 57%），更多的因素集中在 ca、thal_fixed_defect、oldpeak 等蓝色部分。

• P3：预测准确率高达 82%，更多的影响因素在 sel_male=0，thalach=143 等

通过对比不同的患者，我们是可以观察到不同病人之间的预测率和主要影响因素。

dependence_plot

为了理解单个 feature 如何影响模型的输出，我们可以将该 feature 的 SHAP 值与数据集中所有样本的 feature 值进行比较：

多样本可视化探索

下面的图形是针对多患者的预测和影响因素的探索。在 jupyter notebook 的交互式作用下，能够观察不同的特征属性对前 50 个患者的影响：

shap_values = explainer.shap_values(X_train.iloc[:50])shap.force_plot(explainer.expected_value[1],                 shap_values[1],                 X_test.iloc[:50])