Uma primeira olhada em estatística bayesiana e linguagem BUGS

Bem, já vimos que modelos nulos e valores p não são a única forma de os dados fazerem sentido e conseguirmos probabilidades para entender o mundo que nos cerca.
A alguns posts atrás falamos de como teoria da informação pode ajudar também a entendermos o que esta acontecendo. Mas uma terceira abordagem que tem crescido, principalmente agora com a facilidade dos computadores é a bayesiana.
Eu achei que esse artigo da uma boa comparação, mas de forma reduzida, na estatística frequentista, nos estimamos parâmetros usando máxima verosimilhança e comparamos com um modelo nulo, usando a teoria da informação em derivações como o aic, nos podemos comparar modelos. Mas ainda temos a estatística bayesiana, que usa o teorema de bayes para a partir dos dados e do que acreditamos sobre os parâmetros inferir o que deve ser os reais valores dos destes.

stat_schools

Então com estatística bayesiana, tudo parace mais com como a ciência acontece, a gente acha alguma coisa, coleta dados, faz update do que a gente acha usando esses dados coletados e ai interpreta o resultado. No entanto, a parte de fazer o update do que a gente acha mais os dados para produzir a distribuição posterior, pelo menos numericamente eu acho muito complicado.

Uma das vantagens é que existem hoje alguns programas que fazem esse trabalho difícil. Principalmente na genética e evolução, me parece que a estatística bayesiana esta bem mais difundida e tem muitos programa usando ela, apesar de certas criticas relevantes, como essa do M. Henry H. Stevens.

Bem, uma linguagem mais geral para usar a estatística bayesiana é a BUGS.
Bugs é um linguagem igual o R que gera um amostrador para você (de uma olhada no teorema de bayes para saber o que diabos o amostrador amostra), então a ideia é que você decreva o modelo, fale o que você acha (sua distribuição dos parâmetros a priori) e ele gerara a distribuição a posteriori.

O Bugs foi feito em uma versão para windows, o WinBugs, que não esta mais em desenvolvimento mas tem uma versão estavel para download, e o OpenBugs, que é uma versão do WinBugs multiplataforma para Linux, Mac e Windows. Todo mundo pode usar.
Uma derivação que usa a mesma sintaxe mas é desenvolvido a parte por outras pessoas é o JAGS, que funciona igualzinho o OpenBugs que usaremos aqui.

A ultima novidade é o Stan, que usa a mesma sintaxe aproximadamente com outro amostrador, o No-U-Turn sampler, que não tenho ideia do que faz, mas o nome é bonito e segundo os autores é mais rápido.

Porém Bugs e seus derivados, no próprio site são descritos pelos autores como uma caixa de pandora, complicados de entender seu funcionamento, mas podemos resolver um exemplo com eles e ver se o resultado faz sentido. Alias, para todas essas linguagens, podemos conversar com elas direto do R. Com o openbugs temos o pacote R2OpenBugs, que usamos ele do R e fazemos gráficos e outras coisas pelo R.

Particularmente vamos voltar a regressão linear explorada no post sobre AIC.

Então temos um preditor e uma resposta, usamos a formula da regressão, resposta é igual ao intercepto mais a inclinação vezes o preditor.

resposta=a+b*preditor

Vamos primeiro simular os dados onde o “a” (intercepto) vai ser 5 e o “b” (inclinação) vai ser 2.
Temos então trinta amostras usando esse modelo e veremos se usando o OpenBUGS vamos conseguir resgatar esses parâmetros dos dados.

Primeiro olhamos o gráfico para ver como estão os dados.

01

Realizamos uma regressão linear no R para ter resultados para comparar.

> modelo1 summary(modelo1) Call: lm(formula = resposta ~ preditor) Residuals: Min 1Q Median 3Q Max -2.6241 -0.9259 0.1517 0.9380 1.8161 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 4.548 1.136 4.002 0.000418 *** preditor 2.093 0.148 14.147 2.8e-14 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 1.183 on 28 degrees of freedom Multiple R-squared: 0.8773, Adjusted R-squared: 0.8729 F-statistic: 200.1 on 1 and 28 DF, p-value: 2.797e-14

E começamos a análise bayesiana.
Primeiro temos que escrever o modelo na linguagem BUGS. Nele vamos declarar o likelihood, nosso a+bx nesse caso, e nosso prior, o que achamos.
Isso é parte da controvérsia na estatística bayesiana, mas aqui vamos simplesmente partir do pressuposto que não existe relação nenhuma, então vamos começar achando que a=0 e b=0, não existe relação nem intercepto.
E vamos deixar também uma coisa legal do BUGS, que é calcular uma estatística derivada. Vamos ver qual a chance de declínio, dos valores de resposta diminuírem (Que nesse exemplo tem que ser pelo menos muito baixa certo, estamos vendo uma relação positiva da resposta com o preditor).

Então definimos nosso modelo em um arquivo chamado linreg.txt.

sink("linreg.txt") cat(" model { # Priori alpha ~ dnorm(0,0.001) beta ~ dnorm(0,0.001) sigma ~ dunif(0, 100) tau <- 1/ (sigma * sigma) # Verossimilhança for (i in 1:n) { y[i] ~ dnorm(mu[i], tau) mu[i] <- alpha + beta*x[i] } # Quantidade derivadas p.decline <- 1-step(beta) # Probabilidade de declinio } ",fill=TRUE) sink()

Feito isso, arrumamos os dados em uma lista, fazemos uma função para inicializar as cadeias e definimos as características do amostrador.
Dai usamos a função bugs() do pacote R2OpenBugs e temos o resultado da nossa regressão usando estatística bayesiana.

> print(modelo1.bayes, dig = 3) Inference for Bugs model at "linreg.txt", Current: 3 chains, each with 1200 iterations (first 200 discarded) Cumulative: n.sims = 3000 iterations saved mean sd 2.5% 25% 50% 75% 97.5% Rhat n.eff alpha 4.503 1.032 2.411 3.843 4.509 5.205 6.585 1.027 250 beta 2.099 0.135 1.829 2.010 2.098 2.186 2.370 1.018 260 p.decline 0.000 0.000 0.000 0.000 0.000 0.000 0.000 1.000 1 deviance 96.296 2.544 93.370 94.397 95.610 97.480 102.602 1.011 260 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 = Dbar-Dhat) pD = 2.8 and DIC = 99.1 DIC is an estimate of expected predictive error (lower deviance is better).

E o resultado fica assim. Ai o alpha é o “a” ou intercepto, deu um valor muito próximo ao encontrado pela estatística frequentista, com um erro um pouquinho menor. Idem para o beta, que é a inclinação, valor com erro um pouquinho menor. E a chance de declínio é 0, não poderia ser diferente certo, vemos uma forte relação positiva. Mas além disso, vemos a distribuição de como ficou esses parâmetros, alias aqui ha uma grande diferença para a interpretação, na estatística esperamos que as estimativas serão as melhores para aqueles dados, a amostra. Aqui essas estimativas esperamos que representem o valores reais, a “verdade”, mas isso uma discussão mais para frente também.

Continuando na parte prática, temos então como e a distribuição, os quantiles de 2.5%, 25%, 50%, 75% e 95% e o Rhat, que é a convergência da cadeia, a principio basta saber que quanto mais próximo de 1, melhor são nossas estimativas.
Mas é isso, temos nossas estimativas do modelo e interpretamos igual o modelo linear frequentista, a inclinação é 2, então temos uma relação verdadeira (significativa) e positiva, a resposta cresce com o preditor.

Mas como no post anterior, vamos simular dados onde não existe relação, a inclinação é 0 e vejamos.

02

Nas estatística frequentista temos:

> summary(modelo2) Call: lm(formula = resposta2 ~ preditor2) Residuals: Min 1Q Median 3Q Max -5.7710 -3.4117 -0.1875 3.3768 7.2168 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 17.3486 3.9411 4.402 0.000142 *** preditor2 0.3303 0.5099 0.648 0.522397 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 3.902 on 28 degrees of freedom Multiple R-squared: 0.01477, Adjusted R-squared: -0.02042 F-statistic: 0.4196 on 1 and 28 DF, p-value: 0.5224

Agora para análise bayesiana, bem como ja temos a função likelihood escrita, podemos usar as mesmas configurações e função inicializadora das cadeias. Temos que montar os novos dados em uma lista e temos…

> print(modelo1.bayes2, dig = 3) Inference for Bugs model at "linreg.txt", Current: 3 chains, each with 1200 iterations (first 200 discarded) Cumulative: n.sims = 3000 iterations saved mean sd 2.5% 25% 50% 75% 97.5% Rhat n.eff alpha 16.566 3.397 9.721 14.280 16.600 18.82 23.031 1.008 410 beta 0.430 0.439 -0.408 0.133 0.424 0.72 1.322 1.012 240 p.decline 0.156 0.363 0.000 0.000 0.000 0.00 1.000 1.017 220 deviance 167.724 2.432 165.000 166.000 167.100 168.90 173.902 1.008 590 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 = Dbar-Dhat) pD = 2.6 and DIC = 170.3 DIC is an estimate of expected predictive error (lower deviance is better).

Um beta (b ou inclinação) perto de 0, mas se olharmos o erro vemos que uma parte dos valores beta encontrado são zero ou ainda negativos. Vejamos a distribuição, ele começa em -0.4, e termina em 1.3. Vejamos isso num gráfico. Vamos olhar o que achávamos antes de olhar os dados (Prior) e o que achamos depois (Posterior)

03

Vejam como continua quase a mesma coisa, a gente esperava a media em zero e ela continua praticamente la. Para entender melhor como continua a mesma coisa, vejamos como fica esse mesmo gráfico do exemplo 1.

04

Olhem como aqui a distribuição posterior esta para um valor positivo, e muito longe de zero.

Bem isso é, só uma pequena olhada no que se trata estatística bayesiana. Na verdade pode paracer complicado escrever aqueles modelos, mas uma sugestão muito boa no livro do Marc Kéry é que nunca se comece do zero, existem muitos modelos que com poucos ajustes e um copia e cola daqui e dali entre modelos e você consegue montar o que precisa. Claro que existe muita discussão sobre o tema de como analisar dados, mas muito além disso, é importante pelos menos dar uma olhada em como tudo funciona para ler mais claramente a literatura e entender nossos pares.

###############################################################
# Regressão linear com estatitica Bayesina #
###############################################################
set.seed(5)
#exemplo1
#simulação dos dados
preditor<-runif(30,5,10)
resposta<-rnorm(30,5+2*preditor)
 
#grafico
plot(resposta~preditor,pch=19,cex=1.6,xlab="Preditor",ylab="Resposta",
frame.plot=F)
abline(lm(resposta~preditor),lwd=2)
 
#ajuste com estatistica frequentista
modelo1<-lm(resposta~preditor)
summary(modelo1)
 
### Analise Bayesiana usando Openbugs
# Escrevendo o modelo de regressão linear
sink("linreg.txt")
cat("
model {
 
# Priori
alpha ~ dnorm(0,0.001)
beta ~ dnorm(0,0.001)
sigma ~ dunif(0, 100)
tau <- 1/ (sigma * sigma)
 
# Verossimilhança
for (i in 1:n) {
y[i] ~ dnorm(mu[i], tau)
mu[i] <- alpha + beta*x[i]
}
 
# Quantidade derivadas
p.decline <- 1-step(beta) # Probabilidade de declinio
}
",fill=TRUE)
sink()
 
# Junte os dados em uma lista
win.data <- list(x=preditor,y=resposta,n=length(resposta))
 
# Função de inicialização
inits <- function(){ list(alpha=rnorm(1), beta=rnorm(1), sigma = rlnorm(1))}
 
# Escolha os parametros que quer estimar
params <- c("alpha","beta", "p.decline")
 
# Caracteristicas do MCMC
nc = 3 #Numero de cadeias
ni=1200 #Tamanho da cadeira
nb=200 #Numero de simulação que serão descartadas
nt=1 #Thinning rate
 
# Inicie o Amostrador
library(R2OpenBUGS)
modelo1.bayes <-bugs(data = win.data, inits = inits, parameters = params,
model = "linreg.txt",n.thin = nt, n.chains = nc,
n.burnin = nb, n.iter = ni)
 
print(modelo1.bayes, dig = 3)
 
#exemplo2
#simulação, note que eu uso a media e desvio da resposta anterior
preditor2<-runif(30,5,10)
resposta2<-rnorm(30,mean(resposta),sd(resposta))
 
#grafico
plot(resposta2~preditor2,pch=19,cex=1.6,xlab="Preditor",ylab="Resposta",
frame.plot=F)
abline(lm(resposta~preditor),lwd=2)
 
#ajuste com estatistica frequentista
modelo2<-lm(resposta2~preditor2)
summary(modelo2)
 
# Agora é so juntar os dados numa lista
win.data2 <- list(x=preditor2,y=resposta2,n=length(resposta2))
 
# O modelo esta pronto, podemos reutilizar o ja escrito
# As outras funções e configurações podem ser igualmente reutilizadas
#então é so iniciar o amostrador
 
modelo2.bayes <- bugs(data = win.data2, inits = inits, parameters = params,
model = "linreg.txt",n.thin = nt, n.chains = nc,
n.burnin = nb, n.iter = ni)
 
#Resultado
print(modelo1.bayes, dig = 3)
 
#grafico posterior vs prior
#exemplo2
curve(dnorm(x,0.430,0.439),-4,4,lwd=2,ylim=c(0,1),main="Exemplo 2",
xlab="Função densidade")
curve(dnorm(x,0,0.001),add=T,col="red",lwd=2)
legend("topleft",col=c("black","red"),
legend=c("Posterior","Prior"),lty=1,lwd=2)
 
#exemplo 1
curve(dnorm(x,2.099,0.135),-4,4,lwd=2,ylim=c(0,1),main="Exemplo 1",
xlab="Função densidade")
curve(dnorm(x,0,0.001),add=T,col="red",lwd=2)
legend("topleft",col=c("black","red"),
legend=c("Posterior","Prior"),lty=1,lwd=2)

Leave a Reply

Your email address will not be published. Required fields are marked *