Rebuilding summary(lm()) with matrix algebra
An R function that turns a formula into a design matrix and reproduces the summary(lm()) coefficient table with matrix operations, and a test of how much precision the explicit inverse loses next to a QR decomposition.
In October 2020 I wrote an R function that fits linear models using only matrix and scalar operations. It takes a formula and a data frame, builds the matrix of predictors, and returns the coefficient table that summary(lm()) prints. In December I tidied it and put it on my portfolio site next to a hand-written anova().
The algebra
Turn the formula into a design matrix X, with a column of 1s for the intercept and one column per predictor, and pull out the response vector y. The coefficients solve the normal equations: β̂ = (XᵀX)⁻¹Xᵀy.
With n rows and p columns in X, the residuals are e = y − Xβ̂ and σ̂² = eᵀe / (n − p). The divisor n − p is the residual degrees of freedom. The covariance matrix of β̂ is σ̂²(XᵀX)⁻¹, and the standard errors are the square roots of its diagonal. Each t statistic is β̂ⱼ divided by its standard error, and the two-sided p-value comes from a t distribution with n − p degrees of freedom.
For the fit, R² = 1 − eᵀe / Σ(y − ȳ)². Adjusted R² and the overall F statistic follow from the same pieces.
The function
I let R’s own formula machinery build X. model.frame() evaluates the formula against the data, and model.matrix() expands factors into dummy columns, adds interactions and handles transformations, the same way lm() does. Since p counts the columns of X, data frame columns the model doesn’t use have no effect on the degrees of freedom. The rest is the algebra above:
my_lm <- function(formula, data) {
mf <- model.frame(formula, data)
y <- model.response(mf)
X <- model.matrix(attr(mf, "terms"), mf)
n <- nrow(X); p <- ncol(X); df <- n - p
XtX_inv <- solve(crossprod(X))
beta <- drop(XtX_inv %*% crossprod(X, y))
resid <- drop(y - X %*% beta)
sigma2 <- sum(resid^2) / df
se <- sqrt(diag(sigma2 * XtX_inv))
tval <- beta / se
pval <- 2 * pt(abs(tval), df, lower.tail = FALSE)
has_int <- attr(attr(mf, "terms"), "intercept") == 1
tss <- if (has_int) sum((y - mean(y))^2) else sum(y^2)
r2 <- 1 - sum(resid^2) / tss
k <- p - has_int
list(coefficients = cbind(Estimate = beta, `Std. Error` = se,
`t value` = tval, `Pr(>|t|)` = pval),
sigma = sqrt(sigma2), r.squared = r2,
adj.r.squared = 1 - (1 - r2) * (n - has_int) / df,
fstatistic = c(value = (r2 / k) / ((1 - r2) / df), numdf = k, dendf = df))
}The has_int check covers models without an intercept, where lm() measures R² against zero instead of the mean. On mpg ~ wt - 1 both give 0.7197.
Checking it against lm()
The numbers in this post come from running the code in R 4.6.1. On mtcars with mpg ~ wt + hp + factor(cyl), the output matches summary(lm()):
| Term | Estimate | Std. Error | t value |
|---|---|---|---|
| (Intercept) | 35.846 | 2.041 | 17.56 |
| wt | −3.181 | 0.720 | −4.42 |
| hp | −0.0231 | 0.0120 | −1.93 |
| factor(cyl)6 | −3.359 | 1.402 | −2.40 |
| factor(cyl)8 | −3.186 | 2.170 | −1.47 |
The largest absolute difference anywhere in the coefficient table, p-values included, is 5.8e-14. Residual standard error (2.440), R² (0.8572), adjusted R² (0.8361) and the F statistic (40.53 on 4 and 27 degrees of freedom) all agree to within 3e-14. On a gout dataset of 400 people, where race and sex are two-level text columns, su ~ race + sex + age matches to 7e-14.
Inverse versus QR
The remaining differences come from a different algorithm. lm() never forms XᵀX. It takes a QR decomposition of X and solves a triangular system. Forming XᵀX squares the condition number κ, which roughly doubles the number of digits lost to rounding.
To see how much that matters, I fit exact polynomials: 50 evenly spaced points on [0, 1], columns 1, x, x², … up to some degree, every true coefficient equal to 1, and no noise. The table shows the largest coefficient error:
| Degree | κ(X) | Explicit inverse | QR |
|---|---|---|---|
| 2 | 22 | 1.8e-14 | 2.9e-15 |
| 4 | 643 | 9.9e-12 | 3.6e-14 |
| 6 | 2.0e4 | 2.8e-8 | 1.5e-12 |
| 8 | 6.3e5 | 1.0e-4 | 9.1e-11 |
| 10 | 2.0e7 | 0.41 | 1.3e-10 |
| 12 | 6.8e8 | fails | 7.0e-8 |
By degree 10, the explicit inverse misses coefficients that should be 1 by as much as 0.41, and at degree 12 solve() refuses: “system is computationally singular.” QR still gets about seven correct digits, and lm() matched it in every row. solve(XtX, Xty), which skips the explicit inverse, does better but fails the same way at degree 12.
For a few well-scaled predictors, the normal equations are fine, as the mtcars match shows. The version I’d actually use replaces the XtX_inv and beta lines:
qx <- qr(X)
beta <- qr.coef(qx, y)
XtX_inv <- chol2inv(qr.R(qx))chol2inv(qr.R(qx)) builds (XᵀX)⁻¹ from the triangular factor, as summary.lm() does. With that change, the largest difference from lm() on the mtcars model drops from 5.8e-14 to 3.6e-15.