利用SHAP解释Xgboost模型

Xgboost相对于线性模型在进行预测时往往有更好的精度,但是同时也失去了线性模型的可解释性。所以Xgboost通常被认为是黑箱模型。

2017年,Lundberg和Lee的论文提出了SHAP值这一广泛适用的方法用来解释各种模型(分类以及回归),其中最大的受益者莫过于之前难以被理解的黑箱模型,如xgboost和神经网络模型。

本教程中,我们在真实数据集上进行实操,利用SHAP来解释Xgboost模型。

预计学习用时:30分钟。

本教程基于Python 3.6版本、Xgboost 0.82版本以及shap 0.28.5版本。

原创者:东布东 | 修改校对:SofaSofa TeamM |…

http://sofasofa.io/tutorials/shap_xgboost/

1. Feature importance

在这里插入图片描述

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

# 读取数据,目标变量y是球员的身价(万欧元)
data = pd.read_csv('train.csv')

# 获得当时球员年龄
today = pd.to_datetime('2018-01-01')
data['birth_date'] = pd.to_datetime(data['birth_date'])
data['age'] = np.round((today - data['birth_date']).apply(lambda x: x.days) / 365., 1)

# 选择特征,这里只是举例,未必是最佳组合
# 特征依次为身高(厘米)、潜力、速度、射门、传球、带球、防守、体格、国际知名度、年龄
cols = ['height_cm', 'potential', 'pac', 'sho', 'pas', 'dri', 'def', 'phy', 'international_reputation', 'age']

# 训练xgboost回归模型
model = xgb.XGBRegressor(max_depth=4, learning_rate=0.05, n_estimators=150)
model.fit(data[cols], data['y'].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()

在这里插入图片描述
在这里插入图片描述

2. SHAP value

在这里插入图片描述

3. SHAP的Python实现

Python中SHAP值的计算由shap这个package实现,可以通过pip install shap安装。

下面我们针对第1节中训练出的模型model,计算其SHAP值。

引用package并且获得解释器explainer。

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

在这里插入图片描述

shap_values = explainer.shap_values(data[cols])
print(shap_values.shape)

在这里插入图片描述
通过对比发现,我们可以确认基线值就是训练集的目标变量的拟合值的均值。在这里例子中,目标变量是球员的身价(万欧元),也就是球员的平均身价为229万欧元。
在这里插入图片描述

y_base = explainer.expected_value
print(y_base)

data['pred'] = model.predict(data[cols])
print(data['pred'].mean())

在这里插入图片描述

3.1 单个样本的SHAP值

我们可以随机检查其中一位球员身价的预测值以及其特征对预测值的影响。

下面的数据框中第一列是特征名称,第二列是特征的数值,第三列是各个特征在该样本中对应的SHAP值。

# 比如我们挑选数据集中的第30位
j = 30
player_explainer = pd.DataFrame()
player_explainer['feature'] = cols
#某行
player_explainer['feature_value'] = data[cols].iloc[j].values
player_explainer['shap_value'] = shap_values[j]
player_explaine

在这里插入图片描述
我们知道一个样本中各特征SHAP值的和加上基线值应该等于该样本的预测值。

我们可以做如下的验证。

我们知道一个样本中各特征SHAP值的和加上基线值应该等于该样本的预测值。

我们可以做如下的验证。

print('y_base + sum_of_shap_values: %.2f'%(y_base + player_explainer['shap_value'].sum()))
print('y_pred: %.2f'%(data['pred'].iloc[j]))

在这里插入图片描述
shap还提供极其强大的数据可视化功能。下图是对上面数据框的可视化。

指定其中某行进行特征的分析

蓝色表示该特征的贡献是负数,红色则表示该特征的贡献是正数。最长的红色条是潜力值,球员的潜力值很高,而他的身价也因此增加了1092万;最长的蓝色条是年龄,这个球员年龄较小才20岁出头,尚未到职业巅峰,未来也有诸多不确定性,身价也因此降低了180万元

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

在这里插入图片描述

3.2 对特征的总体分析

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

下图中每一行代表一个特征,横坐标为SHAP值。一个点代表一个样本,颜色越红说明特征本身数值越大,颜色越蓝说明特征本身数值越小。

我们可以直观地看出潜力potential是一个很重要的特征,而且基本上是与身价成正相关的。年龄age也会明显影响身价,蓝色点主要集中在SHAP小于0的区域,可见年纪小会降低身价估值,另一方面如果年纪很大,也会降低估值,甚至降低得更明显,因为age这一行最左端的点基本上都是红色的。

shap.summary_plot(shap_values, data[cols])

在这里插入图片描述
我们也可以把一个特征对目标变量影响程度的绝对值的均值作为这个特征的重要性

因为SHAP和feature_importance的计算方法不同,所以我们这里也得到了与第1节不同的重要性排序。

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

在这里插入图片描述

3.3 部分依赖图Partial Dependence Plot

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

SHAP也提供了部分依赖图的功能,与传统的部分依赖图不同的是,这里纵坐标不是目标变量y的数值而是SHAP值。

比如下图中,年纪大概呈现出金字塔分布,也就是24到31岁这个年纪对球员的身价是拉抬作用,小于24以及大于31岁的球员身价则会被年纪所累。
在这里插入图片描述

3.4 对多个变量的交互进行分析

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

在这里插入图片描述

我们也可以用dependence_plot描绘两个变量交互下变量对目标值的影响。

#两个列名  一个potential一个是 international——reputation
shap.dependence_plot('potential', shap_values, data[cols], interaction_index='international_reputation', show=False)

在这里插入图片描述

好的,我会为您提供相关的R语言代码。 首先,我们需要加载所需的包,包括xgboost、ROCR和shap。您可以使用以下命令安装这些包: ```r install.packages("xgboost") install.packages("ROCR") install.packages("shap") ``` 然后,我们可以使用以下代码读取数据集: ```r data <- read.csv("your_data.csv") ``` 接下来,我们需要将因子变量转换为虚拟变量。您可以使用以下代码实现: ```r data <- data.frame(lapply(data, function(x) if(is.factor(x)) as.numeric(as.character(x)) else x)) ``` 然后,我们将数据集分成训练集和测试集。您可以使用以下代码实现: ```r library(caret) set.seed(123) trainIndex <- createDataPartition(data$target, p = .7, list = FALSE, times = 1) train <- data[trainIndex,] test <- data[-trainIndex,] ``` 接下来,我们可以训练xgboost模型并进行预测。您可以使用以下代码实现: ```r library(xgboost) model <- xgboost(data = as.matrix(train[, -1]), label = train$target, nrounds = 100, objective = "binary:logistic", verbose = 0) pred <- predict(model, as.matrix(test[, -1])) ``` 然后,我们可以绘制ROC曲线图来评估模型性能。您可以使用以下代码实现: ```r library(ROCR) predObj <- prediction(pred, test$target) perf <- performance(predObj, measure = "tpr", x.measure = "fpr") plot(perf, colorize = TRUE) ``` 最后,我们可以使用shap包来解释xgboost模型。您可以使用以下代码实现: ```r library(shap) explainer <- shap.explainer(as.matrix(train[, -1]), model) shap_values <- explainer(as.matrix(test[, -1])) plot(shap_values, test[, -1]) ``` 希望这些代码能够帮助您解决问题。
评论 9
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值