之前介绍了SHAP这种模型解释方法,并且给大家详细演示了如何基于DALEX
,用3行代码实现对任意模型的局部(instance-level)SHAP解释,并介绍了结果的可视化,包括默认的方法以及自己通过ggplot2
的画法。
SHAP作为一种经典的变量归因模型解释方法,它和分解解释有着非常紧密的联系,如果你了解了分解解释的原理,那么SHAP解释就非常简单了。其实关于它的内容非常多,目前在R中实现SHAP的方法也有非常多,只用一篇推文是不可能介绍完的哈,后续会持续介绍这个SHAP。
今天给大家介绍一个专门用于SHAP可视化的R包,也是属于DrWhy.AI
生态系统中的一个R包,名字叫shapviz
。
公众号后台回复模型解释即可获取模型解释合集,回复shap即可获取SHAP解释合集。
本文目录:
安装
简介
示例
对接DALEX
# From CRAN
install.packages("shapviz")
# Or the newest version from GitHub:
# install.packages("devtools")
devtools::install_github("ModelOriented/shapviz")
shapviz
主要提供以下几种可视化方法:
sv_importance()
:
变量重要性可视化,可选条形图或者蜂窝图(也就是Python中常见的图)sv_dependence()
、sv_dependence2D()
: 变量依赖图sv_interaction()
: 变量交互效应图sv_waterfall()
: 瀑布图,最常规的SHAP累加归因图sv_force()
: 瀑布图的一种变体如果要使用以上几种可视化方法,需要提供shapviz
对象,这个对象需要两个东西:
S
:SHAP值矩阵X
:包含特征值的数据集为了方便使用,直接支持以下模型的输出:
XGBoost
LightGBM
h2o
kernelshap
fastshap
shapr
treeshap
DALEX
可以看到其中有DALEX
哈,因为DALEX
基本上支持所有的模型,也就是说这个shapviz
支持基本上所有模型。
使用xgboost
建立xgboost模型进行演示。使用diamonds
数据集,共有53940行,10列,回归任务,price
是结果变量,数值型。
library(shapviz)
library(ggplot2)
library(xgboost)
dim(diamonds)
## [1] 53940 10
str(diamonds)
## tibble [53,940 × 10] (S3: tbl_df/tbl/data.frame)
## $ carat : num [1:53940] 0.23 0.21 0.23 0.29 0.31 0.24 0.24 0.26 0.22 0.23 ...
## $ cut : Ord.factor w/ 5 levels "Fair"<"Good"<..: 5 4 2 4 2 3 3 3 1 3 ...
## $ color : Ord.factor w/ 7 levels "D"<"E"<"F"<"G"<..: 2 2 2 6 7 7 6 5 2 5 ...
## $ clarity: Ord.factor w/ 8 levels "I1"<"SI2"<"SI1"<..: 2 3 5 4 2 6 7 3 4 5 ...
## $ depth : num [1:53940] 61.5 59.8 56.9 62.4 63.3 62.8 62.3 61.9 65.1 59.4 ...
## $ table : num [1:53940] 55 61 65 58 58 57 57 55 61 61 ...
## $ price : int [1:53940] 326 326 327 334 335 336 336 337 337 338 ...
## $ x : num [1:53940] 3.95 3.89 4.05 4.2 4.34 3.94 3.95 4.07 3.87 4 ...
## $ y : num [1:53940] 3.98 3.84 4.07 4.23 4.35 3.96 3.98 4.11 3.78 4.05 ...
## $ z : num [1:53940] 2.43 2.31 2.31 2.63 2.75 2.48 2.47 2.53 2.49 2.39 ...
建立xgboost模型:
# 建立xgboost模型
set.seed(3653)
x <- c("carat", "cut", "color", "clarity")
dtrain <- xgb.DMatrix(data.matrix(diamonds[x]), label = diamonds$price)
fit <- xgb.train(params = list(learning_rate = 0.1), data = dtrain, nrounds = 65)
构建shapviz
对象,只使用其中2000个观测:
dia_small <- diamonds[sample(nrow(diamonds), 2000), ]
# 构建shapviz对象
shp <- shapviz(fit, X_pred = data.matrix(dia_small[x]), X = dia_small)
这个shapviz
具体来说支持以下方式创建的对象:
fit
对第一个观测画出瀑布图:这个图的解释请参考之前的推文:分解解释
sv_waterfall(shp, row_id = 1)
瀑布图的另一种呈现方式,就是把几个条形放在一行中,和Python中的展示方式差不多:
sv_force(shp, row_id = 1)
对应的Python的SHAP库中的图形是这种:
看颜值我觉得Python更胜一筹!但是从图形解读的便捷性来看,还是瀑布图更加好理解。
可以聚合多个观测的平均值,比如同时选择color
为D
的样本,此时画出来的SHAP值是平均值:
sv_waterfall(shp, shp$X$color == "D")
我们选择对2000个样本进行SHAP解释,这样对每一个特征,就会得到多个shapley值,对这些值取平均值就可以得到基于SHAP的变量重要性(注意和基于重排的变量重要性区分),这样一来这里的SHAP就变成了dataset-level的模型解释:
sv_importance(shp)
用蜂窝图的形式呈现结果,这个图也是Python的SHAP中常见的shap summary plot:
每一个点都表示一个样本,横坐标是Shapley值,纵坐标是变量名字,颜色深浅表示变量数值大小。这张图可以很清晰的展示每个变量的Shapley值分布情况。比如对于carat
这个变量来说,它的值越小,钻石的价格越低(Shapley值是负的)。
sv_importance(shp, kind = "beeswarm")
还可以展示变量依赖图,非常类似于之前介绍过的CP曲线图。只不过这里的图纵坐标是SHAP值,横坐标依然是预测变量。比如下图就展示了SHAP值随着color
这个变量的改变而变化的情况:
sv_dependence(shp, v = "color")
多个变量同时展示也可以,只不过需要借助patchwork辅助拼图:
library(patchwork)
sv_dependence(shp, v = x) &
theme_gray(base_size = 9) &
ylim(-5000, 15000)
还可以同时展示两个变量对SHAP值的影响,此时的SHAP值是它们的加和:
sv_dependence2D(shp, x = "carat", y = c("clarity", "color"), alpha = 0.5)
如果有交互作用还可以展示主效应和交互效应的SHAP:
shp_i <- shapviz(
fit,
X_pred = data.matrix(dia_small[x]),
X = dia_small, interactions = TRUE
)
# caret主效应以及,它和其他变量的交互效应
sv_dependence(shp_i, v = "carat", color_var = x, interactions = TRUE) &
ylim(-6000, 13000)
还可以用蜂窝图展示所有的主效应和交互效应:
sv_interaction(shp_i) +
theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1))
为了展示,再演示一个SVM模型,使用e1071
包构建一个支持向量机模型:
library(e1071)
set.seed(3653)
diamonds <- diamonds[sample(nrow(diamonds),1000),] #只使用1000行
model_svm <- svm(price ~ carat+cut+color+clarity, data = diamonds)
由于shapviz
不能直接支持e1071
创建的模型对象,所以我们使用DALEX
转换一下。
建立解释器:
library(DALEX)
explain_svm <- DALEX::explain(model = model_svm,
data = diamonds[,1:4],
y = diamonds$price,
label = "randomforest"
)
## Preparation of a new explainer is initiated
## -> model label : randomforest
## -> data : 1000 rows 4 cols
## -> data : tibble converted into a data.frame
## -> target variable : 1000 values
## -> predict function : yhat.svm will be used ( default )
## -> predicted values : No value for predict function target column. ( default )
## -> model_info : package e1071 , ver. 1.7.14 , task regression ( default )
## -> predicted values : numerical, min = -161.2517 , mean = 3783.737 , max = 17316.85
## -> residual function : difference between y and yhat ( default )
## -> residuals : numerical, min = -5942.083 , mean = 40.21553 , max = 8179.3
## A new explainer has been created!
SHAP解释:
set.seed(3653)
shap_svm <- predict_parts(explainer = explain_svm,
new_observation = diamonds[1,1:4],
type = "shap",
B = 15
)
#shap_svm
然后就可以提供给shapviz
用了:
shp_svm <- shapviz(shap_svm)
接下来就可以使用shapviz
包中的各种函数画图了,不做只能使用其中一部分函数,尤其是需要多个观测SHAP的函数,这里是不能用的,因为predict_parts()
只能接受一个观测。
此时只有个别函数可用,比如:
sv_force(shp_svm)
sv_importance(shp_svm)
如果你非要使用多个观测实现SHAP(也就是全局解释的SHAP),还有许多方式其他R包可以实现哈,比如fastshap
、kernelshap
、treeshap
等,慢慢介绍。
联系客服