6. Matrix algebra for regression

Dwaine Studios runs a chain of portrait studios and wants to grow. It already operates in 21 cities, and its managers track two numbers for each: the population aged 16 and younger, the pool of children whose parents buy portraits, and per-capita disposable income, how much money families have to spend. They want to use these two numbers to predict studio sales, so they can rank new cities and open where the forecast is highest. The question is a forecasting one, but underneath it is a modeling one: how do sales depend on both predictors at once?

In Chapter 2 you fit a line with a single predictor by solving two normal equations for the intercept and slope (2.2 Least squares from first principles). That worked because there were only two unknowns. Dwaine has three, one intercept and two slopes, and a realistic study might have twenty. Writing out a separate normal equation for each and solving the tangle by hand is not something anyone wants to do twice. Matrix algebra is the better language: it packs a whole dataset into two symbols, the normal equations into one, and the solution into a formula you write in a single line and compute in a single command.

None of what follows requires you to have seen a matrix before. We build every idea from its definition, test it on a small example you could check with a pencil, then run it on the real 21-city dataset in both R and Python. If you have taken linear algebra, this will read like review with a statistical accent; if not, you will still finish able to fit and explain a multiple regression, the only matrix algebra this course asks of you.

Figure 1 shows that both predictors matter: sales climb as the young population grows, and for a given population the higher-income cities (darker points) tend to sell more. A model that uses both at once is a multiple regression model, and every such model is built, fit, and understood through matrices. By the end you will have built, for the real Dwaine data, the design matrix $\mathbf{X}$, the product $\mathbf{X}'\mathbf{X}$, its inverse, the coefficient estimates, the hat matrix, and every coefficient’s standard error, each from operations you can also do by hand on a small example.

6.1 The data as a vector and a matrix¶

Intuition¶

Start with the response. Dwaine has 21 sales figures, one per city. Stack them into a single column and you have a vector (Definition 6.1), an ordered list of numbers written vertically. Call it $\mathbf{Y}$. A vector is the simplest matrix: many rows, one column.

Now the predictors. Each city carries two numbers, so the predictors form a table with 21 rows and 2 columns, and a table of numbers with rows and columns is a matrix (Definition 6.1). For regression we glue one extra column of all ones onto the front: a bookkeeping trick that lets the intercept ride along as just another coefficient. The result is the design matrix $\mathbf{X}$ (Definition 6.2), with 21 rows and 3 columns. Figure 2 shows how the whole model lines up as stacked arrays.

Reading across a row gives one city’s complete record: a 1, its young population, its disposable income, and, in $\mathbf{Y}$, its sales. Stacking the rows turns 21 separate little equations into the single statement $\mathbf{Y} = \mathbf{X}\boldsymbol{\beta} + \boldsymbol{\varepsilon}$. That compression is the reason to learn matrices: the same three symbols describe a dataset with 21 rows or 21 million.

Formula¶

For a regression with $n$ observations and $p$ parameters (counting the intercept), the pieces are

• $\mathbf{Y}$ is the response vector; row $i$ is the response for case $i$.

• $\mathbf{X}$ is the design matrix; its first column is all ones, and column $k+1$ holds the values of predictor $k$.

• $\boldsymbol{\beta}$ is the parameter vector: the intercept $\beta_0$ followed by the slopes.

• $p$ is the number of regression parameters including the intercept ($p = 3$ for Dwaine), and $n$ is the number of observations ($n = 21$).

For Dwaine, $n = 21$ and $p = 3$: one column of ones, one column of under-16 populations, one column of disposable incomes. The whole model of Chapter 2 generalizes to $\mathbf{Y} = \mathbf{X}\boldsymbol{\beta} + \boldsymbol{\varepsilon}$, which we unpack in 7.1 The model and least squares in matrix form. In words: each observed response equals a weighted sum of that row’s predictors, using the same weights $\boldsymbol{\beta}$ for every row, plus a random error.

R and Python¶

6.2 Transpose and matrix multiplication¶

Intuition¶

We have the data in matrices; now we need to combine them. Two operations do almost all the work. The transpose (Definition 6.3) flips a matrix on its diagonal so rows become columns. Matrix multiplication (Definition 6.4) is the engine that turns the design matrix into the sums of squares and cross-products a regression runs on.

Matrix multiplication is not what a newcomer expects. One entry of the product is built from a whole row of the left matrix and a whole column of the right: line them up, multiply term by term, and add. Applied to $\mathbf{X}'\mathbf{X}$, that rule produces exactly the sums $\sum X_{ik} X_{ij}$ you would otherwise compute one painful sum at a time. Figure 3 shows it on a small example.

Two products, and only two, matter for the fit: $\mathbf{X}'\mathbf{X}$ summarizes the predictors by themselves, and $\mathbf{X}'\mathbf{Y}$ summarizes how each predictor moves with the response. Between them they hold every summary least squares needs, which is why a whole dataset can shrink to a handful of numbers: the fit does not care about individual rows once these sums are known.

Formula¶

Picture it as tipping a matrix over its top-left corner: each row stands up to become a column, so a wide matrix tips into a tall one and back. A column vector, tipped over, becomes a single row, which is why $\mathbf{u}'\mathbf{u}$ multiplies a lying-down copy of a vector by a standing-up copy and collapses to one number.

• The inner dimensions ($k$ and $k$) must match; they are what you sum over.

• The outer dimensions ($m$ and $q$) give the shape of the answer.

• Entry $(i,j)$ is row $i$ of $\mathbf{A}$ dotted with column $j$ of $\mathbf{B}$.

Two products carry the whole load in regression. With $\mathbf{X}$ of size $n \times p$, the transpose $\mathbf{X}'$ is $p \times n$, so

In words: $\mathbf{X}'\mathbf{X}$ collects every sum of squares and cross-product of the predictor columns into a small square table, and $\mathbf{X}'\mathbf{Y}$ collects every predictor-times-response sum into a short column, the raw materials of the normal equations. A rule we will use constantly is that transposing a product reverses the order of its factors.

Proof. We show $(\mathbf{A}\mathbf{B})' = \mathbf{B}'\mathbf{A}'$ by comparing entries. The $(i,j)$ entry of $(\mathbf{A}\mathbf{B})'$ is, by the definition of transpose, the $(j,i)$ entry of $\mathbf{A}\mathbf{B}$, which is $\sum_{\ell} A_{j\ell} B_{\ell i}$. The $(i,j)$ entry of $\mathbf{B}'\mathbf{A}'$ is row $i$ of $\mathbf{B}'$ dotted with column $j$ of $\mathbf{A}'$, that is $\sum_{\ell} (\mathbf{B}')_{i\ell} (\mathbf{A}')_{\ell j} = \sum_{\ell} B_{\ell i} A_{j \ell}$. The two sums have the same terms, so the matrices are equal. $\blacksquare$

R and Python¶

First the small example from Figure 3, so the mechanics are concrete before we turn the engine loose on the real data. In R, %*% is matrix multiplication (a plain * multiplies entry by entry instead); in Python the operator is @.

A <- matrix(c(1, 2, 3,
              4, 5, 6), nrow = 2, byrow = TRUE)
B <- matrix(c(1, 0,
              0, 1,
              1, 1), nrow = 3, byrow = TRUE)
A %*% B

     [,1] [,2]
[1,]    4    5
[2,]   10   11

A = np.array([[1, 2, 3],
              [4, 5, 6]])
B = np.array([[1, 0],
              [0, 1],
              [1, 1]])
print(A @ B)

[[ 4  5]
 [10 11]]

The row-by-column rule is easier to trust once you have watched it happen on the one product a regression actually needs, so build $\mathbf{X}'\mathbf{X}$ yourself before moving on.

6.3 Identity, symmetry, independence, and rank¶

Intuition¶

Before inverting anything, we need vocabulary for the shapes of matrices that show up in regression, and one idea that decides whether a fit is even possible. The identity matrix $\mathbf{I}$ (Definition 6.6) is the matrix version of the number 1: multiplying by it changes nothing. A symmetric matrix (Definition 6.7) equals its own transpose, and $\mathbf{X}'\mathbf{X}$ is always symmetric, which is why its table looks like a mirror across the diagonal.

The idea that decides feasibility is linear independence (Definition 6.8) of the predictor columns. If one predictor column is an exact linear combination of the others, say dispoinc were exactly twice targtpop plus a constant, then the columns carry redundant information, the matrix $\mathbf{X}$ is rank deficient, and $\mathbf{X}'\mathbf{X}$ cannot be inverted. There is then no unique best fit, because the data cannot tell the redundant coefficients apart. This is the matrix face of a problem you will meet again as multicollinearity in 12.1 Multicollinearity and the variance inflation factor.

A picture makes independence concrete. Think of two columns as arrows drawn from the origin. If they point in genuinely different directions, they open up an area between them, and the matrix has full rank. If one is just a stretched copy of the other, both arrows lie on a single line, the area between them is zero, and the matrix is rank deficient. Figure 5 shows both cases side by side.

These names earn their keep: a named property is a fact you can cite instead of recompute. Showing the hat matrix is a projection becomes a one-line argument, “it is symmetric and idempotent,” in place of pages of entry-by-entry checking.

Formula¶

The product $\mathbf{X}'\mathbf{X}$ is always symmetric, because $(\mathbf{X}'\mathbf{X})' = \mathbf{X}'(\mathbf{X}')' = \mathbf{X}'\mathbf{X}$ by the order-reversing rule (Theorem 6.5) of 6.2 Transpose and matrix multiplication.

The central fact of this section connects rank to invertibility.

A tiny example makes rank concrete: $(1, 0)'$ and $(0, 1)'$ are independent, but $(1, 2)'$ and $(2, 4)'$ are dependent because the second is exactly twice the first, so together they span only a line and have rank 1, not 2. In a design matrix, such a dependent column is a predictor that repeats information already present.

R and Python¶

Software reports symmetry, rank, and the determinant (the single number, covered in 6.4 The inverse and the normal equations, whose vanishing signals rank deficiency) directly.

isSymmetric(XtX)
qr(X)$rank
round(det(XtX), 2)

[1] TRUE
[1] 3
[1] 1068302

print(np.allclose(XtX, XtX.T))
print(np.linalg.matrix_rank(X))
print(round(np.linalg.det(XtX), 2))

True
3
1068302.4

The design matrix has rank 3, its full column count, so the two predictors plus the intercept carry three genuinely different pieces of information, and the nonzero determinant confirms that $\mathbf{X}'\mathbf{X}$ is invertible and the Dwaine fit is well posed.

The habit to build: before trusting any multiple regression, ask whether the predictors are genuinely distinct. Software will happily invert a nearly singular matrix and return coefficients with enormous standard errors, the early warning that two predictors carry almost the same information. Rank is the yes-or-no version of that question; Chapter 12 sharpens it into “how close to dependent are they?”

6.4 The inverse and the normal equations¶

Intuition¶

Here is where the machinery pays off. In ordinary algebra, to solve $ax = c$ you divide by $a$, which is the same as multiplying by $a^{-1}$. Matrices have no division, but they have the next best thing: for a square matrix $\mathbf{A}$ with full rank, there is an inverse $\mathbf{A}^{-1}$ (Definition 6.10) that undoes it, meaning $\mathbf{A}^{-1}\mathbf{A} = \mathbf{I}$. To solve a matrix equation $\mathbf{A}\mathbf{x} = \mathbf{c}$, multiply both sides on the left by $\mathbf{A}^{-1}$ and get $\mathbf{x} = \mathbf{A}^{-1}\mathbf{c}$.

The normal equations of least squares, which in Chapter 2 were two scalar equations, collapse in matrix form to the single equation $\mathbf{X}'\mathbf{X}\,\mathbf{b} = \mathbf{X}'\mathbf{Y}$. Chapter 7 derives this from minimizing squared error (7.1 The model and least squares in matrix form); for now take it as the matrix twin of the Chapter 2 normal equations (2.2 Least squares from first principles). Solving it is now a one-liner: multiply by $(\mathbf{X}'\mathbf{X})^{-1}$.

With numbers, $a^{-1}$ exists for every $a$ except zero. With matrices, $\mathbf{A}^{-1}$ exists for every square matrix except the rank-deficient ones, so “can I solve this uniquely,” “is the determinant nonzero,” and “are the columns independent” are three ways of asking one question. For Dwaine, all three say yes.

Formula¶

• The determinant $ad - bc$ measures whether the columns are independent; if it is zero, the inverse formula divides by zero and no inverse exists.

• For larger matrices the determinant is more involved, so we let software compute it, but the rule “invertible if and only if determinant nonzero” stays exact.

This connects back to the picture of rank in Figure 5: the two columns of a two-by-two matrix frame a parallelogram, and the determinant is its area. When the columns point along the same line the parallelogram is squashed flat, its area is zero, and the matrix is singular, so “determinant zero” and “the columns collapse onto one line” are the same event.

Applying the inverse to the normal equations gives the least-squares estimator in one line.

In words: multiply the stored cross-products $\mathbf{X}'\mathbf{Y}$ by the inverse of the cross-product matrix $\mathbf{X}'\mathbf{X}$, and out come all $p$ coefficient estimates at once. This one formula is the whole of least-squares estimation, for any number of predictors.

R and Python¶

First the $2 \times 2$ inverse formula, checked on a small symmetric matrix, so you can see $\mathbf{A}^{-1}\mathbf{A} = \mathbf{I}$ happen. Both languages invert with a single call: solve in R, np.linalg.inv in Python.

M <- matrix(c(2, 1,
              1, 2), nrow = 2, byrow = TRUE)
Minv <- solve(M)
Minv
round(M %*% Minv, 6)

           [,1]       [,2]
[1,]  0.6666667 -0.3333333
[2,] -0.3333333  0.6666667
     [,1] [,2]
[1,]    1    0
[2,]    0    1

M = np.array([[2.0, 1.0],
              [1.0, 2.0]])
Minv = np.linalg.inv(M)
print(Minv)
print(np.round(M @ Minv, 6))

[[ 0.66666667 -0.33333333]
 [-0.33333333  0.66666667]]
[[1. 0.]
 [0. 1.]]

The determinant is $2\cdot 2 - 1 \cdot 1 = 3$, so the inverse is $\tfrac{1}{3}\left(\begin{smallmatrix} 2 & -1 \\ -1 & 2 \end{smallmatrix}\right)$, matching the printout, and $\mathbf{M}^{-1}\mathbf{M}$ returns the identity.

You will almost never invert a matrix bigger than this by hand; working the two-by-two case once is enough to see where the answer comes from and watch the identity appear, so a software inverse is never magic. From here on we let solve and np.linalg.inv do the arithmetic, and the next example turns the inverse loose on the full Dwaine fit.

A determinant of zero is easier to feel than to define, so pull the two columns together yourself and watch the inverse come apart long before they meet.

6.5 Idempotent and projection matrices¶

Intuition¶

The coefficients are done, so the model can now make its own predictions. The question of this section is where those predictions come from and what they look like as a picture. The answer is a single table of numbers that turns the observed sales into the model’s fitted sales in one step, and that step is a shadow: it flattens the data onto the surface the predictors can reach.

With $\mathbf{b}$ in hand, the fitted values are $\widehat{\mathbf{Y}} = \mathbf{X}\mathbf{b}$. Substitute the estimator and something remarkable appears:

Geometrically the hat matrix is a projection (Definition 6.14). Picture $\mathbf{Y}$ as an arrow in $n$-dimensional space; the fitted values a regression can produce all lie in a flat subspace, the span of the columns of $\mathbf{X}$. The hat matrix drops $\mathbf{Y}$ straight down onto that subspace, landing at the closest point, and the residual is the perpendicular leftover. Figure 7 shows this picture, the subject of 7.2 The geometry of least squares.

Two algebraic properties make $\mathbf{H}$ a projection, and they are the workhorses of the diagnostics in Chapter 9. First, $\mathbf{H}$ is symmetric. Second, $\mathbf{H}$ is idempotent: applying it twice is the same as applying it once, $\mathbf{H}\mathbf{H} = \mathbf{H}$. That makes sense for a projection, since once a vector is flattened onto the plane, flattening again does nothing.

The shadow picture is worth holding onto. Your shadow on flat ground keeps your left-right and forward-back position but flattens your height to zero; the hat matrix does the same to $\mathbf{Y}$, flattening away the part of the response that no combination of the predictors could reproduce, which is the residual.

These properties are not decoration: Chapter 9 reads the diagonal of $\mathbf{H}$ to find which cities pull hardest on the fitted surface, Chapter 7 uses its trace to count degrees of freedom, and the perpendicular picture is what makes least squares the closest fit, with no point in the plane nearer to $\mathbf{Y}$ than its own shadow.

Formula¶

The two projections of regression are

with $\widehat{\mathbf{Y}} = \mathbf{H}\mathbf{Y}$ the fitted values and $\mathbf{e} = \mathbf{Y} - \widehat{\mathbf{Y}} = (\mathbf{I} - \mathbf{H})\mathbf{Y}$ the residuals. In words: $\mathbf{H}$ projects onto the space of fittable values, and $\mathbf{I} - \mathbf{H}$ projects onto the leftover space of residuals.

The two defining properties of $\mathbf{H}$ are what make it a projection.

Proof. Write $\mathbf{G} = (\mathbf{X}'\mathbf{X})^{-1}$, which is symmetric because $\mathbf{X}'\mathbf{X}$ is (its inverse inherits symmetry: $\mathbf{G}' = ((\mathbf{X}'\mathbf{X})^{-1})' = ((\mathbf{X}'\mathbf{X})')^{-1} = (\mathbf{X}'\mathbf{X})^{-1} = \mathbf{G}$). For symmetry of $\mathbf{H}$, apply the order-reversing transpose rule of 6.2 Transpose and matrix multiplication twice:

For idempotency, multiply $\mathbf{H}$ by itself and cancel the inner $\mathbf{X}'\mathbf{X}$ against its inverse:

using $(\mathbf{X}'\mathbf{X})\mathbf{G} = (\mathbf{X}'\mathbf{X})(\mathbf{X}'\mathbf{X})^{-1} = \mathbf{I}$. The same two steps show $\mathbf{I} - \mathbf{H}$ is symmetric and idempotent: $(\mathbf{I} - \mathbf{H})(\mathbf{I} - \mathbf{H}) = \mathbf{I} - 2\mathbf{H} + \mathbf{H}\mathbf{H} = \mathbf{I} - 2\mathbf{H} + \mathbf{H} = \mathbf{I} - \mathbf{H}$. $\blacksquare$

R and Python¶

The projection is the one idea in this chapter a flat page cannot really show, so turn the picture around and move $\mathbf{Y}$ with your own finger.

6.6 Quadratic forms and sums of squares¶

Intuition¶

Every “sum of squares” in regression is secretly one matrix expression called a quadratic form (Definition 6.17). The error sum of squares is $\mathrm{SSE} = \sum e_i^2$, and a sum of squared entries of a vector is exactly that vector dotted with itself: $\mathrm{SSE} = \mathbf{e}'\mathbf{e}$. Because $\mathbf{e} = (\mathbf{I} - \mathbf{H})\mathbf{Y}$, a little algebra rewrites SSE using only $\mathbf{Y}$ and a projection matrix in the middle. That middle-matrix pattern, $\mathbf{Y}'\mathbf{A}\mathbf{Y}$, is a quadratic form, and SSTO, SSR, and SSE are all of this shape.

Why care? Two reasons. First, a quadratic form makes the degrees of freedom and the expected value of each sum of squares fall out of the middle matrix, how Chapter 7 explains the ANOVA table. Second, sums of squares can never be negative, and that is guaranteed by a property of the middle matrix called positive semidefiniteness. Figure 9 shows why a positive-definite form is shaped like a bowl that never dips below zero.

The payoff is that this view survives even when the sum no longer looks like a sum of squares. Written as $\mathbf{Y}'(\mathbf{I} - \mathbf{H})\mathbf{Y}$, SSE does not visibly square anything, yet the middle matrix still guarantees it cannot go negative. The same reasoning applies to a variance, a sum of squares in disguise, so the covariance matrix of the next section can never hold a negative variance.

Formula¶

The three regression sums of squares are all quadratic forms in $\mathbf{Y}$, and their middle matrices carry the ANOVA decomposition.

Proof (nonnegativity). Because $\mathbf{I} - \mathbf{H}$ is symmetric and idempotent, $\mathrm{SSE} = \mathbf{Y}'(\mathbf{I} - \mathbf{H})\mathbf{Y} = \mathbf{Y}'(\mathbf{I} - \mathbf{H})'(\mathbf{I} - \mathbf{H})\mathbf{Y} = \mathbf{e}'\mathbf{e} \ge 0$, so SSE is a genuine sum of squares and can never be negative; the same holds for SSR and SSTO. The additive decomposition of the middle matrices is the matrix statement of the ANOVA identity $\mathrm{SSTO} = \mathrm{SSR} + \mathrm{SSE}$ (3.6 The analysis of variance approach). $\blacksquare$

In words: the sums of squares are quadratic forms, and the idempotent middle matrices both guarantee nonnegativity and, in Chapter 7, hand us the degrees of freedom through their traces.

R and Python¶

6.7 Random vectors, expectation, and covariance matrices¶

Intuition¶

So far the matrices held fixed numbers. But $\mathbf{Y}$ is random: rerun the 21 cities under the same conditions and the sales would come out a little different, because of the errors $\boldsymbol{\varepsilon}$. To describe a vector of random variables we need two things: a vector of means and a table of variances and covariances.

The vector of means is just the expectation applied entry by entry. The table is the covariance matrix (Definition 6.19): its diagonal holds the variance of each entry, and its off-diagonals hold the covariance between pairs of entries. For the regression errors, this table is simple: constant variance $\sigma^2$ down the diagonal (every error has the same spread) and zeros off it (distinct errors are uncorrelated), so $\operatorname{Cov}\{\boldsymbol{\varepsilon}\} = \sigma^2\mathbf{I}$. A two-line rule for how means and covariances travel through a linear map then yields the covariance matrix of $\mathbf{b}$, the source of every standard error.

Formula¶

The transformation rules for how these travel through a linear map are the heart of the section.

In words: expectation passes straight through a linear map, and a covariance matrix gets sandwiched between $\mathbf{A}$ and its transpose. Apply these to the regression model, where $E\{\mathbf{Y}\} = \mathbf{X}\boldsymbol{\beta}$ and $\operatorname{Cov}\{\mathbf{Y}\} = \sigma^2\mathbf{I}$, and the mean and covariance of $\mathbf{b}$ follow.

Proof. The estimator is a linear map of $\mathbf{Y}$: with $\mathbf{A} = (\mathbf{X}'\mathbf{X})^{-1}\mathbf{X}'$, we have $\mathbf{b} = \mathbf{A}\mathbf{Y}$. For the mean, use $E\{\mathbf{Y}\} = \mathbf{X}\boldsymbol{\beta}$:

so $\mathbf{b}$ is unbiased for $\boldsymbol{\beta}$. For the covariance, use $\operatorname{Cov}\{\mathbf{Y}\} = \sigma^2\mathbf{I}$ and the sandwich rule:

where the middle $\mathbf{X}'\mathbf{X}$ cancels one inverse. So $\operatorname{Cov}\{\mathbf{b}\} = \sigma^2 (\mathbf{X}'\mathbf{X})^{-1}$: the same inverse we used to compute $\mathbf{b}$ also gives, once scaled by $\sigma^2$, the variance of every coefficient and the covariance of every pair. $\blacksquare$

Replacing the unknown $\sigma^2$ by its estimate MSE gives the estimated covariance matrix $s^2\{\mathbf{b}\} = \mathrm{MSE}\,(\mathbf{X}'\mathbf{X})^{-1}$, whose diagonal square roots are the standard errors $s\{b_k\}$ that software prints.

R and Python¶

6.8 The multivariate normal, in brief¶

Intuition¶

The mean vector and covariance matrix describe a random vector’s center and spread, but not its full shape. For inference we add one assumption, the same one that turned Chapter 2’s estimates into $t$ tests: the errors are normal. Stacked into a vector, normal errors with constant variance and no correlation follow a multivariate normal distribution (Definition 6.22), the several-variable generalization of the bell curve. Its contours are ellipses whose tilt and stretch are read directly off the covariance matrix, as in Figure 10.

What the full distribution adds is how the probability thins out as you move from the center: the same bell-curve falloff you know from one variable, measured along the tilted axes of the covariance ellipse. That is what lets us compute the chance that a coefficient lands within a stated distance of the truth, exactly what a confidence interval reports in the next chapter.

Why this matters: under normal errors the response is multivariate normal, and because $\mathbf{b}$ is a linear map of $\mathbf{Y}$, so is $\mathbf{b}$. That fact, $\mathbf{b} \sim N(\boldsymbol{\beta}, \sigma^2(\mathbf{X}'\mathbf{X})^{-1})$, is the foundation of every confidence interval, $t$ test, and $F$ test ahead.

One caution keeps this honest. Multivariate normality of $\mathbf{b}$ is a consequence of assuming normal errors, not a fact the data hand you for free. If the errors are badly non-normal and the sample is small, that law is only approximate, one reason Chapter 5’s permutation and bootstrap methods exist (5.4 The bootstrap for regression). With a large sample, an averaging effect pulls $\mathbf{b}$ toward normal even when the errors are not.

Formula¶

Applying that property to regression under the normal error model gives the sampling law that drives every later inference.

The response is centered at the regression surface $\mathbf{X}\boldsymbol{\beta}$ with independent, equal-variance coordinates, and the estimator $\mathbf{b}$ is centered at the truth $\boldsymbol{\beta}$ (unbiased) with the covariance matrix we derived.

In words: normal errors make both the data and the estimates normal, so every coefficient is normally distributed with a variance we can read from $(\mathbf{X}'\mathbf{X})^{-1}$, and that is what makes exact inference possible.

R and Python¶

We can demonstrate the covariance-recovery idea by simulation: draw many vectors with a target covariance and confirm the sample covariance matches. Starting from independent standard normals $\mathbf{Z}$, the transformation $\mathbf{V} = \mathbf{Z}\mathbf{L}'$, where $\mathbf{L}$ is a matrix square root of the target $\boldsymbol{\Sigma}$ (the Cholesky factor), produces vectors with covariance $\boldsymbol{\Sigma}$, exactly the sandwich rule at work.

set.seed(4210)
mu <- c(0, 0)
Sigma <- matrix(c(1, 0.8,
                  0.8, 1), nrow = 2)
L <- chol(Sigma)
Z <- matrix(rnorm(2 * 5000), ncol = 2)
V <- Z %*% L
round(cov(V), 3)

      [,1]  [,2]
[1,] 0.967 0.760
[2,] 0.760 0.955

rng = np.random.default_rng(4210)
Sigma = np.array([[1.0, 0.8],
                  [0.8, 1.0]])
L = np.linalg.cholesky(Sigma)
Z = rng.standard_normal((5000, 2))
V = Z @ L.T
print(np.round(np.cov(V, rowvar=False), 3))

[[1.017 0.822]
 [0.822 1.026]]

The sample covariance of the 5000 simulated vectors comes out close to the target $\left(\begin{smallmatrix} 1 & 0.8 \\ 0.8 & 1 \end{smallmatrix}\right)$ in both languages, differing only by sampling noise.

The Cholesky factor can stay a black box for now; the point is that a covariance can be planted on purpose, by passing uncorrelated noise through a fixed matrix and letting the sandwich rule do the rest. The covariance of $\mathbf{b}$ is the same story: error noise $\sigma^2\mathbf{I}$ passed through the fixed matrix $(\mathbf{X}'\mathbf{X})^{-1}\mathbf{X}'$, so the simulation and the regression are two uses of one mechanism.

6.9 Chapter summary¶

You can now speak the matrix language the rest of the book is written in: lay a dataset out as a response vector $\mathbf{Y}$ and a design matrix $\mathbf{X}$, build the products $\mathbf{X}'\mathbf{X}$ and $\mathbf{X}'\mathbf{Y}$ every fit runs on, tell when a matrix is invertible, solve the normal equations for $\mathbf{b} = (\mathbf{X}'\mathbf{X})^{-1}\mathbf{X}'\mathbf{Y}$, build the hat matrix and prove it is a projection, express the sums of squares as quadratic forms, and find the mean and covariance of a random vector. For the real Dwaine data, all of this was computed in both R and numpy: $\mathbf{b} = (-68.86,\ 1.4546,\ 9.3655)'$, $\mathrm{SSE} = 2180.9$ on 18 degrees of freedom, $\mathrm{MSE} = 121.2$, and standard errors $(60.0,\ 0.212,\ 4.06)$. From these tools came the two facts that drive all of multiple regression: $\mathbf{b}$ is unbiased, and its covariance is $\sigma^2(\mathbf{X}'\mathbf{X})^{-1}$.

Every piece fits into one pipeline, drawn in Figure 11: from the data you form two summaries, invert one of them, and everything else, the coefficients, the hat matrix, the residuals, and the standard errors, drops out in a fixed order. If you remember the shape of that pipeline, you can rebuild any single formula in it.

Key results at a glance.

Key terms. matrix, vector, design matrix, transpose, matrix multiplication, conformable, identity matrix, symmetric matrix, linear independence, rank, determinant, inverse, normal equations, hat matrix, idempotent matrix, projection matrix, trace, quadratic form, positive definite, random vector, covariance matrix, multivariate normal.

You should now be able to.

• Represent a regression dataset as a response vector $\mathbf{Y}$ and a design matrix $\mathbf{X}$, and write the model as $\mathbf{Y} = \mathbf{X}\boldsymbol{\beta} + \boldsymbol{\varepsilon}$.

• Compute transposes and matrix products, and build $\mathbf{X}'\mathbf{X}$ and $\mathbf{X}'\mathbf{Y}$ by hand and with software.

• Decide when a square matrix has an inverse using rank and the determinant, and explain what perfect collinearity does to $\mathbf{X}'\mathbf{X}$.

• Solve the normal equations for $\mathbf{b} = (\mathbf{X}'\mathbf{X})^{-1}\mathbf{X}'\mathbf{Y}$.

• Verify that the hat matrix $\mathbf{H}$ is symmetric and idempotent, and explain why fitted values and residuals are projections.

• Express SSE, SSR, and SSTO as quadratic forms, and connect positive definiteness to nonnegative sums of squares.

• Derive the mean vector and covariance matrix of a random vector, and use them to show $\operatorname{Cov}\{\mathbf{b}\} = \sigma^2(\mathbf{X}'\mathbf{X})^{-1}$.

• State the multivariate normal model for $\mathbf{Y}$ and $\mathbf{b}$ and describe how it powers the inference of later chapters.

Where this fits. This chapter is the toolbox that makes the FIT and USE stages of the course workflow (The modeling workflow) work for more than one predictor: the estimator $\mathbf{b} = (\mathbf{X}'\mathbf{X})^{-1}\mathbf{X}'\mathbf{Y}$ fits every multiple regression, and the covariance matrix $\sigma^2(\mathbf{X}'\mathbf{X})^{-1}$ is how every coefficient is later tested. Chapter 2 built these ideas for one predictor with scalar algebra (2.2 Least squares from first principles); we have now generalized them to any number. Next, Chapter 7 takes the same Dwaine $\mathbf{X}'\mathbf{X}$ you computed and proves why these formulas are the right ones: it derives least squares two ways (7.1 The model and least squares in matrix form), reads the degrees of freedom off the hat matrix (7.3 The hat matrix), and proves the Gauss-Markov theorem (7.6 The Gauss-Markov theorem). It is the most abstract chapter in the course, and it is built almost entirely from the matrix facts you assembled here, so the work of this chapter is exactly what makes the next one readable.

6.10 Frequently asked questions¶

Q1. Why do we glue a column of ones onto $\mathbf{X}$? So the intercept can be treated as just another coefficient. Multiplying the ones column by $\beta_0$ gives $\beta_0$ in every row, which is exactly what an intercept does. Without the ones column, the formula $\mathbf{X}\boldsymbol{\beta}$ would have no constant term, and you would be forcing the regression surface through the origin.

Q2. Is $\mathbf{X}'\mathbf{X}$ the same as squaring $\mathbf{X}$? No. $\mathbf{X}$ is usually not square, so $\mathbf{X}^2 = \mathbf{X}\mathbf{X}$ is not even defined. $\mathbf{X}'\mathbf{X}$ multiplies the transpose ($p \times n$) by $\mathbf{X}$ ($n \times p$) to make a $p \times p$ matrix of sums of products. It is the closest matrix analogue of “sum of squares,” which is why it sits at the center of least squares.

Q3. What actually goes wrong when $\mathbf{X}'\mathbf{X}$ is singular? Two predictor columns carry the same information, so the data cannot decide how to split an effect between them. Infinitely many coefficient vectors fit equally well, the inverse does not exist, and software either errors out or silently drops a predictor. This is the extreme end of the multicollinearity you meet in 12.1 Multicollinearity and the variance inflation factor.

Q4. Do I have to invert $\mathbf{X}'\mathbf{X}$ by hand? No. Past $2 \times 2$, let solve or np.linalg.inv do it, and in real work prefer the routines inside lm and statsmodels, which solve the normal equations without forming the inverse. Computing the inverse here just shows that the software does exactly the matrix algebra of this chapter, so it is never a black box.

Q5. Why is the hat matrix called a projection? Because it takes the response vector $\mathbf{Y}$ and drops it perpendicularly onto the flat subspace of all fittable values (the column space of $\mathbf{X}$), landing at the closest point, $\widehat{\mathbf{Y}}$. Projecting a second time changes nothing, which is the algebraic property $\mathbf{H}\mathbf{H} = \mathbf{H}$. The residual is the part of $\mathbf{Y}$ left sticking out, at a right angle to the subspace.

Q6. Where does the normal distribution enter? Only in 6.8 The multivariate normal, in brief, and only for inference. Building $\mathbf{X}$, computing $\mathbf{b}$, the hat matrix, SSE, and the covariance matrix all use no distributional assumption beyond mean zero, constant variance, and uncorrelated errors. Normality is added on top so that $\mathbf{b}$ is multivariate normal and we can build exact $t$ and $F$ procedures, which is Chapter 7’s job.

Q7. Do I have to memorize all these identities, and should the negative entries in the covariance matrix of $\mathbf{b}$ worry me? No on both counts. Hold on to three things: the model $\mathbf{Y} = \mathbf{X}\boldsymbol{\beta} + \boldsymbol{\varepsilon}$, the estimator $\mathbf{b} = (\mathbf{X}'\mathbf{X})^{-1}\mathbf{X}'\mathbf{Y}$, and the covariance $\sigma^2(\mathbf{X}'\mathbf{X})^{-1}$; everything else builds one of those or reads a number out of it, so the identities become things you rederive rather than recall. As for the negative entries, only the off-diagonal covariances can go negative (the diagonal variances never do), and a negative one between two slopes just means that across repeated samples overestimating one tends to accompany underestimating the other. It is information, not an error.

6.11 Practice problems¶

1. (A) Give the dimensions of $\mathbf{X}$, $\mathbf{X}'$, $\mathbf{X}'\mathbf{X}$, $\mathbf{X}'\mathbf{Y}$, $\mathbf{b}$, and $\mathbf{H}$ for the Dwaine model, and say in one phrase what each represents.

2. (A) Explain why the first column of the design matrix is all ones, and what would change in the model if it were removed.

3. (A) State the rule for when two matrices can be multiplied, and use it to explain why $\mathbf{X}\mathbf{X}$ is undefined for the $21 \times 3$ Dwaine design matrix but $\mathbf{X}'\mathbf{X}$ is defined.

4. (A) In words, what is the entry in row 1, column 1 of $\mathbf{X}'\mathbf{X}$ equal to, and why? What are the other entries of the first row?

5. (A) Define an idempotent matrix and a symmetric matrix, and state which of these properties the hat matrix $\mathbf{H}$ has.

6. (A) A colleague writes $\mathbf{X}\mathbf{b} = \mathbf{Y}$ and “cancels $\mathbf{X}$” from $\mathbf{X}'\mathbf{X}\mathbf{b} = \mathbf{X}'\mathbf{Y}$. Explain the two things wrong with this.

7. (A) Explain what it means for the columns of $\mathbf{X}$ to be linearly dependent, and what that does to $\mathbf{X}'\mathbf{X}$ and to the least-squares fit.

8. (A) The covariance matrix of $\mathbf{b}$ has a negative off-diagonal entry between the two slopes. Interpret its sign in one sentence.

9. (B) Prove that $\mathbf{X}'\mathbf{X}$ is symmetric for any matrix $\mathbf{X}$, citing the transpose-of-a-product rule (Theorem 6.5).

10. (B) Prove the transpose-of-a-product rule (Theorem 6.5), $(\mathbf{A}\mathbf{B})' = \mathbf{B}'\mathbf{A}'$, by comparing the $(i,j)$ entries of both sides.

11. (B) Starting from the matrix normal equations $\mathbf{X}'\mathbf{X}\,\mathbf{b} = \mathbf{X}'\mathbf{Y}$, derive $\mathbf{b} = (\mathbf{X}'\mathbf{X})^{-1}\mathbf{X}'\mathbf{Y}$ (Theorem 6.12), stating the condition on $\mathbf{X}$ that the step requires.

12. (B) Prove that the hat matrix $\mathbf{H} = \mathbf{X}(\mathbf{X}'\mathbf{X})^{-1}\mathbf{X}'$ is symmetric and idempotent (Theorem 6.16). Then show $\mathbf{I} - \mathbf{H}$ is idempotent.

13. (B) Show that $\mathbf{H}\mathbf{X} = \mathbf{X}$, and explain geometrically why projecting each column of $\mathbf{X}$ onto the column space of $\mathbf{X}$ leaves it unchanged.

14. (B) Using $\mathbf{e} = (\mathbf{I} - \mathbf{H})\mathbf{Y}$ and $\widehat{\mathbf{Y}} = \mathbf{H}\mathbf{Y}$, prove that $\mathbf{X}'\mathbf{e} = \mathbf{0}$ and that $\widehat{\mathbf{Y}}'\mathbf{e} = 0$ (fitted values are orthogonal to residuals).

15. (B) Prove the covariance transformation rule $\operatorname{Cov}\{\mathbf{A}\mathbf{W}\} = \mathbf{A}\operatorname{Cov}\{\mathbf{W}\}\mathbf{A}'$ (part of Theorem 6.20) for a fixed matrix $\mathbf{A}$, starting from the definition $\operatorname{Cov}\{\mathbf{W}\} = E\{(\mathbf{W} - \boldsymbol{\mu})(\mathbf{W} - \boldsymbol{\mu})'\}$.

16. (B) Use the rules of problem 15 and $E\{\mathbf{Y}\} = \mathbf{X}\boldsymbol{\beta}$, $\operatorname{Cov}\{\mathbf{Y}\} = \sigma^2\mathbf{I}$ to derive $E\{\mathbf{b}\} = \boldsymbol{\beta}$ and $\operatorname{Cov}\{\mathbf{b}\} = \sigma^2(\mathbf{X}'\mathbf{X})^{-1}$ (Theorem 6.21).

17. (B) Show that $\mathrm{SSE} = \mathbf{Y}'(\mathbf{I} - \mathbf{H})\mathbf{Y}$ equals $\mathbf{e}'\mathbf{e}$ (Theorem 6.18), using the symmetry and idempotency of $\mathbf{I} - \mathbf{H}$, and explain why this proves $\mathrm{SSE} \ge 0$.

18. (B) For the $2 \times 2$ matrix $\mathbf{A} = \left(\begin{smallmatrix} a & b \\ b & d \end{smallmatrix}\right)$, write the quadratic form $\mathbf{u}'\mathbf{A}\mathbf{u}$ in terms of $u_1, u_2$, and give a condition on $a, b, d$ that makes it positive definite.

19. (B) The fitted values satisfy $\widehat{\mathbf{Y}} = \mathbf{H}\mathbf{Y}$. Prove that $\operatorname{Cov}\{\widehat{\mathbf{Y}}\} = \sigma^2\mathbf{H}$ and $\operatorname{Cov}\{\mathbf{e}\} = \sigma^2(\mathbf{I} - \mathbf{H})$, using the covariance rule and the properties of $\mathbf{H}$.

20. (C) Read dwaine.csv, build $\mathbf{X}$ and $\mathbf{Y}$, and reproduce $\mathbf{X}'\mathbf{X}$ and $\mathbf{X}'\mathbf{Y}$ in R or Python. Confirm the top-left entry of $\mathbf{X}'\mathbf{X}$ is $n$ and the first entry of $\mathbf{X}'\mathbf{Y}$ is $\sum Y$.

21. (C) Compute $(\mathbf{X}'\mathbf{X})^{-1}$ and $\mathbf{b} = (\mathbf{X}'\mathbf{X})^{-1}\mathbf{X}'\mathbf{Y}$ in software, and check that $\mathbf{b}$ matches lm/statsmodels and that $(\mathbf{X}'\mathbf{X})^{-1}\mathbf{X}'\mathbf{X}$ returns the identity.

22. (C) Build the hat matrix $\mathbf{H}$, verify numerically that $\mathbf{H}\mathbf{H} = \mathbf{H}$ and that $\operatorname{tr}(\mathbf{H}) = 3$, and report the three largest diagonal entries $h_{ii}$ (the leverage values).

23. (C) Compute the residual vector $\mathbf{e} = (\mathbf{I} - \mathbf{H})\mathbf{Y}$ and verify numerically that $\mathbf{X}'\mathbf{e} = \mathbf{0}$ (all three entries near zero) and that $\sum e_i = 0$.

24. (C) Compute $\mathrm{SSE} = \mathbf{e}'\mathbf{e}$, then $\mathrm{MSE}$ and $s = \sqrt{\mathrm{MSE}}$, and confirm they match the residual standard error reported by summary(fit) / fit.summary().

25. (C) Form the estimated covariance matrix $\mathrm{MSE}\,(\mathbf{X}'\mathbf{X})^{-1}$ and report the three standard errors from its diagonal. Confirm they match the software summary, and report the estimated covariance between the two slope estimates.

26. (C) Predict sales for a new city with $\text{targtpop} = 65.4$ and $\text{dispoinc} = 17.6$ by forming the row vector $\mathbf{x}_h = (1, 65.4, 17.6)$ and computing $\mathbf{x}_h \mathbf{b}$. Confirm the result against predict.

27. (C) Add a redundant column to $\mathbf{X}$ equal to $\text{targtpop} + \text{dispoinc}$, and try to compute $(\mathbf{X}'\mathbf{X})^{-1}$. Report what R or Python does (an error, a warning, or a wildly unstable inverse), and connect it to rank deficiency.

28. (C) Center the predictors: replace targtpop and dispoinc by their deviations from their means, refit, and confirm that the two slopes are unchanged while the intercept becomes $\bar{Y}$. Explain, using the normal equations, why centering leaves the slopes alone.

29. (C) Demonstrate the sampling law $\mathbf{b} \sim N(\boldsymbol{\beta}, \sigma^2(\mathbf{X}'\mathbf{X})^{-1})$ (Theorem 6.23, 6.8 The multivariate normal, in brief) by simulation. Treat the fitted $\mathbf{b} = (-68.86,\ 1.4546,\ 9.3655)'$ as the true $\boldsymbol{\beta}$ and $\mathrm{MSE} = 121.16$ as the true $\sigma^2$. Holding the real Dwaine $\mathbf{X}$ fixed, use set.seed(4210) (R) or default_rng(4210) (Python) to generate 5000 response vectors $\mathbf{Y} = \mathbf{X}\boldsymbol{\beta} + \boldsymbol{\varepsilon}$ with $\boldsymbol{\varepsilon} \sim N(\mathbf{0}, \sigma^2\mathbf{I})$, refit each by $(\mathbf{X}'\mathbf{X})^{-1}\mathbf{X}'\mathbf{Y}$, and collect the estimates. (a) Confirm the column means are close to $\boldsymbol{\beta}$ and the sample covariance is close to $\sigma^2(\mathbf{X}'\mathbf{X})^{-1}$. (b) In two sentences, explain why this law lets Chapter 7 attach a $t$ distribution to each coefficient and turn a standard error into a confidence interval.

6.12 Exam practice¶

These five questions match the style of the course exams: each asks you to explain, evaluate, or interpret in full sentences, not to produce a bare number. Write in complete sentences and say which numbers you used. The software output shown was generated from the real dwaine.csv fit. Each model answer shows the depth that earns full marks, followed by one line on what a weak answer misses.

EP 6.1. In the Dwaine fit both estimated slopes are positive: $b_{\text{targtpop}} = 1.4546$ and $b_{\text{dispoinc}} = 9.3655$. A student argues, “since both predictors push sales up, the two slope estimates must be positively correlated across repeated samples.” The estimated covariance matrix of $\mathbf{b}$ is printed below. Evaluate the student’s claim, and explain what the relevant number actually tells you.

          intercept targtpop  dispoinc
intercept 3602.0347   8.7459 -241.4230
targtpop     8.7459    0.0449   -0.6724
dispoinc  -241.4230   -0.6724   16.5158

EP 6.2. Suppose targtpop is re-entered in raw persons instead of thousands, so every value is multiplied by 1000, and the model is refit. Using the output below, state precisely which quantities change and which stay identical, and explain why.

          b            se           t
intercept -68.8570732  60.0169532  -1.1473
targtpop    0.0014546   0.0002118   6.8682
dispoinc    9.3655004   4.0639581   2.3045
SSE = 2180.9274      fitted[1:3] = 187.184, 154.229, 234.396

EP 6.3. The diagonal of the hat matrix holds the leverage values $h_{ii}$. The output below reports their sum, the $2p/n$ rule-of-thumb cutoff, and the five largest values with the two predictor values of each city. Interpret these numbers in context: what do they say about which cities matter most, why do they sum to 3, and would you flag any city as a high-leverage point?

sum of h_ii = 3.0        2p/n cutoff = 0.2857
 city 20:  targtpop 82.7  dispoinc 19.1   h = 0.2788
 city 13:  targtpop 88.4  dispoinc 17.4   h = 0.2390
 city 15:  targtpop 52.5  dispoinc 17.8   h = 0.2095
 city  3:  targtpop 91.3  dispoinc 18.2   h = 0.1737
 city  5:  targtpop 46.9  dispoinc 17.3   h = 0.1620

EP 6.4. A student augments the Dwaine design matrix with a fourth column equal to targtpop + dispoinc, refits, and reports: “Python still returned an inverse and some coefficients, so the augmented model is fine.” The diagnostics below come from the augmented matrix. Explain why the augmented model is not fine, what has gone wrong mathematically, and why R and Python behave differently.

rank of augmented X = 3   (it has 4 columns)
det(X'X)        = -3.7e-05      # against ~1.07e06 for the genuine 3x3 X'X
R  solve(X'X):  Error: system is computationally singular
Python inv(X'X): returns a matrix of enormous (~1e13) entries, no error

EP 6.5. A regional manager wants the model’s predicted mean sales for a new city with targtpop $= 65.4$ and dispoinc $= 17.6$, together with a standard error, so form $\mathbf{x}_h = (1,\ 65.4,\ 17.6)$ and note that the predicted mean is $\widehat{Y}_h = \mathbf{x}_h'\mathbf{b}$. The estimated covariance matrix $s^2\{\mathbf{b}\} = \mathrm{MSE}\,(\mathbf{X}'\mathbf{X})^{-1}$ is reprinted below, and the fit gives $\widehat{Y}_h = 191.10$. A student computes the standard error of $\widehat{Y}_h$ as $\sqrt{65.4^2 \cdot s^2\{b_{\text{targtpop}}\} + 17.6^2 \cdot s^2\{b_{\text{dispoinc}}\} + s^2\{b_0\}} = 94.4$, using only the diagonal variances. Explain why this is wrong, give the correct standard error, and say which way the diagonal-only shortcut errs.

          intercept targtpop  dispoinc
intercept 3602.0347   8.7459 -241.4230
targtpop     8.7459    0.0449   -0.6724
dispoinc  -241.4230   -0.6724   16.5158

x_h'(X'X)^{-1} x_h = 0.063181       MSE = 121.16

6.13 Chapter game¶