r - 具有交互变量的非比例风险 (Cox) 模型的计数过程数据集

Question

我正在尝试运行一个具有时间交互变量的非比例 cox 回归模型，如 Singer 和 Willett 所著的应用纵向数据分析 第 15 章(第 15.3 节)所述。但是我似乎无法得到与这本书一致的答案。

本书中使用的数据和源代码在此 fantastic website 提供。 .不幸的是，最后一章没有提供 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

我尝试按照规定的程序进行操作 here , 和 here ，以及其中列出的来源(例如 survival 包中关于时变协变量的 Thenneau vignette)，当然，当我复制粘贴别人的代码并运行它时，一切正常。但是我正在尝试使用一个数据集从头开始为自己做这件事，我可以将其结果与我的结果进行比较。我就是无法让它发挥作用。

首先我创建了一个事件变量

los$EVENT <- 1 - los$CENSOR

数据集中存在导致问题的重复 ID 号。所以我们要换一个新的身份证号码

los$ID[which(duplicated(los$ID))] <- 842

现在，根据我阅读的内容 here和 here需要拆分数据框，以便对于每个参与者，都有一行指示 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(我还是一名研究生时就开始了)，现在事实证明，最后一章是我迄今为止最大的挑战。我还有三十页要写；任何帮助将不胜感激。

Answer 1

我发现我做错了什么。我创建交互变量时出现了一个愚蠢的错误 TREATINT .而不是

longLOS$TREATINT <- longLOS$EVENT*(longLOS$tstop - 1)

应该是

longLOS$TREATINT <- longLOS$TREAT*(longLOS$tstop - 1)

现在运行模型

summary(modB <- coxph(Surv(tstart, tstop, EVENT) ~ TREAT + TREATINT, data = longLOS))

它不仅有效，而且产生的系数与 Singer 和 Willett 书中报告的系数相匹配。

              coef exp(coef)  se(coef)      z Pr(>|z|)
TREAT     0.706411  2.026705  0.292404  2.416   0.0157
TREATINT -0.020833  0.979383  0.009207 -2.263   0.0237

鉴于我的错误是多么愚蠢，我很想删除整篇文章，但我想我会把它留给像我这样想知道如何在 R 中与时间 Cox 模型进行交互的其他人。