Ошибки сегментации в регрессии Кокса R с точными связями

Я пытаюсь подогнать большую дискретную модель пропорциональных опасностей (~ 100 тыс. строк, ~ 10 тыс. событий). Для этого я использовал coxph(..., method = "exact"), как рекомендовано в документации пакета выживания. документация, в которой указано:

«Точное частичное правдоподобие» эквивалентно условной логистической модели и подходит, когда время представляет собой небольшой набор дискретных значений. Если имеется большое количество ничьих и (старт, стоп) данных о выживании стиля, время вычислений будет чрезмерным.

Были некоторые предупреждения о вычислительной сложности с coxph и большом количестве связей, но согласно документации для clogit в том же пакете:

Однако вычисление точного частичного правдоподобия может быть очень медленным. Если в определенной страте было, скажем, 10 событий из 20 предметов, мы должны сложить знаменатель, включающий все возможные способы выбора 10 из 20, что составляет 20!/(10! 10!) = 184756 членов. Гейл и др. описывают метод быстрой рекурсии, который в значительной степени улучшает это; он был включен в версию 2.36-11 пакета survival.

Так что я не ожидал, что вычислительные проблемы будут слишком плохими. Тем не менее, я столкнулся со многими ошибками сегментации, когда пытался подогнать варианты тривиальной (с одним предиктором) модели Кокса к моему набору данных. Одним из них является «переполнение стека C», что приводит к короткому и приятному (и неинформативному) сообщению:

Error: segfault from C stack overflow
Execution halted

Другая ошибка «память не отображается», которая произошла, когда я случайно перевернул логическое значение «событие», так что у меня было ~ 90 000 событий вместо ~ 10 000:

 *** caught segfault ***
address 0xffffffffac577830, cause 'memory not mapped'

Traceback:
 1: fitter(X, Y, strats, offset, init, control, weights = weights,     method = method, row.names(mf))
 2: coxph(Surv(time, status == EVENT.STATUS) ~ litter, data = data,     method = "exact")
aborting ...

Для справки, код, который я запускаю, просто coxph(Surv(t, d) ~ x, data = data, method = 'exact'). t — столбец целых чисел, d — логическое значение, а x — число с плавающей запятой.

Это известные проблемы? Существуют ли обходные пути?

РЕДАКТИРОВАТЬ: Вот некоторый код, воспроизводящий проблему в наборе данных rats (реплицируется 1000 раз):

library(survival)
print("constructing data")
data <- rats
SIZE <- nrow(rats)
# passes with 100 reps, but fails with 100 on my machine (MacBook Pro, 16g RAM)
REPS <- 1000
# set to 0 for "memory not mapped", 1 for "C stack overflow"
EVENT.STATUS <- 0
data <- data[rep(seq_len(SIZE), REPS), ]
print(summary(data$status == EVENT.STATUS))
print("fitting model")
fit <- coxph(Surv(time, status == EVENT.STATUS) ~ litter,
             data = data, method = "exact")

А вот version:

platform       x86_64-apple-darwin14.0.0
arch           x86_64
os             darwin14.0.0
system         x86_64, darwin14.0.0
status
major          3
minor          1.2
year           2014
month          10
day            31
svn rev        66913
language       R
version.string R version 3.1.2 (2014-10-31)
nickname       Pumpkin Helmet

person Ben Kuhn    schedule 13.03.2015    source источник
comment
можете ли вы воспроизвести любой из наборов данных, включенных в survival ? (используйте `data(package='survival')`, чтобы перечислить их)   -  person jaimedash    schedule 13.03.2015
comment
Segfaults должен предложить вам собрать воспроизводимый пример и отправить его (или их) Терри Терно. Насколько мне известно, он посещает SO нечасто.   -  person IRTFM    schedule 13.03.2015
comment
@BondedDust: спасибо за напоминание. Я также хотел бы знать, есть ли какие-либо возможные обходные пути, пока ошибка не будет исправлена.   -  person Ben Kuhn    schedule 13.03.2015
comment
Я думал об использовании try, но это не сработает с segfaults. Вы уверены, что вам нужно использовать точное? Другими моделями пропорциональных опасностей, которые, возможно, больше подходят для сгруппированных данных, являются модели Пуассона-Глма.   -  person IRTFM    schedule 13.03.2015
comment
Да. Данные о событиях дискретны, поэтому имеют большое количество связей (всего ~ 35 уникальных событий в 100 000 строк). Мне дали понять, что приближения Эфрона и Бреслоу существенно расходятся с точным методом в этом случае.   -  person Ben Kuhn    schedule 13.03.2015
comment
try также не будет полезен, даже если он сработает для segfault, потому что мне нужно фактически обучить модель, а не просто восстановиться после ошибки.   -  person Ben Kuhn    schedule 13.03.2015


Ответы (1)


Я могу сделать модель Пуассона с этим набором данных. (У меня есть большой набор данных, на котором я не хочу рисковать вероятным segfault.)

fit <- glm(  I(status == 0) ~ litter +offset(log(time)), 
               data = data, family=poisson)

> fit

Call:  glm(formula = I(status == 0) ~ litter + offset(log(time)), family = poisson, 
    data = data)

Coefficients:
(Intercept)       litter  
  -4.706485    -0.003883  

Degrees of Freedom: 149999 Total (i.e. Null);  149998 Residual
Null Deviance:      60500 
Residual Deviance: 60150    AIC: 280200

Эта оценка эффекта litter должна быть аналогична той, которую вы получили бы из модели Cox PH.

Если вы хотите увидеть задокументированный «трюк со смещением», обратитесь к классической монографии Бреслоу и Дэя: «Статистические методы в исследованиях рака; Том II — Дизайн и анализ когортных исследований». Они использовали программный пакет GLIM, но его код очень похож на реализацию R glm, поэтому передача концепций должна быть простой. (У меня была возможность ненадолго поработать с Нормом Бреслоу над моей магистерской диссертацией с использованием GLIM. Он был великолепен. Я думаю, что мое предыдущее обучение GLIM сделало изучение R таким легким.)

person IRTFM    schedule 13.03.2015
comment
О, я раньше не видел трюка со смещением Пуассона. Это умно! Но разве эта модель не предполагает постоянной опасности? - person Ben Kuhn; 14.03.2015
comment
Да, это так, но coxph также предполагает постоянную относительную опасность. - person IRTFM; 14.03.2015
comment
Извините, я имел в виду, что, похоже, регрессия Пуассона предполагает, что базовая функция риска постоянна, а не просто постоянная относительная опасность (тогда как регрессия Кокса не делает никаких предположений о форме базовой опасности, хотя он предполагает постоянное отношение рисков между любыми двумя точками данных). - person Ben Kuhn; 14.03.2015
comment
Верно. В этом разница в методах. Вы можете построить модели с кусочно-переменными опасностями, но вам нужно будет разделить время и события на интервалы. - person IRTFM; 14.03.2015
comment
Что ж, профессор Тернау ответил на мой отчет об ошибке, заявив, что такое количество связей было бы вычислительно трудновыполнимым даже с лучшим кодом survival, а ошибка заключалась в том, что он пытался выделить больше памяти, чем было на моей машине. Он рекомендовал мне использовать приближение Эфрона, несмотря на большое количество связей, поэтому я отмечаю этот ответ как правильный по духу (дух не использует точную модель). - person Ben Kuhn; 14.03.2015