图片

图片

图片

图片

图片

图片

图片

图片

SHAP（SHapley Additive exPlanations）是一种模型无关的解释方法，因此它适用于任何模型。前面已经介绍过多次这个SHAP了，在R语言中实现SHAP也是非常简单的。

R语言SHAP模型解释R语言shapviz实现SHAP可视化R语言SHAP模型解释之kernelshap

fastshap基于Rcpp和foreach，所以它的速度非常快（比之前介绍的几种都快），而且fastshap既可以实现局部解释，又能实现全局解释，还可以无缝对接shapviz的可视化。

安装

# Install the latest stable version from CRAN:install.packages("fastshap")# Install the latest development version from GitHub:if (!requireNamespace("remotes")) {  install.packages("remotes")}remotes::install_github("bgreenwell/fastshap")

准备数据和R包

我们使用经过缺失值插补的泰坦尼克号数据集t1，该数据有1309行，6列，其中survived是结果变量，二分类，1代表死亡，2代表存活。

library(fastshap)t1 <- titanic_mice[[1L]]dim(t1)## [1] 1309    6str(t1)## 'data.frame':    1309 obs. of  6 variables:##  $ survived: Factor w/ 2 levels "no","yes": 2 2 1 1 1 2 2 1 2 1 ...##  $ pclass  : int  1 1 1 1 1 1 1 1 1 1 ...##  $ age     : num  29 0.92 2 30 25 48 63 39 53 71 ...##  $ sex     : Factor w/ 2 levels "female","male": 1 2 1 2 1 2 1 2 1 2 ...##  $ sibsp   : int  0 1 1 1 1 0 1 0 2 0 ...##  $ parch   : int  0 2 2 2 2 0 0 0 0 0 ...

把pclass变为有顺序的因子型：

t1$pclass <- as.ordered(t1$pclass)  # makes more sense as an ordered factor

本次使用ranger建立随机森林模型，对于其他模型，fastshap都是支持的。

library(ranger)set.seed(2053)  # for reproducibility(rfo <- ranger(survived ~ ., data = t1, probability = TRUE))## Ranger result## ## Call:##  ranger(survived ~ ., data = t1, probability = TRUE) ## ## Type:                             Probability estimation ## Number of trees:                  500 ## Sample size:                      1309 ## Number of independent variables:  5 ## Mtry:                             2 ## Target node size:                 10 ## Variable importance mode:         none ## Splitrule:                        gini ## OOB prediction error (Brier s.):  0.1337358

局部解释

为了说明如何最简单地使用Shapley值来量化特征贡献，我们需要一个新的观测值来预测。下面我们将为新建一个观测值，给这个观测取个名字叫jack.dawson：

jack.dawson <- data.frame(  #survived = 0L,  # in case you haven't seen the movie  pclass = 3L,     # third-class passenger  age = 20.0,      # twenty years old  sex = factor("male", levels = c("female", "male")),  # male  sibsp = 0L,      # no siblings/spouses aboard  parch = 0L       # no parents/children aboard)

使用fastshap进行模型解释和DALEX一样，首先也是需要建立一个解释器，用来解析模型和数据的各种信息。

在DALEX中，预测函数默认会从模型对象中提取（默认使用predict()），但是在fastshap中需要我们自定义预测函数，用于直接返回数值（回归任务）或者类别概率（分类任务）。

# 自定义预测函数pfun <- function(object, newdata) {  # prediction wrapper  unname(predict(object, data = newdata)$predictions[, "yes"])}# 计算Jack's的生存概率(jack.prob <- pfun(rfo, newdata = jack.dawson))## [1] 0.1314723

如果直接使用以下代码得到的并不是类别概率：

predict(rfo, jack.dawson)## Ranger prediction## ## Type:                             Probability estimation ## Sample size:                      1 ## Number of independent variables:  5

用这个函数可以得到所有t1观测的平均生存概率：

# 计算所有人的平均生存概率(baseline <- mean(pfun(rfo, newdata = t1)))  ## [1] 0.3821045

# Difference between Jack and average(difference <- jack.prob - baseline)## [1] -0.2506322

可以发现Jack的生存可能性低于平均值。下面我们就使用SHAP来解释为什么会这样。

下面就可以建立解释器了，和DALEX一样，也是使用explain()函数，也是需要提供只含预测变量的数据框：

X <- subset(t1, select = -survived)  # 只含预测变量set.seed(2113)# 建立解释器(ex.jack <- explain(rfo, X = X, pred_wrapper = pfun, newdata = jack.dawson))##      pclass          age sex      sibsp parch## [1,]      0 -0.005012306   0 0.02174902     0## attr(,"baseline")## [1] 0## attr(,"class")## [1] "explain" "matrix"  "array"

fastshap包使用高效版本的蒙特卡洛（Monte-Carlo，MC）算法。因此，为了稳定性和准确性，应多次计算特征贡献，并将结果取平均值。为此，只需将nsim参数设置为一个比较大的值即可（默认是1）。下面我们计算Jack的1000个基于Shapley的特征贡献，并获取平均结果：

set.seed(2129)  (ex.jack <- explain(rfo, X = X, pred_wrapper = pfun, newdata = jack.dawson,                    nsim = 1000))##           pclass         age        sex       sibsp       parch## [1,] -0.07554003 -0.01240914 -0.1414107 0.001836116 -0.01103988## attr(,"baseline")## [1] 0## attr(,"class")## [1] "explain" "matrix"  "array"

fastshap使用的MC方法计算的Shapley值的加和不会等于相应预测和基线（即平均预测值）之间的差值。但是借用Python的shap库的技巧，我们可以使用基于回归的调整来校正总和。为此，只需在调用explain()中设置为adjust = TRUE：

set.seed(2133)  (ex.jack.adj <- explain(rfo, X = X, pred_wrapper = pfun, newdata = jack.dawson,                        nsim = 1000, adjust = TRUE))##          pclass         age        sex       sibsp       parch## [1,] -0.0697378 -0.02354202 -0.1485205 0.003980237 -0.01281207## attr(,"baseline")## [1] 0.3821045## attr(,"class")## [1] "explain" "matrix"  "array"

这个结果的加和与difference是不一样的：

sum(ex.jack.adj)  # 和jack.prob - baseline的结果不一样哦## [1] -0.2506322

构建好解释器之后，就可以使用shapviz包进行可视化了。

我们先创建一个简单的瀑布图来可视化这几个特征是如何让Jack的预测生存概率较低的：

library(shapviz)shv <- shapviz(ex.jack.adj, X = jack.dawson, baseline = baseline)sv_waterfall(shv)

图片

显然，Jack是一名男性三等舱乘客，这两个变量对他的较低生存概率贡献最大。

然后是瀑布图的另一种形式：

sv_force(shv)

图片

从图形理解方面来看，这个图虽然和瀑布图表达的意思完全一样，但是不如瀑布图读起来简单。

全局解释

SHAP除了用于局部解释外，还可以用于全局解释。

如果对数据集中所有的观测都进行一遍SHAP解释，然后聚合它们的结果，就可以得到全局的SHAP解释。这个方法和CP聚合成为PDP的方法非常类似。

下面的代码使用1000次MC重复计算训练数据中每个乘客的Shapley解释，并将生成的矩阵强制转换为tibble（以便更好地打印）。设置shap_only=FALSE可以方便shapviz使用。

这个过程很慢哈~

set.seed(2224)  ex.t1 <- explain(rfo, X = X, pred_wrapper = pfun, nsim = 100, adjust = TRUE,                 shap_only = FALSE)tibble::as_tibble(ex.t1$shapley_values)## # A tibble: 1,309 × 5##    pclass      age     sex     sibsp    parch##     <dbl>    <dbl>   <dbl>     <dbl>    <dbl>##  1  0.228  0.00994  0.312   0.0208   -0.00916##  2  0.138  0.331   -0.0717  0.0113    0.0712 ##  3  0.156  0.0213   0.116  -0.0240   -0.0189 ##  4  0.212 -0.0212  -0.183   0.0171    0.00782##  5  0.201 -0.0308   0.282  -0.0149   -0.0369 ##  6  0.168 -0.0393  -0.199  -0.000779 -0.00240##  7  0.177 -0.127    0.346  -0.00829   0.00371##  8  0.153 -0.0644  -0.186  -0.00392  -0.00976##  9  0.234  0.00345  0.296   0.0398    0.00432## 10  0.106 -0.115   -0.212   0.000223 -0.00402## # ℹ 1,299 more rows

可视化全局变量重要性（注意和基于重排的变量重要性区分），这个其实是各个变量的多个Shapley值绝对值的平均值：

shv.global <- shapviz(ex.t1)sv_importance(shv)  

图片

变量重要性图的另一种表现形式，蜂窝图，在Python中被称为shap summary plot：

sv_importance(shv.global, kind = "beeswarm")

图片

变量依赖图，类似于部分依赖图。在这里，我们将看看特征贡献age对其输入值的依赖性：

sv_dependence(shv.global, v = "age")

图片

其他可视化函数可参考关于shapviz的推文：R语言shapviz可视化SHAP

如果你的数据量很大，fastshap还支持使用并行化，借助doParallel包即可实现，大家自己尝试下即可。

本站仅提供存储服务，所有内容均由用户发布，如发现有害或侵权内容，请点击举报。