SHAP（一）：具有 Shapley 值的可解释 AI 简介-慈云数据

SHAP（一）：具有 Shapley 值的可解释 AI 简介

这是用 Shapley 值解释机器学习模型的介绍。沙普利值是合作博弈论中广泛使用的方法，具有理想的特性。本教程旨在帮助您深入了解如何计算和解释基于 Shapley 的机器学习模型解释。我们将采取实用的实践方法，使用“shap”Python 包来逐步解释更复杂的模型。这是一个动态文档，作为“shap” Python 包的介绍。因此，如果您有反馈或贡献，请提出问题或拉取请求，以使本教程变得更好！

大纲

解释线性回归模型
解释广义加性回归模型
解释非加性提升树模型
解释线性逻辑回归模型
解释非加性提升树逻辑回归模型
处理相关输入特征
1.解释线性回归模型

在使用 Shapley 值解释复杂模型之前，了解它们如何适用于简单模型会很有帮助。最简单的模型类型之一是标准线性回归，因此下面我们在[加州住房数据集](https://www.dcc.fc.up.pt/~ltorgo/Regression/cal_housing.html ）。该数据集由 1990 年加利福尼亚州的 20,640 个房屋区块组成，我们的目标是根据 8 个不同的特征预测房价中位数的自然对数：
1. MedInc - 区块组收入中位数
2. HouseAge - 街区组中的房屋年龄中位数
3. AveRooms——每户平均房间数
4. AveBedrms——每户平均卧室数量
  5.人口-区块组人口
5. AveOccup - 家庭成员的平均数量
6. Latitude——块组纬度
7. Longitude——块组经度
```
import sklearn
import shap
# a classic housing price dataset
X, y = shap.datasets.california(n_points=1000)
X100 = shap.utils.sample(X, 100)  # 100 instances for use as the background distribution
# a simple linear model
model = sklearn.linear_model.LinearRegression()
model.fit(X, y)
```
```
LinearRegression()
```
1.1 检查模型系数

理解线性模型的最常见方法是检查为每个特征学习的系数。这些系数告诉我们当我们改变每个输入特征时模型输出会发生多少变化：
```
print("Model coefficients:\n")
for i in range(X.shape[1]):
    print(X.columns[i], "=", model.coef_[i].round(5))
```
```
Model coefficients:
MedInc = 0.45769
HouseAge = 0.01153
AveRooms = -0.12529
AveBedrms = 1.04053
Population = 5e-05
AveOccup = -0.29795
Latitude = -0.41204
Longitude = -0.40125
```
虽然系数非常适合告诉我们当我们改变输入特征的值时会发生什么，但它们本身并不是衡量特征整体重要性的好方法。这是因为每个系数的值取决于输入特征的规模。例如，如果我们以分钟而不是年为单位来测量房屋的年龄，那么 HouseAge 特征的系数将变为 0.0115 / (3652460) = 2.18e-8。显然，房屋建成后的年数并不比分钟数更重要，但其系数值要大得多。这意味着系数的大小不一定能很好地衡量线性模型中特征的重要性。

1.2 使用部分依赖图更完整的图片

要了解模型中特征的重要性，有必要了解更改该特征如何影响模型的输出，以及该特征值的分布。为了将其可视化为线性模型，我们可以构建经典的部分依赖图，并将特征值的分布显示为 x 轴上的直方图：
```
shap.partial_dependence_plot(
    "MedInc",
    model.predict,
    X100,
    ice=False,
    model_expected_value=True,
    feature_expected_value=True,
)
```
上图中的灰色水平线表示模型应用于加州住房数据集时的预期值。垂直灰线表示中位收入特征的平均值。请注意，蓝色部分依赖图线（即我们将中位收入特征固定为给定值时模型输出的平均值）始终穿过两条灰色期望值线的交点。我们可以将该交点视为数据分布的部分依赖图的“中心”。当我们接下来讨论SHAP值时，这种中心化的影响就会变得清晰。

1.3 从部分相关图中读取 SHAP 值

基于 Shapley 值的机器学习模型解释背后的核心思想是使用合作博弈论的公平分配结果，在模型的输入特征中为模型的输出 f ( x ) f(x) f(x) 分配信用。为了将博弈论与机器学习模型联系起来，既需要将模型的输入特征与游戏中的玩家进行匹配，又需要将模型函数与游戏规则进行匹配。由于在博弈论中玩家可以加入或不加入游戏，因此我们需要一种方法来让功能“加入”或“不加入”模型。定义某个特征“加入”模型的含义的最常见方法是，当我们知道该特征的值时，就说该特征已“加入模型”；而当我们不知道该特征的值时，就说该特征尚未加入模型。知道该功能的价值。当只有特征子集 S S S 是模型的一部分时，为了评估现有模型 f f f，我们使用条件期望值公式整合其他特征。这个公式可以有两种形式：

E [ f ( X ) ∣ X S = x S ] E[f(X) \mid X_S = x_S] E[f(X)∣XS=xS]
or
E [ f ( X ) ∣ d o ( X S = x S ) ] E[f(X) \mid do(X_S = x_S)] E[f(X)∣do(XS=xS)]

在第一种形式中，我们知道 S 中特征的值，因为我们“观察”它们。在第二种形式中，我们知道 S 中特征的值，因为我们“设置”了它们。一般来说，第二种形式通常更可取，因为它告诉我们如果我们干预并更改其输入，模型将如何表现，而且因为它更容易计算。在本教程中，我们将完全关注第二种表述。我们还将使用更具体的术语“SHAP 值”来指代应用于机器学习模型的条件期望函数的 Shapley 值。

SHAP 值的计算可能非常复杂（它们通常是 NP 困难的），但线性模型非常简单，我们可以立即从部分依赖图读取 SHAP 值。当我们解释预测 f ( x ) f(x) f(x) 时，特定特征 i i i 的 SHAP 值只是预期模型输出与特征值 x i x_i xi 处的部分依赖图之间的差异：
```
# compute the SHAP values for the linear model
explainer = shap.Explainer(model.predict, X100)
shap_values = explainer(X)
# make a standard partial dependence plot
sample_ind = 20
shap.partial_dependence_plot(
    "MedInc",
    model.predict,
    X100,
    model_expected_value=True,
    feature_expected_value=True,
    ice=False,
    shap_values=shap_values[sample_ind : sample_ind + 1, :],
)
```
经典的部分依赖图和 SHAP 值之间的紧密对应意味着，如果我们在整个数据集中绘制特定特征的 SHAP 值，我们将准确地绘制出该特征的部分依赖图的平均中心版本：
```
shap.plots.scatter(shap_values[:, "MedInc"])
```
1.4 Shapley 值的加性

Shapley 值的基本属性之一是，它们总是总结所有玩家在场时的游戏结果与没有玩家在场时的游戏结果之间的差异。对于机器学习模型，这意味着所有输入特征的 SHAP 值将始终等于基线（预期）模型输出与所解释的预测的当前模型输出之间的差异。最简单的方法是通过瀑布图，该图从我们对房价 E [ f ( X ) ] E[f(X)] E[f(X)] 的背景先验期望开始，然后一次添加一个特征，直到达到当前模型输出 f ( x ) f( x) f(x):
```
# the waterfall_plot shows how we get from shap_values.base_values to model.predict(X)[sample_ind]
shap.plots.waterfall(shap_values[sample_ind], max_display=14)
```
2.解释加性回归模型

线性模型的部分依赖图与 SHAP 值具有如此密切的联系的原因是，模型中的每个特征都是独立于其他每个特征进行处理的（效果只是相加在一起）。我们可以保持这种可加性，同时放宽直线的线性要求。这就产生了众所周知的广义加性模型 (GAM)。虽然有很多方法可以训练这些类型的模型（例如将 XGBoost 模型设置为深度 1），但我们将使用专门为此设计的 InterpretML 可解释的 boosting 机器。
```
# fit a GAM model to the data
import interpret.glassbox
model_ebm = interpret.glassbox.ExplainableBoostingRegressor(interactions=0)
model_ebm.fit(X, y)
# explain the GAM model with SHAP
explainer_ebm = shap.Explainer(model_ebm.predict, X100)
shap_values_ebm = explainer_ebm(X)
# make a standard partial dependence plot with a single SHAP value overlaid
fig, ax = shap.partial_dependence_plot(
    "MedInc",
    model_ebm.predict,
    X100,
    model_expected_value=True,
    feature_expected_value=True,
    show=False,
    ice=False,
    shap_values=shap_values_ebm[sample_ind : sample_ind + 1, :],
)
```
```
shap.plots.scatter(shap_values_ebm[:, "MedInc"])
```
```
# the waterfall_plot shows how we get from explainer.expected_value to model.predict(X)[sample_ind]
shap.plots.waterfall(shap_values_ebm[sample_ind])
```
```
# the waterfall_plot shows how we get from explainer.expected_value to model.predict(X)[sample_ind]
shap.plots.beeswarm(shap_values_ebm)
```
3.解释非加性提升树模型
```
# train XGBoost model
import xgboost
model_xgb = xgboost.XGBRegressor(n_estimators=100, max_depth=2).fit(X, y)
# explain the GAM model with SHAP
explainer_xgb = shap.Explainer(model_xgb, X100)
shap_values_xgb = explainer_xgb(X)
# make a standard partial dependence plot with a single SHAP value overlaid
fig, ax = shap.partial_dependence_plot(
    "MedInc",
    model_xgb.predict,
    X100,
    model_expected_value=True,
    feature_expected_value=True,
    show=False,
    ice=False,
    shap_values=shap_values_xgb[sample_ind : sample_ind + 1, :],
)
```
```
shap.plots.scatter(shap_values_xgb[:, "MedInc"])
```
```
shap.plots.scatter(shap_values_xgb[:, "MedInc"], color=shap_values)
```
4.解释线性逻辑回归模型
```
# a classic adult census dataset price dataset
X_adult, y_adult = shap.datasets.adult()
# a simple linear logistic model
model_adult = sklearn.linear_model.LogisticRegression(max_iter=10000)
model_adult.fit(X_adult, y_adult)
def model_adult_proba(x):
    return model_adult.predict_proba(x)[:, 1]
def model_adult_log_odds(x):
    p = model_adult.predict_log_proba(x)
    return p[:, 1] - p[:, 0]
```
请注意，解释线性逻辑回归模型的概率在输入中不是线性的。
```
# make a standard partial dependence plot
sample_ind = 18
fig, ax = shap.partial_dependence_plot(
    "Capital Gain",
    model_adult_proba,
    X_adult,
    model_expected_value=True,
    feature_expected_value=True,
    show=False,
    ice=False,
)
```
如果我们使用 SHAP 来解释线性逻辑回归模型的概率，我们会看到很强的交互效应。这是因为线性逻辑回归模型在概率空间中不是可加的。
```
# compute the SHAP values for the linear model
background_adult = shap.maskers.Independent(X_adult, max_samples=100)
explainer = shap.Explainer(model_adult_proba, background_adult)
shap_values_adult = explainer(X_adult[:1000])
```
```
Permutation explainer: 1001it [00:58, 14.39it/s]                                  
```
```
shap.plots.scatter(shap_values_adult[:, "Age"])
```
如果我们解释模型的对数赔率输出，我们会看到模型输入和模型输出之间存在完美的线性关系。重要的是要记住您正在解释的模型的单位是什么，并且解释不同的模型输出可能会导致对模型行为的截然不同的看法。
```
# compute the SHAP values for the linear model
explainer_log_odds = shap.Explainer(model_adult_log_odds, background_adult)
shap_values_adult_log_odds = explainer_log_odds(X_adult[:1000])
```
```
Permutation explainer: 1001it [01:01, 13.61it/s]                                  
```
```
shap.plots.scatter(shap_values_adult_log_odds[:, "Age"])
```
```
# make a standard partial dependence plot
sample_ind = 18
fig, ax = shap.partial_dependence_plot(
    "Age",
    model_adult_log_odds,
    X_adult,
    model_expected_value=True,
    feature_expected_value=True,
    show=False,
    ice=False,
)
```
5.解释非加性提升树逻辑回归模型
```
# train XGBoost model
model = xgboost.XGBClassifier(n_estimators=100, max_depth=2).fit(
    X_adult, y_adult * 1, eval_metric="logloss"
)
# compute SHAP values
explainer = shap.Explainer(model, background_adult)
shap_values = explainer(X_adult)
# set a display version of the data to use for plotting (has string values)
shap_values.display_data = shap.datasets.adult(display=True)[0].values
```
```
The use of label encoder in XGBClassifier is deprecated and will be removed in a future release. To remove this warning, do the following: 1) Pass option use_label_encoder=False when constructing XGBClassifier object; and 2) Encode your labels (y) as integers starting with 0, i.e. 0, 1, 2, ..., [num_class - 1].
 98%|===================| 31839/32561 [00:12
```