10. Medidas correctivas y transformaciones (en espanol)

10. Medidas correctivas y transformaciones¶

Un biólogo alinea los pesos del cerebro y del cuerpo de 62 mamíferos terrestres, desde un ratón de 3 gramos hasta un elefante africano de cinco toneladas, y grafica uno contra el otro. La imagen es casi inútil. Sesenta de los puntos quedan aplastados en la esquina inferior izquierda, el elefante y unas pocas ballenas flotan solos en la esquina superior derecha, y cualquier recta que dibujes es arrastrada por los gigantes. La Figure 1 es ese gráfico. La dispersión de los puntos también crece con el tamaño del cuerpo, así que el supuesto ordenado de varianza constante del Capítulo 2 no aparece por ningún lado.

Nada está mal con los datos. Lo que está mal es la escala. Toma el logaritmo de ambos ejes y los mismos 62 puntos se acomodan en una banda limpia, recta y dispersa de forma pareja, como se ve en la Figure 2. La relación siempre fue simple; era simple en términos porcentuales, no en gramos y kilogramos crudos, y la escala logarítmica es donde viven los porcentajes. Este es el movimiento central del capítulo. Cuando los diagnósticos del Capítulo 9 se activan, rara vez necesitas abandonar la regresión. Más a menudo necesitas ajustarla a una versión transformada del problema, o ponderar las observaciones, o reemplazar la función de pérdida, para que los supuestos del modelo se vuelvan ciertos.

El Capítulo 9 te enseñó a detectar problemas: gráficos de residuos curvados, formas de embudo, colas pesadas, una isla o un país que doblan todo el ajuste. Este capítulo es la caja de herramientas de las soluciones. Leerás coeficientes logarítmicos como porcentajes, elegirás una transformación con la verosimilitud de Box-Cox, ponderarás observaciones que llevan ruido desigual, y resistirás un valor atípico terco con la regresión robusta. Al final conocerás un conjunto de datos que ninguna transformación puede salvar, y aprenderás a reconocer cuándo la solución honesta es un modelo completamente distinto.

10.1 Un menú de soluciones¶

Antes de cualquier técnica individual, ayuda ver el menú completo y la lógica que elige entre sus elementos. El Capítulo 9 terminó nombrando problemas, desde gráficos de residuos curvados y dispersión en forma de embudo hasta los casos de alto apalancamiento de 9.3 Influencia: qué puntos cambian realmente el ajuste. Cada problema nombrado tiene una familia de soluciones que le corresponde.

Si el gráfico de residuos se curva, la función de la media está mal: se le está pidiendo a una recta que trace una curva. La solución es doblar el modelo de vuelta a la recta, ya sea transformando el predictor $X$ (dejando el error en paz) o transformando la respuesta $Y$. Si el gráfico de residuos se abre en abanico, la varianza no es constante: el modelo supone un solo $\sigma^2$ pero los datos llevan más ruido en algunos lugares que en otros. La solución es una transformación estabilizadora de la varianza de $Y$ o los mínimos cuadrados ponderados, que le dicen al ajuste que confíe más en los puntos de bajo ruido. Si los residuos tienen colas pesadas o un valor atípico solitario, uno o dos puntos están dominando la pérdida de error cuadrático. La solución es una pérdida que crece más despacio que el cuadrado, que es la regresión robusta. Y si la respuesta es un conteo, una proporción, o una cantidad estrictamente positiva con sesgo estructural, el modelo de errores normales puede ser el marco equivocado por completo, y la solución es un modelo distinto, que es a donde van los Capítulos 13 y 14.

Toda la lógica cabe en una página. La Figure 3 la presenta como un diagrama de decisión: empieza en el gráfico de residuos, lee el patrón, y sigue la flecha hasta la solución. Manténla a tu lado durante el resto del capítulo, que recorre las cuatro columnas de izquierda a derecha.

Las transformaciones del predictor y de la respuesta hacen trabajos distintos, y vale la pena mantenerlas separadas. Transformar $X$ cambia la forma de la curva de la media sin tocar los errores, así que es la herramienta para una relación curvada cuya dispersión ya es pareja. Transformar $Y$ cambia la forma y la dispersión a la vez, así que es la herramienta cuando una curva y un embudo aparecen juntos, exactamente como en los datos de los mamíferos. Una guía aproximada, dibujada en la Figure 4, es la “escalera de potencias”: si un gráfico se dobla en un sentido, baja por la escalera desde $Y$ hacia $\sqrt{Y}$, $\log Y$, $1/Y$; si se dobla en el otro sentido, sube hacia $Y^2$. El método de Box-Cox en 10.3 Box-Cox: dejar que los datos elijan la potencia convierte esta adivinanza en una estimación.

El resto del capítulo toma estas soluciones una a la vez. En todo momento, mantén a la vista el modelo del Capítulo 2: $Y_i = \beta_0 + \beta_1 X_i + \varepsilon_i$ con errores que deben promediar cero, compartir una varianza $\sigma^2$, y estar no correlacionados. Cada solución aquí es un esfuerzo por hacer verdadera una de esas tres condiciones después de que un diagnóstico mostró que era falsa.

10.2 La transformación logarítmica y la lectura de sus coeficientes¶

La transformación logarítmica se gana su propia sección porque es la solución más común y la que más a menudo se malinterpreta. Ahora podemos ajustar la recta de los mamíferos; la pregunta más difícil es qué significa en realidad su pendiente de 0.75, porque un coeficiente sobre una variable en logaritmos no es un cambio en gramos por kilogramo. Es un porcentaje.

Intuición¶

Los logaritmos convierten la multiplicación en suma, así que convierten “crece en un porcentaje” en “crece en una constante”. Una cantidad que se duplica cada vez que algún motor se duplica se ve curvada en ejes ordinarios y recta en ejes logarítmicos. El peso del cerebro se relaciona con el peso del cuerpo de esa manera: entre especies, un cuerpo diez veces más pesado tiende a llevar un cerebro un múltiplo fijo más grande, no un número fijo de gramos más grande. Por eso el gráfico crudo se curvó y el gráfico log-log es recto. El precio de la recta es que la pendiente ahora habla en porcentajes, y tienes que traducir.

Fórmula¶

Dos modelos logarítmicos cubren casi todos los casos. El modelo log-log toma el logaritmo de ambos lados,

y el modelo log-lineal toma el logaritmo solo de la respuesta,

Aquí $\log$ es el logaritmo natural (base $e$) en todo este libro. La pendiente significa algo distinto en cada modelo, y cada significado lleva un nombre (Definición 10.1, Definición 10.2).

En palabras: en la escala log-log la pendiente responde “si $X$ sube 1%, ¿cuántos por ciento sube $Y$ aproximadamente”, y en la escala log-lineal responde “si $X$ sube en una de sus propias unidades, ¿cuántos por ciento sube $Y$ aproximadamente”. Ambas lecturas tienen una forma exacta y una forma aproximada, y la siguiente derivación produce las dos.

Derivación (lecturas porcentuales exacta y aproximada)¶

Demostración. Deshaz el logaritmo sobre la respuesta. El modelo log-log dice $\log Y = \beta_0 + \beta_1 \log X$ para la media, así que exponenciando,

En palabras: la media de $Y$ es una constante por una potencia de $X$. Ahora multiplica $X$ por un factor $c$ (para un aumento del 1%, $c = 1.01$; para una duplicación, $c = 2$). La nueva media dividida entre la vieja es

Así que multiplicar $X$ por $c$ multiplica $Y$ por $c^{\beta_1}$, de forma exacta. El cambio porcentual exacto en $Y$ para una subida del 1% en $X$ es $100\,(1.01^{\beta_1} - 1)$. Para la lectura aproximada, usa $c^{\beta_1} = e^{\beta_1 \log c} \approx 1 + \beta_1 \log c$ cuando $\log c$ es pequeño; con $c = 1.01$, $\log c \approx 0.01$, lo que da un cambio porcentual de cerca de $100 \cdot \beta_1 \cdot 0.01 = \beta_1$ por ciento. Ese es el titular de la elasticidad: en la escala log-log, $\beta_1$ es aproximadamente el cambio porcentual en $Y$ por cambio porcentual en $X$. $\blacksquare$

Demostración (semielasticidad log-lineal). En el modelo log-lineal la media satisface $\log Y = \beta_0 + \beta_1 X$, así que $Y = e^{\beta_0} e^{\beta_1 X}$. Aumenta $X$ en una unidad:

Así que cada paso de una unidad en $X$ multiplica $Y$ por $e^{\beta_1}$, de forma exacta, lo que da un cambio porcentual exacto de $100\,(e^{\beta_1} - 1)$. Como $e^{\beta_1} \approx 1 + \beta_1$ para $\beta_1$ pequeño, el cambio porcentual aproximado es $100\,\beta_1$ por ciento por unidad. $\blacksquare$

La brecha entre lo exacto y lo aproximado se ensancha a medida que el coeficiente crece, dibujada en la Figure 5. Para un coeficiente cercano a $\pm 0.1$ los dos coinciden hasta un error de redondeo, que es la razón por la cual el atajo del “porcentaje por unidad” es seguro para efectos pequeños y engañoso para los grandes. Esta misma maquinaria de $e^{\beta}$ regresa dos veces más adelante: como la razón de momios en la regresión logística (13.3 Leer los coeficientes como razones de momios) y como el crecimiento porcentual por mes en la serie de aerolíneas en logaritmos (15.5 Pronosticar con honestidad: probar en el futuro).

R¶

El ajuste de los mamíferos es un modelo log-log, así que su pendiente es una elasticidad.

mammals <- read.csv("data/mammals.csv")
dim(mammals)
head(mammals, 3)

[1] 62  3
          species  body brain
1      Arctic fox 3.385  44.5
2      Owl monkey 0.480  15.5
3 Mountain beaver 1.350   8.1

La forma más rápida de sentir por qué una elasticidad es un porcentaje es fijar una tú mismo y leer la traducción mientras cambia.

Qué observar. La lectura exacta $100(1.10^{b_1} - 1)$ y el atajo $10\,b_1$ nunca se separan más de una décima de punto en todo el deslizador, que es justo lo que la derivación promete para un coeficiente pequeño. Prueba esto. Baja $b_1$ desde 1 hasta que SCE deje de bajar, llega cerca de 0.75, y lee el porcentaje que 10.2 La transformación logarítmica y la lectura de sus coeficientes te pide reportar.

10.3 Box-Cox: dejar que los datos elijan la potencia¶

Adivinar una transformación a partir de la escalera de potencias funciona, pero se siente arbitrario, y dos analistas pueden estar en desacuerdo. Box y Cox propusieron una solución con principios: pon las potencias candidatas en una sola familia indexada por un número $\lambda$, y luego deja que la máxima verosimilitud elija el $\lambda$ que hace que el modelo normal ajuste mejor. Ahora podemos motivar el método con un conjunto de datos que se curva y se abre a la vez: los clásicos datos cars, donde la distancia de frenado de un auto se mide contra su velocidad.

Intuición¶

Las potencias $\sqrt{Y}$, $\log Y$, $Y$, y $Y^2$ son todas casos especiales de elevar $Y$ a una potencia. Si escribes la transformación como una sola función suave de un parámetro de potencia $\lambda$, entonces elegir una transformación se vuelve elegir un número, y elegir un número es algo que la verosimilitud hace bien. Box-Cox escribe la verosimilitud de los datos como función de $\lambda$, con los coeficientes de regresión y la varianza del error perfilados, y reporta el $\lambda$ que la maximiza, más un intervalo de confianza para que sepas con qué precisión los datos lo fijan.

Fórmula¶

• $\lambda$ es el parámetro de potencia por estimar; $\lambda = 1$ es esencialmente ninguna transformación, $\lambda = 0.5$ es la raíz cuadrada, $\lambda = 0$ es el logaritmo, $\lambda = -1$ es el recíproco.

• El “-1” y la división entre $\lambda$ hacen que la familia sea continua en $\lambda = 0$, donde se vuelve $\log Y$, así que la escalera no tiene huecos.

En palabras: un solo dial $\lambda$ se desliza de forma suave por todos los peldaños de la escalera de potencias, y el logaritmo se sienta exactamente en la muesca de $\lambda = 0$.

Derivación (Box-Cox como verosimilitud perfilada)¶

Demostración. Supón que para la potencia correcta la respuesta transformada sigue el modelo lineal normal,

donde $\mathbf{x}_i'$ es la $i$-ésima fila de la matriz de diseño $\mathbf{X}$ de 7.1 El modelo y los mínimos cuadrados en forma matricial. En palabras: una vez que $Y$ se eleva a la potencia correcta, el resultado obedece el modelo de regresión normal ordinario. Los parámetros son $\lambda$, $\boldsymbol{\beta}$, y $\sigma^2$. El punto sutil es que la verosimilitud debe escribirse para los datos originales $Y_i$, no para los transformados $Y_i^{(\lambda)}$, porque estamos comparando distintos valores de $\lambda$ y cada uno pone $Y^{(\lambda)}$ en una escala distinta. Cambiar variables de $Y_i^{(\lambda)}$ a $Y_i$ trae el jacobiano, el factor de la derivada que introduce cualquier cambio de variables, aquí $\mathrm{d}Y_i^{(\lambda)}/\mathrm{d}Y_i = Y_i^{\lambda - 1}$. Por lo tanto, la log-verosimilitud de los $Y_i$ observados es

En palabras: esta es la log-verosimilitud normal usual para la respuesta transformada, más un término extra que corrige por el estiramiento de la escala de $Y$ (la suma de los log-jacobianos). Ahora perfila $\boldsymbol{\beta}$ y $\sigma^2$: para cualquier $\lambda$ fijo, los valores que maximizan $\ell$ son exactamente el ajuste por mínimos cuadrados de $Y^{(\lambda)}$ sobre $\mathbf{X}$, lo que da una suma de cuadrados de los residuos $\mathrm{RSS}(\lambda)$ y $\hat\sigma^2(\lambda) = \mathrm{RSS}(\lambda)/n$ (el mismo argumento que la estimación de máxima verosimilitud normal de la varianza del Capítulo 2). Sustituyendo estos de vuelta y descartando constantes queda la log-verosimilitud perfilada, una función de $\lambda$ sola:

En palabras: para cada $\lambda$, ajusta la respuesta transformada por mínimos cuadrados ordinarios, lee su suma de cuadrados de los residuos, penaliza por el término jacobiano, y el mejor $\lambda$ es el que maximiza el resultado. El maximizador $\hat\lambda$ es la estimación de Box-Cox. Un intervalo de confianza del $100(1-\alpha)$% es el conjunto de $\lambda$ con $\ell_p(\lambda) \ge \ell_p(\hat\lambda) - \tfrac{1}{2}\chi^2_{1,\,1-\alpha}$, el intervalo de razón de verosimilitud usual. $\blacksquare$

R¶

MASS::boxcox de R calcula $\ell_p(\lambda)$ sobre una malla y lo dibuja. El gráfico en la Figure 8 es la log-verosimilitud perfilada para el modelo de cars, con su pico y su intervalo del 95% marcados.

cars <- read.csv("data/cars.csv")
fit_lin <- lm(dist ~ speed, data = cars)
coef(fit_lin)

(Intercept)       speed
 -17.579095    3.932409

El ajuste de la recta no está obviamente mal, pero sus residuos se abren en abanico a medida que la distancia ajustada crece, el embudo que se ve en la Figure 7. Esa dispersión que se ensancha es la señal de que la respuesta necesita una transformación, y Box-Cox dirá cuál.

La verosimilitud perfil se vuelve más creíble cuando la has visto pelear un rato con el gráfico de residuos, así que toma tú el dial.

Qué observar. La log-verosimilitud perfil y la forma del gráfico de residuos coinciden: $\ell_p$ sube de -135.63 en $\lambda = 1$ hasta su máximo -126.73 en $\lambda = 0.43$ justo cuando el embudo se aplana, y $\mathrm{RSS}(\lambda)$ por sí solo no te habría dicho nada, pues cae hasta 9.56 en el logaritmo. Prueba esto. Encuentra el máximo a mano, y luego sigue bajando más allá de $\lambda = 0.2$ para ver cómo se abre un nuevo valor atípico por abajo, la imagen de transformar de más contra la que advierte 10.3 Box-Cox: dejar que los datos elijan la potencia.

10.4 Mínimos cuadrados ponderados¶

Las transformaciones atacan una media curvada o una varianza que crece con la media. A veces, sin embargo, la varianza es desigual por una razón que no tiene nada que ver con la media: algunas observaciones simplemente se miden con más precisión que otras. Promediar diez lecturas da un punto más confiable que una sola lectura; un instrumento bien calibrado le gana a uno tosco. Cuando conoces la precisión relativa de cada punto, no deberías tirar ese conocimiento tratando cada punto por igual. Los mínimos cuadrados ponderados son cómo lo usas.

Intuición¶

Los mínimos cuadrados ordinarios le dan a cada observación el mismo voto. Si un punto lleva diez veces la varianza de otro, eso es como dejar que un testigo ruidoso declare tan fuerte como uno cuidadoso. Los mínimos cuadrados ponderados bajan el volumen a los puntos ruidosos y lo suben a los precisos, en proporción exactamente inversa a sus varianzas. La Figure 11 muestra la idea: los puntos de bajo ruido se dibujan grandes porque cuentan más.

Fórmula¶

Mantén el modelo lineal pero deja que las varianzas del error difieran:

Aquí $w_i$ se conoce salvo la constante común $\sigma^2$.

En palabras: penaliza cada residuo cuadrado por su peso, para que los puntos precisos contribuyan más al total que el ajuste intenta reducir.

Derivación (WLS a partir de Gauss-Markov)¶

Demostración. El truco es transformar el modelo de varianza desigual en uno de varianza igual y luego citar el resultado de Gauss-Markov que ya demostramos. Multiplica la $i$-ésima observación por $\sqrt{w_i}$. Reuniendo estas en matrices con $\mathbf{W}^{1/2} = \operatorname{diag}(\sqrt{w_1},\dots,\sqrt{w_n})$, define

El modelo transformado es $\mathbf{Y}^{\ast} = \mathbf{X}^{\ast}\boldsymbol{\beta} + \boldsymbol{\varepsilon}^{\ast}$, y sus errores son ahora homocedásticos:

El modelo con estrella satisface cada supuesto de los mínimos cuadrados ordinarios, así que por el teorema matricial de Gauss-Markov (7.6 El teorema de Gauss-Markov) el mejor estimador lineal insesgado es mínimos cuadrados ordinarios sobre los datos con estrella:

Este es el estimador WLS, y como es OLS sobre un modelo que verdaderamente tiene varianza constante, hereda cada buena propiedad: es insesgado y es BLUE, con $\operatorname{Var}\{\mathbf{b}_W\} = \sigma^2 (\mathbf{X}'\mathbf{W}\mathbf{X})^{-1}$. El OLS simple sobre los datos originales sigue siendo insesgado aquí, pero ya no es el mejor: desperdicia la información de precisión que hay en los pesos. $\blacksquare$

La derivación también te dice cómo calcular WLS con nada más que una rutina de OLS: forma $\sqrt{w_i}\,Y_i$ y $\sqrt{w_i}\,\mathbf{x}_i$ y haz la regresión. Lo usamos como verificación más abajo.

R¶

Los datos strongx son un experimento de física: en diez energías de haz, se midió una sección eficaz de dispersión, y cada medición viene con su propia desviación estándar conocida sd. Las desviaciones estándar conocidas son el caso de libro de texto para WLS, con pesos $w_i = 1/\text{sd}_i^2$.

strongx <- read.csv("data/strongx.csv")
strongx[, c("energy", "crossx", "sd")]

   energy crossx sd
1   0.345    367 17
2   0.287    311  9
3   0.251    295  9
4   0.225    268  7
5   0.207    253  7
6   0.186    239  6
7   0.161    220  6
8   0.132    213  6
9   0.084    193  5
10  0.060    192  5

Cuando las varianzas no se conocen, estimas los pesos, por lo general modelando cómo crece la dispersión (por ejemplo, haciendo la regresión de los residuos absolutos sobre un predictor y usando los valores ajustados para construir los pesos). Ese método de dos pasos, iterado, es común y efectivo, y recicla el instinto de modelado de la varianza que el Capítulo 14 reutiliza para manejar la sobredispersión en los modelos de conteo. El caso de strongx es más limpio porque la física nos entregó las varianzas directamente.

Los pesos son más fáciles de creer cuando puedes subirlos poco a poco y ver cómo responde la recta, así que desliza tú mismo entre mínimos cuadrados ordinarios y ponderados.

Qué observar. En $p = 0$ la recta ponderada queda exactamente sobre la ordinaria, que es la afirmación de la derivación de que pesos iguales devuelven los mínimos cuadrados de siempre. Prueba esto. Sube $p$ hasta 1 y observa cómo la pendiente cae de 619.7 a 530.8 mientras el error estándar de $b_0$ baja de 10.08 a 8.08: los números del Ejemplo 10.3, llegando como una rotación que puedes ver. Sigue hasta $p = 2$ y pregúntate si aún confiarías en un ajuste apoyado en tan pocos puntos (10.4 Mínimos cuadrados ponderados).

10.5 Cuando una transformación no basta: los conteos de Galápagos¶

Cada solución hasta ahora supuso que el modelo subyacente era sólido y que solo su escala o su ponderación estaba mal. A veces ese supuesto es falso, y una transformación solo esconde el problema. Los datos de especies de Galápagos del Capítulo 9 son un caso ilustrativo, y seguirlos aquí te muestra cómo distinguir un problema de escala reparable de un modelo que necesita reemplazo.

Recuerda el planteamiento de 9.3 Influencia: qué puntos cambian realmente el ajuste: para 30 islas de Galápagos, el número de especies de plantas se hace la regresión sobre predictores geográficos (área, elevación y distancias), y una isla, Isabela, tiene un valor de apalancamiento enorme porque empequeñece a las demás en área. El ajuste crudo tiene un embudo de residuos, porque las islas con más especies también se dispersan más ampliamente, la firma de los datos de conteo. Una solución natural es una transformación estabilizadora de la varianza. Para conteos, la raíz cuadrada es la elección clásica, ya que un conteo de Poisson con media $\mu$ tiene varianza $\mu$, y $\sqrt{\text{count}}$ tiene una varianza aproximadamente constante.

gala <- read.csv("data/gala.csv")
fit_g_lin  <- lm(Species ~ Area + Elevation + Nearest + Scruz + Adjacent,
                 data = gala)
fit_g_sqrt <- lm(sqrt(Species) ~ Area + Elevation + Nearest + Scruz + Adjacent,
                 data = gala)
round(c(R2_linear = summary(fit_g_lin)$r.squared,
        R2_sqrt   = summary(fit_g_sqrt)$r.squared), 4)

R2_linear   R2_sqrt
   0.7658    0.7827

La raíz cuadrada ayuda un poco, y la Figure 14 muestra que el embudo se alivia. Pero mira más de cerca y las grietas aparecen. Isabela todavía tiene un valor de apalancamiento de 0.97 en el ajuste transformado, esencialmente sin cambio, porque transformar $Y$ no hace nada respecto a una $X$ que está lejos de las demás.

h <- hatvalues(fit_g_sqrt)
names(h) <- gala$island
round(sort(h, decreasing = TRUE)[1:3], 3)

   Isabela Fernandina     Darwin
     0.969      0.950      0.466

La Figure 15 muestra por qué se abre la brecha, usando la curva de elevar al cuadrado que deshace un ajuste de raíz cuadrada. Toma dos predicciones a distancias iguales por encima y por debajo del promedio. Elevar al cuadrado dobla la curva hacia arriba, así que la de arriba sube más de lo que baja la de abajo. Su promedio (sobre la cuerda recta) cae por encima del cuadrado del promedio (sobre la curva), y esa brecha vertical es el sesgo. Cualquier retrotransformación curvada, incluido el logaritmo, hace lo mismo.

Da un paso atrás y cuenta las señales de advertencia. La respuesta es un conteo. Varias islas tienen muy pocas especies (seis islas tienen menos de diez). La varianza crece con la media por construcción. Las predicciones retrotransformadas están sesgadas, y un ajuste lineal simple a los conteos predice un número negativo de especies para las islas más pequeñas, lo cual es imposible. Cada una de estas dice lo mismo: el modelo lineal normal de varianza constante es el marco equivocado, y ninguna potencia de $Y$ lo hará correcto. La solución honesta no es transformar los datos hasta que encajen en el modelo, sino cambiar el modelo por uno construido para conteos. Ese modelo es la regresión de Poisson, y 14.1 Por qué los conteos rompen el modelo lineal regresa a estas exactas islas y las ajusta de manera apropiada, recordando tanto el estatus de alto apalancamiento de Isabela del Capítulo 9 como esta transformación fallida. Por ahora, la lección es de diagnóstico: aprende a distinguir un problema de escala, que una transformación arregla, de un problema de distribución, que no puede.

10.6 Una primera mirada a la regresión robusta¶

La última solución aborda una falla distinta: no una media curvada o una varianza desigual, sino un puñado de puntos con residuos desmesurados que arrastran el ajuste hacia sí mismos. Los mínimos cuadrados son exquisitamente sensibles a tales puntos porque elevan los residuos al cuadrado, así que uno solo grande puede dominar toda la suma. La regresión robusta reemplaza el cuadrado por una pérdida que crece más suavemente, para que ningún punto pueda apoderarse de todo.

La Figure 16 muestra la recompensa sobre una nube simple con un punto lanzado muy lejos. La recta de mínimos cuadrados se inclina hacia arriba para perseguir el punto extraviado; la recta de Huber lo ignora y se queda con la multitud. Esa es toda la promesa de la regresión robusta en una sola imagen, y el resto de la sección explica cómo la pérdida con tope la cumple.

Intuición¶

Imagina la suma de residuos cuadrados como un juego de la cuerda en el que cada punto jala con una fuerza proporcional a su residuo. Un punto que está el doble de lejos jala cuatro veces más fuerte, por el cuadrado. La idea de Huber es ponerle tope a esa escalada: dentro de un rango normal, mantén la pérdida cuadrada y su eficiencia, pero una vez que un residuo es lo bastante grande como para parecer un valor atípico, deja que su jalón crezca solo linealmente, no cuadráticamente. La Figure 17 dibuja el peso resultante que recibe cada punto.

Fórmula¶

La estimación M, donde la M señala un objetivo de estilo de máxima verosimilitud, reemplaza el objetivo de los mínimos cuadrados por

donde $\rho_c$ es la pérdida de Huber.

• $u = (Y_i - \mathbf{x}_i'\boldsymbol{\beta})/s$ es el residuo estandarizado por una estimación de escala resistente $s$.

• $c$ es una constante de ajuste; $c = 1.345$ da cerca del 95% de la eficiencia de OLS cuando los errores realmente son normales, mientras acota la influencia de los valores atípicos.

En palabras: penaliza los residuos ordinarios con el cuadrado familiar, pero penaliza los residuos lejanos solo linealmente, para que un valor atípico sea caro en lugar de catastrófico.

Derivación (IRLS y la limitación honesta)¶

Demostración. Deriva el objetivo e iguálalo a cero. Con $\psi_c = \rho_c'$, las ecuaciones de estimación son $\sum_i \psi_c(u_i)\,\mathbf{x}_i = \mathbf{0}$. Escribe $\psi_c(u) = w(u)\,u$ con el peso

Entonces las ecuaciones de estimación se vuelven $\sum_i w(u_i)\,(Y_i - \mathbf{x}_i'\boldsymbol{\beta})\,\mathbf{x}_i = \mathbf{0}$, que son exactamente las ecuaciones normales de mínimos cuadrados ponderados de 10.4 Mínimos cuadrados ponderados con pesos $w(u_i)$. La trampa es que los pesos dependen de los residuos, que dependen del ajuste. Así que iteramos: empieza desde el ajuste de OLS, calcula los residuos y por tanto los pesos, resuelve el problema de WLS, recalcula los residuos y los pesos, y repite hasta que los coeficientes dejen de moverse. Esto es mínimos cuadrados reponderados iterativamente (IRLS). Un residuo dentro de la banda ($|u| \le c$) mantiene peso 1; un residuo salvaje recibe peso $c/|u|$, que se encoge a medida que crece. $\blacksquare$

La derivación expone la frontera honesta del método. El peso $w(u)$ depende solo del residuo, así que la regresión de Huber protege contra valores atípicos en la dirección de la respuesta, puntos con un gran error vertical. No hace nada respecto a un punto cuyos valores de predictor son extremos, un punto de alto apalancamiento, a menos que ese punto también tenga por casualidad un residuo grande. Un punto de alto apalancamiento a menudo jala la recta hacia sí mismo y por tanto mantiene un residuo pequeño, pasando con peso completo. Protegerse contra eso necesita un estimador de influencia acotada o un estimador MM, que no desarrollamos aquí. Los datos de ahorros hacen concreta la limitación.

R¶

Ajusta el modelo de ahorros del Capítulo 8 de ambas maneras, ordinaria y Huber, sobre los 50 países.

10.7 Resumen del capítulo¶

Ahora puedes responder a un diagnóstico fallido en lugar de solo nombrarlo. Lees un menú de soluciones y haces coincidir cada una con el patrón que la pide: transformar el predictor para una media curvada con dispersión pareja, transformar la respuesta para una curva y un embudo juntos, ponderar para una precisión desigual conocida, usar una pérdida de crecimiento más lento para valores atípicos verticales, y cambiar el modelo cuando la respuesta es un conteo o una proporción. Interpretas los coeficientes logarítmicos como porcentajes en la escala exacta y en la aproximada, eliges una potencia con la verosimilitud perfilada de Box-Cox y su intervalo de confianza para $\lambda$, derivas los mínimos cuadrados ponderados a partir del argumento de Gauss-Markov y los calculas con pesos conocidos, y lees la regresión de Huber como mínimos cuadrados reponderados iterativamente mientras declaras su limitación honesta. Y puedes reconocer, como con los conteos de Galápagos, cuándo ninguna transformación servirá y un modelo distinto es la respuesta correcta.

Resultados clave de un vistazo

Términos clave. Medida correctiva, modelo log-log, modelo log-lineal, elasticidad, semielasticidad, familia de Box-Cox, verosimilitud perfilada, término jacobiano, mínimos cuadrados ponderados, peso, blanqueo, sesgo de retrotransformación, regresión robusta, estimación M, pérdida de Huber, mínimos cuadrados reponderados iterativamente.

Ahora deberías poder

• Elegir una solución (transformar $X$, transformar $Y$, ponderar, o cambiar el modelo) a partir del patrón de diagnóstico que veas.

• Interpretar un coeficiente transformado logarítmicamente en la escala porcentual exacta y aproximada, y derivar ambas lecturas.

• Derivar la log-verosimilitud perfilada de Box-Cox, incluido el término jacobiano, y usarla para elegir una transformación de la respuesta con un intervalo de confianza para $\lambda$.

• Derivar los mínimos cuadrados ponderados a partir del argumento de Gauss-Markov y calcularlos con pesos conocidos.

• Explicar y demostrar la regresión robusta de Huber, y decir contra qué protege y contra qué no.

• Decidir cuándo una transformación es un parche y un modelo distinto es la solución honesta.

Dónde encaja esto. En la columna del flujo de trabajo de El flujo de trabajo del modelado, este capítulo vive en la costura entre CHECK y FIT. El Capítulo 9 (CHECK) te dijo que un supuesto estaba roto; las soluciones aquí te mandan de vuelta a FIT con un modelo mejor especificado, una escala transformada, o un esquema de ponderación, después de lo cual haces CHECK de nuevo y solo entonces USE del resultado. Ese ciclo, diagnosticar luego remediar luego rediagnosticar, es el núcleo de trabajo de la regresión aplicada, y rara vez termina después de una sola pasada. Los hilos continúan: la maquinaria de interpretación logarítmica de 10.2 La transformación logarítmica y la lectura de sus coeficientes regresa como razones de momios en 13.3 Leer los coeficientes como razones de momios y como crecimiento mensual en 15.5 Pronosticar con honestidad: probar en el futuro, el instinto de modelado de la varianza de 10.4 Mínimos cuadrados ponderados regresa como sobredispersión en 14.6 Una familia: el modelo lineal generalizado, y las islas de Galápagos que derrotaron cada transformación aquí reciben su resolución honesta en 14.1 Por qué los conteos rompen el modelo lineal.

10.8 Preguntas frecuentes¶

P1. ¿Debería transformar $X$, transformar $Y$, o ambos? Mira el gráfico de residuos. Una media curvada con dispersión pareja es un problema del predictor: transforma $X$ y deja los errores en paz. Una curva junto con un embudo es un problema de la respuesta: transforma $Y$, lo que arregla ambos a la vez. Cuando solo la varianza se abre en abanico pero la media es recta, prefiere ponderar o una transformación estabilizadora de la varianza de $Y$. Los datos de los mamíferos necesitaron ambas variables en logaritmos porque la relación era multiplicativa en ambas.

P2. ¿Cuándo es seguro el atajo del “porcentaje por unidad”? Cuando el coeficiente es pequeño en magnitud, aproximadamente por debajo de 0.1. Entonces el factor exacto $e^{\beta}$ y la aproximación $1 + \beta$ coinciden hasta una fracción de por ciento (Figure 5). Para un coeficiente como 0.5 o mayor, reporta el porcentaje exacto $100(e^{\beta}-1)$, ya que el atajo puede estar equivocado por varios puntos.

P3. Box-Cox dio $\hat\lambda = 0.43$. ¿Por qué usamos $\lambda = 0.5$ en su lugar? Porque 0.5 está dentro del intervalo de confianza del 95% $[0.23, 0.66]$, así que los datos no lo distinguen del máximo, y la raíz cuadrada es mucho más fácil de interpretar y de justificar que una potencia cruda de 0.43. Redondea a un valor interpretable cercano siempre que el intervalo lo permita.

P4. Mi respuesta tiene ceros o negativos. ¿Aún puedo usar Box-Cox o un logaritmo? No directamente: ambos necesitan $Y > 0$. Las soluciones comunes son un pequeño desplazamiento ($\log(Y + a)$) o, mejor, un modelo construido para el tipo de dato, ya que los ceros en un conteo o una proporción suelen indicar que un modelo lineal generalizado (Capítulos 13 y 14) es la herramienta correcta, no una transformación desplazada.

P5. Si OLS es insesgado incluso bajo varianzas desiguales, ¿para qué molestarse con WLS? OLS sigue siendo insesgado pero deja de ser eficiente: ya no tiene la varianza más pequeña, y sus errores estándar reportados están mal porque suponen un solo $\sigma^2$ común. WLS restaura tanto la eficiencia como los errores estándar honestos usando la precisión conocida de cada punto, como mostró el error estándar más ajustado del intercepto de strongx.

P6. ¿La regresión robusta reemplaza a los diagnósticos? No. Protege contra valores atípicos verticales de forma automática, pero es ciega a los puntos de alto apalancamiento (la lección de Libia), y sus errores estándar son aproximados. Úsala como una herramienta junto a los diagnósticos del Capítulo 9, y siempre investiga por qué un punto es inusual antes de decidir reducir su peso. Un punto con peso reducido a veces es la observación más interesante de los datos.

P7. ¿Cómo obtengo una predicción en la escala original a partir de un modelo transformado? Deshaz la transformación, pero sabe que retrotransformar ingenuamente el valor ajustado está sesgado, como elevar al cuadrado la predicción de Isabela sobrepasó por 46 especies. La media de una función no lineal no es esa función de la media. Para el modelo logarítmico hay una corrección de suavizado estándar; para intervalos de predicción cuidadosos, transforma los extremos del intervalo en lugar de la estimación puntual.

10.9 Problemas de práctica¶

1. (A) Para cada patrón de diagnóstico, nombra la solución que este capítulo recomienda: (i) un gráfico de residuos curvado con dispersión pareja; (ii) una curva y un embudo juntos; (iii) un solo punto con un residuo enorme; (iv) una respuesta de conteo con muchos valores pequeños.

2. (A) Explica la diferencia entre transformar un predictor y transformar la respuesta en términos de lo que cada uno le hace a la función de la media y a la varianza del error.

3. (A) Un modelo log-log tiene pendiente $\beta_1 = 1.3$. Enuncia en palabras qué le hace a $Y$ un aumento del 1% en $X$, y si $Y$ crece más o menos que proporcionalmente.

4. (A) Distingue una elasticidad de una semielasticidad, y di cuál modelo logarítmico produce cada una.

5. (A) En la familia de Box-Cox, ¿a qué transformación corresponde cada uno de $\lambda = 1, 0.5, 0, -1$? ¿Por qué la familia se escribe con el “-1” y la división entre $\lambda$?

6. (A) Explica en una o dos oraciones por qué la log-verosimilitud perfilada de Box-Cox incluye el término jacobiano $(\lambda-1)\sum \log Y_i$, y qué sale mal sin él.

7. (A) Enuncia el criterio de mínimos cuadrados ponderados en palabras, y explica por qué un punto con varianza pequeña debería recibir un peso grande.

8. (A) ¿Por qué los mínimos cuadrados ordinarios siguen siendo insesgados bajo varianzas del error desiguales pero dejan de ser el mejor estimador lineal insesgado?

9. (A) Describe la diferencia entre un valor atípico vertical y un punto de alto apalancamiento, y di contra cuál protege la regresión de Huber.

10. (A) Da dos características de una variable de respuesta que deberían empujarte hacia un modelo distinto en lugar de una transformación, y nombra el capítulo que aporta ese modelo.

11. (B) Partiendo de la relación de la media log-log $Y = e^{\beta_0} X^{\beta_1}$, deriva que multiplicar $X$ por $c$ multiplica $Y$ por $c^{\beta_1}$ de forma exacta (Teorema 10.3), y deriva la lectura aproximada de “$\beta_1$ por ciento por por ciento”.

12. (B) Para el modelo log-lineal, deriva el factor exacto $e^{\beta_1}$ por unidad de $X$ y su aproximación $1 + \beta_1$, y encuentra el valor del coeficiente en el que la aproximación subestima el cambio porcentual verdadero por exactamente un punto porcentual.

13. (B) Deriva la log-verosimilitud perfilada de Box-Cox $\ell_p(\lambda)$ (Teorema 10.6) a partir de la log-verosimilitud completa, mostrando cómo se perfilan $\boldsymbol{\beta}$ y $\sigma^2$ y de dónde viene el término jacobiano.

14. (B) Muestra que la familia de Box-Cox es continua en $\lambda = 0$, es decir, $\lim_{\lambda \to 0}(Y^\lambda - 1)/\lambda = \log Y$.

15. (B) Deriva el estimador de mínimos cuadrados ponderados $\mathbf{b}_W = (\mathbf{X}'\mathbf{W}\mathbf{X})^{-1}\mathbf{X}'\mathbf{W}\mathbf{Y}$ (Teorema 10.8) transformando el modelo con $\mathbf{W}^{1/2}$ y aplicando el teorema de Gauss-Markov. Enuncia la varianza de $\mathbf{b}_W$.

16. (B) Muestra que WLS minimiza $Q_w(\boldsymbol{\beta}) = \sum w_i (Y_i - \mathbf{x}_i'\boldsymbol{\beta})^2$ derivando y obteniendo las ecuaciones normales ponderadas $\mathbf{X}'\mathbf{W}\mathbf{X}\,\mathbf{b}_W = \mathbf{X}'\mathbf{W}\mathbf{Y}$.

17. (B) Para la regresión lineal simple mediante mínimos cuadrados ponderados, deriva la forma cerrada $b_{W,1} = \sum w_i (X_i - \bar{X}_w)(Y_i - \bar{Y}_w) / \sum w_i (X_i - \bar{X}_w)^2$, donde $\bar{X}_w$ y $\bar{Y}_w$ son las medias ponderadas. Define las medias ponderadas.

18. (B) Partiendo de la pérdida de Huber $\rho_c$ (Definición 10.10), calcula $\psi_c = \rho_c'$ y muestra que las ecuaciones de estimación pueden escribirse como ecuaciones normales ponderadas con peso $w(u) = \min(1, c/|u|)$ (Teorema 10.11).

19. (B) Explica, usando el peso $w(u) = \min(1, c/|u|)$, por qué un punto de alto apalancamiento con un residuo pequeño no recibe reducción de peso de la regresión de Huber. Contrasta con un punto de bajo apalancamiento que tiene un residuo grande.

20. (B) Supón que la fila $i$ de un conjunto de datos es el promedio de $n_i$ lecturas independientes cada una de varianza $\tau^2$. Muestra que el peso correcto de WLS es $w_i = n_i$, e identifica la constante común $\sigma^2$ en el modelo $\operatorname{Var}\{\varepsilon_i\} = \sigma^2/w_i$.

21. (C) Ajusta el modelo log-log de los mamíferos en R o Python, reporta $b_0$, $b_1$, y $R^2$, y traduce la pendiente al cambio porcentual exacto en el peso del cerebro para un aumento del 25% en el peso del cuerpo.

22. (C) Usando cars.csv, reproduce la estimación de Box-Cox $\hat\lambda$ y su intervalo del 95%, luego ajusta tanto dist ~ speed como sqrt(dist) ~ speed y compara sus gráficos de residuos contra ajustados. ¿Qué modelo satisface mejor la varianza constante?

23. (C) Usando strongx.csv, ajusta OLS y WLS con pesos $1/\text{sd}^2$, reporta ambas pendientes y ambos errores estándar del intercepto, y confirma los coeficientes de WLS haciendo la regresión de las variables escaladas por $\sqrt{w}$ sin intercepto.

24. (C) En strongx.csv, agrega un término cuadrático (crossx ~ energy + I(energy^2)) al ajuste de WLS y compáralo con el ajuste lineal de WLS. ¿El término de curvatura mejora el ajuste, y qué dice eso sobre el modelo lineal?

25. (C) Usando gala.csv, ajusta Species y sqrt(Species) sobre los cinco predictores geográficos, compara $R^2$ y los gráficos de residuos, y reporta el valor de apalancamiento de Isabela en ambos ajustes. Explica por qué la transformación no reduce el valor de apalancamiento de Isabela.

26. (C) Predice el conteo de especies de Fernandina a partir del modelo de raíz cuadrada de gala, retrotransforma elevando al cuadrado, y compara con su conteo real. Comenta la dirección del sesgo de retrotransformación.

27. (C) Usando savings.csv, ajusta OLS y regresión de Huber de sr sobre los cuatro predictores, lista los cuatro países con los pesos de Huber más pequeños, y cruza cada uno contra su valor de apalancamiento. Identifica un país de alto apalancamiento al que Huber no reduce el peso.

28. (C) Simula datos heterocedásticos (semilla 4210): 40 puntos con $\operatorname{Var}\{\varepsilon_i\}$ proporcional a $X_i^2$, ajusta OLS y WLS con pesos $1/X_i^2$, y compara las dos estimaciones de la pendiente y sus errores estándar a lo largo de 500 réplicas. ¿Cuál es más precisa, y alguna muestra sesgo?

10.10 Práctica de examen¶

Estas cinco preguntas están escritas en el estilo de los exámenes del curso: cada una te pide explicar, evaluar o interpretar en oraciones completas, no solo producir un número. Donde una pregunta muestra salida de software, los números se produjeron en la máquina del curso (R 4.6.0) leyendo los mismos archivos CSV de data/ que usaste todo el semestre; Python con statsmodels da los mismos valores. Lee la salida, luego responde en oraciones completas con unidades. Intenta cada pregunta antes de abrir su respuesta modelo.

EP 10.1 (evaluar una afirmación). Un compañero ajusta un modelo log-lineal de las especies de plantas de Galápagos sobre la elevación de la isla y lee la pendiente como un efecto porcentual.

gala <- read.csv("data/gala.csv")
fit <- lm(log(Species) ~ Elevation, data = gala)
round(coef(fit), 6)

(Intercept)   Elevation
   2.591399    0.002490

El compañero escribe: “El coeficiente 0.00249 significa cerca de $0.25\%$ más especies por metro de elevación, así que una isla 100 metros más alta tiene cerca de $25\%$ más especies”. Explica qué significa el coeficiente en la escala porcentual, luego evalúa ambas lecturas: la de por metro y la de 100 metros. ¿Cuál es correcta, cuál no, y cuál es la cifra correcta para una diferencia de 100 metros?

EP 10.2 (interpretar salida en contexto). Box-Cox corrido sobre la regresión cruda de los mamíferos del peso del cerebro sobre el peso del cuerpo devuelve la salida de abajo.

mammals <- read.csv("data/mammals.csv")
library(MASS)
bc <- boxcox(lm(brain ~ body, data = mammals),
             lambda = seq(-0.5, 0.5, by = 0.01), plotit = FALSE)
lambda_hat <- bc$x[which.max(bc$y)]
ci <- range(bc$x[bc$y > max(bc$y) - 0.5 * qchisq(0.95, 1)])
round(c(lambda_hat = lambda_hat, ci_low = ci[1], ci_high = ci[2]), 3)

lambda_hat     ci_low    ci_high
      0.07      -0.03       0.18

Interpreta $\hat\lambda$ y su intervalo en contexto. ¿Qué transformación de la respuesta recomienda la salida, y por qué? La apertura del capítulo tomó el logaritmo de ambos, el peso del cerebro y el del cuerpo; ¿este resultado de Box-Cox, que transforma solo la respuesta, justifica todo ese movimiento log-log? Explica qué te dice Box-Cox aquí y qué no.

EP 10.3 (qué cambiaría si). Un estudiante ajusta mínimos cuadrados ponderados a los datos de física strongx con pesos $w_i = 1/\text{sd}_i^2$, luego multiplica cada peso por 10 y reajusta, esperando que la recta “confíe diez veces más en los puntos precisos” y se desplace.

strongx <- read.csv("data/strongx.csv")
w <- 1 / strongx$sd^2
fit1 <- lm(crossx ~ energy, data = strongx, weights = w)
fit2 <- lm(crossx ~ energy, data = strongx, weights = 10 * w)
rbind(w = coef(fit1), tenw = coef(fit2))
c(se_energy_w = summary(fit1)$coef["energy", "Std. Error"],
  se_energy_10w = summary(fit2)$coef["energy", "Std. Error"])

     (Intercept)   energy
w       148.4732 530.8354
tenw    148.4732 530.8354
  se_energy_w se_energy_10w
        47.55         47.55

La pendiente, el intercepto, y el error estándar de la pendiente son idénticos hasta el último dígito mostrado. Explica por qué multiplicar cada peso por la misma constante no cambia nada que reportarías, refiriéndote al modelo WLS $\operatorname{Var}\{\varepsilon_i\} = \sigma^2/w_i$. ¿Qué única cantidad sí cambia, y por qué no es un resultado?

EP 10.4 (interpretar salida en contexto). La salida de abajo muestra los valores ajustados de la regresión lineal simple de las especies de plantas de Galápagos sobre los cinco predictores geográficos, para las cinco islas con los valores ajustados más pequeños.

gala <- read.csv("data/gala.csv")
fit_lin <- lm(Species ~ Area + Elevation + Nearest + Scruz + Adjacent, data = gala)
pred <- predict(fit_lin)
head(data.frame(island = gala$island, actual = gala$Species,
                fitted = round(pred, 1))[order(pred), ], 5)

      island actual fitted
   Coamano      2  -36.4
    Darwin     10   -9.0
 Bartolome     31   -7.3
  Gardner1     58   -4.0
  Genovesa     40   -0.5

Interpreta esta salida en el contexto de elegir una solución. ¿Qué es imposible en estas predicciones, por qué sucede, y qué te dice sobre si alguna transformación de la respuesta podría rescatar el modelo lineal normal? Nombra la solución honesta y el capítulo que la aporta.

EP 10.5 (evaluar una afirmación). Un analista reajusta el modelo de ahorros con la regresión robusta de Huber y le dice a un colega: “El ajuste de regresión robusta reduce el peso de cada país inusual de forma automática, así que ya no necesito los diagnósticos del Capítulo 9”. La salida lista los cuatro países de mayor apalancamiento con sus pesos de Huber; para contraste, Zambia tiene un valor de apalancamiento de 0.064 y peso de Huber 0.472.

savings <- read.csv("data/savings.csv")
fit_ols <- lm(sr ~ pop15 + pop75 + dpi + ddpi, data = savings)
fit_hub <- MASS::rlm(sr ~ pop15 + pop75 + dpi + ddpi, data = savings)
h <- hatvalues(fit_ols); wt <- fit_hub$w
data.frame(country = savings$country, leverage_value = round(h, 3),
           huber_weight = round(wt, 3))[order(-h)[1:4], ]

       country leverage_value huber_weight
         Libya          0.531        1.000
 United States          0.334        1.000
         Japan          0.223        0.879
       Ireland          0.212        1.000

Evalúa la afirmación del analista usando estos números. ¿Qué país es el contraejemplo más claro, qué limitación general de la regresión de Huber expone, y qué debería aún hacer el analista?

Juego del capítulo¶