我以前只有 Cox PH 模型及其假设检查的经验。现在,我第一次拥有我的客户数据,其中大多数协变量随时间变化,只有少数是固定值的。我确实将它设置为启动/停止格式,并且我还安装了一些模型,但 cox.zph 显示的 p 值非常小。我的理解是,扩展的 Cox 模型不是 PH 模型。我什至需要对时间相关变量进行 cox.zph 测试吗?或者当包括固定协变量时,我是否需要查看这些变量的大 p 值?
扩展 Cox 模型和 cox.zph
扩展的 Cox 模型在技术上实际上与常规 Cox 模型相同。如果您的数据集被正确构建以适应时间相关的协变量(每个主题多行,开始和结束时间等),cox.ph
那么cox.zph
应该可以很好地处理您的数据。
具有时间相关的协变量并不会改变您应该检查比例假设的事实,在这种情况下,使用 Schoenfeld 残差对转换后的时间使用cox.zph
. 具有非常小的 p 值表明您需要注意与时间相关的系数。
两种主要方法是(1)时间交互和(2)阶跃函数。前者更容易操作和阅读,但如果它不改变 p 值而不是使用后者:
请注意,如果您提供自己的数据会更容易,因此以下内容基于我使用的示例数据
(1) 与时间的互动
在这里,我们对有问题的变量使用与时间的简单交互。请注意,您不需要将时间本身添加到模型中,因为它是基线。
> model.coxph0 <- coxph(Surv(t1, t2, event) ~ female + var2, data = data)
> summary(model.coxph0)
coef exp(coef) se(coef) z Pr(>|z|)
female 0.1699562 1.1852530 0.1605322 1.059 0.290
var2 -0.0002503 0.9997497 0.0004652 -0.538 0.591
检查比例假设违规:
> (viol.cox0<- cox.zph(model.coxph0))
rho chisq p
female 0.0501 1.16 0.2811
var2 0.1020 4.35 0.0370
GLOBAL NA 5.31 0.0704
所以var2
有问题。让我们尝试使用与时间的交互:
> model.coxph0 <- coxph(Surv(t1, t2, event) ~
+ female + var2 + var2:t2, data = data)
> summary(model.coxph0)
coef exp(coef) se(coef) z Pr(>|z|)
female 1.665e-01 1.181e+00 1.605e-01 1.038 0.29948
var2 -1.358e-03 9.986e-01 6.852e-04 -1.982 0.04746 *
var2:t2 5.803e-05 1.000e+00 2.106e-05 2.756 0.00586 **
现在让我们再次检查 zph:
> (viol.cox0<- cox.zph(model.coxph0))
rho chisq p
female 0.0486 1.095 0.295
var2 -0.0250 0.258 0.611
var2:t2 0.0282 0.322 0.570
GLOBAL NA 1.462 0.691
如您所见 - 这就是门票。
(2) 阶跃函数
在这里,我们根据残差的绘制方式创建了一个按时间段划分的模型,并将层添加到特定的问题变量。
> model.coxph1 <- coxph(Surv(t1, t2, event) ~
female + contributions, data = data)
> summary(model.coxph1)
coef exp(coef) se(coef) z Pr(>|z|)
female 1.204e-01 1.128e+00 1.609e-01 0.748 0.454
contributions 2.138e-04 1.000e+00 3.584e-05 5.964 2.46e-09 ***
现在使用 zph:
> (viol.cox1<- cox.zph(model.coxph1))
rho chisq p
female 0.0296 0.41 5.22e-01
contributions 0.2068 21.31 3.91e-06
GLOBAL NA 22.38 1.38e-05
> plot(viol.cox1)
因此,该contributions
系数似乎与时间有关。我尝试与无效的时间互动。所以这里是使用阶跃函数:您首先需要查看图表(上图)并目视检查线条在哪里改变角度。这里似乎是时间法术 8 和 40 左右。所以我们将survSplit
在上述时间使用分组创建数据:
sandbox_data <- survSplit(Surv(t1, t2, event) ~
female +contributions,
data = data, cut = c(8,40), episode = "tgroup", id = "id")
然后运行模型strata
:
> model.coxph2 <- coxph(Surv(t1, t2, event) ~
female + contributions:strata(tgroup), data = sandbox_data)
> summary(model.coxph2)
coef exp(coef) se(coef) z Pr(>|z|)
female 1.249e-01 1.133e+00 1.615e-01 0.774 0.4390
contributions:strata(tgroup)tgroup=1 1.048e-04 1.000e+00 5.380e-05 1.948 0.0514 .
contributions:strata(tgroup)tgroup=2 3.119e-04 1.000e+00 5.825e-05 5.355 8.54e-08 ***
contributions:strata(tgroup)tgroup=3 6.894e-04 1.001e+00 1.179e-04 5.845 5.06e-09 ***
还有中提琴——
> (viol.cox1<- cox.zph(model.coxph1))
rho chisq p
female 0.0410 0.781 0.377
contributions:strata(tgroup)tgroup=1 0.0363 0.826 0.364
contributions:strata(tgroup)tgroup=2 0.0479 0.958 0.328
contributions:strata(tgroup)tgroup=3 0.0172 0.140 0.708
GLOBAL NA 2.956 0.565
这个问题值得在一些帐户上得到更新的答案。首先,随着最新版本的生存包cox.zph()
,该功能发生了很大变化,因此输出看起来不一样可能会造成混淆。其次,当您处理与时间相关的协变量时,可能会有一些隐藏的“陷阱” ,就像在这个问题中一样。第三,尽管其他答案的大部分内容都很好,但其中一种指定时间相关系数的建议方法可能存在严重错误。最后,处理最后一个问题的正确方法使得(至少目前)无法用于检查最终模型中的比例风险(PH)。cox.zph()
多年来,该
cox.zph()
函数使用近似值、缩放 Schoenfeld 残差和(可能转换的)时间之间的相关系数来执行其 PH 测试。如另一个答案所示,该相关系数被报告为“rho”。由于版本 3.0-10 的包,cox.zph()
现在是一个精确的分数测试。不再需要报告“rho”值。对于与时间相关的协变量,因果关系可能存在问题。例如,我最近帮助分析了一些数据,其中将患者对慢性病处方药的使用作为协变量包括在内。随着人们年龄的增长,他们越来越有可能使用这种药物。将该药物作为时间依赖性协变量包括在内是有问题的,因为它可能只是已经存活更长时间的标志。与时间相关的协变量很容易成为更长生存期的代表,这(除了因果关系问题)可能会显示为 PH 问题。我怀疑这可能是这里最初问题中问题的一部分。引用Therneau、Crowson 和 Atkinson的时间相关小插曲:
Cox 模型中时间相关协变量的关键规则很简单,并且与赌博基本相同:你不能展望未来。
- 无论是否存在与时间相关的协变量,时间相关系数都可以帮助解决 PH 问题。作为时间函数的阶跃函数建模系数是另一个答案中提出的方法之一,它是有效的。正如@bandwagoner 在对该问题的评论中指出的那样,另一种提议的方法,即协变量时间交互,可能不是。* 再次引用小插图:
这个错误经常发生,以至于[at]
coxph
例程已更新为打印错误消息以进行此类尝试。问题是上面的代码实际上并没有创建一个时间相关的协变量,而是根据每个主题的协变量值创建一个时间静态值time
;coxph
与我们在调用之外构造变量没有什么不同。这个变量绝对打破了不展望未来的规则,人们很快就会发现循环性:大的值time
似乎预示着长期的生存,因为长期的生存导致 的值很大time
。
- 该
survival
包提供了一种通过用户定义的函数将系数指定为任意时间函数的正确方法tt()
。不幸的是,正如包的 NEWS 文件所说,从版本 3.1-2 开始,“cox.zph 命令现在拒绝带有tt()
术语的模型,在它有一个不正确的计算之前。” 因此,目前看来,时间相关系数的评估将取决于用户定义的tt()
函数与时间相关系数和其他图形评估所看到的时间相关性形式的匹配程度。
* Yuval Spiegler 的回答没有具体说明数据准备的全部性质。如果它是通过Fox 和 Weisbergunfold()
使用的函数来完成的,那么对于每个处于风险的时间,每个人都有一个单独的行。使用这种格式,设计矩阵将包含在任何事件时间处于风险中的所有个人的当前交互值。如果另一个答案使用了以这种方式准备的数据,那么使用显式交互项进行分析就可以了。包时间相关的小插图使用的函数产生的数据stop, start, event
covariate:time
covariate:time
start, stop, event
tmerge()
survival
不会为每个有风险的时间生成单独的行;它将个人的数据分解为具有恒定协变量值的完整时期。使用这种(更短的)数据格式,您必须使用该tt()
功能来正确指定covariate:time
交互。