¿Por qué el Quasipoisson en glm no se trata como un caso especial de Binomio Negativo?

Question

Estoy tratando de ajustar modelos lineales generalizados a algunos conjuntos de datos de conteo que podrían o no estar sobredispersos. Las dos distribuciones canónicas que se aplican aquí son la de Poisson y la del Binomio Negativo (Negbin), con E.V. $ \mu $ y la variación

$Var_P = \mu $

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

que se puede colocar en R usando glm(..,family=poisson) y glm.nb(...) respectivamente. También está el quasipoisson que, a mi entender, es un Poisson ajustado con el mismo E.V. y varianza

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

es decir, cayendo en algún lugar entre Poisson y Negbin. El principal problema de la familia de los cuasipozos es que no existe una probabilidad correspondiente para ella, y por lo tanto no se dispone de muchas pruebas estadísticas y medidas de ajuste extremadamente útiles (AIC, LR, etc.).

Si comparas las variaciones de QP y Negbin, podrías notar que podrías igualarlas poniendo $ \phi = 1 + \frac { \mu }{ \theta }$ . Siguiendo esta lógica, se podría tratar de expresar la distribución de cuasipozos como un caso especial del Negbin:

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

es decir, un Negbin con $ \theta $ linealmente dependiente de $ \mu $ . Traté de verificar esta idea generando una secuencia aleatoria de números según la fórmula anterior y ajustándola con 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 ajustes reproducen los parámetros, y el quasipoisson da una estimación 'razonable' para $ \phi $ . Ahora también podemos definir un valor AIC para el quasipoisson:

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

(Tuve que copiar manualmente el ajuste $ \phi $ valor de summary(glmQP) ya que no pude encontrarlo en el glmQP objeto)

Desde $AIC_{QP} < AIC_{NB}$ esto indicaría que el quasipoisson es, no sorprendentemente, el que mejor se ajusta; así que al menos $AIC_{QP}$ hace lo que debe hacer, y por lo tanto podría ser una definición razonable para el AIC (y por extensión, la probabilidad) de un cuasipunto. Las grandes preguntas que me quedan son entonces

1. ¿Esta idea tiene sentido? ¿Mi verificación se basa en un razonamiento circular?
2. La pregunta principal para cualquiera que "inventa" algo que parece faltar en un tema bien establecido: si esta idea tiene sentido, ¿por qué no está ya implementada en glm ?

Edición: figura añadida

Answer 1

El cuasi-Poisson no es un modelo de máxima verosimilitud (ML) completo, sino un modelo cuasi-ML. Sólo se utiliza la función de estimación (o función de puntuación) del modelo de Poisson para estimar los coeficientes, y luego se emplea una determinada función de varianza para obtener errores estándar adecuados (o más bien una matriz de covarianza completa) para realizar la inferencia. Por lo tanto, glm() no suministra y logLik() o AIC() aquí, etc.

Como bien señalas, un modelo con la misma función de expectativa y varianza puede encajarse en el marco de la binomial negativa (NB) si el size parámetro $\theta_i$ varía junto con la expectativa $\mu_i$ . En la bibliografía, esto se suele denominar parametrización NB1. Véase, por ejemplo, el libro de Cameron y Trivedi (Regression Analysis of Count Data) o "Analysis of Microdata" de Winkelmann y Boes.

Si no hay regresores (sólo un intercepto) la parametrización NB1 y la parametrización NB2 empleadas por MASS 's glm.nb() coinciden. Con los regresores difieren. En la literatura estadística se utiliza con más frecuencia la parametrización NB2, pero algunos paquetes de software también ofrecen la versión NB1. Por ejemplo, en R, se puede utilizar el gamlss paquete para hacer gamlss(y ~ x, family = NBII) . Obsérvese que, de forma algo confusa gamlss utiliza NBI para la parametrización NB2 y NBII para NB1. (Pero la jerga y la terminología no están unificadas en todas las comunidades).

Entonces se podría preguntar, por supuesto, ¿por qué usar cuasi-Poisson si hay NB1 disponible? Sigue habiendo una sutil diferencia: El primero utiliza cuasi-ML y obtiene la estimación a partir de la dispersión de los residuos de desviación al cuadrado (o de Pearson). El segundo utiliza el ML completo. En la práctica, la diferencia no suele ser grande, pero las motivaciones para utilizar uno u otro modelo son ligeramente diferentes.