Um pequeno exemplo de um modelo com Zeros Inflados

Modelos inflados por zeros podem ser muitos interessantes. Vou tentar fazer uma breve introdução com um exemplo simples. 

Vamos supor a seguinte situação, eu gostaria de saber se o número de aranhas presentes numa especie de arbusto A é diferente do número de aranhas comumente encontrado num arbusto B. Utilizando o número magico de coletas que me foi passado na faculdade, eu procurei trinta arbustos da especie A e trinta arbustos da especie B e contei quantas aranhas tinham neles:

#Simulando os dados

set.seed(13)
população<-factor(rep(c("A","B"),each=30))
contagem<-rpois(60,rep(c(2,4),each=30))

#

 Vemos que no Arbusto B, em média temos 4 aranhas, enquanto no arbusto A temos 2 aranhas.

Podemos testar essa afirmação através de um modelo de regressão de poisson.

#Modelo de Poisson

modelo1<-glm(contagem~população,family="poisson")
summary(modelo1)
exp(coef(modelo1))

#

E obtemos como resultado:

Call:
glm(formula = contagem ~ população, family = "poisson")

Deviance Residuals:
    Min       1Q   Median       3Q      Max 
-1.9833  -0.8255   0.0237   0.5868   2.9973 

Coefficients:
            Estimate Std. Error z value Pr(>|z|)   
(Intercept)   0.6763     0.1302   5.195 2.05e-07 ***
populaçãoB    0.6587     0.1604   4.107 4.01e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 90.678  on 59  degrees of freedom
Residual deviance: 72.885  on 58  degrees of freedom
AIC: 235.15

Number of Fisher Scoring iterations: 5

> exp(coef(modelo1))
(Intercept)  populaçãoB
   1.966667    1.932203 

Vejam só, o modelo nos da varias estrelinhas, nossa afirmação foi na mosca, mas mais importante que isso, temos estimativas do tamanho da população, e em média elas ficam bem próxima de realidade, o intercepto que representa a população nos arbustos A tem 1.966.. aranhas, bem próximo do 2 que usamos pra simular os dados, e no arbusto B temos 1.966 + 1.932.., que da praticamente 4, o numero que usamos na simulação. Muito legal, muitos p significativos e condecorados com variasssss estrelinhas.

Mas nem sempre o mundo é tão simples assim. Como nosso amigo Levins propôs uma vez, nem sempre todos os arbustos vão ter aranhas, eles podem não estar ocupados por aranhas. 

Pra gente contar aranhas em um arbusto, tem que ter mais de uma aranha la , ou a gente vai contar zero aranhas, até ai tudo bem pois alguns arbustos podem ser um péssimo lugar pra se viver e nenhuma aranha fica la, ou seja zero aranhas.

Mas podem existir um outro tipo de zero, o zero de aranhas devido a elas nunca terem chegado la, o arbusto pode ser perfeito pra abrigar aranhas, mas nenhuma chegou la. Agora vamos colocar isso nos dados e ver o que acontece com nosso modelo.

#Simulando ocupação

ocupação<-rbinom(60,size=1,prob=0.60)
contagem.ocup<-contagem*ocupação

#Modelo 2 de poisson

 modelo2<-glm(contagem.ocup~população,family="poisson")
summary(modelo2)
exp(coef(modelo2))

#

Call:
glm(formula = contagem.ocup ~ população, family = "poisson")

Deviance Residuals:
    Min       1Q   Median       3Q      Max 
-2.0656  -1.8719  -0.3133   0.6582   4.2838 

Coefficients:
            Estimate Std. Error z value Pr(>|z|)   
(Intercept)   0.4906     0.1429   3.434 0.000594 ***
populaçãoB    0.2671     0.1898   1.407 0.159458   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 160.14  on 59  degrees of freedom
Residual deviance: 158.14  on 58  degrees of freedom
AIC: 264.01

Number of Fisher Scoring iterations: 6

> exp(coef(modelo2))
(Intercept)  populaçãoB
   1.633333    1.306122

É muito triste, mas o fato de nem todos os arbustos estarem ocupados por aranhas nos tira muitas muitas estrelinhas, mas pior que isso, agora concluímos que os arbustos A e B tem o mesmo numero de aranhas, isso é mal, muito mal. Olha a estimativa do intercepto 1.63, usamos 2 na simulação, mas os zeros adicionais puxaram a media do numero de aranhas nos arbustos pra baixo.

Mas o que acontece é que nesse caso não estamos levando em conta que temos 2 tipos de zeros nas coletas, existem zeros de o arbusto ter zero aranhas, mas existem zeros devido ao arbusto não estar ocupado, aranhas nunca chegaram la.

O que deveríamos fazer é tentar fazer um modelo onde esses dois tipos de situação são levados em conta.

Temos de declarar ao modelo que existem zeros devido a uma chance de ocupação que cada arbusto tem, ai considerando os arbustos ocupados que deveríamos fazer as contas para o numero de aranhas no arbusto. Essa é a solução de modelos com zero inflados tenta trazer.

#modelo zero inflado usando o pacote pscl

library(pscl)
modelo3<-zeroinfl(contagem.ocup~população|1)
summary(modelo3)
exp(coef(modelo3)[1:2])
plogis(coef(modelo3)[3])

#

Call:
zeroinfl(formula = contagem.ocup ~ população | 1)

Pearson residuals:
    Min      1Q  Median      3Q     Max
-0.9876 -0.8934 -0.2306  0.4150  3.2704

Count model coefficients (poisson with log link):
            Estimate Std. Error z value Pr(>|z|)   
(Intercept)   0.8296     0.1540   5.388 7.13e-08 ***
populaçãoB    0.5977     0.2019   2.960  0.00307 **

Zero-inflation model coefficients (binomial with logit link):
            Estimate Std. Error z value Pr(>|z|)
(Intercept)  -0.4567     0.2829  -1.615    0.106
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Number of iterations in BFGS optimization: 8
Log-likelihood: -104.9 on 3 Df

> exp(coef(modelo3)[1:2])
count_(Intercept)  count_populaçãoB
         2.292305          1.817889

> plogis(coef(modelo3)[3])
zero_(Intercept)
       0.3877701

A primeira boa noticia é que nossas estrelinhas são devolvidas, temos muitos p significativos. Mas vamos olhar as nossas estimativas, agora estimamos para o arbusto A uma média de 2.29 e para o B 2.29 + 1.81, o que é bem próximo de 4, 2 e 4 foram exatamente os valores que usamos.

Olhamos ainda a estimativa de ocupação, é de 38%, tudo bem que usamos 60%, mas é bastante razoável, e temos que levar em consideração o número de amostras, mas voltamos a concluir certo que arbustos A tem em média 2 aranhas, e arbusto B o dobro.

A parte mais interessante, é que apesar de dois processos estarem acontecendo simultaniamente, um que é o ambiente influenciando no número de aranhas que encontramos e o outro é arbustos estarem ocupados ou não, estamos modelando esses dois processos relativamente bem. Resgatando estimativas razoáveis. Se ignoramos um processo, como ignoramos o processo de ocupação no modelo 2, além de péssimas estimativas do tamanho populacional, podemos concluir erroneamente o que esta acontecendo.

Mas ao pensarmos sobre todos os processos que influenciam as nossas observações, podemos selecionar uma estrategia que leve tudo em conta, todos os processos que temos teoria para subsidiar e gerar boas estimativas.

Nesse caso, pensando que o ambiente influencia o número de aranhas que encontramos, mas encontrar aranhas também depende da ocupação, como prevê a teoria de metapopulações, uma regressão de poisson com zero inflados funciona perfeitamente.  E o preço do modelo errado pode ser varias estrelinhas pior, uma conclusão errada do que está acontecendo.