如何在R的Cox生存模型中为分类协变量加入时变系数

我正在用R的survival包构建Cox比例风险（Cox PH）模型，想给分类变量加入时变系数。可复现的数据准备代码如下：

library(survival)
# Data
stanford <- stanford2
stanford$age_cat <- ifelse(stanford$age > 35, "old", "young")

尝试用tt()函数实现时变系数，但直接用分类变量报错，提示需要哑变量编码：

mod.fail <- coxph(Surv(time, status) ~ tt(age_cat),
             data = stanford,
             tt = function(x, t, ...) x*t)
# Error in x * t : non-numeric argument to binary operator

于是创建了哑变量：

# Create dummy coding of age_cat
stanford$age_cat_d <- ifelse(stanford$age_cat == "old", 1, 0)

但不确定如何正确指定模型，以下两个模型都能运行，但不清楚哪个能实现分类变量效应随时间变化的需求：

# Model 1
mod.t1 <- coxph(Surv(time, status) ~ tt(age_cat_d),
               data = stanford,
               tt = function(x, t, ...) x*t)
# Model 2
mod.t2 <- coxph(Surv(time, status) ~ age_cat_d + tt(age_cat_d),
                data = stanford,
                tt = function(x, t, ...) x*t)

计算t=200时两个模型中年龄分类变量的效应，结果有差异：

# Model 1
coef(mod.t1)[1]*200
# tt(age_cat_d) 
#    0.04425679

# Model 2
coef(mod.t2)[1]+coef(mod.t2)[2]*200
# age_cat_d 
# 0.5424105 

请问上述两个模型中，哪个是为分类协变量实现时变系数的正确方式？相关文档示例多针对连续变量，特此问询。（注：示例仅用于复现，不代表该数据适合构建此类模型）

Model 2是实现分类变量时变系数的正确方式，原因如下：

• Model 1的模型形式为：h(t|x) = h0(t) * exp( β * (age_cat_d * t) )，这意味着当t=0时，分类变量的效应为0，强制效应从0开始随时间线性变化，不符合分类变量的实际逻辑——我们通常需要先考虑t=0时的基线效应，再叠加随时间变化的部分。

• Model 2的模型形式为：h(t|x) = h0(t) * exp( β0 * age_cat_d + β1 * (age_cat_d * t) )，对应的分类变量（old vs young）的系数为β0 + β1*t，其中：

  • β0是t=0时的基线系数
  • β1是系数随时间变化的斜率
    这完全符合“时变系数”的需求：分类变量的效应随时间t线性变化，既保留了初始效应，又体现了效应的时间依赖性。

你计算的t=200时的效应差异也能验证这一点：Model 2的结果包含了基线效应β0，而Model 1仅包含随时间变化的部分，显然Model 2的设定更贴合时变系数的定义。

另外补充：无需手动创建哑变量，也可以直接对因子变量使用tt()，比如在tt函数中处理因子的数值化转换，但你当前的哑变量方式是可行的，Model 2的设定是正确的。

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