基于 shapley-value 的可解释方法及 shap 库实现

一. 为什么需要可解释

几点考虑:

1. 若模型完全黑箱, 会有信任风险, 虽然 Performance 不错, 但在医学诊断等严肃领域, 同样要关心诊断依据.
2. 人类天生的好奇心, 也想知道不同特征到底作了怎样的贡献.
3. 对模型预测的 badcase 作诊断, 增强洞察, 辅助模型与特征的迭代.

模型可解释其实就是想弄懂不同特征到底做了怎样的贡献, 从解释粒度上可以这么分类:

1. over the whole set
   这种好理解, 求一个特征全局视角下的贡献, 可以直接从训练集中把它拎出来, 观察模型收敛后的指标变化.
2. for a particular prediction
   深度模型效果更好的原因之一就是拥有强大的非线性拟合能力, 也就是说同一个特征下的同一个特征值, 会随样本的变化而体现出不同的贡献(受同一个样本内其他邻居特征的影响), 这就导致了仅有全局解释是不够的, 粒度需要细化到具体的单次预测上.
   当能做到了单次预测, 自然也就能统计一个样本集合上的可解释数据了, 此时也拥有了全局视角. 所以本文工作搞单样本的可解释.

形式化定义

形式化的单样本可解释任务就是:
$$g\left(z^{\prime }\right)={\varphi }_0+\sum _{i=1}^M{\varphi }_i{z}_i^{\prime }\text{\hspace{1em}(1)}$$
where z ′ ∈ { 0 , 1 } M z' ∈ \{0, 1\}^M z′∈{0,1}M, M is the number of simplified input features, and ${\varphi }_i\in \mathbb{R}$.
式(1)解读为对出现在待解释样本中的每个特征, 分配一个量化的贡献度.

本文主角是 shap 库, 预备知识是 联盟博弈论的 shapley-value.

二. 联盟博弈论中的 shapley-value

合作博弈论用于多人合作下的收益分配. 主要思想是: 列举出各种不同玩家之间的合作情况, 依据玩家参与与否的边际效应计算贡献.

直接用 shap 库中的公式了. 它是 classic Shapley value equation.
ϕ i = ∑ S ⊆ F \ { i } ∣ S ∣ ! ( ∣ F ∣ − ∣ S ∣ − 1 ) ! ∣ F ∣ ! [ f S ∪ { i } ( x S ∪ { i } ) − f S ( x S ) ] (2) \phi_i=\sum_{S \subseteq F \backslash\{i\}} \frac{|S| !(|F|-|S|-1) !}{|F| !}\left[f_{S \cup\{i\}}\left(x_{S \cup\{i\}}\right)-f_S\left(x_S\right)\right] \tag 2 ϕi​=S⊆F\{i}∑​∣F∣!∣S∣!(∣F∣−∣S∣−1)!​[fS∪{i}​(xS∪{i}​)−fS​(xS​)](2)

其中,

• F, 所有玩家的集合
• i, 待评估贡献的单个玩家
• S, 剩下玩家组成集合的所有可能的子集
• f, 计算合作收益的函数
• $x_S$, 表示一次合作中参与的玩家组合为集合S.

公式解读

式(2)中, 左边的分数部分为系数, 右边的中括号部分为边际收益.
Q: 这个系数(权重)怎么理解呢?
A: 网上资料没看到很好的解读, 说下个人见解吧: 先不管减去一个1这个细节, 这个系数其实就是组合数 $C_{|F|}^{|S|}$ 的倒数.
$\Sigma$ 符号相当于一层 for 循环, 当 S 确定时, 它的出现概率解读为 从 F 中挑 |S| 个特征的种类数, 而当前的 S 又是其中确定的一种, 所以概率就是 $\frac{1}{C_{|F|}^{|S|}}$.

直观例子

参考[6] 中的知乎文章给了一个很好的例子, 这里二次加工下.

甲、乙、丙三人合作经商。不同组合的不同收益列举如下, 问三人合作时如何利益分配?

• 空集
• 甲乙丙单人, 对应三种情况, 都是1万.
• 甲乙, 7万
• 甲丙, 5万
• 乙丙, 4万
• 甲乙丙, 11万

三个元素的所有子集个数为 $2^3=8$, 已经列举在上面了.
代入公式, 见下图(符号有出入, 大体一样)

所以 甲的最终获利为 $\frac{13}{3}$.

计算复杂度

因为要枚举所有玩家的所有不同组合情况, 即集合 F 的所有子集, 这个组合数就是 $2^{|F|}$, 属于 NP-Hard 问题.
怎么去用复杂度更低的计算方法去近似它呢?

1. Shapley sampling values, 见参考 [8].
2. Kernal Shap method, 见参考[4], 下文就作展开.

三. LIME 与 Kernel-Shap

在讲 Kernel-Shap 之前, 先引入预备知识 LIME.

LIME

Lime, Local Interpretable Model-Agnostic Explanations, 是 2017 年提出的一种模型无关的单样本预测解释方法(详见参考[9]). 核心思想是对于要解释的样本
, 把其特征划分为若干个组, 然后对此作局部采样, 得到新的人造样本集合 S. 然后训练新的线性代理模型去学习原模型在样本S附近上的输出.

公式化的表达就是:
$$\xi (x)=\underset{g\in G}{\mathrm{argmin}}L\left(f,g,{\pi }_x\right)+\Omega (g)\text{\hspace{1em}(3)}$$
$g(\cdot )$就是式(1)中的可解释函数, ${\pi }_x$ 是人造样本的权重, $\Omega (\cdot )$限制着g的复杂程度.

图 使用 lime 可解释洞察到模型误把哈士奇(Husky) 预测为狼(Wolf)的原始是只关注到了雪地(snow)这一因素

上图是 lime 论文中的一个素材. 图片的原始特征是像素, 它就通过
选取图像区域 作 super-pixel 的集合映射. 通过可解释分析, 验证了模型学到了预期之外的狼与雪地之间的相关性, 而没有学到因果层面的相关性.

Kernel-Shap

Kernal Shap method 是 LIME + Shapley values 的结合体. 乍一看式(2)与式(3)相去甚远, 但论文中讲到

1. LIME 是一种特征可加性方法(additive feature attribution method);
2. 可加性方法中, 存在满足 {Local accuracy, Missingness, Consistency} 三大特性的唯一解;
3. 而 shap-values 又是符合 LIME 方程约束下的同时具有上句提到的三大特性的唯一解.

就这么一通牵扯, 将二者结合到了一起.
LIME 呢, 可以解读为一种范式而不是确切的一个算法, 因为 loss function $L$ 与 weighting kernel ${\pi }_x$ , 还有 regularization term $\Omega$ 的选取都是缺乏指导的, 而 Kernel Shap 将其作了具化, 满足了上文提到的三大特性.

图. 截取自原论文.

additive 与 三大特性

SHAP 是一类 additive feature attribution (满足可加性的特征归因) 方法. 该类方法更好地满足三大可解释性质:

1. local accuracy
   $$f(x)=g(x^{\prime })={\varphi }_0+\sum _{i=1}^M{\varphi }_i{x}_i^{\prime }\text{\hspace{1em}(1)}$$
   各 feature value 的加和 = 该 sample 的 model output
2. missingness
   $$x_i^{\prime }=0\Rightarrow {\varphi }_i=0\text{\hspace{1em}(2)}$$
   缺失值的 feature attribution value = 0
3. consistency
   当模型有变化, 一个特征变得更重要时, 其 feature attribution value 也不能变小.

四. Shap 库介绍及KernelSHAP实现

SHapley Additive exPlanation, 是一个py 三方库, 依据 合作博弈论 领域 中的 shapley value 思想, 对模型的单个预测作解释.
它把特征比作博弈问题中的玩家, 模型预测比喻玩家合作之后的收益, 于是就顺畅了.

论文[4] 中, 它是这么说的:
“已有的多种方法 LIME, DeepLift 等, 它们之间的联系是啥? 什么情况下, 用其中一种会比另一种更好用?” 难以回答, 而它引入了 shapley 思想作了统一, 既有计算效率上的优化, 又更符合人类直觉.

基于 KernelSHAP 的实现的解释器叫 KernalExplainer.

基本流程

1. 传入待解释样本, 计算待解释的特征个数M
2. 构造不同组合的合成(人造)样本
3. 维护 maskMatrix, 样子见下:
   
4. 计算合成样本的输出
   通过调用 run() 实现.
5. 求解重要性
   这里其实我不太懂. 原公式只是列出了单个特征的重要性求解方法. 但在 solve() 方法实现中, 是批量一次性求解的. 这里面的等价推导shap库没有给出讲解.

源码解读

作了精简, 解读就在注释中.

class Kernel(Explainer):
    def __init__(self, model, data, link=IdentityLink(), **kwargs):
        """
        model: callable对象, 传入特征即可输出结果.
        data: 用来估计数据集上的 期望, 随便传也不会影响一次可解释中的特征重要性
        """
        pass

    def shap_values(self, X, **kwargs):
        """ 主方法 
        X 是待解释的样本, 单个样本解释时, shape 就是 (M,)
        """
        # shape 恢复成了 [1,M]
        data = X.reshape((1, X.shape[0]))
        return self.explain(data, **kwargs)

    def explain(self, incoming_instance:np.ndarray, **kwargs):
        """
        incoming_instance, 是待解释的样本特征
        """
        # 自定义对象, instance.x = incoming_instance
        instance = convert_to_instance(incoming_instance)
        # 哪些特征是可变的, 用下标记录下来, 通常就是 从0到 incoming_instance.shape[1]
        self.varyingInds = self.varying_groups(instance.x)
        self.varyingFeatureGroups = [self.data.groups[i]
                                     for i in self.varyingInds]
        # 即 incoming_instance.shape[1]
        self.M = len(self.varyingFeatureGroups)
        # 原样本(非人造) 的预测结果, 就是待解释的结果
        self.fx = self.model.f(instance.x)
        # 为了控制计算复杂度, 这里限定总的人造样本数
        self.nsamples:int = kwargs.get('nsamples')
        self.allocate()

        # M-1 是因为要把待评估的单个特征摘出来, 除以2 是考虑到成对关系
        num_subset_sizes = np.int(np.ceil((self.M - 1) / 2.0))
        num_paired_subset_sizes = np.int(np.floor((self.M - 1) / 2.0))

        # 与论文 page-6 的 \pi x 权重相对应, 此时还没有 除以 C_M^{subset\_size}
        # weight_vector.shape = (num_subset_sizes,)
        # weight_vector[i] 表示 子集大小为 i+1 时的那些样本的权重之和
        weight_vector = np.array([(self.M - 1.0) / (i * (self.M - i)) for i in range(1, num_subset_sizes + 1)])
        weight_vector[:num_paired_subset_sizes] *= 2
        weight_vector /= np.sum(weight_vector)
        mask = np.zeros(self.M)
        remaining_weight_vector = copy.copy(weight_vector)
        num_samples_left = self.nsamples

        for subset_size in range(1, num_subset_sizes + 1):
            subsets_cnt_of_current_subset_size = binom(self.M, subset_size)
            if num_samples_left * remaining_weight_vector[subset_size - 1] / subsets_cnt_of_current_subset_size >= 1.0 - 1e-8:
                num_full_subsets += 1
                num_samples_left -= subsets_cnt_of_current_subset_size
                for inds_to_be_masked in itertools.combinations(group_inds, subset_size):
                    mask[:] = 0.0
                    mask[np.array(inds_to_be_masked, dtype='int64')] = 1.0
                    self.addsample(instance.x, mask, w)
                    # C_n^k 中的一种情况, 必然与 C_n^{n-k} 的一种情况是 取反的配对关系
                    # 所以可以一次加 俩 样本
                    if subset_size <= num_paired_subset_sizes:
                        mask[:] = np.abs(mask - 1)
                        self.addsample(instance.x, mask, w)
                    else:
                        logger.info(f"对于当前的 subsize={subset_size}, 在给定的总 nsamples={self.nsamples}"
                                    f" 约束下已经不能完全枚举作预测了. 作 跳出 动作")
                        break                        
        
        # 在给定的总样本量约束下, 若 未能全部展开计算, 需要随机采样
        if num_full_subsets != num_subset_sizes:
            # choice() 方法用法, 从 [0, a) 中完成 size 个抽样, 会有重复(有放回抽样)
            # size 为 [1,num_full_subsets) 的 subset 的样本已经添加完毕, 现在只添加后面的
            subset_size_choice = np.random.choice(a=list(range(num_full_subsets+1,
                                                            len(weight_vector)+1)),
                                        # 4倍是为了给去重留 buffer
                                        size=4 * samples_left,
                                        p=remaining_weight_vector)
            ind_set_pos = 0

            while samples_left > 0 and ind_set_pos < len(ind_set):
                mask.fill(0.0)
                # we call np.random.choice once to save time and then just read it here
                subset_size = subset_size_choice[ind_set_pos]
                ind_set_pos += 1
                # 这三行讲 怎么确定 subset_size 下的 subset, 具体哪几个特征其 mask=1
                random_enable_index = np.random.permutation(self.M)
                index_arr = random_enable_index[:subset_size]
                mask[index_arr] = 1.0
                # 虽然这里 w 是 1.0,但不起作用, 因为后面 还要去改 self.kernelWeights[nfixed_samples:]
                self.addsample(instance.x, mask, w=1.0)
            # 就是这里会 覆盖掉上面的 w=1.0
            self.kernelWeights[nfixed_samples:] *= weight_left / self.kernelWeights[nfixed_samples:].sum()


        self.run()
        phi = self.solve()
        return phi

    def allocate(self):
        """ 初始化 mask矩阵, 权重矩阵, 人造样本的预测结果数组 和 人造样本特征.
        """
        self.maskMatrix = np.zeros((self.nsamples, self.M))
        self.kernelWeights = np.zeros(self.nsamples)
        self.y = np.zeros((self.nsamples, self.D))
        self.synth_data = np.tile(self.data.data, (self.nsamples, 1))

    def addsample(self, x, m, w):
        mask = m == 1.0
        evaluation_data = x[0, groups]
        self.synth_data[offset:offset + self.N, groups] = evaluation_data
        self.maskMatrix[self.nsamplesAdded, :] = m
        self.kernelWeights[self.nsamplesAdded] = w
        self.nsamplesAdded += 1
    
    def run(self):
        data = self.synth_data[self.nsamplesRun * self.N:self.nsamplesAdded * self.N, :]
        modelOut = self.model.f(data)
        self.y[self.nsamplesRun * self.N:self.nsamplesAdded * self.N, :] = np.reshape(modelOut, (num_to_run, self.D))

    def solve(self, fraction_evaluated, dim):
        w_aug = np.hstack((self.kernelWeights * (self.M - s), self.kernelWeights * s))
        w_sqrt_aug = np.sqrt(w_aug)
        eyAdj_aug = np.hstack((eyAdj, eyAdj - (self.link.f(self.fx[dim]) - self.link.f(self.fnull[dim]))))
        eyAdj_aug *= w_sqrt_aug
        mask_aug = np.transpose(w_sqrt_aug * np.transpose(np.vstack((self.maskMatrix, self.maskMatrix - 1))))
        coef_ = LassoLarsIC(criterion=c).fit(mask_aug, eyAdj_aug).coef_
        nonzero_inds = np.nonzero(coef_)[0]

        # eliminate one variable with the constraint that all features sum to the output
        eyAdj2 = eyAdj - self.maskMatrix[:, nonzero_inds[-1]] * (
        self.link.f(self.fx[dim]) - self.link.f(self.fnull[dim]))

        # solve a weighted least squares equation to estimate phi
        tmp = np.transpose(np.transpose(etmp) * np.transpose(self.kernelWeights))
        
        w = np.dot(tmp2, np.dot(np.transpose(tmp), eyAdj2))
        phi = np.zeros(self.M)
        phi[nonzero_inds[:-1]] = w
        phi[nonzero_inds[-1]] = (self.link.f(self.fx[dim]) - self.link.f(self.fnull[dim])) - sum(w)

五. shap 库其他 Explainer

LinearExplainer

用于 逻辑回归 模型的可解释, 背后算法是 DeepLIFT algorithm (Deep SHAP) , 官网的例子见 参考[2].
Note that with a linear model the SHAP value for feature i for the prediction $f(x)$ (assuming feature independence) is just
$${\varphi }_i={\beta }_i\cdot (x_i-E[x_i])$$
${\beta }_i$ 是模型学出来的权重, $x_i$ 是当前sample 中特征i的取值, $E[x_i]$ 是数据集中特征i的取值的期望.
可以看出来针对LR, shap 几乎啥都没做.

DeepExplainer

用于 Deep NN 的模型可解释, 原理是 shap 与 ,官网例子见 参考[3].
注意用到了 keras, tensoflow, shap 三个库, 很容易有版本不兼容问题, 导致示例代码不能顺利运行.
类shap.explainers._deep.Deep继承了shap.explainers._explainer.Explainer, 根据model框架不同, 具体干活的又分为 TFDeep 与 PyTorchDeep.

PartitionExplainer

见我的另一篇文章, 参考[10].

参考

1. github, shap
2. 代码例子, Sentiment Analysis with Logistic Regression
3. 代码例子, deepexplainer
4. shap 三方包的论文, a-unified-approach-to-interpreting-model-predictions.pdf
5. 划分树, 联盟划分, owen 计算等, Mutual information-based group explainers with coalition structure
   for machine learning model explanations
6. 知乎文章, 关于Shapley Value（夏普利值）的公式
7. shap库 文档, Brute Force Kernel SHAP
8. Erik Štrumbelj and Igor Kononenko. “Explaining prediction models and individual predictions with feature contributions”. In: Knowledge andd information systems 41.3 (2014), pp. 647–665.
9. Lime paper, “Why Should I Trust You?”: Explaining the Predictions of Any Classifier
10. my blog, shap 中 PartitionExplainer 原理解读与官方demo的调试笔记