La regresión de Poisson es el método de referencia para modelar conteos de eventos y tasas de incidencia en estudios epidemiológicos. Permite estimar la Razón de Tasas de Incidencia (IRR) ajustando simultáneamente por múltiples variables, controlando el tiempo en riesgo mediante un offset.

¿Cuándo usamos la regresión de Poisson?

La distribución de Poisson modela el número de eventos independientes que ocurren en un intervalo de tiempo o espacio fijo, cuando cada evento es raro.

🦟Casos de dengue por municipio por semana epidemiológica
🏥Ingresos hospitalarios por enfermedades respiratorias por mes
💀Muertes por causa específica según grupo etario
🧪Infecciones nosocomiales por cada 1.000 días-paciente
🚑Lesiones de tránsito por 100.000 habitantes-año
👶Defectos al nacer por cada 1.000 nacidos vivos

Supuestos clave

La distribución de Poisson asume que:

  1. Los eventos son independientes entre sí
  2. La variable de respuesta \(Y\) toma valores enteros no negativos: \(Y \in \{0, 1, 2, \ldots\}\)
  3. La media = varianza (\(\mu = \sigma^2\)) — equidispersión
  4. Los eventos son raros en relación con el total de oportunidades

⚠ Sobredispersión: Cuando la varianza observada supera la media (\(\sigma^2 > \mu\)), los errores estándar estarán subestimados y los p-valores serán incorrectos. Explorar siempre la relación media-varianza antes de interpretar resultados.

El modelo matemático

Función de distribución de Poisson

La probabilidad de observar exactamente \(k\) eventos cuando la tasa esperada es \(\mu\):

\[P(Y = k) = \frac{e^{-\mu} \cdot \mu^k}{k!}, \quad k = 0, 1, 2, \ldots\]

El modelo de regresión

El modelo liga el logaritmo de la media a los predictores:

Modelo de Poisson — función de enlace log \[\log(\mu_i) = \beta_0 + \beta_1 X_1 + \beta_2 X_2 + \ldots + \beta_k X_k\]

O equivalentemente en forma exponencial:

\[\mu_i = e^{\beta_0 + \beta_1 X_1 + \ldots + \beta_k X_k}\]

Incorporando el tiempo en riesgo (offset)

Cuando las personas tienen diferente tiempo de seguimiento, o cuando comparamos poblaciones de distinto tamaño, modelamos tasas:

Modelo de tasa con offset \[\log\!\left(\frac{\mu_i}{t_i}\right) = \beta_0 + \beta_1 X_1 + \ldots + \beta_k X_k\] \[\Leftrightarrow \quad \log(\mu_i) = \log(t_i) + \beta_0 + \beta_1 X_1 + \ldots + \beta_k X_k\]

donde \(t_i\) es el tiempo en riesgo (personas-año, días-paciente, tamaño poblacional, etc.)

El offset: log(t_i) se añade como término con coeficiente fijo igual a 1. No se estima — es un ajuste que convierte el conteo en una tasa. En R se especifica como offset(log(tiempo)) dentro de la fórmula.

Interpretación de los coeficientes

Coeficiente Interpretación directa Exponenciado
\(\beta_0\) Log de la tasa base Tasa base: \(e^{\beta_0}\)
\(\beta_j\) Cambio en log-tasa por unidad de \(X_j\) IRR = \(e^{\beta_j}\)
\(\beta_j\) (binaria) Diferencia en log-tasa expuesto vs. no expuesto IRR expuesto vs. referencia

IRR (Incidence Rate Ratio): \(e^{\beta_j}\) estima cuántas veces es mayor la tasa de incidencia en el grupo con \(X_j = 1\) respecto al grupo de referencia, manteniendo todos los demás predictores constantes. Un IRR = 2 significa el doble de tasa; IRR = 0.5 significa la mitad.

Visualizando la distribución de Poisson

# Distribuciones de Poisson con diferentes lambdas
lambdas <- c(1, 3, 7, 15)
colores  <- c("#4ade80", "#fbbf24", "#fb923c", "#f87171")

par(mfrow   = c(1, 4),
    bg       = "#0f1a12",
    col.axis = "#6b8a6e",
    col.lab  = "#d4e4d0",
    col.main = "#fbbf24",
    fg       = "#2d4030",
    mar      = c(4, 3, 3, 1))

for (i in seq_along(lambdas)) {
  lam <- lambdas[i]
  k   <- 0:round(lam * 3 + 5)
  p   <- dpois(k, lam)

  barplot(p,
          names.arg = k,
          col       = paste0(colores[i], "bb"),
          border    = colores[i],
          main      = bquote(lambda == .(lam)),
          xlab      = "Número de eventos (k)",
          ylab      = if (i == 1) "P(Y = k)" else "",
          col.main  = colores[i],
          las       = 1)

  abline(v = lam / max(k) * length(k) * 1.2,
         lty = 2, col = "#ffffff44", lwd = 1)
}

Nota: A medida que \(\lambda\) aumenta, la distribución de Poisson se aproxima a una normal. Para \(\lambda\) pequeños (eventos raros), la distribución es muy asimétrica.

Simulación de un brote — datos de ejemplo

Simulamos un estudio de cohorte que registra casos de enfermedad gastrointestinal aguda (EGA) en 5 municipios durante 12 semanas. Queremos saber si la fuente de agua y la semana epidemiológica predicen la incidencia.

set.seed(2024)

n_municipios <- 5
n_semanas    <- 12

# Marco de datos: municipio × semana
brote <- expand.grid(
  municipio = paste0("Mpio_", LETTERS[1:n_municipios]),
  semana    = 1:n_semanas
) %>%
  as_tibble() %>%
  mutate(
    # Población en riesgo por municipio (personas-semana)
    poblacion    = rep(c(12000, 8500, 15000, 6200, 9800), each = n_semanas),

    # Fuente de agua: 1 = agua no tratada (exposición)
    agua_no_trat = rep(c(1, 0, 0, 1, 0), each = n_semanas),

    # Tendencia temporal y efecto de exposición
    log_mu = log(poblacion / 10000) +          # offset proporcional
             (-2.8) +                           # intercepto (tasa base)
             0.75 * agua_no_trat +              # efecto exposición (IRR ~ 2.1)
             0.04 * semana +                    # tendencia lineal
             rnorm(n(), 0, 0.2),               # variación aleatoria

    # Conteo observado de casos EGA
    casos = rpois(n(), exp(log_mu))
  )

# Primeras filas
head(brote, 10) %>%
  kbl(booktabs = TRUE, digits = 2,
      caption = "Primeras 10 filas — datos de vigilancia EGA simulados") %>%
  kable_styling(bootstrap_options = c("hover", "striped"),
                full_width = FALSE, position = "center") %>%
  row_spec(0, bold = TRUE, color = "#fbbf24", background = "#0a180c") %>%
  row_spec(seq_len(10), color = "#d4e4d0", background = "#0f1a12") %>%
  row_spec(seq(2, 10, 2), background = "#111a12")
Primeras 10 filas — datos de vigilancia EGA simulados
municipio semana poblacion agua_no_trat log_mu casos
Mpio_A 1 12000 1 -1.63 1
Mpio_B 1 12000 1 -1.73 0
Mpio_C 1 12000 1 -1.85 0
Mpio_D 1 12000 1 -1.87 0
Mpio_E 1 12000 1 -1.60 0
Mpio_A 2 12000 1 -1.53 0
Mpio_B 2 12000 1 -1.68 0
Mpio_C 2 12000 1 -1.81 0
Mpio_D 2 12000 1 -2.03 1
Mpio_E 2 12000 1 -2.01 2

Exploración descriptiva

# Tasa por 10.000 hab-semana
brote <- brote %>%
  mutate(
    tasa_10k = casos / poblacion * 10000,
    exposicion = factor(agua_no_trat,
                        labels = c("Agua tratada", "Agua NO tratada"))
  )

# Gráfico 1: Tendencia temporal por municipio
p1 <- ggplot(brote, aes(x = semana, y = tasa_10k,
                         color = exposicion, group = municipio)) +
  geom_line(linewidth = 1.1, alpha = 0.85) +
  geom_point(size = 2, alpha = 0.7) +
  scale_color_manual(values = c("Agua tratada"    = "#4ade80",
                                "Agua NO tratada" = "#f87171")) +
  scale_x_continuous(breaks = 1:12) +
  labs(title    = "Tasa de EGA por semana epidemiológica",
       subtitle = "Por municipio y fuente de agua",
       x = "Semana epidemiológica", y = "Tasa × 10.000 hab",
       color = NULL) +
  theme_minimal(base_size = 12) +
  theme(
    plot.background  = element_rect(fill = "#0f1a12", color = NA),
    panel.background = element_rect(fill = "#172019", color = NA),
    panel.grid.major = element_line(color = "#2d4030"),
    panel.grid.minor = element_blank(),
    plot.title       = element_text(color = "#fbbf24", face = "bold"),
    plot.subtitle    = element_text(color = "#6b8a6e"),
    axis.text        = element_text(color = "#6b8a6e"),
    axis.title       = element_text(color = "#d4e4d0"),
    legend.background = element_rect(fill = "#172019", color = NA),
    legend.text      = element_text(color = "#d4e4d0"),
    strip.text       = element_text(color = "#4ade80", face = "bold")
  )

# Gráfico 2: Distribución de casos
p2 <- ggplot(brote, aes(x = casos, fill = exposicion)) +
  geom_histogram(binwidth = 1, position = "identity",
                 alpha = 0.65, color = NA) +
  scale_fill_manual(values = c("Agua tratada"    = "#4ade80",
                               "Agua NO tratada" = "#f87171")) +
  labs(title = "Distribución de conteos de casos",
       x = "Casos EGA por municipio-semana",
       y = "Frecuencia", fill = NULL) +
  theme_minimal(base_size = 12) +
  theme(
    plot.background  = element_rect(fill = "#0f1a12", color = NA),
    panel.background = element_rect(fill = "#172019", color = NA),
    panel.grid.major = element_line(color = "#2d4030"),
    panel.grid.minor = element_blank(),
    plot.title       = element_text(color = "#fbbf24", face = "bold"),
    axis.text        = element_text(color = "#6b8a6e"),
    axis.title       = element_text(color = "#d4e4d0"),
    legend.background = element_rect(fill = "#172019", color = NA),
    legend.text      = element_text(color = "#d4e4d0")
  )

gridExtra::grid.arrange(p1, p2, ncol = 2)

Verificar supuesto de equidispersión

brote %>%
  group_by(exposicion) %>%
  summarise(
    Media    = mean(casos),
    Varianza = var(casos),
    Razon_V_M = round(var(casos) / mean(casos), 2)
  ) %>%
  kbl(booktabs = TRUE, digits = 2,
      caption = "Verificación de equidispersión por grupo de exposición") %>%
  kable_styling(bootstrap_options = c("hover", "striped"),
                full_width = FALSE, position = "center") %>%
  row_spec(0, bold = TRUE, color = "#4ade80", background = "#0a180c") %>%
  row_spec(1:2, color = "#d4e4d0", background = "#0f1a12") %>%
  row_spec(2, background = "#111a12")
Verificación de equidispersión por grupo de exposición
exposicion Media Varianza Razon_V_M
Agua tratada 0.06 0.05 0.97
Agua NO tratada 0.17 0.23 1.39

Razón varianza/media: Si este valor es cercano a 1, los datos son equidispersos y la regresión de Poisson estándar es adecuada. Si es mayor que 2, considere quasi-Poisson o Binomial Negativa.

Ajustando el modelo de Poisson

Paso a paso

1
Modelo nulo — solo el intercepto con offset. Estima la tasa basal.
2
Modelo simple — una sola exposición. Análogo a la tabla 2×2 en análisis bivariado.
3
Modelo ajustado — múltiples predictores. Controla confusión simultáneamente.
4
Diagnóstico — verificar sobredispersión, residuales, ajuste global.
5
Interpretación — exponenciar coeficientes para obtener IRR con IC 95%.

Modelo 1 — Nulo (solo offset)

mod_nulo <- glm(casos ~ 1 + offset(log(poblacion)),
                data   = brote,
                family = poisson(link = "log"))

summary(mod_nulo)
## 
## Call:
## glm(formula = casos ~ 1 + offset(log(poblacion)), family = poisson(link = "log"), 
##     data = brote)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -11.5425     0.4082  -28.27   <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: 27.678  on 59  degrees of freedom
## Residual deviance: 27.678  on 59  degrees of freedom
## AIC: 40.292
## 
## Number of Fisher Scoring iterations: 6

Tasa basal estimada: exp(coef(mod_nulo)[1]) × 10.000 = 0.1 casos por 10.000 hab-semana (sin ajustar por exposición).

Modelo 2 — Bivariado: agua no tratada

Bivariado offset: log(poblacion)

mod_bivariado <- glm(casos ~ agua_no_trat + offset(log(poblacion)),
                     data   = brote,
                     family = poisson(link = "log"))

summary(mod_bivariado)
## 
## Call:
## glm(formula = casos ~ agua_no_trat + offset(log(poblacion)), 
##     family = poisson(link = "log"), data = brote)
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -12.2051     0.7071 -17.261   <2e-16 ***
## agua_no_trat   1.2973     0.8660   1.498    0.134    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 27.678  on 59  degrees of freedom
## Residual deviance: 25.251  on 58  degrees of freedom
## AIC: 39.864
## 
## Number of Fisher Scoring iterations: 6

Extracción de IRR y IC 95%

# Exponenciar coeficientes → IRR con IC 95%
irr_biv <- tidy(mod_bivariado, exponentiate = TRUE, conf.int = TRUE) %>%
  rename(IRR = estimate, IC_bajo = conf.low, IC_alto = conf.high) %>%
  dplyr::select(term, IRR, IC_bajo, IC_alto, p.value)

irr_biv %>%
  kbl(booktabs = TRUE, digits = 3,
      caption = "IRR — Modelo bivariado: agua no tratada") %>%
  kable_styling(bootstrap_options = c("hover", "striped"),
                full_width = FALSE, position = "center") %>%
  row_spec(0, bold = TRUE, color = "#fbbf24", background = "#0a180c") %>%
  row_spec(1:nrow(irr_biv), color = "#d4e4d0", background = "#0f1a12") %>%
  row_spec(seq(2, nrow(irr_biv), 2), background = "#111a12")
IRR — Modelo bivariado: agua no tratada
term IRR IC_bajo IC_alto p.value
(Intercept) 0.000 0.000 0.000 0.000
agua_no_trat 3.659 0.714 26.398 0.134

Lectura epidemiológica: Los municipios con agua no tratada tienen aproximadamente 3.66 veces la tasa de casos de EGA comparados con los municipios con agua tratada (IC 95%: 0.71 – 26.4).

Modelo 3 — Ajustado: agua + semana

Ajustado agua + tendencia temporal

mod_ajustado <- glm(casos ~ agua_no_trat + semana + offset(log(poblacion)),
                    data   = brote,
                    family = poisson(link = "log"))

summary(mod_ajustado)
## 
## Call:
## glm(formula = casos ~ agua_no_trat + semana + offset(log(poblacion)), 
##     family = poisson(link = "log"), data = brote)
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -10.5863     1.3110  -8.075 6.74e-16 ***
## agua_no_trat   0.4413     1.0996   0.401    0.688    
## semana        -0.2534     0.1957  -1.295    0.195    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 27.678  on 59  degrees of freedom
## Residual deviance: 22.808  on 57  degrees of freedom
## AIC: 39.422
## 
## Number of Fisher Scoring iterations: 6
irr_adj <- tidy(mod_ajustado, exponentiate = TRUE, conf.int = TRUE) %>%
  rename(IRR = estimate, IC_bajo = conf.low, IC_alto = conf.high) %>%
  dplyr::select(term, IRR, IC_bajo, IC_alto, p.value)

irr_adj %>%
  kbl(booktabs = TRUE, digits = 3,
      caption = "IRR — Modelo ajustado: agua + semana epidemiológica") %>%
  kable_styling(bootstrap_options = c("hover", "striped"),
                full_width = FALSE, position = "center") %>%
  row_spec(0, bold = TRUE, color = "#fbbf24", background = "#0a180c") %>%
  row_spec(1:nrow(irr_adj), color = "#d4e4d0", background = "#0f1a12") %>%
  row_spec(seq(2, nrow(irr_adj), 2), background = "#111a12")
IRR — Modelo ajustado: agua + semana epidemiológica
term IRR IC_bajo IC_alto p.value
(Intercept) 0.000 0.000 0.000 0.000
agua_no_trat 1.555 0.165 14.899 0.688
semana 0.776 0.465 1.055 0.195

Resultado ajustado: Controlando la tendencia temporal, el IRR del agua no tratada es 1.55 (IC 95%: 0.16 – 14.9). Cada semana adicional multiplica la tasa en 0.776 (e^-0.253).

Diagnóstico del modelo

Prueba de sobredispersión

# Estadístico de dispersión: deviance / gl debe ser ≈ 1
dispersion_stat <- deviance(mod_ajustado) / df.residual(mod_ajustado)
cat("Estadístico de dispersión (Deviance/gl):", round(dispersion_stat, 3), "\n")
## Estadístico de dispersión (Deviance/gl): 0.4
# Prueba formal con MASS::glm.nb como referencia
# Si el estadístico > 2, considerar quasi-Poisson
if (dispersion_stat > 2) {
  cat("⚠ SOBREDISPERSIÓN detectada — considerar modelo quasi-Poisson\n")
} else {
  cat("✓ Equidispersión aceptable para el modelo de Poisson estándar\n")
}
## ✓ Equidispersión aceptable para el modelo de Poisson estándar

Gráficos de diagnóstico

par(mfrow   = c(1, 3),
    bg       = "#0f1a12",
    col.axis = "#6b8a6e",
    col.lab  = "#d4e4d0",
    col.main = "#fbbf24",
    fg       = "#2d4030",
    mar      = c(4, 4, 3, 1))

# 1. Residuales de Pearson vs valores ajustados
plot(fitted(mod_ajustado),
     residuals(mod_ajustado, type = "pearson"),
     xlab = "Valores ajustados (μ̂)",
     ylab = "Residuales de Pearson",
     main = "Residuales vs Ajustados",
     pch  = 19, col = "#4ade8088", cex = 0.9)
abline(h = 0, lty = 2, col = "#fbbf2488", lwd = 1.5)
abline(h = c(-2, 2), lty = 3, col = "#f8717188", lwd = 1)

# 2. QQ-plot de residuales de deviance
qqnorm(residuals(mod_ajustado, type = "deviance"),
       pch  = 19, col = "#fbbf2488", cex = 0.9,
       main = "QQ-plot — Residuales deviance",
       xlab = "Cuantiles teóricos",
       ylab = "Cuantiles observados")
qqline(residuals(mod_ajustado, type = "deviance"),
       col = "#4ade80", lwd = 1.5, lty = 2)

# 3. Valores ajustados vs observados
plot(brote$casos,
     fitted(mod_ajustado),
     xlab = "Casos observados",
     ylab = "Casos esperados (μ̂)",
     main = "Observado vs Esperado",
     pch  = 19, col = "#7dd3fc88", cex = 0.9)
abline(a = 0, b = 1, col = "#4ade80", lwd = 1.5, lty = 2)

Diagnóstico visual: En el gráfico de residuales de Pearson, los puntos deben distribuirse aleatoriamente alrededor de cero. Patrones en abanico sugieren sobredispersión. En el gráfico observado vs. esperado, los puntos deben aproximarse a la diagonal.

Comparación de modelos (AIC)

modelos_comp <- tibble(
  Modelo     = c("Nulo", "Bivariado (agua)", "Ajustado (agua + semana)"),
  AIC        = c(AIC(mod_nulo), AIC(mod_bivariado), AIC(mod_ajustado)),
  Deviance   = c(deviance(mod_nulo), deviance(mod_bivariado), deviance(mod_ajustado)),
  GL         = c(df.residual(mod_nulo), df.residual(mod_bivariado), df.residual(mod_ajustado))
) %>%
  mutate(Delta_AIC = round(AIC - min(AIC), 1),
         AIC = round(AIC, 1), Deviance = round(Deviance, 1))

modelos_comp %>%
  kbl(booktabs = TRUE,
      caption = "Comparación de modelos — criterio AIC") %>%
  kable_styling(bootstrap_options = c("hover", "striped"),
                full_width = FALSE, position = "center") %>%
  row_spec(0, bold = TRUE, color = "#fbbf24", background = "#0a180c") %>%
  row_spec(1:3, color = "#d4e4d0", background = "#0f1a12") %>%
  row_spec(which.min(modelos_comp$AIC), bold = TRUE,
           color = "#4ade80", background = "#0e220e") %>%
  row_spec(seq(2, 3, 2), background = "#111a12")
Comparación de modelos — criterio AIC
Modelo AIC Deviance GL Delta_AIC
Nulo 40.3 27.7 59 0.9
Bivariado (agua) 39.9 25.3 58 0.4
Ajustado (agua + semana) 39.4 22.8 57 0.0

Criterio AIC: El modelo con menor AIC tiene mejor ajuste penalizado por complejidad. Un \(\Delta\)AIC > 10 indica diferencia sustancial. El modelo resaltado en verde es el seleccionado.

Visualización de resultados

Forest plot de IRR

# Extraer todos los predictores (excluir intercepto)
irr_plot <- irr_adj %>%
  filter(term != "(Intercept)") %>%
  mutate(
    term = recode(term,
      "agua_no_trat" = "Agua NO tratada\n(vs. tratada)",
      "semana"       = "Semana epidemiológica\n(por unidad)"
    )
  )

ggplot(irr_plot, aes(x = IRR, y = fct_rev(factor(term)))) +
  geom_vline(xintercept = 1, linetype = "dashed",
             color = "#fbbf24", linewidth = 0.8, alpha = 0.7) +
  geom_errorbarh(aes(xmin = IC_bajo, xmax = IC_alto),
                 height = 0.15, color = "#6b8a6e", linewidth = 1) +
  geom_point(size = 4, color = "#fb923c") +
  geom_text(aes(label = sprintf("IRR = %.2f\n(%.2f–%.2f)",
                                IRR, IC_bajo, IC_alto)),
            hjust = -0.12, color = "#d4e4d0", size = 3.2) +
  scale_x_log10(limits = c(0.5, 8),
                breaks  = c(0.5, 1, 2, 3, 5)) +
  labs(title    = "Razones de Tasa de Incidencia (IRR)",
       subtitle = "Modelo ajustado con offset log(población)",
       x = "IRR (escala logarítmica)", y = NULL) +
  theme_minimal(base_size = 12) +
  theme(
    plot.background  = element_rect(fill = "#0f1a12", color = NA),
    panel.background = element_rect(fill = "#172019", color = NA),
    panel.grid.major = element_line(color = "#2d4030"),
    panel.grid.minor = element_blank(),
    plot.title       = element_text(color = "#fbbf24", face = "bold", size = 13),
    plot.subtitle    = element_text(color = "#6b8a6e"),
    axis.text        = element_text(color = "#d4e4d0"),
    axis.title       = element_text(color = "#d4e4d0")
  )

Manejo de la sobredispersión

Si el estadístico de dispersión es mayor que 1.5–2, hay dos estrategias principales:

Opción A — Quasi-Poisson

Corrige los errores estándar multiplicándolos por \(\sqrt{\hat{\phi}}\). No cambia los coeficientes, solo los errores e IC.

mod_quasi <- glm(casos ~ agua_no_trat + semana + offset(log(poblacion)),
                 data   = brote,
                 family = quasipoisson(link = "log"))

# Comparar errores estándar
cbind(
  Poisson      = round(coef(summary(mod_ajustado))[, "Std. Error"], 4),
  QuasiPoisson = round(coef(summary(mod_quasi))[, "Std. Error"], 4)
) %>%
  as.data.frame() %>%
  kbl(booktabs = TRUE,
      caption = "Comparación de errores estándar: Poisson vs. Quasi-Poisson") %>%
  kable_styling(bootstrap_options = c("hover", "striped"),
                full_width = FALSE, position = "center") %>%
  row_spec(0, bold = TRUE, color = "#4ade80", background = "#0a180c") %>%
  row_spec(1:3, color = "#d4e4d0", background = "#0f1a12") %>%
  row_spec(2, background = "#111a12")
Comparación de errores estándar: Poisson vs. Quasi-Poisson
Poisson QuasiPoisson
(Intercept) 1.3110 1.1248
agua_no_trat 1.0996 0.9434
semana 0.1957 0.1679

Opción B — Binomial Negativa

Añade un parámetro de dispersión \(\theta\) explícito. Más flexible que quasi-Poisson.

mod_nb <- glm.nb(casos ~ agua_no_trat + semana + offset(log(poblacion)),
                 data = brote)

# IRR del modelo binomial negativo
tidy(mod_nb, exponentiate = TRUE, conf.int = TRUE) %>%
  rename(IRR = estimate, IC_bajo = conf.low, IC_alto = conf.high) %>%
  dplyr::select(term, IRR, IC_bajo, IC_alto, p.value) %>%
  kbl(booktabs = TRUE, digits = 3,
      caption = "IRR — Modelo Binomial Negativo") %>%
  kable_styling(bootstrap_options = c("hover", "striped"),
                full_width = FALSE, position = "center") %>%
  row_spec(0, bold = TRUE, color = "#fbbf24", background = "#0a180c") %>%
  row_spec(1:3, color = "#d4e4d0", background = "#0f1a12") %>%
  row_spec(2, background = "#111a12")
IRR — Modelo Binomial Negativo
term IRR IC_bajo IC_alto p.value
(Intercept) 0.000 0.000 0.000 0.000
agua_no_trat 1.508 0.142 15.299 0.717
semana 0.771 0.449 1.063 0.199

¿Cuándo usar cada uno? La quasi-Poisson es suficiente para corrección de errores estándar con sobredispersión moderada. La Binomial Negativa es preferible cuando la sobredispersión es severa o tiene una interpretación biológica (heterogeneidad individual en el riesgo).

Tabla resumen para el informe epidemiológico

# Tabla de resultados lista para publicar
tabla_final <- bind_rows(
  tidy(mod_bivariado, exponentiate = TRUE, conf.int = TRUE) %>%
    filter(term != "(Intercept)") %>%
    mutate(Modelo = "Bivariado"),
  tidy(mod_ajustado, exponentiate = TRUE, conf.int = TRUE) %>%
    filter(term != "(Intercept)") %>%
    mutate(Modelo = "Ajustado")
) %>%
  mutate(
    Variable = recode(term,
      "agua_no_trat" = "Agua no tratada",
      "semana"       = "Semana epidemiológica"),
    IRR      = round(estimate, 2),
    IC_95    = sprintf("%.2f – %.2f", conf.low, conf.high),
    p_valor  = case_when(
      p.value < 0.001 ~ "< 0.001",
      p.value < 0.01  ~ sprintf("%.3f", p.value),
      TRUE            ~ sprintf("%.2f",  p.value)
    )
  ) %>%
  dplyr::select(Modelo, Variable, IRR, IC_95, p_valor)

tabla_final %>%
  kbl(booktabs = TRUE,
      caption = "Resultados del análisis de Poisson — EGA por municipio") %>%
  kable_styling(bootstrap_options = c("hover", "striped"),
                full_width = FALSE, position = "center") %>%
  row_spec(0, bold = TRUE, color = "#fbbf24", background = "#0a180c") %>%
  row_spec(seq_len(nrow(tabla_final)),
           color = "#d4e4d0", background = "#0f1a12") %>%
  row_spec(seq(2, nrow(tabla_final), 2), background = "#111a12") %>%
  pack_rows("Bivariado", 1, 1, color = "#4ade80",
            background = "#0e1f10", bold = FALSE) %>%
  pack_rows("Ajustado",  2, 3, color = "#fbbf24",
            background = "#1a1500", bold = FALSE)
Resultados del análisis de Poisson — EGA por municipio
Modelo Variable IRR IC_95 p_valor
Bivariado
Bivariado Agua no tratada 3.66 0.71 – 26.40 0.13
Ajustado
Ajustado Agua no tratada 1.55 0.16 – 14.90 0.69
Ajustado Semana epidemiológica 0.78 0.46 – 1.05 0.20

Ejemplo publicado — Melanoma y radiación solar

Este es un ejemplo clásico de la literatura epidemiológica, basado en los datos de Kemp et al. sobre mortalidad por melanoma en 48 estados de Estados Unidos, disponibles en el paquete ISwR de R (Dalgaard, 2008 — Introductory Statistics with R). El estudio examina si la radiación solar ultravioleta explica las diferencias geográficas en tasas de mortalidad por melanoma.

Fuente: Kemp IW, Fraumeni JF Jr., Fears TR (1975). Solar radiation and cancer mortality in the United States. En: Persons at High Risk of Cancer. Academic Press. Reproducido como dataset melanoma en el paquete ISwR.

Descripción del dataset

El dataset contiene datos de mortalidad por melanoma en 48 estados contiguos de EE.UU.:

  • mort — Tasa de mortalidad por melanoma (muertes por millón de habitantes)
  • lat — Latitud del centro del estado (grados norte)
  • lon — Longitud del centro del estado (grados oeste)
  • ocean — Exposición oceánica: 1 = estado costero, 0 = interior
  • UV — Índice de radiación ultravioleta solar
# Datos originales de Kemp et al. (1975) — 48 estados contiguos de EE.UU.
# Fuente: Dalgaard P. Introductory Statistics with R. Springer, 2008. Cap. 13
melanoma <- data.frame(
  state = c("Alabama","Arizona","Arkansas","California","Colorado",
            "Connecticut","Delaware","Florida","Georgia","Idaho",
            "Illinois","Indiana","Iowa","Kansas","Kentucky",
            "Louisiana","Maine","Maryland","Massachusetts","Michigan",
            "Minnesota","Mississippi","Missouri","Montana","Nebraska",
            "Nevada","New Hampshire","New Jersey","New Mexico","New York",
            "North Carolina","North Dakota","Ohio","Oklahoma","Oregon",
            "Pennsylvania","Rhode Island","South Carolina","South Dakota",
            "Tennessee","Texas","Utah","Vermont","Virginia",
            "Washington","West Virginia","Wisconsin","Wyoming"),
  mort  = c(219,160,170,182,149,159,200,217,228,116,
            129,132,116,143,147,190,117,162,143,131,
            109,207,131,129,133,172,125,159,141,152,
            199,115,131,182,136,133,137,178,135,150,
            192,168,124,153,124,152,119,134),
  lat   = c(33.0,34.5,34.7,37.5,39.0,41.8,39.0,28.0,32.5,44.5,
            40.0,40.2,42.0,38.7,37.5,31.2,45.2,39.2,42.3,43.5,
            46.0,32.5,38.5,47.0,41.5,39.5,43.7,40.2,35.0,43.0,
            35.5,47.5,40.2,35.5,44.0,40.8,41.7,33.5,44.5,36.0,
            31.5,39.5,44.0,37.5,47.5,38.5,44.5,43.0),
  lon   = c(86.8,112.0,92.3,119.7,105.5,72.7,75.5,81.5,83.5,114.0,
            89.5,86.2,93.7,98.5,84.8,91.8,69.0,76.7,71.8,84.7,
            94.7,89.8,92.5,110.5,98.5,117.0,71.5,74.5,106.0,75.0,
            79.5,100.5,82.8,97.5,120.5,77.8,71.5,80.8,100.5,86.2,
            98.0,111.5,72.7,78.5,120.5,80.8,90.5,107.5),
  ocean = c(1,0,0,1,0,1,1,1,1,0,
            0,0,0,0,0,1,1,1,1,0,
            0,1,0,0,0,0,1,1,0,1,
            1,0,0,0,1,1,1,1,0,0,
            1,0,1,1,1,0,0,0),
  UV    = c(3.5,5.3,4.3,4.8,4.5,3.0,3.8,5.4,4.8,4.0,
            3.3,3.6,3.1,4.0,3.5,4.7,2.6,3.8,3.0,3.2,
            2.8,4.3,3.7,3.8,3.5,4.8,3.0,3.4,5.0,3.2,
            4.3,3.2,3.4,4.3,3.5,3.3,3.2,4.5,3.5,3.8,
            5.3,4.6,3.0,3.8,3.5,3.3,3.0,4.5)
)

head(melanoma, 8) %>%
  kbl(booktabs = TRUE,
      caption   = "Dataset melanoma — primeras 8 observaciones (Kemp et al., 1975)") %>%
  kable_styling(bootstrap_options = c("hover", "striped"),
                full_width = FALSE, position = "center") %>%
  row_spec(0, bold = TRUE, color = "#fbbf24", background = "#0a180c") %>%
  row_spec(seq_len(8), color = "#d4e4d0", background = "#0f1a12") %>%
  row_spec(seq(2, 8, 2), background = "#111a12")
Dataset melanoma — primeras 8 observaciones (Kemp et al., 1975)
state mort lat lon ocean UV
Alabama 219 33.0 86.8 1 3.5
Arizona 160 34.5 112.0 0 5.3
Arkansas 170 34.7 92.3 0 4.3
California 182 37.5 119.7 1 4.8
Colorado 149 39.0 105.5 0 4.5
Connecticut 159 41.8 72.7 1 3.0
Delaware 200 39.0 75.5 1 3.8
Florida 217 28.0 81.5 1 5.4

Exploración visual

p_uv <- ggplot(melanoma, aes(x = UV, y = mort)) +
  geom_point(color = "#fbbf24", size = 3, alpha = 0.85) +
  geom_smooth(method = "loess", se = TRUE,
              color = "#4ade80", fill = "#4ade8030", linewidth = 1.2) +
  geom_smooth(method = "lm", se = FALSE,
              color = "#f97316", linetype = "dashed", linewidth = 0.9) +
  labs(title    = "Mortalidad por melanoma vs. Radiación UV",
       subtitle = "48 estados de EE.UU. · línea verde = LOESS · naranja = lineal",
       x = "Índice de radiación UV",
       y = "Muertes por millón de hab.") +
  theme_minimal(base_size = 12) +
  theme(
    plot.background  = element_rect(fill = "#0f1a12", color = NA),
    panel.background = element_rect(fill = "#172019", color = NA),
    panel.grid.major = element_line(color = "#2d4030"),
    panel.grid.minor = element_blank(),
    plot.title       = element_text(color = "#fbbf24", face = "bold"),
    plot.subtitle    = element_text(color = "#6b8a6e", size = 9),
    axis.text        = element_text(color = "#6b8a6e"),
    axis.title       = element_text(color = "#d4e4d0")
  )

p_lat <- ggplot(melanoma, aes(x = lat, y = mort,
                               color = factor(ocean),
                               shape = factor(ocean))) +
  geom_point(size = 3, alpha = 0.85) +
  geom_smooth(method = "lm", se = FALSE, linewidth = 1) +
  scale_color_manual(values = c("0" = "#7dd3fc", "1" = "#f97316"),
                     labels = c("Interior", "Costero")) +
  scale_shape_manual(values = c("0" = 16, "1" = 17),
                     labels = c("Interior", "Costero")) +
  labs(title    = "Mortalidad por melanoma vs. Latitud",
       subtitle = "Por exposición oceánica",
       x = "Latitud (grados norte)",
       y = "Muertes por millón de hab.",
       color = NULL, shape = NULL) +
  theme_minimal(base_size = 12) +
  theme(
    plot.background  = element_rect(fill = "#0f1a12", color = NA),
    panel.background = element_rect(fill = "#172019", color = NA),
    panel.grid.major = element_line(color = "#2d4030"),
    panel.grid.minor = element_blank(),
    plot.title       = element_text(color = "#fbbf24", face = "bold"),
    plot.subtitle    = element_text(color = "#6b8a6e", size = 9),
    axis.text        = element_text(color = "#6b8a6e"),
    axis.title       = element_text(color = "#d4e4d0"),
    legend.background = element_rect(fill = "#172019", color = NA),
    legend.text      = element_text(color = "#d4e4d0")
  )

gridExtra::grid.arrange(p_uv, p_lat, ncol = 2)

Patrón visual: A mayor radiación UV (y menor latitud), mayor mortalidad por melanoma. Los estados costeros parecen tener tasas ligeramente diferentes a los del interior con similar latitud — posible modificación de efecto.

Ajuste del modelo de Poisson

En este dataset la variable mort ya es una tasa (por millón de habitantes), no un conteo bruto con tiempo en riesgo conocido, por lo que modelamos directamente el log de la tasa como función de los predictores.

Modelo bivariado UV solamente

# Modelo 1: solo radiación UV
mel_uv <- glm(mort ~ UV,
              data   = melanoma,
              family = poisson(link = "log"))

summary(mel_uv)
## 
## Call:
## glm(formula = mort ~ UV, family = poisson(link = "log"), data = melanoma)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  4.36208    0.06313   69.09   <2e-16 ***
## UV           0.17211    0.01583   10.87   <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: 278.20  on 47  degrees of freedom
## Residual deviance: 162.08  on 46  degrees of freedom
## AIC: 494.83
## 
## Number of Fisher Scoring iterations: 4

Modelo ajustado UV + latitud + costa

# Modelo 2: UV + latitud + exposición oceánica
mel_adj <- glm(mort ~ UV + lat + ocean,
               data   = melanoma,
               family = poisson(link = "log"))

summary(mel_adj)
## 
## Call:
## glm(formula = mort ~ UV + lat + ocean, family = poisson(link = "log"), 
##     data = melanoma)
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  5.819930   0.240108  24.239  < 2e-16 ***
## UV           0.045486   0.024614   1.848   0.0646 .  
## lat         -0.026248   0.003942  -6.658 2.77e-11 ***
## ocean        0.118943   0.025009   4.756 1.98e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 278.205  on 47  degrees of freedom
## Residual deviance:  58.974  on 44  degrees of freedom
## AIC: 395.73
## 
## Number of Fisher Scoring iterations: 4

IRR del modelo ajustado

irr_mel <- tidy(mel_adj, exponentiate = TRUE, conf.int = TRUE) %>%
  rename(IRR = estimate, IC_bajo = conf.low, IC_alto = conf.high) %>%
  dplyr::select(term, IRR, IC_bajo, IC_alto, p.value) %>%
  filter(term != "(Intercept)") %>%
  mutate(
    Variable = recode(term,
      "UV"    = "Radiación UV (por unidad)",
      "lat"   = "Latitud — grados norte (por grado)",
      "ocean" = "Estado costero vs. interior"),
    p_valor = case_when(
      p.value < 0.001 ~ "< 0.001",
      p.value < 0.01  ~ sprintf("%.3f", p.value),
      TRUE            ~ sprintf("%.3f", p.value)
    )
  ) %>%
  dplyr::select(Variable, IRR, IC_bajo, IC_alto, p_valor)

irr_mel %>%
  kbl(booktabs = TRUE, digits = 3,
      caption = "IRR — Modelo de Poisson ajustado: melanoma (Kemp et al.)") %>%
  kable_styling(bootstrap_options = c("hover", "striped"),
                full_width = FALSE, position = "center") %>%
  row_spec(0, bold = TRUE, color = "#fbbf24", background = "#0a180c") %>%
  row_spec(seq_len(nrow(irr_mel)), color = "#d4e4d0", background = "#0f1a12") %>%
  row_spec(seq(2, nrow(irr_mel), 2), background = "#111a12")
IRR — Modelo de Poisson ajustado: melanoma (Kemp et al.)
Variable IRR IC_bajo IC_alto p_valor
Radiación UV (por unidad) 1.047 0.997 1.098 0.065
Latitud — grados norte (por grado) 0.974 0.967 0.982 < 0.001
Estado costero vs. interior 1.126 1.072 1.183 < 0.001

Diagnóstico de sobredispersión

cat("Deviance / gl — UV solo:    ",
    round(deviance(mel_uv)  / df.residual(mel_uv),  3), "\n")
## Deviance / gl — UV solo:     3.523
cat("Deviance / gl — Ajustado:   ",
    round(deviance(mel_adj) / df.residual(mel_adj), 3), "\n")
## Deviance / gl — Ajustado:    1.34

Sobredispersión marcada: El estadístico deviance/gl supera ampliamente 1 en ambos modelos — habitual cuando la variable respuesta es una tasa ya agregada (no un conteo individual). En este contexto se recomienda usar quasi-Poisson para corregir los errores estándar.

# Modelo quasi-Poisson para corregir errores estándar
mel_quasi <- glm(mort ~ UV + lat + ocean,
                 data   = melanoma,
                 family = quasipoisson(link = "log"))

# Comparar IC: Poisson estándar vs. quasi-Poisson
bind_rows(
  tidy(mel_adj,   exponentiate = TRUE, conf.int = TRUE) %>%
    mutate(Modelo = "Poisson"),
  tidy(mel_quasi, exponentiate = TRUE, conf.int = TRUE) %>%
    mutate(Modelo = "Quasi-Poisson")
) %>%
  filter(term != "(Intercept)") %>%
  mutate(
    Variable = recode(term,
      "UV"    = "Radiación UV",
      "lat"   = "Latitud",
      "ocean" = "Costa"),
    IRR  = round(estimate, 3),
    IC   = sprintf("%.3f – %.3f", conf.low, conf.high)
  ) %>%
  dplyr::select(Modelo, Variable, IRR, IC) %>%
  kbl(booktabs = TRUE,
      caption = "Comparación Poisson vs. Quasi-Poisson — IC 95% más amplios con quasi") %>%
  kable_styling(bootstrap_options = c("hover", "striped"),
                full_width = FALSE, position = "center") %>%
  row_spec(0, bold = TRUE, color = "#4ade80", background = "#0a180c") %>%
  row_spec(seq(1, 6), color = "#d4e4d0", background = "#0f1a12") %>%
  row_spec(seq(2, 6, 2), background = "#111a12")
Comparación Poisson vs. Quasi-Poisson — IC 95% más amplios con quasi
Modelo Variable IRR IC
Poisson Radiación UV 1.047 0.997 – 1.098
Poisson Latitud 0.974 0.967 – 0.982
Poisson Costa 1.126 1.072 – 1.183
Quasi-Poisson Radiación UV 1.047 0.989 – 1.107
Quasi-Poisson Latitud 0.974 0.965 – 0.983
Quasi-Poisson Costa 1.126 1.064 – 1.193

Forest plot — modelo ajustado quasi-Poisson

irr_forest <- tidy(mel_quasi, exponentiate = TRUE, conf.int = TRUE) %>%
  filter(term != "(Intercept)") %>%
  mutate(
    Variable = recode(term,
      "UV"    = "Radiación UV\n(por unidad de índice)",
      "lat"   = "Latitud\n(por grado norte)",
      "ocean" = "Estado costero\n(vs. interior)")
  )

ggplot(irr_forest, aes(x = estimate, y = fct_rev(factor(Variable)))) +
  geom_vline(xintercept = 1, linetype = "dashed",
             color = "#fbbf24", linewidth = 0.8, alpha = 0.6) +
  geom_errorbarh(aes(xmin = conf.low, xmax = conf.high),
                 height = 0.15, color = "#6b8a6e", linewidth = 1.1) +
  geom_point(size = 4.5, color = "#f97316") +
  geom_text(aes(label = sprintf("IRR = %.3f\n(%.3f – %.3f)",
                                estimate, conf.low, conf.high)),
            hjust   = -0.1,
            color   = "#d4e4d0",
            size    = 3.1) +
  scale_x_log10(limits = c(0.85, 1.6),
                breaks  = c(0.9, 1.0, 1.1, 1.2, 1.4)) +
  labs(title    = "Mortalidad por melanoma — IRR ajustados",
       subtitle = "Modelo quasi-Poisson · Kemp et al. (1975)",
       x = "IRR (escala log)", y = NULL) +
  theme_minimal(base_size = 12) +
  theme(
    plot.background  = element_rect(fill = "#0f1a12", color = NA),
    panel.background = element_rect(fill = "#172019", color = NA),
    panel.grid.major = element_line(color = "#2d4030"),
    panel.grid.minor = element_blank(),
    plot.title       = element_text(color = "#fbbf24", face = "bold", size = 13),
    plot.subtitle    = element_text(color = "#6b8a6e", size = 9),
    axis.text        = element_text(color = "#d4e4d0"),
    axis.title       = element_text(color = "#d4e4d0")
  )

Interpretación epidemiológica

Radiación UV: Cada unidad adicional en el índice de radiación UV se asocia con un aumento de aproximadamente 4.7% en la tasa de mortalidad por melanoma (IRR = 1.047), ajustando por latitud y exposición oceánica.

Latitud: Por cada grado adicional hacia el norte (menor exposición solar acumulada), la tasa de mortalidad disminuye aproximadamente 2.6% (IRR = 0.974).

Lección metodológica de este ejemplo: La radiación UV y la latitud están altamente correlacionadas (a menor latitud → mayor UV). Al incluir ambas en el modelo, sus efectos se separan: UV captura la exposición acumulada real, mientras que latitud puede absorber otros factores geográficos (dieta, acceso a salud, temperatura). Siempre evalúe colinealidad antes de interpretar coeficientes individuales en modelos ajustados.

✓
Use regresión de Poisson cuando el desenlace es un conteo de eventos raros (casos, muertes, hospitalizaciones).
✓
Siempre incluya el offset log(población) o log(tiempo en riesgo) para modelar tasas en lugar de conteos absolutos.
✓
El IRR es el indicador a reportar. Es el análogo del RR de los estudios de cohorte. Se interpreta igual: veces más (o menos) de tasa en expuestos vs. no expuestos.
✓
Verifique la sobredispersión (deviance/gl ≈ 1). Si hay sobredispersión, use quasi-Poisson o Binomial Negativa para no subestimar los errores estándar ni el intervalo de confianza.
✓
Compare modelos con AIC. Construya el modelo de manera jerárquica: nulo → bivariado → ajustado → con interacciones si hay evidencia biológica de modificación de efecto.
“El objetivo del análisis multivariado no es controlar estadísticamente todas las variables posibles, sino construir un modelo que represente el mecanismo causal que el epidemiólogo tiene en mente.” — Rothman & Greenland