GR’s Website

4.A Modelos para datos de recuento sobredispersos

Utilizamos datos de Long (1990) sobre el número de publicaciones producidas por los bioquímicos doctores para ilustrar la aplicación de los modelos de Poisson, Poisson sobredisperso, binomial negativo y Poisson inflado por cero.

Las variables del conjunto de datos son

artículos en los últimos tres años de doctorado.

codificado uno para las mujeres

codificado uno si está casado

número de hijos menores de seis años

Variable	Descripción
art
fem
mar
kid5
phd	prestigio del Ph.D. programa
ment	artículos por mentor en los últimos tres años

Estos datos también han sido analizados por Long y Freese (2001), y están disponibles en el sitio web de Stata:

> library(foreign)> ab names(ab) "art" "fem" "mar" "kid5" "phd" "ment"> r c(mean=r, var=r, ratio=r/r) mean var ratio 1.692896 3.709742 2.191358 

El número medio de artículos es 1,69 y la varianza es 3,71, un poco más del doble de la media. Los datos están sobredispersos, pero por supuesto no hemos considerado aún ninguna covariable.

Un modelo de Poisson

Ajustemos el modelo utilizado por Long y Freese(2001), un modelo aditivo simple que utiliza los cinco predictores.

> mp summary(mp)Call:glm(formula = art ~ fem + mar + kid5 + phd + ment, family = poisson, data = ab)Deviance Residuals: Min 1Q Median 3Q Max -3.5672 -1.5398 -0.3660 0.5722 5.4467 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 0.304617 0.102981 2.958 0.0031 ** femWomen -0.224594 0.054613 -4.112 3.92e-05 ***marMarried 0.155243 0.061374 2.529 0.0114 * kid5 -0.184883 0.040127 -4.607 4.08e-06 ***phd 0.012823 0.026397 0.486 0.6271 ment 0.025543 0.002006 12.733 
We see that the model obviously doesn't fit the data. The five-percent critical value for a chi-squared with 909 d.f. is
>qchisq(0.95, df.residual(mp)) 980.2518> deviance(mp) 1634.371> pr sum(pr^2) 1662.547
y la desviación y el chi-cuadrado de Pearson están ambos en los 1600.
 4.A.1. Variación extra-Poisson
Ahora suponemos que la varianza es proporcional en lugar de igual a la media, y estimamos el parámetro de escala φ dividiendo el chi-cuadrado de Pearson por su f.d.:
> phi round(c(phi,sqrt(phi)),4) 1.8290 1.3524
Vemos que la varianza es aproximadamente un 83% mayor que la media. Esto significa que debemos ajustar los errores estándar multiplicando por 1,35, la raíz cuadrada de 1,83.
R puede hacer este cálculo por nosotros si utilizamos la familia quasipoisson:
> mq summary(mq)Call:glm(formula = art ~ fem + mar + kid5 + phd + ment, family = quasipoisson, data = ab)Deviance Residuals: Min 1Q Median 3Q Max -3.5672 -1.5398 -0.3660 0.5722 5.4467 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 0.304617 0.139273 2.187 0.028983 * femWomen -0.224594 0.073860 -3.041 0.002427 ** marMarried 0.155243 0.083003 1.870 0.061759 . kid5 -0.184883 0.054268 -3.407 0.000686 ***phd 0.012823 0.035700 0.359 0.719544 ment 0.025543 0.002713 9.415 
The estimates are exactly the same as before, but the standard errorsare about 35% larger. We can verify this fact easily. First we writea useful function to extract standard errors and then use it onour fits:
> se round(data.frame(p=coef(mp), q=coef(mq),+ se.p=se(mp), se.q=se(mq), ratio=se(mq)/se(mp)), 4) p q se.p se.q ratio(Intercept) 0.3046 0.3046 0.1030 0.1393 1.3524femWomen -0.2246 -0.2246 0.0546 0.0739 1.3524marMarried 0.1552 0.1552 0.0614 0.0830 1.3524kid5 -0.1849 -0.1849 0.0401 0.0543 1.3524phd 0.0128 0.0128 0.0264 0.0357 1.3524ment 0.0255 0.0255 0.0020 0.0027 1.3524
Un enfoque alternativo es ajustar un modelo de Poisson y utilizar el estimador robusto o de sándwich de los errores estándar. Esto suele dar resultados muy similares al modelo de Poisson sobredisperso.
 4.A.2 Regresión binomial negativa
Ahora ajustamos un modelo binomial negativo con los mismos predictores.Para ello necesitamos la función glm.nb() del paqueteMASS.
> library(MASS)> mnb summary(mnb)Call:glm.nb(formula = art ~ fem + mar + kid5 + phd + ment, data = ab, init.theta = 2.264387695, link = log)Deviance Residuals: Min 1Q Median 3Q Max -2.1678 -1.3617 -0.2806 0.4476 3.4524 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 0.256144 0.137348 1.865 0.062191 . femWomen -0.216418 0.072636 -2.979 0.002887 ** marMarried 0.150489 0.082097 1.833 0.066791 . kid5 -0.176415 0.052813 -3.340 0.000837 ***phd 0.015271 0.035873 0.426 0.670326 ment 0.029082 0.003214 9.048 1/mnb$theta 0.4416205> -2*(logLik(mp)-logLik(mnb))'log Lik.' 180.196 (df=6)
R's theta es la precisión del efecto multiplicativo aleatorio, y corresponde a 1/σ2 en las notas. La estimación corresponde a una varianza estimada de 0,44 y es altamente significativa.
Para probar la significación de este parámetro se puede pensar en computar el doble de la diferencia de logaritmos de probabilidad entre este modelo y el modelo de Poisson, 180,2, y tratarlo como un chi-cuadrado con una f.d. Sin embargo, la asintótica habitual no se aplica, porque la hipótesis nula está en un límite del espacio de parámetros.
Hay algunos trabajos que muestran que una mejor aproximación es tratar el estadístico como una mezcla 50:50 de cero y un chi-cuadrado con una f.d. Alternativamente, tratar el estadístico como un chi-cuadrado con una f.d. da una prueba conservadora. De cualquier manera, tenemos una evidencia abrumadora de sobredispersión.
Para probar las hipótesis sobre los coeficientes de regresión podemos utilizar las pruebas de Wald o las pruebas de razón de verosimilitud, que son posibles porque hemos hecho supuestos de distribución completa.
También hay una familia negative.binomial para glm y esto se puede utilizar siempre que se den los parámetrostheta. (Esto se basa en el resultado de que la binomial negativa está en la familia glm para la varianza fija). He aquí una comprobación rápida:
> mnbg all(abs(coef(mnb)==coef(mnbg)) 
The standard errors would differ, however, because  allows for the fact that was estimated, whereas  does not.
Unobserved Heterogeneity
R has a function  to compute the density of a gamma distribution with given shape and scale (or its reciprocal the rate). In particular, when the random effect hasvariance  the density is. This can beused to draw the density.
> v = 1/mnb$theta> gd = data.frame(x = seq(0.0, 3, 0.1), g = dgamma(x, shape = 1/v, scale = v))> ggplot(gd, aes(x, g)) + geom_line()> ggsave("gamdenr.png", width=500/72, height=400/72, dpi=72)

 También podemos calcular los cuantiles utilizandoqgamma(p, shape, rate = 1, scale = 1/rate, lower.tail = TRUE). En particular, cuando el efecto aleatorio tiene varianza vlos cuartiles son qgamma((1:3)/4, shap e= 1/v, scale = v).
> qgamma((1:3)/4, shape = 1/v, scale = v) 0.5114167 0.8572697 1.3347651
Los bioquímicos que se encuentran en el Q1 de la distribución de la heterogeneidad no observada publican un 49% menos de trabajos de los que se esperan a partir de sus características observadas, mientras que los que se encuentran en la mediana publican un 14% menos y los que están en el Q3 publican un 33% más de lo esperado.
Comparación de estimaciones y errores estándar
Comparemos las estimaciones de los parámetros y los errores estándar bajo los modelos de Poisson, Poisson sobredispersado y binomial negativo:
> round(data.frame(+ p=coef(mp),q=coef(mq),nb=coef(mnb),+ se.p=se(mp),se.q=se(mq),se.nb=se(mnb)),4) p q nb se.p se.q se.nb(Intercept) 0.3046 0.3046 0.2561 0.1030 0.1393 0.1373femWomen -0.2246 -0.2246 -0.2164 0.0546 0.0739 0.0726marMarried 0.1552 0.1552 0.1505 0.0614 0.0830 0.0821kid5 -0.1849 -0.1849 -0.1764 0.0401 0.0543 0.0528phd 0.0128 0.0128 0.0153 0.0264 0.0357 0.0359ment 0.0255 0.0255 0.0291 0.0020 0.0027 0.0032
Las estimaciones binomiales negativas no son muy diferentes de las basadas en el modelo de Poisson, y ambos conjuntos llevarían a las mismas conclusiones. 
Examinando los errores estándar, vemos que ambos enfoques de la sobredispersión conducen a errores estándar estimados muy similares, y que la regresión de Poisson ordinaria subestima los errores estándar.
Bondad del ajuste
Podemos evaluar la bondad del ajuste del modelo binomial negativo utilizando la desviación
> deviance(mnbg) 1004.281
El modelo binomial negativo se ajusta mejor que el de Poisson, pero todavía tiene una desviación por encima del valor crítico del cinco por ciento de 980,25.
La función de varianza
Los modelos de sobredispersión de Poisson y binomial negativa tienen diferentes funciones de varianza. Una forma de comprobar cuál puede ser más apropiada es crear grupos basados en el predictor lineal, calcular la media y la varianza de cada grupo y, finalmente, trazar la relación media-varianza.
Aquí hay grupos basados en el predictor lineal binomial negativo, creados utilizando cut con cortes en los percentiles 5(5)95 para producir veinte grupos de tamaño aproximadamente igual. Tenga en cuenta que predict() calcula un predictor lineal, a menos que se especifique lo contrario. Para predecir en la escala de la respuestaañada la opción type="response". 
> xb g m v png("c4afig1.png",width=500,height=400)> plot(m, v, xlab="Mean", ylab="Variance", + main="Mean-Variance Relationship")> mtext("Articles Published by Ph.D. Biochemists",padj=-0.5)> x lines(x, x*phi, lty="dashed")> lines(x, x*(1+x/mnb$theta))> legend("topleft", lty=c("dashed","solid"), + legend=c("Q. Poisson","Neg. Binom."), inset=0.05)> dev.off();null device 1 

El gráfico representa la media frente a la varianza y superpone las curvas correspondientes al modelo de Poisson sobredisperso, donde la varianza es φμ, y al modelo binomial negativo, donde la varianza es μ(1+μσ2).
La función de varianza de Poisson hace un trabajo bastante bueno para la mayor parte de los datos, pero no logra capturar las altas varianzas de los académicos más productivos. La función de varianza binomial negativa no es demasiado diferente pero, al ser cuadrática, puede aumentar más rápidamente y hace un mejor trabajo en el extremo superior. Concluimos que el modelo binomial negativo proporciona una mejor descripción de los datos que el modelo de Poisson sobredisperso.
 4.A.3 Modelos de Poisson inflados por ceros
 Un hecho frecuente con los datos de recuento es un exceso de ceros en comparación con lo que se espera bajo un modelo de Poisson. Esto es realmente un problema con nuestros datos: 
> zobs zpoi c(obs=mean(zobs), poi=mean(zpoi)) obs poi 0.3005464 0.2092071 
Vemos que el 30,0% de los científicos de la muestra no publicaron ningún artículo en los últimos tres años de su doctorado, pero el modelo de Poisson predice que sólo el 20,9% no tendría ninguna publicación. Claramente el modelo subestima la probabilidad de recuentos nulos.
Una forma de modelar este tipo de situaciones es suponer que los datos proceden de una mezcla de dos poblaciones, una en la que el recuento es siempre nulo, y otra en la que el recuento tiene una distribución de Poisson con media μ. En este modelo, los recuentos cero pueden proceder de cualquiera de las dos poblaciones, mientras que los recuentos positivos sólo proceden de la segunda.
En el contexto de las publicaciones de los bioquímicos doctorados, podemos imaginar que algunos tenían en mente trabajos en los que las publicaciones no serían importantes, mientras que otros aspiraban a trabajos académicos en los que se esperaba un récord de publicaciones. Los miembros del primer grupo publicarían cero artículos, mientras que los del segundo grupo publicarían 0,1,2,..., un recuento que puede suponerse que tiene una distribución de Poisson.
La distribución del resultado puede entonces modelarse en términos de dos parámetros, π la probabilidad de "siempre cero", y μ, el número medio de publicaciones para los que no están en el grupo "siempre cero". Una forma natural de introducir covariables es modelar el logaritmo de la probabilidad π de siempre cero y el logaritmo de la media μ para los que no están en la clase de siempre cero.
Este tipo de modelo puede ajustarse en R utilizando la función zeroinfl() del paquete pscl. La fórmula del modelo se puede especificar como de costumbre si se van a incluir las mismas variables en ambas ecuaciones. De lo contrario, se pueden proporcionar dos conjuntos de predictores separados por una barra vertical; escriba "zeroinfl" para obtener más información.
Aquí hay un modelo inflado por cero con todas las covariables en ambas ecuaciones:
> library(pscl)> mzip summary(mzip)Call:zeroinfl(formula = art ~ fem + mar + kid5 + phd + ment, data = ab)Pearson residuals: Min 1Q Median 3Q Max -2.3253 -0.8652 -0.2826 0.5404 7.2976 Count model coefficients (poisson with log link): Estimate Std. Error z value Pr(>|z|) (Intercept) 0.640838 0.121307 5.283 1.27e-07 ***femWomen -0.209145 0.063405 -3.299 0.000972 ***marMarried 0.103751 0.071111 1.459 0.144565 kid5 -0.143320 0.047429 -3.022 0.002513 ** phd -0.006166 0.031008 -0.199 0.842378 ment 0.018098 0.002294 7.888 3.07e-15 ***Zero-inflation model coefficients (binomial with logit link): Estimate Std. Error z value Pr(>|z|) (Intercept) -0.577059 0.509386 -1.133 0.25728 femWomen 0.109746 0.280082 0.392 0.69518 marMarried -0.354014 0.317611 -1.115 0.26502 kid5 0.217097 0.196482 1.105 0.26919 phd 0.001274 0.145263 0.009 0.99300 ment -0.134114 0.045243 -2.964 0.00303 **---Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Number of iterations in BFGS optimization: 21 Log-likelihood: -1605 on 12 Df
Mirando la ecuación inflada vemos que el único predictor significativo de estar en la clase "siempre cero" es el número de artículos publicados por el mentor, con cada artículo del mentor asociado con un 12,6% menos de probabilidades de no publicar nunca.
Examinando la ecuación para el número medio o los artículos entre los que no están en la clase siempre cero, encontramos desventajas significativas para las mujeres y los científicos con hijos menores de cinco años, y un gran efecto positivo del número de publicaciones del mentor, con cada artículo asociado a un aumento del 1,8% en el número esperado de publicaciones.
Para verificar que el modelo resuelve el problema del exceso de ceros predecimos π y μ, y calculamos la probabilidad combinada de no publicar. Hay opciones en la función predict() llamadas "zero" y "count" para obtenerlas. También hay una opción "prob" para calcular la densidad predicha, pero esto es excesivo ya que sólo queremos la probabilidad de cero.
> pr mu zip mean(zip) 0.2985679
Así que el modelo resuelve el problema del exceso de ceros, prediciendo que el 29,9% de los bioquímicos no publicarán ningún artículo, mucho más cerca del valor observado del 30,0%.
También hay un modelo binomial negativo inflado a cero, que utiliza una binomial negativa para el recuento en la clase 'no siempre cero'. Este modelo se puede ajustar utilizando zeroinfl() con el parámetro dist="negbin". Los enlaces alternativos para la ecuación de inflexión incluyen el probit, que se puede especificar utilizando link="probit". 
Comparación de modelos con AIC
¡Como sucede, para estos datos la binomial negativa también resuelve el problema! Aquí está la probablidad de los artículos cero en la binomial negativa. Procedemos a partir de los primeros principios, pero también se podría utilizar la densidad binomial negativa incorporada
> munb theta znb # also dnbinom(0, mu=munb, size=theta)> mean(znb) 0.3035957
El modelo predice que el 30,4% de los bioquímicos no publicarían ningún artículo en los últimos tres años de su doctorado, muy cerca del valor observado del 30,0%.
Para elegir entre los modelos binomial negativo y cero inflado tenemos que recurrir a otros criterios. Una forma muy sencilla de comparar modelos con diferentes números de parámetros es calcular el criterio de información de Akaike (AIC), que definimos como
AIC = -2logL + 2p
donde p es el número de parámetros del modelo. El primer término es esencialmente la desviación y el segundo una penalización por el número de parámetros, y se puede calcular para la mayoría de los modelos utilizando la función AIC. También lo obtengo "a mano" para poder verificar el cálculo. Para nuestros datos
> c(nb=AIC(mnb), zip=AIC(mzip)) nb zip 3135.917 3233.546 > mzip$rank aic sapply(list(mnb,mzip), aic) 3133.917 3233.546
 Para este conjunto de datos el modelo binomial negativo es un claro ganador en términos de parsimonia y bondad de ajuste. 
Otros criterios de diagnóstico que podríamos observar son la distribución marginal de los recuentos predichos y observados y las funciones de varianza. 
Modelos truncados en cero y Hurdle
Otros modelos que no hemos cubierto son el Poisson truncado en cero y la binomial negativa, diseñados para datos que no incluyen ceros. Un ejemplo común es la duración de la estancia en un hospital, que es de al menos un día. Un enfoque sensato es ajustar un modelo Poisson o binomial negativo que excluya el cero y reescalar las otras probabilidades para que sumen uno. Hay que tener cuidado al interpretar estos modelos porque μ no es el resultado esperado, sino la media de una distribución subyacente que incluye los ceros. Estos modelos se implementan en la función vglm() del paquete VGAM, utilizando las familias pospoisson y posnegbinomial.
Un enfoque alternativo para el exceso (o la escasez) de ceros es utilizar un proceso de dos etapas, con un modelo logit para distinguir entre los recuentos cero y positivos y, a continuación, un modelo Poisson truncado en cero o binomial negativo para los recuentos positivos. En nuestro ejemplo podríamos utilizar un modelo logit para diferenciar a los que publican de los que no, y luego un modelo Poisson truncado o binomial negativo para el número de artículos de los que publican al menos uno. Estos modelos suelen denominarse modelos de obstáculos. Pueden ajustarse en R utilizando los modelos logit y poisson truncado o binomial negativo por separado y simplemente sumando las probabilidades logarítmicas, o utilizando la función hurdle() del paquete pscl.
Comparando los modelos de obstáculos y los modelos con inflexión cero, encuentro que la distinción entre cero y uno o más es más clara con los modelos de obstáculos, pero la interpretación de la media es más clara con los modelos con inflexión cero.
El siguiente artículo de JSS tiene una discusión útil de todos estos modelos: Modelos de regresión para datos de conteo en R.