6. Álgebra de matrices para la regresión (en espanol) - Regression Modeling and Analysis (MATH 4210)

6. Álgebra de matrices para la regresión¶

Dwaine Studios opera una cadena de estudios de retratos y quiere crecer. Ya trabaja en 21 ciudades, y sus gerentes registran dos números para cada una: la población de 16 años o menos, el grupo de niños cuyos padres compran retratos, y el ingreso disponible per cápita, cuánto dinero tienen las familias para gastar. Quieren usar estos dos números para predecir las ventas de los estudios, de modo que puedan ordenar ciudades nuevas y abrir donde el pronóstico sea más alto. La pregunta es de pronóstico, pero por debajo es de modelado: ¿cómo dependen las ventas de ambos predictores a la vez?

En el Capítulo 2 ajustaste una recta con un solo predictor resolviendo dos ecuaciones normales para el intercepto y la pendiente (2.2 Mínimos cuadrados desde los primeros principios). Eso funcionó porque solo había dos incógnitas. Dwaine tiene tres, un intercepto y dos pendientes, y un estudio realista podría tener veinte. Escribir una ecuación normal separada para cada una y resolver la maraña a mano no es algo que nadie quiera hacer dos veces. El álgebra de matrices es el mejor lenguaje: empaqueta un conjunto de datos completo en dos símbolos, las ecuaciones normales en una sola, y la solución en una fórmula que escribes en una sola línea y calculas con un solo comando.

Nada de lo que sigue requiere que hayas visto una matriz antes. Construimos cada idea a partir de su definición, la probamos en un ejemplo pequeño que podrías verificar con lápiz, y luego la ejecutamos sobre el conjunto de datos real de 21 ciudades tanto en R como en Python. Si has tomado álgebra lineal, esto se leerá como un repaso con acento estadístico; si no, de todos modos terminarás capaz de ajustar y explicar una regresión múltiple, la única álgebra de matrices que este curso te pide.

Un diagrama de dispersión de las ventas de los estudios en miles de dólares en el eje vertical contra la población objetivo menor de 16 años en miles en el eje horizontal, para 21 ciudades. Los puntos suben de abajo a la izquierda hacia arriba a la derecha. Cada punto está sombreado de azul claro a azul oscuro según su ingreso disponible, y los puntos más oscuros (de mayor ingreso) tienden a situarse más arriba para una población dada. — Figure 1:Las ventas de Dwaine Studios contra la población menor de 16 años de una ciudad, con el ingreso disponible indicado por el sombreado. Ambos predictores acompañan mayores ventas, así que un buen pronóstico necesita usarlos juntos, que es exactamente lo que nos permite hacer la maquinaria matricial de este capítulo.

La Figure 1 muestra que ambos predictores importan: las ventas suben a medida que crece la población joven, y para una población dada las ciudades de mayor ingreso (puntos más oscuros) tienden a vender más. Un modelo que usa ambos a la vez es un modelo de regresión múltiple, y todo modelo de este tipo se construye, ajusta y comprende mediante matrices. Al terminar habrás construido, para los datos reales de Dwaine, la matriz de diseño $\mathbf{X}$ , el producto $\mathbf{X}'\mathbf{X}$ , su inversa, las estimaciones de los coeficientes, la matriz sombrero, y el error estándar de cada coeficiente, cada uno a partir de operaciones que también puedes hacer a mano en un ejemplo pequeño.

Esta lección de un vistazo

Qué estamos haciendo: Aprender apenas lo suficiente de álgebra de matrices para escribir y resolver una regresión con varios predictores a la vez, usando los datos de 21 ciudades de Dwaine Studios como ejemplo conductor.
Por qué lo hacemos: Las dos ecuaciones normales del Capítulo 2 (2.2 Mínimos cuadrados desde los primeros principios) funcionaron para un predictor, pero Dwaine tiene dos y un estudio real puede tener veinte; resolverlas predictor por predictor a mano es imposible, así que necesitamos un lenguaje compacto.
Objetivo principal: Representar los datos como $\mathbf{X}$ y $\mathbf{Y}$ , formar $\mathbf{X}'\mathbf{X}$ y $\mathbf{X}'\mathbf{Y}$ , resolver $\mathbf{b} = (\mathbf{X}'\mathbf{X})^{-1}\mathbf{X}'\mathbf{Y}$ , construir la matriz sombrero, y leer la varianza de cada coeficiente en $\sigma^2(\mathbf{X}'\mathbf{X})^{-1}$ , comprobando cada paso primero en un ejemplo pequeño.
Qué cambió respecto de los capítulos anteriores: Hasta ahora, un predictor y álgebra escalar; de aquí en adelante, muchos predictores empaquetados en vectores y matrices, todavía la misma idea de mínimos cuadrados vista a través de una notación nueva. La $\mathbf{X}'\mathbf{X}$ exacta de Dwaine que construimos aquí es el objeto que reutiliza el Capítulo 7 (7.1 El modelo y los mínimos cuadrados en forma matricial).

Objetivos de aprendizaje

Al terminar este capítulo podrás:

Representar un conjunto de datos de regresión como un vector respuesta $\mathbf{Y}$ y una matriz de diseño $\mathbf{X}$ , y escribir el modelo como $\mathbf{Y} = \mathbf{X}\boldsymbol{\beta} + \boldsymbol{\varepsilon}$ .
Calcular transpuestas y productos de matrices, y en particular construir $\mathbf{X}'\mathbf{X}$ y $\mathbf{X}'\mathbf{Y}$ , a mano y con software.
Decidir cuándo una matriz cuadrada tiene inversa, usando el rango y el determinante, y explicar qué le hace la colinealidad perfecta a $\mathbf{X}'\mathbf{X}$ .
Resolver las ecuaciones normales $\mathbf{X}'\mathbf{X}\,\mathbf{b} = \mathbf{X}'\mathbf{Y}$ para el vector de coeficientes $\mathbf{b} = (\mathbf{X}'\mathbf{X})^{-1}\mathbf{X}'\mathbf{Y}$ .
Verificar que la matriz sombrero $\mathbf{H} = \mathbf{X}(\mathbf{X}'\mathbf{X})^{-1}\mathbf{X}'$ es simétrica e idempotente, y explicar por qué los valores ajustados y los residuos son proyecciones.
Expresar SSE, SSR y SSTO como formas cuadráticas, y conectar la definición positiva con sumas de cuadrados no negativas.
Derivar el vector de medias y la matriz de covarianzas de un vector aleatorio, y usarlos para mostrar que $\operatorname{Cov}\{\mathbf{b}\} = \sigma^2 (\mathbf{X}'\mathbf{X})^{-1}$ .
Enunciar el modelo normal multivariante para $\mathbf{Y}$ y $\mathbf{b}$ y describir cómo impulsa la inferencia de los capítulos posteriores.

6.1 Los datos como un vector y una matriz¶

Intuición¶

Empieza por la respuesta. Dwaine tiene 21 cifras de ventas, una por ciudad. Apílalas en una sola columna y tendrás un vector (Definición 6.1), una lista ordenada de números escrita verticalmente. Llámalo $\mathbf{Y}$ . Un vector es la matriz más simple: muchas filas, una columna.

Ahora los predictores. Cada ciudad lleva dos números, así que los predictores forman una tabla con 21 filas y 2 columnas, y una tabla de números con filas y columnas es una matriz (Definición 6.1). Para la regresión pegamos una columna extra de puros unos al frente: un truco de contabilidad que deja que el intercepto viaje como si fuera un coeficiente más. El resultado es la matriz de diseño $\mathbf{X}$ (Definición 6.2), con 21 filas y 3 columnas. La Figure 2 muestra cómo se alinea todo el modelo como arreglos apilados.

Un diagrama que muestra el modelo de regresión como cuatro cajas apiladas. Una caja alta y delgada rotulada Y, 21 por 1, es igual a una caja ancha rotulada X, 21 por 3, con su primera columna delineada y rotulada unos, multiplicada por una caja corta rotulada beta, 3 por 1, más una caja alta y delgada rotulada epsilon, 21 por 1. Una nota dice que las dimensiones internas coinciden: 21 por 3 por 3 por 1 da 21 por 1. — Figure 2:El modelo de regresión en forma matricial. La respuesta Y y los errores son 21 por 1, la matriz de diseño X es 21 por 3 con una columna inicial de unos para el intercepto, y el vector de coeficientes beta es 3 por 1. Las dimensiones internas coinciden, así que el producto está definido.

Leer a lo ancho de una fila da el registro completo de una ciudad: un 1, su población joven, su ingreso disponible, y, en $\mathbf{Y}$ , sus ventas. Apilar las filas convierte 21 pequeñas ecuaciones separadas en el único enunciado $\mathbf{Y} = \mathbf{X}\boldsymbol{\beta} + \boldsymbol{\varepsilon}$ . Esa compresión es la razón para aprender matrices: los mismos tres símbolos describen un conjunto de datos con 21 filas o con 21 millones.

Fórmula¶

Para una regresión con $n$ observaciones y $p$ parámetros (contando el intercepto), las piezas son

\mathbf{Y} = \begin{pmatrix} Y_1 \\ Y_2 \\ \vdots \\ Y_n \end{pmatrix}_{n \times 1}, \qquad \mathbf{X} = \begin{pmatrix} 1 & X_{11} & \cdots & X_{1,p-1} \\ 1 & X_{21} & \cdots & X_{2,p-1} \\ \vdots & \vdots & & \vdots \\ 1 & X_{n1} & \cdots & X_{n,p-1} \end{pmatrix}_{n \times p}, \qquad \boldsymbol{\beta} = \begin{pmatrix} \beta_0 \\ \beta_1 \\ \vdots \\ \beta_{p-1} \end{pmatrix}_{p \times 1}.

$\mathbf{Y}$ es el vector respuesta; la fila $i$ es la respuesta del caso $i$ .
$\mathbf{X}$ es la matriz de diseño; su primera columna es de puros unos, y la columna $k+1$ contiene los valores del predictor $k$ .
$\boldsymbol{\beta}$ es el vector de parámetros: el intercepto $\beta_0$ seguido de las pendientes.
$p$ es el número de parámetros de regresión incluyendo el intercepto ( $p = 3$ para Dwaine), y $n$ es el número de observaciones ( $n = 21$ ).

Para Dwaine, $n = 21$ y $p = 3$ : una columna de unos, una columna de poblaciones menores de 16 años, una columna de ingresos disponibles. Todo el modelo del Capítulo 2 se generaliza a $\mathbf{Y} = \mathbf{X}\boldsymbol{\beta} + \boldsymbol{\varepsilon}$ , que desarrollamos en 7.1 El modelo y los mínimos cuadrados en forma matricial. En palabras: cada respuesta observada es igual a una suma ponderada de los predictores de esa fila, usando los mismos pesos $\boldsymbol{\beta}$ para cada fila, más un error aleatorio.

R y Python¶

Ejemplo 6.1: Construir la matriz de diseño de Dwaine

Pregunta. Lee los datos de Dwaine y arma el vector respuesta $\mathbf{Y}$ y la matriz de diseño $\mathbf{X}$ , con su columna inicial de unos.

Intuición. La respuesta es simplemente la columna de ventas. La matriz de diseño es las dos columnas de predictores con una columna de unos atornillada a la izquierda.

Fórmula. $\mathbf{X}$ es $n \times p = 21 \times 3$ ; la columna 1 es de unos, la columna 2 es targtpop, la columna 3 es dispoinc.

Cálculo.

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

[1] 21  3
  targtpop dispoinc sales
1     68.5     16.7 174.4
2     45.2     16.8 164.4
3     91.3     18.2 244.2

n <- nrow(dwaine)
X <- cbind(1, dwaine$targtpop, dwaine$dispoinc)
colnames(X) <- c("intercept", "targtpop", "dispoinc")
Y <- dwaine$sales
head(X, 3)

     intercept targtpop dispoinc
[1,]         1     68.5     16.7
[2,]         1     45.2     16.8
[3,]         1     91.3     18.2

Ahora en Python, con numpy ensamblando el mismo arreglo a partir del mismo CSV.

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

dwaine = pd.read_csv("data/dwaine.csv")
print(dwaine.shape)
print(dwaine.head(3))

(21, 3)
   targtpop  dispoinc  sales
0      68.5      16.7  174.4
1      45.2      16.8  164.4
2      91.3      18.2  244.2

n = len(dwaine)
X = np.column_stack([np.ones(n), dwaine["targtpop"], dwaine["dispoinc"]])
Y = dwaine["sales"].to_numpy()
print(X[:3])

[[ 1.  68.5 16.7]
 [ 1.  45.2 16.8]
 [ 1.  91.3 18.2]]

Interpretación. La primera ciudad tiene 68,500 residentes jóvenes y $16,700 de ingreso per cápita, y el$ 1 $inicial es el marcador de posición que multiplicará al intercepto. Cada fila de$ \mathbf{X}$ se lee “uno, y luego los predictores de esta ciudad”. El cbind de R y el np.column_stack de Python construyeron el mismo arreglo a partir del mismo archivo, así que los dos lenguajes comparten un solo objeto, y todo el ajuste se vuelve una corta secuencia de operaciones matriciales.

6.2 Transpuesta y multiplicación de matrices¶

Intuición¶

Tenemos los datos en matrices; ahora necesitamos combinarlos. Dos operaciones hacen casi todo el trabajo. La transpuesta (Definición 6.3) voltea una matriz sobre su diagonal, de modo que las filas se vuelven columnas. La multiplicación de matrices (Definición 6.4) es el motor que convierte la matriz de diseño en las sumas de cuadrados y productos cruzados sobre los que corre una regresión.

La multiplicación de matrices no es lo que un principiante espera. Una entrada del producto se construye a partir de toda una fila de la matriz izquierda y toda una columna de la derecha: alinéalas, multiplica término a término, y suma. Aplicada a $\mathbf{X}'\mathbf{X}$ , esa regla produce exactamente las sumas $\sum X_{ik} X_{ij}$ que de otra forma calcularías una penosa suma a la vez. La Figure 3 lo muestra en un ejemplo pequeño.

Un diagrama de multiplicación de matrices. Una matriz A de 2 por 3 con su primera fila sombreada de azul se multiplica por una matriz B de 3 por 2 con su segunda columna sombreada de verde, dando un producto AB de 2 por 2 con la entrada superior derecha resaltada en amarillo. Debajo, la aritmética dice uno por cero más dos por uno más tres por uno igual a cinco, rotulada como fila uno de A punto columna dos de B. — Figure 3:Multiplicación de matrices, entrada por entrada. La entrada resaltada del producto en la fila 1, columna 2 proviene de la fila 1 de A y la columna 2 de B, multiplicadas término a término y sumadas. Cada entrada de un producto es uno de esos productos punto de fila por columna.

Dos productos, y solo dos, importan para el ajuste: $\mathbf{X}'\mathbf{X}$ resume los predictores por sí mismos, y $\mathbf{X}'\mathbf{Y}$ resume cómo se mueve cada predictor con la respuesta. Entre ambos guardan todo el resumen que los mínimos cuadrados necesitan, razón por la cual un conjunto de datos completo puede reducirse a un puñado de números: el ajuste no se interesa por las filas individuales una vez que esas sumas se conocen.

Fórmula¶

Imagínala como volcar una matriz sobre su esquina superior izquierda: cada fila se pone de pie para volverse columna, así que una matriz ancha se vuelca en una alta y de regreso. Un vector columna, volcado, se vuelve una sola fila, razón por la cual $\mathbf{u}'\mathbf{u}$ multiplica una copia acostada de un vector por una copia parada y colapsa a un solo número.

Las dimensiones internas ( $k$ y $k$ ) deben coincidir; son sobre las que se suma.
Las dimensiones externas ( $m$ y $q$ ) dan la forma de la respuesta.
La entrada $(i,j)$ es la fila $i$ de $\mathbf{A}$ punto la columna $j$ de $\mathbf{B}$ .

Dos productos cargan toda la faena en la regresión. Con $\mathbf{X}$ de tamaño $n \times p$ , la transpuesta $\mathbf{X}'$ es $p \times n$ , así que

\mathbf{X}'\mathbf{X} \ \text{es}\ p \times p, \qquad \mathbf{X}'\mathbf{Y}\ \text{es}\ p \times 1.

En palabras: $\mathbf{X}'\mathbf{X}$ reúne toda suma de cuadrados y producto cruzado de las columnas de predictores en una tabla cuadrada pequeña, y $\mathbf{X}'\mathbf{Y}$ reúne toda suma de predictor por respuesta en una columna corta, la materia prima de las ecuaciones normales. Una regla que usaremos constantemente es que transponer un producto invierte el orden de sus factores.

Demostración. Mostramos $(\mathbf{A}\mathbf{B})' = \mathbf{B}'\mathbf{A}'$ comparando entradas. La entrada $(i,j)$ de $(\mathbf{A}\mathbf{B})'$ es, por la definición de transpuesta, la entrada $(j,i)$ de $\mathbf{A}\mathbf{B}$ , que es $\sum_{\ell} A_{j\ell} B_{\ell i}$ . La entrada $(i,j)$ de $\mathbf{B}'\mathbf{A}'$ es la fila $i$ de $\mathbf{B}'$ punto la columna $j$ de $\mathbf{A}'$ , es decir $\sum_{\ell} (\mathbf{B}')_{i\ell} (\mathbf{A}')_{\ell j} = \sum_{\ell} B_{\ell i} A_{j \ell}$ . Las dos sumas tienen los mismos términos, así que las matrices son iguales. $\blacksquare$

R y Python¶

Primero el ejemplo pequeño de la Figure 3, para que la mecánica sea concreta antes de soltar el motor sobre los datos reales. En R, %*% es la multiplicación de matrices (un * simple multiplica entrada por entrada); en Python el operador es @.

A <- matrix(c(1, 2, 3,
              4, 5, 6), nrow = 2, byrow = TRUE)
B <- matrix(c(1, 0,
              0, 1,
              1, 1), nrow = 3, byrow = TRUE)
A %*% B

     [,1] [,2]
[1,]    4    5
[2,]   10   11

A = np.array([[1, 2, 3],
              [4, 5, 6]])
B = np.array([[1, 0],
              [0, 1],
              [1, 1]])
print(A @ B)

[[ 4  5]
 [10 11]]

Ejemplo 6.2: Las matrices

\mathbf{X}'\mathbf{X}

\mathbf{X}'\mathbf{Y}

para Dwaine

Pregunta. Forma $\mathbf{X}'\mathbf{X}$ y $\mathbf{X}'\mathbf{Y}$ para los datos de Dwaine. ¿Qué significan las entradas individuales?

Intuición. $\mathbf{X}'\mathbf{X}$ es una tabla $3 \times 3$ cuyas entradas son sumas de productos de las columnas de predictores; $\mathbf{X}'\mathbf{Y}$ es una columna de longitud 3 que empareja cada columna de predictor con las ventas.

Fórmula. $(\mathbf{X}'\mathbf{X})_{jk} = \sum_i X_{ij} X_{ik}$ y $(\mathbf{X}'\mathbf{Y})_j = \sum_i X_{ij} Y_i$ .

Cálculo.

XtX <- t(X) %*% X
XtX

          intercept targtpop dispoinc
intercept      21.0  1302.40   360.00
targtpop     1302.4 87707.94 22609.19
dispoinc      360.0 22609.19  6190.26

XtY <- t(X) %*% Y
XtY

               [,1]
intercept   3820.00
targtpop  249643.35
dispoinc   66072.75

XtX = X.T @ X
print(XtX)

[[2.100000e+01 1.302400e+03 3.600000e+02]
 [1.302400e+03 8.770794e+04 2.260919e+04]
 [3.600000e+02 2.260919e+04 6.190260e+03]]

XtY = X.T @ Y
print(XtY)

[  3820.   249643.35  66072.75]

Interpretación. La entrada superior izquierda de $\mathbf{X}'\mathbf{X}$ es $21 = \sum 1 \cdot 1 = n$ , el tamaño de muestra. El resto de la primera fila y columna son las sumas de predictores $\sum \text{targtpop} = 1302.4$ y $\sum \text{dispoinc} = 360.0$ . Las entradas de la diagonal 87707.94 y 6190.26 son las sumas de los valores de predictor al cuadrado, y la fuera de diagonal 22609.19 es el producto cruzado $\sum (\text{targtpop})(\text{dispoinc})$ . En $\mathbf{X}'\mathbf{Y}$ , la primera entrada 3820 es $\sum Y$ , y las otras dos son $\sum(\text{targtpop})Y$ y $\sum(\text{dispoinc})Y$ . El Capítulo 7 reutiliza esta misma $\mathbf{X}'\mathbf{X}$ (7.1 El modelo y los mínimos cuadrados en forma matricial).

La regla fila por columna se vuelve creíble cuando la ves ocurrir sobre el único producto que una regresión realmente necesita, así que construye $\mathbf{X}'\mathbf{X}$ tú mismo antes de seguir.

Toca o desliza el dedo sobre cualquier entrada del producto y se resaltarán la fila y la columna que la generaron, con la aritmética escrita debajo.

6.3 Identidad, simetría, independencia y rango¶

Intuición¶

Antes de invertir nada, necesitamos vocabulario para las formas de las matrices que aparecen en la regresión, y una idea que decide si el ajuste es siquiera posible. La matriz identidad $\mathbf{I}$ (Definición 6.6) es la versión matricial del número 1: multiplicar por ella no cambia nada. Una matriz simétrica (Definición 6.7) es igual a su propia transpuesta, y $\mathbf{X}'\mathbf{X}$ siempre es simétrica, razón por la cual su tabla se ve como un espejo a través de la diagonal.

La idea que decide la factibilidad es la independencia lineal (Definición 6.8) de las columnas de predictores. Si una columna de predictor es una combinación lineal exacta de las otras, digamos que dispoinc fuera exactamente el doble de targtpop más una constante, entonces las columnas cargan información redundante, la matriz $\mathbf{X}$ tiene rango deficiente, y $\mathbf{X}'\mathbf{X}$ no se puede invertir. No hay entonces un único mejor ajuste, porque los datos no pueden distinguir los coeficientes redundantes. Esta es la cara matricial de un problema que volverás a encontrar como multicolinealidad en 12.1 La multicolinealidad y el factor de inflación de la varianza.

Una imagen vuelve concreta la independencia. Piensa en dos columnas como flechas dibujadas desde el origen. Si apuntan en direcciones genuinamente distintas, abren un área entre ellas, y la matriz tiene rango completo. Si una es solo una copia estirada de la otra, ambas flechas quedan sobre una sola recta, el área entre ellas es cero, y la matriz tiene rango deficiente. La Figure 5 muestra ambos casos lado a lado.

Dos paneles lado a lado, cada uno mostrando dos flechas dibujadas desde el origen. En el panel izquierdo, rotulado columnas independientes encierran un área, una flecha azul y una flecha verde apuntan en direcciones distintas y el paralelogramo entre ellas está sombreado, marcado área mayor que cero, rango completo, determinante distinto de cero, ajuste único. En el panel derecho, rotulado columnas dependientes colapsan a una recta, la flecha verde es una copia estirada de la flecha azul de modo que ambas quedan sobre una recta punteada, marcado área igual a cero, rango deficiente, determinante cero, sin ajuste único. — Figure 5:El rango visto como geometría. Las columnas independientes se abren y encierran un área, así que el ajuste es único; las columnas dependientes caen sobre una recta y no encierran nada, así que el ajuste no lo es. El área encerrada es exactamente el determinante de la próxima sección, y su colapso a cero es lo que vuelve singular a una matriz.

Estos nombres se ganan su lugar: una propiedad con nombre es un hecho que puedes citar en vez de recalcular. Mostrar que la matriz sombrero es una proyección se vuelve un argumento de una línea, “es simétrica e idempotente”, en lugar de páginas de comprobación entrada por entrada.

Fórmula¶

El producto $\mathbf{X}'\mathbf{X}$ siempre es simétrico, porque $(\mathbf{X}'\mathbf{X})' = \mathbf{X}'(\mathbf{X}')' = \mathbf{X}'\mathbf{X}$ por la regla de inversión del orden (Teorema 6.5) de 6.2 Transpuesta y multiplicación de matrices.

El hecho central de esta sección conecta el rango con la invertibilidad.

Un ejemplo diminuto vuelve concreto el rango: $(1, 0)'$ y $(0, 1)'$ son independientes, pero $(1, 2)'$ y $(2, 4)'$ son dependientes porque el segundo es exactamente el doble del primero, así que juntos abarcan solo una recta y tienen rango 1, no 2. En una matriz de diseño, una columna dependiente así es un predictor que repite información ya presente.

R y Python¶

El software reporta la simetría, el rango y el determinante (el único número, cubierto en 6.4 La inversa y las ecuaciones normales, cuya anulación señala el rango deficiente) directamente.

isSymmetric(XtX)
qr(X)$rank
round(det(XtX), 2)

[1] TRUE
[1] 3
[1] 1068302

print(np.allclose(XtX, XtX.T))
print(np.linalg.matrix_rank(X))
print(round(np.linalg.det(XtX), 2))

True
3
1068302.4

La matriz de diseño tiene rango 3, su cuenta completa de columnas, así que los dos predictores más el intercepto cargan tres piezas de información genuinamente distintas, y el determinante no nulo confirma que $\mathbf{X}'\mathbf{X}$ es invertible y que el ajuste de Dwaine está bien planteado.

El hábito por construir: antes de confiar en cualquier regresión múltiple, pregúntate si los predictores son genuinamente distintos. El software con gusto invertirá una matriz casi singular y devolverá coeficientes con errores estándar enormes, la advertencia temprana de que dos predictores cargan casi la misma información. El rango es la versión de sí o no de esa pregunta; el Capítulo 12 la afina en “¿qué tan cerca de dependientes están?”.

6.4 La inversa y las ecuaciones normales¶

Intuición¶

Aquí es donde la maquinaria rinde frutos. En el álgebra ordinaria, para resolver $ax = c$ divides entre $a$ , lo cual es lo mismo que multiplicar por $a^{-1}$ . Las matrices no tienen división, pero tienen lo mejor que le sigue: para una matriz cuadrada $\mathbf{A}$ de rango completo, existe una inversa $\mathbf{A}^{-1}$ (Definición 6.10) que la deshace, es decir, $\mathbf{A}^{-1}\mathbf{A} = \mathbf{I}$ . Para resolver una ecuación matricial $\mathbf{A}\mathbf{x} = \mathbf{c}$ , multiplica ambos lados por la izquierda por $\mathbf{A}^{-1}$ y obtienes $\mathbf{x} = \mathbf{A}^{-1}\mathbf{c}$ .

Las ecuaciones normales de los mínimos cuadrados, que en el Capítulo 2 eran dos ecuaciones escalares, colapsan en forma matricial a la única ecuación $\mathbf{X}'\mathbf{X}\,\mathbf{b} = \mathbf{X}'\mathbf{Y}$ . El Capítulo 7 la deriva minimizando el error al cuadrado (7.1 El modelo y los mínimos cuadrados en forma matricial); por ahora tómala como la gemela matricial de las ecuaciones normales del Capítulo 2 (2.2 Mínimos cuadrados desde los primeros principios). Resolverla es ahora un solo paso: multiplica por $(\mathbf{X}'\mathbf{X})^{-1}$ .

Con números, $a^{-1}$ existe para todo $a$ excepto el cero. Con matrices, $\mathbf{A}^{-1}$ existe para toda matriz cuadrada excepto las de rango deficiente, así que “¿puedo resolver esto de manera única?”, “¿el determinante es distinto de cero?” y “¿son independientes las columnas?” son tres maneras de hacer una sola pregunta. Para Dwaine, las tres dicen que sí.

Fórmula¶

El determinante $ad - bc$ mide si las columnas son independientes; si es cero, la fórmula de la inversa divide entre cero y no existe inversa.
Para matrices más grandes el determinante es más complicado, así que dejamos que el software lo calcule, pero la regla “invertible si y solo si el determinante es distinto de cero” sigue siendo exacta.

Esto conecta de vuelta con la imagen del rango en Figure 5: las dos columnas de una matriz dos por dos enmarcan un paralelogramo, y el determinante es su área. Cuando las columnas apuntan a lo largo de la misma recta el paralelogramo queda aplastado, su área es cero, y la matriz es singular, así que “determinante cero” y “las columnas colapsan sobre una recta” son el mismo suceso.

Aplicar la inversa a las ecuaciones normales da el estimador de mínimos cuadrados en una sola línea.

En palabras: multiplica los productos cruzados almacenados $\mathbf{X}'\mathbf{Y}$ por la inversa de la matriz de productos cruzados $\mathbf{X}'\mathbf{X}$ , y salen las $p$ estimaciones de coeficientes de una vez. Esta única fórmula es toda la estimación por mínimos cuadrados, para cualquier número de predictores.

R y Python¶

Primero la fórmula de la inversa $2 \times 2$ , comprobada en una matriz simétrica pequeña, para que puedas ver ocurrir $\mathbf{A}^{-1}\mathbf{A} = \mathbf{I}$ . Ambos lenguajes invierten con una sola llamada: solve en R, np.linalg.inv en Python.

M <- matrix(c(2, 1,
              1, 2), nrow = 2, byrow = TRUE)
Minv <- solve(M)
Minv
round(M %*% Minv, 6)

           [,1]       [,2]
[1,]  0.6666667 -0.3333333
[2,] -0.3333333  0.6666667
     [,1] [,2]
[1,]    1    0
[2,]    0    1

M = np.array([[2.0, 1.0],
              [1.0, 2.0]])
Minv = np.linalg.inv(M)
print(Minv)
print(np.round(M @ Minv, 6))

[[ 0.66666667 -0.33333333]
 [-0.33333333  0.66666667]]
[[1. 0.]
 [0. 1.]]

El determinante es $2\cdot 2 - 1 \cdot 1 = 3$ , así que la inversa es $\tfrac{1}{3}\left(\begin{smallmatrix} 2 & -1 \\ -1 & 2 \end{smallmatrix}\right)$ , que concuerda con la impresión, y $\mathbf{M}^{-1}\mathbf{M}$ devuelve la identidad.

Casi nunca invertirás a mano una matriz más grande que esta; trabajar el caso dos por dos una vez basta para ver de dónde viene la respuesta y observar aparecer la identidad, de modo que una inversa por software nunca es magia. De aquí en adelante dejamos que solve y np.linalg.inv hagan la aritmética, y el próximo ejemplo suelta la inversa sobre el ajuste completo de Dwaine.

Ejemplo 6.3: Resolver las ecuaciones normales para Dwaine

Pregunta. Invierte $\mathbf{X}'\mathbf{X}$ y calcula el vector de coeficientes $\mathbf{b} = (\mathbf{X}'\mathbf{X})^{-1}\mathbf{X}'\mathbf{Y}$ . ¿Coinciden los números con lo que reportan lm y statsmodels?

Intuición. Ya almacenamos $\mathbf{X}'\mathbf{X}$ y $\mathbf{X}'\mathbf{Y}$ . Invertir la primera y multiplicar por la segunda resuelve los tres coeficientes de una vez.

Fórmula. $\mathbf{b} = (\mathbf{X}'\mathbf{X})^{-1}\mathbf{X}'\mathbf{Y}$ .

Cálculo.

XtX_inv <- solve(XtX)
round(XtX_inv, 6)
round(XtX %*% XtX_inv, 6)

          intercept  targtpop  dispoinc
intercept 29.728923  0.072183 -1.992553
targtpop   0.072183  0.000370 -0.005550
dispoinc  -1.992553 -0.005550  0.136311
          intercept targtpop dispoinc
intercept         1        0        0
targtpop          0        1        0
dispoinc          0        0        1

b <- XtX_inv %*% XtY
b

               [,1]
intercept -68.85707
targtpop    1.45456
dispoinc    9.36550

fit <- lm(sales ~ targtpop + dispoinc, data = dwaine)
coef(fit)

(Intercept)    targtpop    dispoinc
  -68.85707     1.45456     9.36550

Ahora el mismo cálculo en Python.

XtX_inv = np.linalg.inv(XtX)
print(np.round(XtX_inv, 6))
print(np.round(XtX @ XtX_inv, 6))

[[ 2.9728923e+01  7.2183000e-02 -1.9925530e+00]
 [ 7.2183000e-02  3.7000000e-04 -5.5500000e-03]
 [-1.9925530e+00 -5.5500000e-03  1.3631100e-01]]
[[ 1. -0. -0.]
 [ 0.  1. -0.]
 [-0. -0.  1.]]

b = XtX_inv @ XtY
print(b)

[-68.85707315   1.45455958   9.36550038]

fit = smf.ols("sales ~ targtpop + dispoinc", data=dwaine).fit()
print(fit.params)

Intercept   -68.857073
targtpop      1.454560
dispoinc      9.365500
dtype: float64

Interpretación. El estimador construido a mano da $\mathbf{b} = (-68.86,\ 1.4546,\ 9.3655)'$ , y tanto lm como statsmodels devuelven los números idénticos, lo que confirma que las rutinas internas hacen exactamente este cálculo matricial. El modelo ajustado es $\widehat{\text{sales}} = -68.86 + 1.4546\,(\text{targtpop}) + 9.3655\,(\text{dispoinc})$ . Cada mil residentes jóvenes adicionales suman cerca de 1.45 mil dólares de ventas predichas, manteniendo fijo el ingreso, y cada mil dólares adicionales de ingreso suman cerca de 9.37 mil, manteniendo fija la población. El Capítulo 8 desarrolla esa lectura de “manteniendo fijo” (8.3 Sumas de cuadrados adicionales); aquí, una sola inversa de matriz produjo todo el ajuste.

Inténtalo 6.4

Usando la fórmula $2 \times 2$ , invierte a mano $\mathbf{A} = \left(\begin{smallmatrix} 4 & 2 \\ 1 & 3 \end{smallmatrix}\right)$ . Luego explica por qué $\mathbf{B} = \left(\begin{smallmatrix} 2 & 4 \\ 1 & 2 \end{smallmatrix}\right)$ no tiene inversa.

Solución

Para $\mathbf{A}$ , el determinante es $4 \cdot 3 - 2 \cdot 1 = 10$ , así que $\mathbf{A}^{-1} = \tfrac{1}{10}\left(\begin{smallmatrix} 3 & -2 \\ -1 & 4 \end{smallmatrix}\right) = \left(\begin{smallmatrix} 0.3 & -0.2 \\ -0.1 & 0.4 \end{smallmatrix}\right)$ . Para $\mathbf{B}$ , el determinante es $2 \cdot 2 - 4 \cdot 1 = 0$ , así que no existe inversa. La razón es visible en las columnas: la columna 2 es exactamente el doble de la columna 1, así que las columnas son linealmente dependientes y $\mathbf{B}$ tiene rango deficiente.

Un determinante igual a cero se siente mejor de lo que se define, así que junta tú mismo las dos columnas y observa cómo la inversa se descompone mucho antes de que lleguen a tocarse.

Arrastra la punta de cualquiera de las dos flechas. La caja sombreada tiene área igual al tamaño del determinante, y las entradas de la inversa crecen sin límite conforme esa caja se aplana.

6.5 Matrices idempotentes y de proyección¶

Intuición¶

Los coeficientes están listos, así que el modelo ya puede hacer sus propias predicciones. La pregunta de esta sección es de dónde vienen esas predicciones y cómo lucen como imagen. La respuesta es una sola tabla de números que convierte las ventas observadas en las ventas ajustadas del modelo en un solo paso, y ese paso es una sombra: aplana los datos sobre la superficie que los predictores pueden alcanzar.

Con $\mathbf{b}$ en mano, los valores ajustados son $\widehat{\mathbf{Y}} = \mathbf{X}\mathbf{b}$ . Sustituye el estimador y aparece algo notable:

\widehat{\mathbf{Y}} = \mathbf{X}\mathbf{b} = \mathbf{X}(\mathbf{X}'\mathbf{X})^{-1}\mathbf{X}'\mathbf{Y} = \mathbf{H}\mathbf{Y}, \qquad \mathbf{H} = \mathbf{X}(\mathbf{X}'\mathbf{X})^{-1}\mathbf{X}'.

Geométricamente la matriz sombrero es una proyección (Definición 6.14). Imagina $\mathbf{Y}$ como una flecha en un espacio de $n$ dimensiones; los valores ajustados que una regresión puede producir quedan todos en un subespacio plano, el espacio generado por las columnas de $\mathbf{X}$ . La matriz sombrero deja caer $\mathbf{Y}$ recto hacia abajo sobre ese subespacio, aterrizando en el punto más cercano, y el residuo es lo que sobra en perpendicular. La Figure 7 muestra esta imagen, el tema de 7.2 La geometría de los mínimos cuadrados.

Un boceto tridimensional. Un plano sombreado rotulado espacio columna de X se sitúa cerca del piso. Una flecha azul rotulada Y apunta hacia arriba desde el origen por encima del plano. Una flecha naranja rotulada Y-sombrero igual a H Y queda en el plano, justo debajo de la punta de Y. Una flecha roja rotulada e igual a I menos H por Y conecta la punta de Y-sombrero hacia arriba con la punta de Y, encontrando el plano en ángulo recto marcado con un pequeño cuadrado. — Figure 7:Los mínimos cuadrados como proyección. El vector ajustado Y-sombrero es la sombra de la respuesta Y sobre el espacio columna de X, y el residuo e los une en ángulo recto. La matriz sombrero H realiza la proyección sobre el plano; I menos H produce el residuo perpendicular.

Dos propiedades algebraicas hacen de $\mathbf{H}$ una proyección, y son las bestias de carga de los diagnósticos del Capítulo 9. Primero, $\mathbf{H}$ es simétrica. Segundo, $\mathbf{H}$ es idempotente: aplicarla dos veces es lo mismo que aplicarla una vez, $\mathbf{H}\mathbf{H} = \mathbf{H}$ . Eso tiene sentido para una proyección, pues una vez que un vector queda aplanado sobre el plano, aplanarlo de nuevo no hace nada.

La imagen de la sombra vale la pena retenerla. Tu sombra sobre suelo plano conserva tu posición izquierda-derecha y adelante-atrás pero aplana tu altura a cero; la matriz sombrero le hace lo mismo a $\mathbf{Y}$ , aplanando la parte de la respuesta que ninguna combinación de los predictores podría reproducir, que es el residuo.

Estas propiedades no son adorno: el Capítulo 9 lee la diagonal de $\mathbf{H}$ para encontrar qué ciudades jalan más fuerte sobre la superficie ajustada, el Capítulo 7 usa su traza para contar los grados de libertad, y la imagen perpendicular es lo que hace de los mínimos cuadrados el ajuste más cercano, sin ningún punto en el plano más próximo a $\mathbf{Y}$ que su propia sombra.

Fórmula¶

Las dos proyecciones de la regresión son

\mathbf{H} = \mathbf{X}(\mathbf{X}'\mathbf{X})^{-1}\mathbf{X}', \qquad \mathbf{I} - \mathbf{H},

con $\widehat{\mathbf{Y}} = \mathbf{H}\mathbf{Y}$ los valores ajustados y $\mathbf{e} = \mathbf{Y} - \widehat{\mathbf{Y}} = (\mathbf{I} - \mathbf{H})\mathbf{Y}$ los residuos. En palabras: $\mathbf{H}$ proyecta sobre el espacio de valores ajustables, y $\mathbf{I} - \mathbf{H}$ proyecta sobre el espacio sobrante de los residuos.

Las dos propiedades que definen a $\mathbf{H}$ son las que la vuelven una proyección.

Demostración. Escribe $\mathbf{G} = (\mathbf{X}'\mathbf{X})^{-1}$ , que es simétrica porque $\mathbf{X}'\mathbf{X}$ lo es (su inversa hereda la simetría: $\mathbf{G}' = ((\mathbf{X}'\mathbf{X})^{-1})' = ((\mathbf{X}'\mathbf{X})')^{-1} = (\mathbf{X}'\mathbf{X})^{-1} = \mathbf{G}$ ). Para la simetría de $\mathbf{H}$ , aplica la regla de la transpuesta con inversión del orden de 6.2 Transpuesta y multiplicación de matrices dos veces:

\mathbf{H}' = (\mathbf{X}\mathbf{G}\mathbf{X}')' = (\mathbf{X}')'\mathbf{G}'\mathbf{X}' = \mathbf{X}\mathbf{G}\mathbf{X}' = \mathbf{H}.

Para la idempotencia, multiplica $\mathbf{H}$ por sí misma y cancela la $\mathbf{X}'\mathbf{X}$ interna contra su inversa:

\mathbf{H}\mathbf{H} = \mathbf{X}\mathbf{G}(\mathbf{X}'\mathbf{X})\mathbf{G}\mathbf{X}' = \mathbf{X}\mathbf{G}\mathbf{X}' = \mathbf{H},

usando $(\mathbf{X}'\mathbf{X})\mathbf{G} = (\mathbf{X}'\mathbf{X})(\mathbf{X}'\mathbf{X})^{-1} = \mathbf{I}$ . Los mismos dos pasos muestran que $\mathbf{I} - \mathbf{H}$ es simétrica e idempotente: $(\mathbf{I} - \mathbf{H})(\mathbf{I} - \mathbf{H}) = \mathbf{I} - 2\mathbf{H} + \mathbf{H}\mathbf{H} = \mathbf{I} - 2\mathbf{H} + \mathbf{H} = \mathbf{I} - \mathbf{H}$ . $\blacksquare$

R y Python¶

Ejemplo 6.4: La matriz sombrero en acción

Pregunta. Construye $\mathbf{H}$ para Dwaine, confirma que es simétrica e idempotente con traza 3, y úsala para producir los valores ajustados y los residuos.

Intuición. $\mathbf{H}$ es $21 \times 21$ . Comprobamos numéricamente las dos propiedades que la definen, y luego aplicamos $\mathbf{H}$ y $\mathbf{I} - \mathbf{H}$ a $\mathbf{Y}$ .

Fórmula. $\mathbf{H} = \mathbf{X}(\mathbf{X}'\mathbf{X})^{-1}\mathbf{X}'$ , $\widehat{\mathbf{Y}} = \mathbf{H}\mathbf{Y}$ , $\mathbf{e} = (\mathbf{I} - \mathbf{H})\mathbf{Y}$ .

Cálculo.

H <- X %*% XtX_inv %*% t(X)
c(symmetric = isSymmetric(round(H, 10)),
  idempotent = all(abs(H %*% H - H) < 1e-10))
round(sum(diag(H)), 6)

 symmetric idempotent
      TRUE       TRUE
[1] 3

Yhat <- H %*% Y
e <- (diag(n) - H) %*% Y
round(head(cbind(Y, Yhat, e), 3), 3)

         Y
[1,] 174.4 187.184 -12.784
[2,] 164.4 154.229  10.171
[3,] 244.2 234.396   9.804

H = X @ XtX_inv @ X.T
print(np.allclose(H, H.T), np.allclose(H @ H, H))
print(round(np.trace(H), 6))

True True
3.0

Yhat = H @ Y
e = (np.eye(n) - H) @ Y
print(np.round(np.column_stack([Y, Yhat, e])[:3], 3))

[[174.4   187.184 -12.784]
 [164.4   154.229  10.171]
 [244.2   234.396   9.804]]

Interpretación. La matriz sombrero es simétrica e idempotente, con traza $3 = p$ , exactamente como prometió la derivación. Al aplicarla, las ventas ajustadas de la ciudad 1 son 187.2 contra 174.4 reales, un residuo de -12.8 (el modelo sobrepredice ligeramente). Las entradas de la diagonal de $\mathbf{H}$ , no impresas aquí, son los valores de apalancamiento $h_{ii}$ que el Capítulo 9 convierte en diagnósticos de influencia: uno grande marca una ciudad cuya respuesta tiene una voz desmedida en el ajuste. Todo análisis de residuos del libro parte de esta única matriz.

Inténtalo 6.5

Una matriz de proyección $\mathbf{P}$ satisface $\mathbf{P}\mathbf{P} = \mathbf{P}$ . Muestra que si $\mathbf{P}$ es idempotente entonces también lo es $\mathbf{I} - \mathbf{P}$ , y calcula $(\mathbf{I} - \mathbf{P})\mathbf{P}$ .

Solución

Desarrolla $(\mathbf{I} - \mathbf{P})(\mathbf{I} - \mathbf{P}) = \mathbf{I} - \mathbf{P} - \mathbf{P} + \mathbf{P}\mathbf{P} = \mathbf{I} - 2\mathbf{P} + \mathbf{P} = \mathbf{I} - \mathbf{P}$ , así que $\mathbf{I} - \mathbf{P}$ es idempotente. Y $(\mathbf{I} - \mathbf{P})\mathbf{P} = \mathbf{P} - \mathbf{P}\mathbf{P} = \mathbf{P} - \mathbf{P} = \mathbf{0}$ . Las dos proyecciones se aniquilan mutuamente, que es el álgebra detrás de que los valores ajustados y los residuos sean ortogonales.

La proyección es la única idea de este capítulo que una página plana no puede mostrar de verdad, así que gira la figura y mueve $\mathbf{Y}$ con tu propio dedo.

Arrastra la punta de Y a cualquier lugar del espacio. Yhat permanece sobre el espacio columna de X y el residuo e permanece perpendicular a él, que es justo lo que hace del ajuste el más cercano posible.

6.6 Formas cuadráticas y sumas de cuadrados¶

Intuición¶

Toda “suma de cuadrados” en la regresión es en secreto una expresión matricial llamada forma cuadrática (Definición 6.17). La suma de cuadrados del error es $\mathrm{SSE} = \sum e_i^2$ , y una suma de las entradas al cuadrado de un vector es exactamente ese vector punto consigo mismo: $\mathrm{SSE} = \mathbf{e}'\mathbf{e}$ . Como $\mathbf{e} = (\mathbf{I} - \mathbf{H})\mathbf{Y}$ , un poco de álgebra reescribe SSE usando solo $\mathbf{Y}$ y una matriz de proyección en el medio. Ese patrón de matriz en el medio, $\mathbf{Y}'\mathbf{A}\mathbf{Y}$ , es una forma cuadrática, y SSTO, SSR y SSE tienen todos esta forma.

¿Por qué importa? Por dos razones. Primero, una forma cuadrática hace que los grados de libertad y el valor esperado de cada suma de cuadrados salgan de la matriz del medio, que es como el Capítulo 7 explica la tabla ANOVA. Segundo, las sumas de cuadrados nunca pueden ser negativas, y eso lo garantiza una propiedad de la matriz del medio llamada semidefinición positiva. La Figure 9 muestra por qué una forma definida positiva tiene la forma de un tazón que nunca baja de cero.

Un gráfico de superficie tridimensional de la forma cuadrática u-prima M u para una matriz M dos por dos definida positiva. La superficie es un tazón suave hacia arriba que toca el cero solo en el origen, marcado con un punto rojo rotulado mínimo cero en u igual a cero, y que se eleva en toda dirección alejándose de él. — Figure 9:Una forma cuadrática definida positiva es un tazón. Su valor es cero solo en el origen y estrictamente positivo en toda otra dirección, que es exactamente por qué una suma de cuadrados escrita como tal forma nunca puede ser negativa.

La recompensa es que esta vista sobrevive incluso cuando la suma ya no parece una suma de cuadrados. Escrita como $\mathbf{Y}'(\mathbf{I} - \mathbf{H})\mathbf{Y}$ , SSE no eleva visiblemente nada al cuadrado, y sin embargo la matriz del medio todavía garantiza que no puede volverse negativa. El mismo razonamiento se aplica a una varianza, una suma de cuadrados disfrazada, así que la matriz de covarianzas de la próxima sección nunca puede contener una varianza negativa.

Fórmula¶

Las tres sumas de cuadrados de la regresión son todas formas cuadráticas en $\mathbf{Y}$ , y sus matrices del medio cargan la descomposición ANOVA.

Teorema 6.18: Sumas de cuadrados como formas cuadráticas

Con $\mathbf{J}$ la matriz $n \times n$ de puros unos,

\mathrm{SSE} = \mathbf{Y}'(\mathbf{I} - \mathbf{H})\mathbf{Y}, \qquad \mathrm{SSR} = \mathbf{Y}'\!\left(\mathbf{H} - \tfrac{1}{n}\mathbf{J}\right)\!\mathbf{Y}, \qquad \mathrm{SSTO} = \mathbf{Y}'\!\left(\mathbf{I} - \tfrac{1}{n}\mathbf{J}\right)\!\mathbf{Y},

cada una con una matriz del medio simétrica e idempotente. Las matrices del medio se suman, $(\mathbf{I} - \mathbf{H}) + (\mathbf{H} - \tfrac{1}{n}\mathbf{J}) = \mathbf{I} - \tfrac{1}{n}\mathbf{J}$ , que es el enunciado matricial de $\mathrm{SSTO} = \mathrm{SSR} + \mathrm{SSE}$ , y $\mathrm{SSE} \ge 0$ porque su matriz del medio es semidefinida positiva.

Demostración (no negatividad). Como $\mathbf{I} - \mathbf{H}$ es simétrica e idempotente, $\mathrm{SSE} = \mathbf{Y}'(\mathbf{I} - \mathbf{H})\mathbf{Y} = \mathbf{Y}'(\mathbf{I} - \mathbf{H})'(\mathbf{I} - \mathbf{H})\mathbf{Y} = \mathbf{e}'\mathbf{e} \ge 0$ , así que SSE es una suma de cuadrados genuina y nunca puede ser negativa; lo mismo vale para SSR y SSTO. La descomposición aditiva de las matrices del medio es el enunciado matricial de la identidad ANOVA $\mathrm{SSTO} = \mathrm{SSR} + \mathrm{SSE}$ (3.6 El enfoque del análisis de varianza). $\blacksquare$

En palabras: las sumas de cuadrados son formas cuadráticas, y las matrices del medio idempotentes garantizan la no negatividad y, en el Capítulo 7, nos entregan los grados de libertad a través de sus trazas.

R y Python¶

Ejemplo 6.5: SSE como forma cuadrática

Pregunta. Calcula SSE para Dwaine como la forma cuadrática $\mathbf{Y}'(\mathbf{I} - \mathbf{H})\mathbf{Y}$ , y úsala para estimar la varianza del error $s^2 = \mathrm{MSE}$ .

Intuición. Ya tenemos $\mathbf{e} = (\mathbf{I} - \mathbf{H})\mathbf{Y}$ , así que $\mathrm{SSE} = \mathbf{e}'\mathbf{e}$ es un solo número. Dividir entre $n - p$ grados de libertad da la estimación de la varianza.

Fórmula. $\mathrm{SSE} = \mathbf{e}'\mathbf{e}$ , $\mathrm{MSE} = \mathrm{SSE}/(n - p)$ con $p = 3$ .

Cálculo.

SSE <- as.numeric(t(e) %*% e)
p <- ncol(X)
MSE <- SSE / (n - p)
round(c(SSE = SSE, df = n - p, MSE = MSE), 4)

      SSE        df       MSE
2180.9274   18.0000  121.1626

SSE = float(e.T @ e)
p = X.shape[1]
MSE = SSE / (n - p)
print(round(SSE, 4), n - p, round(MSE, 4))

2180.9274 18 121.1626

Interpretación. La forma cuadrática de los residuos da $\mathrm{SSE} = 2180.9$ sobre $n - p = 21 - 3 = 18$ grados de libertad, así que $\mathrm{MSE} = 121.2$ y el error estándar residual es $s = \sqrt{121.2} \approx 11.0$ mil dólares: aun conociendo la población joven y el ingreso de una ciudad, la predicción de ventas típicamente yerra por cerca de 11 mil dólares. El divisor 18 cuenta las observaciones que quedan después de gastar $p = 3$ en los coeficientes, la versión de regresión múltiple de la regla $n - 2$ del Capítulo 2; el Capítulo 7 lee este $n - p$ en la traza de la matriz sombrero y lo usa para probar que MSE es insesgado (7.3 La matriz sombrero).

6.7 Vectores aleatorios, esperanza y matrices de covarianzas¶

Intuición¶

Hasta ahora las matrices contenían números fijos. Pero $\mathbf{Y}$ es aleatorio: vuelve a correr las 21 ciudades bajo las mismas condiciones y las ventas saldrían un poco distintas, por los errores $\boldsymbol{\varepsilon}$ . Para describir un vector de variables aleatorias necesitamos dos cosas: un vector de medias y una tabla de varianzas y covarianzas.

El vector de medias es simplemente la esperanza aplicada entrada por entrada. La tabla es la matriz de covarianzas (Definición 6.19): su diagonal contiene la varianza de cada entrada, y sus fuera de diagonal contienen la covarianza entre pares de entradas. Para los errores de la regresión, esta tabla es simple: varianza constante $\sigma^2$ por la diagonal (todo error tiene la misma dispersión) y ceros fuera de ella (errores distintos están no correlacionados), así que $\operatorname{Cov}\{\boldsymbol{\varepsilon}\} = \sigma^2\mathbf{I}$ . Una regla de dos líneas para cómo viajan las medias y las covarianzas a través de un mapa lineal entrega entonces la matriz de covarianzas de $\mathbf{b}$ , la fuente de todo error estándar.

Fórmula¶

Definición 6.19: Vector aleatorio, vector de medias y matriz de covarianzas

Un vector aleatorio $\mathbf{W} = (W_1, \dots, W_m)'$ es un vector cuyas entradas son variables aleatorias. Su vector de medias y su matriz de covarianzas son

E\{\mathbf{W}\} = \begin{pmatrix} E\{W_1\} \\ \vdots \\ E\{W_m\} \end{pmatrix}, \qquad \operatorname{Cov}\{\mathbf{W}\} = \begin{pmatrix} \operatorname{Var}\{W_1\} & \cdots & \operatorname{Cov}\{W_1, W_m\} \\ \vdots & \ddots & \vdots \\ \operatorname{Cov}\{W_m, W_1\} & \cdots & \operatorname{Var}\{W_m\} \end{pmatrix},

con las varianzas en la diagonal (siempre no negativas) y una simetría fuera de la diagonal porque $\operatorname{Cov}\{W_i, W_j\} = \operatorname{Cov}\{W_j, W_i\}$ .

Las reglas de transformación de cómo viajan estas a través de un mapa lineal son el corazón de la sección.

En palabras: la esperanza pasa directo a través de un mapa lineal, y una matriz de covarianzas queda emparedada entre $\mathbf{A}$ y su transpuesta. Aplica esto al modelo de regresión, donde $E\{\mathbf{Y}\} = \mathbf{X}\boldsymbol{\beta}$ y $\operatorname{Cov}\{\mathbf{Y}\} = \sigma^2\mathbf{I}$ , y la media y la covarianza de $\mathbf{b}$ se siguen.

Demostración. El estimador es un mapa lineal de $\mathbf{Y}$ : con $\mathbf{A} = (\mathbf{X}'\mathbf{X})^{-1}\mathbf{X}'$ , tenemos $\mathbf{b} = \mathbf{A}\mathbf{Y}$ . Para la media, usa $E\{\mathbf{Y}\} = \mathbf{X}\boldsymbol{\beta}$ :

E\{\mathbf{b}\} = \mathbf{A}\,E\{\mathbf{Y}\} = (\mathbf{X}'\mathbf{X})^{-1}\mathbf{X}'\mathbf{X}\boldsymbol{\beta} = \boldsymbol{\beta},

así que $\mathbf{b}$ es insesgado para $\boldsymbol{\beta}$ . Para la covarianza, usa $\operatorname{Cov}\{\mathbf{Y}\} = \sigma^2\mathbf{I}$ y la regla del emparedado:

\operatorname{Cov}\{\mathbf{b}\} = \mathbf{A}\,(\sigma^2\mathbf{I})\,\mathbf{A}' = \sigma^2 (\mathbf{X}'\mathbf{X})^{-1}\mathbf{X}'\,\mathbf{X}(\mathbf{X}'\mathbf{X})^{-1} = \sigma^2 (\mathbf{X}'\mathbf{X})^{-1},

donde la $\mathbf{X}'\mathbf{X}$ del medio cancela una inversa. Así que $\operatorname{Cov}\{\mathbf{b}\} = \sigma^2 (\mathbf{X}'\mathbf{X})^{-1}$ : la misma inversa que usamos para calcular $\mathbf{b}$ también da, una vez escalada por $\sigma^2$ , la varianza de cada coeficiente y la covarianza de cada par. $\blacksquare$

Reemplazar el desconocido $\sigma^2$ por su estimación MSE da la matriz de covarianzas estimada $s^2\{\mathbf{b}\} = \mathrm{MSE}\,(\mathbf{X}'\mathbf{X})^{-1}$ , cuyas raíces cuadradas de la diagonal son los errores estándar $s\{b_k\}$ que el software imprime.

R y Python¶

Ejemplo 6.6: Errores estándar a partir de la matriz de covarianzas

Pregunta. Forma la matriz de covarianzas estimada $s^2\{\mathbf{b}\} = \mathrm{MSE}\,(\mathbf{X}'\mathbf{X})^{-1}$ para Dwaine, lee los errores estándar de los coeficientes, y compáralos con el resumen del software.

Intuición. Escala la inversa almacenada por MSE; la diagonal contiene las tres varianzas de los coeficientes, y sus raíces cuadradas son los errores estándar.

Fórmula. $s^2\{\mathbf{b}\} = \mathrm{MSE}\,(\mathbf{X}'\mathbf{X})^{-1}$ , y $s\{b_k\} = \sqrt{[s^2\{\mathbf{b}\}]_{kk}}$ .

Cálculo.

cov_b <- MSE * XtX_inv
round(cov_b, 4)
round(sqrt(diag(cov_b)), 4)

          intercept targtpop  dispoinc
intercept 3602.0347   8.7459 -241.4230
targtpop     8.7459   0.0449   -0.6724
dispoinc  -241.4230  -0.6724   16.5158
intercept  targtpop  dispoinc
  60.0170    0.2118    4.0640

round(summary(fit)$coefficients[, 1:2], 4)

            Estimate Std. Error
(Intercept) -68.8571    60.0170
targtpop      1.4546     0.2118
dispoinc      9.3655     4.0640

cov_b = MSE * XtX_inv
print(np.round(cov_b, 4))
print(np.round(np.sqrt(np.diag(cov_b)), 4))

[[ 3.6020347e+03  8.7459000e+00 -2.4142300e+02]
 [ 8.7459000e+00  4.4900000e-02 -6.7240000e-01]
 [-2.4142300e+02 -6.7240000e-01  1.6515800e+01]]
[60.017   0.2118  4.064 ]

print(fit.summary2().tables[1].iloc[:, 0:2].round(4))

             Coef.  Std.Err.
Intercept -68.8571   60.0170
targtpop    1.4546    0.2118
dispoinc    9.3655    4.0640

Interpretación. La diagonal de $s^2\{\mathbf{b}\}$ da las varianzas 3602.0, 0.0449 y 16.52, cuyas raíces cuadradas son los errores estándar 60.02, 0.2118 y 4.064, que concuerdan con ambos resúmenes del software en cada dígito impreso. La fuera de diagonal -0.6724 entre las dos pendientes dice que a lo largo de muestras repetidas una sobreestimación de una tiende a venir con una subestimación de la otra, la consecuencia de que los dos predictores comparten información. El Capítulo 7 usa esta matriz completa, diagonal y fuera de diagonal por igual, para construir cada intervalo de confianza y prueba conjunta.

Inténtalo 6.7

Un vector aleatorio $\mathbf{W}$ tiene matriz de covarianzas $\operatorname{Cov}\{\mathbf{W}\} = \left(\begin{smallmatrix} 4 & 1 \\ 1 & 9 \end{smallmatrix}\right)$ . Sea $\mathbf{A} = (1, 1)$ , una fila $1 \times 2$ , de modo que $\mathbf{A}\mathbf{W} = W_1 + W_2$ . Usa la regla de la covarianza para hallar $\operatorname{Var}\{W_1 + W_2\}$ .

Solución

Por la regla, $\operatorname{Var}\{\mathbf{A}\mathbf{W}\} = \mathbf{A}\operatorname{Cov}\{\mathbf{W}\}\mathbf{A}'$ . Calcula $\mathbf{A}\operatorname{Cov}\{\mathbf{W}\} = (1,1)\left(\begin{smallmatrix} 4 & 1 \\ 1 & 9 \end{smallmatrix}\right) = (5, 10)$ , luego $(5, 10)(1, 1)' = 15$ . Así que $\operatorname{Var}\{W_1 + W_2\} = 15$ . Esto concuerda con la regla escalar $\operatorname{Var}\{W_1\} + \operatorname{Var}\{W_2\} + 2\operatorname{Cov}\{W_1, W_2\} = 4 + 9 + 2(1) = 15$ .

6.8 La normal multivariante, en breve¶

Intuición¶

El vector de medias y la matriz de covarianzas describen el centro y la dispersión de un vector aleatorio, pero no su forma completa. Para la inferencia añadimos un supuesto, el mismo que convirtió las estimaciones del Capítulo 2 en pruebas $t$ : los errores son normales. Apilados en un vector, los errores normales con varianza constante y sin correlación siguen una distribución normal multivariante (Definición 6.22), la generalización a varias variables de la campana. Sus contornos son elipses cuya inclinación y estiramiento se leen directamente en la matriz de covarianzas, como en Figure 10.

Una dispersión de unos mil doscientos puntos que forman una nube elíptica inclinada centrada en el origen, con contornos de densidad elípticos rojos concéntricos dibujados encima. La nube se estira a lo largo de la diagonal de cuarenta y cinco grados porque las dos coordenadas tienen una covarianza positiva de 0.8. — Figure 10:Una distribución normal bivariante. Las dos coordenadas tienen covarianza 0.8, así que la nube y sus contornos elípticos se inclinan a lo largo de la diagonal. La forma de la elipse es exactamente la matriz de covarianzas, razón por la cual la distribución muestral del vector de coeficientes b se describe por su matriz de covarianzas.

Por qué importa esto: bajo errores normales la respuesta es normal multivariante, y como $\mathbf{b}$ es un mapa lineal de $\mathbf{Y}$ , también lo es $\mathbf{b}$ . Ese hecho, $\mathbf{b} \sim N(\boldsymbol{\beta}, \sigma^2(\mathbf{X}'\mathbf{X})^{-1})$ , es el fundamento de todo intervalo de confianza, prueba $t$ y prueba $F$ que vienen.

Lo que la distribución completa añade es cómo se adelgaza la probabilidad a medida que te alejas del centro: la misma caída de campana que conoces de una variable, medida a lo largo de los ejes inclinados de la elipse de covarianza. Eso es lo que nos permite calcular la probabilidad de que un coeficiente caiga dentro de una distancia declarada de la verdad, exactamente lo que reporta un intervalo de confianza en el próximo capítulo.

Una advertencia mantiene esto honesto. La normalidad multivariante de $\mathbf{b}$ es una consecuencia de suponer errores normales, no un hecho que los datos te entregan gratis. Si los errores son muy poco normales y la muestra es pequeña, esa ley es solo aproximada, una razón por la que existen los métodos de permutación y bootstrap del Capítulo 5 (5.4 El bootstrap para la regresión). Con una muestra grande, un efecto de promediado jala a $\mathbf{b}$ hacia la normalidad aun cuando los errores no lo sean.

Fórmula¶

Aplicar esa propiedad a la regresión bajo el modelo de errores normales da la ley muestral que impulsa toda inferencia posterior.

La respuesta está centrada en la superficie de regresión $\mathbf{X}\boldsymbol{\beta}$ con coordenadas independientes y de igual varianza, y el estimador $\mathbf{b}$ está centrado en la verdad $\boldsymbol{\beta}$ (insesgado) con la matriz de covarianzas que derivamos.

En palabras: los errores normales vuelven normales tanto a los datos como a las estimaciones, así que cada coeficiente está distribuido normalmente con una varianza que podemos leer en $(\mathbf{X}'\mathbf{X})^{-1}$ , y eso es lo que hace posible la inferencia exacta.

R y Python¶

Podemos demostrar la idea de recuperación de la covarianza por simulación: extrae muchos vectores con una covarianza objetivo y confirma que la covarianza muestral concuerda. Partiendo de normales estándar independientes $\mathbf{Z}$ , la transformación $\mathbf{V} = \mathbf{Z}\mathbf{L}'$ , donde $\mathbf{L}$ es una raíz cuadrada matricial del objetivo $\boldsymbol{\Sigma}$ (el factor de Cholesky), produce vectores con covarianza $\boldsymbol{\Sigma}$ , exactamente la regla del emparedado en acción.

set.seed(4210)
mu <- c(0, 0)
Sigma <- matrix(c(1, 0.8,
                  0.8, 1), nrow = 2)
L <- chol(Sigma)
Z <- matrix(rnorm(2 * 5000), ncol = 2)
V <- Z %*% L
round(cov(V), 3)

      [,1]  [,2]
[1,] 0.967 0.760
[2,] 0.760 0.955

rng = np.random.default_rng(4210)
Sigma = np.array([[1.0, 0.8],
                  [0.8, 1.0]])
L = np.linalg.cholesky(Sigma)
Z = rng.standard_normal((5000, 2))
V = Z @ L.T
print(np.round(np.cov(V, rowvar=False), 3))

[[1.017 0.822]
 [0.822 1.026]]

La covarianza muestral de los 5000 vectores simulados sale cerca del objetivo $\left(\begin{smallmatrix} 1 & 0.8 \\ 0.8 & 1 \end{smallmatrix}\right)$ en ambos lenguajes, difiriendo solo por el ruido de muestreo.

El factor de Cholesky puede seguir siendo una caja negra por ahora; el punto es que una covarianza se puede plantar a propósito, pasando ruido no correlacionado por una matriz fija y dejando que la regla del emparedado haga el resto. La covarianza de $\mathbf{b}$ es la misma historia: ruido de error $\sigma^2\mathbf{I}$ pasado por la matriz fija $(\mathbf{X}'\mathbf{X})^{-1}\mathbf{X}'$ , así que la simulación y la regresión son dos usos de un mismo mecanismo.

Inténtalo 6.8

Bajo el modelo de errores normales, $\mathbf{b} \sim N(\boldsymbol{\beta}, \sigma^2(\mathbf{X}'\mathbf{X})^{-1})$ . Para el ajuste de Dwaine, ¿qué única matriz almacenada, una vez escalada, da la varianza de la pendiente de dispoinc? ¿Dónde está esa varianza en la impresión del Ejemplo 6.6?

Solución

La matriz de covarianzas de $\mathbf{b}$ es $\sigma^2(\mathbf{X}'\mathbf{X})^{-1}$ , estimada por $\mathrm{MSE}\,(\mathbf{X}'\mathbf{X})^{-1}$ . La varianza de la pendiente de dispoinc es su entrada de la diagonal, la inferior derecha de la matriz de covarianzas estimada, 16.52 en el Ejemplo 6.6. Su raíz cuadrada, $\sqrt{16.52} = 4.064$ , es el error estándar reportado para dispoinc.

6.9 Resumen del capítulo¶

Ahora puedes hablar el lenguaje matricial en el que está escrito el resto del libro: distribuye un conjunto de datos como un vector respuesta $\mathbf{Y}$ y una matriz de diseño $\mathbf{X}$ , construye los productos $\mathbf{X}'\mathbf{X}$ y $\mathbf{X}'\mathbf{Y}$ sobre los que corre cada ajuste, di cuándo una matriz es invertible, resuelve las ecuaciones normales para $\mathbf{b} = (\mathbf{X}'\mathbf{X})^{-1}\mathbf{X}'\mathbf{Y}$ , construye la matriz sombrero y prueba que es una proyección, expresa las sumas de cuadrados como formas cuadráticas, y halla la media y la covarianza de un vector aleatorio. Para los datos reales de Dwaine, todo esto se calculó tanto en R como en numpy: $\mathbf{b} = (-68.86,\ 1.4546,\ 9.3655)'$ , $\mathrm{SSE} = 2180.9$ sobre 18 grados de libertad, $\mathrm{MSE} = 121.2$ , y errores estándar $(60.0,\ 0.212,\ 4.06)$ . De estas herramientas salieron los dos hechos que impulsan toda la regresión múltiple: $\mathbf{b}$ es insesgado, y su covarianza es $\sigma^2(\mathbf{X}'\mathbf{X})^{-1}$ .

Cada pieza encaja en un solo flujo, dibujado en Figure 11: a partir de los datos formas dos resúmenes, inviertes uno de ellos, y todo lo demás, los coeficientes, la matriz sombrero, los residuos y los errores estándar, sale en un orden fijo. Si recuerdas la forma de ese flujo, puedes reconstruir cualquier fórmula individual en él.

Un diagrama de flujo con siete cajas redondeadas conectadas por flechas naranjas. Arriba, una caja dice datos X y Y, sección 6.1. Alimenta dos cajas: X-transpuesta-X y X-transpuesta-Y, los productos cruzados, sección 6.2, y la inversa de X-transpuesta-X, sección 6.4. Ambas alimentan una caja central, b igual a la inversa de X-transpuesta-X por X-transpuesta-Y, los coeficientes, sección 6.4. La caja de coeficientes y la caja de la inversa alimentan luego tres cajas inferiores: la matriz sombrero H con valores ajustados y residuos, sección 6.5; SSE y MSE, sección 6.6; y la matriz de covarianzas de b con errores estándar, sección 6.7. Una leyenda dice: un flujo fijo, forma dos resúmenes, invierte uno, lee todo lo demás. — Figure 11:El capítulo como un solo flujo. Cada cantidad de un ajuste de regresión múltiple, los coeficientes, los valores ajustados, los residuos, la varianza del error y los errores estándar, proviene de los mismos dos resúmenes de productos cruzados y una sola inversa, calculada en un orden fijo.

Resultados clave de un vistazo.

Resultado	Enunciado o fórmula	Válido cuando
Transpuesta de un producto (Teorema 6.5)	$(\mathbf{A}\mathbf{B})' = \mathbf{B}'\mathbf{A}'$	cualesquiera $\mathbf{A}, \mathbf{B}$ conformables
Criterio de invertibilidad (Teorema 6.9)	$\mathbf{X}'\mathbf{X}$ invertible $\Leftrightarrow \mathbf{X}$ tiene rango columna completo	$\mathbf{X}'\mathbf{X}$ cuadrada
Solución de mínimos cuadrados (Teorema 6.12)	$\mathbf{b} = (\mathbf{X}'\mathbf{X})^{-1}\mathbf{X}'\mathbf{Y}$	$\mathbf{X}$ de rango columna completo
La matriz sombrero es una proyección (Teorema 6.16)	$\mathbf{H} = \mathbf{X}(\mathbf{X}'\mathbf{X})^{-1}\mathbf{X}'$ simétrica, idempotente, $\operatorname{tr}\mathbf{H} = p$	$\mathbf{X}$ de rango columna completo
Sumas de cuadrados como formas cuadráticas (Teorema 6.18)	$\mathrm{SSE} = \mathbf{Y}'(\mathbf{I}-\mathbf{H})\mathbf{Y}$ ; $\mathrm{SSTO} = \mathrm{SSR} + \mathrm{SSE}$ ; $\mathrm{SSE} \ge 0$	cualquier respuesta $\mathbf{Y}$
Reglas de transformación lineal (Teorema 6.20)	$E\{\mathbf{A}\mathbf{W}+\mathbf{c}\} = \mathbf{A}E\{\mathbf{W}\}+\mathbf{c}$ ; $\operatorname{Cov}\{\mathbf{A}\mathbf{W}\} = \mathbf{A}\operatorname{Cov}\{\mathbf{W}\}\mathbf{A}'$	$\mathbf{A}, \mathbf{c}$ fijos
Media y covarianza de $\mathbf{b}$ (Teorema 6.21)	$E\{\mathbf{b}\} = \boldsymbol{\beta}$ ; $\operatorname{Cov}\{\mathbf{b}\} = \sigma^2(\mathbf{X}'\mathbf{X})^{-1}$	$E\{\boldsymbol{\varepsilon}\} = \mathbf{0}$ , $\operatorname{Cov}\{\boldsymbol{\varepsilon}\} = \sigma^2\mathbf{I}$
Distribución muestral de $\mathbf{b}$ (Teorema 6.23)	$\mathbf{b} \sim N(\boldsymbol{\beta},\ \sigma^2(\mathbf{X}'\mathbf{X})^{-1})$	errores normales $\boldsymbol{\varepsilon} \sim N(\mathbf{0}, \sigma^2\mathbf{I})$

Términos clave. matriz, vector, matriz de diseño, transpuesta, multiplicación de matrices, conformable, matriz identidad, matriz simétrica, independencia lineal, rango, determinante, inversa, ecuaciones normales, matriz sombrero, matriz idempotente, matriz de proyección, traza, forma cuadrática, definida positiva, vector aleatorio, matriz de covarianzas, normal multivariante.

Ahora deberías poder.

Representar un conjunto de datos de regresión como un vector respuesta $\mathbf{Y}$ y una matriz de diseño $\mathbf{X}$ , y escribir el modelo como $\mathbf{Y} = \mathbf{X}\boldsymbol{\beta} + \boldsymbol{\varepsilon}$ .
Calcular transpuestas y productos de matrices, y construir $\mathbf{X}'\mathbf{X}$ y $\mathbf{X}'\mathbf{Y}$ a mano y con software.
Decidir cuándo una matriz cuadrada tiene inversa usando el rango y el determinante, y explicar qué le hace la colinealidad perfecta a $\mathbf{X}'\mathbf{X}$ .
Resolver las ecuaciones normales para $\mathbf{b} = (\mathbf{X}'\mathbf{X})^{-1}\mathbf{X}'\mathbf{Y}$ .
Verificar que la matriz sombrero $\mathbf{H}$ es simétrica e idempotente, y explicar por qué los valores ajustados y los residuos son proyecciones.
Expresar SSE, SSR y SSTO como formas cuadráticas, y conectar la definición positiva con sumas de cuadrados no negativas.
Derivar el vector de medias y la matriz de covarianzas de un vector aleatorio, y usarlos para mostrar que $\operatorname{Cov}\{\mathbf{b}\} = \sigma^2(\mathbf{X}'\mathbf{X})^{-1}$ .
Enunciar el modelo normal multivariante para $\mathbf{Y}$ y $\mathbf{b}$ y describir cómo impulsa la inferencia de los capítulos posteriores.

Dónde encaja esto. Este capítulo es la caja de herramientas que hace funcionar las etapas AJUSTAR y USAR del flujo de trabajo del curso (El flujo de trabajo del modelado) para más de un predictor: el estimador $\mathbf{b} = (\mathbf{X}'\mathbf{X})^{-1}\mathbf{X}'\mathbf{Y}$ ajusta toda regresión múltiple, y la matriz de covarianzas $\sigma^2(\mathbf{X}'\mathbf{X})^{-1}$ es como se prueba después cada coeficiente. El Capítulo 2 construyó estas ideas para un predictor con álgebra escalar (2.2 Mínimos cuadrados desde los primeros principios); las hemos generalizado ahora a cualquier número. Enseguida, el Capítulo 7 toma la misma $\mathbf{X}'\mathbf{X}$ de Dwaine que calculaste y prueba por qué estas fórmulas son las correctas: deriva los mínimos cuadrados de dos maneras (7.1 El modelo y los mínimos cuadrados en forma matricial), lee los grados de libertad en la matriz sombrero (7.3 La matriz sombrero), y prueba el teorema de Gauss-Markov (7.6 El teorema de Gauss-Markov). Es el capítulo más abstracto del curso, y está construido casi por completo a partir de los hechos matriciales que ensamblaste aquí, así que el trabajo de este capítulo es exactamente lo que hace legible el siguiente.

6.10 Preguntas frecuentes¶

P1. ¿Por qué pegamos una columna de unos a $\mathbf{X}$ ? Para que el intercepto pueda tratarse como un coeficiente más. Multiplicar la columna de unos por $\beta_0$ da $\beta_0$ en cada fila, que es exactamente lo que hace un intercepto. Sin la columna de unos, la fórmula $\mathbf{X}\boldsymbol{\beta}$ no tendría término constante, y estarías forzando la superficie de regresión a pasar por el origen.

P2. ¿Es $\mathbf{X}'\mathbf{X}$ lo mismo que elevar al cuadrado $\mathbf{X}$ ? No. $\mathbf{X}$ por lo general no es cuadrada, así que $\mathbf{X}^2 = \mathbf{X}\mathbf{X}$ ni siquiera está definida. $\mathbf{X}'\mathbf{X}$ multiplica la transpuesta ( $p \times n$ ) por $\mathbf{X}$ ( $n \times p$ ) para hacer una matriz $p \times p$ de sumas de productos. Es el análogo matricial más cercano de “suma de cuadrados”, razón por la cual se sitúa en el centro de los mínimos cuadrados.

P3. ¿Qué sale mal en realidad cuando $\mathbf{X}'\mathbf{X}$ es singular? Dos columnas de predictores cargan la misma información, así que los datos no pueden decidir cómo repartir un efecto entre ellas. Infinitos vectores de coeficientes ajustan igualmente bien, la inversa no existe, y el software o bien da error o bien descarta un predictor en silencio. Este es el extremo de la multicolinealidad que encuentras en 12.1 La multicolinealidad y el factor de inflación de la varianza.

P4. ¿Tengo que invertir $\mathbf{X}'\mathbf{X}$ a mano? No. Más allá de $2 \times 2$ , deja que solve o np.linalg.inv lo hagan, y en el trabajo real prefiere las rutinas dentro de lm y statsmodels, que resuelven las ecuaciones normales sin formar la inversa. Calcular la inversa aquí solo muestra que el software hace exactamente el álgebra de matrices de este capítulo, de modo que nunca es una caja negra.

P5. ¿Por qué se llama proyección a la matriz sombrero? Porque toma el vector respuesta $\mathbf{Y}$ y lo deja caer perpendicularmente sobre el subespacio plano de todos los valores ajustables (el espacio columna de $\mathbf{X}$ ), aterrizando en el punto más cercano, $\widehat{\mathbf{Y}}$ . Proyectar una segunda vez no cambia nada, que es la propiedad algebraica $\mathbf{H}\mathbf{H} = \mathbf{H}$ . El residuo es la parte de $\mathbf{Y}$ que sobresale, en ángulo recto respecto del subespacio.

P6. ¿Dónde entra la distribución normal? Solo en 6.8 La normal multivariante, en breve, y solo para la inferencia. Construir $\mathbf{X}$ , calcular $\mathbf{b}$ , la matriz sombrero, SSE, y la matriz de covarianzas no usan ningún supuesto distribucional más allá de media cero, varianza constante y errores no correlacionados. La normalidad se añade encima para que $\mathbf{b}$ sea normal multivariante y podamos construir procedimientos $t$ y $F$ exactos, que es la tarea del Capítulo 7.

P7. ¿Tengo que memorizar todas estas identidades, y deben preocuparme las entradas negativas en la matriz de covarianzas de $\mathbf{b}$ ? No, en ambos casos. Retén tres cosas: el modelo $\mathbf{Y} = \mathbf{X}\boldsymbol{\beta} + \boldsymbol{\varepsilon}$ , el estimador $\mathbf{b} = (\mathbf{X}'\mathbf{X})^{-1}\mathbf{X}'\mathbf{Y}$ , y la covarianza $\sigma^2(\mathbf{X}'\mathbf{X})^{-1}$ ; todo lo demás construye una de esas o lee un número de ella, así que las identidades se vuelven cosas que rederivas en vez de recordar. En cuanto a las entradas negativas, solo las covarianzas fuera de la diagonal pueden ser negativas (las varianzas de la diagonal nunca lo son), y una negativa entre dos pendientes solo significa que a lo largo de muestras repetidas sobreestimar una tiende a acompañar subestimar la otra. Es información, no un error.

6.11 Problemas de práctica¶

(A) Da las dimensiones de $\mathbf{X}$ , $\mathbf{X}'$ , $\mathbf{X}'\mathbf{X}$ , $\mathbf{X}'\mathbf{Y}$ , $\mathbf{b}$ y $\mathbf{H}$ para el modelo de Dwaine, y di en una frase qué representa cada uno.
(A) Explica por qué la primera columna de la matriz de diseño es de puros unos, y qué cambiaría en el modelo si se eliminara.
(A) Enuncia la regla de cuándo dos matrices se pueden multiplicar, y úsala para explicar por qué $\mathbf{X}\mathbf{X}$ no está definida para la matriz de diseño $21 \times 3$ de Dwaine pero $\mathbf{X}'\mathbf{X}$ sí lo está.
(A) En palabras, ¿a qué es igual la entrada en la fila 1, columna 1 de $\mathbf{X}'\mathbf{X}$ , y por qué? ¿Cuáles son las otras entradas de la primera fila?
(A) Define una matriz idempotente y una matriz simétrica, y di cuáles de estas propiedades tiene la matriz sombrero $\mathbf{H}$ .
(A) Un colega escribe $\mathbf{X}\mathbf{b} = \mathbf{Y}$ y “cancela $\mathbf{X}$ ” de $\mathbf{X}'\mathbf{X}\mathbf{b} = \mathbf{X}'\mathbf{Y}$ . Explica las dos cosas mal en esto.
(A) Explica qué significa que las columnas de $\mathbf{X}$ sean linealmente dependientes, y qué le hace eso a $\mathbf{X}'\mathbf{X}$ y al ajuste de mínimos cuadrados.
(A) La matriz de covarianzas de $\mathbf{b}$ tiene una entrada fuera de diagonal negativa entre las dos pendientes. Interpreta su signo en una oración.
(B) Prueba que $\mathbf{X}'\mathbf{X}$ es simétrica para cualquier matriz $\mathbf{X}$ , citando la regla de la transpuesta de un producto (Teorema 6.5).
(B) Prueba la regla de la transpuesta de un producto (Teorema 6.5), $(\mathbf{A}\mathbf{B})' = \mathbf{B}'\mathbf{A}'$ , comparando las entradas $(i,j)$ de ambos lados.
(B) Partiendo de las ecuaciones normales matriciales $\mathbf{X}'\mathbf{X}\,\mathbf{b} = \mathbf{X}'\mathbf{Y}$ , deriva $\mathbf{b} = (\mathbf{X}'\mathbf{X})^{-1}\mathbf{X}'\mathbf{Y}$ (Teorema 6.12), enunciando la condición sobre $\mathbf{X}$ que requiere el paso.
(B) Prueba que la matriz sombrero $\mathbf{H} = \mathbf{X}(\mathbf{X}'\mathbf{X})^{-1}\mathbf{X}'$ es simétrica e idempotente (Teorema 6.16). Luego muestra que $\mathbf{I} - \mathbf{H}$ es idempotente.
(B) Muestra que $\mathbf{H}\mathbf{X} = \mathbf{X}$ , y explica geométricamente por qué proyectar cada columna de $\mathbf{X}$ sobre el espacio columna de $\mathbf{X}$ la deja sin cambio.
(B) Usando $\mathbf{e} = (\mathbf{I} - \mathbf{H})\mathbf{Y}$ y $\widehat{\mathbf{Y}} = \mathbf{H}\mathbf{Y}$ , prueba que $\mathbf{X}'\mathbf{e} = \mathbf{0}$ y que $\widehat{\mathbf{Y}}'\mathbf{e} = 0$ (los valores ajustados son ortogonales a los residuos).
(B) Prueba la regla de transformación de la covarianza $\operatorname{Cov}\{\mathbf{A}\mathbf{W}\} = \mathbf{A}\operatorname{Cov}\{\mathbf{W}\}\mathbf{A}'$ (parte del Teorema 6.20) para una matriz fija $\mathbf{A}$ , partiendo de la definición $\operatorname{Cov}\{\mathbf{W}\} = E\{(\mathbf{W} - \boldsymbol{\mu})(\mathbf{W} - \boldsymbol{\mu})'\}$ .
(B) Usa las reglas del problema 15 y $E\{\mathbf{Y}\} = \mathbf{X}\boldsymbol{\beta}$ , $\operatorname{Cov}\{\mathbf{Y}\} = \sigma^2\mathbf{I}$ para derivar $E\{\mathbf{b}\} = \boldsymbol{\beta}$ y $\operatorname{Cov}\{\mathbf{b}\} = \sigma^2(\mathbf{X}'\mathbf{X})^{-1}$ (Teorema 6.21).
(B) Muestra que $\mathrm{SSE} = \mathbf{Y}'(\mathbf{I} - \mathbf{H})\mathbf{Y}$ es igual a $\mathbf{e}'\mathbf{e}$ (Teorema 6.18), usando la simetría y la idempotencia de $\mathbf{I} - \mathbf{H}$ , y explica por qué esto prueba $\mathrm{SSE} \ge 0$ .
(B) Para la matriz $2 \times 2$ $\mathbf{A} = \left(\begin{smallmatrix} a & b \\ b & d \end{smallmatrix}\right)$ , escribe la forma cuadrática $\mathbf{u}'\mathbf{A}\mathbf{u}$ en términos de $u_1, u_2$ , y da una condición sobre $a, b, d$ que la vuelva definida positiva.
(B) Los valores ajustados satisfacen $\widehat{\mathbf{Y}} = \mathbf{H}\mathbf{Y}$ . Prueba que $\operatorname{Cov}\{\widehat{\mathbf{Y}}\} = \sigma^2\mathbf{H}$ y $\operatorname{Cov}\{\mathbf{e}\} = \sigma^2(\mathbf{I} - \mathbf{H})$ , usando la regla de la covarianza y las propiedades de $\mathbf{H}$ .
(C) Lee dwaine.csv, construye $\mathbf{X}$ y $\mathbf{Y}$ , y reproduce $\mathbf{X}'\mathbf{X}$ y $\mathbf{X}'\mathbf{Y}$ en R o Python. Confirma que la entrada superior izquierda de $\mathbf{X}'\mathbf{X}$ es $n$ y que la primera entrada de $\mathbf{X}'\mathbf{Y}$ es $\sum Y$ .
(C) Calcula $(\mathbf{X}'\mathbf{X})^{-1}$ y $\mathbf{b} = (\mathbf{X}'\mathbf{X})^{-1}\mathbf{X}'\mathbf{Y}$ en software, y verifica que $\mathbf{b}$ concuerda con lm/statsmodels y que $(\mathbf{X}'\mathbf{X})^{-1}\mathbf{X}'\mathbf{X}$ devuelve la identidad.
(C) Construye la matriz sombrero $\mathbf{H}$ , verifica numéricamente que $\mathbf{H}\mathbf{H} = \mathbf{H}$ y que $\operatorname{tr}(\mathbf{H}) = 3$ , y reporta las tres entradas de la diagonal $h_{ii}$ más grandes (los valores de apalancamiento).
(C) Calcula el vector de residuos $\mathbf{e} = (\mathbf{I} - \mathbf{H})\mathbf{Y}$ y verifica numéricamente que $\mathbf{X}'\mathbf{e} = \mathbf{0}$ (las tres entradas cerca de cero) y que $\sum e_i = 0$ .
(C) Calcula $\mathrm{SSE} = \mathbf{e}'\mathbf{e}$ , luego $\mathrm{MSE}$ y $s = \sqrt{\mathrm{MSE}}$ , y confirma que concuerdan con el error estándar residual reportado por summary(fit) / fit.summary().
(C) Forma la matriz de covarianzas estimada $\mathrm{MSE}\,(\mathbf{X}'\mathbf{X})^{-1}$ y reporta los tres errores estándar de su diagonal. Confirma que concuerdan con el resumen del software, y reporta la covarianza estimada entre las dos estimaciones de pendiente.
(C) Predice las ventas para una ciudad nueva con $\text{targtpop} = 65.4$ y $\text{dispoinc} = 17.6$ formando el vector fila $\mathbf{x}_h = (1, 65.4, 17.6)$ y calculando $\mathbf{x}_h \mathbf{b}$ . Confirma el resultado contra predict.
(C) Añade una columna redundante a $\mathbf{X}$ igual a $\text{targtpop} + \text{dispoinc}$ , e intenta calcular $(\mathbf{X}'\mathbf{X})^{-1}$ . Reporta lo que hace R o Python (un error, una advertencia, o una inversa desmesuradamente inestable), y conéctalo con el rango deficiente.
(C) Centra los predictores: reemplaza targtpop y dispoinc por sus desviaciones respecto de sus medias, reajusta, y confirma que las dos pendientes no cambian mientras que el intercepto se vuelve $\bar{Y}$ . Explica, usando las ecuaciones normales, por qué el centrado deja intactas las pendientes.
(C) Demuestra la ley muestral $\mathbf{b} \sim N(\boldsymbol{\beta}, \sigma^2(\mathbf{X}'\mathbf{X})^{-1})$ (Teorema 6.23, 6.8 La normal multivariante, en breve) por simulación. Trata la $\mathbf{b} = (-68.86,\ 1.4546,\ 9.3655)'$ ajustada como la $\boldsymbol{\beta}$ verdadera y $\mathrm{MSE} = 121.16$ como la $\sigma^2$ verdadera. Manteniendo fija la $\mathbf{X}$ real de Dwaine, usa set.seed(4210) (R) o default_rng(4210) (Python) para generar 5000 vectores respuesta $\mathbf{Y} = \mathbf{X}\boldsymbol{\beta} + \boldsymbol{\varepsilon}$ con $\boldsymbol{\varepsilon} \sim N(\mathbf{0}, \sigma^2\mathbf{I})$ , reajusta cada uno con $(\mathbf{X}'\mathbf{X})^{-1}\mathbf{X}'\mathbf{Y}$ , y recolecta las estimaciones. (a) Confirma que las medias por columna están cerca de $\boldsymbol{\beta}$ y que la covarianza muestral está cerca de $\sigma^2(\mathbf{X}'\mathbf{X})^{-1}$ . (b) En dos oraciones, explica por qué esta ley permite que el Capítulo 7 adjunte una distribución $t$ a cada coeficiente y convierta un error estándar en un intervalo de confianza.

6.12 Práctica de examen¶

Estas cinco preguntas coinciden con el estilo de los exámenes del curso: cada una te pide explicar, evaluar o interpretar en oraciones completas, no producir un número pelado. Escribe en oraciones completas y di qué números usaste. La salida de software mostrada se generó a partir del ajuste real de dwaine.csv. Cada respuesta modelo muestra la profundidad que gana la calificación completa, seguida de una línea sobre lo que una respuesta débil omite.

EP 6.1. En el ajuste de Dwaine ambas pendientes estimadas son positivas: $b_{\text{targtpop}} = 1.4546$ y $b_{\text{dispoinc}} = 9.3655$ . Un estudiante argumenta: “como ambos predictores empujan las ventas hacia arriba, las dos pendientes estimadas deben estar positivamente correlacionadas a lo largo de muestras repetidas”. La matriz de covarianzas estimada de $\mathbf{b}$ está impresa abajo. Evalúa la afirmación del estudiante, y explica qué te dice en realidad el número relevante.

          intercept targtpop  dispoinc
intercept 3602.0347   8.7459 -241.4230
targtpop     8.7459    0.0449   -0.6724
dispoinc  -241.4230   -0.6724   16.5158

Respuesta modelo

La afirmación confunde el signo de los coeficientes con el signo de la covarianza entre sus estimadores, y los dos no tienen relación. El número que resuelve la pregunta es la entrada fuera de diagonal entre las dos pendientes, $s\{b_{\text{targtpop}}, b_{\text{dispoinc}}\} = -0.6724$ , que es negativa. Así que a lo largo de muestras repetidas una sobreestimación de una pendiente tiende a llegar con una subestimación de la otra, exactamente lo opuesto de lo que predice el estudiante. Ese signo negativo viene de la geometría de la matriz de diseño: las dos columnas de predictores están ellas mismas positivamente correlacionadas a lo largo de las 21 ciudades, así que el ajuste no puede subir ambas pendientes a la vez sin contar dos veces la señal compartida, y compensa una contra la otra. Los signos de las pendientes ajustadas, que describen cómo responden las ventas a cada predictor, no llevan información sobre cómo covarían las dos estimaciones. Una respuesta débil solo declara la afirmación verdadera o falsa, o repite que ambas pendientes son positivas, sin leer la entrada -0.6724 ni explicar que es negativa por razones ligadas a los predictores, no a los coeficientes.

EP 6.2. Supón que targtpop se reingresa en personas en bruto en vez de miles, de modo que cada valor se multiplica por 1000, y el modelo se reajusta. Usando la salida de abajo, di con precisión qué cantidades cambian y cuáles quedan idénticas, y explica por qué.

          b            se           t
intercept -68.8570732  60.0169532  -1.1473
targtpop    0.0014546   0.0002118   6.8682
dispoinc    9.3655004   4.0639581   2.3045
SSE = 2180.9274      fitted[1:3] = 187.184, 154.229, 234.396

Respuesta modelo

Multiplicar una columna de predictor por 1000 divide la pendiente de ese predictor y el error estándar de la pendiente por el mismo factor de 1000: la pendiente de targtpop cae de 1.4546 a 0.0014546 y su error estándar de 0.2118 a 0.0002118, así que su cociente, el estadístico $t$ de 6.868, no cambia. Todo lo demás es idéntico al ajuste original: el intercepto -68.857, la pendiente de dispoinc 9.3655, todos los valores ajustados (187.184, 154.229, 234.396, ...), los residuos, $\mathrm{SSE} = 2180.93$ , $\mathrm{MSE}$ , y $R^2$ . La razón es que escalar una columna de $\mathbf{X}$ por una constante distinta de cero no cambia el espacio columna que $\mathbf{X}$ genera, y el ajuste depende de $\mathbf{X}$ solo a través de ese espacio: la matriz sombrero $\mathbf{H}$ , la proyección $\widehat{\mathbf{Y}} = \mathbf{H}\mathbf{Y}$ , y toda suma de cuadrados quedan intactas. Un cambio de unidades por lo tanto reescala el coeficiente afectado y su error estándar al unísono y deja en paz cada valor ajustado y cada inferencia. Una respuesta débil dice solo que “la pendiente cambia” sin notar que su error estándar cambia por el factor idéntico, que el estadístico $t$ , los valores ajustados y SSE son invariantes, y por qué.

EP 6.3. La diagonal de la matriz sombrero contiene los valores de apalancamiento $h_{ii}$ . La salida de abajo reporta su suma, el corte de regla empírica $2p/n$ , y los cinco valores más grandes con los dos valores de predictor de cada ciudad. Interpreta estos números en contexto: ¿qué dicen sobre qué ciudades importan más, por qué suman 3, y marcarías alguna ciudad como punto de apalancamiento alto?

sum of h_ii = 3.0        2p/n cutoff = 0.2857
 city 20:  targtpop 82.7  dispoinc 19.1   h = 0.2788
 city 13:  targtpop 88.4  dispoinc 17.4   h = 0.2390
 city 15:  targtpop 52.5  dispoinc 17.8   h = 0.2095
 city  3:  targtpop 91.3  dispoinc 18.2   h = 0.1737
 city  5:  targtpop 46.9  dispoinc 17.3   h = 0.1620

Respuesta modelo

Cada valor de apalancamiento $h_{ii}$ mide qué tan lejos están los valores de predictor de la ciudad $i$ del centro de la nube de predictores, y es igual al peso que las ventas propias de la ciudad cargan en producir su valor ajustado. La ciudad 20 tiene el mayor valor de apalancamiento ( $h = 0.2788$ ) porque combina una gran población joven (82.7 mil) con el ingreso disponible más alto de los datos (19.1), así que se sitúa en el borde exterior de la nube de predictores y su respuesta tiene la mayor voz en el ajuste local. Los apalancamientos suman 3 porque $\operatorname{tr}(\mathbf{H}) = p = 3$ , el número de parámetros estimados, así que en promedio una ciudad carga $3/21 \approx 0.143$ de su propio valor ajustado; este es el mismo hecho de la traza que produce el divisor $n - p = 18$ para $\mathrm{MSE}$ . Al comparar contra el corte $2p/n = 0.2857$ , ninguna ciudad lo excede, pues aun la ciudad 20 con 0.2788 queda justo por debajo, así que ninguna se marca como punto de apalancamiento inusualmente alto, aunque la ciudad 20 es la que hay que vigilar. Una respuesta débil lee el número más grande sin decir qué mide un valor de apalancamiento, omite que la suma es igual a $p$ porque $\operatorname{tr}(\mathbf{H}) = p$ , o llama a una ciudad de apalancamiento alto por su tamaño en bruto sin compararla con el corte 0.2857.

EP 6.4. Un estudiante aumenta la matriz de diseño de Dwaine con una cuarta columna igual a targtpop + dispoinc, reajusta, y reporta: “Python de todos modos devolvió una inversa y algunos coeficientes, así que el modelo aumentado está bien”. Los diagnósticos de abajo provienen de la matriz aumentada. Explica por qué el modelo aumentado no está bien, qué salió mal matemáticamente, y por qué R y Python se comportan distinto.

rank of augmented X = 3   (it has 4 columns)
det(X'X)        = -3.7e-05      # against ~1.07e06 for the genuine 3x3 X'X
R  solve(X'X):  Error: system is computationally singular
Python inv(X'X): returns a matrix of enormous (~1e13) entries, no error

Respuesta modelo

La nueva columna es una combinación lineal exacta de dos columnas existentes, a saber $\text{targtpop} + \text{dispoinc}$ , así que las cuatro columnas son linealmente dependientes y la matriz de diseño tiene rango 3, no 4. Por el Teorema 6.9, $\mathbf{X}'\mathbf{X}$ es entonces singular y no tiene una inversa genuina: su determinante verdadero es cero, y el $-3.7\times 10^{-5}$ impreso es polvo de punto flotante alrededor de cero, astronómicamente pequeño al lado del determinante de aproximadamente $1.07\times 10^{6}$ de la $\mathbf{X}'\mathbf{X}$ real $3\times 3$ . No hay una única solución de mínimos cuadrados, porque los datos no pueden decidir cómo repartir un efecto entre un predictor y una copia construida de él; infinitos vectores de coeficientes ajustan igualmente bien. R calcula el número de condición recíproco, ve que es esencialmente cero, y se rehúsa con “system is computationally singular”, que es la respuesta honesta. El np.linalg.inv de Python invierte sin comprobar el condicionamiento, así que devuelve una matriz de entradas enormes y sin sentido que son ruido numérico en vez de un ajuste. El estudiante ha confundido una falla numérica silenciosa con un modelo válido: el hecho de que el código “corrió” no es evidencia de que la respuesta signifique algo. Una respuesta débil nota solo que las columnas son colineales sin conectar el rango deficiente con el determinante cero y la ausencia de una única solución, o confía en los números de Python porque no apareció ningún error.

EP 6.5. Un gerente regional quiere las ventas medias predichas por el modelo para una ciudad nueva con targtpop $= 65.4$ y dispoinc $= 17.6$ , junto con un error estándar, así que forma $\mathbf{x}_h = (1,\ 65.4,\ 17.6)$ y nota que la media predicha es $\widehat{Y}_h = \mathbf{x}_h'\mathbf{b}$ . La matriz de covarianzas estimada $s^2\{\mathbf{b}\} = \mathrm{MSE}\,(\mathbf{X}'\mathbf{X})^{-1}$ se reimprime abajo, y el ajuste da $\widehat{Y}_h = 191.10$ . Un estudiante calcula el error estándar de $\widehat{Y}_h$ como $\sqrt{65.4^2 \cdot s^2\{b_{\text{targtpop}}\} + 17.6^2 \cdot s^2\{b_{\text{dispoinc}}\} + s^2\{b_0\}} = 94.4$ , usando solo las varianzas de la diagonal. Explica por qué esto está mal, da el error estándar correcto, y di hacia qué lado yerra el atajo de solo diagonal.

          intercept targtpop  dispoinc
intercept 3602.0347   8.7459 -241.4230
targtpop     8.7459    0.0449   -0.6724
dispoinc  -241.4230   -0.6724   16.5158

x_h'(X'X)^{-1} x_h = 0.063181       MSE = 121.16

Respuesta modelo

La media predicha $\widehat{Y}_h = \mathbf{x}_h'\mathbf{b}$ es una combinación lineal de los tres coeficientes, así que por la regla de la covarianza del Teorema 6.20 su varianza es la forma cuadrática completa $\operatorname{Var}\{\widehat{Y}_h\} = \mathbf{x}_h'\,s^2\{\mathbf{b}\}\,\mathbf{x}_h$ , que incluye toda covarianza fuera de la diagonal, no solo las tres varianzas de la diagonal. El atajo del estudiante descarta los términos cruzados $2 x_{h1} x_{h2}\,s\{b_j, b_k\}$ , y aquí esos términos son grandes y negativos, sobre todo la covarianza -241.42 entre el intercepto y la pendiente de dispoinc, así que ignorarlos exagera desmesuradamente la varianza. Calculándola correctamente, $\operatorname{Var}\{\widehat{Y}_h\} = \mathrm{MSE}\cdot \mathbf{x}_h'(\mathbf{X}'\mathbf{X})^{-1}\mathbf{x}_h = 121.16 \times 0.063181 = 7.655$ , así que el error estándar es $\sqrt{7.655} = 2.77$ mil dólares, no 94.4. El atajo de solo diagonal yerra por un margen enorme hacia el lado alto, porque las covarianzas negativas que los coeficientes genuinamente tienen cancelarían gran parte de la contribución de la diagonal. La lección es que el error estándar de cualquier combinación lineal de $\mathbf{b}$ necesita toda la matriz de covarianzas, ya que los coeficientes están correlacionados. Una respuesta débil calcula $\mathbf{x}_h'\,s^2\{\mathbf{b}\}\,\mathbf{x}_h$ como fórmula sin explicar que los términos fuera de diagonal que faltan son los que vuelven incorrecto el atajo, o no nota que el atajo exagera en vez de subestimar el error estándar.

6.13 Juego del capítulo¶

Chapter summary (in English)

This chapter teaches the matrix algebra that multiple regression needs, assuming no prior linear algebra. It uses the Dwaine Studios data: the sales of 21 portrait studios against two predictors, the population under 16 and per-capita disposable income. The data are laid out as a response vector $\mathbf{Y}$ and a design matrix $\mathbf{X}$ of size $21 \times 3$ , whose first column is all ones for the intercept, and the model is written compactly as $\mathbf{Y} = \mathbf{X}\boldsymbol{\beta} + \boldsymbol{\varepsilon}$ .

We define the transpose and matrix multiplication, whose rule is “row by column,” and with them build the two central pieces: $\mathbf{X}'\mathbf{X}$ , which collects every sum of squares and cross-product, and $\mathbf{X}'\mathbf{Y}$ . A matrix has an inverse only if it has full rank, equivalently if its determinant is nonzero. Solving the normal equations $\mathbf{X}'\mathbf{X}\,\mathbf{b} = \mathbf{X}'\mathbf{Y}$ gives the estimator $\mathbf{b} = (\mathbf{X}'\mathbf{X})^{-1}\mathbf{X}'\mathbf{Y}$ ; for Dwaine, $\mathbf{b} = (-68.86,\ 1.4546,\ 9.3655)$ , identical to lm and statsmodels.

The hat matrix $\mathbf{H} = \mathbf{X}(\mathbf{X}'\mathbf{X})^{-1}\mathbf{X}'$ is symmetric and idempotent: it projects $\mathbf{Y}$ onto the fitted values. The sums of squares SSE, SSR, and SSTO are quadratic forms, and for Dwaine $\mathrm{SSE} = 2180.9$ with $\mathrm{MSE} = 121.2$ on 18 degrees of freedom. For a random vector, the expectation and covariance matrix transform by $E\{\mathbf{A}\mathbf{W}\} = \mathbf{A}E\{\mathbf{W}\}$ and $\operatorname{Cov}\{\mathbf{A}\mathbf{W}\} = \mathbf{A}\operatorname{Cov}\{\mathbf{W}\}\mathbf{A}'$ , from which $\mathbf{b}$ is unbiased and $\operatorname{Cov}\{\mathbf{b}\} = \sigma^2(\mathbf{X}'\mathbf{X})^{-1}$ . Under normal errors, $\mathbf{b}$ follows a multivariate normal distribution, the foundation of all the inference in Chapter 7.