12. Multicolinealidad, selección de variables y validación (en espanol)

12. Multicolinealidad, selección de variables y validación¶

Medir bien la grasa corporal es una molestia. La forma precisa es pesar a una persona bajo el agua y calcular su densidad, lo que requiere un tanque, un técnico y un sujeto cooperativo. La forma barata es una cinta métrica. Así surge una pregunta natural: ¿pueden un puñado de medidas con cinta (cintura, pecho, muslo, muñeca y las demás) predecir el porcentaje de grasa corporal que da el pesaje bajo el agua? Un estudio de 1985 registró ambas cosas para 252 hombres, y ese conjunto de datos, fat.csv, es con el que este capítulo trabaja de principio a fin. (Este es el estudio mayor de 252 hombres con trece medidas, no la tabla de grasa corporal de veinte hombres y tres medidas que conociste para las sumas de cuadrados extra en 8.5 Gráficos de variable agregada; ambos comparten una columna “thigh” pero son datos distintos.)

Hay un plan obvio: hacer la regresión de la grasa corporal sobre las trece medidas y leer los coeficientes. Cuando lo haces, el modelo explica cerca del 75 por ciento de la variación, lo que suena a éxito. Luego miras los coeficientes y algo anda mal. El coeficiente del peso es negativo, como si los hombres más pesados tuvieran menos grasa una vez que fijas las demás medidas. Varios predictores que todos saben que importan salen estadísticamente no significativos. Quita una medida y las otras se tambalean. El problema es que las medidas de cinta miden todas más o menos lo mismo, el tamaño del cuerpo, así que están enredadas entre sí. Figure 1 muestra justo cuán enredadas.

Este capítulo trata de convivir con demasiados predictores, la mayoría redundantes. Tres preguntas lo organizan. Primera, ¿cómo diagnosticamos y describimos el daño que causan los predictores correlacionados, usando el factor de inflación de la varianza? Segunda, cuando hay varios modelos sobre la mesa, ¿cómo elegimos, usando criterios como el $C_p$ de Mallows y la validación cruzada que premian la predicción y no el mero ajuste? Tercera, ¿qué herramientas modernas (el lasso y la regresión ridge) esquivan todo el problema de selección encogiendo los coeficientes en lugar de eliminarlos? Por el camino conocerás la advertencia más importante de la estadística aplicada: un modelo elegido para ajustar tus datos siempre se verá mejor en esos datos de lo que merece.

12.1 La multicolinealidad y el factor de inflación de la varianza¶

Intuición¶

Los capítulos 7 y 8 construyeron una regresión múltiple y leyeron cada coeficiente como una pendiente parcial, el efecto de un predictor con los demás fijos (8.5 Gráficos de variable agregada). Esa lectura suponía calladamente que los predictores llevaban información separada. La viñeta mostró el daño cuando no es así; ahora nombramos la causa, la multicolinealidad (Definición 12.1).

Un poco de ella es inofensiva y normal. Mucha es un problema, porque el modelo no puede distinguir entre los predictores correlacionados. Imagina dos predictores que se mueven casi al unísono, como el peso y el perímetro de la cadera. Los datos muestran qué le pasa a la grasa corporal cuando ambos suben juntos, pero casi no muestran nada sobre qué pasa cuando uno sube y el otro se queda fijo, porque esa combinación apenas ocurre. Así que el modelo no tiene una base firme para repartir el efecto compartido entre ellos. Puede darle al peso un coeficiente positivo grande y a la cadera uno negativo grande, o al revés, y ajustar los datos más o menos igual de bien en cualquier caso. Figure 2 dibuja la escena.

El factor de inflación de la varianza (Definición 12.2) mide exactamente cuánto se infla la varianza.

El modelo completo, y por qué se portan mal sus coeficientes¶

Ajusta primero el modelo completo, para que los síntomas sean concretos. La respuesta es brozek, el porcentaje de grasa corporal por la fórmula de Brozek, y los trece predictores son las medidas de cinta y báscula. Recuerda de 4.2 La correlación y la pendiente de regresión que para un solo predictor, $R^2$ es simplemente la correlación al cuadrado; con muchos predictores $R^2$ sigue subiendo a medida que agregas columnas, ayuden o no, así que un $R^2$ alto por sí solo no prueba nada.

Viste una versión pequeña de esto en 8.5 Gráficos de variable agregada, donde una tabla de grasa corporal de veinte hombres tenía una medida de tríceps y una de muslo tan correlacionadas que sus coeficientes se volvieron inestables. Aquí aparece la misma enfermedad a plena escala, con trece predictores enredados en vez de dos. Un vistazo rápido a las correlaciones entre las medidas de tamaño confirma el diagnóstico.

round(cor(fat[, c("weight", "chest", "abdom", "hip", "thigh")]), 3)

       weight chest abdom   hip thigh
weight  1.000 0.894 0.888 0.941 0.869
chest   0.894 1.000 0.916 0.829 0.730
abdom   0.888 0.916 1.000 0.874 0.767
hip     0.941 0.829 0.874 1.000 0.896
thigh   0.869 0.730 0.767 0.896 1.000

print(fat[["weight", "chest", "abdom", "hip", "thigh"]].corr().round(3))

        weight  chest  abdom    hip  thigh
weight   1.000  0.894  0.888  0.941  0.869
chest    0.894  1.000  0.916  0.829  0.730
abdom    0.888  0.916  1.000  0.874  0.767
hip      0.941  0.829  0.874  1.000  0.896
thigh    0.869  0.730  0.767  0.896  1.000

Cada par correlaciona por encima de 0.7, y el peso con la cadera llega a 0.941. Estas columnas son casi duplicados unas de otras.

Fórmula¶

En palabras: haz la regresión del predictor $k$ sobre los predictores restantes, ve qué tan bien lo predicen, y $\mathrm{VIF}_k$ es uno sobre la fracción que sobra. Si $X_k$ no se relaciona con los demás, $R_k^2 = 0$ y $\mathrm{VIF}_k = 1$. Si los demás explican el 90 por ciento de él, $R_k^2 = 0.9$ y $\mathrm{VIF}_k = 10$. Una regla práctica común señala $\mathrm{VIF}_k > 10$ (equivalentemente $R_k^2 > 0.9$) como grave, aunque el corte es una convención, no una ley.

Vale la pena mirar la forma de esa fórmula, porque explica por qué un poco de colinealidad es inofensiva y mucha es un desastre. Figure 3 grafica el VIF contra $R_k^2$. La curva es casi plana mientras un predictor es en su mayoría propio, luego se dobla bruscamente hacia arriba a medida que los demás se acercan a reproducirlo. El peso, del cual las otras doce medidas explican el 97 por ciento, se sitúa muy arriba en la parte empinada.

El nombre es literal. El VIF es el factor por el cual se multiplica la varianza de $b_k$, en comparación con la varianza que tendría si $X_k$ no estuviera correlacionado con los demás predictores:

En palabras: la varianza de un coeficiente es el usual $\sigma^2 / S_{kk}$ (pequeño cuando el predictor está muy disperso) por el factor de inflación. La colinealidad entra solo a través de ese último factor.

Deducción (el VIF a partir de la regresión auxiliar)¶

Demostración. Fija un predictor $X_k$. Por la construcción de variable agregada de 8.5 Gráficos de variable agregada, el coeficiente de mínimos cuadrados $b_k$ es igual a la pendiente de una regresión simple de $Y$ sobre la parte de $X_k$ que los demás predictores no explican. Haz concreta esa parte: ejecuta la regresión auxiliar de $X_k$ sobre todos los demás predictores, y llama a sus residuos $\tilde{X}_{ik}$. Estos residuos son $X_k$ con la información de los demás predictores removida.

El resultado de variable agregada dice que $b_k = \sum_i \tilde{X}_{ik} Y_i / \sum_i \tilde{X}_{ik}^2$, una pendiente de regresión simple con los $\tilde{X}_{ik}$ como predictor. De 2.5 Comportamiento muestral y el teorema de Gauss-Markov, una pendiente de regresión simple $\sum k_i Y_i$ con pesos $k_i = \tilde{X}_{ik} / \sum_j \tilde{X}_{jk}^2$ tiene varianza $\sigma^2 \sum_i k_i^2 = \sigma^2 / \sum_i \tilde{X}_{ik}^2$. Así que

Ahora identifica el denominador. La regresión auxiliar tiene suma total de cuadrados $\sum_i (X_{ik} - \bar{X}_k)^2 = S_{kk}$ y suma de cuadrados del error $\sum_i \tilde{X}_{ik}^2$. Por la definición de $R^2$ aplicada a esa regresión auxiliar,

Sustituye:

El primer factor $\sigma^2 / S_{kk}$ es lo que sería la varianza si $X_k$ estuviera libre de los demás ($R_k^2 = 0$). El segundo factor $1/(1 - R_k^2) = \mathrm{VIF}_k$ es el precio de la colinealidad. Cuando los demás predictores casi reproducen $X_k$, $R_k^2 \to 1$ y la varianza explota. $\blacksquare$

R y Python¶

El paquete car de R calcula todos los VIF a la vez. La fórmula de arriba te permite comprobar cualquiera de ellos a mano a partir de una sola regresión auxiliar.

Remedios¶

Hay cuatro respuestas honestas a los VIF altos. Primera, no hacer nada, si solo te importa la predicción dentro del rango de los datos: la colinealidad no sesga $\hat{Y}$ ni ensancha mucho su intervalo, así que un modelo predictivo puede llevar predictores redundantes sin peligro. Segunda, quitar predictores: si el peso, la cadera y el abdomen dicen casi lo mismo, conserva el que puedas medir mejor y descarta el resto. Tercera, combinarlos en un solo índice (un promedio, o una componente principal) que capture la señal de tamaño compartida. Cuarta, encoger los coeficientes con regresión ridge (Sección 12.5), que tolera la colinealidad por diseño. Lo que no debes hacer es leer un coeficiente colineal como un efecto causal real y actuar según su signo.

Leer sobre varianza inflada es una cosa; aquí está la misma enfermedad con una perilla, para que decidas por ti mismo cuánta correlación es demasiada.

12.2 Elegir entre modelos: criterios de selección¶

Intuición¶

El modelo completo tiene predictores redundantes; un modelo más pequeño podría predecir igual de bien y ser más fácil de confiar. ¿Pero cuál modelo más pequeño? Con trece predictores hay $2^{13} = 8192$ subconjuntos posibles. Necesitamos una puntuación que los clasifique. La puntuación obvia, $R^2$, es inútil para esto, porque $R^2$ nunca disminuye cuando agregas un predictor: el modelo más grande siempre gana, incluso si los predictores agregados son ruido. Recuerda de 4.2 La correlación y la pendiente de regresión que $R^2$ es la correlación al cuadrado, una medida de ajuste, no de predicción. Un buen criterio de selección debe premiar el ajuste cobrando a la vez una penalización por cada parámetro extra, de modo que un predictor gane su lugar solo si se paga a sí mismo.

Cuatro criterios hacen esto, cada uno en su propia moneda. El $R^2$ ajustado descuenta $R^2$ por los grados de libertad gastados. El $C_p$ de Mallows estima el error total de predicción y lo compara con el número de parámetros. AIC y BIC intercambian bondad de ajuste contra una penalización por tamaño de la teoría de la información. PRESS predice cada caso a partir de un modelo ajustado sin él. Suelen coincidir en el mensaje general (quita el peso muerto) mientras discrepan en el tamaño exacto, lo cual es en sí una lección: no hay un único modelo correcto.

Fórmula¶

Sea un modelo candidato con $p$ parámetros (predictores más intercepto), suma de cuadrados del error $\mathrm{SSE}_p = \sum_i (Y_i - \hat{Y}_i)^2$, y $\mathrm{SSTO} = \sum_i (Y_i - \bar{Y})^2$. Sea $\hat{\sigma}^2$ el cuadrado medio del error del modelo más grande (completo), nuestra mejor estimación de la varianza del ruido $\sigma^2$. Los cuatro criterios son

• $R^2_{a,p}$ ($R^2$ ajustado): como $R^2$ pero con $\mathrm{SSE}$ y $\mathrm{SSTO}$ divididos por sus grados de libertad. Agregar un predictor inútil lo baja. Más grande es mejor.

• $C_p$ (Mallows): más pequeño es mejor; un modelo sin sesgo tiene $C_p \approx p$, así que buscas modelos con $C_p$ bajo que también se sitúen cerca de la línea $C_p = p$.

• $\mathrm{AIC}_p$, $\mathrm{BIC}_p$: término de ajuste más penalización; más pequeño es mejor. La penalización de BIC $p \ln n$ es más pesada que la de AIC $2p$ siempre que $n > 7$, así que BIC prefiere modelos más pequeños.

Esa última afirmación es fácil de ver en una imagen. Figure 6 dibuja las dos penalizaciones por tamaño como líneas rectas. Con $n = 252$ hombres, BIC cobra alrededor de $\ln 252 \approx 5.5$ por cada parámetro extra mientras que AIC cobra solo 2, así que la línea de BIC sube casi tres veces más empinada. Un predictor tiene que mejorar el ajuste más para pagar su peaje BIC más pesado, que es exactamente por qué BIC aterriza en modelos más pequeños.

El quinto criterio, PRESS, recibe su propia subsección porque mide algo distinto: el error honesto fuera de la muestra.

donde $\hat{Y}_{(i)}$ es la predicción para el caso $i$ del modelo reajustado con el caso $i$ removido. En palabras: oculta cada observación, predícela a partir del resto, y suma los fallos al cuadrado. Más pequeño es mejor, y a diferencia de $\mathrm{SSE}$, PRESS no puede reducirse solo agregando predictores, porque un predictor basura ayuda al ajuste pero perjudica la predicción del caso dejado fuera.

Deducción (el $C_p$ de Mallows)¶

Demostración. Queremos un criterio que sea pequeño cuando los valores ajustados $\hat{Y}_i$ de un modelo estén cerca de las medias verdaderas $\mu_i = E\{Y_i\}$, a través de los $n$ casos. Define el objetivo

el error cuadrático medio total de los valores ajustados, estandarizado por $\sigma^2$. Divide cada término en varianza y sesgo al cuadrado con la identidad $E\{W^2\} = \operatorname{Var}\{W\} + (E\{W\})^2$ aplicada a $W = \hat{Y}_i - \mu_i$:

Para el término de varianza, los valores ajustados son $\hat{\mathbf{Y}} = \mathbf{H}\mathbf{Y}$ con la matriz sombrero $\mathbf{H}$ del modelo candidato, y $\sum_i \operatorname{Var}\{\hat{Y}_i\} = \sigma^2\,\mathrm{tr}(\mathbf{H}) = p\sigma^2$, usando $\mathrm{tr}(\mathbf{H}) = p$ de 7.3 La matriz sombrero. Así que

Ahora conecta $B$ con algo que podamos calcular, la suma de cuadrados del error esperada. Escribe $\mathbf{Y} = \boldsymbol{\mu} + \boldsymbol{\varepsilon}$ y $\mathrm{SSE}_p = \|(\mathbf{I} - \mathbf{H})\mathbf{Y}\|^2$. Entonces

porque $(\mathbf{I}-\mathbf{H})\boldsymbol{\mu}$ tiene $i$-ésima entrada $\mu_i - E\{\hat{Y}_i\}$ así que su longitud al cuadrado es $B$, y $\mathbf{I}-\mathbf{H}$ es idempotente con traza $n - p$ (de nuevo 7.3 La matriz sombrero). Despeja el sesgo: $B = E\{\mathrm{SSE}_p\} - (n - p)\sigma^2$. Sustituye en $\Gamma_p$:

Reemplaza el desconocido $E\{\mathrm{SSE}_p\}$ por el observado $\mathrm{SSE}_p$ y el desconocido $\sigma^2$ por el $\hat{\sigma}^2$ del modelo completo, y tienes el estadístico de Mallows:

Salen dos consecuencias. Si el modelo candidato no tiene sesgo ($B = 0$, es decir, contiene todos los predictores que importan), entonces $E\{\mathrm{SSE}_p\} = (n - p)\sigma^2$ y $\Gamma_p = (n - p) - (n - 2p) = p$. Así que los modelos sin sesgo tienen $C_p \approx p$: buscas modelos donde $C_p$ sea a la vez pequeño y cercano a $p$. Un modelo que omite un predictor importante tiene $B > 0$, lo que empuja $C_p$ bien por encima de $p$. $\blacksquare$

R y Python¶

El paquete leaps de R busca los 2¹³ subconjuntos e informa el mejor modelo de cada tamaño.

Deducir y comprobar PRESS¶

PRESS parece caro: reajustar el modelo $n$ veces, una por cada caso eliminado, son 252 ajustes por modelo aquí. Una identidad hermosa lo hace gratis. Puedes obtener cada predicción eliminada a partir de un solo ajuste con todos los datos, usando el residuo ordinario y el valor de apalancamiento.

Demostración. Escribe la matriz de diseño $\mathbf{X}$ con $i$-ésima fila $\mathbf{x}_i'$, la estimación con todos los datos $\mathbf{b} = (\mathbf{X}'\mathbf{X})^{-1}\mathbf{X}'\mathbf{Y}$, el valor ajustado $\hat{Y}_i = \mathbf{x}_i'\mathbf{b}$, el residuo $e_i = Y_i - \hat{Y}_i$, y el valor de apalancamiento $h_{ii} = \mathbf{x}_i'(\mathbf{X}'\mathbf{X})^{-1}\mathbf{x}_i$ de 7.3 La matriz sombrero (la diagonal de la matriz sombrero). Sea $\mathbf{b}_{(i)}$ la estimación a partir de los datos con el caso $i$ removido, y $\hat{Y}_{(i)} = \mathbf{x}_i'\mathbf{b}_{(i)}$ la predicción eliminada. Mostramos

Eliminar el caso $i$ cambia la matriz de productos cruzados y el vector a $\mathbf{X}'\mathbf{X} - \mathbf{x}_i\mathbf{x}_i'$ y $\mathbf{X}'\mathbf{Y} - \mathbf{x}_i Y_i$. La identidad de actualización del inverso de Sherman-Morrison (deducida y comprobada en 9.3 Influencia: qué puntos cambian realmente el ajuste) da

Escribe $\mathbf{u} = (\mathbf{X}'\mathbf{X})^{-1}\mathbf{x}_i$, de modo que $\mathbf{x}_i'\mathbf{u} = h_{ii}$ y $\mathbf{x}_i'(\mathbf{X}'\mathbf{X})^{-1}\mathbf{X}'\mathbf{Y} = \mathbf{x}_i'\mathbf{b} = \hat{Y}_i$. Multiplica el inverso actualizado por $\mathbf{X}'\mathbf{Y} - \mathbf{x}_i Y_i$:

donde el numerador colapsó a $\hat{Y}_i - Y_i = -e_i$. Ahora forma la predicción eliminada:

Así que cada residuo eliminado es simplemente el residuo ordinario dividido por $1 - h_{ii}$, y

Los puntos de alto apalancamiento (grande $h_{ii}$) ven sus residuos magnificados, lo cual es correcto: esos son los casos sobre los que el modelo se apoya, y los predice peor cuando se remueven. $\blacksquare$

Antes de aceptar el veredicto de un solo criterio, recorre tú mismo la escalera completa de tamaños de modelo y observa cómo las cuatro puntuaciones tiran en direcciones distintas.

12.3 Los peligros de la selección automática¶

Intuición¶

Es tentador automatizar la búsqueda: dejar que la computadora agregue y quite predictores por sus valores $p$ hasta que nada mejore. Esto es la selección paso a paso (Definición 12.7), y está disponible en una línea en todo paquete estadístico. También es la fuente de más mala ciencia que casi cualquier otro procedimiento rutinario, y deberías entender exactamente por qué antes de usarla.

Dos cosas salen mal. La primera es el sobreajuste (Definición 12.8): un procedimiento que persigue el mejor ajuste en tu muestra se aferrará a accidentes de esa muestra, patrones que son ruido y no se repetirán. La segunda, más sutil y peor, es la inferencia posterior a la selección (Definición 12.9): una vez que has elegido un modelo porque ajustó bien, los valores $p$, los intervalos de confianza y las pruebas $F$ que el software imprime para ese modelo están mal. Se dedujeron suponiendo que el modelo estaba fijo de antemano, antes de mirar los datos. Elegir el modelo mirando los datos rompe esa suposición, y la significancia reportada se vuelve ficción.

Una demostración con puro ruido¶

La forma más limpia de ver el desastre es ejecutar la selección sobre datos sin señal alguna. Toma una respuesta y cuarenta predictores que sean todos ruido aleatorio independiente, sin relación con la respuesta ni entre sí. Una prueba honesta casi nunca debería encontrar nada. La selección encuentra bastante.

set.seed(4210)
noise <- data.frame(y = rnorm(100), matrix(rnorm(100 * 40), 100, 40))
names(noise)[-1] <- paste0("x", 1:40)
pvals <- sapply(1:40, function(j)
  summary(lm(noise$y ~ noise[[paste0("x", j)]]))$coefficients[2, 4])
best5 <- order(pvals)[1:5]
selected <- lm(y ~ ., data = noise[, c("y", paste0("x", best5))])
fstat <- summary(selected)$fstatistic
c(R2 = summary(selected)$r.squared,
  overall_F_pvalue = as.numeric(pf(fstat[1], fstat[2], fstat[3], lower.tail = FALSE)))

              R2 overall_F_pvalue
      0.13497399       0.01663716

Elegimos los cinco predictores con los valores $p$ individuales más pequeños de cuarenta columnas de puro ruido, ajustamos un modelo solo sobre esos cinco, y obtuvimos un valor $p$ de la prueba $F$ general de 0.017. Por la regla habitual, ese modelo es “significativo al nivel 0.05”, y su $R^2 = 0.135$ se ve respetable. Todo ello es una ilusión: no hay relación alguna en los datos. El paso de selección fabricó la significancia escogiendo a dedo, y la prueba $F$, ciega al escogido a dedo, la reportó como real.

Esto no es una casualidad única de la semilla. Figure 10 repite todo el experimento 2000 veces. Bajo una prueba honesta el valor $p$ de la prueba $F$ general debería ser uniforme en $[0, 1]$, así que solo el 5 por ciento debería caer por debajo de 0.05. Tras la selección, el 94 por ciento lo hace.

Qué hacer en su lugar¶

Nada de esto significa que la selección de variables esté prohibida. Significa tres cosas. Primera, prefiere los métodos basados en criterios de la Sección 12.2 (compara modelos completos por $C_p$, BIC, o error validado cruzadamente) sobre la búsqueda paso a paso guiada por valores $p$, e informa que buscaste. Segunda, nunca cites los valores $p$ de un modelo seleccionado como si el modelo hubiera estado fijo de antemano; recuerda de 8.4 La prueba lineal general que la distribución de la prueba $F$ supone que la comparación se eligió antes de ver los datos. Tercera, y más confiable de todo, juzga el modelo final con datos con los que no fue elegido. Eso es validación, el tema de la siguiente sección, y es el único árbitro honesto para un modelo que construiste buscando.

12.4 Validación: entrenar, probar y validar cruzadamente¶

Intuición¶

Los residuos de un modelo ajustado te dicen qué tan bien ajusta los datos sobre los que fue construido, lo cual siempre es halagador. Para saber qué tan bien predice un modelo, debes probarlo con datos que nunca ha visto. La idea hace eco del pensamiento fuera de la muestra detrás del bootstrap en 5.4 El bootstrap para la regresión y de los intervalos de predicción de 3.5 Intervalo de predicción para una observación nueva: las estimaciones honestas del error vienen de tratar algunos datos como genuinamente nuevos.

La versión más simple divide los datos una vez: un conjunto de entrenamiento (Definición 12.10) para ajustar el modelo, y un conjunto de prueba, apartado, para medir su error. La división debe hacerse antes de todo modelado, y el conjunto de prueba nunca debe influir en una sola decisión. Una versión más eficiente, la validación cruzada de $k$ particiones (Definición 12.11), divide los datos en $k$ particiones iguales, luego rota: cada partición toma su turno como conjunto de prueba mientras las otras $k - 1$ particiones entrenan el modelo. Promediar los $k$ errores de prueba usa cada caso para probar exactamente una vez, así que no desperdicia datos. Figure 11 muestra la rotación.

Medimos el error con la raíz del error cuadrático medio (Definición 12.12), el tamaño típico de un fallo de predicción, en las unidades de la respuesta (puntos porcentuales de grasa corporal). Una más pequeña significa que el modelo está más cerca, en promedio, de la verdad para un hombre nuevo.

La división única entrenamiento/prueba¶

Divide los 252 hombres con una semilla fija para que el análisis se reproduzca, ajusta sobre el conjunto de entrenamiento, y compara el error de entrenamiento con el error de prueba.

Validación cruzada de $k$ particiones, escrita a mano¶

Una sola división es ruidosa: un conjunto de prueba afortunado o desafortunado puede engañar. La validación cruzada promedia sobre $k$ divisiones. Vale la pena escribir el bucle tú mismo una vez, para que el mecanismo no guarde misterio. Asigna a cada caso un número de partición, luego haz el bucle: entrena con las otras particiones, predice la partición apartada, recolecta los errores al cuadrado.

El Ejemplo 12.5 reportó la RMSE de prueba de una sola división. Este widget repite esa división mil veces, para que veas todas las demás respuestas que igual de fácil te pudo haber dado.

12.5 Un vistazo a la contracción: ridge y lasso¶

Intuición¶

La selección es un instrumento tosco: un predictor está o dentro (coeficiente completo) o fuera (cero). La contracción (Definición 12.13) ofrece una alternativa más suave. En vez de eliminar predictores, los conserva todos pero jala sus coeficientes hacia cero, cambiando un poco de sesgo por una gran caída de varianza. Cuando los predictores son colineales, los mínimos cuadrados ordinarios pueden producir coeficientes enormes que se cancelan entre sí (el peso negativo del Ejemplo 12.1). La contracción se niega a dejar que los coeficientes crezcan tanto, lo que los estabiliza.

Fórmula¶

Los mínimos cuadrados ordinarios minimizan la suma de residuos al cuadrado. Ridge (Definición 12.14) y lasso (Definición 12.15) agregan una penalización sobre el tamaño de los coeficientes.

• El primer término es la pérdida usual de mínimos cuadrados; quiere un buen ajuste.

• El segundo término es la penalización; quiere coeficientes pequeños. La suma corre sobre los $q$ coeficientes de pendiente (no el intercepto), y los predictores se estandarizan primero para que la penalización los trate de forma justa.

• $\lambda \ge 0$ es el parámetro de ajuste que fija la fuerza de la penalización. En $\lambda = 0$ ambos se reducen a mínimos cuadrados; a medida que $\lambda \to \infty$ todos los coeficientes son aplastados a cero. Elegimos $\lambda$ por validación cruzada.

Ridge usa la penalización al cuadrado ($\sum \beta_k^2$, la norma $L_2$); lasso usa la penalización de valor absoluto ($\sum |\beta_k|$, la norma $L_1$). La diferencia suena menor y no lo es. Ridge encoge cada coeficiente hacia cero pero rara vez exactamente a cero. Lasso puede fijar coeficientes exactamente a cero, así que hace selección de variables y contracción a la vez. La razón es geométrica, mostrada en Figure 15: la región de restricción del lasso es un rombo con esquinas sobre los ejes, y el contorno de mejor ajuste tiende a tocar una esquina, donde un coeficiente es cero. La región de restricción de ridge es un disco suave sin esquinas.

R y Python¶

El paquete de R glmnet es la herramienta estándar. Su cv.glmnet ajusta toda la trayectoria de valores de $\lambda$ y elige uno por validación cruzada, devolviendo lambda.min (el minimizador del error) y lambda.1se (el mayor $\lambda$ dentro de un error estándar del mínimo, una elección deliberadamente más dispersa). Los predictores se estandarizan automáticamente.

El resumen sencillo de la contracción¶

Ridge y lasso responden a los problemas de colinealidad y selección en un solo movimiento. Ridge conserva cada predictor pero encoge los coeficientes de modo que dos predictores colineales no puedan estallar uno contra otro; es la herramienta natural cuando crees que muchos predictores contribuyen cada uno un poco. Lasso encoge y selecciona a la vez, poniendo en cero los coeficientes que no necesita, así que es la herramienta natural cuando crees que solo unos pocos predictores importan y quieres que el modelo los encuentre. Ambos reemplazan la pregunta discreta e inestable “¿cuáles predictores están dentro?” con un dial suave y ajustable $\lambda$, fijado por validación cruzada. Por eso la contracción, no la selección paso a paso, es el estándar moderno cuando el objetivo es predecir.

12.6 Resumen del capítulo¶

Ahora puedes manejar una regresión con más predictores de los que puedes confiar. Diagnosticaste la multicolinealidad con el factor de inflación de la varianza y lo dedujiste de una regresión auxiliar; comparaste modelos candidatos con $R^2$ ajustado, $C_p$ de Mallows, AIC, BIC y PRESS, deduciendo $C_p$ y demostrando el atajo del residuo eliminado que hace gratis el dejar uno fuera; viste, en una demostración de puro ruido, por qué la selección paso a paso sobreajusta y por qué los valores $p$ posteriores a la selección mienten; escribiste un bucle de validación cruzada de $k$ particiones a mano y usaste un conjunto de prueba apartado para juzgar un modelo honestamente; y describiste y ejecutaste ridge y lasso, leyendo la contracción como un intercambio sesgo-varianza hecho por una penalización ajustable. En los datos de grasa corporal el candidato de ocho predictores venció al modelo completo en cada medida honesta: menor $C_p$, menor PRESS, y menor error validado cruzadamente.

Figure 17 pone todo el capítulo en una sola página: diagnostica la redundancia, deja que tu objetivo decida el remedio, elige un modelo por una puntuación con mentalidad de predicción, y confirma al ganador con datos que nunca vio.

Resultados clave de un vistazo

Términos clave. Multicolinealidad, factor de inflación de la varianza, regresión auxiliar, $R^2$ ajustado, $C_p$ de Mallows, AIC, BIC, PRESS, residuo eliminado, selección paso a paso, sobreajuste, inferencia posterior a la selección, conjuntos de entrenamiento y de prueba, validación cruzada de $k$ particiones, raíz del error cuadrático medio, equilibrio sesgo-varianza, contracción, regresión ridge, lasso, parámetro de ajuste $\lambda$.

Ahora deberías ser capaz de

• Diagnosticar la multicolinealidad con el VIF y deducirlo de una regresión auxiliar (Teorema 12.4).

• Explicar qué daña y qué no daña la multicolinealidad, y enumerar los remedios estándar.

• Calcular y comparar $R^2$ ajustado, $C_p$ de Mallows, AIC, BIC y PRESS.

• Deducir el $C_p$ de Mallows (Teorema 12.5) y demostrar el atajo del residuo eliminado de PRESS (Teorema 12.6).

• Explicar, sin rodeos, por qué la selección paso a paso y la inferencia posterior a la selección engañan.

• Implementar validación cruzada de $k$ particiones a mano en R y Python.

• Describir la contracción de ridge y lasso y demostrarlas en los datos de grasa corporal.

Dónde encaja esto. Este capítulo son las etapas AJUSTAR y COMPROBAR de El flujo de trabajo del modelado hechas juntas y repetidamente: ajustas muchos modelos, compruebas cada uno por colinealidad y error de predicción, e iteras hacia uno que puedas defender. Los capítulos 8 a 11 te dieron una sola regresión múltiple y las herramientas para interpretarla y diagnosticarla; este capítulo enfrentó la situación más difícil de muchos modelos competidores y te enseñó a elegir entre ellos por la predicción en vez del ajuste. La mentalidad de validación pasa directo a 13. Regresión logística, donde el mismo pensamiento de datos apartados se vuelve las métricas de clasificación (exactitud, ROC) que juzgan un modelo logístico, y las ideas de contracción y de criterios reaparecen dondequiera que un modelo tenga más parámetros de los que los datos pueden sostener cómodamente.

12.7 Preguntas frecuentes¶

P1. ¿La multicolinealidad sesga mis coeficientes? No. Los mínimos cuadrados siguen siendo insesgados bajo colinealidad: en promedio los coeficientes siguen siendo correctos. Lo que la colinealidad hace es inflar su varianza, así que las estimaciones de cualquier muestra individual pueden estar tremendamente desviadas, con signos cambiados y errores estándar enormes. El promedio está bien; la estimación individual no es confiable.

P2. Si la colinealidad no perjudica la predicción, ¿por qué preocuparse por ella? Depende de tu objetivo. Si solo quieres predecir la grasa corporal para hombres como los de los datos, un modelo colineal predice bien y puedes ignorar los VIF. Si quieres interpretar los coeficientes (para decir cuál medida importa y por cuánto), la colinealidad vuelve esos coeficientes sin sentido, y debes arreglarla quitando, combinando o encogiendo predictores.

P3. ¿Qué criterio de selección debo usar? Para interpretación, BIC y su penalización más pesada te dan un modelo austero y defendible. Para predicción, el error validado cruzadamente es la respuesta más directa, y el $C_p$ de Mallows es un sustituto rápido y bien motivado. No agonices por diferencias pequeñas: informa que los criterios apuntaron a un rango de modelos (aquí, de cuatro a ocho predictores) y elige dentro de ese rango por una razón que puedas enunciar.

P4. ¿Por qué el $C_p$ del modelo completo siempre es igual al número de sus parámetros? Por construcción. $C_p = \mathrm{SSE}_p/\hat{\sigma}^2 - (n - 2p)$, y para el modelo completo $\hat{\sigma}^2 = \mathrm{SSE}_{\text{full}}/(n - p_{\text{full}})$, así que $\mathrm{SSE}_{\text{full}}/\hat{\sigma}^2 = n - p_{\text{full}}$ y $C_p = (n - p_{\text{full}})

• (n - 2p_{\text{full}}) = p_{\text{full}}$. El modelo completo es el punto de referencia, así que su$C_p$aterriza exactamente sobre la línea$C_p = p$.

P5. ¿Un PRESS más bajo siempre es mejor? Como comparación de modelos sobre los mismos datos, sí: un PRESS más bajo significa mejor predicción dejando uno fuera. Pero PRESS sigue siendo una estimación dentro de la muestra en el sentido de que usa los mismos 252 hombres para construir y para probar (solo que nunca el mismo hombre para ambos a la vez), así que puede ser optimista si todo el conjunto de datos es inusual. Un conjunto de prueba realmente apartado, o validación cruzada de $k$ particiones repetida con semillas distintas, es una comprobación más fuerte.

P6. ¿Puedo simplemente confiar en el modelo que devuelve la función step() de mi software? No. La selección paso a paso sobreajusta y sus valores $p$ reportados son inválidos (Sección 12.3). Si la usas, trata su salida como un candidato a validar con datos apartados, nunca como un modelo terminado, y nunca cites la significancia de sus coeficientes como si el modelo se hubiera elegido de antemano.

P7. El lasso y el paso a paso ambos quitan predictores. ¿Cuál es la diferencia? El paso a paso toma decisiones duras e inestables de dentro o fuera por valor $p$ y reporta inferencia deshonesta. El lasso toma la decisión de forma continua a través de una sola penalización $\lambda$ elegida por validación cruzada, encogiendo los coeficientes suavemente y apagándolos solo a medida que crece la penalización. El lasso es más estable, tiene una justificación coherente de sesgo-varianza, y está construido para predecir. En caso de duda, prefiérelo.

P8. ¿Por qué estandarizar los predictores antes de ridge o lasso? La penalización $\sum \beta_k^2$ o $\sum |\beta_k|$ suma coeficientes medidos en las unidades propias de cada predictor. Sin estandarizar, un predictor medido en unidades pequeñas (y por eso con un coeficiente grande) sería penalizado más que un predictor idéntico en unidades grandes. Estandarizar pone cada predictor en la misma escala para que la penalización sea justa. Tanto glmnet como el ajuste escalado de scikit-learn aquí hacen esto.

12.8 Problemas de práctica¶

1. (A) En una sola oración cada uno, di qué le hace la multicolinealidad al sesgo de $b_k$, a la varianza de $b_k$, y a las predicciones $\hat{Y}$ hechas dentro del rango de los datos. Luego nombra dos de los remedios estándar para un predictor con VIF alto.

2. (A) Un predictor tiene $\mathrm{VIF} = 5$. ¿Cuál es el $R_k^2$ auxiliar, y cuántas veces más grande es $\operatorname{Var}\{b_k\}$ que lo que sería si el predictor no estuviera correlacionado con los demás?

3. (A) Explica por qué $R^2$ no puede usarse para elegir entre un modelo con 5 predictores y uno con 8 predictores, pero el $R^2$ ajustado sí.

4. (A) Enuncia la regla práctica que conecta el $C_p$ de un “buen” modelo con su número de parámetros $p$, y di qué indica un $C_p$ muy por encima de $p$.

5. (A) ¿Por qué BIC tiende a seleccionar modelos más pequeños que AIC? Señala el término específico que difiere.

6. (A) En palabras sencillas, ¿qué mide el estadístico PRESS que el $\mathrm{SSE}$ ordinario no?

7. (A) Un estudiante selecciona predictores por regresión paso a paso, luego reporta los valores $p$ del modelo como evidencia de que es correcto. Nombra la falacia y explícala en dos oraciones.

8. (A) Describe la diferencia entre lo que ridge y lasso le hacen a un coeficiente, y di cuál de los dos realiza selección de variables.

9. (A) La RMSE de entrenamiento de un modelo es 3.5 y su RMSE de validación cruzada de 10 particiones es 4.6. ¿Qué te dice la brecha, y cuál número estima el desempeño con datos nuevos?

10. (A) Explica por qué los predictores se estandarizan antes de aplicar una penalización ridge o lasso.

11. (B) Deduce $\operatorname{Var}\{b_k\} = \sigma^2/(S_{kk}(1 - R_k^2))$ (Teorema 12.4) partiendo del hecho de variable agregada de que $b_k$ es la pendiente de $Y$ sobre los residuos de $X_k$ regresado sobre los demás predictores. Identifica dónde aparece $\mathrm{VIF}_k$.

12. (B) Muestra que la suma de cuadrados del error de la regresión auxiliar es igual a $S_{kk}(1 - R_k^2)$, y úsalo para explicar por qué un predictor casi determinado por los demás tiene un denominador casi nulo en la varianza de su coeficiente.

13. (B) Deduce el $C_p$ de Mallows (Teorema 12.5) a partir del error cuadrático medio total estandarizado $\Gamma_p$, incluyendo los pasos $\sum_i \operatorname{Var}\{\hat{Y}_i\} = p\sigma^2$ y $E\{\mathrm{SSE}_p\} = B + (n - p)\sigma^2$. Concluye que un modelo sin sesgo tiene $C_p \approx p$.

14. (B) Demuestra el atajo del residuo eliminado de PRESS $Y_i - \hat{Y}_{(i)} = e_i/(1 - h_{ii})$ (Teorema 12.6), enunciando con claridad dónde se usa la actualización del inverso de Sherman-Morrison.

15. (B) Muestra que para el modelo completo, $C_p$ es igual a su número de parámetros exactamente. (Parte de la definición de $C_p$ y del hecho de que $\hat{\sigma}^2$ es el MSE del modelo completo.)

16. (B) Un caso de alto apalancamiento tiene $h_{ii} = 0.9$ y residuo ordinario $e_i = 2$. Calcula su residuo eliminado, y explica por qué la predicción dejando uno fuera magnifica los residuos de los casos de alto apalancamiento.

17. (B) Explica, usando la descomposición sesgo-varianza del error de predicción, por qué un modelo puede tener a la vez menor error de entrenamiento y mayor error de prueba que un modelo más pequeño. ¿Cuál término de la descomposición sube a medida que el modelo crece?

18. (B) Muestra que a medida que $\lambda \to \infty$ el coeficiente ridge de un solo predictor estandarizado $\hat{\beta}_{\text{ridge}} = \sum x_i Y_i / (\sum x_i^2 + \lambda)$ tiende a cero, y que en $\lambda = 0$ es igual a la pendiente de mínimos cuadrados. (Usa la solución ridge de un predictor, que puedes deducir en el Problema 19.)

19. (B) Deduce el estimador ridge de un predictor del Problema 18 igualando a cero la derivada de $\sum (Y_i - \beta x_i)^2 + \lambda \beta^2$, e interpreta el denominador como “señal más penalización”.

20. (C) Ajusta el modelo completo y reproduce los VIF. Identifica cada predictor con $\mathrm{VIF}

    10$ y, para uno de ellos, confirma el VIF ejecutando su regresión auxiliar a mano.

21. (C) Ejecuta la búsqueda exhaustiva del mejor subconjunto (R leaps, o el bucle de Python) e informa el mejor modelo de cada tamaño con su $C_p$. Enuncia los tamaños elegidos por $C_p$, $R^2$ ajustado y BIC.

22. (C) Para el mejor subconjunto de cuatro predictores (peso, abdomen, antebrazo, muñeca), calcula $C_p$, AIC, BIC y PRESS. Compara cada uno con el candidato de ocho predictores y di cuál modelo reportarías y por qué.

23. (C) Escribe una función que calcule PRESS a partir de un modelo ajustado usando el atajo, y verifícala contra un bucle de fuerza bruta dejando uno fuera sobre el modelo candidato. Informa ambos valores.

24. (C) Reproduce la demostración de puro ruido posterior a la selección con una semilla distinta. Informa el $R^2$ del modelo seleccionado y el valor $p$ de la prueba $F$ general, y explica por qué son engañosos.

25. (C) Realiza una división entrenamiento/prueba 70/30 con semilla. Ajusta los modelos completo y candidato sobre el conjunto de entrenamiento e informa ambas RMSE del conjunto de prueba. ¿Cuál modelo predice mejor, y coincide el orden con el ajuste del conjunto de entrenamiento?

26. (C) Escribe un bucle de validación cruzada de 5 particiones a mano y compara el modelo completo, el candidato de ocho predictores, y el modelo de cuatro predictores (peso, abdomen, antebrazo, muñeca) por CV RMSE. Ordénalos.

27. (C) Ajusta un lasso con cv.glmnet (R) o LassoCV (Python) a todo el conjunto de datos. Informa cuáles predictores tienen coeficientes distintos de cero en la penalización validada cruzadamente, y compara ese conjunto con el modelo $C_p$ de ocho predictores.

28. (C) Ajusta la regresión ridge a través de una malla de valores de $\lambda$ y grafica el coeficiente de weight contra $\log\lambda$. Describe cómo cambia el coeficiente, y conecta su comportamiento con el VIF de 33.5 que encontraste para weight.

29. (C) Usando el modelo candidato, grafica los residuos eliminados $e_i/(1 - h_{ii})$ contra los residuos ordinarios $e_i$. ¿Cuáles casos se mueven más, y qué propiedad de esos casos lo explica?

30. (C) Toma el 25 por ciento de los hombres con los mayores apalancamientos y reajusta el modelo candidato sin ellos. Informa cuánto cambian los coeficientes y PRESS, y discute si el modelo se está apoyando en unos pocos casos.

12.9 Práctica de examen¶

Estas cinco preguntas están escritas al estilo de los exámenes del curso: cada una te pide explicar tu razonamiento en oraciones completas, no solo producir un número. Donde se muestra salida, se produjo en la máquina del curso a partir de fat.csv; lee los números que necesites del impreso y di cuáles usaste. Trabaja cada una antes de abrir su respuesta modelo.

EP 12.1 (conceptos, evaluar una afirmación). El modelo completo de grasa corporal del Ejemplo 12.1 fue reajustado y se imprimieron el coeficiente de weight y tres factores de inflación de la varianza.

full <- lm(brozek ~ age + weight + height + neck + chest + abdom + hip +
             thigh + knee + ankle + biceps + forearm + wrist, data = fat)
round(coef(summary(full))["weight", ], 4)
round(car::vif(full)[c("weight", "hip", "abdom")], 2)

  Estimate Std. Error    t value   Pr(>|t|)
   -0.0803     0.0496    -1.6198     0.1066
weight    hip  abdom
 33.51  14.80  11.77

Un estudiante lee el coeficiente negativo de weight y concluye: “manteniendo fijas las otras doce medidas, los hombres más pesados tienen menos grasa corporal, así que deberíamos decirles a los clientes que ganar peso baja su porcentaje de grasa corporal”. Evalúa esta afirmación. Usa el factor de inflación de la varianza en tu respuesta, y enuncia con claridad qué puede y qué no puede sostener el coeficiente.

EP 12.2 (interpretar salida en contexto). Se ejecutó una búsqueda exhaustiva del mejor subconjunto sobre una versión de seis predictores del modelo, brozek ~ age + weight + height + abdom + hip + thigh, informando el mejor subconjunto de cada tamaño.

rs <- regsubsets(brozek ~ age + weight + height + abdom + hip + thigh,
                 data = fat, nvmax = 6)
ss <- summary(rs)
rbind(Cp = round(ss$cp, 2), adjR2 = round(ss$adjr2, 4), BIC = round(ss$bic, 1))
c(min_Cp = which.min(ss$cp), max_adjR2 = which.max(ss$adjr2), min_BIC = which.min(ss$bic))

        [,1]    [,2]    [,3]    [,4]    [,5]    [,6]
Cp     55.07    6.29    3.88    4.40    5.13    7.00
adjR2 0.6608  0.7165  0.7203  0.7208  0.7212  0.7202
BIC  -262.40 -303.10 -302.00 -298.00 -293.70 -288.30
   min_Cp max_adjR2   min_BIC
        3         5         2

El mejor modelo de tamaño 2 es weight + abdom, el mejor modelo de tamaño 3 es weight + abdom + thigh, y el mejor modelo de tamaño 5 es weight + height + abdom + hip + thigh. Interpreta esta salida. Enuncia cuál tamaño de modelo selecciona cada uno de los tres criterios, explica en oraciones completas por qué el $C_p$ de Mallows y BIC aterrizan en tamaños distintos, y di cuál modelo reportarías a un colega que quiere interpretar los coeficientes, con tu razón.

EP 12.3 (qué cambiaría si). Para el modelo de cuatro predictores brozek ~ weight + abdom + forearm + wrist, la suma de cuadrados del error es $\mathrm{SSE}_p = 3994.31$ con $n = 252$ y $p = 5$ parámetros. La estimación de la varianza del ruido del modelo completo es $\hat{\sigma}^2 = 15.90$.

cand4 <- lm(brozek ~ weight + abdom + forearm + wrist, data = fat)
c(SSE = sum(resid(cand4)^2), sigma2_full = summary(full)$sigma^2)

       SSE sigma2_full
  3994.311      15.904

Primero calcula $C_p$ para este modelo. Luego responde: ¿qué cambiaría si estimaras $\sigma^2$ no del modelo completo sino del modelo de solo abdomen, cuyo cuadrado medio residual es 20.38? Calcula el $C_p$ resultante y explica, en oraciones completas, por qué $C_p$ se vuelve sin sentido con esa elección.

EP 12.4 (evaluar una afirmación, interpretar salida). Tres modelos se compararon por su RMSE de entrenamiento y por la RMSE validada cruzadamente de 5 particiones sobre todos los datos de grasa corporal: el modelo completo de trece predictores, el candidato de ocho predictores (edad, peso, cuello, abdomen, cadera, muslo, antebrazo, muñeca), y un modelo de dos predictores (abdomen, peso).

set.seed(4210)
folds <- sample(rep(1:5, length.out = nrow(fat)))
# training RMSE and 5-fold CV RMSE for each model

full   train RMSE 3.876   5-fold CV RMSE 4.097
cand8  train RMSE 3.893   5-fold CV RMSE 4.003
two    train RMSE 4.103   5-fold CV RMSE 4.183

Un estudiante argumenta: “el modelo completo tiene la RMSE de entrenamiento más baja (3.876), así que predice mejor; reporta el modelo completo”. Evalúa este argumento. Di cuál número estima cómo le irá al modelo con un hombre nuevo y por qué, explica por qué el orden se invierte entre las dos columnas, y enuncia cuál modelo reportarías.

EP 12.5 (qué cambiaría si). Se ajustó un lasso a todos los datos de grasa corporal con la penalización $\lambda$ elegida por validación cruzada, y los coeficientes se leyeron en dos penalizaciones, lambda.min (el minimizador del error) y el más disperso lambda.1se.

X <- as.matrix(fat[, preds]); Y <- fat$brozek
set.seed(4210)
cvl <- cv.glmnet(X, Y, alpha = 1)
round(cbind(min = coef(cvl, s = "lambda.min")[, 1],
            onese = coef(cvl, s = "lambda.1se")[, 1]), 3)

                min   onese
(Intercept) -13.683  -6.733
age           0.055   0.052
weight       -0.065   0.000
height       -0.081  -0.149
neck         -0.403  -0.137
chest         0.000   0.000
abdom         0.830   0.635
hip          -0.147   0.000
thigh         0.171   0.000
knee          0.000   0.000
ankle         0.086   0.000
biceps        0.114   0.000
forearm       0.387   0.112
wrist        -1.437  -1.265

En lambda.min once predictores tienen coeficientes distintos de cero; en lambda.1se solo seis. Explica qué le pasa a los coeficientes a medida que $\lambda$ aumenta de lambda.min a lambda.1se, por qué el lasso fija algunos coeficientes exactamente en cero mientras que la regresión ridge no lo haría, y qué cambiaría si los predictores no se estandarizaran antes de ajustar.

Juego del capítulo¶