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

Question

我正在尝试使用R中的"rpart"包来构建生存树,我希望使用这棵树然后对其他观察进行预测.

我知道有很多涉及rpart和预测的SO问题; 但是,我无法找到任何解决(我认为)特定于使用带有"Surv"对象的rpart的问题.

我的特殊问题涉及解释"预测"功能的结果.一个例子很有帮助:

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]是.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绘制每个分区的非参数曲线.那是黑线.我还绘制了与插入(我认为是)'rate'参数(我认为是)生存指数公式的结果对应的行.

我知道非参数和参数拟合不一定相同,但这似乎不止于此:似乎我需要扩展我的X变量或其他东西.

基本上,我似乎并不理解rpart/survival在引擎盖下使用的公式.任何人都可以帮助我从(1)rpart模型得到(2)任意观察的生存方程式吗？

Answer 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()基础设施,根据(使用logrank得分)或完全参数生存树中的条件推断来考虑非参数生存树partykit.