Comenzaré diciendo gracias por leer esto, y cualquier comentarista debería sentirse libre de criticar cualquier aspecto de esta especificación de modelado o flujo de trabajo que considere ingenuo o desinformado.
Estoy simulando dos grupos a partir de este proceso de curva doble sigmoidea:
$$ y = \frac{y_{\text{mid}}}{1 + e^{-(a \cdot \text{time} + b)}} + \frac{y_{\text{max}}}{1 + e^{-(c \cdot \text{time} + d)}}, \quad 0 \leq y \leq 100, \quad \text{time} \in [0, 300] $$
Usando el siguiente código en brms:
sw_dbl_log <- function(a, b, c, d, time){
((50)/(1+exp(-a*Time + b)))+((50)/(1+exp(-c*Time + d)))
}
nn = 300
Raw<- tibble(n = 1:(nn*2),
Group = rep(c("Control", "Experimental"), each = nn),
a = c(rnorm(nn, 0.1, 0.01),
rnorm(nn, 0.1, 0.01)),
b = c(rnorm(nn, 5, 0.01),
rnorm(nn, 5, 0.01)),
c = c(rnorm(nn, 0.1, 0.01),
rnorm(nn, 0.1, 0.002)),
d = c(rnorm(nn, 5, 0.01),
rnorm(nn, 25, 0.01))
) %>%
tidyr::expand(nesting(n,Group, a, b, c, d),
Time = seq(0, 300, by = 10)) %>%
mutate(y = sw_dbl_log(a, b, c, d, Time))Como se esperaba, esto produce distribuciones gaussianas para los parámetros de cada grupo (a:d), con la serie temporal del grupo de Control que se asemeja a una sola sigmoide con una fase de crecimiento. La serie temporal del grupo Experimental tiene dos fases de crecimiento, como se pretendía. Estoy aprovechando esta especificación porque cada parámetro da un resultado claramente interpretable, es decir, a y c representan tasas de crecimiento y b representa el punto de inflexión de cada una de las fases de crecimiento.
Mi procedimiento de ajuste es el siguiente:
formula <- bf(y ~ (50/(1+exp(-a*Time + b))) + (50/(1+exp(-c*Time + d))),
a + b + c + d ~ 0 + Group,
set_rescor(FALSE),
nl = TRUE)
prior1 <- c(
prior(normal(0.1, 0.1), nlpar = "a", lb = 0),
prior(normal(5, 2), nlpar = "b", lb = 0, ub = 10),
prior(normal(0.1, 0.1), nlpar = "c",lb = 0),
prior(normal(25, 5), nlpar = "d", coef = "GroupExperimental"),
prior(normal(5, 2), nlpar = "d", coef = "GroupControl")
)
fit <- brm(formula,
data = Raw,
prior = prior1,
control = list(adapt_delta = 0.95, max_treedepth = 15),
iter = 10000,
warmup = 5000,
cores = 3,
chains = 3)He omitido el ajuste previo y las comprobaciones predictivas del código anterior, pero la inspección visual se ve bien.
El problema que tengo actualmente es que el ajuste del modelo es pobre, es decir, los números de Rhat y ESS son demasiado bajos, y se debe principalmente a la bimodalidad en la distribución posterior de los parámetros para uno de los grupos, a veces Control y a veces Experimental, que hace que las cadenas no converjan. Por lo tanto, para algunos ajustes, cualquiera o todas las distribuciones posteriores de a, b, c, d por Grupo (por ejemplo, "b_a_GroupExperimental") serán bimodales.
No puedo determinar por qué el modelo se ajusta mal y causa esta bimodalidad. Mi primera suposición es que los priors están causando problemas, pero, según mi entendimiento ingenuo, los priors están especificados bastante ajustados alrededor de las distribuciones originales de simulación, por lo que no puedo ver por qué podrían estar causando el problema. Debo admitir que esta conclusión podría exponer mi ignorancia sobre el tema, lo que está alimentando mi incapacidad para solucionar los problemas de ajuste.
¿Debo cambiar los priors, el modelo, especificar valores iniciales o simplemente ejecutar más iteraciones? Cualquier ayuda sería muy apreciada.
Editar: Como actualización, aquí están los posteriores de cada parámetro por Grupo.
Además, si cambio los valores de simulación del grupo Control para b y d, y así hacer que la serie temporal represente una sigmoide doble en lugar de una única sigmoide, sus posteriores retornan como se esperaba. Por ejemplo, si $$ a_{C} = 0.1, b_{C} = 5, c_{C} = 0.1, d_{C} = 20 $$ , esos valores coinciden con la media de los posteriores y las cadenas se mezclan como se espera. Sin embargo, los del grupo Experimental siguen siendo bimodales y las cadenas no se mezclan.
