R语言：随机生存森林和SHAP分析

随机生存森林算法作为随机森林的拓展算法，主要应用于生存资料的分析。SHAP作为外部解释算法之一，被广泛应用于机器学习模型的可解释性分析中。现尝试在R语言中，采用SHAP算法对训练好的随机生存森林模型进行可解释性分析。

首先是训练随机生存森林模型。

library(randomForestSRC)
data("veteran")
rfsrc_fit <- rfsrc(
  Surv(time,status)~., #公式
  ntree = 100,         # 树的棵树 
  nsplit = 5,          # 最小节点数
  importance = TRUE,   #变量重要性
  tree.err=TRUE,       # 误差
  data=veteran         # 数据集
)
plot(rfsrc_fit)#绘制误差曲线和特征重要性

然后，使用SHAP算法对模型进行可解释性分析。

但是，在实际应用时发现，在R语言中没有随机生存森林对应的SHAP算法包（有XGBoost对应的包），因此我们首先使用XGBoost模型生成一个储存SHAP特征值的“躯壳”，随后计算随机生存森林模型中各个特征的SHAP值，最后使用XGBoost对应的SHAP分析包，生成我们需要的SHAP瀑布图或蜂窝图等。

首先利用XGBoost生成“躯壳”。

X_train <- data.matrix(iris[, -1])
dtrain <- xgboost::xgb.DMatrix(X_train, label = iris[, 1], nthread = 1)
fit <- xgboost::xgb.train(data = dtrain, nrounds = 10, nthread = 1)
x <- shapviz(fit, X_pred = X_train)#数据“壳子”
#class(x)
#"shapviz"

通过class对x进行解析，发现x的格式是shapviz。在x中存在S、X、baseline和S_inter四个子变量。X：变量值，S：变量对应的SHAP值，baseline：基准值，S_inter：并没有发现有具体的作用。

生成躯壳之后，计算随机生存森林模型中不同特征对应的SHAP值。计算SHAP值可能会需要较多的时间。

#####计算shap值
#定义预测函数
p_function <- function(model, data)
{
  preds <- predict(model, newdata = data,type = "response")
  return(preds$predicted)  
}
#构建解释器
rf_exp <- DALEX::explain(rfsrc_fit,
                         data = veteran[,-3:-4],
                         y=veteran$status==1,
                         predict_function = p_function)
#计算变量的shap值，可以按照自己的需要设置数据条目数量
n <- nrow(veteran)#调整n值，设置自己需要的数据量
shapvalue <- as.data.frame(matrix(nrow=6,ncol=n+1))#生成存储shap值的空数据框
for(i in 1:n){
  shap_rf <- shap(rf_exp,veteran[i,])
  #对数据进行整合
  write.csv(shap_rf,'1234.csv',row.names = F)
  shapvalue1 <- read.csv('1234.csv')
  shapvalue[,1] <- shapvalue1[1:6,3]
  shapvalue[,i+1] <-shapvalue1[1:6,2]
  print(i)
}
#对数据的格式进行调整（这部分可以简化）
write.csv(t(shapvalue),'shapvalue1.csv',row.names = F)
shapvalue1 <- read.csv('shapvalue1.csv',header = T)
colnames(shapvalue1) <- shapvalue1[1,]#设置列名
shapvalue1 <- shapvalue1[-1,]#将第一行数据删除
shapvalue1[,1]
shapvalue1=as.data.frame(lapply(shapvalue1,as.numeric))#将数据调整为num类型

随后，按照驱壳中数据要求，将随机生存森林中的数据输入。经过网络检索发现基线值一般是模型预测值的均值。

x$S <- as.matrix(shapvalue1)#将计算出的shap值导入x$S
x$X <- veteran[1:nrow(x$S),-3:-4]#将shap值对应的变量值导入x$X
x$baseline <- mean(rfsrc_fit$predicted[1:nrow(x$S)])

最后就可以进行可视化。

在瀑布图中我们以第20个样本数据为例进行可视化展示。

#shap可视化
sv_importance(x)
sv_importance(x, kind = "no")
sv_importance(x, kind = "beeswarm")

row <- 20
sv_waterfall(x,row_id = row,fill_colors = c("#f7d13d", "#a52c60"))
sv_force(x,row_id = row)

总结一下，在R语言中没有随机生存森林对应的SHAP分析包，我们通过使用XGBoost模型的SHAP分析包，借助“躯壳”，实现对随机生存森林SHAP分析的可视化。