R语言机器学习与临床预测模型69--机器学习模型解释利器：SHAP

R小盐准备介绍R语言机器学习与预测模型的学习笔记， 快来收藏关注【科研私家菜】

01 机器学习的可解释性

对于集成学习方法，效果虽好，但一直无法解决可解释性的问题。我们知道一个xgboost或lightgbm模型，是由N棵树组成，所以对于特定的一个样本，我们无法知道这个样本的特征值是如何影响最终结果。虽说“不管白猫黑猫，抓住耗子的就是好猫”，但在具体任务中，我们还是希望能够获得样本每个特征与其结果之间的关系，特别是针对模型误分的那些样本，如果能够从特征和结果的角度进行分析，对于提高模型效果或是分析异常样本，是非常有帮助的。但是，其可解释性相对困难。
对于集成树模型来说，当做分类任务时，模型输出的是一个概率值。前文提到，SHAP是SHapley Additive exPlanations的缩写，即沙普利加和解释，因此SHAP实际是将输出值归因到每一个特征的shapely值上，换句话说，就是计算每一个特征的shapley值，依此来衡量特征对最终输出值的影响。

其原理及推到公式不再赘述。。。

02 SHAP的R语言实现

SHAP（SHapley Additive exPlanations）

library(tidyverse)
library(xgboost)
library(caret)
library(dplyr)
source("shap.R")

bike <- read.csv("../shap-values-master/bike.csv",header = T)


bike_2=select(bike, -days_since_2011, -cnt, -yr)

bike_dmy = dummyVars(" ~ .", data = bike_2, fullRank=T)
bike_x = predict(bike_dmy, newdata = bike_2)

## Create the xgboost model
model_bike = xgboost(data = bike_x, 
                   nround = 10, 
                   objective="reg:linear",
                   label= bike$cnt)  


## Calculate shap values
shap_result_bike = shap.score.rank(xgb_model = model_bike, 
                              X_train =bike_x,
                              shap_approx = F
                              )

# `shap_approx` comes from `approxcontrib` from xgboost documentation. 
# Faster but less accurate if true. Read more: help(xgboost)

## Plot var importance based on SHAP
var_importance(shap_result_bike, top_n=15)

## Prepare data for top N variables
shap_long_bike = shap.prep(shap = shap_result_bike,
                           X_train = bike_x , 
                           top_n = 10
                           )

## Plot shap overall metrics
plot.shap.summary(data_long = shap_long_bike)


## 
xgb.plot.shap(data = bike_x, # input data
              model = model_bike, # xgboost model
              features = names(shap_result_bike$mean_shap_score[1:10]), # only top 10 var
              n_col = 3, # layout option
              plot_loess = T # add red line to plot
              )

效果如下：

03 SHAP R语言示例

data("iris")
X1 = as.matrix(iris[,-5])
mod1 = xgboost::xgboost(
  data = X1, label = iris$Species, gamma = 0, eta = 1,
  lambda = 0, nrounds = 1, verbose = FALSE)

# shap.values(model, X_dataset) returns the SHAP
# data matrix and ranked features by mean|SHAP|
shap_values <- shap.values(xgb_model = mod1, X_train = X1)
shap_values$mean_shap_score
shap_values_iris <- shap_values$shap_score

# shap.prep() returns the long-format SHAP data from either model or
shap_long_iris <- shap.prep(xgb_model = mod1, X_train = X1)
# is the same as: using given shap_contrib
shap_long_iris <- shap.prep(shap_contrib = shap_values_iris, X_train = X1)

# **SHAP summary plot**
shap.plot.summary(shap_long_iris, scientific = TRUE)
shap.plot.summary(shap_long_iris, x_bound  = 1.5, dilute = 10)

# Alternatives options to make the same plot:
# option 1: from the xgboost model
shap.plot.summary.wrap1(mod1, X = as.matrix(iris[,-5]), top_n = 3)

# option 2: supply a self-made SHAP values dataset
# (e.g. sometimes as output from cross-validation)
shap.plot.summary.wrap2(shap_score = shap_values_iris, X = X1, top_n = 3)

效果如下：

关注R小盐，关注科研私家菜(VX_GZH: SciPrivate)，有问题请联系R小盐。让我们一起来学习 R语言机器学习与临床预测模型