Análisis jerárquico
ejemplo con datos simulados
objetivos:
- simular datos con estructura jerárquica
- realizar análisis de complete pooling y no-pooling
- comparar los análisis anteriores con un modelo jerárquico usando JAGS
Suponiendo que estudiamos 10 sub-poblaciones de una región y que seguimos la supervivencia de un número variable de individuos en cada sub-población. Vamos a simular datos con estructura jerárquica. La idea es que todas las sub-poblaciones son parte de una misma "metapoblación", son de la misma especie, etc. Esto hace que esperemos que los parámetros de supervivencia sean parecidos entre las sub-poblaciones, pero que tengan cierto grado de diferencia. Para modelar la variabilidad entre sub-poblaciones podemos usar una distribución Beta, por ejemplo con parámetros \(a_s = 2\) y \(b_s = 10\). Si repasamos el Bestiario vemos que esta combinación de parámetros resulta en un valor esperado de \(\frac{a_s}{(a_s+b_s)}\) = 0.1666667 y una varianza de \(\frac{(a_s \times b_s)}{(a_s + b_s)^2 \times (a_s + b_s + 1)}\) = 0.0106838. Entonces, cada sub-población tendrá su propia tasa de supervivencia. Para simular los datos hacemos:
set.seed(1234)
n <- 10 # subgrupos
m <- c(30, 28, 20, 30, 30, 26, 29, 5, 3, 27) # individuos muestraeados por grupo
a_s <- 2
b_s <- 10
theta <- rbeta(n, a_s, b_s) # generamos n tasas de mortalidad
y <- rbinom(n, size = m, prob = theta) # simulamos número de muertes por grupo
op <- par(mar = c(5, 6, 4, 5) + 0.1, mgp = c(3.5, 1, 0), cex.lab = 1.5, font.lab = 2,
cex.axis = 1.3, las = 1, bty = "n")
plot(table(y), xlab = "número de muertes", ylab = "frecuencia")
par(op)Podemos modelar estos datos asumiendo que: 1. todas las sub-poblaciones tienen igual mortalidad (complete pooling) 2. cada sub-población tiene una tasa de mortalidad independiente de las otras (no pooling) 3. cada sub-población es diferente pero parecida a las otras (partial pooling)
Vamos por partes:
(1) complete pooling
En este caso, como el modelo de datos (cuántos individuos mueren) es una Binomial, podemos usar una Beta como previa conjugada. Si queremos usar previas no-informativas para la tasa de mortalidad hacemos:alpha <- 1
beta <- 1
x <- seq(from = 0, to = 1, by = 0.001)
pos_theta1 <- dbeta(x, alpha + sum(y), beta + sum((m - y)))
op <- par(mar = c(5, 6, 4, 5) + 0.1, mgp = c(3.5, 1, 0), cex.lab = 1.5, font.lab = 2,
cex.axis = 1.3, las = 1, bty = "n")
plot(x, pos_theta1, type = "l", lty = 2, lwd = 3, ylab = "densidad", xlab = "tasa de mortalidad")
par(op)Preguntas:
1- ¿Cuál es el valor esperado de la posterior de la tasa de mortalidad asumiendo "complete pooling"? 2- ¿Qué tan diferente o parecida es esa posterior a la media de la supervivencia?
(2) no-pooling
De nuevo usando previas no informativas:
pos_theta2 <- matrix(0, n, length(x))
op <- par(mar = c(5, 6, 4, 5) + 0.1, mgp = c(3.5, 1, 0), cex.lab = 1.5, font.lab = 2,
cex.axis = 1.3, las = 1, bty = "n")
plot(x, pos_theta1, type = "l", lty = 2, lwd = 3, ylab = "densidad", xlab = "tasa de mortalidad")
for (i in 1:10) {
pos_theta2[i, ] <- dbeta(x, alpha + y[i], beta + m[i] - y[i])
lines(x, pos_theta2[i, ], lwd = 3, col = "darkgrey")
}
par(op)Preguntas:
3- ¿Qué se puede decir de las distintas posteriores?
4- ¿Cómo estimaríamos una tasa de mortalidad promedio para la meta-población a partir de estas posteriores?
(3) partial pooling
Finalmente, podemos hacer un análisis jerárquico (partial pooling), reconociendo explícitamente que cada sub-población es potencialmente distinta de las otras, pero sin embargo, todas son parte de la misma meta-población y por ende las observaciones en las sub-poblaciones no son del todo independientes.
Just Another Gibs Sampler
Para este análisis vamos a usar JAGS. JAGS es un software que programa cadenas de Markov chain Monte Carto (MCMC) para modelos bayesianos
JAGS es un sucesor de BUGS, que es Bayesian inference Using Gibbs Sampling (Lunn, Jackson, Best, Thomas, & Spiegelhalter, 2013; Lunn, Thomas, Best, & Spiegelhalter, 2000). JAGS es muy parecido a BUGS pero tiene algunas funciones extra y a veces es más rápido. Además, JAGS se puede usar en Windows, Mac y Linux.
Para usar JAGS (o BUGS) primero tenemos que definir el modelo usando un lenguaje sintético y guardarlo como un archivo de texto. Podemos escribir el modelo en un editor de texto o, como en el ejemplo de abajo, usar la función cat de R que nos permite escribir y guardar texto desde un script. En este caso, llamaremos al archivo "hier.bug":
cat(file="hier.bug",
"
model
{
for( i in 1 : n ) {
y[i] ~ dbin(theta[i], m[i])
theta[i] ~ dbeta(a,b)
}
a ~ dnorm(0,0.001)T(0,)
b ~ dnorm(0,0.001)T(0,)
mean_pobl <- a/(a+b)
}
")Luego definimos la lista de datos que vamos a pasarle a JAGS, una función para los valores iniciales de las cadenas Markovianas, y la lista de parámetros que queremos guardar.
data <- list("y", "m", "n")
inits <- function() list(a = runif(1, 1, 5), b = runif(1, 5, 20))
params <- c("a", "b", "theta", "mean_pobl")Ahora solo queda cargar los paquetes de R que sirven para comunicarse con JAGS, definir las iteraciones, número de cadenas a correr y cuántos valores descartar como burn-in.
library(jagsUI)
ni <- 5000
nc <- 3
nt <- 4
nb <- 2500
hier.sim <- jags(data, inits, params, model.file = "hier.bug", n.chains = nc,
n.iter = ni, n.burnin = nb, n.thin = nt)##
## Processing function input.......
##
## Done.
##
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph Size: 37
##
## Initializing model
##
## Adaptive phase, 100 iterations x 3 chains
## If no progress bar appears JAGS has decided not to adapt
##
##
## Burn-in phase, 2500 iterations x 3 chains
##
##
## Sampling from joint posterior, 2500 iterations x 3 chains
##
##
## Calculating statistics.......
##
## Done.
print(hier.sim)## JAGS output for model 'hier.bug', generated by jagsUI.
## Estimates based on 3 chains of 5000 iterations,
## burn-in = 2500 iterations and thin rate = 4,
## yielding 1875 total samples from the joint posterior.
## MCMC ran for 0.023 minutes at time 2015-09-22 14:03:40.
##
## mean sd 2.5% 50% 97.5% overlap0 f Rhat n.eff
## a 3.962 2.831 0.718 3.367 11.285 FALSE 1 1.076 49
## b 27.184 17.012 5.712 24.015 69.386 FALSE 1 1.072 45
## theta[1] 0.059 0.038 0.003 0.054 0.148 FALSE 1 1.017 120
## theta[2] 0.097 0.043 0.030 0.093 0.196 FALSE 1 1.003 701
## theta[3] 0.070 0.044 0.005 0.066 0.173 FALSE 1 1.007 308
## theta[4] 0.221 0.063 0.117 0.215 0.363 FALSE 1 1.006 352
## theta[5] 0.185 0.057 0.090 0.179 0.305 FALSE 1 1.006 338
## theta[6] 0.162 0.053 0.074 0.157 0.278 FALSE 1 1.000 1875
## theta[7] 0.152 0.051 0.071 0.147 0.268 FALSE 1 1.000 1875
## theta[8] 0.105 0.063 0.010 0.097 0.252 FALSE 1 1.005 365
## theta[9] 0.112 0.066 0.010 0.106 0.268 FALSE 1 1.001 1585
## theta[10] 0.100 0.044 0.030 0.094 0.198 FALSE 1 1.004 496
## mean_pobl 0.127 0.033 0.069 0.125 0.199 FALSE 1 1.001 1387
## deviance 33.528 5.169 24.354 33.271 44.425 FALSE 1 1.018 115
##
## Successful convergence based on Rhat values (all < 1.1).
## Rhat is the potential scale reduction factor (at convergence, Rhat=1).
## For each parameter, n.eff is a crude measure of effective sample size.
##
## overlap0 checks if 0 falls in the parameter's 95% credible interval.
## f is the proportion of the posterior with the same sign as the mean;
## i.e., our confidence that the parameter is positive or negative.
##
## DIC info: (pD = var(deviance)/2)
## pD = 13.1 and DIC = 46.668
## DIC is an estimate of expected predictive error (lower is better).
a.sim <- hier.sim$sims.list$a
b.sim <- hier.sim$sims.list$b
mean(hier.sim$sims.list$mean_pobl) # comparar con el valor usado para simular los datos ## [1] 0.1272574
alpha = 1
beta = 1
x = seq(from = 0, to = 1, by = 0.001)
op = par(mfrow = c(1, 2), mar = c(4, 4, 1, 1) + 0.1, las = 1, cex.lab = 1.5,
cex.axis = 1.3)
pos_theta1 = dbeta(x, alpha + sum(y), beta + sum((m - y)))
plot(x, pos_theta1, type = "l", lty = 2, lwd = 3, ylab = "densidad", xlab = "tasa de mortalidad",
ylim = c(0, 20))
pos_theta2 = matrix(0, 10, length(x))
for (i in 1:10) {
pos_theta2[i, ] = dbeta(x, alpha + y[i], beta + m[i] - y[i])
lines(x, pos_theta2[i, ], lwd = 3, col = "darkgrey")
}
plot(x, dbeta(x, mean(hier.sim$sims.list$a), mean(hier.sim$sims.list$b)), lwd = 2,
col = 2, xlab = "tasa de mortalidad", ylab = "", type = "l", ylim = c(0,
20))
lines(x, pos_theta1, type = "l", lty = 2, lwd = 3)
for (i in 1:10) {
lines(density(hier.sim$sims.list$theta[, i]), col = "blue", lwd = 2)
}
lines(density(hier.sim$sims.list$mean_pobl), lwd = 3, col = 2)
par(op)