我正在使用 R 中的生存包构建 Cox PH 模型,并希望包含我的分类变量的时间相关系数。可重复的数据设置:
library(survival)
# Data
stanford <- stanford2
stanford$age_cat <- ifelse(stanford$age > 35, "old", "young")
从时间相关的小插图这里开始工作,用于生存包,我需要使用
tt()
函数。尝试 1 显示我需要虚拟编码。
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)
下面是我认为我们应该如何估计每个模型中 time=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
那么,上述任一模型是否是实现年龄类别的时间相关系数的正确方法?链接的小插图中的示例(以及我发现的其他使用
tt()
的指南)重点关注连续变量的时间相关系数。 (注意:上面的例子只是为了重现性;我并不是说我们应该为给定的数据创建这样一个与时间相关的模型)
[1]:https://cran.r-project.org/web/packages/survival/vignettes/timedep.pdf
由于
tt()
声明了时变系数的变换,无论您的协变量是连续的还是离散的,当您从时变系数 Cox 模型中删除“主项”时,这是一个理解您所拟合的模型的问题以及如何解释参数估计。
回答这个问题的最简单方法可能是浏览不同的模型规范(通过语法)并解释它们在做什么。
library(survival)
# Data
stanford <- stanford2
stanford$age_b <- ifelse(stanford$age > 35, 1, 0) # add binary covariate
# Function computing time functional form
myfun <- function(x, t, ...){ x * log(t + 20)}
年龄具有一种时不变效应。
coxph(Surv(time, status) ~ age, data = stanford)
#> Call:
#> coxph(formula = Surv(time, status) ~ age, data = stanford)
#>
#> coef exp(coef) se(coef) z p
#> age 0.02917 1.02960 0.01064 2.741 0.00613
#>
#> Likelihood ratio test=8.27 on 1 df, p=0.004034
#> n= 184, number of events= 113
年龄的影响现在也随着时间的推移而变化。年龄的总影响被分解为时不变项(
age
上的系数)和时变项(tt(age)
上的系数)。根据 tt()
使用的函数,本例中年龄的总影响为 -.007 + .007*log(t+20)。这种解释在时变系数插图中提供。
coxph(Surv(time, status) ~ age + tt(age),
tt = myfun,
data = stanford)
#> Call:
#> coxph(formula = Surv(time, status) ~ age + tt(age), data = stanford,
#> tt = myfun)
#>
#> coef exp(coef) se(coef) z p
#> age -0.007256 0.992770 0.042434 -0.171 0.864
#> tt(age) 0.007182 1.007208 0.008190 0.877 0.381
#>
#> Likelihood ratio test=9.04 on 2 df, p=0.01086
#> n= 184, number of events= 113
与模型 2 类似,我们让年龄的影响随时间变化。然而,我们不再单独估计时变分量和时不变分量。相反,我们直接估计年龄的总体影响,它会随着时间的推移而变化。年龄的总影响是 0.006*log(t+20)。
coxph(Surv(time, status) ~ tt(age),
tt = myfun,
data = stanford)
#> Call:
#> coxph(formula = Surv(time, status) ~ tt(age), data = stanford,
#> tt = myfun)
#>
#> coef exp(coef) se(coef) z p
#> tt(age) 0.005829 1.005846 0.002046 2.849 0.00439
#>
#> Likelihood ratio test=9.02 on 1 df, p=0.002677
#> n= 184, number of events= 113
现在让我们尝试用二元协变量而不是连续协变量来拟合这些模型。系数估计值发生变化,但它们仍然代表与时间相关的相同概念。
与模型 1 相同连续:年龄具有一种时不变效应。现在,该效应不再是连续年龄变化 1 个单位的效应,而是年老而不是年轻的效应。
coxph(Surv(time, status) ~ age_b, data = stanford)
#> Call:
#> coxph(formula = Surv(time, status) ~ age_b, data = stanford)
#>
#> coef exp(coef) se(coef) z p
#> age_b 0.2721 1.3128 0.2304 1.181 0.238
#>
#> Likelihood ratio test=1.47 on 1 df, p=0.2258
#> n= 184, number of events= 113
### Model 2 Binary: binary covariate, adding time-varying coefficient
Same as Model 2 Continuous: the effect of age now also varies across time. The total effect of age is decomposed into a time-invariant term (the coefficient on age) and a time-varying term (the coefficient on tt(age)). The total effect of age in this example is .025 + .050*log(t+20) based on the function used for `tt()`. That is the effect of being old rather than young.
```r
coxph(Surv(time, status) ~ age_b + tt(age_b),
tt = myfun,
data = stanford)
#> Call:
#> coxph(formula = Surv(time, status) ~ age_b + tt(age_b), data = stanford,
#> tt = myfun)
#>
#> coef exp(coef) se(coef) z p
#> age_b 0.02475 1.02506 0.92143 0.027 0.979
#> tt(age_b) 0.04680 1.04791 0.16956 0.276 0.783
#>
#> Likelihood ratio test=1.54 on 2 df, p=0.4621
#> n= 184, number of events= 113
我们现在再次估计年老与年轻的总时变效应,而不是将总效应分解为时变和时不变分量。
coxph(Surv(time, status) ~ tt(age_b),
tt = myfun,
data = stanford)
#> Call:
#> coxph(formula = Surv(time, status) ~ tt(age_b), data = stanford,
#> tt = myfun)
#>
#> coef exp(coef) se(coef) z p
#> tt(age_b) 0.05121 1.05255 0.04239 1.208 0.227
#>
#> Likelihood ratio test=1.54 on 1 df, p=0.2142
#> n= 184, number of events= 113