O Banco de dados encheu :(

Ae galera, infelizmente o blog ta hospedado no biz.nf, que como a maioria dos hosts gratuitos, no plano gratuito so permitem 10 megas de banco de dados.

Eu tenho 9.9 megas atualmente e num cabe mais figura 🙁

Quando eu conseguir algum dinheiro eu pago um host para ter mais espaço e continuar postando. Até la da para deixar a pagina com os posts que já tem.

Qualquer coisa eu estou disponível por e-mail.

GLM com distribuição de Poisson e zeros inflados

Opa, nesse post aqui do ano passado, falamos sobre overdisperion em um glm usando distribuição de poisson. Vimos como ele pode ajudar a levar em conta talvez fatores que não tínhamos em mente enquanto coletávamos dados. Nesse caso a gente contabilizava uma dispersão extra, para mais ou para menos para modelos gerais lineares usando a distribuição de poisson.
Agora a gente vai avaliar um outro caso similar, que tenta resolver o problema de maior dispersão dos dados que esperamos. Agora especificamente devido a muito zeros.

Acho que ter muito zeros é um problema comum em ecologia, principalmente em dados de censo ou abundância. Mas como a gente ja viu a muito tempo atras aqui, modelos com zeros-inflados podem ser uma boa para atacar essa questão.

As vezes eu fico perdido quando e porque as pessoas inventam as coisas. Por exemplo a distribuição de poisson teve o primeiro uso avaliando as contagens de mortes por coices de cavalos no exercito, depois ficou mais comum devido ao uso na engenharia de produção, para prever defeitos, etc.

E aqui a gente usa para conta bixinhos e plantinha. Modelos com zeros inflados eu não tenho ideia de onde ou para que foi criado inicialmente, nem o uso mais comum, mas cai como uma luva para a biologia. Porque a ideia é que podemos ter vários zeros de contagens, mas esses zeros podem ser por dois motivos. Por exemplo se estamos contando aranhas nos arbustos, podem haver arbustos com zero, com uma, com duas aranhas, e assim vai. Ela vão dispersando, e chegando nos arbusto ao acaso, tudo bem.

Mas alguns arbustos podem não ter aranhas não porque não chegou nenhuma aranha, mas porque ele esta num lugar que não é habitável, por algum motivo. Na natureza ocorrem coisas “malucas”, vamos supor que existam lagartas que geram uma nevoa que impede a colonização de aranhas, alias isso a gente não precisa imaginar, existe mesmo hehehe, olhe aqui. Certo mas nesse exemplo então, devem haver zeros devido a presença dessa lagarta, que não deixa as aranhas chegarem, e devem haver zeros que são devido ao fato de nenhuma aranha ter chegado. So que isso faz com que tenhamos uma dispersão maior dos dados.

Lembra-se quando falamos que na distribuição de poisson, só temos um parâmetro, o \lambda. Acontece que essa distribuição vem de um processo, onde eventos independentes acontecem, as aranhas chegando. Mas esses eventos deixam de ser independentes, quanto tem uma coisa fazendo eles não acontecerem. Logo ferrou tudo.
Mas se essas lagartas estão la, então a gente pode supor que existe uma chance \psi delas não estarem no arbusto, por outro lado 1-\psi de ter uma lagartinha no arbusto.
Podemos escrever isso da seguinte forma:

{p(y_i=0)} = \psi+(1-\psi)e^{-\lambda}

Com isso liberamos o processo de poisson desses zeros, e podemos desconsiderar os zeros que não são devido a aranhas não simplesmente terem chegado la.

{p(y_i=n_i)} = {(1-\psi)} \cdot { \lambda^{n_i} e^{-\lambda}\over{n_i!}} para n_i>0

Agora isso pode parecer complicado, mas não é, y_i é simplesmente nossa iesima amostra, que vai ter seu valor de \lambda_i, mas temos a chance \psi de termos um zero ali, que zera tudo. Veja que na segunda equação temos um 1-\psi multiplicando o processo de poison, a função de massa de probabilidade de poisson, no caso extremo, que \psi=1, será 1-1=0, ou seja, se multiplicarmos tudo por zero, a probabilidade sempre é zero, e teremos um pmf sempre zero pro processo de poisson.

Para não confundir muito vamos fazer um exemplo com dados simulados para ver como funciona. Vamos comparar a abundância de aranhas em arbusto com e sem flores, e vamos supor que essas lagartas tem uma chance de 30% de estarem nesses arbustos. Simulamos os dados e obtemos o seguinte:

01

Legal, agora precisamos fazer nosso modelo. Vamos primeiro fazer o modelo bayesiano logo, para entender melhor, e depois olhamos com mais algumas considerações o modelo frequentista.

model {
# Priors
 psi ~ dunif(0,1)
 alpha ~ dnorm(0,0.001)
 beta ~ dnorm(0,0.001)
 
# Likelihood
 for (i in 1:n) {
    w[i] ~ dbern(psi)
    C[i] ~ dpois(eff.lambda[i])
    eff.lambda[i] <- w[i]*lambda[i]
    log(lambda[i]) <- alpha + beta *x[i]
 }

Aqui temos, como no caso do overdispersion, dois processos estocásticos, mas la tinha uma processo de poisson e um processo gaussiano(normal) e aqui temos um processo de poisson e um binomial.
Veja que aqui temos exatamente o que descrevemos la em cima. Temos o processo binomial.

Esse aqui:
{p(y_i=0)} = \psi+(1-\psi)e^{-\lambda}
é essa linha aqui

w[i] ~ dbern(psi)

Dai vemos nosso valor de lambda

log(lambda[i]) <- alpha + beta *x[i]

vemos se ele é um zero binomial ou de poisson

eff.lambda[i] <- w[i]*lambda[i]

e simulamos o processo de poisson.

C[i] ~ dpois(eff.lambda[i])

Veja que essas três linhas é o que colocamos aqui:

{p(y_i=n_i)} = {(1-\psi)} \cdot { \lambda^{n_i} e^{-\lambda}\over{n_i!}} para n_i>0

Aqui estamos usando um prior não informativo, na verdade um prior de zero para tudo. Mas com o modelo pronto, agora a gente separa os dados, faz uma função para inicializar as cadeias, configura as cadeias e roda o modelo. Um detalhe interessante, é que como temos dois processos estocásticos, as coisas começam a ficar mais complicadas e precisamos de uma amostra maior da distribuição posterior para conseguir estabilizar as cadeias. Lembre-se que estamos otimizando dois processos estocásticos, então nada é tão simples como visto aqui.

Mas vamos olhar o resultado:

> print(out, dig = 3) Inference for Bugs model at "ZIP.txt", Current: 3 chains, each with 50000 iterations (first 10000 discarded), n.thin = 4 Cumulative: n.sims = 120000 iterations saved mean sd 2.5% 25% 50% 75% 97.5% Rhat n.eff alpha 0.474 0.200 0.072 0.341 0.478 0.611 0.854 1.001 77000 beta 1.113 0.222 0.688 0.962 1.110 1.262 1.556 1.001 75000 psi 0.722 0.068 0.584 0.676 0.724 0.769 0.848 1.001 120000 deviance 162.205 8.901 146.700 155.700 161.450 167.600 181.600 1.001 51000 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 = 39.6 and DIC = 201.8 DIC is an estimate of expected predictive error (lower deviance is better).

Bem, conseguimos estimar o valor de zeros inflados, com uma média de 72%, veja o valor de psi, bem próximo ao que usamos na simulação, de 70%. E ficamos com os parâmetros lambda para os arbustos com e sem flores.

Podemos olhar melhor o que conseguimos num gráfico.

01

Temos como ficou a distribuição do psi, e os valores de lambda ficaram bem condizentes com o nosso primeiro gráfico, veja que talvez um valor um pouquinho acima da mediana do gráfico, já que estamos calculando ele excluindo os zeros por causa das lagartas. Podemos concluir aqui agora que arbustos com flores tem mais aranhas que arbustos sem flores.

Certo, agora vamos fazer a mesma coisa usando estatística frequentista. O pacote pslc tem a função zeroinfl, que faz o mesmo que fizemos com o bugs, usando máxima verosimilhança.

Primeiro ajustamos o modelo e vemos como ficam as estimativas:

> summary(modelo.z) Call: zeroinfl(formula = C ~ x | 1, dist = "poisson") Pearson residuals: Min 1Q Median 3Q Max -1.2347 -0.9056 -0.1420 0.6217 2.2187 Count model coefficients (poisson with log link): Estimate Std. Error z value Pr(>|z|) (Intercept) 0.4898 0.2015 2.431 0.0151 * xCom Flor 1.1036 0.2231 4.946 7.59e-07 *** Zero-inflation model coefficients (binomial with logit link): Estimate Std. Error z value Pr(>|z|) (Intercept) -0.9776 0.3535 -2.765 0.00569 ** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Number of iterations in BFGS optimization: 8 Log-likelihood: -106.1 on 3 Df

Os valores de alpha e beta ficam parecidos, são o primeiro intercept e o xCom Flor. Depois temos o Zero-inflation model coefficients, temos somente um parametro, o nosso psi, mas que precisa ser convertido pela equação logística para vermos em termos de porcentagem. Mas ele deu o mesmo valor. E não gostaríamos que fosse diferente certo?

Agora para comparar o que acontece quando não levamos em conta os zeros a mais, a inflação de zeros, vamos fazer um glm de poisson.

> summary(modelo.glm) Call: glm(formula = C ~ x, family = "poisson") Deviance Residuals: Min 1Q Median 3Q Max -2.6957 -1.5275 -0.1582 0.6995 2.7414 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 0.1542 0.1690 0.912 0.362 xCom Flor 1.1360 0.1943 5.847 5e-09 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for poisson family taken to be 1) Null deviance: 171.28 on 59 degrees of freedom Residual deviance: 131.37 on 58 degrees of freedom AIC: 253.59 Number of Fisher Scoring iterations: 5

Veja como o valor de intercept e xCom Flor ficaram muito menores. Isso porque os zeros estão puxando eles para baixo, ou seja, acabamos por subestimar os parâmetros que estimamos, e vamos assim predizer erroneamente que devemos ver menos aranhas nos arbustos.

Mas podemos ver como o modelo com zeros inflados tem um valor menor de máxima verossimilhança (melhor ajuste).

> logLik(modelo.z) 'log Lik.' -106.0822 (df=3) > logLik(modelo.glm) 'log Lik.' -124.7955 (df=2)

O que acontece, que ao usar o parâmetro psi, temos um custo menor, a função final explica melhor os dados, mas temos um parâmetro a mais, podemos ponderar isso usando aic, e re-comparar o ajuste.

> extractAIC(modelo.z) [1] 3.0000 218.1645 > extractAIC(modelo.glm) [1] 2.0000 253.5911

E vemos novamente que mesmo com o peso que o AIC da ao número de parâmetros do modelo, temos um valor menor de AIC para o modelo com zeros inflados. Podemos construir um intervalo de confiança aqui para ver se eles são realmente diferentes, e se deveríamos pegar o mais complexo (com zeros inflados) em detrimento do outro, mas eu não sei fazer exatamente esse teste, então vamos olhar outra coisa antes de se preocupar com isso.

Por último ainda damos uma olhada nos resíduos desses dois modelos.

01

Agora olha que legal, veja primeiro que os resíduos de forma geral, estão mais espalhados no glm que no modelo de zero inflados. Outra coisa, veja como os resíduos do índice de 30 a 60, tem alguns valores muitos negativos no glm, que são as barras maiores no histograma. Isso é porque são os zeros que não vem do processo de poisson, mas que estão ali zuando o modelo.
Não é a toa que sempre temos que olhar os resíduos dos modelos que ajustamos, conseguimos sempre encontrar através da analise de resíduos muitos problemas, mas podemos já começar a pensar nas soluções olhando para eles.

Bem esse é o modelo de zero inflados, agora melhor destrinchado que da primeira vez que olhamos para ele. E é legal que ele cai como uma luva para a biologia. É daqui que é derivado aqueles modelos que contabilizam os erros de detecção complexos como vemos aqui.
Para entender aqueles modelos mais complexos, ajuda muito pensar e entender isso aqui.

Outra coisa que achei muito interessante pensando sobre overdispersion e modelos com zeros-inflados, é o fato de como eles podem ajudar a resolver críticas que a gente recebe em manuscritos. Quem nunca enviou um manuscrito para alguma revista e não recebeu uma crítica do tipo. “Você deveria ter coletado tal coisa”, “Você deveria ter (insira aqui uma coisa que não da mais para fazer)”. Poxa podemos ainda tratar certas incertezas com esse tipo de abordagem. Ao invés de jogar tudo fora, sempre tem uma solução (talvez não exatamente sempre, mas quase sempre). E é nessas horas que é bom saber pelo menos que esse tipo de abordagem existe, nunca se sabe quando ela pode salvar a vida. Esse tipo de modelo com zeros inflados pode também ser usado na regressão logística com distribuição binomial.

Bem é isso. o script vai estar la no repositório recologia, se eu escrevi alguma bobeira, algo errado, deixe um comentário corrigindo, ou mande um e-mail. e é isso ai 🙂

Referencia:
Marc Kéry (2010) Introduction to WinBUGS for Ecologists. Academic Press, Burlington.

### Modelo com  Zeros Inflados
set.seed(123)
### Gerando os dados
psi <- 0.7
n.amostras <- 30
x <- gl(n = 2, k = n.amostras, labels = c("Sem Flor", "Com Flor"))
w <- rbinom(n = 2*n.amostras, size = 1, prob = psi)
lambda <- exp(0.69 +(0.92*(as.numeric(x)-1)) )
C <- rpois(n = 2*n.amostras, lambda = w * lambda)
cbind(x, w, C)
 
#Figura 01
boxplot(C~x,frame=F,ylab="Abundância de aranhas",xlab="Arbusto")
 
### Analise Bayesiana
# Definindo o modelo
sink("ZIP.txt")
cat("
model {
# Priors
 psi ~ dunif(0,1)
 alpha ~ dnorm(0,0.001)
 beta ~ dnorm(0,0.001)
 
# Likelihood
 for (i in 1:n) {
    w[i] ~ dbern(psi)
    C[i] ~ dpois(eff.lambda[i])
    eff.lambda[i] <- w[i]*lambda[i]
    log(lambda[i]) <- alpha + beta *x[i]
 }
 
# Quantidade derivada (Para comparar com o resultado frequentista)
 R.lpsi <- logit(1-psi)
}
",fill=TRUE)
sink()
 
# Juntando os dados
bugs.data <- list(C = C, x = as.numeric(x)-1, n = length(x))
 
# Função para gerar o inicio da cadeia
inits <- function(){ list(alpha=rlnorm(1), beta=rlnorm(1), w = rep(1, 2*n.amostras))}
 
# Parametros a Estimar
params <- c("alpha","beta","psi","R.lpsi")
 
# Caracteristicas do MCMC
nc <- 3					# Número de cadeias
ni <- 50000				# Número de amostras da distribuição posterior
nb <- 10000				# Numero de amostras a ser descartadas do inicio da cadeia como burn-in
nt <- 4					# Thinning rate
 
# Iniciando o Gibbs Sampler
library(R2OpenBUGS)
out <- bugs(data=bugs.data, inits=inits, parameters.to.save=params, model.file="ZIP.txt",
n.thin=nt,n.chains=nc,n.burnin=nb, n.iter=ni)
 
print(out, dig = 3)
 
#figura 2
layout(matrix(c(1,1,2,3),ncol=2,nrow=2,byrow=T))
hist(out$sims.list$psi,breaks=seq(0,1,by=0.025),freq=F,xlab="Probabilidade",main="Chance do arbusto ser habitável"
     ,ylab="Frequência")
abline(v=out$summary[3,1],col="red",lty=2,lwd=2)
lines(x=c(out$summary[3,3],out$summary[3,3]), y=c(0,1),lwd=4,col="red")
lines(x=c(out$summary[3,7],out$summary[3,7]), y=c(0,1),lwd=4,col="red")
legend("topleft",lty=c(2,1),lwd=c(2,4),col="red",legend=c("Média","Intervalo de confiança de 95%"),bty="n")
 
 
hist(exp(out$sims.list$alpha),freq=F,breaks=seq(0,8,by=1),main=expression(paste(lambda," Sem flor")),
     xlab="Lambda Esperado",ylab="Frequência",xaxt="n")
axis(1,at=0:8)
hist(exp(out$sims.list$alpha+out$sims.list$beta),freq=F,breaks=seq(0,8,by=1),main=expression(paste(lambda," Com flor")),
     xlab="Lambda Esperado",ylab="Frequência",xaxt="n")
axis(1,at=0:8)
 
### Analise Frequentista
#install.packages("pscl")
library(pscl)
modelo.z <- zeroinfl(C ~ x | 1, dist = "poisson")
summary(modelo.z)
 
modelo.glm <- glm(C ~ x , family = "poisson")
summary(modelo.glm)
 
logLik(modelo.z)
logLik(modelo.glm)
 
extractAIC(modelo.z)
extractAIC(modelo.glm)
 
#Figura 3
layout(matrix(c(1,2,3,4),ncol=2,nrow=2,byrow=T))
 
menor<-min(c(resid(modelo.z),resid(modelo.glm)))
maior<-max(c(resid(modelo.z),resid(modelo.glm)))
 
plot(resid(modelo.z),frame=F,main="Modelo com Zeros-Inflados",ylim=c(menor,maior))
abline(h=0,lty=2,lwd=2)
hist(resid(modelo.z),main="Modelo com Zeros-Inflados",breaks=seq(menor,maior+0.5,by=0.25),freq=F)
lines(density(resid(modelo.z)),lwd=3,lty=3,col="red")
 
plot(resid(modelo.glm),frame=F,main="GLM",ylim=c(menor,maior))
abline(h=0,lty=2,lwd=2)
hist(resid(modelo.glm),main="GLM",breaks=seq(menor,maior+0.5,by=0.25),freq=F)
lines(density(resid(modelo.glm)),lwd=3,lty=3,col="red")

Como diabos funciona um Gibbs Sampler?

Vamos dar uma olhada aqui num algorítimo de Gibbs Sampler. Bem o Openbugs, seu irmão winbugs e o Jags todos usam Gibbs Sampler para ajustar nossos modelos. O que esses programas fazem é fazer um sampler parecido com o que a gente vai ver aqui, so que ele monta todo o código sozinho, somente com o modelo que a gente manda pra ele e dados.
Definitivamente um grande corte de caminho, já que nos permite um bocado de abstração, e se focar no modelo que nos interessa, e em questões numéricas.

Mas nem por isso não vamos deixar de pelo menos entender um pouco o que está acontecendo, ter um pouco de intuição em como funciona o Gibbs Sampler. Apesar de ser um caso especial do Metropolis-Hasting, que a gente viu aqui, aqui e aqui, ele é um pouco mais complicado de entender. Enquanto no Metropolis-Hasting a gente tinha uma chance de aceitar cada novo valor para a distribuição posterior, ou repetir o anterior, aqui a gente sempre aceita o novo valor para a distribuição posterior.

Mas vamos olhar um exemplo funcionando. Primeiro vamos criar uns dados de exemplo. Vamos simular 1000 valores com média -2 e desvio padrão de 2.

Então ficamos com os seguintes dados:

01

Então vamos lembrar o que queremos fazer, queremos inferir a distribuição posterior, saber quais parâmetros realmente geraram os dados a partir dos dados e de um pressuposto que temos, a nossa distribui a priori, nosso prior.

Então queremos:

p(\mu,\sigma \mid Dados) \propto p(Dados \mid \mu,\sigma) \cdot p(\mu,\sigma)

O que são essas coisas?

Queremos saber distribuição posterior p(\mu,\sigma \mid dados), ou seja, quais são os valores de média e desvio padrão que geraram os nossos dados.

Sendo que nos temos os valores de média e desvio padrão que mais provavelmente gerarão os dados que estamos olhando p(Dados \mid \mu,\sigma) , o likelihood, máxima verossimilhança, que vimos até como conseguir ela aqui.

E temos também um presuposto, algo que achamos que deve ser, dado a literatura, os nossa intuição, o que temos p(\mu,\sigma).

Então é isso. Como aqui é um exemplo vamos usar um Prior bem ruim, chutando que a média deveria ser 5 e o desvio 3.162278. Veja que o desvio ta no prior como o tau.prior, e como a variância é 10, então colocar o tau.prior como 0.1 é a mesma coisa que 1/10, ou seja estamos colocando a variância como 10, e a raiz de 10 é 3 e alguma coisa.

# Valores do Prior  #
media.prior<-5
tau.prior<-0.1

Além disso setamos valores iniciais para começar o Gibbs Sampler.

# Valores iniciais #
tau<-1
mtau<-n/2

Criamos também vetores para receber os resultados

nsim<-5000
cadeia.media<-rep(0,nsim)
cadeia.variancia<-rep(0,nsim)

E temos todos os ingredientes, temos os dados, o nosso prior, valores iniciais para iniciar, e agora vamos ao Gibbs Sampler propriamente dito.

for (i in 1:nsim) {
 
  v<-1/(tau*n + tau.prior)
  m<-v*(tau*sum(dados)+tau.prior*media.prior)
  cadeia.media[i]<-rnorm(1,m,sqrt(v))
 
  tau<-rgamma(1,mtau,sum((dados-cadeia.media[i])^2)/2)
  cadeia.variancia[i]<-1/tau
 
}

E acabou, 5 linha de código, e se olharmos o resultado:

O valor estimado para média é de:

> mean(cadeia.media[1501:nsim]) [1] -1.999777

Com uma precisão:

> sd(cadeia.media[1501:nsim]) [1] 0.06566074

E no vaso de desvio padrão:

> sqrt(mean(cadeia.variancia[1501:nsim])) [1] 2.039191

Com uma precisão:

> sd(cadeia.variancia[1501:nsim]) [1] 0.1877103

Realmente funciona bem, já que os dados foram gerados com média -2 e desvio padrão 2, os distribuição posterior está bem próximo, englobando os valores reais.

Mas a pergunta é, o que aconteceu durante aquele loop que gerou esses resultados, o que o Gibbs Sampler fez?

Realmente é difícil entender aquele monte de multiplicação e somas, pelo menos pra mim é complicado. Mas vamos olhar um pouco melhor as variáveis, primeiro vamos ver como é a relação entre elas no grafo.

Vamos olhar ele e o script ao mesmo tempo.

01

v<-1/(tau*n + tau.prior)
m<-v*(tau*sum(dados)+tau.prior*media.prior)

No grafo, as variáveis são vértices, e as arestas aponta a direção em que uma variável influencia outra, se o vértice não recebe nenhuma aresta, ele não é influenciado por ninguém.

Veja que começamos gerando valores para v e m, para gerar esses valores, usamos os dados, os priores e o valor de tau. Com exceção de tau, todos são valores “estáticos”. Eles não mudam, veja que no grafo, só tem setinhas saindo deles a indo para os valores que eles influenciam, mas eles não recebem nenhuma setinha.

cadeia.media[i]<-rnorm(1,m,sqrt(v))

Agora esses dois valores são usados para gerar uma amostra da distribuição posterior, que vai ser guardado nesse vetor chamado cadeia.media.

tau<-rgamma(1,mtau,sum((dados-cadeia.media[i])^2)/2)

Depois usamos ele, mais um parâmetro externo chamado mtau, que depois vamos comentar, mas que serve para regularmos os tamanhos dos passos que vamos dar na caminhada aleatória do Gibbs Sampler, mas foco na hora que usamos a distribuição gamma para gerar um novo valor de tau, usamos o valor recém adicionado a distribuição posterior, usamos ele para calcular os mínimos quadrados, que é o que? Uma função de custo, se o nosso valor novo é “bom”, perto do real, a soma aqui sera algo mais próximo de zero, quanto pior o novo valor adicionado, mais sera o valor aqui, mas por enquanto vamos continuar olhando o grafo e se ater ao fato que ele é usado para estabelecer o novo valor de tau.

cadeia.variancia[i]<-1/tau

Usamos tau também para estabelecer o novo valor de variância, nosso outro parâmetro da distribuição normal. E assim acabou uma iteração, e agora começamos de novo.

Mas o detalhe é olhar esse ciclo que é formado aqui.

01

Um valor da distribuição posterior é usado para gerar o próximo valor de tau, que é usado para gerar os próximos valores de v e m que são usados para gerar o próximo valor da distribuição posterior, que é usado para gerar o próximo valor de tau, e assim por quantas iterações quisermos, para o tamanho da cadeia de mcmc que definimos.

Agora pensando um pouco mais nos valores, vamos olhar os primeiros 50 valores de tau, v e m.

01

Veja que quando o valor da média para o próximo valor que vamos sortear sobe, a valor de tau desce, quando tau desce o valor de m sobe, e assim vai.
Se a gente olhar com mais carinho essa linha

tau<-rgamma(1,mtau,sum((dados-cadeia.media[i])^2)/2)

Pense que se estamos com um valor muito bom para o parâmetro média, essa soma dará um valor muito próximo de zero, e a implicação disso ao se gerar um número aleatório para tau é o seguinte:

01

Ou seja, próximo de zero, deixamos tau ter a possibilidade de ser grande, já se o valor é grande, pegamos um valor muito ruim para a média, deixamos o valor sempre pequeno, próximo de zero.

Agora ao olhar a linha que gera o valor que será usado para propor a próxima amostra para a distribuição posterior da média a gente entende isso.

m<-v*(tau*sum(dados)+tau.prior*media.prior)

Veja que se, ao extremo, tau for zero, ele tira a influencia dos dados, do likelihood e deixa os valores de prior comandando o que deve ser a próxima amostra da média. O contrario, se tau é um valor grande, ele aumenta a influencia do likelihood no que será o a próxima amostra da média. E lembre-se que o valor do posterior é algo entre a máxima verossimilhança e o nosso prior. E assim estamos caminhando, entre esses valores.
Outra coisa legal a se notar é que temos uma multiplicação pelo n, e na linha anterior

v<-1/(tau*n + tau.prior)

Veja que o tau é multiplicado pelo n, e n é o tamanho da amostra, e daqui sai a nossa influencia da quantidade de dados, no que será a distribuição posterior, lembra quando comentamos que quanto mais dados, menor a importância relativa do prior, em algum post que não lembro mais hehe?

Mas é isso ai, assim funciona o Gibbs sampler, sempre caminhando sem parar, diferente do Metropolis-Hastings.

Podemos olhar que o resultado e ver como a convergência foi rápida

01

Não é a toa, usamos um valor gigantesco de amostra, 1000 amostras.
E veja como os valores de média e desvio que conseguimos como a distribuição posterior ficam próximos ao valor real, mais precisamente, o intervalo de confiança de 95% engloba o valor real de média e desvio padrão usado para produzir os dados.

01

E esse é o Gibbs Sampler, que tanto falamos. E o que o OpenBugs faz quando a gente chama ele? Ele constrói um algorítimo desse tipo, mas com o modelo que definimos no model. Mas a idéia geral é por aqui.

É um bocado mais complicado que o Metropolis-Hasting, principalmente até captar o que ele está fazendo com os dados. Eu precisei ficar rodando ele com menos dados, e mudando vários dos parâmetros para conseguir entender mais ou menos o que está acontecendo e espero não ter entendido nada errado, por favor me corrijam se eu estou falando alguma abobrinha.

Outra coisa, se olharem no script original, do Brian Neelon, verão que ele tem 2 parâmetros, a e b, que deixa em zero. Eu suprimi esses parâmetros aqui, mas veja que eles influenciam nos dois números que servem de parâmetro para a função rgamma que faz o novo valor de tau a cada iteração. Isso possibilita o controle do tamanhos dos passos que damos, e vai influenciar na distribuição posterior. Esse tipo de coisa pode ser otimizada para valores que dão soluções melhores, mas por enquanto, como mal entendemos o Gibbs Sampler, não vamos mexer neles. Mas vendo o algorítimo eu acho que ajuda a desmistificar um pouco e nos deixar mais confiantes com a estatística bayesiana, ja que não tem nada muito milagroso aqui.

Bem vamos ficar por aqui, que esse post ja está meio grande, mas o script está no repositório recologia e no topo dele existe um link para o script original da onde eu copie o Gibbs Sampler mostrado aqui.

# Original:  http://people.duke.edu/~neelo003/r/Normal%20Mean%20Gibbs.r
# Gibbs sampler para uma distribuição normal
###################################
set.seed(250)
 
########
# Data #
########
n<-1000
media.real<--2
desvio.padrao.real<-2
dados<-rnorm(n,media.real,desvio.padrao.real)
 
hist(dados,main="Histograma",ylab="Frequência",xlab="Valores")
abline(v=mean(dados),lty=2,col="red",lwd=2)
 
#####################
# Valores do Prior  #
#####################
media.prior<-5
tau.prior<-0.1
 
####################
# Valores iniciais #
####################
tau<-1
mtau<-n/2
 
####################################
# Vetor para guardar os resultados #
####################################
 
nsim<-5000
cadeia.media<-rep(0,nsim)
cadeia.variancia<-rep(0,nsim)
 
#################
# Gibbs Sampler #
#################
vc<-rep(0,nsim)
mc<-rep(0,nsim)
tauc<-rep(0,nsim)
 
 
for (i in 1:nsim) {
 
  v<-1/(tau*n + tau.prior)
  m<-v*(tau*sum(dados)+tau.prior*media.prior)
  cadeia.media[i]<-rnorm(1,m,sqrt(v))
 
  tau<-rgamma(1,mtau,sum((dados-cadeia.media[i])^2)/2)
  cadeia.variancia[i]<-1/tau
 
}
 
library(igraph)
grafo<-graph.formula(n-+v,tau-+v,tau.prior-+v,v-+m,dados-+m,tau.prior-+m,media.prior-+m,m-+cadeia.media,v-+cadeia.media,
                   cadeia.media-+tau,dados-+tau,tau-+cadeia.variancia,tau-+m)
tkplot(grafo,vertex.size=60,vertex.color="gray",vertex.label.color="black",edge.arrow.size=2,edge.color="black",
       vertex.label.cex=1.2)
 
grafo<-graph.formula(m-+cadeia.media,v-+cadeia.media,cadeia.media-+tau,tau-+cadeia.variancia,tau-+m,tau-+v)
tkplot(grafo,vertex.size=60,vertex.color="gray",vertex.label.color="black",edge.arrow.size=2,edge.color="black",
       vertex.label.cex=1.2)
 
layout(matrix(c(1,2,3,4),nrow=2,ncol=2,byrow=T))
curve(dgamma(x,shape=1,rate=0.01),0,10,main="Rate = 0.01",frame=F)
curve(dgamma(x,shape=1,rate=0.1),0,10,main="Rate = 0.1",frame=F)
curve(dgamma(x,shape=1,rate=1),0,10,main="Rate = 1",frame=F)
curve(dgamma(x,shape=1,rate=10),0,10,main="Rate = 10",frame=F)
 
##########################
# Distribuição Posterior #
##########################
mean(cadeia.media[1501:nsim])
sd(cadeia.media[1501:nsim])
 
sqrt(mean(cadeia.variancia[1501:nsim]))
sd(cadeia.variancia[1501:nsim])
 
 
#Evolução da cadeia
plot(1:nsim,cadeia.media,type="l",col="red",frame=F)
abline(h=mean(cadeia.media),lwd=2,lty=2)
abline(v=1500,lwd=4,lty=3,col="gray")
text(750,-1.76,"Burn-in")
 
#Distribuição posterior
plot(cadeia.media[1:500],sqrt(cadeia.variancia[1:500]),type="b",frame=F,xlab="Média",ylab="Desvio Padrão",cex=0.5,pch=19)
library(MASS)
contornos<-kde2d(cadeia.media[1:500],sqrt(cadeia.variancia[1:500]))
contour(contornos,add=T,lwd=4,nlevels = 6,col=rev(heat.colors(6)))

Encontrando a máxima verossimilhança analiticamente

Em 2012, a gente comentou o que é e como encontrar a máxima verossimilhança de uma distribuição binomial aqui, agora vamos olhar com mais carinho o que fizemos la e tentar entender um pouco melhor como esse negócio funciona de vedadinha.

Para uma única observação da distribuição binomial, a verossimilhança (ou likelihood) de k ocorrências de um total de N tem uma probabilidade p (e.g.k arbustos com aranhas de um total de N arbustos ou k passarinhos parasitados de um total N de passarinhos observados) é:

\mbox{Prob}(k|p,N)=\binom{N}{k} p^k (1-p)^{N-k}

Pra quem não sabe o que \binom{N}{k} é, de uma olhada em combinação.

Mas continuando, agora se a gente tem n observações, cada uma com o mesmo número N de arbustos (no caso das aranhas, ou pode ser n observações com N passarinhos em cada amostra, ou o que diabos for do seu interesse), e o número de arbustos habitados por aranhas da i^{\mbox{th}} observação é k_i, então a verosimilhança (ou likelihood) vai ser:

{\cal L} = {\prod_{i=1}^n} {\binom{N}{k_i} p^{k_i} (1-p)^{N-k_i}}
Quem nunca viu esse simbolo {\prod_{i=1}^n}, ele é a mesma coisa que somatório, mas ao invés de somar a gente multiplica as coisas, ou seja ele é o simbolo do produtório, lembre-se que a chance de 2 coisas acontecerem é uma coisa vezes a outra, lembre-se dos “e” ou “ou” da estatística.

Nos já comentamos o porque as pessoas transformam a verosimilhança em log (log-likelihood), que segundos as regras de log, uma multiplicação vira uma soma, além de ter que lidar com números com mais casas decimais, permitindo maior precisão, mas no caso do log-likelihood, não sei como escreve em português, se log-verossimilhança ou verossimilhança logarítmica, se alguém puder me corrigir seria legal, mas o que temos é:

 L=\sum_{i=1}^n (\log \binom{N}{k_i} + k_i \log p + (N-k_i) \log (1-p))

No R a gente escreve isso ai em cima assim:

sum(dbinom(k,size=N,prob=p,log=TRUE))

Agora a gente vê formulas pra caramba em ciência, faz muitas contas que simplesmente sabemos que é. E dificilmente descobrimos o porque é, pegue por exemplo quando a gente calcula a prevalência, o total de casos ocorrido em todas as observações pelo total de casos avaliado em todo o experimento de algum parasita.

Agora vamos olhar algo legal relativo a equação ali de cima, partindo de uma abordagem analítica (leia calculo ^^)

Nesse caso nos podemos resolver o problema analiticamente, diferenciando e equação a partir do p e igualando a derivada a zero. Existe também esse pequeno post aqui sobre derivadas no R. Então a gente não precisa necessariamente saber fazer essa conta, mas é legal o que se desenrola daqui.

Primeiro lembre-se como se parece o gráfico da verossimilhança de uma distribuição binomial.

.

Veja que ele pode começar em baixo, vai chegar uma hora que sobe, numa barriguinha e depois vai descer.

Se você sair na calculando derivadas e desenhando retas tangente assim.

01

O topo da barriguinha do gráfico da verossimilhança da distribuição binomial vai ter uma reta paralela ao eixo x, ou seja a inclinação vai ser zero, ou seja com derivada zero.

Então suponha que \hat p é a máxima verossimilhança da estimativa do parâmetro p que satisfaça:

 \frac{dL}{dp} =  \frac{d\sum_{i=1}^n \left(  \log \binom{N}{k_i} + k_i \log p + (N-k_i) \log (1-p) \right)}{dp}  =  0.

Como a derivativa das soma é igual a soma das derivativas, temos:

  \sum_{i=1}^n \frac{d \log \binom{N}{k_i}}{dp} +  \sum_{i=1}^n \frac{d k_i \log p}{dp} +  \sum_{i=1}^n \frac{d (N-k_i) \log (1-p)}{dp} = 0

O termo \log \binom{N}{k_i} é uma constante em relação ao p, logo sua derivativa é zero, e assim o primeiro termo desaparece.

Como k_i e (N-k_i) são fatores constantes, eles saem das derivativas e a equação então se torna:

  \sum_{i=1}^n k_i \frac{d \log p}{dp} +  \sum_{i=1}^n (N-k_i) \frac{d  \log (1-p)}{dp} = 0.

A derivativa de \log p é 1/p, então a regra da cadeia diz que a derivada de

\log (1-p) é

d(\log (1-p))/d(1-p) \cdot d(1-p)/dp = -1/(1-p).

Vamos escrever o valor particular de p que queremos como \hat p.

Então…

 \frac{1}{\hat p} \sum_{i=1}^n k_i - \frac{1}{1-\hat p} \sum_{i=1}^n (N-k_i)  =  0

Jogamos aquele pedacinho negativo pro outro lado, ou como é o correto que aprendi no Khan Academy, somamos o que está depois do igual dos dois lados.
 \frac{1}{\hat p} \sum_{i=1}^n k_i  =  \frac{1}{1-\hat p} \sum_{i=1}^n (N-k_i)

Dai a gente troca os 1-\hat p e \hat p que estão embaixo do 1 de lado…
 (1-\hat p) \sum_{i=1}^n k_i  =  \hat p \sum_{i=1}^n (N-k_i)

Depois aplica a distributiva ali no 1-\hat p,
 \sum_{i=1}^n k_i - \hat p \sum_{i=1}^n k_i =  \hat p \sum_{i=1}^n (N-k_i)

Passa pro outro lado
 \sum_{i=1}^n k_i =  \hat p \sum_{i=1}^n (N-k_i) + \hat p \sum_{i=1}^n k_i

E coloca o \hat p em evidência…
 \sum_{i=1}^n k_i  =  \hat p \left( \sum_{i=1}^n k_i + \sum_{i=1}^n (N-k_i) \right)

Agora note que aqui,\sum_{i=1}^n k_i + \sum_{i=1}^n (N-k_i) , a gente ta somando os casos que aconteceram, k, com o total de casos N, menos os casos que acontecerem k, poxa, isso simplismente N.

 \sum_{i=1}^n k_i = \hat p \sum_{i=1}^n N

Sem desistir de entender o que está acontencedo, a somatoria \hat p \sum_{i=1}^n N é o numero de observações com N amostras, ou  \hat p n N , depois vamos ainda passar tudo que multiplica o \hat p pro outro lado dividindo….
 \sum_{i=1}^n k_i  =  \hat p n N

E acabamos com…
 \hat p  =  \frac{\sum_{i=1}^n k_i}{nN}

Então a estimativa da máxima verissimilhança \hat p, é apenas o fração do total ocorrências de algo, \sum k_i, dividido por pela soma de todas as observações, nN.

Por exemplo o total de arbustos com aranhas pelo total de arbustos vistoriados, ou o total de passarinhos parasitados dividido pelo total de passarinhos verificados.

Parece um bocado de esforço para provar o obvio, que a melhor estimativa do parasitismo é meramente a frequência de casos de parasitismos, mas a gente começa a ver a relação daquela continha que sempre fazemos, com as curvas que likelihood que gostamos de desenhar.

Outras distribuições, como a distribuição do poisson se comportam de forma similiar. Se a gente diferenciar a verissimilhança (likelihood) ou a verossimilhança logaritimizada (log-likelihood) e resolver para a estimativa da máxima verossimilhança, nos conseguimos a resposta.
Com um pouco mais de esforço (ja que a distribuição normal tem dois parâmetros, média e desvio padrão), podemos fazer o mesmo para média e desvio padrão da distribuição normal (também conhecida com Gaussiana).

Para a distribuição de poisson, a estimativa do parâmetro de frequência \hat \lambda é igual a média das contagens observadas por amostra.

Para a distribuição normal, com dois parâmetros \mu e \sigma^2, nos temos que calcular as derivativas parciais da verossimilhança para os dois parâmetros, e resolver as duas equações. (\partial L/\partial \mu =\partial L/\partial \sigma^2 = 0).

E como a gente disse, sem conta nenhuma, a resposta obvia daquilo que a gente calculou a vida inteira.
\hat \mu=\bar x (a estimativa é a média observada nas amostras) e a variância é \hat{\sigma^2}=\sum (x_i-\bar x)^2/n (a estimativa da variância é a variância das amostras)

Um pequeno detalhe, é que a máxima verossimilhança é uma estimativa tendenciosa da variância, dividindo a soma dos quadrados \sum (x_i-\bar x)^2 por n-1 ao invés de n é o usado comumente. Mas a explicação fica para a próxima hehe.

Para outras distribuições simples como a binomial negativa, e mais alguns problemas complexos, não há uma solução analítica simples assim, como a que vimos aqui. E temos que aplicar soluções numéricas para encontrar a máxima verossimilhança para parâmetros das distribuições e modelos que temos interesse. Para praticamente tudo, existem programas prontos, mas o ponto legal do calculo visto aqui é criar alguma intuição do que está acontecendo em casos simples, e que as estimativas fazem sentido quando pensamos nas distribuições, ou quando olhamos um histograma.

Bem eu ainda não tenho os culhões para deduzir tudo isso sozinho, essa conta tem no livro Ecological models and Data in R, que é um dos livros mais legais que ja li, apesar de não entender 3/4 do que está escrito, foi o primeiro livro que vi tentando explicar estatística e modelagem, de uma perspectiva mais biológica que bolinhas pretas e brancas em potinhos, alias eu sou um grande fa do autor, o Ben Bolker, uma das imagens do cabeçalho do blog é a foto dele, se prestar atenção nesse nome, direto vai encontrar ele espalhado pelos pacotes do R. Além de tudo acho um link legal para ver onde o calculo entra na bio-estatística, e como a gente usa calculo muito mais rotineiramente que pensamos, gostamos dele ou não. E ainda bem que a distribuição binomial tem somente um parâmetro, para facilitar entender a derivação, que é bem complicada, mas com algum esforço da para entender.

Bem ficamos por aqui hoje. Feliz ano novo a todos que passarem por aqui, e que tenhamos um 2014 melhor que 2013.

Referência:

Bolker, B. 2008 – Ecological Models and Data in R 408pp Princeton University Press