library(glmmTMB)
library(performance)
library(gtsummary)
library(janitor)
library(tidyverse)Dispersión y exceso de ceros
Este material es parte de la Unidad 5 del Curso de Epidemiología - Nivel Avanzado del Instituto Nacional de Epidemiología “Dr. Juan H. Jara” - ANLIS
Dispersión y exceso de ceros by Tamara Ricardo is licensed by CC BY-NC 4.0 
Introducción
En modelos de regresión Poisson se asume equidispersión, es decir que la varianza de la variable dependiente es igual a su media: \[ var(Y) = \mu \]
Sin embargo, en la práctica, esta suposición frecuentemente no se cumple y la varianza puede ser mayor que la media (sobredispersión) o menor que la media (subdispersión). Además, en algunos casos, los datos pueden contener un exceso de ceros que el modelo Poisson no ajusta adecuadamente. Estos fenómenos pueden llevar a que las inferencias que hagamos sean incorrectas, por lo que deben controlarse al ajustar los modelos.
Sobredispersión
La sobredispersión es común en datos de conteo y puede ser causada por varios factores, como la heterogeneidad no observada, correlación entre observaciones o una alta incidencia de valores extremos. La sobredispersión puede inflar los errores estándar, lo que lleva a una subestimación de la significancia estadística. Matemáticamente, se expresa como:
\[ Var(Y) > \mu \]
Para detectar si un modelo tiene sobredispersión, se puede calcular el ratio de dispersión, que es el cociente entre la devianza del modelo y los grados de libertad. Un ratio significativamente mayor que 1 indica sobredispersión.
En R, la función check_overdispersion() del paquete performance evalúa el ratio de dispersión en un modelo. En el ejemplo de regresión Poisson el cociente del modelo final era cercano a 1, lo cual indicaba que no existía evidencia de sobredispersión y la distribución seleccionada ajustaba bien los datos.
Veamos ahora un ejemplo de datos con sobredispersión a partir de los casos de dengue reportados en los departamentos de la Costa Atlántica de la provincia de Buenos Aires (PBA)1 entre las semanas epidemiológicas 1 y 20 de 20242. Los mismos se encuentran en el archivo “dengue_costa.txt”.
Comenzaremos por cargar los paquetes necesarios para el análisis:
Ahora cargamos los datos y exploramos su estructura:
# Carga datos
datos <- read_delim("dengue_costa.txt")
# Explora datos
glimpse(datos)Rows: 1,047
Columns: 6
$ provincia <chr> "Buenos Aires", "Buenos Aires", "Buenos Aires", "Buenos…
$ departamento <chr> "BAHIA BLANCA", "BAHIA BLANCA", "BAHIA BLANCA", "BAHIA …
$ periodo <dbl> 2024, 2024, 2024, 2024, 2024, 2024, 2024, 2024, 2024, 2…
$ semana_epi <dbl> 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 4, 4, 4, 4…
$ grupo_edad_cat <chr> "0-9 años", "10-19 años", "20-34 años", "35-65 años", "…
$ casos <dbl> 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0…Las variables de interés para el modelo son:
semana_epi: semana epidemiológica de reporte de los casos.grupo_edad_cat: grupo etario.casos: número de casos reportados por departamento, semana epidemiológica y grupo etario. Es la variable dependiente.
Ajustamos un modelo Poisson con todas las variables de interés:
fit_pois <- glm(casos ~ semana_epi + grupo_edad_cat,
data = datos,
family = poisson)
# Salida del modelo
summary(fit_pois)
Call:
glm(formula = casos ~ semana_epi + grupo_edad_cat, family = poisson,
data = datos)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -3.092154 0.279663 -11.057 < 2e-16 ***
semana_epi 0.044249 0.008993 4.920 8.65e-07 ***
grupo_edad_cat10-19 años 1.172838 0.294392 3.984 6.78e-05 ***
grupo_edad_cat20-34 años 2.195273 0.271054 8.099 5.54e-16 ***
grupo_edad_cat35-65 años 2.312218 0.269545 8.578 < 2e-16 ***
grupo_edad_cat65+ años 0.730874 0.314523 2.324 0.0201 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for poisson family taken to be 1)
Null deviance: 1436.5 on 1046 degrees of freedom
Residual deviance: 1187.9 on 1041 degrees of freedom
AIC: 1735.1
Number of Fisher Scoring iterations: 6El modelo ajustado no incluye un término offset, ya que este se utiliza cuando es necesario ajustar la tasa de eventos por una variable de tiempo o tamaño poblacional. En este ejemplo, podemos asumir que la población de cada departamento se mantiene constante en el periodo estudiado. Además, dado que una misma persona puede infectarse más de una vez con el virus del dengue, no es necesario incorporar un término de ajuste adicional.
Chequeamos si existe sobredispersión:
check_overdispersion(fit_pois)# Overdispersion test
dispersion ratio = 1.955
Pearson's Chi-Squared = 2035.204
p-value = < 0.001El cociente de dispersión es de 1.96, por lo que el modelo Poisson no es el más adecuado para representar los datos.
Subdispersión
La subdispersión, aunque menos común, puede ocurrir en situaciones donde los conteos están más controlados o restringidos. Este fenómeno puede hacer que los errores estándar se subestimen, resultando en una sobreestimación de la significancia estadística. Matemáticamente se expresa como:
\[ Var(Y) < \mu \]
Para detectar si un modelo tiene subdispersión también se usa el ratio de dispersión. Donde un ratio significativamente menor que 1 indica subdispersión.
Control de la sobredispersión y subdispersión
Modelos quasi-Poisson
Una forma común de manejar la sobredispersión en los datos es mediante el uso de modelos quasi-Poisson. Estos modelos ajustan los errores estándar sin modificar la media predicha, incorporando un parámetro adicional que permite que la varianza sea distinta a la media. Aunque los modelos quasi-Poisson son útiles para corregir la sobredispersión, no abordan las causas subyacentes de esta ni son adecuados para datos con subdispersión o un exceso de ceros. Además, la estimación de parámetros en estos modelos tiende a ser menos eficiente, ya que asumen que la sobredispersión es constante en todas las observaciones, lo que puede llevar a inferencias menos robustas. Debido a estas limitaciones, no profundizaremos en el uso de los modelos quasi-Poisson.
Distribución binomial negativa
Dentro de los modelos lineales generalizados (GLM), la distribución binomial negativa es especialmente útil para manejar la sobredispersión en datos discretos. Esta distribución aporta flexibilidad, mejora las inferencias estadísticas con estimaciones más precisas, permite procesar grandes volúmenes de datos y se puede combinar con componentes de inflación de ceros.
La familia binomial negativa no está integrada en la función glm() del paquete stats. No obstante, existen varios paquetes que permiten su implementación. En el contexto del curso, usaremos el paquete glmmTMB (Brooks et al. 2017) para explorar alternativas a los modelos Poisson, dado su alto grado de versatilidad. La función base glmmTMB() permite ajustar desde regresiones lineales y modelos lineales generalizados (GLMs) hasta modelos mucho más complejos, lo que lo convierte en una herramienta poderosa para manejar una amplia variedad de estructuras de datos.
A modo de ejemplo, ajustaremos el modelo Poisson anterior con glmmTMB():
fit_pois_tmb <- glmmTMB(casos ~ semana_epi + grupo_edad_cat,
data = datos,
family = poisson)
# Salida modelo
summary(fit_pois_tmb) Family: poisson ( log )
Formula: casos ~ semana_epi + grupo_edad_cat
Data: datos
AIC BIC logLik deviance df.resid
1735.1 1764.8 -861.5 1723.1 1041
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -3.092153 0.279663 -11.057 < 2e-16 ***
semana_epi 0.044249 0.008993 4.920 8.64e-07 ***
grupo_edad_cat10-19 años 1.172837 0.294392 3.984 6.78e-05 ***
grupo_edad_cat20-34 años 2.195272 0.271054 8.099 5.54e-16 ***
grupo_edad_cat35-65 años 2.312217 0.269545 8.578 < 2e-16 ***
grupo_edad_cat65+ años 0.730874 0.314523 2.324 0.0201 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Al comparar los resultados del modelo ajustado con glmmTMB y el modelo ajustado con el paquete stats, podemos observar que, aunque no presentan exactamente la misma información, los coeficientes y el AIC son prácticamente idénticos.
Ahora, ajustaremos el modelo de regresión binomial negativa utilizando el argumento family = nbinom2, que tiene como función de enlace por defecto el logaritmo (link = "log"):
# Modelo binomial negativo
fit_nb <- glmmTMB(casos ~ semana_epi + grupo_edad_cat,
data = datos,
family = nbinom2)
# Salida del modelo
summary(fit_nb) Family: nbinom2 ( log )
Formula: casos ~ semana_epi + grupo_edad_cat
Data: datos
AIC BIC logLik deviance df.resid
1538.4 1573.0 -762.2 1524.4 1040
Dispersion parameter for nbinom2 family (): 0.467
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -3.27950 0.33522 -9.783 < 2e-16 ***
semana_epi 0.06159 0.01570 3.923 8.73e-05 ***
grupo_edad_cat10-19 años 1.16315 0.33067 3.518 0.000435 ***
grupo_edad_cat20-34 años 2.18318 0.30879 7.070 1.55e-12 ***
grupo_edad_cat35-65 años 2.30095 0.30722 7.490 6.90e-14 ***
grupo_edad_cat65+ años 0.71124 0.34949 2.035 0.041842 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Si comparamos los coeficientes de este modelo con el de Poisson, observamos que son similares, pero los intervalos de confianza varían ligeramente debido al ajuste binomial negativo:
Characteristic |
Poisson |
Binomial negativa |
||||
|---|---|---|---|---|---|---|
exp(Beta) |
95% CI 1 |
p-value |
exp(Beta) |
95% CI 1 |
p-value |
|
| semana_epi | 1.05 | 1.03, 1.06 | <0.001 | 1.06 | 1.03, 1.10 | <0.001 |
| grupo_edad_cat | ||||||
| grupo_edad_cat10-19 años | 3.23 | 1.81, 5.75 | <0.001 | 3.20 | 1.67, 6.12 | <0.001 |
| grupo_edad_cat20-34 años | 8.98 | 5.28, 15.3 | <0.001 | 8.87 | 4.85, 16.3 | <0.001 |
| grupo_edad_cat35-65 años | 10.1 | 5.95, 17.1 | <0.001 | 9.98 | 5.47, 18.2 | <0.001 |
| grupo_edad_cat65+ años | 2.08 | 1.12, 3.85 | 0.020 | 2.04 | 1.03, 4.04 | 0.042 |
| 1
CI = Confidence Interval |
||||||
La regresión binomial negativa tiene como limitación no ser adecuada para datos con subdispersión. Además, son modelos computacionalmente complejos y su ajuste depende de la estructura específica de los datos, lo que puede afectar la estabilidad de las estimaciones.
Distribución Conway-Maxwell Poisson
Otra alternativa a la regresión de Poisson, que permite controlar tanto sobredispersión como subdispersión y exceso de ceros es el modelo Conway-Maxwell Poisson (COM-Poisson). El mismo ofrece una generalización del modelo Poisson estándar, introduciendo un parámetro adicional que permite ajustar de forma independiente la media y la varianza.
Podemos ajustar estos modelos con la función glmmTMB() especificando el argumento family = compois(), que tiene como función de enlace por defecto el logaritmo (link = "log"):
# Modelo COMPOIS
fit_compois <- glmmTMB(casos ~ semana_epi + grupo_edad_cat,
data = datos,
family = compois)
# Salida del modelo
summary(fit_compois)Una desventaja de estos modelos es que son computacionalmente complejos y están implementados en pocos software estadísticos y paquetes de R. Por otro lado, los parámetros son menos intuitivos de interpretar, y puede presentar dificultades en la estimación de parámetros, especialmente en casos de datos limitados o de mala calidad.
Exceso de ceros (zero-inflation)
El exceso de ceros ocurre cuando se observan más ceros en los datos de lo que el modelo Poisson predice. Este fenómeno puede ser un indicativo de que el modelo Poisson no es adecuado para los datos.
Para identificar un exceso de ceros, se puede utilizar la función check_zeroinflation() del paquete performance:
check_zeroinflation(fit_pois_tmb)# Check for zero-inflation
Observed zeros: 821
Predicted zeros: 739
Ratio: 0.90Los resultados del test indican que los datos presentan un exceso de ceros. Para explorar el origen de este fenómeno, comenzamos tabulando los datos por semana epidemiológica:
datos |>
group_by(semana_epi) |>
summarise(casos = sum(casos))# A tibble: 20 × 2
semana_epi casos
<dbl> <dbl>
1 1 4
2 2 2
3 3 4
4 4 6
5 5 3
6 6 6
7 7 9
8 8 11
9 9 20
10 10 48
11 11 51
12 12 59
13 13 63
14 14 49
15 15 31
16 16 11
17 17 12
18 18 9
19 19 11
20 20 1Observamos que se reportaron casos de dengue en todas las semanas epidemiológicas. Ahora, veamos la distribución de los casos por departamento:
datos |>
group_by(departamento) |>
summarise(casos = sum(casos))# A tibble: 10 × 2
departamento casos
<chr> <dbl>
1 BAHIA BLANCA 117
2 CORONEL DORREGO 2
3 GENERAL ALVARADO 6
4 GENERAL PUEYRREDON 75
5 LA COSTA 150
6 LOBERIA 2
7 MAR CHIQUITA 8
8 NECOCHEA 14
9 PINAMAR 28
10 TRES ARROYOS 8En este caso, todos los departamentos reportaron casos de dengue durante el período de estudio. Procedamos a cruzar los datos por semana epidemiológica y departamento:
datos |>
group_by(semana_epi, departamento) |>
summarise(casos = sum(casos)) |>
print(n = 20) # Muestra las primeras 20 observaciones# A tibble: 200 × 3
# Groups: semana_epi [20]
semana_epi departamento casos
<dbl> <chr> <dbl>
1 1 BAHIA BLANCA 1
2 1 CORONEL DORREGO 0
3 1 GENERAL ALVARADO 0
4 1 GENERAL PUEYRREDON 2
5 1 LA COSTA 1
6 1 LOBERIA 0
7 1 MAR CHIQUITA 0
8 1 NECOCHEA 0
9 1 PINAMAR 0
10 1 TRES ARROYOS 0
11 2 BAHIA BLANCA 0
12 2 CORONEL DORREGO 0
13 2 GENERAL ALVARADO 0
14 2 GENERAL PUEYRREDON 1
15 2 LA COSTA 0
16 2 LOBERIA 0
17 2 MAR CHIQUITA 0
18 2 NECOCHEA 0
19 2 PINAMAR 1
20 2 TRES ARROYOS 0
# ℹ 180 more rowsEn esta tabla, se observa que no todos los departamentos en La Costa reportaron casos de dengue cada semana epidemiológica. Además, dado que los datos están separados por grupo etario, es razonable suponer que tampoco se reportaron casos para cada grupo en cada departamento y semana epidemiológica.
Para obtener una visión más detallada, generamos una tabla de frecuencia para los reportes de cero casos:
tabyl(datos$casos) |>
adorn_pct_formatting() datos$casos n percent
0 821 78.4%
1 136 13.0%
2 44 4.2%
3 23 2.2%
4 13 1.2%
5 5 0.5%
6 2 0.2%
7 1 0.1%
9 1 0.1%
12 1 0.1%La tabla muestra que el 78.4% de los datos corresponden a observaciones con cero reportes de casos de dengue. Esto es esperado en el contexto de datos de vigilancia epidemiológica, donde es común encontrar muchas observaciones sin eventos.
Este exceso de ceros se puede controlar usando modelos zero-inflated, que tendrán distribución Poisson o binomial negativa según exista o no sobredispersión de los datos.
Añadiendo el argumento zi = ~ 1, podemos controlar el exceso de ceros en el modelo de regresión binomial negativa ajustada anteriormente. Indicando que la inflación en ceros sigue un modelo de intercepto único (sin predictores adicionales):
# Modelo zero-inflated
fit_zi_nb <- glmmTMB(casos ~ semana_epi + grupo_edad_cat,
zi = ~ 1,
data = datos,
family = nbinom2)
# Salida del modelo
summary(fit_zi_nb) Family: nbinom2 ( log )
Formula: casos ~ semana_epi + grupo_edad_cat
Zero inflation: ~1
Data: datos
AIC BIC logLik deviance df.resid
1540.4 1580.0 -762.2 1524.4 1039
Dispersion parameter for nbinom2 family (): 0.467
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -3.27950 0.33522 -9.783 < 2e-16 ***
semana_epi 0.06159 0.01570 3.923 8.73e-05 ***
grupo_edad_cat10-19 años 1.16315 0.33067 3.518 0.000435 ***
grupo_edad_cat20-34 años 2.18318 0.30879 7.070 1.55e-12 ***
grupo_edad_cat35-65 años 2.30095 0.30722 7.490 6.90e-14 ***
grupo_edad_cat65+ años 0.71124 0.34949 2.035 0.041842 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Zero-inflation model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -16.69 3025.54 -0.006 0.996La salida de estos modelos se divide en tres partes que incluyen:
Conditional model: coeficientes, error estándar, valor z y significancia para el modelo de conteo binomial negativo.Zero-inflation model: coeficiente para la probabilidad de observar un cero adicional, con su respectivo error estándar, valor z y significancia. El coeficiente del intercepto no es significativo, lo que indica que la probabilidad de observar un cero adicional no está significativamente influenciada por las variables en el modelo y podría ser necesario incluir un predictor enzipara mejorar el ajuste del modelo.Estadísticas del modelo: AIC, BIC, logLik, devianza y grados de libertad residuales.
En resumen
Existen diversas soluciones para tratar la sobredispersión, subdispersión y exceso de ceros en modelos de conteo. En este documento, nos enfocamos en las estrategias más comunes que pueden ser abordadas mediante modelos lineales generalizados utilizando el paquete glmmTMB.
| Problema | Modelo |
|---|---|
| Sobredispersión | glmmTMB(formula, family = nbinom2) |
| Subdispersión | glmmTMB(formula, family = compois) |
| Sobredispersión + Exceso de Ceros | glmmTMB(formula, family = nbinom2, zi = ~1) |
Referencias
Agresti, Alan. 2015. “Chapter 75: Models for Count Data.” In Foundations of Linear and Generalized Linear Models, 228–67. Wiley Series in Probability and Statistics. Hoboken, New Jersey: John Wiley & Sons, Inc.
Bolker, et al., Ben. 2024. “GLMM FAQ: Model Extensions.” https://bbolker.github.io/mixedmodels-misc/glmmFAQ.html#.
Brooks, Mollie E., Kasper Kristensen, Koen J. van, Arni Magnusson, Casper W. Berg, Anders Nielsen, Hans J. Skaug, Martin Maechler, and Benjamin M. Bolker. 2017. “glmmTMB Balances Speed and Flexibility Among Packages for Zero-Inflated Generalized Linear Mixed Modeling” 9. https://doi.org/10.32614/RJ-2017-066.
Firke, Sam. 2023. “Janitor: Simple Tools for Examining and Cleaning Dirty Data.” https://CRAN.R-project.org/package=janitor.
R Core Team. 2024. “R: A Language and Environment for Statistical Computing.” https://www.R-project.org/.
Salinas-Rodríguez, Aarón, Betty Manrique-Espinoza, and Sandra G Sosa-Rubí. 2009. “Análisis Estadístico Para Datos de Conteo: Aplicaciones Para El Uso de Los Servicios de Salud.” Salud Pública de México 51: 397–406.
Sjoberg, Daniel D., Karissa Whiting, Michael Curry, Jessica A. Lavery, and Joseph Larmarange. 2021. “Reproducible Summary Tables with the Gtsummary Package” 13: 570–80. https://doi.org/10.32614/RJ-2021-053.
Wickham, Hadley, Mara Averick, Jennifer Bryan, Winston Chang, Lucy D’Agostino McGowan, Romain François, Garrett Grolemund, et al. 2019. “Welcome to the Tidyverse” 4: 1686. https://doi.org/10.21105/joss.01686.
Footnotes
Fuente: https://es.wikipedia.org/wiki/Localidades_balnearias_del_mar_Argentino↩︎
Fuente: Vigilancia de las enfermedades por virus del Dengue y Zika. Datos Abiertos del Ministerio de Salud. Disponible en: http://datos.salud.gob.ar/dataset/vigilancia-de-dengue-y-zika↩︎