Modelo de Poisson com Superdispersão

Quase-verossimilhança e Modelo Binomial Negativo com aplicações em R

Author

Prof. Dr. Dennison Carvalho - Baseado em Dobson & Barnett (2018)

Sobre esta apostila. O modelo de Poisson assume que média e variância são iguais. Na prática, dados de contagem quase sempre têm variância maior do que a média — fenômeno chamado superdispersão. Esta apostila apresenta o diagnóstico de superdispersão e duas soluções: o modelo quasi-Poisson (quasi-verossimilhança) e o modelo de Binomial Negativa. A ênfase é na prática em R com interpretação detalhada dos resultados.

1 Objetivos

Esta apostila apresenta, de forma prática e didática, três abordagens para dados de contagem:

  1. Modelo de Poisson;
  2. Modelo de Poisson com superdispersão, via quase-verossimilhança;
  3. Modelo Binomial Negativo.

A ênfase será na aplicação em R, na interpretação dos resultados e na avaliação da qualidade do ajuste.

2 Contexto geral

Muitos problemas em Estatística envolvem uma variável resposta que representa uma contagem. Exemplos:

  • número de acidentes em uma rodovia;
  • número de casos de uma doença;
  • número de defeitos em uma peça;
  • número de espécies observadas em uma área;
  • número de publicações de um pesquisador;
  • número de sinistros em seguros.

Quando a resposta é uma contagem, um primeiro modelo natural é o modelo de Poisson.

3 Modelo de Poisson

Considere uma variável aleatória de contagem:

\[ Y_i = \text{número de eventos observados na unidade } i. \]

No modelo de Poisson, assumimos:

\[ Y_i \sim \text{Poisson}(\mu_i), \]

em que:

\[ E(Y_i)=\mu_i \]

e

\[ \operatorname{Var}(Y_i)=\mu_i. \]

Essa igualdade é muito importante.

Importante

No modelo de Poisson clássico, a média e a variância são iguais:

\[ E(Y_i)=\operatorname{Var}(Y_i)=\mu_i. \]

3.1 Regressão de Poisson

Na regressão de Poisson, relacionamos a média esperada \(\mu_i\) com variáveis explicativas. A ligação mais usada é a ligação log:

\[ \log(\mu_i)=\beta_0+\beta_1x_{i1}+\cdots+\beta_px_{ip}. \]

Aplicando exponencial:

\[ \mu_i=\exp(\beta_0+\beta_1x_{i1}+\cdots+\beta_px_{ip}). \]

3.2 Por que usar ligação log?

A ligação log é usada porque:

  1. garante que \(\mu_i>0\);
  2. é a ligação canônica da Poisson;
  3. permite interpretação multiplicativa dos efeitos;
  4. é adequada quando os efeitos são proporcionais, e não apenas aditivos.

Por exemplo, se o modelo é:

\[ \log(\mu_i)=\beta_0+\beta_1x_i, \]

então:

\[ \mu_i=\exp(\beta_0)\exp(\beta_1x_i). \]

Quando \(x\) aumenta uma unidade, a média esperada é multiplicada por:

\[ \exp(\beta_1). \]

4 O problema da superdispersão

A principal limitação do modelo de Poisson é assumir que:

\[ \operatorname{Var}(Y_i)=\mu_i. \]

Na prática, muitas vezes observamos:

\[ \operatorname{Var}(Y_i)>\mu_i. \]

Esse fenômeno é chamado de superdispersão.

4.1 O que é superdispersão?

Existe superdispersão quando a variabilidade dos dados é maior do que a prevista pelo modelo de Poisson.

Em termos simples:

O modelo de Poisson espera uma certa variabilidade, mas os dados variam mais do que ele consegue explicar.

4.2 Consequências da superdispersão

Se ignoramos a superdispersão:

  • os erros-padrão podem ficar subestimados;
  • os testes podem parecer significativos demais;
  • os intervalos de confiança ficam estreitos demais;
  • podemos concluir que uma variável é importante quando, na verdade, a evidência não é tão forte.
Cuidado

A superdispersão não costuma alterar muito os coeficientes estimados no modelo de Poisson, mas pode alterar bastante os erros-padrão, os valores-p e as conclusões inferenciais.

5 Exemplo prático em R

Vamos construir um exemplo simples. Suponha que estamos estudando o número de atendimentos em diferentes unidades de saúde.

A variável resposta será:

\[ Y_i = \text{número de atendimentos na unidade } i. \]

As variáveis explicativas serão:

  • tempo: tempo de funcionamento da unidade, em anos;
  • tipo: tipo da unidade, A ou B.

Para fins didáticos, vamos simular dados com superdispersão.

Code 
set.seed(123)

n <- 120

tempo <- runif(n, min = 1, max = 10)
tipo <- factor(sample(c("A", "B"), n, replace = TRUE))

# Média verdadeira usada para simular os dados
eta <- 1.2 + 0.12 * tempo + 0.45 * (tipo == "B")
mu <- exp(eta)

# Gerando dados com superdispersão usando Binomial Negativa
# Quanto menor o size, maior a superdispersão.
y <- rnbinom(n, size = 3, mu = mu)

dados <- data.frame(y, tempo, tipo)

head(dados)
   y    tempo tipo
1  3 3.588198    A
2 16 8.094746    B
3 14 4.680792    A
4  5 8.947157    A
5 12 9.464206    A
6 13 1.410008    B

6 Análise exploratória

Antes de ajustar modelos, devemos olhar os dados.

Code 
summary(dados)
       y              tempo       tipo  
 Min.   : 0.000   Min.   :1.006   A:53  
 1st Qu.: 4.000   1st Qu.:3.451   B:67  
 Median : 8.000   Median :5.323         
 Mean   : 9.108   Mean   :5.597         
 3rd Qu.:13.000   3rd Qu.:7.880         
 Max.   :46.000   Max.   :9.948

6.1 Gráfico da contagem contra o tempo

Code 
plot(dados$tempo, dados$y,
     pch = 19,
     xlab = "Tempo de funcionamento",
     ylab = "Número de atendimentos",
     main = "Contagem observada versus tempo")

6.2 Boxplot por tipo de unidade

Code 
boxplot(y ~ tipo,
        data = dados,
        xlab = "Tipo de unidade",
        ylab = "Número de atendimentos",
        main = "Contagem por tipo de unidade")

6.3 Média e variância observadas

Code 
mean(dados$y)
[1] 9.108333
Code 
var(dados$y)
[1] 49.91254

Se a variância amostral for muito maior que a média amostral, isso já sugere possível superdispersão. Essa comparação é apenas exploratória, pois a média também muda com as covariáveis.

7 Ajuste do modelo de Poisson

Vamos começar com o modelo de Poisson usual:

\[ Y_i \sim \text{Poisson}(\mu_i) \]

com

\[ \log(\mu_i)=\beta_0+\beta_1\text{tempo}_i+\beta_2\text{tipoB}_i. \]

Code 
mod_pois <- glm(y ~ tempo + tipo,
                family = poisson(link = "log"),
                data = dados)

summary(mod_pois)

Call:
glm(formula = y ~ tempo + tipo, family = poisson(link = "log"), 
    data = dados)

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  1.25736    0.09507   13.22  < 2e-16 ***
tempo        0.12384    0.01220   10.15  < 2e-16 ***
tipoB        0.35424    0.06292    5.63 1.81e-08 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 584.05  on 119  degrees of freedom
Residual deviance: 453.89  on 117  degrees of freedom
AIC: 908.7

Number of Fisher Scoring iterations: 5

7.1 Interpretação dos coeficientes

No modelo com ligação log, os coeficientes estão na escala do logaritmo da média.

Para interpretar na escala original, usamos a exponencial dos coeficientes:

Code 
coef(mod_pois)
(Intercept)       tempo       tipoB 
  1.2573632   0.1238385   0.3542391
Code 
exp(coef(mod_pois))
(Intercept)       tempo       tipoB 
   3.516138    1.131833    1.425096

7.1.1 Intercepto

O intercepto representa o logaritmo da média esperada quando:

  • tempo = 0;
  • tipo = A.

Como tempo = 0 pode estar fora da faixa observada, o intercepto nem sempre tem interpretação prática direta.

7.1.2 Coeficiente de tempo

Se \(\hat\beta_1\) é o coeficiente de tempo, então:

\[ \exp(\hat\beta_1) \]

é o fator multiplicativo na média esperada para cada aumento de uma unidade em tempo, mantendo o tipo de unidade fixo.

Exemplo de interpretação:

Para cada ano adicional de funcionamento, o número médio esperado de atendimentos é multiplicado por \(\exp(\hat\beta_1)\), mantendo o tipo de unidade constante.

7.1.3 Coeficiente de tipoB

Se \(\hat\beta_2\) é o coeficiente de tipoB, então:

\[ \exp(\hat\beta_2) \]

é a razão entre a média esperada do tipo B e a média esperada do tipo A, mantendo tempo constante.

Exemplo de interpretação:

Unidades do tipo B apresentam, em média, \(\exp(\hat\beta_2)\) vezes o número esperado de atendimentos das unidades do tipo A, para o mesmo tempo de funcionamento.

8 Avaliando a superdispersão

Uma forma simples de avaliar superdispersão é calcular:

\[ \hat\phi_D = \frac{D}{gl}, \]

em que:

  • \(D\) é a deviance residual;
  • \(gl\) são os graus de liberdade residuais.

Também podemos usar:

\[ \hat\phi_P = \frac{X^2}{gl}, \]

em que \(X^2\) é a soma dos resíduos de Pearson ao quadrado.

Code 
# Deviance residual dividida pelos graus de liberdade
phi_deviance <- deviance(mod_pois) / df.residual(mod_pois)
phi_deviance
[1] 3.879426
Code 
# Estatística de Pearson dividida pelos graus de liberdade
res_pearson <- residuals(mod_pois, type = "pearson")
phi_pearson <- sum(res_pearson^2) / df.residual(mod_pois)
phi_pearson
[1] 3.9918

8.1 Interpretação

Regra prática:

  • valor próximo de 1: sem forte evidência de superdispersão;
  • valor muito maior que 1: possível superdispersão;
  • valor menor que 1: possível subdispersão.
Dica

Na prática, valores como 1,5, 2 ou maiores já merecem atenção. Não existe uma fronteira universal, mas quanto maior o valor, maior a evidência de superdispersão.

9 Teste aproximado de bondade do ajuste

Podemos calcular um p-valor aproximado usando a deviance residual:

Code 
p_dev <- 1 - pchisq(deviance(mod_pois), df = df.residual(mod_pois))
p_dev
[1] 0

E outro usando a estatística de Pearson:

Code 
X2 <- sum(residuals(mod_pois, type = "pearson")^2)
p_pearson <- 1 - pchisq(X2, df = df.residual(mod_pois))
p_pearson
[1] 0

Se os p-valores forem muito pequenos, isso sugere que o modelo de Poisson pode não estar ajustando adequadamente os dados.

10 Gráficos de diagnóstico

10.1 Resíduos de Pearson versus valores ajustados

Code 
plot(fitted(mod_pois), residuals(mod_pois, type = "pearson"),
     pch = 19,
     xlab = "Valores ajustados",
     ylab = "Resíduos de Pearson",
     main = "Poisson: resíduos de Pearson vs ajustados")
abline(h = 0, lty = 2)

10.2 Resíduos deviance versus valores ajustados

Code 
plot(fitted(mod_pois), residuals(mod_pois, type = "deviance"),
     pch = 19,
     xlab = "Valores ajustados",
     ylab = "Resíduos deviance",
     main = "Poisson: resíduos deviance vs ajustados")
abline(h = 0, lty = 2)

10.3 Valores observados versus ajustados

Code 
plot(fitted(mod_pois), dados$y,
     pch = 19,
     xlab = "Valores ajustados",
     ylab = "Valores observados",
     main = "Poisson: observados vs ajustados")
abline(0, 1, lty = 2)

11 Modelo de quase-verossimilhança

Quando há superdispersão, uma solução simples é usar o modelo quasi-Poisson.

No modelo quasi-Poisson, mantemos a média:

\[ E(Y_i)=\mu_i, \]

mas relaxamos a variância:

\[ \operatorname{Var}(Y_i)=\phi\mu_i. \]

Aqui, \(\phi\) é o parâmetro de dispersão.

  • Se \(\phi=1\), voltamos à Poisson usual;
  • Se \(\phi>1\), há superdispersão;
  • Se \(\phi<1\), há subdispersão.

11.1 Demonstração simples da consequência nos erros-padrão

No modelo de Poisson, aproximadamente:

\[ \operatorname{Var}(\hat\beta) \approx (X^TWX)^{-1}. \]

No modelo quasi-Poisson:

\[ \operatorname{Var}(\hat\beta) \approx \phi (X^TWX)^{-1}. \]

Assim, os erros-padrão são multiplicados por:

\[ \sqrt{\phi}. \]

Importante

O quasi-Poisson geralmente mantém os mesmos coeficientes do modelo de Poisson, mas aumenta os erros-padrão quando há superdispersão. Isso torna os testes mais cautelosos.

11.2 Ajuste no R

Code 
mod_quasi <- glm(y ~ tempo + tipo,
                 family = quasipoisson(link = "log"),
                 data = dados)

summary(mod_quasi)

Call:
glm(formula = y ~ tempo + tipo, family = quasipoisson(link = "log"), 
    data = dados)

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  1.25736    0.18995   6.619 1.14e-09 ***
tempo        0.12384    0.02437   5.082 1.43e-06 ***
tipoB        0.35424    0.12572   2.818  0.00568 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for quasipoisson family taken to be 3.9918)

    Null deviance: 584.05  on 119  degrees of freedom
Residual deviance: 453.89  on 117  degrees of freedom
AIC: NA

Number of Fisher Scoring iterations: 5

11.3 Comparando Poisson e quasi-Poisson

Code 
summary(mod_pois)$coefficients
             Estimate Std. Error   z value     Pr(>|z|)
(Intercept) 1.2573632 0.09507397 13.225104 6.285022e-40
tempo       0.1238385 0.01219727 10.152972 3.214229e-24
tipoB       0.3542391 0.06292490  5.629554 1.806767e-08
Code 
summary(mod_quasi)$coefficients
             Estimate Std. Error  t value     Pr(>|t|)
(Intercept) 1.2573632 0.18995295 6.619340 1.144654e-09
tempo       0.1238385 0.02436952 5.081697 1.432714e-06
tipoB       0.3542391 0.12572074 2.817666 5.680702e-03

Observe que os coeficientes estimados tendem a ser iguais ou muito próximos. A diferença principal aparece nos erros-padrão, estatísticas de teste e valores-p.

Code 
# Parâmetro de dispersão estimado no quasi-Poisson
summary(mod_quasi)$dispersion
[1] 3.9918

11.4 Interpretação prática

Se o parâmetro de dispersão estimado for, por exemplo, \(\hat\phi=3\), isso significa que a variância dos dados é aproximadamente três vezes maior do que a esperada pela Poisson.

Nesse caso, usar o quasi-Poisson evita que os erros-padrão fiquem artificialmente pequenos.

12 Modelo Binomial Negativo

Outra forma muito usada para lidar com superdispersão é o modelo Binomial Negativo.

Uma parametrização comum é:

\[ E(Y_i)=\mu_i \]

e

\[ \operatorname{Var}(Y_i)=\mu_i+\frac{\mu_i^2}{\theta}. \]

Aqui, \(\theta\) é um parâmetro de dispersão.

12.1 Interpretação da variância

A variância da Binomial Negativa é maior que a média:

\[ \operatorname{Var}(Y_i)=\mu_i+\frac{\mu_i^2}{\theta}. \]

Como o termo

\[ \frac{\mu_i^2}{\theta} \]

é positivo, temos:

\[ \operatorname{Var}(Y_i)>\mu_i. \]

Ou seja, a Binomial Negativa acomoda naturalmente superdispersão.

12.2 Demonstração didática: por que a Binomial Negativa gera superdispersão?

Imagine que, condicionalmente a uma taxa \(\Lambda_i\), temos:

\[ Y_i|\Lambda_i \sim \text{Poisson}(\Lambda_i). \]

Mas agora suponha que a própria taxa \(\Lambda_i\) varie entre unidades. Então, pela lei da variância total:

\[ \operatorname{Var}(Y_i)=E[\operatorname{Var}(Y_i|\Lambda_i)] + \operatorname{Var}[E(Y_i|\Lambda_i)]. \]

Como, para Poisson:

\[ E(Y_i|\Lambda_i)=\Lambda_i \]

and

\[ \operatorname{Var}(Y_i|\Lambda_i)=\Lambda_i, \]

segue que:

\[ \operatorname{Var}(Y_i)=E(\Lambda_i)+\operatorname{Var}(\Lambda_i). \]

Se \(E(\Lambda_i)=\mu_i\), então:

\[ \operatorname{Var}(Y_i)=\mu_i+\operatorname{Var}(\Lambda_i). \]

Como \(\operatorname{Var}(\Lambda_i)>0\), a variância fica maior que a média. Essa é a ideia central da superdispersão.

13 Ajuste da Binomial Negativa no R

Para ajustar a Binomial Negativa, usaremos a função glm.nb() do pacote MASS.

Code 
library(MASS)

mod_nb <- glm.nb(y ~ tempo + tipo,
                 data = dados)

summary(mod_nb)

Call:
glm.nb(formula = y ~ tempo + tipo, data = dados, init.theta = 3.150441707, 
    link = log)

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  1.24747    0.17147   7.275 3.46e-13 ***
tempo        0.12379    0.02373   5.216 1.82e-07 ***
tipoB        0.37104    0.12280   3.022  0.00251 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for Negative Binomial(3.1504) family taken to be 1)

    Null deviance: 163.71  on 119  degrees of freedom
Residual deviance: 129.57  on 117  degrees of freedom
AIC: 736.25

Number of Fisher Scoring iterations: 1

              Theta:  3.150 
          Std. Err.:  0.573 

 2 x log-likelihood:  -728.248

13.1 Coeficientes exponenciados

Code 
exp(coef(mod_nb))
(Intercept)       tempo       tipoB 
   3.481533    1.131783    1.449245

A interpretação dos coeficientes é semelhante à regressão de Poisson com ligação log.

  • \(\exp(\hat\beta_1)\): fator multiplicativo na média esperada para cada aumento de uma unidade em tempo;
  • \(\exp(\hat\beta_2)\): razão entre as médias esperadas dos grupos, mantendo as demais variáveis constantes.

13.2 Parâmetro theta

Code 
mod_nb$theta
[1] 3.150442

No modelo Binomial Negativo, o parâmetro theta controla a superdispersão. Valores menores de theta indicam maior superdispersão.

14 Comparação entre os modelos

14.1 Coeficientes

Code 
coef(mod_pois)
(Intercept)       tempo       tipoB 
  1.2573632   0.1238385   0.3542391
Code 
coef(mod_quasi)
(Intercept)       tempo       tipoB 
  1.2573632   0.1238385   0.3542391
Code 
coef(mod_nb)
(Intercept)       tempo       tipoB 
  1.2474726   0.1237940   0.3710425

14.2 Coeficientes exponenciados

Code 
exp(coef(mod_pois))
(Intercept)       tempo       tipoB 
   3.516138    1.131833    1.425096
Code 
exp(coef(mod_quasi))
(Intercept)       tempo       tipoB 
   3.516138    1.131833    1.425096
Code 
exp(coef(mod_nb))
(Intercept)       tempo       tipoB 
   3.481533    1.131783    1.449245

14.3 AIC

O AIC pode ser usado para comparar Poisson e Binomial Negativa, mas não é apropriado para quasi-Poisson, pois modelos de quase-verossimilhança não têm uma verossimilhança completa no mesmo sentido.

Code 
AIC(mod_pois, mod_nb)
         df      AIC
mod_pois  3 908.6995
mod_nb    4 736.2479
Cuidado

Não use AIC para comparar diretamente modelos quasi-Poisson com Poisson ou Binomial Negativa, pois o quasi-Poisson é baseado em quase-verossimilhança.

15 Comparando valores ajustados

Code 
dados$fit_pois <- fitted(mod_pois)
dados$fit_quasi <- fitted(mod_quasi)
dados$fit_nb <- fitted(mod_nb)

head(dados)
   y    tempo tipo  fit_pois fit_quasi    fit_nb
1  3 3.588198    A  5.483372  5.483372  5.428539
2 16 8.094746    B 13.654150 13.654150 13.743918
3 14 4.680792    A  6.277837  6.277837  6.214758
4  5 8.947157    A 10.647931 10.647931 10.538943
5 12 9.464206    A 11.352026 11.352026 11.235572
6 13 1.410008    B  5.966829  5.966829  6.007843
Code 
plot(dados$y, dados$fit_pois,
     pch = 19,
     col = "blue",
     xlab = "Observado",
     ylab = "Ajustado",
     main = "Valores observados e ajustados")
points(dados$y, 
       dados$fit_nb, 
       pch = 17,
       col = "red4")
abline(0, 1, lty = 2)
legend("topright",
       legend = c("Poisson", "Binomial Negativa"),
       pch = c(19, 17),
       bty = "n",
       col = c("blue","red4"))

16 Predição com os modelos

Vamos prever a média esperada para uma unidade:

  • tempo = 5 anos;
  • tipo = B.
Code 
novo <- data.frame(tempo = 5,
                   tipo = factor("B", levels = levels(dados$tipo)))

predict(mod_pois, newdata = novo, type = "response")
      1 
9.30726
Code 
predict(mod_quasi, newdata = novo, type = "response")
      1 
9.30726
Code 
predict(mod_nb, newdata = novo, type = "response")
       1 
9.369739

A opção type = "response" devolve a predição na escala da média, isto é, no número esperado de eventos.

17 Intervalos de confiança para os efeitos

17.1 Poisson

Code 
confint(mod_pois)
                2.5 %    97.5 %
(Intercept) 1.0687853 1.4415155
tempo       0.1000183 0.1478396
tipoB       0.2315278 0.4782675
Code 
exp(confint(mod_pois))
               2.5 %   97.5 %
(Intercept) 2.911840 4.227097
tempo       1.105191 1.159327
tipoB       1.260524 1.613277

17.2 Quasi-Poisson

Para quasi-Poisson, podemos construir intervalos usando os erros-padrão do resumo do modelo.

Code 
est <- coef(mod_quasi)
se <- summary(mod_quasi)$coefficients[, "Std. Error"]

IC_quasi <- cbind(
  LI = est - 1.96 * se,
  LS = est + 1.96 * se
)

IC_quasi
                    LI        LS
(Intercept) 0.88505539 1.6296709
tempo       0.07607425 0.1716028
tipoB       0.10782644 0.6006517
Code 
exp(IC_quasi)
                  LI       LS
(Intercept) 2.423119 5.102196
tempo       1.079043 1.187206
tipoB       1.113854 1.823307

17.3 Binomial Negativa

Code 
confint(mod_nb)
                 2.5 %    97.5 %
(Intercept) 0.91173235 1.5859358
tempo       0.07781264 0.1700241
tipoB       0.12954394 0.6121600
Code 
exp(confint(mod_nb))
               2.5 %   97.5 %
(Intercept) 2.488630 4.883859
tempo       1.080920 1.185333
tipoB       1.138309 1.844411

18 Quando usar cada abordagem?

18.1 Poisson

Use quando:

  • a resposta é contagem;
  • a variância é aproximadamente igual à média, após considerar as covariáveis;
  • não há evidência forte de superdispersão.

18.2 Quasi-Poisson

Use quando:

  • há superdispersão;
  • seu foco principal é corrigir erros-padrão e testes;
  • você não precisa comparar modelos por AIC.

18.3 Binomial Negativa

Use quando:

  • há superdispersão importante;
  • você quer um modelo probabilístico completo;
  • deseja comparar modelos por AIC;
  • há heterogeneidade não explicada entre unidades.

19 Resumo final

Modelo Média Variância Uso principal
Poisson \(\mu\) \(\mu\) Contagens sem superdispersão
Quasi-Poisson \(\mu\) \(\phi\mu\) Corrigir erros-padrão
Binomial Negativa \(\mu\) \(\mu+\mu^2/\theta\) Modelar superdispersão explicitamente

20 Script final compacto

O bloco abaixo reúne a parte principal da análise.

Code 
# Dados simulados
set.seed(123)
n <- 120
tempo <- runif(n, 1, 10)
tipo <- factor(sample(c("A", "B"), n, replace = TRUE))
eta <- 1.2 + 0.12 * tempo + 0.45 * (tipo == "B")
mu <- exp(eta)
y <- rnbinom(n, size = 3, mu = mu)
dados <- data.frame(y, tempo, tipo)

# Poisson
mod_pois <- glm(y ~ tempo + tipo,
                family = poisson(link = "log"),
                data = dados)

# Diagnóstico de superdispersão
phi_dev <- deviance(mod_pois) / df.residual(mod_pois)
phi_pear <- sum(residuals(mod_pois, type = "pearson")^2) / df.residual(mod_pois)

phi_dev
[1] 3.879426
Code 
phi_pear
[1] 3.9918
Code 
# Quasi-Poisson
mod_quasi <- glm(y ~ tempo + tipo,
                 family = quasipoisson(link = "log"),
                 data = dados)

# Binomial Negativa
library(MASS)
mod_nb <- glm.nb(y ~ tempo + tipo, data = dados)

# Comparação dos resultados
summary(mod_pois)

Call:
glm(formula = y ~ tempo + tipo, family = poisson(link = "log"), 
    data = dados)

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  1.25736    0.09507   13.22  < 2e-16 ***
tempo        0.12384    0.01220   10.15  < 2e-16 ***
tipoB        0.35424    0.06292    5.63 1.81e-08 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 584.05  on 119  degrees of freedom
Residual deviance: 453.89  on 117  degrees of freedom
AIC: 908.7

Number of Fisher Scoring iterations: 5
Code 
summary(mod_quasi)

Call:
glm(formula = y ~ tempo + tipo, family = quasipoisson(link = "log"), 
    data = dados)

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  1.25736    0.18995   6.619 1.14e-09 ***
tempo        0.12384    0.02437   5.082 1.43e-06 ***
tipoB        0.35424    0.12572   2.818  0.00568 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for quasipoisson family taken to be 3.9918)

    Null deviance: 584.05  on 119  degrees of freedom
Residual deviance: 453.89  on 117  degrees of freedom
AIC: NA

Number of Fisher Scoring iterations: 5
Code 
summary(mod_nb)

Call:
glm.nb(formula = y ~ tempo + tipo, data = dados, init.theta = 3.150441707, 
    link = log)

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  1.24747    0.17147   7.275 3.46e-13 ***
tempo        0.12379    0.02373   5.216 1.82e-07 ***
tipoB        0.37104    0.12280   3.022  0.00251 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for Negative Binomial(3.1504) family taken to be 1)

    Null deviance: 163.71  on 119  degrees of freedom
Residual deviance: 129.57  on 117  degrees of freedom
AIC: 736.25

Number of Fisher Scoring iterations: 1

              Theta:  3.150 
          Std. Err.:  0.573 

 2 x log-likelihood:  -728.248
Code 
# Efeitos multiplicativos
exp(coef(mod_pois))
(Intercept)       tempo       tipoB 
   3.516138    1.131833    1.425096
Code 
exp(coef(mod_quasi))
(Intercept)       tempo       tipoB 
   3.516138    1.131833    1.425096
Code 
exp(coef(mod_nb))
(Intercept)       tempo       tipoB 
   3.481533    1.131783    1.449245
Code 
# Comparação por AIC apenas entre modelos com verossimilhança
AIC(mod_pois, mod_nb)
         df      AIC
mod_pois  3 908.6995
mod_nb    4 736.2479

21 Conclusão

Dados de contagem são frequentemente analisados com modelos de Poisson. Entretanto, a suposição de igualdade entre média e variância pode ser restritiva. Quando a variância observada é maior do que a prevista pelo modelo, temos superdispersão.

Nessa situação, o quasi-Poisson é uma alternativa simples para corrigir erros-padrão e inferências. Já a Binomial Negativa oferece um modelo probabilístico mais flexível, com uma estrutura de variância que cresce mais rapidamente que a média.

Na prática, uma boa estratégia é:

  1. ajustar Poisson;
  2. verificar superdispersão;
  3. se houver superdispersão, ajustar quasi-Poisson e Binomial Negativa;
  4. comparar interpretações, erros-padrão, resíduos e, quando possível, AIC.