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

Como precuela de una pregunta sobre modelos lineales mixtos en R, y para compartir como referencia para principiantes / intermedios aficionados a la estadística, decidí publicar como un "estilo de preguntas y respuestas" independiente los pasos involucrados en el cálculo "manual" del coeficientes y valores predichos de una regresión lineal simple.

El ejemplo es con el conjunto de datos R incorporado mtcars, y se configuraría como millas por galón consumido por un vehículo que actúa como la variable independiente, retrocediendo sobre el peso del automóvil (variable continua) y el número de cilindros como un factor con tres niveles (4, 6 u 8) sin interacciones.

EDITAR: Si está interesado en esta pregunta, definitivamente encontrará una respuesta detallada y satisfactoria en esta publicación de Matthew Drury fuera de CV .

Nota : publiqué una versión ampliada de esta respuesta en mi sitio web .

¿Podría considerar publicar una respuesta similar con el motor R real expuesto?

¡Seguro! Bajamos por la madriguera del conejo.

La primera capa es lmla interfaz expuesta al programador R. Puede ver la fuente de esto simplemente escribiendo lmen la consola R. La mayor parte (como la mayoría de la mayoría del código de nivel de producción) está ocupado revisando entradas, configurando atributos de objeto y arrojando errores; pero esta línea sobresale

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

lm.fites otra función R, puedes llamarla tú mismo. Si lmbien funciona convenientemente con fórmulas y marcos de datos, lm.fitdesea matrices, por lo que ese es un nivel de abstracción eliminado. Verificando la fuente lm.fit, más trabajo ocupado y la siguiente línea realmente interesante

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

Ahora estamos llegando a alguna parte. .Calles la forma en que R llama al código C. Hay una función C, C_Cdqrls en la fuente R en alguna parte, y necesitamos encontrarla. Aqui esta .

Mirando la función C, de nuevo, encontramos en su mayoría comprobación de límites, limpieza de errores y trabajo ocupado. Pero esta linea 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 idioma, R ha llamado 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 finalmente estamos en el punto donde podemos hacer algo de álgebra lineal, y en realidad resolver el sistema de ecuaciones. Este es el tipo de cosas en las que fortran es realmente bueno, lo que explica por qué pasamos por tantas capas para llegar aquí.

El comentario también explica qué 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.

Entonces fortran va a resolver el sistema encontrando la descomposición .$QR$

Lo primero que sucede, 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 dqrdc2en nuestra matriz de entrada x. ¿Qué es esto?

 c     dqrfit uses the linpack routines dqrdc and dqrsl.

Así que finalmente hemos llegado a linpack . Linpack es una biblioteca de álgebra lineal fortran que existe desde los años 70. El álgebra lineal más grave eventualmente encuentra su camino hacia 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 realiza el trabajo real. Me llevaría un buen día completo descubrir qué está haciendo este código, es tan bajo como es posible. Pero genéricamente, tenemos una matriz y queremos factorizarla en un producto X = Q R donde Q es una matriz ortogonal y R es una matriz triangular superior. Es algo inteligente, porque una vez que tienes Q y R puedes resolver las ecuaciones lineales para la regresión$X$$X = QR$$Q$$R$$Q$$R$

muy facilmente. En efecto

entonces todo el sistema se convierte

pero es triangular superior y tiene el mismo rango que X t X , por lo que mientras nuestro problema esté bien planteado, es un rango completo, y también podríamos resolver el sistema reducido$R$$X^t X$

Pero aquí está lo asombroso. es triangular superior, por lo que la última ecuación lineal aquí es justa , por lo que la resolución de β n es trivial. Luego puede subir las filas, una por una, y sustituirlas en las β que ya conoce, obteniendo cada vez una ecuación lineal variable simple para resolver. Entonces, una vez que tienes Q y R , todo se derrumba a lo que se llama sustitución hacia atrás , lo cual es fácil. Puede leer sobre esto con más detalle aquí , donde se elabora un pequeño ejemplo explícito.$R$constant * beta_n = constant$\beta_n$$\beta$$Q$$R$

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 probarse a sí mismo que los resultados en R con un ejemplo simple se pueden alcanzar siguiendo el álgebra lineal de proyecciones en el espacio de la columna y el concepto de errores perpendiculares (producto de puntos), ilustrado en diferentes publicaciones , y bien explicado por el Dr. Strang en Álgebra lineal y sus aplicaciones , y fácilmente accesible aquí .

$\small \beta$

$\small D1$$\small D2$$\small X$

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

$\small [*]$lm

Continuando entonces, para calcular los coeficientes ($\small\beta$$\small ProjMatrix = \small (X^{T}X)^{-1}X^{T}$$\small [ProjMatrix] \,[y]\, =\, [RegrCoef's]$$\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éntico a: coef(lm(mpg ~ wt + as.factor(cyl)-1)) .

$\small Hat Matrix = \small X(X^{T}X)^{-1}X^{T}$

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

$\hat{y}$$\small X(X^{T}X)^{-1}X^{T}\,y$y_hat <- HAT %*% mpg

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