注意：我已经在我的网站上发布了这个答案的扩展版本。

您会考虑在实际暴露的 R 引擎的情况下发布类似的答案吗？

当然！我们去兔子洞。

第一层是lm暴露给 R 程序员的接口。lm您只需在 R 控制台上键入即可查看源代码。其中大部分（就像大多数生产级代码一样）忙于检查输入、设置对象属性和抛出错误；但是这条线很突出

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

lm.fit是另一个R函数，你可以自己调用。虽然lm可以方便地使用公式和数据框，lm.fit但需要矩阵，因此删除了一层抽象。检查源代码lm.fit，更繁忙的工作，以及以下非常有趣的行

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

现在我们正在取得进展。 .Call是 R 调用 C 代码的方式。在某个地方的 R 源代码中有一个 C 函数 C_Cdqrls，我们需要找到它。 在这里。

再次查看 C 函数，我们发现主要是边界检查、错误清理和繁忙的工作。但这条线不同

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

所以现在我们在使用我们的第三种语言，R 调用了 C，它调用了 fortran。 这是fortran代码。

第一条评论说明了一切

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

（有趣的是，这个例程的名称似乎在某些时候被更改了，但有人忘记更新评论）。所以我们终于可以做一些线性代数，并实际求解方程组。这是 fortran 真正擅长的事情，这就解释了为什么我们要经过这么多层才能到达这里。

该评论还解释了代码将要做什么

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

因此，fortran 将通过查找分解来解决系统问题。$QR$

发生的第一件事，也是迄今为止最重要的，是

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

dqrdc2这会在我们的输入矩阵上调用 fortran 函数x。这是什么？

 c     dqrfit uses the linpack routines dqrdc and dqrsl.

所以我们终于做到了linpack。Linpack 是一个自 70 年代以来一直存在的 fortran 线性代数库。大多数严肃的线性代数最终都找到了 linpack 的方式。在我们的例子中，我们使用函数dqrdc2

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

这是完成实际工作的地方。我要花一整天的时间才能弄清楚这段代码在做什么，它的级别很低。但一般来说，我们有一个矩阵，我们想将其分解为乘积，其中是正交矩阵，是上三角矩阵。这是一件很聪明的事情，因为一旦有了和，您就可以求解线性方程以进行回归$X$$X = QR$$Q$$R$$Q$$R$

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

非常简单地。确实

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

所以整个系统变成

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

但是具有相同的秩，所以只要我们的问题是良定的，它就是满秩的，我们还不如解决约化系统$R$$X^t X$

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

但这是很棒的事情。 是上三角函数，所以这里的最后一个线性方程就是，所以求解很简单。然后，您可以逐行向上，并代入您已经知道的 s，每次都得到一个简单的单变量线性方程来求解。所以，一旦你有了和，整个事情就会崩溃到所谓的向后替换，这很容易。您可以在此处更详细地阅读此内容，其中一个明确的小示例已完全解决。$R$constant * beta_n = constant$\beta_n$$\beta$$Q$$R$

R 中的实际逐步计算在 Matthew Drury 在同一线程中的回答中得到了精美的描述。在这个答案中，我想通过一个简单的例子向自己证明 R 中的结果可以通过投影到列空间的线性代数和垂直（点积）误差概念来达到，在不同的帖子中进行了说明，并由 Strang 博士在线性代数及其应用中进行了很好的解释，并且在此处易于访问。

为了估计回归中的系数，$\small \beta$

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

用和表示值为 [0,1] 的虚拟变量，我们首先需要在设计矩阵中包含($\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没有“-1”的情况下，它需要第一列只有一个，但是如果没有这个截距列，我们将得出相同的结果。

然后继续计算系数（$\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

等同于：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