Sesión 26 (A) — Laboratorio integrador II: palanca, distancia de Cook y el diccionario geometría ↔ probabilidad

La actividad 2 (html) de la práctica A de la sesión 25 dejó estimado el modelo del Boston original: el logaritmo del precio mediano sobre rm, lstat, dis y la indicadora chas. Esta práctica pregunta cuánto depende ese ajuste de barrios concretos. La actividad 6 (html) de la práctica A de la sesión 14 hizo la misma pregunta con catorce viviendas: las quitó una a una y miró cómo se movían los coeficientes. Esta práctica usa dos medidas que resumen esa comprobación para las \(506\) observaciones a la vez: la palanca de cada barrio, que es geometría de la lección 10, y la distancia de Cook, que combina la palanca con el residuo.

La práctica cierra con un texto sin código: el diccionario que traduce la geometría de \(\mathbb R^n\) de la primera mitad del curso al lenguaje de la probabilidad de la segunda. Es una analogía explícita, no una demostración, y no es materia de examen.

Objetivos

1. Calcular la palanca de cada barrio con leverage y leerla como la longitud al cuadrado de una proyección: cuánto pesa cada fila en su propio valor ajustado.
2. Calcular la distancia de Cook a partir de la palanca y el residuo, localizar el barrio más influyente y comprobar qué cambia al quitarlo.
3. Reconocer un grupo de observaciones influyentes, los \(16\) barrios con el precio censurado en \(50\), y ver qué cambia sin ellos.
4. Leer el diccionario geometría ↔ probabilidad del curso.

Comandos nuevos de esta práctica. Dos:

leverage

tras estimar un modelo con ols, calcula para cada observación su palanca y dos medidas más que esta práctica no usa, influ y dffits (la segunda, en Belsley, Kuh y Welsch); con --save guarda las tres como series, lever, influ y dffits, y con --quiet no imprime la tabla de las \(506\) filas.

imaxc

dada una matriz, devuelve la fila en que cada columna alcanza su máximo; con {lever}, la serie convertida en matriz de una columna, da el número de la observación con la palanca máxima.

Los demás ya se han usado: open, series, ols y sus accesores ($uhat, $sigma, $coeff, $stderr, $nobs), summary --simple, smpl --restrict y smpl full, printf, y las funciones sum, min, max, mean y median.

Importe boston.csv, estime el modelo con chas (Modelo -> Mínimos cuadrados ordinarios) y pida la palanca de cada observación (en la ventana del modelo, Análisis -> Observaciones influyentes).

El valor ajustado es la proyección de \(\boldsymbol y\) sobre el subespacio generado por las columnas de \(\boldsymbol{\mathsf X}\) (lección 10), y la componente \(i\) de esa proyección es una combinación lineal de todas las componentes de \(\boldsymbol y\). El peso de \(y_i\) en su propio valor ajustado \(\hat y_i\) es la palanca \(h_i\): la longitud al cuadrado de la proyección sobre ese subespacio del vector canónico \(\boldsymbol e_i\), el que tiene un \(1\) en la fila \(i\) y ceros en las demás. Es un número entre \(0\) y \(1\), y las \(n\) palancas suman \(k\), la dimensión del subespacio; su media es \(k/n\). Una palanca grande dice que la fila \(i\) está lejos de las demás en el espacio de los regresores: el ajuste pasa cerca de ella porque pocas filas más ocupan esa zona de los regresores.

leverage --quiet --save guarda las tres series. summary lever --simple da sus estadísticos; sum(lever) comprueba que suman \(k=5\) y median(lever) da la mediana, que summary --simple no imprime. El printf compara con la convención de los manuales, \(2k/n\), y cuenta cuántos barrios la superan. imaxc({lever}) devuelve la fila de la palanca máxima, y el último printf imprime los datos de ese barrio.

o bien teclee en línea de comandos:

open ../../data/boston.csv --quiet
series lmedv = log(medv)

ols lmedv const rm lstat dis chas --quiet
scalar k = 5
leverage --quiet --save
summary lever --simple
printf "Suma de las palancas: %.2f ; media k/n = %.4f ; mediana %.4f ; convencion 2k/n = %.4f ; barrios por encima: %d\n", sum(lever), k/$nobs, median(lever), 2*k/$nobs, sum(lever > 2*k/$nobs)
scalar i = imaxc({lever})
printf "Palanca maxima %.4f en el barrio %d: medv %.1f  rm %.2f  lstat %.2f  dis %.2f  chas %d\n", lever[i], i, medv[i], rm[i], lstat[i], dis[i], chas[i]

Estadísticos principales, usando las observaciones 1 - 506
para la variable 'lever' (506 observaciones válidas)

  Media                        0,0098814
  Mínimo                       0,0021583
  Máximo                        0,075518
  Desviación típica            0,0094198
  Observaciones ausentes               0

Suma de las palancas: 5,00 ; media k/n = 0,0099 ; mediana 0,0061 ; convencion 2k/n = 0,0198 ; barrios por encima: 62
Palanca maxima 0,0755 en el barrio 366: medv 27,5  rm 3,56  lstat 7,12  dis 1,61  chas 0

Una palanca grande dice que la fila puede mover el ajuste; si además su residuo es grande, lo mueve. La distancia de Cook combina las dos cosas: \[D_i=\frac{\hat e_i^2\,h_i}{k\,\mathfrak s^2\,(1-h_i)^2},\] con \(\hat e_i\) el residuo, \(h_i\) la palanca y \(\mathfrak s^2\) la cuasivarianza de los residuos. Mide cuánto se desplaza el vector de valores ajustados al quitar la observación \(i\), en unidades de \(\mathfrak s^2\) y por dimensión. Gretl no la imprime, pero se calcula en una línea con las series que leverage guardó.

$uhat es la serie de residuos del último modelo y $sigma la desviación típica de la regresión, \(\mathfrak s\), de modo que $sigma^2 es \(\mathfrak s^2\). La serie cook aplica la fórmula fila a fila. El printf cuenta cuántos barrios superan la convención \(4/n\) y localiza el máximo; el siguiente imprime el residuo y la distancia de Cook del barrio con la palanca máxima, el \(366\) de la Actividad 1. Después, smpl cook < max(cook) --restrict deja fuera ese barrio, se vuelve a estimar y smpl full recupera la muestra entera.

series e = $uhat
scalar s2 = $sigma^2
series cook = e^2 * lever / (k * s2 * (1 - lever)^2)
printf "Convencion 4/n = %.4f ; barrios por encima: %d\n", 4/$nobs, sum(cook > 4/$nobs)
scalar j = imaxc({cook})
printf "Cook maxima %.3f en el barrio %d: medv %.1f  rm %.2f  lstat %.2f  dis %.2f  chas %d  palanca %.4f  residuo %.3f\n", cook[j], j, medv[j], rm[j], lstat[j], dis[j], chas[j], lever[j], e[j]
printf "Barrio %d, el de palanca maxima: residuo %.3f  Cook %.4f ; barrios con Cook mayor: %d\n", i, e[i], cook[i], sum(cook > cook[i])
printf "Con todos:        chas %.4f  rm %.4f  lstat %.4f  dis %.4f\n", $coeff(chas), $coeff(rm), $coeff(lstat), $coeff(dis)
smpl cook < max(cook) --restrict
ols lmedv const rm lstat dis chas --quiet
printf "Sin el barrio %d: chas %.4f  rm %.4f  lstat %.4f  dis %.4f  (n = %d)\n", j, $coeff(chas), $coeff(rm), $coeff(lstat), $coeff(dis), $nobs
smpl full

Convencion 4/n = 0,0079 ; barrios por encima: 40
Cook maxima 0,113 en el barrio 365: medv 21,9  rm 8,78  lstat 5,29  dis 1,90  chas 1  palanca 0,0542  residuo -0,699
Barrio 366, el de palanca maxima: residuo 0,387  Cook 0,0506 ; barrios con Cook mayor: 6
Con todos:        chas 0,1597  rm 0,1203  lstat -0,0394  dis -0,0062
Sin el barrio 365: chas 0,1785  rm 0,1292  lstat -0,0391  dis -0,0068  (n = 505)

En el Boston original el precio mediano está censurado: el valor máximo, \(50\), aparece \(16\) veces, y en el fichero significa ``\(50\) o más''. Esos \(16\) barrios tienen un precio real desconocido y un residuo del mismo signo en todos ellos. Compare su distancia de Cook con la de los demás y estime el modelo sin ellos.

sum(medv = 50)= cuenta los barrios censurados, porque la comparación vale \(1\) cuando es cierta. Las dos medias de cook se calculan como la media de grupo de la lección 16 (``Una indicadora desplaza el intercepto: el río Charles''): la media del producto de cook por la indicadora del grupo, dividida por la media de la indicadora. El primer smpl deja solo los \(16\) barrios para imprimir el menor y el mayor de sus residuos, sus medias de rm, lstat y palanca, cuántos tienen palanca mayor que la media \(k/n\) y cuántos superan la convención \(4/n\) de la distancia de Cook; smpl medv < 50 --restrict los deja fuera.

printf "Barrios con medv = 50: %d\n", sum(medv == 50)
printf "Cook media de esos barrios: %.4f ; de los demas: %.4f\n", mean(cook * (medv == 50)) / mean(medv == 50), mean(cook * (medv < 50)) / mean(medv < 50)
smpl medv == 50 --restrict
printf "Residuos de esos barrios: minimo %.3f ; maximo %.3f\n", min(e), max(e)
printf "Medias de esos barrios: rm %.2f  lstat %.2f  palanca %.4f ; con palanca mayor que k/n: %d ; con Cook mayor que 4/n: %d\n", mean(rm), mean(lstat), mean(lever), sum(lever > k/506), sum(cook > 4/506)
smpl full
smpl medv < 50 --restrict
ols lmedv const rm lstat dis chas
printf "Sin los 16: chas %.4f (ee %.4f)  rm %.4f  lstat %.4f  dis %.4f\n", $coeff(chas), $stderr(chas), $coeff(rm), $coeff(lstat), $coeff(dis)
smpl full

Barrios con medv = 50: 16
Cook media de esos barrios: 0,0161 ; de los demas: 0,0026
Residuos de esos barrios: minimo 0,099 ; maximo 0,757
Medias de esos barrios: rm 7,48  lstat 4,36  palanca 0,0279 ; con palanca mayor que k/n: 15 ; con Cook mayor que 4/n: 8

Modelo 1: MCO, usando las observaciones 1-490
Variable dependiente: lmedv

             coeficiente   Desv. típica   Estadístico t   valor p 
  ----------------------------------------------------------------
  const       2,79539       0,144073         19,40        1,72e-62 ***
  rm          0,107222      0,0193339         5,546       4,82e-08 ***
  lstat      -0,0369454     0,00204660      -18,05        4,23e-56 ***
  dis         0,00286834    0,00563626        0,5089      0,6110  
  chas        0,123365      0,0422265         2,922       0,0036   ***

Media de la vble. dep.  3,005859   D.T. de la vble. dep.   0,382794
Suma de cuad. residuos  23,33702   D.T. de la regresión    0,219357
R-cuadrado              0,674309   R-cuadrado corregido    0,671623
F(4, 485)               251,0351   Valor p (de F)          1,2e-116
Log-verosimilitud       50,58940   Criterio de Akaike     -91,17880
Criterio de Schwarz    -70,20677   Crit. de Hannan-Quinn  -82,94233

Log-verosimilitud para medv = -1422,28

Sin considerar la constante, el valor p más alto fue el de la variable 8 (dis)

Sin los 16: chas 0,1234 (ee 0,0422)  rm 0,1072  lstat -0,0369  dis 0,0029

El curso ha presentado la misma construcción dos veces. Las lecciones 3 a 8 la presentaron con vectores de \(\mathbb R^n\): la media es el coeficiente de la proyección sobre la constante, la regresión es una proyección sobre el subespacio de los regresores, el residuo es ortogonal a ese subespacio y el \(R^2\) es el cuadrado de un coseno. Las lecciones 9 a 12 la presentaron con variables aleatorias: el valor esperado es el coeficiente de la proyección sobre la constante, la esperanza condicional es una proyección sobre las funciones de los regresores, la perturbación es ortogonal a esas funciones. Las dos construcciones coinciden porque los dos espacios tienen un producto escalar y una proyección, y la tabla recoge la traducción, ampliando la tabla-diccionario con la que cerró la lección 9. Es una analogía explícita, que el apéndice de geometría desarrolla en sus secciones 1 a 4, desde el apéndice (sección 1) (html), y que aquí no se demuestra.

La columna de la izquierda es lo que Gretl calcula; la de la derecha es lo que se quiere conocer. El puente entre las dos es el muestreo aleatorio de la lección 11: los datos son una columna de la matriz de copias, y cada fórmula de la izquierda, evaluada sobre la muestra aleatoria, es un estadístico que estima la cantidad de la derecha.¹

1. Las palancas suman \(5\) y su media es \(0{,}0099\). ¿Por qué suman el número de columnas de \(\boldsymbol{\mathsf X}\) y no el de filas?
2. El barrio \(366\) tiene la palanca máxima y el \(365\) la distancia de Cook máxima. ¿Por qué no coinciden?
3. ¿Qué cambia al quitar el barrio \(365\) y qué se concluye de ello? ¿Es una razón para quitarlo?
4. Los \(16\) barrios con medv igual a \(50\) tienen una distancia de Cook media seis veces mayor que los demás. ¿Qué tienen en común que explique esa influencia?
5. En la tabla del diccionario, la media de los datos y el valor esperado de una variable aleatoria ocupan la misma fila. ¿En qué sentido son el mismo objeto y en qué sentido no lo son?

• Wooldridge, J. M. (2020), Introductory Econometrics, cap. 9, sección 9.5 (observaciones influyentes).
• Belsley, D. A., Kuh, E. y Welsch, R. E. (1980), Regression Diagnostics, Wiley: la palanca, DFFITS y las convenciones \(2k/n\) y \(2\sqrt{k/n}\).
• Cook, R. D. (1977), ``Detection of influential observation in linear regression'', Technometrics 19, 15–18.
• Apéndice de geometría del curso, secciones 1 a 4: el espacio de las variables aleatorias, las copias idénticas e independientes, la matriz de copias y el puente entre los dos espacios.

1. La palanca \(h_i\) es la longitud al cuadrado de la proyección de \(\boldsymbol e_i\) sobre el subespacio de los regresores. Los \(n\) vectores canónicos forman una base ortonormal de \(\mathbb R^n\), y la suma de las longitudes al cuadrado de sus proyecciones sobre un subespacio es la dimensión de ese subespacio, \(k\), porque con una base ortonormal \(\boldsymbol b_1,\ldots,\boldsymbol b_k\) del subespacio, \[\sum_{i=1}^n\Vert P\boldsymbol e_i\Vert_e^2=\sum_{i=1}^n\sum_{j=1}^k\langle\boldsymbol e_i,\boldsymbol b_j\rangle_e^2=\sum_{j=1}^k\sum_{i=1}^n b_{ji}^2=\sum_{j=1}^k\Vert\boldsymbol b_j\Vert_e^2=k;\] esa suma es la traza de la matriz de proyección.² Si el subespacio fuese todo \(\mathbb R^n\), cada palanca valdría \(1\) y el ajuste reproduciría los datos.
2. La palanca solo depende de los regresores, y la distancia de Cook multiplica la palanca por el residuo al cuadrado. El barrio \(366\) está lejos de los demás en rm, y su residuo, \(0{,}39\), no es pequeño: su distancia de Cook es \(0{,}051\), la séptima del fichero. El \(365\) tiene menos palanca, \(0{,}054\) frente a \(0{,}076\), pero su residuo, \(-0{,}70\), es casi el doble en valor absoluto y, al cuadrado, más del triple: por eso su distancia de Cook, \(0{,}113\), es la mayor.
3. El coeficiente de chas pasa de \(0{,}160\) a \(0{,}179\): la diferencia entre barrios junto al río y sin río depende en un \(12\) % de un solo barrio. Esa dependencia no es una razón para quitarlo. La distancia de Cook señala qué observación mirar; quitarla exige una razón sobre los datos, como un error de registro. Si el barrio es real, se queda, y el informe dice de qué depende el resultado.
4. Su precio registrado, \(50\), es una cota inferior del precio real, y sus residuos tienen todos el mismo signo: son positivos, el ajuste queda por debajo de \(50\) en los \(16\), y más por debajo aún del precio verdadero. Que los \(16\) queden por debajo es sistemático, no azar. Además son barrios con rm alto y lstat bajo, con medias de \(7{,}48\) habitaciones y un \(4{,}36\) % de población de estatus bajo, y \(15\) de los \(16\) tienen palanca mayor que la media \(k/n\). Un grupo de observaciones puede ser influyente aunque ninguna lo sea mucho por separado.
5. Los dos son el coeficiente de la proyección de un objeto sobre la constante, en un espacio con producto escalar: \(\mu_{\boldsymbol x}=\langle\boldsymbol x,\boldsymbol 1\rangle_s\) en \(\mathbb R^n\) y \(\mathrm E(X)=\langle X,\mathit 1\rangle_P\) en el espacio de las variables aleatorias. No son el mismo objeto: la media es un número calculado con los datos observados, y el valor esperado es un número de la población, que no se observa. La media muestral, la fórmula de la media evaluada sobre la muestra aleatoria, es la variable aleatoria que estima el valor esperado; que lo estime sin sesgo es la lección 11.