Представление параметрической модели выживания в форме "Процесс подсчета" в JAGS
Эта проблема
Я пытаюсь построить модель выживания в JAGS, которая учитывает изменяющиеся во времени ковариаты. Я хотел бы, чтобы это была параметрическая модель - например, предполагая, что выживание следует распределению Вейбулла (но я бы хотел, чтобы опасность изменялась, поэтому экспонента слишком проста). Так что это, по сути, байесовская версия того, что можно сделать в пакете flexsurv, который допускает изменяющиеся во времени ковариаты в параметрических моделях.
Поэтому я хочу иметь возможность вводить данные в форме "процесса подсчета", где у каждого субъекта есть несколько строк, каждая из которых соответствует временному интервалу, в котором их ковариаты оставались постоянными (как описано в этом PDF-документе или здесь. Это (start, stop] формулировка, flexurv позволяют пакеты survival или flexurv.
К сожалению, каждое объяснение того, как проводить анализ выживаемости в JAGS, похоже, предполагает один ряд на субъекта.
Я попытался использовать этот более простой подход и распространить его на формат процесса подсчета, но модель неправильно оценивает распределение.
Неудачная попытка:
Вот пример. Сначала мы генерируем некоторые данные:
library('dplyr')
library('survival')
## Make the Data: -----
set.seed(3)
n_sub <- 1000
current_date <- 365*2
true_shape <- 2
true_scale <- 365
dat <- data_frame(person = 1:n_sub,
true_duration = rweibull(n = n_sub, shape = true_shape, scale = true_scale),
person_start_time = runif(n_sub, min= 0, max= true_scale*2),
person_censored = (person_start_time + true_duration) > current_date,
person_duration = ifelse(person_censored, current_date - person_start_time, true_duration)
)
person person_start_time person_censored person_duration
(int) (dbl) (lgl) (dbl)
1 1 11.81416 FALSE 487.4553
2 2 114.20900 FALSE 168.7674
3 3 75.34220 FALSE 356.6298
4 4 339.98225 FALSE 385.5119
5 5 389.23357 FALSE 259.9791
6 6 253.71067 FALSE 259.0032
7 7 419.52305 TRUE 310.4770
Затем мы разбиваем данные на 2 наблюдения по каждому предмету. Я просто делю каждый предмет на время = 300 (если они не успели на время = 300, в котором они получают только одно наблюдение).
## Split into multiple observations per person: --------
cens_point <- 300 # <----- try changing to 0 for no split; if so, model correctly estimates
dat_split <- dat %>%
group_by(person) %>%
do(data_frame(
split = ifelse(.$person_duration > cens_point, cens_point, .$person_duration),
START = c(0, split[1]),
END = c(split[1], .$person_duration),
TINTERVAL = c(split[1], .$person_duration - split[1]),
CENS = c(ifelse(.$person_duration > cens_point, 1, .$person_censored), .$person_censored), # <— edited original post here due to bug; but problem still present when fixing bug
TINTERVAL_CENS = ifelse(CENS, NA, TINTERVAL),
END_CENS = ifelse(CENS, NA, END)
)) %>%
filter(TINTERVAL != 0)
person split START END TINTERVAL CENS TINTERVAL_CENS
(int) (dbl) (dbl) (dbl) (dbl) (dbl) (dbl)
1 1 300.0000 0 300.0000 300.00000 1 NA
2 1 300.0000 300 487.4553 187.45530 0 187.45530
3 2 168.7674 0 168.7674 168.76738 1 NA
4 3 300.0000 0 300.0000 300.00000 1 NA
5 3 300.0000 300 356.6298 56.62979 0 56.62979
6 4 300.0000 0 300.0000 300.00000 1 NA
Теперь мы можем настроить модель JAGS.
## Set-Up JAGS Model -------
dat_jags <- as.list(dat_split)
dat_jags$N <- length(dat_jags$TINTERVAL)
inits <- replicate(n = 2, simplify = FALSE, expr = {
list(TINTERVAL_CENS = with(dat_jags, ifelse(CENS, TINTERVAL + 1, NA)),
END_CENS = with(dat_jags, ifelse(CENS, END + 1, NA)) )
})
model_string <-
"
model {
# set priors on reparameterized version, as suggested
# here: https://sourceforge.net/p/mcmc-jags/discussion/610036/thread/d5249e71/?limit=25#8c3b
log_a ~ dnorm(0, .001)
log(a) <- log_a
log_b ~ dnorm(0, .001)
log(b) <- log_b
nu <- a
lambda <- (1/b)^a
for (i in 1:N) {
# Estimate Subject-Durations:
CENS[i] ~ dinterval(TINTERVAL_CENS[i], TINTERVAL[i])
TINTERVAL_CENS[i] ~ dweibull( nu, lambda )
}
}
"
library('runjags')
param_monitors <- c('a', 'b', 'nu', 'lambda')
fit_jags <- run.jags(model = model_string,
burnin = 1000, sample = 1000,
monitor = param_monitors,
n.chains = 2, data = dat_jags, inits = inits)
# estimates:
fit_jags
# actual:
c(a=true_shape, b=true_scale)
В зависимости от того, где находится точка разделения, модель оценивает очень разные параметры для базового распределения. Он получает правильные параметры, только если данные не разбиты на форму подсчета. Кажется, что это не способ форматирования данных для такого рода проблемы.
Если я пропускаю предположение, и моя проблема в меньшей степени связана с JAGS и больше связана с тем, как я ее формулирую, предложения приветствуются. Я мог бы отчаиваться, что изменяющиеся во времени ковариаты не могут использоваться в параметрических моделях выживания (и могут использоваться только в моделях, таких как модель Кокса, которая предполагает постоянные опасности и которая фактически не оценивает лежащее в основе распределение) - однако, как Как я упоминал выше, пакет flexsurvreg в R действительно учитывает формулировку (start, stop] в параметрических моделях.
И если кто-нибудь знает, как построить такую модель на другом языке (например, STAN вместо JAGS), это тоже будет оценено.
В любом случае, я был озадачен этим более недели, поэтому любая помощь или предложения очень ценятся.
РЕДАКТИРОВАТЬ:
Крис Джексон дает несколько полезных советов по электронной почте:
Я думаю, что конструкция T() для усечения в JAGS необходима здесь. По существу, для каждого периода (t [i], t [i + 1]), где человек жив, но ковариата постоянна, время выживания сокращается влево в начале периода и, возможно, также подвергается цензуре вправо в конец. Так что вы бы написали что-то вроде
y[i] ~ dweib(shape, scale[i])T(t[i], )
Я попытался реализовать это предложение следующим образом:
model {
# same as before
log_a ~ dnorm(0, .01)
log(a) <- log_a
log_b ~ dnorm(0, .01)
log(b) <- log_b
nu <- a
lambda <- (1/b)^a
for (i in 1:N) {
# modified to include left-truncation
CENS[i] ~ dinterval(END_CENS[i], END[i])
END_CENS[i] ~ dweibull( nu, lambda )T(START[i],)
}
}
К сожалению, это не совсем подходит. Со старым кодом модель в основном правильно понимала параметр масштаба, но очень плохо справлялась с параметром формы. С этим новым кодом он становится очень близким к правильному параметру формы, но постоянно переоценивает параметр масштаба. Я заметил, что степень переоценки коррелирует с тем, насколько поздно наступает точка разделения. Если точка разделения ранняя (cens_point = 50), переоценка на самом деле отсутствует; если уже поздно (cens_point = 350), их много.
Я подумал, что, возможно, проблема может быть связана с "двойным подсчетом" наблюдений: если мы видим цензурированное наблюдение в момент времени t = 300, то от этого же человека, наблюдение без цензуры в момент времени t = 400, мне кажется интуитивным, что этот человек вносит две точки данных в наш вывод о параметрах Вейбулла, тогда как на самом деле они должны просто вносить одну точку. Поэтому я попытался включить случайный эффект для каждого человека; однако, это полностью провалилось, с огромными оценками (в диапазоне 50-90) для параметра nu. Я не уверен, почему это так, но, возможно, это вопрос для отдельного поста. Поскольку я не знаю, связаны ли проблемы, вы можете найти код для всего этого поста, включая код JAGS для этой модели, здесь.