4  Modelos para Datos de Recuento

4.1 Cuando el resultado es un número entero no negativo

Muchas variables de interés en economía y ciencias sociales no son continuas ni binarias: son enteros no negativos. El número de visitas al médico en un año, los accidentes de tráfico en un tramo de carretera, las publicaciones académicas de un investigador o los delitos registrados en un municipio son ejemplos de datos de recuento. Comparten tres rasgos que hacen inapropiada la regresión lineal ordinaria.

En primer lugar, solo pueden tomar valores en \(\{0, 1, 2, \ldots\}\): nunca serán negativos ni fraccionarios. En segundo lugar, su distribución es asimétricamente positiva: hay muchos valores pequeños y pocos muy elevados. En tercer lugar, la varianza suele crecer con la media — fenómeno conocido como sobredispersión — de modo que los errores estándar de MCO son incorrectos incluso si el modelo lineal captura bien la relación en la media. Los modelos de Poisson y Binomial Negativa están diseñados exactamente para respetar estas características.

4.2 El modelo de regresión de Poisson

4.2.1 La distribución de Poisson

La distribución de Poisson con parámetro \(\lambda > 0\) asigna a cada entero \(k \geq 0\) la probabilidad:

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

Su rasgo definitorio es la equidispersión: \(E[Y] = \text{Var}[Y] = \lambda\). La distribución está completamente caracterizada por un único parámetro, de modo que media y varianza son iguales. Cuando los datos reales muestran varianza superior a la media — lo que es lo habitual — hablamos de sobredispersión y la Poisson infravalora la variabilidad.

par(mfrow=c(1,3), mar=c(4.5,4,2.2,0.5), cex.main=0.88)
for (lam in c(1, 5, 15)) {
  k_max <- qpois(0.999, lam)
  k_seq <- 0:k_max
  bp    <- barplot(dpois(k_seq, lam), col=gray(0.65), border="white",
                   names.arg=as.character(k_seq),
                   cex.names=ifelse(lam>10, 0.55, 0.75),
                   ylab="P(Y = k)", main=bquote(Poisson(lambda==.(lam))),
                   ylim=c(0, max(dpois(k_seq,lam))*1.2))
  abline(h=0)
  text(max(bp)*0.6, max(dpois(k_seq,lam))*1.05,
       sprintf("E = Var = %d", lam), cex=0.75, font=2)
}
par(mfrow=c(1,1))

Figure 4.1: Distribución de Poisson para tres valores del parámetro lambda. A medida que lambda crece, la distribución se desplaza a la derecha y se vuelve más simétrica (se aproxima a la Normal por el TCL). En todos los casos E[Y] = Var[Y] = lambda.

Menos común, pero posible es la subdispersión: \(\text{Var}(Y_i) < E[Y_i]\). Este fenómeno ocurre cuando los datos están fuertemente acotados o cuando hay mecanismos reguladores.

4.2.2 El modelo de regresión de Poisson

Para modelar datos de recuento con covariables necesitamos que el parámetro \(\lambda_i\) dependa de \(\mathbf{x}_i\). La función de enlace natural es el logaritmo, que garantiza \(\lambda_i > 0\):

\[\lambda_i = E[Y_i \mid \mathbf{x}_i] = e^{\mathbf{x}_i'\boldsymbol{\beta}}\]

o equivalentemente, en escala log-lineal: \(\ln \lambda_i = \mathbf{x}_i'\boldsymbol{\beta}\). Esta especificación hace que los efectos de las covariables sean multiplicativos sobre la tasa esperada de ocurrencias: un incremento unitario en \(x_j\) multiplica \(\lambda\) por el factor \(e^{\beta_j}\).

El parámetro \(\lambda\) es la intensidad o tasa del proceso. Si \(\lambda = 2\), en promedio ocurren 2 eventos por unidad de tiempo. Si \(\lambda = 0.5\), en promedio ocurre medio evento (o equivalentemente, 1 evento cada 2 unidades de tiempo).

Visualmente, cuando \(\lambda\) es pequeño (ej. \(\lambda = 1\)), la distribución está concentrada en valores bajos (muchos ceros). Cuando \(\lambda\) es grande (ej. \(\lambda = 10\)), la distribución se extiende más y parece más simétrica.

4.2.3 Ejemplos de cálculo

Supongamos que el número de accidentes laborales por mes en una fábrica sigue una distribución Poisson con \(\lambda = 2.5\).

Probabilidad de exactamente 3 accidentes: \[P(Y = 3) = \frac{e^{-2.5} \cdot 2.5^3}{3!} = \frac{0.0821 \cdot 15.625}{6} = 0.2138\]

Es decir, hay aproximadamente un 21.4% de probabilidad de observar exactamente 3 accidentes.

4.2.4 Función de log-verosimilitud

Con \(n\) observaciones independientes, la log-verosimilitud es:

\[\ell(\boldsymbol{\beta}) = \sum_{i=1}^n \left[ y_i \cdot \mathbf{x}_i'\boldsymbol{\beta} - e^{\mathbf{x}_i'\boldsymbol{\beta}} - \ln(y_i!)\right]\]

El tercer sumando, \(\ln(y_i!)\), no depende de \(\boldsymbol{\beta}\) y no afecta a la maximización. Las condiciones de primer orden son \(\sum_i \mathbf{x}_i(y_i - \hat{\lambda}_i) = \mathbf{0}\), que expresan que los residuos son ortogonales a las covariables — exactamente igual que en MCO, pero la solución es no lineal y requiere métodos numéricos (Newton-Raphson o IRLS).

4.2.5 Tasas de incidencia (IRR) y efectos marginales

Los coeficientes del modelo Poisson tienen dos interpretaciones complementarias, ninguna de ellas directa en términos de cambios en \(E[y]\):

Tasa de incidencia (IRR): El cociente de tasas esperadas cuando \(x_j\) aumenta en una unidad, manteniendo el resto constante, es:

\[\text{IRR}_j = \frac{E[Y \mid x_j+1]}{E[Y \mid x_j]} = \frac{e^{\mathbf{x}'\boldsymbol{\beta}+\beta_j}}{e^{\mathbf{x}'\boldsymbol{\beta}}} = e^{\beta_j}\]

Un \(\text{IRR} > 1\) indica que la covariable aumenta la tasa de recuento; un \(\text{IRR} < 1\), que la reduce. Es el equivalente del odds ratio del Logit, pero para recuentos.

Efecto marginal (AME): El efecto sobre el recuento esperado es no constante:

\[\frac{\partial E[Y \mid \mathbf{x}]}{\partial x_j} = \lambda_i \cdot \beta_j\]

El AME promedia este efecto sobre la muestra: \(\text{AME}_j = \bar\lambda \cdot \beta_j = \frac{1}{n}\sum_i \hat\lambda_i \cdot \beta_j\).

4.3 La sobredispersión y el modelo Binomial Negativa

4.3.1 El problema de la sobredispersión

En la práctica, los datos de recuento presentan casi siempre varianza mayor que la media. Esto tiene dos consecuencias graves si se ignora. La primera es que los errores estándar de Poisson están sesgados hacia abajo, produciendo estadísticos \(z\) inflados y \(p\)-valores demasiado pequeños: el modelo declara significativas variables que en realidad no lo son. La segunda es que las predicciones de la cola derecha — los recuentos muy elevados — son sistemáticamente incorrectas.

set.seed(2031)
par(mfrow=c(1,2), mar=c(4.5,4.5,1.8,0.5), cex.main=0.88)

# Panel 1: media vs varianza por grupos
grupos <- cut(rpois(500, exp(2 + 0.3*rnorm(500))), breaks=5)
lam_g  <- tapply(rpois(500, exp(2+0.3*rnorm(500))), grupos, mean)
# Simular sobredispersión: NB
y_nb   <- MASS::rnegbin(500, mu=exp(2+0.3*rnorm(500)), theta=2)
var_g  <- tapply(y_nb, cut(y_nb, breaks=5), var, na.rm=TRUE)
med_g  <- tapply(y_nb, cut(y_nb, breaks=5), mean, na.rm=TRUE)
med_g  <- med_g[!is.na(med_g)]; var_g <- var_g[!is.na(var_g)]

plot(med_g, var_g, pch=16, col=gray(0.2), cex=1.5,
     xlab="Media del grupo", ylab="Varianza del grupo",
     main="Media vs Varianza")
abline(0,1,  lwd=2, lty=1)  # línea Poisson: Var=Media
abline(lm(var_g ~ med_g), lwd=2, lty=2)  # línea datos reales
legend("topleft", legend=c("Poisson (Var=Media)","Datos reales"),
       lty=c(1,2), lwd=2, bty="n", cex=0.85)

# Panel 2: Poisson vs NB en la distribución
mu_ex  <- 5
k_seq  <- 0:25
plot(k_seq-0.15, dpois(k_seq, mu_ex), type="h", lwd=3, lty=1,
     xlab="k", ylab="P(Y=k)", main="Poisson vs Binomial Negativa")
lines(k_seq+0.15, dnbinom(k_seq, mu=mu_ex, size=2), type="h",
      lwd=3, lty=2, col=gray(0.5))
legend("topright", legend=c("Poisson(5)","NB(mu=5, theta=2)"),
       lty=c(1,2), lwd=3, col=c("black",gray(0.5)), bty="n", cex=0.85)
par(mfrow=c(1,1))

Figure 4.2: La sobredispersión en datos de recuento. Panel izquierdo: bajo equidispersión (Poisson) todos los puntos caen sobre la línea de 45°. Con sobredispersión, la varianza supera a la media en todos los grupos. Panel derecho: la distribución Binomial Negativa (línea discontinua) asigna más probabilidad a la cola derecha que Poisson (línea sólida), ajustando mejor los recuentos extremos.

El modelo quasi-Poisson surge cuando el modelo Poisson es demasiado restrictivo, ya que impone que la varianza sea igual a la media, algo que rara vez se cumple en datos reales (donde suele haber sobredispersión). El quasi-Poisson relaja este supuesto permitiendo que la varianza sea proporcional a la media mediante un parámetro de dispersión. Esto tiene una consecuencia importante: los coeficientes estimados son exactamente los mismos que en Poisson, pero los errores estándar se ajustan (normalmente aumentan), corrigiendo así la inferencia. En otras palabras, el quasi-Poisson no cambia el modelo de la media, sino que corrige cómo medimos la incertidumbre cuando los datos son más variables de lo que Poisson permite.

4.3.2 El test formal de sobredispersión

El test de Cameron y Trivedi (1990), implementado en AER::dispersiontest(), contrasta: \[H_0: \text{Var}(y) = E(y) \quad \text{vs.} \quad H_1: \text{Var}(y) = E(y) + \alpha \cdot E(y)^c\] donde \(c=2\) corresponde a la especificación Binomial Negativa. Un \(\alpha > 0\) significativo indica sobredispersión y justifica abandonar el modelo Poisson.

4.3.3 El modelo Binomial Negativa (NB2)

La distribución Binomial Negativa generaliza la Poisson añadiendo un parámetro de dispersión \(\theta > 0\):

\[P(Y=k \mid \mu, \theta) = \binom{k + \theta - 1}{k}\left(\frac{\theta}{\theta+\mu}\right)^\theta\left(\frac{\mu}{\theta+\mu}\right)^k\]

Sus momentos son \(E[Y] = \mu\) y \(\text{Var}[Y] = \mu + \mu^2/\theta\). El ratio varianza-media es:

\[\text{VMR} = 1 + \frac{\mu}{\theta}\]

Cuando \(\theta \to \infty\), VMR \(\to 1\) y la NB converge a la Poisson: el modelo Poisson es un caso límite de la NB. Cuanto más pequeño es \(\theta\), mayor es la sobredispersión. En el modelo de regresión NB2, \(\mu_i = e^{\mathbf{x}_i'\boldsymbol{\beta}}\) y \(\theta\) se estima simultáneamente con \(\boldsymbol{\beta}\) por máxima verosimilitud. En R: MASS::glm.nb(y ~ x, data=df).

4.4 Modelo ilustrativo: visitas al médico

Para el análisis completo utilizamos el dataset de visitas al médico, que ya conocemos del Capítulo 1: (T01_CP03_visitas_medico.RData, \(n=600\)). La variable visitas presenta un ratio VMR = 10.3, lo que anticipa una sobredispersión severa y hace de este dataset un caso ilustrativo ideal.

4.4.1 Análisis exploratorio y diagnóstico previo

library(kableExtra); library(MASS); library(AER)
load("data/T01_CP03_visitas_medico.RData")

par(mfrow=c(1,3), mar=c(4.5,4.5,1.8,0.5), cex.main=0.88)

# Panel 1: observado vs Poisson
mu_v  <- mean(medico$visitas)
k_lim <- min(max(medico$visitas), 25)
obs_f <- table(factor(pmin(medico$visitas,25), levels=0:25))/nrow(medico)
plot(0:25-0.2, as.numeric(obs_f), type="h", lwd=3, lty=1,
     xlab="N.º visitas", ylab="Proporción", main="Obs. vs Poisson")
lines(0:25+0.2, dpois(0:25, mu_v), type="h", lwd=2.5, lty=2, col=gray(0.55))
legend("topright", legend=c("Observado","Poisson"),
       lty=c(1,2), lwd=c(3,2.5), col=c("black",gray(0.55)), bty="n", cex=0.8)

# Panel 2: media vs varianza por grupos de edad
medico$edad_g <- cut(medico$edad, breaks=c(17,30,45,60,80), labels=c("18-30","31-45","46-60","61+"))
m_g <- tapply(medico$visitas, medico$edad_g, mean)
v_g <- tapply(medico$visitas, medico$edad_g, var)
plot(m_g, v_g, pch=16, cex=1.8, col=gray(0.2),
     xlab="Media por grupo edad", ylab="Varianza por grupo edad",
     main="Media vs Varianza")
abline(0,1, lwd=2, lty=1)
abline(lm(v_g~m_g), lwd=2, lty=2)
text(m_g, v_g+1.5, names(m_g), cex=0.75)
legend("topleft", legend=c("Poisson (Var=Med)","Datos reales"),
       lty=c(1,2), lwd=2, bty="n", cex=0.8)

# Panel 3: medias por grupos
grupos_v <- c("Sin crónico","Con crónico","Sin seg.priv.","Con seg.priv.")
medias_v <- c(mean(medico$visitas[medico$cronico==0]),
              mean(medico$visitas[medico$cronico==1]),
              mean(medico$visitas[medico$seguro_privado==0]),
              mean(medico$visitas[medico$seguro_privado==1]))
bp_v <- barplot(medias_v, col=gray(c(0.80,0.40,0.80,0.40)), border="white",
        names.arg=grupos_v, cex.names=0.72,
        ylab="Media visitas/año", main="Medias por grupo", las=2)
abline(h=0)
text(bp_v, medias_v+0.3, sprintf("%.1f",medias_v), cex=0.8, font=2)
par(mfrow=c(1,1))

Figure 4.3: Diagnóstico previo para el modelo de recuento. Panel izquierdo: distribución observada vs Poisson teórica. La distribución real tiene cola más larga que Poisson, confirmando sobredispersión. Panel central: media y varianza por grupos de edad — la varianza crece más que la media (patrón NB). Panel derecho: diferencias de media de visitas entre grupos de cronicidad y seguro privado.

4.4.2 Estimación: Poisson, Quasi-Poisson y Binomial Negativa

fml_v <- visitas ~ edad + cronico + ingreso + seguro_privado + limitacion_funcional + sexo

pois_v <- glm(fml_v, data=medico, family=poisson(link="log"))
qpois_v<- glm(fml_v, data=medico, family=quasipoisson(link="log"))
nb_v   <- glm.nb(fml_v, data=medico)

vars_v <- c("edad","cronico","ingreso","seguro_privado","limitacion_funcional","sexo")
labs_v <- c("Edad","Crónico (=1)","Ingreso","Seg. privado (=1)","Limitación func.","Mujer (=1)")

fmt_coef_count <- function(mod, family="poisson") {
  ct <- coef(summary(mod))
  pname <- if (family=="nb") "Pr(>|z|)" else "Pr(>|z|)"
  b  <- ct[vars_v,"Estimate"]
  pv <- ct[vars_v, ncol(ct)]
  sig <- ifelse(pv<0.001,"***",ifelse(pv<0.01,"**",ifelse(pv<0.05,"*","")))
  paste0(sprintf("%.3f",b), sig)
}

df_est_v <- data.frame(
  Variable = labs_v,
  Poisson  = fmt_coef_count(pois_v),
  `EE Pois.` = sprintf("%.3f", coef(summary(pois_v))[vars_v,"Std. Error"]),
  `EE Q-Pois.` = sprintf("%.3f", coef(summary(qpois_v))[vars_v,"Std. Error"]),
  NB       = fmt_coef_count(nb_v, "nb"),
  `EE NB` = sprintf("%.3f", coef(summary(nb_v))[vars_v,"Std. Error"]),
  check.names=FALSE
)
kbl(df_est_v, booktabs=TRUE, row.names=FALSE, align=c("l","r","r","r","r","r"),
    col.names=c("Variable","Poisson","EE Pois.","EE Q-Pois.","NB","EE NB")) |>
  kable_styling(latex_options=c("hold_position","scale_down"),
                bootstrap_options=c("condensed","striped"), full_width=FALSE) |>
  footnote(general=sprintf("***p<.001, **p<.01, *p<.05. n=%d. VMR observada=%.1f.",
                           nrow(medico), var(medico$visitas)/mean(medico$visitas)),
           general_title="Nota:", footnote_as_chunk=FALSE)

Table 4.1: Estimación comparada Poisson, Quasi-Poisson y Binomial Negativa para el número de visitas al médico. Los coeficientes son idénticos en Poisson y Quasi-Poisson, pero los errores estándar del Quasi-Poisson son mayores (corrigen la sobredispersión). La NB da coeficientes similares pero con errores estándar aún diferentes.

Variable	Poisson	EE Pois.	EE Q-Pois.	NB	EE NB
Edad	0.016***	0.001	0.002	0.019***	0.002
Crónico (=1)	1.151***	0.035	0.068	1.162***	0.073
Ingreso	-0.007***	0.001	0.001	-0.007***	0.001
Seg. privado (=1)	0.218***	0.031	0.060	0.290***	0.069
Limitación func.	0.399***	0.032	0.063	0.452***	0.081
Mujer (=1)	0.196***	0.029	0.058	0.213**	0.065
Nota:
***p<.001, **p<.01, *p<.05. n=600. VMR observada=10.3.

La Table 4.1 ilustra el efecto de la sobredispersión sobre la inferencia. Los coeficientes de Poisson y Quasi-Poisson son idénticos — ambos maximizan la misma ecuación de estimación —, pero los errores estándar del Quasi-Poisson son un factor \(\sqrt{\hat\phi}\) mayores, donde \(\hat\phi\) es el parámetro de dispersión estimado (3.83). Esto hace que algunos coeficientes que parecen significativos con Poisson dejen de serlo con el Quasi-Poisson.

4.4.3 Test de sobredispersión y comparación de modelos

disp_test <- dispersiontest(pois_v, trafo=2)

df_comp <- data.frame(
  Modelo   = c("Poisson","Quasi-Poisson","Binomial Negativa"),
  LogLik   = c(sprintf("%.1f", as.numeric(logLik(pois_v))),
               "—",
               sprintf("%.1f", as.numeric(logLik(nb_v)))),
  AIC      = c(sprintf("%.1f", AIC(pois_v)), "—", sprintf("%.1f", AIC(nb_v))),
  `Phi/theta` = c("1 (fijo)",
                  sprintf("%.2f (phi)", summary(qpois_v)$dispersion),
                  sprintf("%.3f (theta)", nb_v$theta)),
  check.names=FALSE
)
kbl(df_comp, booktabs=TRUE, row.names=FALSE, align=c("l","r","r","r"),
    col.names=c("Modelo","Log-Verosim.","AIC","Dispersión")) |>
  kable_styling(latex_options=c("hold_position","scale_down"),
                bootstrap_options=c("condensed","striped"), full_width=FALSE) |>
  footnote(general=sprintf("Test Cameron-Trivedi: alpha=%.3f, z=%.2f, p<0.001. NB preferida.",
                           disp_test$estimate, disp_test$statistic),
           general_title="Nota:", footnote_as_chunk=FALSE)

Table 4.2: Test de sobredispersión (Cameron-Trivedi) y comparación de modelos Poisson vs NB. El test rechaza con claridad la equidispersión. La NB mejora sustancialmente el ajuste (AIC mucho menor) y el parámetro theta estimado es pequeño, confirmando la presencia de sobredispersión.

Modelo	Log-Verosim.	AIC	Dispersión
Poisson	-2074.1	4162.2	1 (fijo)
Quasi-Poisson	—	—	3.83 (phi)
Binomial Negativa	-1632.9	3281.7	2.360 (theta)
Nota:
Test Cameron-Trivedi: alpha=0.254, z=9.18, p<0.001. NB preferida.

4.4.4 IRR y efectos marginales del modelo NB

irr_v  <- exp(coef(nb_v)[vars_v])
ci_nb  <- suppressMessages(confint(nb_v, level=0.95))
irr_lo <- exp(ci_nb[vars_v,1])
irr_hi <- exp(ci_nb[vars_v,2])

lambda_hat <- fitted(nb_v)
ame_nb     <- mean(lambda_hat) * coef(nb_v)[vars_v]

df_irr <- data.frame(
  Variable = labs_v,
  IRR      = sprintf("%.3f", irr_v),
  `IC inf.95%` = sprintf("%.3f", irr_lo),
  `IC sup.95%` = sprintf("%.3f", irr_hi),
  `AME (vis./año)` = sprintf("%.3f", ame_nb),
  check.names=FALSE
)
kbl(df_irr, booktabs=TRUE, row.names=FALSE, align=c("l","r","r","r","r"),
    col.names=c("Variable","IRR","IC 95% inf.","IC 95% sup.","AME")) |>
  kable_styling(latex_options=c("hold_position","scale_down"),
                bootstrap_options=c("condensed","striped"), full_width=FALSE) |>
  footnote(general="IC al 95% con perfil de verosimilitud.",
           general_title="Nota:", footnote_as_chunk=FALSE)

Table 4.3: Tasas de incidencia (IRR) y efectos marginales AME del modelo Binomial Negativa. El IRR del cronico indica que padecer una enfermedad crónica multiplica el número esperado de visitas por 2.49. El AME expresa el incremento en visitas en unidades absolutas.

Variable	IRR	IC 95% inf.	IC 95% sup.	AME
Edad	1.019	1.016	1.023	0.152
Crónico (=1)	3.196	2.777	3.681	9.302
Ingreso	0.993	0.990	0.995	-0.058
Seg. privado (=1)	1.336	1.166	1.531	2.322
Limitación func.	1.571	1.340	1.846	3.619
Mujer (=1)	1.238	1.089	1.407	1.708
Nota:
IC al 95% con perfil de verosimilitud.

La Table 4.3 sintetiza las dos interpretaciones clave. El IRR del cronico es 3.20: padecer una enfermedad crónica multiplica el número esperado de visitas por 3.20 — es decir, un 220% más. En términos del AME, esto equivale a 9.3 visitas adicionales al año respecto a un individuo sin enfermedad crónica. El seguro privado aumenta las visitas en un factor de 1.34 (AME = 2.32 visitas más).

4.5 Extensiones: modelos inflados en cero

Hay situaciones en las que la proporción de ceros en los datos supera incluso lo que predice el modelo NB. Esto ocurre cuando el proceso generador de datos tiene dos regímenes: un régimen que siempre produce cero (individuos que por construcción no pueden tener el evento) y un régimen que produce recuentos según un Poisson o NB. Los modelos inflados en cero (ZIP y ZINB) modelizan explícitamente estos dos regímenes con una mezcla:

\[P(Y=0) = \pi + (1-\pi)\,g(0) \qquad P(Y=k>0) = (1-\pi)\,g(k)\]

donde \(\pi\) es la probabilidad de pertenecer al régimen de ceros estructurales y \(g(\cdot)\) es la distribución de recuento (Poisson o NB) del régimen activo. En R, pscl::zeroinfl() estima estos modelos. Su uso está justificado cuando el test de Vuong y el gráfico de frecuencias muestran un exceso de ceros que NB no puede absorber. El desarrollo detallado de estos modelos queda fuera del alcance de este capítulo introductorio.

4.6 Casos prácticos

4.6.1 Caso práctico 1: Delitos por municipio

Contexto. El dataset CP01 recoge 400 municipios con variables socioeconómicas y de seguridad. La variable delitos tiene media 8.8 y VMR = 6.8, indicando sobredispersión moderada donde la NB supera al Poisson.

load("data/T04_CP01_criminalidad.RData")
fml_c <- delitos ~ densidad_pob + tasa_desempleo + ingresos_medios + policia_per_capita + zona_urbana
nb_c  <- glm.nb(fml_c, data=crimen)
p_c   <- glm(fml_c, data=crimen, family=poisson)
dt_c  <- dispersiontest(p_c, trafo=2)

vars_c <- c("densidad_pob","tasa_desempleo","ingresos_medios","policia_per_capita","zona_urbana")
labs_c <- c("Densidad pob.","Tasa desempleo (%)","Ingresos medios (€m)","Policías/cápita","Zona urbana")

ct_c   <- coef(summary(nb_c))
irr_c  <- exp(coef(nb_c)[vars_c])
lam_c  <- fitted(nb_c)
ame_c  <- mean(lam_c) * coef(nb_c)[vars_c]

df_c1 <- data.frame(
  Variable = labs_c,
  Coef = paste0(sprintf("%.4f",coef(nb_c)[vars_c]),
         ifelse(ct_c[vars_c,4]<0.001,"***",ifelse(ct_c[vars_c,4]<0.01,"**",
                ifelse(ct_c[vars_c,4]<0.05,"*","")))),
  IRR  = sprintf("%.3f", irr_c),
  AME  = sprintf("%.3f", ame_c),
  check.names=FALSE
)
kbl(df_c1, booktabs=TRUE, row.names=FALSE, align=c("l","r","r","r"),
    col.names=c("Variable","Coef. NB","IRR","AME")) |>
  kable_styling(latex_options=c("hold_position","scale_down"),
                bootstrap_options=c("condensed","striped"), full_width=FALSE) |>
  footnote(general=sprintf("n=%d | VMR=%.1f | theta=%.3f | AIC NB=%.1f vs AIC Pois=%.1f.",
                           nrow(crimen), var(crimen$delitos)/mean(crimen$delitos),
                           nb_c$theta, AIC(nb_c), AIC(p_c)),
           general_title="Nota:", footnote_as_chunk=FALSE)

Table 4.4: Modelo NB para delitos municipales (CP01). El desempleo y la densidad de población aumentan la criminalidad; ingresos más altos y mayor presencia policial la reducen. VMR=6.8 justifica la NB sobre Poisson.

Variable	Coef. NB	IRR	AME
Densidad pob.	0.0006***	1.001	0.005
Tasa desempleo (%)	0.0525***	1.054	0.463
Ingresos medios (€m)	-0.0285***	0.972	-0.251
Policías/cápita	-0.0994***	0.905	-0.877
Zona urbana	0.4057***	1.500	3.579
Nota:
n=400 | VMR=6.8 | theta=3.530 | AIC NB=2358.8 vs AIC Pois=2853.5.

El alumno puede reproducir este análisis con el siguiente script:

T04_CP01_mECO_Poisson_Criminalidad.R

4.6.2 Caso práctico 2: Accidentes de tráfico por tramo

Contexto. El dataset CP02 recoge 500 tramos de carretera. La variable accidentes tiene media 4 y VMR = 3.2. Dado el 15% de tramos sin ningún accidente, se aplica la NB y se compara con el Poisson.

par(mfrow=c(1,2), mar=c(4.5,4.5,1.8,0.5), cex.main=0.88)

tab_ac <- table(factor(pmin(accid$accidentes,15), levels=0:15))
bp_ac  <- barplot(tab_ac, col=gray(0.65), border="white",
                  names.arg=c(as.character(0:14),"15+"),
                  xlab="N.º accidentes", ylab="Frecuencia",
                  main="Distribución accidentes")
abline(h=0)

boxplot(accidentes ~ iluminacion, data=accid,
        names=c("Sin ilum.","Con ilum."),
        col=gray(c(0.80,0.45)), border=gray(0.2),
        xlab="", ylab="N.º accidentes", main="Accidentes por iluminación")
par(mfrow=c(1,1))

Figure 4.4: Accidentes de tráfico (CP02). Panel izquierdo: la distribución muestra fuerte concentración en 0 y cola derecha larga. Panel derecho: la iluminación del tramo reduce los accidentes de forma muy apreciable.

fml_a  <- accidentes ~ trafico_diario + velocidad_max + num_curvas + dias_lluvia + iluminacion + longitud_km
p_a    <- glm(fml_a, data=accid, family=poisson)
nb_a   <- glm.nb(fml_a, data=accid)
dt_a   <- dispersiontest(p_a, trafo=2)

vars_a <- c("trafico_diario","velocidad_max","num_curvas","dias_lluvia","iluminacion","longitud_km")
labs_a <- c("Tráfico diario","Vel. máx. (km/h)","N.º curvas","Días lluvia/año","Iluminación (=1)","Longitud (km)")

ct_a  <- coef(summary(nb_a))
irr_a <- exp(coef(nb_a)[vars_a])
ame_a <- mean(fitted(nb_a)) * coef(nb_a)[vars_a]

df_a2 <- data.frame(
  Variable = labs_a,
  `Pois.` = sprintf("%.4f",coef(p_a)[vars_a]),
  `NB` = paste0(sprintf("%.4f",coef(nb_a)[vars_a]),
         ifelse(ct_a[vars_a,4]<0.001,"***",ifelse(ct_a[vars_a,4]<0.01,"**",
                ifelse(ct_a[vars_a,4]<0.05,"*","")))),
  IRR = sprintf("%.3f", irr_a),
  check.names=FALSE
)
kbl(df_a2, booktabs=TRUE, row.names=FALSE, align=c("l","r","r","r"),
    col.names=c("Variable","Poisson","NB","IRR")) |>
  kable_styling(latex_options=c("hold_position","scale_down"),
                bootstrap_options=c("condensed","striped"), full_width=FALSE) |>
  footnote(general=sprintf("n=%d | VMR=%.1f | Cameron-Trivedi p<.001 | AIC NB=%.1f vs Pois=%.1f.",
                           nrow(accid), var(accid$accidentes)/mean(accid$accidentes),
                           AIC(nb_a), AIC(p_a)),
           general_title="Nota:", footnote_as_chunk=FALSE)

Table 4.5: Modelos Poisson y NB para accidentes de tráfico (CP02). La iluminación y la velocidad máxima son los predictores más potentes. El test Cameron-Trivedi confirma sobredispersión y la NB es claramente preferida por AIC.

Variable	Poisson	NB	IRR
Tráfico diario	0.0000	0.0000	1.000
Vel. máx. (km/h)	0.0061	0.0069***	1.007
N.º curvas	0.1141	0.1255***	1.134
Días lluvia/año	0.0038	0.0043**	1.004
Iluminación (=1)	-0.2247	-0.2405**	0.786
Longitud (km)	0.0290	0.0328***	1.033
Nota:
n=500 | VMR=3.2 | Cameron-Trivedi p<.001 | AIC NB=2427.7 vs Pois=2802.4.

El alumno puede reproducir este análisis con el siguiente script:

T04_CP02_mECO_NB_Accidentes.R

4.6.3 Caso práctico 3: Publicaciones académicas

Contexto. El dataset CP03 recoge 600 investigadores. La variable publicaciones tiene media 8 y VMR = 26.9 — sobredispersión muy severa donde la NB es imprescindible. Los predictores son experiencia, financiación, doctorado, h-index previo, sexo y número de colaboradores.

fml_i  <- publicaciones ~ experiencia + financiacion + doctorado + h_index_previo + sexo + colaboradores
p_i    <- glm(fml_i, data=invest, family=poisson)
nb_i   <- glm.nb(fml_i, data=invest)

vars_i <- c("experiencia","financiacion","doctorado","h_index_previo","sexo","colaboradores")
labs_i <- c("Experiencia (años)","Financiación (€m)","Doctorado (=1)","H-index previo","Mujer (=1)","Colaboradores")

ct_i  <- coef(summary(nb_i)); ct_p <- coef(summary(p_i))
irr_i <- exp(coef(nb_i)[vars_i])
ame_i <- mean(fitted(nb_i)) * coef(nb_i)[vars_i]

sig_mark <- function(pv) ifelse(pv<0.001,"***",ifelse(pv<0.01,"**",ifelse(pv<0.05,"*","")))
df_i3 <- data.frame(
  Variable = labs_i,
  `Pois. EE`  = sprintf("%.3f",ct_p[vars_i,"Std. Error"]),
  `NB coef.`  = paste0(sprintf("%.3f",coef(nb_i)[vars_i]), sig_mark(ct_i[vars_i,4])),
  `NB EE`     = sprintf("%.3f",ct_i[vars_i,"Std. Error"]),
  IRR         = sprintf("%.3f",irr_i),
  check.names=FALSE
)
kbl(df_i3, booktabs=TRUE, row.names=FALSE, align=c("l","r","r","r","r"),
    col.names=c("Variable","EE Poisson","Coef. NB","EE NB","IRR")) |>
  kable_styling(latex_options=c("hold_position","scale_down"),
                bootstrap_options=c("condensed","striped"), full_width=FALSE) |>
  footnote(general=sprintf("n=%d | VMR=%.1f | theta NB=%.3f. EE Poisson << EE NB: inflación de significatividad.",
                           nrow(invest), var(invest$publicaciones)/mean(invest$publicaciones), nb_i$theta),
           general_title="Nota:", footnote_as_chunk=FALSE)

Table 4.6: Modelo NB para publicaciones académicas (CP03). El doctorado y el h-index previo son los predictores más potentes. Con VMR=26.9, el Poisson infraestima gravemente los errores estándar, declarando significativas variables que con NB no lo son.

Variable	EE Poisson	Coef. NB	EE NB	IRR
Experiencia (años)	0.002	0.055***	0.006	1.057
Financiación (€m)	0.000	0.006***	0.001	1.006
Doctorado (=1)	0.037	0.592***	0.135	1.807
H-index previo	0.007	0.076*	0.031	1.079
Mujer (=1)	0.029	0.158	0.125	1.172
Colaboradores	0.008	0.153***	0.036	1.166
Nota:
n=600 | VMR=26.9 | theta NB=0.482. EE Poisson << EE NB: inflación de significatividad.

La Table 4.6 pone de manifiesto de forma contundente el riesgo de ignorar la sobredispersión. Los errores estándar de Poisson son sistemáticamente menores que los de NB: por ejemplo, el error estándar de experiencia pasa de 0.002 (Poisson) a 0.006 (NB), una diferencia que puede cambiar las conclusiones sobre significatividad. El alumno puede reproducir este análisis con el siguiente script:

T04_CP03_mECO_NB_Publicaciones.R

4.7 Guía de comandos R para este capítulo

4.7.1 glm(…, family = poisson) — Regresión de Poisson

La regresión de Poisson se estima con glm() especificando family = poisson(link = "log"). La función de enlace logarítmica garantiza que la tasa esperada \(\lambda_i = E[y_i|\mathbf{x}_i]\) sea siempre positiva. Los coeficientes estimados son efectos sobre el log de la tasa, no sobre la tasa directamente.

pois <- glm(visitas ~ edad + cronico + ingreso + seguro_privado,
            data = medico, family = poisson(link = "log"))
summary(pois)

Los parámetros principales son: la fórmula especifica la variable de recuento y las covariables; family = poisson(link = "log") activa el modelo Poisson con enlace logarítmico (es el enlace canónico y el más habitual); data es el dataframe. summary(pois) muestra los coeficientes, errores estándar, estadísticos \(z\) y \(p\)-valores.

4.7.2 AER::dispersiontest() — Test de sobredispersión de Cameron-Trivedi

Antes de concluir que el modelo Poisson es adecuado, es obligatorio contrastar formalmente si existe sobredispersión. dispersiontest() del paquete AER aplica el test de Cameron y Trivedi (1990), que contrasta \(H_0: \text{Var}(y) = E(y)\) frente a \(H_1: \text{Var}(y) = E(y) + \alpha \cdot E(y)^c\).

library(AER)
dt <- dispersiontest(pois, trafo = 2)
dt$estimate    # estimación de alpha
dt$statistic   # estadístico z del test
dt$p.value     # p-valor

El argumento trafo = 2 especifica \(c = 2\), que corresponde a la especificación Binomial Negativa (NB2) más habitual. Un \(\hat{\alpha} > 0\) con \(p\)-valor pequeño confirma sobredispersión y recomienda abandonar el Poisson en favor de la NB o el Quasi-Poisson.

4.7.3 MASS::glm.nb() — Regresión Binomial Negativa

Cuando el test de sobredispersión confirma que la varianza supera a la media, el modelo Binomial Negativa (NB2) es el sucesor natural del Poisson. glm.nb() del paquete MASS estima por MLE tanto los coeficientes \(\boldsymbol{\beta}\) como el parámetro de dispersión \(\theta\) (\(\text{VMR} = 1 + \mu/\theta\)).

library(MASS)
nb <- glm.nb(visitas ~ edad + cronico + ingreso + seguro_privado,
             data = medico)
nb$theta        # parámetro de dispersión estimado
AIC(pois, nb)   # comparar con Poisson por AIC

La sintaxis es idéntica a glm() salvo que no hay argumento family: glm.nb() siempre usa la distribución Binomial Negativa. nb$theta devuelve el parámetro \(\theta\) estimado: valores pequeños (< 1) indican sobredispersión severa. AIC() compara ambos modelos en escala penalizada por complejidad.

4.7.4 exp(coef()) — Tasas de Incidencia (IRR)

Los coeficientes de los modelos Poisson y NB son efectos sobre el logaritmo de la tasa esperada \(\ln \lambda\). Para obtener las tasas de incidencia (IRR), que tienen una interpretación multiplicativa directa sobre el recuento esperado, se aplica la función exponencial.

exp(coef(nb))                       # IRR puntuales
exp(confint(nb, level = 0.95))      # IC al 95% de los IRR

Un IRR de 2.5 para cronico significa que padecer una enfermedad crónica multiplica el número esperado de visitas por 2.5 (es decir, un 150% más que los no crónicos, manteniendo el resto constante). Los IRR son equivalentes a los odds ratios del Logit, pero para recuentos.

4.7.5 mean(fitted(nb)) * coef(nb) — Efectos marginales AME

Los IRR miden efectos multiplicativos; los efectos marginales AME miden el cambio absoluto en el recuento esperado cuando una covariable aumenta en una unidad. El AME es el producto del coeficiente por la media de los recuentos predichos.

# AME del modelo NB sobre E[y|x]
ame_nb <- mean(fitted(nb)) * coef(nb)["cronico"]

fitted(nb) devuelve los recuentos predichos \(\hat{\lambda}_i = e^{\mathbf{x}_i'\hat{\boldsymbol{\beta}}}\) para cada observación. El AME es, por tanto, la media de esas predicciones multiplicada por el coeficiente. Para el AME de las variables continuas, la interpretación es: un aumento unitario en \(x_j\) aumenta el recuento esperado en ame_nb unidades.

4.7.6 AIC() — Comparación de modelos por criterio de información

Para comparar modelos con distinto número de parámetros (como Poisson vs NB) se usa el criterio de información de Akaike (AIC), que penaliza la complejidad del modelo. El modelo con menor AIC es preferido.

AIC(pois)         # AIC del modelo Poisson
AIC(nb)           # AIC del modelo NB
AIC(pois) - AIC(nb)   # diferencia: positivo → NB mejor

\(\text{AIC} = -2\ell(\hat{\boldsymbol{\beta}}) + 2k\), donde \(\ell\) es la log-verosimilitud y \(k\) el número de parámetros. El Quasi-Poisson no tiene AIC porque no es un modelo MLE completo; en ese caso, la comparación debe hacerse mediante el test F o simplemente observando si el parámetro de dispersión \(\hat{\phi}\) es claramente distinto de 1.

Los scripts de este capítulo están disponibles en la carpeta scripts/ del repositorio del manual, en:

github.com/carlanta/MicroEconometrics

Cada script cubre el ciclo completo de análisis: EDA con diagnóstico VMR, test de sobredispersión Cameron-Trivedi, estimación Poisson y NB, comparación de errores estándar, IRR, AME y comparación por AIC.

Cameron, A. Colin, and Pravin K. Trivedi. 1990. “Regression-Based Tests for Overdispersion in the Poisson Model.” Journal of Econometrics 46: 347–64.