Fellipe Gomes

12 minute read

Inferência bayesiana

Imagem da Internet

Quando estamos falando de Inferência nosso objetivo normalmente é tentar verificar alguma informação sobre uma quantidade desconhecida.

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.

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:

suppressMessages(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 modelo

O 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ções

Implementando 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)
## module glm loaded
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 31
##    Unobserved stochastic nodes: 3
##    Total graph size: 166
## 
## Initializing model
print(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
## beta[1]    5.042   0.434   4.200   4.758   5.041   5.330   5.901 1.001
## beta[2]   -0.477   0.077  -0.629  -0.528  -0.477  -0.426  -0.326 1.001
## sigma      6.436   0.887   4.992   5.808   6.343   6.952   8.441 1.001
## tau        0.025   0.007   0.014   0.021   0.025   0.030   0.040 1.001
## deviance 201.885   2.653 198.890 199.946 201.204 203.075 208.666 1.001
##          n.eff
## beta[1]  15000
## beta[2]  15000
## sigma    15000
## tau      15000
## deviance 15000
## 
## 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.5 and DIC = 205.4
## 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:

suppressMessages(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.04204 0.433538 3.540e-03      3.551e-03
## beta[2]   -0.47726 0.077375 6.318e-04      6.224e-04
## deviance 201.88535 2.652675 2.166e-02      2.159e-02
## sigma      6.43598 0.887477 7.246e-03      7.247e-03
## tau        0.02548 0.006758 5.518e-05      5.518e-05
## 
## 2. Quantiles for each variable:
## 
##               2.5%       25%       50%       75%     97.5%
## beta[1]    4.20009   4.75805   5.04070   5.32982   5.90077
## beta[2]   -0.62947  -0.52847  -0.47701  -0.42621  -0.32594
## deviance 198.89013 199.94648 201.20413 203.07457 208.66556
## sigma      4.99192   5.80761   6.34326   6.95162   8.44128
## tau        0.01403   0.02069   0.02485   0.02965   0.04013

O 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.1443  -0.9566   0.4801   0.4322  -0.8687 
## 
## 
## [[2]]
## 
## Fraction in 1st window = 0.1
## Fraction in 2nd window = 0.5 
## 
##  beta[1]  beta[2] deviance    sigma      tau 
##   1.1847  -0.9427   1.5877   1.0883  -0.7937 
## 
## 
## [[3]]
## 
## Fraction in 1st window = 0.1
## Fraction in 2nd window = 0.5 
## 
##  beta[1]  beta[2] deviance    sigma      tau 
##   0.1536  -0.3203   0.5631  -0.1977   0.2429
geweke.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       36800 3746          9.82     
##  beta[2]  20       38030 3746         10.20     
##  deviance 20       37410 3746          9.99     
##  sigma    20       38030 3746         10.20     
##  tau      20       38660 3746         10.30     
## 
## 
## [[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       38660 3746         10.30     
##  beta[2]  20       37410 3746          9.99     
##  deviance 20       37410 3746          9.99     
##  sigma    20       36800 3746          9.82     
##  tau      20       38030 3746         10.20     
## 
## 
## [[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       39300 3746         10.50     
##  deviance 20       37410 3746          9.99     
##  sigma    20       35610 3746          9.51     
##  tau      30       40620 3746         10.80
heidel.diag(bayes.mod.fit.mcmc)
## [[1]]
##                                        
##          Stationarity start     p-value
##          test         iteration        
## beta[1]  passed       1         0.741  
## beta[2]  passed       1         0.754  
## deviance passed       1         0.223  
## sigma    passed       1         0.345  
## tau      passed       1         0.275  
##                                      
##          Halfwidth Mean     Halfwidth
##          test                        
## beta[1]  passed      5.0424 0.011659 
## beta[2]  passed     -0.4775 0.002054 
## deviance passed    201.8770 0.071447 
## sigma    passed      6.4365 0.024810 
## tau      passed      0.0255 0.000188 
## 
## [[2]]
##                                        
##          Stationarity start     p-value
##          test         iteration        
## beta[1]  passed         1       0.489  
## beta[2]  passed         1       0.465  
## deviance passed       501       0.099  
## sigma    passed         1       0.498  
## tau      passed         1       0.814  
##                                      
##          Halfwidth Mean     Halfwidth
##          test                        
## beta[1]  passed      5.0389 0.012539 
## beta[2]  passed     -0.4764 0.002150 
## deviance passed    201.8989 0.076619 
## sigma    passed      6.4324 0.024676 
## tau      passed      0.0255 0.000189 
## 
## [[3]]
##                                        
##          Stationarity start     p-value
##          test         iteration        
## beta[1]  passed       1         0.830  
## beta[2]  passed       1         0.773  
## deviance passed       1         0.212  
## sigma    passed       1         0.215  
## tau      passed       1         0.333  
##                                      
##          Halfwidth Mean     Halfwidth
##          test                        
## beta[1]  passed      5.0448 0.011953 
## beta[2]  passed     -0.4779 0.002133 
## deviance passed    201.8638 0.073117 
## sigma    passed      6.4390 0.024315 
## tau      passed      0.0254 0.000186

Diagnostico 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)
## Carregando pacotes exigidos: boa
superdiag(bayes.mod.fit.mcmc, burnin = 100)
## Number of chains = 3 
## Number of iterations = 5000 per chain before discarding the burn-in period
## The burn-in period = 100 per chain
## Sample size in total = 14700 
## 
## ********** The Geweke diagnostic: **********
## Z-scores:
##                       chain1     chain 2    chain 3
## beta[1]            0.8026201 -1.53071053  1.0931724
## beta[2]           -0.6583132  1.56659356 -1.1504195
## deviance           0.7675920  1.76928077 -2.4736755
## sigma              1.2001822  0.34672962 -0.5374265
## tau               -1.5392509 -0.07493664  0.2399271
## Window From Start  0.1000000  0.51269000  0.2127100
## Window From Stop   0.5000000  0.19937000  0.0123800
## 
## ********** The Gelman-Rubin diagnostic: **********
## Potential scale reduction factors:
## 
##          Point est. Upper C.I.
## beta[1]           1          1
## beta[2]           1          1
## deviance          1          1
## sigma             1          1
## tau               1          1
## 
## Multivariate psrf
## 
## 1
## 
## ********** The Heidelberger-Welch diagnostic: **********
## 
## Chain 1, epsilon=0.1, alpha=0.05                                       
##          Stationarity start     p-value
##          test         iteration        
## beta[1]  passed       1         0.792  
## beta[2]  passed       1         0.752  
## deviance passed       1         0.167  
## sigma    passed       1         0.208  
## tau      passed       1         0.168  
##                                      
##          Halfwidth Mean     Halfwidth
##          test                        
## beta[1]  passed      5.0420 0.01176  
## beta[2]  passed     -0.4775 0.00207  
## deviance passed    201.8802 0.07231  
## sigma    passed      6.4387 0.02511  
## tau      passed      0.0255 0.00019  
## 
## Chain 2, epsilon=0.037, alpha=0.05                                       
##          Stationarity start     p-value
##          test         iteration        
## beta[1]  passed         1       0.375  
## beta[2]  passed         1       0.389  
## deviance passed       491       0.113  
## sigma    passed         1       0.425  
## tau      passed         1       0.782  
##                                      
##          Halfwidth Mean     Halfwidth
##          test                        
## beta[1]  passed      5.0375 0.012597 
## beta[2]  passed     -0.4762 0.002168 
## deviance passed    201.8963 0.077470 
## sigma    passed      6.4330 0.024955 
## tau      passed      0.0255 0.000191 
## 
## Chain 3, epsilon=0.198, alpha=0.1                                       
##          Stationarity start     p-value
##          test         iteration        
## beta[1]  passed       1         0.817  
## beta[2]  passed       1         0.759  
## deviance passed       1         0.207  
## sigma    passed       1         0.197  
## tau      passed       1         0.318  
##                                      
##          Halfwidth Mean     Halfwidth
##          test                        
## beta[1]  passed      5.0448 0.012066 
## beta[2]  passed     -0.4779 0.002153 
## deviance passed    201.8627 0.073956 
## sigma    passed      6.4390 0.024598 
## tau      passed      0.0254 0.000188 
## 
## ********** The Raftery-Lewis diagnostic: **********
## 
## Chain 1, converge.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]  2        3635  3746         0.970     
##  beta[2]  2        3821  3746         1.020     
##  deviance 2        3696  3746         0.987     
##  sigma    2        3821  3746         1.020     
##  tau      2        3856  3746         1.030     
## 
## 
## Chain 2, converge.eps = 0.005
## Quantile (q) = 0.001
## Accuracy (r) = +/- 0.005
## Probability (s) = 0.95 
##                                                 
##           Burn-in  Total Lower bound  Dependence
##           (M)      (N)   (Nmin)       factor (I)
##  beta[1]  1        158   154          1.03      
##  beta[2]  1        158   154          1.03      
##  deviance 1        158   154          1.03      
##  sigma    1        158   154          1.03      
##  tau      1        158   154          1.03      
## 
## 
## Chain 3, converge.eps = 5e-04
## Quantile (q) = 0.05
## Accuracy (r) = +/- 0.005
## Probability (s) = 0.95 
## 
## You need a sample size of at least 7299 with these values of q, r and s

Para 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.

comments powered by Disqus