Modelos para proporciones: GLM binomiales

Dunn & Smyth (2018), Capítulo 9

Nelson Muriel

Universidad Iberoamericana CDMX

Plan de la sesión

Bloque	Tema	Libro	Min
0	El GLM binomial y formas de especificar la respuesta	§9.2	15
1	Funciones de enlace y distribuciones de tolerancia	§9.3–9.4	25
2	Odds, razones de momios y logit	§9.5	20
3	Dosis efectiva mediana, ED50	§9.6	15
	Descanso		10
4	Log-log complementario: ensayos de dilución	§9.7	15
5	Sobredispersión	§9.8	30
6	Wald falla (Hauck–Donner) y datos binarios	§9.9–9.10	20
7	Estudio de caso: insecticidas	§9.11	25
8	Cierre y ejercicios	§9.13	5

Datos: paquete GLMsData. El código del libro se moderniza con ggplot2, broom, car::deltaMethod() y alternativas actuales (DHARMa, glmmTMB, logistf).

Bloque 0 — El GLM binomial

Paquetería

library(GLMsData)
library(ggplot2)
library(broom)
library(car)
library(statmod)
theme_set(theme_minimal(base_size = 14))

data(turbines); data(germ); data(mammary); data(nminer); data(deposit)

La distribución binomial como EDM (§9.2)

\(y = (\text{éxitos})/m\) con \(m\) ensayos independientes, \(y = 0, 1/m, \dots, 1\):

\[ f(y;\mu,m) = \binom{m}{my}\mu^{my}(1-\mu)^{m(1-y)}, \qquad \phi = 1 \]

Pesos previos: \(w = m\)

La distribución binomial como EDM (p 2)

Desviación unitaria: \[ d(y,\mu) = 2\left\{ y\log\frac{y}{\mu} + (1-y)\log\frac{1-y}{1-\mu} \right\} \]
Desviación residual: \(D = \sum_i m_i\, d(y_i,\widehat{\mu_i}) \;\overset{\text{a}}{\sim}\; \chi^2_{n-p'}\)
La aproximación de punto de silla es adecuada si \(\min\{m_i y_i\} \ge 3\) y \(\min\{m_i(1-y_i)\} \ge 3\)

Ejemplo 9.1: turbinas con fisuras

turbines$prop <- turbines$Fissures / turbines$Turbines
turbines$caso <- seq_len(nrow(turbines))     # numeración de la Tabla 9.1 (no viene en los datos)
turbines

   Hours Turbines Fissures       prop caso
1    400       39        0 0.00000000    1
2   1000       53        4 0.07547170    2
3   1400       33        2 0.06060606    3
4   1800       73        7 0.09589041    4
5   2200       30        5 0.16666667    5
6   2600       39        9 0.23076923    6
7   3000       42        9 0.21428571    7
8   3400       13        6 0.46153846    8
9   3800       34       22 0.64705882    9
10  4200       40       21 0.52500000   10
11  4600       36       21 0.58333333   11

ggplot(turbines, aes(Hours, prop, size = Turbines)) +
  geom_point(shape = 1) +
  geom_text(aes(label = caso), nudge_y = 0.04, size = 3.5, show.legend = FALSE) +
  labs(x = "Horas de uso", y = "Proporción con fisuras", size = "Turbinas (m)")

Tres formas de especificar la respuesta

# (1) proporciones + pesos (m)
tur_m1 <- glm(Fissures / Turbines ~ Hours, family = binomial,
              weights = Turbines, data = turbines)
# (2) matriz de éxitos y fracasos   <-- la recomendada
tur_m2 <- glm(cbind(Fissures, Turbines - Fissures) ~ Hours, family = binomial,
              data = turbines)
# (3) una fila por individuo (Bernoulli): hay que expandir los datos
tur_largo <- with(turbines, data.frame(
  Hours  = rep(Hours, Turbines),
  fisura = unlist(Map(function(f, m) rep(c(1, 0), c(f, m - f)), Fissures, Turbines))
))
nrow(tur_largo)

[1] 432

tur_m3 <- glm(fisura ~ Hours, family = binomial, data = tur_largo)

rbind(m1 = coef(tur_m1), m2 = coef(tur_m2), m3 = coef(tur_m3))

   (Intercept)        Hours
m1   -3.923597 0.0009992372
m2   -3.923597 0.0009992372
m3   -3.923597 0.0009992372

Mismos coeficientes, pero la forma Bernoulli cambia la desviación y los gl (volveremos en el Bloque 6):

sapply(list(m1 = tur_m1, m2 = tur_m2, m3 = tur_m3),
       function(m) c(deviance = deviance(m), df = df.residual(m)))

               m1       m2       m3
deviance 10.33147 10.33147 379.0749
df        9.00000  9.00000 430.0000

La desviación “a mano”

unit_dev <- function(y, mu) {
  t1 <- ifelse(y == 0, 0, y * log(y / mu))
  t2 <- ifelse(y == 1, 0, (1 - y) * log((1 - y) / (1 - mu)))
  2 * (t1 + t2)
}
y <- with(turbines, Fissures / Turbines); m <- turbines$Turbines
c(manual = sum(m * unit_dev(y, fitted(tur_m2))), glm = deviance(tur_m2))

  manual      glm 
10.33147 10.33147

¿Se cumple la regla de la aproximación de punto de silla?

c(min_exitos = min(turbines$Fissures), min_fracasos = min(turbines$Turbines - turbines$Fissures))

  min_exitos min_fracasos 
           0            7

Hay conteos menores que 3: conviene cautela al interpretar \(\chi^2\).

Bloque 1 — Enlaces y tolerancia

Funciones de enlace para el GLM binomial (§9.3)

Se necesita \(0 < \mu < 1\):

Enlace	\(\eta = g(\mu)\)	En R
Logit (canónico, por defecto)	\(\log\{\mu/(1-\mu)\}\)	`"logit"`
Probit	\(\Phi^{-1}(\mu)\)	`"probit"`
Log-log complementario	\(\log\{-\log(1-\mu)\}\)	`"cloglog"`
Cauchit, log		`"cauchit"`, `"log"`

Logit y probit son simétricos en torno a \(\mu = 0.5\); cloglog no.

Las curvas inversas (Fig. 9.2)

eta <- seq(-5, 5, length.out = 200)
enlaces <- do.call(rbind, lapply(c("logit", "probit", "cloglog", "cauchit"), 
  function(l)
  data.frame(eta, mu = binomial(link = l)$linkinv(eta), enlace = l)))

ggplot(enlaces, aes(eta, mu, color = enlace, linetype = enlace)) +
  geom_hline(yintercept = 0.5, color = "grey75") +
  geom_vline(xintercept = 0, color = "grey75") +
  geom_line(linewidth = 1) +
  labs(x = expression("Predictor lineal " * eta), y = expression("Proporción esperada " * mu))

Ejemplo 9.2: tres enlaces sobre las turbinas

tur_logit  <- tur_m2
tur_probit <- update(tur_logit, family = binomial(link = "probit"))
tur_cll    <- update(tur_logit, family = binomial(link = "cloglog"))

mods <- list(Logit = tur_logit, Probit = tur_probit, "Comp. log-log" = tur_cll)
cbind(as.data.frame(t(sapply(mods, coef))),
      Deviance = sapply(mods, deviance),
      AIC      = sapply(mods, AIC))

              (Intercept)        Hours  Deviance      AIC
Logit           -3.923597 0.0009992372 10.331466 49.80819
Probit          -2.275807 0.0005783211  9.814837 49.29156
Comp. log-log   -3.603280 0.0008104936 12.227914 51.70463

Los coeficientes son distintos (escalas distintas), pero el ajuste es casi el mismo:

nuevas <- data.frame(Hours = seq(0, 5000, length.out = 100))
curvas <- do.call(rbind, lapply(names(mods), function(nm)
  data.frame(enlace = nm, nuevas, mu = predict(mods[[nm]], nuevas, type = "response"))))

ggplot() +
  geom_point(data = turbines, aes(Hours, prop), size = 2) +
  geom_line(data = curvas, aes(Hours, mu, color = enlace, linetype = enlace), linewidth = 1) +
  labs(x = "Horas de uso", y = "Proporción con fisuras", color = NULL, linetype = NULL)

Interpretación por umbral (§9.4)

Idea: Cada turbina tiene una tolerancia \(t_i \sim N(\tau_i, \sigma^2)\) con \(\tau_i = \beta_0' + \beta_1' x_i\), y desarrolla fisuras si \(t_i \le T\):

\[ \begin{split} \mu_i =& \Pr(t_i \le T) = \Phi\!\left(\frac{T - \tau_i}{\sigma}\right) = \Phi(\beta_0 + \beta_1 x_i),\\ \beta_0 =& \frac{T-\beta_0'}{\sigma},\qquad \beta_1 = -\frac{\beta_1'}{\sigma} \end{split} \]

Intepretación por umbral, enlaces

Enlace	Distribución de tolerancia
Logit	Logística
Probit	Normal
Cloglog	Valor extremo (Gumbel)
Cauchit	Cauchy

Verificación por simulación

set.seed(1)
n <- 20000; T_umbral <- 2000; b0p <- 5000; b1p <- -0.8; sigma <- 1500
x   <- runif(n, 0, 5000)
tol <- rnorm(n, mean = b0p + b1p * x, sd = sigma)
y_s <- as.numeric(tol <= T_umbral)

fit_tol <- glm(y_s ~ x, family = binomial(link = "probit"))
rbind(estimado = coef(fit_tol),
      teorico  = c((T_umbral - b0p) / sigma, -b1p / sigma))

         (Intercept)            x
estimado   -1.989722 0.0005306401
teorico    -2.000000 0.0005333333

Con tolerancia logística se recupera el logit:

s <- 800
tol_l <- rlogis(n, location = b0p + b1p * x, scale = s)
fit_tol_l <- glm(as.numeric(tol_l <= T_umbral) ~ x, family = binomial(link = "logit"))
rbind(estimado = coef(fit_tol_l),
      teorico  = c((T_umbral - b0p) / s, -b1p / s))

         (Intercept)            x
estimado   -3.661032 0.0009751574
teorico    -3.750000 0.0010000000

Bloque 2 — Odds y razones de momios

Odds y el enlace logit (§9.5)

Si \(A\) ocurre con probabilidad \(\mu\): \(\text{odds}(A) = \mu/(1-\mu)\). Por ejemplo, \(\mu = 0.6 \Rightarrow\) odds \(= 1.5\).

\[ \log(\text{odds}) = \beta_0 + \beta_1 x \quad\Longleftrightarrow\quad \text{odds} = \exp(\beta_0)\,\{\exp(\beta_1)\}^x \]

Si \(x\) sube una unidad, los log-odds suben \(\beta_1\)
Equivalente: los odds se multiplican por \(\exp(\beta_1)\)
Con \(x \in \{0, 1\}\): \(\exp(\beta_1)\) es la razón de momios (odds ratio, OR)

Ejemplo 9.3: turbinas

coef(tur_logit)

  (Intercept)         Hours 
-3.9235965551  0.0009992372

exp(coef(tur_logit)["Hours"])            # factor por cada hora

Hours 
1.001

exp(1000 * coef(tur_logit)["Hours"])     # factor por cada 1000 horas

   Hours 
2.716209

El modelo supone log-odds lineales en Hours. Lo comprobamos con el logit empírico (sumando 0.5 para poder calcularlo cuando \(y = 0\)):

turbines$elogit <- with(turbines, log((Fissures + 0.5) / (Turbines - Fissures + 0.5)))
turbines$eta    <- predict(tur_logit)

ggplot(turbines, aes(Hours)) +
  geom_point(aes(y = elogit)) +
  geom_line(aes(y = eta), color = "darkred") +
  labs(x = "Horas de uso", y = "Log-odds")
ggplot(turbines, aes(Hours)) +
  geom_point(aes(y = exp(elogit))) +
  geom_line(aes(y = exp(eta)), color = "darkred") +
  labs(x = "Horas de uso", y = "Odds")

Ejemplo 9.4: germinación de semillas

Dos tipos de semilla (O. aegyptiaca 73 y 75) y dos extractos de raíz (frijol y pepino).

str(germ)

'data.frame':   21 obs. of  4 variables:
 $ Germ   : int  10 23 23 26 17 5 53 55 32 46 ...
 $ Total  : int  39 62 81 51 39 6 74 72 51 79 ...
 $ Extract: Factor w/ 2 levels "Bean","Cucumber": 1 1 1 1 1 2 2 2 2 2 ...
 $ Seeds  : Factor w/ 2 levels "OA73","OA75": 2 2 2 2 2 2 2 2 2 2 ...

germ$Prop <- germ$Germ / germ$Total

rbind(data.frame(factor = "Extracto", nivel = germ$Extract, Prop = germ$Prop),
      data.frame(factor = "Semillas", nivel = germ$Seeds,   Prop = germ$Prop)) |>
  ggplot(aes(nivel, Prop)) +
  geom_boxplot() +
  facet_wrap(~ factor, scales = "free_x") +
  labs(x = NULL, y = "Germ/Total")

Razones de momios con intervalos

gm_m1 <- glm(cbind(Germ, Total - Germ) ~ Seeds + Extract, family = binomial, data = germ)
tidy(gm_m1)

# A tibble: 3 × 5
  term            estimate std.error statistic  p.value
  <chr>              <dbl>     <dbl>     <dbl>    <dbl>
1 (Intercept)       -0.700     0.151     -4.65 3.36e- 6
2 SeedsOA75          0.270     0.155      1.75 8.04e- 2
3 ExtractCucumber    1.06      0.144      7.38 1.55e-13

tidy(gm_m1, exponentiate = TRUE, conf.int = TRUE)

# A tibble: 3 × 7
  term            estimate std.error statistic  p.value conf.low conf.high
  <chr>              <dbl>     <dbl>     <dbl>    <dbl>    <dbl>     <dbl>
1 (Intercept)        0.496     0.151     -4.65 3.36e- 6    0.368     0.665
2 SeedsOA75          1.31      0.155      1.75 8.04e- 2    0.968     1.78 
3 ExtractCucumber    2.90      0.144      7.38 1.55e-13    2.19      3.85

Con extracto de pepino, los odds de germinar son \(\approx 2.9\) veces los del extracto de frijol (con el mismo tipo de semilla).

Wald vs. perfil de verosimilitud:

rbind(Wald = exp(confint.default(gm_m1))["ExtractCucumber", ],
      Perfil = exp(confint(gm_m1))["ExtractCucumber", ])

          2.5 %   97.5 %
Wald   2.186048 3.847426
Perfil 2.189641 3.854706

Bloque 3 — ED50

Dosis efectiva mediana (§9.6)

En experimentos dosis–respuesta, ED50 es la dosis con \(\mu = 0.5\). Con \(g(\mu) = \beta_0 + \beta_1 d\):

\[ \mathrm{ED}(\rho) = \frac{g(\rho) - \beta_0}{\beta_1}, \qquad \mathrm{ED50}_{\text{logit}} = -\frac{\beta_0}{\beta_1} \]

Es no lineal en \(\beta\): el error estándar se obtiene con el método delta, con gradiente \(\left(-1/\beta_1,\ \beta_0/\beta_1^2\right)\).

ED50 de las turbinas

A mano con el método delta:

b <- coef(tur_logit); V <- vcov(tur_logit)
ed50_hat <- unname(-b[1] / b[2])
grad     <- c(-1 / b[2], b[1] / b[2]^2)
se_ed50  <- sqrt(drop(t(grad) %*% V %*% grad))
c(ED50 = ed50_hat, SE = se_ed50)

     ED50        SE 
3926.5916  158.0381

Con MASS::dose.p() y con car::deltaMethod() (la misma herramienta de la sesión sobre pruebas de hipótesis):

MASS::dose.p(tur_logit)

             Dose       SE
p = 0.5: 3926.592 158.0381

deltaMethod(tur_logit, "-`(Intercept)` / Hours")

                 Estimate      SE   2.5 % 97.5 %
-Intercept/Hours  3926.59  158.04 3616.84 4236.3

Otros niveles y otros enlaces

MASS::dose.p(tur_logit, p = c(0.1, 0.5, 0.9))

             Dose       SE
p = 0.1: 1727.690 203.7928
p = 0.5: 3926.592 158.0381
p = 0.9: 6125.493 366.9071

sapply(mods, function(m) {
  d <- MASS::dose.p(m)
  c(ED50 = as.numeric(d), SE = as.numeric(attr(d, "SE")))
})

         Logit    Probit Comp. log-log
ED50 3926.5916 3935.1971     3993.5748
SE    158.0381  167.6132      137.8352

IC de Wald al 95% (logit):

ed50_hat + c(-1, 1) * qnorm(0.975) * se_ed50

[1] 3616.843 4236.341

Correr las turbinas unas 3900 horas produce fisuras en alrededor de la mitad de ellas, con cualquiera de los tres enlaces.

Descanso

10 minutos

Bloque 4 — Log-log complementario

Ensayos de dilución (§9.7)

Un cultivo con “dosis” \(d_i\) da un resultado positivo si contiene al menos una célula activa. Si el número de células sigue una Poisson con media \(\lambda d_i\):

\[ 1 - \mu_i = e^{-\lambda d_i} \;\Longrightarrow\; \log\{-\log(1-\mu_i)\} = \log\lambda + \log d_i \]

GLM binomial con enlace cloglog y offset \(\log d_i\):

\[ g(\mu_i) = \beta_0 + \log d_i, \qquad \hat\lambda = \exp(\hat\beta_0), \qquad 1/\hat\lambda = \text{células por cultivo positivo} \]

Prueba del modelo de un impacto (single-hit): ajustar \(g(\mu_i) = \beta_0 + \beta_1\log d_i\) y contrastar \(H_0:\beta_1 = 1\).

Ejemplo 9.6: células madre mamarias

mammary

  N.Cells N.Assays N.Outgrowths
1      15       38            3
2      40        6            6
3      60       17           13
4      90        8            6
5     125       12            9

fit_mm <- glm(cbind(N.Outgrowths, N.Assays - N.Outgrowths) ~ offset(log(N.Cells)),
              family = binomial(link = "cloglog"), data = mammary)
summary(fit_mm)$coefficients

             Estimate Std. Error   z value      Pr(>|z|)
(Intercept) -4.163625  0.1744346 -23.86925 6.391454e-126

1 / exp(coef(fit_mm)[1])

(Intercept) 
   64.30418

Una célula madre por cada \(\approx 64\) células. Intervalos de confianza (se invierten, porque \(1/e^{x}\) es decreciente):

est <- summary(fit_mm)$coefficients
ic_wald   <- 1 / exp(est[1, "Estimate"] + c(1, -1) * qnorm(0.975) * est[1, "Std. Error"])
ic_perfil <- rev(1 / exp(confint(fit_mm)))
rbind(Wald = ic_wald, Perfil = ic_perfil)

         97.5 %    2.5 %
Wald   45.68359 90.51452
Perfil 46.27686 92.08723

¿Se cumple el modelo de un impacto?

fit_mm1 <- glm(cbind(N.Outgrowths, N.Assays - N.Outgrowths) ~ log(N.Cells),
               family = binomial(link = "cloglog"), data = mammary)
anova(fit_mm, fit_mm1, test = "Chisq")

Analysis of Deviance Table

Model 1: cbind(N.Outgrowths, N.Assays - N.Outgrowths) ~ offset(log(N.Cells))
Model 2: cbind(N.Outgrowths, N.Assays - N.Outgrowths) ~ log(N.Cells)
  Resid. Df Resid. Dev Df Deviance Pr(>Chi)
1         4     16.852                     
2         3     16.205  1   0.6468   0.4213

No hay evidencia contra \(H_0: \beta_1 = 1\): se mantiene el modelo de un impacto.

Bloque 5 — Sobredispersión

Qué es y por qué importa (§9.8)

Binomial: \(\mathrm{var}[y] = \mu(1-\mu)/m\). En la práctica puede haber más variación.

Las probabilidades \(p_i\) varían aun con las mismas covariables (modelo jerárquico, beta-binomial): \[ \mathrm{var}[y_i] = \phi_i\,\mu_i(1-\mu_i)/m_i, \qquad \phi_i = 1 + (m_i - 1)\rho \]
Los \(m_i\) casos están correlacionados (conglomerados), con correlación \(\rho\)

Consecuencias: errores estándar subestimados, pruebas “demasiado significativas”, modelos demasiado complejos.

No puede haber sobredispersión con datos binarios (\(m = 1\)).

Detección: desviación y Pearson vs. gl

Primero, las turbinas (Ejemplo 9.7):

gof <- function(m) {
  c(gl = df.residual(m),
    desviacion = deviance(m),
    pearson = sum(resid(m, type = "pearson")^2))
}
gof(tur_logit)

        gl desviacion    pearson 
  9.000000  10.331466   9.250839

pchisq(gof(tur_logit)["desviacion"], gof(tur_logit)["gl"], lower.tail = FALSE)

desviacion 
 0.3243236

Sin evidencia de sobredispersión: las turbinas operan más o menos de manera independiente.

Ejemplo 9.8: germinación con interacción

Las semillas se plantan juntas en parcelas: es plausible que estén correlacionadas.

gm_int <- glm(cbind(Germ, Total - Germ) ~ Extract * Seeds, family = binomial, data = germ)
anova(gm_int, test = "Chisq")

Analysis of Deviance Table

Model: binomial, link: logit

Response: cbind(Germ, Total - Germ)

Terms added sequentially (first to last)

              Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
NULL                             20     98.719              
Extract        1   55.969        19     42.751 7.364e-14 ***
Seeds          1    3.065        18     39.686   0.08000 .  
Extract:Seeds  1    6.408        17     33.278   0.01136 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

gof(gm_int)

        gl desviacion    pearson 
  17.00000   33.27779   31.65114

pchisq(deviance(gm_int), df.residual(gm_int), lower.tail = FALSE)

[1] 0.01039184

Desviación \(33.3\) y Pearson \(31.7\) con 17 gl: evidencia de sobredispersión aun con el modelo máximo.

Descartar otras causas: residuales cuantílicos

Para GLM binomiales se recomiendan los residuales cuantílicos (statmod::qresid):

set.seed(1)
qres <- qresid(gm_int)
ggplot(data.frame(ajustado = fitted(gm_int), qres), aes(ajustado, qres)) +
  geom_hline(yintercept = 0, color = "grey60") +
  geom_point() + geom_smooth(se = FALSE, color = "darkred") +
  labs(x = "Valor ajustado", y = "Residual cuantílico")

Sin residuales grandes, sin observaciones influyentes (diseño balanceado): la causa es sobredispersión.

Alternativa moderna: residuales simulados con DHARMa

sim <- DHARMa::simulateResiduals(gm_int, plot = FALSE)
DHARMa::testDispersion(sim)


    DHARMa nonparametric dispersion test via sd of residuals fitted vs.
    simulated

data:  simulationOutput
dispersion = 1.0185, p-value = 0.816
alternative hypothesis: two.sided

plot(sim)

(Se evalúa sólo si DHARMa está instalado.)

Modelo cuasi-binomial

\(V(\mu) = \mu(1-\mu)\) pero con dispersión libre \(\phi\) (estimada con Pearson):

gm_od <- update(gm_int, family = quasibinomial)
anova(gm_od, test = "F")                 # F porque phi se estima

Analysis of Deviance Table

Model: quasibinomial, link: logit

Response: cbind(Germ, Total - Germ)

Terms added sequentially (first to last)

              Df Deviance Resid. Df Resid. Dev       F    Pr(>F)    
NULL                             20     98.719                      
Extract        1   55.969        19     42.751 30.0610 4.043e-05 ***
Seeds          1    3.065        18     39.686  1.6462   0.21669    
Extract:Seeds  1    6.408        17     33.278  3.4418   0.08099 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Con binomial había interacción significativa; con cuasi-binomial sólo Extract lo es (la interacción queda marginal).

Qué cambia: los errores estándar

phi <- summary(gm_od)$dispersion
c(phi_Pearson = phi,
  chequeo = sum(resid(gm_int, type = "pearson")^2) / df.residual(gm_int))

phi_Pearson     chequeo 
   1.861832    1.861832

data.frame(Estimate = coef(summary(gm_int))[, "Estimate"],
           Binom.SE = coef(summary(gm_int))[, "Std. Error"],
           Quasi.SE = coef(summary(gm_od))[, "Std. Error"],
           Razon    = coef(summary(gm_od))[, "Std. Error"] / coef(summary(gm_int))[, "Std. Error"])

                            Estimate  Binom.SE  Quasi.SE   Razon
(Intercept)               -0.4122448 0.1841784 0.2513095 1.36449
ExtractCucumber            0.5400782 0.2498130 0.3408672 1.36449
SeedsOA75                 -0.1459269 0.2231659 0.3045076 1.36449
ExtractCucumber:SeedsOA75  0.7781037 0.3064332 0.4181249 1.36449

sqrt(phi)

[1] 1.36449

Mismas estimaciones; los EE se multiplican por \(\sqrt{\hat\phi}\). El cuasi-binomial no es un modelo de probabilidad: no hay AIC.

Alternativa basada en verosimilitud: beta-binomial

gm_bb <- glmmTMB::glmmTMB(cbind(Germ, Total - Germ) ~ Extract * Seeds,
                          family = glmmTMB::betabinomial, data = germ)
summary(gm_bb)$coefficients$cond

                            Estimate Std. Error    z value   Pr(>|z|)
(Intercept)               -0.4445600  0.2182557 -2.0368773 0.04166234
ExtractCucumber            0.5221435  0.2968227  1.7591089 0.07855901
SeedsOA75                 -0.0973906  0.2736686 -0.3558706 0.72193751
ExtractCucumber:SeedsOA75  0.7979248  0.3779548  2.1111644 0.03475819
attr(,"ddf")
[1] "asymptotic"

c(AIC_binomial = AIC(gm_int), AIC_betabinomial = AIC(gm_bb))

    AIC_binomial AIC_betabinomial 
        117.8740         117.5335

La beta-binomial sí es un modelo de probabilidad completo (permite AIC, simulación, predicción de la distribución).

¿Cuánto importa ignorarlo? Simulación

Sin efecto verdadero de \(x\), ¿con qué frecuencia rechazamos \(H_0\) al 5%? Datos beta-binomiales con \(a = b = 5\) (\(\rho \approx 0.09\)):

sim_una <- function(n = 20, m = 30, a = 5, b = 5) {
  p <- rbeta(n, a, b)
  x <- rep(0:1, each = n / 2)
  y <- rbinom(n, m, p)
  d <- data.frame(y, m, x)
  f0 <- glm(cbind(y, m - y) ~ 1, family = binomial, data = d)
  f1 <- glm(cbind(y, m - y) ~ x, family = binomial, data = d)
  q0 <- glm(cbind(y, m - y) ~ 1, family = quasibinomial, data = d)
  q1 <- glm(cbind(y, m - y) ~ x, family = quasibinomial, data = d)
  c(binomial = anova(f0, f1, test = "Chisq")[2, "Pr(>Chi)"],
    cuasi    = anova(q0, q1, test = "F")[2, "Pr(>F)"])
}
set.seed(2026)
pvals <- replicate(500, sim_una())
rowMeans(pvals < 0.05)

binomial    cuasi 
   0.268    0.032

La binomial rechaza muy por encima del 5% nominal; la cuasi-binomial se acerca al 5%.

Bloque 6 — Wald falla y datos binarios

Casi-separación y Hauck–Donner (§9.9)

Si en algún nivel de un factor todos los \(y_i\) son 0 (o todos 1), \(\hat\mu \to 0\) ó \(1\), así que \(\hat\eta \to \pm\infty\) y algún coeficiente tiende a \(\pm\infty\).

Los errores estándar crecen más rápido que los coeficientes
El estadístico de Wald \(\hat\beta/\mathrm{se}(\hat\beta) \to 0\), sea cual sea la significancia real (efecto Hauck–Donner)
El cociente de verosimilitud y la prueba score siguen siendo confiables (invariantes a reparametrizaciones)

Ejemplo 9.9: mineros ruidosos y eucaliptos

nminer$Eucs15 <- nminer$Eucs > 15
table(Miners = nminer$Miners, Eucs15 = nminer$Eucs15)

      Eucs15
Miners FALSE TRUE
     0    14    0
     1     6   11

ggplot(nminer, aes(Eucs, Miners)) +
  geom_point(size = 2) +
  geom_vline(xintercept = 15.5, color = "grey50") +
  labs(x = "Eucaliptos", y = "Presencia de mineros ruidosos")

Con más de 15 eucaliptos, siempre hay mineros ruidosos.

Wald vs. LR vs. score

nm1 <- glm(Miners ~ Eucs15, family = binomial, data = nminer)
nm0 <- glm(Miners ~ 1, family = binomial, data = nminer)
tidy(nm1)                                  # Wald: EE enorme, p ~ 0.995

# A tibble: 2 × 5
  term        estimate std.error statistic p.value
  <chr>          <dbl>     <dbl>     <dbl>   <dbl>
1 (Intercept)   -0.847     0.488  -1.74     0.0825
2 Eucs15TRUE    20.4    3242.      0.00630  0.995

anova(nm1, test = "Chisq")                 # razón de verosimilitud

Analysis of Deviance Table

Model: binomial, link: logit

Response: Miners

Terms added sequentially (first to last)

       Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
NULL                      30     42.684              
Eucs15  1    18.25        29     24.435 1.937e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

anova(nm0, nm1, test = "Rao")              # prueba score

Analysis of Deviance Table

Model 1: Miners ~ 1
Model 2: Miners ~ Eucs15
  Resid. Df Resid. Dev Df Deviance    Rao  Pr(>Chi)    
1        30     42.684                                 
2        29     24.435  1    18.25 14.041 0.0001788 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

glm.scoretest(nm0, nminer$Eucs15)          # versión de statmod (estadístico z)

[1] 3.747173

tapply(fitted(nm1), nminer$Eucs15, mean)   # mu_hat = 1 cuando Eucs15 = TRUE

FALSE  TRUE 
  0.3   1.0

El Wald dice “no significativo” (\(P \approx 0.995\)); LR y score dicen “altamente significativo”. Confiar en LR/score.

Una solución moderna: corrección de Firth

La regresión logística con sesgo reducido (Firth) produce estimadores finitos aun con separación:

lf <- logistf::logistf(Miners ~ Eucs15, data = nminer)
summary(lf)

logistf::logistf(formula = Miners ~ Eucs15, data = nminer)

Model fitted by Penalized ML
Coefficients:
                  coef  se(coef) lower 0.95 upper 0.95     Chisq            p
(Intercept) -0.8023465 0.4720293  -1.800369 0.08472562  3.125978 7.705361e-02
Eucs15TRUE   3.9378407 1.5197922   1.661929 8.84865570 15.331482 9.020075e-05
            method
(Intercept)      2
Eucs15TRUE       2

Method: 1-Wald, 2-Profile penalized log-likelihood, 3-None

Likelihood ratio test=15.33148 on 1 df, p=9.020075e-05, n=31
Wald test = 7.600098 on 1 df, p = 0.005836512

(Se evalúa sólo si logistf está instalado.)

Datos binarios: no hay bondad de ajuste (§9.10)

Con \(m_i = 1\) para todo \(i\), la desviación y Pearson dependen sólo de los valores ajustados:

mu_hat <- fitted(nm1); y_b <- nminer$Miners
c(desviacion = deviance(nm1),
  a_mano = -2 * sum(ifelse(y_b == 1, log(mu_hat), log(1 - mu_hat))))

desviacion     a_mano 
  24.43457   24.43457

c(gl = df.residual(nm1), desviacion = deviance(nm1))

        gl desviacion 
  29.00000   24.43457

Desviación \(<\) gl no indica subdispersión: no tiene significado.

Para binarios: pruebas LR / score (con \(p'\) mucho menor que \(n\)), y evaluar calibración y discriminación (ROC/AUC).

Bloque 7 — Estudio de caso: insecticidas

Los datos (§9.11)

Lotes de insectos expuestos a depósitos (mg) de tres insecticidas; proporción muerta tras la exposición.

str(deposit)

'data.frame':   18 obs. of  4 variables:
 $ Killed     : int  3 5 19 19 24 35 2 14 20 27 ...
 $ Number     : int  50 49 47 38 29 50 50 49 50 50 ...
 $ Insecticide: Factor w/ 3 levels "A","B","C": 1 1 1 1 1 1 2 2 2 2 ...
 $ Deposit    : num  2 2.64 3.48 4.59 6.06 8 2 2.64 3.48 4.59 ...

deposit$Prop   <- deposit$Killed / deposit$Number
deposit$logDep <- log(deposit$Deposit)
deposit$elogit <- with(deposit, log((Killed + 0.5) / (Number - Killed + 0.5)))

ggplot(deposit, aes(Deposit, Prop, color = Insecticide, shape = Insecticide)) +
  geom_point(size = 3) + geom_line() +
  labs(x = "Depósito (mg)", y = "Proporción de insectos muertos")

A y B parecen similares; C parece distinto.

Modelo 1: `Deposit + Insecticide`

ins_m1 <- glm(cbind(Killed, Number - Killed) ~ Deposit + Insecticide,
              family = binomial, data = deposit)
coef(ins_m1)

 (Intercept)      Deposit InsecticideB InsecticideC 
  -3.2213638    0.6316762    0.3695267    2.6880162

MASS::dose.p(ins_m1, c(1, 2))        # ED50 del insecticida A

             Dose        SE
p = 0.5: 5.099708 0.2468085

Para B y C, el intercepto no está en un solo parámetro: reajustamos sin constante (un intercepto por insecticida):

ins_m1A <- update(ins_m1, . ~ . - 1)
coef(ins_m1A)

     Deposit InsecticideA InsecticideB InsecticideC 
   0.6316762   -3.2213638   -2.8518371   -0.5333477

cbind(A = MASS::dose.p(ins_m1A, c(2, 1)),
      B = MASS::dose.p(ins_m1A, c(3, 1)),
      C = MASS::dose.p(ins_m1A, c(4, 1)))

                A        B         C
p = 0.5: 5.099708 4.514714 0.8443372

Diagnóstico del modelo 1

set.seed(1)
diag1 <- data.frame(ajustado = fitted(ins_m1), qres = qresid(ins_m1),
                    Insecticide = deposit$Insecticide)
ggplot(diag1, aes(ajustado, qres, color = Insecticide)) +
  geom_hline(yintercept = 0, color = "grey60") +
  geom_point() + geom_line() +
  labs(x = "Valor ajustado", y = "Residual cuantílico", title = "Patrón sistemático")
ggplot(deposit, aes(Deposit, elogit, color = Insecticide)) +
  geom_point() + geom_line() +
  labs(x = "Depósito (mg)", y = "Logit empírico", title = "Log-odds no lineales")

El modelo 1 es inadecuado: sobrestima/subestima sistemáticamente en dosis bajas y altas, y los log-odds no son lineales en Deposit.

Modelo 2: `log(Deposit)`

En dosis–respuesta es habitual usar el logaritmo de la dosis:

ins_m2 <- glm(cbind(Killed, Number - Killed) ~ logDep + Insecticide - 1,
              family = binomial, data = deposit)
ED50_m2_log <- cbind(A = MASS::dose.p(ins_m2, c(2, 1)),
                     B = MASS::dose.p(ins_m2, c(3, 1)),
                     C = MASS::dose.p(ins_m2, c(4, 1)))
exp(ED50_m2_log)                     # ED50 en mg

                A        B        C
p = 0.5: 4.688232 4.154625 1.753202

Mejor, pero aún subestima las proporciones en dosis bajas y altas: los log-odds no son constantes en \(\log d\), sino quizá cuadráticos.

Modelo 3: cuadrático en \(\log d\)

ins_m3 <- glm(cbind(Killed, Number - Killed) ~ poly(logDep, 2) + Insecticide,
              family = binomial, data = deposit)
anova(ins_m2, ins_m3, test = "Chisq")

Analysis of Deviance Table

Model 1: cbind(Killed, Number - Killed) ~ logDep + Insecticide - 1
Model 2: cbind(Killed, Number - Killed) ~ poly(logDep, 2) + Insecticide
  Resid. Df Resid. Dev Df Deviance Pr(>Chi)   
1        14     23.385                        
2        13     15.090  1   8.2949 0.003976 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

c(deviance = deviance(ins_m3), gl = df.residual(ins_m3))

deviance       gl 
15.09036 13.00000

AIC(ins_m1, ins_m2, ins_m3)

       df       AIC
ins_m1  4 118.22335
ins_m2  4  93.58276
ins_m3  5  87.28790

Mejora significativa, y la desviación (\(\approx 15\) con 13 gl) no muestra sobredispersión. Pero:

c(min(deposit$Killed), min(deposit$Number - deposit$Killed))

[1] 2 0

Hay conteos menores que 3: la aproximación de punto de silla es dudosa y la conclusión no es del todo confiable.

Curvas ajustadas de los tres modelos

malla <- expand.grid(Deposit = seq(min(deposit$Deposit), max(deposit$Deposit), length.out = 100),
                     Insecticide = levels(deposit$Insecticide))
malla$logDep <- log(malla$Deposit)
curvas_ins <- rbind(
  data.frame(modelo = "m1: Deposit",         malla, mu = predict(ins_m1, malla, type = "response")),
  data.frame(modelo = "m2: log(Deposit)",    malla, mu = predict(ins_m2, malla, type = "response")),
  data.frame(modelo = "m3: poly(logDep, 2)", malla, mu = predict(ins_m3, malla, type = "response")))

ggplot() +
  geom_point(data = deposit, aes(Deposit, Prop, color = Insecticide), size = 2) +
  geom_line(data = curvas_ins, aes(Deposit, mu, color = Insecticide), linewidth = 0.9) +
  facet_wrap(~ modelo) +
  labs(x = "Depósito (mg)", y = "Proporción de insectos muertos")

ED50 con el modelo cuadrático (Problema 9.3)

dose.p() no sirve con el término cuadrático, pero se despeja a mano. Para el insecticida \(j\) con intercepto \(a_j\):

\[ a_j + \beta_1 u + \beta_2 u^2 = 0, \quad u = \log(\text{dosis}) \;\Rightarrow\; u^* = \frac{-\beta_1 + \sqrt{\beta_1^2 - 4\beta_2 a_j}}{2\beta_2} \]

(la raíz con signo \(+\) es la que tiene logit creciente). Hay que reajustar con potencias crudas, porque con poly() los coeficientes son ortogonales y no tienen esta forma:

ins_m3r <- glm(cbind(Killed, Number - Killed) ~ Insecticide + logDep + I(logDep^2),
               family = binomial, data = deposit)
all.equal(fitted(ins_m3), fitted(ins_m3r))     # mismo modelo, otra parametrización

[1] TRUE

EE por método delta

ed50_cuad <- function(a_expr) {
  g <- sprintf("exp((-logDep + sqrt(logDep^2 - 4*`I(logDep^2)`*(%s))) / (2*`I(logDep^2)`))", a_expr)
  deltaMethod(ins_m3r, g)
}
rbind(A = ed50_cuad("`(Intercept)`"),
      B = ed50_cuad("`(Intercept)` + InsecticideB"),
      C = ed50_cuad("`(Intercept)` + InsecticideC"))

  Estimate       SE    2.5 % 97.5 %
A 4.300614 0.275063 3.761500 4.8397
B 3.839975 0.203386 3.441345 4.2386
C 1.968594 0.092697 1.786912 2.1503

ED50 en mg, con errores estándar delta. Si el estimado cae por debajo de 2 mg (la dosis mínima observada), es una extrapolación.

Cierre

Ejercicios (Problemas 9.2 y 9.5)

9.2. Demuestre que \(\mathrm{ED50} = -\beta_0/\beta_1\) (probit), \(\{\log(\log 2) - \beta_0\}/\beta_1\) (cloglog) y \((\log 0.5 - \beta_0)/\beta_1\) (enlace log). Verifíquelo con dose.p() en las turbinas.

9.5. Con la Tabla 9.6 del libro (hipertensión y apnea del sueño): prueba de Wald, IC al 95% y razones de momios para cada coeficiente:

tabla96 <- data.frame(variable = c("Intercepto", "Edad", "Sexo", "IMC", "IAH"),
                      est = c(-6.949, 0.805, 0.161, 0.332, 0.116),
                      ee  = c(0.377, 0.0444, 0.113, 0.0393, 0.0204))
tabla96$z     <- tabla96$est / tabla96$ee
tabla96$p     <- 2 * pnorm(-abs(tabla96$z))
tabla96$OR    <- exp(tabla96$est)
tabla96$OR_li <- exp(tabla96$est - qnorm(0.975) * tabla96$ee)
tabla96$OR_ls <- exp(tabla96$est + qnorm(0.975) * tabla96$ee)
round(tabla96[, -1], 3)

     est    ee       z     p    OR OR_li OR_ls
1 -6.949 0.377 -18.432 0.000 0.001 0.000 0.002
2  0.805 0.044  18.131 0.000 2.237 2.050 2.440
3  0.161 0.113   1.425 0.154 1.175 0.941 1.466
4  0.332 0.039   8.448 0.000 1.394 1.290 1.505
5  0.116 0.020   5.686 0.000 1.123 1.079 1.169

¿Qué variables parecen importantes? Interprete cada razón de momios. Para el inciso de predicción, codifique las covariables dicotomizadas según el enunciado del libro y use plogis(eta).

Mensajes clave

GLM binomial: respuesta = proporción con pesos \(m\); usar cbind(éxitos, fracasos)
Logit (odds), probit (tolerancia normal), cloglog (asimétrico, diluciones): ajustes parecidos, coeficientes en escalas distintas
\(\exp(\beta)\) en el logit es una razón de momios
\(\mathrm{ED}(\rho) = (g(\rho) - \beta_0)/\beta_1\) con error estándar por método delta
Sobredispersión: comparar desviación y Pearson con los gl; cuasi-binomial (EE \(\times\sqrt\phi\), pruebas F) o beta-binomial
Con casi-separación Wald falla: usar LR o score (o Firth)
Con datos binarios no hay bondad de ajuste basada en la desviación

Referencias

Dunn, P. K., Smyth, G. K. (2018). Generalized Linear Models With Examples in R. Springer. Capítulo 9.
Dunn, P. K., Smyth, G. K. (2017). GLMsData: Generalized Linear Model Data Sets. R package.
Hauck, W. W., Donner, A. (1977). Wald’s test as applied to hypotheses in logit analysis. JASA, 72, 851–853.
Firth, D. (1993). Bias reduction of maximum likelihood estimates. Biometrika, 80, 27–38.
Hartig, F. (2022). DHARMa: Residual Diagnostics for Hierarchical (Multi-Level / Mixed) Regression Models. R package.
Brooks, M. E. et al. (2017). glmmTMB balances speed and flexibility among packages for zero-inflated generalized linear mixed modeling. The R Journal, 9(2), 378–400.