Наборы данных процесса подсчета для моделей непропорциональной опасности (Кокса) с переменными взаимодействия

Я пытаюсь запустить модель непропорциональной регрессии Кокса с переменной взаимодействия со временем, как описано в главе 15 (раздел 15.3) Прикладной лонгитюдный анализ данных Зингера и Виллетта. Однако я не могу получить ответы, согласующиеся с книгой.

Данные, используемые в этой книге, и исходный код предоставлены в этом фантастический сайт. К сожалению, в последней главе отсутствует код R, а предоставленный набор данных для R для примера, обсуждаемого в тексте, является неполным и дает неправильные ответы для простейшей модели (которую я, делать, знаю, как запустить). Вместо этого, чтобы получить полный набор данных для этого примера, нужно щелкнуть ссылку «Загрузить» в столбце «SAS» (который имеет правильный набор данных), а затем, после установки пакета haven (который позволяет читать в сторонних форматах данных) ), прочтите соответствующий набор данных через:

haven::read_sas("alda/lengthofstay.sas7bdat")

В этом наборе данных указывается продолжительность пребывания участников (переменная ID) (переменная DAYS) при стационарном лечении в больнице. Переменная цензуры - CENSOR. Исследователи выдвинули гипотезу, что два разных типа лечения (бинарная переменная TREAT) могут предсказать дифференциальные значения риска отказа от лечения. Кроме того, они ожидали, что разница в рисках между группами нет будет постоянной с течением времени, поэтому требуется создание условия взаимодействия. Я могу заставить работать простую модель основного эффекта, возвращая те же коэффициенты опасности, которые указаны в книге (именно так я в конечном итоге обнаружил, что файл .csv, поставляемый с кодом R, был неполным).

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), и, конечно же, когда я копирую чужой код и запускаю его. все работает нормально. Но я пытаюсь сделать это для себя с нуля с помощью набора данных, результаты которого я могу сравнить со своими. И я просто не могу заставить его работать.

сначала я создал переменную EVENT

los$EVENT <- 1 - los$CENSOR

в наборе данных есть повторяющийся идентификационный номер, который вызывает проблемы. Поэтому нам нужно изменить его на новый идентификационный номер.

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.

Что я делаю неправильно? Я медленно прорабатывал Зингера и Уиллетта в течение почти трех лет (я начал, еще будучи аспирантом), и теперь последняя глава оказывается, безусловно, моей самой большой проблемой. У меня осталось тридцать страниц; любая помощь будет невероятно оценена.