本文介绍了模型可解释性的重要性，特别是金融和风险管理领域。文章探讨了基于统计学的特征选择，如WOE和IV值在分类问题中的应用，以及回归问题中特征的相关性分析。此外，文章详细阐述了基于博弈论的SHAP值计算方法，用于评估模型中每个特征的边际贡献，强调了SHAP在模型解释中的作用。最后，通过示例展示了如何使用SHAP进行特征重要性可视化，包括单样本预测值解释和特征重要度输出。
前言
随着人工智能技术的发展，我们可以使用更加复杂和高级的模型来解决复杂的业务问题。然而，这些复杂的模型也带来了可解释性的问题，特别是在金融和风险管理领域。因此，我们需要确保我们的模型结果是可解释和可追溯的。这意味着我们需要开发和应用透明、可理解的模型，确保模型结果的准确性和有效性。
在模型可解释方面的研究主要集中在以下几点：
- 可解释性标准的制定（忽略）：提出了一些新的度量标准，如透明度指数和解释性指数等，以便更好地比较和选择不同的模型，便于选出精度尚可且具有可解释性能力的论文
- 自适应解释性模型的研究：根据不同的数据和应用场景，自适应地调整模型参数和结构，以提高模型的可解释性和实用性。（近期论文：AutoLIME: Automated Local Interpretable Model Explanations for Predictive Models）
- 因果推断的研究（各种因果图实现等）：有因才有果，利用贝叶斯网络或uplift等模型，研究因子之间的联系，从而找到影响模型结果的几个根因。
- 模型特征重要性方面的研究（shap等）：找到对模型结果影响重要的特征，组织特征表现，对模型结果进行可视化的解释。（一般业务上会优先考虑该方法）
模型可解释性是一个很大很泛的命题，本文主要介绍特征重要性在工程上的使用以及一些理论基础。
特征重要性
在建模过程中，特征重要性的思考分为：
- 基于统计学的特征选择
- 基于模型的特征选择
- 基于博弈论的特征选择
基于统计学的特征选择
在搭建模型中，样本特征的“独立性”越好，不同目标间样本的特征差异性越大，其特征往往能够对模型产生比较正面的效果。
-
分类问题的特征选择
建模目标为分类问题，特别是二分类问题时，WOE和IV值的计算常出现在特征选择中，WOE和IV值的计算主要是为了衡量：在特征分箱中，评估正样本与负样本间的差异，从而评估自变量对分类结果的影响。-
woe值计算
WOE：着重于对比不同自变量值之间的影响大小，可以将WOE值降序后进行特征选择。woe的计算公式如下：
w o e = l n ( g o o d n u m s n u m s ) − l n ( b a d n u m s n u m s ) = l n ( 正样本占比 负样本占比 ) woe = ln(\frac{good_{nums}}{nums})-ln(\frac{bad_{nums}}{nums})=ln(\frac{正样本占比}{负样本占比}) woe=ln(numsgoodnums)−ln(numsbadnums)=ln(负样本占比正样本占比)
woe值越大，表示自变量对目标变量的影响越大。 -
iv值计算
IV方法更加侧重于特征对标签的预测能力评估，IV是一种用于评估自变量的预测能力的方法。其计算公式为:
I V = ∑ i = 0 p ( g o o d n u m s i n u m s i ) ∗ w o e i IV = \sum_{i=0}^p(\frac{good_{nums_{i}}}{nums_{i}})*woe_{i} IV=i=0∑p(numsigoodnumsi)∗woei
p：特征的分箱数。IV值越高，表示自变量对目标变量的预测能力越强，对分类结果的拟合效果越好。通常，IV值小于0.02的特征可以认为是无用的，0.02-0.1的特征为较弱的特征，0.1-0.3的特征为中等强度特征，于大0.3的特征为强特征。在特征选择中，我们通常会选择IV值较高的特征作为模型的输入特征，以提高模型的预测能力和泛化能力。
-
-
回归问题的特征选择
对于回归模型，在进行特征选择时，主要是利用自变量和目标变量间的相关性、以及自变量间的相关性进行：- 冗余特征的剔除
- 相关性弱的特征剔除
相关系数的计算方法多种多样，常用的计算方法有皮尔逊相关系数、斯皮尔曼等级相关系数以及余弦相似度等计算方法，其中皮尔逊相关系数的计算方法如下：
r x y = s u m ( ( x − m e a n ( x ) ) ∗ ( y − m e a n ( y ) ) ) / s q r t ( s u m ( ( x − m e a n ( x ) ) 2 ) ∗ s u m ( ( y − m e a n ( y ) ) 2 ) ) r_{xy} = sum((x - mean(x)) * (y - mean(y))) / sqrt(sum((x - mean(x))^2) * sum((y - mean(y))^2)) rxy=sum((x−mean(x))∗(y−mean(y)))/sqrt(sum((x−mean(x))2)∗sum((y−mean(y))2))
一般而言，皮尔逊系数的取值范围在-1到1之间，取值为1表示两个变量呈完全正相关关系，取值为-1表示两个变量呈完全负相关关系，取值为0表示两个变量之间没有线性相关关系。一般来说，取值在0.8-1之间的皮尔逊系数表示两个变量之间呈强正相关关系，取值在0.6-0.8之间表示两个变量之间呈中等正相关关系，取值在0.4-0.6之间表示两个变量之间呈弱正相关关系，取值在***-0.4-0.4之间表示两个变量之间没有线性相关关系***，取值在-0.6–0.4之间表示两个变量之间呈弱负相关关系，取值在-0.8–0.6之间表示两个量变之间呈中等负相关关系，取值在-1–0.8之间表示两个变量之间呈强负相关关系。因此，当皮尔逊系数的取值越接近于1或-1时，说明两个变量之间的相关性越强。 -
无监督类型的特征选择
后期详细介绍，某种意义上来看无监督学习（IF、Kmeans以及深度聚类等方法）对特征的选择更加敏感。简单来说：特征需要让样本具有差异性。思想：- 模型方法：PCA以及SVD等奇异值分析的降维方法。这类方法简单，且具有不错的效果，但是模型难以保留其可解释性。
- 统计学的方法：基于频次等信息的计算
基于模型的特征选择（略）
一本基于LR以及树模型等进行特征选择，方法：基于模型训练后利用Feature_important进行特征选择。
案例可参考：https://blog.csdn.net/qq_42363032/article/details/113695733
基于博弈论的特征选择（建模后）
shap作为“一种通用的解释性AI框架”在建模后的模型可解释性和特征筛选上发挥着越来越重要的作用。
本节将从以下几个方面来介绍Shap：
shap产生的基本思想
SHAP的思想源自于博弈论，旨在通过计算每个特征对预测结果的边际贡献来刻画其对于模型的重要性。
具体来说，SHAP使用Shapley值的概念，即在一个合作博弈中，对于每一个参与者，其对于合作所产生的价值的贡献，来计算每个特征的边际效应。例如：
假设有一个公司，有三个员工甲、乙和丙。当甲、乙以及丙一起工作时，为公司带来了120元的效益；乙和丙一起工作时，为公司带来了80元的效益；甲和丙一起工作为公司带来了60的效益；甲乙一起工作为公司带来50的效应；甲单独工作时为公司带来20的效应；乙单独工作时为公司带来60的效应；丙单独工作时为公司带来60的效应。为了计算甲的边际效应，我们可以使用SHAP的方法进行计算：
首先，我们计算空联盟（甲乙丙均不工作）时，甲单独工作带来的收益为20元。
接着，我们计算乙丙工作时，甲的边际贡献为120 - 80 = 40元。
然后，我们计算乙工作时，甲的边际贡献为50 - 60 = -10元。
最后，我们计算丙工作时，甲的边际贡献为60 - 60 = 0元。
通过加权计算这些边际贡献，我们可以得到甲的Shapley值为（40 + 20 - 10 + 0）/ 4 = 12.5元。同样的方式，我们也可以计算出乙和丙的Shapley值。
为什么不直接采用控制变量法来得出甲能够带来40元的效应呢？原因有两点（个人理解）：一是特征之间可能存在相互作用，控制变量法可能无法完全消除这些影响；二是概率和边缘概率的关系可能会影响结果的计算。因此，SHAP提供了一种更加有效、精确的计算特征重要性的方法。
shap的公式推导
参考一般线性模型中，计算特征的各个贡献度，本质是计算每个特征对应的权重：
f
(
x
)
=
β
0
+
β
1
∗
x
1
+
β
2
∗
x
2
+
.
.
.
+
β
p
∗
x
p
\\f(x)=\beta_0+\beta_1*x_1+\beta_2*x_2+...+\beta_p*x_p
f(x)=β0+β1∗x1+β2∗x2+...+βp∗xp
但是对于非线性模型，是很难之际计算出
β
\beta
β的值，所以在shapely值的计算中，其实采用了遍历的方法进行计算。实际在shapely的计算中，会定义特征
i
i
i 对预测结果的贡献如下：
ψ
i
=
β
i
∗
x
i
−
β
i
∗
E
(
X
i
)
\psi_i = \beta_i*x_i-\beta_i*E(X_i)
ψi=βi∗xi−βi∗E(Xi)
其中
E
(
β
i
∗
X
i
)
E(\beta_i*X_i)
E(βi∗Xi)认为是特征
i
i
i 对模型结果的平均影响估计，则用样本的特征
i
i
i的贡献减去平均贡献，即是该样本中特征
i
i
i 对模型的贡献程度。故该样本的所有特征贡献可得到计算：
∑
i
=
1
p
ψ
i
(
f
)
=
∑
i
=
1
p
(
β
i
∗
x
i
−
β
i
∗
E
(
X
i
)
)
=
f
(
x
)
−
E
(
f
(
x
)
)
\sum_{i=1}^p\psi_i(f)=\sum_{i=1}^p(\beta_i*x_i-\beta_i*E(X_i))=f(x)-E(f(x))
i=1∑pψi(f)=i=1∑p(βi∗xi−βi∗E(Xi))=f(x)−E(f(x))
即：数据点
x
x
x 的贡献之和就等于预测结果减去模型的平均预测值，这也就解释了模型对该样本的模型预测结果，为何相对模型预测的平均值偏大或偏小，利用
ψ
i
\psi_i
ψi 也可以知道哪个特征对该结果有正向作用，何种特征具有负向作用。
ψ
i
\psi_i
ψi的计算公式，根据上文的定义，可通俗的理解为：计算特征
i
i
i 在于不在时，面对不同的特征组合时对模型的预测结果的“增益” ，故其计算公式可以修改为：
ψ
i
(
v
a
l
)
=
∑
S
∈
x
i
,
x
2
,
.
.
.
,
x
p
/
x
j
∣
S
∣
!
(
p
−
∣
S
∣
−
1
)
!
p
!
(
v
a
l
(
S
∪
x
i
)
−
v
a
l
(
S
)
)
\psi_i(val) = \sum_{S\in{x_i,x_2,...,x_p}/{x_j}}\frac{|S|!(p-|S|-1)!}{p!}(val(S\cup x_{i})-val(S))
ψi(val)=S∈xi,x2,...,xp/xj∑p!∣S∣!(p−∣S∣−1)!(val(S∪xi)−val(S))
其中，
(
v
a
l
(
S
∪
x
i
)
−
v
a
l
(
S
)
)
(val(S\cup x_{i})-val(S))
(val(S∪xi)−val(S))表示在联盟中，加上特征
i
i
i ，和不加入特征
i
i
i 时，特征
i
i
i 产生的贡献，前面的
∣
S
∣
!
(
p
−
∣
S
∣
−
1
)
!
p
!
\frac{|S|!(p-|S|-1)!}{p!}
p!∣S∣!(p−∣S∣−1)!为权重，p为联盟中的特征数，S为选择的子集。
从上文公式可知，随着模型特征数和样本数的增加，shapley值在计算的过程中，其时间复杂度会大大增加，为了简化计算，在shapley值的计算过程中，引入了树模型的方法来加速计算（TreeExplainer等），使用Tree Shap的计算方法，理论上能够将模型的时间复杂度由
O
(
T
L
2
M
)
O(TL2^M)
O(TL2M) 降低到
O
(
T
L
D
2
)
O(TLD^2)
O(TLD2)，其中：T为树的数量，L是所有树的最大节点数量，M是特征的数量，D是树的最大深度。其Tree Shap的具体计算过程可以参考论文：https://arxiv.org/pdf/1706.06060.pdf【git源码】
shap在模型可解释上的实际应用
现阶段Python中的shap包已经能够实现各种各样的特征重要性可视化功能，这在工作中经常使用。本节将通过案例形式介绍部分常用的可视化功能。
具体来说，我们将展示如何使用shap来进行以下可视化：
- 单样本预测值解释性可视化
- shap特征重要度输出和可视化
- shap值的摘要图
- shap依赖图——单特征和预测值关系的可视化展示
此处用公开数据集Adult Census Income数据集为例，模型采用lightgbm模型，介绍特征可视化的使用。
模型搭建
coding is:
import pandas as pd
import numpy as np
import lightgbm as lgb
#from pycaret.regression import *
from sklearn.preprocessing import LabelEncoder
import joblib
class Estimator:
def __init__(self,path,category_name=None,label_name=None):
self.data = self._read_data(path=path)
self.category_name = category_name
self.label_name = label_name
def _read_data(self,path):
data = pd.read_csv(path)
print(data)
print(list(data.columns))
return data
# 数据类型转换
def transform_feature(self):
'''
:desc: 简单建模使用labelencoder进行编码
:return: dataframe
'''
for cate in self.category_name:
encoder = LabelEncoder()
self.data[cate] = encoder.fit_transform(self.data[cate])
self.data['label'] = self.data[self.label_name]
self.data.drop([self.label_name],axis=1,inplace=True)
# 补充缺失值---》直接补0(简化计算）
self.data = self.data.fillna(0)
return self.data
def _data_split(self,data,params=0.2):
from sklearn.model_selection import train_test_split
Y = data['label']
X = data.drop(['label'],axis=1)
X_train, X_test, y_train, y_test = train_test_split(X, Y, test_size=params, random_state=666)
from sklearn.preprocessing import StandardScaler
sc = StandardScaler()
X_train = sc.fit_transform(X_train)
X_test = sc.transform(X_test)
return X_train,y_train,X_test,y_test
def build_model(self,grid_search=None,is_save=True):
'''
:desc : 构建模型 默认使用默认的参数
:param grid_search:若不为空，使用grid_search遍历搜索参数
:return:
'''
self.transform_feature()
X_train, y_train,X_test, y_test = self._data_split(data=self.data)
if grid_search is None:
model = lgb.LGBMClassifier(
boosting='gbdt',
objective='binary',
metric='auc'
)
model.fit(X_train, y_train)
else:
model = None
if is_save == True:
joblib.dump(model, save_path)
#clf = joblib.load('dota_model.pkl')
# 评估模型
predcit = model.predict(X_test)
# 准召率计算
from sklearn import metrics
accuracy = metrics.accuracy_score(list(y_test),list(predcit))
recall = metrics.recall_score(list(y_test),list(predcit))
# auc计算
auc = metrics.roc_auc_score(list(y_test),list(predcit))
print("============= the model of accuracy is {} ===========".format(accuracy))
print("============= the model of recall is {} ===========".format(recall))
print("============= the model of auc is {} ============".format(auc))
if __name__ == '__main__':
path = data_path
category_name = [
'workclass','education','marital.status','occupation',
'relationship','race','sex','native.country','income'
]
label_name = 'income'
work = Estimator(path,category_name=category_name,label_name=label_name)
work.build_model()
训练结果：
shap解释模型
在使用shap进行模型解释的时候，需要注意验证shap对模型解释的有效性，并非所有的模型均可用shap来进行解释，比如在使用lstm的时候，若采用GPU加速，形成cudnnLSTM模型，shap是无法对该模型进行解释的。根本原因在于，使用shap进行解释时，无法满足shap的**可加性**原则。
f
(
x
)
=
f
b
a
s
e
+
∑
i
=
1
p
(
ψ
i
(
f
)
)
f(x)=f_{base}+\sum_{i=1}^p(\psi_i(f))
f(x)=fbase+i=1∑p(ψi(f))
其中，
f
(
x
)
f(x)
f(x)为模型预测的结果，
f
b
a
s
e
f_{base}
fbase为模型的期望输出结果，
ψ
i
(
f
)
\psi_i(f)
ψi(f)为每个特征为模型贡献的shap值。当shap对模型具有可解释性时，其可加性间的偏差不高于0.04.即：
a
b
s
(
f
(
x
)
−
(
f
b
a
s
e
+
∑
i
=
1
p
(
ψ
i
(
f
)
)
)
)
<
0.04
abs(f(x)-(f_{base}+\sum_{i=1}^p(\psi_i(f))))<0.04
abs(f(x)−(fbase+i=1∑p(ψi(f))))<0.04
coding 如下：
'''
desc:利用shap进行模型可解释性
'''
import pandas as pd
import numpy as np
import lightgbm as lgb
#from pycaret.regression import *
import random
import seaborn as sns
import shap
from tqdm import tqdm
from sklearn.preprocessing import LabelEncoder
import joblib
class Estimator:
def __init__(self,path,category_name=None,label_name=None):
self.data = self._read_data(path=path)
self.category_name = category_name
self.label_name = label_name
def _read_data(self,path):
data = pd.read_csv(path)
print(data)
print(list(data.columns))
return data
# 数据类型转换
def transform_feature(self):
'''
:desc: 简单建模使用labelencoder进行编码
:return: dataframe
'''
for cate in self.category_name:
encoder = LabelEncoder()
self.data[cate] = encoder.fit_transform(self.data[cate])
self.data['label'] = self.data[self.label_name]
self.data.drop([self.label_name],axis=1,inplace=True)
# 补充缺失值---》直接补0(简化计算）
self.data = self.data.fillna(0)
return self.data
def _data_split(self,data,params=0.2):
from sklearn.model_selection import train_test_split
Y = data['label']
X = data.drop(['label'],axis=1)
X_train, X_test, y_train, y_test = train_test_split(X, Y, test_size=params, random_state=666)
from sklearn.preprocessing import StandardScaler
sc = StandardScaler()
X_train = sc.fit_transform(X_train)
X_test = sc.transform(X_test)
return X_train,y_train,X_test,y_test
def build_model(self,grid_search=None,is_save=True):
'''
:desc : 构建模型 默认使用默认的参数
:param grid_search:若不为空，使用grid_search遍历搜索参数
:return:
'''
self.transform_feature()
X_train, y_train,X_test, y_test = self._data_split(data=self.data)
if grid_search is None:
model = lgb.LGBMClassifier(
boosting='gbdt',
objective='binary',
metric='auc'
)
model.fit(X_train, y_train)
else:
model = None
if is_save == True:
joblib.dump(model, r'F:\博客相关\shap相关\coding\model\lightgbm.pkl')
#clf = joblib.load('dota_model.pkl')
# 评估模型
predcit = model.predict(X_test)
# 准召率计算
from sklearn import metrics
accuracy = metrics.accuracy_score(list(y_test),list(predcit))
recall = metrics.recall_score(list(y_test),list(predcit))
# auc计算
auc = metrics.roc_auc_score(list(y_test),list(predcit))
print("============= the model of accuracy is {} ===========".format(accuracy))
print("============= the model of recall is {} ===========".format(recall))
print("============= the model of auc is {} ============".format(auc))
columns = list(self.data.columns)
print(columns)
columns.remove('label')
X_train = pd.DataFrame(X_train)
X_train.columns = columns
X_test = pd.DataFrame(X_test)
X_test.columns = columns
return X_train, y_train,X_test, y_test, model
def tree_explainer(self,model,data,sample_index = None,check_type = False,label='label'):
'''
:desc: 使用shap的treeExplanier进行解析
:param sample: 单条样本
:return: 单条样本基于shap的模型可解释性
'''
explainer = shap.TreeExplainer(model)
columns = list(data.columns)
try:
columns.remove(label)
except:
pass
#shap_values = explainer.shap_values(data)
shap_values = []
for i in tqdm(range(len(data))):
#print(pd.DataFrame(data.loc[i]).T)
shap_value = explainer.shap_values(pd.DataFrame(data.loc[i]).T)
shap_values.append(shap_value[0][0])
shap_values = pd.DataFrame(shap_values)
print(shap_values)
shap_values.columns = columns
# 若样本的feature太长，使用
# single sample plot
if(sample_index is None):
sample_index = random.randint(0,len(data))
# check shap 是否有效==>利用可加性验证模型
y_base = explainer.expected_value
if(check_type == True):
error = 0
for j in tqdm(range(len(data))):
y = self.model.predict(data.iloc[j])[label]
shap_sum = sum(shap_values[j])
error += abs(y-(y_base+shap_sum))
if(error/len(data)>0.05):
print("======================== 可加性不满足 shap解释不成立 ==========================")
return shap_values,shap_values
# 瀑布图（单样本）
shap.plots._waterfall.waterfall_legacy(y_base[0]
,np.array(shap_values.iloc[sample_index]),
shap_values.iloc[0])
# 力图(作用和瀑布图差不多）
#shap.plots.force(shap_values[sample_index])
# 决策图（n个样本shap值的可视化表达）
shap.decision_plot(y_base[0],np.array(shap_values[0:10]),features=columns)
# 特征重要性（使用模型的shap的绝对均值）
# 分组后计算
#cluster = shap.utils.hclust(data.iloc[:,:-1],data.iloc[:,-1])
#shap.plots.bar(shap_values,cluster=cluster) ====shap_values的计算方式需要改变
shap.summary_plot(np.array(shap_values),columns,plot_type="bar")
# 常见的蜂窝图
#shap.plots.beeswarm(shap_values)
shap.summary_plot(np.array(shap_values),columns,plot_size=0.3)
# 单变量散点图
# shap.plots.scatter(shap_values[:,"xxx"])
# 特征组合对模型的重要性
#shap_interaction_values = explainer.shap_interaction_values(data.iloc[:,-1])
#shap.summary_plot(shap_interaction_values,data.iloc[:,-1],max_display=10,plot_type='compact_dot')
feature_shap = pd.DataFrame(shap_values)
feature_shap.columns = columns
return feature_shap
if __name__ == '__main__':
sns.set(font='SimHei', font_scale=1.5)
path = r'F:\博客相关\shap相关\数据\Adult Census Income\adult.csv'
category_name = [
'workclass','education','marital.status','occupation',
'relationship','race','sex','native.country','income'
]
label_name = 'income'
work = Estimator(path,category_name=category_name,label_name=label_name)
X_train, y_train,X_test, y_test, model = work.build_model(is_save=False)
# shap特征重要性可视化
sns.set(font='Microsoft YaHei, SimHei') # 设置画图中的中文为黑体
feature_shap = work.tree_explainer(model,X_test)
# 加载模型
#model = joblib.load(r'F:\博客相关\shap相关\coding\model\lightgbm.pkl')
单样本预测值解释性可视化
-
单样本预测值解释性可视化
- 样本瀑布图
代码：
shap.plots._waterfall.waterfall_legacy(y_base[0] ,np.array(shap_values.iloc[sample_index])
其实用shap.plots.waterfall()展示更加简单。在样本较大的情况下，一般不太使用shap.plots.waterfall()，原因在于shap_values的两种求解方案。（在文末会提到）
图片解释：图片代表每个样本最终预测结果偏离期望值（均值）时，各个特征带来的边界效应。比如：age对该样本产生的影响最大，相对于期望值，其甙类-0.5的边际效应。
- 样本瀑布图
-
shap特征重要度输出和可视化
shap特征重要性的输出一般有两种形式，一种是将shap值进行输出，然后自己进行归一化处理后，将特征维度上的shap绝对值进行求和后进行输出。该方法的灵活性更高，比如在上一例子中，可以对不同性别的用户分别求解特征重要性。还有一种是直接用shap的自带函数：
shap.summary_plot(np.array(shap_values),columns,plot_type="bar") #cluster = shap.utils.hclust(data.iloc[:,:-1],data.iloc[:,-1]) #shap.plots.bar(shap_values,cluster=cluster)
其他关于shap的模型可视化可见：
https://zhuanlan.zhihu.com/p/441302127
备注
shap values的计算方式有两种：
shap_valuer = explainer.shap_values(X) # 为一个二维numpy
shap_values2 = explainer(X) # 除了shap值还有base值等
虽然大多画shap图的资料更倾向于第二种方案（因为画图更加简单），但是当data的样本量过大时，结合tqdm来进行shap values计算更容易看到shap值计算的进度，后续取出shap值后也更加方便分析。
扫一扫
1128
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
