如Singer 和 Willett的Applied Longitudinal Data Analaysis第 15 章(第 15.3 节)所述,我正在尝试运行一个非比例 cox 回归模型,该模型具有与时间的相互作用变量。但是,我似乎无法得到与这本书一致的答案。
本书中使用的数据和源代码都在这个奇妙的网站上提供。不幸的是,最后一章没有提供 R 代码,并且为文本中讨论的示例提供的 R 数据集不完整,并且为最简单的模型(我知道如何运行)提供了不正确的答案。相反,要获取此示例的完整数据集,必须单击“SAS”列中的“下载”链接(具有正确的数据集),然后在安装haven
包后(允许读取外国数据格式) ,通过以下方式读取相关数据集:
haven::read_sas("alda/lengthofstay.sas7bdat")
该数据集表示参与者在医院住院治疗中的(变量)住院ID
时间(变量)。DAYS
审查变量是CENSOR
。研究人员假设两种不同类型的治疗(二元变量TREAT
)可以预测检查出治疗风险的差异值。此外,他们预计组间的危险差异不会随着时间的推移而保持不变,因此需要创建一个交互项。我可以让简单的主效应模型工作,返回书中报告的相同危险系数(这就是我最终发现 R 代码提供的 .csv 文件不完整的原因)。
summary(modA <- coxph(Surv(DAYS,1-CENSOR) ~ TREAT, data = los))
coef exp(coef) se(coef) z Pr(>|z|)
TREAT 0.1457 1.1568 0.1541 0.945 0.345
我尝试遵循此处和此处列出的程序以及其中列出的来源(例如,Therneau vignette on time-variable covariates on the survival
package),当然,当我复制粘贴其他人的代码并运行它时一切正常。但是我正在尝试使用一个数据集从头开始为自己做这件事,我可以将其结果与我的结果进行比较。而我就是无法让它发挥作用。
首先我创建了一个 EVENT 变量
los$EVENT <- 1 - los$CENSOR
数据集中存在导致问题的重复 ID 号。所以我们必须将其更改为新的ID号
los$ID[which(duplicated(los$ID))] <- 842
现在,根据我在此处和此处阅读 的内容,需要拆分数据框,以便对于每个参与者,当任何其他参与者经历事件EVENT
时,在他们的事件(或审查)时间之前的每个时间点都有一行指示状态。因此,我们需要创建一个包含所有唯一事件时间的向量,然后在这些事件时间上拆分数据集
cutPoints <- sort(unique(los$DAYS[los$EVENT == 1]))
# now split the dataset
longLOS <- survSplit(Surv(DAYS,EVENT)~ ., data = los, cut = cutPoints)
# and (just because I'm anal) rename the interval upper bound column (formerly "DAYS")
names(longLOS)[5] <- "tstop"
当我查看这个数据集时,它似乎就是我所追求的,(1)每个参与者的行数与数据集中其他人经历事件的事件时间之前的间隔一样多,(2)两列表示每个间隔的下限和上限,以及 (3) 事件列,当受访者没有经历事件时,所有行的所有行都为 0,当他们经历事件或被审查时,最后一行为 1。
接下来,我创建了与时间的交互变量,从“间隔上限”列中减去 1,以便主效应TREAT
代表住院第一天的治疗效果。
longLOS$TREATINT <- longLOS$EVENT*(longLOS$tstop - 1)
并运行模型
summary(modB <- coxph(Surv(tstart, tstop, EVENT) ~ TREAT + TREATINT, data = longLOS))
但它不起作用!我收到了(相当无用的)错误消息
Error in fitter(X, Y, strats, offset, init, control, weights = weights, :
routine failed due to numeric overflow.This should never happen. Please contact the author.
我究竟做错了什么?近三年来,我一直在慢慢地通过 Singer 和 Willett 工作(我还是研究生时开始的),现在最后一章被证明是迄今为止我最大的挑战。我还有三十页要写;任何帮助将不胜感激。