Da Poisson à binomial negativa

Você ajusta um MLG de Poisson a dados de contagem, olha o diagnóstico e descobre que a variância dos dados é três, cinco, dez vezes a média. O manual diz: “há sobredispersão, use a binomial negativa”. Este artigo explica por quê — não como receita, mas como consequência de um raciocínio. A ideia central é escrever a contagem em duas etapas: primeiro a natureza sorteia uma taxa \(Z\) para cada unidade; depois, dada a taxa, a contagem é Poisson. Com duas ferramentas — a esperança iterada \(\mathbb E[\mathbb E[Y \mid Z]]\) e a variância total — e uma integral, a binomial negativa aparece sozinha. No fim, ajustamos os dois modelos a um mesmo conjunto de dados e medimos, por simulação, quanto custa ignorar a sobredispersão. Todos os números saem de um script reprodutível.

§1O ponto de partida: o MLG de Poisson

Dados de contagem — consultas médicas por ano, sinistros por apólice, defeitos por peça, gols por jogo — são inteiros não negativos. O modelo mais simples para eles é a Poisson:

Num MLG de Poisson, cada observação \(i\) tem sua própria média, ligada às covariáveis pelo logaritmo:

A Poisson tem uma propriedade que define tudo o que vem a seguir: a variância é igual à média,

Isso se chama equidispersão. É uma restrição forte: o modelo não tem nenhum parâmetro livre para a variância. Uma vez escolhida a média, a variância está decidida. Na prática, dados reais quase sempre têm variância maior que a média. Esse é o fenômeno da sobredispersão.

§2Como a sobredispersão aparece

Se o modelo (2) estiver certo, os resíduos de Pearson \(r_i = (y_i - \hat\mu_i)/\sqrt{\hat\mu_i}\) têm variância perto de 1, e a soma dos seus quadrados fica perto dos graus de liberdade. Daí a estatística de dispersão:

em que \(p\) é o número de coeficientes. Sob o modelo Poisson, \(\hat\phi \approx 1\). Um \(\hat\phi\) bem acima de 1 indica sobredispersão (a razão deviance/graus de liberdade conta a mesma história). No exemplo que vamos analisar na §11 — 500 contagens com uma covariável — o ajuste de Poisson dá \(\hat\phi = 3{,}20\): a variância dos dados é cerca de três vezes o que a Poisson admite.

Um gráfico ajuda mais que um número. Ordene as observações pela média ajustada, separe em grupos e, em cada grupo, compare a média com a variância amostral. Sob Poisson, os pontos ficariam em torno da reta \(\operatorname{Var} = \mu\):

Os pontos estão todos acima da reta, e cada vez mais acima conforme a média cresce. A curva cheia,\(\mu + \mu^2/\theta\), é a variância da binomial negativa — chegaremos a ela na §9. Antes, a pergunta importante: por que os dados se comportam assim?

§3De onde vem a sobredispersão: heterogeneidade não observada

A Poisson surge naturalmente quando eventos raros acontecem independentemente, todos com a mesma taxa. É essa última condição que costuma falhar. Pense em consultas médicas por ano: duas pessoas com a mesma idade, o mesmo sexo e a mesma renda (as covariáveis do modelo) não têm a mesma propensão a ir ao médico. Uma tem uma doença crônica, a outra mora longe do posto, uma terceira é mais ansiosa. Nada disso está nos dados. O modelo atribui às duas a mesma \(\mu_i\), mas as taxas verdadeiras são diferentes.

A forma de formalizar isso é dar a cada unidade uma taxa aleatória. Chame de \(Z\) a taxa que a natureza sorteia para uma unidade. O modelo passa a ter dois andares:

Condicionalmente a \(Z\), tudo continua sendo Poisson: a variância é igual à média. Mas nós não observamos \(Z\); observamos só \(Y\). A distribuição que os dados seguem é a marginal de \(Y\) — a Poisson “misturada” sobre todas as taxas possíveis. Toda a pergunta se reduz a: quais são a média, a variância e a distribuição dessa marginal?

Sobre a Gama, uma convenção: usamos forma \(\alpha\) e taxa \(\beta\), com densidade

(Alguns livros e o numpy usam a escala \(1/\beta\) em vez da taxa. As contas são as mesmas, trocando \(\beta\) por \(1/\text{escala}\).) Por que Gama e não outra distribuição? Por enquanto, aceite que ela é flexível (assume valores positivos, com assimetria ajustável) — a §8 mostra a razão mais profunda.

§4Ferramenta 1: a esperança iterada, \(\mathbb E[\mathbb E[Y \mid Z]]\)

A primeira ferramenta é a lei da esperança total (ou esperança iterada):

Leia de dentro para fora. \(\mathbb E[Y \mid Z]\) é a média de \(Y\) se soubéssemos a taxa \(Z\). Como \(Z\) é aleatória, \(\mathbb E[Y \mid Z]\) é uma variável aleatória — uma função de \(Z\). A esperança de fora tira a média dessa função sobre todos os valores possíveis de \(Z\). Em palavras: a média geral é a média das médias condicionais.

Demonstração (para \(Y\) discreta e \(Z\) contínua, nosso caso). Escreva \(p(y \mid z) = \Pr(Y = y \mid Z = z)\). Então \(\mathbb E[Y \mid Z = z] = \sum_y y\, p(y \mid z)\) e

Na segunda linha trocamos a soma com a integral (pode, porque todos os termos são não negativos — teorema de Tonelli). O passo seguinte usa a lei da probabilidade total: \(\Pr(Y = y) = \int p(y \mid z) f_Z(z)\,dz\). Guarde essa integral — ela vai reaparecer na §7 para dar a distribuição inteira de \(Y\).

§5Ferramenta 2: a lei da variância total

A média sozinha não revela sobredispersão; precisamos da variância. A lei da variância total diz:

Demonstração. Parta de \(\operatorname{Var}(Y) = \mathbb E[Y^2] - (\mathbb E[Y])^2\) e aplique (7) aos dois termos — a esperança iterada vale para qualquer função de \(Y\), inclusive \(Y^2\):

Na segunda linha usamos, condicionalmente a \(Z\), a identidade \(\mathbb E[Y^2 \mid Z] = \operatorname{Var}(Y \mid Z) + (\mathbb E[Y \mid Z])^2\). Na terceira, a chave é reconhecer que a expressão entre chaves é a variância da variável aleatória \(\mathbb E[Y \mid Z]\).

A interpretação é o coração do argumento. A variância de \(Y\) tem duas fontes: o ruído dentro de cada unidade (mesmo sabendo a taxa, a contagem varia — é o ruído Poisson) e a diferença entre as unidades (as taxas variam de uma unidade para outra). O MLG de Poisson só enxerga a primeira.

§6Aplicando as duas leis: a sobredispersão aparece sozinha

No modelo (5), dado \(Z\), \(Y\) é Poisson com média \(Z\). Então

A média. Pela esperança iterada (7):

A variância. Pela variância total (9):

Compare com a média:

Esse é o resultado principal do artigo, e vale a pena lê-lo devagar. O excesso de variância é exatamente a variância das taxas. Se todas as unidades tivessem a mesma taxa (\(\operatorname{Var}(Z) = 0\)), voltaríamos à equidispersão. Qualquer heterogeneidade não observada, por menor que seja, empurra a variância para cima da média. Repare também que (14) não usou nada específico da Gama além de \(\operatorname{Var}(Z)\): toda mistura de Poissons é sobredispersa. A Gama só vai entrar para darmos nome à distribuição resultante.

Dá para conferir isso sem acreditar em nenhuma conta: basta simular o experimento em duas etapas. Com \(\alpha = 2\) e \(\beta = 0{,}5\), as fórmulas preveem \(\mathbb E[Y] = 4\) e \(\operatorname{Var}(Y) = 4 + 8 = 12\):

import numpy as np
rng = np.random.default_rng(20261109)
alpha, beta = 2.0, 0.5
z = rng.gamma(shape=alpha, scale=1/beta, size=1_000_000)   # 1ª etapa: uma taxa por unidade
y = rng.poisson(z)                                          # 2ª etapa: Poisson dada a taxa
print(y.mean(), y.var())                                    # ≈ 4 e ≈ 12

\(\alpha = 2\), \(\beta = 0{,}5\) (taxa), \(p = \beta/(1+\beta) = 1/3\). A variância de \(Y\) é a da Poisson (4) mais a variância de \(Z\) (8): a heterogeneidade triplica a variância sem mexer na média.

§7A distribuição marginal: a integral que dá a binomial negativa

Média e variância já explicam a sobredispersão, mas para ajustar um modelo por máxima verossimilhança precisamos da distribuição de \(Y\). Ela vem da lei da probabilidade total — a mesma integral que apareceu dentro da demonstração (8):

Passo 1 — tire da integral o que não depende de \(z\) e junte as potências:

Passo 2 — reconheça o núcleo de uma Gama. O integrando é, a menos da constante, a densidade de uma \(\text{Gama}(y + \alpha,\; \beta + 1)\). Como toda densidade integra 1,

Passo 3 — substitua e reorganize:

Passo 4 — dê nome. Defina \(p = \beta/(\beta + 1)\), de modo que \(1 - p = 1/(\beta+1)\), e \(r = \alpha\). Então

Essa é a binomial negativa \(\text{BN}(r, p)\). Quando \(r\) é inteiro, \(\Gamma(y + r)/(\Gamma(r)\,y!) = \binom{y + r - 1}{y}\) e (19) é a distribuição clássica do número de fracassos antes do \(r\)-ésimo sucesso em ensaios de Bernoulli com probabilidade \(p\). Aqui \(r = \alpha\) pode ser qualquer real positivo; a fórmula com funções Gama é a extensão natural. O nome “negativa” vem da série binomial com expoente negativo, \((1 - q)^{-r} = \sum_y \binom{y + r - 1}{y} q^y\), que mostra de uma vez que (19) soma 1.

Conferindo os momentos com as fórmulas da BN, \(\mathbb E[Y] = r(1-p)/p\) e \(\operatorname{Var}(Y) = r(1-p)/p^2\):

exatamente (12) e (13). Os dois caminhos — momentos pela esperança iterada e distribuição pela integral — chegam ao mesmo lugar. E a simulação concorda, ponto a ponto:

Um atalho: a função geradora de probabilidades

A esperança iterada dá também um caminho sem integrar densidades. A função geradora de probabilidades de \(Y\) é \(G_Y(s) = \mathbb E[s^{Y}]\). Condicionando em \(Z\) e usando a da Poisson, \(\mathbb E[s^{Y} \mid Z] = e^{Z(s-1)}\):

em que \(M_Z(t) = (\beta/(\beta - t))^{\alpha}\) é a geradora de momentos da Gama. O último termo é a geradora da \(\text{BN}(\alpha, p)\). Mesmo resultado, agora com \(\mathbb E[\mathbb E[\,\cdot \mid Z]]\) fazendo todo o trabalho.

§8Por que justamente a Gama?

A §6 mostrou que qualquer distribuição para \(Z\) gera sobredispersão. A Gama é escolhida porque ela é a conjugada da Poisson: o produto “Poisson × Gama” dentro da integral (15) é outra Gama, e a integral tem forma fechada. Com uma lognormal para \(Z\), por exemplo, o modelo Poisson-lognormal é perfeitamente razoável, mas \(\Pr(Y = y)\) só sai por integração numérica.

A conjugação tem um bônus. Pela regra de Bayes, a distribuição da taxa de uma unidade depois de observar sua contagem é

com média

A melhor estimativa da taxa de uma unidade é uma média ponderada entre o que ela mostrou (\(y\)) e a média da população (\(\alpha/\beta\)). É o princípio da credibilidade em seguros: um segurado com muitos sinistros paga mais, mas não proporcionalmente ao seu histórico curto.

§9Do modelo de mistura ao MLG binomial negativo

Até aqui havia uma única população. Num MLG cada observação tem sua média \(\mu_i = \exp(\mathbf x_i^\top \boldsymbol\beta)\). Para acomodar isso, escreva a taxa da unidade \(i\) como a média do modelo vezes um efeito multiplicativo não observado:

Forma e taxa iguais a \(\theta\) fazem \(\mathbb E[\varepsilon_i] = 1\). Isso não é arbitrário: se \(\varepsilon_i\) tivesse média \(c \ne 1\), \(c\) seria absorvido pelo intercepto \(\beta_0\) e não daria para separá-los. Fixar a média em 1 garante que \(\mu_i\) continua sendo a média de \(Y_i\). Note que \(Z_i = \mu_i \varepsilon_i \sim \text{Gama}(\theta, \theta/\mu_i)\): é o modelo (5) com \(\alpha = \theta\) e \(\beta = \theta/\mu_i\). Substituindo em (12), (13) e (19):

A conta da variância, de novo pela lei da variância total: \(\mathbb E[\mu_i\varepsilon_i] + \operatorname{Var}(\mu_i\varepsilon_i) = \mu_i + \mu_i^2/\theta\). O termo \(\mu_i^2/\theta\) é o excesso; ele cresce com o quadrado da média, por isso a curva cheia da figura da §2 se afasta da reta. Essa versão se chama NB2 (variância quadrática na média) e é a padrão em quase todo software.

O parâmetro \(\theta\) controla a heterogeneidade. Muitos textos usam \(\kappa = 1/\theta\) (o alpha do statsmodels; o R reporta theta). \(\theta\) pequeno = taxas muito diferentes entre unidades = muita sobredispersão. E quando \(\theta \to \infty\), \(\operatorname{Var}(\varepsilon_i) \to 0\), todas as unidades têm a mesma taxa e deveríamos recuperar a Poisson. De fato, reescrevendo (25):

A Poisson é um caso-limite da binomial negativa. Essa relação de “modelos encaixados” vai nos dar um teste na §11.

§10A binomial negativa é um MLG?

Um MLG exige que a distribuição pertença à família exponencial. Com \(\theta\) fixo, a BN pertence. Tome o logaritmo de (25):

Com o parâmetro canônico \(\vartheta = \log(\mu/(\mu+\theta))\), temos \(1 - e^{\vartheta} = \theta/(\mu + \theta)\), e a forma é \(y\vartheta - b(\vartheta) + c(y)\) com \(b(\vartheta) = -\theta \log(1 - e^{\vartheta})\). As derivadas de \(b\) devolvem média e variância, como deve ser numa família exponencial:

A ligação canônica seria \(\log(\mu/(\mu+\theta))\), que depende de \(\theta\) e é pouco interpretável; na prática usa-se a ligação log, \(\log\mu_i = \mathbf x_i^\top\boldsymbol\beta\), a mesma da Poisson. Assim os coeficientes têm a mesma leitura nos dois modelos: \(e^{\beta_j}\) é a razão de médias para um aumento de uma unidade em \(x_j\).

Estimação

Para \(\theta\) fixo, \(\boldsymbol\beta\) sai pelo algoritmo usual dos MLGs, mínimos quadrados reponderados iterativamente (IRLS). Com ligação log, \(d\mu/d\eta = \mu\), e os pesos e as equações de escore ficam

Compare com a Poisson, cujo escore é \(\sum_i (y_i - \mu_i)\,\mathbf x_i = \mathbf 0\) e cujos pesos são \(w_i = \mu_i\). A BN dá menos peso às observações de média alta, justamente as mais ruidosas. Para \(\theta\), maximiza-se a log-verossimilhança com \(\boldsymbol\beta\) fixo e alterna-se entre os dois passos até convergir (é o que faz o MASS::glm.nb do R). Como \(\mathbb E[\partial^2\ell/\partial\boldsymbol\beta\,\partial\theta] = 0\), as estimativas de \(\boldsymbol\beta\) e \(\theta\) são assintoticamente independentes, e a matriz de covariância de \(\hat{\boldsymbol\beta}\) é \((\mathbf X^\top \mathbf W \mathbf X)^{-1}\) com os pesos (29).

O que exatamente a Poisson erra

A equação de escore da Poisson tem esperança zero sempre que a média estiver certa, \(\mathbb E[Y_i] = \mu_i\) — não importa a variância. Por isso \(\hat{\boldsymbol\beta}_{\text{Poisson}}\) continua consistente sob sobredispersão. O erro está na incerteza. O MLG de Poisson reporta \((\mathbf X^\top \operatorname{diag}(\mu_i)\,\mathbf X)^{-1}\), que pressupõe \(\operatorname{Var}(Y_i) = \mu_i\). A covariância verdadeira do estimador é o “sanduíche”

com \(\mathbf M = \operatorname{diag}(\mu_i)\), maior que a reportada. Há três caminhos para corrigir:

1. Quasi-Poisson: supõe \(\operatorname{Var}(Y_i) = \phi\,\mu_i\) e multiplica os erros-padrão por \(\sqrt{\hat\phi}\). Mantém as estimativas da Poisson. Simples, mas a forma da variância (proporcional à média) pode estar errada, e não há verossimilhança (sem AIC, sem probabilidades previstas).
2. Erros-padrão robustos (sanduíche): estimam (30) sem supor forma para a variância. Bons com \(n\) grande.
3. Binomial negativa: modela a variância explicitamente com (25). Tem verossimilhança completa, dá probabilidades para cada contagem e é mais eficiente quando o modelo está certo.

§11Um exemplo completo

Simulamos \(n = 500\) contagens exatamente pelo caminho de (24): uma covariável \(x_i \sim \mathcal N(0, 1)\), \(\log\mu_i = 1 + 0{,}5\,x_i\), um efeito \(\varepsilon_i \sim \text{Gama}(1{,}5;\ 1{,}5)\) para cada unidade e \(Y_i \mid \varepsilon_i \sim \text{Poisson}(\mu_i\varepsilon_i)\). Como a simulação é nossa, sabemos a resposta: \(\beta_0 = 1\), \(\beta_1 = 0{,}5\), \(\theta = 1{,}5\). Depois ajustamos Poisson, quasi-Poisson e binomial negativa:

n = 500 contagens simuladas com \(\log\mu_i = 1 + 0{,}5\,x_i\) e \(\theta = 1{,}5\). Entre parênteses, erros-padrão. Média amostral de \(y\): 2,89; variância: 11,49; zeros: 124. A quasi-Poisson tem as mesmas estimativas da Poisson, só multiplica os erros-padrão por \(\sqrt{\hat\phi}\).

Leia a tabela por linhas:

• As estimativas pontuais quase coincidem (\(\hat\beta_1 = 0{,}463\) na Poisson, \(0{,}478\) na BN). Era o esperado: as duas equações de escore são não viesadas quando a média está certa.
• Os erros-padrão não. O da Poisson para \(\beta_1\) é \(0{,}028\); o da BN, \(0{,}053\) — quase o dobro. O intervalo de 95% da Poisson, \([0{,}407;\ 0{,}518]\), contém o valor verdadeiro 0,5, mas tem pouco mais da metade da largura do intervalo da BN, \([0{,}375;\ 0{,}581]\): promete uma precisão que os dados não têm. Nesta amostra deu certo; a §12 mostra quantas vezes falha.
• O diagnóstico de dispersão vai de 3,20 na Poisson para 1,02 na BN: depois de modelar a heterogeneidade, os resíduos de Pearson voltam a ter variância 1.
• \(\hat\theta = 1{,}29\) (erro-padrão 0,13), compatível com o verdadeiro 1,5.
• AIC cai de 2594,9 para 2146,9 — uma diferença enorme por um único parâmetro extra.

Teste da razão de verossimilhanças: \(H_0\colon\) Poisson

Como a Poisson é o limite \(\theta \to \infty\) da BN (26), podemos testar \(H_0\colon \kappa = 1/\theta = 0\) pela razão de verossimilhanças:

Um detalhe que muitos esquecem: \(\kappa = 0\) está na fronteira do espaço paramétrico (\(\kappa\) não pode ser negativo). Nesse caso a distribuição nula de RV não é \(\chi^2_1\), e sim uma mistura meio a meio de uma massa em zero e de \(\chi^2_1\). Na prática: divida por 2 o p-valor da \(\chi^2_1\). Aqui isso não muda nada — o p-valor fica abaixo de \(10^{-99}\) —, mas em casos limítrofes muda.

As frequências

O teste mais concreto é comparar quantas observações de cada valor o modelo previa com quantas apareceram:

A Poisson erra exatamente onde a teoria diz que erraria. A mistura espalha massa para os dois lados: unidades com taxa baixa geram zeros demais, unidades com taxa alta geram valores altos demais. A Poisson, presa a \(\operatorname{Var} = \mu\), não consegue fazer as duas coisas e concentra tudo no meio. Muito do que se chama “excesso de zeros” é só isso — sobredispersão — e a BN resolve sem precisar de um modelo inflacionado de zeros.

§12Quanto custa ignorar a sobredispersão

Um exemplo é uma amostra. Para medir o comportamento dos métodos, repetimos o experimento 2.000 vezes para cada nível de heterogeneidade: geramos dados com \(n = 200\), ajustamos os três modelos e verificamos se o intervalo de Wald de 95% para \(\beta_1\) contém o valor verdadeiro 0,5. Um bom método acerta 95% das vezes.

2.000 conjuntos de dados por linha, n = 200, \(\log\mu = 1 + 0{,}5x\). Cada célula: cobertura do IC de Wald de 95% para \(\beta_1\) · largura média. Erro de Monte Carlo da cobertura: ± 0,5 p.p. A última coluna é o desvio-padrão de \(\hat\beta_1\) entre os conjuntos: o que um bom erro-padrão deveria estimar.

Quatro conclusões:

1. Com dados Poisson (\(\theta = \infty\)), os três acertam. A BN não custa nada aqui: \(\hat\theta\) simplesmente vai para valores enormes e o modelo vira Poisson. Usar BN “por precaução” é barato.
2. A Poisson piora rapidamente com a heterogeneidade. Com \(\theta = 2\) (taxas com desvio-padrão de 71% em torno da média), o intervalo de “95%” cobre 73%; com \(\theta = 0{,}5\), menos da metade das vezes. Um teste de nível 5% rejeitaria \(H_0\) verdadeira mais de 50% das vezes.
3. A quasi-Poisson ajuda muito, mas não totalmente. Ela corrige o tamanho médio da variância, não a forma: supõe \(\operatorname{Var} \propto \mu\) quando a verdadeira é \(\mu + \mu^2/\theta\). Com \(\theta = 0{,}5\), cobre 88,6%.
4. A BN fica perto de 95% em todos os casos e, além disso, é mais eficiente: o desvio-padrão real de \(\hat\beta_1\) é menor com a BN (0,116 contra 0,132 da Poisson com \(\theta = 0{,}5\)) porque os pesos (29) dão menos importância às observações mais ruidosas.

§13Na prática: Python e R

import statsmodels.api as sm
import statsmodels.formula.api as smf

pois = smf.glm("y ~ x", data=df, family=sm.families.Poisson()).fit()
phi = pois.pearson_chi2 / pois.df_resid                 # dispersão (4): perto de 1?

quasi = smf.glm("y ~ x", data=df, family=sm.families.Poisson()).fit(scale="X2")   # quasi-Poisson
robust = smf.glm("y ~ x", data=df, family=sm.families.Poisson()).fit(cov_type="HC0")  # sanduíche (30)

nb = smf.negativebinomial("y ~ x", data=df).fit()      # NB2; reporta alpha = 1/theta
lr = 2 * (nb.llf - pois.llf)                            # p-valor: chi2(1) dividido por 2
print(nb.summary())

pois  <- glm(y ~ x, family = poisson, data = d)
sum(residuals(pois, type = "pearson")^2) / df.residual(pois)   # dispersão

quasi <- glm(y ~ x, family = quasipoisson, data = d)

library(MASS)
nb <- glm.nb(y ~ x, data = d)          # reporta theta (= 1/alpha)
lr <- 2 * (logLik(nb) - logLik(pois))
pchisq(lr, df = 1, lower.tail = FALSE) / 2   # teste na fronteira

Atenção às parametrizações: o R reporta \(\theta\); o statsmodels reporta \(\alpha = 1/\theta\); alguns textos chamam \(1/\theta\) de \(k\) ou \(\kappa\). É sempre o mesmo modelo (25).

§14Cuidados

• Primeiro acerte a média. Uma covariável importante esquecida, uma interação ausente ou uma relação não linear também produzem \(\hat\phi \gt 1\). A sobredispersão pode ser sintoma de um preditor linear errado; a BN “absorve” o problema, mas os coeficientes continuam errados.
• Excesso de zeros de verdade. Se há um processo que gera zeros estruturais (quem nunca fuma tem zero cigarros, sempre), nem a BN dá conta; modelos inflacionados de zeros ou hurdle separam os dois processos. Compare as frequências de zero, como na §11, antes de concluir.
• Dependência. Observações agrupadas (alunos na mesma escola, medidas repetidas no mesmo paciente) geram sobredispersão que não é heterogeneidade individual. Aí o caminho é um modelo misto (efeito aleatório por grupo) ou GEE.
• NB1 × NB2. Se a variância cresce linearmente com a média, \(\operatorname{Var} = \mu(1 + \delta)\), a NB1 ou a quasi-Poisson descrevem melhor. Um gráfico como o da §2 indica qual forma é mais compatível com os dados.
• Subdispersão (\(\hat\phi \lt 1\)) existe, mas é rara; a mistura nunca a produz, por (14). Quando aparece, pense em processos regulares (contagens com teto, eventos que se inibem) e em modelos como Conway–Maxwell–Poisson ou binomial.

§15Resumo

1. O MLG de Poisson impõe \(\operatorname{Var}(Y_i) = \mu_i\). Dados reais costumam ter variância maior; o diagnóstico é \(\hat\phi = X^2/(n-p) \gt 1\).
2. A causa mais comum é heterogeneidade não observada: unidades com as mesmas covariáveis têm taxas diferentes. Modelamos isso em dois andares: \(Z \sim \text{Gama}(\alpha, \beta)\) e \(Y \mid Z \sim \text{Poisson}(Z)\).
3. A esperança iterada dá \(\mathbb E[Y] = \mathbb E[\mathbb E[Y \mid Z]] = \mathbb E[Z]\). A variância total dá \(\operatorname{Var}(Y) = \mathbb E[Z] + \operatorname{Var}(Z)\). Logo \(\operatorname{Var}(Y) - \mathbb E[Y] = \operatorname{Var}(Z) \gt 0\): o excesso de variância é a variância das taxas.
4. Integrando a Poisson contra a Gama, a marginal é a binomial negativa \(\text{BN}(r = \alpha,\; p = \beta/(1+\beta))\). A Gama é escolhida por ser conjugada: a integral tem forma fechada.
5. Para o MLG, \(Y_i \mid \varepsilon_i \sim \text{Poisson}(\mu_i\varepsilon_i)\) com \(\varepsilon_i \sim \text{Gama}(\theta, \theta)\), o que dá \(\mathbb E[Y_i] = \mu_i\) e \(\operatorname{Var}(Y_i) = \mu_i + \mu_i^2/\theta\). Quando \(\theta \to \infty\), volta-se à Poisson.
6. Com \(\theta\) fixo a BN é família exponencial, então é um MLG. Usa-se ligação log, os coeficientes têm a mesma leitura da Poisson, e \(\theta\) é estimado por máxima verossimilhança.
7. Ignorar a sobredispersão não vicia as estimativas, mas encolhe os erros-padrão: na simulação, o IC de “95%” da Poisson cobriu 48% das vezes com \(\theta = 0{,}5\). A BN ficou perto de 95% em todos os cenários e não custou nada quando os dados eram Poisson.

§16Para ir além

• M. Greenwood e G. U. Yule (1920). An inquiry into the nature of frequency distributions representative of multiple happenings… Journal of the Royal Statistical Society, 83(2). A origem da mistura Poisson–Gama, com dados de acidentes.
• P. McCullagh e J. A. Nelder (1989). Generalized Linear Models, 2ª ed., cap. 6. Chapman & Hall. Quasi-verossimilhança e sobredispersão.
• J. F. Lawless (1987). Negative binomial and mixed Poisson regression. Canadian Journal of Statistics, 15(3).
• A. C. Cameron e P. K. Trivedi (2013). Regression Analysis of Count Data, 2ª ed. Cambridge. A referência completa sobre modelos de contagem.
• J. M. Hilbe (2011). Negative Binomial Regression, 2ª ed. Cambridge. Prático, com muitos exemplos e as variantes NB1, NB2 e outras.
• A. Agresti (2015). Foundations of Linear and Generalized Linear Models, cap. 7. Wiley.
• W. N. Venables e B. D. Ripley (2002). Modern Applied Statistics with S, 4ª ed., seção 7.4. Springer. A base do glm.nb.