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


SHAP的理解与应用

SHAP有两个核心，分别是shap values和shap interaction values，在官方的应用中，主要有三种，分别是force plot、summary plot和dependence plot，这三种应用都是对shap values和shap interaction values进行处理后得到的。下面会介绍SHAP的官方示例，以及我个人对SHAP的理解和应用。

1. SHAP官方示例

首先简单介绍下shap values和shap interaction values到底是什么东西。

根据官方的示例运行下面代码读取数据并得到shap values和shap interaction values：

这里的y是个人年收入超过5万美元的可能性，是一个二分类模型。

1.1 shap values介绍

SHAP将模型的预测值解释为每个输入特征的归因值之和，其中归因值就是shap values，根据公式

但是lightgbm实际得到的shap values与模型预测值的关系并不是​

# 如果explainer.expected_value得到两个值，则应该为explainer.expected_value[1]

print('解释模型的常数:', explainer.expected_value)

print('训练样本预测值的log odds ratio的均值:', np.log(model.predict_proba(X_train)[:, 1]/ (1 - model.predict_proba(X_train)[:, 1])).mean())

print('常数与归因值之和:', explainer.expected_value + shap_values[0].sum())

print('预测值:', model.predict_proba(X.iloc[0:1])[:, 1][0])

print('预测值的log odds ratio:', np.log(model.predict_proba(X.iloc[0:1])[:, 1][0] / (1 - model.predict_proba(X.iloc[0:1])[:, 1][0])))

解释模型的常数: -2.359295492472404

训练样本预测值的log odds ratio的均值: -2.3592954924724046

常数与归因值之和: -4.342723316250527

预测值: 0.012834215330737286

预测值的log odds ratio: -4.342723316250522

根据上面的输出，可以看到：

1. 解释模型的常数 = 训练样本模型预测值的log odds ratio的均值
2. 常数与归因值之和 = 模型预测值的log odds ratio

1.2 shap interaction values介绍

shap interaction values则是特征俩俩之间的交互归因值，用于捕捉成对的相互作用效果，与shap values的关系为

可以与​

由于shap interaction values得到的是相互作用的交互归因值，假设有N个样本M个特征时，shap values的维度是N×M，而shap interaction values的维度是N×M×M，也就是说一个样本的一个特征，shap valus由一个归因值​

同样用第一个样本进行验证：

1.3 force plot

force plot是针对单个样本预测的解释，它可以将shap values可视化为force，每个特征值都是一个增加或减少预测的force，预测从基线开始，基线是解释模型的常数，每个归因值是一个箭头，增加（正值）或减少（负值）预测。

同样是第一个样本，其force plot如下：

红色的为正贡献，蓝色为负贡献，对于第一个样本的模型预测值0.0128，由上图可以解释为特征Education-Num=13的正贡献最大，其次是Age=39，但是Capital Gain=2174的负贡献很大，其次是Relationship=0，意思是没有家庭，在所有特征的综合影响下，该样本个人年收入超过5万美元的可能性只有0.0128。

这个图是直接由shap values绘成的，可以比较第一个样本的shap values具体数值：

sample_0_shap = pd.DataFrame(X.iloc[0,:])

sample_0_shap.rename(columns={0: 'feature_value'}, inplace=True)

sample_0_shap['shap_value'] = shap_values[0]

sample_0_shap.sort_values('shap_value', ascending=False)

所以就算不绘制force plot，直接获取样本的shap values，就可以知道每个特征值是如何贡献得到模型预测值的。

1.4 summary plot

summary plot是针对全部样本预测的解释，有两种图，一种是取每个特征的shap values的平均绝对值来获得标准条形图，这个其实就是全局重要度，另一种是通过散点简单绘制每个样本的每个特征的shap values，通过颜色可以看到特征值大小与预测影响之间的关系，同时展示其特征值分布。两种图分别如下：

shap.summary_plot(shap_values, X, plot_type="bar")

shap.summary_plot(shap_values, X)

两个图都可以看到Relationship全局重要度是最高的，其次是Age。第一个图可以看到各个特征重要度的相对关系，虽然Capital Gain是第三，但是重要度只有Relationship的60%，而第二个图由颜色深浅则可以看到Relationship和Age都是值越大，个人年收入超过5万美元的可能性越大。

其实如果要查看特征值大小与预测影响之间的关系的话，第二种图还是不够清楚，所以这里主要讲第一种图，第一种图其实就是对shap values按照特征维度聚合计算平均绝对值，也就是

feature_importance = pd.DataFrame()

feature_importance['feature'] = X.columns

feature_importance['importance'] = np.abs(shap_values).mean(0)

feature_importance.sort_values('importance', ascending=False)

可以看到和条形图的关系一模一样，所以获取到shap values后，直接对shap values按照特征维度聚合计算平均绝对值，就可以得到summary plot的结果。

1.5 dependence plot

如果要看特征值大小与预测影响之间的关系使用dependence plot更合适，dependence plot清楚地展示了单个特征是如何影响模型的预测结果的，dependence plot同样有多种使用方式，一种是查看某个特征是如何影响到模型预测结果的，另一种是一个特征是如何和另一个特征交互影响到模型预测结果的。

首先先看某个特征是如何影响到模型预测结果的，这里以Age为例子，运行下面代码：

shap.dependence_plot('Age', shap_values, X, interaction_index=None)

可以看到Age越大，个人年收入超过5万美元的可能性越大，但是在65岁后出现波动，80岁后可能性下降。

这个图的绘制也容易理解，是由Age的特征值与Age的shap值进行对应，分别为散点图的x与y即可绘制得到，运行下面代码，可以得到与上面一摸一样的图：

plt.figure(figsize=(7.5, 5))

plt.scatter(X['Age'], shap_values[:, 0], s=10, alpha=1)

接下来看一个特征是如何和另一个特征交互影响到模型预测结果的，这里以Age和Education-Num为例子，运行下面代码：

shap.dependence_plot('Age', shap_values, X, interaction_index='Education-Num')

虽然散点还是由Age绘制，但是Age的预测其实是有其他特征相互作用的，散点图垂直分散就是由相互作用效应驱动，所以用Education-Num进行着色以突出显示可能的相互作用，由上图可以看到，在50岁左右的时候，Education-Num更大，个人年收入超过5万美元的可能性越大。

我个人使用SHAP主要还是看中了单个特征与模型预测的归因关系，但是由上图可以看到用shap values绘制归因关系是有其他特征的相互作用使图垂直分散的，那怎么才能消除这种相互作用？这时候就可以使用shap interaction values，前面介绍了​

shap.dependence_plot(('Age', 'Age'), shap_interaction_values, X, interaction_index=None)

可以看到新的归因关系图因为去除了其他特征的相互作用，垂直分散就没那么厉害了，而这个图的绘制也很简单， 因为shap interaction values的维度是N×M×M，我们只需用[:, feature_index, feature_index]切片就可以得到​，用下面代码可以得到上图的结果：

plt.figure(figsize=(7.5, 5))

plt.scatter(X['Age'], shap_interaction_values[:, 0, 0], s=10, alpha=1)

这里还有一个使用经验，如果只是用shap interaction values绘制归因关系的话，一开始explainer.shap_interaction_values(X)返回shap interaction values时不需要使用全量的X，对X进行抽样即可，毕竟shap interaction values的计算时间还是较长的，而且最后绘制归因关系的散点图时也不需要全量数据，用部分数据就可以代表归因关系了。

现在force plot、summary plot和dependence plot是如何绘制的就介绍完了，我个人使用喜欢直接对shap_values和shap_interaction_values进行处理，自己去绘制全局重要度和归因关系，一个原因是我喜欢把图封装到html中进行展示，比如下图，在html可以选择特征去展示：

另一个原因是并不是所有模型都支持类别型变量的，需要将特征进行one-hot编码，这时候就需要自己去处理shap_values和shap_interaction_values，除非将类别下变量进行label编码当作数值型去处理，然后再绘制归因关系时特征值使用原始值。

2 个人使用示例

前面说了并不是所有模型都支持类别型特征，需要将特征进行one-hot编码，而且在我的使用中，就算是支持类别特征的LightGBM，SHAP在获取shap interaction values时也会报错。

如果有注意数据集，X_display其实就是X的原始值，里面Relationship、Workclass等值其实是类别特征，我们用X_display去建模，在explainer.shap_interaction_values(X_display)时会报错。

上面代码在最后一部会返回错误ValueError: could not convert string to float: ’ State-gov’，原因是TreeEnsemble类把input_dtype属性赋了np.float64，并且对X的特征做了个X.astype(self.model.input_dtype)转换，这里引发了错误，而且在后面使用values属性时也不行，因为有一个把trees属性赋None的操作，这也是会在后面引发错误的。

这个错误不知道是不是SHAP的版本问题，但是为了让对类别特征one-hot编码的数据也使用SHAP，有必要对shap_values和shap interaction values进行另外的处理。

2.1 类别特征的处理

这里介绍下如果对类别特征one-hot编码了，如何去还原shap values和shap interaction values。

举一个例子，比如特征Sex，分别有值Male和Female，对该特征进行one-hot后有特征Sex_Male和Sex_Female，用one-hot后的数据建模并获取shap values，此时shap values在Sex特征处应该有两列归因值，分别是

shap interaction values同理，由于我们使用shap interaction values绘制归因图是为了去除其他特征的相互作用，我们只要对

有了上面的理论基础，就可以进行shap values和shap interaction values的处理了，这里对方法进行了封装：

def onehot_pipeline(model, X_train, y_train, char_cols=None, num_fillna=None, char_fillna=None):

    传入带有参数的模型,封装成类别特征one-hot的pipline

    ————————————————————————————————————

        X_train:训练集的特征,pd.DataFrame格式

        char_cols:类别特征的列表,如不传入自动根据数据类型获取

        num_fillna:数值特征的缺失填充值,可支持不填充

        char_fillna:类别特征的缺失填充值,可支持不传入,但模型会自动填充null用于one-hot

        pipeline:封装好mapper和model的pipeline,并训练完成

    from sklearn_pandas import DataFrameMapper

    from sklearn.preprocessing import OneHotEncoder

    from sklearn2pmml.decoration import ContinuousDomain, CategoricalDomain

    from sklearn.pipeline import Pipeline

        col_types = X_train.dtypes

        char_cols = list(col_types[col_types.apply(lambda x: 'int' not in str(x) and 'float' not in str(x))].index)

    num_cols = list(set(X_train.columns) - set(char_cols))

    if not isinstance(char_fillna, str):

        char_fillna = 'null'

    mapper = DataFrameMapper(

        [(num_cols, ContinuousDomain(missing_value_replacement=num_fillna, with_data=False))] +

        [([char_col], [CategoricalDomain(missing_value_replacement=char_fillna, invalid_value_treatment='as_is'),

                       OneHotEncoder(handle_unknown='ignore')])

         for char_col in char_cols]

    pipeline = Pipeline(steps=[('mapper', mapper), ('model', model)])

    pipeline.fit(X_train, y_train)

def pipeline_shap(pipeline, X_train, y_train, interaction=False, sample=None):

    获取由onehot_pipeline返回的pipeline的shap值

    ————————————————————————————————————

        pipeline:onehot_pipeline的返回对象

        X_train:训练集的特征,pd.DataFrame格式

        interaction:是否返回shap interaction values

        sample:抽样数int或抽样比例float,不传入则不抽样

        feature_values:如传入sample则是抽样后的X_train,否则为X_train

        shap_values:pd.DataFrame格式shap values,如interaction传入True,则为shap interaction values

    if isinstance(sample, int):

        feature_values = X_train.sample(n=sample)

    elif isinstance(sample, float):

        feature_values = X_train.sample(frac=sample)

        feature_values = X_train

    mapper = pipeline.steps[0][1]

    model = pipeline._final_estimator

    sort_cols, onehot_cols = [], []

    for i in mapper.features:

        sort_cols += i[0]

        if 'OneHot' in str(i[1]):

            onehot_cols += i[0]

    feature_values = feature_values[sort_cols]

    X_train_mapper = mapper.transform(X_train)

    feature_values_mapper = mapper.transform(feature_values)

    model.fit(X_train_mapper, y_train)

    shap_values = pd.DataFrame(index=feature_values.index, columns=feature_values.columns)

    explainer = shap.TreeExplainer(model)

        mapper_shap_values = explainer.shap_interaction_values(feature_values_mapper)

        for col in sort_cols:

            if col in onehot_cols:

                col_index_span = len(X_train[col].unique())

                shap_values[col] = mapper_shap_values[

                    :, col_index: col_index + col_index_span, col_index: col_index + col_index_span

                ].sum(2).sum(1)

                col_index += col_index_span

                shap_values[col] = mapper_shap_values[:, col_index, col_index]

                col_index += 1

        mapper_shap_values = explainer.shap_values(feature_values_mapper)

        if len(mapper_shap_values) == 2:

            mapper_shap_values = mapper_shap_values[1]

        for col in sort_cols:

            if col in onehot_cols:

                col_index_span = len(X_train[col].unique())

                shap_values[col] = mapper_shap_values[

                    :, col_index: col_index + col_index_span

                col_index += col_index_span

                shap_values[col] = mapper_shap_values[:, col_index]

                col_index += 1

    return feature_values, shap_values

运行下面代码，查看特征重要度：

可以看到和前面将类别特征当作数值型特征去处理的特征重要度差异较大，one-hot的特征重要度是Marital Status最重要，而之前重要度最高的Relationship在这里只有第8，主要原因是之前将Relationship当作数值特征去处理了，由于shap values还是支持类别特征的，我们可以看下LightGBM直接处理类别特征的结果：

feature_importance = pd.DataFrame()

feature_importance['feature'] = X_display.columns

feature_importance['importance'] = np.abs(explainer.shap_values(X_display)[1]).mean(0)

feature_importance.sort_values('importance', ascending=False)

这个得到的特征重要度则和one-hot处理的比较相近，都是Marital Status最重要，但还有些区别，主要还是one-hot处理和LightGBM直接处理类别特征的差异。

最后看下归因关系：

plt.figure(figsize=(7.5, 5))

plt.scatter(feature_values['Age'], shap_interaction_values['Age'], s=10, alpha=1)

欢迎关注《python金融风控评分卡模型和数据分析(加强版)》，学习更多相关内容。