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)))

Overdispersion em um modelo para contagens (Poisson)

Ola, após quase um mês enrolado, o blog está de volta. O ultimo host nos baniu rapidinho hehehe e ainda no meio das provas.
Mas agora vamos ser no que da, achei um novo host sugerido no Reddit, o biz.nf. Espero que pelo menos no fim de ano o blog fique ativo.

Mas chega de blablabla. Vamos continuar onde paramos no último post quando falamos de GLM usando distribuição de Poisson.
Um dos problemas, não é bem um problema, mais um característica da distribuição de poisson é que ele so tem um parâmetro, o \lambda (lambda), ele determina onde está a média e a variação é igual também a esse valor. Daqui a alguns post espero escrever sobre como é o processo de poisson, mas essas características, como falamos muitas vezes, fazem todo o sentido se pensar que quanto contamos poucas coisas, erramos por pouquinho, agora quando contamos muito erramos por muito. Porém as vezes , quando estamos pegando dados reais, muitas vezes a variação pode ser maior ou menor que isso, por exemplo se existe algum fator que não estamos levando em conta, isso faz com que a variação seja diferente de \lambda, apesar de média estar certa. Isso vai levar a gente a fazer previsões erradas, ja que se parar para pensar um momento, se a variação a mais ou a menos não for levada em conta, a gente pode ter confiança demais num valor, ou de menos.

Vamos tentar ilustrar com um exemplo.
Imagine que visitamos várias fazendas e reservas, e sempre parávamos em alguns locais e contávamos quantos coelhos víamos e ficamos com a suspeita de que coelhos são mais comuns onde a terra está arada, foi trabalhada do que em locais nativos. Então aqui esta um box-plot de 20 contagens para cada um níveis que estamos considerando (2 níveis, arado e nativo), 20 para cada, da um total de 40 amostras.

Uma tabelinha assim:

Local Contagem 1 Nativo 0 2 Nativo 2 3 Nativo 3 . . . 39 Arado 2 40 Arado 2

O que nos da um gráfico assim:

Coelhos

Certo. A primeira coisa que vem a cabeça é um glm com distribuição de Poisson, mas agora podemos ou não pensar que há uma maior variação nos dados.

Bem vamos fazer os dois usando o comando GLM do R e ver o que acontece.

> summary(glm.fit.no.OD) Call: glm(formula = C.OD ~ x, family = poisson) Deviance Residuals: Min 1Q Median 3Q Max -2.00000 -0.78339 -0.02871 0.65787 2.27670 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 0.6931 0.1581 4.384 1.17e-05 *** xArado 0.4220 0.2035 2.074 0.0381 * --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for poisson family taken to be 1) Null deviance: 54.580 on 39 degrees of freedom Residual deviance: 50.182 on 38 degrees of freedom AIC: 154.28 Number of Fisher Scoring iterations: 5

Esse primeiro caso ta considerando o mesmo que fizemos no post teste t de poisson. Vemos as estimativas de \lambda e um valor p significativo, legal.

> summary(glm.fit.with.OD) Call: glm(formula = C.OD ~ x, family = quasipoisson) Deviance Residuals: Min 1Q Median 3Q Max -2.00000 -0.78339 -0.02871 0.65787 2.27670 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 0.6931 0.1757 3.945 0.000332 *** xArado 0.4220 0.2261 1.867 0.069681 . --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for quasipoisson family taken to be 1.234685) Null deviance: 54.580 on 39 degrees of freedom Residual deviance: 50.182 on 38 degrees of freedom AIC: NA Number of Fisher Scoring iterations: 5

E agora, olhando o número vermelhinho ali destacado, o valor p não é mais significativo.
Mas o que mudou? Se a gente olhar as estimativas de \lambda aqui, são iguais ao modelo inicial. Mas veja que os erros(“Std. Error”) do valor de \lambda são maiores, agora estamos levando em conta uma gama maior de valores, por isso a menor certeza da diferença, e a diferença nos valores p. Se pensar em valor p como a chance de observamos novamente um valor tão extremo, claro que uma variação maior diminui nossas certezas.
Mas veja também que temos um outro valor diferente ainda, o “Dispersion parameter”. Agora é maior que 1. Não é difícil chutar que ele ta maior porque estamos levando em conta a maior variação nos dados.

Bem podemos repetir essa análise usando uma análise bayesiana com o Openbugs para melhorar nossa intuição sobre o assunto.

Então vamos olhar como fica o modelinho na linguagem bugs.

sink("overdispersion.txt")
cat("
model {
# Priors
 alpha ~ dnorm(0,0.001)
 beta ~ dnorm(0,0.001)
 sigma ~ dunif(0, 10)
 tau <- 1 / (sigma * sigma)
 
# Likelihood
 for (i in 1:n) {
    C.OD[i] ~ dpois(lambda[i])
    log(lambda[i]) <- alpha + beta *x[i] + eps[i]
    eps[i] ~ dnorm(0, tau)
 }
}
",fill=TRUE)
sink()

Veja que o modelo é igual a todos que temos feito, a não ser por uma adição, que é o parâmetro de dispersão eps, que esta contabilizando a “overdispersion”.

Muito cuido na leitura a partir daqui, que eu posso ter entendido as coisas erradas, mas vamos indo.

Bem, para diferenciar os lugares temos:

C.OD[i] ~ dpois(lambda[i])
log(lambda[i]) <- alpha + beta *x[i]

Ou seja, cada local depende dos seus parâmetros locais, que é o \alpha(alpha), e o \beta(beta), lembrando que o x vai ser 0 para locais arado e 1 para nativo. Com isso construímos o \lambda.
Mas ai adicionamos a maior variação no \lambda

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

Além disso tudo, temos mais alguma coisa. Que coisa? A coisa que influência a quantidade de coelho e que não conseguimos controlar, ou seja a quantidade de coelhos observada não é um processo de Poisson perfeito, temos outras coisinhas, ou coisonas acontecendo durantes as observações que não controlamos. Talvez a quantidade de alimento disponível, presença de predadores, a ação desses predadores. O mundo é mais complexo que estamos descrevendo nesse modelinho. Muitas coisas acontecendo ao mesmo tempo. So que não temos como contabilizar tudo isso, não há tempo viável para ficar vasculhando o local para disponibilidade de alimento, nem ficar vários dias la olhando o que aparece de predador, como eles agem, como são todas as interações. E tem outra, os dados ja estão coletados, não vamos pular do barco e não analisar os dados so porque não existem outras 5 milhões de variáveis disponíveis, tempo que sempre tentar fazer o melhor possível com o que temos a mão, e não reclamar do que não temos.

Mas não é por causa disso que precisamos “mentir” para o nosso modelo, podemos dizer, modelo esta acontecendo alguma coisa (essa coisa é o eps ali), que sempre acontece mas eu não tehho controle, leve isso em consideração ok? E ficamos assim:

log(lambda[i]) <- alpha + beta *x[i] + eps[i]
eps[i] ~ dnorm(0, tau)

Veja que temos um erro de distribuição normal, com média zero e um desvio padrão sigma, e a variância que chamamos de tau. É isso que tamos fazendo. Se for um processo perfeito, esperamos que esse valor seja próximo de zero, se não esperamos que esse valor contabilize essa maior, ou menor dispersão dos dados.

Mas depois disso separamos os dados numa lista no formato certo para mandar para o bugs, Fazemos uma função para iniciar valores para os parâmetros de interesse. configuramos como queremos a cadeia e pimba, temos o resultado.

> print(out, dig = 3) Inference for Bugs model at "overdispersion.txt", Current: 3 chains, each with 3000 iterations (first 1000 discarded), n.thin = 5 Cumulative: n.sims = 6000 iterations saved mean sd 2.5% 25% 50% 75% 97.5% Rhat n.eff alpha 0.625 0.185 0.224 0.506 0.635 0.755 0.956 1.001 6000 beta 0.431 0.236 -0.022 0.274 0.427 0.586 0.916 1.001 6000 sigma 0.308 0.168 0.029 0.181 0.300 0.420 0.655 1.001 5300 deviance 144.011 7.084 129.600 139.000 144.600 149.600 155.700 1.001 6000 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 = 10.3 and DIC = 154.3 DIC is an estimate of expected predictive error (lower deviance is better).

E aqui estamos, as estimativas do \lambda são bem próximas daquelas usando o comando GLM. A conclusão também, assim como o valor p não era significativo, se a gente olhar o \beta e sua distribuição, vemos que o intervalo de 95% vai de -0.022 a 0.916. Ou seja o zero, ou zero de diferença está dentro daquilo que consideramos como verdadeiro, então estamos inclinados a que tanto o locais arados como nativos tem a mesma quantidade de coelhinhos, tem um \lambda igual.

Mas agora como pensar nesse desvio, nesse “Overdispersion”?

Eu não tenho muita certeza desses gráficos, por causa do kernel que usei, mas é algo mais ou menos assim, vamos colocar esses valores num gráfico e olhar para locais nativos e arados.

Arado

Nativo

Bem veja a diferença entre como é o ajuste com e sem overdispersion, veja que o nos dois casos, apesar de não dar para ver direito para locais arado, a distribuição azul, que é quando consideramos o “overdispersion” é mais abertinha, que a linha vermelha, considerando somente a variação esperada no processo de Poisson, pode parecer pouca diferença, mas veja que fez uma grande diferença nas nossas conclusões, sendo em um caso com valor p menor que 0.05 e no outro não, ao abrir as perninhas, a distribuição fica mais sobrepostas, uma na outra, entre arado e nativo, fazendo com que começamos a levar mais em consideração a afirmação que elas são iaguais, e que a diferença que estamos vendo é ao acaso. Mas a idéia é essa, estamos assumindo que existe algo a mais na variação, correlações que não levamos em conta, algo a mais, e estamos deixando isso aparecer no modelo, e somente então tomando nossa conclusão.

Para confirmar que estamos no caminho certo, podemos gerar um conjunto de dados sem Overdispersion e ver o que acontece. Usaremos os mesmo parâmetros alpha e beta, so vamos tirar a maior variação no lambda e usar o mesmo modelo que usamos no caso acima.
Primeiro geramos os dados e vemos como fica.

exemplo2

Ai fazemos tudo denovo, juntar dados e tal e rodamos o openbugs denovo para esse novo conjunto de dados.

> print(out2, dig = 3) Inference for Bugs model at "overdispersion.txt", Current: 3 chains, each with 3000 iterations (first 1000 discarded), n.thin = 5 Cumulative: n.sims = 6000 iterations saved mean sd 2.5% 25% 50% 75% 97.5% Rhat n.eff alpha 0.335 0.113 0.111 0.261 0.335 0.413 0.554 1.001 3500 beta 0.604 0.139 0.325 0.513 0.604 0.696 0.881 1.001 6000 sigma 0.171 0.101 0.009 0.095 0.161 0.238 0.387 1.120 49 deviance 397.125 8.296 377.100 392.600 399.000 403.100 408.900 1.006 410 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 = 10.4 and DIC = 407.5 DIC is an estimate of expected predictive error (lower deviance is better).

Agora vamos olhar o que aconteceu com o valor de sigma. Veja agora que ele é baixinho, comparado com o caso anterior. E ele vai de 0.009 a 0.387. Bem temos praticamente um zero dentro da amostra, ja que 0.009 está bem próximo de zero. Esse valor na verdade é 0.00852687, mas ta arredondado para 3 casas depois do zero, a gente sempre imprime o resultado print(out2, dig = 3), 3 dígitos de precisão.

Mas parece que essa idéa funciona, e é bem interessante para produzir modelos melhores. E se o sigma for algo perto de zero, a gente pode excluir ele do modelo e fazer um modelo mais simples. E assim começamos a tentar ver um mundo um pouco maior além dos glm convencionais. Tudo isso tem pronto no R como vimos ao usar family=quasipoison no R, então uma vez entendendo, vale a pena considerar a possibilidade de uso, dependendo do que estamos considerando.

Num próximo post, vamos estender essa ideia para os modelos com zero-inflados, que ja fizemos um post aqui a milhões de anos atras, mas agora vamos voltar nele com estatística bayesiana.

Bem ficamos por aqui, legal ter o blog devolta, espero que esse host não chute minha bunda também como os anteriores, deixei até uma propagandinha ali em cima na lateral do blog pra ver se faço algum dinheiro para pagar um servidor hehe.

O script vai estar la no repositório recologia do github além de aqui embaixo, e é isso ai, até a próxima.

### Gerando os dados
set.seed(123)
n.locais <- 20
x <- gl(n = 2, k = n.locais, labels = c("Nativo", "Arado"))
eps <- rnorm(2*n.locais, mean = 0, sd = 0.5)
lambda.OD <- exp(0.60 +(0.35*(as.numeric(x)-1) + eps) )
 
C.OD <- rpois(n = 2*n.locais, lambda = lambda.OD)
 
data.frame(Local=x,Contagem=C.OD)
 
#GRafico 1
boxplot(C.OD ~ x, col = "grey", xlab = "Uso da terra",frame=F,
ylab = "Contagem de coelhos", las = 1, ylim = c(0, max(C.OD)))
 
 
### Analise usando o R
glm.fit.no.OD <- glm(C.OD ~ x, family = poisson)
glm.fit.with.OD <- glm(C.OD ~ x, family = quasipoisson)
summary(glm.fit.no.OD)
summary(glm.fit.with.OD)
 
anova(glm.fit.no.OD, test = "Chisq")
anova(glm.fit.with.OD, test = "F")
 
### Analise com OpenBugs
 
sink("overdispersion.txt")
cat("
model {
# Priors
 alpha ~ dnorm(0,0.001)
 beta ~ dnorm(0,0.001)
 sigma ~ dunif(0, 10)
 tau <- 1 / (sigma * sigma)
 
# Likelihood
 for (i in 1:n) {
    C.OD[i] ~ dpois(lambda[i])
    log(lambda[i]) <- alpha + beta *x[i] + eps[i]
    eps[i] ~ dnorm(0, tau)
 }
}
",fill=TRUE)
sink()
 
# Juntando os dados
bugs.data <- list(C.OD = C.OD, x = as.numeric(x)-1, n = length(x))
 
# Função de inicialização
inits <- function(){ list(alpha=rlnorm(1), beta=rlnorm(1), sigma = rlnorm(1))}
 
# Parametros a serem guardados na saida.
params <- c("alpha", "beta", "sigma")
 
# Configurando o MCMC
nc <- 3		# Número de cadeias
ni <- 3000	# Tamanho da cadeia
nb <- 1000	# burn-in, quanto do inicio da cadeia vai ser jogado fora
nt <- 5		# Thinning rate (A cada uma amostra que pegamos, jogamos 5 fora, para evitar auto-correlação)
 
# Iniciando o Gibbs sampling
library(R2OpenBUGS)
out <- bugs(data=bugs.data, inits=inits, parameters.to.save=params,model.file="overdispersion.txt", n.thin=nt,
            n.chains=nc,n.burnin=nb, n.iter=ni)
print(out, dig = 3)
 
 
#O que realmente estamos estimando.
 
#Grafico 2
hist((C.OD[which(x=='Arado')]),freq=F,main="Arado",xlab="Histograma de Contagens",ylab="Frequência da contagem")
lines(density(rpois(1000,out$summary[1,1]),adjust = 4),lwd=2,col="red",lty=2)
valores<-rep(out$summary[1,1],1000)+rnorm(1000,0,out$summary[3,1])
valores[valores<0]<-0
lines(density(rpois(1000,valores),adjust = 3),lwd=2,col="blue",lty=2)
legend("topright",col=c("red","blue"),lty=2,lwd=2,legend=c("Sem overdispersion","Com overdispersion"),bty="n")
 
#Grafico 3
hist((C.OD[which(x=='Nativo')]),freq=F,main="Nativo",xlab="Histograma de Contagens",ylab="Frequêcia da contagem")
lines(density(rpois(1000,out$summary[1,1]+out$summary[2,1]),adjust = 3),lwd=2,col="red",lty=2)
valores<-rep(out$summary[1,1]+out$summary[2,1],1000)+rnorm(1000,0,out$summary[3,1])
valores[valores<0]<-0
lines(density(rpois(1000,valores),adjust = 3),lwd=2,col="blue",lty=2)
legend("topright",col=c("red","blue"),lty=2,lwd=2,legend=c("Sem overdispersion","Com overdispersion"),bty="n")
 
 
#Exemplo 2
n.locais <- 60
x <- gl(n = 2, k = n.locais, labels = c("Nativo", "Arado"))
lambda.OD <- exp(0.60 +(0.35*(as.numeric(x)-1)) )
 
C.OD <- rpois(n = 2*n.locais, lambda = lambda.OD)
 
boxplot(C.OD ~ x, col = "grey", xlab = "Uso da terra",frame=F,
ylab = "Contagem de coelhos", las = 1, ylim = c(0, max(C.OD)))
 
bugs.data2 <- list(C.OD = C.OD, x = as.numeric(x)-1, n = length(x))
 
params <- c("alpha", "beta", "sigma")
 
out2 <- bugs(data=bugs.data2, inits=inits, parameters.to.save=params,model.file="overdispersion.txt", n.thin=nt,
            n.chains=nc,n.burnin=nb, n.iter=ni)
print(out2, dig = 3)

Ancova Binomial Bayesiana avaliando Prevalências

Esses nomes de posts são de certa forma uma sacanagem. Mas normalmente é sempre assim que a gente vê as coisas.
Análises normalmente estão “empacotadas” em nomes, teste t, diferença de médias, com mesma variação, sem mesma variação, anova, ancova, anova de medidas repetidas , regressão linear e etc.
Eu entendo que capítulos de livros ou posts precisam de títulos, índices, mas a gente tem que tomar sempre cuidado para não focar demais num nome e não olhar o todo, se você focar apenas em um ponto você pode perder todo o esplendor do por do Sol, como diria Bruce Lee em Operação Dragão (Isso que é citação heim).

Após o blablabla inicial, vamos voltar ao GLM binomial que estávamos fazendo aqui.

Certo, no caso anterior estávamos contando e vendo a proporção, ou a palavra mais usada para parasitismos que é prevalência, como ela muda entre duas espécies.

Apesar de estarmos simulando dados, vamos pensar em algo legal. Um exemplo que eu acho legal é o Pardal ou Passer domesticus, esse que infesta as cidades e a gente vê um bocado. A distribuição geográfica era originalmente, a maior parte da Europa e da Ásia, mas ele foi introduzido na América, África, Austrália e Nova Zelândia, e olha que esses dois últimos são ilhas e eles não chegaram voando eu acho. Mas por algum motivo sua população tem decaído na europa. Inclusive o jornal The Independent oferece 5 mil euros ou libras, não sei reconhecer o simbolo do dinheiro que eles usam, mas então aqui tem a reportagem que eles abrem essa oportunidade. Ja houveram três entradas sérias para ganhar o prêmio, onde alguns caras publicaram artigos em revistas científicas e colocaram suas teorias para explicar o que esta acontecendo. Mas até agora nada. Bem apesar de a gente não ligar muito pro Pardal, são 5 mil euros no bolso, pode ser interessante tentar entender o que ta acontecendo. Mas como a gente estava pensando em parasitismos e prevalências, vamos continuar pensando em exemplos assim, uma ideia interessante é ver se o mesmo parasita observado anteriormente, o Plasmodium que causa a Malaria Aviaria não pode ter um dedo nesse caso. Como assim?

Uma teoria interessante é essa:

immunityinvasion

Espécies tem sucesso ao invadir novos ambientes, porque elas trazem seus parasitinhas, e ja vivem bem com eles, mas ai esses parasitas podem atacar novas espécies ao chegarem, ai os parasitas rebentam todo mundo do local enquanto a espécie nova procria e aumenta sua população tranquila, lembra portugueses chegando ao Brasil e os índios morrendo não? Mas essa é a idea.

Acontece que os pardais estão se ferrando no seu local de origem. Tudo é bem complicado de se encaixar. Mas suponha que tenhamos os seguintes dados:

> dados status local peso 1 1 Mata Atlantica 79 2 1 Mata Atlantica 25 3 1 Mata Atlantica 68 . . . 277 1 Mata Atlantica 23 278 1 Mata Atlantica 19 279 1 Amazonia 62 280 1 Amazonia 28 281 1 Amazonia 38 . . . 496 1 Amazonia 20 497 1 Amazonia 50 498 1 Amazonia 15 499 1 Cerrado 50 500 1 Cerrado 34 . . . 726 1 Cerrado 53 727 1 Cerrado 26 728 1 Cerrado 78

A gente tem pego alguns pardais aqui no Brasil, cada linha é um passarinho diferente, e vendo toda essa discussão, a gente resolveu dar uma olhadinha aqui, será que a prevalência dos pardais varia aqui dentro do Brasil, o que a gente ve no cerrado é igual na amazônia ou na mata atlântica?

Então temos três locais, e nesses locais a gente foi colocando redes pra capturar passarinhos, ai pegávamos eles, tirávamos sangue de um por um, pesavamos e soltavamos, depois a gente pega o sangue, joga num frasco, poe um primer que amplifica a sequência do plasmodium, faz uma pcr e olha se amplificou, se amplificou ta doente, status é 1, ta doente, tem o parasita, caso contrario é zero.

Esses são os dados resultantes desse processo todo. A primeira coisa a notar, é que não temos muito como controlar o número de amostras, o número de passarinhos em cada caso, claro que precisamos de um minimo, se pegar apenas um ou nenhum pardal num da para querer ver muita coisa, mas as populações variam, a chance de captura depende de muitas coisas, então o número de pardais que teremos vai variar entre locais, não tem muito jeito, daria para homogeniezar isso jogando fora amostras, mas ai nos teríamos que se basear no pior caso, no local que temos menos amostras para ver quantas amostras usar, e convenhamos, faz pouco sentido ou nenhum jogar informação fora. Mas para nossa satisfação é possível usar a distribuição binomial para descrever a prevalência, ou chance e não teremos problemas com a variação do número de passarinhos por amostra.

Agora aqui temos uma diferença forte em relação ao que estávamos fazendo no post anterior. Cada individuo tem uma característica somente sua, o peso, então não da para simplesmente ver a contagem total de passarinhos com o parasita e sem de acordo com o local.

Bem antes de começar a pensar em como fazer o modelo, vamos olhar melhor os dados.

01

Primeiro vamos olhar quantos passarinhos estão parasitados e não parasitados de forma geral, usando todos os passarinhos.
Aqui parece que temos meio a meio, mais ou menos.

Mas vamos olhar se isso continua do mesmo jeito por local.

02

Nossa, agora virou uma zona, na Amazônia parece que é meio a meio, agora no cerrado parece que tem muito mais indivíduos parasitados que não parasitados, e pra piorar tudo parece que a mata atlântica apresenta um padrão oposto.
Mas ainda tem o peso dos passarinhos na jogada.

Colocando num gráfico vemos

03

Bem como so tem 0 ou 1, os pontinhos ficam no teto do gráfico ou no chão. Mas a gente pode fazer uma retinha chamada Loess, que segue a seguinte ideia, pensa que a gente pega os 5 primeiro pontinhos, sendo que eles estão organizados do menor pro maior peso, que podem ser 0 ou 1, somamos e dividimos por 5, ai temos um pontinho, uma prevalência desses 5 pontinhos, ai a gente faz isso com os próximos 5 pontinhos e fazemos mais um pontinho, e depois denovo e denovo, e juntamos uma linha nesses pontinhos, essa é a linha vermelha, que nesse caso ajuda a visualizar como a prevalência esta mudando de acordo com o peso. E vemos uma mudança nos três locais, mais mais forte no cerrado aparentemente, menos nos outros dois casos, mas mudança de acordo com o peso. Não é difícil imaginar motivos que ocasionariam tal tendência, maior peso pode significar indivíduos mais velhos, então que ficaram mais tempo expostos ao ambiente com o parasita, existem ainda a possibilidade do sistema imuno piorar com a idade, muitas possíveis causas, mas antes de cogitar isso, temos que verificar que isso esta acontecendo, e parece que é o caso.

Então parece que podemos tentar pensar da seguinte forma:

04

O status, o estado relativo a presença do parasita depende do peso e do local, mas essas variáveis interagem, uma influência a outra, como sugerido no segundo e terceiro gráfico, que mostra padrões diferentes de acordo com o local e vice versa.

Então fazemos nosso modelinho para mandar para o OpenBugs.

model {
# Priors
 for (i in 1:n.local) {
    alpha[i] ~ dnorm(0, 0.01)  # Interceptos (Um pra cada local)
    beta[i] ~ dnorm(0, 0.01)   # Inclinações (Um pra cada local)
 }
 
# Likelihood
 for (i in 1:n) {
    status[i] ~ dbin(p[i], 1)
    logit(p[i]) <- alpha[local[i]] + beta[local[i]]* peso[i]
 }

Olha, primeiro estabelecemos nossos prior, vamos assumir nenhum efeito e zero parasitismos com esses prior. Veja que essas distribuições do jeito que estão são um pico sobre o zero.
Outra coisa aqui é que precisamos de três equações da reta, uma para cada local, então temos três interceptos e três inclinações, três a+bx.

Agora o likelihood fica assim, o status vem de uma distribuição binomial, dai o uso do status[i] ~ dbin(p[i], 1), como é 0 ou 1, temos o dbin com a probabilidade do indivíduo, e para 1 indivíduo, os dois parâmetros dentro do dbin.
Dai a segunda linha é a parte determinística do modelo, que é a+bx, aqui alpha[local[i]] + beta[local[i]]* peso[i], agora o alpha e o beta, o a e o b dependem do local que o passarinho esta e do seu peso. E acabou o modelo. Como estamos trabalhando com a distribuição binomial, ou seja um GLM, temos que transformar essa probabilidade com a função logística, veja que tem uma função pronta no bugs pra isso para facilitar a vida, isso é a parte escrita logit(p[i]).

Ai antes de fechar a chave do modelo calculamos mais alguns parametros de interesse.

# Recuperando os efeitos baseados no caso base "Amazonia"
 a.effe2 <- alpha[2] - alpha[1]		# Intercepto Amazonia vs Cerrado
 a.effe3 <- alpha[3] - alpha[1]		# Intercepto Amazonia vs Mata Atlantica
 b.effe2 <- beta[2] - beta[1]		# Inclinação Amazonia vs Cerrado
 b.effe3 <- beta[3] - beta[1]		# Inclinação Amazonia vs Mata Atlantica
 
# Teste Personalizado (Respondendo a minha pergunta)
 test1 <- alpha[2] - alpha[1]		# Diferença entre Cerrado a Amazonia
 test2 <- alpha[2] - alpha[3]		# Diferença entre Cerrado e Mata Atlantica
}

Antes de fechar a chave do modelo, recuperamos os efeitos baseados num intercepto geral, para comparar com os valores usando estatística frequentista, usando o comando glm do R e a última parte eu faço meus testes que tenho interesse, que é ver se a prevalência muda do cerrado para a Amazônia e Mata Atlântica, afinal antes de querer saber porque, eu tenho que saber que ha mudança, se não houver mudança a pergunta porque não faz sentido, que é o meu interesse a principio. Isso seria algo como comparações entre grupos, mas eu não vou sair comparando tudo com tudo que nem louco, eu quero saber isso a principio e não tentar todas as comparações possíveis.

Pronto, agora juntamos os dados numa lista.

bugs.data <- list(status = dados$status, local = as.numeric(dados$local), n.local = 3, n = nrow(dados),peso=dados$peso)

Lembrando que mandamos os locais como números, então temos que lembrar depois quem é o 1, quem é 2 e quem é três, como isso vem de um fator do R, a não ser que declaremos uma ordem, ela vai ser a alfabética ou lexografica, então é Amazônia, que começa com A é 1, Cerrado é 2 e Mata atlântica é 3. Simples.

Outro detalhe é que o peso foi Uniformizado, meramente isso é diminuir cada valor da média e dividir pelo desvio padrão, isso nos deixa com os valores variando em unidades de desvio padrão e com média zero. As tendencias vão continuar as mesmas com ou sem essa transformação, mas podem haver problemas se o bugs tiver que lidar com números muitos grandes, por isso ta assim, talvez não haja diferença aqui, ja que os pesos so tem duas casas decimais, mas transformamos de qualquer forma.

Fazemos uma função para gerar valores para o inicio da cadeia ao acaso:

> # Função para iniciar a cadeia > inits <- function(){ list(alpha = rlnorm(3, 3, 1), beta = rlnorm(3, 2, 1))} # Note log-normal inits > #Ela so gera valores ao acaso para iniciar as cadeias, veja so > inits() $alpha [1] 6.996297 22.373664 6.808293 $beta [1] 12.366757 6.475181 5.738183

Olha o que isso faz, gera 3 interceptos e 3 inclinações para iniciar a cadeia, esses valores tem que ser razoáveis, mas como vamos queimar o início da cadeia não são muito relevantes, mas como no caso em que calculamos o pi, o começo sempre é ruim por isso é bom descartar, mas temos que começar de algum lugar, senão não da.

Depois disso falamos também quais parâmetros vamos guardar e configuramos o MCMC.

Mandamos tudo do R para o OpenBugs e recebemos de volta resultado.

> print(out,digits=3) Inference for Bugs model at "ancovabin.txt", Current: 3 chains, each with 1500 iterations (first 500 discarded), n.thin = 5 Cumulative: n.sims = 3000 iterations saved mean sd 2.5% 25% 50% 75% 97.5% Rhat n.eff alpha[1] -0.054 0.136 -0.325 -0.141 -0.051 0.036 0.216 1.001 3000 alpha[2] 2.205 0.274 1.731 2.011 2.194 2.372 2.789 1.001 3000 alpha[3] -0.521 0.128 -0.765 -0.609 -0.522 -0.435 -0.272 1.001 3000 beta[1] 0.286 0.137 0.017 0.191 0.283 0.378 0.553 1.003 870 beta[2] 1.526 0.288 0.982 1.330 1.519 1.712 2.130 1.001 3000 beta[3] 0.392 0.126 0.146 0.304 0.392 0.478 0.638 1.001 3000 a.effe2 2.258 0.303 1.695 2.053 2.242 2.456 2.878 1.001 3000 a.effe3 -0.468 0.187 -0.845 -0.595 -0.466 -0.340 -0.101 1.001 3000 b.effe2 1.240 0.321 0.637 1.024 1.234 1.452 1.888 1.003 3000 b.effe3 0.107 0.187 -0.258 -0.018 0.108 0.234 0.469 1.003 860 test1 2.258 0.303 1.695 2.053 2.242 2.456 2.878 1.001 3000 test2 2.726 0.302 2.164 2.523 2.713 2.913 3.363 1.001 2400 deviance 817.706 3.492 812.900 815.100 817.100 819.600 826.202 1.001 3000 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 = 6.0 and DIC = 823.7 DIC is an estimate of expected predictive error (lower deviance is better).

E temos ai nosso resultado, o intercepto e inclinação de cada caso, cada local, parâmetros para desenhar uma função logística para cada quadro da figura 3, so pensar em uma função logística desenhada em cada quadradinho ali.
Outra coisa é que recuperamos os coeficientes baseados no intercepto geral, o caso base.

Se verificarmos, para ter certeza que a análise funciona, vemos que recuperamos os valores originais usados para fazer a simulação dos dados.

> beta.vec [1] 0.0 2.0 -0.5 0.5 1.0 0.0

alpha[1]=-0.054 é o intercepto básico, o primeiro valor ali em cima, e a.effe2=2.258 e a.effe3=-0.468 que são as diferença de intercepto, estimativas muito boas, não estamos concluindo nada errado até aqui.
Dai temos as inclinações beta[1]=0.286 e as interações b.effe2=1.240 e b.effe3=0.107, que são razoáveis, apesar de termos subestimado um pouco talvez a inclinação base, ja que a média deu 0.286, mas o valor real usado na simulação de 0.5 está dentro do intervalo de confiança de 95% que vai de 0.017 até 0.553, so olhar ali o beta[1].

Usando o comando glm do R podemos ainda fazer a mesma analise num framework frequentista e temos:

> summary(glm(status ~ local * peso, family = binomial,data=dados)) Call: glm(formula = status ~ local * peso, family = binomial, data = dados) Deviance Residuals: Min 1Q Median 3Q Max -3.0142 -1.0129 0.2111 1.0985 1.7024 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -0.04842 0.13629 -0.355 0.722375 localCerrado 2.21404 0.30130 7.348 2.01e-13 *** localMata Atlantica -0.46942 0.18620 -2.521 0.011699 * peso 0.27947 0.13866 2.015 0.043858 * localCerrado:peso 1.20495 0.31480 3.828 0.000129 *** localMata Atlantica:peso 0.11294 0.18557 0.609 0.542798 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: 999.97 on 727 degrees of freedom Residual deviance: 811.67 on 722 degrees of freedom AIC: 823.67 Number of Fisher Scoring iterations: 6

Estimativas muito semelhantes dos parâmetros, aqui teríamos que fazer as contas para extrair os parâmetros de cada local, mas como a gente viu no modelo bayesiano, são contas de mais e menos, então não é muito difícil. E temos que chegar a mesma conclusão de qualquer forma, senão algo estaria errado.

Mas após verificar os resíduos e validar o modelo, damos uma olhada no resultado dos nossos testes específicos.

05

E vemos uma diferença significativa, ou seja é pouco provável que não haja diferença, veja como a distribuição das diferenças esta muito longe do zero, ou seja podemos afirmar que existe uma diferença entre a prevalência do cerrado e os outros locais. Agora porque? Porque será que isso deveria acontecer? Eu agora começaria a sequenciar os Plasmodiuns e veria se eles são muito diferentes entre os locais, isso poderia ajudar a pensar se o que esta acontecendo tem haver com o ambiente, ou com a adaptação do parasita as espécies, algo assim, mas são dados simulados, a gente só queria fazer a análise pra ver como era 🙂

Bem é isso, analise de covariância binomial, covariância porque tem fatores contínuos (peso) e fatores discretos (locais). Binomial que os dados que tinha era contagem de 0 e 1, e é um caso de modelo geral linar (glm) ja que eu lianerizo uma relação que não é linear com a função logística, pra estimar as coisas como se fossem uma reta, usando a distribuição binomial, então glm binomial.

O script inteiro vai estar no repositório Recologia no Github e aqui em baixo. Até a próxima.

###############################################################
#    ANCOVA Binomial                                          #
###############################################################
#Gerando dados para um exemplo
set.seed(51)
quantidade<-round(runif(3, 200, 300))
peso<-round(runif(sum(quantidade),10,80))
hist(peso)
peso<-(peso-mean(peso))/sd(peso)
hist(peso)
 
#Preditores
dados<-data.frame(local=factor(rep(c("Mata Atlantica", "Amazonia", "Cerrado"),quantidade)),peso)
dados
 
Xmat <- model.matrix(~ local*peso,data=dados)   #Matriz de Efeitos
print(Xmat, dig = 2)
beta.vec <- c(0, 2, -0.5, 0.5, 1, 0)            #Valor dos parametros
 
lin.pred <- Xmat[,] %*% beta.vec	        # Valor do preditor Linear
exp.p <- exp(lin.pred) / (1 + exp(lin.pred))    # Chance esperada
exp.p
 
#Resposta
status <- rbinom(n = nrow(dados), size = 1, prob = exp.p) # Adicionando erro binomial
status				                          # Checando o que foi criado
 
#Juntando tudo num data.frame
dados<-data.frame(status,dados)
 
##################################################################
 
#Olhando os dados
dados
 
#Grafico dos Diagnósticos
#jpeg("01.jpg")
barplot(table(dados$status),ylim=c(0,500),xlab="Diagnóstico",ylab="Número de passarinhos")
#dev.off()
#Grafico dos diagnósticos por locais
#jpeg("02.jpg")
barplot(table(dados$status,dados$local),ylim=c(0,300),xlab="Diagnóstico",ylab="Número de passarinhos",beside=T,
        legend.text = c("Não parasitado","Parasitado"))
#dev.off()
#Grafico dos diagnósticos por locais e peso.
library(lattice)
 
panel.smoother <- function(x, y) {
  panel.xyplot(x, y,pch=19,col="black")             # Plota os pontos
  panel.loess(x, y,span=0.6,lwd=2,lty=2,col="red")  # Mostra a linha que muda de acordo com a média
}
 
#jpeg("03.jpg")
xyplot(status~peso|local,data=dados,panel=panel.smoother,layout = c(1, 3))
#dev.off()
 
#Grafo do modelo
library(igraph)
#jpeg("04.jpg")
plot(graph.formula( Local -+ Status , Peso -+ Status , Local ++ Peso ),vertex.size=25,edge.color="black",
     vertex.color="white")
#dev.off()
 
### Analise usando OpenBugs
# Definindo um modelo
sink("ancovabin.txt")
cat("
model {
 
# Priors
 for (i in 1:n.local) {
    alpha[i] ~ dnorm(0, 0.01)  # Interceptos (Um pra cada local)
    beta[i] ~ dnorm(0, 0.01)   # Inclinações (Um pra cada local)
 }
 
# Likelihood
 for (i in 1:n) {
    status[i] ~ dbin(p[i], 1)
    logit(p[i]) <- alpha[local[i]] + beta[local[i]]* peso[i]
 }
 
# Recuperando os efeitos baseados no caso base "Amazonia"
 a.effe2 <- alpha[2] - alpha[1]		# Intercepto Amazonia vs Cerrado
 a.effe3 <- alpha[3] - alpha[1]		# Intercept Amazonia vs Mata Atlantica
 b.effe2 <- beta[2] - beta[1]		# Inclinação Amazonia vs Cerrado
 b.effe3 <- beta[3] - beta[1]		# Inclinação Amazonia vs Mata Atlantica
 
# Teste Personalizado (Respondendo a minha pergunta)
 test1 <- alpha[2] - alpha[1]		# Diferença entre Cerrado a Amazonia
 test2 <- alpha[2] - alpha[3]		# Diferença entre Cerrado e Mata Atlantica
}
",fill=TRUE)
sink()
 
# Juntando os dados numa lista (Veja que o local a gente manda como um número apenas,
#vc tem que lembrar que número é o que)
bugs.data <- list(status = dados$status, local = as.numeric(dados$local), n.local = 3, n = nrow(dados),peso=dados$peso)
 
# Função para iniciar a cadeia
inits <- function(){ list(alpha = rlnorm(3, 3, 1), beta = rlnorm(3, 2, 1))} # Note log-normal inits
#Ela so gera valores ao acaso para iniciar as cadeias, veja so
inits()
 
# Parametros que queremos salvar os resultados
params <- c("alpha", "beta","a.effe2","a.effe3","b.effe2","b.effe3","test1","test2")
 
# Configurações do MCMC
ni <- 1500
nb <- 500
nt <- 5
nc <- 3
 
# Iniciando o Gibbs Sampler
library(R2OpenBUGS)
out <-bugs(data = bugs.data, inits = inits, parameters.to.save = params, model.file = "ancovabin.txt",
           n.thin = nt, n.chains = nc, n.burnin = nb, n.iter = ni)
 
print(out,digits=3)
beta.vec
summary(glm(status ~ local * peso, family = binomial,data=dados))
 
str(out)
 
#jpeg("05.jpg")
par(mfrow=c(2,1))
hist(out$sims.list$test1,xlim=c(0,max(out$sims.list$test2)),main="Diferença entre Cerrado a Amazonia",
     xlab="Tamanho da Diferença")
lines(c(0,0),c(0,400),lty=2,lwd=2,col="red")
 
hist(out$sims.list$test2,xlim=c(0,max(out$sims.list$test2)),main="Diferença entre Cerrado e Mata Atlantica",
     xlab="Tamanho da Diferença")
lines(c(0,0),c(0,400),lty=2,lwd=2,col="red")
#dev.off()

Teste t Binomial – GLM com prevalências.

A gente ja usou um bocado estatística bayesiana usando Openbugs. Os modelos feitos no Openbugs são facilmente replicados no Winbugs, Jags e com um pouquinho de dedicação no Stan. Mas a maioria dos posts anteriores era apenas usando a distribuição normal para descrever os dados.

Mas as vezes não é possível descrever o mundo assim. Um exemplo é quando a gente trabalha com algo como a prevalência de um parasita.
Vamos pegar um exemplo aqui, passarinhos sendo parasitados por Haemoproteus columbae.

Haemoproteus é um genero de protozoa que parasita passarinhos, repteis e anfíbios. O seu nome é derivado do grego, “Haima” = Sangue e
“Proteus” = Deus dos mares que tem o poder de assumir varias formas diferentes. Dai da para sacar que que ele vive no sangue e tem mais de uma forma.

vol102(3)_foto_capa2

Então temos passarinhos, e eles podem ou não estar contaminado. Mas como descobrimos?

Bem existem dois jeitos, ou a gente faz uma lâmina do sangue, prepara ela e olha no microscópio e fica procurando o Haemoproteus ou a gente prepara o sangue, faz uma pcr colocando primers que replicam o DNA do haemoproteus mas não do passarinho hospedeiro (devem vir varias células, logo vários DNAS diferentes num pouquinho de sangue) e se houver bastante DNA amplificado no final, bingo, o parasita tem que estar ali.

Mas no final das contas o que vamos falar é o seguinte, vamos ter sangue de vários hospedeiros, e nos dois casos vamos apenas falar se o indivíduo, o passarinho esta ou não parasitado pelo Haemoproteus.

Individuo rep.dados.Sp..each...15. Parasita 1 1 Columbina talpacoti Sim 2 2 Columbina talpacoti Sim 3 3 Columbina talpacoti Sim 4 4 Columbina talpacoti Sim 5 5 Columbina talpacoti Não . . . 25 10 Scardafella squammata Sim 26 11 Scardafella squammata Não 27 12 Scardafella squammata Não 28 13 Scardafella squammata Sim 29 14 Scardafella squammata Não 30 15 Scardafella squammata Não . . .

Acho que da para entender como estão os dados, temos uma identidade para cada indivíduo, cada amostra, ou seja um passarinho é uma amostra. Cada amostra pertence a uma espécie, e sabemos se ele esta ou não parasitada, temos algo como Sim ou Não, parasitado ou não parasitado, como quiser chamar.

Bem como são sim ou não, normalmente as pessoas resumem essa informação em contagens de quantos passarinhos estão parasitados e quantos não estão parasitados.

> tabela Parasita Sp Sim Não Columbina talpacoti 32 30 Scardafella squammata 10 47

Aqui mantemos a informação de quantos passarinhos foram amostrados (soma de todos as células da matriz), quantos de cada espécie, quantos parasitados ou não por espécie, notem que com essa matriz podemos reconstruir a planilha la em cima, isso se não nos importamos com a identidade de cada passarinho para outras informações, mas temos aqui a mesma informação que la em cima.

Virou uma tabelinha assim, a primeira coisa que vem a cabeça é teste de chi-quadrado para ver se temos um desvio nas proporções. Como assim? Uai se o parasita ocorre numa frequência constante, independente da espécie, então esperaríamos uma proporção igual entre as linhas e colunas, agora se algo acontece, uma espécie é mais parasitada ou, sei la, a prevalência é baixinha nas espécies, ou a gente pode observar se o desvio das contagens é devido ao acaso ou há razões para desconfiarmos que algo a mais acontece.

Então realizamos o teste de chi-quadrado:

> teste.chi <- chisq.test(tabela) > teste.chi Pearson's Chi-squared test with Yates' continuity correction data: tabela X-squared = 13.6387, df = 1, p-value = 0.0002216 > teste.chi$observed # Contagens observadas Parasita Sp Sim Não Columbina talpacoti 32 30 Scardafella squammata 10 47 > teste.chi$expected # Contagens esperada sobre a hipotese nula Parasita Sp Sim Não Columbina talpacoti 21.88235 40.11765 Scardafella squammata 20.11765 36.88235

E temos o clássico p muito menor que 0.05, e mais um monte de valor. Bem, como não fornecemos as probabilidades para o esperado, ele calcula automaticamente, mas podemos ver o que esperávamos no resultado olhando o $expected do modelo ajustado. Bem eu ja fiquei perdido nesses números. Quando a gente olha por exemplo o resultado de cruzamentos de Drosophila melanogaster em relação a algumas mutações, é fácil estabelecer mentalmente o que é o esperado, 1:1 em caso do sexo, 3:1 no caso de uma mutação dominante e assim vai, mas essa pergunta pode ser mais complicada para alguns casos. Então a pergunta, “o que você espera de resultado” é sempre uma pergunta interessante. Mas de qualquer forma aqui a gente vê que algo acontece, agora vamos olhar mais a frente.

Outra forma de olhar os dados é calcular a prevalência, que é, para uma espécie nesse exemplo, o quanto de casos de parasitados tivemos do total de dados.

> dados Sp N Prevalence simul 1 Columbina talpacoti 62 0.516 32 2 Scardafella squammata 57 0.193 10

E temos a prevalência, que é mais ou menos a chance de um indivíduo estar parasitado. E Para nossa surpresa, pelo menos pra minha, dividir o número de indivíduos parasitados pelo total de indivíduos amostrados da a máxima verosimilhança (Maximum likelihood) da distribuição Binomial, mas no próximo post vamos olhar isso.
Mas antes de continuar, vamos ver o que estamos querendo falar ao dizer que 0.516, ou 52% é a prevalência do Haemoproteus na espécie Columbina talpacoti

01

Bem, são dados de contagem aqui, mas temos que ter em mente que essa contagem tem um limite, que no primeiro caso é 62 e no segundo 57. Deixamos riscado no gráfico esse limite.
Agora para o gráfico 1, onde esta a bolinha mais alta? Oras, no 32, o maximum likelihood, mas existe um erro inerente a qualquer amostragem, mas para essa mesma prevalência, ainda seria razoável pensar nela como o valor real de prevalência, a real chance de um passarinho ser contaminado mesmo com 31 passarinhos doentes. Agora pense no valor de 30 ou 29, ja estamos ficando com metade da altura que tem o maximum likelihood, aqui talvez não seja muito legal pensar em 51.6% como prevalência, e pensar que a prevalência é um pouco mais baixa. O mesmo vale para outra especie.

Mas então temos crenças diferentes para diferentes valores de prevalência para cada uma das especies. Mas além de estabelecer quanto é a prevalência de cada espécie, podemos tentar estimar a diferença entre essas prevalência. Que seria o que? Algo como o quanto deve ser a diferença entre elas.

So que aqui tem um esquema, meramente falando em diferença, tudo fica meio confuso, a gente não esta meramente diminuindo 0.516 de 0.193. A gente quer saber a distribuição, assim como essas que estimamos para cada especie, a distribuição das diferenças, porque não vamos simplismente falar um valor, ele pode estar errado, a gente tem que ver todas as possibilidades de diferenças.

E a gente pode fazer isso como similar a forma como fizemos o teste t aqui, so que para usar a distribuição binomial, precisamos de uma função de ligação, que nesse caso é a função logistica.

Vamos ver primeiro como é a função na linguagem Bugs:

sink("Binomial.t.test.txt")
cat("
model {
 
# Priors
 alpha ~ dnorm(0,0.01)
 beta ~ dnorm(0,0.01)
 
# Likelihood
 for (i in 1:n) {
    C[i] ~ dbin(p[i], N[i])
    logit(p[i]) <- alpha + beta *species[i]
 }
# Derived quantities
 prev.C_talpacoti <- exp(alpha) / (1 + exp(alpha))
 prev.S_squammata <- exp(alpha + beta) / (1 + exp(alpha + beta))
 prev.Diff <- prev.C_talpacoti - prev.S_squammata	# Teste
}
",fill=TRUE)
sink()

A gente ja fez quase isso aqui, sem usar MCMC. Mas agora a diferenças que a gente tem, ao invés de apenas 1 parâmetro, um valor de probabilidade (que aqui a gente le prevalência), a gente tem dois, uma para uma espécie e outro para outra espécie.

Então fazemos aqui temos que fazer isso para as duas espécies.

Então o likelihood tem um loop que vai de 1:2, aqui o n é dois quando a gente preenche os dados.

# Likelihood
 for (i in 1:n) {
    "Meu modelo aqui"
 }

Dentro desse loop a gente vai por nosso modelo.
E nosso modelo é como?
Tem uma parte estocástica e uma parte determinística.

Como é a parte estocástica?

C[i] ~ dbin(p[i], N[i])

Oras, a parte estocástica é uma probabilidade, que depende do N, o número de amostras, é uma distribuição binomial. Mas esse p é determinado pela nossa função determinística, e como é ela?

logit(p[i]) <- alpha + beta *species[i]

Igual foi em todo modelo linear, a+bx, intercepto + inclinação * espécie. Vezes a espécie porque eu vou dar um valor para cada espécie. Aqui que ta o macete, eu vou atribuir a probabilidade de uma espécie como alpha, e o da outra como alpha main beta. so que beta é multiplicado por um valor, que é a hora que eu mando junto dos dados o 0 ou 1, ou seja, vamos supor um exemplo simples, eu tenho espécies A e B, A tem prevalência de 20 e B tem prevalência de 30
Nesse caso, para esse modelo, o intercepto é 20, a inclinação é 10 e o beta da especie A é 0 e o beta da especie B é 1.
Assim ficamos com o seguinte:

Para a espécie A A = interceopto + inclinação * especie A = 20 + 10 * 0 A = 20

Como qualquer número multiplicado por zero é zero, esse 10 some para a especie 1, e ela é meramente o valor do intercepto.
Ja pra a espécie B, a espécie 1, temos o seguinte:

Para a espécie B B = interceopto + inclinação * especie B = 20 + 10 * 1 B = 30

Oras, 10 vezes 1 é 10 mesmo, que somado com o 20 da 30, a prevalência da segunda espécie.
Isso é só pra facilitar a visualização do que ta dentro do likelihood ali, mas observe que o valor de p(de prevalência) esta dentro da função logit, a função logística, usado para generalizar para linearizar a relação e assim usar a distribuição binomial. Por isso o a+bx da equação vai sair com um número esquisito de olhar, que não é porcentagem, não é um número entre 0 e 1, entre 0% e 100%. Para olhar ele depois temos que destransformar, e fazemos isso usando as quantidades derivadas no Openbugs.

# Derived quantities
 prev.C_talpacoti <- exp(alpha) / (1 + exp(alpha))
 prev.S_squammata <- exp(alpha + beta) / (1 + exp(alpha + beta))
 prev.Diff <- prev.C_talpacoti - prev.S_squammata	# Teste

Agora olha essa equação ai e lembra da formula da equação logística, poxa aqui a gente so ta destransformando a continha que falamos acima, veja a formula no wikipedia, é exatamente isso a formula da função logística, so lembrar que o inverso de log é exponencial. So que além de acessar os valores de prevalência das duas espécies, olhamos também a diferença entre eles, e ja teremos uma saída direto.

Agora vamos rodar o modelo e ver o que acontece.

> out <- bugs(data=bugs.dados, inits=inits, parameters.to.save=params,model.file="Binomial.t.test.txt", + n.thin=nt, n.chains=nc, n.burnin=nb,n.iter=ni) > print(out, dig = 3) Inference for Bugs model at "Binomial.t.test.txt", Current: 3 chains, each with 1200 iterations (first 200 discarded), n.thin = 2 Cumulative: n.sims = 3000 iterations saved mean sd 2.5% 25% 50% 75% 97.5% Rhat n.eff alpha 0.065 0.253 -0.432 -0.113 0.068 0.243 0.539 1.001 3000 beta -1.636 0.424 -2.452 -1.929 -1.637 -1.342 -0.807 1.001 2400 prev.C_talpacoti 0.516 0.062 0.394 0.472 0.517 0.561 0.632 1.001 3000 prev.S_squammata 0.177 0.049 0.095 0.142 0.174 0.208 0.285 1.001 3000 prev.Diff 0.339 0.079 0.180 0.284 0.342 0.395 0.481 1.003 2800 deviance 10.488 1.855 8.609 9.127 9.914 11.250 15.331 1.001 2100 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 = 1.9 and DIC = 12.4 DIC is an estimate of expected predictive error (lower deviance is better).

Primeiro olhamos o Rhat para ver a convergência dos valores, senão nem adianta olhar o resultado, como todos estão ok a gente procede para olhar os resultados, se você não lembra o que é o Rhat, pode dar uma conferida aqui.

Mas então o alpha e beta, que é o a+bx da equação da reta ta ae, mas são números que estão transformados pela função logística. Vamos olhar direto a prevalência das espécies, que ja destransformamos os dados diretamente no bugs, e la temos, dado que simulamos dados com essas mesmas características, recuperamos muito bem os valores de prevalência, como recuperamos bem, vamos olhar a diferença, e a diferença de prevalências tem média de 0.339 e desvio de 0.079. Ta mas mais legal que isso vamos olhar qual o intervalo de confiança de 95% dele, é algo que vai de 0.180 a 0.481, ou seja dentro desse intervalo não tem nenhum valor de 0, ou negativo, isso quer dizer que no minimo existe uma diferença de 0.18 e no máximo 0.481 para as prevalência, então se pegarmos por exemplo a média de C. talpacoti, que é 0.516 de prevalência, nesse caso o S. squammata tem no minimo 0.516-0.180 e no maximo 0.516-0.481.
Essa conta é esse teste aqui:

prev.Diff <- prev.C_talpacoti - prev.S_squammata	# Teste

Mas fazemos ele para cada valor do MCMC, por isso temos a cadeiazinha dele e assim uma distribuição.
Vamos olhar esses três valores das quantidades derivadas num gráfico.

02

O primeiro e o segundo, são as distribuições de qual a prevalência, e nossa confiança no valor. O risco vermelho é a media da distribuições, bem próximo das distribuições que tinhamos originalmente certo?

So que aqui também calculamos a distribuição da diferença, e vemos que ela é quase que certamente maior que zero, a gente não vê a distribuição nem encostando no zero, e a diferença ser zero é o mesmo que não haver diferença, se você é zero cm mais alto que eu, eu posso afirmar que temos a mesma altura, mas se você é no minimo três centímetros mais alto que eu, eu não posso afirmar que temos a mesma altura, simples assim.

Então nesse caso podemos afirmar que “Existe uma diferença na prevalência de Haemoproteus nessas duas espécies”, a mais ainda, podemos afirmar que “Existe uma diferença, e ela deve ser algo entre 18% e 34% menos prevalente em S. squammata do que em C. talpacoti“, agora porque isso ocorre, precisamos ler mais do que eu li para escrever esse post.

Bem, para constar, podemos fazer o mesmo teste usando estatística frequentista assim.

> modelo.glm<-glm(cbind(dados$simul, dados$N - dados$simul) ~ dados$Sp, family = binomial) > summary(modelo.glm) Call: glm(formula = cbind(dados$simul, dados$N - dados$simul) ~ dados$Sp, family = binomial) Deviance Residuals: [1] 0 0 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 0.06454 0.25413 0.254 0.799529 dados$SpScardafella squammata -1.61210 0.43111 -3.739 0.000184 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: 1.5693e+01 on 1 degrees of freedom Residual deviance: -3.5527e-15 on 0 degrees of freedom AIC: 12.551 Number of Fisher Scoring iterations: 3

Notem que temos os mesmos valores para alpha e beta praticamente, não esperamos nada diferente, já que o prior que usamos foi pouco informativo, mas como fizemos aqui a analise de três formas diferentes, chi-quadrado, glm num framework bayesiano e glm num framework frequentistca, e tudo convergiu para a mesma resposta, existe uma diferença entre a espécies quanto a sua prevalência de Haemoproteus, devemos estar no caminho certo.

Bem acho que acabei falando um monte de coisa aqui, mas acho que agora com o tempo da para ir destrinchando o glm melhor nos próximos posts. Mas glms são muito legais, ja que da para usar, nesse mesmo esquema, muitos tipos de distribuições, cobrindo uma gama muito maior de tipos de dados que com os modelo lineares mas usando todas as estrategias que ja tinhamos para cobrir os modelos lineares, aqueles em que os dados tem uma distribuição normal.

Referencias
Adriano E.A. & Cordeiro N.S. 2001 – Prevalence and Intensity of Haemoproteus columbae in Three Species of Wild Doves from Brazil – Mem. Inst. Oswaldo Cruz 96(2) : 175-178

Marc Kéry 2010 – Introduction to WinBUGS for Ecologists.

set.seed(123)
### Gerando dados
dados<-data.frame(Sp=c("Columbina talpacoti","Scardafella squammata"),N=c(62,57),
                  Prevalence=c(0.516,0.193))
 
dados$simul<-0
 
for(i in 1:2) {
    dados[i,"simul"]<-rbinom(1, dados[i,"N"], prob = dados[i,"Prevalence"])
    }
 
layout(matrix(1:2,ncol=1,nrow=2))
for(i in 1:2) {
    plot(dbinom(x=0:dados$N[i],size=dados$N[i],prob=dados$Prevalence[i])~c(0:dados$N[i]),frame=F,pch=19,
         main=paste(dados$Sp[i]," N=",dados$N[i],sep=""),ylab="Frequência",xlim=c(0,70),xlab="",ylim=c(0,0.15))
    legend("topright",legend=paste("Prevalência =",dados$Prevalence[i]*100,"%"),bty="n")
    lines(c(dados$N[i],dados$N[i]),c(0.05,0),col="blue",lty=2,lwd=3)
    text(dados$N[i]+5,0.05,paste("Limite N =",dados$N[i]),cex=0.7)
    }
 
tabela<- as.table(cbind(dados$simul, dados$N - dados$simul))
dimnames(tabela) <- list(Sp = c("Columbina talpacoti","Scardafella squammata"),Parasita = c("Sim","Não"))
tabela
 
teste.chi <- chisq.test(tabela)
 
str(teste.chi)
teste.chi$observed   # Contagens observadas
teste.chi$expected   # Contagens esperada sobre a hipotese nula
teste.chi$residuals  # Residuo de Pearson
teste.chi$stdres     # Residos "standardized"
 
### Análise usando Bugs
# Definindo Modelo
 
sink("Binomial.t.test.txt")
cat("
model {
 
# Priors
 alpha ~ dnorm(0,0.01)
 beta ~ dnorm(0,0.01)
 
# Likelihood
 for (i in 1:n) {
    C[i] ~ dbin(p[i], N[i]) 
    logit(p[i]) <- alpha + beta *species[i]
 }
# Derived quantities
 prev.C_talpacoti <- exp(alpha) / (1 + exp(alpha))
 prev.S_squammata <- exp(alpha + beta) / (1 + exp(alpha + beta))
 prev.Diff <- prev.C_talpacoti - prev.S_squammata	# Teste
}
",fill=TRUE)
sink()
 
# Junte os dados numa lista
bugs.dados <- list(C = dados$simul, N = dados$N, species = c(0,1), n = length(dados$Sp))
 
# Função para iniciar os dados
inits <- function(){ list(alpha=rlnorm(1), beta=rlnorm(1))}
 
# Parametros para acompanhar
params <- c("alpha", "beta", "prev.C_talpacoti", "prev.S_squammata", "prev.Diff")
 
# Configurações do MCMC
nc <- 3
ni <- 1200
nb <- 200
nt <- 2
 
# Iniciando o Gibbs sampling
library(R2OpenBUGS)
out <- bugs(data=bugs.dados, inits=inits, parameters.to.save=params,model.file="Binomial.t.test.txt",
            n.thin=nt, n.chains=nc, n.burnin=nb,n.iter=ni)
 
print(out, dig = 3)
 
out$summary
 
#Figura 2
str(out)
par(mfrow = c(3,1))
hist(out$sims.list$prev.C_talpacoti, col = "grey", xlim = c(0,1), xlab="",
	main = "Prevalência de C. talpacoti", breaks = 30,ylim=c(0,300))
abline(v = out$mean$prev.C_talpacoti, lwd = 3, col = "red")
hist(out$sims.list$prev.S_squammata, col = "grey", xlim = c(0,1), xlab= "",
	main = "Prevalência de S. squammata", breaks = 30,ylim=c(0,300))
abline(v = out$mean$prev.S_squammata, lwd = 3, col = "red")
hist(out$sims.list$prev.Diff, col = "grey", xlim = c(0,1), xlab = "",
	main = "Diferença nas Prevalências", breaks = 30,ylim=c(0,300))
abline(v = 0, lwd = 3, col = "red")
 
### Analise usando GLM
modelo.glm<-glm(cbind(dados$simul, dados$N - dados$simul) ~ dados$Sp, family = binomial)
summary(modelo.glm)