当 R 中的生存分析中违反比例假设时，如何对协变量与时间的相互作用进行建模

在 R 中，当比例检验（使用 coxph）显示违反了 Cox 模型中的比例假设时，合并协变量和时间之间的交互项的最佳方法是什么？我知道您可以使用分层或与时间项交互，我对后者感兴趣。我无法在互联网上找到明确的解释以及如何执行此操作的示例。在使用 Rossi 数据集的最常见示例中，Fox 建议这样做：

coxph(formula = Surv(start, stop, arrest.time) ~ fin + age + age:stop + prio, data = Rossi.2)

使用age:stop 建模与age:start 建模之间有区别吗？公式必须使用这种格式吗？如果我使用带有两个参数格式的 Surv，以下内容是否也有意义？

coxph(formula = Surv(week, arrest) ~ fin + age + age:week + prio, data = Rossi)

或者您必须拆分数据集并使用 Surv(start,stop,event) 方法？ 另外，还有时间变换方法，所以，

coxph(formula = Surv(week, arrest) ~ fin + age + tt(age) + prio, data = Rossi, tt=function(x,t,...) x*t)

我知道有些人更喜欢模型log(t)代替t这里。但其中哪一种是模拟与时间相互作用的正确方法呢？这些都涉及相同/不同的基础统计模型吗？最后，都是建模（对于交互项）：h(t) = h0(t)exp(b*X*t)?

这本质上是一个由三部分组成的问题：

1. 如何估计随时间变化的影响？
2. 使用不同规格的时变效果有什么区别survival::coxph功能
3. 如何确定时间变化的形状，即线性、对数……

下面我将尝试使用资深数据示例来回答这些问题，该示例在 4.2 节中进行了介绍。关于时间相关协变量和时间相关系数的小插图 https://cran.microsoft.com/web/packages/survival/vignettes/timedep.pdf（也称为时变效应）survival包裹：

library(dplyr)
library(survival)

data("veteran", package = "survival")
veteran <- veteran %>%
  mutate(
    trt   = 1L * (trt == 2),
    prior = 1L * (prior == 10))
head(veteran)
#>   trt celltype time status karno diagtime age prior
#> 1   0 squamous   72      1    60        7  69     0
#> 2   0 squamous  411      1    70        5  64     1
#> 3   0 squamous  228      1    60        3  38     0
#> 4   0 squamous  126      1    60        9  63     1
#> 5   0 squamous  118      1    70       11  65     1
#> 6   0 squamous   10      1    20        5  49     0

1. 如何估计时变效应

有不同的流行方法和实现，例如survival::coxph, timereg::aalen或在适当的数据转换后使用 GAM（见下文）。

尽管具体方法及其实现有所不同，但总体思路是创建一个长格式数据集，其中

• 后续活动被划分为多个时间间隔
• 对于每个主题，除了最后一个（如果是事件）之外的所有间隔中的状态均为 0
• 时间变量在每个时间间隔内更新

然后，时间（或时间的变换，例如 log(t)）只是一个协变量，并且可以通过感兴趣的协变量与时间（变换后的）协变量之间的相互作用来估计随时间变化的效应。

如果时间变化的函数形式已知，则可以使用tt()方法：

cph_tt <- coxph(
      formula = Surv(time, status) ~ trt + prior + karno + tt(karno),
      data    = veteran,
      tt      = function(x, t, ...) x * log(t + 20))

2、不同规格的时变效果使用时有什么区别survival::coxph功能

没有区别。我假设tt()函数只是通过转换为长格式进行估计的捷径。您可以使用以下代码验证这两种方法是否等效：

转换为长格式

veteran_long <- survSplit(Surv(time, status)~., data = veteran, id = "id",
  cut = unique(veteran$time)) %>%
  mutate(log_time = log(time + 20))
head(veteran_long) %>% select(id, trt, age, tstart, time, log_time, status)
#>   id trt age tstart time log_time status
#> 1  1   0  69      0    1 3.044522      0
#> 2  1   0  69      1    2 3.091042      0
#> 3  1   0  69      2    3 3.135494      0
#> 4  1   0  69      3    4 3.178054      0
#> 5  1   0  69      4    7 3.295837      0
#> 6  1   0  69      7    8 3.332205      0

cph_long <- coxph(formula = Surv(tstart, time, status)~
    trt + prior + karno + karno:log_time, data = veteran_long)

## models are equivalent, just different specification
cbind(coef(cph_long), coef(cph_tt))
#>                       [,1]        [,2]
#> trt             0.01647766  0.01647766
#> prior          -0.09317362 -0.09317362
#> karno          -0.12466229 -0.12466229
#> karno:log_time  0.02130957  0.02130957

3. 如何确定时间变化的形状？

如前所述，时变效应只是协变量的相互作用x和时间t，因此时变效应可以有不同的规格，相当于标准回归模型中的交互作用，例如

• x*t：线性协变量效应，线性时变效应
• f(x)*t：非线性协变量效应，线性时变效应
• f(t)*x：线性协变量效应，非线性时变（对于分类 x）这本质上代表了分层基线风险
• f(x, t)：非线性、非线性时变效应

在每种情况下，效果的函数形式f可以根据数据估计或预先指定（例如f(t)*x = karno * log(t + 20)多于）。

在大多数情况下，您更愿意估计f从数据来看。据我所知，对此类影响（惩罚性）估计的支持仅限于survival包裹。但是，您可以使用mgcv::gam估计上面指定的任何影响（在适当的数据转换之后）。下面给出一个例子并显示效果karno随着时间的推移趋向于 0，无论后续开始时的卡诺夫斯基分数如何（参见here https://adibender.github.io/pammtools/articles/tveffects.html详细信息以及第 4.2 节here https://arxiv.org/pdf/1806.01042.pdf):

library(pammtools)
# data transformation
ped <- as_ped(veteran, Surv(time, status)~., max_time = 400)
# model
pam <- mgcv::gam(ped_status ~ s(tend) + trt + prior + te(tend, karno, k = 10),
  data = ped, family = poisson(), offset = offset, method = "REML")
p_2d <- gg_tensor(pam)
p_slice <- gg_slice(ped, pam, "karno", tend = unique(tend), karno = c(20, 50, 80), reference = list(karno = 60))
gridExtra::grid.arrange(p_2d, p_slice, nrow = 1)