解释性差？机器学习解释利器：SHAP原理及实战_shap可解释性

武汉加油 热干面，你要好起来啊！

本文首发于公众号：AI小老弟，全文约5000字，阅读时长5-10分钟

导读
本文首先介绍了机器学习解释包SHAP原理和计算方法，然后基于kaggle竞赛Home Credit数据构建用户违约预测的二分类模型，实战演练了SHAP的几个常用功能。

针对结构化的数据以及分类任务，集成模型往往会有较好的效果，如XGBOOST的诞生，不仅风靡各大数据竞赛，也在工程中得到了广泛的应用。

对于集成学习方法，效果虽好，但一直无法解决可解释性的问题。我们知道一个xgboost或lightgbm模型，是由N棵树组成，所以对于特定的一个样本，我们无法知道这个样本的特征值是如何影响最终结果。虽说“不管白猫黑猫，抓住耗子的就是好猫”，但在具体任务中，我们还是希望能够获得样本每个特征与其结果之间的关系，特别是针对模型误分的那些样本，如果能够从特征和结果的角度进行分析，对于提高模型效果或是分析异常样本，是非常有帮助的，此外，针对某些特定的场景，如风控等，我们甚至可以分析拒绝样本的具体原因，提供给运营部门，完善业务流程。

最重要的是，在做PPT的时候，也会让自己的报告合理并且高大上起来，赢得老板赞赏，走向人生巅峰（皮）。

所以今天，就要介绍一个来自于博弈论的方法–SHAP（SHapley Additive exPlanations），解决上面的问题。

在介绍SHAP之前，首先思考这样一个问题：

_

小明，小军，小强（是的，就是他们小学课本三人组），组队参加王者农药大赛，大赛设定哪个队先拿100个人头，可以获得一万元奖金。终于在他们三个的通力配合下，赢得了比赛，获得一万元奖金，但到了分钱阶段，分歧出现了，因为他们三个人的水平、角色不一，小强个人实力最强善用高输出角色，光是他自己就拿了大半的人头；但若是按照人头数分，也不合适，因为前期小强有几次差点挂掉，多亏队友及时治疗，另外有好多人头也是靠攻速抢到的。那么，应该怎么分配这一万块钱，才是最公平的呢？

_

机智（厚脸皮）的小老弟登场了，这个问题，我们可以从贡献出发：

小强自己可以干掉35个人，小军可以干掉15个人，小明可以干掉10个人，显然，他们每个人都无法独立完成目标。

小强和小军联手可以干掉70个人，小强和小明可以干掉60个人，小军和小明联手可以干掉40个人，同样，无法完成任务。

而只有当三个人联手时候，才可以干掉100个人，达成目标。

（可以理解为小军小强小明有辅助有输出有坦克）

考虑先后顺序，他们三个人一共有6种组队方式。

我们定义边际贡献：假定初始组合为小强一个人，贡献35，加入小军后，他们俩的组合贡献为70，则小军的边际贡献就是70-35=35，再加入小明，三个人的组合贡献为100，则小明的边际贡献为100-70=30；

根据组合次序计算全部的组合及组合中每个人的边际贡献，如表格所示，则可以求得每个人的贡献度。

因此，最合理的分配方案就是小强小军小明分别获得4916、2917、2167。

大家可能奇怪上面的故事和本文主题有啥关系，实际上，上面的问题就是通俗的解释了SHAP方法的核心：shapley value(沙普利值)的简化版计算方法。

沙普利值是博弈论大师劳埃德·沙普利提出来的一种针对合作博弈的解决方案，对公式推导等感兴趣的，可以参考文献【1】，这里就不展开细讲了。

SHAP方法几乎可以给所有机器学习、深度学习提供一个解释的方案，包括树模型、线性模型以及神经网络模型。我们重点关注树模型，研究SHAP是如何评价树模型中的特征对于结果的贡献度。主要参考论文为【2】【3】【4】。

_

对实战更感兴趣的朋友可以直接拖到后面。

_

对于集成树模型来说，当做分类任务时，模型输出的是一个概率值。前文提到，SHAP是SHapley Additive exPlanations的缩写，即沙普利加和解释，因此SHAP实际是将输出值归因到每一个特征的shapely值上，换句话说，就是计算每一个特征的shapley值，依此来衡量特征对最终输出值的影响。用公式表示：

其中，其中g是解释模型，M是输入特征的数目， z表示相应特征是否存在（1或0），这里的存在是针对如图像和文本数据（如文本中，将词one-hot后，某个句子中并不会出现所有词）；Φ是每个特征的归因值（Shapley值）， Φ_0是一个常数。由于树模型的输入是结构化数据，对于样本x，所有的特征都是存在的，因此公式可以写为：

公式中的难点是求j，针对如果是线性模型，

样本有p个特征，β是特征权重，定义第j个特征预测的贡献是

其中E(β_j X__j)是特征j的平均影响估计（即全部样本该特征的均值）。将一个样本的所有特征的贡献相加，则：

也就是说样本 x 的所有特征的贡献之和等于预测值减去平均预测值。当然上述情况是在模型为一个线性模型且特征之间相互独立，在树模型中这种方法显然不成立，因此，对于某个特征j，需要针对所有可能的特征组合（包括不同顺序）计算shapley值，然后加权求和，即：

其中 S是模型中使用的特征的子集，x是要解释的样本的特征值的向量，p为特征数量 ， val(S)指在特征组合S下的模型输出值，解释一下权重，一共p个特征，则在考虑顺序的情况下这p个特征共有p!种组合，固定了某个特征j，则剩余的有(p-| S |-1)!S!种组合。

从上面看，这个方法主要存在两个问题：

1、计算非常耗时，这是因为子集S的数量太多，如果要精确计算每个S下的模型的输出，计算量是非常庞大的。

2、Shapley值有可能会造成误导，这是因为特征值的Shapley值不是从模型训练中删除特征后的预测值之差。 Shapley值的解释是：给定当前的一组特征值，特征值对实际预测值与平均预测值之差的贡献就是估计的Shapley值。

针对这两个问题，Lundberg提出了TreeSHAP，这是SHAP的一种变体，适用于基于树的机器学习模型，例如决策树，随机森林和GBDT。 TreeSHAP速度很快，可以计算精确的Shapley值，并且在特征间存在相关性时正确估计Shapley值。

首先简单介绍下计算原理，假定只有一棵树，针对某个样本x和特征子集S。如果考虑全部特征，即子集S为全部特征集合，则叶子节点的预测值就是该特征下的模型的输出值；如果子集为空，则使用所有叶子节点的预测值的加权平均作为输出值；如果子集包含部分特征， 则忽略那些由于去掉部分特征后所无法触达的叶子节点，从剩下的叶子节点中，再取加权（权重，指的是叶子节点的大小，也就是叶子节点上的样本的数量）取平均。这里存在一个问题就是计算一个特征的shapley值的时候，需要把这一过程应用到所有的可能的特征集S上，TreeSHAP方法的解决思路是一次性把所有可能的子集S沿树展开，对于每一个决策节点，我们必须关注子集的大小，而这取决于在父节点以及分割特征上的子集。举例来说，一棵树的首个分割决策点为特征x3，所有包含特征x3的特征子集将划到一个节点上；而不包含x3的特征子集，将会划到降低权重后的节点上。但是不同的子集大小不同，意味着权重也就不同。算法必须持续关注每一个节点上所有子集的权重（大小），而这也增加了算法的复杂度。针对这一问题，TreeSHAP也有解决办法， 下面将继续介绍。

首先在O(TL2^M)时间内直接估计SHAP值

这里不考虑计算复杂度；计算某特征的SHAP值，我们一共有M个特征，则可能的S子集有2^{M个，时间复杂度为O(2}M)，而对于每一个S子集，需遍历所有树的所有节点去计算最后每个叶子节点的预测期望值，因此每个子集S在所有子树中的时间复杂度是O(TL)，综合起来，针对所有的特征子集S计算的时间复杂度是 O(TL2^M)。

如图，其中，其中v是节点值的向量，而内部节点的值是为internal，向量a和b是每个内部节点的左右节点索引，向量t为每个内部节点的阈值，d是在内部节点中分割的特征的索引向量，向量r表示每个节点的覆盖的样本数量（即有多少数据样本落在该节点中），权重w表示与子集S匹配的训练样本落在节点的比例。

算法1具体流程：从树的根节点开始执行程序G(j, w)，即计算G(1, 1)，如果不是内部节点，即为叶子节点，则直接计算（叶子节点预测值 * 匹配集合S的训练样本落在新节点的比例）；如果当前节点是内部节点，判断该内部节点是否属于子集S，如果属于子集S则按照当前节点的阈值分割，在新的节点计算G，如果不属于子集S则根据左右节点的样本比例，计算两个新节点的G的和。

在O(TLD^2)时间内直接估计SHAP值

算法1的时间复杂度较高，TreeSHAP提出了一种优化算法，将时间复杂度从指数降低为多项式。具体来说，时间复杂度降为 O(TLD2)，空间复杂度为O(D2 + M)，同时在树为平衡树的时候，深度D变成logL。其中L是树中的最大叶子数量，M是特征数量。

上面算法中的m，是已经分裂的特定特征的路径，它包含四个属性：

d：特征索引；

z：流经该分支的“0”（即该特征没在特征子集S中）路径的分数；

o：流经这个分支的“1”路径(即该特征在特征子集S中)的分数；

w：用于保存给定基数集合的比例。

具体算法解释，这里就不展开了，可以参考论文【2】。

至此，TreeSHAP介绍完毕。

实际上，SHAP不止有TreeSHAP方法，针对深度学习等，SHAP目前也有一些解决框架，SHAP方法的优良特性，将为模型的解释性上提供合理的方案，辅助人们更好的去认识机器学习和深度学习。

_

实战演示部分开始

_

首先开始SHAP的安装和应用，常规操作

pip install shap
OR
conda install shap  #针对anaconda用户
1
2
3

我们在jupyter notebook 里调用，可视化操作方便。

本例数据来自kaggle竞赛Home Credit Default Risk，也就是网贷风险用户识别，为一个二分类的模型。小老弟的机器性能太渣，所以只能选取20个特征，10000个样本进行测试，模型为Lightgbm，模型参数并没有做太多调整。需要注意一个地方，数据中有类别变量，在Lightgbm中可以首先将其转化为int型，然后直接指定类别变量，在利用SHAP展示的时候，可以准备两份数据，一份是原始的数据，另一份是将类别变量转化为int数字的数据，如下面调用SHAP的代码中，data_model为转化后的数据，而data为原始数据。

首先看一下LGB自带的特征重要性，这里选择importance_type=‘gain’，也就是通过模型训练过程中，分割节点时特征的总“gain”来计算特征的重要性排序。

下面来使用SHAP

explainer = shap.TreeExplainer(model) # 初始化解释器
shap.initjs() #初始化JS
shap_values = explainer.shap_values(data_model[use_cols]) #计算每个样本的每个特征的SHAP值
1
2
3

接下来就可以查看某个样本的SHAP情况，我们选择三个样本，其中两个预测正确的样本，即正例预测为正例，反例预测为反例，另外一个为高风险未识别样本。

查看三个样本的shaple值，其中红色代表增加shaple值，蓝色代表降低，也即红色代表增加样本风险（下同）。

从上之下依次是低风险，高风险，未识别三个特定样本的shap值可视化图。

shap.force_plot(explainer.expected_value[1], shap_values[1][3860,:], data[use_cols].iloc[3860,:]) #3860为样本在数据集中的索引
1

上图中有一个base value和output value，base value，也就是所有样本的平均预测水平，SHAP对于xgboost和lightgbm取的是对数几率转换，也就是说，

可以做一个验证：

explainer.expected_value[1],  np.log(data['PREDICT']/ (1 - data['PREDICT'])).mean()
#output:(-2.583467142969475, -2.583400390403734)
1
2

总体来说并不影响我们使用SHAP，继续看图中高风险样本，output value 是要高于base value的，也就是说这个样本的风险高于整体水平，而高的0.51-（-2.58）=3.09则是由不同特征贡献的。

shap_values[1][6578,:].sum()
# output:3.094375265491124
1
2

我们再看一个异常样本，即标签为1，但模型预测概率很低的样本，相当于未识别出的坏样本，对比三个样本的特征shap值，从上至下依次是低风险样本，高风险样本和未识别样本。

shap.bar_plot(shap_values=shap_values[1][3860,:],feature_names=use_cols)
1

可以看到，未识别样本的各特征贡献上与低风险样本类似，这也是造成模型误判的原因。

再来看概括图，即 summary plot，该图是对全部样本全部特征的shaple值进行求和，可以反映出特征重要性及每个特征对样本正负预测的贡献。

shap.summary_plot(shap_values, data[use_cols])
1

第二种summary_plot图，是把所有的样本点都呈现在图中，如图，此时颜色代表特征值的大小，而横坐标为shap值的大小，从图中可以看到 days_credit这一特征，值越小，shap值越大，换句话来说就是days_credit越大，风险越高。

shap.summary_plot(shap_values[0], data[use_cols])
1

进一步，如果我们想看单个特征与预测结果之间的关系，那么可以使用dependence plot，dependence_plot有两种模式，一种是单个特征与模型预测结果的关系，另一种是两个特征交互作用下的与模型预测结果的关系。

首先是DAYS_CREDIT单特征：

shap.dependence_plot('DAYS_CREDIT', shap_values[1], data_model[use_cols], display_features=data[use_cols],interaction_index=None)
1

如图，横坐标是特征的值，纵坐标是对应的shap值，整体趋势也是随着days_credit增大，shap值增大。

我们再来看与NAME_INCOME_TYPE特征交互影响下的图。

shap.dependence_plot('DAYS_CREDIT', shap_values[1], data_model[use_cols], display_features=data[use_cols],interaction_index='NAME_INCOME_TYPE')
1

如图，主体仍然是DAYS_CREDIT的分布，但是加入了NAME_INCOME_TYPE特征，右坐标轴代表不同的NAME_INCOME_TYPE的颜色，可以看到，在不同的DAYS_CREDIT区间上，由不同的NAME_INCOME_TYPE带来的shap值也不同，如在较小的DAYS_CREDIT上，NAME_INCOME_TYPE为working有更低的shap值，而在中间的的DAYS_CREDIT上，working有更高的shap值。

我们也可以转换一下坐标轴，以name为主体。

shap.dependence_plot('NAME_INCOME_TYPE', shap_values[1], data_model[use_cols], display_features=data[use_cols],interaction_index='DAYS_CREDIT')
1

可以看到，收入类型为working的shap值整体要高不少。

至此，SHAP的基本使用介绍完毕，还有其他一些功能，大家也可以踊跃探索。

参考文献

[1] Shapley, L.S… (1953). A value for n-persons games. Annals of Mathematics Studies. 28. 307-318.[2] Lundberg, Scott & Lee, Su-In. (2017). Consistent feature attribution for tree ensembles.

[3] Scott M. Lundberg and Su-In Lee. 2017. A unified approach to interpreting model predictions.

[4] Marco Tulio Ribeiro, Sameer Singh, and Carlos Guestrin. 2016. “Why Should I Trust You?”: Explaining the Predictions of Any Classifier.
转载：https://zhuanlan.zhihu.com/p/106320452
欢迎关注作者主页，学习更多相关内容

在公众号「python风控模型」里回复关键字：学习资料