具有交互变量的非比例风险（Cox）模型的计数过程数据集

Question

我正在尝试运行一个非比例风险回归模型，其中包含一个与时间交互作用的变量，正如Singer和Willett在《应用纵向数据分析》第15章（第15.3节）中所述。然而，我似乎无法得到与书中相符的答案。

本书和源代码使用的数据可以在这个fantastic website上找到。不幸的是，最后一章没有提供R代码，并且在文本中讨论的示例的R所提供的数据集是不完整的，并且对于最简单的模型（我知道如何运行）提供了错误的答案。相反，要获取此示例的完整数据集，必须单击“SAS”列中的“Download”链接（具有正确数据集），然后在安装允许读取外部数据格式的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

我试图按照这里和这里以及其中列出的来源（例如，survival包中有关时间变量协变量的Therneau小册子），当然，当我复制和运行别人的代码时，一切都很顺利。但是我正在尝试使用一个数据集从头开始做到这一点，并将其结果与我的结果进行比较。但我就是无法让它起作用。

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

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 的书籍，现在最后一章证明是迄今为止我最大的挑战。我还有三十页要读完，任何帮助都将不胜感激。

- llewmills

Answer 1

我找出了我的错误。当我创建交互变量TREATINT时，犯了一个愚蠢的错误。应该是：

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

而不是：

longLOS$TREATINT <- longLOS$EVENT*(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模型交互的其他人。

- llewmills