Análisis Bayesianos

Los análisis Bayesianos son similares a los que vimos en máxima verosimilitud en el sentido de 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 datos nuevos que todavía no hemos colectado. 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)}}{ p(\boldsymbol{y}) } = \frac{p(\boldsymbol{y} \lvert \boldsymbol{\theta}) p(\boldsymbol{\theta)}}{\int p(\boldsymbol{y} \lvert \boldsymbol{\theta)} p(\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 que define la previa. 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 generar datos con la computadora a partir de un modelo (vamos a simular los datos). Imaginen que queremos estudiar la remoción de frutos en \(30\) plantas. En cada una de las plantas marcamos \(20\) frutos y contamos cuántos son removidos por dispersores luego de un tiempo fijo. Si suponemos que un buen modelo para este tipo de datos es una distribución Binomial con una probabilidad de éxito fija hacemos:

set.seed(1234)
nobs <- 30    # número de observaciones (plantas)
frutos <- rep(20, nobs) # frutos disponibles
p_rem <- 0.2  # probabilidad de remoción por fruto
removidos <- rbinom(nobs, size = frutos, prob = p_rem)

El modelo de datos (cuántos frutos son removidos) en este casi es una Binomial con número de pruebas (la cantidad de frutos disponibles) conocido. Para hacer un análisis Bayesiano de estos datos, tenemos que definir una previa para la probabilidad de éxito (p_rem). Esa previa tiene que tomar valores continuos entre \(0\) y \(1\). Una opción sería una distribución uniforme con esos límites, pero si usamos una distribución Beta, 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

alpha_p <- alpha + sum(removidos)
beta_p  <- beta  + sum(frutos - removidos)

Para obtener el valor esperado de una distribución Beta hacermos

alpha_p / (alpha_p + beta_p)
## [1] 0.1877076

Eso nos da un estimador puntual de la probabilidad de remoción por fruto p_rem. Para tener una medida de incertidumbre alrededor de este valor, podemos ver los cuantiles de la posterior

qbeta(c(0.025, 0.975), alpha_p, beta_p)
## [1] 0.1575462 0.2198340

También podemos graficar la distribución posterior y compararla con la previa para ver cuánto aprendimos haciendo el análisis de datos.

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, 12, "posterior")

par(op)

¿Cómo sabemos si estas estimaciones tienen sentido para nuestros datos? En este caso, la pregunta es trivial porque conocemos cómo se generaron los datos, pero cuando trabajamos con datos de verdad, el modelo de datos es un supuesto y tenemos que ver si ese supuesto tiene sentido.

Una opción para contestar esta pregunta es hacer simulaciones a partir de la posterior y compararlas con los datos.

nreps <- 10000
vals  <- 0:20 # posibles valores de remoción
res <- matrix(NA, nreps, length(vals) - 1) # matriz para resultados
p_sim <- rbeta(nreps, alpha_p, beta_p) # muestra aleatoria de la posterior

for(i in 1:nreps){
  tmp <- rbinom(nobs, frutos, p_sim[i])
  res[i, ] <- hist(tmp, right = FALSE, breaks = vals, plot = FALSE)$density
}

plot(table(removidos)/nobs, xlim = c(0, 10), ylim = c(0,0.6),
     ylab = "frecuencia", type = "p", pch = 19)
library(coda)
ci <- HPDinterval(as.mcmc(res))
lines(0:19, ci[,2])
lines(0:19, ci[,1])

En el gráfico podemos ver que las observaciones (puntos), están dentro del intervalo de credibilidad para muestras de una Binomial con probabilidad de éxito dada por la posterior de p_rem.

Veamos ahora un análisis con datos del mundo real. Vamos a usar los datos de remoción de frutos de quintral por el monito del monte que vimos en el ejercicio anterior.

url <- "https://github.com/jmmorales/cursoMyD/raw/master/Data/quintral.txt"
quintral <- read.table(url, header = TRUE)

# previas
alpha <- 1
beta  <- 1

alpha_p <- alpha + sum(quintral$Removidos)
beta_p  <- beta  + sum(quintral$Frutos - quintral$Removidos)

op <- par(cex.lab = 1.5 , font.lab = 2, cex.axis = 1.3, las = 1, bty = "n")
curve(dbeta(x, alpha_p, beta_p), lwd = 2, ylab = "Densidad de probabilidad", xlab = "Probabilidad de remoción")
curve(dbeta(x, 1, 1), lwd = 2, col = "gray", add = TRUE)
text(0.2, 2.5, "previa")
text(0.65, 20, "posterior")

par(op)

Vemos que la posterior es bastante acotada, implicando que tendríamos bastante seguridad respecto a que la probabilidad de remoción de un fruto de quintral es cercana al \(60 \%\).

Veamos qué pasa si simulamos datos con este modelo y esta posterior y contrastamos con las observaciones.

nreps <- 10000
nobs <- length(quintral$Removidos)

vals  <- 0:60 # posibles valores de remoción
res <- matrix(NA, nreps, length(vals) - 1) # matriz para resultados
p_sim <- rbeta(nreps, alpha_p, beta_p) # muestra aleatoria de la posterior

for(i in 1:nreps){
  tmp <- rbinom(nobs, quintral$Frutos, p_sim[i])
  res[i, ] <- hist(tmp, right = FALSE, breaks = vals, plot = FALSE)$density
}

plot(table(quintral$Removidos)/nobs, xlim = c(0, 60), ylim = c(0,0.2),
     ylab = "frecuencia", xlab = "removidos", type = "p", pch = 19)
library(coda)
ci <- HPDinterval(as.mcmc(res))
lines(0:59, ci[,2])
lines(0:59, ci[,1])

Pregunta: ¿Qué pueden decir sobre la probabilidad de observar estos datos bajo el modelo Binomial que ajustamos?

En general, trabajamos con modelos de datos que no tienen solución analítica para las posteriores y usamos algúm método de tipo Markov Chain Monte Carlo para obtener muestras de las posteriores. Veamos como hacer eso para este modelo con JAGS. Primero, repitamos el proceso de generar datos y calcular analíticamente las posteriores para poder comparar con los resultados de JAGS

set.seed(123)
nobs <- 30    # número de observaciones (plantas)
frutos <- rep(20, nobs) # frutos disponibles
p_rem <- 0.2  # probabilidad de remoción por fruto
removidos <- rbinom(nobs, size = frutos, prob = p_rem)
# previas
alpha <- 1
beta  <- 1

alpha_p <- alpha + sum(removidos)
beta_p  <- beta  + sum(frutos - removidos)

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 = nobs, 
               frutos = frutos,
               removidos = removidos) 
inits  <- function() list(theta = runif(1, 0, 1))
params <- c("theta")

Finalmente, definimos cuántas iteraciones queremos ejecutar por cada 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 corremos

Para 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: 30
##    Unobserved stochastic nodes: 1
##    Total graph size: 64
## 
## Initializing model
## 
## Adaptive phase..... 
## Adaptive phase complete 
##  
## 
##  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). 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 de sacar cualquier conclusión a partir de 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,
## adaptation = 100 iterations (sufficient),
## burn-in = 500 iterations and thin rate = 1,
## yielding 1500 total samples from the joint posterior. 
## MCMC ran for 0 minutes at time 2021-12-06 20:11:17.
## 
##             mean    sd    2.5%     50%   97.5% overlap0 f  Rhat n.eff
## theta      0.222 0.017   0.190   0.222   0.257    FALSE 1 1.002  1432
## deviance 121.724 1.427 120.733 121.167 125.925    FALSE 1 1.006   712
## 
## 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 = 122.74 
## DIC is an estimate of expected predictive error (lower is better).

Pregunta: ¿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.214 0.21 0.205 0.217 0.214 ...
##   ..$ deviance: num [1:1500] 121 121 122 121 121 ...
##  $ mean       :List of 2
##   ..$ theta   : num 0.222
##   ..$ deviance: num 122
##  $ sd         :List of 2
##   ..$ theta   : num 0.0169
##   ..$ deviance: num 1.43
##  $ q2.5       :List of 2
##   ..$ theta   : num 0.19
##   ..$ deviance: num 121
##  $ q25        :List of 2
##   ..$ theta   : num 0.211
##   ..$ deviance: num 121
##  $ q50        :List of 2
##   ..$ theta   : num 0.222
##   ..$ deviance: num 121
##  $ q75        :List of 2
##   ..$ theta   : num 0.233
##   ..$ deviance: num 122
##  $ q97.5      :List of 2
##   ..$ theta   : num 0.257
##   ..$ deviance: num 126
##  $ 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 1432
##   ..$ deviance: num 712
##  $ pD         : num 1.02
##  $ DIC        : num 123
##  $ summary    : num [1:2, 1:11] 0.2224 121.7236 0.0169 1.4265 0.1904 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:2] "theta" "deviance"
##   .. ..$ : chr [1:11] "mean" "sd" "2.5%" "25%" ...
##  $ samples    :List of 3
##   ..$ : 'mcmc' num [1:500, 1:2] 0.214 0.21 0.205 0.217 0.214 ...
##   .. ..- attr(*, "dimnames")=List of 2
##   .. .. ..$ : NULL
##   .. .. ..$ : chr [1:2] "theta" "deviance"
##   .. ..- attr(*, "mcpar")= num [1:3] 501 1000 1
##   ..$ : 'mcmc' num [1:500, 1:2] 0.2 0.25 0.191 0.223 0.203 ...
##   .. ..- attr(*, "dimnames")=List of 2
##   .. .. ..$ : NULL
##   .. .. ..$ : chr [1:2] "theta" "deviance"
##   .. ..- attr(*, "mcpar")= num [1:3] 501 1000 1
##   ..$ : 'mcmc' num [1:500, 1:2] 0.236 0.223 0.237 0.186 0.229 ...
##   .. ..- 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 9
##   ..$ n.chains        : num 3
##   ..$ n.adapt         : num 100
##   ..$ sufficient.adapt: logi TRUE
##   ..$ 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] 121.667 0.206
##   .. .. ..- attr(*, "names")= chr [1:2] "deviance" "theta"
##   .. ..$ : Named num [1:2] 121.861 0.204
##   .. .. ..- attr(*, "names")= chr [1:2] "deviance" "theta"
##   .. ..$ : Named num [1:2] 120.784 0.218
##   .. .. ..- attr(*, "names")= chr [1:2] "deviance" "theta"
##   ..$ elapsed.mins    : num 0
##  $ run.date   : POSIXct[1:1], format: "2021-12-06 20:11:17"
##  $ parallel   : logi FALSE
##  $ bugs.format: logi FALSE
##  $ calc.DIC   : logi TRUE
##  - 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"))

Pregunta: ¿Qué les parece importante ver en este tipo de gráficos?

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, alpha_p, beta_p), add = TRUE, lwd = 2, lty = 2)

par(op)

Pregunta: ¿Qué pueden decir de las diferencias entre estas dos curvas?

Si queremos por ejemplo 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

y 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.91

Pregunta: ¿Cómo harían esos cáculos en base a los resultados teóricos?

También podemos calcular el intervalo de credibilidad:

library(coda)
HPDinterval(as.mcmc(m.sim$sims.list$theta))
##          lower     upper
## var1 0.1878296 0.2542638
## attr(,"Probability")
## [1] 0.95

A lo largo de otros prácticos iremos viendo cómo aprovechar la capacidad de generar muestras de las posteriores de nuestros modelos.

Opcional: Ajustando modelos con Stan y brms

Stan

Veamos como hacer este mismo análisis en Stan.

El proceso general es parecido a lo que hacemos con JAGS, pero tenemos que definir más detalles en el documento del modelo. Para complicar un poco la cosa, los nombres que se usan para las distribuciones son diferentes a los de JAGS o R . Por suerte, existe un manual muy detallado sobre Stan que lo pueden encontrar aquí.

En el bloque data definimos las cantidades conocidas y qué formato tienen. Luego, definimos en el bloque de parameters a las cantidades no observadas que queremos cuantificar. Finalmente. el bloque model contiene las previas y el modelo de datos.

cat(file = "rem.stan", 
"
data { 
  int<lower=1> N;             // total number of observations 
  int<lower=0> removidos[N];  // response variable 
  int<lower=1> frutos[N];
}

parameters {
  real<lower=0,upper=1> theta;
}

model{
  theta ~ beta(1,1);
  removidos ~ binomial(frutos, theta);
}
"
)

Para llamar a Stan desde R, usamos el paquete rstan. Para poder correr cadenas en paralelo usamos las opciones rstan_options(auto_write = TRUE) y options(mc.cores = parallel::detectCores()). Luego definimos los datos que vamos a pasarle a Stan y ejecutamos stan definiendo dónde está el modelo, los datos, cuántas iteraciones, el “thining” y cuántas cadenas queremos correr.

library(rstan)
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())

rem_dat <- list(N = nobs,
               removidos = removidos,
               frutos = frutos)
fit <- stan(file = 'rem.stan', data = rem_dat,
            iter = 1000, thin = 1, chains = 3)

Para ver los resultados, podemos hacer un print del objeto de salida de stan

print(fit)
## Inference for Stan model: rem.
## 3 chains, each with iter=1000; warmup=500; thin=1; 
## post-warmup draws per chain=500, total post-warmup draws=1500.
## 
##          mean se_mean   sd    2.5%     25%     50%     75%   97.5% n_eff Rhat
## theta    0.22    0.00 0.02    0.19    0.21    0.22    0.23    0.26   541 1.01
## lp__  -319.66    0.03 0.67 -321.66 -319.84 -319.39 -319.21 -319.16   612 1.01
## 
## Samples were drawn using NUTS(diag_e) at Mon Dec  6 20:11:24 2021.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at 
## convergence, Rhat=1).

Comparemos con los resultados teóricos

psam <- extract(fit, pars = c("theta"))
op <- par(cex.lab = 1.5 , font.lab = 2, cex.axis = 1.3, las = 1, bty = "n")
plot(density(psam$theta ), lwd = 2, main = "", ylim = c(0, 35), xlab = expression(theta))
curve(dbeta(x, alpha_p, beta_p), add = TRUE, lwd = 2, lty = 2)

par(op)

También podemos usar el paquete brms que sirve como un atajo para correr modelos lineales en Stan. Pero brms es para modelos de tipo regresión, así que tenemos que plantear el modelo como una regresión simple donde queremos estimar el intercepto

library(brms)

datos <- data.frame(frutos, removidos)

m1 <- brm(removidos | trials(frutos) ~ 1,
  data = datos, 
  family = binomial(),
  iter = 1000,
  chains = 3)

Podemos hacer un print para ver los resultados

print(m1)
##  Family: binomial 
##   Links: mu = logit 
## Formula: removidos | trials(frutos) ~ 1 
##    Data: datos (Number of observations: 30) 
##   Draws: 3 chains, each with iter = 1000; warmup = 500; thin = 1;
##          total post-warmup draws = 1500
## 
## Population-Level Effects: 
##           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept    -1.25      0.09    -1.44    -1.07 1.00      701      841
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

Valores negativos?! Sí, porque brms usó un “logit link” para ajustar la regresión. Si queremos valores en la escala original, tenemos que transformarlos (usando plogis). Primero recuperemos los valores de las muestras de la posterior, y luego hagamos la comparación entre la posterior y los resultados teóricos

psims <- posterior_samples(m1, "^b")
op <- par(cex.lab = 1.5 , font.lab = 2, cex.axis = 1.3, las = 1, bty = "n")
plot(density(plogis(psims$b_Intercept)), lwd = 2, main = "", ylim = c(0, 35), xlab = expression(theta))
curve(dbeta(x, alpha_p, beta_p), add = TRUE, lwd = 2, lty = 2)

par(op)

Juan Manuel Morales. 26 de Septiembre del 2015. última actualización: 2021-12-06