¿Cómo calculo la distribución predictiva posterior en WinBUGS?

Question

Me gustaría trabajar en la distribución predictiva posterior del siguiente modelo multinivel (efectos mixtos) en WinBUGS. El ejemplo es tomado de un ejemplo que encontré en línea con fines ilustrativos.

Yi ~ Binomial(pi,1) where

logit(pi) = b0+ b1 log(income) + b2 distance + b3 dropout + b4 college + uj(i)

Se dan priors no informativos para los efectos fijos, asumiendo bk ~ Normal(0, 0.000001). El segundo parámetro es la precisión (el recíproco de la varianza), por lo que la varianza es un millón. Asumimos que

uj(i) ~ N(0, t)

donde la precisión t tiene una prior gamma con parámetros 0.001 y 0.001, por lo que la media es uno y la varianza es 1000.

En WinBUGS:

model {
# N observaciones
for (i in 1:N) {
    hospital[i] ~ dbern(p[i])
    logit(p[i]) <- bcons + blonginc*loginc[i] + bdistance*distance[i] + 
        bdropout*dropout[i] + bcollege*college[i] + u[group[i]] 
}
# M grupos
for (j in 1:M) {
    u[j] ~ dnorm(0,tau)
}
# Priors
bcons     ~ dnorm(0.0,1.0E-6)
bloginc   ~ dnorm(0.0,1.0E-6)
bdistance ~ dnorm(0.0,1.0E-6)
bdropout  ~ dnorm(0.0,1.0E-6)
bcollege  ~ dnorm(0.0,1.0E-6)
# Hyperprior
tau ~ dgamma(0.001,0.001)
}

Puedo ajustar este modelo usando stan_glmer que automáticamente me da el promedio de muestra de la distribución predictiva posterior mean_ppd. Me gustaría hacer lo mismo en WinBUGS pero no sé cómo

Answer 1

Para monitorear un PPD, simplemente escribe la distribución de la variable que deseas observar. Esa variable puede ser un nodo observado, o una función del nodo observado y otras variables en tu modelo.

Por ejemplo, digamos que deseas ver si el nodo "hospital[1]" se ajusta bien. Necesitas una línea adicional de código como esta:

hospital_1_ppd ~ dbern(p[1])

, y agrega 'hospital_1_ppd' a la lista de variables que estás monitoreando. Puedes averiguar el resto.

Espero que ayude.