| SOCR ≫ | DSPA ≫ | DSPA3 Topics ≫ |
DSPA3: Data Science and Predictive Analytics (UMich HS650)
Chapter 3: Linear Algebra, Matrix Computing, and Regression
Chapter 3: Linear Algebra, Matrix Computing, and Regression
SOCR/MCAIM (Ivo Dinov)
SOCR/MCAIM (Ivo Dinov)
August 25, 2026
library(ggplot2)
library(dplyr)
library(tidyr)
library(patchwork)
library(plotly) # interactive figures and ALL 3-D graphics
library(Matrix)
library(DT)How this chapter uses graphics
Every two-dimensional figure is drawn with
ggplot2and rendered statically. Immediately after each one, the equivalentplot_ly()code appears in a chunk markedeval=FALSE, echo=TRUE, visible in these notes, ready to paste into a live session for an interactive demonstration.Every three-dimensional figure is drawn with
plot_ly()and evaluated, because rotation carries information a fixed projection destroys. In this chapter that includes the residual-sum-of-squares surface over parameter space, the covariance ellipsoid with its principal axes, the geometry of orthogonal projection, and the fitted regression plane.
1 Learning objectives
After completing this chapter you will be able to:
- Construct and manipulate matrices in
R, and state the cost of each operation in floating-point operations. - Define rank, null space, and the four fundamental subspaces, and explain why \(\operatorname{rank}(X)=p\) is the condition least squares actually requires.
- Derive the orthogonal projection \(\hat{Y}=HY\) and prove that \(H\) is symmetric and idempotent with \(\operatorname{tr}(H)=p\).
- Compute and compare the LU, Cholesky, QR, eigen, and singular value decompositions, and choose the right one for a given task.
- Quantify numerical conditioning, and explain why forming \(X^\top X\) costs you half your significant digits.
- Derive the OLS estimator three ways, calculus, geometry, and maximum likelihood, and state the Gauss–Markov theorem precisely.
- Derive \(\operatorname{Var}(\hat\beta)=\sigma^2(X^\top X)^{-1}\) and use it to construct every standard error, \(t\)-statistic, and \(F\)-test in a regression summary.
- Diagnose a fitted model using leverage, standardized residuals, and Cook’s distance, and explain each in terms of the hat matrix.
- Distinguish multicollinearity from dependence, compute VIFs, and explain why spurious correlation is guaranteed when \(k \gg n\).
- Fit and interpret regression trees, model trees, and Bayesian additive regression trees, and state the computational cost of each.
Estimated time: 9–12 hours including exercises.
Prerequisites: Chapter 1 (reproducibility
conventions, dspa_read()) and Chapter 2 (covariance, correlation, EDA).
Readers wanting a refresher on calculus, series, and complex numbers may
find the mathematical
foundations module useful.
2 PART I: LINEAR ALGEBRA AND MATRIX COMPUTING
Linear algebra is the study of linear maps between vector spaces. Its practical importance rests on a single observation: almost every process is locally linear. A smooth function \(f\) near a point \(x_0\) satisfies
\[f(x) = f(x_0) + \underbrace{Df(x_0)}_{\text{a matrix}}(x-x_0) + O\!\left(\|x-x_0\|^2\right),\]
so the first-order behaviour of any differentiable system is a matrix. That is why arbitrarily complicated problems collapse into two canonical forms:
\[\underbrace{A\mathbf{x}=\mathbf{b}}_{\text{solve a system}} \qquad\text{and}\qquad \underbrace{A\mathbf{v}=\lambda\mathbf{v}}_{\text{find invariant directions}}.\]
Everything in this chapter is an elaboration of those two equations.
3 Building and indexing matrices
3.1 Construction
#> [,1] [,2] [,3]
#> [1,] 1 3 5
#> [2,] 2 4 6#> [,1] [,2] [,3]
#> [1,] 1 2 3
#> [2,] 4 5 6#> [,1] [,2] [,3] [,4]
#> [1,] -0.591 -0.934 -0.828 0.012
#> [2,] 0.027 1.324 -0.348 -0.223
#> [3,] -1.517 0.625 -1.538 0.888
#> [4,] -1.363 -0.046 -0.256 -0.592
#> [5,] 1.178 -1.004 -1.150 -0.656R stores matrices in column-major
order, a matrix is a vector with a dim attribute,
laid out one column after another. This is the same convention as
Fortran and LAPACK, and the reason column operations
(colSums, X[, j]) touch contiguous memory and
are faster than the corresponding row operations on large matrices.
#> [1] 1 2 3 4 5 6#> $dim
#> [1] 2 3diag() is overloaded three ways, keyed on the type of
its argument:
#> [,1] [,2] [,3]
#> [1,] 1 0 0
#> [2,] 0 2 0
#> [3,] 0 0 3#> [1] 1 4#> [,1] [,2] [,3] [,4]
#> [1,] 1 0 0 0
#> [2,] 0 1 0 0
#> [3,] 0 0 1 0
#> [4,] 0 0 0 1c1 <- 1:5
m4 <- cbind(m3, c1) # column concatenation; picks up the name "c1"
r1 <- 1:4
m5 <- rbind(m3, r1)
dimnames(m5) <- list(NULL, NULL) # strip inherited names
round(m5, 3)#> [,1] [,2] [,3] [,4]
#> [1,] -0.591 -0.934 -0.828 0.012
#> [2,] 0.027 1.324 -0.348 -0.223
#> [3,] -1.517 0.625 -1.538 0.888
#> [4,] -1.363 -0.046 -0.256 -0.592
#> [5,] 1.178 -1.004 -1.150 -0.656
#> [6,] 1.000 2.000 3.000 4.0003.2 Indexing
#> [,1] [,2] [,3] [,4]
#> [1,] 1 4 7 10
#> [2,] 2 5 8 11
#> [3,] 3 6 9 12#> [1] 4#> [1] 1 4 7 10#> [,1] [,2] [,3] [,4]
#> [1,] 1 4 7 10#> [,1] [,2]
#> [1,] 4 7
#> [2,] 5 8
#> [3,] 6 9Common misconception: “
A[i, ]gives me a row of a matrix.” It gives a vector.Rsilently drops the dimension when a subscript selects a single row or column, sodim(m6[1, ])isNULLand any subsequent%*%will treat it as a column vector. In code that must work for both a single row and several, always writedrop = FALSE. This is one of the most common sources of “non-conformable arguments” errors.
4 Matrix algebra
4.1 Addition, subtraction, and scalar recycling
Addition and subtraction are elementwise and require identical dimensions.
#> [,1] [,2] [,3]
#> [1,] 3 7 11
#> [2,] 5 9 13#> [,1] [,2] [,3]
#> [1,] 1 1 1
#> [2,] 1 1 1#> [,1] [,2] [,3]
#> [1,] 1 3 5
#> [2,] 2 4 6Cost: \(O(mn)\), one operation per entry, no arithmetic can avoid it.
4.2 Elementwise product versus matrix product
R uses * for the Hadamard
(elementwise) product and %*% for the
matrix product. They are entirely different operations
that happen to share notation in loose writing.
\[(A\circ B)_{ij}=A_{ij}B_{ij} \qquad\text{versus}\qquad (AB)_{ij}=\sum_{\ell=1}^{n}A_{i\ell}B_{\ell j}.\]
The Hadamard product needs \(A\) and \(B\) to have the same shape. The matrix product needs the inner dimensions to match: \(P_{m\times k}=L_{m\times n}R_{n\times k}\).
#> [,1] [,2] [,3]
#> [1,] 2 12 30
#> [2,] 6 20 42#> [1] 2 3#> [1] 3 2#> [,1] [,2]
#> [1,] 52 88
#> [2,] 64 109Complexity. The definition requires \(mnk\) multiply–add pairs, so the textbook
algorithm is \(O(mnk)\), or \(O(n^3)\) for square matrices. This is not
optimal: Strassen’s
algorithm (1969) achieves \(O(n^{\log_2
7})=O(n^{2.807})\) by trading one multiplication for several
additions, and successive refinements have pushed the theoretical
exponent below \(2.372\). In practice
R calls a tuned BLAS (dgemm), which uses the
\(O(n^3)\) algorithm with cache
blocking, asymptotically worse but far faster at realistic sizes,
because memory traffic dominates flop count.
set.seed(3)
A <- matrix(rnorm(500 * 500), 500)
B <- matrix(rnorm(500 * 500), 500)
system.time(A %*% B)[["elapsed"]] # BLAS dgemm#> [1] 0.034.3 Inner and outer products
Two vectors combine in two fundamentally different ways.
\[\underbrace{\mathbf{u}^\top\mathbf{v}=\sum_i u_iv_i}_{\text{inner product: a SCALAR}} \qquad\qquad \underbrace{\mathbf{u}\mathbf{v}^\top=\big(u_iv_j\big)}_{\text{outer product: a MATRIX of rank 1}}\]
u <- c(1, 2, 3, 4, 5)
v <- c(4, 5, 6, 7, 8)
u %*% v # 1x1: R silently makes the left operand a ROW vector#> [,1]
#> [1,] 100#> [,1]
#> [1,] 100#> [1] 100#> [,1] [,2] [,3] [,4] [,5]
#> [1,] 4 5 6 7 8
#> [2,] 8 10 12 14 16
#> [3,] 12 15 18 21 24
#> [4,] 16 20 24 28 32
#> [5,] 20 25 30 35 40#> [,1] [,2] [,3] [,4] [,5]
#> [1,] 4 5 6 7 8
#> [2,] 8 10 12 14 16
#> [3,] 12 15 18 21 24
#> [4,] 16 20 24 28 32
#> [5,] 20 25 30 35 40#> [1] 1Common misconception: “multiplying two vectors gives the outer product.” It depends entirely on orientation. \(\mathbf{u}^\top\mathbf{v}\) is a scalar measuring alignment; \(\mathbf{u}\mathbf{v}^\top\) is an \(n\times n\) matrix of rank one.
Rblurs this because a bare vector has no orientation, inu %*% vit is treated as a row on the left and a column on the right. Writecrossprod(u, v)for the inner product andtcrossprod(u, v)for the outer product, and the ambiguity disappears.
Rank-one outer products are the atoms of matrix analysis: the singular value decomposition (§3.6.5) writes any matrix as a weighted sum of them.
4.4 Transpose,
crossprod, and why it matters
#> [,1] [,2] [,3]
#> [1,] 2 4 6
#> [2,] 3 5 7#> [,1] [,2]
#> [1,] 2 3
#> [2,] 4 5
#> [3,] 6 7#> [1] TRUEThe combination \(A^\top B\) is so
ubiquitous in statistics that R provides a dedicated
routine:
set.seed(5)
X <- matrix(rnorm(2000 * 200), 2000)
t1 <- system.time(r1 <- t(X) %*% X)[["elapsed"]]
t2 <- system.time(r2 <- crossprod(X))[["elapsed"]]
c(explicit_transpose = t1, crossprod = t2, identical = all.equal(r1, r2))#> explicit_transpose crossprod identical
#> 0.01 0.02 1.00crossprod(X) never materializes \(X^\top\), it calls BLAS dsyrk,
which exploits the symmetry of the result and computes only the lower
triangle, roughly halving the work. Use crossprod(A, B) for
\(A^\top B\) and
tcrossprod(A, B) for \(AB^\top\).
5 Vector spaces, rank, and the four subspaces
Everything about the solvability of \(A\mathbf{x}=\mathbf{b}\) and the identifiability of a regression model is a statement about subspaces.
5.1 Span, independence, basis
Given vectors \(\mathbf{a}_1,\dots,\mathbf{a}_n\in\mathbb{R}^m\), their span is the set of all linear combinations,
\[\operatorname{span}\{\mathbf{a}_1,\dots,\mathbf{a}_n\}=\left\{\sum_{j=1}^n c_j\mathbf{a}_j : c_j\in\mathbb{R}\right\}.\]
They are linearly independent if \(\sum_j c_j\mathbf{a}_j=\mathbf{0}\) forces every \(c_j=0\). A basis for a subspace is an independent spanning set; every basis of a given subspace has the same size, called its dimension.
5.2 The four fundamental subspaces
For \(A\in\mathbb{R}^{m\times n}\):
| Subspace | Definition | Lives in | Dimension |
|---|---|---|---|
| Column space, aka Image space, \(\operatorname{im}(A)=\operatorname{col}(A)\) | \(\{A\mathbf{x}:\mathbf{x}\in\mathbb{R}^n\}\) | \(\mathbb{R}^m\) | \(r\) |
| Null space, Kernel space, \(\operatorname{ker(A)}\mathcal{N}(A)\) | \(\{\mathbf{x}:A\mathbf{x}=\mathbf{0}\}\) | \(\mathbb{R}^n\) | \(n-r\) |
| Row space \(\operatorname{col}(A^\top)\) | \(\{A^\top\mathbf{y}:\mathbf{y}\in\mathbb{R}^m\}\) | \(\mathbb{R}^n\) | \(r\) |
| Left null space \(\mathcal{N}(A^\top)\) | \(\{\mathbf{y}:A^\top\mathbf{y}=\mathbf{0}\}\) | \(\mathbb{R}^m\) | \(m-r\) |
Here \(r=\operatorname{rank}(A)\), the number of linearly independent columns, which equals the number of linearly independent rows, a fact that is not obvious and is worth pausing on.
Rank–nullity theorem. \[\operatorname{rank}(A)+\dim\mathcal{N}(A)=n.\]
Orthogonality relations. \(\mathcal{N}(A)\perp\operatorname{col}(A^\top)\) and \(\mathcal{N}(A^\top)\perp\operatorname{col}(A)\). The second is the one that generates least squares (§3.9): the residual lives in the left null space, orthogonal to everything the model can produce.
#> [,1] [,2] [,3]
#> [1,] 1 2 3
#> [2,] 2 4 6
#> [3,] 1 1 1#> [1] 2# A basis for the null space, from the SVD's trailing right singular vectors
sv <- svd(A)
r <- sum(sv$d > max(dim(A)) * .Machine$double.eps * max(sv$d))
null_basis <- sv$v[, (r + 1):ncol(A), drop = FALSE]
round(null_basis, 4)#> [,1]
#> [1,] 0.4082
#> [2,] -0.8165
#> [3,] 0.4082#> [,1]
#> [1,] 0
#> [2,] 0
#> [3,] 0Why this matters for regression. A design matrix
\(X_{n\times p}\) produces a unique OLS
solution iff \(\operatorname{rank}(X)=p\), i.e. \(\mathcal{N}(X)=\{\mathbf{0}\}\). If a
column is an exact linear combination of others, a dummy-variable trap,
a redundant total, \(n<p\), then
infinitely many \(\hat\beta\) give
identical fitted values, and no data can distinguish them.
lm() responds by reporting NA for the aliased
coefficients rather than failing.
set.seed(7)
d <- data.frame(x1 = rnorm(50), x2 = rnorm(50))
d$x3 <- 2 * d$x1 - 3 * d$x2 # exactly collinear by construction
d$y <- d$x1 + d$x2 + rnorm(50)
coef(lm(y ~ x1 + x2 + x3, data = d)) # x3 is aliased -> NA#> (Intercept) x1 x2 x3
#> 0.173273 1.081103 0.808563 NA6 Norms, inner products, and orthogonal projection
6.1 Norms
A norm assigns a length. The \(\ell_p\) family on \(\mathbb{R}^n\):
\[\|\mathbf{x}\|_p=\left(\sum_{i=1}^n|x_i|^p\right)^{1/p},\qquad \|\mathbf{x}\|_1=\sum_i|x_i|,\quad \|\mathbf{x}\|_2=\sqrt{\sum_i x_i^2},\quad \|\mathbf{x}\|_\infty=\max_i|x_i|.\]
The choice of norm is a modelling decision with consequences you will meet repeatedly: \(\ell_2\) gives ridge regression and OLS, \(\ell_1\) gives LASSO and median regression, and the difference between them is the reason \(\ell_1\) produces exact zeros.
The inner product \(\langle\mathbf{u},\mathbf{v}\rangle=\mathbf{u}^\top\mathbf{v}\) induces \(\|\mathbf{x}\|_2=\sqrt{\langle\mathbf{x},\mathbf{x}\rangle}\) and satisfies the Cauchy–Schwarz inequality
\[|\langle\mathbf{u},\mathbf{v}\rangle|\le\|\mathbf{u}\|_2\,\|\mathbf{v}\|_2,\]
with equality iff \(\mathbf{u}\) and \(\mathbf{v}\) are parallel. Applying this to centred data vectors is exactly the statement \(|\rho_{xy}|\le 1\), the correlation coefficient is the cosine of the angle between centred variables:
\[\rho_{xy}=\frac{\langle \mathbf{x}-\bar x\mathbf{1},\; \mathbf{y}-\bar y\mathbf{1}\rangle} {\|\mathbf{x}-\bar x\mathbf{1}\|_2\,\|\mathbf{y}-\bar y\mathbf{1}\|_2}=\cos\theta .\]
set.seed(13)
x <- rnorm(100); y <- 0.6 * x + rnorm(100, sd = 0.8)
xc <- x - mean(x); yc <- y - mean(y)
c(cor = cor(x, y),
cosine_angle = sum(xc * yc) / (sqrt(sum(xc^2)) * sqrt(sum(yc^2))),
angle_degrees = acos(cor(x, y)) * 180 / pi)#> cor cosine_angle angle_degrees
#> 0.518186 0.518186 58.7893516.2 Orthogonal projection
Let \(S=\operatorname{col}(X)\subseteq\mathbb{R}^n\) with \(X\) of full column rank \(p\). The orthogonal projection of \(\mathbf{y}\) onto \(S\) is the unique point \(\hat{\mathbf{y}}\in S\) minimizing \(\|\mathbf{y}-\hat{\mathbf{y}}\|_2\). Writing \(\hat{\mathbf{y}}=X\mathbf{b}\) and requiring the residual to be orthogonal to every column of \(X\),
\[X^\top(\mathbf{y}-X\mathbf{b})=\mathbf{0} \;\Longrightarrow\; \mathbf{b}=(X^\top X)^{-1}X^\top\mathbf{y} \;\Longrightarrow\; \hat{\mathbf{y}}=\underbrace{X(X^\top X)^{-1}X^\top}_{H}\mathbf{y}.\]
\(H\) is the hat matrix, it puts the hat on \(\mathbf{y}\).
Proposition. \(H\) is symmetric, idempotent, and \(\operatorname{tr}(H)=p\). Likewise \(M=I-H\) is symmetric idempotent with \(\operatorname{tr}(M)=n-p\), and \(HM=0\).
Proof. Symmetry: \(H^\top=X\big((X^\top X)^{-1}\big)^\top X^\top=H\) since \(X^\top X\) is symmetric. Idempotence: \(H^2=X(X^\top X)^{-1}\underbrace{X^\top X(X^\top X)^{-1}}_{I}X^\top=H\). Trace: using \(\operatorname{tr}(AB)=\operatorname{tr}(BA)\), \(\operatorname{tr}(H)=\operatorname{tr}\!\big((X^\top X)^{-1}X^\top X\big)=\operatorname{tr}(I_p)=p\). \(\blacksquare\)
These three facts carry an enormous amount of downstream weight: the \(p\) in “degrees of freedom”, the \(n-p\) divisor in \(\hat\sigma^2\), the diagonal entries \(h_{ii}\) that define leverage, and the \(F\)-test’s numerator and denominator degrees of freedom are all consequences.
set.seed(19)
n <- 60; p <- 3
Xd <- cbind(1, matrix(rnorm(n * (p - 1)), n))
H <- Xd %*% solve(crossprod(Xd)) %*% t(Xd)
c(symmetric = all.equal(H, t(H)),
idempotent = all.equal(H %*% H, H),
trace_H = round(sum(diag(H)), 10),
p = p,
trace_I_H = round(sum(diag(diag(n) - H)), 10),
n_minus_p = n - p)#> symmetric idempotent trace_H p trace_I_H n_minus_p
#> 1 1 3 3 57 57Projection is a geometric statement, so it deserves a geometric picture. This one is interactive, the perpendicularity of the residual to the plane is the entire content of least squares, and it is only convincing when you can rotate it.
set.seed(23)
# A plane through the origin spanned by two vectors in R^3
Xg <- cbind(c(1, 0, 0.3), c(0, 1, 0.5))
Hg <- Xg %*% solve(crossprod(Xg)) %*% t(Xg)
yv <- c(0.9, 0.4, 2.2) # a point off the plane
yhat <- as.vector(Hg %*% yv) # its projection
resid <- yv - yhat
gs <- seq(-1.5, 1.5, length.out = 25)
grid <- expand.grid(a = gs, b = gs)
px <- matrix(Xg[1, 1] * grid$a + Xg[1, 2] * grid$b, 25)
py <- matrix(Xg[2, 1] * grid$a + Xg[2, 2] * grid$b, 25)
pz <- matrix(Xg[3, 1] * grid$a + Xg[3, 2] * grid$b, 25)
plot_ly() |>
add_surface(x = px, y = py, z = pz, opacity = 0.45, showscale = FALSE,
colorscale = list(c(0, "#BBD5EE"), c(1, "#BBD5EE"))) |>
add_trace(x = c(0, yv[1]), y = c(0, yv[2]), z = c(0, yv[3]),
type = "scatter3d", mode = "lines+markers", name = "y",
line = list(width = 6, color = "black"),
marker = list(size = 4, color = "black")) |>
add_trace(x = c(0, yhat[1]), y = c(0, yhat[2]), z = c(0, yhat[3]),
type = "scatter3d", mode = "lines+markers", name = "Hy (fitted)",
line = list(width = 6, color = "#1F77B4"),
marker = list(size = 4, color = "#1F77B4")) |>
add_trace(x = c(yhat[1], yv[1]), y = c(yhat[2], yv[2]), z = c(yhat[3], yv[3]),
type = "scatter3d", mode = "lines", name = "residual (I-H)y",
line = list(width = 6, color = "firebrick", dash = "dash")) |>
layout(title = "Least squares is orthogonal projection onto col(X)",
scene = list(xaxis = list(title = "e1"), yaxis = list(title = "e2"),
zaxis = list(title = "e3"), aspectmode = "cube"))c(residual_dot_col1 = round(sum(resid * Xg[, 1]), 12),
residual_dot_col2 = round(sum(resid * Xg[, 2]), 12),
pythagoras = all.equal(sum(yv^2), sum(yhat^2) + sum(resid^2)))#> residual_dot_col1 residual_dot_col2 pythagoras
#> 0 0 1The residual is numerically orthogonal to both spanning vectors, and \(\|\mathbf{y}\|^2=\|\hat{\mathbf{y}}\|^2+\|\mathbf{y}-\hat{\mathbf{y}}\|^2\). That Pythagorean identity, applied to centred data, is the ANOVA decomposition \(\mathrm{SST}=\mathrm{SSR}+\mathrm{SSE}\) (§3.10).
7 Inverses, pseudo-inverses, and solving systems
7.1 When does an inverse exist?
For a square \(A\in\mathbb{R}^{n\times n}\), a two-sided inverse \(A^{-1}\) satisfying \(AA^{-1}=A^{-1}A=I_n\) exists iff \(A\) is nonsingular, i.e. \(\operatorname{rank}(A)=n\), equivalently \(\det(A)\ne 0\), equivalently \(\mathcal{N}(A)=\{\mathbf{0}\}\). When it exists it is unique.
For a rectangular \(A\in\mathbb{R}^{m\times n}\) one-sided inverses may exist:
- a left inverse \(L\) with \(LA=I_n\) exists iff \(\operatorname{rank}(A)=n\) (full column rank), and then \(L=(A^\top A)^{-1}A^\top\) works;
- a right inverse \(R\) with \(AR=I_m\) exists iff \(\operatorname{rank}(A)=m\) (full row rank), and then \(R=A^\top(AA^\top)^{-1}\) works.
Common misconception: “\(A^{-1}\) means the same thing for rectangular matrices.” One-sided inverses are not unique when \(m\ne n\), there are infinitely many left inverses of a tall full-column-rank matrix. Writing \(A^{-1}\) for one of them invites the reader to cancel it on either side, which is invalid. Use \(A^{+}\) (the pseudo-inverse) when you need a canonical choice.
For \(2\times 2\) the closed form is worth memorizing:
\[A=\begin{pmatrix}a&b\\c&d\end{pmatrix},\qquad A^{-1}=\frac{1}{ad-bc}\begin{pmatrix}d&-b\\-c&a\end{pmatrix},\qquad ad-bc\ne 0 .\]
#> [,1] [,2]
#> [1,] 1 3
#> [2,] 2 4#> [1] -2#> [,1] [,2]
#> [1,] -2 1.5
#> [2,] 1 -0.5#> [,1] [,2]
#> [1,] 1 0
#> [2,] 0 17.2 The Moore–Penrose pseudo-inverse
For any \(A\in\mathbb{R}^{m\times n}\), singular or rectangular, there is a unique \(A^{+}\in\mathbb{R}^{n\times m}\) satisfying the four Penrose conditions
\[AA^{+}A=A,\qquad A^{+}AA^{+}=A^{+},\qquad (AA^{+})^\top=AA^{+},\qquad (A^{+}A)^\top=A^{+}A .\]
It is computed from the SVD (§3.6.5) as \(A^{+}=V\Sigma^{+}U^\top\), where \(\Sigma^{+}\) inverts the nonzero singular values and transposes. When \(A\) is invertible, \(A^{+}=A^{-1}\); when \(A\) has full column rank, \(A^{+}=(A^\top A)^{-1}A^\top\).
\(A^{+}\mathbf{b}\) is the minimum-norm least-squares solution: among all \(\mathbf{x}\) minimizing \(\|A\mathbf{x}-\mathbf{b}\|_2\), it is the one with smallest \(\|\mathbf{x}\|_2\). That property is what makes it the right default for rank-deficient problems.
Asing <- matrix(c(1, 2, 3,
2, 4, 6,
1, 1, 1), nrow = 3, byrow = TRUE) # rank 2
Ap <- MASS::ginv(Asing)
c(penrose_1 = all.equal(Asing %*% Ap %*% Asing, Asing),
penrose_2 = all.equal(Ap %*% Asing %*% Ap, Ap),
penrose_3 = all.equal(Asing %*% Ap, t(Asing %*% Ap)),
penrose_4 = all.equal(Ap %*% Asing, t(Ap %*% Asing)))#> penrose_1 penrose_2 penrose_3 penrose_4
#> TRUE TRUE TRUE TRUE#> [1] "try-error"7.3 Solving \(A\mathbf{x}=\mathbf{b}\)
Consider the system
\[\begin{aligned} a + b + 2c &= 6\\ 3a - 2b + c &= 2\\ 2a + b - c &= 3 \end{aligned} \qquad\Longleftrightarrow\qquad \underbrace{\begin{pmatrix}1&1&2\\3&-2&1\\2&1&-1\end{pmatrix}}_{A} \underbrace{\begin{pmatrix}a\\b\\c\end{pmatrix}}_{\mathbf{x}} =\underbrace{\begin{pmatrix}6\\2\\3\end{pmatrix}}_{\mathbf{b}} .\]
This generalizes the scalar case exactly. To solve \(2x-3=5\), move the constant to the right (\(2x=8\)) and multiply by the reciprocal (\(x=4\)). The matrix version replaces “reciprocal” with “inverse”:
A <- matrix(c(1, 1, 2,
3, -2, 1,
2, 1, -1), nrow = 3, byrow = TRUE)
b <- c(6, 2, 3)
x_solve <- solve(A, b) # PREFERRED: LU factorization, never forms A^{-1}
x_inv <- solve(A) %*% b # forms the explicit inverse first
rbind(`solve(A, b)` = x_solve, `solve(A) %*% b` = as.vector(x_inv))#> [,1] [,2] [,3]
#> solve(A, b) 1.35 1.75 1.45
#> solve(A) %*% b 1.35 1.75 1.45#> [,1]
#> [1,] 0
#> [2,] 0
#> [3,] 0Common misconception: “to solve \(A\mathbf{x}=\mathbf{b}\), compute \(A^{-1}\) and multiply.” Mathematically fine; computationally wrong.
solve(A, b)performs an LU factorization and two triangular solves, about \(\tfrac{2}{3}n^3\) flops.solve(A) %*% bcomputes the full inverse (\(\approx 2n^3\) flops, three times the work), then multiplies. It is also less accurate: the explicit inverse accumulates rounding error in every one of its \(n^2\) entries, and the subsequent multiplication propagates all of them. The rule in numerical linear algebra is: never compute an inverse you only intend to multiply by something.
Lower-Upper (LU) factorization, aka LU decomposition, takes a matrix \(A\) and breaks it down into two special, easier-to-work-with matrices
- \(L\) (Lower Triangular Matrix), where all entries above the main diagonal are zero.
- \(U\) (Upper Triangular Matrix), where all entries below the main diagonal are zero.
Then, the original matrix equation \(A \cdot x = b\) is transformed into \[(L \cdot U) \cdot x = b\]
Instead of solving the complex system all at once, LU decomposition splits \(A\) into \(L\) and \(U\) supporting computational solving the problem in two simpler, highly efficient steps
- Forward Substitution (\(Ly = b\)): Because \(L\) is lower triangular, the computer can easily solve for a temporary vector \(y\) from the top down.
- Back Substitution (\(Ux = y\)): Because \(U\) is upper triangular, the computer can easily solve for your final answer \(x\) from the bottom up.
The main benefit is gain in speed, as doing full Gaussian elimination on a large matrix is computationally expensive. By factoring \(A\) into \(L\) and \(U\) once, we reuse that factorization to solve \(Ax = b\) very quickly for many different vectors \(b\).
set.seed(31)
n <- 700
Ab <- matrix(rnorm(n * n), n); bb <- rnorm(n)
t_solve <- system.time(x1 <- solve(Ab, bb))[["elapsed"]]
t_inv <- system.time(x2 <- solve(Ab) %*% bb)[["elapsed"]]
c(solve_A_b = t_solve, invert_then_multiply = t_inv,
speedup = round(t_inv / t_solve, 2),
residual_solve = norm(Ab %*% x1 - bb, "2"),
residual_inv = norm(Ab %*% x2 - bb, "2"))#> solve_A_b invert_then_multiply speedup
#> 5.00000e-02 1.50000e-01 3.00000e+00
#> residual_solve residual_inv
#> 3.69212e-12 7.59754e-12And a word on Cramer’s rule. It expresses \(x_j=\det(A_j)/\det(A)\), where \(A_j\) replaces column \(j\) of \(A\) by \(\mathbf{b}\). It is a beautiful theoretical result and a terrible algorithm: cofactor expansion of a determinant is \(O(n!)\), and even computing the \(n+1\) determinants by LU costs \(O(n^4)\) against \(O(n^3)\) for solving directly. It is also numerically unstable. Use it to prove things, never to compute them.
7.4 Reference: core matrix operations
| Expression | Meaning | Cost |
|---|---|---|
t(X) |
transpose | \(O(mn)\) (often free, a view) |
diag(x) |
diagonal matrix / extract diagonal / identity | \(O(n)\) / \(O(n)\) / \(O(n^2)\) |
A %*% B |
matrix product | \(O(mnk)\) |
crossprod(A, B) |
\(A^\top B\) without forming \(A^\top\) | \(O(mnk)\), halved when \(A=B\) |
tcrossprod(A, B) |
\(AB^\top\) | \(O(mnk)\) |
solve(A, b) |
solve \(A\mathbf{x}=\mathbf{b}\) via LU | \(\tfrac{2}{3}n^3\) |
solve(A) |
explicit inverse | \(2n^3\) |
qr.solve(A, b) |
least squares via QR | \(2mn^2-\tfrac{2}{3}n^3\) |
chol(A) |
Cholesky, \(A\) symmetric positive definite | \(\tfrac{1}{3}n^3\) |
eigen(A) |
eigendecomposition | \(O(n^3)\), \(\approx 10n^3\) symmetric |
svd(A) |
singular value decomposition | \(O(mn^2)\), \(\approx 14mn^2\) thin |
rowSums / colSums / rowMeans
/ colMeans |
margins | \(O(mn)\) |
8 Matrix decompositions
A decomposition rewrites \(A\) as a product of matrices with special structure, triangular, orthogonal, diagonal. The point is always the same: the factored form makes a hard problem easy.
8.1 LU
Gaussian elimination with partial pivoting factors \(PA=LU\) with \(P\) a permutation, \(L\) unit lower-triangular, \(U\) upper-triangular. Then \(A\mathbf{x}=\mathbf{b}\) becomes two triangular solves, each \(O(n^2)\):
\[L\mathbf{z}=P\mathbf{b}\ \ (\text{forward}),\qquad U\mathbf{x}=\mathbf{z}\ \ (\text{back}).\]
Factorization costs \(\tfrac{2}{3}n^3\); each subsequent right-hand side costs only \(2n^2\). This is why you factor once and solve many times.
#> [,1] [,2] [,3]
#> [1,] 1.0000 0.0000 0
#> [2,] 0.6667 1.0000 0
#> [3,] 0.3333 0.7143 1#> [,1] [,2] [,3]
#> [1,] 3 -2.0000 1.0000
#> [2,] 0 2.3333 -1.6667
#> [3,] 0 0.0000 2.8571#> [1] TRUE8.2 Cholesky decomposition
If \(A\) is symmetric positive definite (\(\mathbf{x}^\top A\mathbf{x}>0\) for all \(\mathbf{x}\ne 0\)), then \(A=R^\top R\) with \(R\) upper-triangular and positive diagonal. Cholesky decomposition exploits symmetry to halve LU’s work: \(\tfrac{1}{3}n^3\) flops, no pivoting needed, and it is unconditionally stable.
Cholesky doubles as a positive-definiteness test: it
fails exactly when the matrix is not positive definite. Every covariance
matrix in this book is symmetric positive semi-definite, which is why
chol() appears constantly in simulation and Gaussian
likelihood code.
set.seed(37)
Z <- matrix(rnorm(200 * 4), 200)
S <- crossprod(Z) / 199 # a sample covariance: symmetric PSD
R <- chol(S)
round(R, 4)#> [,1] [,2] [,3] [,4]
#> [1,] 1.0258 0.0004 0.0135 -0.0586
#> [2,] 0.0000 1.0786 0.0152 -0.1092
#> [3,] 0.0000 0.0000 0.9738 -0.0830
#> [4,] 0.0000 0.0000 0.0000 0.9674#> [1] TRUE# Simulating correlated Gaussians is one Cholesky away
Sigma <- matrix(c(4, 2, 2, 3), 2)
L <- t(chol(Sigma))
set.seed(41)
sim <- t(L %*% matrix(rnorm(2 * 10000), 2))
round(cov(sim), 3); Sigma#> [,1] [,2]
#> [1,] 3.879 1.883
#> [2,] 1.883 2.887#> [,1] [,2]
#> [1,] 4 2
#> [2,] 2 38.3 QR decomposition
QR factorization (or QR decomposition) is a fundamental mathematical technique in linear algebra that breaks down any matrix \(A\) into two specific matrices:
\[A = QR,\] where
- \(Q\) (Orthogonal Matrix) is a matrix whose columns are orthonormal (meaning they are perpendicular to each other and each have a length of \(1\)). Geometrically, multiplying by \(Q\) represents a rotation or reflection, which preserves lengths and angles.
- \(R\) (Upper Triangular Matrix) is a matrix where all the entries below the main diagonal are zero. Multiplying by \(R\) represents scaling and shearing.
While LU factorization works for square matrices, QR factorization is especially powerful for non-square, rectangular matrices, \(n \ge p\), where there are more rows than columns.
Solving Least Squares Problems (Linear Regression): To find the line of best fit for a dataset with more data points (\(n\)) than variables (\(p\)), the standard normal equation approach (\(X^\top X b = X^\top y\)) can be numerically unstable. Using QR factorization (\(X = QR\)) makes solving the least squares problem numerically stable and robust.
Unlike simple methods, like standard Gram-Schmidt orthogonalization, proper QR algorithms, using Householder reflections or Givens rotations, avoid catastrophic round-off computational errors, ensuring the columns of \(Q\) remain strictly orthogonal. The iterative QR algorithm is one of the most famous methods used by software, e.g., NumPy, Julia, to find all the eigenvalues of a matrix simultaneously.
Any \(X\in\mathbb{R}^{n\times p}\)
with \(n\ge p\) factors as \(X=QR\) with \(Q\in\mathbb{R}^{n\times p}\) having
orthonormal columns (\(Q^\top Q=I_p\))
and \(R\in\mathbb{R}^{p\times p}\)
upper-triangular. R uses Householder
reflections, which cost \(2np^2-\tfrac{2}{3}p^3\) and are backward
stable. (Classical Gram–Schmidt produces the same factors in exact
arithmetic but loses orthogonality catastrophically in floating point;
modified Gram–Schmidt is better but Householder is better still.)
QR is the workhorse of least squares. Substituting \(X=QR\):
\[\|\mathbf{y}-X\boldsymbol\beta\|_2^2=\|\mathbf{y}-QR\boldsymbol\beta\|_2^2 =\|Q^\top\mathbf{y}-R\boldsymbol\beta\|_2^2+\|(I-QQ^\top)\mathbf{y}\|_2^2 ,\]
because multiplying by an orthogonal matrix preserves \(\ell_2\) length. The second term does not involve \(\boldsymbol\beta\), so the minimizer solves the triangular system
\[R\hat{\boldsymbol\beta}=Q^\top\mathbf{y},\]
by back-substitution in \(O(p^2)\). Crucially, \(X^\top X\) is never formed.
set.seed(43)
Xq <- cbind(1, matrix(rnorm(100 * 3), 100))
yq <- Xq %*% c(2, -1, 0.5, 3) + rnorm(100)
qrX <- qr(Xq)
Q <- qr.Q(qrX); Rm <- qr.R(qrX)
c(orthonormal = all.equal(crossprod(Q), diag(ncol(Xq))),
reconstructs = all.equal(Q %*% Rm, Xq, check.attributes = FALSE),
R_upper_triangular = all(abs(Rm[lower.tri(Rm)]) < 1e-12))#> orthonormal reconstructs R_upper_triangular
#> TRUE TRUE TRUE#> [1] 1.921448 -0.949995 0.585666 3.017408#> [1] 1.921448 -0.949995 0.585666 3.0174088.4 Eigen-spectra, Eigen-decomposition
A nonzero \(\mathbf{v}\) with \(A\mathbf{v}=\lambda\mathbf{v}\) is an eigenvector; \(\lambda\) is the corresponding eigenvalue. The eigenvectors are the directions the transformation leaves invariant up to scaling, the axes along which \(A\) merely stretches, compresses, or flips.
The historical thread runs from Euler’s work on rotational motion through Lagrange’s inertia matrices (principal axes) to Hilbert, who fixed the prefix eigen, meaning “own” or “characteristic”.
Collecting eigenpairs into \(V=[\mathbf{v}_1\cdots\mathbf{v}_n]\) and \(\Lambda=\operatorname{diag}(\lambda_1,\dots,\lambda_n)\) gives \(AV=V\Lambda\); if \(V\) is invertible, \(A=V\Lambda V^{-1}\).
Spectral theorem. If \(A\) is real symmetric, its eigenvalues are real, its eigenvectors can be chosen orthonormal, and \(A=Q\Lambda Q^\top\) with \(Q^\top Q=I\). Moreover \(A\) is positive definite iff every \(\lambda_j>0\).
#> [1] 5 2#> [,1] [,2]
#> [1,] 0.7071 -0.4472
#> [2,] 0.7071 0.8944# The correct verification is A V = V Lambda (equivalently A v_j = lambda_j v_j)
V <- e$vectors; L <- diag(e$values)
round(Ae %*% V - V %*% L, 12)#> [,1] [,2]
#> [1,] 0 0
#> [2,] 0 0#> [,1] [,2]
#> [1,] 4 1
#> [2,] 2 3Common misconception: “check eigenvectors with
(values * diag(n) - A) %*% vectors.” That expression evaluates \((\Lambda-A)V\), which is not zero in general, the eigen-relation is \(AV=V\Lambda\), and \(\Lambda\) must multiply from the right. The two agree only when \(A\) is a scalar multiple of the identity, where every vector is an eigenvector and the check is vacuous. Verify one pair at a time, or use \(AV-V\Lambda\).
For a symmetric matrix the picture is cleaner and geometrically meaningful:
Sy <- matrix(c(4, 1.5, 1.5, 2), 2)
es <- eigen(Sy, symmetric = TRUE)
c(eigenvalues = es$values,
orthonormal = all.equal(crossprod(es$vectors), diag(2)),
reconstructs = all.equal(es$vectors %*% diag(es$values) %*% t(es$vectors), Sy),
positive_definite = all(es$values > 0))#> eigenvalues1 eigenvalues2 orthonormal reconstructs
#> 4.80278 1.19722 1.00000 1.00000
#> positive_definite
#> 1.00000A sample covariance matrix is symmetric positive semi-definite, so its eigendecomposition has a direct interpretation: the eigenvectors are the principal axes of the data ellipsoid and the eigenvalues are the variances along those axes. That is the whole content of principal component analysis (Chapter 4), and it is best seen in three dimensions, rotating:
mlb <- dspa_read(
"https://umich.instructure.com/files/330381/download?download_frd=1",
"01a_data.txt", reader = read.table, header = TRUE)
M <- as.matrix(mlb[, c("Height", "Weight", "Age")])
ctr <- colMeans(M); Sg <- cov(M)
eg <- eigen(Sg, symmetric = TRUE)
# 1-sd ellipsoid {x : (x-mu)' S^{-1} (x-mu) = 1} = mu + Q L^{1/2} u, |u| = 1
th <- seq(0, 2 * pi, length.out = 60)
ph <- seq(0, pi, length.out = 30)
U <- rbind(as.vector(outer(cos(th), sin(ph))),
as.vector(outer(sin(th), sin(ph))),
as.vector(outer(rep(1, 60), cos(ph))))
E <- eg$vectors %*% diag(sqrt(eg$values)) %*% U + ctr
ex <- matrix(E[1, ], 60); ey <- matrix(E[2, ], 60); ez <- matrix(E[3, ], 60)
p3 <- plot_ly() |>
add_surface(x = ex, y = ey, z = ez, opacity = 0.30, showscale = FALSE,
colorscale = list(c(0, "#9EC5E8"), c(1, "#9EC5E8"))) |>
add_trace(x = M[, 1], y = M[, 2], z = M[, 3], type = "scatter3d",
mode = "markers", name = "MLB players",
marker = list(size = 2, opacity = 0.35, color = "grey30"))
for (j in 1:3) {
ax <- ctr + eg$vectors[, j] * sqrt(eg$values[j])
p3 <- add_trace(p3, x = c(ctr[1], ax[1]), y = c(ctr[2], ax[2]),
z = c(ctr[3], ax[3]), type = "scatter3d", mode = "lines",
name = sprintf("PC%d (var %.1f)", j, eg$values[j]),
line = list(width = 8))
}
p3 |> layout(title = "Covariance ellipsoid and its principal axes (eigenvectors)",
scene = list(xaxis = list(title = "Height (in)"),
yaxis = list(title = "Weight (lb)"),
zaxis = list(title = "Age (yr)")))data.frame(component = paste0("PC", 1:3),
eigenvalue = round(eg$values, 3),
prop_variance = round(eg$values / sum(eg$values), 4),
cumulative = round(cumsum(eg$values) / sum(eg$values), 4))8.5 Singular value decomposition (SVD)
The SVD is the most general and most useful decomposition. Every matrix \(A\in\mathbb{R}^{m\times n}\) factors as
\[A=U\Sigma V^\top,\qquad U^\top U=I,\quad V^\top V=I,\quad \Sigma=\operatorname{diag}(\sigma_1\ge\sigma_2\ge\cdots\ge\sigma_r>0).\]
No squareness, no invertibility, no symmetry required. Equivalently,
\[A=\sum_{j=1}^{r}\sigma_j\,\mathbf{u}_j\mathbf{v}_j^\top\]
a sum of \(r\) rank-one outer products (§3.2.3), each weighted by its singular value.
Relationships worth knowing: \(\sigma_j=\sqrt{\lambda_j(A^\top A)}\); the columns of \(V\) are eigenvectors of \(A^\top A\) and of \(U\) eigenvectors of \(AA^\top\); \(\operatorname{rank}(A)=\#\{\sigma_j>0\}\); \(\|A\|_2=\sigma_1\); and the first \(r\) columns of \(U\) span \(\operatorname{col}(A)\) while the last \(n-r\) columns of \(V\) span \(\mathcal{N}(A)\), the SVD hands you all four fundamental subspaces at once.
Eckart–Young–Mirsky theorem. Among all matrices of rank at most \(k\), the truncation \(A_k=\sum_{j\le k}\sigma_j\mathbf{u}_j\mathbf{v}_j^\top\) minimizes both \(\|A-A_k\|_2=\sigma_{k+1}\) and \(\|A-A_k\|_F=\big(\sum_{j>k}\sigma_j^2\big)^{1/2}\).
This is the theoretical basis of PCA, latent semantic analysis, image compression, recommender systems, and denoising. Cost: \(\approx 14mn^2\) for the thin SVD when \(m\ge n\).
set.seed(47)
Am <- matrix(rnorm(6 * 4), 6)
s <- svd(Am)
c(reconstructs = all.equal(s$u %*% diag(s$d) %*% t(s$v), Am),
U_orthonormal = all.equal(crossprod(s$u), diag(4)),
V_orthonormal = all.equal(crossprod(s$v), diag(4)),
spectral_norm = all.equal(max(s$d), norm(Am, "2")),
sigma_vs_eigen = all.equal(s$d^2, eigen(crossprod(Am))$values))#> reconstructs U_orthonormal V_orthonormal spectral_norm sigma_vs_eigen
#> TRUE TRUE TRUE TRUE TRUELow-rank approximation, demonstrated on real high-dimensional data. The GSE5859 microarray study measured expression for 8,793 genes across 208 lymphoblastoid samples from three HapMap populations (CEU, CHB, JPT).
gene <- dspa_read("https://umich.instructure.com/files/2001417/download?download_frd=1",
"exprs_GSE5859.csv", header = TRUE)
info <- dspa_read("https://umich.instructure.com/files/2001418/download?download_frd=1",
"exprs_MetaData_GSE5859.csv", header = TRUE)
dim(gene); dim(info)#> [1] 8793 209#> [1] 208 4G <- as.matrix(gene[1:200, 2:101]) # 200 genes x 100 samples
G <- sweep(G, 1, rowMeans(G)) # centre each gene
sg <- svd(G)
energy <- cumsum(sg$d^2) / sum(sg$d^2)
ranks <- c(1, 2, 5, 10, 25, 50)
approx_err <- sapply(ranks, function(k) {
Gk <- sg$u[, 1:k, drop = FALSE] %*% diag(sg$d[1:k], k, k) %*% t(sg$v[, 1:k, drop = FALSE])
c(frobenius_rel = norm(G - Gk, "F") / norm(G, "F"),
predicted = sqrt(sum(sg$d[-(1:k)]^2)) / sqrt(sum(sg$d^2)))
})
data.frame(rank = ranks, t(round(approx_err, 5)),
variance_captured = round(energy[ranks], 4))Observed and Eckart–Young-predicted errors agree to five decimals, the theorem is not an approximation, it is exact.
scree <- data.frame(k = 1:40, sigma = sg$d[1:40], cum = energy[1:40])
p_a <- ggplot(scree, aes(k, sigma)) +
geom_col(fill = "steelblue") +
labs(title = "Singular value spectrum", x = "index j",
y = expression(sigma[j])) + theme_dspa(10)
p_b <- ggplot(scree, aes(k, cum)) +
geom_line(linewidth = 0.9, colour = "firebrick") + geom_point(size = 1.3) +
geom_hline(yintercept = 0.9, linetype = "dashed", colour = "grey40") +
scale_y_continuous(labels = scales::percent) +
labs(title = "Cumulative variance captured", x = "rank k", y = NULL) +
theme_dspa(10)
p_a | p_b
# --- Interactive equivalent ------------------------------------------------
plot_ly(x = 1:40, y = sg$d[1:40], type = "bar", name = "singular values") |>
add_lines(x = 1:40, y = energy[1:40] * max(sg$d), yaxis = "y2",
name = "cumulative variance") |>
layout(title = "Singular value spectrum (GSE5859 submatrix)",
xaxis = list(title = "index j"),
yaxis = list(title = "sigma_j"),
yaxis2 = list(overlaying = "y", side = "right", tickformat = ".0%",
range = c(0, 1), title = "cumulative"),
legend = list(orientation = "h"))The reconstruction quality is a surface over (gene, sample, expression), so this one is interactive:
mk <- function(k) {
sg$u[, 1:k, drop = FALSE] %*% diag(sg$d[1:k], k, k) %*% t(sg$v[, 1:k, drop = FALSE])
}
plot_ly() |>
add_surface(z = G, showscale = FALSE, opacity = 0.95,
colorscale = "Viridis", name = "full rank") |>
add_surface(z = mk(3) + 6, showscale = FALSE, opacity = 0.95,
colorscale = "Viridis", name = "rank 3 (offset +6)") |>
layout(title = "Gene-expression matrix: full rank (bottom) vs. rank-3 approximation (top)",
scene = list(xaxis = list(title = "Sample"),
yaxis = list(title = "Gene"),
zaxis = list(title = "Centred expression")))Three of one hundred components reproduce the dominant structure. The residual is not noise in any trivial sense, it is where the gene-specific signal lives, but the rank-3 surface shows how much of the matrix is driven by a handful of sample-level factors (batch, population, processing date).
9 Conditioning and numerical stability
Floating-point arithmetic is exact only for a measure-zero set of numbers. Whether that matters depends on the condition number of your problem.
For an invertible \(A\), the 2-norm condition number is
\[\kappa_2(A)=\|A\|_2\,\|A^{-1}\|_2=\frac{\sigma_{\max}(A)}{\sigma_{\min}(A)}\ \ge 1 .\]
It bounds how much a relative perturbation of the input can be amplified in the output: if \(A\mathbf{x}=\mathbf{b}\) and \(A(\mathbf{x}+\delta\mathbf{x})=\mathbf{b}+\delta\mathbf{b}\), then
\[\frac{\|\delta\mathbf{x}\|}{\|\mathbf{x}\|}\le\kappa_2(A)\,\frac{\|\delta\mathbf{b}\|}{\|\mathbf{b}\|} .\]
Rule of thumb. In IEEE double precision the machine epsilon is \(\epsilon\approx 2.2\times10^{-16}\), i.e. about 16 significant decimal digits. A problem with \(\kappa\approx 10^{t}\) loses roughly \(t\) of them.
# Hilbert matrix is a classic pathological example
hilbert <- function(n) outer(1:n, 1:n, function(i, j) 1 / (i + j - 1))
data.frame(n = 3:9,
kappa = sapply(3:9, \(n) kappa(hilbert(n), exact = TRUE)),
digits_lost = round(log10(sapply(3:9, \(n) kappa(hilbert(n), exact = TRUE))), 1))Hilbert matrix Example: The Hilbert matrix is notoriously ill-conditioned because its columns are nearly collinear as \(n\) grows, making it nearly rank-deficient in finite precision.
The Hilbert matrix is the classic pathological example, by \(n=9\) its condition number exceeds \(10^{11}\), so a linear solve retains only about five correct digits.
n <- 9
Hn <- hilbert(n)
x_true <- rep(1, n)
b_h <- Hn %*% x_true
x_hat <- solve(Hn, b_h)
c(max_abs_error = max(abs(x_hat - x_true)),
kappa = kappa(Hn, exact = TRUE),
residual_norm = norm(Hn %*% x_hat - b_h, "2"))#> max_abs_error kappa residual_norm
#> 1.98434e-05 4.93154e+11 3.68219e-16Note the problem: the residual is tiny while the error is large. A small residual proves the algorithm was backward stable. It says nothing about forward accuracy when the problem is ill-conditioned. Reporting only a residual is how bad numerics goes undetected.
9.1 The result that governs least squares
Proposition. For \(X\) of full column rank, \[\kappa_2(X^\top X)=\kappa_2(X)^2 .\]
Proof. The singular values of \(X^\top X\) are the squares of those of \(X\) (§3.6.5), so \(\kappa_2(X^\top X)=\sigma_{\max}^2/\sigma_{\min}^2=\kappa_2(X)^2\). \(\blacksquare\)
Consequence. Solving least squares by forming and inverting \(X^\top X\), the “normal equations” route, doubles the number of digits you lose. If \(\kappa_2(X)=10^{6}\), which happens routinely with polynomial terms or near-duplicate predictors, then \(\kappa_2(X^\top X)=10^{12}\) and only about four significant digits survive. QR works with \(X\) directly and loses only \(\log_{10}\kappa_2(X)\) digits.
This is the single most important computational fact in the chapter, and §3.11 measures it.
10 PART II: LINEAR REGRESSION
11 The linear model
11.1 Notation
Let \(Y_1,\dots,Y_n\) be an observed response and \(x_{i1},\dots,x_{ip-1}\) the covariates for case \(i\). The linear model is
\[Y_i=\beta_0+\beta_1x_{i1}+\cdots+\beta_{p-1}x_{i,p-1}+\varepsilon_i, \qquad i=1,\dots,n,\]
or, stacking cases,
\[\underbrace{\mathbf{Y}}_{n\times1}=\underbrace{X}_{n\times p}\underbrace{\boldsymbol\beta}_{p\times1}+\underbrace{\boldsymbol\varepsilon}_{n\times1}, \qquad X=\begin{pmatrix}1&x_{11}&\cdots&x_{1,p-1}\\ 1&x_{21}&\cdots&x_{2,p-1}\\ \vdots&\vdots&\ddots&\vdots\\ 1&x_{n1}&\cdots&x_{n,p-1}\end{pmatrix}.\]
The leading column of ones supplies the intercept. R
represents vectors as \(n\times 1\)
columns by default.
“Linear” refers to linearity in \(\boldsymbol\beta\), not in \(x\). The model \(Y=\beta_0+\beta_1x+\beta_2x^2+\varepsilon\) is a linear model, the design matrix simply has a column \(x^2\). So is \(Y=\beta_0+\beta_1\log x+\beta_2\sin x+\varepsilon\). What is not linear is \(Y=\beta_0 e^{\beta_1 x}+\varepsilon\), where the parameter sits inside a nonlinear function. This distinction is what lets a “linear” model capture curvature (§3.16).
11.2 Sample statistics as matrix operations
Let \(\mathbf{1}\) be the \(n\times1\) vector of ones. Then
\[\bar Y=\frac{1}{n}\mathbf{1}^\top\mathbf{Y}, \qquad \widehat{\operatorname{Var}}(Y)=\frac{1}{n-1}\big(\mathbf{Y}-\bar Y\mathbf{1}\big)^\top\big(\mathbf{Y}-\bar Y\mathbf{1}\big).\]
y <- mlb$Height
n <- length(y)
Y <- matrix(y, n, 1)
one <- matrix(1, n, 1)
c(matrix_mean = as.numeric(crossprod(one, Y) / n),
base_mean = mean(y))#> matrix_mean base_mean
#> 73.6973 73.6973#> matrix_var base_var
#> 5.3168 5.3168Centring is itself a projection: \(\mathbf{Y}-\bar Y\mathbf{1}=(I-P_{\mathbf 1})\mathbf{Y}\) where \(P_{\mathbf 1}=\tfrac{1}{n}\mathbf{1}\mathbf{1}^\top\) is the hat matrix of an intercept-only model. So “subtract the mean” and “regress on a constant and keep the residual” are literally the same operation.
Weighted averages generalize the same way. For \(A=(a_1,\dots,a_n)\),
\[A X=\left(\sum_i a_i x_{i1},\ \dots,\ \sum_i a_i x_{ip}\right),\]
which lets a single matrix product compute any set of linear contrasts at once.
G1 <- as.matrix(gene[, -1])
colmeans_apply <- colMeans(G1)
colmeans_matrix <- as.vector(crossprod(rep(1 / nrow(G1), nrow(G1)), G1))
c(agree = all.equal(unname(colmeans_apply), colmeans_matrix),
first_five = round(colmeans_matrix[1:5], 4))#> agree first_five1 first_five2 first_five3 first_five4 first_five5
#> 1.0000 5.7040 5.7218 5.7263 5.7436 5.8355ggplot(data.frame(m = colmeans_matrix), aes(m)) +
geom_histogram(bins = nclass.FD(colmeans_matrix),
fill = "steelblue", colour = "white") +
labs(title = "Average gene expression per sample (GSE5859)",
subtitle = "Roughly symmetric and unimodal, as array normalization intends",
x = "Sample mean expression", y = "Count") +
theme_dspa()
# --- Interactive equivalent ------------------------------------------------
h <- hist(colmeans_matrix, plot = FALSE)
plot_ly(x = h$mids, y = h$counts, type = "bar", name = "Column Averages") |>
layout(title = "Average Gene Expression Histogram",
xaxis = list(title = "Column means"),
yaxis = list(title = "Frequency"),
legend = list(orientation = "h"))A single matrix product can compute a mean and a group contrast simultaneously. Choosing weights \(a_i=1/n_M\) for males and \(a_i=-1/n_F\) for females makes the second output column the male-minus-female difference for each gene:
\[X\begin{pmatrix}1/p & a_1\\ 1/p & a_2\\ \vdots&\vdots\\ 1/p & a_p\end{pmatrix} =\begin{pmatrix}\bar X_1 & \text{diff}_1\\ \vdots&\vdots\\ \bar X_N & \text{diff}_N\end{pmatrix}.\]
sex_map <- info[, c("filename", "sex")]
rownames(sex_map) <- sex_map$filename
sex_map <- sex_map[colnames(G1), ] # align to the expression columns
stopifnot(identical(rownames(sex_map), colnames(G1)))
nF <- sum(sex_map$sex == "F"); nM <- sum(sex_map$sex == "M")
c(females = nF, males = nM)#> females males
#> 86 122W <- cbind(mean = rep(1 / ncol(G1), ncol(G1)),
male_minus_female = ifelse(sex_map$sex == "F", -1 / nF, 1 / nM))
contrasts_mat <- G1 %*% W
round(head(contrasts_mat, 6), 4)#> mean male_minus_female
#> [1,] 6.3833 -0.0032
#> [2,] 7.0916 -0.0313
#> [3,] 5.4770 0.0648
#> [4,] 7.5840 -0.0013
#> [5,] 3.1977 0.0153
#> [6,] 7.3382 0.078412 Ordinary least squares
12.1 Three derivations of the same estimator
Define the residual sum of squares
\[f(\boldsymbol\beta)=\|\mathbf{Y}-X\boldsymbol\beta\|_2^2=(\mathbf{Y}-X\boldsymbol\beta)^\top(\mathbf{Y}-X\boldsymbol\beta).\]
Derivation 1: calculus. Expanding, \(f(\boldsymbol\beta)=\mathbf{Y}^\top\mathbf{Y}-2\boldsymbol\beta^\top X^\top\mathbf{Y}+\boldsymbol\beta^\top X^\top X\boldsymbol\beta\). Using \(\nabla_{\boldsymbol\beta}(\mathbf{a}^\top\boldsymbol\beta)=\mathbf{a}\) and \(\nabla_{\boldsymbol\beta}(\boldsymbol\beta^\top A\boldsymbol\beta)=2A\boldsymbol\beta\) for symmetric \(A\),
\[\nabla f(\boldsymbol\beta)=-2X^\top\mathbf{Y}+2X^\top X\boldsymbol\beta \;\stackrel{!}{=}\;\mathbf{0} \quad\Longrightarrow\quad \underbrace{X^\top X\hat{\boldsymbol\beta}=X^\top\mathbf{Y}}_{\text{normal equations}} .\]
The Hessian is \(\nabla^2 f=2X^\top X\), positive definite whenever \(\operatorname{rank}(X)=p\), so this critical point is the unique global minimum, \(f\) is a convex quadratic. Hence
\[\boxed{\;\hat{\boldsymbol\beta}=(X^\top X)^{-1}X^\top\mathbf{Y}\;}\]
Derivation 2: geometry. \(X\boldsymbol\beta\) ranges over \(\operatorname{col}(X)\). The closest point of a subspace to \(\mathbf{Y}\) is its orthogonal projection (§3.4.2), characterized by residual orthogonality \(X^\top(\mathbf{Y}-X\hat{\boldsymbol\beta})=0\), the normal equations again, obtained with no calculus at all.
Derivation 3: maximum likelihood. If \(\boldsymbol\varepsilon\sim N_n(\mathbf{0},\sigma^2I)\), the log-likelihood is
\[\ell(\boldsymbol\beta,\sigma^2)=-\frac{n}{2}\log(2\pi\sigma^2)-\frac{1}{2\sigma^2}\|\mathbf{Y}-X\boldsymbol\beta\|_2^2 ,\]
so maximizing over \(\boldsymbol\beta\) is minimizing the RSS. Under Gaussian errors, OLS is the MLE.
12.2 Why squared loss?
Three reasons, none of which is “the minima of \(f\) and \(f^2\) coincide”, a claim that is false in general, since \(\tfrac{d}{dx}f^2=2ff'\) vanishes wherever \(f=0\) as well as wherever \(f'=0\).
- Smoothness and a closed form. \(\sum_i\varepsilon_i^2\) is differentiable everywhere and quadratic, giving a unique analytic solution. \(\sum_i|\varepsilon_i|\) is not differentiable at zero and requires linear programming.
- Distributional match. Squared loss is the Gaussian log-likelihood (Derivation 3). Absolute loss corresponds to a Laplace error model.
- Optimality. Gauss–Markov (§3.9.3) says OLS is best-in-class under second-moment assumptions alone.
The tradeoff is robustness: squaring gives an outlier at \(5\sigma\) twenty-five times the influence of one at \(1\sigma\). Median regression (\(L_1\)) has a breakdown point above zero; OLS does not.
set.seed(53)
xs <- 1:40
ys <- 2 + 0.8 * xs + rnorm(40, sd = 1.5)
ys[38] <- ys[38] + 28 # one contaminated observation
ols <- lm(ys ~ xs)
lad <- quantreg::rq(ys ~ xs, tau = 0.5)
ggplot(data.frame(xs, ys), aes(xs, ys)) +
geom_point(size = 1.8, colour = "grey30") +
geom_point(data = data.frame(xs = 38, ys = ys[38]), colour = "firebrick", size = 3.5) +
geom_abline(aes(intercept = coef(ols)[1], slope = coef(ols)[2],
colour = "OLS (L2)"), linewidth = 1) +
geom_abline(aes(intercept = coef(lad)[1], slope = coef(lad)[2],
colour = "Median regression (L1)"), linewidth = 1) +
scale_colour_manual(values = c("OLS (L2)" = "steelblue",
"Median regression (L1)" = "darkgreen")) +
labs(title = "One outlier is enough to tilt the least-squares line",
subtitle = sprintf("OLS slope %.3f | L1 slope %.3f | true slope 0.800",
coef(ols)[2], coef(lad)[2]),
x = "x", y = "y", colour = NULL) +
theme_dspa()
12.3 The loss surface
\(f(\boldsymbol\beta)\) is a convex quadratic bowl in parameter space, and OLS finds its bottom. With two parameters this is literally a surface, so it is worth rotating:
ha <- dspa_read("https://umich.instructure.com/files/1644953/download?download_frd=1",
"CaseStudy12_AdultsHeartAttack_Data.csv", stringsAsFactors = FALSE)
ha$CHARGES <- suppressWarnings(as.numeric(ha$CHARGES))
ha <- ha[complete.cases(ha), ]
Xh <- cbind(1, ha$LOS); yh <- ha$CHARGES
fit_h <- lm(CHARGES ~ LOS, data = ha)
bh <- unname(coef(fit_h))
g0 <- seq(bh[1] - 3000, bh[1] + 3000, length.out = 70)
g1 <- seq(bh[2] - 300, bh[2] + 300, length.out = 70)
RSS <- outer(g0, g1, Vectorize(function(a, b) sum((yh - a - b * ha$LOS)^2)))
plot_ly() |>
add_surface(x = g1, y = g0, z = RSS, opacity = 0.9, showscale = FALSE,
colorscale = "Viridis") |>
add_trace(x = bh[2], y = bh[1], z = sum(residuals(fit_h)^2),
type = "scatter3d", mode = "markers", name = "OLS minimum",
marker = list(size = 7, color = "red")) |>
layout(title = "Residual sum of squares over parameter space",
scene = list(xaxis = list(title = "slope beta1"),
yaxis = list(title = "intercept beta0"),
zaxis = list(title = "RSS"),
camera = list(eye = list(x = 1.6, y = -1.6, z = 0.9))))The surface is an elliptic paraboloid with a single minimum. The elongation and tilt of its contours are governed by \(X^\top X\): a nearly singular \(X^\top X\) gives a long narrow valley, which is exactly what multicollinearity looks like geometrically (§3.13) and why the estimates become unstable.
12.4 The Gauss–Markov theorem
Let’s state the assumptions precisely, because which conclusions survive depends on which assumption fails.
| Assumption | Needed for | |
|---|---|---|
| A1 | \(\mathbf{Y}=X\boldsymbol\beta+\boldsymbol\varepsilon\) (linear in \(\boldsymbol\beta\)) | Unbiasedness |
| A2 | \(E[\boldsymbol\varepsilon\mid X]=\mathbf{0}\) | Unbiasedness |
| A3 | \(\operatorname{Var}(\boldsymbol\varepsilon\mid X)=\sigma^2I_n\) (homoscedastic, uncorrelated) | Efficiency, valid SEs |
| A4 | \(\operatorname{rank}(X)=p\) | Existence and uniqueness |
| A5 | \(\boldsymbol\varepsilon\mid X\sim N(\mathbf{0},\sigma^2I)\) | Exact \(t\) and \(F\) distributions |
Gauss–Markov theorem. Under A1–A4, \(\hat{\boldsymbol\beta}_{\text{OLS}}\) is the Best Linear Unbiased Estimator: for any other estimator \(\tilde{\boldsymbol\beta}=C\mathbf{Y}\) that is linear in \(\mathbf{Y}\) and unbiased, \(\operatorname{Var}(\tilde{\boldsymbol\beta})-\operatorname{Var}(\hat{\boldsymbol\beta})\) is positive semi-definite. Normality (A5) is not required.
Proof sketch. Write \(C=(X^\top X)^{-1}X^\top+D\). Unbiasedness for all \(\boldsymbol\beta\) forces \(DX=0\). Then \(\operatorname{Var}(C\mathbf{Y})=\sigma^2\big[(X^\top X)^{-1}+DD^\top\big]\), and \(DD^\top\succeq 0\). \(\blacksquare\)
Read the qualifiers carefully, each one is a doorway to a later chapter. “Best” means minimum variance among linear unbiased estimators. Ridge and LASSO are biased and can have strictly smaller mean squared error (Chapter 11). Under heavy-tailed errors, nonlinear estimators beat OLS. And when A3 fails, generalized least squares does better.
13 Inference
13.1 The sampling distribution of \(\hat{\boldsymbol\beta}\)
Substituting \(\mathbf{Y}=X\boldsymbol\beta+\boldsymbol\varepsilon\):
\[\hat{\boldsymbol\beta}=(X^\top X)^{-1}X^\top(X\boldsymbol\beta+\boldsymbol\varepsilon) =\boldsymbol\beta+(X^\top X)^{-1}X^\top\boldsymbol\varepsilon .\]
Hence, using \(\operatorname{Var}(A\boldsymbol\varepsilon)=A\operatorname{Var}(\boldsymbol\varepsilon)A^\top\),
\[E[\hat{\boldsymbol\beta}]=\boldsymbol\beta,\qquad \boxed{\;\operatorname{Var}(\hat{\boldsymbol\beta})=\sigma^2(X^\top X)^{-1}\;}\]
Every standard error in every regression table is \(\hat\sigma\sqrt{[(X^\top X)^{-1}]_{jj}}\).
The error variance is estimated by
\[\hat\sigma^2=\frac{\mathrm{RSS}}{n-p}=\frac{\|(I-H)\mathbf{Y}\|_2^2}{n-p},\]
and the divisor is \(n-p\) rather than \(n\) precisely because \(\operatorname{tr}(I-H)=n-p\) (§3.4.2), the residual vector lives in an \((n-p)\)-dimensional subspace, so it has that many degrees of freedom. Then \(E[\hat\sigma^2]=\sigma^2\).
Under A5, \(\hat{\boldsymbol\beta}\sim N_p\big(\boldsymbol\beta,\sigma^2(X^\top X)^{-1}\big)\), \((n-p)\hat\sigma^2/\sigma^2\sim\chi^2_{n-p}\), the two are independent, and
\[t_j=\frac{\hat\beta_j-\beta_j^{0}}{\hat\sigma\sqrt{[(X^\top X)^{-1}]_{jj}}}\sim t_{n-p}.\]
Xm <- model.matrix(~ LOS, data = ha)
XtXi <- solve(crossprod(Xm))
beta_hat <- XtXi %*% crossprod(Xm, ha$CHARGES)
res <- ha$CHARGES - Xm %*% beta_hat
n_h <- nrow(Xm); p_h <- ncol(Xm)
sigma2 <- sum(res^2) / (n_h - p_h)
se <- sqrt(sigma2 * diag(XtXi))
tval <- as.vector(beta_hat) / se
pval <- 2 * pt(abs(tval), df = n_h - p_h, lower.tail = FALSE)
manual <- data.frame(Estimate = as.vector(beta_hat), `Std. Error` = se,
`t value` = tval, `Pr(>|t|)` = pval, check.names = FALSE)
round(manual, 5)#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) 4582.700 399.6377 11.46714 0.00000
#> LOS 212.287 69.5324 3.05306 0.0026913.2 Decomposing variance
The Pythagorean identity from §3.4.2, applied after centring:
\[\underbrace{\|\mathbf{Y}-\bar Y\mathbf{1}\|^2}_{\mathrm{SST}} =\underbrace{\|\hat{\mathbf{Y}}-\bar Y\mathbf{1}\|^2}_{\mathrm{SSR}} +\underbrace{\|\mathbf{Y}-\hat{\mathbf{Y}}\|^2}_{\mathrm{SSE}} .\]
\[R^2=\frac{\mathrm{SSR}}{\mathrm{SST}}=1-\frac{\mathrm{SSE}}{\mathrm{SST}}, \qquad R^2_{\text{adj}}=1-\frac{\mathrm{SSE}/(n-p)}{\mathrm{SST}/(n-1)} .\]
Common misconception: “a higher \(R^2\) means a better model.” \(R^2\) is non-decreasing in the number of predictors, adding a column of pure noise can never lower it, because the larger column space contains the smaller one, so the projection can only get closer to \(\mathbf{Y}\). Add \(n-1\) random predictors and \(R^2\) reaches 1 exactly. Adjusted \(R^2\) penalizes \(p\) and can decrease. Neither is a threshold: what counts as a good \(R^2\) depends entirely on the field, a value of 0.02 is meaningful in asset pricing and a value of 0.95 is suspicious in survey research.
set.seed(59)
n_i <- 60
y_i <- rnorm(n_i)
Xnoise <- matrix(rnorm(n_i * 50), n_i) # 50 columns of PURE NOISE
r2 <- sapply(1:50, function(k) {
m <- lm(y_i ~ Xnoise[, 1:k])
c(R2 = summary(m)$r.squared, adjR2 = summary(m)$adj.r.squared)
})
data.frame(k = 1:50, t(r2)) |>
pivot_longer(-k, names_to = "measure", values_to = "value") |>
ggplot(aes(k, value, colour = measure)) +
geom_line(linewidth = 1) +
geom_hline(yintercept = 0, linetype = "dashed", colour = "grey50") +
scale_colour_manual(values = c(R2 = "firebrick", adjR2 = "steelblue"),
labels = c(R2 = expression(R^2),
adjR2 = expression(Adjusted~R^2))) +
labs(title = "Fifty predictors of pure noise, one response of pure noise",
subtitle = "R-squared climbs toward 1 regardless; adjusted R-squared does not",
x = "Number of noise predictors", y = NULL, colour = NULL) +
theme_dspa()
The \(F\)-test for the whole model compares the fitted model against the intercept-only model:
\[F=\frac{\mathrm{SSR}/(p-1)}{\mathrm{SSE}/(n-p)}\sim F_{p-1,\,n-p}\quad\text{under }H_0:\beta_1=\cdots=\beta_{p-1}=0 .\]
SST <- sum((ha$CHARGES - mean(ha$CHARGES))^2)
SSE <- sum(residuals(fit_h)^2)
SSR <- SST - SSE
Fstat <- (SSR / (p_h - 1)) / (SSE / (n_h - p_h))
c(R2 = SSR / SST,
R2_from_lm = summary(fit_h)$r.squared,
F = Fstat,
F_from_lm = summary(fit_h)$fstatistic[["value"]],
p_value = pf(Fstat, p_h - 1, n_h - p_h, lower.tail = FALSE))#> R2 R2_from_lm F F_from_lm p_value
#> 0.06001243 0.06001243 9.32120238 9.32120238 0.0026924114 Computing OLS: three routes, one answer
The formula \(\hat{\boldsymbol\beta}=(X^\top X)^{-1}X^\top\mathbf{Y}\) tells you what to compute, not how. There are three standard routes with very different numerical behaviour.
| Route | What it does | Flops | Error scales as |
|---|---|---|---|
| Normal equations | form \(X^\top X\), Cholesky-solve | \(np^2+\tfrac{1}{3}p^3\) | \(\kappa_2(X)^2\) |
| QR (Householder) | \(X=QR\), back-solve \(R\hat\beta=Q^\top\mathbf{Y}\) | \(2np^2-\tfrac{2}{3}p^3\) | \(\kappa_2(X)\) |
| SVD | \(X=U\Sigma V^\top\), \(\hat\beta=V\Sigma^{+}U^\top\mathbf{Y}\) | \(\approx 14np^2\) | \(\kappa_2(X)\), handles rank deficiency |
The normal equations are roughly twice as fast as
QR. They are also the only route whose error grows with the
square of the condition number (§3.7.1). lm() uses QR,
deliberately, and accepts the factor of two.
solve_normal <- function(X, y) as.vector(solve(crossprod(X), crossprod(X, y)))
solve_qr <- function(X, y) { q <- qr(X); as.vector(backsolve(qr.R(q), crossprod(qr.Q(q), y))) }
solve_svd <- function(X, y) {
s <- svd(X); d <- s$d; keep <- d > max(dim(X)) * .Machine$double.eps * max(d)
as.vector(s$v[, keep, drop = FALSE] %*% ((crossprod(s$u[, keep, drop = FALSE], y)) / d[keep]))
}Now build a design whose conditioning we control, with a known true coefficient vector, so the error is measurable rather than inferred.
make_design <- function(n, p, kappa, seed = 67) {
set.seed(seed)
U <- qr.Q(qr(matrix(rnorm(n * p), n)))
V <- qr.Q(qr(matrix(rnorm(p * p), p)))
d <- 10^seq(0, -log10(kappa), length.out = p) # singular values spanning kappa
U %*% diag(d) %*% t(V)
}
# experiment <- do.call(rbind, lapply(10^c(2, 4, 6, 8, 10), function(kap) {
# Xk <- make_design(400, 8, kap)
# beta_true <- rep(1, 8)
# yk <- Xk %*% beta_true # NOISELESS: any error is numerical
# rel <- function(b) max(abs(b - beta_true)) / max(abs(beta_true))
# data.frame(kappa_X = kap,
# kappa_XtX = kappa(crossprod(Xk), exact = TRUE),
# normal_eqns = rel(solve_normal(Xk, yk)),
# QR = rel(solve_qr(Xk, yk)),
# SVD = rel(solve_svd(Xk, yk)),
# lm_fit = rel(unname(lm.fit(Xk, yk)$coefficients)))
# }))
# experiment |> mutate(across(everything(), \(z) signif(z, 3)))
## Tis revision prevents the loop from crashing when $\kappa$ gets too large,
## using tryCatch() inside the lapply loop, to gracefully record an error or
## assign NA for the normal equations method at extreme condition numbers
## while completing successfully the QR and SVD methods
experiment <- do.call(rbind, lapply(10^c(2, 4, 6, 8, 10), function(kap) {
Xk <- make_design(400, 8, kap)
beta_true <- rep(1, 8)
yk <- Xk %*% beta_true # NOISELESS
rel <- function(b) max(abs(b - beta_true)) / max(abs(beta_true))
# Safely handle normal equations when matrix becomes computationally singular
normal_res <- tryCatch(
rel(solve_normal(Xk, yk)),
error = function(e) NA_real_
)
data.frame(
kappa_X = kap,
kappa_XtX = kappa(crossprod(Xk), exact = TRUE),
normal_eqns = normal_res,
QR = rel(solve_qr(Xk, yk)),
SVD = rel(solve_svd(Xk, yk)),
lm_fit = rel(unname(lm.fit(Xk, yk)$coefficients))
)
}))
#
experiment <- experiment |> mutate(across(everything(), \(z) signif(z, 3)))
print(experiment)#> kappa_X kappa_XtX normal_eqns QR SVD lm_fit
#> 1 1e+02 1.00e+04 2.89e-13 1.29e-14 8.44e-15 3.33e-15
#> 2 1e+04 1.00e+08 1.90e-08 2.70e-13 4.08e-13 3.13e-13
#> 3 1e+06 1.00e+12 1.91e-04 4.94e-11 4.43e-11 3.98e-12
#> 4 1e+08 6.29e+15 NA 1.70e+00 1.88e-09 NA
#> 5 1e+10 2.13e+16 NA 1.04e+00 2.89e-07 NAexperiment |>
dplyr::select(kappa_X, normal_eqns, QR, SVD) |>
pivot_longer(-kappa_X, names_to = "method", values_to = "error") |>
mutate(error = pmax(error, 1e-17)) |>
ggplot(aes(kappa_X, error, colour = method)) +
geom_line(linewidth = 1) + geom_point(size = 2.2) +
scale_x_log10(labels = scales::label_log()) +
scale_y_log10(labels = scales::label_log()) +
scale_colour_manual(values = c(normal_eqns = "firebrick", QR = "steelblue",
SVD = "darkgreen")) +
labs(title = "Accuracy of three least-squares solvers versus conditioning",
subtitle = "Noiseless data, so every deviation from beta = 1 is numerical error",
x = expression(kappa[2](X)), y = "Relative error in beta-hat",
colour = NULL) +
theme_dspa()
# --- Interactive equivalent ------------------------------------------------
plot_ly(experiment, x = ~kappa_X, y = ~normal_eqns, type = "scatter",
mode = "lines+markers", name = "Normal equations") |>
add_trace(y = ~QR, name = "QR") |>
add_trace(y = ~SVD, name = "SVD") |>
layout(title = "Least-squares solver accuracy vs. condition number",
xaxis = list(type = "log", title = "kappa(X)"),
yaxis = list(type = "log", title = "relative error"),
legend = list(orientation = "h"))The normal-equations curve rises at twice the slope of the others on a log-log scale, the \(\kappa^2\) penalty, measured. By \(\kappa_2(X)=10^{8}\) the normal equations have no correct digits left, while QR and SVD are still accurate to about eight.
set.seed(71)
Xt <- cbind(1, matrix(rnorm(20000 * 30), 20000)); yt <- rnorm(20000)
timings <- sapply(list(normal_eqns = \() solve_normal(Xt, yt),
QR = \() solve_qr(Xt, yt),
SVD = \() solve_svd(Xt, yt),
lm.fit = \() lm.fit(Xt, yt),
lm = \() lm(yt ~ Xt - 1)),
\(f) system.time(f())[["elapsed"]])
round(timings, 4)#> normal_eqns QR SVD lm.fit lm
#> 0.00 0.03 0.03 0.01 0.02Common misconception: “the manual matrix formula is faster than
lm(), so use it.” It is faster, for two reasons that are both bad reasons to prefer it. First,lm()does far more work, it builds a model frame, expands factors into contrasts, and computes standard errors, fitted values, residuals, and the entiresummary()apparatus.lm.fit()is the like-for-like comparison. Second, the speed advantage of the normal equations comes precisely from skipping the factorization that protects accuracy. In production, uselm()for modelling andqr.solve()orlm.fit()when you need only the coefficients.
Xb <- cbind(1, mlb$Height)
yb <- mlb$Weight
rbind(normal_equations = solve_normal(Xb, yb),
QR = solve_qr(Xb, yb),
lm = unname(coef(lm(Weight ~ Height, data = mlb))))#> [,1] [,2]
#> normal_equations -154.224 4.82977
#> QR -154.224 4.82977
#> lm -154.224 4.82977ggplot(mlb, aes(Height, Weight)) +
geom_point(alpha = 0.18, size = 1.1, colour = "grey25") +
geom_smooth(method = "lm", formula = y ~ x, colour = "firebrick",
fill = "grey70", linewidth = 1) +
labs(title = "MLB players: weight regressed on height",
subtitle = sprintf("Weight-hat = %.1f + %.2f x Height | r = %.3f | n = %d",
coef(lm(Weight ~ Height, mlb))[1],
coef(lm(Weight ~ Height, mlb))[2],
cor(mlb$Height, mlb$Weight), nrow(mlb)),
x = "Height (in)", y = "Weight (lb)") +
theme_dspa()
# --- Interactive equivalent ------------------------------------------------
fitted_manual <- Xb %*% solve_qr(Xb, yb)
plot_ly(x = ~mlb$Height) |>
add_markers(y = ~mlb$Weight, name = "Data") |>
add_lines(x = ~mlb$Height, y = ~fitted_manual[, 1],
name = "Manual least squares (QR)") |>
add_lines(x = ~mlb$Height, y = ~fitted(lm(Weight ~ Height, mlb)),
name = "lm(Weight ~ Height)",
line = list(width = 4, dash = "dash")) |>
layout(title = "Baseball Players: Linear Model of Weight vs. Height",
xaxis = list(title = "Height (in)"),
yaxis = list(title = "Weight (lb)"),
legend = list(orientation = "h"))15 Covariance as a matrix computation
For \(X_{n\times k}\) with column means \(\bar{\mathbf{x}}\),
\[\Sigma=\frac{1}{n-1}\big(X-\mathbf{1}\bar{\mathbf{x}}^\top\big)^\top\big(X-\mathbf{1}\bar{\mathbf{x}}^\top\big), \qquad \Sigma_{ij}=\frac{1}{n-1}\sum_{m=1}^{n}(x_{mi}-\bar x_i)(x_{mj}-\bar x_j).\]
\(\Sigma\) is symmetric positive semi-definite, with variances on the diagonal and covariances off it.
Xc <- matrix(c(4.0, 4.2, 3.9, 4.3, 4.1,
2.0, 2.1, 2.0, 2.1, 2.2,
0.60, 0.59, 0.58, 0.62, 0.63), ncol = 3)
Xcentred <- sweep(Xc, 2, colMeans(Xc)) # clearer than rep()/matrix() gymnastics
S_manual <- crossprod(Xcentred) / (nrow(Xc) - 1)
round(S_manual, 6)#> [,1] [,2] [,3]
#> [1,] 0.02500 0.00750 0.00175
#> [2,] 0.00750 0.00700 0.00135
#> [3,] 0.00175 0.00135 0.00043#> [,1] [,2] [,3]
#> [1,] 0.02500 0.00750 0.00175
#> [2,] 0.00750 0.00700 0.00135
#> [3,] 0.00175 0.00135 0.00043#> [1] TRUECorrelation rescales covariance by the standard deviations:
\[\rho_{ij}=\frac{\Sigma_{ij}}{\sqrt{\Sigma_{ii}\Sigma_{jj}}}, \qquad\text{equivalently}\qquad P=D^{-1/2}\Sigma D^{-1/2},\ \ D=\operatorname{diag}(\Sigma).\]
#> [,1] [,2] [,3]
#> [1,] 1.0000 0.5303 0.1578
#> [2,] 0.5303 1.0000 -0.0737
#> [3,] 0.1578 -0.0737 1.0000#> Weight Height Age
#> Weight 1.0000 0.5303 0.1578
#> Height 0.5303 1.0000 -0.0737
#> Age 0.1578 -0.0737 1.0000Common misconception: “\(\operatorname{cov}(x,x)=1\).” \(\operatorname{cov}(x,x)=\operatorname{var}(x)\), the diagonal of a covariance matrix holds variances, in squared units of each variable. It is the correlation that satisfies \(\operatorname{cor}(x,x)=1\), because dividing by \(\sqrt{\operatorname{var}(x)\operatorname{var}(x)}\) normalizes it. Both matrices are symmetric; only the correlation matrix has unit diagonal.
Interpreting the magnitude of a correlation requires context, not a lookup table. The bands \(|\rho|<0.1\) negligible, \(0.1\)–\(0.3\) weak, \(0.3\)–\(0.5\) moderate, \(>0.5\) strong come from Cohen’s conventions for behavioural science; they are a starting point for discussion, not a property of the number. A correlation of \(0.2\) between a biomarker and a clinical outcome may be highly actionable; a correlation of \(0.7\) between two instruments measuring the same construct is disappointing.
16 Regression diagnostics
Every diagnostic in this section is a function of the hat matrix \(H=X(X^\top X)^{-1}X^\top\). Knowing that turns a checklist into a set of derivations.
16.1 Assumptions, and what breaks when each fails
| Assumption | Diagnostic | Consequence of violation | Remedy |
|---|---|---|---|
| Linearity in \(\boldsymbol\beta\) | Residuals vs. fitted | Biased \(\hat{\boldsymbol\beta}\) | Transform; add terms; nonlinear model |
| \(E[\varepsilon\mid X]=0\) | Residual patterns; theory | Biased \(\hat{\boldsymbol\beta}\) | Omitted variables; instruments |
| Homoscedasticity | Scale–location plot | SEs wrong, \(\hat{\boldsymbol\beta}\) still unbiased | Robust (HC) SEs; WLS; transform |
| Uncorrelated errors | Durbin–Watson; ACF of residuals | SEs wrong | GLS; cluster-robust SEs; time-series model |
| Normality of \(\varepsilon\) | Q-Q of standardized residuals | Exact \(t\)/\(F\) invalid (asymptotics often save you) | Bootstrap; robust regression |
| \(\operatorname{rank}(X)=p\) | VIF; \(\kappa_2(X)\) | No unique solution | Drop aliased columns; regularize |
Note the asymmetry: heteroscedasticity and error correlation break the standard errors but leave \(\hat{\boldsymbol\beta}\) unbiased. Nonlinearity and omitted variables break the estimates themselves. The second failure is far more serious and far less often checked.
16.2 Leverage
The diagonal entry \(h_{ii}=\mathbf{x}_i^\top(X^\top X)^{-1}\mathbf{x}_i\) is the leverage of case \(i\), the sensitivity of its own fitted value to its own response, \(h_{ii}=\partial\hat y_i/\partial y_i\). From §3.4.2:
\[0\le h_{ii}\le 1,\qquad \sum_{i=1}^{n}h_{ii}=\operatorname{tr}(H)=p, \qquad\text{so}\qquad \bar h=\frac{p}{n}.\]
Leverage depends only on \(X\), never on \(\mathbf{Y}\): it measures how unusual a case’s predictor values are. The conventional flag is \(h_{ii}>2p/n\).
16.3 Residuals, standardized and studentized
Since \(\hat{\boldsymbol\varepsilon}=(I-H)\mathbf{Y}\) and \(I-H\) is symmetric idempotent,
\[\operatorname{Var}(\hat\varepsilon_i)=\sigma^2(1-h_{ii}).\]
Raw residuals therefore have unequal variances, high-leverage points get artificially small residuals, which is exactly backwards for outlier detection. Standardizing fixes it:
\[r_i=\frac{\hat\varepsilon_i}{\hat\sigma\sqrt{1-h_{ii}}} \qquad\text{(internally studentized)},\]
\[t_i=\frac{\hat\varepsilon_i}{\hat\sigma_{(i)}\sqrt{1-h_{ii}}}\sim t_{n-p-1} \qquad\text{(externally studentized)},\]
where \(\hat\sigma_{(i)}\) omits case \(i\). Only \(t_i\) has an exact null distribution, which is what makes it usable as a test.
16.4 Cook’s distance
Cook’s distance measures how far all fitted values move when case \(i\) is deleted:
\[D_i=\frac{\sum_{j}(\hat y_j-\hat y_{j(i)})^2}{p\,\hat\sigma^2} =\frac{r_i^2}{p}\cdot\frac{h_{ii}}{1-h_{ii}} .\]
The factorization is the useful part: influence = outlyingness \(\times\) leverage. A point is influential only if it is both poorly fitted and extreme in predictor space. High leverage alone is harmless; a large residual at low leverage moves nothing. Conventional flags: \(D_i>4/n\), or \(D_i>1\).
mlb2 <- mlb |> dplyr::select(-Name) |>
mutate(Team = factor(Team), Position = factor(Position))
fit_mlb <- lm(Weight ~ Height + Age, data = mlb2)
diag_df <- data.frame(
fitted = fitted(fit_mlb),
resid = residuals(fit_mlb),
std_res = rstandard(fit_mlb),
stu_res = rstudent(fit_mlb),
leverage = hatvalues(fit_mlb),
cooks = cooks.distance(fit_mlb),
idx = seq_len(nrow(mlb2))
)
p_lm <- length(coef(fit_mlb)); n_lm <- nrow(mlb2)
thresh <- c(leverage = 2 * p_lm / n_lm, cooks = 4 / n_lm)
round(thresh, 5)#> leverage cooks
#> 0.00580 0.00387# Identify by THRESHOLD, never by hard-coded index
diag_df |>
filter(cooks > thresh[["cooks"]]) |>
arrange(desc(cooks)) |> head(6) |>
mutate(across(where(is.numeric), \(z) round(z, 4)))p1 <- ggplot(diag_df, aes(fitted, resid)) +
geom_point(alpha = 0.25, size = 1) +
geom_hline(yintercept = 0, colour = "grey50") +
geom_smooth(method = "loess", formula = y ~ x, se = FALSE,
colour = "firebrick", linewidth = 0.8) +
labs(title = "Residuals vs. fitted", subtitle = "Curvature signals a missing term",
x = "Fitted", y = "Residual") + theme_dspa(9)
p2 <- ggplot(diag_df, aes(sample = std_res)) +
stat_qq(size = 0.8, alpha = 0.4) + stat_qq_line(colour = "firebrick") +
labs(title = "Normal Q-Q of STANDARDIZED residuals",
subtitle = "Theoretical normal quantiles on the x axis",
x = "Theoretical quantiles", y = "Standardized residuals") + theme_dspa(9)
p3 <- ggplot(diag_df, aes(fitted, sqrt(abs(std_res)))) +
geom_point(alpha = 0.25, size = 1) +
geom_smooth(method = "loess", formula = y ~ x, se = FALSE,
colour = "firebrick", linewidth = 0.8) +
labs(title = "Scale-location", subtitle = "A trend indicates heteroscedasticity",
x = "Fitted", y = expression(sqrt(abs(standardized~residual)))) + theme_dspa(9)
p4 <- ggplot(diag_df, aes(leverage, std_res)) +
geom_point(aes(size = cooks), alpha = 0.3) +
geom_hline(yintercept = c(-2, 2), linetype = "dashed", colour = "grey50") +
geom_vline(xintercept = thresh[["leverage"]], linetype = "dashed",
colour = "firebrick") +
scale_size_continuous(range = c(0.4, 4)) +
labs(title = "Residuals vs. leverage",
subtitle = "Point size is Cook's D = outlyingness x leverage",
x = expression(h[ii]), y = "Standardized residual", size = "Cook's D") +
theme_dspa(9)
(p1 | p2) / (p3 | p4)
# --- Interactive equivalents ----------------------------------------------
plot_ly(diag_df, x = ~fitted, y = ~resid, type = "scatter", mode = "markers",
text = ~paste0("case ", idx), marker = list(opacity = 0.4)) |>
layout(title = "Fitted values vs. model residuals",
xaxis = list(title = "Fitted"), yaxis = list(title = "Residuals"))
qq <- qqnorm(diag_df$std_res, plot.it = FALSE)
plot_ly(x = qq$x, y = qq$y, type = "scatter", mode = "markers",
name = "Standardized residuals") |>
add_lines(x = range(qq$x), y = range(qq$x), name = "y = x",
line = list(color = "red", width = 3)) |>
layout(title = "Normal Q-Q plot",
xaxis = list(title = "Theoretical quantiles"),
yaxis = list(title = "Standardized residuals"))# Half-normal plot of leverages, built inline (no extra package dependency)
hv <- sort(diag_df$leverage)
nh <- length(hv)
hn <- data.frame(theoretical = qnorm((nh + seq_len(nh)) / (2 * nh + 1)),
leverage = hv,
label = order(diag_df$leverage))
top2 <- tail(hn, 2)
ggplot(hn, aes(theoretical, leverage)) +
geom_point(size = 1, alpha = 0.5) +
geom_text(data = top2, aes(label = label), hjust = 1.3, size = 3,
colour = "firebrick") +
labs(title = "Half-normal plot of leverages",
subtitle = "The two largest are labelled by case index",
x = "Half-normal quantiles", y = expression(h[ii])) +
theme_dspa()
17 Multicollinearity
17.1 What OLS actually requires
Common misconception: “multiple regression assumes the predictors are independent of each other.” It does not. The requirement is \(\operatorname{rank}(X)=p\), no column may be an exact linear combination of the others. Correlated predictors are entirely permissible. Under A1–A4 the estimator remains unbiased no matter how correlated the columns are; what changes is \(\operatorname{Var}(\hat{\boldsymbol\beta})=\sigma^2(X^\top X)^{-1}\), which grows as \(X^\top X\) approaches singularity.
Believing the false version leads to dropping informative predictors because they correlate with others, trading a variance problem you could have reported for an omitted-variable bias you cannot.
Geometrically: when two columns are nearly parallel, the RSS bowl of §3.9.4 becomes a long narrow valley. Every point along the valley floor fits almost equally well, so the estimates shift substantially under small perturbations of the data, while the fitted values stay stable. That is the signature of multicollinearity: unstable coefficients, stable predictions.
sim_collinear <- function(rho, n = 200, reps = 400, seed = 79) {
set.seed(seed)
Sig <- matrix(c(1, rho, rho, 1), 2)
L <- t(chol(Sig))
out <- t(replicate(reps, {
Z <- t(L %*% matrix(rnorm(2 * n), 2))
y <- 1 + 2 * Z[, 1] + 3 * Z[, 2] + rnorm(n)
b <- coef(lm(y ~ Z[, 1] + Z[, 2]))
c(b1 = b[2], b2 = b[3], fit_sd = sd(fitted(lm(y ~ Z[, 1] + Z[, 2]))))
}))
data.frame(rho = rho,
mean_b1 = mean(out[, 1]), sd_b1 = sd(out[, 1]),
mean_b2 = mean(out[, 2]), sd_b2 = sd(out[, 2]),
VIF = 1 / (1 - rho^2))
}
do.call(rbind, lapply(c(0, 0.5, 0.9, 0.99), sim_collinear)) |>
mutate(across(where(is.numeric), \(z) round(z, 4)))The means stay at the true values 2 and 3 for every \(\rho\), unbiased. The standard deviations inflate by exactly \(\sqrt{\mathrm{VIF}}\).
17.2 The variance inflation factor (VIF)
With a single predictor, the minimum attainable variance of \(\hat\beta_\ell\) is
\[\operatorname{Var}(\hat\beta_\ell)_{\min}=\frac{\sigma^2}{\sum_{i=1}^{n}(x_{i\ell}-\bar x_\ell)^2}.\]
In a multiple regression, the actual variance is
\[\operatorname{Var}(\hat\beta_\ell)=\frac{\sigma^2}{\sum_{i=1}^{n}(x_{i\ell}-\bar x_\ell)^2}\cdot\frac{1}{1-R_\ell^2},\]
where \(R_\ell^2\) is the \(R^2\) from regressing \(x_\ell\) on the other \(p-2\) predictors. The ratio is the variance inflation factor
\[\boxed{\;\mathrm{VIF}_\ell=\frac{\operatorname{Var}(\hat\beta_\ell)}{\operatorname{Var}(\hat\beta_\ell)_{\min}}=\frac{1}{1-R_\ell^2}\;}\]
so the standard error is inflated by \(\sqrt{\mathrm{VIF}_\ell}\). Conventional flags are \(\mathrm{VIF}>4\) (worth noting) and \(\mathrm{VIF}>10\) (worth acting on), but these are conventions, not tests, and their practical meaning depends on how large \(\sigma^2\) and \(n\) are.
#> Height Age
#> 1.00546 1.00546# VIF from first principles, to make the definition concrete
R2_height <- summary(lm(Height ~ Age, data = mlb2))$r.squared
c(vif_height_manual = 1 / (1 - R2_height))#> vif_height_manual
#> 1.00546Remedies, in rough order of preference. Accept it and report wider intervals (the estimates are still unbiased); centre variables before forming polynomial or interaction terms, which often removes most of the problem; regularize with ridge, which trades a little bias for a large variance reduction (Chapter 11); combine collinear predictors into a single index or principal component (Chapter 4). Dropping a predictor is a last resort, because it converts a variance problem into a bias problem.
17.3 Spurious correlation when \(k \gg n\)
In high dimensions, strong sample correlations arise with certainty among variables that are independent in the population. For \(k\) IID standard Gaussian columns of length \(n\), extreme-value theory gives
\[\max_{i\ne j}\big|\hat\rho_{ij}\big|\;\approx\;2\sqrt{\frac{\log k}{n}}\,,\]
which tends to 1 as \(k\) grows with \(n\) fixed.
max_cor_sim <- function(n, k, reps = 20, seed = 83) {
set.seed(seed)
mean(replicate(reps, {
Z <- matrix(rnorm(n * k), n)
R <- cor(Z)
max(abs(R[upper.tri(R)]))
}))
}
grid <- expand.grid(n = c(50, 100), k = c(10, 50, 200, 1000))
grid$observed <- mapply(max_cor_sim, grid$n, grid$k)
grid$predicted <- 2 * sqrt(log(grid$k) / grid$n)
grid |> mutate(across(c(observed, predicted), \(z) round(z, 3)))ggplot(grid, aes(k, observed, colour = factor(n))) +
geom_line(linewidth = 1) + geom_point(size = 2.2) +
geom_line(aes(y = pmin(predicted, 1)), linetype = "dashed", linewidth = 0.7) +
scale_x_log10() +
labs(title = "Maximum sample correlation among INDEPENDENT Gaussian variables",
subtitle = "Solid: simulated. Dashed: 2 sqrt(log k / n)",
x = "Number of variables k (log scale)",
y = "Max |correlation|", colour = "n") +
theme_dspa()
At \(n=50\), \(k=1000\), modest by genomics standards, the largest pairwise correlation among completely independent variables averages above 0.5. Any procedure that screens features by marginal correlation will select noise.
What to do when \(k>n\). VIF is undefined (\(X^\top X\) is singular), so:
- Dimension reduction, PCA, PLS, ICA, or factor analysis reduce to \(k'<n\) components (Chapter 4).
- Regularization, ridge and LASSO give unique solutions even when \(X^\top X\) is singular; ridge does so by adding \(\lambda I\), which shifts every eigenvalue away from zero.
- Sure Independence Screening (Fan & Lv, 2008) ranks predictors by marginal correlation and reduces \(k\) to \(O(n)\) before applying a careful method. Nonparametric and model-free variants relax the linearity of the screening step.
- Stability selection, repeat selection over bootstrap resamples and keep what recurs, controlling the expected number of false selections.
18 Case study: MLB players
18.1 The data
The MLB
dataset records 1,034 players: Name, Team,
Position, Height (in), Weight
(lb), Age (yr).
#> 'data.frame': 1034 obs. of 5 variables:
#> $ Team : Factor w/ 30 levels "ANA","ARZ","ATL",..: 4 4 4 4 4 4 4 4 4 4 ...
#> $ Position: Factor w/ 9 levels "Catcher","Designated_Hitter",..: 1 1 1 3 3 6 7 9 9 4 ...
#> $ Height : int 74 74 72 72 73 69 69 71 76 71 ...
#> $ Weight : int 180 215 210 210 188 176 209 200 231 180 ...
#> $ Age : num 23 34.7 30.8 35.4 35.7 ...#> Weight Height Age
#> Min. :150 Min. :67.0 Min. :20.9
#> 1st Qu.:187 1st Qu.:72.0 1st Qu.:25.4
#> Median :200 Median :74.0 Median :27.9
#> Mean :202 Mean :73.7 Mean :28.7
#> 3rd Qu.:215 3rd Qu.:75.0 3rd Qu.:31.2
#> Max. :290 Max. :83.0 Max. :48.5#>
#> Catcher Designated_Hitter First_Baseman Outfielder
#> 76 18 55 194
#> Relief_Pitcher Second_Baseman Shortstop Starting_Pitcher
#> 315 58 52 221
#> Third_Baseman
#> 45ggplot(mlb2, aes(Weight)) +
geom_histogram(bins = nclass.FD(mlb2$Weight), fill = "steelblue",
colour = "white") +
geom_vline(xintercept = c(mean(mlb2$Weight), median(mlb2$Weight)),
colour = c("firebrick", "darkgreen"), linetype = c(1, 2),
linewidth = 0.9) +
labs(title = "Distribution of player weight",
subtitle = "Red solid: mean. Green dashed: median. Mean > median indicates right skew",
x = "Weight (lb)", y = "Count") +
theme_dspa()
# --- Interactive equivalent ------------------------------------------------
plot_ly(x = mlb2$Weight, type = "histogram", name = "Weight") |>
layout(title = "Baseball Players' Weight Histogram", bargap = 0.1,
xaxis = list(title = "Weight (lb)"),
yaxis = list(title = "Frequency"))# GGally::ggpairs(mlb2, columns = c("Height", "Weight", "Age"),
# lower = list(continuous = GGally::wrap("points", size = 0.4,
# alpha = 0.25)),
# upper = list(continuous = GGally::wrap("cor", size = 3.4))) +
# labs(title = "MLB pairs plot") + theme_dspa(9)
GGally::ggpairs(
mlb2,
columns = c("Height", "Weight", "Age"),
title = "MLB pairs plot", # Add title here
lower = list(continuous = GGally::wrap("points", size = 0.4, alpha = 0.25)),
upper = list(continuous = GGally::wrap("cor", size = 3.4))
)
Because a scatterplot matrix invites brushing, selecting points in one panel to see where they land in the others, the interactive version carries information the static one cannot, so it is evaluated:
plot_ly(mlb2) |>
add_trace(type = "splom",
dimensions = list(list(label = "Height", values = ~Height),
list(label = "Weight", values = ~Weight),
list(label = "Age", values = ~Age)),
text = ~Position,
marker = list(color = as.integer(mlb2$Position), size = 4,
opacity = 0.55,
line = list(width = 0.5, color = "rgb(230,230,230)"))) |>
layout(title = "MLB pairs plot (drag to brush; colour = position)",
hovermode = "closest", dragmode = "select",
plot_bgcolor = "rgba(240,240,240,0.95)")#> Weight Height Age
#> Weight 1.0000 0.5303 0.1578
#> Height 0.5303 1.0000 -0.0737
#> Age 0.1578 -0.0737 1.0000Height and weight are strongly positively associated; age relates only weakly to either. No pair is near-collinear, so no identifiability problem arises.
18.2 Fitting
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) -164.9995 19.3828 -8.5127 0.0000
#> TeamARZ 7.1881 4.2590 1.6877 0.0918
#> TeamATL -1.5631 3.9757 -0.3932 0.6943
#> TeamBAL -5.3128 4.0193 -1.3218 0.1865
#> TeamBOS -0.2838 4.0034 -0.0709 0.9435
#> TeamCHC 0.4026 3.9949 0.1008 0.9197
#> TeamCIN 2.1051 3.9934 0.5271 0.5982
#> TeamCLE -1.3160 4.0356 -0.3261 0.7444c(n = nobs(fit_all), coefficients = length(coef(fit_all)),
R2 = summary(fit_all)$r.squared,
adj_R2 = summary(fit_all)$adj.r.squared,
sigma_hat = summary(fit_all)$sigma)#> n coefficients R2 adj_R2 sigma_hat
#> 1034.000000 40.000000 0.385804 0.361706 16.777432Factor predictors are expanded by treatment
contrasts: a \(k\)-level
factor contributes \(k-1\) indicator
columns, and each coefficient is a difference from the reference level.
That is why Team alone adds about thirty coefficients.
A caution on reading the stars. With thirty-plus coefficients tested at \(\alpha=0.05\), roughly one or two will be “significant” under a complete null. Individual \(p\)-values in a wide model are not independent evidence about individual predictors; use the \(F\)-test for the factor as a whole (
anova()), or control the false discovery rate (Chapter 11).
The Team term is tested with a single \(F\)-statistic on its full set of degrees of
freedom, which is the right question, “does team membership explain
weight?”, rather than thirty separate ones.
19 Case study: heart-attack hospitalization
#> 'data.frame': 148 obs. of 8 variables:
#> $ Patient : int 1 2 3 4 5 6 7 8 9 10 ...
#> $ DIAGNOSIS: int 41041 41041 41091 41081 41091 41091 41091 41091 41041 41041 ...
#> $ SEX : chr "F" "F" "F" "F" ...
#> $ DRG : int 122 122 122 122 122 121 121 121 121 123 ...
#> $ DIED : int 0 0 0 0 0 0 0 0 0 1 ...
#> $ CHARGES : num 4752 3941 3657 1481 1681 ...
#> $ LOS : int 10 6 5 2 1 9 15 15 2 1 ...
#> $ AGE : int 79 34 76 80 55 84 84 70 76 65 ...#> n complete
#> 148 148ggplot(ha, aes(LOS, CHARGES)) +
geom_point(alpha = 0.55, colour = "steelblue", size = 1.8) +
geom_smooth(method = "lm", formula = y ~ x, colour = "firebrick",
fill = "grey75") +
geom_point(aes(x = mean(LOS), y = mean(CHARGES)), colour = "darkblue",
size = 5, shape = 21, fill = "gold", stroke = 1.2) +
labs(title = "Hospital charges versus length of stay",
subtitle = sprintf("CHARGES-hat = %.2f + %.2f x LOS | r = %.3f | the fitted line passes through (mean LOS, mean CHARGES)",
coef(fit_h)[1], coef(fit_h)[2],
cor(ha$LOS, ha$CHARGES)),
x = "Length of stay (days)", y = "Charges ($)") +
theme_dspa()
# --- Interactive equivalent ------------------------------------------------
plot_ly(ha, x = ~LOS, y = ~CHARGES, type = "scatter", mode = "markers",
name = "Data") |>
add_trace(x = ~mean(LOS), y = ~mean(CHARGES), type = "scatter",
mode = "markers", name = "(mean LOS, mean CHARGES)",
marker = list(size = 20, color = "blue",
line = list(color = "yellow", width = 2))) |>
add_lines(x = ~LOS, y = fitted(fit_h), name = "Linear model") |>
layout(title = paste0("lm(CHARGES ~ LOS), Cor(LOS, CHARGES) = ",
round(cor(ha$LOS, ha$CHARGES), 3)))That the least-squares line passes exactly through \((\bar x,\bar y)\) is not a coincidence: it follows from the intercept \(\hat a=\bar y-\hat b\bar x\), which is itself forced by the first normal equation \(\mathbf{1}^\top(\mathbf{y}-X\hat{\boldsymbol\beta})=0\), the residuals of any model with an intercept sum to zero.
19.1 Simple regression from the covariance
For a single predictor the closed form is
\[\hat b=\frac{\sum_i(x_i-\bar x)(y_i-\bar y)}{\sum_i(x_i-\bar x)^2} =\frac{\widehat{\operatorname{Cov}}(x,y)}{\widehat{\operatorname{Var}}(x)}, \qquad \hat a=\bar y-\hat b\bar x .\]
b_cov <- cov(ha$LOS, ha$CHARGES) / var(ha$LOS)
a_cov <- mean(ha$CHARGES) - b_cov * mean(ha$LOS)
rbind(from_covariance = c(intercept = a_cov, slope = b_cov),
from_lm = unname(coef(fit_h)))#> intercept slope
#> from_covariance 4582.7 212.287
#> from_lm 4582.7 212.287A patient staying 10 days has predicted charges \(\hat a+10\hat b\):
#> fit lwr upr
#> 1 6705.57 -3.22108 13414.4Report the prediction interval, not just the point estimate. It is far wider than the confidence interval for the mean response, because it must absorb the irreducible \(\sigma^2\) of an individual outcome in addition to uncertainty about the line:
\[\operatorname{Var}(\hat y_{\text{new}})=\sigma^2\Big(\underbrace{1}_{\text{new observation}}+\underbrace{\mathbf{x}_0^\top(X^\top X)^{-1}\mathbf{x}_0}_{\text{uncertainty in }\hat{\boldsymbol\beta}}\Big).\]
newx <- data.frame(LOS = seq(min(ha$LOS), max(ha$LOS), length.out = 120))
ci <- as.data.frame(predict(fit_h, newx, interval = "confidence"))
pi <- as.data.frame(predict(fit_h, newx, interval = "prediction"))
bands <- cbind(newx, fit = ci$fit, ci_lo = ci$lwr, ci_hi = ci$upr,
pi_lo = pi$lwr, pi_hi = pi$upr)
ggplot(bands, aes(LOS)) +
geom_ribbon(aes(ymin = pi_lo, ymax = pi_hi, fill = "95% prediction"), alpha = 0.25) +
geom_ribbon(aes(ymin = ci_lo, ymax = ci_hi, fill = "95% confidence"), alpha = 0.45) +
geom_line(aes(y = fit), colour = "firebrick", linewidth = 1) +
geom_point(data = ha, aes(LOS, CHARGES), alpha = 0.45, size = 1.5,
inherit.aes = FALSE) +
scale_fill_manual(values = c("95% prediction" = "steelblue",
"95% confidence" = "orange")) +
labs(title = "Confidence band for the mean vs. prediction band for an individual",
x = "Length of stay (days)", y = "Charges ($)", fill = NULL) +
theme_dspa()
19.2 Multiple predictors and a de novo estimator
# A small OLS routine that fails informatively on rank-deficient designs
reg <- function(y, x) {
X <- cbind(Intercept = 1, as.matrix(x))
r <- Matrix::rankMatrix(X)[1]
if (r < ncol(X))
stop(sprintf("Design is rank deficient: rank %d < %d columns.", r, ncol(X)),
call. = FALSE)
qrX <- qr(X) # QR, not the normal equations
drop(backsolve(qr.R(qrX), crossprod(qr.Q(qrX), y)))
}
rbind(reg_manual = reg(ha$CHARGES, ha[, c("LOS", "AGE")]),
lm = unname(coef(lm(CHARGES ~ LOS + AGE, data = ha))))#> [,1] [,2] [,3]
#> reg_manual 7280.55 259.674 -43.6768
#> lm 7280.55 259.674 -43.6768#> Error : Design is rank deficient: rank 2 < 3 columns.20 Baseball offense: multiple regression in 3-D
The Lahman database records
team-season statistics. We model runs scored (R) from walks
(BB), singles, and home runs (HR) for full
162-game seasons before 2002.
Teams <- dspa_read("https://umich.instructure.com/files/2798317/download?download_frd=1",
"teamsData.csv", header = TRUE)
dat <- Teams |>
filter(G == 162, yearID < 2002) |>
mutate(Singles = H - X2B - X3B - HR) |>
dplyr::select(name, R, Singles, HR, BB)
dim(dat); head(dat, 4)#> [1] 663 5fit_bb <- lm(R ~ BB, data = dat)
fit_bh <- lm(R ~ BB + HR, data = dat)
fit_all3 <- lm(R ~ BB + Singles + HR, data = dat)
data.frame(
model = c("R ~ BB", "R ~ BB + HR", "R ~ BB + Singles + HR"),
R2 = c(summary(fit_bb)$r.squared, summary(fit_bh)$r.squared,
summary(fit_all3)$r.squared),
adj_R2 = c(summary(fit_bb)$adj.r.squared, summary(fit_bh)$adj.r.squared,
summary(fit_all3)$adj.r.squared),
sigma = c(summary(fit_bb)$sigma, summary(fit_bh)$sigma,
summary(fit_all3)$sigma),
AIC = c(AIC(fit_bb), AIC(fit_bh), AIC(fit_all3))
) |> mutate(across(where(is.numeric), \(z) round(z, 3)))ggplot(dat, aes(BB, R)) +
geom_point(alpha = 0.35, size = 1.4, colour = "steelblue") +
geom_smooth(method = "lm", formula = y ~ x, colour = "firebrick") +
labs(title = "Runs scored versus walks drawn",
subtitle = sprintf("Team-seasons of 162 games before 2002; R-squared = %.3f",
summary(fit_bb)$r.squared),
x = "Walks by batters (BB)", y = "Runs scored (R)") +
theme_dspa()
# --- Interactive equivalent ------------------------------------------------
plot_ly(x = ~dat$BB) |>
add_markers(y = ~dat$R, name = "Data") |>
add_lines(x = ~dat$BB, y = ~fitted(fit_bb), name = "lm(R ~ BB)",
line = list(width = 4)) |>
layout(title = "Scatter plot / regression for baseball data",
xaxis = list(title = "(BB) Walks by batters"),
yaxis = list(title = "(R) Runs scored"),
legend = list(orientation = "h"))With two predictors the fitted object is a plane in
three dimensions, and a plane must be rotated to be understood. Note
that the surface below is generated from fit_bh, the same
model whose coefficients are reported, and is evaluated over the
observed range of the data.
cf <- coef(fit_bh)
ax_bb <- seq(min(dat$BB), max(dat$BB), length.out = 40)
ax_hr <- seq(min(dat$HR), max(dat$HR), length.out = 40)
# z[i, j] must have nrow = length(y) = length(ax_hr), ncol = length(x) = length(ax_bb)
zplane <- outer(ax_hr, ax_bb, function(hr, bb) cf[["(Intercept)"]] +
cf[["BB"]] * bb + cf[["HR"]] * hr)
plot_ly() |>
add_markers(data = dat, x = ~BB, y = ~HR, z = ~R, text = ~name,
type = "scatter3d", mode = "markers",
marker = list(size = 3, opacity = 0.6, color = "#1F77B4"),
name = "Team-seasons") |>
add_surface(x = ax_bb, y = ax_hr, z = zplane, opacity = 0.55,
showscale = FALSE, colorscale = list(c(0, "#E8A24C"), c(1, "#E8A24C")),
name = "Fitted plane") |>
layout(title = sprintf("Fitted plane R = %.1f + %.3f BB + %.3f HR (R-squared = %.3f)",
cf[1], cf[2], cf[3], summary(fit_bh)$r.squared),
scene = list(xaxis = list(title = "(BB) Walks"),
yaxis = list(title = "(HR) Home runs"),
zaxis = list(title = "(R) Runs scored")))# Residuals as vertical segments from each point to the plane make the
# geometry of "least squares" explicit in three dimensions.
sel <- sample(nrow(dat), 90)
fitted_bh <- fitted(fit_bh)
p <- plot_ly() |>
add_surface(x = ax_bb, y = ax_hr, z = zplane, opacity = 0.35,
showscale = FALSE,
colorscale = list(c(0, "#CCCCCC"), c(1, "#CCCCCC"))) |>
add_markers(x = dat$BB[sel], y = dat$HR[sel], z = dat$R[sel],
type = "scatter3d", mode = "markers", name = "Observed",
marker = list(size = 3.5, color = "#1F77B4"))
for (i in sel[1:45])
p <- add_trace(p, x = c(dat$BB[i], dat$BB[i]), y = c(dat$HR[i], dat$HR[i]),
z = c(dat$R[i], fitted_bh[i]), type = "scatter3d",
mode = "lines", showlegend = FALSE,
line = list(color = "firebrick", width = 2))
p |> layout(title = "Residuals are vertical distances to the fitted plane",
scene = list(xaxis = list(title = "BB"), yaxis = list(title = "HR"),
zaxis = list(title = "R")))21 Expanding the model
21.1 Polynomial terms
A quadratic in Age is still a linear
model, the design matrix simply gains an Age²
column.
Always centre before squaring. Raw
AgeandAge²correlate at about 0.99 over a narrow age range, which makes \(X^\top X\) nearly singular, the VIF enormous, and the individual coefficients uninterpretable. Centring makes the columns nearly orthogonal without changing the fit at all.
mlb3 <- mlb2 |> mutate(age_c = Age - mean(Age), age_c2 = age_c^2,
age_raw2 = Age^2)
c(cor_raw = cor(mlb3$Age, mlb3$age_raw2),
cor_centred = cor(mlb3$age_c, mlb3$age_c2))#> cor_raw cor_centred
#> 0.995536 0.538704fit_raw <- lm(Weight ~ Height + Age + age_raw2, data = mlb3)
fit_ctr <- lm(Weight ~ Height + age_c + age_c2, data = mlb3)
rbind(raw = car::vif(fit_raw),
centred = car::vif(fit_ctr)) |> round(2)#> Height Age age_raw2
#> raw 1.01 112.34 112.30
#> centred 1.01 1.42 1.41c(identical_fit = all.equal(fitted(fit_raw), fitted(fit_ctr)),
R2_linear = summary(lm(Weight ~ Height + Age, data = mlb3))$r.squared,
R2_quadratic = summary(fit_ctr)$r.squared)#> identical_fit R2_linear R2_quadratic
#> 1.000000 0.320225 0.323683Identical fitted values, wildly different VIFs. The quadratic term buys almost nothing here, but the technique matters.
ggplot(mlb3, aes(Age, Weight)) +
geom_point(alpha = 0.15, size = 1) +
geom_smooth(aes(colour = "Linear"), method = "lm", formula = y ~ x, se = FALSE) +
geom_smooth(aes(colour = "Quadratic"), method = "lm",
formula = y ~ poly(x, 2), se = FALSE) +
geom_smooth(aes(colour = "LOESS"), method = "loess", formula = y ~ x, se = FALSE) +
scale_colour_manual(values = c(Linear = "firebrick", Quadratic = "steelblue",
LOESS = "darkgreen")) +
labs(title = "Is the age-weight relationship curved?",
subtitle = "poly(x, 2) generates ORTHOGONAL polynomials, sidestepping the collinearity entirely",
x = "Age (yr)", y = "Weight (lb)", colour = NULL) +
theme_dspa()
poly(x, 2) returns orthogonal polynomial contrasts, so
it is the cleanest solution of all, no centring required, VIF exactly
1.
21.2 Binary indicators
Thresholding a numeric predictor creates a dummy variable whose coefficient has a clean interpretation as a difference of conditional expectations:
\[\beta_{\text{age30}}=E[\,\text{Weight}\mid \text{age30}=1,\ \text{others fixed}\,]-E[\,\text{Weight}\mid \text{age30}=0,\ \text{others fixed}\,].\]
mlb3$age30 <- as.integer(mlb3$Age >= 30)
fit_ind <- lm(Weight ~ Height + Age + age30, data = mlb3)
round(coef(summary(fit_ind)), 4)#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) -186.6102 18.4297 -10.1255 0.0000
#> Height 4.9645 0.2345 21.1747 0.0000
#> Age 0.7544 0.2175 3.4678 0.0005
#> age30 2.3251 1.9896 1.1686 0.2428Dichotomizing a continuous predictor discards information and reduces power: it replaces a graded relationship with a step. Do it when the threshold has substantive meaning (a clinical cutoff, a policy eligibility rule), not to simplify a plot.
21.3 Interactions
An interaction lets one predictor’s effect depend on another:
\[Y=\beta_0+\beta_1x_1+\beta_2x_2+\beta_3x_1x_2+\varepsilon \quad\Longrightarrow\quad \frac{\partial E[Y]}{\partial x_1}=\beta_1+\beta_3x_2 .\]
fit_int <- lm(Weight ~ Team + Height + Age * Position, data = mlb3)
anova(lm(Weight ~ Team + Height + Age + Position, data = mlb3), fit_int)c(R2_without = summary(lm(Weight ~ Team + Height + Age + Position, data = mlb3))$r.squared,
R2_with = summary(fit_int)$r.squared,
adjR2_without = summary(lm(Weight ~ Team + Height + Age + Position, data = mlb3))$adj.r.squared,
adjR2_with = summary(fit_int)$adj.r.squared)#> R2_without R2_with adjR2_without adjR2_with
#> 0.385804 0.392186 0.361706 0.363213The nested \(F\)-test is the right tool: it asks whether the whole block of interaction terms improves the fit more than chance, rather than scanning individual \(p\)-values.
Hierarchy principle. If an interaction \(x_1x_2\) is in the model, keep both main effects even if they are individually “insignificant”. Removing them makes the model non-invariant to shifts of origin, rescaling \(x_1\to x_1+c\) would change the fit, which is not a property any sensible model should have.
22 Model selection, and its limits
Stepwise procedures explore the model space greedily using an information criterion:
\[\mathrm{AIC}=-2\hat\ell+2p,\qquad \mathrm{BIC}=-2\hat\ell+p\log n .\]
In R, step(fit, k = 2) gives AIC and
step(fit, k = log(n)) gives BIC. Direction matters, and
forward selection needs a scope, starting
from a full model with direction = "forward" leaves nothing
to add.
full <- lm(Weight ~ Height + Age + Position, data = mlb3)
null <- lm(Weight ~ 1, data = mlb3)
back <- stats::step(full, direction = "backward", trace = 0)
fwd <- stats::step(null, scope = list(lower = ~1, upper = formula(full)),
direction = "forward", trace = 0)
both <- stats::step(null, scope = list(lower = ~1, upper = formula(full)),
direction = "both", trace = 0)
bic <- stats::step(full, direction = "backward", k = log(nobs(full)), trace = 0)data.frame(
search = c("backward (AIC)", "forward (AIC)", "both (AIC)", "backward (BIC)"),
terms = sapply(list(back, fwd, both, bic),
\(m) paste(attr(terms(m), "term.labels"), collapse = " + ")),
AIC = round(sapply(list(back, fwd, both, bic), AIC), 1),
adj_R2 = round(sapply(list(back, fwd, both, bic), \(m) summary(m)$adj.r.squared), 4)
)Common misconception: “stepwise selection finds the best model.” It does not, in three distinct senses.
It is greedy. Exhaustive search over \(p\) predictors requires \(2^p\) model fits, \(2^{20}\approx 10^6\), \(2^{50}\approx 10^{15}\). Stepwise examines \(O(p^2)\) of them and can never recover a pair of predictors that is jointly informative but individually weak.
The reported inference is invalid. \(p\)-values, confidence intervals, and \(R^2\) from the selected model ignore the search itself. They are computed as if the model had been specified in advance, so they are systematically optimistic, sometimes dramatically.
It is unstable. Small perturbations of the data change which variables survive.
set.seed(89)
sel_counts <- table(unlist(replicate(150, {
b <- mlb3[sample(nrow(mlb3), replace = TRUE), ]
m <- stats::step(lm(Weight ~ Height + Age + Position + age30, data = b),
direction = "backward", trace = 0)
attr(terms(m), "term.labels")
}, simplify = FALSE)))
data.frame(term = names(sel_counts),
selected_in = as.integer(sel_counts),
pct_of_150 = round(100 * as.integer(sel_counts) / 150, 1)) |>
arrange(desc(selected_in))Terms selected in only a fraction of bootstrap resamples are not reliably in the model, and that fraction, stability selection, is far more informative than the single answer stepwise returns. Regularized selection (Chapter 11) addresses all three problems more directly, and post-selection inference methods exist for the second.
23 PART III: TREE-BASED REGRESSION
Linear models impose a global functional form. Trees impose none: they recursively partition the predictor space into rectangles and fit a simple model inside each. The price of that flexibility is instability and a loss of smoothness; the payoff is automatic detection of interactions and nonlinearities.
24 Regression trees
24.1 The splitting criterion
A regression tree partitions \(\mathbb{R}^{p}\) into disjoint regions \(R_1,\dots,R_M\) and predicts a constant in each:
\[\hat f(\mathbf{x})=\sum_{m=1}^{M}c_m\,\mathbb{1}\{\mathbf{x}\in R_m\}.\]
For squared-error loss the optimal constant is the region mean, \(\hat c_m=\bar y_{R_m}\), so growing the tree means choosing splits that minimize within-region sum of squares. Finding the globally optimal partition is NP-hard, so CART proceeds greedily: at each node, choose the variable \(j\) and cut point \(s\) maximizing the variance reduction
\[\Delta(j,s)=\underbrace{\sum_{i\in T}(y_i-\bar y_T)^2}_{\mathrm{SSE}_T} -\Big(\underbrace{\sum_{i\in T_L}(y_i-\bar y_{T_L})^2}_{\mathrm{SSE}_L} +\underbrace{\sum_{i\in T_R}(y_i-\bar y_{T_R})^2}_{\mathrm{SSE}_R}\Big),\]
which is always \(\ge 0\), splitting can never increase within-node SSE. Normalizing by \(n_T\) gives the equivalent statement in terms of variance:
\[\Delta_{\text{var}}=\operatorname{Var}(T)-\frac{n_L}{n_T}\operatorname{Var}(T_L)-\frac{n_R}{n_T}\operatorname{Var}(T_R).\]
A closely related heuristic, the standard deviation reduction, replaces variances with standard deviations:
\[\mathrm{SDR}=\operatorname{sd}(T)-\sum_{m}\frac{n_m}{n_T}\operatorname{sd}(T_m).\]
SDR is the criterion used by M5 and is intuitive to read, but it is not the same objective as SSE minimization, \(\operatorname{sd}\) is a concave function of variance, so the two can rank splits differently. CART minimizes SSE directly.
orig <- c(1, 2, 3, 3, 4, 5, 6, 6, 7, 8)
split_stats <- function(L, R, full) {
n <- length(full)
sse <- function(v) sum((v - mean(v))^2)
c(SSE_reduction = sse(full) - (sse(L) + sse(R)),
SDR = sd(full) - (length(L) / n * sd(L) + length(R) / n * sd(R)),
pred_left = mean(L), pred_right = mean(R))
}
rbind(
`Split A: {1,2,3} | {3,4,5,6,6,7,8}` =
split_stats(c(1, 2, 3), c(3, 4, 5, 6, 6, 7, 8), orig),
`Split B: {1,2,3,3,4,5} | {6,6,7,8}` =
split_stats(c(1, 2, 3, 3, 4, 5), c(6, 6, 7, 8), orig)
) |> round(4)#> SSE_reduction SDR pred_left pred_right
#> Split A: {1,2,3} | {3,4,5,6,6,7,8} 26.7857 0.7703 2 5.5714
#> Split B: {1,2,3,3,4,5} | {6,6,7,8} 33.7500 1.0415 3 6.7500Both criteria prefer Split B, which produces more homogeneous children. Terminal-node predictions are the region means, 3.0 and 6.75.
24.2 Cost and control
At each node the greedy search evaluates every variable at every candidate cut point. Sorting each variable once costs \(O(n\log n)\), after which the SSE for all \(n-1\) cuts is obtained in a single pass using running sums. Per node the cost is \(O(pn)\) after sorting, so building a balanced tree of depth \(O(\log n)\) costs
\[O\big(pn\log n\big).\]
Trees left to grow will fit the training data perfectly and
generalize badly. Two controls: pre-pruning via
minsplit, minbucket, and
maxdepth; and post-pruning via
cost-complexity, which minimizes
\[C_\alpha(T)=\sum_{m=1}^{|T|}\mathrm{SSE}_m+\alpha|T|\]
over subtrees, with \(\alpha\)
chosen by cross-validation. rpart computes the whole
cost-complexity path and reports the cross-validated error for each,
which is what printcp() shows.
library(rpart)
library(rpart.plot)
set.seed(1234)
idx <- sample(seq_len(nrow(mlb2)), size = floor(0.75 * nrow(mlb2)))
mlb_train <- mlb2[idx, ]
mlb_test <- mlb2[-idx, ]
mlb_tree <- rpart(Weight ~ Height + Age, data = mlb_train)
mlb_tree#> n= 775
#>
#> node), split, n, deviance, yval
#> * denotes terminal node
#>
#> 1) root 775 341952.00 201.858
#> 2) Height< 74.5 500 164839.00 195.120
#> 4) Height< 72.5 239 73317.90 189.757
#> 8) Height< 70.5 63 15833.70 182.063 *
#> 9) Height>=70.5 176 52420.00 192.511 *
#> 5) Height>=72.5 261 78353.80 200.031
#> 10) Age< 27.28 106 30455.10 194.274 *
#> 11) Age>=27.28 155 41982.80 203.968 *
#> 3) Height>=74.5 275 113139.00 214.109
#> 6) Age< 30.015 201 86834.90 211.159
#> 12) Height< 75.5 91 35635.80 205.846 *
#> 13) Height>=75.5 110 46505.20 215.555
#> 26) Age< 24.795 23 7365.22 204.348 *
#> 27) Age>=24.795 87 35487.70 218.517 *
#> 7) Age>=30.015 74 19803.90 222.122 *#>
#> Regression tree:
#> rpart(formula = Weight ~ Height + Age, data = mlb_train)
#>
#> Variables actually used in tree construction:
#> [1] Age Height
#>
#> Root node error: 341952/775 = 441.2
#>
#> n= 775
#>
#> CP nsplit rel error xerror xstd
#> 1 0.18709 0 1.0000 1.0028 0.05526
#> 2 0.03851 1 0.8129 0.8170 0.04465
#> 3 0.01901 2 0.7744 0.7964 0.04469
#> 4 0.01730 3 0.7554 0.8161 0.04739
#> 5 0.01481 4 0.7381 0.8079 0.04798
#> 6 0.01373 5 0.7233 0.7999 0.04752
#> 7 0.01068 6 0.7096 0.7887 0.04727
#> 8 0.01000 7 0.6989 0.7789 0.04738cptab <- as.data.frame(mlb_tree$cptable)
ggplot(cptab, aes(nsplit, xerror)) +
geom_line(linewidth = 0.9, colour = "steelblue") +
geom_point(size = 2) +
geom_errorbar(aes(ymin = xerror - xstd, ymax = xerror + xstd), width = 0.15,
colour = "grey45") +
geom_hline(yintercept = min(cptab$xerror) + cptab$xstd[which.min(cptab$xerror)],
linetype = "dashed", colour = "firebrick") +
labs(title = "Cost-complexity pruning path",
subtitle = "Dashed line: the 1-SE rule threshold; choose the smallest tree beneath it",
x = "Number of splits", y = "Cross-validated relative error") +
theme_dspa()
rpart.plot(mlb_tree, digits = 4, fallen.leaves = TRUE, type = 3, extra = 101,
box.palette = "BuGn",
main = "Regression tree: MLB weight from height and age")
Because a tree with two numeric predictors is a piecewise-constant surface, the partition is best appreciated in three dimensions:
gh <- seq(min(mlb2$Height), max(mlb2$Height), length.out = 60)
ga <- seq(min(mlb2$Age), max(mlb2$Age), length.out = 60)
gridp <- expand.grid(Height = gh, Age = ga)
zt <- matrix(predict(mlb_tree, gridp), nrow = 60, ncol = 60) # rows = Height
plot_ly() |>
add_surface(x = ga, y = gh, z = zt, opacity = 0.85, showscale = FALSE,
colorscale = "Viridis") |>
add_markers(x = mlb2$Age, y = mlb2$Height, z = mlb2$Weight,
type = "scatter3d", mode = "markers", name = "Players",
marker = list(size = 1.6, opacity = 0.3, color = "black")) |>
layout(title = "A regression tree is a piecewise-constant surface",
scene = list(xaxis = list(title = "Age (yr)"),
yaxis = list(title = "Height (in)"),
zaxis = list(title = "Predicted weight (lb)")))The flat plateaus and vertical cliffs are the tree’s defining feature. They make it robust to monotone transformations and able to represent interactions automatically, and they make it a poor model for a genuinely smooth relationship, which is exactly what height-to-weight is.
24.3 Evaluating a regression tree
\[\mathrm{MAE}=\frac{1}{n}\sum_{i=1}^{n}\big|\hat y_i-y_i\big|, \qquad \mathrm{RMSE}=\sqrt{\frac{1}{n}\sum_{i=1}^{n}(\hat y_i-y_i)^2}.\]
MAE is in the units of the response and is robust; RMSE penalizes large errors more heavily and is the quantity least squares minimizes. Report both, and always compare against a baseline.
MAE <- function(obs, pred) mean(abs(obs - pred))
RMSE <- function(obs, pred) sqrt(mean((obs - pred)^2))
pred_tree <- predict(mlb_tree, mlb_test)
pred_lin <- predict(lm(Weight ~ Height + Age, data = mlb_train), mlb_test)
pred_base <- rep(mean(mlb_train$Weight), nrow(mlb_test))
data.frame(
model = c("Baseline (training mean)", "Linear model", "Regression tree"),
MAE = c(MAE(mlb_test$Weight, pred_base), MAE(mlb_test$Weight, pred_lin),
MAE(mlb_test$Weight, pred_tree)),
RMSE = c(RMSE(mlb_test$Weight, pred_base), RMSE(mlb_test$Weight, pred_lin),
RMSE(mlb_test$Weight, pred_tree)),
cor = c(NA, cor(pred_lin, mlb_test$Weight), cor(pred_tree, mlb_test$Weight))
) |> mutate(across(where(is.numeric), \(z) round(z, 4)))cmp <- bind_rows(
data.frame(obs = mlb_test$Weight, pred = pred_lin, model = "Linear model"),
data.frame(obs = mlb_test$Weight, pred = pred_tree, model = "Regression tree"))
ggplot(cmp, aes(pred, obs)) +
geom_point(alpha = 0.35, size = 1.3, colour = "steelblue") +
geom_abline(slope = 1, intercept = 0, colour = "firebrick", linetype = "dashed") +
facet_wrap(~ model) +
labs(title = "Predicted versus observed weight on held-out data",
subtitle = "The tree's vertical stripes are its finitely many terminal-node predictions",
x = "Predicted (lb)", y = "Observed (lb)") +
theme_dspa(10)
The tree can emit only as many distinct predictions as it has leaves, which is why its panel shows discrete vertical stripes. It cannot extrapolate beyond the range of the training response at all, the maximum prediction is the largest leaf mean, never the largest observed value.
25 Model trees
A model tree replaces each terminal node’s constant with a linear model fitted to the cases in that node. The result is a piecewise-linear surface: it inherits the tree’s ability to find regions while recovering smoothness and extrapolation within them.
Quinlan’s M5 introduced the idea; Cubist is its direct
successor by the same author, implemented in C with no Java dependency,
and it adds rule-based simplification and optional nearest-neighbour
correction of predictions.
library(Cubist)
cub <- cubist(x = mlb_train[, c("Height", "Age")], y = mlb_train$Weight)
summary(cub)$output |> strsplit("\n") |> _[[1]] |> head(28) |> cat(sep = "\n")#>
#> Cubist [Release 2.07 GPL Edition] Tue Aug 25 16:23:49 2026
#> ---------------------------------
#>
#> Target attribute `outcome'
#>
#> Read 775 cases (3 attributes) from undefined.data
#>
#> Model:
#>
#> Rule 1: [775 cases, mean 201.9, range 150 to 290, est err 13.9]
#>
#> outcome = -179.6 + 4.8 Height + 0.95 Age
#>
#>
#> Evaluation on training data (775 cases):
#>
#> Average |error| 15.0
#> Relative |error| 0.90
#> Correlation coefficient 0.46
#>
#>
#> Attribute usage:
#> Conds Model
#>
#> 100% Height
#> 100% Agepred_cub <- predict(cub, mlb_test[, c("Height", "Age")])
data.frame(
model = c("Baseline", "Linear model", "Regression tree", "Model tree (Cubist)"),
MAE = c(MAE(mlb_test$Weight, pred_base), MAE(mlb_test$Weight, pred_lin),
MAE(mlb_test$Weight, pred_tree), MAE(mlb_test$Weight, pred_cub)),
RMSE = c(RMSE(mlb_test$Weight, pred_base), RMSE(mlb_test$Weight, pred_lin),
RMSE(mlb_test$Weight, pred_tree), RMSE(mlb_test$Weight, pred_cub)),
cor = c(NA, cor(pred_lin, mlb_test$Weight), cor(pred_tree, mlb_test$Weight),
cor(pred_cub, mlb_test$Weight))
) |> mutate(across(where(is.numeric), \(z) round(z, 4)))# --- The original M5 implementation, for reference -------------------------
# RWeka wraps Weka's Java implementation. It requires a working rJava and a
# Java runtime; Cubist above is the same algorithm family without that
# dependency. Do NOT set WEKA_HOME unless you know your installation path --
# overwriting it will break a working install.
library(RWeka)
mlb_m5 <- M5P(Weight ~ Height + Age, data = mlb_train)
summary(mlb_m5)
pred_m5 <- predict(mlb_m5, mlb_test)
c(cor = cor(pred_m5, mlb_test$Weight), MAE = MAE(mlb_test$Weight, pred_m5))26 Bayesian additive regression trees
26.1 The model
BART (Chipman, George & McCulloch, 2010) is a sum-of-trees model with a regularizing prior:
\[y_i=\sum_{j=1}^{m}g\big(\mathbf{x}_i\mid T_j,M_j\big)+\varepsilon_i, \qquad \varepsilon_i\stackrel{iid}{\sim}N(0,\sigma^2),\]
where \(T_j\) is the \(j\)-th tree structure and \(M_j=\{\mu_{1j},\dots,\mu_{b_jj}\}\) its terminal-node values. Typically \(m=200\) trees, each deliberately weak.
The distinction from boosting matters. Boosting fits trees sequentially to residuals and has no probability model. BART places a prior over the whole ensemble and samples from the posterior, so uncertainty quantification comes free and correctly, every prediction arrives as a full posterior distribution rather than a point.
The prior factorizes as
\[p\big(\sigma,\{(T_j,M_j)\}_{j=1}^m\big)=p(\sigma)\prod_{j=1}^{m}\left[p(T_j)\prod_{k}p\big(\mu_{kj}\mid T_j\big)\right],\]
with four regularizing components:
- Tree depth. A node at depth \(d\) splits with probability \(\alpha(1+d)^{-\beta}\); defaults \(\alpha=0.95,\beta=2\) give split probabilities \(0.95, 0.24, 0.11, 0.06,\dots\), so trees are strongly discouraged from growing deep. This is the main regularizer, it keeps each tree weak so that the ensemble, not any member, does the work.
- Splitting variable. Uniform over available predictors.
- Split point. Uniform over the observed values of the chosen predictor.
- Leaf values and noise. \(\mu_{kj}\sim N(0,\sigma_\mu^2)\) with \(\sigma_\mu=0.5/(\kappa\sqrt m)\) after rescaling \(y\) to \([-0.5,0.5]\), so each tree contributes only \(O(1/\sqrt m)\); and \(\sigma^2\sim\nu\lambda/\chi^2_\nu\), a conjugate scaled-inverse-\(\chi^2\) calibrated from the data.
26.2 Fitting by MCMC
BART is fitted by Bayesian backfitting Gibbs sampling. Holding all other trees fixed, define the partial residual
\[R_{j}=\mathbf{y}-\sum_{j'\ne j}g\big(\mathbf{x}\mid T_{j'},M_{j'}\big),\]
then draw \((T_j,M_j)\mid R_j,\sigma\) via a Metropolis–Hastings step over local tree moves (grow, prune, change, swap), and finally \(\sigma\mid\text{all}\) from its conjugate posterior. Repeat over \(j=1,\dots,m\) and iterate.
Cost is \(O\big(\text{iterations}\times m\times n\big)\), every sweep touches every tree and every observation. This is the main practical constraint.
Three steps summarize the algorithm: initialize the prior on \((f,\sigma)\); run the chain whose stationary distribution is the posterior \(p(f,\sigma\mid\text{data})\); and treat the retained draws as a sample from the full posterior, so that pointwise means give predictions and pointwise quantiles give credible intervals.
26.3 One-dimensional demonstration
library(BART)
f_true <- function(x) sin(x) * x^3
sig <- 0.2
set.seed(1234)
n_b <- 300
xb <- sort(2 * runif(n_b) - 1)
yb <- f_true(xb) + sig * rnorm(n_b)
xtest <- seq(-1, 1, by = 0.02)
bart1 <- wbart(x.train = data.frame(x = xb), y.train = yb,
x.test = data.frame(x = xtest),
nskip = 300, ndpost = 1000, printevery = 100000L)#> *****Into main of wbart
#> *****Data:
#> data:n,p,np: 300, 1, 101
#> y1,yn: 0.572918, 1.073126
#> x1,x[n*p]: -0.993435, 0.997482
#> xp1,xp[np*p]: -1.000000, 1.000000
#> *****Number of Trees: 200
#> *****Number of Cut Points: 100 ... 100
#> *****burn and ndpost: 300, 1000
#> *****Prior:beta,alpha,tau,nu,lambda: 2.000000,0.950000,0.032889,3.000000,0.018662
#> *****sigma: 0.309524
#> *****w (weights): 1.000000 ... 1.000000
#> *****Dirichlet:sparse,theta,omega,a,b,rho,augment: 0,0,1,0.5,1,1,0
#> *****nkeeptrain,nkeeptest,nkeeptestme,nkeeptreedraws: 1000,1000,1000,1000
#> *****printevery: 100000
#> *****skiptr,skipte,skipteme,skiptreedraws: 1,1,1,1
#>
#> MCMC
#> done 0 (out of 1300)
#> time: 2s
#> check counts
#> trcnt,tecnt,temecnt,treedrawscnt: 1000,1000,1000,1000#> [1] 1000 101qm <- apply(bart1$yhat.test, 2, quantile, probs = c(0.025, 0.975))
band <- data.frame(x = xtest, fit = colMeans(bart1$yhat.test),
lo = qm[1, ], hi = qm[2, ], truth = f_true(xtest))
ggplot() +
geom_point(data = data.frame(x = xb, y = yb), aes(x, y),
colour = "grey55", size = 1, alpha = 0.7) +
geom_ribbon(data = band, aes(x, ymin = lo, ymax = hi, fill = "95% credible band"),
alpha = 0.30) +
geom_line(data = band, aes(x, truth, colour = "True signal"),
linewidth = 1, linetype = "dashed") +
geom_line(data = band, aes(x, fit, colour = "Posterior mean"), linewidth = 1) +
scale_fill_manual(values = c("95% credible band" = "steelblue")) +
scale_colour_manual(values = c("True signal" = "black",
"Posterior mean" = "firebrick")) +
labs(title = expression(paste("BART recovers ", f(x) == sin(x) %.% x^3, " with n = 300")),
x = "x", y = "y", colour = NULL, fill = NULL) +
theme_dspa()
# --- Interactive equivalent ------------------------------------------------
plot_ly(x = xb, y = yb, type = "scatter", mode = "markers", name = "Data") |>
add_trace(x = xtest, y = f_true(xtest), mode = "lines", name = "True signal") |>
add_trace(x = xtest, y = colMeans(bart1$yhat.test), mode = "lines",
name = "Posterior mean") |>
add_trace(x = xtest, y = qm[1, ], mode = "lines", name = "Lower 2.5%") |>
add_trace(x = xtest, y = qm[2, ], mode = "lines", name = "Upper 97.5%") |>
layout(title = "BART Model (n = 300)",
xaxis = list(title = "X"), yaxis = list(title = "Y"),
legend = list(orientation = "h"))More data should tighten the bands. It does, at the predicted rate:
set.seed(1234)
n2 <- 3000
x2 <- sort(2 * runif(n2) - 1)
y2 <- f_true(x2) + sig * rnorm(n2)
bart2 <- wbart(x.train = data.frame(x = x2), y.train = y2,
x.test = data.frame(x = xtest),
nskip = 300, ndpost = 1000, printevery = 100000L)#> *****Into main of wbart
#> *****Data:
#> data:n,p,np: 3000, 1, 101
#> y1,yn: 1.021931, 0.626872
#> x1,x[n*p]: -0.999316, 0.998606
#> xp1,xp[np*p]: -1.000000, 1.000000
#> *****Number of Trees: 200
#> *****Number of Cut Points: 100 ... 100
#> *****burn and ndpost: 300, 1000
#> *****Prior:beta,alpha,tau,nu,lambda: 2.000000,0.950000,0.035536,3.000000,0.017828
#> *****sigma: 0.302529
#> *****w (weights): 1.000000 ... 1.000000
#> *****Dirichlet:sparse,theta,omega,a,b,rho,augment: 0,0,1,0.5,1,1,0
#> *****nkeeptrain,nkeeptest,nkeeptestme,nkeeptreedraws: 1000,1000,1000,1000
#> *****printevery: 100000
#> *****skiptr,skipte,skipteme,skiptreedraws: 1,1,1,1
#>
#> MCMC
#> done 0 (out of 1300)
#> time: 8s
#> check counts
#> trcnt,tecnt,temecnt,treedrawscnt: 1000,1000,1000,1000q2 <- apply(bart2$yhat.test, 2, quantile, probs = c(0.025, 0.975))
c(mean_band_width_n300 = mean(qm[2, ] - qm[1, ]),
mean_band_width_n3000 = mean(q2[2, ] - q2[1, ]),
ratio = mean(qm[2, ] - qm[1, ]) / mean(q2[2, ] - q2[1, ]),
sqrt_10 = sqrt(10))#> mean_band_width_n300 mean_band_width_n3000 ratio
#> 0.255149 0.118263 2.157468
#> sqrt_10
#> 3.162278The band width shrinks by a factor close to \(\sqrt{10}\) when \(n\) increases tenfold, the \(O(n^{-1/2})\) rate expected of a well-behaved posterior.
bind_rows(
data.frame(x = xtest, lo = qm[1, ], hi = qm[2, ], n = "n = 300"),
data.frame(x = xtest, lo = q2[1, ], hi = q2[2, ], n = "n = 3,000")) |>
ggplot(aes(x)) +
geom_ribbon(aes(ymin = lo, ymax = hi, fill = n), alpha = 0.4) +
geom_line(data = data.frame(x = xtest, y = f_true(xtest)), aes(x, y),
linewidth = 0.9, colour = "black") +
scale_fill_manual(values = c("n = 300" = "firebrick", "n = 3,000" = "steelblue")) +
labs(title = "Posterior credible bands narrow as the sample grows",
x = "x", y = "y", fill = NULL) +
theme_dspa()
26.4 Higher-dimensional simulation
set.seed(1234)
n_h2 <- 2000; p_h2 <- 20
beta_h <- 3 * (1:p_h2) / p_h2
Xh2 <- matrix(rnorm(n_h2 * p_h2), ncol = p_h2)
yh2 <- as.double(10 + Xh2 %*% beta_h + rnorm(n_h2))
n_test <- 5000
Xp <- matrix(rnorm(n_test * p_h2), ncol = p_h2)
y_test_true <- as.double(10 + Xp %*% beta_h)
t_bart <- system.time(
bart_md <- wbart(x.train = as.data.frame(Xh2), y.train = yh2,
x.test = as.data.frame(Xp),
nskip = 200, ndpost = 400, printevery = 100000L)
)[["elapsed"]]#> *****Into main of wbart
#> *****Data:
#> data:n,p,np: 2000, 20, 5000
#> y1,yn: 6.218160, 6.111044
#> x1,x[n*p]: -1.207066, -0.943031
#> xp1,xp[np*p]: -0.561920, 0.628735
#> *****Number of Trees: 200
#> *****Number of Cut Points: 100 ... 100
#> *****burn and ndpost: 200, 400
#> *****Prior:beta,alpha,tau,nu,lambda: 2.000000,0.950000,0.994832,3.000000,0.209999
#> *****sigma: 1.038301
#> *****w (weights): 1.000000 ... 1.000000
#> *****Dirichlet:sparse,theta,omega,a,b,rho,augment: 0,0,1,0.5,1,20,0
#> *****nkeeptrain,nkeeptest,nkeeptestme,nkeeptreedraws: 400,400,400,400
#> *****printevery: 100000
#> *****skiptr,skipte,skipteme,skiptreedraws: 1,1,1,1
#>
#> MCMC
#> done 0 (out of 600)
#> time: 8s
#> check counts
#> trcnt,tecnt,temecnt,treedrawscnt: 400,400,400,400lm_md <- lm(y ~ ., data = data.frame(Xh2, y = yh2))
pred_lm <- predict(lm_md, data.frame(Xp))
c(bart_seconds = round(t_bart, 2),
RMSE_BART = RMSE(y_test_true, bart_md$yhat.test.mean),
RMSE_LM = RMSE(y_test_true, pred_lm))#> bart_seconds RMSE_BART RMSE_LM
#> 7.5200000 1.0797239 0.0950987ggplot(data.frame(lm = pred_lm, bart = bart_md$yhat.test.mean),
aes(lm, bart)) +
geom_point(alpha = 0.12, size = 0.7, colour = "steelblue") +
geom_abline(slope = 1, intercept = 0, colour = "firebrick", linewidth = 1) +
labs(title = "BART versus the linear model on a truly linear process",
subtitle = "The generating mechanism IS linear, so OLS is optimal and BART pays for its flexibility",
x = "Linear model prediction", y = "BART posterior mean") +
theme_dspa()
This is an instructive negative result. When the data-generating process is exactly linear, OLS is the MLE and BART’s extra flexibility buys nothing, it costs both accuracy and several orders of magnitude in compute. Model flexibility is not free. Reach for BART when you suspect interactions or nonlinearity you cannot specify, not by default.
26.5 Heart-attack charges: BART versus LASSO
ha_model <- ha |>
mutate(gender = as.integer(SEX == "F")) |>
dplyr::select(-any_of(c("Patient", "DIAGNOSIS", "SEX")))
names(ha_model)#> [1] "DRG" "DIED" "CHARGES" "LOS" "AGE" "gender"x_all <- as.matrix(ha_model[, setdiff(names(ha_model), "CHARGES")])
y_all <- ha_model$CHARGES
c(n = nrow(x_all), p = ncol(x_all))#> n p
#> 148 5Standardization must use training statistics only. Scaling the test set by its own mean and standard deviation would place the fitted model and the new data in different coordinate systems, silently corrupting every prediction.
library(glmnet)
nd <- 10
n_a <- length(y_all)
ntrain <- floor(0.8 * n_a)
res <- data.frame(split = 1:nd, RMSE_BART = NA_real_, RMSE_LASSO = NA_real_)
pred_store <- list()
for (i in seq_len(nd)) {
set.seed(1234 * i)
tr <- sample(n_a, ntrain)
y_tr <- y_all[tr]; y_te <- y_all[-tr]
x_tr <- x_all[tr, , drop = FALSE]; x_te <- x_all[-tr, , drop = FALSE]
# Standardize with TRAINING statistics, then apply to test
ctr <- colMeans(x_tr)
scl <- apply(x_tr, 2, sd); scl[scl == 0] <- 1
x_tr_s <- scale(x_tr, center = ctr, scale = scl)
x_te_s <- scale(x_te, center = ctr, scale = scl)
m_bart <- wbart(x_tr, y_tr, x_te, nskip = 200, ndpost = 400,
printevery = 100000L)
cvl <- cv.glmnet(x_tr_s, y_tr, family = "gaussian")
m_lasso <- glmnet(x_tr_s, y_tr, family = "gaussian", lambda = cvl$lambda.min)
p_b <- m_bart$yhat.test.mean
p_l <- as.vector(predict(m_lasso, x_te_s, s = cvl$lambda.min))
res$RMSE_BART[i] <- RMSE(y_te, p_b)
res$RMSE_LASSO[i] <- RMSE(y_te, p_l)
pred_store[[i]] <- data.frame(bart = p_b, lasso = p_l, obs = y_te)
}#> *****Into main of wbart
#> *****Data:
#> data:n,p,np: 118, 5, 30
#> y1,yn: 2105.135593, -2641.864407
#> x1,x[n*p]: 121.000000, 0.000000
#> xp1,xp[np*p]: 122.000000, 1.000000
#> *****Number of Trees: 200
#> *****Number of Cut Points: 2 ... 1
#> *****burn and ndpost: 200, 400
#> *****Prior:beta,alpha,tau,nu,lambda: 2.000000,0.950000,290.550176,3.000000,2228477.133405
#> *****sigma: 3382.354604
#> *****w (weights): 1.000000 ... 1.000000
#> *****Dirichlet:sparse,theta,omega,a,b,rho,augment: 0,0,1,0.5,1,5,0
#> *****nkeeptrain,nkeeptest,nkeeptestme,nkeeptreedraws: 400,400,400,400
#> *****printevery: 100000
#> *****skiptr,skipte,skipteme,skiptreedraws: 1,1,1,1
#>
#> MCMC
#> done 0 (out of 600)
#> time: 1s
#> check counts
#> trcnt,tecnt,temecnt,treedrawscnt: 400,400,400,400
#> *****Into main of wbart
#> *****Data:
#> data:n,p,np: 118, 5, 30
#> y1,yn: -3870.449153, 5404.550847
#> x1,x[n*p]: 122.000000, 1.000000
#> xp1,xp[np*p]: 122.000000, 0.000000
#> *****Number of Trees: 200
#> *****Number of Cut Points: 2 ... 1
#> *****burn and ndpost: 200, 400
#> *****Prior:beta,alpha,tau,nu,lambda: 2.000000,0.950000,290.550176,3.000000,2261599.844262
#> *****sigma: 3407.398506
#> *****w (weights): 1.000000 ... 1.000000
#> *****Dirichlet:sparse,theta,omega,a,b,rho,augment: 0,0,1,0.5,1,5,0
#> *****nkeeptrain,nkeeptest,nkeeptestme,nkeeptreedraws: 400,400,400,400
#> *****printevery: 100000
#> *****skiptr,skipte,skipteme,skiptreedraws: 1,1,1,1
#>
#> MCMC
#> done 0 (out of 600)
#> time: 1s
#> check counts
#> trcnt,tecnt,temecnt,treedrawscnt: 400,400,400,400
#> *****Into main of wbart
#> *****Data:
#> data:n,p,np: 118, 5, 30
#> y1,yn: 2659.855932, -1749.144068
#> x1,x[n*p]: 121.000000, 1.000000
#> xp1,xp[np*p]: 122.000000, 0.000000
#> *****Number of Trees: 200
#> *****Number of Cut Points: 2 ... 1
#> *****burn and ndpost: 200, 400
#> *****Prior:beta,alpha,tau,nu,lambda: 2.000000,0.950000,290.550176,3.000000,2182151.820220
#> *****sigma: 3347.013986
#> *****w (weights): 1.000000 ... 1.000000
#> *****Dirichlet:sparse,theta,omega,a,b,rho,augment: 0,0,1,0.5,1,5,0
#> *****nkeeptrain,nkeeptest,nkeeptestme,nkeeptreedraws: 400,400,400,400
#> *****printevery: 100000
#> *****skiptr,skipte,skipteme,skiptreedraws: 1,1,1,1
#>
#> MCMC
#> done 0 (out of 600)
#> time: 0s
#> check counts
#> trcnt,tecnt,temecnt,treedrawscnt: 400,400,400,400
#> *****Into main of wbart
#> *****Data:
#> data:n,p,np: 118, 5, 30
#> y1,yn: 714.000000, -4103.000000
#> x1,x[n*p]: 121.000000, 1.000000
#> xp1,xp[np*p]: 122.000000, 1.000000
#> *****Number of Trees: 200
#> *****Number of Cut Points: 2 ... 1
#> *****burn and ndpost: 200, 400
#> *****Prior:beta,alpha,tau,nu,lambda: 2.000000,0.950000,290.550176,3.000000,2310883.132530
#> *****sigma: 3444.324312
#> *****w (weights): 1.000000 ... 1.000000
#> *****Dirichlet:sparse,theta,omega,a,b,rho,augment: 0,0,1,0.5,1,5,0
#> *****nkeeptrain,nkeeptest,nkeeptestme,nkeeptreedraws: 400,400,400,400
#> *****printevery: 100000
#> *****skiptr,skipte,skipteme,skiptreedraws: 1,1,1,1
#>
#> MCMC
#> done 0 (out of 600)
#> time: 1s
#> check counts
#> trcnt,tecnt,temecnt,treedrawscnt: 400,400,400,400
#> *****Into main of wbart
#> *****Data:
#> data:n,p,np: 118, 5, 30
#> y1,yn: 3229.305085, 1884.305085
#> x1,x[n*p]: 122.000000, 1.000000
#> xp1,xp[np*p]: 122.000000, 0.000000
#> *****Number of Trees: 200
#> *****Number of Cut Points: 2 ... 1
#> *****burn and ndpost: 200, 400
#> *****Prior:beta,alpha,tau,nu,lambda: 2.000000,0.950000,290.550176,3.000000,2169347.152131
#> *****sigma: 3337.179552
#> *****w (weights): 1.000000 ... 1.000000
#> *****Dirichlet:sparse,theta,omega,a,b,rho,augment: 0,0,1,0.5,1,5,0
#> *****nkeeptrain,nkeeptest,nkeeptestme,nkeeptreedraws: 400,400,400,400
#> *****printevery: 100000
#> *****skiptr,skipte,skipteme,skiptreedraws: 1,1,1,1
#>
#> MCMC
#> done 0 (out of 600)
#> time: 1s
#> check counts
#> trcnt,tecnt,temecnt,treedrawscnt: 400,400,400,400
#> *****Into main of wbart
#> *****Data:
#> data:n,p,np: 118, 5, 30
#> y1,yn: 2115.474576, 1998.474576
#> x1,x[n*p]: 121.000000, 1.000000
#> xp1,xp[np*p]: 122.000000, 1.000000
#> *****Number of Trees: 200
#> *****Number of Cut Points: 2 ... 1
#> *****burn and ndpost: 200, 400
#> *****Prior:beta,alpha,tau,nu,lambda: 2.000000,0.950000,290.550176,3.000000,2209447.645085
#> *****sigma: 3367.882283
#> *****w (weights): 1.000000 ... 1.000000
#> *****Dirichlet:sparse,theta,omega,a,b,rho,augment: 0,0,1,0.5,1,5,0
#> *****nkeeptrain,nkeeptest,nkeeptestme,nkeeptreedraws: 400,400,400,400
#> *****printevery: 100000
#> *****skiptr,skipte,skipteme,skiptreedraws: 1,1,1,1
#>
#> MCMC
#> done 0 (out of 600)
#> time: 0s
#> check counts
#> trcnt,tecnt,temecnt,treedrawscnt: 400,400,400,400
#> *****Into main of wbart
#> *****Data:
#> data:n,p,np: 118, 5, 30
#> y1,yn: -1255.525424, -3624.525424
#> x1,x[n*p]: 123.000000, 0.000000
#> xp1,xp[np*p]: 122.000000, 0.000000
#> *****Number of Trees: 200
#> *****Number of Cut Points: 2 ... 1
#> *****burn and ndpost: 200, 400
#> *****Prior:beta,alpha,tau,nu,lambda: 2.000000,0.950000,290.550176,3.000000,2177441.668674
#> *****sigma: 3343.399788
#> *****w (weights): 1.000000 ... 1.000000
#> *****Dirichlet:sparse,theta,omega,a,b,rho,augment: 0,0,1,0.5,1,5,0
#> *****nkeeptrain,nkeeptest,nkeeptestme,nkeeptreedraws: 400,400,400,400
#> *****printevery: 100000
#> *****skiptr,skipte,skipteme,skiptreedraws: 1,1,1,1
#>
#> MCMC
#> done 0 (out of 600)
#> time: 1s
#> check counts
#> trcnt,tecnt,temecnt,treedrawscnt: 400,400,400,400
#> *****Into main of wbart
#> *****Data:
#> data:n,p,np: 118, 5, 30
#> y1,yn: 3795.457627, -388.542373
#> x1,x[n*p]: 121.000000, 0.000000
#> xp1,xp[np*p]: 122.000000, 1.000000
#> *****Number of Trees: 200
#> *****Number of Cut Points: 2 ... 1
#> *****burn and ndpost: 200, 400
#> *****Prior:beta,alpha,tau,nu,lambda: 2.000000,0.950000,288.534922,3.000000,1972180.926129
#> *****sigma: 3181.913893
#> *****w (weights): 1.000000 ... 1.000000
#> *****Dirichlet:sparse,theta,omega,a,b,rho,augment: 0,0,1,0.5,1,5,0
#> *****nkeeptrain,nkeeptest,nkeeptestme,nkeeptreedraws: 400,400,400,400
#> *****printevery: 100000
#> *****skiptr,skipte,skipteme,skiptreedraws: 1,1,1,1
#>
#> MCMC
#> done 0 (out of 600)
#> time: 1s
#> check counts
#> trcnt,tecnt,temecnt,treedrawscnt: 400,400,400,400
#> *****Into main of wbart
#> *****Data:
#> data:n,p,np: 118, 5, 30
#> y1,yn: -3943.779661, -524.779661
#> x1,x[n*p]: 123.000000, 0.000000
#> xp1,xp[np*p]: 122.000000, 1.000000
#> *****Number of Trees: 200
#> *****Number of Cut Points: 2 ... 1
#> *****burn and ndpost: 200, 400
#> *****Prior:beta,alpha,tau,nu,lambda: 2.000000,0.950000,290.550176,3.000000,2184631.918632
#> *****sigma: 3348.915451
#> *****w (weights): 1.000000 ... 1.000000
#> *****Dirichlet:sparse,theta,omega,a,b,rho,augment: 0,0,1,0.5,1,5,0
#> *****nkeeptrain,nkeeptest,nkeeptestme,nkeeptreedraws: 400,400,400,400
#> *****printevery: 100000
#> *****skiptr,skipte,skipteme,skiptreedraws: 1,1,1,1
#>
#> MCMC
#> done 0 (out of 600)
#> time: 0s
#> check counts
#> trcnt,tecnt,temecnt,treedrawscnt: 400,400,400,400
#> *****Into main of wbart
#> *****Data:
#> data:n,p,np: 118, 5, 30
#> y1,yn: 5245.610169, 4149.610169
#> x1,x[n*p]: 121.000000, 0.000000
#> xp1,xp[np*p]: 122.000000, 1.000000
#> *****Number of Trees: 200
#> *****Number of Cut Points: 2 ... 1
#> *****burn and ndpost: 200, 400
#> *****Prior:beta,alpha,tau,nu,lambda: 2.000000,0.950000,289.984491,3.000000,2281761.566623
#> *****sigma: 3422.552953
#> *****w (weights): 1.000000 ... 1.000000
#> *****Dirichlet:sparse,theta,omega,a,b,rho,augment: 0,0,1,0.5,1,5,0
#> *****nkeeptrain,nkeeptest,nkeeptestme,nkeeptreedraws: 400,400,400,400
#> *****printevery: 100000
#> *****skiptr,skipte,skipteme,skiptreedraws: 1,1,1,1
#>
#> MCMC
#> done 0 (out of 600)
#> time: 1s
#> check counts
#> trcnt,tecnt,temecnt,treedrawscnt: 400,400,400,400#> RMSE_BART RMSE_LASSO
#> 3308.0 3292.5preds <- bind_rows(pred_store)
p_rmse <- ggplot(res, aes(RMSE_BART, RMSE_LASSO)) +
geom_point(size = 2.6, colour = "steelblue") +
geom_abline(slope = 1, intercept = 0, colour = "firebrick", linetype = "dashed") +
labs(title = "Out-of-sample RMSE across 10 splits",
subtitle = "Points below the line favour LASSO",
x = "RMSE (BART)", y = "RMSE (LASSO)") + theme_dspa(10)
p_pred <- ggplot(preds, aes(bart, lasso)) +
geom_point(alpha = 0.3, size = 1.2, colour = "grey30") +
geom_abline(slope = 1, intercept = 0, colour = "firebrick", linewidth = 0.9) +
geom_smooth(method = "lm", formula = y ~ x, colour = "steelblue",
linetype = "dashed", se = FALSE) +
labs(title = "Held-out predictions agree closely",
subtitle = "Red: perfect agreement. Blue dashed: OLS of LASSO on BART",
x = "BART prediction ($)", y = "LASSO prediction ($)") + theme_dspa(10)
p_rmse | p_pred
# --- Interactive equivalents ----------------------------------------------
plot_ly() |>
add_markers(x = res$RMSE_BART, y = res$RMSE_LASSO, name = "Splits") |>
add_trace(x = range(res$RMSE_BART), y = range(res$RMSE_BART),
type = "scatter", mode = "lines", name = "Equal RMSE",
line = list(width = 3)) |>
layout(title = "Out-of-sample RMSE: BART vs. LASSO",
xaxis = list(title = "RMSE (BART)"),
yaxis = list(title = "RMSE (LASSO)"),
legend = list(orientation = "h"))
plot_ly() |>
add_markers(x = preds$bart, y = preds$lasso, name = "Predictions") |>
add_trace(x = range(preds$bart), y = range(preds$bart), type = "scatter",
mode = "lines", name = "Ideal agreement", line = list(width = 3)) |>
layout(title = "Scatter plot predictions (BART vs. LASSO)",
xaxis = list(title = "BART predictions"),
yaxis = list(title = "LASSO predictions"),
legend = list(orientation = "h"))26.6 Posterior uncertainty, case by case
The reason to fit BART rather than a boosted ensemble is that every prediction comes with a distribution.
#> *****Into main of wbart
#> *****Data:
#> data:n,p,np: 148, 5, 0
#> y1,yn: -712.837838, 2893.162162
#> x1,x[n*p]: 122.000000, 1.000000
#> *****Number of Trees: 200
#> *****Number of Cut Points: 2 ... 1
#> *****burn and ndpost: 500, 2000
#> *****Prior:beta,alpha,tau,nu,lambda: 2.000000,0.950000,290.550176,3.000000,2178098.906176
#> *****sigma: 3343.904335
#> *****w (weights): 1.000000 ... 1.000000
#> *****Dirichlet:sparse,theta,omega,a,b,rho,augment: 0,0,1,0.5,1,5,0
#> *****nkeeptrain,nkeeptest,nkeeptestme,nkeeptreedraws: 2000,2000,2000,2000
#> *****printevery: 100000
#> *****skiptr,skipte,skipteme,skiptreedraws: 1,1,1,1
#>
#> MCMC
#> done 0 (out of 2500)
#> time: 3s
#> check counts
#> trcnt,tecnt,temecnt,treedrawscnt: 2000,0,0,2000#> [1] 2000 148sig_df <- data.frame(draw = seq_along(bart_long$sigma), sigma = bart_long$sigma)
ggplot(sig_df, aes(draw, sigma)) +
geom_line(alpha = 0.6, colour = "steelblue", linewidth = 0.3) +
geom_hline(yintercept = mean(bart_long$sigma), colour = "firebrick",
linewidth = 0.9) +
labs(title = expression(paste("Posterior draws of ", sigma)),
subtitle = "A stationary trace with no drift indicates the chain has converged",
x = "Draw (post burn-in)", y = expression(sigma)) +
theme_dspa()
# --- Interactive equivalent ------------------------------------------------
plot_ly(x = seq_along(bart_long$sigma), y = bart_long$sigma,
type = "scatter", mode = "markers", name = "sigma draws") |>
layout(title = "Post burn-in draws of sigma",
xaxis = list(title = "Number of posterior draws returned"),
yaxis = list(title = "sigma"))ggplot(data.frame(pred = bart_long$yhat.train.mean, obs = y_all),
aes(pred, obs)) +
geom_point(size = 2, alpha = 0.6, colour = "steelblue") +
geom_abline(slope = 1, intercept = 0, colour = "firebrick",
linetype = "dashed", linewidth = 0.9) +
labs(title = sprintf("Observed vs. BART-predicted charges (r = %.3f)",
cor(bart_long$yhat.train.mean, y_all)),
x = "BART posterior mean ($)", y = "Observed charges ($)") +
theme_dspa()
ord <- order(bart_long$yhat.train.mean)
sel_cases <- ord[round(seq(1, length(ord), length.out = 40))]
post <- as.data.frame(bart_long$yhat.train[, sel_cases, drop = FALSE])
names(post) <- paste0("case_", sel_cases)
post_long <- post |>
pivot_longer(everything(), names_to = "case", values_to = "draw") |>
mutate(case = factor(case, levels = paste0("case_", sel_cases)))
obs_pts <- data.frame(case = factor(paste0("case_", sel_cases),
levels = paste0("case_", sel_cases)),
observed = y_all[sel_cases])
ggplot(post_long, aes(case, draw)) +
geom_boxplot(outlier.shape = NA, fill = "grey88", colour = "grey40",
linewidth = 0.3) +
geom_point(data = obs_pts, aes(case, observed), colour = "darkgreen",
size = 2.4) +
labs(title = "Posterior predictive spread for 40 cases, ordered by fitted value",
subtitle = "Green points are the observed charges; many lie outside the posterior box",
x = "Case (ordered by posterior mean)", y = "Predicted charges ($)") +
theme_dspa(9) + theme(axis.text.x = element_blank())
# --- Interactive equivalent ------------------------------------------------
plot_ly() |>
add_trace(data = post_long, x = ~case, y = ~draw, type = "box",
color = ~case, showlegend = FALSE) |>
add_trace(data = obs_pts, x = ~case, y = ~observed, type = "scatter",
mode = "markers", name = "Observed charge",
marker = list(size = 14, color = "green",
line = list(color = "yellow", width = 2))) |>
layout(title = "Posterior predictive ranges with observed charges highlighted",
xaxis = list(title = "Cases"),
yaxis = list(title = "BART hospitalization charge prediction range"),
showlegend = FALSE)Read this figure honestly. Many observed charges fall
outside their posterior box, meaning the available
covariates (DRG, DIED, LOS,
AGE, gender) do not determine hospital charges. The right
conclusion is not that BART failed but that the information
content of the features is limited, and a model that reports
its own uncertainty is telling you so, which a point prediction never
would.
27 Computational complexity summary
Throughout, \(n\) = observations (rows), \(p\) = predictors (columns), \(m\) = trees, \(T\) = iterations. Flop counts are leading-order for the dense, real case.
| Operation | Time | Memory | Note |
|---|---|---|---|
| \(A+B\), \(A\circ B\) | \(O(mn)\) | \(O(mn)\) | Memory-bound, not compute-bound |
| \(AB\) (textbook) | \(O(mnk)\) | \(O(mk)\) | BLAS dgemm with cache blocking |
| \(AB\) (Strassen) | \(O(n^{2.807})\) | \(O(n^2)\) | Wins only at very large \(n\) |
\(A^\top A\) via
crossprod |
\(\tfrac{1}{2}np^2\) | \(O(p^2)\) | dsyrk exploits symmetry |
| LU factorization | \(\tfrac{2}{3}n^3\) | \(O(n^2)\) | Each extra right-hand side: \(2n^2\) |
| Cholesky | \(\tfrac{1}{3}n^3\) | \(O(n^2)\) | SPD only; no pivoting; also a PD test |
| QR (Householder) | \(2np^2-\tfrac{2}{3}p^3\) | \(O(np)\) | Backward stable; what lm() uses |
| Eigen (symmetric) | \(\approx 10n^3\) | \(O(n^2)\) | Tridiagonalize, then QR iteration |
| Eigen (general) | \(\approx 25n^3\) | \(O(n^2)\) | Hessenberg + shifted QR |
| SVD (thin, \(n\ge p\)) | \(\approx 14np^2\) | \(O(np)\) | Most expensive, most informative |
Explicit inverse solve(A) |
\(2n^3\) | \(O(n^2)\) | Three times LU; avoid |
| Cramer’s rule | \(O(n^4)\)–\(O(n!)\) | — | Never as an algorithm |
| OLS via normal equations | \(np^2+\tfrac{1}{3}p^3\) | \(O(p^2)\) | Fastest; error \(\propto\kappa_2(X)^2\) |
| OLS via QR | \(2np^2-\tfrac{2}{3}p^3\) | \(O(np)\) | \(\approx2\times\) slower; error \(\propto\kappa_2(X)\) |
| OLS via SVD | \(\approx 14np^2\) | \(O(np)\) | Slowest; handles rank deficiency |
| Hat matrix \(H\) (explicit) | \(O(n^2p)\) | \(\mathbf{O(n^2)}\) | Use hatvalues(), it needs only the diagonal |
hatvalues() (diagonal only) |
\(O(np)\) | \(O(n)\) | From the QR object |
| VIF (all \(p\)) | \(O(np^2)\) | \(O(p^2)\) | \(p\) auxiliary regressions, or one inverse |
| Stepwise selection | \(O(p^2)\times\) fit cost | — | Greedy; exhaustive is \(2^p\) |
| Regression tree (CART) | \(O(pn\log n)\) | \(O(n)\) | Sort once per variable, then running sums |
| Cost-complexity pruning + CV | \(O(k\cdot pn\log n)\) | \(O(n)\) | \(k\)-fold |
| Model tree (M5 / Cubist) | \(O(pn\log n+Lp^3)\) | \(O(Lp)\) | \(L\) leaves, one regression each |
| BART | \(O(T\,m\,n)\) | \(O(T_{\text{keep}}\,n)\) | Every sweep touches every tree |
Three entries deserve emphasis.
Never form \(H\)
explicitly. It is \(n\times
n\); at \(n=10^5\) that is 80 GB
in double precision. Leverage needs only its diagonal, and
hatvalues() extracts that from the QR factorization in
\(O(np)\).
The normal-equations speedup is real but purchased with digits. Roughly a factor of two, in exchange for squaring the condition number (§3.11).
BART’s cost is linear in \(n\) but multiplied by \(T\times m\). With \(m=200\) trees and \(T=1{,}200\) sweeps that is 240,000 tree-updates per fit, which is why BART sections dominate the runtime of this chapter.
28 Common pitfalls
| # | Pitfall | Consequence | Fix |
|---|---|---|---|
| 1 | A[i, ] where a matrix was expected |
Silent dimension drop → non-conformable arguments | drop = FALSE |
| 2 | Confusing * with %*% |
Wrong answer or a cryptic error | Hadamard vs. matrix product |
| 3 | u %*% v for an inner product |
Orientation-dependent; brittle | crossprod(u, v) / tcrossprod(u, v) |
| 4 | solve(A) %*% b |
3× the work, worse accuracy | solve(A, b) |
| 5 | Cramer’s rule or the adjugate as an algorithm | \(O(n^4)\) at best, unstable | LU, Cholesky, or QR |
| 6 | Least squares via \((X^\top X)^{-1}\) | Loses \(2\log_{10}\kappa\) digits | QR (qr.solve, lm.fit) |
| 7 | Benchmarking a coefficient formula against lm() |
Not like-for-like | Compare against lm.fit() |
| 8 | Verifying eigenvectors with \((\Lambda-A)V\) | Only correct when \(A\propto I\) | Check \(AV=V\Lambda\) |
| 9 | Believing OLS requires independent predictors | Drops informative variables; induces bias | Only \(\operatorname{rank}(X)=p\) is required |
| 10 | Writing \(\operatorname{cov}(x,x)=1\) | Confuses covariance with correlation | \(\operatorname{cov}(x,x)=\operatorname{var}(x)\) |
| 11 | Comparing models by \(R^2\) | Non-decreasing in \(p\); noise raises it | Adjusted \(R^2\), AIC/BIC, or held-out error |
| 12 | Fixed \(R^2\) or correlation thresholds | Field-dependent conventions, not facts | Report effect size with context |
| 13 | Raw x and x² in the same model |
\(\rho\approx0.99\); VIF explodes | Centre first, or use poly() |
| 14 | Dropping the main effect but keeping the interaction | Model not invariant to shifts of origin | Hierarchy principle |
| 15 | Q-Q plot of fitted values against residuals | Not a normality diagnostic | stat_qq() on standardized residuals |
| 16 | Hard-coding “influential observation 65” | Silently wrong after any data or seed change | Identify by threshold, report by name |
| 17 | Reading stars in a 30-coefficient model | Multiplicity: some will be “significant” by chance | Block \(F\)-test; FDR control |
| 18 | Forward step() without a scope |
Returns the starting model | scope = list(lower = ~1, upper = ...) |
| 19 | Reporting inference from a selected model | \(p\)-values ignore the search | Bootstrap stability; post-selection inference |
| 20 | Standardizing test data with its own mean and SD | Train and test in different coordinates | Compute centre/scale on training data only |
| 21 | Judging numerics by residual size alone | Small residual ≠ accurate solution | Check \(\kappa\); a residual measures backward error |
| 22 | Reaching for a flexible model by default | Costs accuracy when the truth is simple | Baseline first; flexibility must earn its place |
29 Practice problems
29.1 Problem 1: Rank, null space, and aliasing
Build a design matrix in which one column is an exact linear
combination of two others. Predict what lm() will report,
verify, and exhibit a basis for \(\mathcal{N}(X)\).
Solution
set.seed(101)
n1 <- 40
X1 <- cbind(1, rnorm(n1), rnorm(n1))
X1 <- cbind(X1, 3 * X1[, 2] - 2 * X1[, 3]) # column 4 is redundant
colnames(X1) <- c("int", "a", "b", "c")
y1 <- X1[, 2] + X1[, 3] + rnorm(n1)
c(columns = ncol(X1), rank = Matrix::rankMatrix(X1)[1])#> columns rank
#> 4 3#> (Intercept) X1[, -1]a X1[, -1]b X1[, -1]c
#> 0.192328 1.018158 1.110664 NAsv1 <- svd(X1)
tol <- max(dim(X1)) * .Machine$double.eps * max(sv1$d)
nullvec <- sv1$v[, sv1$d <= tol, drop = FALSE]
round(nullvec / nullvec[which.max(abs(nullvec))], 4)#> [,1]
#> [1,] 0.0000
#> [2,] 1.0000
#> [3,] -0.6667
#> [4,] -0.3333#> [,1]
#> [1,] 0
#> [2,] 0
#> [3,] 0
#> [4,] 0
#> [5,] 0
#> [6,] 0
#> [7,] 0
#> [8,] 0
#> [9,] 0
#> [10,] 0
#> [11,] 0
#> [12,] 0
#> [13,] 0
#> [14,] 0
#> [15,] 0
#> [16,] 0
#> [17,] 0
#> [18,] 0
#> [19,] 0
#> [20,] 0
#> [21,] 0
#> [22,] 0
#> [23,] 0
#> [24,] 0
#> [25,] 0
#> [26,] 0
#> [27,] 0
#> [28,] 0
#> [29,] 0
#> [30,] 0
#> [31,] 0
#> [32,] 0
#> [33,] 0
#> [34,] 0
#> [35,] 0
#> [36,] 0
#> [37,] 0
#> [38,] 0
#> [39,] 0
#> [40,] 0lm() reports NA rather
than an arbitrary choice.
29.2 Problem 2: Prove and verify \(\operatorname{tr}(H)=p\)
Prove that the trace of the hat matrix equals the number of columns of \(X\), and verify numerically for three different designs.
Solution
Proof. Using cyclicity of the trace, \[\operatorname{tr}(H)=\operatorname{tr}\!\big(X(X^\top X)^{-1}X^\top\big) =\operatorname{tr}\!\big((X^\top X)^{-1}X^\top X\big)=\operatorname{tr}(I_p)=p .\]
check_trace <- function(n, p, seed) {
set.seed(seed)
X <- cbind(1, matrix(rnorm(n * (p - 1)), n))
c(n = n, p = p, trace_H = sum(hatvalues(lm(rnorm(n) ~ X - 1))),
mean_leverage = mean(hatvalues(lm(rnorm(n) ~ X - 1))), p_over_n = p / n)
}
rbind(check_trace(50, 3, 1), check_trace(200, 8, 2), check_trace(1000, 15, 3)) |>
round(6)#> n p trace_H mean_leverage p_over_n
#> [1,] 50 3 3 0.060 0.060
#> [2,] 200 8 8 0.040 0.040
#> [3,] 1000 15 15 0.015 0.01529.3 Problem 3: Measure the \(\kappa^2\) penalty
Construct a design with \(\kappa_2(X)\approx 10^{7}\) and a known \(\boldsymbol\beta\). Solve by normal equations and by QR, and report the digits lost by each.
Solution
Xk <- make_design(300, 6, kappa = 1e7, seed = 202)
beta0 <- seq_len(6)
yk <- Xk %*% beta0
b_ne <- solve_normal(Xk, yk)
b_qr <- solve_qr(Xk, yk)
digits <- \(b) -log10(max(abs(b - beta0)) / max(abs(beta0)))
c(kappa_X = signif(kappa(Xk, exact = TRUE), 3),
kappa_XtX = signif(kappa(crossprod(Xk), exact = TRUE), 3),
correct_digits_normal_eqns = round(digits(b_ne), 1),
correct_digits_QR = round(digits(b_qr), 1))#> kappa_X kappa_XtX
#> 1.00e+07 1.00e+14
#> correct_digits_normal_eqns correct_digits_QR
#> 2.90e+00 1.01e+0129.4 Problem 4: Multicollinearity inflates variance without inducing bias
Simulate two predictors with correlation \(\rho\), fit OLS over many replicates, and confirm that \(\operatorname{sd}(\hat\beta_1)\) scales as \(\sqrt{\mathrm{VIF}}\) while \(E[\hat\beta_1]\) stays at the truth.
Solution
vif_study <- function(rho, n = 150, reps = 800, seed = 303) {
set.seed(seed)
L <- t(chol(matrix(c(1, rho, rho, 1), 2)))
b1 <- replicate(reps, {
Z <- t(L %*% matrix(rnorm(2 * n), 2))
y <- 1 + 2 * Z[, 1] + 3 * Z[, 2] + rnorm(n)
coef(lm(y ~ Z))[2]
})
data.frame(rho = rho, mean_b1 = mean(b1), sd_b1 = sd(b1),
VIF = 1 / (1 - rho^2))
}
vs <- do.call(rbind, lapply(c(0, 0.3, 0.6, 0.8, 0.9, 0.95, 0.99), vif_study))
vs$sd_ratio <- vs$sd_b1 / vs$sd_b1[1]
vs$sqrt_VIF <- sqrt(vs$VIF)
round(vs, 4)ggplot(vs, aes(sqrt_VIF, sd_ratio)) +
geom_abline(slope = 1, intercept = 0, colour = "firebrick", linetype = "dashed") +
geom_point(size = 2.6, colour = "steelblue") + geom_line(colour = "steelblue") +
labs(title = "Observed SE inflation versus the theoretical sqrt(VIF)",
x = expression(sqrt(VIF)), y = expression(sd(hat(beta)[1])/sd[rho==0])) +
theme_dspa()
mean_b1 stays at
2.0 throughout. Collinearity is a precision problem,
never an accuracy problem.
29.5 Problem 5: Eckart–Young by hand
Verify that the rank-\(k\) SVD truncation attains the theoretical Frobenius error \(\big(\sum_{j>k}\sigma_j^2\big)^{1/2}\), and confirm no other rank-\(k\) matrix does better.
Solution
set.seed(404)
Ay <- matrix(rnorm(60 * 30), 60)
sy <- svd(Ay)
check_k <- function(k) {
Ak <- sy$u[, 1:k] %*% diag(sy$d[1:k], k, k) %*% t(sy$v[, 1:k])
# A competitor: keep k components chosen at random rather than the top k
set.seed(k); pick <- sample(length(sy$d), k)
Ar <- sy$u[, pick] %*% diag(sy$d[pick], k, k) %*% t(sy$v[, pick])
c(k = k,
svd_error = norm(Ay - Ak, "F"),
theory = sqrt(sum(sy$d[-(1:k)]^2)),
spectral = norm(Ay - Ak, "2"),
sigma_k1 = sy$d[k + 1],
random_k_error = norm(Ay - Ar, "F"))
}
do.call(rbind, lapply(c(1, 5, 10, 20), check_k)) |> round(5)#> k svd_error theory spectral sigma_k1 random_k_error
#> [1,] 1 39.4804 39.4804 11.96676 11.96676 41.2248
#> [2,] 5 32.5819 32.5819 9.98429 9.98429 37.8444
#> [3,] 10 25.3422 25.3422 8.19112 8.19112 33.8201
#> [4,] 20 12.8286 12.8286 5.13535 5.13535 20.5050svd_error matches theory and
spectral matches \(\sigma_{k+1}\) exactly, as
Eckart–Young–Mirsky guarantees. The randomly chosen rank-\(k\) competitor is always worse, the top
\(k\) components are optimal, not
merely convenient.
29.6 Problem 6: Robust standard errors
Simulate heteroscedastic errors, show that OLS remains unbiased but its nominal standard errors are wrong, and repair them with a sandwich estimator.
Solution
set.seed(505)
n6 <- 400
x6 <- runif(n6, 1, 10)
y6 <- 2 + 1.5 * x6 + rnorm(n6, sd = 0.4 * x6) # SD grows with x
m6 <- lm(y6 ~ x6)
se_ols <- coef(summary(m6))[, "Std. Error"]
se_robust <- sqrt(diag(sandwich::vcovHC(m6, type = "HC3")))
# Ground truth: the empirical SD of the estimator over many replicates
true_sd <- apply(replicate(2000, {
xx <- runif(n6, 1, 10)
yy <- 2 + 1.5 * xx + rnorm(n6, sd = 0.4 * xx)
coef(lm(yy ~ xx))
}), 1, sd)
rbind(nominal_OLS = se_ols, robust_HC3 = se_robust, empirical_truth = true_sd) |>
round(5)#> (Intercept) x6
#> nominal_OLS 0.26888 0.04540
#> robust_HC3 0.20645 0.05015
#> empirical_truth 0.20796 0.05020sandwich::vcovHC() with
lmtest::coeftest() whenever the scale–location plot trends.
29.7 Problem 7: Trees cannot extrapolate
Fit a regression tree and a linear model on \(x\in[0,5]\) and predict on \(x\in[5,10]\). Explain the difference.
Solution
set.seed(606)
xtr <- runif(200, 0, 5); ytr <- 3 + 2 * xtr + rnorm(200, sd = 0.6)
xte <- seq(0, 10, by = 0.05)
t7 <- rpart(y ~ x, data = data.frame(x = xtr, y = ytr))
l7 <- lm(y ~ x, data = data.frame(x = xtr, y = ytr))
pr <- data.frame(x = xte,
tree = predict(t7, data.frame(x = xte)),
linear = predict(l7, data.frame(x = xte)),
truth = 3 + 2 * xte) |>
pivot_longer(-x, names_to = "model", values_to = "pred")
ggplot(pr, aes(x, pred, colour = model)) +
geom_line(linewidth = 1) +
geom_point(data = data.frame(x = xtr, y = ytr), aes(x, y),
inherit.aes = FALSE, alpha = 0.2, size = 0.8) +
geom_vline(xintercept = 5, linetype = "dashed", colour = "grey40") +
scale_colour_manual(values = c(tree = "firebrick", linear = "steelblue",
truth = "black")) +
labs(title = "Beyond the training range, the tree flatlines",
subtitle = "Training data end at x = 5 (dashed line)",
x = "x", y = "Prediction", colour = NULL) +
theme_dspa()
c(max_tree_prediction = max(predict(t7, data.frame(x = xte))),
max_training_y = max(ytr),
truth_at_x10 = 23)#> max_tree_prediction max_training_y truth_at_x10
#> 12.5876 13.5098 23.0000max(ytr). The linear model extrapolates correctly here
because the true relationship is linear; a parametric
form that is wrong will extrapolate confidently and wrongly, which is
its own failure mode.
29.8 Problem 8: Ridge as an eigenvalue shift
Show that ridge regression, \(\hat\beta_\lambda=(X^\top X+\lambda I)^{-1}X^\top y\), shifts every eigenvalue of \(X^\top X\) upward by \(\lambda\), and demonstrate the effect on the condition number and on estimator variance.
Solution
If \(X^\top X=Q\Lambda Q^\top\) then \(X^\top X+\lambda I=Q(\Lambda+\lambda I)Q^\top\), so every eigenvalue becomes \(\lambda_j+\lambda\) and
\[\kappa_2(X^\top X+\lambda I)=\frac{\lambda_{\max}+\lambda}{\lambda_{\min}+\lambda}\;<\;\frac{\lambda_{\max}}{\lambda_{\min}} .\]
The floor at \(\lambda\) is what makes the problem solvable even when \(X^\top X\) is singular, including when \(p>n\).
set.seed(707)
Xr <- make_design(120, 10, kappa = 1e6, seed = 707)
btrue <- rep(1, 10)
yr <- Xr %*% btrue + rnorm(120, sd = 0.05)
XtX <- crossprod(Xr)
ev <- eigen(XtX, symmetric = TRUE)$values
lams <- 10^seq(-10, 0, length.out = 40)
ridge <- do.call(rbind, lapply(lams, function(l) {
bl <- solve(XtX + l * diag(10), crossprod(Xr, yr))
c(lambda = l,
kappa = (max(ev) + l) / (min(ev) + l),
bias2 = sum((colMeans(t(bl)) - btrue)^2),
mse = sum((bl - btrue)^2))
}))
ridge <- as.data.frame(ridge)
ggplot(ridge, aes(lambda, mse)) +
geom_line(linewidth = 1, colour = "steelblue") +
geom_point(data = ridge[which.min(ridge$mse), ], colour = "firebrick",
size = 3) +
scale_x_log10() + scale_y_log10() +
labs(title = "Ridge trades bias for variance and wins on total error",
subtitle = sprintf("Optimal lambda = %.2e reduces MSE by a factor of %.0f vs. OLS",
ridge$lambda[which.min(ridge$mse)],
ridge$mse[1] / min(ridge$mse)),
x = expression(lambda), y = expression(group("||", hat(beta)-beta, "||")^2)) +
theme_dspa()
30 Checkpoint
- A colleague reports that their hand-coded
solve(t(X) %*% X) %*% t(X) %*% yis twice as fast aslm(). What two things would you tell them? - A design has \(\kappa_2(X)=10^{5}\). Roughly how many correct digits survive a QR solve, and how many survive a normal-equations solve?
- In a fitted model \(\sum_i h_{ii}=7\). What does that tell you about \(X\)?
- Adding a predictor raised \(R^2\) from 0.61 to 0.63 but lowered adjusted \(R^2\). What happened, and which do you report?
- Two predictors correlate at 0.97. Is \(\hat\beta\) biased? What is affected, and by how much?
- Why can a regression tree never predict a value larger than the largest response in its training data, and when does that matter?
Answers
- First, the comparison is not like-for-like,
lm()also builds a model frame, expands factors, and computes standard errors, residuals, and the full summary; the fair comparison islm.fit(). Second, and more importantly, the speed comes from skipping the QR factorization that protects accuracy: forming \(X^\top X\) squares the condition number, so the fast method loses twice as many digits. Useqr.solve()orlm.fit()when only coefficients are needed. - Double precision carries about 16 digits. QR loses \(\log_{10}\kappa=5\), leaving about 11. Normal equations lose \(\log_{10}\kappa^2=10\), leaving about 6.
- \(\operatorname{tr}(H)=p\), so \(X\) has 7 columns, 7 estimated coefficients including the intercept. It also fixes mean leverage at \(7/n\) and residual degrees of freedom at \(n-7\).
- A predictor was added that explains less variance than the degree of freedom it consumed. \(R^2\) cannot decrease when a column is added, because the larger column space contains the smaller one. Adjusted \(R^2\) penalizes \(p\) and can decrease. Report adjusted \(R^2\), or better, held-out error or AIC/BIC.
- Not biased, unbiasedness requires only A1, A2, and A4. What inflates is the variance: \(\mathrm{VIF}=1/(1-0.97^2)\approx 16.9\), so each standard error is about \(\sqrt{16.9}\approx 4.1\) times larger than it would be with orthogonal predictors. Confidence intervals widen by the same factor; fitted values are essentially unaffected.
- Each terminal node predicts the mean of the training responses falling in it, and a mean is bounded by the values it averages. This matters for any extrapolation task, forecasting growth, dose–response beyond the tested range, pricing outside the observed band, where a tree will confidently flatline. It is also why tree ensembles are poor at trend extrapolation and why hybrid model trees (linear fits in the leaves) exist.
31 Summary
Linear algebra
- Rank, not independence, is the condition that makes a model identifiable; the four fundamental subspaces organize every solvability question.
- Least squares is orthogonal projection. The hat matrix \(H=X(X^\top X)^{-1}X^\top\) is symmetric idempotent with \(\operatorname{tr}(H)=p\), and that single fact produces degrees of freedom, leverage, and the ANOVA decomposition.
- Decompositions turn hard problems into easy ones: LU for square solves, Cholesky for symmetric positive definite systems, QR for least squares, eigen for symmetric structure, SVD for everything else.
- The condition number governs how many digits survive. \(\kappa_2(X^\top X)=\kappa_2(X)^2\) is why
lm()uses QR and accepts being twice as slow.
Regression
- OLS follows from calculus, geometry, or Gaussian maximum likelihood, three routes to one estimator. Gauss–Markov says it is best among linear unbiased estimators, a qualifier that leaves room for ridge and LASSO.
- \(\operatorname{Var}(\hat{\boldsymbol\beta})=\sigma^2(X^\top X)^{-1}\) generates every standard error, \(t\)-statistic, and confidence interval.
- Correlated predictors inflate variance by \(\mathrm{VIF}\) but never bias \(\hat{\boldsymbol\beta}\). When \(k\gg n\), spurious correlation is guaranteed at scale \(2\sqrt{\log k/n}\).
- Diagnostics are functions of \(H\): leverage is \(h_{ii}\), standardized residuals divide by \(\sqrt{1-h_{ii}}\), and Cook’s distance factors as outlyingness \(\times\) leverage.
- \(R^2\) cannot decrease when predictors are added. Stepwise selection is greedy, unstable, and invalidates the inference it reports.
Trees
- CART partitions greedily by variance reduction at \(O(pn\log n)\), producing a piecewise-constant surface that cannot extrapolate.
- Model trees restore local linearity inside each region.
- BART places a prior over an ensemble of deliberately weak trees and samples the posterior, so every prediction arrives with calibrated uncertainty, at a cost of \(O(Tmn)\), and with no advantage when the truth is linear.
Where these threads continue
| Thread | Continues in |
|---|---|
| SVD, eigendecomposition, PCA | Dimensionality reduction |
| Classification with the same design-matrix machinery | Supervised classification |
| Kernel methods; trees inside ensembles | Black-box methods |
| Cross-validation, ROC, calibration | Model assessment |
| Ridge, LASSO, elastic net, FDR control | Feature selection |
| Mixed models, GEE, autocorrelated errors | Longitudinal analysis |
| Gradient descent, Newton, conjugate gradient | Function optimization |
32 Chapter roadmap
- Chapter 1:
Foundations. R toolchain, reproducibility conventions,
dspa_read(), simulation. - Chapter 2: Data quality and exploratory visual analytics. Centrality, dispersion, density estimation, missingness, contingency tables, the grammar of graphics.
- Dimensionality reduction. PCA, ICA, factor analysis, t-SNE, UMAP, all built on the SVD of §3.6.5.
- Supervised classification. kNN, naive Bayes, decision trees, logistic regression.
- Black-box methods. Neural networks, SVM, random forests, gradient boosting.
- Model assessment, validation, improvement. Cross-validation, ROC/AUC, calibration, hyperparameter tuning.
- Variable importance and feature selection. Ridge, LASSO, elastic net, stability selection, FDR control.
- Longitudinal and time-series analysis. Mixed models, GEE, ARIMA, state-space models.
- Function optimization. Gradient descent, Newton and quasi-Newton methods, EM, Bayesian optimization.
33 Session information
#> R version 4.3.3 (2024-02-29 ucrt)
#> Platform: x86_64-w64-mingw32/x64 (64-bit)
#> Running under: Windows 11 x64 (build 26200)
#>
#> Matrix products: default
#>
#>
#> locale:
#> [1] LC_COLLATE=English_United States.utf8
#> [2] LC_CTYPE=English_United States.utf8
#> [3] LC_MONETARY=English_United States.utf8
#> [4] LC_NUMERIC=C
#> [5] LC_TIME=English_United States.utf8
#>
#> time zone: America/New_York
#> tzcode source: internal
#>
#> attached base packages:
#> [1] stats graphics grDevices utils datasets methods base
#>
#> other attached packages:
#> [1] Cubist_0.4.4 lattice_0.22-6 rpart.plot_3.1.2 rpart_4.1.23
#> [5] DT_0.33 Matrix_1.6-5 plotly_4.12.0 patchwork_1.3.0
#> [9] tidyr_1.3.1 dplyr_1.1.4 ggplot2_4.0.1
#>
#> loaded via a namespace (and not attached):
#> [1] sandwich_3.1-0 sass_0.4.9 generics_0.1.3 stringi_1.8.4
#> [5] digest_0.6.37 magrittr_2.0.3 evaluate_1.0.3 grid_4.3.3
#> [9] RColorBrewer_1.1-3 fastmap_1.2.0 plyr_1.8.9 jsonlite_1.8.9
#> [13] survival_3.7-0 GGally_2.2.1 httr_1.4.7 mgcv_1.9-1
#> [17] purrr_1.0.2 crosstalk_1.2.1 viridisLite_0.4.2 scales_1.4.0
#> [21] lazyeval_0.2.2 jquerylib_0.1.4 abind_1.4-5 cli_3.6.3
#> [25] rlang_1.1.5 splines_4.3.3 withr_3.0.2 cachem_1.1.0
#> [29] yaml_2.3.10 otel_0.2.0 tools_4.3.3 reshape2_1.4.4
#> [33] SparseM_1.84 MatrixModels_0.5-3 ggstats_0.6.0 vctrs_0.6.5
#> [37] R6_2.6.1 zoo_1.8-12 lifecycle_1.0.5 stringr_1.5.1
#> [41] car_3.1-2 htmlwidgets_1.6.4 MASS_7.3-60.0.1 pkgconfig_2.0.3
#> [45] pillar_1.10.1 bslib_0.9.0 gtable_0.3.6 Rcpp_1.0.14
#> [49] data.table_1.16.4 glue_1.8.0 xfun_0.52 tibble_3.2.1
#> [53] tidyselect_1.2.1 rstudioapi_0.18.0 knitr_1.51 farver_2.1.2
#> [57] nlme_3.1-165 htmltools_0.5.8.1 carData_3.0-5 rmarkdown_2.31
#> [61] labeling_0.4.3 compiler_4.3.3 quantreg_5.98 S7_0.2.1