专栏名称: 走天涯徐小洋地理数据科学

一个爱生活的地理土博，分享GIS、遥感、空间分析、R语言、景观生态等地理数据科学实操教程、经典文献、数据资源

HF.099 | 黑箱模型的事后可解释性-SHAP框架

走天涯徐小洋地理数据科学 · 公众号 · · 2025-03-03 15:39

正文

请到「今天看啥」查看全文

SHAP （SHapley Additive exPlanations）是一种基于博弈论的 模型解释方法 ，能够量化每个特征对机器学习模型预测的贡献，帮助理解 黑盒模型 的 内部决策逻辑 。其作为可解释机器学习领域的核心工具，近年来在全球水文领域的学术研究和实际应用中产生了深远影响。

论文速递1：最新研究《Compounding effects in flood drivers challenge estimates of extreme river floods》创新性引入 SHAP模型 ，首次量化降水、地形、土地利用等多因子的非线性协同效应，揭示传统方法低估的复合致灾链通过全局特征贡献排序+局部样本归因，SHAP不仅定位关键驱动因子（如短时强降雨贡献度达52%），更动态拆解“极端事件预测逻辑”，为跨区域洪水风险评估提供“可解释-可干预”的决策地图！

图1.说明了在一个流域内的两个洪水事件样本的复合效应的重要性。（A和B），两个洪水样本的输入特征（彩色条）和模型输出（橙色点）。(C)特征之间的成对交互效应（包括对角线中特征的主效应）；(D)：(B)中洪水样特征间的两两交互作用效应，复合效应的重要性为22.9%。

引文链接：Jiang, S., Tarasova, L., Yu, G., & Zscheischler, J. (2024). Compounding effects in flood drivers challenge estimates of extreme river floods. Science Advances , 10(13), eadl4005.

论文速递 2：最新研究《How Interpretable ML Can Benefit Process Understanding in the Geosciences》通过SHAP模型首次实现地球科学多因子耦合作用的动态归因，从冰川消融到岩石风化，拆解非线性机制与时空异质性，为复杂地质过程建模提供“可量化-可验证”的解释框架。文中创新性融合SHAP依赖图+地质过程链，不仅揭示关键驱动因子（如温度突变的SHAP贡献度达64%），更通过局部样本归因追溯极端事件的因果路径，推动地学模型从“经验驱动”向“数据-机理双驱动”跨越！

图2：将可解释的机器学习应用于地球科学过程理解的工作流程和示例

引文链接： Jiang, S., Sweet, L. B., Blougouras, G., Brenning, A., Li, W., Reichstein, M., ... & Zscheischler, J. (2024). How interpretable machine learning can benefit process understanding in the geosciences. Earth's Future , 12(7), e2024EF004540.

SHAP概念与理论基础

SHAP 将博弈论的Shapley值与机器学习结合，提出了加性特征归因方法（Additive Feature Attribution），统一了LIME、DeepLIFT等6种传统解释方法的理论基础，解决了不同方法结果不一致的问题。例如，通过公式化证明，SHAP是唯一满足局部准确性（local accuracy）、缺失特征归零（missingness）和一致性（consistency）三大公理的解释模型。

简单举个例子：多个工人合作完成工件时，Shapley值通过计算所有可能的合作组合中每个参与者的 边际贡献 ，确定其应得报酬。SHAP将此理论迁移到机器学习中，将特征视为“参与者”，模型预测结果视为“总收益”，特征贡献通过Shapley值量化。

SHAP通过构建加性解释模型，将预测值分解为每个特征的线性组合：

其中， ϕ ₀ 是基线预测（所有特征取平均时的预测值）， ϕ _i 表示第ii个特征的SHAP值（贡献值），正值表示正向影响，负值反之。

SHAP核心算法

1.Shapley值计算原理：每个特征的SHAP值通过计算其在所有特征组合中的边际贡献加权平均得到。

其中，SS为不含特征ii的子集，MM为总特征数，f(S)f(S)表示仅使用子集SS的模型预测值。

2.近似估计算法：直接计算Shapley值的复杂度为 O(2 ^M ) 针对高维数据需采用近似方法：

（1）随机采样法：随机生成特征排列，构造包含/不包含目标特征的虚拟样本，计算预测差异并取平均；

（2）蒙特卡洛积分：通过大量迭代逼近真实值，适用于非树模型。

3.TreeSHAP——树模型的高效计算：针对决策树类模型（如XGBoost、LightGBM），TreeSHAP将复杂度从 O(TL2 ^M ) 优化至 O(TLD ² ) （TT为树数量，LL为叶子节点数，DD为树深度），通过动态规划递归计算特征路径概率，避免穷举所有子集。

SHAP代码实现

一、Python实现（基于XGBoost与SHAP库）

1. 环境安装与数据准备

# 安装核心库!pip install xgboost shap pandas scikit-learn 
# 导入模块 import xgboost as xgb import shap from sklearn.datasets  import make_classification from sklearn.model_selection  import train_test_split 
# 生成模拟数据（二分类任务）X, y = make_classification(n_samples=1000, n_features=10, random_state=42)X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, random_state=42)

2. 模型训练与评估

# 构建XGBoost模型 params = {'objective':'binary:logistic', 'eval_metric':'logloss', 'max_depth':3}dtrain = xgb.DMatrix(X_train, label=y_train)model = xgb.train(params,  dtrain, num_boost_round=100)
# 预测与性能评估 dtest = xgb.DMatrix(X_test)y_pred = model.predict(dtest)

3. SHAP值计算与可视化

# 创建SHAP解释器explainer = shap.TreeExplainer(model)shap_values = explainer.shap_values(X_test) 
# 全局解释：特征重要性条形图shap.summary_plot(shap_values,  X_test, feature_names=[f'Feature_{i}' for i in range(10)])
# 局部解释：单个样本决策力图 shap.force_plot(explainer.expected_value,  shap_values[0](), X_test[0](), feature_names=feature_names)
# 交互作用分析：依赖图shap.dependence_plot("Feature_0",  shap_values, X_test, interaction_index="Feature_1")

二、R语言实现（通过shapviz包）

1. 环境配置与数据加载

# 安装包install.packages(c("xgboost",  "shapviz", "ggplot2"))





    
# 加载库library(xgboost)library(shapviz)
# 加载数据 data(iris)X 5]  # 特征矩阵 y as.numeric(iris$Species  == "versicolor")  # 二分类标签

2. 模型训练与SHAP分析

# 训练XGBoost模型 model   data = as.matrix(X),   label = y,  nrounds = 100,  objective = "binary:logistic")
# 计算SHAP值shap as.matrix(X)) 
# 可视化：蜂群图 sv_importance(shap, kind = "bee") +   ggtitle("Global Feature Importance")

SHAP相关示例

示例1：LightGBM-SHAP（R语言实现）：

本示例介绍了如何在R中利用SHAP值解释LightGBM模型的预测结果。

LightGBM是一个基于Gradient Boosting的高效、分布式、可扩展和并行的开源库，是一种高效的梯度提升决策树算法，但其黑盒性质使得理解模型变得困难。LightGBM的核心特点是它使用了一种基于分区的决策树学习算法，这种算法可以有效地处理大规模数据集和高维特征。为了提高模型的可解释性，这里我们使用SHAP进行模型解释和提高可解释性，并提供相应代码示例：

library(shapviz)library(ggplot2)install.packages("lightgbm") library(lightgbm)
#数据加载及处理df "MPE.csv")# df vfactor "MPE","Gender","Cough")df[vfactor] df[vfactor], factor)
set.seed(2)dftrain df[,-1]),label = df[,1])lgb "binary",                               learning_rate=0.3),                 data=dftrain,                 nrounds=1000,                 verbose = -1)
shp df[,-1]))sv_waterfall(shp,row_id = 2,max_display = 10)

sv_force(shp,row_id = 2)

sv_importance(shp,fill="#0085FF")

sv_importance(shp, kind = "beeswarm")

sv_importance(shp,               kind = "both",              show_numbers = TRUE,               bee_width = 




    
0.2)

sv_dependence(shp,              v="ADAPE",              color_var = NULL)

示例2： NHANES I 生存 模型（Python实现）

该研究基于美国国家健康与营养调查一期(NHANES I)的随访死亡率数据构建Cox比例风险模型，通过整合NHANES I流行病学随访研究的纵向追踪信息，旨在验证SHAP解释框架在XGBoost生存模型中的可解释性优势。研究发现，SHAP值能够将树集成模型的复杂决策过程解构为具有临床可读性的特征贡献度图谱，其解释粒度可达到与传统线性模型的可解释性层级。

（1）数据导入：

import shapimport xgboostfrom sklearn.model_selection import train_test_splitimport matplotlib.pylab as plX,y = shap.datasets.nhanesi()X_display,y_display = shap.datasets.nhanesi(display=True) # human readable feature values
xgb_full = xgboost.DMatrix(X, label=y)

（2）创建 XGBoost 数据对象：

# create a train/test splitX_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)

（3）训练 XGBoost 模型

# use validation set to choose # of treesparams = {    "eta": 0.002,    "max_depth": 3,    "objective": "survival:cox",    "subsample": 0.5}model_train = xgboost.train(params, xgb_train, 5000, evals = [(xgb_test, "test")], verbose_eval=1000)# train final model on the full data setparams = {    "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)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

（4）检查性能及解释模型对整个数据集的预测

# see how well we can order people by survivalc_statistic_harrell(model_train.predict(xgb_test,  iteration_range=(0, 5000)), y_test)

（5）SHAP图绘制代码：

shap.summary_plot(shap_values, X)




    
shap.dependence_plot("age", shap_values, X)shap.dependence_plot("sex_isFemale", shap_values, X, display_features=X_display)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("cholesterol", shap_values, X, show=False)pl.xlim(100,400)pl.show()shap.dependence_plot("pulse_pressure", shap_values, X)shap.dependence_plot("red_blood_cells", shap_values, X)

shap_interaction_values = shap.TreeExplainer(model).shap_interaction_values(X.iloc[:2000,:])shap.summary_plot(shap_interaction_values, X.iloc[:2000,:])shap.dependence_plot(    ("age", "age"),    shap_interaction_values, X.iloc[:2000,:],    display_features=X_display.iloc[:2000,:])shap.dependence_plot(    ("age", "sex_isFemale"),    shap_interaction_values, X.iloc[:2000,:],    display_features=X_display.iloc[:2000,:])shap.dependence_plot(    ("age", "systolic_blood_pressure"),    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", "bmi"),    shap_interaction_values, X.iloc[:2000,:],    display_features=X_display.iloc[:2000,:])shap.dependence_plot(    ("age", "serum_protein"),    shap_interaction_values, X.iloc[:2000,:],    display_features=X_display.iloc[:2000,:])

（6） SHAP图汇总：

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

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

（7）SHAP依赖性图：

SHAP 摘要图提供了每个特征的一般概述，而 SHAP 依赖性图显示了模型输出如何随特征值而变化。请注意，每个点都是一个人，单个特征值的垂直离散是模型中交互效应的结果。SHAP 汇总图的行是将 SHAP 依赖性图的点投影到 y 轴上，然后按特征本身重新着色的结果。

下面我们给出了部分 NHANES I 特征的 SHAP 依赖性图，揭示了有趣但意料之中的趋势。

（8）SHAP 交互作用值汇总图

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

（9）SHAP 交互作用值依赖性图

对 SHAP 交互作用值 a 运行依赖图，可以让我们分别观察主效应和交互作用效应。绘制涉及年龄的交互效应。这些效应捕获了原始 SHAP 图中存在但主效应图中缺少的所有垂直离散。

总结： LightGBM与XGBoost的主要区别在于它们的决策树学习算法。LightGBM使用了基于分区的决策树学习算法，而XGBoost使用了基于页面的决策树学习算法。此外，LightGBM使用了Histogram Based Method来估计损失函数，这种方法可以减少内存占用和计算量。

一图胜千言！水文图绘改版后致力于分享水文相关的精美图表，为读者提供作图思路和经验，帮助大家制作更漂亮丰富的图表。同时欢迎留言咨询绘图难点，我们会针对性地分享相关绘制经验。另外也期待读者踊跃来稿，分享更好的构图思维和技巧，稿件可发送至邮箱[email protected], 或者联系微信17339888901投稿。