shap计算交互效应

字数统计: 2.4k   |   阅读时长≈ 10 分钟

本笔记展示了如何计算一个非常简单函数的SHAP交互值。我们从一个简单的线性函数开始，然后添加一个交互项，观察它如何改变SHAP值和SHAP交互值。

首先导入需要用到的包

1
2
3
4
5
import numpy as np
import xgboost
from sklearn.linear_model import LinearRegression

import shap

没有互作的例子

模拟一些二元数据和一个包含交互项的线性结果。

请注意，我们使 X 中的各个特征完全独立，以便轻松求解精确的 SHAP 值。

这里 X 和 y 分别是输入和输出数据，无互作效应

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
N = 2_000
X = np.zeros((N, 5))

X[:1_000, 0] = 1

X[:500, 1] = 1
X[1_000:1_500, 1] = 1

X[:250, 2] = 1
X[500:750, 2] = 1
X[1_000:1_250, 2] = 1
X[1_500:1_750, 2] = 1

# mean-center the data
X[:, 0:3] -= 0.5

y = 2 * X[:, 0] - 3 * X[:, 1]

我们看到这些变量确实是独立的

1
2
3
4
5
6
7
8
In [3]: np.cov(X.T)
Out[3]:
array([[0.25012506, 0.        , 0.        , 0.        , 0.        ],
       [0.        , 0.25012506, 0.        , 0.        , 0.        ],
       [0.        , 0.        , 0.25012506, 0.        , 0.        ],
       [0.        , 0.        , 0.        , 0.        , 0.        ],
       [0.        , 0.        , 0.        , 0.        , 0.        ]])

并且每列均值均为0

1
2
In [4]: X.mean(axis=0)
Out[4]: array([0., 0., 0., 0., 0.])

构建一个简单的 xgboost 模型

1
2
3
4
5
6
# train a model with single tree
Xd = xgboost.DMatrix(X, label=y)
model = xgboost.train({"eta": 1, "max_depth": 3, "base_score": 0, "lambda": 0}, Xd, 1)

print("Model error =", np.linalg.norm(y - model.predict(Xd)))
print(model.get_dump(with_stats=True)[0])

计算 shap 值

1
2
3
4
5
6
7
8
pred = model.predict(Xd, output_margin=True)

explainer = shap.TreeExplainer(model)
explanation = explainer(Xd)

shap_values = explanation.values
# make sure the SHAP values add up to marginal predictions
np.abs(shap_values.sum(axis=1) + explanation.base_values - pred).max()

如果我们绘制蜂群图，就会发现只有特征0和1对输出有影响，且它们的影响仅有两种可能的幅度（分别为1/-1和1.5/-1.5）。

训练一个线性模型

1
2
3
4
5
In [10]: lr = LinearRegression()
    ...: lr.fit(X, y)
    ...: lr_pred = lr.predict(X)
    ...: lr.coef_.round(2)
Out[10]: array([ 2., -3., -0.,  0.,  0.])

确保计算出的SHAP值与真实的SHAP值一致（对于任何线性回归，我们都可以直接计算真实的SHAP值）

1
2
3
In [11]: main_effect_shap_values = lr.coef_ * (X - X.mean(0))
    ...: np.linalg.norm(shap_values - main_effect_shap_values)
Out[11]: np.float64(1.6542433490447965e-13)

这个极小的数值证明了：在没有特征交互的情况下，树模型的 SHAP 值精确退化为线性模型的主效应，SHAP 方法具有一致性。

计算shap互作效应

请注意，当不存在交互作用时，SHAP交互值仅为对角矩阵，其对角线上的元素为SHAP值。

1
2
shap_interaction_values = explainer.shap_interaction_values(Xd)
shap_interaction_values[0]

第一行是计算交互效应

  • explainer 是之前创建的 shap.TreeExplainer(model)。
  • shap_interaction_values(Xd) 是 TreeExplainer 的一个方法，专门用于计算 SHAP 交互值矩阵。
  • 返回值的形状是 (n_samples, n_features, n_features)：
    • 第一个维度：样本数。
    • 第二、三个维度：特征 × 特征的方阵，每个样本一个矩阵。
  • 计算复杂度比普通 SHAP 值高，因为要枚举所有特征对，但 TreeExplainer 对树模型有精确的多项式时间算法。

第二行输出第一个样本的结果

1
2
3
4
5
6
Out[14]:
array([[ 1. ,  0. ,  0. ,  0. ,  0. ],
       [ 0. , -1.5,  0. ,  0. ,  0. ],
       [ 0. ,  0. ,  0. ,  0. ,  0. ],
       [ 0. ,  0. ,  0. ,  0. ,  0. ],
       [ 0. ,  0. ,  0. ,  0. ,  0. ]], dtype=float32)
  • 取出第一个样本的交互值矩阵，形状为 (n_features, n_features)。
  • 这个矩阵的结构：
    • 对角线元素 [i, i]：特征 i 的主效应（main effect）。
    • 非对角线元素 [i, j]（i≠j）：特征 i 和特征 j 之间的交互效应（）。
  • 矩阵是对称的，即 [i, j] = [j, i]，因为交互是双向的。

让我们确保SHAP交互值的总和等于边际预测：这里就是检验每个样本的交互值矩阵中所有元素加在一起 + 基准值 ，是否精确等于模型的原始预测值。

1
np.abs(shap_interaction_values.sum((1, 2)) + explainer.expected_value - pred).max()

并确保SHAP交互作用值的主要效应与线性模型中的效应相匹配：

1
2
3
4
5
total = 0
for i in range(N):
    for j in range(5):
        total += np.abs(shap_interaction_values[i, j, j] - main_effect_shap_values[i, j])
total
1
Out[17]: np.float64(1.0533118387982904e-11)

额外：详细解释交互效应

假设特征 1 和特征 2 之间存在交互。在 shap_interaction_values 矩阵中：

1
2
3
       特征1   特征2
特征1 [ Φ[1,1]  Φ[1,2] ]
特征2 [ Φ[2,1]  Φ[2,2] ]

由于矩阵对称，Φ[1,2] = Φ[2,1]。

这里有两个概念：

名称 表达式 含义
单个特征分到的交互贡献 Φ[1,2]（或 Φ[2,1]） 特征 1 从交互中分到的部分，也是特征 2 分到的部分
两个特征之间的总交互效应 Φ[1,2] + Φ[2,1] = 2 * Φ[1,2] 交互整体对预测值的额外贡献

所以：

非对角线元素 Φ[1,2] 是“总交互效应”的一半，但它是“特征 1 分到的交互贡献”的全部。

不矛盾，只是参照物不同。

额外： shap_values 和 shap_interaction_values 矩阵的关系

每个特征的 shap_values 和 shap_interaction_values 的关系如下

对于每个样本、每个特征 i： \[ \boxed{\text{shap_values}[i] = \text{shap_interaction_values}[i,i] + \sum_{j \neq i} \text{shap_interaction_values}[i,j]} \] 用文字说就是：

普通 SHAP 值 = 主效应（对角线） + 该特征与其他所有特征的交互贡献之和

所以 shap_values 就是 shap_interaction_values 矩阵第 i 行的行和（也是第 i 列的列和，因为矩阵对称）。

有互作的例子

创建数据如下

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
N = 2_000
X = np.zeros((N, 5))
X[:1_000, 0] = 1

X[:500, 1] = 1
X[1_000:1_500, 1] = 1

X[:250, 2] = 1
X[500:750, 2] = 1
X[1_000:1_250, 2] = 1
X[1_500:1_750, 2] = 1

X[:125, 3] = 1
X[250:375, 3] = 1
X[500:625, 3] = 1
X[750:875, 3] = 1
X[1_000:1_125, 3] = 1
X[1_250:1_375, 3] = 1
X[1_500:1_625, 3] = 1
X[1_750:1_875, 3] = 1

# we can't exactly mean center the data or XGBoost has trouble finding the splits
X[:, :4] -= 0.4999

# interaction of features is implemented as the multiplication of the features. Note that any other function of the
#  features would also work, but is harder to interpret (e.g. sin(x1*x2)).
y = 2 * X[:, 0] - 3 * X[:, 1] + 2 * X[:, 1] * X[:, 2]

查看每列均值是否为0

1
2
In [19]: X.mean(axis=0)
Out[19]: array([1.e-04, 1.e-04, 1.e-04, 1.e-04, 0.e+00])

训练一个简单的 xgboost 模型

1
2
3
4
5
# train a model with single tree
Xd = xgboost.DMatrix(X, label=y)
model = xgboost.train({"eta": 1, "max_depth": 4, "base_score": 0, "lambda": 0}, Xd, 1)
print("Model error =", np.linalg.norm(y - model.predict(Xd)))
print(model.get_dump(with_stats=True)[0])

计算 shap 值

1
2
3
4
5
6
7
8
pred = model.predict(Xd, output_margin=True)

explainer = shap.TreeExplainer(model)
explanation = explainer(Xd)

shap_values = explanation.values
# make sure the SHAP values add up to marginal predictions
np.abs(shap_values.sum(axis=1) + explanation.base_values - pred).max()
1
Out[21]: np.float32(4.7683716e-07)

如果我们构建一个蜂群图，就会发现现在只有特征3和4无关紧要，而特征1由于交互作用，可能有四种不同的效应大小。

1
shap.plots.beeswarm(explanation)

训练一个线性模型

1
2
3
4
lr = LinearRegression()
lr.fit(X, y)
lr_pred = lr.predict(X)
lr.coef_.round(2)
1
Out[23]: array([ 2., -3.,  0., -0.,  0.])

请注意，SHAP 值不再与主效应一致，因为它们现在包含了交互效应。

1
2
main_effect_shap_values = lr.coef_ * (X - X.mean(0))
np.linalg.norm(shap_values - main_effect_shap_values)
1
Out[24]: np.float64(15.811387829626835)

计算shap互作效应

同样代码如下

1
2
shap_interaction_values = explainer.shap_interaction_values(Xd)
shap_interaction_values[0].round(2)
1
2
3
4
5
6
7
Out[26]:
array([[ 1.  ,  0.  ,  0.  ,  0.  ,  0.  ],
       [ 0.  , -1.5 ,  0.25,  0.  ,  0.  ],
       [ 0.  ,  0.25,  0.  ,  0.  ,  0.  ],
       [ 0.  ,  0.  ,  0.  ,  0.  ,  0.  ],
       [ 0.  ,  0.  ,  0.  ,  0.  ,  0.  ]], dtype=float32)

我看了一下不算互作效应时的 shap 值如下，可以看到确实是： shap_values 就是 shap_interaction_values 矩阵第 i 行的行和

1
2
3
4
In [28]: shap_values[0]
Out[28]:
array([ 1.        , -1.2499001 ,  0.25010002,  0.        ,  0.        ],
      dtype=float32)

确保SHAP交互值之和等于边际预测

1
np.abs(shap_interaction_values.sum((1, 2)) + explainer.expected_value - pred).max()
1
Out[29]: np.float32(4.7683716e-07)

当存在交互作用时，主效应不再与 SHAP 值匹配，但它们确实与线性模型的主效应在 SHAP 交互值矩阵的对角线上匹配。

1
2
3
4
5
total = 0
for i in range(N):
    for j in range(5):
        total += np.abs(shap_interaction_values[i, j, j] - main_effect_shap_values[i, j])
total
1
Out[30]: np.float64(0.0005347490392160024)

绘制SHAP 依赖图（dependence plot） ：作用是可视化单个特征对模型预测的影响，并揭示该特征是否与其他特征存在交互作用。

1
shap.dependence_plot(0, shap_values, X)
  • ind：要分析的特征。可以是特征索引（整数）或特征名称（字符串）。这里传 0，表示分析第 0 个特征（即 X 的第一列）。
  • shap_values：SHAP 值矩阵，形状为 (n_samples, n_features)。通常由 explainer.shap_values(X) 或 explanation.values 得到。
  • features：原始特征数据矩阵 X，形状为 (n_samples, n_features)。函数会从中提取第 ind 列作为 x 轴，并自动寻找与第 ind 个特征交互最强的另一个特征用于着色。

调用后会生成一张散点图：

  • x 轴：第 ind 个特征的特征值（来自 X[:, ind]）。
  • y 轴：该特征的 SHAP 值（来自 shap_values[:, ind]），表示该特征对每个样本预测值的贡献。
  • 每个点：一个样本。
  • 颜色：默认情况下，函数会自动选择与第 ind 个特征交互作用最强的另一个特征，用它的值给点上色。颜色深浅反映那个特征的大小。
  • 趋势线：一条平滑曲线，展示 SHAP 值随特征值变化的平均趋势。

这里我们为特征0绘制依赖关系图，就会发现它只需要两个值，且这些值完全取决于该特征的取值。因此，它们位于一条直线上（因为特征0没有与其他特征发生交互，所以其效果完全由该特征的取值决定）。

相反，如果我们为特征2构建一个依赖关系图，就会发现它有4个可能的取值，而这些取值并非完全由特征2的值决定，它们还取决于特征3的值。依赖关系图中的这种垂直分布反映了非线性交互的影响。

1
shap.dependence_plot(2, shap_values, X)

参考资料

  1. https://shap.readthedocs.io/en/latest/example_notebooks/tabular_examples/tree_based_models/Basic%20SHAP%20Interaction%20Value%20Example%20in%20XGBoost.html
  • 版权声明： 本博客所有文章除特别声明外，著作权归作者所有。转载请注明出处！
  • Copyrights © 2019-2026 Vincere Zhou
  • 访问人数: | 浏览次数:

请我喝杯茶吧~

支付宝
微信