1

我正在尝试使用该survival包使用左截断数据来拟合生存模型,但是我不确定语法是否正确。

假设我们正在测量受雇年龄 ( age) 和工作类型 ( parttime) 对公共卫生诊所医生就业时间的影响。医生是否退出或被审查由censor变量表示(0 表示退出,1 表示审查)。这种行为是在 18 个月的窗口中测量的。退出或审查的时间由两个变量表示,entry(开始时间)和exit(停止时间)表示医生在诊所受雇的年数。如果医生在窗口“打开”后开始就业,他们的entry时间设置为 0。如果他们在窗口“打开”之前开始就业,他们的entry时间表示当窗口“打开”时他们已经在该职位上工作了多长时间,他们exit时间是从他们最初被雇用时起,他们要么辞职,要么被“关闭”窗口审查。我们还假设就业时间和就业时间之间存在双向交互作用ageexit)。

这是玩具数据集。它比普通数据集小得多,因此在survival给定数据结构的情况下,估计本身并不像包含的语法和变量(使用 R 中的包)是否正确那么重要。玩具数据与 Singer 和 Willet 的Applied Longitudinal Data Analysis第 15 章中讨论的数据集具有完全相同的结构。我试图匹配他们报告的结果,但没有成功。网上没有很多关于如何在 R 中对左截断数据进行生存分析的明确信息,以及为本书提供代码的网站(这里) 不提供相关章节的 R 代码。在 R 中建模时变协变量和交互效应的方法非常复杂,我只是想知道我是否遗漏了一些重要的东西。

这是玩具数据

id <- 1:40
entry <- c(2.3,2.5,2.5,1.2,3.5,3.1,2.5,2.5,1.5,2.5,1.4,1.6,3.5,1.5,2.5,2.5,3.5,2.5,2.5,0.5,rep(0,20))
exit <- c(5.0,5.2,5.2,3.9,4.0,3.6,4.0,3.0,4.2,4.0,2.9,4.3,6.2,4.2,3.0,3.9,4.1,4.0,3.0,2.0,0.2,1.2,0.6,1.9,1.7,1.1,0.2,2.2,0.8,1.9,1.2,2.3,2.2,0.2,1.7,1.0,0.6,0.2,1.1,1.3)
censor <- c(1,1,1,1,0,0,0,0,1,0,0,1,1,1,0,0,0,0,0,0,rep(1,20))
parttime <- c(1,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,1,1,0,0,0,0,0,0,0,0,0)
age <- c(34,28,29,38,33,33,32,28,40,30,29,34,31,33,28,29,29,31,29,29,30,37,33,38,34,37,37,40,29,38 ,49,32,30,27,35,34,35,30,35,34)

doctors <- data.frame(id,entry,exit,censor,parttime,age)

现在为模型。

coxph(Surv(entry, exit, 1-censor) ~ parttime + age + age:exit, data = doctors)

给定数据结构和我们想知道的内容,这是指定模型的正确方法吗?此处的答案表明它是正确的,但我不确定是否正确指定了交互变量。

4

1 回答 1

1

通常情况下,直到我在 SO 上发布有关某个问题的问题后,我才会自己解决问题。如果与时间预测器有交互,我们需要将数据集转换为计数过程,人周期格式(即长格式)。这是因为每个参与者都需要一个时间间隔来跟踪他们在事件发生在数据集中其他任何人身上的每个时间点的事件状态,直到他们退出研究为止。

首先让我们创建一个事件变量

doctors$event <- 1 - doctors$censor

在我们运行 cox 模型之前,我们需要使用包survSplit中的函数survival。为此,我们需要制作事件发生时所有时间点的向量

cutPoints <- order(unique(doctors$exit[doctors$event == 1]))

现在我们可以将它传递给survSplit函数来创建一个新的数据集......

docNew <- survSplit(Surv(entry, exit, event)~.,
                   data = doctors,
                   cut = cutPoints,
                   end = "exit")

...然后我们在上面运行我们的模型

coxph(Surv(entry,exit,event) ~ parttime + age + age:exit, data = docNew)

瞧!

于 2018-05-01T20:25:11.433 回答