R语言绘制SHAP图

SHAP图是使用SHAP值生成的图形，用于展示机器学习模型预测结果中各个特征的重要性及其影响。要说提高模型结果的可解释性，SHAP图绝对是一把好手。这次，我们借助R语言的tidymodels、shapviz包来绘制XGBoost模型预测结果的SHAP图。

图1 模型预测与对应的SHAP图

基本R包

tidymodels

tidymodels 是一个R语言的元包（meta-package），它集合了一系列用于建模和预测的包，旨在提供一个统一的、面向数据科学的机器学习工作流程。这个集合包括了数据预处理（如recipes）、模型训练（如parsnip）、模型评估（如yardstick）、模型选择（如tune和workflows）等各个方面。tidymodels 的设计哲学是强调整洁的数据和代码，与tidyverse系列包保持一致的语法和风格，使得数据科学家能够以一种连贯和可重复的方式处理复杂的机器学习问题。

shapviz

shapviz 是一个R语言的包，专门用于可视化SHAP（SHapley Additive exPlanations）值。SHAP值是一种用于解释机器学习模型预测结果的工具，它通过将模型的预测归因于每个特征的影响来工作。shapviz 提供了一系列函数，帮助用户生成易于理解的图形，以展示哪些特征对模型的预测有重要影响，以及这些影响的方向和大小。这对于提高模型的可解释性和透明度非常有用。

SHAP图

原理

SHAP（SHapley Additive exPlanations）值基于博弈论中的Shapley值，用于估计每个特征对模型预测结果的贡献。SHAP值的核心思想是通过考虑所有可能的特征组合，计算每个特征在不同组合中的边际贡献。具体来说，它考虑了特征在加入或移除时对模型预测结果的影响，并据此分配贡献度。

用处

SHAP图提供了一种直观的方式来理解机器学习模型的预测。通过查看SHAP图，我们可以快速识别哪些特征对模型的预测结果有重要影响，以及这些影响是正面的还是负面的。这对于理解模型的决策过程、识别潜在的偏差或错误、以及改进模型的可解释性和透明度非常有帮助。

PimaIndiansDiabetes数据集

我们这次使用R中mlbench包中的PimaIndiansDiabetes数据集，它是一个经典的糖尿病相关数据集，该数据集最初来自国家糖尿病/消化/肾脏疾病研究所，包含了Pima印第安人的医疗记录信息。目标是基于数据集中包含的某些诊断测量来预测患者是否患有糖尿病。这些实例是从一个较大的数据库中根据特定条件选择出来的，特别是这里的所有患者都是Pima印第安至少21岁的女性。

该数据集一共有768个样本，8个特征变量和一个目标变量，具体如下：

• 特征变量：
  1. pregnant：怀孕次数
  2. glucose：口服葡萄糖耐量试验中2小时的血糖浓度
  3. pressure：舒张压（mmHg）
  4. triceps：三头肌皮肤褶皱厚度（mm）
  5. insulin：2小时血清胰岛素（mu U/ml）
  6. mass：体重指数（体重（kg）/身高（m）^2）
  7. pedigree：糖尿病谱系功能
  8. age：年龄（岁）
• 目标变量：
  1. diabetes：是否患有糖尿病

代码

下面是绘制SHAP图的具体代码，由于代码注释比较详细，就不进行文字描述了。

# 加载必要的库
library(dplyr)       # 数据处理
library(mlbench)     # 加载数据集
library(shapviz)     # SHAP值计算与可视化
library(tidyverse)   # 数据操作
library(tidymodels)  # 包含了rsample等相关包

# 替换为你自己的工作目录
setwd('C:/Users/Administrator/Desktop/检验')

# 加载数据
data(PimaIndiansDiabetes)
data= PimaIndiansDiabetes
# 设置因变量
target<- 'diabetes'
# 转变为因子水平
factor_columns<- c('pregnant', 'diabetes')
data<- data%>%
  mutate(across(all_of(factor_columns), as.factor))
# 创建一个函数来转换符合条件的列  如果因子类型只有两个水平，就将其转变为0，1取值的数值变量
data_cols<- names(data)
for(col in data_cols){
  if(col!= target&& is.factor(data[[col]])&& n_distinct(data[[col]])== 2){
    # 转换并替换原始列
    data[[col]]<- ifelse(data[[col]]== levels(data[[col]])[1], 0, 1)
  }
}

# 分割数据集为训练集和测试集
set.seed(4321)
data_split<- initial_split(data, prop= 0.8, strata= diabetes)
train_data<- training(data_split)
test_data<- testing(data_split)

# 创建数据预处理的流程
# 模型要预测target, 用其他变量进行预测
data_recipe<- recipe(as.formula(paste(target, "~ .")), data = train_data) %>%
  step_unknown(all_nominal_predictors()) %>%  # step_unknown()步骤来处理所有名义型（分类）预测变量中的未知（或新出现）类别；将未知类别编码为“未知”
  step_dummy(all_nominal_predictors()) %>%  # 将分类变量转换为虚拟变量
  step_zv(all_predictors()) %>%  # step_zv()步骤来识别并移除（或处理）所有预测变量中的零方差变量。零方差变量是指那些在所有观测中取值都相同的变量，它们不提供任何有关目标变量的信息，因此可以安全地从数据中移除。这个步骤有助于减少模型的复杂度，提高模型的泛化能力。
  step_normalize(all_numeric_predictors())  # 标准化数值型变量
# 预处理步骤应用于指定的训练数据集，train_data
trained_recipe<- prep(data_recipe, training= train_data)
# 使用训练好的recipe对训练和测试数据进行预处理
train_data_prepped<- bake(trained_recipe, new_data= train_data)
test_data_prepped<- bake(trained_recipe, new_data= test_data)

# 定义提升树模型
boosted_tree_model<- boost_tree(
  mode= "classification",
  trees= 1000,
  mtry= tune(),
  min_n= tune(),
  tree_depth= tune(),
  learn_rate= tune(),
  loss_reduction= tune(),
  sample_size= tune()
)%>%
  set_engine('xgboost')

# 快速生成数据框对应的统计信息
skimr::skim(train_data_prepped)
skimr::skim(test_data_prepped)

# 从测试数据中删除标签列
test_data_prepped<- test_data_prepped%>% select(-all_of(target))
# 将所有因子型变量转换为数值型并转换为矩阵
test_data_matrix<- as.matrix(test_data_prepped%>% mutate(across(where(is.factor), as.numeric)))
# 创建workflow
# 首先，通过调用workflow()函数（不带任何参数）来创建一个空的工作流程对象，并使用管道操作符%>%来链式调用后续的函数。
# 然后，使用add_model()函数将指定的模型（在这个例子中是boosted_tree_model）添加到工作流程中。
# 最后，使用add_recipe()函数将之前定义的数据预处理配方（data_recipe）添加到工作流程中。
workflow<-
  workflow()%>%
  add_model(boosted_tree_model)%>%
  add_recipe(data_recipe)

# 调优模型
tune_res<- tune_grid(
  workflow,
  resamples= vfold_cv(train_data, v= 5),
  grid= 10
)
# 根据accuracy指标，选择最佳参数
best_params<- select_best(tune_res, metric= "accuracy")
# 训练最终模型
final_workflow<- finalize_workflow(workflow, best_params)
# 这里，final_workflow是一个之前已经定义并可能包含了模型规格（boosted_tree_model）和数据预处理配方（data_recipe）的工作流程对象。
final_model<- fit(final_workflow, data= train_data)

# 从最终模型中提取 xgboost 模型对象
xgboost_model<- extract_fit_engine(final_model)
# 计算SHAP值
shap_values<- shapviz(xgboost_model, X_pred= test_data_matrix)
# 绘制 SHAP 值图
# kind: 指定要显示的图表类型。 both 表示同时显示条形图和散点图。bar 仅显示条形图，beeswarm 仅显示散点图。
# SHAP value不显示，设置show_numbers= False
sv_importance(shap_values, kind= "both", show_numbers= TRUE, color_bar_title= 'Feature value')

相应的SHAP图，如图2所示。

图2 XGBoost预测结果对应的SHAP图