Por que o quase-Poisson no GLM não é tratado como um caso especial de binômio negativo?

Estou tentando ajustar modelos lineares generalizados a alguns conjuntos de dados de contagem que podem ou não estar superdispersos. As duas distribuições canônicas que se aplicam aqui são o Poisson e o Binomial Negativo (Negbin), com EV e variância$\mu$

$Var_P = \mu$

$Var_{NB} = \mu + \frac{\mu^2}{\theta}$

que pode ser montado em R usando glm(..,family=poisson)e glm.nb(...), respectivamente. Há também a quasipoissonfamília, que no meu entender é um Poisson ajustado com o mesmo VE e variância

$Var_{QP} = \phi\mu$ ,

ou seja, caindo em algum lugar entre Poisson e Negbin. O principal problema com a família quasipoisson é que não há uma probabilidade correspondente para isso e, portanto, muitos testes estatísticos extremamente úteis e medidas de ajuste (AIC, LR etc.) não estão disponíveis.

Se você comparar as variações de QP e Negbin, poderá perceber que pode equipará-las colocando . Continuando com essa lógica, você pode tentar expressar a distribuição de quaseipoisson como um caso especial do Negbin:$\phi = 1 + \frac{\mu}{\theta}$

$QP\,(\mu,\phi) = NB\,(\mu,\theta = \frac{\mu}{\phi-1})$ ,

isto é, um Negbin com linearmente dependente de . Tentei verificar essa ideia gerando uma sequência aleatória de números de acordo com a fórmula acima e ajustando-a com :$\theta$$\mu$glm

#fix parameters

phi = 3
a = 1/50
b = 3
x = 1:100

#generating points according to an exp-linear curve
#this way the default log-link recovers the same parameters for comparison

mu = exp(a*x+b) 
y = rnbinom(n = length(mu), mu = mu, size = mu/(phi-1)) #random negbin generator

#fit a generalized linear model y = f(x)  
glmQP = glm(y~x, family=quasipoisson) #quasipoisson
glmNB = glm.nb(y~x) #negative binomial

> glmQP

Call:  glm(formula = y ~ x, family = quasipoisson)

Coefficients:
(Intercept)            x  
    3.11257      0.01854  
(Dispersion parameter for quasipoisson family taken to be 3.613573)

Degrees of Freedom: 99 Total (i.e. Null);  98 Residual
Null Deviance:      2097 
Residual Deviance: 356.8    AIC: NA

> glmNB

Call:  glm.nb(formula = y ~ x, init.theta = 23.36389741, link = log)

Coefficients:
(Intercept)            x  
    3.10182      0.01873  

Degrees of Freedom: 99 Total (i.e. Null);  98 Residual
Null Deviance:      578.1 
Residual Deviance: 107.8    AIC: 824.7

Ambos os ajustes reproduzem os parâmetros, e o quasipoisson fornece uma estimativa 'razoável' para . Agora também podemos definir um valor AIC para o quasipoisson:$\phi$

df = 3 # three model parameters: a,b, and phi
phi.fit = 3.613573 #fitted phi value copied from summary(glmQP)
mu.fit = glmQP$fitted.values 

#dnbinom = negbin density, log=T returns log probabilities
AIC = 2*df - 2*sum(dnbinom(y, mu=mu.fit, size = mu.fit/(phi.fit - 1), log=T))
> AIC
[1] 819.329

(Eu tive que copiar manualmente o valor ajustado , pois não o encontrei no objeto)$\phi$summary(glmQP)glmQP

Uma vez que , isso indicaria que o quasipoisson é, sem surpresa, o melhor ajuste; portanto, pelo menos, faz o que deve fazer e, portanto, pode ser uma definição razoável para o AIC (e, por extensão, probabilidade) de um quasipoisson. As grandes questões que me restam são então A I C Q $AIC_{QP} < AIC_{NB}$$AIC_{QP}$

1. Essa ideia faz sentido? Minha verificação é baseada em raciocínio circular?
2. A principal questão para quem "inventa" algo que parece estar faltando em um tópico bem estabelecido: se essa idéia faz sentido, por que ela já não está implementada glm?

Editar: figura adicionada

O quase-Poisson não é um modelo de máxima verossimilhança máxima (ML), mas um modelo de quase-ML. Você apenas usa a função de estimativa (ou função de pontuação) do modelo de Poisson para estimar os coeficientes e, em seguida, emprega uma certa função de variação para obter erros padrão adequados (ou melhor, uma matriz de covariância completa) para realizar inferência. Portanto, glm()não fornece e logLik()ou AIC()aqui etc.

Conforme apontado corretamente, um modelo com a mesma função de expectativa e variação pode ser incorporado na estrutura do binômio negativo (NB) se o sizeparâmetro variar junto com a expectativa . Na literatura, isso é tipicamente chamado de parametrização NB1. Veja, por exemplo, o livro de Cameron & Trivedi (Regression Analysis of Count Data) ou "Analysis of Microdata" de Winkelmann & Boes.μ $\theta_i$$\mu_i$

Se não houver regressores (apenas uma interceptação) a parametrização NB1 ea parametrização NB2 empregado por MASS's glm.nb()coincidem. Com regressores eles diferem. Na literatura estatística, a parametrização do NB2 é usada com mais frequência, mas alguns pacotes de software também oferecem a versão NB1. Por exemplo, em R, você pode usar o gamlsspacote para fazer gamlss(y ~ x, family = NBII). Observe que um pouco confuso é gamlssusado NBIpara a parametrização do NB2 e NBIIpara o NB1. (Mas o jargão e a terminologia não são unificados em todas as comunidades.)

Então você pode perguntar, é claro, por que usar quase-Poisson se há NB1 disponível? Ainda existe uma diferença sutil: o primeiro usa quase-ML e obtém a estimativa da dispersão a partir dos resíduos do desvio ao quadrado (ou Pearson). Este último usa ML completo. Na prática, a diferença geralmente não é grande, mas as motivações para o uso de qualquer modelo são ligeiramente diferentes.