4

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

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

使用年龄：停止与年龄：开始建模之间有区别吗？公式是否必须使用这种格式？如果我使用具有两个参数格式的 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)：？

4

1 回答 1

8

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

  1. 如何估计时变效应？
  2. survival::coxph不同规格的时变效果使用函数有什么区别
  3. 如何确定时变的形状，即线性、对数、...

我将尝试在下面使用经验丰富的数据示例来回答这些问题，该示例在包中时间相关协变量和时间相关系数（也称为时变效应）的小插图的第 4.2 节中有特色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))

survival::coxph2.不同规格的时变效果使用函数有什么区别

没有区别。我假设该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随着时间的推移，无论后续开始时的 Karnofsky 分数如何，随着时间的推移，效果趋于 0（有关详细信息，请参见此处以及此处的第 4.2 节）：

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)

于 2018-12-09T01:00:58.237 回答