Inferência bayesiana
Quando estamos falando de Inferência nosso objetivo normalmente é tentar verificar alguma informação sobre uma quantidade desconhecida.
TL;DR
- Na inferência bayesiana, o parâmetro é tratado como aleatório e sua "crença" é representada por uma distribuição de probabilidade.
- O pacote
R2jagsajusta modelos bayesianos no R chamando o JAGS (amostrador de Gibbs).- O
mcmcplotse osuperdiagdiagnosticam a convergência das cadeias MCMC.
Para isso devemos utilizar toda informação disponível, seja ela objetiva ou subjetiva (isto é, vinda de umam amostra ou de algum conhecimento préveo ou intuitivo)
Segundo o ponto de vista Bayesiano essa informação subjetiva também será incorporada na análise graças ao teorema de bayes.
Foi o Statistical Rethinking, do McElreath, que me fez realmente entender bayesiana como um jeito de raciocinar e não só como fórmula pra decorar. Recomendo bastante pra quem tá começando a mexer com isso agora.
Como no ponto de vista Bayesiano atribuímos aleatoriedade ao parâmetro, nossa “crença” será representada por uma distribuição de probabilidade (ou modelo probabilístico)
Teorema de bayes: \[ p(\theta|x)=\frac{p(x,\theta)}{p(x)}=\frac{p(x|\theta)p(\theta)}{p(x)} \]
onde:
- \(p(x|\theta)\): função de verossimilhança (modelo)
- \(p(\theta)\): distribuição a priori
- \(p(x)\): distribuição marginal de \(x\).
A estimação muitas vezes envolve o cálculo de integrais nada simples analiticamente porém, alguns algorítimos como o amostrador de Gibbs pode relizar aproximações muito relevantes.
Modelo linear bayesiano
Para entender como funciona o modelo bayesiano, primeiramente vamos começar com algo bem simples, suponha:
\[ Y_i \sim N(\mu_i,\tau) \] onde \(\mu\) é definido como \(\mu_i= X \mathbf{\beta}\).
Incialmente vamos considerar que não existe relação nenhuma, então utilizaremos a priori:
\[ \beta \sim N(0,\tau_{\beta}) \]
onde \(\tau\) é conhecido.
Nem sempre é uma tarefa simples determinar a distribuição posteri de um modelo bayesiano e é neste ponto que o pacote jagsserá bastante útil (existem outras alternativas como o WinBugs, OpenBugs, Stan, mas aqui resolvi trazer apenas o jags por possuir vantagens bem interessantes.)
Jags
O pacote R2jags é exatamente o que seu nome significa: “Just Another Gibbs Sampler”. Possui as mesmas funcionalidades do nosso querido OpenBugs possibilitando também que seja utilizado inteiramente dentro do ambiente R.
Assim como o OpenBugs, ele também trabalha chamando o software oficial que precisa ser baixado no site.
Para começar a utilizar basta baixar o pacote e acessá-lo na biblioteca:
library(R2jags)Declarando o modelo
A base de dados que será utilizada para ajustar o modelo será a base nativa do R chamada trees:
X<-trees[,1:2] #Matriz de variáveis explanatórias
Y<- trees[,3] #Vetor da variável resposta
p <- ncol(X) #p é o número de parâmetros do modelo (nesse caso é o número de colunas)
n <- nrow(X) #n é o número de observações do modeloO modelo deve estar declarado e salvo em um arquivo .txt (ou mesmo um outro arquivo .r) da seguinte maneira:
### Declarando o modelo Bayesiano
sink("linreg.txt")
cat("
model {
# Prioris
for(j in 1:p)
{
beta[j] ~ dnorm(mu.beta, tau.beta)
}
sigma ~ dunif(0, 100)
tau <- 1/ (sigma * sigma)
# Verossimilhança
for (i in 1:n) {
y[i] ~ dnorm(mu[i], tau)
mu[i] <- inprod(X[i,], beta)
}
}
",fill=TRUE)
sink()Uma vez que o modelo esta declarado, é a hora de nomear os parametros da função que fará o ajuste do modelo
#Parametros da Priori
mu.beta <- 0
tau.beta <- 0.001
#Set Working Directory
wd <- getwd()
# Junte os dados em uma lista
win.data <- list(X=X,y=Y,p=p,n=n,mu.beta=mu.beta,tau.beta=tau.beta)
# Função de inicialização
inits <- function(){ list(beta=rnorm(p), sigma = rlnorm(1))}
# Os parametros que desejamos estimar
params <- c("beta","sigma","tau")
# Caracteristicas do MCMC
n.burnin <- 500 #Número de iterações que serão descartadas
n.thin <- 10 #para economizar memória e tempo de computação se n.iter for grande
n.post <- 5000
n.chains <- 3 #Número de cadeias
n.iter <- n.burnin + n.thin*n.post #Número de iteraçõesImplementando o modelo
Após ter em mãos todos esses resultados, já podemos ajustar o modelo com o comando jags(), veja:
bayes.mod.fit <-jags(data = win.data,
inits = inits,
parameters = params,
model.file = "linreg.txt", # O arquivo "linreg.txt" deve estar no mesmo diretório
n.iter = n.iter,
n.thin=n.thin,
n.burnin=n.burnin,
n.chains=n.chains,
working.directory=wd,DIC = T)## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 31
## Unobserved stochastic nodes: 3
## Total graph size: 166
##
## Initializing modelprint(bayes.mod.fit, dig = 3)## Inference for Bugs model at "linreg.txt", fit using jags,
## 3 chains, each with 50500 iterations (first 500 discarded), n.thin = 10
## n.sims = 15000 iterations saved
## mu.vect sd.vect 2.5% 25% 50% 75% 97.5% Rhat n.eff
## beta[1] 5.045 0.435 4.183 4.757 5.043 5.324 5.916 1.001 15000
## beta[2] -0.478 0.078 -0.633 -0.527 -0.477 -0.427 -0.324 1.001 15000
## sigma 6.448 0.904 4.995 5.805 6.335 6.970 8.502 1.001 15000
## tau 0.025 0.007 0.014 0.021 0.025 0.030 0.040 1.001 15000
## deviance 201.924 2.682 198.881 199.970 201.244 203.149 208.856 1.001 7200
##
## For each parameter, n.eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor (at convergence, Rhat=1).
##
## DIC info (using the rule, pD = var(deviance)/2)
## pD = 3.6 and DIC = 205.5
## DIC is an estimate of expected predictive error (lower deviance is better).Com os resultados em mãos podemos avaliar o ajuste do modelo, o jags nos fornece os intervalos de credibilidade e o Rhat, que é a convergência da cadeia, a princípio vamos apenas considerar o fato de que quanto mais próximo de 1, melhor são as estimativas.
Não vou me extender neste post com a interpretação do modelo pois o objetivo esta sendo mostrar a funcionalidade do jags em conjunto com o R.
Diagnósticos do modelo com mcmcplots
Para o diagnóstico do modelo podemos utilizar o pacote mcmcplots que fornece de maneira bem agradável os resultados gerados pelo amostrador, primeiramente vamos carregar o pacote:
library(mcmcplots)Em seguida precisar informar para o R que o resultado do algorítimo se trata de um objeto mcmc, portanto:
bayes.mod.fit.mcmc <- as.mcmc(bayes.mod.fit)
summary(bayes.mod.fit.mcmc)##
## Iterations = 1:49991
## Thinning interval = 10
## Number of chains = 3
## Sample size per chain = 5000
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## beta[1] 5.04490 0.435344 3.555e-03 3.555e-03
## beta[2] -0.47754 0.077588 6.335e-04 6.335e-04
## deviance 201.92383 2.682384 2.190e-02 2.144e-02
## sigma 6.44763 0.903646 7.378e-03 7.359e-03
## tau 0.02542 0.006784 5.539e-05 5.524e-05
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## beta[1] 4.18250 4.75721 5.04333 5.32437 5.91642
## beta[2] -0.63255 -0.52732 -0.47726 -0.42674 -0.32376
## deviance 198.88143 199.97019 201.24393 203.14881 208.85648
## sigma 4.99470 5.80492 6.33492 6.96990 8.50193
## tau 0.01383 0.02058 0.02492 0.02968 0.04008O pacote nos fornece alguns tipos de gráficos para diagnóstico
caterplot(bayes.mod.fit.mcmc) #Observando todas as estimativas
caterplot(bayes.mod.fit.mcmc,parms = params) #Observando as estimativas de todos os parâmetros menos o desvio
denplot(bayes.mod.fit.mcmc) #Densidade das estimativas de cada cadeia
traplot(bayes.mod.fit.mcmc,greek = T) #Avaliando a convergência
E por fim, para diagnósticos rápidos, pode produzir arquivos html com traço, densidade e autocorrelação.
O comando traça tudo em uma página e os arquivos serão exibidos em seu navegador de internet padrão.
mcmcplot(bayes.mod.fit.mcmc)Vai retornar um relatório resumido para todos os parâmetros como nesta imagem da internet como:
Como o objetivo do post é trazer a funcionalidade do pacote, vou apenas deixar ilustrado quais são algumas das funções mais comumente utilizadas para avaliar estatísticamente o desempenho dos modelos.
Diagnosticos estatísticos do modelo:
#Mais diagnosticos:
gelman.plot(bayes.mod.fit.mcmc)
geweke.diag(bayes.mod.fit.mcmc)## [[1]]
##
## Fraction in 1st window = 0.1
## Fraction in 2nd window = 0.5
##
## beta[1] beta[2] deviance sigma tau
## -1.6717 1.1790 -0.4485 0.1854 -0.6815
##
##
## [[2]]
##
## Fraction in 1st window = 0.1
## Fraction in 2nd window = 0.5
##
## beta[1] beta[2] deviance sigma tau
## 0.37278 -0.36960 -0.24342 -0.08007 0.30725
##
##
## [[3]]
##
## Fraction in 1st window = 0.1
## Fraction in 2nd window = 0.5
##
## beta[1] beta[2] deviance sigma tau
## -0.15725 0.19911 -0.08445 -0.34043 0.35357geweke.plot(bayes.mod.fit.mcmc)
raftery.diag(bayes.mod.fit.mcmc)## [[1]]
##
## Quantile (q) = 0.025
## Accuracy (r) = +/- 0.005
## Probability (s) = 0.95
##
## Burn-in Total Lower bound Dependence
## (M) (N) (Nmin) factor (I)
## beta[1] 20 39950 3746 10.70
## beta[2] 20 36200 3746 9.66
## deviance 20 37410 3746 9.99
## sigma 20 38030 3746 10.20
## tau 20 36800 3746 9.82
##
##
## [[2]]
##
## Quantile (q) = 0.025
## Accuracy (r) = +/- 0.005
## Probability (s) = 0.95
##
## Burn-in Total Lower bound Dependence
## (M) (N) (Nmin) factor (I)
## beta[1] 20 38030 3746 10.20
## beta[2] 20 36800 3746 9.82
## deviance 20 37410 3746 9.99
## sigma 20 37410 3746 9.99
## tau 20 35610 3746 9.51
##
##
## [[3]]
##
## Quantile (q) = 0.025
## Accuracy (r) = +/- 0.005
## Probability (s) = 0.95
##
## Burn-in Total Lower bound Dependence
## (M) (N) (Nmin) factor (I)
## beta[1] 20 37410 3746 9.99
## beta[2] 20 38030 3746 10.20
## deviance 20 37410 3746 9.99
## sigma 30 40620 3746 10.80
## tau 20 39300 3746 10.50heidel.diag(bayes.mod.fit.mcmc)## [[1]]
##
## Stationarity start p-value
## test iteration
## beta[1] passed 1 0.292
## beta[2] passed 1 0.455
## deviance passed 1 0.733
## sigma passed 1 0.881
## tau passed 1 0.816
##
## Halfwidth Mean Halfwidth
## test
## beta[1] passed 5.0481 0.012089
## beta[2] passed -0.4780 0.002155
## deviance passed 201.8829 0.073069
## sigma passed 6.4367 0.024544
## tau passed 0.0255 0.000187
##
## [[2]]
##
## Stationarity start p-value
## test iteration
## beta[1] passed 1 0.246
## beta[2] passed 1 0.249
## deviance passed 1 0.967
## sigma passed 1 0.950
## tau passed 1 0.770
##
## Halfwidth Mean Halfwidth
## test
## beta[1] passed 5.0386 0.011955
## beta[2] passed -0.4765 0.002134
## deviance passed 201.9023 0.068414
## sigma passed 6.4571 0.025014
## tau passed 0.0253 0.000188
##
## [[3]]
##
## Stationarity start p-value
## test iteration
## beta[1] passed 1 0.657
## beta[2] passed 1 0.690
## deviance passed 1 0.544
## sigma passed 1 0.813
## tau passed 1 0.873
##
## Halfwidth Mean Halfwidth
## test
## beta[1] passed 5.0480 0.012156
## beta[2] passed -0.4781 0.002163
## deviance passed 201.9863 0.076685
## sigma passed 6.4491 0.025385
## tau passed 0.0254 0.000188Diagnostico de convergencia rapida: superdiag
Uma função muito conveniente para analisar representações numéricas de diagnósticos em um ajuste é o pacote superdiag de Tsai, Gill e Rapkin, 2012 que trás uma série de estatísticas para avaliar o desempenho dos ajustes do modelo.
library(superdiag)
superdiag(bayes.mod.fit.mcmc, burnin = 100)## Number of chains = 3
## Number of iterations = 5000 per chain before discarding the burn-in period
## Burn-in period = 100 per chain
## Sample size in total = 14703
##
## ****************** The Geweke diagnostic: ******************
## Windows:
## chain 1 chain 2 chain 3
## From start 0.1 0.5420 0.2999
## From stop 0.5 0.3511 0.6893
##
## Z-scores:
## chain 1 chain 2 chain 3
## beta[1] -1.85586 0.3331 -1.66699
## beta[2] 1.57605 -0.2271 1.53584
## deviance 0.02463 0.3356 -1.14324
## sigma -0.15363 -0.8820 -0.33962
## tau -0.09745 0.9937 0.01232
##
## *************** The Gelman-Rubin diagnostic: ***************
## Potential scale reduction factors:
## Point est. Upper C.I.
## beta[1] 1.0001 1.001
## beta[2] 1.0000 1.000
## deviance 1.0009 1.002
## sigma 1.0002 1.001
## tau 0.9999 1.000
##
## Multivariate psrf: 1.0005
##
## ************* The Heidelberger-Welch diagnostic ************
## Chain 1:
## epsilon=0.1, alpha=0.05
## Stationarity start p-value
## test iteration
## beta[1] passed 1 0.1576
## beta[2] passed 1 0.2864
## deviance passed 1 0.8399
## sigma passed 1 0.8207
## tau passed 1 0.7405
##
## Halfwidth Mean Halfwidth
## test
## beta[1] passed 5.04671 0.012211
## beta[2] passed -0.47775 0.002177
## deviance passed 201.89097 0.074094
## sigma passed 6.43566 0.024772
## tau passed 0.02549 0.000189
##
## Chain 2:
## epsilon=0.079, alpha=0.1
## Stationarity start p-value
## test iteration
## beta[1] passed 1 0.3032
## beta[2] passed 1 0.3259
## deviance passed 1 0.9562
## sigma passed 1 0.7462
## tau passed 1 0.5362
##
## Halfwidth Mean Halfwidth
## test
## beta[1] passed 5.03850 0.0120853
## beta[2] passed -0.47646 0.0021574
## deviance passed 201.90084 0.0693125
## sigma passed 6.45467 0.0252168
## tau passed 0.02536 0.0001894
##
## Chain 3:
## epsilon=0.054, alpha=0.005
## Stationarity start p-value
## test iteration
## beta[1] passed 1 0.5489
## beta[2] passed 1 0.5665
## deviance passed 1 0.5038
## sigma passed 1 0.8038
## tau passed 1 0.8898
##
## Halfwidth Mean Halfwidth
## test
## beta[1] passed 5.04719 0.0122925
## beta[2] passed -0.47794 0.0021858
## deviance passed 201.98956 0.0775537
## sigma passed 6.44893 0.0256817
## tau passed 0.02544 0.0001937
##
## *************** The Raftery-Lewis diagnostic ***************
## Chain 1:
## Convergence eps = 0.001
## Quantile (q) = 0.025
## Accuracy (r) = +/- 0.005
## Probability (s) = 0.95
##
## Burn-in Total Lower bound Dependence
## (M) (N) (Nmin) factor (I)
## beta[1] 30 40170 3746 10.70
## beta[2] 20 36340 3746 9.70
## deviance 20 38200 3746 10.20
## sigma 20 38200 3746 10.20
## tau 20 36950 3746 9.86
##
## Chain 2:
## Convergence eps = 5e-04
## Quantile (q) = 0.25
## Accuracy (r) = +/- 0.001
## Probability (s) = 0.99
##
## You need a sample size of at least 1244044 with these values of q, r and s
##
## Chain 3:
## Convergence eps = 0.005
## Quantile (q) = 0.25
## Accuracy (r) = +/- 5e-04
## Probability (s) = 0.999
##
## You need a sample size of at least 8120675 with these values of q, r and s
##
## ************* The Hellinger distance diagnostic ************
## Between chains:
## Min Max
## beta[1] 0.01735 0.02915
## beta[2] 0.02015 0.02620
## deviance 0.03155 0.03413
## sigma 0.01858 0.02731
## tau 0.01538 0.02810
##
## Within chain 1:
## 980 1960 2940 3920
## beta[1] 0.05231 0.03952 0.04017 0.04259
## beta[2] 0.04261 0.05034 0.04320 0.04782
## deviance 0.05880 0.04060 0.06297 0.04311
## sigma 0.03871 0.03667 0.06465 0.04285
## tau 0.03668 0.03996 0.03633 0.04083
##
## Within chain 2:
## 980 1960 2940 3920
## beta[1] 0.03098 0.04075 0.04281 0.03887
## beta[2] 0.03050 0.03770 0.03887 0.04216
## deviance 0.04541 0.03992 0.03390 0.04730
## sigma 0.04660 0.03876 0.03090 0.02866
## tau 0.03648 0.03773 0.02967 0.03589
##
## Within chain 3:
## 980 1960 2940 3920
## beta[1] 0.03356 0.03988 0.03146 0.02986
## beta[2] 0.03425 0.04729 0.03175 0.03219
## deviance 0.05894 0.03553 0.05018 0.04509
## sigma 0.04392 0.04245 0.03858 0.03760
## tau 0.04089 0.03458 0.04512 0.03047Para finalizar, outra função que pode ser útil pata atualizando o modelo, se necessário - por exemplo, se não houver convergência ou pouca convergencia:
bayes.mod.fit.upd <- update(bayes.mod.fit, n.iter=1000)
bayes.mod.fit.upd <- autojags(bayes.mod.fit)Muito a estudar
Assim como toda a Estatística, inferência bayesiana não funciona se a teoria não for aplicada corretamente. É uma ferramenta muito poderosa e necessita ser usada com cautela pois demanda bastante o uso de metodologias estatísticas.
Como dizia o tio Ben: “grandes poderes trazem grandes responsabilidades” então vamos tomar cuidado com os resultados que encontramos. Para outras formas de avaliar e apresentar o ajuste de modelos, veja os posts sobre análise multivariada e sobre tabelas para apresentar resultados. Para ver os cálculos por trás de um modelo bayesiano implementados do zero (sem o JAGS), veja este outro post.
Atualização (2026)
Nota acrescentada em 2026. O JAGS continua funcionando e é uma ótima porta de entrada para MCMC, mas quem está começando um projeto novo hoje frequentemente prefere o Stan (via rstan ou o pacote brms, que gera o código Stan automaticamente a partir de uma sintaxe parecida com lm()). A lógica de priori, verossimilhança e diagnóstico de convergência deste post vale para qualquer um dos dois.
Perguntas frequentes
Qual a diferença entre estatística bayesiana e frequentista?
Na abordagem frequentista, o parâmetro é fixo e desconhecido, e a incerteza está na amostra. Na bayesiana, o parâmetro é tratado como uma variável aleatória com uma distribuição de probabilidade (a priori), atualizada pelos dados observados via o teorema de Bayes para gerar a distribuição posterior.
Como sei se minha cadeia MCMC convergiu?
O indicador mais usado é o Rhat (fator de redução de escala potencial): valores próximos de 1 indicam convergência. Complementarmente, o pacote mcmcplots permite inspecionar visualmente o traço das cadeias, e testes como Geweke e Heidelberger-Welch (do pacote superdiag) formalizam esse diagnóstico.
Preciso instalar algo além do R para usar JAGS?
Sim. O JAGS é um software externo (não é só um pacote R) e precisa ser baixado e instalado separadamente antes de carregar o R2jags no R, que funciona como uma interface para chamá-lo.
Referencias
- Uma primeira olhada em estatística bayesiana e linguagem BUGS por Augusto Ribas - blog Recologia
- KRUSCHKE, J. K. Doing Bayesian Data Analysis: A Tutorial with R, JAGS, and Stan. 2. ed. Academic Press / Elsevier, 2014.
