使用 R 包 'survival' 和 'rms' 校准 Cox PH 模型：时间单位混淆
Calibrating a Cox PH model with R packages 'survival' and 'rms': time unit confusion
我用 R 包“rms”构建了一个 Cox 比例风险模型，并试图对其进行交叉验证。将数据拆分为训练集和测试集是我想做的，但我是生存分析的新手，除了 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 年）- 查看模型 $surv.summary，我可以看到 20 年中每一年的生存率和 'no. at risk' 数字。所以这是有道理的。但是调用 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, B = 40, 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 选项。 rms 相关的 Hmisc 包中的 units() 和 label() 函数可能非常有用，但是如果你不重新 运行 和重置 datadist() 在使用它们之后，我怀疑下游的软件会变得混乱。如果您在一个地方指定一个单位，它可能会期望在另一个地方使用相同的单位。
但是，这些命令不进行任何转换。默认的假设是时间单位是“天”，所以这就是输出中默认打印的内容。如果将“单位”更改为“年”，打印输出将显示“年”而不是“日”，但基础计算不会改变。
因此，虽然 calibrate() 首先声称在“1 天”进行计算，但事实并非如此；那只是它的默认打印单位。它仍然在 time = 1 处进行了校准。这么早的校准可能不是你想要的。
如果原始 cph() 调用中的 time.inc 设置与 calibrate() 调用中的 u 设置不匹配，我依稀记得会遇到一些问题。我通常的做法是知道我想要校准的时间点（例如，某些类型的癌症数据的 3 年生存期）并将其用于这两种设置。尝试一下玩具数据集，看看如何让它为您所用。
最后，calibrate() 最好与 plot() 一起使用以显示校准曲线（理想的、建模的、乐观校正的 bootstrap）。如果您尝试 print() calibrate 对象，可能会出现故障。标准图上显示的值是正确的。
我用 R 包“rms”构建了一个 Cox 比例风险模型，并试图对其进行交叉验证。将数据拆分为训练集和测试集是我想做的，但我是生存分析的新手，除了 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 年）- 查看模型 $surv.summary，我可以看到 20 年中每一年的生存率和 'no. at risk' 数字。所以这是有道理的。但是调用 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, B = 40, 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 选项。 rms 相关的 Hmisc 包中的 units() 和 label() 函数可能非常有用，但是如果你不重新 运行 和重置 datadist() 在使用它们之后，我怀疑下游的软件会变得混乱。如果您在一个地方指定一个单位，它可能会期望在另一个地方使用相同的单位。
但是，这些命令不进行任何转换。默认的假设是时间单位是“天”，所以这就是输出中默认打印的内容。如果将“单位”更改为“年”，打印输出将显示“年”而不是“日”，但基础计算不会改变。
因此，虽然 calibrate() 首先声称在“1 天”进行计算，但事实并非如此；那只是它的默认打印单位。它仍然在 time = 1 处进行了校准。这么早的校准可能不是你想要的。
如果原始 cph() 调用中的 time.inc 设置与 calibrate() 调用中的 u 设置不匹配，我依稀记得会遇到一些问题。我通常的做法是知道我想要校准的时间点（例如，某些类型的癌症数据的 3 年生存期）并将其用于这两种设置。尝试一下玩具数据集，看看如何让它为您所用。
最后，calibrate() 最好与 plot() 一起使用以显示校准曲线（理想的、建模的、乐观校正的 bootstrap）。如果您尝试 print() calibrate 对象，可能会出现故障。标准图上显示的值是正确的。