如何在R语言的Cox比例风险模型中处理时变变量？

核心结论 #

你不能直接用 coxph(Surv(time, event) ~ tr1 * tr2 + tr1 * tr3 + tr2 * tr3, data) 拟合，因为这个公式把tr1/tr2/tr3当作固定协变量处理，完全没体现它们的时间依赖性。要处理这类基于时间节点划分的多阶段时变变量，必须先将宽格式数据转换为计数过程格式（长数据），再用Surv(start, stop, event)的形式建模。

问题本质 #

你的tr1/tr2/tr3是划分时间阶段的节点，每个个体的生存历程被这些节点分割成多个风险区间，不同区间内的暴露状态（或风险特征）不同。比如：

• 个体6的生存时间7.52介于tr2=7.33和tr3=9.11之间，其历程分为：[0,4.63)（对应tr1阶段）、[4.63,7.33)（对应tr2阶段）、[7.33,7.52]（事件发生）
• 个体1的生存时间远大于tr3，历程分为4个区间：[0,11.6)、[11.6,63.85)、[63.85,66.51)、[66.51,125.81]

具体实现步骤 #

1. 数据预处理（宽转长）

先将原始宽格式数据转换为计数过程的长格式，每个区间对应一行数据，包含区间起始/结束时间、事件指示、该区间的时变协变量状态：

library(survival)
library(dplyr)

# 构造原始数据框
tr1 <- c(11.6, 19.04, NA, NA, 10.39, 4.63)
tr2 <- c(63.85, 20.29, NA, NA, 11.64, 7.33)
tr3 <- c(66.51, 22.95, NA, NA, 16.05, 9.11)
if_event <- c(1, 1, 0, 0, 1, 1)
sur_time <- c(125.81, 30.23, 59.27, 161.17, 51.36, 7.52)
income <- c(103844, 57246, 83056, 38380, 37518, 900)
population <- c(10000, 20000, 40000, 88000, 90000, 102034)

dat <- data.frame(tr1, tr2, tr3, if_event, sur_time, income, population)

# 定义函数生成每个个体的计数过程行
create_counting_process <- function(row) {
  # 收集所有时间节点并去重排序
  times <- c(0, row$tr1, row$tr2, row$tr3, row$sur_time) %>%
    na.omit() %>%
    unique() %>%
    sort()
  
  # 生成区间的start和stop
  start <- times[-length(times)]
  stop <- times[-1]
  
  # 事件指示：仅最后一个区间可能发生事件
  event <- rep(0, length(start))
  event[length(event)] <- row$if_event
  
  # 生成时变暴露变量：exp1=1表示当前区间在tr1之前，以此类推
  exp1 <- as.integer(stop <= row$tr1)
  exp2 <- as.integer(stop <= row$tr2)
  exp3 <- as.integer(stop <= row$tr3)
  
  # 处理NA情况：无对应时间节点则暴露为0
  exp1[is.na(row$tr1)] <- 0
  exp2[is.na(row$tr2)] <- 0
  exp3[is.na(row$tr3)] <- 0
  
  # 合并为数据框
  data.frame(
    id = row.names(row),
    start, stop, event,
    exp1, exp2, exp3,
    income = row$income,
    population = row$population
  )
}

# 应用到所有个体，得到长格式数据
dat_long <- dat %>%
  rowwise() %>%
  do(create_counting_process(.)) %>%
  ungroup()

2. 拟合Cox模型

用Surv(start, stop, event)定义生存对象，即可正确纳入时变变量：

• 主效应模型：拟合每个阶段暴露的独立效应

model_main <- coxph(
  Surv(start, stop, event) ~ exp1 + exp2 + exp3 + income + population,
  data = dat_long
)
summary(model_main)

• 交互效应模型：如果你需要拟合阶段之间的交互作用（对应你原本想做的tr1*tr2等），可以直接添加交互项：

model_interact <- coxph(
  Surv(start, stop, event) ~ exp1*exp2 + exp1*exp3 + exp2*exp3 + income + population,
  data = dat_long
)
summary(model_interact)

• 阶段分类模型：如果想直接对比不同阶段的风险差异，可以生成阶段变量建模：

# 先在长数据中添加阶段变量
dat_long <- dat_long %>%
  mutate(stage = case_when(
    exp1 == 1 & exp2 == 1 & exp3 == 1 ~ 1,
    exp1 == 0 & exp2 == 1 & exp3 == 1 ~ 2,
    exp1 == 0 & exp2 == 0 & exp3 == 1 ~ 3,
    TRUE ~ 4
  ))

# 拟合阶段效应模型
model_stage <- coxph(
  Surv(start, stop, event) ~ factor(stage) + income + population,
  data = dat_long
)
summary(model_stage)

关键注意点 #

• 原始公式的错误在于：将时变节点当作固定数值，忽略了时间区间内的风险变化，无法反映个体在不同阶段的暴露状态。
• 对于tr1/tr2/tr3全为NA的个体（如样本3、4），直接生成[0, sur_time]单个区间，暴露变量全为0即可。

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