如何在R的rfsrc竞争风险模型中结合SHAP分析变量交互方向？

问题背景 #

使用randomForestSRC包处理含生存时间与状态的生物信息学数据，构建了包含3种结局的竞争风险模型，已完成手动嵌套交叉验证并得到VIP值，但需要明确变量间交互作用的方向。尝试用SHAP生成可视化图表时，现有代码无法适配竞争风险模型的多结局输出逻辑。

适配竞争风险模型的SHAP实现方案 #

核心是修改预测函数，使其正确提取rfsrc竞争风险模型的多结局预测结果，再结合kernelshap和shapviz完成可视化。

步骤1：加载依赖包 #

library(randomForestSRC)
library(kernelshap)
library(shapviz)

步骤2：构建竞争风险模型（示例用模拟数据） #

用mtcars模拟含3类结局的生存数据，mpg作为生存时间，status作为竞争风险事件状态（0=截尾，1/2/3对应三类结局）：

set.seed(123)
df <- mtcars
# 模拟三类结局的竞争风险状态
df$status <- sample(c(0,1,2,3), nrow(df), replace = TRUE)
# 构建竞争风险模型
cr_mod <- rfsrc(Surv(mpg, status) ~ ., data = df, importance = TRUE, ntree = 200)

步骤3：定义适配竞争风险的预测函数 #

竞争风险模型的预测结果存储在predict()返回的predicted列表中，每个元素对应一类结局的累积风险预测，需将其转换为矩阵格式：

pred_fun_cr <- function(model, new_data) {
  pred_obj <- predict(model, newdata = new_data)
  # 提取每类结局的预测值，返回行=样本、列=结局的矩阵
  do.call(cbind, pred_obj$predicted)
}

# 测试预测函数输出格式
head(pred_fun_cr(cr_mod, df))

步骤4：计算SHAP值 #

背景数据集bg_X建议使用训练集的随机子集（如50个样本），大幅提升计算效率：

# 筛选特征列（排除生存时间和状态列）
feat_cols <- setdiff(colnames(df), c("mpg", "status"))

set.seed(456)
shap_cr <- kernelshap(
  model = cr_mod,
  X = df[feat_cols],
  bg_X = df[sample(nrow(df), 50), feat_cols],
  pred_fun = pred_fun_cr
)

步骤5：SHAP可视化分析 #

给SHAP结果命名对应三类结局，生成变量重要性图和依赖图（依赖图可直观展示变量交互方向）：

# 打包SHAP结果并指定结局名称
sv_cr <- shapviz(shap_cr)
names(sv_cr) <- c("结局1", "结局2", "结局3")

# 1. 分结局的SHAP变量重要性图
sv_importance(sv_cr, facet_ncol = 1)

# 2. 单个变量的SHAP依赖图（展示交互方向）
# 示例：查看变量cyl对结局1的影响，用disp作为交互参考变量
sv_dependence(sv_cr, "cyl", outcome = "结局1", color_var = "disp")

关键说明 #

• 竞争风险模型的预测输出是多结局累积风险，预测函数必须返回矩阵格式，每列对应一个结局
• 背景集bg_X避免使用全量数据，随机采样可显著降低计算耗时
• sv_dependence的color_var参数可直观展示变量交互：若目标变量的SHAP值随参考变量的变化趋势明显，趋势方向即为交互作用的方向

内容的提问来源于stack exchange，提问作者陈德赛