Modelos para datos de recuento excesivamente dispersos.

Usamos datos de Long (1990) sobre el número de publicaciones producidas por Ph.D. bioquímicos para ilustrar la aplicación de modelos de Poisson, Poisson excesivamente disperso, binomial negativa y Poisson inflado cero.

Las variables en el conjunto de datos son

Variable Descripción
art artículos de arte en los últimos tres años de doctorado.
fem codificado uno para mujeres
mar codificado uno si está casado
kid5 número de niños menores de seis años
phd prestigio de 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 <- read.dta("http://www.stata-press.com/data/lf2/couart2.dta")
head(ab)
names(ab)
## [1] "art"  "fem"  "mar"  "kid5" "phd"  "ment"
r <- c(mean(ab$art), var(ab$art))
r
## [1] 1.692896 3.709742
c(mean=r[1], var=r[2], ratio=r[2]/r[1])
##     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 muy dispersos, pero, por supuesto, todavía no hemos considerado 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 <- glm(art~fem+mar+kid5+phd+ment, family=poisson, data=ab)
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  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 1817.4  on 914  degrees of freedom
## Residual deviance: 1634.4  on 909  degrees of freedom
## AIC: 3314.1
## 
## Number of Fisher Scoring iterations: 5

Vemos que el modelo obviamente no se ajusta a los datos. El valor crítico del cinco por ciento para un chi-cuadrado con 909 d.f. es

qchisq(0.95, df.residual(mp))
## [1] 980.2518
deviance(mp)
## [1] 1634.371
pr <- residuals(mp,"pearson")
head(pr)
##          1          2          3          4          5          6 
## -1.3986202 -1.1385809 -1.1510584 -1.2012566 -1.4860318 -0.9814832
sum(pr^2)
## [1] 1662.547

y la desviación y el chi-cuadrado de Pearson están ambos en el 1600.

Variación extra-Poisson.

Ahora asumimos que la varianza es proporcional en lugar de igual a la media, y estimamos el parámetro de escala \(\phi\) dividiendo la chi-cuadrado de Pearson por su d.f.

 phi <- sum(pr^2)/df.residual(mp)
phi
## [1] 1.828984
round(c(phi,sqrt(phi)),4)
## [1] 1.8290 1.3524

Vemos que la varianza es aproximadamente un \(83\%\) mayor que la media(por el 1.83 del phi). 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 usamos la familia de quasipoisson:

mq <- glm(art~fem+mar+kid5+phd+ment, family=quasipoisson, data=ab)
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  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 1.829006)
## 
##     Null deviance: 1817.4  on 914  degrees of freedom
## Residual deviance: 1634.4  on 909  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 5

Las estimaciones son exactamente las mismas que antes, pero los errores estándar son aproximadamente un \(35\%\) mayores. Podemos verificar este hecho fácilmente. Primero escribimos una función útil para extraer errores estándar y luego la usamos en nuestros ajustes:

se <- function(model) sqrt(diag(vcov(model)))
round(data.frame(p=coef(mp), q=coef(mq),
se.p=se(mp), se.q=se(mq), ratio=se(mq)/se(mp)), 4)

Un enfoque alternativo es ajustar un modelo de Poisson y utilizar el estimador robusto o sándwich de los errores estándar. Por lo general, esto da resultados muy similares al modelo de Poisson demasiado disperso.

Regresión binomial negativa.

Ahora ajustamos un modelo binomial negativo con los mismos predictores. Para hacer esto, necesitamos la función glm.nb () en el paquete MASS.

library(MASS)

mnb <- glm.nb(art~fem+mar+kid5+phd+ment, data=ab)

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  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Negative Binomial(2.2644) family taken to be 1)
## 
##     Null deviance: 1109.0  on 914  degrees of freedom
## Residual deviance: 1004.3  on 909  degrees of freedom
## AIC: 3135.9
## 
## Number of Fisher Scoring iterations: 1
## 
## 
##               Theta:  2.264 
##           Std. Err.:  0.271 
## 
##  2 x log-likelihood:  -3121.917
1/mnb$theta
## [1] 0.4416205
-2*(logLik(mp)-logLik(mnb))
## 'log Lik.' 180.196 (df=6)

Theta de R es la precisión del efecto aleatorio multiplicativo, y corresponde a \(\frac{1}{\sigma^2}\) en las notas. La estimación corresponde a una varianza estimada de 0.44 y es muy significativa.

Para probar la significancia de este parámetro, puede pensar en calcular el doble de la diferencia en las verosimilitudes logarítmicas entre este modelo y el modelo de Poisson, 180,2, y tratarlo como un chi cuadrado con un d.f. Sin embargo, las asintóticas habituales no se aplican porque la hipótesis nula se encuentra 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 un d.f. Alternativamente, tratar el estadístico como un chi cuadrado con un d.f. da una prueba conservadora. De cualquier manera, tenemos pruebas abrumadoras de dispersión excesiva.

Para probar hipótesis sobre los coeficientes de regresión, podemos usar las pruebas de Wald o las pruebas de razón de verosimilitud, que son posibles porque hemos hecho supuestos de distribución completos.

También existe la familia negative.binomial de glm , se puede utilizar siempre que se proporcionen los parámetros theta. (Esto se basa en el resultado de que el binomio negativo está en la familia glm para la varianza fija). Aquí hay una comprobación rápida:

 mnbg <- glm(art~fem+mar+kid5+phd+ment, 
 family=negative.binomial(mnb$theta), data=ab)
summary(mnbg)
## 
## Call:
## glm(formula = art ~ fem + mar + kid5 + phd + ment, family = negative.binomial(mnb$theta), 
##     data = ab)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.1678  -1.3617  -0.2806   0.4476   3.4524  
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.256149   0.140008   1.830  0.06765 .  
## femWomen    -0.216420   0.074043  -2.923  0.00355 ** 
## marMarried   0.150490   0.083686   1.798  0.07247 .  
## kid5        -0.176416   0.053836  -3.277  0.00109 ** 
## phd          0.015271   0.036568   0.418  0.67633    
## ment         0.029082   0.003277   8.876  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Negative Binomial(2.2644) family taken to be 1.03911)
## 
##     Null deviance: 1109.0  on 914  degrees of freedom
## Residual deviance: 1004.3  on 909  degrees of freedom
## AIC: 3133.9
## 
## Number of Fisher Scoring iterations: 6
all(abs(coef(mnb)==coef(mnbg)) < 1e-12)
## [1] TRUE

Sin embargo, los errores estándar diferirían porque glm.nb admite el hecho de que se estimó theta, mientras que glm no.

Heterogeneidad no observada.

R tiene una función dgamma (x, forma, tasa = 1, escala = 1 / tasa) para calcular la densidad de una distribución gamma con forma y escala dadas (o su recíproca la tasa). En particular, cuando el efecto aleatorio tiene varianza v, la densidad es dgamma (x, forma = 1 / v, escala = v). Esto se puede usar para dibujar la densidad.

#library(ggplot2)
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 cuantiles usando qgamma (p, shape, rate = 1, scale = 1 / rate, lower.tail = TRUE). En particular, cuando el efecto aleatorio tiene varianza v, los cuartiles son qgamma ((1: 3) / 4, forma e = 1 / v, escala = v).

qgamma((1:3)/4, shape = 1/v, scale = v)
## [1] 0.5114167 0.8572697 1.3347651

Los bioquímicos en el primer trimestre de la distribución de heterogeneidad no observada publican un \(49\%\) menos de artículos de lo esperado a partir de sus características observadas, mientras que los de la mediana publican un 14% menos y los del tercer trimestre publican un \(33\%\) más de lo esperado.

Comparación de estimaciones y errores estándar.

Comparemos estimaciones de parámetros y errores estándar bajo los modelos de Poisson, Poisson excesivamente disperso y binomio 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)

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.

Al observar los errores estándar, vemos que ambos enfoques de sobredispersión conducen a errores estándar estimados muy similares y que la regresión de Poisson ordinaria subestima los errores estándar.

Bondad de ajuste.

Podemos evaluar la bondad de ajuste del modelo binomial negativo utilizando la desviación

deviance(mnbg)
## [1] 1004.281

El modelo binomial negativo se ajusta mejor que el de Poisson, pero aún 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 Poisson excesivamente dispersos y binomiales negativos tienen diferentes funciones de varianza. Una forma de comprobar cuál puede ser más apropiado es crear grupos basados en el predictor lineal, calcular la media y la varianza para cada grupo y, finalmente, graficar la relación media-varianza.

A continuación, se muestran grupos basados en el predictor lineal binomial negativo, creado utilizando cortes con roturas en los percentiles 5 (5) 95 para producir veinte grupos de aproximadamente el mismo tamaño. Tenga en cuenta que predict () calcula un predictor lineal, a menos que se especifique lo contrario. Para predecir en la escala de la respuesta agregue la opción tipo = “respuesta”.

xb <- predict(mnb) 

g <- cut(xb, breaks=quantile(xb,seq(0,100,5)/100))

m <- tapply(ab$art, g, mean)

v <- tapply(ab$art, g, var)
plot(m, v, xlab="Mean", ylab="Variance", 
 main="Mean-Variance Relationship")

mtext("Articles Published by Ph.D. Biochemists",padj=-0.5)

x <- seq(0.63,3.37,0.02)

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)

El gráfico traza la media frente a la varianza y superpone las curvas correspondientes al modelo de Poisson demasiado disperso, donde la varianza es \(\phi\mu\), y el modelo binomial negativo, donde la varianza es \(\mu(1+\mu\sigma^2)\).

La función de varianza de Poisson hace un buen trabajo para la mayor parte de los datos, pero no capta las altas varianzas de los estudiosos más productivos. La función de varianza binomial negativa no es muy diferente pero, al ser cuadrática, puede aumentar más rápido 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 excesivamente disperso.

Modelos de Poisson inflados con cero.

Una ocurrencia frecuente con los datos de recuento es un exceso de ceros en comparación con lo que se espera en un modelo de Poisson. En realidad, esto es un problema con nuestros datos:

zobs <- ab$art == 0
head(zobs)
## [1] TRUE TRUE TRUE TRUE TRUE TRUE
 zpoi <- exp(-exp(predict(mp))) # or dpois(0,exp(predict(mp)))
head(dpois(0,exp(predict(mp))))
##         1         2         3         4         5         6 
## 0.1414034 0.2735238 0.2658201 0.2362139 0.1098883 0.3816279
 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 publicó ningún artículo en los últimos tres años de su doctorado, pero el modelo de Poisson predice que solo el \(20,9\%\) no tendría publicaciones. Claramente, el modelo subestima la probabilidad de conteos cero.

Una forma de modelar este tipo de situación es asumir que los datos provienen de una mezcla de dos poblaciones, una donde el conteo es siempre cero y otra donde el conteo tiene una distribución de Poisson con media \(\mu\)ras que los conteos positivos provienen solo de la segunda.

En el contexto de las publicaciones de Ph.D. bioquímicos, podemos imaginar que algunos tenían en mente trabajos donde las publicaciones no serían importantes, mientras que otros apuntaban a trabajos académicos donde se esperaba un récord de publicaciones. Los miembros del primer grupo publicarían cero artículos, mientras que los miembros del segundo grupo publicarían 0,1,2, …, un recuento que se puede suponer que tiene una distribución de Poisson.

La distribución del resultado puede luego modelarse en términos de dos parámetros, \(\pi\) la probabilidad de ‘siempre cero’ y \(\mu\), el número medio de publicaciones para aquellos que no están en el grupo ‘siempre cero’. Una forma natural de introducir covariables es modelar el logit de la probabilidad \(\pi\) de siempre cero y el logaritmo de la media \(\mu\) para aquellos que no están en la clase siempre cero.

Este tipo de modelo se puede ajustar en R usando la función zeroinfl () en el 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 a cero con todas las covariables en ambas ecuaciones:

library(pscl)
## Classes and Methods for R developed in the
## Political Science Computational Laboratory
## Department of Political Science
## Stanford University
## Simon Jackman
## hurdle and zeroinfl functions by Achim Zeileis
mzip <- zeroinfl(art~fem+mar+kid5+phd+ment, data=ab)

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.640839   0.121307   5.283 1.27e-07 ***
## femWomen    -0.209144   0.063405  -3.299 0.000972 ***
## marMarried   0.103750   0.071111   1.459 0.144567    
## kid5        -0.143320   0.047429  -3.022 0.002513 ** 
## phd         -0.006166   0.031008  -0.199 0.842376    
## 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.577060   0.509386  -1.133  0.25728   
## femWomen     0.109752   0.280082   0.392  0.69517   
## marMarried  -0.354018   0.317611  -1.115  0.26501   
## kid5         0.217095   0.196483   1.105  0.26920   
## phd          0.001275   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: 19 
## Log-likelihood: -1605 on 12 Df

Al observar la ecuación inflar, 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.

Al observar la ecuación para el número medio de artículos entre los que no están en la clase siempre cero, encontramos desventajas significativas para las mujeres y los científicos con niños menores de cinco años, y un gran efecto positivo del número de publicaciones del mentor, con cada artículo asociado. con 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 \(\pi\) y \(\mu\), y calculamos la probabilidad combinada de que no haya publicaciones. Hay opciones en la función de predict() llamadas “cero” y “contar” para obtenerlas. También hay una opción “prob” para calcular la densidad predicha, pero esto es excesivo ya que solo queremos la probabilidad de cero.

pr <- predict(mzip,type="zero")  # π

mu <- predict(mzip,type="count") # μ

zip <- pr + (1-pr)*exp(-mu) # also predict(mzip,type="prob")[,1]

mean(zip)
## [1] 0.2985679

También hay un modelo binomial negativo inflado con cero, que utiliza un binomio negativo para el recuento en la clase “no siempre cero”. Este modelo se puede ajustar usando zeroinfl () con el parámetro dist = “negbin”. Los enlaces alternativos para la ecuación inflar incluyen el probit, que se puede especificar usando link = “probit”

Comparación de modelos con AIC.

Da la casualidad de que para estos datos el binomio negativo también resuelve el problema. Aquí está la probabilidad de cero artículos en el binomio negativo. Partimos de los primeros principios, pero también se podría utilizar la densidad binomial negativa incorporada

munb <- exp(predict(mnb))

theta <- mnb$theta

znb <- (theta/(munb+theta))^theta

# also dnbinom(0, mu=munb, size=theta)
mean(znb)
## [1] 0.3035957

El modelo predice que \(30,4\%\) de los bioquímicos no publicaría ningún artículo en los últimos tres años de su doctorado, muy cerca del valor observado de \(30,0\%\)

Para elegir entre el binomio negativo y los modelos inflados cero debemos 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 = -2\log L + 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 usando la función AIc. También lo obtengo “a mano” para que podamos verificar el cálculo. Para nuestros datos

c(nb=AIC(mnb), zip=AIC(mzip))
##       nb      zip 
## 3135.917 3233.546
mzip$rank <- length(coef(mzip)) # add a rank component

aic <- function(model) -2*logLik(model) + 2*model$rank

sapply(list(mnb,mzip), aic)
## [1] 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 considerar son la distribución marginal de los recuentos previstos y observados y las funciones de varianza.

Modelos sin truncamiento y con obstáculos(Zero-Truncated and Hurdle Models)

Otros modelos que no hemos cubierto son el de Poisson truncado en cero y el binomio negativo, diseñados para datos que no incluyen ceros. Un ejemplo común es la duración de la estadía en un hospital, que es de al menos un día. Un enfoque sensato es ajustar un modelo binomial negativo o de Poisson que excluye cero y reescala las otras probabilidades para sumar uno. Se debe 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 () en el paquete VGAM, utilizando las familias pospoisson y posnegbinomial.

Un enfoque alternativo al exceso (o escasez) de ceros es utilizar un proceso de dos etapas, con un modelo logit para distinguir entre recuentos cero y positivos y luego un modelo de 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 truncado de Poisson o binomial negativo para el número de artículos de los que publican al menos uno. Estos modelos a menudo se denominan modelos de obstáculos. Se pueden ajustar en R usando los modelos logit separados y Poisson truncado en cero o modelos binomiales negativos y simplemente agregando las verosimilitudes log, o usando la función hurdle () en el paquete pscl.

Al comparar los modelos con obstáculos y los modelos inflados con cero, encuentro que la distinción entre cero y uno o más es más clara con los modelos con obstáculos, pero la interpretación de la media es más clara con los modelos inflados con cero.