Regresión por mínimos cuadrados Cálculo de álgebra lineal paso a paso

Question

Como precuela a una pregunta sobre modelos lineales mixtos en R, y para compartir como referencia para los aficionados a la estadística principiantes/intermedios, he decidido publicar como "Q&A-style" independiente los pasos implicados en el cálculo "manual" de los coeficientes y valores predichos de una regresión lineal simple.

El ejemplo es con el conjunto de datos incorporado en R, mtcars y se configuraría como millas por galón consumidas por un vehículo actuando como variable independiente, regresionada sobre el peso del coche (variable continua), y el número de cilindros como factor con tres niveles (4, 6 u 8) sin interacciones.

EDIT: Si le interesa esta cuestión, encontrará sin duda una respuesta detallada y satisfactoria en esto post de Matthew Drury fuera de CV .

Answer 1

Nota : He publicado una versión ampliada de esta respuesta en mi sitio web .

¿Sería tan amable de publicar una respuesta similar con el motor R real expuesto?

Claro que sí. Vamos a la madriguera del conejo.

La primera capa es lm la interfaz expuesta al programador de R. Puedes ver el código fuente de esto simplemente escribiendo lm en la consola R. La mayor parte (como la mayor parte del código de producción) está ocupada comprobando entradas, estableciendo atributos de objetos y lanzando errores; pero esta línea llama la atención

lm.fit(x, y, offset = offset, singular.ok = singular.ok, 
                ...)

lm.fit es otra función de R, puedes llamarla tú mismo. Mientras que lm trabaja cómodamente con fórmulas y marcos de datos, lm.fit quiere matrices, así que eso es un nivel de abstracción eliminado. Comprobación de la fuente de lm.fit más trabajo, y la siguiente línea realmente interesante

z <- .Call(C_Cdqrls, x, y, tol, FALSE)

Ahora estamos llegando a alguna parte. .Call es la forma que tiene R de llamar al código C. Hay una función C, C_Cdqrls en el código fuente de R en alguna parte, y tenemos que encontrarla. Aquí está .

Observando la función C, de nuevo, encontramos principalmente comprobación de límites, limpieza de errores y trabajo ocupado. Pero esta línea es diferente

F77_CALL(dqrls)(REAL(qr), &n, &p, REAL(y), &ny, &rtol,
        REAL(coefficients), REAL(residuals), REAL(effects),
        &rank, INTEGER(pivot), REAL(qraux), work);

Así que ahora estamos en nuestro tercer lenguaje, R ha llamado a C que está llamando a fortran. Aquí está el código fortran .

El primer comentario lo dice todo

c     dqrfit is a subroutine to compute least squares solutions
c     to the system
c
c     (1)               x * b = y

(curiosamente, parece que el nombre de esta rutina se cambió en algún momento, pero alguien olvidó actualizar el comentario). Así que por fin estamos en el punto donde podemos hacer un poco de álgebra lineal, y realmente resolver el sistema de ecuaciones. Este es el tipo de cosas para las que Fortran es realmente bueno, lo que explica por qué hemos pasado por tantas capas para llegar hasta aquí.

El comentario también explica lo que va a hacer el código

c     on return
c
c        x      contains the output array from dqrdc2.
c               namely the qr decomposition of x stored in
c               compact form.

Así que fortran va a resolver el sistema encontrando el $QR$ descomposición.

Lo primero que ocurre, y con mucho lo más importante, es

call dqrdc2(x,n,n,p,tol,k,qraux,jpvt,work)

Esto llama a la función fortran dqrdc2 en nuestra matriz de entrada x . ¿Qué es esto?

 c     dqrfit uses the linpack routines dqrdc and dqrsl.

Por fin hemos llegado a linpack . Linpack es una biblioteca de álgebra lineal fortran que existe desde los años 70. La mayoría del álgebra lineal seria eventualmente encuentra su camino a linpack. En nuestro caso, estamos utilizando la función dqrdc2

c     dqrdc2 uses householder transformations to compute the qr
c     factorization of an n by p matrix x.

Aquí es donde se hace el trabajo de verdad. Me llevaría un día entero entender lo que hace este código, es de muy bajo nivel. Pero genéricamente, tenemos una matriz $X$ y queremos factorizarlo en un producto $X = QR$ donde $Q$ es una matriz ortogonal y $R$ es una matriz triangular superior. Esto es algo inteligente, porque una vez que se tiene $Q$ y $R$ puede resolver las ecuaciones lineales de regresión

$$ X^t X \beta = X^t Y $$

muy fácilmente. De hecho

$$ X^t X = R^t Q^t Q R = R^t R $$

por lo que todo el sistema se convierte en

$$ R^t R \beta = R^t Q^t y $$

pero $R$ es triangular superior y tiene el mismo rango que $X^t X$ por lo que si nuestro problema está bien planteado, es de rango completo, y podemos resolver el sistema reducido.

$$ R \beta = Q^t y $$

Pero esto es lo increíble. $R$ es triangular superior, por lo que la última ecuación lineal aquí es simplemente constant * beta_n = constant por lo que resolviendo para $\beta_n$ es trivial. A continuación, puede subir las filas, una por una, y sustituir en el $\beta$ s que ya conoces, cada vez obteniendo una simple ecuación lineal de una variable para resolver. Así, una vez que tengas $Q$ y $R$ todo se reduce a lo que se llama sustitución regresiva que es fácil. Puede leer más información al respecto aquí en el que se elabora un pequeño ejemplo explícito.

Answer 2

Los cálculos reales paso a paso en R están bellamente descritos en la respuesta de Matthew Drury en este mismo hilo. En esta respuesta quiero caminar a través del proceso de demostrar a uno mismo que los resultados en R con un ejemplo sencillo se puede llegar siguiendo el álgebra lineal de las proyecciones sobre el espacio de columnas, y perpendicular (producto punto) errores concepto, ilustrado en diferentes mensajes, y muy bien explicado por el Dr. Strang en Álgebra lineal y sus aplicaciones y de fácil acceso aquí .

Para estimar los coeficientes $\small \beta$ en la regresión,

$$\small mpg = intercept\,(cyl=4) + \beta_1\,*\,weight + D1\,* intercept\,(cyl=6) + D2\, * intercept\,(cyl=8)\,\,\,\,[*]$$

con $\small D1$ y $\small D2$ que representan variables ficticias con valores [0,1], primero tendríamos que incluir en el matriz de diseño ( $\small X$ ) la codificación ficticia del número de cilindros, como sigue:

attach(mtcars)    
x1 <- wt

    x2 <- cyl; x2[x2==4] <- 1; x2[!x2==1] <-0

    x3 <- cyl; x3[x3==6] <- 1; x3[!x3==1] <-0

    x4 <- cyl; x4[x4==8] <- 1; x4[!x4==1] <-0

    X <- cbind(x1, x2, x3, x4)
    colnames(X) <-c('wt','4cyl', '6cyl', '8cyl')

head(X)
        wt 4cyl 6cyl 8cyl
[1,] 2.620    0    1    0
[2,] 2.875    0    1    0
[3,] 2.320    1    0    0
[4,] 3.215    0    1    0
[5,] 3.440    0    0    1
[6,] 3.460    0    1    0

Si la matriz de diseño tuviera que ser estrictamente paralela a la ecuación $\small [*]$ (arriba), donde el primer intercepto corresponde a coches de cuatro cilindros, como en el lm sin un `-1', requeriría una primera columna de sólo unos, pero obtendremos los mismos resultados sin esta columna de intercepción.

Continuando entonces, para calcular los coeficientes ( $\small\beta$ ) necesitamos el matriz de proyección - proyectamos el vector de valores de la variable independiente sobre el espacio de columnas de los vectores que constituyen el matriz de diseño . El álgebra lineal es $\small ProjMatrix = \small (X^{T}X)^{-1}X^{T}$ que se multiplica por el vector de la variable independiente: $\small [ProjMatrix] \,[y]\, =\, [RegrCoef's]$ o $\small (X^{T}X)^{-1}X^{T}\,y = \beta$ :

X_tr_X_inv <- solve(t(X) %*% X)    
Proj_M <- X_tr_X_inv %*% t(X)
Proj_M %*% mpg

          [,1]
wt   -3.205613
4cyl 33.990794
6cyl 29.735212
8cyl 27.919934

Idéntica a: coef(lm(mpg ~ wt + as.factor(cyl)-1)) .

Por último, para calcular los valores previstos, necesitaremos el matriz de sombreros que se define como, $\small Hat Matrix = \small X(X^{T}X)^{-1}X^{T}$ . Esto se calcula fácilmente como:

HAT <- X %*% X_tr_X_inv %*% t(X)

Y la estimación ( $\hat{y}$ ) como $\small X(X^{T}X)^{-1}X^{T}\,y$ en este caso: y_hat <- HAT %*% mpg que da valores idénticos a:

cyl <- as.factor(cyl); OLS <- lm(mpg ~ wt + cyl); predict(OLS) :

y_hat <- as.numeric(y_hat)
predicted <- as.numeric(predict(OLS))
all.equal(y_hat,predicted)
[1] TRUE