Cuando se trata de modelizar recuentos (es decir, números enteros mayores o iguales a 0), solemos empezar con la regresión de Poisson. Se trata de un modelo lineal generalizado en el que se supone que una respuesta tiene una distribución de Poisson condicionada a una suma ponderada de predictores. Por ejemplo, podríamos modelar el número de conmociones cerebrales documentadas de los mariscales de campo de la NFL como una función de los golpes jugados y los años totales de experiencia de su línea ofensiva. Sin embargo, un posible inconveniente de la regresión de Poisson es que puede no describir con precisión la variabilidad de los recuentos.

Una distribución de Poisson está parametrizada por \(\lambda\), que resulta ser tanto su media como su varianza. Aunque es cómodo de recordar, no suele ser realista. Una distribución de recuentos suele tener una varianza que no es igual a su media. Cuando vemos que esto ocurre con los datos que suponemos (o esperamos) que se distribuyen en Poisson, decimos que tenemos infra o sobredispersión, dependiendo de si la varianza es menor o mayor que la media. Realizar una regresión de Poisson sobre datos de recuento que presentan este comportamiento da como resultado un modelo que no se ajusta bien.

Un enfoque que aborda este problema es la regresión binomial negativa. La distribución binomial negativa, al igual que la distribución de Poisson, describe las probabilidades de ocurrencia de números enteros mayores o iguales a 0. A diferencia de la distribución de Poisson, la varianza y la media no son equivalentes. Esto sugiere que podría servir como una aproximación útil para modelar recuentos con una variabilidad diferente a su media. La varianza de una distribución binomial negativa es una función de su media y tiene un parámetro adicional, k, llamado parámetro de dispersión. Digamos que nuestro recuento es la variable aleatoria Y de una distribución binomial negativa, entonces la varianza de Y es

$$ var(Y) = \mu + \mu^{2}/k $$

A medida que el parámetro de dispersión se hace más grande, la varianza converge al mismo valor que la media, y la binomial negativa se convierte en una distribución de Poisson.

Para ilustrar la distribución binomial negativa, vamos a trabajar con algunos datos del libro, Categorical Data Analysis, de Alan Agresti (2002). Los datos se presentan en la tabla 13.6 de la sección 13.4.3. Los datos proceden de una encuesta realizada a 1308 personas en la que se les preguntó cuántas víctimas de homicidio conocen. Las variables son resp, el número de víctimas que el encuestado conoce, y race, la raza del encuestado (blanco o negro). ¿Ayuda la raza a explicar el número de víctimas de homicidio que conoce una persona? Primero hay que introducir los datos en R:

> # Table 13.6> # Agresti, p. 561> black <- c(119,16,12,7,3,2,0)> white <- c(1070,60,14,4,0,0,1)> resp <- c(rep(0:6,times=black), rep(0:6,times=white))> race <- factor(c(rep("black", sum(black)), rep("white", sum(white))),+ levels = c("white","black"))> victim <- data.frame(resp, race)

Antes de pasar a la modelización, exploremos los datos. En primer lugar, observamos que la mayoría de los encuestados son blancos:

> table(race)racewhite black 1149 159 

Los negros tienen una media más alta que los blancos:

> with(victim, tapply(resp, race, mean)) white black 0.09225413 0.52201258 

Para cada raza la varianza de la muestra es aproximadamente el doble de la media. Parece que tenemos sobredispersión.

> with(victim, tapply(resp, race, var)) white black 0.1552448 1.1498288 

Por último, observamos la distribución de los recuentos por raza.

> table(resp, race) raceresp white black 0 1070 119 1 60 16 2 14 12 3 4 7 4 0 3 5 0 2 6 1 0

Pasamos al ajuste del modelo. Primero probamos la regresión de Poisson utilizando la función glm() y mostramos una parte del resumen de salida.

> # Poisson model> pGLM <- glm(resp ~ race, data=victim, family = poisson)> summary(pGLM)Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -2.38321 0.09713 -24.54 <2e-16 ***raceblack 1.73314 0.14657 11.82 <2e-16 ***

La raza es muy significativa. Parece que los negros son mucho más propensos a conocer a alguien que fue víctima de un homicidio. Pero, ¿qué significa el coeficiente 1,73? En este modelo simple con un predictor dicotómico, es la diferencia en los recuentos log esperados. Si exponenciamos el coeficiente obtenemos una relación de medias muestrales:

> exp(coef(pGLM))raceblack 5.658419 > # same thing> mean(victim$resp)/mean(victim$resp) 5.658419

De hecho, si hacemos una predicción con este modelo y exponenciamos los resultados, obtenemos las medias muestrales:

> exp(predict(pGLM, newdata = data.frame(race=c("white","black")))) 1 2 0.09225413 0.52201258 > # same thing> with(victim, tapply(resp, race, mean)) white black 0.09225413 0.52201258 

Esto dice que el recuento de víctimas conocidas para los blancos se distribuye como una Poisson con media y varianza iguales a 0.09, mientras que el recuento de víctimas conocidas para los negros se distribuye como una Poisson con media y varianza iguales a 0,52. Ya sabemos por nuestro análisis exploratorio que las varianzas observadas eran mucho mayores, por lo que no deberíamos estar muy satisfechos con las varianzas estimadas del modelo. Si examinamos los recuentos ajustados, veremos aún más pruebas de la falta de ajuste:

> # fitted counts for Poisson GLM:> fmeans <- exp(predict(pGLM, newdata = data.frame(race = c("white","black"))))> fmeans 1 2 0.09225413 0.52201258 > fittedW <- dpois(0:6,lambda = fmeans) * sum(victim$race=="white") > fittedB <- dpois(0:6,lambda = fmeans) * sum(victim$race=="black") > data.frame(Response=0:6,BlackObs=black, BlackFit=round(fittedB,1), + WhiteObs=white, WhiteFit=round(fittedW,1)) Response BlackObs BlackFit WhiteObs WhiteFit1 0 119 94.3 1070 1047.72 1 16 49.2 60 96.73 2 12 12.9 14 4.54 3 7 2.2 4 0.15 4 3 0.3 0 0.06 5 2 0.0 0 0.07 6 0 0.0 1 0.0

Antes guardamos las medias predichas en un objeto llamado fmeans. A continuación, generamos recuentos ajustados utilizando la función dpois junto con las medias estimadas para predecir la probabilidad de obtener de 0 a 6. A continuación, multiplicamos esas probabilidades por el número de encuestados para obtener recuentos ajustados. Por último, combinamos todo en un marco de datos para comparar fácilmente los valores observados y los ajustados. Dos de las cosas más dramáticas a tener en cuenta es que estamos infraajustando los recuentos de 0 y sobreajustando los recuentos de 1.

Podemos utilizar un rootogram para visualizar el ajuste de un modelo de regresión de conteo. La función rootogram() del paquete countreg lo facilita.

> countreg::rootogram(pGLM)

La línea curva roja es el ajuste teórico de Poisson. “De cada punto de la línea roja cuelga una barra, cuya altura representa la diferencia entre los recuentos esperados y los observados. Una barra que cuelga por debajo de 0 indica un ajuste insuficiente. Una barra que cuelga por encima de 0 indica un exceso de ajuste. Los recuentos se han transformado con una transformación de raíz cuadrada para evitar que los recuentos más pequeños queden ocultos y abrumados por los recuentos más grandes. Vemos una gran cantidad de infraajuste para los recuentos 2 y superiores y un sobreajuste masivo para el recuento 1.

Ahora vamos a intentar ajustar un modelo binomial negativo. Observamos que la variabilidad de los recuentos es mayor para ambas razas. Parece que la distribución binomial negativa se aproximaría mejor a la distribución de los recuentos.

Para ajustar un modelo binomial negativo en R recurrimos a la función glm.nb() del paquete MASS (un paquete que viene instalado con R). De nuevo sólo mostramos parte de la salida resumida:

> # Fit a negative binomial model> library(MASS)> nbGLM <- glm.nb(resp ~ race, data=victim)> summary(nbGLM)Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -2.3832 0.1172 -20.335 < 2e-16 ***raceblack 1.7331 0.2385 7.268 3.66e-13 *** Theta: 0.2023 

Nota primero que los coeficientes son los mismos que antes. Una vez más podemos exponer el coeficiente de la carrera para obtener una relación de medias muestrales y hacer predicciones para obtener las medias muestrales originales.

> # fitted counts for Negative Binomial GLM:> fmeans <- exp(predict(pGLM, newdata = data.frame(race = c("white","black"))))> fmeans # same as pGLM 1 2 0.09225413 0.52201258 

Pero observe que el error estándar del coeficiente racial es mayor, lo que indica más incertidumbre en nuestra estimación (0,24 frente a 0,15). Esto tiene sentido dada la variabilidad observada en nuestros recuentos. Obsérvese también la estimación de Theta. Es nuestro parámetro de dispersión. Podemos acceder a él directamente desde nuestro objeto modelo de la siguiente manera:

> nbGLM$theta 0.2023119

Y podemos utilizarlo para obtener las varianzas estimadas para los recuentos:

> fmeans + fmeans^2 * (1/nbGLM$theta) 1 2 0.134322 1.868928 

Estas están mucho más cerca de las varianzas observadas que las dadas por el modelo de Poisson.

Una vez más, visualizamos el ajuste utilizando un rootograma:

> countreg::rootogram(nbGLM)

Este parece mucho mejor que el rootograma del modelo de Poisson. Hay un ligero infraajuste/sobreajuste para los recuentos 1 a 3, pero por lo demás se ve bastante bien.

Para obtener más información sobre nuestro modelo binomial negativo, vamos a utilizar sus parámetros para simular datos y comparar los datos simulados con los datos observados. A continuación cargamos el paquete magrittr para acceder al operador %>% que nos permite “encadenar” funciones.

Primero tabulamos los recuentos y creamos un gráfico de barras para los participantes blancos y negros, respectivamente. A continuación, utilizamos los parámetros del modelo para simular los datos de una distribución binomial negativa. Los dos parámetros son mu y tamaño (es decir, parámetro de dispersión). Observe que utilizamos la función coef() para extraer los coeficientes adecuados para cada raza. Para los blancos es sólo el intercepto, para los negros es el intercepto y la pendiente (por lo que los sumamos). Luego exponenciamos para convertir de la escala logarítmica a la escala original. Simulamos el mismo número de observaciones que tenemos en nuestros datos originales.

library(magrittr) # for %>% operatorop <- par(mfrow=c(1,2))set.seed(1)victim$resp %>% `)) %>% table() %>% barplot(main = "Simulated White")victim$resp %>% `[`(victim$race=="black") %>% table() %>% barplot(main = "Observed Black")rnbinom(n = 159, size = nbGLM$theta, mu = exp(sum(coef(nbGLM)))) %>% table() %>% barplot(main = "Simulated Black")par(op)

Los datos simulados son muy similares a los datos observados, lo que nos da de nuevo confianza para elegir la regresión binomial negativa para modelar estos datos. Estos gráficos también demuestran la naturaleza condicional de nuestro modelo. La distribución binomial negativa de los recuentos depende, o está condicionada, por la raza. Cada raza tiene una media diferente pero un parámetro de dispersión común.

• Agresti, Alan (2002), Categorical Data Analysis, Wiley.
• Kleiber, Christian & Zeileis, Achim (2016): Visualizing Count Data Regressions Using Rootograms, The American Statistician, DOI: 10.1080/00031305.2016.1173590

Para preguntas o aclaraciones sobre este artículo, póngase en contacto con el StatLab de la Biblioteca de la UVA: [email protected]

Vea toda la colección de artículos del StatLab de la Biblioteca de la UVA.

Clay Ford
Consultor de Investigación Estadística
Biblioteca de la Universidad de Virginia
.