使用 R 包“survival”和“rms”校准 Cox PH 模型：时间单位混淆

如何解决使用 R 包“survival”和“rms”校准 Cox PH 模型：时间单位混淆

我使用 R 包“rms”构建了 Cox Proportional Hazards 模型，并尝试对其进行交叉验证。将数据分成训练集和测试集是我想做的，但我是生存分析的新手，除了 rms::calibrate 之外，在文献中找不到任何东西。我无法让它工作。

代码如下：

# using the package 'survival',I make a survival object with
# follow-up time (2000 to 2020) and status (event=1,survival/censoring 0)
surv2 <- Surv(grid3@data$def_mean,grid3@data$status)
d <- datadist(grid3@data) # stores distribution summaries for potential variables??
options(datadist = "d") # seems to help cph() refer to variables
model <- cph(surv2 ~ cost_mean + elev_mean + popn_mean + cop99 + PAs_mean,data = grid3@data,x = TRUE,y = TRUE,surv = TRUE,time.inc=1)
modrms <- rms::calibrate(model,B = 40,u = 1)

'time.inc' 是时间增量（1 年） - 查看 model$surv.summary，我可以看到生存和“不”。 20 年中每一年的风险数据。所以这是有道理的。但是调用 rms::calibrate 我得到的第一条消息是 Using Cox survival estimates at 1 Days...并查看我得到的校准：

> summary(attr(modrms,"predicted"))
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
      1       1       1       1       1       1

...所以看起来模型已经校准了 1 天？当然，每个人都活了下来 (1=100%) ...使用 rms::calibrate(model,u = 20) 也会发生同样的事情。

我再次尝试开始：

units(grid3@data$def_mean) <- "year"
surv3 <- Surv(grid3@data$def_mean,grid3@data$status)

...但这给了我一个错误！

Error in Ops.units(time,origin) : 
  both operands of the expression should be "units" objects

我不知道接下来该尝试什么。如果我能用 2000-10 年的数据建立一个模型，用它来预测 2010-20 年，并查看预测与实际，那不是很好吗？但我一直坚持校准，而且文档假设了比我更多的统计专业知识（大学统计数据加上努力提高我的数学）。

这是数据结构（不知道如何使这个可重现）：

> str(grid3@data)
'data.frame':   36918 obs. of  7 variables:
 $ def_mean  : num  20 20 20 20 20 20 20 20 20 20 ...
 $ status    : num  0 0 0 0 0 0 0 0 0 0 ...
 $ elev_mean : num  -0.664 -0.664 -0.664 -0.664 -0.664 ...
 $ popn_mean : num  -0.1658 0.0664 -0.1484 0.0601 -0.0381 ...
 $ cost_mean : num  1.53 1.48 1.43 1.66 1.6 ...
 $ PAs_mean  : num  0 0 0 0 0 0 0 0 0 0 ...
 $ cop99     : Factor w/ 12 levels "10","20","30",..: 5 5 5 5 5 5 5 5 5 5 ...

> summary(grid3@data$def_mean)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   3.20   20.00   20.00   19.59   20.00   20.00

> table(grid3@data$status)

    0     1 
34696  2222

解决方法

rms 软件包在某些方面可能具有陡峭的学习曲线。您已经了解到有时会被忽略的一件事：使用 datadist() 来总结预测变量的重要性，然后设置 datadist 选项（使用 datadist对象），以便汇总函数具有合理的默认显示选择。

关于第二个错误，我想知道您是否在更改时间单位后没有重新运行datadist()命令并重置datadist选项。 units() 相关的 label() 包中的 rms 和 Hmisc 函数可能非常有用，但如果您在之后不重新运行和重置 datadist()使用它们我怀疑事情会让下游的软件感到困惑。如果你在一个地方指定一个单位，它可能会在另一个地方期待相同的单位。

然而，这些命令不做任何转换。默认假设是时间单位是“天”，因此这就是默认情况下在输出中打印的内容。如果您将“单位”更改为“年”，打印输出将显示“年”而不是“日”，但基础计算不会改变。

因此，尽管 calibrate() 最初声称以“1 天”计算，但事实并非如此；那只是它的默认打印单位。它仍然在 time = 1 处进行校准。这么早的校准可能不是您想要的。

我依稀记得如果原始 time.inc 调用中的 cph() 设置与 u 调用中的 calibrate() 设置不匹配，会遇到一些问题。我通常的做法是知道我想要校准的时间点（例如，某些类型癌症数据的 3 年生存率）并将其用于这两种设置。玩一下玩具数据集，看看如何让它适合你。

最后，calibrate() 最好与 plot() 一起使用来显示校准曲线（理想的、建模的、由引导程序校正的乐观）。如果您尝试 print() calibrate 对象，可能会出现故障。标准图上显示的值是正确的。