Cox版本的Boruta+SHAP分析（心力衰竭数据集）
Boruta算法是变量筛选的有力工具，而SHAP分析是观察预测变量与结局变量间关系的不错的方法，在传统的分析方法的基础上提供了一个全新的视角。Boruta算法+SHAP分析，正在逐渐成为一种流行的分析策略。
COX分析是医学分析中最重要的一种类型，分类结局中纳入了时间因素，较单纯的二分类结局能够提供更多的信息，比如某变量对于结局的贡献随时间是如何变化的？Boruta+SHAP的分析策略也可以用在Cox分析中，并且展示了一些cox分析独特的特点。
1. Boruta算法进行变量筛选代码
这部分还是使用Boruta包，如图中仅作稍许变化即可。
library(Boruta)
library(survival)
set.seed(1)
Y = Surv(data$time, data$DEATH_EVENT)
boruta_obj<-Boruta(y=Y,x=data[,c(1:11)],doTrace=0,ntree=500,pValue=0.01)
print(TentativeRoughFix(boruta_obj))#分两类
print(boruta_obj)#三类，包含怀疑的数据
2. SHAP分析代码
参考treeshap包的代码，一键运行可以到和鲸社区相关项目查看。
library(tidyverse)
library(treeshap)
#先转成data.table格式
data <- data.table::data.table(data[,c('age','ejection_fraction','serum_creatinine','serum_sodium','time','DEATH_EVENT')])
surv_cols <- c("DEATH_EVENT", "time")
feature_cols <- colnames(data)
#这里比较少见的函数的操作，属于data.table内容
train_x <- model.matrix(
~ -1 + .,
data[, .SD, .SDcols = setdiff(feature_cols, surv_cols[1:2])]
)
train_y <- survival::Surv(
event = (data[, get("DEATH_EVENT")] %>%
as.character() %>%
as.integer()),
time = data[, get("time")],
type = "right"
)
rf <- ranger::ranger(
x = train_x,
y = train_y,
data = data,
max.depth = 10,
num.trees = 10
)
unified_model_risk <- ranger_surv.unify(rf, train_x, type = "risk")#"risk"type计算累积风险的sum，
shaps <- treeshap(unified_model_risk, train_x)
# compute shaps for 3 selected months of followup points
unified_model_surv <- ranger_surv.unify(rf, train_x, type = "survival", time = c(36, 60, 120))
shaps_surv <- treeshap(unified_model_surv, train_x)
3. 心力衰竭数据集分析的部分结果
- 四个变量对于心力衰竭的结局有所贡献，最重要的是ejection_fraction，接下来我们单独考查这一个指标。
- ejection_fraction与心力衰竭的关系，30以下是有风险的，越低风险越大。
- 年龄与心力衰竭的关系
分析的结果展现了一些心力衰竭相关因素的变化，还可以指定时间点，对相关的因素进行更加细致的考察。但是，因为使用的数据集是公开的数据集，我们并不确定来源的可靠性，所以并不确定所展示相关关系的准确性，而只做技术上的讨论。
