如何在R中为randomForestSRC多变量随机森林模型绘制SHAP值？

一、基于现有randomForestSRC模型计算SHAP值 #

randomForestSRC本身没有内置SHAP计算工具，但可以借助通用SHAP包，针对每个因变量单独处理：

1. 使用fastshap包（快速蒙特卡洛近似） #

通过自定义预测函数，为每个目标变量单独计算SHAP值：

library(randomForestSRC)
library(fastshap)
library(ggplot2)
library(reshape2)

# 构建多因变量模型
mod1 <- rfsrc(Multivar(mpg, wt) ~ ., data = mtcars, importance = TRUE)

# 定义针对mpg的预测函数
predict_mpg <- function(object, newdata) {
  predict(object, newdata = newdata, outcome.target = "mpg")$predicted
}

# 计算mpg的SHAP值
shap_mpg <- fastshap::explain(mod1, X = mtcars[, -c(1, 6)], 
                              pred_wrapper = predict_mpg, nsim = 10)

# 绘制mpg的SHAP汇总图
ggplot(melt(shap_mpg), aes(x = value, y = Var1)) +
  geom_boxplot(fill = "#377EB8", alpha = 0.7) +
  labs(title = "SHAP Summary Plot for mpg", x = "SHAP Value", y = "Feature") +
  theme_minimal()

# 同理处理wt的SHAP值
predict_wt <- function(object, newdata) {
  predict(object, newdata = newdata, outcome.target = "wt")$predicted
}

shap_wt <- fastshap::explain(mod1, X = mtcars[, -c(1, 6)], 
                             pred_wrapper = predict_wt, nsim = 10)

ggplot(melt(shap_wt), aes(x = value, y = Var1)) +
  geom_boxplot(fill = "#E41A1C", alpha = 0.7) +
  labs(title = "SHAP Summary Plot for wt", x = "SHAP Value", y = "Feature") +
  theme_minimal()

2. 使用shapr包（精确SHAP值计算） #

适合小数据集，计算精确SHAP值：

library(shapr)
library(reshape2)
library(ggplot2)

# 定义返回所有因变量预测值的函数
predict_multivar <- function(object, newdata) {
  predict(object, newdata = newdata)$predicted
}

# 初始化SHAP解释器
explainer <- shapr(mtcars[, -c(1, 6)], mod1)

# 计算所有因变量的SHAP值
shap_values <- explain(explainer, X = mtcars[, -c(1, 6)], 
                       prediction_zero = colMeans(mtcars[, c(1, 6)]),
                       predict_model = predict_multivar)

# 提取mpg的SHAP值并绘图
shap_mpg_df <- shap_values$shap_values[, , "mpg"]
ggplot(melt(shap_mpg_df), aes(x = value, y = Var1)) +
  geom_boxplot(fill = "#377EB8", alpha = 0.7) +
  labs(title = "SHAP Summary Plot for mpg", x = "SHAP Value", y = "Feature") +
  theme_minimal()

# 提取wt的SHAP值并绘图
shap_wt_df <- shap_values$shap_values[, , "wt"]
ggplot(melt(shap_wt_df), aes(x = value, y = Var1)) +
  geom_boxplot(fill = "#E41A1C", alpha = 0.7) +
  labs(title = "SHAP Summary Plot for wt", x = "SHAP Value", y = "Feature") +
  theme_minimal()

二、换用支持多因变量+SHAP的替代包 #

1. mlr3生态（一站式多输出工作流） #

借助mlr3multioutput和mlr3shap实现多因变量模型的SHAP可视化：

library(mlr3)
library(mlr3multioutput)
library(mlr3learners)
library(mlr3shap)

# 创建多输出回归任务
task <- TaskMultioutput$new(id = "mtcars", backend = mtcars, target = c("mpg", "wt"))

# 初始化多输出随机森林学习器
learner <- lrn("regr.ranger_multioutput")

# 训练模型
learner$train(task)

# 计算并绘制mpg的SHAP汇总图
shap_result <- learner$predict(task)$explain(task, method = "shap")
shap_result$filter(target = "mpg")$plot(type = "summary") +
  ggtitle("SHAP Summary Plot for mpg")

# 计算并绘制wt的SHAP汇总图
shap_result$filter(target = "wt")$plot(type = "summary") +
  ggtitle("SHAP Summary Plot for wt")

2. xgboost多输出梯度提升树（替代随机森林） #

如果可接受梯度提升树，xgboost支持多输出任务，结合Python的shap库（通过reticulate调用）实现可视化：

library(xgboost)
library(reticulate)
shap <- import("shap")

# 准备数据
X <- as.matrix(mtcars[, -c(1, 6)])
y <- as.matrix(mtcars[, c(1, 6)])

# 构建多输出xgboost模型
xgb_mod <- xgboost(data = X, label = y, nrounds = 100, 
                   objective = "reg:squarederror")

# 计算SHAP值
explainer <- shap$TreeExplainer(xgb_mod)
shap_values <- explainer$shap_values(X)

# 绘制mpg的SHAP汇总图
shap$summary_plot(shap_values[[1]], X, feature_names = colnames(X))

# 绘制wt的SHAP汇总图
shap$summary_plot(shap_values[[2]], X, feature_names = colnames(X))

关键注意事项 #

• 多因变量模型的SHAP值需按每个目标变量单独计算和可视化，因每个因变量的特征贡献逻辑独立。
• fastshap的nsim参数控制模拟次数，值越大结果越精确但计算越慢；shapr适合小数据集的精确计算。

内容的提问来源于stack exchange，提问作者Silvia