含时变协变量的Cox回归模型正确性验证及结果异常原因咨询

我旨在研究接受干细胞移植的血液病患者中非复发死亡率（TRM）的潜在风险因素。首先针对各变量构建单变量Cox回归模型以单独分析其效应，随后将p值≤0.20的变量纳入多变量Cox回归模型。

通过分析变量残差验证比例风险假设时，发现慢性GVHD（EICRc）变量的残差呈现明显规律分布，不满足假设，因此决定将其作为时变协变量构建Cox回归模型。

所用R语言代码及慢性GVHD残差图如下：

install.packages("survival")
library(survival)

# EICRc: chronic GVHD variable 
# TRM: non-relapse mortality (a dummy variable where 1 represents the event)
# T_seguimiento: time until the event (in months)

# here the univariate cox regression model, with the EICRc variable as non-time-dependent variable.
TRM_cGVHD <- coxph(Surv(T_seguimiento, TRM) ~ EICRc, data = dat)
test.ph <- cox.zph(TRM_cGVHD)
test.ph

# I check for the proportional hazard assumption:
plot(test.ph)

（残差图显示残差呈明显规律分布，不满足比例风险假设）

基于此，我使用survival包的tt()函数构建含时变协变量的模型，代码及结果如下：

##I use the tt() function of the survival package:

dat$EICRc <- relevel(as.factor(dat$EICRc), ref = "no")
TRM_cGVHD <- coxph(Surv(T_seguimiento, TRM) ~ tt(EICRc), data = dat)  
summary(TRM_cGVHD)

Call:
coxph(formula = Surv(T_seguimiento, TRM) ~ EICRc, data = dat)

  n= 226, number of events= 40 
   (1 observation deleted due to missingness)

           coef exp(coef) se(coef)     z Pr(>|z|)   
EICRcsi -0.8917    0.4100   0.3416 -2.61  0.00905 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

        exp(coef) exp(-coef) lower .95 upper .95
EICRcsi      0.41      2.439    0.2099    0.8008

Concordance= 0.651  (se = 0.028 )
Likelihood ratio test= 7.14  on 1 df,   p=0.008
Wald test            = 6.81  on 1 df,   p=0.009
Score (logrank) test = 7.19  on 1 df,   p=0.007

但预期中慢性GVHD应为非复发死亡率的风险因素，而模型结果显示其为保护因素（HR=0.41，p<0.01）。现咨询：该时变协变量Cox回归模型的构建是否正确？结果异常的可能原因有哪些？

一、时变协变量模型构建是否正确？ #

你当前的代码没有正确实现时变协变量的Cox模型。tt()函数的用法需要自定义时间变换逻辑，仅写tt(EICRc)不会自动将变量转换为时变形式，从输出的公式coxph(formula = Surv(T_seguimiento, TRM) ~ EICRc, data = dat)也能看出，模型实际还是按普通非时变协变量拟合的，根本没用到时变逻辑。

正确实现时变协变量的方式分两种，取决于你的数据结构：

1. 宽转长数据结构：如果数据是宽格式（每个患者一行），需先转换为长格式，新增时间区间列，明确每个区间内EICRc的状态。例如，若患者在移植后第6个月发生慢性GVHD，那么0-6个月区间EICRc为"no"，6个月之后为"si"。
2. 使用tt()自定义函数：必须为tt()指定时间变换函数，示例如下：

# 定义时间变换函数，假设EICRc的效应随时间线性变化
tt_func <- function(x, t, ...) {
  x * t  # 协变量与时间的交互项，代表效应随时间变化
}

TRM_cGVHD <- coxph(Surv(T_seguimiento, TRM) ~ tt(EICRc), 
                   data = dat, 
                   tt = tt_func)  # 传入自定义的tt函数

如果是分段时变（比如某个时间点后效应改变），也可以在函数里设置条件判断。

二、结果异常的可能原因 #

即使模型构建正确，出现HR与预期相反的情况，可能有以下原因：

• 变量编码或参考组设置错误：你将EICRc的参考组设为"no"，结果中EICRcsi的HR=0.41代表"si"组相对于"no"组TRM风险更低。需确认"si"是否确实代表发生慢性GVHD，有没有把编码含义搞反。
• 时间依赖关系的方向特殊：慢性GVHD的风险效应可能随时间变化——早期发生GVHD可能是风险因素，但后期发生反而可能有保护作用（比如提示移植物抗宿主效应同时控制了潜在疾病）。残差图显示不满足比例风险假设，说明效应随时间改变，整体拟合后可能出现反向HR。
• 混杂因素未控制：当前是单变量模型，未调整年龄、移植类型、预处理方案等混杂因素。这些因素可能与EICRc和TRM都相关，单变量分析结果会被混杂偏移，导致HR方向反转。
• 事件数较少：样本中仅40个TRM事件，小样本量可能导致估计结果不稳定，甚至出现与预期相反的效应。
• 数据缺失或错误：有1个观察因缺失被删除，需确认缺失的是否为关键数据；同时检查T_seguimiento的时间单位、TRM的事件定义是否正确（比如是否把复发误标记为TRM，或随访时间计算错误）。

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