Data Augmentation para dados de Captura e Recaptura

Acho que um tipo de dados muito comum em ecologia são os de Captura e Recaptura. Desde aqueles feitos com armadilhas, onde elas são vistoriadas periodicamente num grid, até as armadilhas de câmeras, que tiram fotos por sensor de movimento, redes capturando passarinhos, peixes, borboletas, talvez mais algumas formas de conseguir dados assim.

A ideia é que capturando indivíduos de alguma espécie, e recapturando depois (basta dar um jeito de identificar o um indivíduo) a gente pode calcular uma serie parâmetros populacionais que temos muito interesse, taxas de crescimento, sobrevivência, migração e assim vai.

Mas esses modelo sempre foram meio misticos pra mim, um monte de formula fechada que vem em um programa, a gente usa sem saber o que é, e usa da forma fechadinha. Mas se a gente conseguir ver esse tipo de modelo como um glm, as portas se abrem para um monte de possibilidades. Nesse post aqui a gente viu como estimar a colonização e extinção, usando o pacote unmarked e um modelo Bugs equivalente.

No modelo final, a gente vê que coisas como colonização, extinção são na verdade medidas derivadas, e que a resposta propriamente dita é encontrar ou não um individuo em uma amostra, observar ele ou não, a captura. Sendo que indivíduos podem ser recapturados. Resumindo, tentamos fazer um modelo assim.

logit(p_{i,j})=\alpha_i+\beta_j+\delta \cdot F_{i,j}+ \gamma \cdot X_{i,j}

E o \alpha_i \sim Normal(\mu_{\alpha},\sigma^2_{\alpha})

Vamos entender isso tudo, primeiro aqui, a gente ta lidando com dois índices, o i e o j, aqui o i diz respeito aos indivíduos e j ao evento de captura. Por exemplo pegando um exemplo em que a gente arma uma rede para pegar passarinhos, o i é referente a identidade do passarinho, o número da anilha dele, e o j é o dia que a rede for armada, no final isso resulta numa matriz de dados onde os indivíduos são linhas e os dias de rede armada as colunas, e dentro dessa matriz a gente marca as captura com um e zero para quando não houve captura, sendo que uma linha pode ter mais de uma captura, a gente pode recapturar um mesmo passarinho.

Agora seguindo esse exemplo, indivíduos, passarinhos podem diferir em capturabilidade, a gente pode tratar isso, tratar o fato de alguns passarinhos serem mais bobinhos e outros não assumindo que eles vem de uma distribuição normal para essa variabilidade, teremos então alguns bobinhos, alguns espertalhões mas a maioria vai ter uma determinada capturabilidade, isso é o \alpha_i da formula, veja que ele usa o índice i de individuo, e dentro desse \alpha podemos modelar outras coisas dos indivíduos, sei la, sexo por exemplo, se machos são mais curiosos, derrepente sejam mais capturados.

Agora a chance de uma captura não depende só de um individuo, por exemplo um dia ensolarado pode ter mais passarinho voando que um dia de chuva, ou dia depois da chuva, ou seja, o dia da amostra pode influenciar na capturabilidade. Isso ai a gente vei colocar dentro do \beta_j, lembrando que j é o índice dos eventos de captura. Então a gente pode modelar as coisas desse tipo aqui dentro. Lembrando que tanto no \alpha quanto o \beta, a gente ta colocando uma letrinha, mas cabe uma formulinha aqui dentro, com intercepto e inclinação de por exemplo, dados que você coleta dos passarinhos, tipo a morfometria, se você acha que ela tem a ver com a chance do passarinho ser capturado, podemos tentar observar isso.

Agora as outras duas coisas, o \delta \cdot F_{i,j} tem haver com o famoso ditado, gato escaldado tem medo de água, ou o contrario. Apos primeira comida, a chance de captura pode ser alterada, tanto para menos quanto pra mais, por exemplo, se ser capturado é muito estressante, o individuo pode começar a ficar mais esperto e difícil de capturar, ou se a gente ta lidando com uma gaiola com alguma isca, é fácil pensar que passa pela cabeça do individuo a frase “boca livre”, e ele pode começar a ser capturado com mais facilidade por procurar uma boca livre, ou sei la, veja que esse parâmetro depende tanto do individuo i como da captura j, já que esse efeito a partir da segunda captura, se alguém é capturado apenas uma vez, não tem como apresentar esse efeito, por isso leva os dois índices.

O último parâmetro que tem ali é o \gamma \cdot X_{i,j}, que ai são as questões ambientais, como é o ambiente, ou a disponibilidade de algo especial. Isso a gente enquadra ai nesse \gamma, veja que ele também depende de i e j.

Certo então com o modelo ali em cima a gente conseguiria estimar uma probabilidade de captura, mas qual o interesse nisso?
A gente se interessa em quantidades mais palpáveis, o tamanho da população, sobrevivência, imigração, mas essas quantidade pode ser derivada dessa chance de captura, sendo que com o modelo acima a gente tem a chance de considerar praticamente tudo que é do nosso interesse.

Agora vamos simular uns dados aqui para ver melhor como isso acontece. Vamos supor uma população de 100 indivíduos, ai vamos supor também uma chance de captura de 50% e que coletamos por três vezes, em três eventos de coleta.

Beleza, se tudo que colocamos no modelo la em cima for constante, teríamos o seguinte.

> N = 100 > p = 0.5 > T = 3 > Real<-array(NA, dim = c(N, T)) > for (j in 1:T){ + Real[,j] <- rbinom(n = N, size = 1, prob = p) + } > Real [,1] [,2] [,3] [1,] 0 1 1 [2,] 1 0 1 [3,] 0 1 0 [4,] 0 0 1 . . . [95,] 0 1 1 [96,] 1 0 1 [97,] 0 0 0 [98,] 1 0 1 [99,] 1 1 1 [100,] 0 1 1

Certo, colocamos nossos parâmetros e meramente sorteamos de uma distribuição binomial a com a chance de captura que definimos se capturamos ou não. Temos uma matriz de 100 linhas, que são 100 indivíduos e 3 colunas, que são os três eventos de coleta, agora da uma olhada o que aconteceu na linha 97, a gente nunca capturou esse individuo nos três eventos.

Se a gente somar as linhas dessa matriz, a gente vai ver que

> detectado <- apply(Real, 1, max) > sum(detectado) [1] 86

vimos somente 86 indivíduos dos 100 que existe, existem 14 sortudos que nunca caíram na rede, simples assim.
Agora a pergunta é a seguinte, existem 100 indivíduos, mas a gente não ve uma matriz como a acima, esse valor de 100 esta oculto pela natureza, aqui a gente sabe porque inventamos ele, mas com dados reais a gente não sabe.
Os dados reais que a gente ve vão ser assim:

> dados<-array(NA, dim = c(N, T)) > dados <- Real[detectado == 1,] > dados [,1] [,2] [,3] [1,] 0 1 1 [2,] 1 0 1 [3,] 0 1 0 . . . [83,] 1 0 1 [84,] 1 0 1 [85,] 1 1 1 [86,] 0 1 1

So vão até a linha 86, porque as 14 linhas que so tem zero a gente não viu, não ve e nunca verá. Nossas observações são sempre imperfeitas. Agora como faz para chegar aquele valor de 100 indivíduos, se a gente so tem nossas observações acima?

Existiam outros estimadores, presentas em programas como o Capture e o Mark, eles estimavam tudo e depois faziam uma conta para estimar o tamanho da população, aquele valor de 100. Agora quando a gente ta usando estatística bayesiana, com MCMC e tal, existe um problema que a gente vai de iteração a iteração, e o vetor de N pode mudar de iteração a iteração, ou seja não cabe nos dados que a gente manda para o BUGS. Mas como problemas existem para serem solucionados, dois caras chamados Tanner e Wong inventaram uma técnica chamada de data augmentation, que o Royle e companhia trouxeram pros modelos de captura e recaptura usando MCMC.

Basicamente duas coisas são feitas, adicionamos um numero arbitrário de zeros a matriz de dados analisamos uma versão reparametrizada dos dados originais. Essencialmente convertendo um modelo de população fechada para um modelo de ocupação.

Então a gente está adicionando um grande numero de observações em potencial que não foram observadas. Os dados aumentados terão um número de linhas M e o mesmo numero de colunas T das outras observações, se o tamanho da população deveria ser o número de linhas, que deveria ser 100, e só temos 86 linhas, ao aumentarmos os dados teremos uma quantidade linhas agora muito muito maior que o tamanho da população, que é o que nos interessa. Nesses dados a gente ajusta um tipo de modelo de zero inflado se o tamanho da população fosse conhecido. Para isso a gente adiciona um indicador binario que indica se uma linha da matriz de dados aumentada representa um individuo real ou um que não existe na pratica. Esse indicador é dado por um prior da distribuição de Bernoulli, ou Binomial, e esse parâmetro, que podemos chamar de \omega (omega) é estimável dos dados. Uma probabilidade de inclusão. Com essa probabilidade de inclusão a gente tem o tamanho da população saindo dos dados.

Então criamos nosso modelo na linguagem bugs

# Modelo na linguagem Bugs
sink("modelo.txt")
cat("
    model {
        # Priors
        omega ~ dunif(0, 1)
        p ~ dunif(0, 1)
        # Likelihood
        for (i in 1:M){
            z[i] ~ dbern(omega)
            # Indicador de inclusão
            for (j in 1:T){
                dadosaug[i,j] ~ dbern(p.eff[i,j])
                p.eff[i,j] <- z[i] * p # So pode ser detectado se z=1
                } #j
            } #i
        # Quantidades derivadas
        N <- sum(z[])
        }
    ",fill = TRUE)
sink()

Vamos observar alguma coisas, primeiro que temos dois priores

omega ~ dunif(0, 1)
p ~ dunif(0, 1)

Que são as chances daquela linha ser um indivíduo e a chance do indivíduo ter sido capturado.

Veja então que no modelo de verossimilhança, a gente tem uma chance do individuo ser existir

for (i in 1:M){
    z[i] ~ dbern(omega)

E a chance dele ter sido capturado num dado evento pra esse individuo

for (j in 1:T){
    dadosaug[i,j] ~ dbern(p.eff[i,j])
    p.eff[i,j] <- z[i] * p # So pode ser detectado se z=1

Mas veza que ele é multiplicado pelo valor de z, que vai ser zero ou um, com uma chance omega, que é a chance dele ser um individuo real.

N <- sum(z[])

Ou seja a soma do vetor de z vai ser exatamente a estimativa do tamanho populacional.

Agora é só juntar os dados, gerar valores iniciais, configurar as cadeias e rodar o esse modelo. E temos o seguinte resultado.

> print(out, dig = 3) Inference for Bugs model at "modelo.txt", Current: 3 chains, each with 2500 iterations (first 500 discarded), n.thin = 2 Cumulative: n.sims = 6000 iterations saved mean sd 2.5% 25% 50% 75% 97.5% Rhat n.eff N 95.451 4.114 89.000 92.000 95.000 98.000 105.000 1.001 2700 p 0.542 0.037 0.466 0.517 0.542 0.567 0.613 1.001 4100 omega 0.405 0.036 0.337 0.380 0.404 0.429 0.478 1.001 4000 deviance 395.385 18.993 363.200 383.300 393.500 407.200 436.600 1.001 2800 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 = 180.3 and DIC = 575.7 DIC is an estimate of expected predictive error (lower deviance is better).

Sempre temos que fazer aquelas validações, que tivemos uma convergência dos parâmetros, mas podemos ver que funciona bem.

Primeiro vamos olhar o N que é o tamanho populacional

mean sd 2.5% 25% 50% 75% 97.5% Rhat n.eff N 95.451 4.114 89.000 92.000 95.000 98.000 105.000 1.001 2700

Veja que temos uma média de 95 para o parâmetro, que é perto de 100, o valor real do tamanho populacional, com um desvio padrão de 4.
Olhando o intervalo de confiança vemos que ele vai de 89 a 105, o que inclui o valor real de 100, ou seja, deu certo, esse negocio funciona mesmo para estimar o tamanho populacional.

E estimamos bem também a chance de captura, que na realidade era de 0.5 (50%), que esta dentro do intervalo de confiança.

mean sd 2.5% 25% 50% 75% 97.5% Rhat n.eff p 0.542 0.037 0.466 0.517 0.542 0.567 0.613 1.001 4100

Então vamos tentar fazer um gráfico para ver melhor essa saída.

01

Certo, nossa estimativa cobre o valor real, dentro do intervalo de confiança de 95%, não sei nem se posso falar intervalo de confiança, mas é o 95% ali das saídas, as vezes esses termos tem pequenos pressupostos que a gente não sabe, até escrever errado.
Mas agora podemos ver como o \omega (omega) faz, vamos olhar o valor dele.

mean sd 2.5% 25% 50% 75% 97.5% Rhat n.eff omega 0.405 0.036 0.337 0.380 0.404 0.429 0.478 1.001 4000

Média de 0.405 com desvio de 0.036, esse 40.5% é o número de linhas da matriz de dados aumentada que mandamos que são provavelmente realmente indivíduos. Como a gente tinha 86 observações, e aumentamos a matriz em 150 linhas, temos um total de 236 linhas. E dessas 40.5% consideramos como indivíduos.

> 0.405*(86+150) [1] 95.58

As outras não devem ser observações, apenas linhas. Mas o aumento deu a possibilidade de considerar esses individuos que não observamos, mas deviam estar la.
Mas uma pergunta razoavel é o quanto devemos aumentar a matriz. Porque a gente aumentou 150 aqui e deu certo, mas será que qualquer número funciona?
Bem sabemos que o tamanho da população vai ser igual ao número das observações, ou maior. Se aumentamos pouco a matriz, por exemplo, somente em 5 linhas, podemos ver o que acontece.

> print(out5, dig = 3) Inference for Bugs model at "modelo.txt", Current: 3 chains, each with 2500 iterations (first 500 discarded), n.thin = 2 Cumulative: n.sims = 6000 iterations saved mean sd 2.5% 25% 50% 75% 97.5% Rhat n.eff N 90.004 1.093 87.000 89.000 90.000 91.000 91.000 1.001 6000 p 0.574 0.031 0.512 0.553 0.574 0.595 0.636 1.001 5300 omega 0.978 0.019 0.930 0.970 0.984 0.992 0.999 1.002 6000 deviance 369.337 5.821 354.500 365.200 369.700 373.700 376.900 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 = var(deviance)/2) pD = 16.9 and DIC = 386.3 DIC is an estimate of expected predictive error (lower deviance is better).

Veja que o valor de omega ficou espremido, grudado no 100%, olhando o grafico isso fica mais facil de visualizar.

02

Veja que a distribuição é bem truncada, como se tivesse uma barreira ali segurando a distribuição de ir mais para direita, ou seja a gente precisava de mais “espaço”. Então se a gente ver uma distribuição truncada assim, sabemos que temos que aumentar o tamanho da matriz mais um pouco. Com 150 deu legal, podemos olhar denovo isso.

> print(out150, dig = 3) Inference for Bugs model at "modelo.txt", Current: 3 chains, each with 2500 iterations (first 500 discarded), n.thin = 2 Cumulative: n.sims = 6000 iterations saved mean sd 2.5% 25% 50% 75% 97.5% Rhat n.eff N 95.551 4.222 89.000 93.000 95.000 98.000 105.000 1.001 6000 p 0.542 0.038 0.467 0.517 0.543 0.568 0.615 1.002 1700 omega 0.406 0.037 0.339 0.381 0.405 0.430 0.481 1.001 5800 deviance 395.845 19.407 363.300 383.300 393.500 407.400 438.710 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 = var(deviance)/2) pD = 188.3 and DIC = 584.2 DIC is an estimate of expected predictive error (lower deviance is better).

03

Aqui a distribuição ficou legal, veja que o omega não ficou espremido, truncado nem para a esquerda nem para a direita. Então está tudo ok. Mas e se a gente aumentar demais, tem problema? Será que o precavido aqui vai se ferrar?

> print(out1500, dig = 3) Inference for Bugs model at "modelo.txt", Current: 3 chains, each with 2500 iterations (first 500 discarded), n.thin = 2 Cumulative: n.sims = 6000 iterations saved mean sd 2.5% 25% 50% 75% 97.5% Rhat n.eff N 95.501 4.175 89.000 92.000 95.000 98.000 105.000 1.001 6000 p 0.542 0.037 0.468 0.517 0.542 0.568 0.614 1.001 2800 omega 0.061 0.006 0.049 0.056 0.061 0.065 0.074 1.001 3400 deviance 395.610 19.188 363.300 383.300 393.600 407.100 437.600 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 = var(deviance)/2) pD = 184.1 and DIC = 579.7 DIC is an estimate of expected predictive error (lower deviance is better).

04

E a resposta é não, não tem problema nenhum aumentar muito. O que acontece é que o omega vai ficar baixinho, mas não ficando truncado ta tudo certo. Veja que aumentando em 1500 linhas a matriz, temos exatamente a mesma estimativa de tamanho populacional, unica coisa é que aqui tivemos uma matriz gigante para processar, e aqui o 6.1% de 1586 linhas que é o total de indivíduos.

Eu achei essa solução muito mirabolante. Primeiro porque a gente cai em um glm, o que torna um modelo de captura e recaptura algo bem próximo dos glm que comumente olhamos e essa sacada de aumentar a matriz é muito foda. Melhor sobrar que faltar. Agora aqui conseguimos ter uma noção de que data augumentation funciona, mas em um modelo simples.

No inicio do post colocamos o seguinte modelo.
logit(p_{i,j})=\alpha_i+\beta_j+\delta \cdot F_{i,j}+ \gamma \cdot X_{i,j}
E o \alpha_i \sim Normal(\mu_{\alpha},\sigma^2_{\alpha})

Mas o que ajustamos foi algo como

logit(p_{i,j})=\alpha

Somente um parâmetro, so lembrar do modelo, que so tem uma linha aproximando nossas observações.

dadosaug[i,j] ~ dbern(p.eff[i,j])

Agora temos que começar a expandir esse parametro para o modelo completo. Sendo que não somos obrigados a tentar ajustar tudo, podemos usar somente aquilo que nos interessa, ou que temos dados coletados.

Bem é isso ai, hora de ir dormir. 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 boa noite.

Referencia:
Kéry e Shaub 2011 – Bayesian population analysis using WinBUGS. A hierarchical perspective – Academic Press 554pp

#Gerando dados
N = 100
p = 0.5
T = 3
Real<-array(NA, dim = c(N, T))
for (j in 1:T){
    Real[,j] <- rbinom(n = N, size = 1, prob = p)
    }
#A visão da natureza
Real
 
detectado <- apply(Real, 1, max)
sum(detectado)
 
dados<-array(NA, dim = c(N, T))
dados <- Real[detectado == 1,]
#A nossa visão
dados
 
# Aumentando os dados com 150 individuos em potencial
nz <- 150
dadosaug <- rbind(dados, array(0, dim = c(nz,T)))
 
# Modelo na linguagem Bugs
sink("modelo.txt")
cat("
    model {
        # Priors
        omega ~ dunif(0, 1)
        p ~ dunif(0, 1)
        # Likelihood
        for (i in 1:M){
            z[i] ~ dbern(omega)
            # Indicador de inclusão
            for (j in 1:T){
                dadosaug[i,j] ~ dbern(p.eff[i,j])
                p.eff[i,j] <- z[i] * p # So pode ser detectado se z=1
                } #j
            } #i
        # Quantidades derivadas
        N <- sum(z[])
        }
    ",fill = TRUE)
sink()
 
# Juntando os dados para mandar para o bugs
bugs.data <- list(dadosaug = dadosaug, M = nrow(dadosaug), T = ncol(dadosaug))
# Função geradora de valores iniciais
inits <- function() list(z = rep(1, nrow(dadosaug)), p = runif(1, 0, 1))
# Parametros a serem monitorados
params <- c("N", "p", "omega")
 
# Configurando o MCMC
ni <- 2500
nt <- 2
nb <- 500
nc <- 3
# Chamando o openbugs do R
library(R2OpenBUGS)
out <- bugs(data=bugs.data,inits=inits,parameters=params,model="modelo.txt",
            n.chains = nc,n.thin = nt, n.iter = ni, n.burnin = nb)
 
 
# Distribuições posteriores
print(out, dig = 3)
 
str(out)
out$summary
 
#Figura 1
hist(out$sims.list$N,nclass=50,col="gray",main="",xlab ="Tamanho Populacional",las = 1, xlim = c(80, 150),
     ylab="Frequência")
abline(v = sum(detectado),lwd = 3)
abline(v = 100,lwd=3,col="red")
abline(v = out$summary[1,c(3,7)],lwd = 1,col="blue",lty=2)
legend("topright",legend=c("Detectados","Tamanho real","Intervalo de 95% estimado"),lwd=c(3,3,1),lty=c(1,1,2),
       col=c("black","red","blue"))
 
 
nz <- 5
dadosaug <- rbind(dados, array(0, dim = c(nz,T)))
bugs.data <- list(dadosaug = dadosaug, M = nrow(dadosaug), T = ncol(dadosaug))
 
out5 <- bugs(data=bugs.data,inits=inits,parameters=params,model="modelo.txt",
            n.chains = nc,n.thin = nt, n.iter = ni, n.burnin = nb)
print(out5, dig = 3)
 
#Figura 2
hist(out5$sims.list$N,nclass=50,col="gray",main="Dados aumentados em 5",xlab ="Tamanho Populacional",las = 1, xlim = c(80, 150),
     ylab="Frequência")
abline(v = sum(detectado),lwd = 3)
abline(v = 100,lwd=3,col="red")
abline(v = out5$summary[1,c(3,7)],lwd = 1,col="blue",lty=2)
legend("topright",legend=c("Detectados","Tamanho real","Intervalo de 95% estimado"),lwd=c(3,3,1),lty=c(1,1,2),
       col=c("black","red","blue"))
 
 
nz <- 150
dadosaug <- rbind(dados, array(0, dim = c(nz,T)))
bugs.data <- list(dadosaug = dadosaug, M = nrow(dadosaug), T = ncol(dadosaug))
 
 
out150 <-bugs(data=bugs.data,inits=inits,parameters=params,model="modelo.txt",
            n.chains = nc,n.thin = nt, n.iter = ni, n.burnin = nb)
print(out150, dig = 3)
 
#Figura 3
hist(out150$sims.list$N,nclass=50,col="gray",main="Dados aumentados em 150",xlab ="Tamanho Populacional",las = 1, xlim = c(80, 150),
     ylab="Frequência")
abline(v = sum(detectado),lwd = 3)
abline(v = 100,lwd=3,col="red")
abline(v = out150$summary[1,c(3,7)],lwd = 1,col="blue",lty=2)
legend("topright",legend=c("Detectados","Tamanho real","Intervalo de 95% estimado"),lwd=c(3,3,1),lty=c(1,1,2),
       col=c("black","red","blue"))
 
 
nz <- 1500
dadosaug <- rbind(dados, array(0, dim = c(nz,T)))
bugs.data <- list(dadosaug = dadosaug, M = nrow(dadosaug), T = ncol(dadosaug))
 
out1500 <- bugs(data=bugs.data,inits=inits,parameters=params,model="modelo.txt",
            n.chains = nc,n.thin = nt, n.iter = ni, n.burnin = nb)
print(out1500, dig = 3)
 
#Figura 4
hist(out1500$sims.list$N,nclass=50,col="gray",main="Dados aumentados em 1500",xlab ="Tamanho Populacional",las = 1, xlim = c(80, 150),
     ylab="Frequência")
abline(v = sum(detectado),lwd = 3)
abline(v = 100,lwd=3,col="red")
abline(v = out1500$summary[1,c(3,7)],lwd = 1,col="blue",lty=2)
legend("topright",legend=c("Detectados","Tamanho real","Intervalo de 95% estimado"),lwd=c(3,3,1),lty=c(1,1,2),
       col=c("black","red","blue"))

Usando Offset em modelos lineares

A gente falou aqui sobre modelos com Overdispersion e aqui sobre modelos com Zeros-Inflados.

Usando a linguagem Bugs, por um lado a gente pode sofre um pouco por ser obrigado a definir o modelo, mas por outro lado, a gente ganha no melhor entendimento do que está acontecendo. Overdispersion é um exemplo, que somente depois que eu vi o modelo no Bugs que eu fui entender o que diabos era aquele número que saio no glm quando a gente usava a family=quasipoisson.

Agora ta na hora de entender o que significa offset, aquela opção que você pode usar nos modelos lineares da função lm, nos modelos gerais lineares na função glm, afinal, em quase todo lugar.

Se a gente olhar a documentação da função com ?lm, a gente vai ver a descrição o argumento offset:

“this can be used to specify an a priori known component to be included in the linear predictor during fitting. This should be NULL or a numeric vector of length equal to the number of cases. One or more offset terms can be included in the formula instead or as well, and if more than one are specified their sum is used. See model.offset.”

No meu inglês nórdico, podemos ler:

“Este pode ser usado para especificar um componente conhecido a priori para ser incluído como preditor linear no ajuste. Este deve ser NULL ou um vetor numérico de tamanho igual ao número de casos. Um ou mais termos de offset podem ser incluídos na formula ao invés ou tal como, e se mais de um for especificado, sua soma será usada, olhe model.offset.”

Se você olhar a função model.offset, vai cair nas funções que transformam formula em matrizes, álgebra linear e não da mais para entender nada.

Mas vamos pensar de maneira mais simples. No Glm de Poisson, como estamos olhando, nos assumimos que as contagens esperadas são adequadamente descritas pelos efeitos das covariáveis. Porém, frequentemente, nos temos algo como uma “janela de contagem” que não é constante. Por exemplo, áreas de estudos que não tem o mesmo tamanho ou amostras que tem duração diferentes, a gente pode não determinar antes o tempo que vai olhar uma lagoa, simplesmente vamos olhando e contando quantos sapinhos vemos la. Para adicionar essa variação no modelo de Poisson, por exemplo se temos contagens de coelhos em vários campos de tamanho diferentes, nos podemos definir a área como offset, mais especificamente, como o modelo de Poisson a gente usa função de ligação logarítmica, a gente define o log da área de cada amostras. Simplificando, a gente modela a densidade como a resposta, e não a contagem per se.

Vamos fazer nosso exemplo aqui, vamos pensar em coelhinhos. Estamos contando coelhinhos em pastos e plantações, vamos fazer varias amostras, mais ou menos assim:

01

Cada pontinho é uma amostra. Mas cada uma dessas amostras são de areas cercados, com curva de nível, varias coisas no campo que dividem elas, então não são contínuos iguais, temos tamanhos diferentes para cada amostras. Podemos representar esses tamanhos aqui com o tamanho da bolinhas sendo equivalente ao tamanho da amostra.

02

No modelo de Poisson a gente tem:

C_i \sim Poisson(\lambda_i)

Cada amostra vem de uma distribuição de Poisson com média \lambda.
Agora como cada amostra foi feita num pequeno pasto de área A, nosso modelo então é na verdade:

C_i \sim Poisson(A_i \cdot \lambda_i)

Dessa forma, o preditor linear fica assim:

log(A_i \cdot \lambda_i)

E pelas regras de logaritmo, a gente pode separar essa multiplicação dentro do log em dois logs.

log(A_i) + log(\lambda_i)

Agora esse \lambda pode vir de um preditor linear, por exemplo com a variável x sendo o tipo de pasto que coletamos. Então

log(A_i \cdot \lambda_i) = log(A_i) + \alpha + \beta \cdot x_i

Veja que ficou assim porque podemos fazer aquela separação da área e do \lambda, por causa da multiplicação no log.

De forma mais simples, estamos pensando nesse log da área com um coeficiente 1 multiplicando ele, ou seja, a área fazendo o efeito dela proporcionalmente, mas isso pode ser mudado, usando um coeficiente para a área.

log(A_i \cdot \lambda_i) = \beta_0 \cdot log(A_i) + \alpha + \beta \cdot x_i

Agora é so passar isso para a linguagem Bugs da seguinte forma.

 for (i in 1:n) {
    C[i] ~ dpois(lambda[i])
    log(lambda[i]) <- 1 * logA[i] + alpha + beta *x[i]
 }

São duas linhas, a linha que tem a ~ é a parte estocástica, e a linha com <- é a parte determinística, que está sendo somado, além do preditor linear, o offset da área, ja logaritimizado, a gente não ta mandando pro bugs a área real, a gente já calcula o log no R e só mando o valor do log. Então é juntar os dados, configurar o MCMC, fazer a função de inicialização das cadeias e rodar o modelo. E para frisar, veja que a gente ta mandando o termo logA que é o log(A).

# Juntando os dados numa lista para o OpenBugs
bugs.data <- list(C = C, x = as.numeric(x)-1, logA = log(A), n = length(x))

Rodamos o modelo, mesmo com cadeias pequenas ele converge rápido, é um modelo relativamente simples e vemos o seguinte.

> print(out, dig = 3) Inference for Bugs model at "Offset.txt", Current: 3 chains, each with 1100 iterations (first 100 discarded), n.thin = 2 Cumulative: n.sims = 3000 iterations saved mean sd 2.5% 25% 50% 75% 97.5% Rhat n.eff alpha 0.579 0.126 0.323 0.493 0.583 0.669 0.809 1.003 3000 beta 0.966 0.148 0.685 0.861 0.966 1.064 1.255 1.001 3000 deviance 97.335 1.974 95.430 95.960 96.745 98.080 102.902 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 = 2.0 and DIC = 99.3 DIC is an estimate of expected predictive error (lower deviance is better).

Todos os parâmetros com um Rhat perto de 1, e vemos a diferença entre os dois tipos de locais que observamos.
Lembrando que beta aqui é a diferença entre os dois locais, isso pela forma como escrevemos o preditor linear, de uma olhada aqui se não entendeu essa parte. Mas vemos que essa diferença tem média 0.966, mas mais importante que isso o intervalo de confiança de 95% dessa distribuição vai de 0.685 a 1.255, ou seja , não engloba 0 (nenhuma) diferença, alias, com 0.148 de desvio padrão para esse parâmetro, o valor de zero diferença está muitos desvios padrão de distancia da média, temos muita certeza da diferença.

Agora a titulo de curiosidade vamos explorar aqui o que acontece ao ignorarmos o offset, ou seja construir o modelo de forma errada.

> glm.fit.sem.offset <- glm(C ~ x, family = poisson) > summary(glm.fit.sem.offset) Call: glm(formula = C ~ x, family = poisson) Deviance Residuals: Min 1Q Median 3Q Max -1.4298 -0.7816 -0.1419 0.8385 1.7360 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 1.8245 0.1270 14.367 < 2e-16 *** xPlantação 0.9355 0.1499 6.242 4.31e-10 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for poisson family taken to be 1) Null deviance: 60.948 on 19 degrees of freedom Residual deviance: 17.615 on 18 degrees of freedom AIC: 103.45 Number of Fisher Scoring iterations: 4

Veja que no modelo bayesiano com offset, o alpha ou intercepto, como quiser chamar, ficou 0.579, sem o offset aqui ele ficou como 1.8245, ou seja muito maior, super-estimamos muito a contagem de coelhos, esse valor vai dar quantos coelhos esperamos achar no pasto e nas plantações, lembrando que nas plantações o lambda vai ser alpha mais beta, ou aqui intercept + xPlantação, então por mais que a diferença até que fica parecida, temos os valores de lambda bem super-estimado.

Mas não confunda com o problema ser modelo bayesiano ou não, se ajustarmos o modelo frequentista temos praticamente os mesmo parametros do modelo bayesiana, com o mesmo erro.

> glm.fit.com.offset <- glm(C ~ x, family = poisson, offset = log(A)) > summary(glm.fit.com.offset) Call: glm(formula = C ~ x, family = poisson, offset = log(A)) Deviance Residuals: Min 1Q Median 3Q Max -1.68068 -0.44061 -0.06695 0.31836 1.69544 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 0.5811 0.1270 4.575 4.75e-06 *** xPlantação 0.9658 0.1499 6.445 1.16e-10 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for poisson family taken to be 1) Null deviance: 59.836 on 19 degrees of freedom Residual deviance: 13.535 on 18 degrees of freedom AIC: 99.376 Number of Fisher Scoring iterations: 4

Apesar da grande diferença nas estimativas dos parâmetros, o uso ou não de offset, não é facil detectar que ele esta faltando, temos que ter pensado isso a priori. Se a gente olhar os resíduos de ambos os modelos não vemos nenhuma tendencia, pelo menos eu não vejo.

03

Temos um desvio residual menor no modelo empregando o offset, 13.535 com offset contra 17.615 sem offset, mas a gente nunca saberia disso se não fizer o ajuste com o offset, a novamente, não pegaria a necessidade disso olhando somente um plot dos resíduos sem offset.

Bem, é isso, mais um conceito legal para se entender, e facilitar a vida ao ler artigos e quem saber resolver os problemas de alguém. O repositório recologia esta la. E é isso, até o próximo post.

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

### Gerando dados
set.seed(1234)
n.locais <- 10
A <- runif(n = 2*n.locais, 2,5)		# Area de 2 a 5 km2 
x <- gl(n = 2, k = n.locais, labels = c("Pasto", "Plantação"))
preditor.linear <- log(A) + 0.69 +(0.92*(as.numeric(x)-1))
lambda <- exp(preditor.linear)
C <- rpois(n = 2*n.locais, lambda = lambda)
 
#Figura 1
stripchart(C~x,vertical=T,pch=19,method="stack",offset=0.8,at=c(1.2,2.2),frame=F,xlab="Contagem",ylab="Local")
 
#Figura 2
stripchart(C~x,vertical=T,pch=19,method="jitter",offset=0.8,at=c(1,2),frame=F,xlab="Contagem",ylab="Local",cex=A/2)
 
### Análise com OpenBugs
# Definindo o modelo
sink("Offset.txt")
cat("
model {
# Priors
 alpha ~ dnorm(0,0.001)
 beta ~ dnorm(0,0.001)
 
# Likelihood
 for (i in 1:n) {
    C[i] ~ dpois(lambda[i])
    log(lambda[i]) <- 1 * logA[i] + alpha + beta *x[i]
 }
}
",fill=TRUE)
sink()
 
# Juntando os dados numa lista para o OpenBugs
bugs.data <- list(C = C, x = as.numeric(x)-1, logA = log(A), n = length(x))
 
# Função para iniciar as cadeias com valores ao acaso
inits <- function(){ list(alpha=rlnorm(1), beta=rlnorm(1))}
 
# Parametros a estimar
params <- c("alpha","beta")
 
# MCMC settings
nc <- 3    # Número de cadeias
ni <- 1100 # Tamanho de cada cadeia
nb <- 100  # Número de amostras a descartar como burn-in
nt <- 2    # Thinning rate
 
# Inicie o Gibbs-sampler
library(R2OpenBUGS)
 
out <- bugs(data=bugs.data, inits=inits, parameters.to.save=params, model.file="Offset.txt",
            n.thin=nt,n.chains=nc,n.burnin=nb, n.iter=ni)
 
#Resultado
print(out, dig = 3)
 
### Analise usando R
 
glm.fit.sem.offset <- glm(C ~ x, family = poisson)
glm.fit.com.offset <- glm(C ~ x, family = poisson, offset = log(A))
summary(glm.fit.sem.offset)
summary(glm.fit.com.offset)
 
#figura 3
par(mfrow=c(2,1))
plot(resid(glm.fit.sem.offset),ylim=c(-2.5,2),main="Sem Offset",pch=19,ylab="Resíduos",xlab="Amostras",cex=1.3)
plot(resid(glm.fit.com.offset),ylim=c(-2.5,2),main="Com Offset",pch=19,ylab="Resíduos",xlab="Amostras",cex=1.3)

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)