模型解释工具：SHAP --- 论文实战
模型解释工具：SHAP --- 论文实战
模型可解释性指对模型内部机制的理解以及对模型结果的理解，即模型中的那些特征最重要；哪个特征造成效果的好与不好。在众多模型中，可解释性最强的当然是线性模型；但是由于线性模型在特征提取方面存在天然的弱势，以及过于依赖原始的特征工程，在稍微复杂一些的场景中，基本见不到线性模型的影子。现在效果比较好的几类模型，SVR、 neural network 模型可获得较好的效果，但他们同时也被称为“黑盒”模型，原因是这类模型不能描述特征重要性。所以，打开“黑盒”模型就成了一个值得研究的课题。
关于模型解释性，除了线性模型和决策树这种天生就有很好解释性的模型以外，sklean中有很多模型都有importance这一接口，可以查看特征的重要性。其实这已经含沙射影地体现了模型解释性的理念。只不过传统的 importance 的计算方法其实有很多争议，且并不总是一致。
模型可解释性工具 SHAP
Shapley value最早由加州大学洛杉矶分校（UCLA）的教授 Lloyd Shapley 提出, 主要是用来解决合作博弈论中的分配均衡问题。Lloyd Shapley 是 2012 年的诺贝尔经济学奖获得者，也是博弈论领域的无冕之王。
SHAP，作为一种经典的事后解释框架，可以对每一个样本中的每一个特征变量，计算出其重要性值，达到解释的效果。该值在 SHAP 中被专门称为Shapley Value。因此 Shapley Value 是 SHAP 方法的核心所在，理解好该值背后的含义将大大有助于我们理解 SHAP 的思想。
我们已有的数据集中会包含很多特征变量，从博弈论的角度，可以把每一个特征变量当成一个玩家。用该数据集去训练模型得到的预测结果， 可以看成众多玩家合作完成一个项目的收益。Shapley value，通过考虑各个玩家做出的贡献，来公平地分配合作的收益。
支持SHAP的模型
SHAP 可以用来解释很多模型。目前支持的模型有：
Tree Explainer（树模型）：Tree Explainer 首先根据树结构将待解释样本划分为一个叶子节点，并计算该叶子节点对预测结果的贡献值。然后，从叶子节点开始向上遍历树结构，计算每个节点对预测结果的贡献差异，并根据特征的排序顺序依次累加每个特征的贡献，得到其 Shapley value。
Kernel Explainer（核函数）：Kernel Explainer 首先通过随机抽样获取一组样本，并在这些样本上拟合一个全局模型。然后，在每个样本上，根据每个特征的取值和全局模型的预测值，构建一个新的局部模型。最后，通过核函数对每个局部模型进行加权平均，得到基于全局模型的 Shapley value 近似值。
Deep Explainer（神经网络模型）：DeepExplainer 会首先通过正向传播计算出模型在待解释样本上的预测值和激活值。然后，通过反向传播计算各层梯度，将梯度乘以激活值，并根据 Shapley value 对每个神经元进行权重分配，得到每个特征的重要性贡献。
Gradient Explainer（梯度信息）：GradientExplainer 首先选择一个参考点作为基准，并在该点处计算模型输出的梯度。然后，在每个样本上，通过对输入特征进行微小扰动来计算新的输入数据点，并计算模型输出的梯度差异。根据梯度的加权平均值，可以得到每个特征的 Shapley value 近似值。
Linear Explainer（线性模型）：LinearExplainer 首先通过拟合一个线性模型来计算每个特征对预测结果的贡献值。然后，根据特征的排序顺序依次累加每个特征的贡献，得到其 Shapley value。
SHAP解释
导入Python库
import pandas as pd #python科学计算库
import numpy as np #Python的一个开源数据分析处理库。
import matplotlib.pyplot as plt #常用Python画图工具
from xgboost import XGBRegressor # 导入 XGBRegressor 模型
from sklearn.model_selection import train_test_split # 数据划分模块
from sklearn.preprocessing import StandardScaler # 标准化模块
from sklearn.metrics import mean_squared_error, r2_score #误差函数MSE,误差函数R^2,
from sklearn.model_selection import GridSearchCV #超参数网格搜索
import shap # 导入SHAP模型解释工具
导入数据并标准化
data = pd.read_excel('D:/复现/trainset_loop6.xlsx') #读取xlsx格式数据
# date = pd.read_csv('D:/复现/trainset_loop6.csv') #读取csv格式数据
print(data.isnull().sum()) #检查数据中是否存在缺失值
print(data.shape) #检查维度
print(data.columns) #数据的标签
data = data.drop(["PN","AN"], axis = 1) #axis = 1表示对列进行处理，0表示对行
Y, X = data['Eads'] , data.drop(['Eads'] , axis = 1) #对Y、X分别赋值
columns = X.columns
scaler = StandardScaler()
X = scaler.fit_transform(X)
X_train,X_test,y_train,y_test = train_test_split(X , Y , test_size=0.2 , random_state=42)
XGBoost模型训练
model = XGBRegressor() # 模型实例化。
n_estimators = [50,100,150]
max_depth = [2,3,4]
reg_alpha = [0.1,0.5,1.0]
reg_lambda = [0.1,0.5,1.0]
grid_search = GridSearchCV(estimator = model,cv = 5, param_grid={'max_depth': max_depth,"n_estimators":n_estimators,"reg_alpha":reg_alpha,"reg_lambda":reg_lambda}, scoring='neg_mean_squared_error')
grid_search.fit(X_train, y_train)
best_model = grid_search.best_estimator_
best_max_depth = grid_search.best_params_['max_depth']
best_n_estimators = grid_search.best_params_['n_estimators']
best_reg_alpha = grid_search.best_params_['reg_alpha']
best_reg_lambda = grid_search.best_params_['reg_lambda']
print(f'Best max_depth:{best_max_depth:.3f}',"\n",f'Best n_estimators:{best_n_estimators:.3f}', # ,"\n",用于换行
"\n",f'Best reg_alpha:{best_reg_alpha:.3f}',"\n",f'Best reg_lambda:{best_reg_lambda:.3f}')
SHAP解释模型
explainer = shap.TreeExplainer(best_model) # 传入训练好的模型。
shap_values = explainer.shap_values(X_test) # 这里拿验证数据集进行呈现。输入X_train拿训练数据集进行呈现。
shap.summary_plot(shap_values = shap_values, # 一个数组，其中包含了样本的 Shapley 值，大小为 [n_samples, n_features]
features = X_test , # 一个数组，其中包含了样本的特征矩阵。这需要与 shap_values 中的样本数量一致。
feature_names = columns, # 特征名称列表，长度应该和 features 的列数相同。
max_display = 25, # 可选参数，用于指定要显示的最多特征数量。默认情况下，将显示所有特征。
plot_type = "dot", # 可选参数，用于指定图形类型。可以设置为 ‘dot’ 或 ‘bar’。
)
shap.summary_plot（）的超参数：
- shap_values ：一个数组，其中包含了样本的Shapley 值，大小为 [n_samples, n_features]
- features ：一个数组，其中包含了X_test所有值。需要与 shap_values 的样本数量一致。
- feature_names：特征名称列表，长度应该和 features 的列数相同。
- max_display（default=all）：用于指定要显示的最多特征数量。
- plot_type: 用于指定图形类型。可以设置为 ‘dot’ 或 ‘bar’。
蜂巢与特征重要性图结合
fig, ax1 = plt.subplots()
shap.summary_plot(shap_values, X_test, feature_names=columns,max_display = 25, plot_type="dot", show=False, color_bar=True)
plt.gca().set_position([0.2, 0.2, 0.65, 0.65]) # 调整图表位置，留出右侧空间放热度条
ax1 = plt.gca()
ax2 = ax1.twiny()
shap.summary_plot(shap_values, X_test,feature_names=columns, max_display = 25,plot_type="bar", show=False)
plt.gca().set_position([0.2, 0.2, 0.65, 0.65]) # 调整图表位置，与蜂巢图对齐
ax2.axhline(y=25, color='gray', linestyle='-', linewidth=1) # 注意y值应该对应顶部
bars = ax2.patches # 获取所有的柱状图对象
for bar in bars:
bar.set_alpha(0.3)# 设置透明度
ax1.set_xlabel('Shapley Value Contribution (Bee Swarm)', fontsize=12)
ax2.set_xlabel('Mean Shapley Value (Feature Importance)', fontsize=12)
ax2.xaxis.set_label_position('top') # 将标签移动到顶部
ax2.xaxis.tick_top() # 将刻度也移动到顶部
ax1.set_ylabel('Features', fontsize=12)
plt.tight_layout()
plt.savefig("SHAP_combined_with_top_line_corrected.jpg", format='jpg', bbox_inches='tight')
plt.show()
参考链接
- Welcome to the SHAP documentation — SHAP latest documentation
- 机器学习模型可解释性
- Exercise: Advanced Uses of SHAP Values | Kaggle