使用 R 中“rpart”包中的生存树来预测新的观察结果

Posted 2023-02-19

【中文标题】使用 R 中“rpart”包中的生存树来预测新的观察结果

【英文标题】：Using a survival tree from the 'rpart' package in R to predict new observations

【发布时间】：2015-08-22 10:00:06

【问题描述】：

我正在尝试使用 R 中的“rpart”包来构建生存树，我希望使用这棵树来预测其他观察结果。

我知道有很多关于 rpart 和预测的 SO 问题；但是，我找不到任何解决（我认为）特定于将 rpart 与“Surv”对象一起使用的问题。

我的特殊问题涉及解释“预测”功能的结果。一个例子很有帮助：

library(rpart)
library(OIsurv)

# Make Data:
set.seed(4)
dat = data.frame(X1 = sample(x = c(1,2,3,4,5), size = 1000, replace=T))
dat$t = rexp(1000, rate=dat$X1)
dat$t = dat$t / max(dat$t)
dat$e = rbinom(n = 1000, size = 1, prob = 1-dat$t )

# Survival Fit:
sfit = survfit(Surv(t, event = e) ~ 1, data=dat)
plot(sfit)

# Tree Fit:
tfit = rpart(formula = Surv(t, event = e) ~ X1 , data = dat, control=rpart.control(minsplit=30, cp=0.01))
plot(tfit); text(tfit)

# Survival Fit, Broken by Node in Tree:
dat$node = as.factor(tfit$where)
plot( survfit(Surv(dat$t, event = dat$e)~dat$node) )

到目前为止一切顺利。我对这里发生的事情的理解是，rpart 正试图将指数生存曲线拟合到我的数据子集。基于这种理解，我相信当我调用predict(tfit) 时，对于每个观察，我都会得到一个与该观察的指数曲线参数相对应的数字。因此，例如，如果 predict(fit)[1] 是 0.46，那么这意味着对于我的原始数据集中的第一次观察，曲线由公式 P(s) = exp(−λt) 给出，其中 λ=.46。

这似乎正是我想要的。对于每个观察（或任何新观察），我可以获得该观察在给定时间点内存活/死亡的预测概率。 （编辑：我意识到这可能是一个误解——这些曲线没有给出生/死的概率，而是给出了在一个区间内存活的概率。不过，这不会改变下面描述的问题。）

但是，当我尝试使用指数公式时...

# Predict:
# an attempt to use the rates extracted from the tree to
# capture the survival curve formula in each tree node.
rates = unique(predict(tfit))
for (rate in rates) 
  grid= seq(0,1,length.out = 100)
  lines(x= grid, y= exp(-rate*(grid)), col=2)

我在这里所做的是以与生存树相同的方式拆分数据集，然后使用survfit 为每个分区绘制非参数曲线。那是黑线。我还绘制了与将（我认为是）“速率”参数插入（我认为是）生存指数公式的结果相对应的线条。

我知道非参数拟合和参数拟合不一定相同，但这似乎不止于此：似乎我需要缩放我的 X 变量或其他东西。

基本上，我似乎不明白 rpart/survival 在幕后使用的公式。谁能帮我从 (1) rpart 模型到 (2) 任意观察的生存方程？

【参考方案1】：

生存数据在内部以指数方式缩放，因此根节点中的预测速率始终固定为1.000。然后，predict() 方法报告的预测总是与根节点中的生存相关，即高或低某个因素。有关详细信息，请参阅vignette("longintro", package = "rpart") 中的第 8.4 节。无论如何，您报告的 Kaplan-Meier 曲线与 rpart 小插图中报告的完全一致。

如果您想直接获得树中的 Kaplan-Meier 曲线图并获得预测的中位生存时间，您可以将 rpart 树强制转换为 constparty 树，由 partykit 包提供：

library("partykit")
(tfit2 <- as.party(tfit))
## Model formula:
## Surv(t, event = e) ~ X1
## 
## Fitted party:
## [1] root
## |   [2] X1 < 2.5
## |   |   [3] X1 < 1.5: 0.192 (n = 213)
## |   |   [4] X1 >= 1.5: 0.082 (n = 213)
## |   [5] X1 >= 2.5: 0.037 (n = 574)
## 
## Number of inner nodes:    2
## Number of terminal nodes: 3
##
plot(tfit2)

打印输出显示中位生存时间和相应的 Kaplan-Meier 曲线的可视化。也可以通过predict() 方法将type 参数分别设置为"response" 和"prob" 来获得。

predict(tfit2, type = "response")[1]
##          5 
## 0.03671885 
predict(tfit2, type = "prob")[[1]]
## Call: survfit(formula = y ~ 1, weights = w, subset = w > 0)
## 
##  records    n.max  n.start   events   median  0.95LCL  0.95UCL 
## 574.0000 574.0000 574.0000 542.0000   0.0367   0.0323   0.0408 

作为rpart 生存树的替代方案，您还可以考虑基于ctree() 中的条件推理的非参数生存树（使用对数秩分数）或使用通用mob() 基础设施的完全参数生存树partykit 包。

【讨论】：

感谢详细的回复！不过，我的目标是在任何时间点为任何实例获取 P(alive)。这似乎应该给我更多的信息，而不仅仅是提取与每个实例的树节点相关的中位生存时间。我能够做到这一点的唯一方法是使用“pec”包中的 predictSurvProb 函数，但这个函数有点问题，我也希望从生存曲线计算生存概率会更有效他们自己，而不是依赖这个函数。

是的，Kaplan-Meier 函数是幸存者函数 S(t) 的（非参数）估计量，即在某个时间仍然活着的概率吨。 Kaplan-Meier 函数既可以使用survfit() 手动计算，也可以像您一样使用基于$where 的因子 - 或通过带有type = "prob" 的partykit。如果您想在每个叶子中拟合参数模型（例如，指数或 Weibull），您可以使用 survreg() 而不是 survfit()。

对不起，我没有完全关注：您能否编辑您的帖子以提供实际代码，为给定的 t 和给定的实例提供 S(t)？例如，给定一个 rpart 对象 tfit 和一个实例 dat[1,]，以及一个时间 dat[1,'t']，我应该使用什么代码来获取那个实例和那个 t 的 S(t)？

我不明白您为什么要编辑我的答案。上面显示的代码 sn-p predict(tfit2, type = "prob")[[1]] 提取拟合的 survfit 对象用于第一次观察。从中您可以提取您喜欢的所有“通常”数量。例如，查看对象的summary()，它显示了完整的 Kaplan-Meier 曲线坐标以及一些附加信息。

但这确实是一个关于survfit 和survival 的问题，其中还有有用的书籍、教程等。但我认为如果你这样做：km1 <- predict(tfit2, type = "prob")[[1]] 然后@987654353 @你应该看到你需要的一切。你可以很容易地从中得到分位数，例如quantile(km1, c(0.2, 0.5, 0.8))，它给出了 S(t) 分别为 0.8、0.5 和 0.2 的时间。或者，如果您想要一个功能，您可以执行 km1f <- approxfun(km1$time, km1$surv) 然后 km1f(c(0.011, 0.037, 0.094)) 等。

【参考方案2】：

@Achim Zeileis 的回答很有帮助，但似乎没有回答确切的@jwdink 问题。我将其理解为“如果 RPart 树按照最佳指数生存拟合进行分裂，那么这些拟合的绝对值的 Lambda 是多少，因此我们可以使用这些指数生存函数进行预测”。 RPart 摘要确实显示了估计的比率，但只是在假设整个人口的比率为 1 的相对条件下。要克服，可以拟合指数 survreg，从那里获取参考的 lambda，然后将 RPart 预测比率乘以该数字（见下面的代码）。

也就是说，这不是如何从树中预测 RPart 中的存活率。我没有直接在 RPart 中找到生存预测函数，但是正如 Achim 上面指出的那样，partykit 使用 Kaplan-Meier 估计，即来自最终叶的非参数生存。我认为在生存随机森林树中也是如此，在最终的叶子中使用了 K-M 曲线。

此问题中的模拟数据使用指数分布，因此 KM 和指数生存曲线在设计上是相似的，但是对于不同的模拟或实际分布，通过 RPart 树估计指数率并在最终叶子中使用 KM 曲线（同一棵树的）会给出不同的存活率。

sfit = survfit(Surv(t, event = e) ~ 1, data=dat)
tfit = rpart(formula = Surv(t, event = e) ~ X1 , data = dat, control=rpart.control(minsplit=30, cp=0.01))
plot(tfit); text(tfit)

# Survival Fit, Broken by Node in Tree:
dat$node = as.factor(tfit$where)
table(dat$node)
s0 = survreg(Surv(t,e)~ 1, data =  dat, dist = "exponential") #-0.6175
e0 = exp(-summary(s0)$coefficients[1]); e0 #1.854
rates = unique(predict(tfit))
#1) plot K-M curves by node (black):
plot( survfit(Surv(dat$t, event = dat$e)~dat$node) )

#2) plot exponential survival with rates = e0 * RPart rates (red):
for (rate in rates) 
  grid= seq(0,1,length.out = 100)
  lines(x= grid, y= exp(-e0*rate*(grid)), col=2)

#3) plot partykit survival curves based on RPart tree (green)
library(partykit)
tfit2 <- as.party(tfit)
col_n = 1
for (node in names(table(dat$node)))
  predict_curve = predict(tfit2, newdata = dat[dat$node == node, ], type = "prob")  
  surv_esitmated = approxfun(predict_curve[[1]]$time, predict_curve[[1]]$surv)
  lines(x= grid, y= surv_esitmated(grid), col = 2+col_n)
  col_n=+1