14. Regresión de Poisson y la idea del MLG (en espanol) - Regression Modeling and Analysis (MATH 4210)

14. Regresión de Poisson y la idea del MLG¶

En 1973 dos biólogos, Michael Johnson y Peter Raven, publicaron un conteo. Para cada una de 30 islas del archipiélago de las Galápagos registraron cuántas especies de plantas crecían allí. Junto a cada conteo anotaron un puñado de hechos geográficos: el área de la isla, su elevación máxima, su distancia a la isla vecina más cercana, su distancia a la gran isla de Santa Cruz, y el área de la isla adyacente más próxima. La pregunta era vieja y sencilla. ¿Por qué algunas islas albergan cientos de especies y otras solo una docena?

Ya te has topado con este conjunto de datos dos veces. En el Capítulo 9 fue una advertencia sobre diagnósticos. Ajusta una recta por mínimos cuadrados ordinarios, y una isla, Isabela, la más grande de la cadena, queda tan lejos en el espacio de predictores que dobla el ajuste entero hacia ella misma. En el Capítulo 10 probaste el remedio estándar, transformar la respuesta, y ayudó con las gráficas de residuos pero nunca terminó de sentirse bien. Algo en los datos se resistía al modelo lineal, y te prometimos una respuesta honesta más adelante. Más adelante es ahora.

El problema se ve en una sola imagen. Figure 1 grafica las especies contra el área de la isla con la recta de mínimos cuadrados trazada por encima. El ajuste hace algo imposible. Su banda de predicción tiene un único ancho fijo, así que empujada hacia las islas pequeñas la banda baja por debajo de cero y promete conteos negativos de especies. Y un ancho fijo no puede ajustar ambos extremos: las islas pequeñas se agrupan apretadas cerca de cero mientras las islas grandes se dispersan por cientos. Un conteo de especies no puede ser negativo, no puede ser fraccionario, y su dispersión crece con su tamaño. El modelo lineal no supone nada de eso. Este capítulo construye el modelo que sí lo hace.

Dos diagramas de dispersión lado a lado de las 30 islas Galápagos. El panel izquierdo grafica el conteo de especies contra el área de la isla con una recta de mínimos cuadrados ordinarios; la mayoría de las islas se apiñan en la esquina inferior izquierda cerca de área cero y pocas especies, un punto (Isabela) queda muy a la derecha, y la recta ajustada es arrastrada hacia arriba por él mientras pasa por encima de la mayoría de los puntos apiñados. El panel derecho grafica los mismos conteos de especies contra el área en un eje horizontal logarítmico, desplegando las islas apiñadas y mostrando una nube creciente y curvada cuya dispersión vertical se ensancha a medida que crecen los conteos de especies. — Figure 1:Los conteos de especies de las Galápagos contra el área de la isla. Izquierda: una recta de mínimos cuadrados ordinarios es distorsionada por la única isla enorme y supone una dispersión constante que los datos no tienen. Derecha: en un eje de área logarítmica las islas se despliegan y la dispersión claramente se ensancha a medida que crecen los conteos, la firma de datos de conteo que un modelo lineal no puede captar.

Al terminar este capítulo los conteos de especies tendrán un modelo que respeta lo que son, el mismo modelo manejará tasas de daño de barcos con exposiciones muy distintas, y todo ello, junto con las rectas del Capítulo 2 y las curvas logísticas del Capítulo 13, se colapsará en una sola familia con una única idea de estimación. Esa familia es el modelo lineal generalizado, y es la recompensa de toda esta parte del curso.

Esta lección de un vistazo

Qué estamos haciendo: Construir la regresión de Poisson para respuestas de conteo, luego ensamblar la familia del modelo lineal generalizado (familia más enlace) que contiene la regresión lineal, la logística y la de Poisson como una sola idea, y resolver la pregunta de las especies de las Galápagos.
Por qué lo hacemos: Los conteos de especies no pueden ser negativos ni fraccionarios y su dispersión crece con su media, así que el modelo lineal, e incluso las transformaciones del Capítulo 10, no pueden ajustarlos de forma honesta.
Objetivo principal: Enunciar el modelo log-lineal de Poisson, interpretar sus coeficientes como razones de tasas, diagnosticar la sobredispersión y aplicar la corrección cuasi-Poisson, y ensamblar el marco del MLG en una sola tabla.
Qué cambió respecto a los capítulos anteriores: El Capítulo 13 rompió por primera vez con la respuesta normal para datos de sí o no; este capítulo hace lo mismo para datos de cuántos, usando la misma maquinaria de máxima verosimilitud y devianza (13.2 Ajuste por máxima verosimilitud), y cierra el hilo de las Galápagos que quedó abierto en los Capítulos 9 y 10.

Objetivos de aprendizaje

Al terminar este capítulo serás capaz de:

Explicar por qué los mínimos cuadrados ordinarios son la herramienta equivocada para respuestas de conteo, nombrando los dos supuestos del modelo que violan.
Enunciar el modelo de regresión log-lineal de Poisson e interpretar sus coeficientes como razones de tasas, en las unidades del problema.
Derivar la log-verosimilitud de Poisson y sus ecuaciones de puntaje, y conectarlas con la maquinaria de máxima verosimilitud del Capítulo 13.
Usar un desplazamiento (offset) para modelar una tasa en vez de un conteo crudo, deduciendo por qué el logaritmo de la exposición entra en el predictor lineal con un coeficiente fijo de uno.
Diagnosticar la sobredispersión a partir del estadístico de Pearson y aplicar la corrección cuasi-Poisson, y decir cuándo un modelo binomial negativo es el mejor arreglo.
Resolver la pregunta de las especies de las Galápagos con un modelo de conteo, y explicar por qué tiene éxito donde las transformaciones del Capítulo 10 se quedaron cortas.
Ensamblar el marco del modelo lineal generalizado en una sola tabla (familia más enlace) que contiene la regresión lineal, la logística y la de Poisson como tres instancias de una única idea.

14.1 Por qué los conteos rompen el modelo lineal¶

Intuición¶

Un conteo es un número entero de cosas que ocurrieron: especies en una isla, incidentes de daño en un barco, clientes en una hora, convulsiones en una semana. Los conteos tienen tres rasgos tercos. Nunca son negativos. Son números enteros, no fracciones. Y su variabilidad crece con su promedio: una isla que promedia 5 especies varía en una especie o dos, mientras que una isla que promedia 200 varía en docenas. El modelo lineal de los Capítulos 2 al 8 supone una respuesta que puede vagar por cualquier punto de la recta numérica, toma cualquier valor real, y se dispersa con la misma $\sigma^2$ en cada nivel de los predictores. Cada uno de esos supuestos es falso para un conteo.

El modelo de probabilidad natural para un conteo es la distribución de Poisson (Poisson distribution) (Definición 14.1). Describe el número de eventos que ocurren cuando cada evento es raro, los eventos ocurren de forma independiente, y llegan a una tasa promedio estable. Una distribución de Poisson tiene un solo parámetro, su media $\mu$ , y una propiedad notable: su varianza es igual a su media. Figure 2 muestra tres distribuciones de Poisson. A medida que la media sube de 1 a 10, la distribución se desliza a la derecha, se dispersa, y se vuelve más simétrica, y en cada paso su dispersión queda atada a su centro.

Tres gráficas de barras de distribuciones de probabilidad de Poisson con medias 1, 4 y 10, dibujadas en tres colores. En media 1 las barras se amontonan en 0 y 1 y están fuertemente sesgadas a la derecha. En media 4 el pico se mueve a 3 y 4 y la forma está menos sesgada. En media 10 las barras forman un montículo ancho y casi simétrico centrado cerca de 10. Una nota en cada panel indica que la varianza es igual a la media. — Figure 2:Tres distribuciones de Poisson. Toda la distribución queda fijada por un solo número, la media, y su dispersión siempre es igual a esa media. Los conteos de media pequeña se amontonan cerca de cero y están fuertemente sesgados a la derecha; los conteos de media grande se dispersan ampliamente y se ven casi normales.

Ese vínculo media-igual-a-varianza es exactamente el patrón que mostró Figure 1 y exactamente lo que los mínimos cuadrados de varianza constante no pueden representar. Figure 3 lo hace explícito: el modelo lineal traza una línea de varianza plana, mientras que los datos de conteo cabalgan la diagonal donde la varianza sigue a la media.

Una gráfica con la respuesta media en el eje horizontal y la varianza en el eje vertical. Una línea horizontal plana etiquetada varianza constante sigma cuadrado representa el supuesto de mínimos cuadrados ordinarios. Una línea creciente a 45 grados etiquetada varianza igual a la media representa el supuesto de Poisson. Puntos dispersos que representan islas Galápagos agrupadas suben desde la parte inferior izquierda hacia la superior derecha, siguiendo la diagonal en vez de la línea plana. — Figure 3:El supuesto que falla. Los mínimos cuadrados ordinarios suponen varianza constante (línea plana); el modelo de Poisson supone que la varianza crece con la media (diagonal). Los conteos agrupados de las Galápagos suben por la diagonal, así que el modelo de varianza plana es simplemente la forma equivocada.

Fórmula¶

Un conteo aleatorio $Y$ sigue una distribución de Poisson con media $\mu > 0$ cuando

P(Y = y) = \frac{e^{-\mu}\,\mu^{y}}{y!}, \qquad y = 0, 1, 2, \dots

$Y$ es el conteo, un número entero no negativo.
$\mu$ es la media de la distribución, el conteo promedio.
$y!$ es el factorial de $y$ , que normaliza las probabilidades para que sumen uno.

En palabras: la probabilidad de ver exactamente $y$ eventos queda fijada por completo por la tasa promedio $\mu$ . La distribución acarrea dos hechos que usaremos constantemente,

E\{Y\} = \mu, \qquad \operatorname{Var}\{Y\} = \mu .

En palabras: para un conteo de Poisson la media y la varianza son el mismo número. Esta única ecuación es lo que separa la regresión de conteo de todo lo anterior, y es el hilo que seguimos hasta el final del capítulo.

Inténtalo 14.1

Un estudiante ajusta mínimos cuadrados ordinarios a los conteos de especies de las Galápagos y reporta un intervalo de predicción del 95% para una isla nueva y diminuta que va de -18 a 46 especies. Nombra las dos cosas que están mal con este intervalo, y ata cada una a un rasgo de los datos de conteo.

Solución

Primero, el intervalo incluye conteos negativos de especies, que son imposibles; un conteo está acotado por debajo por cero, pero el modelo de mínimos cuadrados pone probabilidad positiva sobre cualquier valor real. Segundo, el intervalo tiene un ancho fijo determinado por la única estimación $s$ , ignorando que una isla pequeña (media baja) debería tener un intervalo estrecho y una isla grande (media alta) uno ancho. Ambas fallas se remontan a la misma fuente: los mínimos cuadrados suponen una respuesta no acotada con varianza constante, mientras que un conteo es no negativo con varianza igual a su media.

Antes de construir un modelo para conteos, dedica un minuto a la distribución misma, porque el deslizador de abajo es el único parámetro que tiene una Poisson.

Mueve el conteo medio mu y observa cómo la distribución se desliza a la derecha, se ensancha y se aplana a la vez: con una Poisson no puedes mover el centro sin mover la dispersión, porque la varianza es la media. Volver a la Sección 14.1.

14.2 El modelo log-lineal de Poisson y las razones de tasas¶

Intuición¶

Queremos que el conteo medio $\mu_i$ dependa de los predictores, la misma ambición que la regresión ordinaria, pero con dos reparaciones. La media debe mantenerse positiva sin importar lo que hagan los predictores, y el modelo debería describir cómo los predictores multiplican el conteo en vez de sumarle, porque los procesos de conteo suelen ser multiplicativos: duplicar el área de una isla no suma un número fijo de especies, escala el número esperado.

Ambas reparaciones vienen de un solo movimiento: modelar el logaritmo de la media como una función lineal de los predictores. Como el logaritmo puede ser cualquier número real mientras que su inverso, la exponencial, siempre es positivo, un modelo lineal en la escala logarítmica puede variar libremente mientras la media que produce se mantiene positiva. Y como una suma en la escala logarítmica es un producto en la escala original, los coeficientes se vuelven multiplicadores. Figure 5 muestra el enlace en acción: una recta en el predictor lineal $\eta$ se convierte en una curva en la media $\mu$ que se pega a cero para $\eta$ negativo y sube empinada para $\eta$ positivo, sin bajar ni una vez por debajo de cero.

Dos paneles apilados que comparten un eje horizontal etiquetado predictor lineal eta. El panel superior muestra una recta creciente, el predictor lineal en sí. El panel inferior muestra la media mu igual a exp de eta como una curva que se mantiene apenas por encima de cero para eta negativo y sube empinada para eta positivo, sin cruzar nunca por debajo de cero. Una línea guía discontinua marca que eta igual a cero se aplica a mu igual a uno. — Figure 5:El enlace logarítmico. Una recta en el predictor lineal eta (arriba) se convierte en una curva exponencial en la media mu (abajo). Como mu es igual a exp de eta, la media siempre es positiva, y pasos iguales en eta multiplican la media por un factor constante.

Fórmula¶

El modelo de regresión log-lineal de Poisson dice que para cada caso $i$ ,

Y_i \sim \text{Poisson}(\mu_i), \qquad \log \mu_i = \eta_i = \beta_0 + \beta_1 X_{i1} + \dots + \beta_{p-1} X_{i,p-1} .

$Y_i$ es el conteo observado para el caso $i$ , que se supone Poisson dados los predictores.
$\mu_i = E\{Y_i\}$ es su media.
$\eta_i$ es el predictor lineal, la familiar combinación de recta de los predictores.
$\log$ es aquí la función de enlace: aplica la media positiva sobre toda la recta real donde vive el predictor lineal.
$\beta_0, \dots, \beta_{p-1}$ son los coeficientes de regresión, estimados por $b_0, \dots, b_{p-1}$ .

Despejar la media del enlace da la forma multiplicativa

\mu_i = e^{\eta_i} = e^{\beta_0}\, e^{\beta_1 X_{i1}} \cdots e^{\beta_{p-1} X_{i,p-1}} .

En palabras: el conteo medio es un producto de factores, uno por predictor. Por esto los coeficientes se leen como razones de tasas (rate ratios). Aumenta $X_{ij}$ en una unidad, manteniendo fijos los demás, y $\eta_i$ sube en $\beta_j$ , así que $\mu_i$ se multiplica por $e^{\beta_j}$ :

\frac{\mu(X_j + 1)}{\mu(X_j)} = e^{\beta_j} .

En palabras: $e^{\beta_j}$ es el factor por el cual cambia el conteo esperado ante un aumento de una unidad en $X_j$ . Un coeficiente de cero significa una razón de uno, sin efecto; un coeficiente positivo significa una razón mayor que uno, el conteo crece; un coeficiente negativo significa una razón menor que uno, el conteo se encoge. Este es el primo, para datos de conteo, de la razón de momios $e^{\beta_j}$ que conociste para la regresión logística en 13.3 Leer los coeficientes como razones de momios, y de la lectura de cambio porcentual de un coeficiente en escala logarítmica de 10.2 La transformación logarítmica y la lectura de sus coeficientes.

Esta es la imagen que hay que tener en la cabeza. En la escala logarítmica el modelo es una recta ordinaria, así que pasos iguales en un predictor suman la misma cantidad $\beta_j$ cada vez. Exponencia, y ese sumar parejo se convierte en multiplicar parejo: cada paso multiplica el conteo por el mismo factor $e^{\beta_j}$ . Figure 6 muestra ambos lados a la vez, usando el efecto de la elevación del ajuste de más abajo, donde cien metros extra multiplican el conteo esperado de especies por cerca de 1.42.

Dos gráficas de líneas lado a lado que comparten un eje horizontal de elevación en cientos de metros, de 0 a 4. El panel izquierdo grafica el logaritmo del conteo medio y es una recta creciente; etiquetas naranjas igualmente espaciadas que dicen más beta marcan que cada paso hacia arriba suma la misma cantidad. El panel derecho grafica el conteo medio en sí y es una curva que sube cada vez más empinada; etiquetas naranjas que dicen por 1.42 marcan que cada paso hacia arriba multiplica el conteo por el mismo factor, así que las brechas entre puntos crecen. — Figure 6:Un coeficiente es una razón de tasas. En la escala logarítmica (izquierda) cada paso de elevación suma la misma cantidad al predictor lineal, una recta. En la escala de conteo (derecha) ese mismo paso multiplica el conteo esperado por un factor fijo de cerca de 1.42, así que el conteo sube cada vez más rápido. Sumar en la escala logarítmica es multiplicar en la escala de conteo.

Derivación (la log-verosimilitud de Poisson y sus ecuaciones de puntaje)¶

Los mínimos cuadrados no estimarán este modelo: la respuesta no es normal y la media no es lineal. Como en la regresión logística, estimamos por máxima verosimilitud (maximum likelihood), el principio de 13.2 Ajuste por máxima verosimilitud que elige los coeficientes que hacen más probables los conteos observados.

Demostración. Los conteos $Y_1, \dots, Y_n$ son independientes dados los predictores, así que la verosimilitud es el producto de sus probabilidades de Poisson, y la log-verosimilitud es la suma de los logaritmos:

\ell(\boldsymbol\beta) = \sum_{i=1}^n \Big[ y_i \log \mu_i - \mu_i - \log(y_i!) \Big] .

Sustituye el modelo $\log \mu_i = \mathbf{x}_i'\boldsymbol\beta$ , de modo que $\mu_i = e^{\mathbf{x}_i'\boldsymbol\beta}$ , donde $\mathbf{x}_i = (1, X_{i1}, \dots, X_{i,p-1})'$ es el vector de predictores para el caso $i$ :

\ell(\boldsymbol\beta) = \sum_{i=1}^n \Big[ y_i\, \mathbf{x}_i'\boldsymbol\beta - e^{\mathbf{x}_i'\boldsymbol\beta} - \log(y_i!) \Big] .

Para maximizar, deriva con respecto al vector de coeficientes e iguala el resultado a cero. El último término no involucra a $\boldsymbol\beta$ . Para los dos primeros, usando $\partial (\mathbf{x}_i'\boldsymbol\beta)/\partial \boldsymbol\beta = \mathbf{x}_i$ y $\partial e^{\mathbf{x}_i'\boldsymbol\beta}/\partial \boldsymbol\beta = e^{\mathbf{x}_i'\boldsymbol\beta}\mathbf{x}_i = \mu_i \mathbf{x}_i$ ,

\frac{\partial \ell}{\partial \boldsymbol\beta} = \sum_{i=1}^n \Big( y_i \mathbf{x}_i - \mu_i \mathbf{x}_i \Big) = \sum_{i=1}^n (y_i - \mu_i)\,\mathbf{x}_i = \mathbf{0} .

Estas son las ecuaciones de puntaje (score equations) de la regresión de Poisson. Lee la primera componente (la del intercepto, donde $\mathbf{x}_i$ tiene un 1): dice $\sum_i (y_i - \hat\mu_i) = 0$ , las medias ajustadas reproducen el conteo total. Cada otra componente dice $\sum_i (y_i - \hat\mu_i) X_{ij} = 0$ , los residuos $y_i - \hat\mu_i$ son ortogonales a cada predictor. Estos son los análogos de las ecuaciones normales de 2.2 Mínimos cuadrados desde los primeros principios, con la media ajustada $\hat\mu_i$ desempeñando el papel del valor ajustado. La diferencia es que $\hat\mu_i = e^{\mathbf{x}_i'\mathbf{b}}$ es no lineal en los coeficientes, así que a diferencia de los mínimos cuadrados no hay solución en forma cerrada. Las ecuaciones se resuelven por los mismos mínimos cuadrados reponderados iterativamente que ajustan la regresión logística, esbozados en 13.2 Ajuste por máxima verosimilitud; el software lo ejecuta en un puñado de pasos. $\blacksquare$

Fíjate en lo poco que cambió respecto a la regresión logística. Allá las ecuaciones de puntaje eran $\sum_i (y_i - \hat\pi_i)\mathbf{x}_i = \mathbf{0}$ con $\hat\pi_i$ una probabilidad; aquí son $\sum_i (y_i - \hat\mu_i)\mathbf{x}_i = \mathbf{0}$ con $\hat\mu_i$ un conteo medio. Misma forma, distinto modelo de la media. La Sección 14.6 muestra que esto no es un accidente.

R¶

La función de ajuste es glm, la misma que se usa para la regresión logística, con family = poisson. El enlace por defecto para la familia de Poisson es el logarítmico, así que rara vez lo nombras.

gala <- read.csv("data/gala.csv")
pois_fit <- glm(Species ~ Area + Elevation + Nearest + Scruz + Adjacent,
                family = poisson, data = gala)
summary(pois_fit)

Call:
glm(formula = Species ~ Area + Elevation + Nearest + Scruz +
    Adjacent, family = poisson, data = gala)

Coefficients:
              Estimate Std. Error z value Pr(>|z|)
(Intercept)  3.155e+00  5.175e-02  60.963  < 2e-16 ***
Area        -5.799e-04  2.627e-05 -22.074  < 2e-16 ***
Elevation    3.541e-03  8.741e-05  40.507  < 2e-16 ***
Nearest      8.826e-03  1.821e-03   4.846 1.26e-06 ***
Scruz       -5.709e-03  6.256e-04  -9.126  < 2e-16 ***
Adjacent    -6.630e-04  2.933e-05 -22.608  < 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: 3510.73  on 29  degrees of freedom
Residual deviance:  716.85  on 24  degrees of freedom
AIC: 889.68

Number of Fisher Scoring iterations: 5

Ejemplo 14.1: Leer los coeficientes de las Galápagos como razones de tasas

Pregunta. En el ajuste de Poisson, ¿qué dice el coeficiente de elevación sobre los conteos de especies, en unidades sencillas que usaría un biólogo?

Intuición. En la escala logarítmica el coeficiente es una pendiente aditiva; exponéncialo para convertirlo en un multiplicador del conteo esperado por unidad de elevación.

Fórmula. La razón de tasas para un aumento de una unidad en un predictor es $e^{b_j}$ ; para un aumento de $c$ unidades es $e^{c\,b_j}$ .

Cálculo.

b <- coef(pois_fit)
c(per_meter = exp(b["Elevation"]),
  per_100m  = exp(100 * b["Elevation"]))

per_meter.Elevation  per_100m.Elevation
           1.003547            1.424840

import numpy as np
import pandas as pd
import statsmodels.formula.api as smf
import statsmodels.api as sm

gala = pd.read_csv("data/gala.csv")
pois_fit = smf.glm("Species ~ Area + Elevation + Nearest + Scruz + Adjacent",
                   data=gala, family=sm.families.Poisson()).fit()
b = pois_fit.params
print(np.exp(b["Elevation"]), np.exp(100 * b["Elevation"]))

1.0035468693157312 1.4248398207117439

Interpretación. Cada metro adicional de elevación multiplica el conteo esperado de especies por cerca de 1.0035, un aumento de aproximadamente un tercio de por ciento por metro. Eso se siente más fácil sobre un paso más grande: cada 100 metros extra de elevación multiplican el conteo esperado por cerca de 1.42, un aumento del 42%, manteniendo fijos los demás predictores. El coeficiente es una razón de tasas, no un número aditivo de especies, que es exactamente el comportamiento multiplicativo que exigen los datos de conteo.

Inténtalo 14.2

El coeficiente de Adjacent es $b = -6.630 \times 10^{-4}$ por kilómetro cuadrado. Calcula la razón de tasas para un aumento de 100 km $^2$ en el área de la isla adyacente, y di en palabras qué significa para el conteo esperado de especies.

Solución

La razón de tasas para un aumento de 100 unidades es $e^{100 b} = e^{100 \times (-0.0006630)} = e^{-0.06630} = 0.9358$ . Cada 100 km $^2$ extra de área de isla adyacente multiplican el conteo esperado de especies por cerca de 0.936, una caída de aproximadamente 6.4%, manteniendo fijos los demás predictores. Un vecino grande se asocia con menos especies, consistente con una isla adyacente grande que atrae colonizadores o compite por ellos.

Guiar el enlace logarítmico a mano es la forma más rápida de sentir lo que hace la máxima verosimilitud, así que toma tú los dos coeficientes y deja que la log-verosimilitud te avise cuando estés cerca.

Mueve b0 y b1 para ajustar a mano los conteos de especies de las Galápagos contra la elevación. La curva exp(b0 + b1 x) nunca puede bajar de cero, la razón de tasas que se muestra es exp(b1), el multiplicador por cada 100 metros, y la banda ámbar es la dispersión de Poisson dentro de la cual se supone que deben quedar los conteos. Volver a la Sección 14.2.

14.3 Desplazamientos: modelar una tasa en vez de un conteo¶

Intuición¶

Los conteos crudos solo son comparables cuando vienen de igualdad de oportunidad. Un buque de carga en el mar durante 40,000 meses tiene mucha más probabilidad de sufrir daño por oleaje que uno en el mar durante 100 meses, así que comparar sus conteos crudos de incidentes carece de sentido. Lo que queremos es la tasa: incidentes por mes de servicio. La exposición, meses de servicio, difiere de barco a barco y debe integrarse al modelo, no ignorarse.

El arreglo es elegante. Modela la tasa como log-lineal, luego reescríbela como un modelo del conteo con una pieza fija añadida al predictor lineal. Esa pieza fija, el logaritmo de la exposición, se llama desplazamiento (offset) (Definición 14.4): un predictor cuyo coeficiente no se estima sino que se fija en uno, porque sabemos exactamente cómo debería escalar el conteo con la exposición. Duplica los meses en el mar y, en igualdad de condiciones, esperas el doble de incidentes. Figure 8 muestra por qué los conteos crudos engañan: el tipo de barco B tiene por mucho el mayor número de incidentes, pero también por mucho el mayor servicio; por mes de servicio su tasa de daño es ordinaria.

Dos gráficas de barras de los cinco tipos de barco A a E. La gráfica izquierda muestra el total de incidentes de daño; el tipo B se eleva por encima de los demás con bastante más de cien incidentes mientras el resto son mucho más pequeños. La gráfica derecha muestra incidentes por mil meses de servicio; ahora el tipo B es poco notable e intermedio, mientras los tipos A y E muestran las tasas de daño más altas, invirtiendo el orden de la gráfica izquierda. — Figure 8:Por qué importa la exposición. Izquierda: los incidentes de daño crudos hacen que el tipo de barco B parezca por mucho el peor. Derecha: una vez que dividimos por meses de servicio para obtener una tasa, el tipo B es promedio y los tipos A y E son los preocupantes. El desplazamiento es lo que le permite al modelo comparar tasas en vez de conteos crudos.

Fórmula¶

Sea $t_i$ la exposición para el caso $i$ (meses de servicio, personas-año, área inspeccionada) y sea $\lambda_i$ la tasa de eventos por unidad de exposición. La media del conteo es la tasa por la exposición, $\mu_i = \lambda_i t_i$ . Modela el logaritmo de la tasa como lineal:

\log \lambda_i = \beta_0 + \beta_1 X_{i1} + \dots + \beta_{p-1} X_{i,p-1} .

Derivación (el desplazamiento). Como $\mu_i = \lambda_i t_i$ , toma logaritmos:

\log \mu_i = \log \lambda_i + \log t_i = \underbrace{\beta_0 + \beta_1 X_{i1} + \dots + \beta_{p-1} X_{i,p-1}}_{\text{estimado}} + \underbrace{\log t_i}_{\text{desplazamiento fijo}} .

En palabras: modelar el logaritmo de la tasa es lo mismo que modelar el logaritmo del conteo con $\log t_i$ añadido al predictor lineal como un término con coeficiente conocido exactamente uno. Ese término extra es el desplazamiento. Desplaza el predictor lineal de cada caso por su logaritmo de exposición sin costar un parámetro, así que los coeficientes estimados describen la tasa $\lambda_i$ , y $e^{\beta_j}$ es ahora una razón de tasas por unidad de exposición. Si erróneamente quitaras el desplazamiento, el modelo explicaría los incidentes en parte a través de cualquier predictor que resulte correlacionado con el tiempo de servicio, confundiendo la tasa con la exposición. $\blacksquare$

R¶

Los barcos con cero servicio registrado no aportan exposición ni oportunidad de incidentes, así que $\log t_i$ está indefinido para ellos; primero descartamos esas filas. El desplazamiento entra a través del argumento offset.

ships <- read.csv("data/ships.csv")
ships <- ships[ships$service > 0, ]
ships$year   <- factor(ships$year)
ships$period <- factor(ships$period)
ship_fit <- glm(incidents ~ type + year + period,
                family = poisson, offset = log(service), data = ships)
round(summary(ship_fit)$coefficients, 4)

            Estimate Std. Error  z value Pr(>|z|)
(Intercept)  -6.4059     0.2174 -29.4600   0.0000
typeB        -0.5433     0.1776  -3.0595   0.0022
typeC        -0.6874     0.3290  -2.0891   0.0367
typeD        -0.0760     0.2906  -0.2614   0.7938
typeE         0.3256     0.2359   1.3803   0.1675
year65        0.6971     0.1496   4.6587   0.0000
year70        0.8184     0.1698   4.8207   0.0000
year75        0.4534     0.2332   1.9446   0.0518
period75      0.3845     0.1183   3.2507   0.0012

Ejemplo 14.2: Razones de tasas de daño de los barcos

Pregunta. Manteniendo fijos el año de construcción y el periodo de operación, ¿cómo se compara la tasa de daño por oleaje del barco tipo B con la del tipo A de referencia? ¿Y cuánto cambia la tasa operar en el periodo posterior (1975 a 1979)?

Intuición. Los coeficientes están en la escala del logaritmo de la tasa, así que exponéncialos para obtener razones de tasas. El tipo A es el nivel de referencia, así que su razón de tasas es uno por construcción.

Fórmula. Razón de tasas para un nivel de factor $= e^{b}$ relativa al nivel de referencia.

Cálculo.

round(exp(coef(ship_fit)), 3)

(Intercept)       typeB       typeC       typeD       typeE      year65
      0.002       0.581       0.503       0.927       1.385       2.008
     year70      year75    period75
      2.267       1.574       1.469

ships = pd.read_csv("data/ships.csv")
ships = ships[ships["service"] > 0].copy()
ship_fit = smf.glm("incidents ~ C(type) + C(year) + C(period)", data=ships,
                   family=sm.families.Poisson(),
                   offset=np.log(ships["service"])).fit()
print(np.exp(ship_fit.params).round(3))

Intercept          0.002
C(type)[T.B]       0.581
C(type)[T.C]       0.503
C(type)[T.D]       0.927
C(type)[T.E]       1.385
C(year)[T.65]      2.008
C(year)[T.70]      2.267
C(year)[T.75]      1.574
C(period)[T.75]    1.469
dtype: float64

Interpretación. La razón de tasas del tipo B es 0.581: tras ajustar por exposición, época de construcción y periodo de operación, los barcos tipo B sufren daño a cerca del 58% de la tasa del tipo A, no la más alta como sugerían los conteos crudos en Figure 8. La razón de tasas de period75 es 1.469: los barcos operados en el periodo posterior tuvieron cerca de 47% más incidentes de daño por mes de servicio que en el periodo anterior, manteniendo fijos el tipo y el año de construcción. El intercepto exponencia a 0.0017 incidentes por mes para un barco de referencia, por lo cual la escala práctica es incidentes por mil meses.

Inténtalo 14.3

Un analista tiene datos de crimen de la ciudad: crimes (número de crímenes reportados) y population para cada ciudad, y quiere tasas de crimen por persona. Escribe la fórmula de R y el argumento extra para un modelo de Poisson de la tasa de crimen sobre un predictor unemployment. ¿Por qué no modelar simplemente crimes / population directamente con mínimos cuadrados?

Solución

El modelo es glm(crimes ~ unemployment, family = poisson, offset = log(population), data = ...). El desplazamiento log(population) hace que los coeficientes describan la tasa por persona. Modelar crimes / population con mínimos cuadrados es peor por dos razones: la razón no es normal ni de varianza constante (una tasa de un pueblo pequeño es mucho más ruidosa que una de una ciudad grande, pero los mínimos cuadrados las ponderarían por igual), y los valores ajustados podrían volverse negativos. El enfoque del desplazamiento mantiene el conteo como una respuesta de Poisson, así que la varianza escala correctamente con la cantidad de exposición detrás de cada ciudad.

Un desplazamiento se cree más fácil cuando ya probaste las alternativas, así que pon el exponente de exposición bajo tu dedo y mira qué afirma cada ajuste sobre los barcos.

Pon el exponente de exposición en 0 y el modelo ignora del todo los meses en el mar, prediciendo el mismo puñado de incidentes para un barco con 45 meses de servicio y para uno con 44,882. Ponlo en 1 y duplicar el servicio duplica el conteo esperado, que es exactamente lo que fija el desplazamiento log t. Volver a la Sección 14.3.

14.4 Sobredispersión y la corrección cuasi-Poisson¶

Intuición¶

El modelo de Poisson hace una promesa fuerte: la varianza es igual a la media. Los datos de conteo reales a menudo rompen esa promesa dispersándose más de lo que Poisson permite, una condición llamada sobredispersión (overdispersion) (Definición 14.5). Sus causas son ordinarias: predictores que no mediste, agrupamiento de eventos, diferencias de individuo a individuo en la tasa subyacente. Sea cual sea la causa, el síntoma es el mismo, los conteos se dispersan más ancho que sus medias ajustadas.

La sobredispersión rara vez cambia mucho las estimaciones de los coeficientes, pero arruina sus errores estándar. Un ajuste de Poisson que supone varianza igual a la media, cuando la varianza verdadera es varias veces mayor, reporta errores estándar varias veces demasiado pequeños, convirtiendo el ruido en falsa significancia. Así que la sobredispersión no es una molestia para encogerse de hombros; es la diferencia entre un análisis honesto y uno engañoso. Figure 10 muestra el diagnóstico para el ajuste de las Galápagos: los residuos al cuadrado, que deberían rondar la media ajustada bajo un verdadero Poisson, quedan muy por encima de la línea de varianza igual a la media.

Un diagrama de dispersión para el ajuste de Poisson de las Galápagos con la media ajustada en el eje horizontal y la desviación al cuadrado en el eje vertical, ambos en escalas logarítmicas. Una línea de referencia diagonal marca donde la desviación al cuadrado es igual a la media ajustada, la promesa de Poisson. Casi toda isla queda bastante por encima de la línea, muchas por uno o dos órdenes de magnitud, mostrando que los conteos varían mucho más de lo que Poisson permite. — Figure 10:Diagnosticar la sobredispersión en el ajuste de las Galápagos. Bajo un verdadero Poisson las desviaciones al cuadrado se dispersarían alrededor de la diagonal (varianza igual a la media). En cambio casi toda isla queda muy por encima de ella, así que los datos varían mucho más de lo que el modelo de Poisson afirma, y los errores estándar del modelo no pueden confiarse tal como están.

Fórmula¶

Mide la dispersión con el estadístico de Pearson (Pearson statistic) (Definición 14.6) dividido entre sus grados de libertad. Define el residuo de Pearson para el caso $i$ y la dispersión estimada $\hat\phi$ :

r_i^{P} = \frac{Y_i - \hat\mu_i}{\sqrt{\hat\mu_i}}, \qquad \hat\phi = \frac{1}{n-p}\sum_{i=1}^n \big(r_i^{P}\big)^2 = \frac{X^2}{n-p} .

$r_i^{P}$ estandariza cada residuo por la desviación estándar de Poisson $\sqrt{\hat\mu_i}$ .
$X^2 = \sum_i (r_i^{P})^2$ es el estadístico ji-cuadrado de Pearson.
$n - p$ son los grados de libertad residuales, casos menos parámetros estimados.

En palabras: $\hat\phi$ es el residuo estandarizado al cuadrado promedio. Bajo un modelo de Poisson correcto debería estar cerca de 1. Un valor bastante por encima de 1 señala sobredispersión; un valor de $\hat\phi$ alrededor de 2 significa que la varianza verdadera es aproximadamente el doble de la media, y así sucesivamente.

La corrección cuasi-Poisson mantiene intactos el modelo de la media de Poisson y sus estimaciones de coeficientes, porque las ecuaciones de puntaje $\sum_i (y_i - \mu_i)\mathbf{x}_i = \mathbf{0}$ no involucran la varianza. Solo reescala la varianza, suponiendo

\operatorname{Var}\{Y_i\} = \phi\, \mu_i

para una constante desconocida $\phi > 1$ , estimada por $\hat\phi$ . Cada error estándar se multiplica entonces por $\sqrt{\hat\phi}$ :

s_{\text{quasi}}\{b_j\} = \sqrt{\hat\phi}\; s_{\text{Poisson}}\{b_j\} .

En palabras: infla cada error estándar de Poisson por la raíz cuadrada de la dispersión estimada. Esta es la misma idea que viste en 10.4 Mínimos cuadrados ponderados, donde modelar la varianza no constante cambió los pesos y los errores estándar mientras dejaba intacta la historia de la media. Aquí se permite que la varianza sea un múltiplo constante de la media en vez de exactamente la media.

R¶

glm reporta el supuesto de dispersión crudo, así que calcula $\hat\phi$ directamente a partir de los residuos de Pearson, luego reajusta con family = quasipoisson, que hace el reescalado por ti.

phi_hat <- sum(residuals(pois_fit, type = "pearson")^2) / df.residual(pois_fit)
phi_hat

[1] 31.74914

Ejemplo 14.3: Diagnosticar y corregir la sobredispersión en las Galápagos

Pregunta. ¿Está sobredispersado el ajuste de Poisson de las Galápagos, y qué le pasa al error estándar del coeficiente de elevación una vez que lo tomamos en cuenta?

Intuición. Calcula $\hat\phi$ ; si está muy por encima de 1, los errores estándar de Poisson son demasiado pequeños por un factor de $\sqrt{\hat\phi}$ . Reajusta como cuasi-Poisson y compara.

Fórmula. $\hat\phi = X^2/(n-p)$ , y $s_{\text{quasi}} = \sqrt{\hat\phi}\, s_{\text{Poisson}}$ .

Cálculo.

quasi_fit <- glm(Species ~ Area + Elevation + Nearest + Scruz + Adjacent,
                 family = quasipoisson, data = gala)
round(summary(quasi_fit)$coefficients, 5)

             Estimate Std. Error  t value Pr(>|t|)
(Intercept)  3.15481    0.29159 10.81933  0.00000
Area        -0.00058    0.00015 -3.91750  0.00065
Elevation    0.00354    0.00049  7.18891  0.00000
Nearest      0.00883    0.01026  0.86001  0.39829
Scruz       -0.00571    0.00353 -1.61963  0.11838
Adjacent    -0.00066    0.00017 -4.01229  0.00051

quasi_fit = smf.glm("Species ~ Area + Elevation + Nearest + Scruz + Adjacent",
                    data=gala, family=sm.families.Poisson()).fit(scale="X2")
print(quasi_fit.scale)
print(quasi_fit.bse["Elevation"])

31.749135323558495
0.0004925072249019893

Interpretación. La dispersión es $\hat\phi = 31.7$ , no el 1 que exige un verdadero Poisson: los conteos de especies varían cerca de 32 veces más de lo que Poisson permite. El error estándar de la elevación crece de $8.74 \times 10^{-5}$ bajo Poisson a $4.93 \times 10^{-4}$ bajo cuasi-Poisson, un factor de $\sqrt{31.7} \approx 5.6$ . El coeficiente en sí no cambia, se queda en 0.00354, pero ahora se conoce cerca de 5.6 veces con menos precisión, y los predictores débiles Nearest y Scruz, significativos bajo el Poisson ingenuo, ya no son significativos. El análisis cuasi-Poisson es el honesto.

Figure 11 lleva el punto a casa a lo largo de los cinco predictores: dividir cada estadístico de prueba entre $\sqrt{\hat\phi}$ jala dos de ellos por debajo del umbral aproximado de significancia que el Poisson ingenuo les hacía superar con holgura.

Una gráfica de barras agrupadas con los cinco predictores de las Galápagos Area, Elevation, Nearest, Scruz y Adjacent en el eje horizontal. Para cada predictor dos barras muestran el valor absoluto del estadístico de prueba, una barra de Poisson más alta y una barra cuasi-Poisson mucho más corta, cerca de 5.6 veces más pequeña. Una línea horizontal discontinua en 2 marca el umbral aproximado de significancia. Bajo Poisson las cinco barras se elevan muy por encima de la línea; bajo cuasi-Poisson, Nearest y Scruz caen por debajo de la línea mientras Area, Elevation y Adjacent se quedan por encima. — Figure 11:El costo de ignorar la sobredispersión. Los estadísticos de prueba de Poisson ingenuo (barras altas) están todos muy por encima del corte aproximado de significancia de 2 (línea discontinua). Dividir entre la raíz cuadrada de la dispersión da los estadísticos cuasi-Poisson honestos (barras cortas), y dos predictores, Nearest y Scruz, caen por debajo del corte.

Una mención al binomial negativo¶

El cuasi-Poisson es un parche rápido y ligero en supuestos: no nombra una distribución de probabilidad, solo reescala varianzas. Cuando quieres una verosimilitud genuina, y por lo tanto AIC, pruebas de razón de verosimilitud, e intervalos de predicción honestos, la elección estándar es el modelo binomial negativo (negative binomial). Añade un parámetro $\theta$ y permite que la varianza crezca cuadráticamente, $\operatorname{Var}\{Y\} = \mu + \mu^2/\theta$ , lo que ajusta mejor la sobredispersión fuerte que la varianza estrictamente lineal del cuasi-Poisson. En R es MASS::glm.nb; los datos de las Galápagos dan $\hat\theta = 1.67$ , lo bastante pequeño para confirmar sobredispersión severa. Para este curso, diagnostica con $\hat\phi$ , echa mano del cuasi-Poisson como el arreglo de todos los días, y ten presente que el binomial negativo es el siguiente paso cuando necesitas una verosimilitud completa.

Inténtalo 14.4

Un ajuste de Poisson reporta un error estándar de la pendiente de 0.04 y una dispersión de Pearson de $\hat\phi = 9$ . ¿Cuál es el error estándar cuasi-Poisson honesto, y por qué factor se encoge el estadístico de prueba basado en $z$ ?

Solución

El error estándar cuasi-Poisson es $\sqrt{9}\times 0.04 = 3 \times 0.04 = 0.12$ . Como el estadístico de prueba es el coeficiente dividido entre su error estándar, triplicar el error estándar divide el estadístico entre 3. Un predictor que parecía fuertemente significativo con un $z$ de, digamos, 6 bajo Poisson cae a un $t$ de 2 bajo cuasi-Poisson, una conclusión completamente distinta. Por esto exactamente la sobredispersión debe revisarse antes de cualquier afirmación de significancia.

La dispersión es solo un número hasta que la ves mover un valor p, así que arrástrala tú y mira qué predictores de las Galápagos sobreviven a la corrección.

Sube la dispersión hasta que la línea ámbar cuasi-Poisson por fin supere la nube de islas, cerca de 32, y observa cómo los estadísticos de prueba de Nearest y Scruz bajan de 2 mientras cada error estándar se infla por la raíz cuadrada de phi. Volver a la Sección 14.4.

14.5 El hilo de las Galápagos, resuelto¶

Ahora podemos cerrar la historia que ha corrido desde el Capítulo 9. Recuerda las dos visitas anteriores. El Capítulo 9 hizo de los datos de las Galápagos su lección sobre influencia: Isabela, la isla más grande, es el punto de alto apalancamiento al que apunta 9.1 Los valores de apalancamiento y la matriz sombrero, uno que un ajuste de mínimos cuadrados ordinarios no puede evitar perseguir, porque queda solo en el borde más lejano del espacio de predictores. En 10.2 La transformación logarítmica y la lectura de sus coeficientes y 10.3 Box-Cox: dejar que los datos elijan la potencia probaste el remedio clásico, transformar la respuesta, tomando logaritmos o una potencia de Box-Cox para domar el sesgo y la dispersión creciente. Mejoró las gráficas de residuos, pero nunca quitó la incomodidad, porque el problema real nunca fue la escala de una respuesta continua. La respuesta era un conteo, y ninguna transformación de un conteo lo convierte en una medición normal de varianza constante.

El modelo de Poisson corrige el defecto real. Integra la no negatividad, la estructura multiplicativa, y el vínculo media-varianza que los conteos de especies genuinamente tienen, y una vez que tomamos en cuenta la sobredispersión sus errores estándar son honestos. Figure 13 pone las dos filosofías lado a lado: la curva de log-transformar-luego-mínimos-cuadrados del Capítulo 10, y la curva de la media cuasi-Poisson de este capítulo, ambas contra los conteos crudos en un eje de área logarítmica. Coinciden en el medio y se separan en los extremos, donde el modelo de conteo mantiene las predicciones positivas y deja que la dispersión se ensanche como debe.

Un diagrama de dispersión de los conteos de especies de las Galápagos contra el área de la isla en un eje horizontal logarítmico, con dos curvas ajustadas superpuestas. Una curva discontinua etiquetada mínimos cuadrados con log-transformación del Capítulo 10 y una curva sólida etiquetada Poisson del Capítulo 14 ambas suben con el área y siguen los datos de cerca por el medio. En las islas más pequeñas la curva de mínimos cuadrados sobre logaritmos se dobla de forma distinta mientras la curva de Poisson se mantiene positiva y se pega a los conteos bajos; la dispersión de los datos es visiblemente más apretada en los conteos pequeños y más ancha en los conteos grandes. — Figure 13:La resolución. La curva de transformar-y-ajustar del Capítulo 10 y el modelo de conteo del Capítulo 14 coinciden donde los datos son densos y divergen en los extremos. El modelo de Poisson está construido para lo que la respuesta es, un conteo, así que mantiene las predicciones positivas y coincide con la forma en que la dispersión crece con la media.

La lección es más grande que un archipiélago. Las transformaciones del Capítulo 10 son la herramienta correcta cuando una respuesta continua está en una escala incómoda. Cuando la respuesta es una clase de cosa fundamentalmente distinta, un conteo, un sí o no, una tasa, no doblas la respuesta para que encaje en el modelo lineal; cambias el modelo. Ese cambio es el tema de la última sección, y es la razón por la que la regresión logística y la de Poisson pertenecen al mismo curso que los mínimos cuadrados.

Así que la primera decisión en cualquier análisis no es qué predictores usar; es qué clase de cosa es la respuesta. Figure 14 convierte esa decisión en un diagrama de flujo breve: nombra el tipo de respuesta, sigue la flecha, y aterrizas en el modelo correspondiente. La siguiente sección muestra por qué los cuatro destinos son en realidad una sola máquina.

Un diagrama de flujo de decisión. Una caja en la parte superior pregunta qué clase de cosa es la respuesta. Cuatro flechas se abren hacia abajo a cuatro cajas de tipo de respuesta: continua y aproximadamente simétrica; continua, positiva y sesgada; un conteo 0, 1, 2, 3 y así sucesivamente; y sí o no o una proporción. Cada caja de respuesta apunta hacia abajo a su modelo: regresión lineal con una distribución Normal y enlace identidad; una transformación del Capítulo 10 o un MLG gamma; regresión de Poisson con enlace logarítmico; y regresión logística con una distribución Binomial y enlace logit. Debajo de la caja de Poisson cuelgan dos cajas más: una que dice si la respuesta es una tasa, añade un desplazamiento de log t por exposición, y una que dice si los conteos están demasiado dispersos, revisa phi-sombrero y usa cuasi-Poisson. — Figure 14:Diagnostica la respuesta, luego elige el modelo. Comienza arriba, decide qué clase de cantidad es la respuesta, y el diagrama de flujo te dirige al modelo correspondiente. Los conteos llevan a la regresión de Poisson, con un desplazamiento cuando la respuesta es en realidad una tasa y una corrección cuasi-Poisson cuando los conteos están sobredispersos.

14.6 Una familia: el modelo lineal generalizado¶

Intuición¶

Da un paso atrás y mira lo que comparten los últimos tres capítulos. La regresión ordinaria modela una media continua como una recta. La regresión logística modela una probabilidad pasando una recta a través del enlace logit. La regresión de Poisson modela un conteo pasando una recta a través del enlace logarítmico. En cada caso hay un predictor lineal $\eta = \mathbf{x}'\boldsymbol\beta$ haciendo el mismo trabajo; lo que cambia son solo dos elecciones: qué distribución de probabilidad sigue la respuesta, y qué función de enlace conecta su media con el predictor lineal. Fija esas dos elecciones y todo lo demás, la estimación por máxima verosimilitud, la devianza, el flujo de trabajo, es maquinaria compartida. Ese es el modelo lineal generalizado (generalized linear model), o MLG (Definición 14.8).

Fórmula¶

Un modelo lineal generalizado tiene tres partes:

\underbrace{Y_i \sim \text{una distribución de la familia exponencial}}_{\text{componente aleatorio}}, \qquad \underbrace{\eta_i = \mathbf{x}_i'\boldsymbol\beta}_{\text{predictor lineal}}, \qquad \underbrace{g(\mu_i) = \eta_i}_{\text{enlace}} .

El componente aleatorio nombra la distribución de la respuesta (normal, binomial, Poisson, y otras), que fija cómo depende la varianza de la media.
El predictor lineal $\eta_i$ es la misma combinación lineal de predictores usada durante todo el libro.
La función de enlace $g$ aplica la media $\mu_i$ sobre la escala del predictor lineal; su inversa la aplica de vuelta, $\mu_i = g^{-1}(\eta_i)$ .

En palabras: elige una distribución y un enlace, y has especificado un modelo; el mismo motor de estimación los ajusta todos. Los tres modelos de esta parte son tres filas de una sola tabla.

Modelo	Respuesta	Distribución	Enlace $g(\mu)$	Modelo de la media	$\operatorname{Var}\{Y\}$	Capítulo
Lineal	continua	Normal	identidad $\mu$	$\mu = \mathbf{x}'\boldsymbol\beta$	$\sigma^2$ (constante)	2 a 8
Logística	binaria / proporción	Binomial	logit $\log\frac{\mu}{1-\mu}$	$\mu = \dfrac{e^{\mathbf{x}'\boldsymbol\beta}}{1+e^{\mathbf{x}'\boldsymbol\beta}}$	$\mu(1-\mu)$	13
Poisson	conteo	Poisson	log $\log\mu$	$\mu = e^{\mathbf{x}'\boldsymbol\beta}$	$\mu$	14

Figure 15 dibuja la misma idea como un diagrama: un predictor lineal alimentando tres enlaces hacia tres tipos de respuesta.

Un diagrama esquemático. A la izquierda, una sola caja etiquetada predictor lineal eta igual a x-transpuesta beta. Una flecha se abre hacia tres cajas de enlace apiladas verticalmente: identidad, logit y log. Cada caja de enlace apunta a una caja de respuesta a la derecha: identidad apunta a una respuesta continua con una distribución normal y varianza constante, logit apunta a una respuesta binaria con una distribución binomial y varianza mu por uno menos mu, log apunta a una respuesta de conteo con una distribución de Poisson y varianza igual a mu. Una franja de leyenda debajo dice un predictor lineal, un motor de estimación, tres familias. — Figure 15:El modelo lineal generalizado en una imagen. Un solo predictor lineal alimenta tres enlaces distintos hacia tres tipos de respuesta distintos. Cambia el enlace y la familia y te mueves entre la regresión lineal, la logística y la de Poisson sin cambiar la maquinaria de estimación subyacente.

R¶

La unificación no es solo conceptual; es literalmente cómo está construido el software. La función glm ajusta cada fila de la tabla, y elegir family = gaussian(link = "identity") reproduce los mínimos cuadrados ordinarios exactamente.

Ejemplo 14.4: Los mínimos cuadrados son un MLG

Pregunta. ¿Ajustar los datos de Toluca con glm y la familia gaussiana con enlace identidad da los mismos coeficientes que el ajuste lm del Capítulo 2?

Intuición. La regresión ordinaria es la fila del MLG con una respuesta normal y un enlace identidad, así que los dos ajustes deberían coincidir en cada dígito.

Fórmula. Modelo lineal como MLG: respuesta Normal, $g(\mu) = \mu$ , $\mu = \mathbf{x}'\boldsymbol\beta$ .

Cálculo.

toluca <- read.csv("data/toluca.csv")
lm_fit  <- lm(hours ~ lotsize, data = toluca)
glm_fit <- glm(hours ~ lotsize, family = gaussian(link = "identity"), data = toluca)
rbind(lm = coef(lm_fit), glm = coef(glm_fit))

    (Intercept)  lotsize
lm     62.36586 3.570202
glm    62.36586 3.570202

toluca = pd.read_csv("data/toluca.csv")
lm_fit  = smf.ols("hours ~ lotsize", data=toluca).fit()
glm_fit = smf.glm("hours ~ lotsize", data=toluca,
                  family=sm.families.Gaussian()).fit()
print(lm_fit.params.values)
print(glm_fit.params.values)

[62.36585859  3.57020202]
[62.36585859  3.57020202]

Interpretación. Los dos ajustes coinciden exactamente, $b_0 = 62.37$ y $b_1 = 3.5702$ , los mismos números de Toluca de 2.2 Mínimos cuadrados desde los primeros principios. Los mínimos cuadrados ordinarios no son una técnica separada que vive fuera del mundo del MLG; son el MLG que obtienes al elegir una respuesta normal y un enlace identidad. La regresión logística es la elección binomial-y-logit, la regresión de Poisson es la elección conteo-y-log, y un solo motor de estimación ajusta las tres.

Inténtalo 14.5

Para cada respuesta de abajo, nombra la familia y el enlace del MLG que elegirías, y di qué hace la varianza a medida que crece la media: (a) el número de erratas por página en un manuscrito; (b) si un préstamo cae en incumplimiento o no; (c) el precio de venta de una casa en dólares.

Solución

(a) Erratas por página es un conteo: familia de Poisson, enlace logarítmico, varianza igual a la media (crece linealmente con ella). (b) Incumplimiento o no es binario: familia binomial, enlace logit, varianza $\mu(1-\mu)$ , mayor cerca de $\mu = 0.5$ y encogiéndose hacia los extremos. (c) El precio de venta es continuo (aunque a menudo sesgado, así que una transformación logarítmica o un MLG gamma pueden ayudar): la elección base es la familia gaussiana con enlace identidad y varianza constante $\sigma^2$ . La primera pregunta es siempre qué clase de cosa es la respuesta.

14.7 Resumen del capítulo¶

Ahora puedes modelar conteos. Este capítulo comenzó con por qué los mínimos cuadrados ordinarios son la herramienta equivocada para un conteo, cuyo comportamiento no negativo, de número entero, y de media-igual-a-varianza viola la respuesta no acotada de varianza constante que el modelo lineal supone. La reparación fue el modelo log-lineal de Poisson, que modela el logaritmo de la media como un predictor lineal para que la media ajustada se mantenga positiva y cada coeficiente se lea como una razón de tasas $e^{\beta_j}$ . Maximizar la log-verosimilitud de Poisson dio ecuaciones de puntaje de exactamente la misma forma que las ecuaciones normales del Capítulo 2. Un desplazamiento, el logaritmo de la exposición entrando con coeficiente fijo uno, convirtió el modelo de conteo en un modelo de tasa para los datos de barcos. La dispersión de Pearson $\hat\phi$ diagnosticó la sobredispersión severa en los conteos de las Galápagos, y la corrección cuasi-Poisson infló los errores estándar a tamaños honestos mientras dejaba en paz las estimaciones. Eso finalmente cerró el hilo de las Galápagos de los Capítulos 9 y 10: el problema nunca fue la escala de la respuesta, era que la respuesta es un conteo. Y todo ello, junto con los modelos lineal y logístico, son tres filas de una sola tabla del modelo lineal generalizado.

Resultados clave de un vistazo

Resultado	Enunciado o fórmula	Válido cuando
Distribución de Poisson (Def 14.1)	$P(Y=y) = e^{-\mu}\mu^y / y!$ , $\ E\{Y\}=\operatorname{Var}\{Y\}=\mu$	respuesta de conteo, eventos raros e independientes a una tasa estable
Modelo log-lineal de Poisson (Def 14.2)	$\log \mu_i = \mathbf{x}_i'\boldsymbol\beta$ , así que $\mu_i = e^{\mathbf{x}_i'\boldsymbol\beta}$	la media debe ser positiva y multiplicativa en los predictores
Razón de tasas	$\mu(X_j+1)/\mu(X_j) = e^{\beta_j}$	interpretar un coeficiente log-lineal
Ecuaciones de puntaje de Poisson (Thm 14.3)	$\sum_i (y_i - \hat\mu_i)\mathbf{x}_i = \mathbf{0}$	ajuste por máxima verosimilitud del modelo log-lineal
Desplazamiento (Def 14.4)	$\log\mu_i = \mathbf{x}_i'\boldsymbol\beta + \log t_i$	modelar una tasa con exposición $t_i$
Dispersión de Pearson (Def 14.6)	$\hat\phi = X^2/(n-p)$	revisar el supuesto de varianza-igual-a-media
Error estándar cuasi-Poisson (Def 14.7)	$s_{\text{quasi}} = \sqrt{\hat\phi}\, s_{\text{Poisson}}$	sobredispersión, $\hat\phi > 1$
Modelo lineal generalizado (Def 14.8)	familia $+$ enlace $g(\mu)=\eta$	cualquier respuesta de la familia exponencial

Términos clave. distribución de Poisson, modelo de regresión log-lineal de Poisson, predictor lineal, función de enlace, razón de tasas, ecuaciones de puntaje, desplazamiento, exposición, sobredispersión, dispersión de Pearson, cuasi-Poisson, binomial negativo, modelo lineal generalizado.

Ahora deberías ser capaz de

Explicar por qué los mínimos cuadrados ordinarios son la herramienta equivocada para respuestas de conteo, nombrando los dos supuestos del modelo que violan.
Enunciar el modelo log-lineal de Poisson e interpretar sus coeficientes como razones de tasas, en las unidades del problema.
Derivar la log-verosimilitud de Poisson y sus ecuaciones de puntaje, y conectarlas con la maquinaria de máxima verosimilitud del Capítulo 13.
Usar un desplazamiento para modelar una tasa, y explicar por qué el logaritmo de la exposición entra en el predictor lineal con un coeficiente fijo de uno.
Diagnosticar la sobredispersión a partir de la dispersión de Pearson y aplicar la corrección cuasi-Poisson, y decir cuándo un modelo binomial negativo es el mejor arreglo.
Resolver la pregunta de las especies de las Galápagos con un modelo de conteo, y explicar por qué tiene éxito donde las transformaciones del Capítulo 10 se quedaron cortas.
Ensamblar el marco del MLG en una sola tabla (familia más enlace) que contiene la regresión lineal, la logística y la de Poisson como instancias de una única idea.

Dónde encaja esto. En la columna vertebral del flujo de trabajo de El flujo de trabajo del modelado este capítulo vive sobre todo en AJUSTAR y REVISAR, con una fuerte recompensa en USAR. El modelo log-lineal de Poisson es un paso de AJUSTAR para una respuesta que las herramientas de AJUSTAR de los capítulos anteriores no podían manejar; el diagnóstico de sobredispersión es un paso de REVISAR, el análogo para datos de conteo de las revisiones de varianza de residuos del Capítulo 9 y del modelado de varianza de 10.4 Mínimos cuadrados ponderados; y leer los coeficientes como razones de tasas es el paso de USAR, el primo para datos de conteo de las razones de momios de 13.3 Leer los coeficientes como razones de momios. Mirando atrás, este capítulo completa el arco que comenzó cuando el Capítulo 13 rompió por primera vez con la respuesta normal: la regresión logística manejó el sí o no, la regresión de Poisson maneja el cuántos, y la tabla del MLG muestra ambas como extensiones de los mínimos cuadrados que conoces desde el Capítulo 2. Mirando adelante, el mismo hábito de “elige una familia y un enlace” se lleva a los modelos gamma para datos positivos sesgados, los modelos multinomiales para varias categorías, y los modelos mixtos y jerárquicos que un segundo curso de regresión aborda.

14.8 Preguntas frecuentes¶

P1. ¿No puedo simplemente log-transformar los conteos y correr mínimos cuadrados ordinarios? Puedes, y la gente lo hizo durante décadas, pero tiene costos reales. Los conteos de cero hacen que $\log 0$ sea indefinido, forzando parches incómodos de “sumar 0.5”. Las predicciones retransformadas describen una mediana, no una media, así que no suman correctamente. Y el supuesto de varianza constante sigue siendo erróneo incluso después de tomar logaritmos. La regresión de Poisson esquiva los tres modelando el conteo directamente en sus propios términos.

P2. ¿Por qué el enlace logarítmico específicamente? ¿Podría usar uno distinto? El enlace logarítmico es la elección natural porque hace multiplicativos los coeficientes (razones de tasas), mantiene la media positiva, y es el enlace “canónico” que hace más limpia la matemática de Poisson. Puedes usar otros (el enlace identidad da un modelo de conteo aditivo, ocasionalmente útil), pero el enlace logarítmico es el que viene por defecto por buenas razones y es lo que significa “regresión de Poisson” a menos que se diga otra cosa.

P3. Mi dispersión es $\hat\phi = 0.8$ , por debajo de 1. ¿Es subdispersión, y la corrijo? Los valores modestamente por debajo de 1 son comunes y por lo general inofensivos; hacen los errores estándar de Poisson ligeramente conservadores (demasiado grandes), así que tus pruebas son, en todo caso, demasiado cautelosas. La subdispersión genuina (de, digamos, conteos que son más regulares que aleatorios, como huevos en un cartón lleno) existe pero es rara. Preocúpate por $\hat\phi$ bastante por encima de 1, no modestamente por debajo.

P4. ¿El cuasi-Poisson cambia mis estimaciones de coeficientes? No. Los coeficientes resuelven las ecuaciones de puntaje $\sum_i (y_i - \mu_i)\mathbf{x}_i = \mathbf{0}$ , que nunca mencionan la varianza, así que el cuasi-Poisson los deja idénticos al ajuste de Poisson. Solo cambian los errores estándar, y por lo tanto los valores p y los intervalos de confianza. Ese es todo el punto: las estimaciones estaban bien, la incertidumbre estaba subestimada.

P5. ¿Cuándo debería usar el binomial negativo en vez del cuasi-Poisson? Usa el cuasi-Poisson para una corrección rápida y libre de distribución a los errores estándar. Echa mano del binomial negativo cuando necesites una verosimilitud genuina: para comparar modelos por AIC, para correr pruebas de razón de verosimilitud, o para formar intervalos de predicción para conteos nuevos. El binomial negativo también ajusta mejor la sobredispersión muy fuerte porque su varianza crece cuadráticamente en la media.

P6. ¿Por qué descartar los barcos con cero servicio no sesgó el análisis? Un barco con cero meses de servicio no tuvo oportunidad de sufrir daño, así que no acarrea información sobre la tasa de daño, y $\log 0$ está indefinido de todas formas. Quitar filas con cero exposición es estándar y correcto; no estás descartando datos informativos, estás descartando filas que no tuvieron exposición que modelar.

P7. El Poisson y el cuasi-Poisson dieron la misma devianza residual. ¿No debería la corrección cambiar el ajuste? La devianza mide qué tan lejos quedan las medias ajustadas de los datos, y el cuasi-Poisson no mueve las medias ajustadas, así que la devianza no cambia. Lo que cambia es cómo juzgamos la devianza y los coeficientes contra el azar: el cuasi-Poisson divide entre la dispersión, así que la misma devianza ahora cuenta una historia más cautelosa.

14.9 Problemas de práctica¶

(A) Enumera los tres rasgos de los datos de conteo que el modelo lineal ordinario no respeta, y di qué supuesto del modelo viola cada uno.
(A) Enuncia la media y la varianza de una variable aleatoria de Poisson. ¿Por qué este único hecho hace de los mínimos cuadrados de varianza constante el modelo equivocado para conteos?
(A) En un modelo log-lineal de Poisson un predictor tiene coeficiente $b = 0.5$ . Enuncia su razón de tasas e interprétala en palabras.
(A) Explica en una o dos oraciones qué es un desplazamiento y por qué su coeficiente se fija en uno en vez de estimarse.
(A) Un ajuste de Poisson y su reajuste cuasi-Poisson reportan coeficientes idénticos pero errores estándar distintos. Explica por qué.
(A) Define la sobredispersión y nombra dos causas ordinarias de ella en datos de conteo reales.
(A) Da la familia y la función de enlace para cada fila del MLG: regresión lineal, regresión logística, regresión de Poisson.
(A) Un colega calcula tasas de crimen como crimes / population y ajusta mínimos cuadrados ordinarios. Da dos razones por las que un modelo de Poisson basado en un desplazamiento es mejor.
(A) ¿Por qué un modelo de Poisson nunca puede predecir un conteo negativo, cualesquiera que sean los valores de los predictores? Señala la pieza del modelo que lo garantiza.
(B) Partiendo de la fmp de Poisson $P(Y=y) = e^{-\mu}\mu^y/y!$ , escribe la log-verosimilitud para conteos independientes $y_1, \dots, y_n$ con $\log\mu_i = \mathbf{x}_i'\boldsymbol\beta$ , y deriva las ecuaciones de puntaje $\sum_i (y_i - \mu_i)\mathbf{x}_i = \mathbf{0}$ (Teorema 14.3).
(B) A partir de las ecuaciones de puntaje del problema 10, muestra que la componente del intercepto obliga a $\sum_i \hat\mu_i = \sum_i y_i$ : las medias ajustadas reproducen el conteo total observado.
(B) Muestra que las ecuaciones de puntaje de Poisson tienen la misma forma $\sum_i (y_i - \hat m_i)\mathbf{x}_i = \mathbf{0}$ que las ecuaciones de puntaje logísticas del Capítulo 13 y las ecuaciones normales del Capítulo 2, identificando qué es $\hat m_i$ en cada caso.
(B) Deriva el resultado del desplazamiento: partiendo de $\mu_i = \lambda_i t_i$ y $\log\lambda_i = \mathbf{x}_i'\boldsymbol\beta$ , muestra que $\log\mu_i = \mathbf{x}_i'\boldsymbol\beta + \log t_i$ , y explica por qué $\log t_i$ entra con coeficiente uno.
(B) Muestra que en un modelo log-lineal de Poisson con un solo predictor binario $X \in \{0, 1\}$ , el coeficiente $\beta_1$ es igual a $\log(\mu_1/\mu_0)$ , el logaritmo de la razón de conteos medios entre los dos grupos. Por lo tanto $e^{\beta_1}$ es esa razón.
(B) El modelo cuasi-Poisson supone $\operatorname{Var}\{Y_i\} = \phi\mu_i$ . Explica por qué este reescalado multiplica cada error estándar de coeficiente por $\sqrt\phi$ pero deja las estimaciones de los coeficientes sin cambios, refiriéndote a las ecuaciones de puntaje.
(B) Escribe el estadístico de Pearson $X^2 = \sum_i (y_i - \hat\mu_i)^2/\hat\mu_i$ y explica por qué $E\{X^2\} \approx n - p$ bajo un modelo de Poisson correcto, de modo que $\hat\phi = X^2/(n-p)$ debería estar cerca de 1.
(B) Para dos medias de Poisson $\mu_A$ y $\mu_B$ con exposiciones $t_A, t_B$ entradas como desplazamientos, muestra que la razón de tasas ajustada $\lambda_B/\lambda_A = e^{\beta}$ no depende de las exposiciones, así que la comparación es una comparación justa por exposición.
(B) Explica por qué se usa la máxima verosimilitud, no los mínimos cuadrados, para ajustar el modelo de Poisson, y describe en palabras cómo los mínimos cuadrados reponderados iterativamente alcanzan las estimaciones (puedes citar el esbozo del Capítulo 13).
(C) Ajusta el modelo de Poisson Species ~ Area + Elevation + Nearest + Scruz + Adjacent a gala.csv en R o Python. Reporta el coeficiente de Area y su razón de tasas para un aumento de 100 km $^2$ .
(C) Para el ajuste de Poisson de las Galápagos, calcula la dispersión de Pearson $\hat\phi$ a partir de los residuos de Pearson y confirma que coincide con el valor 31.7 del capítulo. Enuncia qué dice sobre el modelo.
(C) Reajusta el modelo de las Galápagos con family = quasipoisson (R) o la opción scale="X2" (Python). Reporta cuáles predictores son significativos al nivel del 5% bajo cuasi-Poisson pero eran significativos bajo Poisson ingenuo, y explica el cambio.
(C) Ajusta el modelo de daño de los barcos incidents ~ type + year + period con offset = log(service) en las filas de servicio positivo. Reporta la razón de tasas de period75 e interprétala.
(C) Reajusta el modelo de barcos sin el desplazamiento (quita el argumento offset). Compara los coeficientes de type con el modelo con desplazamiento y explica por qué cambian, conectando el cambio con la confusión de tasa y exposición.
(C) Calcula la dispersión de Pearson $\hat\phi$ de los barcos, reajusta como cuasi-Poisson, y reporta cómo cambia el error estándar del coeficiente typeB. ¿Están los datos de barcos tan sobredispersos como los de las Galápagos?
(C) Ajusta el modelo de Toluca hours ~ lotsize de dos maneras: lm y glm con family = gaussian(link = "identity"). Muestra que los coeficientes coinciden, y explica qué dice esto sobre la relación entre los mínimos cuadrados y el MLG.
(C) Usando el ajuste de Poisson de las Galápagos, predice el conteo esperado de especies para una isla con Area = 40, Elevation = 200, Nearest = 1, Scruz = 30, Adjacent = 5. Haz la predicción en la escala del predictor lineal y luego exponencia, y confirma que el software da el mismo número con type = "response".
(C) Ajusta un modelo binomial negativo a los datos de las Galápagos con MASS::glm.nb (R) o NegativeBinomial de statsmodels (Python) y reporta el $\theta$ estimado. Compara su coeficiente de Elevation con los ajustes de Poisson y cuasi-Poisson.
(C) Haz una gráfica diagnóstica de los residuos de Pearson al cuadrado contra las medias ajustadas para el ajuste de Poisson de las Galápagos, con la línea de referencia de varianza igual a la media, y describe qué muestra sobre la sobredispersión.
(A) El Capítulo 10 probó transformaciones logarítmica y de Box-Cox sobre los conteos de especies de las Galápagos, y ayudaron con las gráficas de residuos pero nunca satisficieron del todo. Explica por qué el modelo de Poisson tiene éxito donde esas transformaciones se quedaron cortas. Nombra al menos dos defectos específicos de transformar un conteo que el modelo de conteo evita.

14.10 Práctica de examen¶

PE 14.1. La regresión de Poisson hace dos desviaciones de los mínimos cuadrados ordinarios: modela el logaritmo de la media en vez de la media misma, y estima los coeficientes por máxima verosimilitud en vez de minimizando una suma de cuadrados. Explica por qué cada desviación nos es impuesta por la naturaleza de los datos de conteo. Ata el enlace logarítmico a una propiedad específica de los conteos, y ata la elección de la máxima verosimilitud a otra.

Respuesta modelo

Las dos desviaciones responden a dos fallas distintas de los mínimos cuadrados. El enlace logarítmico es impuesto por el hecho de que una media de conteo debe ser positiva: las especies, los incidentes y las llegadas nunca pueden promediar un número negativo. Si modeláramos la media directamente como $\mathbf{x}'\boldsymbol\beta$ , algunos valores de los predictores empujarían la media ajustada por debajo de cero, un conteo imposible. Modelar $\log\mu_i = \mathbf{x}'\boldsymbol\beta$ en su lugar permite que el predictor lineal recorra toda la recta real mientras la media $\mu_i = e^{\mathbf{x}'\boldsymbol\beta}$ se mantiene estrictamente positiva, y como bono hace multiplicativos los coeficientes, lo que coincide con la forma en que los procesos de conteo realmente escalan.

La máxima verosimilitud es impuesta por una propiedad distinta: la varianza de un conteo no es constante sino igual a su media, así que la respuesta no es ni normal ni de varianza constante. Los mínimos cuadrados son el estimador correcto cuando los errores son normales con varianza constante, porque bajo esas condiciones minimizar la suma de cuadrados es la solución de máxima verosimilitud y es eficiente por Gauss-Markov. Un conteo rompe ambas condiciones, así que los mínimos cuadrados pierden su justificación. La máxima verosimilitud en cambio usa la probabilidad de Poisson real de los conteos observados, ponderando cada caso por cuánta información acarrea, y entrega las ecuaciones de puntaje $\sum_i (y_i - \hat\mu_i)\mathbf{x}_i = \mathbf{0}$ que definen el ajuste.

Una respuesta débil nombra solo el punto de “los conteos no pueden ser negativos” y se detiene, sin conectar la elección de la máxima verosimilitud con la falla de los supuestos de normalidad y varianza constante que hacen óptimos los mínimos cuadrados en primer lugar.

PE 14.2. Un biólogo ajusta un modelo de Poisson de dos predictores a los datos de las Galápagos y te entrega esta salida. Interpreta el coeficiente de Elevation como una razón de tasas para un aumento de 100 metros, di qué significa el intercepto exponenciado y si describe alguna isla real, y enuncia si confiarías en los diminutos errores estándar tal como están.

gala <- read.csv("data/gala.csv")
f <- glm(Species ~ Elevation + Nearest, family = poisson, data = gala)
summary(f)

Coefficients:
             Estimate Std. Error z value Pr(>|z|)
(Intercept) 3.574e+00  3.819e-02  93.603  < 2e-16 ***
Elevation   1.475e-03  3.341e-05  44.151  < 2e-16 ***
Nearest     6.559e-03  1.444e-03   4.541 5.59e-06 ***

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 3510.7  on 29  degrees of freedom
Residual deviance: 1806.4  on 27  degrees of freedom

Respuesta modelo

El coeficiente de elevación es 0.001475 en la escala logarítmica, así que su razón de tasas para un aumento de 100 metros es $e^{100 \times 0.001475} = e^{0.1475} = 1.159$ . En contexto: manteniendo fija la distancia a la isla más cercana, cada 100 metros adicionales de elevación multiplican el número esperado de especies de plantas por cerca de 1.16, un aumento del 16%. Este es un multiplicador, no un conteo aditivo; sería incorrecto decir “16 especies más”.

El intercepto es 3.574 en la escala logarítmica, así que $e^{3.574} = 35.7$ . Literalmente es el conteo esperado de especies para una isla con elevación cero y distancia cero a la isla más cercana. Ninguna isla real tiene elevación cero, así que esto es una extrapolación fuera de los datos y no debería leerse como una predicción para ninguna isla real; es una línea base matemática por la que pasa el ajuste.

No confiaría en los errores estándar tal como están impresos. Los valores $z$ de 44 y 94 son inverosímilmente grandes, y la devianza residual de 1806 en 27 grados de libertad está muy por encima de sus grados de libertad, una señal clara de sobredispersión severa (el ajuste completo de las Galápagos tenía $\hat\phi \approx 32$ ). Los errores estándar impresos suponen varianza igual a la media; bajo sobredispersión son varias veces demasiado pequeños, así que las estrellas de significancia no pueden creerse hasta que la dispersión se estime y se aplique una corrección cuasi-Poisson.

Una respuesta débil lee el efecto de la elevación de forma aditiva o simplemente repite las estrellas de significancia sin notar que los enormes valores $z$ y la devianza residual inflada apuntan a una sobredispersión que hace no confiables esos errores estándar.

PE 14.3. En el modelo de daño de los barcos, ajustado con un desplazamiento de $\log(\text{service})$ , un estudiante ve el coeficiente period75 $= 0.385$ y escribe: “Los barcos operados en el periodo posterior sufren 0.385 incidentes de daño más que los del periodo anterior”. Dos cosas distintas están mal con esa oración. Identifica y corrige ambas, luego escribe la afirmación que el estudiante debió haber hecho.

Respuesta modelo

El primer error es la escala. El coeficiente 0.385 vive en la escala del logaritmo de la media, donde los efectos son aditivos, pero en la escala de conteo los efectos son multiplicativos. Para interpretarlo exponencias: $e^{0.385} = 1.469$ . Así que el periodo posterior no suma 0.385 incidentes; multiplica el conteo esperado por cerca de 1.47, una razón de tasas, que es un aumento del 47%, no un número fijo de incidentes extra.

El segundo error es el desplazamiento. Como el modelo lleva $\log(\text{service})$ como un desplazamiento, sus coeficientes describen una tasa de daño por mes de servicio, no un conteo crudo de incidentes. Así que la comparación es entre tasas de daño por unidad de exposición, no entre totales; un barco ajetreado y un barco ocioso en el mismo periodo tienen la misma tasa pero conteos crudos muy distintos. Ignorar el desplazamiento convierte una afirmación sobre tasas en una afirmación falsa sobre conteos.

La oración correcta: “Manteniendo fijos el tipo de barco y el año de construcción, los barcos que operan en el periodo posterior tuvieron cerca de 47% más incidentes de daño por mes de servicio que los barcos del periodo anterior”.

Una respuesta débil arregla solo la exponenciación y reporta “47% más incidentes”, olvidando que el desplazamiento hace de la cantidad una tasa por mes de servicio en vez de un conteo crudo de incidentes.

PE 14.4. Los dos ajustes de abajo son el mismo modelo de barcos, una vez con el desplazamiento de exposición y una vez sin él. Los coeficientes de type cambian de forma pronunciada. Explica qué le pasa al coeficiente del tipo B cuando se quita el desplazamiento y por qué, usando las cifras de exposición, y enuncia cuál de los dos ajustes responde a la pregunta “¿qué tipo de barco es el más peligroso de operar?”.

ships <- read.csv("data/ships.csv"); ships <- ships[ships$service > 0, ]
with_offset <- glm(incidents ~ type + year + period, family = poisson,
                   offset = log(service), data = ships)
no_offset   <- glm(incidents ~ type + year + period, family = poisson, data = ships)

type coefficients WITH offset:    typeB -0.5433  typeC -0.6874  typeD -0.0760  typeE  0.3256
type coefficients NO offset:      typeB  1.7957  typeC -1.2528  typeD -0.9045  typeE -0.1463

mean months of service by type:   A 1356   B 19760   C 885   D 635   E 855
total incidents by type:          A   42   B   253   C  12   D  17   E  32

Respuesta modelo

Con el desplazamiento, typeB $= -0.5433$ , así que su razón de tasas es $e^{-0.5433} = 0.581$ : por mes de servicio, los barcos tipo B sufren daño a cerca del 58% de la tasa base del tipo A, haciéndolos más seguros que el tipo A, no más peligrosos. Quita el desplazamiento y el coeficiente se voltea a +1.7957, una razón de $e^{1.7957} = 6.0$ , que hace que el tipo B parezca seis veces peor que el tipo A.

Las cifras de exposición explican el vuelco. Los barcos tipo B acumulan una media de $19{,}760$ meses de servicio, cerca de quince veces la media del tipo A de $1{,}356$ , y amasan 253 incidentes totales, muchos más que cualquier otro tipo. Sin el desplazamiento el modelo está explicando conteos crudos, así que le acredita al tipo B un coeficiente enorme simplemente porque esos barcos están en el mar muchísimo más tiempo y por lo tanto tienen más oportunidades de ser dañados. El desplazamiento $\log(\text{service})$ divide fuera esa exposición, convirtiendo la respuesta en una tasa por mes; una vez que se toma en cuenta la exposición, la tasa de daño del tipo B está en realidad por debajo de la base del tipo A. El volumen de incidentes y la tasa de incidentes son preguntas distintas, y el desplazamiento es lo que las separa.

El ajuste con desplazamiento es el que responde a “qué tipo es el más peligroso de operar”, porque el peligro es una tasa por unidad de exposición, no un total crudo. Por ese ajuste el tipo E, con el único coeficiente positivo, tiene la tasa de daño más alta; el ajuste sin desplazamiento solo te dice qué tipo resultó registrar más incidentes, lo que es sobre todo una afirmación sobre el tamaño de la flota y el tiempo en el mar.

Una respuesta débil nota el cambio de signo pero no lo conecta con el tiempo de servicio muy grande del tipo B, sin ver que el modelo sin desplazamiento está confundiendo la exposición con la tasa.

PE 14.5. Para el ajuste completo de Poisson de las Galápagos con cinco predictores la dispersión de Pearson es $\hat\phi = 31.7$ . (a) Explica en contexto qué dice este único número sobre los conteos de especies y sobre los errores estándar de Poisson. (b) Un colega reajusta los mismos datos con un modelo binomial negativo y reporta $\hat\theta = 1.675$ y un coeficiente de Elevation de 0.003855, contra el valor de Poisson 0.003541. Interpreta $\hat\theta$ , y explica por qué los dos coeficientes de elevación son casi idénticos aunque sus errores estándar no lo son.

Respuesta modelo

(a) Un modelo de Poisson promete que la varianza de un conteo es igual a su media, lo que haría $\hat\phi \approx 1$ . El valor $\hat\phi = 31.7$ dice que los conteos de especies de las Galápagos se dispersan cerca de 32 veces más ampliamente de lo que Poisson permite, sobredispersión severa, muy probablemente por rasgos de las islas que los cinco predictores no captan, agrupamiento espacial de especies, y diferencias genuinas de isla a isla en la tasa de colonización subyacente. La consecuencia es específica: la sobredispersión apenas mueve las estimaciones de los coeficientes, pero hace los errores estándar de Poisson demasiado pequeños por un factor de cerca de $\sqrt{31.7} \approx 5.6$ . Así que los valores p de Poisson ingenuo son demasiado optimistas, y los predictores que parecen fuertemente significativos pueden no sobrevivir a una corrección honesta; las estimaciones puntuales, sin embargo, siguen siendo confiables.

(b) El parámetro $\hat\theta = 1.675$ es el parámetro de dispersión extra del binomial negativo, que fija la varianza a $\mu + \mu^2/\theta$ . Un $\theta$ pequeño hace grande ese término cuadrático, así que $\hat\theta = 1.675$ señala sobredispersión fuerte y confirma la historia que contó $\hat\phi$ . Los dos coeficientes de elevación, 0.003855 y 0.003541, coinciden a dos dígitos porque las correcciones de sobredispersión cambian cómo se modela la varianza, no cómo se modela la media. Ambos ajustes describen esencialmente la misma relación log-lineal de la media y resuelven casi las mismas ecuaciones de ajuste de la media; lo que el binomial negativo cambia es el ancho de los errores estándar, ampliándolos respecto al Poisson ingenuo para reflejar la dispersión verdadera. La media es estable; solo se mueve la incertidumbre declarada.

Una respuesta débil trata $\hat\phi = 31.7$ como evidencia de que todo el modelo está mal y debe descartarse, o lee un $\hat\theta$ grande como más dispersión, sin ver que la sobredispersión infla los errores estándar mientras deja las estimaciones de los coeficientes esencialmente sin cambios.

Juego del capítulo¶

Chapter summary (in English)

This chapter covers Poisson regression for count responses, using two datasets: plant species on 30 Galapagos islands and damage incidents on cargo ships. A count is never negative, is a whole number, and its variability grows with its average; the ordinary linear model assumes the opposite, so it fails. The Poisson distribution describes counts with a single quantity, its mean $\mu$ , and a key property: the variance equals the mean.

The Poisson log-linear model proposes $\log\mu_i = \mathbf{x}_i'\boldsymbol\beta$ : the log of the mean is linear in the predictors, so the mean is never negative and the coefficients are read multiplicatively: $e^{\beta_j}$ is a rate ratio, the count-data cousin of the odds ratio of logistic regression. We estimate by maximum likelihood; deriving the log-likelihood gives the score equations $\sum_i (y_i - \mu_i)\mathbf{x}_i = \mathbf{0}$ , the same form as the normal equations of Chapter 2.

To model a rate instead of a raw count, we add the log of the exposure as an offset, with fixed coefficient one. For the ships, the exposure is months of service; type B, which looked worst by raw count, turns out average once we adjust for exposure.

Overdispersion happens when counts vary more than Poisson allows. It is diagnosed with the Pearson dispersion $\hat\phi = X^2/(n-p)$ , which should sit near 1; for the Galapagos it is 31.7. The quasi-Poisson correction multiplies each standard error by $\sqrt{\hat\phi}$ without changing the coefficients; the negative binomial is the full-likelihood alternative. This finally resolves the Galapagos thread from Chapters 9 and 10: the problem was never the scale, it was that the response is a count. Finally, the generalized linear model unifies everything: a distribution plus a link function $g(\mu)=\eta$ produce linear, logistic, and Poisson regression as three rows of one table, with a single estimation engine.