sobre - principios de conteo pdf
Representación del modelo de supervivencia paramétrico en forma de ''Proceso de conteo'' en JAGS (0)
El problema
Estoy tratando de construir un modelo de supervivencia en JAGS que permita covariables variables en el tiempo. Me gustaría que fuera un modelo paramétrico, por ejemplo, asumiendo que la supervivencia sigue la distribución de Weibull (pero me gustaría permitir que el peligro varíe, por lo que la exponencial es demasiado simple). Entonces, esta es esencialmente una versión bayesiana de lo que se puede hacer en el paquete flexsurv , que permite covariables variables en el tiempo en modelos paramétricos.
Por lo tanto, quiero poder ingresar los datos en una forma de ''proceso de conteo'', donde cada sujeto tiene múltiples filas, cada una correspondiente a un intervalo de tiempo en el que sus covariables permanecieron constantes (como se describe en este pdf o here . Esta es la formulación (start, stop] que permiten los paquetes de survival o flexurv .
Desafortunadamente, cada explicación de cómo realizar el análisis de supervivencia en JAGS parece asumir una fila por sujeto.
Intenté adoptar este enfoque más simple y extenderlo al formato del proceso de conteo, pero el modelo no estima correctamente la distribución.
Un intento fallido:
Aquí hay un ejemplo. Primero, generamos algunos datos:
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
Luego, dividimos los datos en 2 observaciones por materia. Solo estoy dividiendo cada tema en el tiempo = 300 (a menos que no hayan llegado al tiempo = 300, en el que obtienen solo una observación).
## 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
Ahora podemos configurar el modelo 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)
Dependiendo de dónde esté el punto de división, el modelo estima parámetros muy diferentes para la distribución subyacente. Solo obtiene los parámetros correctos si los datos no se dividen en el formulario de proceso de conteo. Parece que esta no es la forma de formatear los datos para este tipo de problema.
Si me falta una suposición y mi problema está menos relacionado con JAGS y más relacionado con cómo estoy formulando el problema, las sugerencias son muy bienvenidas. Podría estar desesperado porque las covariables variables en el tiempo no pueden usarse en modelos de supervivencia paramétricos (y solo pueden usarse en modelos como el modelo cox, que asume riesgos constantes y que en realidad no estiman la distribución subyacente), sin embargo, como Como mencioné anteriormente, el paquete flexsurvreg en R admite la formulación (start, stop] en modelos paramétricos.
Y si alguien sabe cómo construir un modelo como este en otro idioma (por ejemplo, STAN en lugar de JAGS), eso también sería apreciado.
De cualquier manera, he estado perplejo por más de una semana, por lo que cualquier ayuda o sugerencia es muy apreciada.
EDITAR:
Chris Jackson proporciona algunos consejos útiles por correo electrónico:
Creo que la construcción T () para el truncamiento en JAGS es necesaria aquí. Esencialmente para cada período (t [i], t [i + 1]) donde una persona está viva pero la covariable es constante, el tiempo de supervivencia se trunca a la izquierda al comienzo del período, y posiblemente también se censura a la derecha en el fin. Así que escribirías algo como
y[i] ~ dweib(shape, scale[i])T(t[i], )
Intenté implementar esta sugerencia de la siguiente manera:
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],)
}
}
Desafortunadamente, esto no acaba de hacer el truco. Con el código anterior, el modelo estaba en su mayoría acertando el parámetro de escala, pero haciendo un trabajo muy malo en el parámetro de forma. Con este nuevo código, se acerca mucho al parámetro de forma correcta, pero sobreestima constantemente el parámetro de escala. Me he dado cuenta de que el grado de sobreestimación se correlaciona con la tardanza del punto de división. Si el punto de división es temprano ( cens_point = 50 ), no hay realmente una sobreestimación; si es tarde ( cens_point = 350 ), hay mucho.
Pensé que tal vez el problema podría estar relacionado con el "doble conteo" de las observaciones: si vemos una observación censurada en t = 300, entonces, desde esa misma persona, una observación sin censura en t = 400, me parece intuitivo que esta persona está aportando dos puntos de datos a nuestra inferencia sobre los parámetros de Weibull cuando en realidad deberían estar contribuyendo solo un punto. Por lo tanto, traté de incorporar un efecto aleatorio para cada persona; sin embargo, esto fracasó completamente, con estimaciones enormes (en el rango de 50-90) para el parámetro nu . No estoy seguro de por qué, pero tal vez sea una pregunta para una publicación separada. Dado que no sé si los problemas están relacionados, puede encontrar el código para toda esta publicación, incluido el código JAGS para ese modelo, here .