Análisis Bayesianos
Previas Conjugadas
Simulemos datos de remoción de frutos en 10 plantas con 100 frutos cada una usando una Binomial:
set.seed(1234)
nobs <- 10 # número de observaciones
frutos <- 100 # frutos disponibles
p.rem <- 0.2 # probabilidad de remoción por fruto
removidos <- rbinom(nobs, size = frutos, prob = p.rem)La idea principal detrás de los análisis Bayesianos es poder decir algo acerca de las variables aleatorias no observadas en base a las que sí observamos. Para llegar a ese objetivo, usamos reglas de probabilidades. Podemos usar la regla de Bayes en términos de los parámetros (θ) y los datos (y): p(θ|y) = p(y|θ) × p(θ)/p(y). En este caso, como el "modelo de datos" es una Binomial con una probabilidad de éxito fija, logramos una solución analítica si usamos como previa para θ una distribución Beta (ver el Bestiario).
La distribución Beta tiene dos parámetros (α y β) que se interpretan como la cantidad éxitos y fracasos observados previamente + 1. Entonces, si consideramos que hasta ahora no hemos visto nada (cero éxitos y cero fracasos), usamos como previa una Beta con α = 1 y β = 1. Esta opción corresponde a una previa "chata":
alpha <- 1
beta <- 1Para obetener la posterior de la probabilidad de remoción por fruto, actualizamos los parámetros de la Beta en base a la cantidad de frutos removidos (éxitos) y la cantidad de frutos que no fueron removidos (fracasos). Podemos ver por ejemplo, cómo queda la posetior luego de ver un sólo dato y luego de ver todos los datos:
xvec <- seq(0, 1, by = 0.001)
op <- par(mfrow = c(1, 2))
plot(xvec, dbeta(xvec, alpha + removidos[1], beta + frutos - removidos[1]),
type = "l", lwd = 2, col = "gray", xlab = expression(paste("valor de ",
theta)), ylab = "Densidad de Probabilidad", main = "n=1")
lines(xvec, dbeta(xvec, alpha, beta))
abline(v = removidos[1]/frutos, lwd = 2, lty = 2)
text(0.3, 5, "posterior", col = "darkgray", adj = 0)
text(0.25, 7.5, "dato")
text(0.6, 1.5, "previa")
plot(xvec, dbeta(xvec, alpha + sum(removidos), beta + sum(frutos - removidos)),
type = "l", lwd = 2, col = "gray", xlab = expression(paste("valor de ",
theta)), ylab = "Densidad de Probabilidad", main = "n=10")
lines(xvec, dbeta(xvec, alpha, beta))
abline(v = sum(removidos)/(frutos * nobs), lwd = 2, lty = 2)
text(0.3, 10, "posterior", col = "darkgray", adj = 0)
text(0.25, 25, "dato")
text(0.6, 3, "previa")
Una vez que tenemos la posterior, podemos definir distitas cosas. Por ejemplo el valor esperado de tasa de remoción por fruto es (alpha+sum(removidos)) / (alpha+sum(removidos) + beta+sum(frutos-removidos)) = 0.1936128. La probabilidad de que la tasa de remoción por fruto sea menor al 16% es
pbeta(0.16, alpha + sum(removidos), beta + sum(frutos - removidos))## [1] 0.002478181
1 - pbeta(0.2, alpha + sum(removidos), beta + sum(frutos - removidos))## [1] 0.3002499
qbeta(c(0.025, 0.975), alpha + sum(removidos), beta + sum(frutos - removidos))## [1] 0.1697490 0.2186349
Just Another Gibbs Sampler
Vamos a hacer el mismo análisis pero con JAGS. JAGS es un software que programa cadenas de Markov chain Monte Carto (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. (2012). JAGS version 3.3.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{
theta ~ dbeta(1, 1)
for (i in 1:nobs) {
removidos[i] ~ dbin(theta, frutos)
}
}"
)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 una lista de 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 # si descartamos 'thin' algunos valores
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 Size: 14
##
## 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.
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.001 minutes at time 2015-09-30 07:55:50.
##
## mean sd 2.5% 50% 97.5% overlap0 f Rhat n.eff
## theta 0.193 0.013 0.170 0.193 0.218 FALSE 1 1.002 1500
## deviance 55.942 1.446 54.925 55.429 60.396 FALSE 1 1.013 505
##
## 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 and DIC = 56.984
## DIC is an estimate of expected predictive error (lower is better).
jagsUI::traceplot(m.sim, parameters = c("theta"))
Podemos graficar las posterior aproximada y compararla con la analítica
plot(density(m.sim$sims.list$theta), lwd = 2, main = "", ylim = c(0, 35))
lines(xvec, dbeta(xvec, alpha + sum(removidos), beta + sum(frutos - removidos)),
lwd = 2, col = "gray")
Y podemos comparar las estimaciones acerca de la posterior con las que calculamos analíticamente. Por ejemplo, el valor esperado de tasa de remoción por fruto es
m.sim$mean$theta## [1] 0.1934056
# o bien
mean(m.sim$sims.list$theta)## [1] 0.1934056
# que es parecido a
(alpha + sum(removidos))/(alpha + sum(removidos) + beta + sum(frutos - removidos))## [1] 0.1936128
Si queremos ver la probabilidad de que la tasa de remoción por fruto sea menor al 15% hacemos
length(which(m.sim$sims.list$theta < 0.16))/length(m.sim$sims.list$theta)## [1] 0.002666667
length(which(m.sim$sims.list$theta > 0.2))/length(m.sim$sims.list$theta)## [1] 0.312
Con JAGS (o BUGS) podemos ajustar modelos con distinto grado de complejidad. Vamos a ir de a poco... Para empezar, veamos una regresión Binomial. Para corroborar que todo funciona, vamos a simular datos. Veamos por ejemplo un caso en el que las plantas tienen distinto tamaño de cosecha, y que los consumidores se ven atraidos a estas plantas.
nplants <- 20
lambdas <- seq(50, 150, length = nplants)
tam <- rpois(nplants, lambdas)
b0 <- -5
b1 <- 0.05
p.rem <- plogis(b0 + b1 * tam)
removidos <- rbinom(nplants, size = tam, prob = p.rem)Para ajustar el modelo tenemos que modificar el modelo bugs
cat(file = "remr.bug",
"
model{
b0 ~ dnorm(0, 0.001)
b1 ~ dnorm(0, 0.001)
for (i in 1:nplants) {
removidos[i] ~ dbin(prem[i], tam[i])
logit(prem[i]) <- b0 + b1*tam[i]
}
}"
)Como antes, tenemos que armar una lista de datos para pasarle a jags, una función para los valores iniciales y la lista de parámetros que queremos tener en cuenta.
m.data <- list("nplants", "tam", "removidos")
inits <- function() list(b0 = rnorm(1, -5, 0.1), b1 = rnorm(1, 0, 0.1))
params <- c("b0", "b1")
ni <- 1000 # número de iteraciones
nt <- 1 # si descartamos 'thin' algunos valores
nb <- 500 # cuantas iteraciones usamos de 'burn in'
nc <- 3 # y cuantas cadenas corremos
library(jagsUI)
mr.sim <- jags(data = m.data, inits = inits, parameters.to.save = params, model.file = "remr.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 Size: 105
##
## 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.
print(mr.sim)## JAGS output for model 'remr.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 2015-09-30 07:55:50.
##
## mean sd 2.5% 50% 97.5% overlap0 f Rhat n.eff
## b0 -5.148 0.254 -5.695 -5.150 -4.636 FALSE 1 1.007 368
## b1 0.052 0.002 0.047 0.051 0.056 FALSE 1 1.004 573
## deviance 119.182 1.941 117.381 118.598 124.347 FALSE 1 1.018 601
##
## 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.9 and DIC = 121.063
## DIC is an estimate of expected predictive error (lower is better).
Como siempre, tenemos que ver que todas las cadenas hayan converjido y que tengamos un tamaño de muestra de las posteriores más o menos razonables antes de hacer inferencias. Este ejemplo usa previas no informativas para los parámetros de la regresión, de manera que tenemos los mismos peligros que vimos en el ejemplo de pruebas de potencia. La diferencia es que podríamos usar previas más informativas para evitar errores tipo M y tipo S.
sessionInfo()## R version 3.2.2 (2015-08-14)
## Platform: x86_64-apple-darwin13.4.0 (64-bit)
## Running under: OS X 10.10.5 (Yosemite)
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] jagsUI_1.3.7 lattice_0.20-33 knitr_1.11
##
## loaded via a namespace (and not attached):
## [1] digest_0.6.8 mime_0.4 grid_3.2.2
## [4] knitrBootstrap_1.0.0 formatR_1.2 magrittr_1.5
## [7] coda_0.17-1 evaluate_0.7.2 stringi_0.5-5
## [10] rjags_3-15 rmarkdown_0.8 tools_3.2.2
## [13] stringr_1.0.0 markdown_0.7.7 parallel_3.2.2
## [16] yaml_2.1.13 htmltools_0.2.6
Juan Manuel Morales. 26 de Septiembre del 2015. Última actualización: 2015-09-30