NHANES I 生存模型

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

import matplotlib.pylab as pl
import numpy as np
import xgboost
from sklearn.model_selection import train_test_split

import shap

创建 XGBoost 数据对象

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

X, y = shap.datasets.nhanesi()

# human readable feature values
X_display, y_display = shap.datasets.nhanesi(display=True)

xgb_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)

训练 XGBoost 模型

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.67918
[1000]  test-cox-nloglik:7.02985
[2000]  test-cox-nloglik:6.97516
[3000]  test-cox-nloglik:6.96240
[4000]  test-cox-nloglik:6.96217
[5000]  test-cox-nloglik:6.96558
[6000]  test-cox-nloglik:6.96995
[7000]  test-cox-nloglik:6.97310
[8000]  test-cox-nloglik:6.97721
[9000]  test-cox-nloglik:6.98116
[9999]  test-cox-nloglik:6.98511

# 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:9.28404
[1000]  test-cox-nloglik:8.60868
[2000]  test-cox-nloglik:8.53134
[3000]  test-cox-nloglik:8.49490
[4000]  test-cox-nloglik:8.47122
[4999]  test-cox-nloglik:8.45328

检查性能

C-统计量衡量我们按生存时间对人们进行排序的能力（1.0 表示完美的排序）。

def c_statistic_harrell(pred, labels):
    total = 0
    matches = 0
    for i in range(len(labels)):
        for j in range(len(labels)):
            if labels[j] > 0 and abs(labels[i]) > labels[j]:
                total += 1
                if pred[j] > pred[i]:
                    matches += 1
    return matches / total


# see how well we can order people by survival
c_statistic_harrell(model_train.predict(xgb_test), y_test)

0.817035332310394

解释模型在整个数据集上的预测

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

SHAP 摘要图

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

这个总结图取代了典型的特征重要性条形图。它告诉我们哪些特征最重要，以及它们在整个数据集上的影响范围。颜色帮助我们可视化特征值的变化如何影响风险的变化（例如，高白细胞计数导致高死亡风险）。

shap.summary_plot(shap_values, X)

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_isFemale", 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_blood_pressure", shap_values, X, show=False)
pl.xlim(80, 225)
pl.show()

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("sedimentation_rate", shap_values, X)

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

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

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

计算 SHAP 交互值

更多细节请参阅 Tree SHAP 论文，但简而言之，SHAP 交互值是 SHAP 值对高阶交互的推广。在 XGBoost 的后期版本（>=1.0.0）中，通过 pred_interactions 标志实现了成对交互的快速精确计算。使用此标志，XGBoost 会为每个预测返回一个矩阵，其中主效应在对角线上，交互效应在非对角线上。

主要效果类似于你在线性模型中得到的 SHAP 值，而交互效果则捕捉所有高阶交互并将它们分配给成对的交互项。请注意，整个交互矩阵的总和是模型当前输出与预期输出之间的差异，因此非对角线上的交互效果被分成两半（因为每对有两个）。

在绘制交互效应时，SHAP 包会自动将非对角线值乘以二以获得完整的交互效应。

# we only take 300 people in order to run quicker
number_patients = 300

shap_interaction_values = shap.TreeExplainer(model).shap_interaction_values(
    X.iloc[:number_patients, :]
)

SHAP 交互值总结图

SHAP 交互值矩阵的摘要图绘制了一个摘要图矩阵，对角线上是主效应，非对角线上是交互效应。

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

SHAP 交互值依赖图

在SHAP交互值上运行依赖图，使我们能够分别观察主要效应和交互效应。

下面我们绘制了年龄的主要效应以及年龄的一些交互效应。将年龄的主要效应图与之前的年龄SHAP值图进行比较是有信息量的。主要效应图没有垂直分散，因为所有的交互效应都被捕捉到了非对角项中。

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

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

shap.dependence_plot(
    ("age", "sex_isFemale"),
    shap_interaction_values,
    X.iloc[:number_patients, :],
    display_features=X_display.iloc[:number_patients, :],
)

shap.dependence_plot(
    ("age", "systolic_blood_pressure"),
    shap_interaction_values,
    X.iloc[:number_patients, :],
    display_features=X_display.iloc[:number_patients, :],
)

shap.dependence_plot(
    ("age", "white_blood_cells"),
    shap_interaction_values,
    X.iloc[:number_patients, :],
    display_features=X_display.iloc[:number_patients, :],
)

shap.dependence_plot(
    ("age", "bmi"),
    shap_interaction_values,
    X.iloc[:number_patients, :],
    display_features=X_display.iloc[:number_patients, :],
)

现在我们展示几个关于收缩压的例子。

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

shap.dependence_plot(
    ("systolic_blood_pressure", "age"),
    shap_interaction_values,
    X.iloc[:number_patients, :],
    display_features=X_display.iloc[:number_patients, :],
)

shap.dependence_plot(
    ("systolic_blood_pressure", "age"),
    shap_interaction_values,
    X.iloc[:number_patients, :],
    display_features=X_display.iloc[:number_patients, :],
)

我们将所有患者之间的交互汇总，以绘制哪些特征交互最多。请注意，我们不分析特征如何交互（高特征A + 低特征B 导致结果C，等等）。浅色显示强烈的交互效应。

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