¿Por qué el cuasi-Poisson en GLM no se trata como un caso especial de binomio negativo?

Estoy tratando de ajustar modelos lineales generalizados a algunos conjuntos de datos de conteo que pueden o no estar dispersos. Las dos distribuciones canónicas que se aplican aquí son Poisson y Binomial Negativo (Negbin), con EV y varianza.$\mu$

$Var_P = \mu$

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

que pueden ajustarse en R usando glm(..,family=poisson)y glm.nb(...), respectivamente. También está la quasipoissonfamilia, que en mi opinión es un Poisson ajustado con el mismo EV y varianza.

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

es decir, caer en algún lugar entre Poisson y Negbin. El principal problema con la familia de cuasipoisson es que no hay una probabilidad correspondiente para ello y, por lo tanto, muchas pruebas estadísticas y medidas de ajuste extremadamente útiles (AIC, LR, etcétera) no están disponibles.

Si compara las variaciones QP y Negbin, podría notar que podría igualarlas poniendo . Continuando con esta lógica, podría intentar expresar la distribución de cuasipoisson como un caso especial de Negbin:$\phi = 1 + \frac{\mu}{\theta}$

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

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

(Tuve que copiar manualmente el equipada valor a partir , ya que no pude encontrar en el objeto)$\phi$summary(glmQP)glmQP

$AIC_{QP} < AIC_{NB}$$AIC_{QP}$

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

Editar: figura agregada

El cuasi-Poisson no es un modelo de máxima verosimilitud total (ML) sino un modelo cuasi-ML. Simplemente use la función de estimación (o función de puntuación) del modelo de Poisson para estimar los coeficientes, y luego emplee 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.

size$\theta_i$$\mu_i$

Si no hay regresores (solo una intersección), 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, la parametrización NB2 se usa con mayor frecuencia, pero algunos paquetes de software también ofrecen la versión NB1. Por ejemplo, en R, puede usar el gamlsspaquete para hacer gamlss(y ~ x, family = NBII). Tenga en cuenta que gamlssutiliza algo confuso NBIpara la parametrización NB2 y NBIIpara NB1. (Pero la jerga y la terminología no están unificadas en todas las comunidades).

Entonces podría preguntar, por supuesto, ¿por qué usar cuasi-Poisson si hay NB1 disponible? Todavía hay una sutil diferencia: el primero usa cuasi-ML y obtiene la estimación de la dispersión a partir de los residuos de desviación cuadrática (o Pearson). Este último usa ML completo. En la práctica, la diferencia a menudo no es grande, pero las motivaciones para usar cualquiera de los modelos son ligeramente diferentes.