# --- Instalar (solo los que falten) ------------------------------------
paquetes <- c("tidyverse", # ggplot2, dplyr, tidyr, purrr, tibble
"GLMsData", # datos del libro: danishlc, kstones, pock, hcrabs ...
"statmod", # qresid(): residuos cuantílicos aleatorizados
"broom", # tidy(): tablas de coeficientes
"patchwork", # combinar gráficas
"MASS", # glm.nb() (se usa siempre como MASS::)
"sandwich", "lmtest") # errores estándar robustos (opcional)
faltan <- setdiff(paquetes, rownames(installed.packages()))
if (length(faltan) > 0) install.packages(faltan)
# --- Cargar -------------------------------------------------------------
library(tidyverse)
library(GLMsData)
library(statmod)
library(broom)
library(patchwork)
theme_set(theme_minimal(base_size = 14))
col1 <- "#2BA8E0"; col2 <- "#C45BC4"; col3 <- "#E8A33D"
set.seed(2026)Modelos para conteos: GLM Poisson y binomial negativa
Dunn & Smyth (2018), Capítulo 10 — con R
Datos: cáncer de pulmón en cuatro ciudades danesas
data(danishlc)
dlc <- as_tibble(danishlc) |>
mutate(Age = factor(Age, levels = unique(as.character(Age))), # orden natural de la tabla (40-54, ..., >74)
Tasa = Cases / Pop * 1000, # casos por 1000
AgeNum = c(40, 55, 60, 65, 70, 75)[as.integer(Age)]) # límite inferior de la clase
dlc# A tibble: 24 × 6
Cases Pop Age City Tasa AgeNum
<int> <int> <fct> <fct> <dbl> <dbl>
1 11 3059 40-54 Fredericia 3.60 40
2 11 800 55-59 Fredericia 13.8 55
3 11 710 60-64 Fredericia 15.5 60
4 10 581 65-69 Fredericia 17.2 65
5 11 509 70-74 Fredericia 21.6 70
6 10 605 >74 Fredericia 16.5 75
7 13 2879 40-54 Horsens 4.52 40
8 6 1083 55-59 Horsens 5.54 55
9 15 923 60-64 Horsens 16.3 60
10 10 834 65-69 Horsens 12.0 65
# ℹ 14 more rowsggplot(dlc, aes(Age, Tasa, colour = City, group = City, shape = City)) +
geom_line(linewidth = 1) + geom_point(size = 3) +
labs(x = "Grupo de edad", y = "Casos por 1000", colour = NULL, shape = NULL)
La tasa crece con la edad y cae en el grupo \(>74\); no hay patrón claro por ciudad.
Diagnósticos: factor edad vs. cuadrático
diag_glm(m2, "Edad como factor (m2)") / diag_glm(m4, "Edad cuadrática (m4)")
m4 muestra observaciones más influyentes y es menos simple: preferimos m2.
Paradoja de Simpson
tasas <- function(d, por) d |>
summarise(n = sum(Counts), .by = all_of(c("Method", "Outcome", por))) |>
mutate(prop = n / sum(n), .by = all_of(c("Method", por))) |>
filter(Outcome == "Success")
bind_rows(tasas(ks, character(0)) |> mutate(Size = "Todos (colapsado)"), tasas(ks, "Size")) |>
mutate(Size = factor(Size, levels = c("Small", "Large", "Todos (colapsado)"))) |>
ggplot(aes(Method, prop, fill = Method)) +
geom_col(width = .6) + geom_text(aes(label = scales::percent(prop, accuracy = .1)), vjust = -.4) +
facet_wrap(~Size) + scale_fill_manual(values = c(col1, col3), guide = "none") +
scale_y_continuous(labels = scales::percent, limits = c(0, 1)) +
labs(x = "Método", y = "Tasa de éxito")
Colapsado, el método B parece mejor (83 % vs. 78 %); por tamaño, el método A gana en ambos estratos. A se usa más con cálculos grandes (más difíciles).
¿Qué significa la triple interacción?
dyouth |>
summarise(n = sum(Obs), .by = c(Age, Gender, Depression)) |>
mutate(prop = n / sum(n), .by = c(Age, Gender)) |>
filter(Depression == "H") |>
ggplot(aes(Age, prop, colour = Gender, group = Gender)) +
geom_line(linewidth = 1.1) + geom_point(size = 3) +
scale_colour_manual(values = c(F = col2, M = col1)) +
scale_y_continuous(labels = scales::percent) +
labs(x = "Edad", y = "% con depresión alta", colour = "Género")
El patrón por edad es distinto en mujeres y varones. La interacción Age:Group solo refleja distinta composición por edad entre grupos; se conserva para que las pruebas con Depression queden ajustadas por esa demografía.
Binomial negativa como mezcla Poisson–gamma
\[ y_i\mid\lambda_i\sim\text{Pois}(\lambda_i),\quad\lambda_i\sim\text{Gamma}(\text{media }\mu_i,\ \text{CV}^2=\psi) \ \Longrightarrow\ E[y_i]=\mu_i,\ \ \text{var}[y_i]=\mu_i+\psi\mu_i^2=\mu_i+\frac{\mu_i^2}{k},\ \ k=1/\psi \]
Para \(k\) fijo es una familia exponencial (desviación unitaria \(2[y\log\frac y\mu-(y+k)\log\frac{y+k}{\mu+k}]\), \(\phi=1\)). En R \(k\) se estima junto con \(\beta\) con MASS::glm.nb() (ahí se llama theta).
mu <- 10; k <- 2
lambda <- rgamma(1e5, shape = k, rate = k / mu) # media mu, CV² = 1/k
y_mezcla <- rpois(1e5, lambda)
c(media = mean(y_mezcla), varianza = var(y_mezcla), var_teorica = mu + mu^2 / k) media varianza var_teorica
9.96351 59.63891 60.00000 tibble(y = y_mezcla) |> count(y) |> mutate(frec = n / sum(n)) |>
ggplot(aes(y, frec)) + geom_col(fill = col1, alpha = .6) +
geom_line(data = tibble(y = 0:50), aes(y, dnbinom(y, mu = mu, size = k)), colour = col2, linewidth = 1) +
geom_line(data = tibble(y = 0:50), aes(y, dpois(y, mu)), colour = col3, linewidth = 1, linetype = 2) +
labs(x = "y", y = "Frecuencia", subtitle = "Barras: Poisson–gamma simulada · morado: NB teórica · naranja: Poisson")
Datos pock: marcas en membranas
Marcas (pocks) contadas a distintas diluciones de un medio viral. Usamos \(\log_2(\text{Dilución})\).
data(pock)
pk <- as_tibble(pock) |> mutate(ld = log2(Dilution))
resumen <- pk |> summarise(media = mean(Count), varianza = var(Count), .by = ld) |>
mutate(razon = varianza / media) |> arrange(ld)
resumen# A tibble: 5 × 4
ld media varianza razon
<dbl> <dbl> <dbl> <dbl>
1 0 187. 1781. 9.54
2 1 105. 667. 6.37
3 2 50.8 360. 7.09
4 3 27.1 195. 7.20
5 4 9.1 17.7 1.94p1 <- ggplot(pk, aes(ld, Count)) + geom_jitter(width = .08, height = 0, alpha = .7) +
labs(x = "log2(dilución)", y = "Conteo de marcas")
p2 <- ggplot(resumen, aes(log(media), log(varianza))) + geom_point(size = 3) +
geom_smooth(formula = y ~ x, method = "lm", se = FALSE, colour = col1) +
labs(x = "log(media del grupo)", y = "log(varianza del grupo)")
p1 + p2
coef(lm(log(varianza) ~ log(media), data = resumen)) # pendiente ≈ 1.4: V(μ) ≈ μ^1.5(Intercept) log(media)
0.02861162 1.44318666 La varianza es mucho mayor que la media y la razón crece con la media: entre Poisson (\(V=\mu\)) y gamma (\(V=\mu^2\)).
Interpretación y ajuste
nuevo <- tibble(ld = seq(0, 4, length.out = 100))
pred <- \(m, nm, crit, ...) {
p <- predict(m, newdata = nuevo, se.fit = TRUE, ...)
tibble(modelo = nm, ld = nuevo$ld, mu = exp(p$fit),
lo = exp(p$fit - crit * p$se.fit), hi = exp(p$fit + crit * p$se.fit))
}
curvas <- bind_rows(
pred(pk_pois, "Poisson", qnorm(.975)),
pred(pk_qp, "Cuasi-Poisson", qt(.975, df.residual(pk_qp))),
pred(pk_nb, "Bin. negativa", qnorm(.975), dispersion = 1)) |>
mutate(modelo = factor(modelo, levels = c("Poisson", "Cuasi-Poisson", "Bin. negativa")))
ggplot() +
geom_ribbon(data = curvas, aes(ld, ymin = lo, ymax = hi), alpha = .25, fill = col1) +
geom_line(data = curvas, aes(ld, mu), colour = col1, linewidth = 1) +
geom_point(data = pk, aes(ld, Count), alpha = .6) +
facet_wrap(~modelo) + labs(x = "log2(dilución)", y = "Conteo de marcas", subtitle = "IC del 95 % para μ")
Cada dilución al doble multiplica el conteo esperado por \(e^{-0.7}\approx0.5\): el número de marcas es proporcional a la concentración del virus. Poisson da bandas engañosamente estrechas.
Diagnósticos de los modelos para pock
diag_glm(pk_nb, "Binomial negativa", cuantil = TRUE) /
diag_glm(pk_qp, "Cuasi-Poisson (desviación estandarizada)", cuantil = FALSE)
La NB subestima un poco la varianza de conteos pequeños; la cuasi-Poisson, la de conteos grandes (distinta función de varianza).
(Opcional) Prueba de sobredispersión con DHARMa
sim <- DHARMa::simulateResiduals(pk_pois, n = 500)
DHARMa::testDispersion(sim)
DHARMa nonparametric dispersion test via sd of residuals fitted vs.
simulated
data: simulationOutput
dispersion = 7.4117, p-value < 2.2e-16
alternative hypothesis: two.sidedExploración
hc <- hc |> mutate(logSat = log(Sat + 1))
g1 <- ggplot(hc, aes(Wt, Sat)) + geom_jitter(height = .2, alpha = .6)
g2 <- ggplot(hc, aes(log(Wt), logSat)) + geom_jitter(height = .05, alpha = .6) + labs(y = "log(Sat+1)")
g3 <- ggplot(hc, aes(Col, logSat)) + geom_boxplot(fill = col1, alpha = .5) + labs(y = "log(Sat+1)")
g4 <- ggplot(hc, aes(Width, Sat)) + geom_jitter(height = .2, alpha = .6)
g5 <- ggplot(hc, aes(log(Width), logSat)) + geom_jitter(height = .05, alpha = .6) + labs(y = "log(Sat+1)")
g6 <- ggplot(hc, aes(Spine, logSat)) + geom_boxplot(fill = col2, alpha = .5) + labs(y = "log(Sat+1)")
(g1 + g2 + g3) / (g4 + g5 + g6)
Más satélites con hembras más grandes, más claras y sin daño en la espina.
Las covariables están relacionadas
coef(lm(log(Wt) ~ log(Width), data = hc)) # pendiente ≈ 2.6 (≈ 3 si fueran todas de la misma forma)(Intercept) log(Width)
-0.6004161 2.5629298 ggplot(hc, aes(log(Width), log(Wt), colour = Col)) + geom_point(alpha = .7) +
scale_colour_manual(values = c("#F4E3A1", col3, "#8C5A2B", "#2D1B0E")) +
labs(colour = "Color")
El peso resume el tamaño: las hembras más claras y sin daño tienden a pesar más. Color y espina podrían ser efectos aparentes vía el peso.
Cuasi-Poisson vs. binomial negativa: ajuste e incertidumbre
nuevoW <- tibble(Wt = seq(min(hc$Wt), max(hc$Wt), length.out = 100))
pr <- \(m, nm, crit, ...) {
p <- predict(m, newdata = nuevoW, se.fit = TRUE, ...)
tibble(modelo = nm, lWt = log(nuevoW$Wt), mu = exp(p$fit),
lo = exp(p$fit - crit * p$se.fit), hi = exp(p$fit + crit * p$se.fit))
}
curv <- bind_rows(pr(cr_m2, "Cuasi-Poisson", qt(.975, df.residual(cr_m2))),
pr(cr_nb, "Bin. negativa", qnorm(.975), dispersion = 1))
ggplot() +
geom_point(data = hc, aes(log(Wt), Sat), alpha = .5) +
geom_ribbon(data = curv, aes(lWt, ymin = lo, ymax = hi, fill = modelo), alpha = .25) +
geom_line(data = curv, aes(lWt, mu, colour = modelo, linetype = modelo), linewidth = 1) +
scale_colour_manual(values = c(col1, col2)) + scale_fill_manual(values = c(col1, col2)) +
labs(x = "log(peso)", y = "Número de satélites", colour = NULL, fill = NULL, linetype = NULL)
Las diferencias aparecen en los pesos grandes (sistemático) y en la forma de las bandas (aleatorio).