特征重要性之shap value

20221118补充

关于分类问题的base，回归是均值，分类是：base_values + sum(shap_values[0].values) = ln(p/1-p)

下面直接看代码：

# -*- coding: utf-8 -*-
"""
Created on Fri Nov 18 10:07:54 2022

@author: xxx
"""

# 引入数据
from sklearn import datasets
import numpy as np

data = pd.read_csv('F:/give_me_some_credit/cs-training.csv')
data.pop('id')
y = data.SeriousDlqin2yrs
x = data[['RevolvingUtilizationOfUnsecuredLines', 'age',
       'NumberOfTime30-59DaysPastDueNotWorse', 'DebtRatio', 'MonthlyIncome',
       'NumberOfOpenCreditLinesAndLoans', 'NumberOfTimes90DaysLate',
       'NumberRealEstateLoansOrLines', 'NumberOfTime60-89DaysPastDueNotWorse',
       'NumberOfDependents']]



import xgboost as xgb
from sklearn.model_selection import train_test_split
params={
        'eta':0.2,                        #特征权重，取值范围0~1，通常最后设置eta为0.01~0.2
        'max_depth':3,                    #树的深度，通常取值3-10，过大容易过拟合，过小欠拟合  230
        'min_child_weight':230,             #最小样本的权重，调大参数可以繁殖过拟合
        'gamma':0.4,                      #控制是否后剪枝，越大越保守，一般0.1、 0.2的样子
        'subsample':0.8,                  #随机取样比例
        'colsample_bytree':0.8 ,          #默认为1，取值0~1，对特征随机采集比例
        'lambda':0.8,
        'alpha':0.6,
        'n_estimators':500,
        'booster':'gbtree',               #迭代树
        'objective':'binary:logistic',    #逻辑回归，输出为概率
        'nthread':6,                      #设置最大的进程量，若不设置则会使用全部资源
        'scale_pos_weight':3,             #默认为0,1可以处理类别不平衡

        'seed':1234,                      #随机树种子
        'silent':1,                       #0表示输出结果
        'eval_metric':'auc'               #评分指标
        }
model = xgb.XGBClassifier(**params) #使用时主要区别在这里，其实接口形式的和其他的模型用法基本一样
model.fit(x, y)  
y_predicted = model.predict_proba(x)
y_predicted.mean()
p = y_predicted[:,1]
  

方法一：
import shap
explainer = shap.TreeExplainer(model)
shap_values = explainer.shap_values(x)  # 传入特征矩阵X，计算SHAP值
y_base = explainer.expected_value


方法二：
shap_values = model.get_booster().predict(xgb.DMatrix(x), pred_contribs=True) #这一步将每一行数据转成shap构成的形式，最后一列是shap的base值，其余的是每个变量的shap的值。
#shap_df = pd.DataFrame(np.abs(shap_values[:, :-1]), columns=cols) #这一步是先取绝对值然后再讲array转成pf
#shap_importance = shap_df.mean().sort_values(ascending=False).reset_index()  #最后计算每个特征的shap值的均值
#shap_importance.columns = ['Feature', 'Shap_Importance']

根据公式验证：
base_values + sum(shap_values[0].values) = ln(p/1-p)

 

 

基于模型刷选特征方法有：排列重要性、shap value、null importance

下面来说一下 shap value

一、shap value的原理

在SHAP被广泛使用之前，我们通常用feature importance或者partial dependence plot来解释xgboost。Feature importance可以直观地反映出特征的重要性，看出哪些特征对最终的模型影响较大。但是无法判断特征与最终预测结果的关系是如何的，是正相关、负相关还是其他更复杂的相关性？因此就引起来SHAP。

SHAP的名称来源于SHapley Additive exPlanation。Shapley value起源于合作博弈论。比如说甲乙丙丁四个工人一起打工，甲和乙完成了价值100元的工件，甲、乙、丙完成了价值120元的工件，乙、丙、丁完成了价值150元的工件，甲、丁完成了价值90元的工件，那么该如何公平、合理地分配这四个人的工钱呢？Shapley提出了一个合理的计算方法（有兴趣地可以查看原论文），我们称每个参与者分配到的数额为Shapley value。

SHAP是由Shapley value启发的可加性解释模型。对于每个预测样本，模型都产生一个预测值，SHAP value就是该样本中每个特征所分配到的数值。 假设第i个样本为xi，第i个样本的第j个特征为xi,j，模型对第i个样本的预测值为yi，整个模型的基线（通常是所有样本的目标变量的均值）为ybase，那么SHAP value服从以下等式。

其中f(xi,1)为xi,j的SHAP值。直观上看，f(xi,1)就是第i个样本中第1个特征对最终预测值yi的贡献值，当f(xi,1)>0，说明该特征提升了预测值，也正向作用；反之，说明该特征使得预测值降低，有反作用。SHAP value最大的优势是SHAP能对于反映出每一个样本中的特征的影响力，而且还表现出影响的正负性

二、shap的实现

前期简单建模，使用波士顿房价数据

# -*- coding: utf-8 -*-
"""
Created on Sun Sep 26 15:51:26 2021

@author: chenguimei
"""



# 加载模块
import xgboost as xgb
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt; plt.style.use('seaborn')

#波士顿房价数据集
from  sklearn.datasets import load_boston
boston=load_boston()
boston.data
boston.target
boston.feature_names
boston_df=pd.DataFrame(boston.data,columns=boston.feature_names)
boston_df['target'] = boston.target 

cols = [i for i in boston_df.columns[:-1]]

# 训练xgboost回归模型
model = xgb.XGBRegressor(max_depth=4, learning_rate=0.05, n_estimators=150)
model.fit(boston_df[cols], boston_df['target'].values)


# 获取feature importance
plt.figure(figsize=(15, 5))
plt.bar(range(len(cols)), model.feature_importances_)
plt.xticks(range(len(cols)), cols, rotation=-45, fontsize=14)
plt.title('Feature importance', fontsize=14)
plt.show()

 将每一行每个特征的shap value值放在一个df中，也就是将上面的公式实例出来

import shap
# model是在第1节中训练的模型
explainer = shap.TreeExplainer(model)
shap_values = explainer.shap_values(boston_df[cols])

 

（1）单个样本（一行数据）的SHAP值

随机检查其中一个样本的每个特征对预测值的影响

#第一列是特征名称，第二列是特征的数值，第三列是各个特征在该样本中对应的SHAP值。
# 比如我们挑选数据集中的第30位
j = 30
player_explainer = pd.DataFrame()
player_explainer['feature'] = cols
player_explainer['feature_value'] = boston_df[cols].iloc[j].values
player_explainer['shap_value'] = shap_values[j]
player_explainer['base'] = model.predict(boston_df[cols]).mean() #就是预测的分数的均值
player_explainer['sum'] = player_explainer['shap_value'].sum() #特征的shap和
player_explainer['base+sum'] = player_explainer['base']+player_explainer['sum']
player_explainer

shap还提供极其强大的数据可视化功能

shap.initjs()
shap.force_plot(explainer.expected_value, shap_values[j], boston_df[cols].iloc[j])

 蓝色表示该特征的贡献是负数，红色则表示该特征的贡献是正数，黑色加粗的是该样本最后的预测值。最长的蓝色条是LSTAT，该特征的含义是人口密度，密度低，房价自然也低。

 

（2）全部特征的分析

除了能对单个样本的SHAP值进行可视化之外，还能对特征进行整体的可视化

shap.summary_plot(shap_values, boston_df[cols])

 图中每一行代表一个特征，横坐标为SHAP值。一个点代表一个样本，颜色越红说明特征本身数值越大，颜色越蓝说明特征本身数值越小，可以看出LSTAT越大，房价越小，和房价成反比关系

 

也可以把一个特征对目标变量影响程度的绝对值的均值作为这个特征的重要性

shap.summary_plot(shap_values, boston_df[cols], plot_type="bar")

 上面这个图怎么得来的，我就演示一下LSTAT变量的值

#LSTAT位于最后一个，因此我们只需要提取最后一列
pd.DataFrame(shap_values).iloc[:,-1].apply(lambda x:abs(x)).mean()  #输出 3.7574333926117998

 

 (3)部分依赖图Partial Dependence Plot

就是shap值和原值的散点图，可以看出趋势，是单调还是U型，等等

shap.dependence_plot('LSTAT', shap_values, boston_df[cols], interaction_index=None, show=False)

 

 （4）对多个变量的交互进行分析

我们也可以多个变量的交互作用进行分析。一种方式是采用summary_plot描绘出散点图

shap interaction values则是特征俩俩之间的交互归因值，用于捕捉成对的相互作用效果，由于shap interaction values得到的是相互作用的交互归因值，假设有N个样本M个特征时，shap values的维度是N×M，而shap interaction values的维度是N×M×M，也就是说一个样本的一个特征shap valus由一个归因值对应，而shap interaction values由一系列交互归因值对应，并且交互归因值的和等于归因值。

我是这样理解：

shap valus = base+x1_shap+x2_shap+...+xn_shap

shap interaction values = base + x1_shap+x2_shap+...+xn_shap+x2*x3*x4...xn_shap_interaction_values+....

shap_interaction_values = shap.TreeExplainer(model).shap_interaction_values(boston_df[cols])
shap.summary_plot(shap_interaction_values, boston_df[cols], max_display=4)

 看这个图不是很明白，我们看具体数据，就很清晰明了

print('第一个特征的shap values:', shap_values[0][0])
print('第一个特征的shap interaction values:', shap_interaction_values[0][0])
print('第一个特征的shap interaction values和:', shap_interaction_values[0][0].sum())

'''
第一个特征的shap values: -0.38626033
第一个特征的shap interaction values: [ 8.65335166e-02  5.80601394e-04  5.75962476e-03 -1.57106481e-03
 -7.26313889e-02 -6.15637302e-02 -7.20938668e-02 -6.41713059e-03
 -4.72193956e-03 -3.00363600e-02 -1.32852495e-02  1.28373504e-05
 -2.16826200e-01]
第一个特征的shap interaction values和: -0.38626036
'''

 

 我们也可以用dependence_plot描绘两个变量交互下变量对目标值的影响，这样看更加明白

for i in ['LSTAT','RM','DIS','PTRATIO']:
    shap.dependence_plot('LSTAT', shap_values, boston_df[cols], interaction_index=i, show=False)

 

 

 

 

 

 从上面的图片我们可以看出，'RM','DIS','PTRATIO'，这个颜色代表着值得大小，越红代表越大，越蓝代表越小，比如说LSTAT和RM，可以看出，在LSTAT越小的时候，RM的颜色基本是红色的，说明此时RM的值是比较大的，而当LSTAT越大是，RM的颜色基本是蓝色，说明RM的值是越来越小，因此可以推断出LSTAT和RM的交互系数应该是一份负值。

for i in ['RM','LSTAT']:
    shap.dependence_plot('RM', shap_values, boston_df[cols], interaction_index=i, show=False)

我们看具体的LSTAT和RM的交互系数分布

 

2022年6月8号补充

xgb自带的计算shap值

# 加载模块
import xgboost as xgb
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt; plt.style.use('seaborn')

#波士顿房价数据集
from  sklearn.datasets import load_boston
boston=load_boston()
boston.data
boston.target
boston.feature_names
boston_df=pd.DataFrame(boston.data,columns=boston.feature_names)
boston_df['target'] = boston.target 


cols = [i for i in boston_df.columns[:-1]]

# 训练xgboost回归模型
model = xgb.XGBRegressor(max_depth=4, learning_rate=0.05, n_estimators=150)
model.fit(boston_df[cols], boston_df['target'].values)



shap_values = model.get_booster().predict(xgb.DMatrix(boston_df[cols]), pred_contribs=True) #这一步将每一行数据转成shap构成的形式，最后一列是shap的base值，其余的是每个变量的shap的值。
shap_df = pd.DataFrame(np.abs(shap_values[:, :-1]), columns=cols) #这一步是先取绝对值然后再讲array转成pf
shap_importance = shap_df.mean().sort_values(ascending=False).reset_index()  #最后计算每个特征的shap值的均值
shap_importance.columns = ['Feature', 'Shap_Importance']

我们来看看细节，重点看上面代码的注释，然后再联系到开头说shap值的定义时，一切就不言而喻了。