9. Diagnósticos del modelo (en espanol)

9. Diagnósticos del modelo¶

En 1973 dos biólogos, Michael Johnson y Peter Raven, contaron el número de especies de plantas en 30 islas del archipiélago de las Galápagos y registraron unos cuantos datos geográficos de cada isla: su área, su elevación máxima, qué tan lejos está de la isla más cercana, qué tan lejos de la isla central de Santa Cruz, y el área de la isla vecina más próxima. La pregunta natural es si la geografía predice la riqueza biológica. Las islas más grandes, más altas y más aisladas podrían albergar más especies. Una regresión múltiple del conteo de especies sobre esos cinco predictores es el primer paso obvio, y parece un éxito: explica alrededor del 77 por ciento de la variación en los conteos de especies.

Luego graficas los residuos, y el modelo se desmorona. La Figure 1 muestra los residuos contra los conteos ajustados. Tres problemas saltan a la vista de inmediato. La dispersión de los residuos crece a medida que crece el conteo ajustado, un embudo que rompe el supuesto de varianza constante. Una isla, Isabela, queda lejos del resto. Y para varias islas pequeñas el modelo predice un número negativo de especies, algo imposible. Un ajuste que reporta un $R^2$ alto puede aun así estar equivocado en todo lo que importa.

Este capítulo trata la cuarta etapa del flujo de trabajo de modelado, CHECK (verificar). Has aprendido a PREGUNTAR (ASK) una pregunta, EXPLORAR (EXPLORE) los datos, y AJUSTAR (FIT) un modelo en forma matricial (7.1 El modelo y los mínimos cuadrados en forma matricial) con toda su maquinaria de inferencia (8.4 La prueba lineal general). Ajustar es ya la parte fácil; cualquier computadora lo hará, y cualquier estudiante también. La destreza escasa es el criterio: decidir si un modelo ajustado merece confianza, encontrar las observaciones que en silencio llevan la voz cantante, y nombrar con exactitud qué supuesto rompe un modelo. Eso es lo que hacen los diagnósticos. Al final del capítulo tendrás un conjunto de herramientas que convierte “el ajuste se ve bien” en una lista de comprobación que puedes defender.

Trabajamos dos conjuntos de datos reales en paralelo. Los datos de especies de las Galápagos (gala) son el modelo que falla en todas las pruebas. Los datos de ahorro entre países (savings), que conociste por primera vez en el capítulo de correlación, son el modelo cuyos supuestos en su mayoría se cumplen pero cuyo ajuste está gobernado por un solo país. En el camino, los datos de producción de Toluca regresan una última vez para responder una pregunta que no pudimos zanjar antes: ¿es de verdad una recta la forma correcta?

Todo en este capítulo responde una de dos preguntas sencillas. ¿Algún punto de datos es problemático, ya sea porque ocupa un lugar extraño o porque en silencio dirige el ajuste? ¿O el modelo rompe una promesa que hizo sobre los errores tomados en conjunto? La Figure 2 es el mapa. La rama izquierda persigue los casos problemáticos (Secciones 9.1 a 9.3); la rama derecha verifica los supuestos (Secciones 9.4 a 9.6); cada casilla nombra la herramienta y la sección que la responde. Vuelve a esta imagen cada vez que pierdas el hilo.

9.1 Los valores de apalancamiento y la matriz sombrero¶

Intuición¶

Antes de juzgar un residuo tenemos que saber cuánta libertad le dio el modelo a cada punto. Imagina la recta ajustada como un balancín. Un punto que está cerca del centro de los valores del predictor puede subir o bajar bastante sin desviar mucho la recta; la multitud a su alrededor la mantiene en su sitio. Un punto muy al borde del rango del predictor es distinto: como es la única observación por allá, la recta tiene que seguirlo, así que el valor ajustado en ese punto es jalado hacia la propia respuesta del punto. Ese jalón lo mide el valor de apalancamiento $h_{ii}$ (Definición 9.1). Un punto de alto apalancamiento es aquel cuyos valores del predictor son inusuales, de modo que el ajuste se ve obligado a prestarle atención especial sin importar su respuesta.

El valor de apalancamiento es una afirmación sobre los predictores únicamente. Todavía no dice nada sobre si un punto está bien o mal ajustado. Mide potencial, la capacidad de influir en la recta, y dedicaremos las próximas dos secciones a convertir el potencial en un veredicto.

Fórmula¶

Recuerda la matriz sombrero de 7.3 La matriz sombrero. Apilando las observaciones en la matriz de diseño $\mathbf{X}$ (una matriz $n \times p$ cuya primera columna es de unos, con $p$ el número de parámetros incluyendo el intercepto), los valores ajustados son $\hat{\mathbf{Y}} = \mathbf{H}\mathbf{Y}$ donde

• $\mathbf{H}$ es la matriz sombrero $n \times n$, la proyección ortogonal sobre el espacio de columnas de $\mathbf{X}$; “le pone el sombrero a $\mathbf{Y}$”.

• $\mathbf{X}'$ es la transpuesta de $\mathbf{X}$, y $(\mathbf{X}'\mathbf{X})^{-1}$ es la inversa construida en 6.4 La inversa y las ecuaciones normales.

Derivación (propiedades de los valores de apalancamiento)¶

Demostración. En 7.3 La matriz sombrero demostramos que $\mathbf{H}$ es simétrica ($\mathbf{H}' = \mathbf{H}$) e idempotente ($\mathbf{H}\mathbf{H} = \mathbf{H}$), con traza igual a $p$. Tres hechos se siguen directamente.

Primero, como $\mathbf{H}$ es simétrica e idempotente, $\mathbf{H} = \mathbf{H}\mathbf{H}$ leída en la diagonal da

Reordenando, $h_{ii} - h_{ii}^2 = \sum_{j \neq i} h_{ij}^2 \ge 0$, así que $h_{ii}(1 - h_{ii}) \ge 0$. Como también $h_{ii} = \sum_j h_{ij}^2 \ge 0$, obtenemos

En palabras: un valor de apalancamiento queda atrapado entre 0 y 1, y un valor cercano a 1 significa que el ajuste en ese punto está determinado casi por completo por la propia respuesta de ese punto.

Segundo, los valores de apalancamiento suman la traza de $\mathbf{H}$:

Tercero, dividiendo entre $n$, el valor de apalancamiento promedio es $p/n$. Un punto no lleva su parte justa del espacio de predictores si su valor de apalancamiento está muy por encima de ese promedio, lo que motiva la regla práctica común: examina cualquier caso con

El factor 2 es una convención, no una ley; marca aproximadamente los puntos que merecen una segunda mirada. $\blacksquare$

Para la regresión simple el valor de apalancamiento tiene una forma cerrada limpia que muestra con claridad lo que mide.

Demostración. Con un predictor, $\mathbf{x}_i' = (1, X_i)$ y, como se calculó en 7.3 La matriz sombrero,

Desarrollando la forma cuadrática $\mathbf{x}_i'(\mathbf{X}'\mathbf{X})^{-1}\mathbf{x}_i$ y simplificando con $\sum X_i = n\bar X$ se obtiene

En palabras: el valor de apalancamiento en la regresión simple es una base $1/n$ compartida por todos, más una penalización que crece con la distancia al cuadrado de $X_i$ respecto al valor medio del predictor. La Figure 3 dibuja esta parábola sobre los tamaños de lote de Toluca. $\blacksquare$

R y Python¶

Un gráfico de índice deja obvias las dos que sobresalen. La Figure 4 grafica el valor de apalancamiento de cada isla contra su posición en los datos, con la línea $2p/n$ trazada: Isabela y Fernandina se elevan por encima de las otras 28 islas, y todo lo demás cae en una banda inofensiva.

La forma cerrada $h_{ii} = 1/n + (X_i - \bar X)^2/S_{xx}$ vale la pena comprobarla contra el software una vez, para que confíes en que hatvalues y la matriz sombrero coinciden.

toluca <- read.csv("data/toluca.csv")
tfit <- lm(hours ~ lotsize, data = toluca)
x <- toluca$lotsize
h_closed <- 1 / nrow(toluca) + (x - mean(x))^2 / sum((x - mean(x))^2)
max(abs(hatvalues(tfit) - h_closed))

[1] 2.775558e-17

toluca = pd.read_csv("data/toluca.csv")
tfit = smf.ols("hours ~ lotsize", data=toluca).fit()
x = toluca["lotsize"].to_numpy()
h_closed = 1 / len(toluca) + (x - x.mean()) ** 2 / np.sum((x - x.mean()) ** 2)
print(np.max(np.abs(tfit.get_influence().hat_matrix_diag - h_closed)))

4.163336342344337e-17

Ambas coinciden hasta el polvo de la aritmética de punto flotante, lo que confirma la fórmula.

El Teorema 9.2 dice que los valores de apalancamiento son un presupuesto fijo que siempre suma $p$, y esa afirmación resulta mucho más creíble después de ver cómo se gasta ese presupuesto.

Qué observar. Ningún valor de la respuesta entra en esta imagen, porque $h_{ii}$ se construye solo con los predictores. Prueba esto. Empuja la corrida nueva hasta un tamaño de lote de 200 y confirma que su propio apalancamiento sube por encima de 0.47 mientras el de cada corrida original baja, con la suma fijada en $p = 2$ todo el tiempo (9.1 Los valores de apalancamiento y la matriz sombrero).

9.2 Cuatro variantes del residuo¶

Intuición¶

El valor de apalancamiento midió qué puntos tenían el poder de distorsionar el ajuste. La siguiente pregunta es a qué puntos el ajuste de hecho no acierta. Para eso leemos el residuo $e_i = Y_i - \hat Y_i$, la sobra del modelo: la parte de $Y_i$ que el ajuste no capturó. Es tentador cazar valores atípicos revisando los residuos crudos en busca de valores grandes. Es una trampa, porque los residuos crudos no están en igualdad de condiciones. Un punto de alto apalancamiento jala la recta hacia sí mismo, lo que encoge su propio residuo. Así que los mismos puntos con mayor capacidad de distorsionar el ajuste tienden a tener residuos crudos engañosamente pequeños. Para comparar residuos de forma justa tenemos que ponerlos en una escala común, y eso significa dividir entre la desviación estándar correcta, que depende del valor de apalancamiento.

Fórmula¶

El factor de escala $\sqrt{1 - h_{ii}}$ no es adorno. Proviene de una fórmula exacta de la varianza de los residuos.

Derivación (la varianza de un residuo)¶

Demostración. Escribe el vector de residuos como una transformación lineal de las respuestas. Como $\hat{\mathbf{Y}} = \mathbf{H}\mathbf{Y}$,

Bajo el modelo $\operatorname{Var}\{\mathbf{Y}\} = \sigma^2 \mathbf{I}$ (los errores no están correlacionados y tienen varianza común, de 6.7 Vectores aleatorios, esperanza y matrices de covarianzas). La matriz $\mathbf{I} - \mathbf{H}$ es simétrica e idempotente, igual que $\mathbf{H}$. Aplicando la regla de covarianza para una transformación lineal,

Leyendo la $i$-ésima entrada diagonal,

En palabras: el residuo en un punto de alto apalancamiento tiene menor varianza, porque el ajuste se dobla para alcanzarlo. Por eso justamente un residuo crudo subestima el problema en los puntos de alto apalancamiento, y por eso dividir entre $\sqrt{1 - h_{ii}}$ restablece una escala común. La Figure 6 muestra el encogimiento y su corrección. $\blacksquare$

Derivación (el residuo eliminado studentizado y la prueba de valor atípico)¶

El residuo estandarizado $r_i$ tiene una falla sutil: si el caso $i$ es de verdad un valor atípico, infla $s$, la misma cantidad que está en su propio denominador, y se enmascara a sí mismo. El arreglo es estimar la desviación estándar del error a partir de los otros casos. Sea $s_{(i)}$ la raíz del error cuadrático medio del ajuste con el caso $i$ eliminado. Una identidad algebraica útil, que evita reajustar $n$ veces, es

En palabras: esto construye la varianza del error de dejar uno fuera a partir de cantidades que el ajuste completo ya te da, restando la propia contribución de error al cuadrado del caso $i$, de modo que nunca reajustas de verdad.

Demostración. Considera el residuo eliminado $d_i = Y_i - \hat Y_{i(i)}$, la brecha entre $Y_i$ y la predicción en $\mathbf{x}_i$ de un modelo que nunca vio el caso $i$. Como esa predicción no usa $Y_i$, los dos son independientes. Una identidad que probamos a partir de la actualización de dejar uno fuera en 9.3 Influencia: qué puntos cambian realmente el ajuste da $d_i = e_i / (1 - h_{ii})$, con varianza

Estimando $\sigma^2$ con el $s_{(i)}^2$ independiente y dividiendo $d_i$ entre su error estándar estimado se obtiene

bajo el modelo, cuando el caso $i$ no es un valor atípico. Como el numerador y el denominador son independientes, este es un cociente $t$ exacto. $\blacksquare$

Para probar si el caso más extremo es un valor atípico genuino miras $\max_i |t_i|$. Como implícitamente estás realizando $n$ pruebas, usa un punto de corte de Bonferroni: marca el caso $i$ cuando $|t_i| > t(1 - \alpha/(2n);\, n - p - 1)$. La corrección mantiene la probabilidad de una sola falsa alarma entre todos los $n$ casos por debajo de $\alpha$.

9.3 Influencia: qué puntos cambian realmente el ajuste¶

Intuición¶

Un valor de apalancamiento alto es potencial y un residuo grande es un síntoma; la influencia es lo que de verdad nos importa, la cantidad en que cambiaría el ajuste si un punto no estuviera (Definición 9.7). Un punto puede ser de alto apalancamiento y aun así tener poca influencia, si resulta que cae justo sobre la tendencia que marcan los demás puntos. Y un punto con un residuo grande pero un valor de apalancamiento ordinario puede jalar la recta solo un poco. La influencia es el producto de los dos ingredientes: una fila de predictores inusual y una respuesta que el modelo ajusta mal.

La Figure 7 lo ilustra con cuatro pequeñas imágenes. En cada una dejamos caer un solo punto rojo en la misma nube y reajustamos. Empuja el punto muy lejos en $X$ pero mantenlo sobre la tendencia, y la recta no se mueve: alto apalancamiento, ninguna influencia. Dale un error vertical grande pero un $X$ ordinario, y la recta apenas se inclina: un valor atípico simple. Solo cuando el punto está a la vez lejos en $X$ y fuera de la tendencia la recta se balancea para perseguirlo. Ese último panel es influencia, y requiere ambos ingredientes al mismo tiempo.

Fórmula¶

La segunda forma hace explícita la receta: la influencia es igual a una parte de residuo ($r_i^2$) por un factor $h_{ii}/(1 - h_{ii})$ que crece con el valor de apalancamiento. Ambos deben ser considerables para que $D_i$ sea grande. Una guía aproximada trata $D_i$ cercano o por encima de 1 como digno de atención seria, y cualquier caso muy por encima del resto como digno de nombrar.

Derivación (la actualización de dejar uno fuera y la distancia de Cook)¶

Cada medida de influencia descansa sobre un hecho algebraico: puedes calcular el ajuste sin el caso $i$ a partir del ajuste completo, sin necesidad de reajustar.

Demostración. Eliminar el caso $i$ quita la fila $\mathbf{x}_i'$ de $\mathbf{X}$ y el valor $Y_i$ de $\mathbf{Y}$, así que

La fórmula de Sherman-Morrison da la inversa de una actualización de rango uno por resta: para una $\mathbf{A}$ invertible y un vector $\mathbf{u}$ con $\mathbf{u}'\mathbf{A}^{-1}\mathbf{u} \neq 1$,

que puedes verificar multiplicando el lado derecho por $\mathbf{A} - \mathbf{u}\mathbf{u}'$ y agrupando términos para obtener $\mathbf{I}$. Fija $\mathbf{A} = \mathbf{X}'\mathbf{X}$ y $\mathbf{u} = \mathbf{x}_i$, de modo que $\mathbf{u}'\mathbf{A}^{-1}\mathbf{u} = h_{ii}$. Multiplicando la inversa actualizada por $\mathbf{X}_{(i)}'\mathbf{Y}_{(i)}$ y simplificando (usando $\mathbf{x}_i'\mathbf{b} = \hat Y_i$ y $e_i = Y_i - \hat Y_i$) se colapsa al resultado compacto

En palabras: descartar el caso $i$ desplaza el vector de coeficientes en un múltiplo de $(\mathbf{X}'\mathbf{X})^{-1}\mathbf{x}_i$, escalado por el residuo del caso e inflado por $1/(1 - h_{ii})$ cuando el valor de apalancamiento es alto. Evaluando el efecto del desplazamiento en $\mathbf{x}_i$ mismo se obtiene $\hat Y_i - \hat Y_{i(i)} = h_{ii} e_i/(1 - h_{ii})$, así que el residuo eliminado es $d_i = Y_i - \hat Y_{i(i)} = e_i + h_{ii}e_i/(1-h_{ii}) = e_i/(1 - h_{ii})$, la identidad prometida en 9.2 Cuatro variantes del residuo. $\blacksquare$

Demostración. El numerador de $D_i$ es una forma cuadrática en ese desplazamiento, porque $\hat{\mathbf{Y}} - \hat{\mathbf{Y}}_{(i)} = \mathbf{X}(\mathbf{b} - \mathbf{b}_{(i)})$:

Sustituye la fórmula de eliminación. El $\mathbf{X}'\mathbf{X}$ del medio cancela una inversa a cada lado, dejando $\mathbf{x}_i'(\mathbf{X}'\mathbf{X})^{-1}\mathbf{x}_i = h_{ii}$:

Divide entre $p\,\mathrm{MSE}$ y reconoce $r_i^2 = e_i^2/(\mathrm{MSE}(1 - h_{ii}))$:

Esta es la separación prometida de $D_i$ en un residuo estandarizado al cuadrado y un término en el valor de apalancamiento $h_{ii}$. $\blacksquare$

R y Python¶

El gráfico de influencia de la Figure 9 combina toda la historia: el valor de apalancamiento en un eje, el residuo eliminado studentizado en el otro, y la distancia de Cook como tamaño de la burbuja. Separa los papeles que un resumen de un solo número mezcla.

Nombrar un punto influyente no es lo mismo que eliminarlo. El flujo de trabajo honesto es reajustar sin el caso y reportar cuánto se mueve la respuesta, para que un lector pueda juzgar. La Figure 10 y los números detrás de ella hacen exactamente eso.

Los cuatro paneles de la Figure 7 son imágenes fijas de ese argumento. En el widget de abajo tú mueves el caso solitario y observas como responde $D_i$.

Qué observar. El caso solitario empieza con el mayor valor de apalancamiento de los datos y una distancia de Cook prácticamente nula, que es el panel superior derecho de la Figure 7. Prueba esto. Arrástralo hacia abajo hasta cerca de 6.5, luego reinicia y arrastra el caso en $X = 5$ esa misma distancia: el mismo error vertical, y casi nada de influencia (9.3 Influencia: qué puntos cambian realmente el ajuste).

9.4 Leer los gráficos de diagnóstico¶

Intuición¶

Los números formales son afilados, pero la forma más rápida de atrapar un modelo roto es mirarlo. Tres gráficos hacen la mayor parte del trabajo, y cada uno apunta a un supuesto. El gráfico de residuos contra valores ajustados verifica la forma de la función media y la constancia de la varianza. El gráfico cuantil-cuantil normal verifica la forma de la distribución del error. El gráfico escala-ubicación verifica de nuevo la varianza constante, de forma más sensible. Aprendiste a dibujar el primero de ellos allá en 2.3 Valores ajustados y las propiedades de los residuos; ahora puedes leer los tres como un conjunto.

La Figure 12 es una guía de campo del gráfico de residuos contra valores ajustados. Cuatro formas se repiten: una banda horizontal sin forma (sano), un embudo (la varianza crece con la media), una curvatura (la función media tiene la forma equivocada), y un punto solitario descarriado (un valor atípico). Entrena el ojo con estas y la mayoría de los diagnósticos se vuelven reconocimiento de patrones.

R y Python¶

Ya viste el gráfico de residuos contra valores ajustados de las Galápagos en la Figure 1: un embudo de libro de texto. El gráfico cuantil-cuantil normal de la Figure 13 cuenta la mitad de la historia sobre la normalidad. En R es una sola línea, qqnorm(residuals(gfit)); qqline(residuals(gfit)), y statsmodels de Python ofrece sm.qqplot(gfit.resid, line="s"). Si los errores fueran normales los puntos abrazarían la recta de referencia; aquí la cola superior se aparta e Isabela se hunde muy por debajo, el eco visual de su residuo studentizado de -5.33.

El gráfico escala-ubicación de la Figure 14 grafica $\sqrt{|r_i|}$ contra el valor ajustado. Una nube plana significa varianza constante; una tendencia creciente significa que la dispersión crece con la media. La nube de las Galápagos sube de manera constante, confirmando el embudo desde un segundo ángulo.

9.5 Pruebas formales para los supuestos¶

Intuición¶

Los gráficos persuaden; las pruebas ponen un número sobre la sospecha. Úsalos juntos, y apóyate en los gráficos cuando los dos discrepen, porque una prueba puede pasar por alto un patrón que un gráfico deja obvio, y una prueba puede marcar una arruga demasiado pequeña para importar. Hay una prueba para cada supuesto: varianza constante, errores normales e independencia.

Para la varianza constante, dos pruebas son estándar. La prueba de Breusch-Pagan regresa los residuos al cuadrado sobre los predictores y pregunta si esa regresión explica algo; si la dispersión depende de los predictores, lo hace. La prueba de Brown-Forsythe divide los datos en grupos, compara la dispersión de los residuos entre grupos usando medianas (de modo que resiste valores atípicos), y es buena elección cuando sospechas que la varianza cambia con un predictor en particular.

Para la normalidad, la prueba de Shapiro-Wilk compara los residuos ordenados con lo que producirían datos normales; un valor p pequeño dice que los residuos no son normales. Un pariente cercano es la prueba de correlación para la normalidad, que calcula la correlación entre los residuos ordenados y sus puntuaciones normales (los puntos del gráfico cuantil-cuantil) y rechaza la normalidad cuando esa correlación queda demasiado por debajo de 1.

Para la independencia en datos ordenados en el tiempo, la prueba de Durbin-Watson verifica si los errores consecutivos están correlacionados.

Fórmula¶

El estadístico de Durbin-Watson es

donde $\hat\rho$ es la autocorrelación de rezago 1 estimada de los residuos, la correlación entre cada residuo y el anterior a él en orden temporal.

• $D$ va de 0 a 4. $D \approx 2$ significa ausencia de autocorrelación ($\hat\rho \approx 0$).

• $D < 2$ señala autocorrelación positiva (los errores se mueven en rachas, $\hat\rho > 0$).

• $D > 2$ señala autocorrelación negativa (los errores alternan de signo).

En palabras: Durbin-Watson compara cuánto difieren los residuos vecinos con qué tan grandes son los residuos en conjunto; cuando los vecinos son parecidos, el numerador es pequeño y $D$ cae por debajo de 2. La distribución nula exacta depende de la matriz de diseño, así que las tablas clásicas dan dos cotas $d_L$ y $d_U$; el software moderno devuelve en su lugar un valor p exacto. Este diagnóstico es la puerta de entrada a 15.2 Por que el tiempo rompe los minimos cuadrados ordinarios, donde los errores ordenados en el tiempo son todo el tema y donde se construyen modelos que arreglan el problema en lugar de solo detectarlo.

R y Python¶

Cuando la normalidad falla, como pasa con gala, una respuesta honesta es dejar de apoyarse en la maquinaria de $t$ y $F$ de la teoría normal y cambiar a la inferencia por remuestreo de 5.4 El bootstrap para la regresión, que no supone una forma normal del error. El arreglo más profundo, cambiar el modelo para que sus errores se comporten, es el tema del Capítulo 10 y, para datos de conteo como los totales de especies, del Capítulo 14.

La prueba de Durbin-Watson necesita datos en un orden con sentido, así que no aplica de verdad a los países de savings, que no tienen una secuencia natural. Para verla funcionar simulamos dos series ordenadas en el tiempo: una con errores independientes, otra con errores que se acarrean de un paso al siguiente.

9.6 La prueba F de falta de ajuste¶

Intuición¶

Cada diagnóstico hasta ahora pregunta si los errores se comportan mal. La prueba de falta de ajuste hace una pregunta distinta: ¿es correcta del todo la forma de la función de regresión? Supón que la media verdadera de $Y$ se curva, pero ajustas una recta. Entonces, incluso con datos perfectos, la recta no acertará las medias verdaderas. El problema es que la SSE ordinaria no puede distinguir “el modelo está equivocado” de “los datos tienen ruido”, porque ambas cosas la inflan. La salida es la repetición. Si tienes varias observaciones en el mismo valor del predictor, la dispersión entre ellas mide el ruido puro, sin ayuda de ningún modelo. Eso da una vara de error irreducible contra la cual comparar los errores del modelo.

Los datos de Toluca están hechos para esto. Las 25 corridas usaron solo 11 tamaños de lote distintos, así que la mayoría de los tamaños de lote se repitió. La Figure 16 muestra la idea: en cada tamaño de lote las corridas se dispersan alrededor de su propia media de grupo (error puro), y las medias de grupo se dispersan alrededor de la recta ajustada (falta de ajuste).

Fórmula¶

Derivación (la falta de ajuste como prueba lineal general)¶

Demostración. Esta es la prueba lineal general de 8.4 La prueba lineal general con un par específico de modelos. El modelo completo pone una media libre $\mu_j$ en cada uno de los $c$ niveles distintos del predictor: $Y_{jk} = \mu_j + \varepsilon_{jk}$. Su valor ajustado por mínimos cuadrados en el nivel $j$ es la media de grupo $\bar Y_j$, así que su suma de cuadrados del error es exactamente $\mathrm{SSPE} = \sum_j \sum_k (Y_{jk} - \bar Y_j)^2$, con $n - c$ grados de libertad porque estima $c$ medias. El modelo reducido es la recta $Y_{jk} = \beta_0 + \beta_1 X_j + \varepsilon_{jk}$, cuya suma de cuadrados del error es la $\mathrm{SSE}$ ordinaria con $n - 2$ grados de libertad.

El estadístico de la prueba lineal general es

Que la diferencia $\mathrm{SSE} - \mathrm{SSPE}$ sea igual a $\sum_j n_j(\bar Y_j - \hat Y_j)^2$ se sigue de desarrollar $Y_{jk} - \hat Y_j = (Y_{jk} - \bar Y_j) + (\bar Y_j - \hat Y_j)$, elevar al cuadrado, y sumar: el término cruzado se anula porque $\sum_k (Y_{jk} - \bar Y_j) = 0$ dentro de cada grupo. Bajo el modelo reducido las dos medias de cuadrados estiman ambas $\sigma^2$, así que $F^\ast$ está cerca de 1; cuando las medias verdaderas se apartan de una recta, solo $\mathrm{MSLF}$ se infla, ya que $E\{\mathrm{MSPE}\} = \sigma^2$ siempre. $\blacksquare$

R y Python¶

Vale la pena desarmar con las propias manos la separación de la Figure 16, porque las dos sumas de cuadrados responden a movimientos completamente distintos.

Qué observar. Mover todo un grupo de repeticiones a la vez deja intacto el error puro y le carga el movimiento entero a la falta de ajuste. Prueba esto. Sube unas 6 unidades las tres corridas en $X = 6$ y observa como el valor p cae de 1.0000 a 0.017; luego sube solo una de ellas y observa como la prueba se queda callada otra vez (Figure 16).

9.7 Resumen del capítulo¶

Los diagnósticos convierten “el ajuste se ve bien” en una lista de comprobación defendible. Ahora puedes derivar los valores de apalancamiento a partir de la matriz sombrero, acotarlos entre 0 y 1, mostrar que suman $p$, y calcularlos a mano en la regresión simple; colocar un residuo en cuatro escalas y explicar por qué el residuo crudo engaña en los puntos de alto apalancamiento; probar si un caso es atípico con el residuo eliminado studentizado y un punto de corte de Bonferroni; derivar la distancia de Cook a partir del cambio en los coeficientes al dejar uno fuera y leerla, junto con DFFITS y DFBETAS, para encontrar los puntos que de verdad mueven el ajuste; leer los gráficos de residuos, cuantil-cuantil y escala-ubicación y aplicar las pruebas de Breusch-Pagan, Brown-Forsythe, Shapiro-Wilk y Durbin-Watson; y derivar y aplicar la prueba F de falta de ajuste con error puro. El ajuste de las Galápagos falló en todos los frentes mientras que el de savings estaba gobernado por un solo país, y la recta de Toluca quedó vindicada al fin.

Resultados clave de un vistazo

Términos clave. valor de apalancamiento, punto de alto apalancamiento, residuo crudo, residuo semiestudentizado, residuo estandarizado, residuo eliminado studentizado, residuo eliminado, valor atípico, observación influyente, distancia de Cook, DFFITS, DFBETAS, prueba de Breusch-Pagan, prueba de Brown-Forsythe, prueba de Shapiro-Wilk, prueba de correlación para la normalidad, prueba de Durbin-Watson, autocorrelación, error puro, falta de ajuste.

Ahora deberías poder

• Derivar las propiedades de los valores de apalancamiento a partir de la matriz sombrero y calcularlos a mano en la regresión simple.

• Distinguir los residuos crudo, semiestudentizado, estandarizado y eliminado studentizado, y explicar por qué el residuo crudo es una escala engañosa.

• Probar si un solo caso es atípico con el residuo eliminado studentizado y un punto de corte de Bonferroni.

• Derivar la distancia de Cook a partir del cambio en los coeficientes al dejar uno fuera, y calcular DFFITS y DFBETAS.

• Separar un punto de alto apalancamiento de uno influyente, y diagnosticar por nombre los ajustes de savings y gala.

• Leer los gráficos de residuos, cuantil-cuantil y escala-ubicación, y aplicar las pruebas de Breusch-Pagan, Brown-Forsythe, Shapiro-Wilk y Durbin-Watson.

• Derivar y aplicar la prueba F de falta de ajuste con error puro a partir de valores repetidos del predictor.

Dónde encaja esto. Este capítulo es la etapa CHECK del flujo de trabajo de El flujo de trabajo del modelado: después de AJUSTAR un modelo debes averiguar si merece confianza antes de USARLO para predecir o decidir. Los diagnósticos miran hacia atrás, a la maquinaria matricial de 7.3 La matriz sombrero, que construyó la matriz sombrero cuya diagonal contiene los valores de apalancamiento y cuyas propiedades alimentan cada medida de influencia, y a la prueba lineal general de 8.4 La prueba lineal general, que es la prueba de falta de ajuste disfrazada. Miran hacia adelante, al Capítulo 10, que toma los dos ajustes rotos de este capítulo, el embudo de las Galápagos y sus errores no normales, e intenta repararlos con transformaciones y mínimos cuadrados ponderados, continuando por nombre los casos de savings y gala. El hilo de gala corre aún más lejos: el Capítulo 14 da la resolución honesta con la regresión de Poisson, la herramienta correcta para una respuesta de conteo que nunca puede ser negativa. Un diagnóstico aquí solo se introduce, no se usa por completo: el estadístico de Durbin-Watson se vuelve todo el tema de 15.2 Por que el tiempo rompe los minimos cuadrados ordinarios, donde los errores ordenados en el tiempo se modelan y corrigen en lugar de solo detectarse.

9.8 Preguntas frecuentes¶

P1. ¿Es un punto de alto apalancamiento algo malo? No por sí mismo. Un valor de apalancamiento alto significa una fila de predictores inusual, que es potencial para influir en el ajuste, no prueba de daño. Un punto de alto apalancamiento que cae justo sobre la tendencia que marcan los demás apenas mueve la recta. Un valor de apalancamiento alto se vuelve un problema solo cuando se empareja con un residuo grande, y esa combinación es lo que mide la distancia de Cook.

P2. ¿Debería eliminar un punto influyente? Casi nunca de forma automática. Un punto influyente es una señal para investigar, no basura para descartar. Revísalo primero por un error de captura de datos. Si es real, el reporte honesto muestra el ajuste con y sin él y explica la diferencia, como hicimos con Libia. Eliminar datos reales para obtener una respuesta más limpia es como las análisis salen mal.

P3. ¿Cuál es la diferencia entre un valor atípico y un punto influyente? Un valor atípico tiene un residuo grande: el modelo lo ajusta mal. Un punto influyente cambia el ajuste cuando se elimina. Pueden ocurrir por separado. Zambia está cerca de ser un valor atípico con un valor de apalancamiento ordinario, así que se ajusta mal pero no dirige el modelo. Libia es influyente sin un residuo extremo, porque su valor de apalancamiento es muy alto.

P4. Mi gráfico de residuos se ve bien pero Shapiro-Wilk rechaza la normalidad. ¿Quién gana? Lee el tamaño de muestra. En muestras grandes las pruebas formales rechazan desviaciones diminutas e inofensivas de la normalidad; en muestras pequeñas pasan por alto las reales (el valor p de Breusch-Pagan de las Galápagos de 0.08 es un ejemplo de ese fallo). Usa el gráfico para juzgar si la desviación es lo bastante grande como para importar, y recuerda que las pruebas $t$ y $F$ son bastante tolerantes con la no normalidad leve.

P5. ¿Por qué dividir entre $s_{(i)}$ en lugar de $s$ para la prueba de valor atípico? Si el caso $i$ es de verdad un valor atípico, incluirlo infla el $s$ ordinario, que está en el denominador de $r_i$ y esconde el mismo caso que estás probando. Estimar la desviación estándar del error a partir de los otros casos, $s_{(i)}$, quita ese autoenmascaramiento y da una distribución $t$ exacta para $t_i$.

P6. ¿Aplica Durbin-Watson a cualquier regresión? No. Prueba la correlación entre errores consecutivos, así que solo significa algo cuando las filas tienen un orden real, por lo general el tiempo. Aplicarlo a datos transversales como los países de savings no tiene sentido, porque “consecutivo” es solo el orden alfabético. El Capítulo 15 lo usa donde corresponde.

P7. La distancia de Cook no tiene valor p. ¿Qué tan grande es demasiado grande? La distancia de Cook es una escala descriptiva, no una prueba formal. La guía común marca $D_i$ cercano o por encima de 1, pero la lectura más útil es relativa: busca casos que se separen mucho del resto en un gráfico de índice. Libia con 0.27 está por debajo de 1 pero claramente separada del grupo, lo que basta para justificar una mirada más de cerca.

P8. Si un modelo falla un diagnóstico, ¿es inútil? No. Un diagnóstico te dice qué supuesto se rompió y por tanto de qué conclusiones desconfiar. El ajuste de las Galápagos aún describe una asociación real entre la geografía y la riqueza de especies; lo que no puede sostener es un intervalo de predicción de teoría normal o una afirmación sobre el conteo exacto de una isla pequeña. Conocer los límites de un modelo es lo que te permite usar la parte que se sostiene.

9.9 Problemas de práctica¶

1. (A) En una oración cada uno, define el valor de apalancamiento, el valor atípico y el punto influyente, y di cuáles dos se combinan para formar el tercero.

2. (A) En la tarea de alguien sale un valor de apalancamiento $h_{ii} = 1.3$. Explica cómo sabes que está mal sin ver los datos.

3. (A) ¿Por qué un punto de alto apalancamiento tiende a tener un residuo crudo pequeño? Remítete a la fórmula de la varianza de $e_i$.

4. (A) Enuncia la regla $2p/n$ y calcula el punto de corte para el modelo de savings. ¿Qué único hecho sobre la matriz sombrero hace que $p/n$ sea el valor de apalancamiento promedio?

5. (A) Explica en palabras sencillas por qué el residuo eliminado studentizado usa $s_{(i)}$ en lugar de $s$.

6. (A) Da el patrón de residuos contra valores ajustados que esperarías para (i) varianza no constante, (ii) un término cuadrático faltante, (iii) un solo valor atípico grave.

7. (A) Un colega elimina cada punto con $|t_i| > 2$ y reporta el ajuste limpio. Da dos razones por las que esto es mala práctica.

8. (A) El estadístico de Durbin-Watson de savings es 1.93. Explica por qué esto no es evidencia de que los errores de savings sean independientes.

9. (A) Interpreta el DFBETAS de Libia de -1.02 para el coeficiente de ddpi para un lector que nunca ha visto el estadístico.

10. (B) Usando la simetría y la idempotencia de $\mathbf{H}$, demuestra que $0 \le h_{ii} \le 1$ y que $\sum_i h_{ii} = p$ (Teorema 9.2).

11. (B) Deriva la fórmula del valor de apalancamiento en regresión simple $h_{ii} = 1/n + (X_i - \bar X)^2/S_{xx}$ a partir de $h_{ii} = \mathbf{x}_i'(\mathbf{X}'\mathbf{X})^{-1}\mathbf{x}_i$ (Teorema 9.3).

12. (B) Demuestra $\operatorname{Var}\{\mathbf{e}\} = \sigma^2(\mathbf{I} - \mathbf{H})$, y deduce tanto $\operatorname{Var}\{e_i\}$ como $\operatorname{Cov}\{e_i, e_j\}$ (Teorema 9.5).

13. (B) Verifica la fórmula de Sherman-Morrison para $(\mathbf{A} - \mathbf{u}\mathbf{u}')^{-1}$ multiplicándola por $\mathbf{A} - \mathbf{u}\mathbf{u}'$ y simplificando hasta $\mathbf{I}$.

14. (B) Partiendo de la fórmula de eliminación $\mathbf{b} - \mathbf{b}_{(i)} = (\mathbf{X}'\mathbf{X})^{-1}\mathbf{x}_i e_i/(1 - h_{ii})$ (Teorema 9.10), deriva la forma cerrada $D_i = (r_i^2/p)\,h_{ii}/(1 - h_{ii})$ (Teorema 9.11).

15. (B) Muestra que el residuo eliminado satisface $d_i = e_i/(1 - h_{ii})$ y que $\operatorname{Var}\{d_i\} = \sigma^2/(1 - h_{ii})$.

16. (B) Deriva la descomposición de falta de ajuste $\mathrm{SSE} = \mathrm{SSPE} + \mathrm{SSLF}$ desarrollando $Y_{jk} - \hat Y_j$ y mostrando que el término cruzado se anula (Teorema 9.13).

17. (B) Explica por qué $E\{\mathrm{MSPE}\} = \sigma^2$ sea o no correcto el modelo lineal, pero $E\{\mathrm{MSLF}\} > \sigma^2$ cuando no lo es. Identifica los modelos completo y reducido detrás de la prueba F de falta de ajuste.

18. (B) Muestra algebraicamente que $\mathrm{DFFITS}_i = t_i\sqrt{h_{ii}/(1 - h_{ii})}$ partiendo del cambio en el único valor ajustado $\hat Y_i - \hat Y_{i(i)}$ y la fórmula de eliminación.

19. (B) Usando $D = 2(1 - \hat\rho)$, encuentra el valor de $\hat\rho$ que implica un estadístico de Durbin-Watson de 0.67, y explica qué significa esa autocorrelación para los errores.

20. (C) Ajusta el modelo gala. Reporta los valores de apalancamiento por encima de $2p/n$, y confirma que los valores de apalancamiento suman $p$.

21. (C) Para el modelo de savings, calcula las cuatro escalas del residuo para los Estados Unidos y Japón, y comenta cómo la studentización cambia la imagen para cada uno.

22. (C) Aplica la prueba de valor atípico de Bonferroni al modelo de savings. Reporta el residuo eliminado studentizado más extremo, el punto de corte, y tu decisión.

23. (C) Calcula la distancia de Cook para el modelo gala, haz un gráfico de índice, y confirma el valor de Isabela con la fórmula $(r_i^2/p)\,h_{ii}/(1-h_{ii})$.

24. (C) Reajusta el modelo gala sin Isabela y reporta cómo cambian los coeficientes de Elevation y Area. Comenta si las conclusiones son frágiles.

25. (C) Haz los gráficos de residuos contra valores ajustados, cuantil-cuantil normal y escala-ubicación para el modelo de savings, y argumenta a partir de ellos que sus supuestos son aceptables.

26. (C) Aplica las pruebas de Breusch-Pagan y Shapiro-Wilk al modelo gala, luego a un modelo que use $\log(\text{Species})$ como respuesta. ¿Ayuda el logaritmo? (Esto adelanta el Capítulo 10.) Aplica también la prueba de Brown-Forsythe (Levene modificada) para la varianza constante sobre el modelo gala crudo dividiendo los casos en la mediana del valor ajustado y aplicando una prueba $t$ de dos muestras a las desviaciones absolutas de los residuos respecto a sus medianas de grupo; ¿coincide con Breusch-Pagan?

27. (C) Usando toluca.csv, aplica la prueba F de falta de ajuste a mano: construye las medias de grupo, calcula SSPE y SSLF, forma $F^\ast$, y compáralo con anova.

28. (C) Simula una regresión con errores positivamente autocorrelacionados (semilla 4210, $\rho = 0.8$, $n = 40$) y una con errores independientes, y compara sus estadísticos de Durbin-Watson y valores p.

29. (B) Un punto tiene $h_{ii} = 0.8$ y cae exactamente sobre la recta ajustada ($e_i = 0$). Calcula su distancia de Cook y DFFITS, y explica por qué un punto de alto apalancamiento puede tener influencia cero.

30. (C) Para el modelo gala, produce el gráfico de influencia (valor de apalancamiento contra residuo eliminado studentizado, tamaño de burbuja distancia de Cook) e identifica cada isla que sea inusual en al menos un eje.

9.10 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 reportar un número. Donde una pregunta muestra salida de software, se produjo en la máquina del curso a partir de los mismos archivos de data/ que usaste todo el semestre; lees los números del impreso y los interpretas. Trabaja cada una antes de abrir la respuesta modelo.

EP 9.1. El modelo de savings sr ~ pop15 + pop75 + dpi + ddpi ($n = 50$, $p = 5$) se ajusta, y las escalas del residuo para Chile y Zambia, los dos países con los residuos crudos más grandes, se imprimen abajo junto con el punto de corte de valor atípico de Bonferroni.

savings <- read.csv("data/savings.csv")
sfit <- lm(sr ~ pop15 + pop75 + dpi + ddpi, data = savings)
n <- nrow(savings); p <- length(coef(sfit))
i <- which(savings$country %in% c("Chile", "Zambia"))
data.frame(country = savings$country[i],
           hii = round(hatvalues(sfit)[i], 3),
           raw_e = round(residuals(sfit)[i], 3),
           standardized = round(rstandard(sfit)[i], 3),
           stud_deleted = round(rstudent(sfit)[i], 3))
qt(1 - 0.05 / (2 * n), n - p - 1)

 country   hii  raw_e standardized stud_deleted
   Chile 0.037 -8.242       -2.209       -2.313
  Zambia 0.064  9.751        2.651        2.854
[1] 3.5258

Explica por qué el residuo de Zambia crece de 2.564 en la escala semiestudentizada (su residuo crudo 9.751 dividido entre $s = 3.803$) a 2.854 una vez que se studentiza y elimina, y di qué característica de Zambia en la tabla impulsa ese aumento. Luego indica si la prueba de valor atípico marca a Zambia como un valor atípico genuino, nombra el número contra el que lo comparas, y explica en una oración de qué te protege la corrección de Bonferroni.

EP 9.2. El mismo ajuste de savings da el resumen de influencia de abajo para sus tres países más comentados.

 country   hii rstudent cooksD dffits
   Japan 0.223    1.603  0.143  0.860
  Zambia 0.064    2.854  0.097  0.748
   Libya 0.531   -1.089  0.268 -1.160
   ( 2p/n = 0.200,   DFFITS cutoff 2*sqrt(p/n) = 0.632 )

Un estudiante escribe: “Zambia tiene el residuo studentizado más grande de cualquier país, así que Zambia es el único punto que más distorsiona el ajuste de savings, y reajustaría sin él”. Evalúa esta afirmación. Indica si es correcta, corrige el razonamiento usando los números de la tabla, y nombra el país que de verdad influye más en el ajuste, con la medida que lo zanja.

EP 9.3. Para el modelo de las Galápagos Species ~ Area + Elevation + Nearest + Scruz + Adjacent ($n = 30$, $p = 6$), un evaluador aplica las pruebas de Breusch-Pagan y Shapiro-Wilk e inspecciona los valores ajustados.

gala <- read.csv("data/gala.csv")
gfit <- lm(Species ~ Area + Elevation + Nearest + Scruz + Adjacent, data = gala)
library(lmtest)
bptest(gfit)
shapiro.test(residuals(gfit))
range(fitted(gfit)); sum(fitted(gfit) < 0)

	studentized Breusch-Pagan test
data:  gfit
BP = 9.7959, df = 5, p-value = 0.08123

	Shapiro-Wilk normality test
data:  residuals(gfit)
W = 0.91351, p-value = 0.01826

[1] -36.38392 386.40356
[1] 5

Un estudiante concluye: “El valor p de Breusch-Pagan es 0.081, que está por encima de 0.05, así que el supuesto de varianza constante se cumple y el modelo de las Galápagos está bien”. Evalúa esta conclusión. Explica qué muestra y qué no muestra la prueba aquí, aporta la otra evidencia del impreso, e indica la regla general sobre pruebas contra gráficos que ilustra este caso.

EP 9.4. Continuando con el ajuste de las Galápagos, un analista teme que Isabela, la isla más grande, esté dirigiendo todo el modelo. El ajuste se reajusta con Isabela eliminada y se comparan los coeficientes.

iso <- which(gala$island == "Isabela")
c(h = hatvalues(gfit)[iso], rstudent = rstudent(gfit)[iso],
  cooksD = cooks.distance(gfit)[iso])
gfit2 <- lm(Species ~ Area + Elevation + Nearest + Scruz + Adjacent,
            data = gala[-iso, ])
round(rbind(all = coef(gfit), no_Isabela = coef(gfit2)), 4)
c(R2_all = summary(gfit)$r.squared, R2_no_Isabela = summary(gfit2)$r.squared)

       h.16 rstudent.16   cooksD.16
     0.9685     -5.3337     68.0764
           (Intercept)    Area Elevation Nearest   Scruz Adjacent
all             7.0682 -0.0239    0.3195  0.0091 -0.2405  -0.0748
no_Isabela     22.5861  0.2957    0.1404 -0.2552 -0.0901  -0.0650
   R2_all R2_no_Isabela
    0.766         0.871

Interpreta esta salida en contexto. Di qué cambiaría en las conclusiones reportadas sobre Area y Elevation si se descartara Isabela, explica por qué una sola isla puede hacer esto usando sus tres números de diagnóstico, e indica cuál es la forma honesta de reportar un resultado tan frágil.

EP 9.5. Los datos de Toluca registran las hours de trabajo contra el lotsize para 25 corridas de producción que usaron solo 11 tamaños de lote distintos. La prueba F de falta de ajuste compara el ajuste de recta con un modelo que tiene una media libre en cada tamaño de lote.

toluca <- read.csv("data/toluca.csv")
reduced <- lm(hours ~ lotsize, data = toluca)
full    <- lm(hours ~ factor(lotsize), data = toluca)
anova(reduced, full)

Analysis of Variance Table

Model 1: hours ~ lotsize
Model 2: hours ~ factor(lotsize)
  Res.Df   RSS Df Sum of Sq      F Pr(>F)
1     23 54825
2     14 37581  9     17245 0.7138 0.6893

Interpreta esta salida en contexto: identifica qué número impreso es la suma de cuadrados del error puro y cuál es la suma de cuadrados de la falta de ajuste, enuncia la conclusión de la prueba en lenguaje sencillo, y explica qué compara la prueba para llegar a ella. Luego responde: ¿qué cambiaría en tu capacidad de realizar esta prueba si cada una de las 25 corridas hubiera usado un tamaño de lote distinto?

Juego del capítulo¶