SHAP（四）：NHANES I 生存模型

这是一个 Cox 比例风险模型，基于来自 NHANES I 的数据以及来自 NHANES I 流行病学随访研究。它旨在说明 SHAP 值如何能够以传统上仅由线性模型提供的清晰度解释 XGBoost 模型。我们在数据中看到有趣的非线性模式，这表明了这种方法的潜力。请记住，我们尚未对数据进行检查以校准当前的实验室测试，因此您不应将结果视为可操作的医学见解，而应将其视为概念证明。

请注意，对 Cox 损失和 SHAP 交互效果的支持最近才合并，因此您需要最新的 XGBoost 主版本才能运行此笔记本。

import matplotlib.pylab as pl
import xgboost
from sklearn.model_selection import train_test_splitimport shap

1.创建 XGBoost 数据对象

这使用了 SHAP 数据集模块中可用的 NHANES I 数据的预处理子集。

X, y = shap.datasets.nhanesi()
X_display, y_display = shap.datasets.nhanesi(display=True
)  # human readable feature valuesxgb_full = xgboost.DMatrix(X, label=y)# create a train/test split
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=7)
xgb_train = xgboost.DMatrix(X_train, label=y_train)
xgb_test = xgboost.DMatrix(X_test, label=y_test)

2.训练 XGBoost 模型

# use validation set to choose # of trees
params = {"eta": 0.002, "max_depth": 3, "objective": "survival:cox", "subsample": 0.5}
model_train = xgboost.train(params, xgb_train, 10000, evals=[(xgb_test, "test")], verbose_eval=1000
)

[0]	test-cox-nloglik:7.26952
[1000]	test-cox-nloglik:6.55767
[2000]	test-cox-nloglik:6.48836
[3000]	test-cox-nloglik:6.47129
[4000]	test-cox-nloglik:6.46786
[5000]	test-cox-nloglik:6.46583
[6000]	test-cox-nloglik:6.46623
[7000]	test-cox-nloglik:6.46841
[8000]	test-cox-nloglik:6.46972
[9000]	test-cox-nloglik:6.47175
[9999]	test-cox-nloglik:6.47396

# train final model on the full data set
params = {"eta": 0.002, "max_depth": 3, "objective": "survival:cox", "subsample": 0.5}
model = xgboost.train(params, xgb_full, 5000, evals=[(xgb_full, "test")], verbose_eval=1000
)

[0]	test-cox-nloglik:8.88073
[1000]	test-cox-nloglik:8.17142
[2000]	test-cox-nloglik:8.08556
[3000]	test-cox-nloglik:8.04853
[4000]	test-cox-nloglik:8.0248
[4999]	test-cox-nloglik:8.00511

3.检查性能

C 统计量衡量我们如何根据人们的生存时间对他们进行排序（1.0 是完美排序）。

def c_statistic_harrell(pred, labels):total = 0matches = 0for i in range(len(labels)):for j in range(len(labels)):if labels[j] > 0 and abs(labels[i]) > labels[j]:total += 1if pred[j] > pred[i]:matches += 1return matches / total# see how well we can order people by survival
c_statistic_harrell(model_train.predict(xgb_test, ntree_limit=5000), y_test)

0.835090082176807

4.解释模型对整个数据集的预测

shap_values = shap.TreeExplainer(model).shap_values(X)

4.1 SHAP 摘要图

XGBoost 的 SHAP 值解释了模型的边际输出，即 Cox 比例风险模型的死亡对数几率的变化。我们可以从下面看到，根据模型，死亡的主要危险因素是年老。死亡风险的下一个最有力的指标是男性。

该摘要图取代了特征重要性的典型条形图。它告诉我们哪些特征是最重要的，以及它们对数据集的影响范围。颜色使我们能够匹配特征值的变化如何影响风险的变化（例如高白细胞计数导致高死亡风险）。

shap.summary_plot(shap_values, X)

在这里插入图片描述

4.2 SHAP 相关图

SHAP 摘要图给出了每个特征的总体概述，而 SHAP 依赖图显示了模型输出如何随特征值变化。请注意，每个点都是一个人，单个特征值的垂直分散是由模型中的交互效应产生的。自动选择用于着色的功能来突出显示可能驱动这些交互的因素。稍后我们将了解如何使用 SHAP 交互值检查模型中是否确实存在交互。请注意，SHAP 汇总图的行是将 SHAP 相关图的点投影到 y 轴上，然后由特征本身重新着色得到的。

下面我们给出了每个 NHANES I 特征的 SHAP 依赖图，揭示了有趣但预期的趋势。请记住，其中一些值的校准可能与现代实验室测试不同，因此得出结论时要小心。

# we pass "Age" instead of an index because dependence_plot() will find it in X's column names for us
# Systolic BP was automatically chosen for coloring based on a potential interaction to check that
# the interaction is really in the model see SHAP interaction values below
shap.dependence_plot("Age", shap_values, X)

在这里插入图片描述

# we pass display_features so we get text display values for sex
shap.dependence_plot("Sex", shap_values, X, display_features=X_display)

在这里插入图片描述

# setting show=False allows us to continue customizing the matplotlib plot before displaying it
shap.dependence_plot("Systolic BP", shap_values, X, show=False)
pl.xlim(80, 225)
pl.show()

在这里插入图片描述

shap.dependence_plot("Poverty index", shap_values, X)

在这里插入图片描述

shap.dependence_plot("White blood cells", shap_values, X, display_features=X_display, show=False
)
pl.xlim(2, 15)
pl.show()

在这里插入图片描述

shap.dependence_plot("BMI", shap_values, X, display_features=X_display, show=False)
pl.xlim(15, 50)
pl.show()

在这里插入图片描述

shap.dependence_plot("Serum magnesium", shap_values, X, show=False)
pl.xlim(1.2, 2.2)
pl.show()

在这里插入图片描述

shap.dependence_plot("Sedimentation rate", shap_values, X)

在这里插入图片描述

shap.dependence_plot("Serum protein", shap_values, X)

在这里插入图片描述

shap.dependence_plot("Serum cholesterol", shap_values, X, show=False)
pl.xlim(100, 400)
pl.show()

请添加图片描述

shap.dependence_plot("Pulse pressure", shap_values, X)

在这里插入图片描述

shap.dependence_plot("Serum iron", shap_values, X, display_features=X_display)

在这里插入图片描述

shap.dependence_plot("TS", shap_values, X)

在这里插入图片描述

shap.dependence_plot("Red blood cells", shap_values, X)

在这里插入图片描述

5.计算 SHAP 交互值

有关更多详细信息，请参阅 Tree SHAP 论文，但简单地说，SHAP 交互值是 SHAP 值对更高阶交互的推广。最新版本的 XGBoost 中使用 pred_interactions 标志实现了成对交互的快速精确计算。使用此标志，XGBoost 为每个预测返回一个矩阵，其中主效应位于对角线上，交互效应位于非对角线上。主效应类似于线性模型获得的 SHAP 值，交互效应捕获所有高阶交互，并将它们划分为成对交互项。请注意，整个交互矩阵的总和是模型当前输出与预期输出之间的差，因此非对角线上的交互效应被分成两半（因为每个都有两个）。绘制交互效果时，SHAP 包会自动将非对角线值乘以 2，以获得完整的交互效果。

# takes a couple minutes since SHAP interaction values take a factor of 2 * # features
# more time than SHAP values to compute, since this is just an example we only explain
# the first 2,000 people in order to run quicker
shap_interaction_values = shap.TreeExplainer(model).shap_interaction_values(X.iloc[:2000, :]
)

5.1 SHAP 交互值汇总图

SHAP 交互值矩阵的汇总图绘制了汇总图矩阵，其中对角线上有主效应，对角线外有交互效应。

shap.summary_plot(shap_interaction_values, X.iloc[:2000, :])

在这里插入图片描述

5.2 SHAP 交互值依赖图

对 SHAP 交互值 a 运行依赖图可以让我们分别观察主效应和交互效应。

下面我们绘制了年龄的主要影响以及年龄的一些交互影响。将年龄的主效应图与早期的年龄 SHAP 值图进行比较可以提供丰富的信息。主效应图没有垂直分散，因为相互作用效应全部以非对角线项捕获。

shap.dependence_plot(("Age", "Age"),shap_interaction_values,X.iloc[:2000, :],display_features=X_display.iloc[:2000, :],
)

在这里插入图片描述

现在我们绘制涉及年龄的交互效应。这些效应捕获了原始 SHAP 图中存在但上面的主效应图中缺失的所有垂直色散。下图涉及年龄和性别，显示基于性别的死亡风险差距因年龄而异，并在 60 岁时达到峰值。

shap.dependence_plot(("Age", "Sex"),shap_interaction_values,X.iloc[:2000, :],display_features=X_display.iloc[:2000, :],
)

在这里插入图片描述

shap.dependence_plot(("Age", "Systolic BP"),shap_interaction_values,X.iloc[:2000, :],display_features=X_display.iloc[:2000, :],
)

在这里插入图片描述

shap.dependence_plot(("Age", "White blood cells"),shap_interaction_values,X.iloc[:2000, :],display_features=X_display.iloc[:2000, :],
)

在这里插入图片描述

shap.dependence_plot(("Age", "Poverty index"),shap_interaction_values,X.iloc[:2000, :],display_features=X_display.iloc[:2000, :],
)

在这里插入图片描述

shap.dependence_plot(("Age", "BMI"),shap_interaction_values,X.iloc[:2000, :],display_features=X_display.iloc[:2000, :],
)

在这里插入图片描述

shap.dependence_plot(("Age", "Serum magnesium"),shap_interaction_values,X.iloc[:2000, :],display_features=X_display.iloc[:2000, :],
)

在这里插入图片描述

Now we show a couple examples with systolic blood pressure.

shap.dependence_plot(("Systolic BP", "Systolic BP"),shap_interaction_values,X.iloc[:2000, :],display_features=X_display.iloc[:2000, :],
)

在这里插入图片描述

shap.dependence_plot(("Systolic BP", "Age"),shap_interaction_values,X.iloc[:2000, :],display_features=X_display.iloc[:2000, :],
)

在这里插入图片描述

shap.dependence_plot(("Systolic BP", "Age"),shap_interaction_values,X.iloc[:2000, :],display_features=X_display.iloc[:2000, :],
)

在这里插入图片描述

import matplotlib.pylab as pl
import numpy as np

tmp = np.abs(shap_interaction_values).sum(0)
for i in range(tmp.shape[0]):tmp[i, i] = 0
inds = np.argsort(-tmp.sum(0))[:50]
tmp2 = tmp[inds, :][:, inds]
pl.figure(figsize=(12, 12))
pl.imshow(tmp2)
pl.yticks(range(tmp2.shape[0]), X.columns[inds], rotation=50.4, horizontalalignment="right"
)
pl.xticks(range(tmp2.shape[0]), X.columns[inds], rotation=50.4, horizontalalignment="left"
)
pl.gca().xaxis.tick_top()
pl.show()