Superdispersão em GLM de Poisson

Este post era inicialmente um email que mandei para o grupo de emails do laboratório. E agora transformei num post*, incluindo algumas contribuições do meu colega Luiz Magnago (gratidão!).

Escrevo porque é uma assunto importante relacionado a análise de dados
e que não vejo muito por aí (e que eu mesmo ignorava até uns meses
atrás… Shame on me, shame! rs) E bom, se eu tiver entendido certo o
que li no Zuur 2009, AKA o Livro dos Pinguins, o que acontece é o
seguinte:

É de conhecimento relativamente comum que, quando trabalhamos com
dados de contagem, uma solução é usar GLM com distribuição de Poisson.
Dados de contagem tendem a seguir a distribuição de Poisson, e muitas
vezes usar ela é melhor do que por exemplo tentar normalizar os dados
usando logaritmo, etc.

O que frequentemente esquecemos é que a distribuição de Poisson também
tem as suas premissas! Sim, não é qualquer dado de contagem que vai
seguir ela necessariamente. Ao contrário da distribuição normal, que
tem um parâmetro pra média e outro pra variância, a distribuição de
Poisson tem um parâmetro só – assumindo assim que a média é igual à
variância.

Mas e se não for? É bem comum termos dados em que a variância é na
verdade maior do que a média (visualize um gráfico de dispersão em
que, quanto mais longe no eixo X, mais espalhados são os valores de Y
sem que a média aumente proporcionalmente). Isso é o que chamamos de
superdispersão – os dados são mais dispersos do que seria esperado
numa distribuição de Poisson.

Como verificar isso? Bom, se eu entendi certo (reli a mesma página
umas 15 vezes pra garantir, rs), depois de rodarmos um GLM de Poisson
no R, rodamos o comando summary sobre os seus resultados. Entre outras
coisas, o summary vai mostrar algo como “Residual deviance: 46.4 on 40
degrees of freedom”. Se não houver superdispersão, o residual deviance
vai ser aproximadamente igual ao número de graus de liberdade – neste
caso estamos bem.

Mas pode aparecer algo como “Residual deviance: 1641.7 on 40 degrees
of freedom” – e isso indicaria superdispersão. Se eu entendi certo, se
o residual deviance for tipo duas vezes maior do que os graus de
liberdade, já estamos com problemas. Em análises que fiz, já vi ele
ser dez ou até cem vezes maior – o que pode trazer sérios problemas às
nossas inferências.

Luiz explica um pouco mais:

No summary do glm sai o tamanho do parâmetro de dispersão que deve ser respeitado para cada distribuição. Por exemplo, “parâmetro de dispersão deve ser igual a 1”. Então essa divisão do Residual deviance pelo degrees of freedom tem que dar próximo de 1. O problema é o quão próximo? A prova real da ausência de overdispersion pode ser apenas vista quando plotada a regressão na forma de y=a+(b*x)+e e não no formato normal y=a+(b*x). Sem o erro (e) adicionado não podemos ver o quanto a overdispersion realmente prejudicou nosso modelo. Modelos com Poisson e Binomial são os mais comprometidos por esse problema.

E que problemas seriam esses? Bom, o Livro dos Pinguins em um momento
mostra resultados de uma seleção de modelos mostrando vários efeitos
significativos; e depois diz que esses resultados devem todos ser
jogados no lixo por causa da superdispersão. Pois a superdispersão
aumenta o erro do tipo I, faz com que variáveis explanatórias que não
têm efeito real apareçam como sendo altamente significativas. Num GAM,
superdispersão pode resultar em padrões muito complexos, com muitos
graus de liberdade efeitos, que não fazem sentido; ao corrigirmos pra
superdispersão, variáveis que tinham um p de 0.00001 de repente passam
a um p de 0.03, e outras que eram altamente significativas perdem toda
a sua significância.

E ninguém quer chegar à conclusão de que uma variável é importante e
descobrir que isso foi tudo um artefato ou um erro estatístico, né?
:-)

Como lidar com isso? Bom, Zuur dá duas soluções relativamente fáceis:
1) incluir um parâmetro para a variância, uma correção no modelo –
isso substituindo family=poisson por family=quasipoisson. Não é uma
distribuição chamada quasipoisson, é apenas uma correção pra
variância. A desvantagem é que não tem AIC associado. 2) usar a
família binomial negativa. Prum glm, usamos library(MASS); glm.nb(y~x)
(não precisa incluir família porque o glm.nb é só pra binomial
negativa); prum gam ajustado no mgcv, usamos family=nb. Binomial
negativa é uma distribuição em que a variância pode ser maior do que a
média (portanto, mais complexa do que a de Poisson). Isso dá um valor
de AIC; a desvantagem? Às vezes os modelos não convergem. E existem
outras possibilidades – podemos ter zeros demais no modelo, e existem
pra isso os zero-inflated models, os quais não me atrevo a tentar
explicar [ainda].

Se não estão convencidas/os, reparem no seguinte. Fiz uma simulação
(script abaixo): simulei dados de contagem usando um processo
qualquer, rodei um GLM de Poisson e deu um p=1.31e-07 – pra 42 pontos
sem nenhuma relação com a variável explanatória! Usando binomial
negativa, deu p=0.425.

library(MASS)
set.seed(17)
a <- abs(round(rnorm(42, 50, 100)))
b <- runif(42)
mod.pois <- glm(a~b, family=poisson)
summary(mod.pois)
mod.nb <- glm.nb(a~b)
summary(mod.nb)

E aqui um exemplo com dados reais, fornecido por Luiz. São dados referentes a vento e dano nas árvores. Ajustando um GLM de Poisson, um quasipoisson e um binomial negativo, os resultados numéricos são os seguintes:

Parâmetros Poisson Quasipoisson Binomial negativa
Intercepto 2.0835 2.0835 1.8082
Erro do intercepto 0.1343 0.4022 0.3541
Inclinação 0.0203 0.0203 0.0272
Erro da inclinação 0.003 0.0089 0.0085
p (inclinação) 6E-12 0.0315 0.00144
Null deviance 252.32 252.32 34.734
Residual deviance 202.16 202.16 27.099
Graus de liberdade 22 22 22
AIC 309.53 NA 185.01
Dispersion parameter 1 8.96 (2.1799) 1

 

Reparem que os parâmetros (intercepto e inclinação) são exatamente os mesmos no Poisson e no quasipoisson. O que muda são os erros associados – e, mudando os erros, a significância também muda. Então, enquanto no Poisson temos um p de 6E-12, ou seja, p=0.000000000006 (6E-12 significa 0 seguido de ponto seguido de onze zeros seguido de 6), no quasipoisson o p é de 0.0315. Dependendo de como interpretamos o p-valor, a interpretação muda, pois o quasipoisson dá evidência relativamente fraca dos efeitos do vento (às vezes usamos p crítico de 0.01, ou 0.005, ou até menor, por exemplo pra corrigir pra testes múltiplos). E residual deviance de 202 sobre 22 graus de liberdade é indicativo de superdispersão. O binomial negativo dá um resultado intermediário, com p=0.001. O modelo binomial negativo diz que o parâmetro de dispersão é igual a 1, mas o parâmetro da distribuição binomial negativa é de 2.18. (Se alguém quiser explicar melhor isso aí, saibam que este blog adora comentários!). O modelo binomial negativo tem um AIC bem menor que o modelo Poisson, portanto é muito mais parcimonioso.

Mas o mais interessante não é isso. O mais interessante é a visualização gráfica.

Plotando os dados e a curva ajustada, temos o seguinte para os três modelos:

overdisp1.png

Sem grandes diferenças, né? Inclusive os dois primeiros vão ser idênticos porque os parâmetros são idênticos.

Mas, plotando na forma y=a+(b*x)+e – ou seja, plotando valores preditos + o erro (para o erro estou usando resíduos padronizados, ou Pearson residuals, mas o gráfico usando deviance residuals foi essencialmente o mesmo), temos:

overdisp2.png

Isso seria a “prova real” da superdispersão. Aquele pontinho aí em cima não é legal. :-)

Script pra fazer esses gráficos pra um dos modelos:

m1<-glm(danos~vento, family=poisson, data=dados1)
summary(m1)

#Plotando os dados brutos e a curva esperada
plot(danos~vento, xlab="Vento (%)", ylab="Copas com danos")
curve(exp(coef(m1)[1]+coef(m1)[2]*x),add=T,lty=1)

#Plotando o esperado mais o erro
predict1<-coef(m1)[1]+(coef(m1)[2]*vento)
res1<-resid(m1, type="pearson")
est1<-exp(predict1+res1)

plot(est1~vento, xlab="Vento (%)", ylab="Copas com danos")
curve(exp(coef(m1)[1]+coef(m1)[2]*x),add=T,lty=1)

Bom, em resumo é isso. Tudo que é ouro não reluz, nem todos que
caminham perdidos estão, o velho que é forte não definha, e nem tudo
que é contagem é Poisson. :-)

* Essa transformação de email em post com adição das figuras etc deve ter demorado umas duas horas. Acham que blogar é fácil? :-)

 

13 pensamentos sobre “Superdispersão em GLM de Poisson

  1. Oi Pav,

    obrigado pelo post e pelas reflexões! De fato, esses são alguns elementos que geralmente consideramos como “asterisco” (aquilo que está em letras pequenas e podemos deixar para ver depois – e frequentemente nunca voltamos para ver isso), mas que na verdade pode ter sérias consequências sobre as inferências e interpretações que fizermos delas (como você mostrou).

    Nesse contexto, me lembrei de um curso dado pelos professores Paulo Inácio Prado e João Baptista, na USP, sobre modelagem estatística utilizando a abordagem da máxima verossimilhança (http://cmq.esalq.usp.br/BIE5781/doku.php) – eu sou fã deles, e recomendo muito o curso!

    Tem um tópico em que eles explicam, por texto, exemplos e código, a aplicação de GLMs Poisson, com um comentário sobre o uso de um GLM com binomial negativa. Acho que vale bastante a pena a leitura, como complemento ao post: http://cmq.esalq.usp.br/BIE5781/doku.php?id=05-binomial-poisson:05-binomial-poisson
    (tem um jeito mais bonito de compartilhar link por aqui?)

    Aliás, como eles se propõe ao exercício de escrever as funções de (log) verossimilhança (negativa) de maneira explícita, isso pode ser um ponto a mais para eu recomendar a leitura, uma vez que isso deixa mais claro os pormenores por trás da análise, que geralmente ficam obscuros nos mistérios da função ‘glm’ do R.

    Abraços

    Curtir

  2. Pingback: Machismo estatístico (tradução) – Mais Um Blog de Ecologia e Estatística

  3. Pingback: Interpretando resultados – Mais Um Blog de Ecologia e Estatística

  4. Pingback: Fazendo pesquisa em quarentena – Mais Um Blog de Ecologia e Estatística

  5. Prezado Pavel Dodonov,

    Como sugerido no blog, você abriu a oportunidade de fazermos comentários ou discussão sobre o parâmetro de dispersão da Binomial negativa que é igual a 1 e no resultado ele indicou o valor 2.18.

    Primeiramente parabéns ao blog. Muito legal conduzir suas perspectivas sobre o assunto de forma descontraída, muito bem colocado seu ponto de vista.

    A função glm.nb(y~x) retorna como output o seguinte resultado “Dispersion parameter for Negative Binomial(2.18) family taken to be 1”. Esse valor (2.18) é a estimativa do parâmetro de dispersão da família exponencial de distribuições de probabilidade que normalmente (para um modelo da Binomial Negativa, sem a superdispersão) é igual a 1 (assim como a distribuição Bernoulli e Poisson também). Entretanto, por vários motivos (que aqui no seu blog você indica alguns) há a ocorrência da superdispersão e, portanto, irá diferir de 1. Como houve uma parametrização da variância (método de quasiverossimilhança) muitas vezes (se essa parametrização for linear, como é normalmente) podemos concluir, nesse caso, que a dispersão é maior 2.18 vezes do que o valor esperado.

    Portanto, se você fatorar a distribuição Binomial Negativa na forma da família exponencial, você irá concluir que: o parâmetro de dispersão phi=1, o parâmetro canônico (ou natural) theta= ln(um/(um+k)), a função do parâmetro b(theta)=-kln(1-e^{theta}) e a função c(y,phi)= ln( (gama(k+y))/(gama(y)y!) ), e consequentemente entenderá o porquê que ele teoricamente é igual a 1, no entanto foi estimado 2.18 no modelo proposto.

    Atenciosamente Carlos Antônio Zarzar (carloszarzar_@hotmail.com)

    Curtir

  6. Pingback: Interpretando os resultados de um GLM – Mais Um Blog de Ecologia e Estatística

  7. Quando utilizo o quase-poisson, como avalio o ajuste do modelo? Estou tentando entender um pouco mais sobre o chi2gof cells (numlist) no Stata. Saberia explicar?

    Curtir

  8. Oi Pavel!! Que post lindo e como sempre, muito esclarecedor. Parabéns pela disposição de continuar traduzindo a estatística e isso faz toda diferença pra nós estudantes = )
    Grande abraço!!

    Curtir

  9. PREZADO PAVEL DODONOV,

    Excelente e acessível, na forma de trajetória da descoberta, o que facilitou muito o entendimento. Obrigado.

    Rodei um modelo poisson bem simples ACIDENTES = b0 + b1PREVENCAO, em que PREVENCAO é uma variável dummy para a existência de um programa de prevenção e, ao contrário de você, deparei-me com SUBDISPERSÃO, enquanto o output alertava:
    “(Dispersion parameter for poisson family taken to be 1)”

    Fator de dispersão:
    mod1$deviance/mod1$df.residual = 0.664928. (Sub dispersão)

    Coefficients MOD1:
    Estimate Std. Error z value Pr(>|z|)
    (Intercept) 1.31283 0.04735 27.73 <2e-16 (Olho nos erros)
    progCom -2.00597 0.14408 -13.92 |z|)
    (Intercept) 1.31283 0.04735 27.73 <2e-16 *** (Erros maiores agora)
    progCom -2.00597 0.14408 -13.92 <2e-16 ***

    – Então a subdispersão aumenta a chance do ERRO TIPO II?
    – Nesse caso, ainda seria aplicável a distribuição binomial negativa?

    Obrigado
    ALESSANDRO
    alessandro.correa4@gmail.com

    Curtir

    • Oi Alessandro,
      Muito obrigado pelo comentário!
      Então… O livro de Zuur et al. fala que subdispersão é algo raro e eu nem me atentei muito a isso. Percebi pela primeira vez ano passado, durante uma aula de estatística, em que realizei um GLM de Poisson que não deu nem um pouco significativo mas o teste de Kruskal-Wallis deu significativo, e eu acho que pelos gráficos havia indícios de haver diferença.

      Éis o que Zuur fala:
      “If phi <1, we have underdispersion. The latter means that the variance of the response variable is smaller than you would expect from a Poisson distribution. Reasons for underdispersion are the modl is fitting a couple of outliers rather too well or there are too many explanatory variables or interactions in the model (overfitting). If this is not the case, then the consensus is not to correct for underdispersion. Models that take underdispersion into account are discussed in Chapter 7 oif Hilbe 2007 Negative Binomial Regression. Cambridge.

      E éis o que o livro Hilbe 2013 Methods of Statistical Model Estimation fala sobre underdispersion:
      "Under-dispersed Poisson data may be handled by scaling, by use of a hurdle model, by generalized Poisson or by using a generalized negative binomial model."

      Hurdle model é um tipo de modelo para inflação por zeros, nos quais presença-ausência é ajustada como um modelo binomial e a abundância onde há presença é ajustada como um zero-truncated Poisson.

      Então, minha sugestão seria… Para um modelo simples como o seu, com uma única variável explanatória categórica, use um teste não-paramétrico como Mann-Whitney ou Kruskal-Wallis, cujas propriedades são bem conhecidas e que no geral têm se mostrado confiáveis (até onde sei). Para coisas mais complexas, veja essas possibilidades aí; elas estão acima no meu atual nível de conhecimento. :-)

      Abraço!

      Curtir

Deixe um comentário

Preencha os seus dados abaixo ou clique em um ícone para log in:

Logotipo do WordPress.com

Você está comentando utilizando sua conta WordPress.com. Sair /  Alterar )

Foto do Google

Você está comentando utilizando sua conta Google. Sair /  Alterar )

Imagem do Twitter

Você está comentando utilizando sua conta Twitter. Sair /  Alterar )

Foto do Facebook

Você está comentando utilizando sua conta Facebook. Sair /  Alterar )

Conectando a %s