如何计算 randomForestSRC 中的中位生存期

How to calculate median survival in randomForestSRC

我尝试计算 randomForestSRC 中的中位生存期。

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

我想获得生存函数的中位数(veteran.pred$survival== 0.5 的值 veteran.pred$time.interest),例如,对于第一行

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 & 0.48 如何计算中位生存期?

我认为 time.interest 不是用于该目的的正确项目,因为它的长度比验证集中的案例数长得多。 (也就是说,我不确定它到底应该告诉你什么。)如果你看一下 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是veteran[-train , ]的行数,所以您需要使用预测列表中的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 行(因此它是原始数据)。如果你在未经调整的分析上查看传统 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 个案例将为中位数提供非常不稳定的估计方法,尤其是在存在分类预测变量时。我还认为随机福雷斯特范式应该能够从整个数据集中得出有用的预测,而无需使用数据的交叉验证类型拆分。