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