SHAP (SHapley Additive exPlanati)及DALEX预测单样本变量情况和shapviz可视化学习

SHAP (SHapley Additive exPlanati)及DALEX预测单样本变量情况和shapviz可视化学习

SHapley Additive exPlanati（SHAP） 是一种广泛用于解释机器学习模型预测的方法。由于目前所构建的模型都是非常复杂的，从而难以去解释其中的结果，因此就提出了一个预测解释框架——SHAP。

SHAP的优势

每种模型的行为和预测结果需要不同的解释方法。一些模型由于其结构和工作原理较为简单直接，因此可以视为“可自解释”的。例如，决策树、回归树、基于规则的模型，以及基于最近邻的方法，它们的决策过程较为透明，用户可从模型的结构直接理解其决策逻辑。此外，线性回归和其他加性模型可以通过展示每个输入变量的边际效应来进行解释，其中“边际效应”指的是在其他变量保持不变的情况下，一个变量改变一个单位对模型所造成的平均影响。这种方法有助于清晰地展示每个变量对模型预测结果的贡献，是解释和可视化线性回归及相关模型行为的常用手段。 通过这种方法，研究者或数据分析师能够更好地理解数据特征如何影响模型预测，进而作出更有效的决策或优化模型。

在机器学习中，尤其是对于如随机森林和梯度提升决策树等复杂的“黑箱”模型，虽然它们能在许多场景下取得卓越的预测效果，但这些模型的内部工作机制和预测逻辑往往不易于解释。SHAP是一种用于模型解释的工具，它通过为每个输入特征分配一个“归因值”来量化该特征对模型预测结果的贡献。 SHAP基于博弈论中的Shapley值，确保了解释的数学一致性和公平性。通过SHAP框架，研究者可以解释每个输入特征是如何影响模型输出的，从而使得通常被视为黑箱的模型变得更加透明和可解释。

因此，笔者个人的感觉是SHAP框架其实可以去解释任何简单和复杂，线性和非线性模型中的变量贡献，但是简单的模型没必要使用SHAP去解释，作为模型的构建者来说，我们必须要记住一句话：“The best explanation of a simple model is the model itself; it perfectly represents itself and is easy to understand”。

SHAP 的原理

SHAP的核心思想来源于博弈论中的Shapley值，它是一种用来衡量玩家在合作博弈中的贡献的方法。SHAP将模型的每一个特征看作博弈中的一个“玩家”，其贡献可以理解为模型输出变化的部分。

1. 假设我们有一个合作博弈，多个玩家一起合作产生一个总收益。Shapley值用于公平地分配这些收益，衡量每个玩家对最终结果的平均边际贡献。
2. 在SHAP中，玩家是特征变量，合作博弈的总收益是模型的预测值。
3. SHAP值是通过考虑所有特征组合的边际贡献来计算的，即考虑在有或没有某个特征的情况下，对模型预测输出的影响。

目前能够实现SHAP分析的工具有很多，本次我们来学习一下DALEX和shapviz这两个R包。

DALEX（Descriptive mAchine Learning EXplanati） 是一个专门用于解释机器学习模型行为的R包。它提供了一种统一的框架，可以帮助用户理解模型的预测结果以及模型的整体行为，尤其是对于复杂的黑箱模型（如随机森林、梯度提升决策树和神经网络）。

主要功能：

局部和全局解释：1. 局部解释（Local explanati）：针对单个预测结果，分析特定特征对预测值的贡献。2. 全局解释（Global explanati）：分析整个模型的行为，例如特征重要性或交互关系。

支持多种解释方法：1.特征重要性（Feature Importance)；2. 部分依赖图（Partial Dependence Plots, PDPs)；. 累积局部效应图（Accumulated Local Effects, ALEs)； 4. Shapley值（Shapley Values)；4. 模型诊断工具，如残差分析

相关拓展包：1. DALEXtra：扩展了DALEX，支持与更多的机器学习框架（如H2O、Keras、PyTorch）无缝集成。2. iBreakDown：专注于局部解释，如Shapley值或LIME。

主要函数

1. explain()： 用于包装一个模型，使其可以与DALEX的解释函数一起使用。它接受模型、数据和目标变量作为输入。
2. model_parts()： 计算全局特征重要性，显示每个特征对模型输出的总体影响。
3. model_profile()： 绘制特征依赖图（Partial Dependence Plots）或累积局部效应图（ALEs）。
4. model_performance()： 用于诊断模型性能，例如计算残差或绘制残差分布。
5. predict_parts()： 用于单个预测的局部解释，例如LIME或SHAP值。

shapviz 一个用于可视化和解释 SHAP (SHapley Additive exPlanati) 值的工具包。它旨在为用户提供简单而高效的可视化功能，帮助更直观地解释机器学习模型的预测

1. 局部解释（Local Explanation）：显示单个样本的特征对预测结果的具体贡献，例如：瀑布图（Waterfall Plot）：展示从基线值到最终预测值的累积贡献。力图（Force Plot）：类似瀑布图的另一种表示方法，更注重特征贡献方向。
2. 全局解释（Global Explanation）：分析特征在整个数据集中的影响，例如：特征重要性图（Feature Importance）：展示哪些特征对预测最重要。特征依赖图（Dependence Plot）：探索特定特征的值如何影响预测结果。交互作用图（Interaction Plot）：展示两个特征之间的交互效应。

分析步骤

1.导入

代码语言：javascript代码运行次数：0运行复制

rm(list = ls())
load("TCGA_HSC_practice.Rdata")

library(DALEX)
library(shapviz)

2.数据预处理

代码语言：javascript代码运行次数：0运行复制

# prepare model
colnames(meta)
head(meta)[1:,1:8]
#                                ID OS   gender  race site neoadjuvant age
# TCGA-CR-774-01A TCGA-CR-774-01A  0 1.000000 FEMALE WHITE   CR          o  67
# TCGA-CV-A45V-01A TCGA-CV-A45V-01A  1 1.066667 FEMALE WHITE   CV          o  87
# TCGA-CV-7102-01A TCGA-CV-7102-01A  1 1.866667 FEMALE WHITE   CV          o  76
meta$OS <- factor(meta$OS)

# 提取数据，去除A
# 要注意在建立模型的时候变量可以使非数字
# shap解释变量的时候一定要把变量内容设定为数字
dat <- (meta[,c(2,8,11,12,17,41:4)])
head(dat)
#                  OS age day_completion month_completion lymph_count T  M
# TCGA-CV-7102-01A  1  76              6               10          24  1 0
# TCGA-MT-A67D-01A  0  55             21                5          2 2 0 0
# TCGA-P-A6T4-01A  1  54             21               10          14 4 1 0
# TCGA-CV-694-01A  1  66             1                9          19  2 0
# TCGA-D6-6824-01A  0  61              7                9          25 4 1 0
# TCGA-C-A642-01A  1  57             17                6         110  2 0
dat[,-1] <- data.frame(lapply(dat[,-1], ))
dat$age <- ifelse(dat$age > 60, 1,0)
dat$day_completion <- ifelse(dat$day_completion > 15, 1,0)
dat$month_completion <- ifelse(dat$month_completion > 6, 1,0)
dat$lymph_count <- ifelse(dat$lymph_count > 20, 1,0)
dat$`T` <- ifelse(dat$`T` > 2, 1,0)
dat$`` <- ifelse(dat$`` > 1, 1,0)
dat$`M` <- ifelse(dat$`M` > 0, 1,0)
str(dat)

# 使用randomForest/ranger去构建随机森林树模型
library(randomForest)
model <- randomForest(OS ~.,data =dat)

library(ranger)
model <- ranger(OS ~  age + day_completion + month_completion + lymph_count + `T` + `` +`M`, 
                data = dat,classification = TRUE,
                probability = TRUE) # 设置输出概率
model
# Ranger result
# 
# Call:
#  ranger(OS ~ age + day_completion + month_completion + lymph_count +      T +  + M, data = dat, classification = TRUE, probability = TRUE) 
# 
# Type:                             Probability estimation 
# umber of trees:                  500 
# Sample size:                      69 
# umber of independent variables:  7 
# Mtry:                             2 
# Target node size:                 10 
# Variable importance mode:         none 
# Splitrule:                        gini 
# OOB prediction error (Brier s.):  0.2189129 

.使用DALEX建立解释器

代码语言：javascript代码运行次数：0运行复制

explainer <- explain(model, 
                     data = dat[,-1], 
                     y = dat$OS==1,
                     label = "Random forest",
                     colorize = FALSE)
# Preparation of a new explainer is initiated
#   -> model label       :  Random forest 
#   -> data              :  69  rows  7  cols 
#   -> target variable   :  69  values 
#   -> predict function  :  yhat.ranger  will be used (  default  )
#   -> predicted values  :  o value for predict function target column. (  default  )
#   -> model_info        :  package ranger , ver. 0.16.0 , task classification (  default  ) 
#   -> model_info        :  Model info detected classification task but 'y' is a logical . Converted to numeric.  (  OTE  )
#   -> predicted values  :  numerical, min =  0.061755 , mean =  0.45077 , max =  0.882942  
#   -> residual function :  difference between y and yhat (  default  )
#   -> residuals         :  numerical, min =  -0.56149 , mean =  -0.0017059 , max =  0.6008169  
#   A new explainer has been created!   

4.shap分析

代码语言：javascript代码运行次数：0运行复制

# 创建一个空列表来存储每个样本的 SHAP 值
shap_values_list <- list()

# 循环遍历每一行数据并计算 SHAP 值
for (i in 1:nrow(dat)) {
  # 提取当前样本的特征值
  new_observation <- dat[i, -1] # 去掉目标列
  names(new_observation) <- colnames(dat)[-1] # 确保列名一致
  
  # 计算 SHAP 值并存储
  shap_values_list[[i]] <- predict_parts(
    explainer = explainer,
    new_observation = new_observation,
    type = "shap",
    B = 25 # 选择多少个排列组合
  )
  
  # 打印当前进度
  cat("Processed observation:", i, "\n")
}

# 查看 SHAP 值列表的结构
str(shap_values_list)

4.1 shap分析-全局(model_parts不支持shap分析!）

代码语言：javascript代码运行次数：0运行复制

shap_total <- model_parts(
    explainer = explainer,
    new_observation = dat[,-1],
    type = "variable_importance", # difference,ratio, raw
    B = 25 # 选择多少个排列组合
  )
shap_total

5.DALEX可视化

代码语言：javascript代码运行次数：0运行复制

plot(shap_values_list[1])

6.对接shapviz可视化

代码语言：javascript代码运行次数：0运行复制

# 初始化一个空列表，用于存储 shapviz 结果
shapviz_res_list <- list()
# 遍历 shap_values_list
for (i in seq_along(shap_values_list)) {
  # 对第 i 个 SHAP 值计算 shapviz
  shapviz_res_list[[i]] <- shapviz(shap_values_list[[i]])
  # 打印进度
  cat("Processed shapviz for sample:", i, "\n")
}
# 检查生成的 shapviz 对象列表
print(shapviz_res_list)

sv_waterfall(shapviz_res_list[[2]])
sv_force(shapviz_res_list[[2]])
sv_importance(shapviz_res_list[[2]])

SHAP结果中最基础的图形结果解释：描述了SHAP值如何给每个特征分配改变期望模型预测的贡献，这种改变是在考虑该特征的情况下产生的。它们解释了如何从基础值 E[f(z)]（如果我们不知道任何特征时的预测值）到当前输出 f(x) 的转变过程。这个图示展示了单一的排序。然而，当模型是非线性的或者输入特征不是独立的时候，特征添加到"期望事件"？(expectation matters)中的顺序很重要，SHAP值是通过平均所有可能排序的φi值得出的。这意味着SHAP值不仅反映了每个特征对预测输出的平均影响，还考虑了特征之间的相互作用和依赖关系，从而提供了一种更全面和精确的特征重要性度量方式。

以下是shapviz可视化的结果，基础值是0.45，通过一个个参数影响的叠加之后最后的值为0.168

但目前DALEX和shapviz的组合拳方式似乎只能展示单个样本的shap结果~

参考资料：

1. github: /
2. DrWhy github: .html
3. shapviz github：
4. Explaining prediction models and individual predicti with feature contributi. Knowl Inf Syst (2014) 41:647–665
5. A Unified Approach to Interpreting Model Predicti. Advances in eural Information Processing Systems 0 (2017)
6. 医学和生信笔记：

注：若对内容有疑惑或者有发现明确错误的朋友，请联系后台(欢迎交流)。更多内容可关注：生信方舟

- ED -