Quasi-verossimilhança em MLG

Para ajustar um MLG você escolhe uma distribuição — Poisson, binomial, Gama — e maximiza a verossimilhança. Mas, olhando com atenção para as equações que o algoritmo resolve, a distribuição quase não aparece: só entram a média \(\mu_i\) e a função de variância \(V(\mu_i)\). Robert Wedderburn percebeu isso em 1974 e tirou a consequência: se só a média e a variância importam, dá para especificar só a média e a variância e jogar fora o resto da distribuição. O resultado é a quasi-verossimilhança, a base da quasi-Poisson, da quasi-binomial e de modelos para dados que não seguem distribuição conhecida nenhuma. Este artigo conta de onde ela veio, demonstra por que funciona, mostra os estimadores, as propriedades algébricas que a tornam tão conveniente e — com três exemplos e quatro simulações — quando ela é a ferramenta certa e quando não é. Todos os números saem de um script reprodutível.

§1O que um MLG realmente usa da distribuição

Num MLG, a resposta pertence à família exponencial:

\(\theta\) é o parâmetro canônico, \(\phi\) o de dispersão e \(V(\mu) = b''(\theta)\), escrita como função de \(\mu\), é a função de variância: \(V(\mu) = \mu\) na Poisson, \(\mu(1-\mu)\) na binomial (para proporções), \(\mu^2\) na Gama, \(1\) na normal. Agora derive a log-verossimilhança de uma observação em relação a um coeficiente \(\beta_j\), pela regra da cadeia \(\beta_j \to \eta \to \mu \to \theta\):

(O termo do meio vem de \(d\mu/d\theta = b''(\theta) = V(\mu)\).) Somando nas observações, as equações de escore que definem o estimador de máxima verossimilhança são

Leia (3) com cuidado. Não há \(b(\cdot)\), não há \(c(y,\phi)\), não há forma da densidade. Só aparecem a média \(\mu_i(\boldsymbol\beta)\) e a variância \(V(\mu_i)\). Duas distribuições com a mesma média e a mesma função de variância dariam exatamente o mesmo \(\hat{\boldsymbol\beta}\). Essa é a observação que dá origem a tudo o que segue.

§2A origem: Wedderburn e a mancha-da-folha

Os MLGs foram propostos por John Nelder e Robert Wedderburn em 1972. Dois anos depois, Wedderburn — um jovem estatístico da Rothamsted Experimental Station — publicou na Biometrika “Quasi-likelihood functions, generalized linear models, and the Gauss–Newton method”. O problema concreto era um conjunto de dados sobre a mancha-da-folha (Rhynchosporium secalis) em cevada: a porcentagem da área foliar atingida, para 10 variedades em 9 locais.

Esses dados são proporções, mas não são binomiais: não há “\(m\) ensaios” e “\(y\) sucessos”, há uma área contínua. Não existe distribuição natural para eles. Mesmo assim, dá para dizer bastante coisa sem escolher uma — primeiro sobre a média, depois sobre a variância.

A média. O experimento cruza dois fatores, variedade \(i\) e local \(j\), com uma observação por combinação. O modelo natural é

com 1 + 9 + 8 = 18 parâmetros para 90 observações. A ligação logit faz duas coisas. Garante que \(\mu_{ij}\) fique entre 0 e 1 quaisquer que sejam \(\alpha_i\) e \(\beta_j\). E dá aos efeitos uma leitura multiplicativa: \(e^{\alpha_i - \alpha_{i'}}\) é a razão entre as chances de infecção de duas variedades, a mesma em todos os locais. Uma variedade resistente reduz a chance pelo mesmo fator num local seco e num local úmido — é a hipótese de ausência de interação, e com uma observação por célula ela nem poderia ser testada.

A variância: o que é obrigatório. Uma proporção vive em \([0, 1]\), e isso já limita a variância. Como \(Y^2 \le Y\) nesse intervalo, \(\mathbb E[Y^2] \le \mathbb E[Y] = \mu\) e

Qualquer variável em \([0,1]\) obedece a essa desigualdade, binomial ou não. Consequência: a variância tem de ir a zero quando a proporção média vai a 0 ou a 1. Um local onde a doença quase não aparece não pode ter proporções muito espalhadas. Mas isso ainda deixa muitas formas possíveis, e duas são candidatas naturais: \(\mu(1-\mu)\), a forma binomial, e \(\mu^2(1-\mu)^2\).

A variância: qual forma. A diferença entre elas aparece quando se olha o coeficiente de variação, o desvio-padrão relativo à média, perto de zero:

Com a forma binomial, o coeficiente de variação cresce como \(1/\sqrt\mu\): uma variedade que perde em média 1% da folha teria flutuações relativas mais de três vezes maiores que uma que perde 10% (\(\sqrt{99}/\sqrt{9} \approx 3{,}3\)). Com \(\mu^2(1-\mu)^2\), o coeficiente de variação é aproximadamente constante para proporções pequenas: uma doença em 1% da folha flutua, em proporção, tanto quanto uma em 10%. É o comportamento esperado de uma epidemia, que cresce de forma multiplicativa. E o mesmo vale, por simetria, perto de 1, trocando a área doente pela área sadia \(1-\mu\).

Há uma segunda leitura, que liga a variância à ligação. Pelo método delta, \(\operatorname{Var}(\operatorname{logit} Y) \approx \operatorname{Var}(Y)\,\big/\,[\mu(1-\mu)]^2\), pois \(d\operatorname{logit}(\mu)/d\mu = 1/[\mu(1-\mu)]\). Então

A variância de Wedderburn é exatamente a que torna \(\operatorname{logit} Y\) homocedástico: ela traduz, na escala original, o modelo que um agrônomo ajustaria transformando os dados, \(\operatorname{logit} y_{ij} = \alpha_i + \beta_j + \varepsilon_{ij}\) com erros de variância constante. (A forma binomial faz o mesmo para a transformação arco-seno da raiz.) A vantagem de não transformar é que o modelo continua sendo para \(\mathbb E[Y]\), e não para \(\mathbb E[\operatorname{logit} Y]\), e que proporções iguais a 0 — comuns em dados de doença — não quebram nada.

Wedderburn ajustou o modelo com a variância binomial, viu que ela não descrevia bem os resíduos e propôs

Não se conhece nenhuma distribuição da família exponencial com essa função de variância. Curiosamente, a parte que depende de \(\mu\) até tem a forma certa: integrando \(d\theta/d\mu = 1/V(\mu)\) obtém-se o “parâmetro canônico” \(\theta = 2\operatorname{logit}\mu - 1/\mu + 1/(1-\mu)\), e a função \(Q\) da §3.2 se escreve como \(y\theta - b(\theta)\), com \(b = \operatorname{logit}\mu + 1/(1-\mu)\). O que falta é um termo \(c(y,\phi)\) que faça \(\exp\{\cdot\}\) integrar 1 — ou seja, uma densidade de verdade. Wedderburn não precisou dele. A pergunta do artigo era: dá para estimar \(\boldsymbol\beta\) resolvendo (3) com uma \(V\) qualquer, mesmo sem verossimilhança por trás? E o estimador resultante tem boas propriedades? A resposta, sim e sim, é o conteúdo das §3–§5. A §8 refaz o exemplo com dados simulados no mesmo formato.

A ideia foi completada por outros:

• Peter McCullagh (1983), nos Annals of Statistics, deu a teoria assintótica: consistência, normalidade e a matriz de covariância do estimador, inclusive para observações dependentes.
• Kung-Yee Liang e Scott Zeger (1986) estenderam a ideia para dados longitudinais e agrupados: as GEE (equações de estimação generalizadas) são equações de quasi-verossimilhança multivariadas, com uma matriz de correlação “de trabalho”.
• John Nelder e Daryl Pregibon (1987) propuseram a quasi-verossimilhança estendida, que permite comparar funções de variância diferentes (§4.6).
• O livro de McCullagh e Nelder (1989), cap. 9, consolidou o método — é lá que a maioria dos estatísticos o aprendeu, e é lá que estão os dados originais da mancha-da-folha.

Wedderburn morreu em 1975, um ano depois do artigo, com menos de 30 anos. A quasi-Poisson e a quasi-binomial que hoje estão em qualquer software — family = quasipoisson no R — são consequência direta dele.

§3A fundamentação matemática

3.1 O quasi-escore

Suponha só duas coisas sobre uma observação \(Y\):

com \(V\) uma função conhecida e \(\phi \gt 0\) possivelmente desconhecido. Nada sobre a forma da distribuição. Defina o quasi-escore imitando o termo de (2):

Um escore de verossimilhança verdadeiro, \(\partial\ell/\partial\mu\), tem três propriedades que sustentam toda a teoria assintótica da máxima verossimilhança. O quasi-escore tem as mesmas três, usando apenas (5).

Propriedade 1 — média zero.

Propriedade 2 — a variância. Como \(\mu\) e \(V(\mu)\) são constantes,

Propriedade 3 — a identidade da informação. Derive \(U\) em relação a \(\mu\) (regra do quociente):

O termo com \(V'(\mu)\) some porque multiplica \(Y - \mu\), que tem média zero. A igualdade “\(-\mathbb E[\text{derivada do escore}] = \operatorname{Var}(\text{escore})\)” é a identidade de Bartlett que, na verossimilhança, faz a informação de Fisher ser ao mesmo tempo a curvatura esperada e a variância do escore. É ela que dá \(\operatorname{Cov}(\hat{\boldsymbol\beta}) \approx \mathcal I^{-1}\). Como (7)–(9) valem, todos os argumentos assintóticos da máxima verossimilhança que só usam essas três propriedades valem também para o quasi-escore.

3.2 A quasi-verossimilhança

Se \(U\) se comporta como uma derivada de log-verossimilhança, integre-o:

Esta é a quasi-(log-)verossimilhança. Por construção \(\partial Q/\partial\mu = U\) e, pela regra da cadeia, \(\partial Q/\partial\beta_j = \sum_i U_i\,\partial\mu_i/\partial\beta_j\): maximizar \(Q\) é resolver (3). O limite inferior \(y\) só fixa a constante — faz \(Q(y;y) = 0\) — e não afeta nada que dependa de \(\boldsymbol\beta\). Para cada função de variância, a integral dá (com \(\phi = 1\) e a menos de termos que só dependem de \(y\)):

Quando existe uma distribuição da família exponencial com aquela \(V\), \(Q\) coincide com a sua log-verossimilhança (a menos de termos em \(y\)). A última linha é a de Wedderburn: uma “verossimilhança” sem distribuição.

Exemplo: conferindo a linha de Wedderburn. Derive cada termo de \(Q = (2y-1)\log\frac{\mu}{1-\mu} - \frac{y}{\mu} - \frac{1-y}{1-\mu}\):

Expandindo o numerador: \((2y\mu - 2y\mu^2 - \mu + \mu^2) + (y - 2y\mu + y\mu^2) - (\mu^2 - y\mu^2) = y - \mu\). Logo \(\partial Q/\partial\mu = (y-\mu)/[\mu^2(1-\mu)^2]\), como deveria. É assim que se verifica qualquer linha da tabela.

3.3 A quasi-deviance

Como na verossimilhança, a medida de ajuste é a distância ao modelo saturado (\(\mu_i = y_i\)):

Multiplicar por \(\phi\) tira o parâmetro de dispersão: \(D\) só depende de \(V\). Para \(V = \mu\) ela é a deviance da Poisson, \(2\sum[y\log(y/\mu) - (y - \mu)]\), e assim por diante. Ela serve para comparar modelos encaixados com a mesma \(V\) (§4.5) — e não para comparar funções de variância diferentes, porque cada \(V\) define uma escala diferente.

§4Os estimadores

4.1 Os coeficientes: as equações de quasi-escore

Em notação matricial, com \(\mathbf D = \partial\boldsymbol\mu/\partial\boldsymbol\beta^\top\) (matriz \(n\times p\), \(D_{ij} = \partial\mu_i/\partial\beta_j\)) e \(\mathbf V = \operatorname{diag}\{V(\mu_i)\}\), as equações (3) ficam

\(\hat{\boldsymbol\beta}\) é a raiz de (13). Repare que \(\phi\) multiplica a equação inteira e pode ser cancelado: \(\hat{\boldsymbol\beta}\) não depende de \(\phi\).

4.2 O algoritmo: escore de Fisher = IRLS = Gauss–Newton

O método de escore de Fisher troca a derivada de \(\mathbf U\) pela sua esperança, \(-\mathbb E[\partial\mathbf U/\partial\boldsymbol\beta^\top] = \mathbf D^\top\mathbf V^{-1}\mathbf D/\phi\) (a versão matricial de (9)), e itera

com tudo avaliado em \(\boldsymbol\beta^{(r)}\). Mais uma vez \(\phi\) cancela. Para ver que isso é o IRLS dos MLGs, escreva \(\mathbf D = \boldsymbol\Delta\mathbf X\), com \(\boldsymbol\Delta = \operatorname{diag}(d\mu_i/d\eta_i)\), e defina os pesos e a resposta de trabalho

Então \(\mathbf D^\top\mathbf V^{-1}\mathbf D = \mathbf X^\top\mathbf W\mathbf X\) e \(\mathbf D^\top\mathbf V^{-1}(\mathbf y-\boldsymbol\mu) = \mathbf X^\top\mathbf W\boldsymbol\Delta^{-1}(\mathbf y-\boldsymbol\mu)\), e (14) vira

uma regressão de mínimos quadrados ponderados de \(\mathbf z\) em \(\mathbf X\), repetida até convergir. É o mesmo algoritmo que ajusta qualquer MLG; a única coisa que muda de um modelo quasi para outro é a função \(V\) dentro dos pesos. O “Gauss–Newton” do título de Wedderburn vem daqui: (14) é exatamente o passo de Gauss–Newton para minimizar \(\sum_i (y_i - \mu_i)^2/V_i\) com os \(V_i\) congelados na iteração atual — mínimos quadrados não lineares ponderados.

Exemplo. No conjunto de dados de insetos da §6 (quatro inseticidas, ligação log, \(V(\mu) = \mu\)), as iterações de (16) a partir de um chute ingênuo são:

Ponto de partida: \(\mu_i^{(0)} = (y_i + \bar y)/2 + 0{,}1\). Na terceira iteração os coeficientes já estão certos até a 5ª casa; as duas últimas só confirmam a convergência. Nenhuma linha depende de \(\phi\).

4.3 A dispersão \(\phi\): o estimador de Pearson

Por (5), \(\mathbb E[(Y_i - \mu_i)^2/V(\mu_i)] = \phi\) para toda observação. Isso sugere o estimador de momentos

a estatística \(X^2\) de Pearson dividida pelos graus de liberdade. O \(n - p\) no lugar de \(n\) compensa os \(p\) parâmetros estimados, exatamente como em \(s^2 = \text{SQR}/(n-p)\) na regressão linear (que é o caso \(V = 1\)). Uma alternativa é \(D/(n-p)\), a deviance dividida pelos graus de liberdade, mas ela é viesada quando as contagens são pequenas ou as proporções estão perto de 0 e 1; \(\hat\phi\) de Pearson é consistente sempre que (5) vale, e é o padrão no R e no statsmodels.

4.4 A covariância de \(\hat{\boldsymbol\beta}\)

Expanda (13) em torno do valor verdadeiro: \(\mathbf 0 = \mathbf U(\hat{\boldsymbol\beta}) \approx \mathbf U(\boldsymbol\beta) - \mathcal I(\hat{\boldsymbol\beta} - \boldsymbol\beta)\), com \(\mathcal I = \mathbf D^\top\mathbf V^{-1}\mathbf D/\phi\). Logo \(\hat{\boldsymbol\beta} - \boldsymbol\beta \approx \mathcal I^{-1}\mathbf U(\boldsymbol\beta)\) e, como as observações são independentes,

O passo do meio usa \(\operatorname{Cov}(\mathbf U) = \mathcal I\), a versão matricial de (8)–(9). Na prática, \(\phi\) é trocado por \(\hat\phi\) de (17). Resultado: os erros-padrão do modelo quasi são os do MLG “nominal” multiplicados por \(\sqrt{\hat\phi}\). Na quasi-Poisson, os da Poisson vezes \(\sqrt{\hat\phi}\); na quasi-binomial, os da binomial vezes \(\sqrt{\hat\phi}\).

E se a função de variância estiver errada — se a variância verdadeira for \(\boldsymbol\Sigma \ne \phi\mathbf V\)? Então \(\operatorname{Cov}(\mathbf U) \ne \mathcal I\) e a covariância correta é o sanduíche

estimada trocando \(\boldsymbol\Sigma\) por \(\operatorname{diag}\{(y_i - \hat\mu_i)^2\}\). Ela não supõe forma nenhuma para a variância, mas é ruidosa com \(n\) pequeno (a §9 mede isso).

4.5 Testes: \(t\) e \(F\) no lugar de \(z\) e \(\chi^2\)

Como \(\phi\) é estimado, a inferência segue o modelo da regressão linear, não o da verossimilhança:

• Um coeficiente: \(t = \hat\beta_j/\text{e.p.}(\hat\beta_j)\), comparado com a \(t_{n-p}\).
• Modelos encaixados (o menor com \(p_0\) parâmetros, o maior com \(p_1\)):

com \(\hat\phi_1\) estimado no modelo maior. Se \(\phi\) fosse conhecido e igual a 1, \((D_0 - D_1)\) seria o teste da razão de verossimilhanças com \(\chi^2_{p_1-p_0}\). Dividir por \(\hat\phi\) “desconta” a sobredispersão, e usar \(F\) em vez de \(\chi^2\) desconta a incerteza de ter estimado \(\hat\phi\).

4.6 Quando é preciso comparar funções de variância

\(Q\) não serve para escolher entre \(V(\mu) = \mu(1-\mu)\) e \(V(\mu) = \mu^2(1-\mu)^2\): as duas estão em escalas diferentes. Nelder e Pregibon (1987) resolveram isso com a quasi-verossimilhança estendida

em que \(d_i\) é a contribuição da observação \(i\) para a deviance (12). O termo \(\log V(y_i)\) é o que permite comparar \(V\)s diferentes: \(Q^+\) é a aproximação de ponto de sela da log-verossimilhança, e para a normal e a normal inversa é exata. Maximizando em \(\phi\) dá \(\hat\phi = D/n\). Para seleção de modelos com sobredispersão em ecologia usa-se também o QAIC de Lebreton et al. (1992), \(\text{QAIC} = -2\ell/\hat c + 2k\), com \(\hat c = \hat\phi\) e \(k\) incluindo o \(\phi\).

§5As vantagens algébricas

As propriedades a seguir são a razão de a quasi-verossimilhança ser tão usada. Cada uma é uma conta curta.

5.1 Consistência exige só a média certa

Calcule a esperança de (13) no valor verdadeiro de \(\boldsymbol\beta\):

A conta não usou a variância. Uma equação de estimação com esperança zero no valor verdadeiro (uma equação de estimação não viesada) produz, sob condições de regularidade, um estimador consistente. Logo \(\hat{\boldsymbol\beta}\) é consistente mesmo que \(V\) esteja errada. A variância errada estraga os erros-padrão (que (19) conserta) e a eficiência (§5.3), mas não aponta o estimador para o lugar errado.

5.2 \(\phi\) se separa de \(\boldsymbol\beta\)

Já vimos que \(\phi\) cancela em (13) e (14). Algebricamente, isso quer dizer que o problema se decompõe em dois: primeiro estima-se \(\boldsymbol\beta\) como se \(\phi = 1\); depois, \(\phi\) por (17); por fim, a covariância (18) é a nominal multiplicada por \(\hat\phi\). Exemplo: na §6, a quasi-Poisson devolve exatamente os mesmos \(\hat\beta\) que a Poisson (na tabela, uma única coluna de estimativas serve às duas). Por isso trocar Poisson por quasi-Poisson nunca muda uma previsão; muda só a incerteza.

5.3 Otimalidade: o “Gauss–Markov” das equações de estimação

Por que justamente \(\mathbf D^\top\mathbf V^{-1}(\mathbf y-\boldsymbol\mu)\)? Considere qualquer equação de estimação linear em \(\mathbf y\),

com \(\mathbf A\) uma matriz \(n\times p\) qualquer (que pode depender de \(\boldsymbol\beta\)). Pelo mesmo argumento de (18), a raiz tem \(\hat{\boldsymbol\beta} - \boldsymbol\beta \approx (\mathbf A^\top\mathbf D)^{-1}\mathbf A^\top(\mathbf y-\boldsymbol\mu)\) e covariância assintótica

A escolha de quasi-verossimilhança é \(\mathbf A^\star = \boldsymbol\Sigma^{-1}\mathbf D\), que dá \(\mathbf C(\mathbf A^\star) = (\mathbf D^\top\boldsymbol\Sigma^{-1}\mathbf D)^{-1}\). Afirmação: \(\mathbf C(\mathbf A) - \mathbf C(\mathbf A^\star)\) é semidefinida positiva para todo \(\mathbf A\). Demonstração. Defina

\(\mathbf M\boldsymbol\Sigma\mathbf M^\top\) é semidefinida positiva (é uma covariância). Expandindo, os dois termos cruzados são iguais a

e o último termo também é \(\mathbf C(\mathbf A^\star)\). Logo

Entre todos os estimadores definidos por equações lineares em \(\mathbf y\), o de quasi-verossimilhança tem a menor variância assintótica (Godambe, 1960; McCullagh, 1983). É a generalização do teorema de Gauss–Markov: no caso \(V = 1\) e ligação identidade, \(\mathbf A^\star \propto \mathbf X\) e (23) são as equações normais dos mínimos quadrados. Duas consequências práticas: (i) acertar a forma de \(V\) compra eficiência — a §9 mede quanto; (ii) a garantia é só dentro da classe linear: se a distribuição verdadeira for conhecida e não for da família exponencial, a máxima verossimilhança completa pode ser mais eficiente.

5.4 A família exponencial é um caso particular

Se existe uma distribuição da família exponencial com função de variância \(V\), então \(Q\) é a sua log-verossimilhança (a menos de termos em \(y\)) — basta comparar (2) com (6) e integrar. Nesse caso o estimador de quasi-verossimilhança é o de máxima verossimilhança, com toda a sua eficiência. A quasi-verossimilhança não é um substituto aproximado do MLG: ela é a teoria geral, e o MLG clássico é o caso em que existe uma distribuição por trás.

5.5 Mais liberdade com o mesmo algoritmo

• Qualquer \(V\) positiva, mesmo sem distribuição (como (4)), cabe no mesmo IRLS (16).
• \(\phi\) livre nos dois sentidos: \(\hat\phi \lt 1\) (subdispersão) é tão admissível quanto \(\hat\phi \gt 1\). Na Poisson e na binomial, \(\phi = 1\) é imposto.
• Respostas fora do suporte “oficial”: a quasi-Poisson aceita \(y\) contínuo e não negativo (gastos, concentrações), e a quasi-binomial aceita proporções contínuas — \(Q\) só exige que a integral (10) exista.

§6Exemplo 1 — quasi-Poisson: insetos por parcela

Quatro inseticidas (A, B, C, D) são aplicados em 12 parcelas cada, e conta-se o número de insetos vivos por parcela. As médias verdadeiras são 15, 12, 4 e 6, mas as parcelas diferem entre si (microclima, borda, vegetação), e os dados foram gerados com \(\operatorname{Var}(Y) = 3\mu\) — o caso exato que a quasi-Poisson descreve, com \(\phi = 3\). O primeiro diagnóstico não precisa de modelo: com grupos, compare média e variância.

12 parcelas por inseticida. Os dados foram gerados com \(\operatorname{Var}(Y) = 3\mu\). Sob Poisson, a última coluna ficaria perto de 1; aqui ela fica entre 2 e 6, sem tendência clara com a média — o padrão que a quasi-Poisson supõe.

O modelo é \(\log\mu_{i} = \beta_0 + \beta_1[\text{B}] + \beta_2[\text{C}] + \beta_3[\text{D}]\). Ajustando Poisson e quasi-Poisson:

n = 48, \(\hat\phi = X^2/(n-p) = 147,1/44 = 3,34\) (verdadeiro: 3). A estimativa é a mesma nas três colunas; só muda o erro-padrão. Quasi-Poisson = Poisson × \(\sqrt{\hat\phi} = 1,83\). O intervalo quasi usa o quantil \(t_{44}\) = 2,015.

Três leituras:

• \(\hat\phi = 3{,}34\), perto do verdadeiro 3. Os erros-padrão da Poisson estão \(\sqrt{3{,}34} = 1{,}83\) vezes pequenos demais.
• O intervalo para B − A passa de \([-0{,}272;\ 0{,}170]\) para \([-0{,}467;\ 0{,}365]\). Em razão de taxas, \(e^{\hat\beta_1} = 0{,}95\): os dois inseticidas são indistinguíveis nos dois modelos, mas a Poisson promete uma precisão quase duas vezes maior do que os dados têm.
• O teste global de efeito do inseticida continua significativo — \(F_{3,44} = 7{,}18\), p = 0,0005 — mas a evidência é mais de dez ordens de grandeza menos extrema que a da Poisson (\(\chi^2_3 = 72{,}0\), p \(\approx 10^{-15}\)). Com efeitos menores, é nessa distância que moram os falsos positivos.

A coluna sanduíche fica entre as duas para alguns coeficientes e acima da quasi para C − A. Com 12 parcelas por grupo, ela estima uma variância separada por grupo e é bem mais ruidosa que a quasi-Poisson, que junta todas as 48 observações num único \(\hat\phi\).

§7Exemplo 2 — quasi-binomial: germinação de sementes

Dois tipos de semente são postos para germinar em placas de 30 sementes, com água ou com um extrato vegetal; oito placas por combinação. A pergunta é se o extrato age diferente nos dois tipos — a interação. Sementes da mesma placa compartilham umidade e temperatura, então germinam de forma correlacionada: o número de germinadas por placa não é binomial. Se a probabilidade de cada placa varia em torno da média com correlação intraplaca \(\rho\), a variância é

ou seja, a variância binomial vezes \(\phi = 1 + (m-1)\rho\). Com \(m\) igual em todas as placas, isso é exatamente a quasi-binomial. Aqui \(\rho = 0{,}06\) dá \(\phi = 2{,}74\), e a interação verdadeira é zero.

8 placas de 30 sementes por combinação; probabilidades de placa sorteadas de uma Beta com correlação intraplaca \(\rho = 0,06\), o que dá \(\phi = 1 + (m-1)\rho = 2,74\). A interação verdadeira é zero. Entre parênteses, erros-padrão.

Este conjunto de dados é instrutivo justamente por ser “azarado”. A combinação tipo 1 · extrato germinou 15% contra 31% esperados, o que produziu uma interação estimada de 1,2 onde a verdadeira é zero, e \(\hat\phi = 1{,}54\), bem abaixo do verdadeiro 2,74. A binomial rejeita com p \(\approx 3\times 10^{-5}\); a quasi-binomial também rejeita, com p = 0,002. Uma amostra só não diz qual método é calibrado. Para isso, repetimos o mesmo desenho 2.000 vezes, sempre sem interação, e contamos as rejeições a 5%:

• binomial (\(\chi^2_1\)): rejeita 22,6% das vezes — mais de quatro vezes o nível nominal;
• quasi-binomial (\(F_{1,28}\)): rejeita 4,4% — dentro do erro de Monte Carlo de 5%.

Nessas 2.000 repetições, \(\hat\phi\) teve média 2,77 (verdadeiro 2,74), com 80% dos valores entre 1,92 e 3,68: a nossa amostra estava nos 10% mais baixos. A quasi-binomial é calibrada no agregado, e a variabilidade de \(\hat\phi\) com 28 graus de liberdade é o preço de não saber \(\phi\) — o \(F\) já o leva em conta.

§8Exemplo 3 — a variância de Wedderburn: proporções contínuas

De volta a 1974. Simulamos dados no formato da mancha-da-folha: 10 variedades × 9 locais, a proporção da área foliar atingida com \(\operatorname{logit}\mu_{ij} = \alpha_i + \beta_j\) — de proporções quase nulas a quase 1 — e uma distribuição Beta cuja variância é exatamente \(\mu^2(1-\mu)^2\). Ajustamos o mesmo modelo para a média com duas funções de variância: a binomial \(\mu(1-\mu)\) (a quasi-binomial, que é o que muita gente faria) e a de Wedderburn (4).

As cinco primeiras linhas: média do resíduo padronizado ao quadrado em cada quinto das observações, ordenadas pelo valor ajustado. Se a função de variância estiver certa, todas ficam perto de 1. \(Q^+\) maior = melhor; ela pode comparar funções de variância diferentes, coisa que a quasi-deviance não pode.

• O padrão dos resíduos. Com a variância binomial, a média do resíduo ao quadrado sobe de 0,68 no quinto mais baixo para 1,55 no quarto quinto: a variância suposta é grande demais para as proporções pequenas, em relação às do meio. Com a de Wedderburn, os cinco quintos ficam entre 0,84 e 1,14. É o mesmo diagnóstico que Wedderburn fez em 1974. (O padrão é menos dramático do que a razão teórica \(\mu(1-\mu)\) sugeriria, porque a própria \(V\) errada distorce o ajuste: com pesos \(\propto\mu(1-\mu)\), as observações extremas quase não pesam e acabam mal ajustadas.)
• \(\hat\phi\) não compara modelos. 0,158 contra 1,113 não significa que o primeiro ajusta “melhor”; as duas \(V\) estão em escalas diferentes. Quem compara é \(Q^+\): 121,8 contra 97,1, a favor de \(\mu^2(1-\mu)^2\).
• Eficiência. O contraste variedade 10 − variedade 1 (verdadeiro 3,5) sai 4,00 com e.p. 0,50 na variância certa, e 4,17 com e.p. 0,68 na binomial. O sanduíche conta a mesma história (0,35 contra 0,46): mesmo com erros-padrão robustos, a variância errada desperdiça informação — é o teorema da §5.3 em ação.

§9Quando os modelos quasi são bons

Para medir em vez de opinar, quatro estudos de simulação com 2.000 conjuntos de dados cada. S1 e S3 geram contagens com \(\operatorname{Var} = 3\mu\) (a hipótese da quasi-Poisson é verdadeira); S2 gera dados contínuos positivos com \(\operatorname{Var} = 0{,}25\,\mu^2\) e os ajusta com a \(V\) errada (\(\mu\)) e a certa (\(\mu^2\)).

2000 conjuntos de dados por cenário; \(\log\mu = 1 + 0{,}5x\) (S1, S3) ou \(1 + 1{,}5x\) (S2), \(x \sim U(0, 2)\). Cada célula da primeira tabela: proporção de intervalos de 95% para \(\beta_1\) que contêm o valor verdadeiro. Erro de Monte Carlo: ± 0,5 p.p.

A quasi-verossimilhança é uma boa escolha quando:

1. A sobredispersão é proporcional à variância nominal. Em S1, a Poisson cobre 73,8% e a quasi-Poisson, 94,5%. É o caso de contagens agregadas de unidades homogêneas, de proporções com tamanho de grupo constante (28), e de qualquer situação em que comparar variância e média (como na tabela da §6) mostre uma razão aproximadamente constante.
2. O interesse está em \(\boldsymbol\beta\), não na distribuição. Efeitos, razões de taxas, razões de chances, testes de hipótese: tudo isso só precisa da média e de erros-padrão honestos.
3. Não existe distribuição razoável. Proporções contínuas, índices, contagens com dependência fraca e desconhecida — a mancha-da-folha é o exemplo-tipo.
4. A amostra é pequena e a \(V\) é plausível. Em S3, com \(n = 12\), a quasi-Poisson com quantil \(t\) cobre 95,2%, enquanto o sanduíche cai para 85,8%: um único parâmetro \(\phi\) é muito mais estável que uma variância por observação. Repare na diferença entre quantil normal (92,0%) e \(t\) (95,2%): com \(n\) pequeno, usar \(t\) e \(F\) importa.
5. Você quer robustez sem complicar. O ajuste é o do MLG nominal, e \(\hat{\boldsymbol\beta}\) continua consistente mesmo que a \(V\) esteja errada (5.1).

§10Quando os modelos quasi são ruins

1. Quando a forma da variância está errada. Em S2, usar \(V = \mu\) para dados com variância \(\propto\mu^2\) mantém \(\hat\beta_1\) sem viés, mas o desvio-padrão real do estimador sobe de 0,124 para 0,162 — a variância é 1,7 vez maior, ou seja, é como jogar fora cerca de 40% dos dados. E o erro-padrão reportado (0,148 em média) subestima o real (0,162): a cobertura cai para 93,3%. Multiplicar por \(\sqrt{\hat\phi}\) corrige o tamanho médio da variância, não a sua forma. O mesmo vale para contagens com \(\operatorname{Var} = \mu + \mu^2/\theta\): a quasi-Poisson dá peso \(\propto\mu\) justamente às observações mais ruidosas — o artigo Da Poisson à binomial negativa mostra a cobertura dela caindo para 88,6% nesse caso.
2. Quando você precisa de probabilidades. Não há distribuição, então não há \(\Pr(Y = 0)\), quantis preditivos, intervalos de predição ou simulação de novos dados. Para prever “quantas apólices terão zero sinistros”, é preciso um modelo completo (binomial negativa, beta-binomial).
3. Quando você precisa de verossimilhança. Sem \(\ell\) não há AIC, BIC, razão de verossimilhanças entre famílias diferentes, nem inferência bayesiana. \(Q^+\) e o QAIC são remendos úteis, não substitutos completos.
4. Com dados binários (0/1). Para \(Y \in \{0, 1\}\), \(\operatorname{Var}(Y) = \mu(1-\mu)\) obrigatoriamente: não existe sobredispersão em uma Bernoulli, e \(\phi\) não é identificável. Simulamos 2.000 conjuntos de \(n = 200\) respostas binárias com uma forte heterogeneidade não observada (um efeito aleatório \(\mathcal N(0;\ 1{,}5^2)\) na escala logit) e ajustamos a quasi-binomial: \(\hat\phi\) teve média 1,009, com 95% dos valores entre 0,979 e 1,046. A heterogeneidade existe, mas \(\hat\phi\) não a vê. Ela só aparece agrupando (\(m_i \gt 1\)) ou com um modelo misto.
5. Quando o problema é a média, não a variância. Uma covariável esquecida, uma interação ausente ou uma relação não linear inflam \(\hat\phi\). A quasi “absorve” o sintoma com erros-padrão maiores, mas os coeficientes continuam viesados. Acerte \(\mathbb E[Y]\) primeiro — é a única hipótese da qual a consistência depende (5.1).
6. Com dependência. Medidas repetidas no mesmo indivíduo e unidades agrupadas violam a independência que (18) usa. A extensão natural é a GEE, que é quasi-verossimilhança com uma matriz de correlação de trabalho — ou um modelo misto.
7. Com estruturas que a média e a variância não capturam. Excesso de zeros estrutural, truncamento, mistura de populações: aí a forma da distribuição é o fenômeno, e um modelo para ela (inflacionado de zeros, hurdle) é necessário.

§11Na prática: R e Python

# R — quasi-Poisson e quasi-binomial
qp <- glm(insetos ~ inseticida, family = quasipoisson, data = d)
summary(qp)$dispersion                      # phi-hat de Pearson (17)
anova(update(qp, . ~ 1), qp, test = "F")    # teste F (20)

qb <- glm(cbind(germ, 30 - germ) ~ semente * extrato, family = quasibinomial, data = g)

# variância de Wedderburn: família quasi com V e deviance próprias
wedd <- quasi(link = "logit", variance = "mu(1-mu)")
wedd$variance   <- function(mu) (mu * (1 - mu))^2
wedd$dev.resids <- function(y, mu, wt)
  wt * 2 * ((2*y - 1) * (qlogis(y) - qlogis(mu)) - 2 + y/mu + (1 - y)/(1 - mu))
fw <- glm(prop ~ local + variedade, family = wedd, data = folha)

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

# quasi-Poisson: o ajuste da Poisson com phi estimado por Pearson
qp = smf.glm("insetos ~ inseticida", data=d, family=sm.families.Poisson()).fit(scale="X2")
print(qp.scale)                       # phi-hat (17)
sand = smf.glm("insetos ~ inseticida", data=d, family=sm.families.Poisson()).fit(cov_type="HC0")  # (19)

# quasi-binomial com proporções e número de sementes como peso
qb = smf.glm("germ + falha ~ semente * extrato", data=g,          # duas colunas: sucessos e fracassos
             family=sm.families.Binomial()).fit(scale="X2")

Para uma \(V\) sem família pronta, o script deste artigo implementa o IRLS (16), o \(\hat\phi\) (17), a deviance (12), o sanduíche (19) e \(Q^+\) (21) em umas 60 linhas de numpy.

§12Resumo

1. As equações de escore de um MLG só usam a média \(\mu_i(\boldsymbol\beta)\) e a função de variância \(V(\mu_i)\). Wedderburn (1974) tirou a conclusão: especifique só \(\mathbb E[Y] = \mu\) e \(\operatorname{Var}(Y) = \phi V(\mu)\).
2. O quasi-escore \(U = (Y-\mu)/(\phi V(\mu))\) tem média zero, variância \(1/(\phi V)\) e satisfaz a identidade da informação — as três propriedades que sustentam a teoria assintótica da verossimilhança. Sua integral \(Q\) é a quasi-verossimilhança.
3. \(\hat{\boldsymbol\beta}\) resolve \(\mathbf D^\top\mathbf V^{-1}(\mathbf y-\boldsymbol\mu) = \mathbf 0\) pelo mesmo IRLS dos MLGs. \(\phi\) cancela; é estimado depois por \(X^2/(n-p)\); a covariância é \(\hat\phi(\mathbf X^\top\mathbf W\mathbf X)^{-1}\). Testes usam \(t\) e \(F\).
4. Vantagens: consistência exigindo só a média certa; separação entre \(\boldsymbol\beta\) e \(\phi\); otimalidade entre as equações de estimação lineares (um Gauss–Markov generalizado); a família exponencial como caso particular; qualquer \(V\), inclusive sem distribuição, e \(\phi \lt 1\).
5. Exemplos: quasi-Poisson (\(\hat\phi = 3{,}34\), erros-padrão 1,8 vez maiores que os da Poisson), quasi-binomial (erro tipo I de 22,6% → 4,4%) e a variância de Wedderburn (resíduos sem padrão, \(Q^+\) maior, contraste mais preciso).
6. São bons quando a variância é proporcional à nominal, o interesse está nos coeficientes, não há distribuição natural ou a amostra é pequena. São ruins quando a forma de \(V\) está errada (perda de eficiência e cobertura), quando se precisa de probabilidades ou verossimilhança, com dados binários, com dependência ou quando o problema está na média.

§13Para ir além

• R. W. M. Wedderburn (1974). Quasi-likelihood functions, generalized linear models, and the Gauss–Newton method. Biometrika, 61(3), 439–447. O artigo original, curto e legível.
• J. A. Nelder e R. W. M. Wedderburn (1972). Generalized linear models. Journal of the Royal Statistical Society A, 135(3), 370–384.
• P. McCullagh (1983). Quasi-likelihood functions. Annals of Statistics, 11(1), 59–67. A teoria assintótica.
• P. McCullagh e J. A. Nelder (1989). Generalized Linear Models, 2ª ed., caps. 9 e 10. Chapman & Hall. Inclui os dados da mancha-da-folha.
• J. A. Nelder e D. Pregibon (1987). An extended quasi-likelihood function. Biometrika, 74(2), 221–232.
• K.-Y. Liang e S. L. Zeger (1986). Longitudinal data analysis using generalized linear models. Biometrika, 73(1), 13–22. As GEE.
• V. P. Godambe (1960). An optimum property of regular maximum likelihood estimation. Annals of Mathematical Statistics, 31(4), 1208–1211. A origem da teoria de equações de estimação ótimas.
• D. A. Williams (1982). Extra-binomial variation in logistic linear models. Applied Statistics, 31(2), 144–148.
• J.-D. Lebreton, K. P. Burnham, J. Clobert e D. R. Anderson (1992). Modeling survival and testing biological hypotheses using marked animals. Ecological Monographs, 62(1), 67–118. Origem do QAIC.
• A. Agresti (2015). Foundations of Linear and Generalized Linear Models, cap. 8. Wiley.