如何在随机森林SRC中计算中位生存率

Question

我试图计算随机森林SRC的中位生存率。

library(randomForestSRC)
data(veteran, package = "randomForestSRC")
train <- sample(1:nrow(veteran), round(nrow(veteran) * 0.80))
veteran.grow <- rfsrc(Surv(time, status) ~ ., veteran[train, ], ntree = 100)
veteran.pred <- predict(veteran.grow, veteran[-train , ])
print(veteran.grow)
print(veteran.pred)

我得到了生存功能

veteran.pred$survival

我想得到生存函数的中值(值< code >老兵. pred $生存== 0.5，例如，第一行)

a = as.data.frame(veteran.pred$time.interest)
b= as.data.frame(veteran.pred$survival[1,])
df =cbind(a, b)
df <-rename(df, `time` =`veteran.pred$time.interest`)
df$`veteran.pred$survival[1, ]` = round(df$`veteran.pred$survival[15, ]`, 2)
subset(df, df$`veteran.pred$survival[1, ]`== 0.5 )$time

问题是，生存函数很少有意义。在我们的例子中，精确到0.5

[1] 1.00 1.00 1.00 0.95 0.95 0.95 0.95 0.93 0.93 0.93 0.93 0.93 0.93 0.93 0.90 0.90 0.90 0.87 0.87 0.87
[21] 0.87 0.87 0.86 0.76 0.76 0.67 0.67 0.67 0.67 0.67 0.67 0.66 0.66 0.66 0.66 0.66 0.65 0.64 0.57 0.57
[41] 0.57 0.57 0.57 0.57 0.57 0.48 0.45 0.39 0.39 0.32 0.32 0.29 0.29 0.29 0.29 0.26 0.26 0.25 0.25 0.25
[61] 0.25 0.25 0.25 0.25 0.20 0.20 0.20 0.20 0.20 0.20 0.18 0.18 0.18 0.16 0.16 0.16 0.14 0.14 0.14 0.11
[81] 0.07 0.07

所以我们有最接近的值0.57

Answer 1

我不认为< code>time.interest是用于该目的的正确项目，因为它的长度比验证集中的案例数要长得多。(也就是说，我不确定它到底想告诉你什么。)如果您查看< code>str(veteran.pred)的输出，您会在顶部看到:

str(veteran.pred)
List of 32
 $ call          : language generic.predict.rfsrc(object = object, newdata = newdata, outcome.target = outcome.target,      importance = impo| __truncated__ ...
 $ family        : chr "surv"
 $ n             : int 27
 $ ntree         : num 100
 $ yvar          :'data.frame': 27 obs. of  2 variables:
  ..$ time  : int [1:27] 100 384 123 22 21 139 31 51 54 132 ...
  ..$ status: int [1:27] 0 1 0 1 1 1 1 1 1 1 ...
 $ yvar.names    : chr [1:2] "time" "status"
 $ xvar          :'data.frame': 27 obs. of  6 variables:
  ..$ trt     : int [1:27] 1 1 1 1 1 1 1 1 1 1 ...
  ..$ celltype: int [1:27] 1 2 2 2 2 2 2 2 2 3 ...
  ..$ karno   : int [1:27] 70 60 40 60 40 80 75 60 70 80 ...
  ..$ diagtime: int [1:27] 6 9 3 4 2 2 3 1 1 5 ...
  ..$ age     : int [1:27] 70 42 55 68 55 64 65 67 67 50 ...
  ..$ prior   : int [1:27] 0 0 0 0 10 0 0 0 0 0 ...
 $ xvar.names    : chr [1:6] "trt" "celltype" "karno" "diagtime" ...
 # --- snipped

我认为，由于27是< code >老兵[-train，]的行数，您需要使用预测列表中的< code>yvar项:

 str(veteran.pred$yvar)
#'data.frame':  27 obs. of  2 variables:
# $ time  : int  100 384 123 22 21 139 31 51 54 132 ...
# $ status: int  0 1 0 1 1 1 1 1 1 1 ...
?survfit
 survfit(Surv(time,status)~1  , data=veteran.pred$yvar)
#Call: survfit(formula = Surv(time, status) ~ 1, data = veteran.pred$yvar)
#
      n  events  median 0.95LCL 0.95UCL 
     27      24      54      49     139 
plot( survfit(Surv(time,status)~1  , data=veteran.pred$yvar) )

我对建议这一程序持严重保留意见。请注意，在该列表的林节点内有另一个名为 yvar 的项目，它有 110 行（因此它是原始数据）。IOf你看看传统KM曲线在未经调整的镇痛上的结果，你会得到：

survfit(Surv(time,status)~1  , data=veteran.pred$forest$yvar) 
Call: survfit(formula = Surv(time, status) ~ 1, data = veteran.pred$forest$yvar)

      n  events  median 0.95LCL 0.95UCL 
    110     104      87      53     111 

我认为来自80:20 CV策略的27个案例将给出一个非常不稳定的中位数估计方法，特别是当存在分类预测时。我还认为，随机forrest范式应该能够从整个数据集中得出有用的预测，而无需使用交叉验证类型的数据拆分。