如何基于R语言带时变协变量的Cox模型生成用户级生存预测？

问题背景 #

我们有一个右删失生存分析场景：从t=0开始按固定时间间隔追踪用户，包含静态协变量X1、X2和时变协变量X3（每个时间间隔测量一次，对应区间删失数据）。

数据准备 #

模拟宽表数据 #

library(survival)
library(MASS)
library(data.table)

dt.wide <-
  data.table("id" = c(1, 2, 3, 4, 5, 6),
             "time" = c(5, 4, 3, 9, 1, 8),
             "status" = c(1, 1, 1, 0, 0, 1),
             "X1" = c(52, 59, 54, 60, 39, 26),
             "X2" = c(1, 1, 1, 0, 0, 1),
             "X3.1" = c(21, 17, 17, 23, 26, 19),
             "X3.2" = c(18, 13, 15, 19, NA, 14),
             "X3.3" = c(16, 22, 17, 17, NA, 20),
             "X3.4" = c(19, 13, NA, 20, NA, 15),
             "X3.5" = c(25, NA, NA, 19, NA, 20),
             "X3.6" = c(NA, NA, NA, 19, NA, 20),
             "X3.7" = c(NA, NA, NA, 16, NA, 22),
             "X3.8" = c(NA, NA, NA, 13, NA, 20),
             "X3.9" = c(NA, NA, NA, 23, NA, NA))

转换为长表格式（适配时变协变量） #

将宽表转换为单用户多行的区间删失长表，满足Cox模型对时变协变量的处理要求：

dt.long <-
  data.table("id"    = c(1, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 4, 4, 4, 4, 4, 4,
                         4, 4, 4, 5, 6, 6, 6, 6, 6, 6, 6, 6),
             "start" = c(0, 1, 2, 3, 4, 0, 1, 2, 3, 0, 1, 2, 0, 1, 2, 3, 4, 5,
                         6, 7, 8, 0, 0, 1, 2, 3, 4, 5, 6, 7),
             "stop"  = c(1, 2, 3, 4, 5, 1, 2, 3, 4, 1, 2, 3, 1, 2, 3, 4, 5, 6,
                         7, 8, 9, 1, 1, 2, 3, 4, 5, 6, 7, 8),
             "event" = c(0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0,
                         0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1),
             "X1"    = c(52, 52, 52, 52, 52, 59, 59, 59, 59, 54, 54, 54, 60,
                         60, 60, 60, 60, 60, 60, 60, 60, 39, 26, 26, 26, 26,
                         26, 26, 26, 26),
             "X2"    = c(1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0,
                         0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1),
             "X3"    = c(21, 18, 16, 19, 25, 17, 13, 22, 13, 17, 15, 17, 23,
                         19, 17, 20, 19, 19, 16, 13, 23, 26, 19, 14, 20, 15,
                         20, 20, 22, 20))

拟合Cox比例风险模型 #

cph.model <- coxph(Surv(start, stop, event) ~ X1 + X2 + X3, data = dt.long)

当前核心问题 #

直接使用survfit()为中期用户预测时存在两个关键缺陷：

1. 未考虑用户已存活至当前区间的事实：修改新用户的start/stop值，预测结果无变化。
2. 未识别时变协变量的区间特异性：不同区间的新用户（如6-7区间 vs 3-4区间）使用相同X3值时，预测结果完全一致，忽略了X3仅适用于当前区间的特性。

例如，定义6-7区间的新用户：

nd <- data.table("id"    = 7,
                 "start" = 6,
                 "stop"  = 7,
                 "X1"    = 30,
                 "X2"    = 1,
                 "X3"    = 20)

运行预测后得到的生存概率矩阵为：

[,1] [,2]    [,3]    [,4]    [,5]    [,6]    [,7]    [,8]    [,9]
[1,]   1    1 0.84815 0.69705 0.53001 0.53001 0.53001 0.27879 0.27879

将新用户改为3-4区间后，预测结果完全相同，不符合业务逻辑。

解决方案步骤 #

1. 构建新用户的完整区间长表 #

对于处于current_start到current_stop区间的新用户，生成从当前区间到后续所有可能区间的长表数据：

• 静态协变量X1、X2保持不变。
• 时变协变量X3：当前区间使用已知值，后续区间若无观测值，可根据业务场景选择（如保持当前值、使用训练集均值）。

以6-7区间的新用户为例：

# 假设后续区间X3保持当前值20
nd_long <- data.table(
  id = 7,
  start = c(6,7,8),  # 从当前区间开始到训练数据的最大观测区间
  stop = c(7,8,9),
  X1 = 30,
  X2 = 1,
  X3 = 20  # 当前及后续区间的X3取值
)

2. 计算条件生存概率 #

先预测从t=0开始的完整生存曲线，再基于用户已存活到当前区间的事实，计算条件生存概率：

# 预测完整生存曲线
sf_full <- survfit(cph.model, newdata = nd_long)
# 获取当前区间结束时间（7）对应的生存概率
S_current <- sf_full$surv[sf_full$time == 7]
# 计算后续时间点的条件生存概率（已存活到7的前提下）
conditional_surv <- sf_full$surv[sf_full$time >=7] / S_current
# 对应后续时间点
conditional_times <- sf_full$time[sf_full$time >=7]

3. 计算剩余期望生存期 #

剩余期望生存期是后续条件生存曲线下的面积，用梯形法计算：

# 计算相邻时间点的间隔
time_diff <- diff(conditional_times)
# 计算各区间的梯形面积
area_parts <- (conditional_surv[-length(conditional_surv)] + conditional_surv[-1])/2 * time_diff
# 处理最后一个区间的剩余面积（假设生存概率保持到最大观测区间结束）
last_area <- conditional_surv[length(conditional_surv)] * (max(dt.long$stop) - max(conditional_times))
# 总剩余期望生存期
remaining_expectation <- sum(area_parts) + last_area

关键注意事项 #

• 时变协变量的后续取值：需根据业务逻辑合理假设，这会直接影响预测结果的准确性。
• 条件生存概率的必要性：必须基于用户已存活到当前区间的事实进行归一化，否则预测结果是从t=0开始的无条件生存期，不符合中期用户的需求。

内容的提问来源于stack exchange，提问作者George