Análisis Bayesianos
Los análisis Bayesianos son similares a los que vimos en máxima verosimilitud en el sentido que dependen explícitamente de modelos probabilísticos para los datos. Es decir, tenemos que definir un modelo de datos. La gran diferencia es que con Bayes, podemos obtener distribuciones de probabilidades para todas las cantidades no observadas, incluyendo parámetros, valores perdidos o nuevas posibles observaciones. De esta manera, los análisis Bayesianos nos permiten cuantificar incertidumbre y armar modelos realistas que tienen en cuenta por ejemplo observaciones imperfectas.
Como vimos en la teórica, la regla de Bayes planteada en términos de datos y parámetros es:
\[ p(\boldsymbol{\theta} \lvert \boldsymbol{y}) = \frac{p(\boldsymbol{y} \lvert \boldsymbol{\theta}) p(\boldsymbol{\theta)}}{\int p(\boldsymbol{y} \lvert \boldsymbol{\theta)} d \boldsymbol{\theta} } \]
Es decir que la probabilidad posterior de los parámetros \(\boldsymbol{\theta}\) dado que observamos los datos \(\mathbf{y}\) es igual al likelihood multiplicado por las previas y dividido por la probabilidad total de los datos. La función de likelihood nos da la probabilidad de observar los datos condicional al valor de los parámetros \(p(\mathbf{y} \lvert \boldsymbol{\theta})\). La previa de los parámetros \(p(\boldsymbol{\theta})\) refleja los posibles valores de los parámetros de acuerdo con nuestras “creencias” previas, o los resultados de estudios anteriores, o lo que nos parece que tiene sentido para el sistema de estudio (en definitiva, en base a información previa). Finalmente, la probabilidad total de los datos se obtiene integrando la función de lilkelihood sobre los posibles valores de los parámetros. Como veremos más adelante, los análisis Bayesianos combinados con métodos numéricos permiten analizar modelos con muchos parámetros, niveles de variabilidad y variables “ocultas” pero primero vamos a empezar por casos simples donde podemos calcular las posteriores directamente.
Para entender bien cómo es todo el proceso, vamos a simular los datos. Supongamos que queremos estudiar la remoción de frutos en \(10\) plantas con \(100\) frutos cada una y suponemos que un buen modelo para este tipo de datos es una distribución Binomial con una probabilidad de éxito fija:
set.seed(1234)
nobs <- 10 # número de observaciones (plantas)
frutos <- rep(100, nobs) # frutos disponibles
p.rem <- 0.2 # probabilidad de remoción por fruto
removidos <- rbinom(nobs, size = frutos, prob = p.rem)En este caso, como el modelo de datos (cuántos frutos son removidos) es una Binomial, si usamos una distribución Beta como previa para el parámetro de probabilidad de éxito (p.rem) es posible obtener un resultado analítico para la posterior. En este caso, la posterior es otra distribución Beta pero con sus parámetros actualizados en base a las observaciones. Se dice entonces que la distribución Beta es la conjugada de la Binomial. Si la previa de la tasa de remoción por fruto es una distribución Beta con parámetros \(\alpha\) y \(\beta\), actualizamos los valores de \(\alpha\) y \(\beta\) en base a la cantidad de éxitos y fracasos obervados. La posterior de la tasa de remoción por fruto es entonces una Beta con \(\alpha = \sum y\), \(\beta = \sum (n-y)\) donde \(y\) representa a los frutos removidos de los \(n\) disponibles. Veamos como hacer esto en R.
# previas
alpha <- 1
beta <- 1
pos.alpha <- alpha + sum(removidos)
pos.beta <- beta + sum(frutos - removidos)
# valor esperado de la posterior
pos.alpha/(pos.alpha + pos.beta)## [1] 0.1936128# quantiles de la posterior
qbeta(c(0.025, 0.975), pos.alpha, pos.beta)## [1] 0.1697490 0.2186349# ahora un gráfico
op <- par(cex.lab = 1.5, font.lab = 2, cex.axis = 1.3, las = 1, bty = "n")
curve(dbeta(x, alpha + sum(removidos), beta + sum(frutos - removidos)), lwd = 2,
ylab = "Densidad de probabilidad", xlab = "Probabilidad de remoción")
curve(dbeta(x, 1, 1), lwd = 2, col = "gray", add = TRUE)
text(0.6, 2.5, "previa")
text(0.35, 25, "posterior")
par(op)Ahora veamos como ajustar este modelo con JAGS
Just Another Gibbs Sampler
JAGS es un software que programa cadenas de Markov Chain Monte Carlo (MCMC) para modelos bayesianos (Plummer, M. (2003). JAGS: A program for analysis of Bayesian graphical models using Gibbs sampling. In Proceedings of the 3rd international workshop on distributed statistical computing (dsc 2003), Vienna, Austria.ISSN 1609-395X. Plummer, M. (2015). JAGS version 4.0 user manual).
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.
cat(file = "rem.bug", "
model{
# likelihood
for (i in 1:nobs) {
removidos[i] ~ dbin(theta, frutos[i])
}
# previa
theta ~ dbeta(1, 1)
}")Además, tenemos que armar una lista con los datos que vamos a pasarle a JAGS, una función para generar los valores iniciales de las cadenas Markovianas, y definir los parámetros que queremos guardar.
m.data <- list("nobs", "frutos", "removidos")
inits <- function() list(theta = runif(1, 0, 1))
params <- c("theta")Finalmente, definimos cuántas iteraciones queremos correr por cadena, cuántos valores vamos a ir descartado en una secuencia (“thin”), qué largo tiene el “burn in” y cuántas cadenas queremos correr.
ni <- 1000 # número de iteraciones
nt <- 1 # tasa de 'thining'
nb <- 500 # cuantas iteraciones usamos de 'burn in'
nc <- 3 # y cuantas cadenas corremosPara llamar a JAGS desde R usamos el paquete jagsUI y la función jags de ese paquete
library(jagsUI)
m.sim <- jags(data = m.data, inits = inits, parameters.to.save = params, model.file = "rem.bug",
n.chains = nc, n.thin = nt, n.iter = ni, n.burnin = nb)##
## Processing function input.......
##
## Done.
##
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 10
## Unobserved stochastic nodes: 1
## Total graph size: 24
##
## Initializing model
##
## Adaptive phase, 100 iterations x 3 chains
## If no progress bar appears JAGS has decided not to adapt
##
##
## Burn-in phase, 500 iterations x 3 chains
##
##
## Sampling from joint posterior, 500 iterations x 3 chains
##
##
## Calculating statistics.......
##
## Done.Para ver los resultados de este análisis podemos pedirle a R que “imprima” la salida de la función jags usando print(m.sim) en este caso. Vemos que se reporta el nombre del modelo que se usó, cuántas cadenas se simularon y otros detalles como el tiempo de ejecución. Después aparece una tabla con la lista de parámetros que le pedimos que registre (en este caso uno solo) y la devianza. En la tabla aparece la media, desvío y cuantiles estimados a partir de las cadenas Markovianas. También aparece una columna (overlap0) que nos dice si la posterior incluye al cero o no, y otra (f) que nos dice qué fracción de la posterior es del mismo signo que la media. Finalmente, aparecen dos columnas con información importante. Rhat estima si las cadenas convergieron a una distribución estable y n.eff estima el número efectivo de muestras de la posterior que surgen de las cadenas. Antes que hacer nada con la salida de JAGS, tenemos que chequear que las cadenas hayan convergido (Rhat \(\leq 1.1\)). También, ver que el tamaño efectivo de la muestra de la posterior (n.eff) sea suficiente.
print(m.sim)## JAGS output for model 'rem.bug', generated by jagsUI.
## Estimates based on 3 chains of 1000 iterations,
## burn-in = 500 iterations and thin rate = 1,
## yielding 1500 total samples from the joint posterior.
## MCMC ran for 0.002 minutes at time 2017-06-19 21:38:16.
##
## mean sd 2.5% 50% 97.5% overlap0 f Rhat n.eff
## theta 0.193 0.013 0.169 0.193 0.220 FALSE 1 1.005 401
## deviance 56.010 1.529 54.925 55.425 60.334 FALSE 1 1.006 1500
##
## 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 = 1.2 and DIC = 57.179
## DIC is an estimate of expected predictive error (lower is better).¿Qué usamos para comparar el valor estimado de la probabilidad de remoción que estima JAGS con el resultado analítico que vimos antes?
En el objeto de salida de jags quedan guardadas varias cosas. Para ver qué cosas, podemos hacer
str(m.sim)## List of 24
## $ sims.list :List of 2
## ..$ theta : num [1:1500] 0.173 0.181 0.183 0.211 0.179 ...
## ..$ deviance: num [1:1500] 57.5 56 55.6 56.9 56.2 ...
## $ mean :List of 2
## ..$ theta : num 0.193
## ..$ deviance: num 56
## $ sd :List of 2
## ..$ theta : num 0.013
## ..$ deviance: num 1.53
## $ q2.5 :List of 2
## ..$ theta : num 0.169
## ..$ deviance: num 54.9
## $ q25 :List of 2
## ..$ theta : num 0.184
## ..$ deviance: num 55
## $ q50 :List of 2
## ..$ theta : num 0.193
## ..$ deviance: num 55.4
## $ q75 :List of 2
## ..$ theta : num 0.202
## ..$ deviance: num 56.4
## $ q97.5 :List of 2
## ..$ theta : num 0.22
## ..$ deviance: num 60.3
## $ overlap0 :List of 2
## ..$ theta : logi FALSE
## ..$ deviance: logi FALSE
## $ f :List of 2
## ..$ theta : num 1
## ..$ deviance: num 1
## $ Rhat :List of 2
## ..$ theta : num 1
## ..$ deviance: num 1.01
## $ n.eff :List of 2
## ..$ theta : num 401
## ..$ deviance: num 1500
## $ pD : num 1.17
## $ DIC : num 57.2
## $ summary : num [1:2, 1:11] 0.193 56.01 0.013 1.529 0.169 ...
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : chr [1:2] "theta" "deviance"
## .. ..$ : chr [1:11] "mean" "sd" "2.5%" "25%" ...
## $ samples :List of 3
## ..$ : mcmc [1:500, 1:2] 0.173 0.181 0.183 0.211 0.179 ...
## .. ..- attr(*, "dimnames")=List of 2
## .. .. ..$ : NULL
## .. .. ..$ : chr [1:2] "theta" "deviance"
## .. ..- attr(*, "mcpar")= num [1:3] 501 1000 1
## ..$ : mcmc [1:500, 1:2] 0.19 0.174 0.194 0.186 0.173 ...
## .. ..- attr(*, "dimnames")=List of 2
## .. .. ..$ : NULL
## .. .. ..$ : chr [1:2] "theta" "deviance"
## .. ..- attr(*, "mcpar")= num [1:3] 501 1000 1
## ..$ : mcmc [1:500, 1:2] 0.185 0.176 0.196 0.179 0.199 ...
## .. ..- attr(*, "dimnames")=List of 2
## .. .. ..$ : NULL
## .. .. ..$ : chr [1:2] "theta" "deviance"
## .. ..- attr(*, "mcpar")= num [1:3] 501 1000 1
## ..- attr(*, "class")= chr "mcmc.list"
## $ modfile : chr "rem.bug"
## $ model :List of 8
## ..$ ptr :function ()
## ..$ data :function ()
## ..$ model :function ()
## ..$ state :function (internal = FALSE)
## ..$ nchain :function ()
## ..$ iter :function ()
## ..$ sync :function ()
## ..$ recompile:function ()
## ..- attr(*, "class")= chr "jags"
## $ parameters : chr [1:2] "theta" "deviance"
## $ mcmc.info :List of 8
## ..$ n.chains : num 3
## ..$ n.adapt : num 100
## ..$ n.iter : num 1000
## ..$ n.burnin : num 500
## ..$ n.thin : num 1
## ..$ n.samples : num 1500
## ..$ end.values :List of 3
## .. ..$ : Named num [1:2] 58.294 0.171
## .. .. ..- attr(*, "names")= chr [1:2] "deviance" "theta"
## .. ..$ : Named num [1:2] 54.924 0.193
## .. .. ..- attr(*, "names")= chr [1:2] "deviance" "theta"
## .. ..$ : Named num [1:2] 55.155 0.187
## .. .. ..- attr(*, "names")= chr [1:2] "deviance" "theta"
## ..$ elapsed.mins: num 0.002
## $ run.date : POSIXct[1:1], format: "2017-06-19 21:38:16"
## $ random.seed: int 1497919096
## $ parallel : logi FALSE
## $ bugs.format: logi FALSE
## - attr(*, "class")= chr "jagsUI"Podemos hacer un gráfico de las cadenas Markovianas del parámetro de probabilidad de remoción
jagsUI::traceplot(m.sim, parameters = c("theta"))
Podemos graficar la posterior:
op <- par(cex.lab = 1.5, font.lab = 2, cex.axis = 1.3, las = 1, bty = "n")
plot(density(m.sim$sims.list$theta), lwd = 2, main = "", ylim = c(0, 35), xlab = expression(theta))
par(op)y comparar con el resultado teórico:
op <- par(cex.lab = 1.5, font.lab = 2, cex.axis = 1.3, las = 1, bty = "n")
plot(density(m.sim$sims.list$theta), lwd = 2, main = "", ylim = c(0, 35), xlab = expression(theta))
curve(dbeta(x, pos.alpha, pos.beta), add = TRUE, lwd = 2, lty = 2)
par(op)Si queremos estimar la probabilidad de que la tasa de remoción por fruto sea menor al \(16\)% hacemos:
length(which(m.sim$sims.list$theta < 0.16))/length(m.sim$sims.list$theta)## [1] 0.002y la probabilidad de que la tasa de remoción por fruto sea mayor que 20%:
length(which(m.sim$sims.list$theta > 0.2))/length(m.sim$sims.list$theta)## [1] 0.2906667También podemos calcular el intervalo de credibilidad:
library(coda)
HPDinterval(as.mcmc(m.sim$sims.list$theta))## lower upper
## var1 0.1671327 0.2177921
## attr(,"Probability")
## [1] 0.95A lo largo de otros prácticos iremos viendo cómo aprovechar la capacidad de generar muestras de las posteriores de nuestros modelos.
Juan Manuel Morales. 26 de Septiembre del 2015. última actualización: 2017-06-19