7. The general linear model

We return to the Dwaine Studios cities of 6. Matrix algebra for regression. The company runs a chain of portrait studios in mid-sized cities, the kind of place families visit for graduation and holiday photos, and it wants a rule for guessing how much a new studio will sell before it signs a lease. Two facts about a city are easy to look up in advance: the number of young people (population aged 16 or younger, in thousands) and the per-capita disposable income (in thousands of dollars). The question is whether those two numbers together predict annual studio sales well enough to guide the expansion.

In Chapter 6 you already assembled this fit by hand: you built the design matrix, formed $\mathbf{X}'\mathbf{X}$ and $\mathbf{X}'\mathbf{Y}$, inverted the small matrix, and solved for the coefficients $\mathbf{b}$. What that arithmetic could not tell you is why the recipe $\mathbf{b} = (\mathbf{X}'\mathbf{X})^{-1}\mathbf{X}'\mathbf{Y}$ is the right one, or how far to trust the numbers it returns. This chapter proves why, and turns the fit into inference. Figure 1 shows the raw picture: sales climb with each predictor on its own, but neither plot alone tells you what the two predictors do together, and that “together” is the whole point of multiple regression.

Simple linear regression, from Chapters 2 and 3, handles one predictor. Real questions almost always come with several. We could try to rebuild the whole theory predictor by predictor, but that path gets ugly fast: two predictors, then three, each with its own normal equations and variance formulas. Instead we do what 6. Matrix algebra for regression set up. We write the model once, in matrix form, so that one predictor or ten predictors look identical on the page. Every result from Chapter 2, the least-squares estimates, their unbiasedness, their variances, the Gauss-Markov optimality, the $n-2$ divisor, reappears here in a single clean form that works for any number of predictors at once. This chapter is where regression stops being a collection of special cases and becomes one general method.

This is the most abstract chapter in the course, and worth a word on how to read it. The geometry in 7.2 The geometry of least squares is the idea to hold onto: once you picture fitting as dropping a perpendicular from the data onto the space of possible fits, every proof that follows reduces to one of two facts about the hat matrix. Hold that picture and the algebra becomes bookkeeping.

7.1 The model and least squares in matrix form¶

Intuition¶

In Chapter 2 a single predictor gave one slope and one intercept. Add a second predictor and you have two slopes plus an intercept; the mean response is now a tilted plane instead of a line, as Figure 2 shows for Dwaine. With more predictors it is a flat sheet in a space you cannot picture. Writing separate symbols for every coefficient and separate sums for every normal equation would bury the idea under bookkeeping. The fix is to stack the observations into vectors and matrices and let one line of matrix algebra carry all the bookkeeping for us.

The systematic part of each observation is a weighted sum of that city’s predictor values, the same weights ($\beta_0, \beta_1, \dots$) reused for every city. Stacking the cities vertically turns “the same weights applied to every row” into a single matrix-times-vector product. Figure 3 shows the four objects and their shapes.

Formula¶

The general linear model (Definition 7.1) collects all $n$ observations into one matrix equation.

In words: each observation’s response equals a weighted sum of that observation’s predictor values plus a random error, and the same weights are reused for every row. The four objects are:

• $\mathbf{Y}$ is the $n \times 1$ vector of responses, one entry per observation.

• $\mathbf{X}$ is the $n \times p$ design matrix (Definition 7.2). Its first column is all ones (for the intercept); each remaining column holds one predictor’s values. Row $i$ is city $i$’s data.

• $\boldsymbol{\beta}$ is the $p \times 1$ vector of unknown parameters $(\beta_0, \beta_1, \dots, \beta_{p-1})'$.

• $\boldsymbol{\varepsilon}$ is the $n \times 1$ vector of random errors.

• $p$ is the number of parameters, counting the intercept; with $k$ predictors, $p = k + 1$. Here $n = 21$ and $p = 3$.

The error assumptions of Chapter 2 become two compact statements about the whole vector at once:

In words: every error averages to zero, and the covariance matrix is $\sigma^2$ times the identity. That one symbol packs both assumptions at once: each error has the same variance $\sigma^2$ (down the diagonal), and distinct errors are uncorrelated (off the diagonal). Because $\mathbf{X}$ and $\boldsymbol{\beta}$ are fixed, $E\{\mathbf{Y}\} = \mathbf{X}\boldsymbol{\beta}$ and $\operatorname{Cov}\{\mathbf{Y}\} = \sigma^2 \mathbf{I}$ as well. We assume throughout that $\mathbf{X}$ has full column rank $p$: no predictor is a copy of another, or an exact linear combination of the others. That is exactly what makes $\mathbf{X}'\mathbf{X}$ invertible (6.4 The inverse and the normal equations).

To fit the model we reuse the least-squares idea from 2.2 Least squares from first principles: pick the coefficient vector $\mathbf{b}$ that makes the fitted values closest to the data in total squared error. Writing a candidate vector as $\mathbf{b}$, the sum of squared errors is

where $\mathbf{x}_i'$ is row $i$ of $\mathbf{X}$. In words: $Q$ is the squared length of the residual vector, the same quantity Chapter 2 minimized, now written as one dot product. The minimizer is the least-squares estimator (Theorem 7.3)

In words: gather the cross-products of the predictors ($\mathbf{X}'\mathbf{X}$) and of the predictors with the response ($\mathbf{X}'\mathbf{Y}$), and multiply by the inverse of the first. We derive this next, then check it on Dwaine.

Derivation by calculus¶

Proof. Expand $Q$ using the transpose rules of 6.2 Transpose and matrix multiplication. Since $\mathbf{b}'\mathbf{X}'\mathbf{Y}$ is a scalar, it equals its own transpose $\mathbf{Y}'\mathbf{X}\mathbf{b}$, so the two middle terms combine:

We minimize by taking the gradient with respect to $\mathbf{b}$ and setting it to zero. The gradient is just the vector of partial derivatives, one per entry of $\mathbf{b}$, and setting it to zero finds the bottom of the bowl. Two vector-derivative rules do the work: for a constant vector $\mathbf{a}$, $\nabla_{\mathbf{b}}(\mathbf{a}'\mathbf{b}) = \mathbf{a}$; and for a symmetric matrix $\mathbf{M}$, $\nabla_{\mathbf{b}}(\mathbf{b}'\mathbf{M}\mathbf{b}) = 2\mathbf{M}\mathbf{b}$. These are the vector twins of the scalar rules $\frac{d}{db}(ab) = a$ and $\frac{d}{db}(mb^2) = 2mb$ you already know. The matrix $\mathbf{X}'\mathbf{X}$ is symmetric, so

Setting the gradient to zero and dividing by 2 gives the normal equations (Definition 7.4)

Because $\mathbf{X}$ has full column rank, $\mathbf{X}'\mathbf{X}$ is invertible (6.4 The inverse and the normal equations), so we may multiply on the left by $(\mathbf{X}'\mathbf{X})^{-1}$:

This stationary point is the global minimum, not a saddle or maximum, because the matrix of second derivatives of $Q$ is $2\mathbf{X}'\mathbf{X}$, which is positive definite; a quadratic with a positive-definite Hessian is strictly convex and has a single lowest point. $\blacksquare$

These are the same normal equations as Chapter 2, now stated once for any number of predictors. In 2.2 Least squares from first principles the two scalar normal equations were $\sum e_i = 0$ and $\sum X_i e_i = 0$; the single matrix equation $\mathbf{X}'\mathbf{X}\mathbf{b} = \mathbf{X}'\mathbf{Y}$ contains both of those and one more for each extra predictor, as we will see in 7.7 Simple regression as a special case.

R¶

Chapter 6 built $\mathbf{X}'\mathbf{X}$ and $\mathbf{X}'\mathbf{Y}$ for Dwaine by hand. Here we let the computer form them and solve.

dwaine <- read.csv("data/dwaine.csv")
head(dwaine, 3)
X <- cbind(1, dwaine$targtpop, dwaine$dispoinc)
Y <- dwaine$sales
dim(X)

  targtpop dispoinc sales
1     68.5     16.7 174.4
2     45.2     16.8 164.4
3     91.3     18.2 244.2
[1] 21  3

7.2 The geometry of least squares¶

Intuition¶

The calculus derivation finds the answer but hides why it is the answer. There is a picture that makes least squares obvious, and it is worth carrying for the rest of the course. Think of the response $\mathbf{Y}$ as a single arrow (a point) in $n$-dimensional space, one axis per observation. As the coefficient vector $\mathbf{b}$ ranges over all possible values, the fitted vector $\hat{\mathbf{Y}} = \mathbf{X}\mathbf{b}$ traces out a flat subspace: every linear combination of the columns of $\mathbf{X}$. This subspace is the column space of $\mathbf{X}$ (Definition 7.5), and it holds every prediction the model can possibly make.

The data point $\mathbf{Y}$ almost never sits exactly in that subspace, because real data has noise. So the fitting problem is: which point in the column space is closest to $\mathbf{Y}$? The closest point on a flat surface to an outside point is the foot of the perpendicular dropped from the point to the surface. That foot is $\hat{\mathbf{Y}}$, and the perpendicular segment is the residual vector $\mathbf{e} = \mathbf{Y} - \hat{\mathbf{Y}}$. Figure 4 draws exactly this for the smallest case that fits on a page, three observations and two columns.

Formula¶

The geometric condition is that the residual is orthogonal to every column of $\mathbf{X}$:

In words: each column of $\mathbf{X}$, dotted with the residual vector, gives zero, so the residual points in a direction the model cannot represent. Expanding $\mathbf{X}'(\mathbf{Y} - \mathbf{X}\mathbf{b}) = \mathbf{0}$ gives back the normal equations $\mathbf{X}'\mathbf{X}\mathbf{b} = \mathbf{X}'\mathbf{Y}$, so the geometry and the calculus lead to the identical solution.

Derivation by projection¶

Derivation (least squares as orthogonal projection). Let $\mathbf{b}$ solve the normal equations, so $\mathbf{X}'\mathbf{e} = \mathbf{0}$ with $\mathbf{e} = \mathbf{Y} - \mathbf{X}\mathbf{b}$. Take any other candidate coefficient vector $\mathbf{b}^{\ast}$, with fitted vector $\mathbf{X}\mathbf{b}^{\ast}$. Split the gap from $\mathbf{Y}$ to that candidate by routing through $\hat{\mathbf{Y}} = \mathbf{X}\mathbf{b}$:

The second piece $\mathbf{g} = \mathbf{X}(\mathbf{b} - \mathbf{b}^{\ast})$ lives in the column space. Its squared length adds without a cross term, because the cross term vanishes by orthogonality:

Therefore, by the Pythagorean theorem in $n$ dimensions,

with equality only when $\|\mathbf{g}\|^2 = 0$, that is $\mathbf{X}\mathbf{b}^{\ast} = \mathbf{X}\mathbf{b}$. Since $\mathbf{X}$ has full column rank, this forces $\mathbf{b}^{\ast} = \mathbf{b}$. So the normal-equations solution gives the unique smallest sum of squared errors, and $\hat{\mathbf{Y}}$ is the closest point of the column space to $\mathbf{Y}$. $\blacksquare$

This is the same algebraic identity as the Chapter 2 proof in 2.2 Least squares from first principles, now read as geometry: the cross term vanishes there because $\sum e_i = 0$ and $\sum X_i e_i = 0$, which is precisely $\mathbf{X}'\mathbf{e} = \mathbf{0}$ written out. The projection view (in 6.5 Idempotent and projection matrices) is the one to remember, because diagnostics in Chapter 9 are all about which observations sit in awkward corners of the column space.

R and Python¶

Orthogonality is easy to state and easy to doubt, so the simulation below hands you the data and invites you to break it.

7.3 The hat matrix¶

Intuition¶

Section 7.2 pictured the fit as a projection of $\mathbf{Y}$ onto the column space. That projection is carried out by a single matrix, and giving it a name pays off for the rest of the chapter. The fitted vector is a fixed linear function of the data: plug the formula for $\mathbf{b}$ into $\hat{\mathbf{Y}} = \mathbf{X}\mathbf{b}$ and you get $\hat{\mathbf{Y}} = \mathbf{X}(\mathbf{X}'\mathbf{X})^{-1}\mathbf{X}'\mathbf{Y}$. The matrix sitting in front of $\mathbf{Y}$ turns the raw responses into fitted values, so it earns the name hat matrix (Definition 7.6): it puts the hat on $\mathbf{Y}$. It is the single most useful matrix in the rest of the course. Its diagonal entries are the leverages that flag unusual predictor values (Chapter 9), and its algebra makes the variance results below fall out in a line.

Formula¶

The hat matrix and the fitted and residual vectors are

• $\mathbf{H}$ is $n \times n$ and depends only on the predictors, not on $\mathbf{Y}$.

• $\hat{\mathbf{Y}} = \mathbf{H}\mathbf{Y}$: each fitted value is a fixed weighted average of all the responses, with weights read off a row of $\mathbf{H}$.

• $\mathbf{e} = (\mathbf{I} - \mathbf{H})\mathbf{Y}$: the residual is what is left after projecting, so $\mathbf{I} - \mathbf{H}$ is itself a projection, onto the space perpendicular to the columns of $\mathbf{X}$.

Figure 7 shows the Toluca hat matrix as a grid of weights; the diagonal entries, outlined, are the leverages $h_{ii}$.

Derivation (properties of H)¶

Proof. Write $\mathbf{A} = (\mathbf{X}'\mathbf{X})^{-1}$, which is symmetric because $\mathbf{X}'\mathbf{X}$ is (the inverse of a symmetric matrix is symmetric, by 6.4 The inverse and the normal equations).

Symmetric. Using $(\mathbf{BCD})' = \mathbf{D}'\mathbf{C}'\mathbf{B}'$ from 6.2 Transpose and matrix multiplication,

Idempotent. Projecting twice equals projecting once. The middle $\mathbf{X}'\mathbf{X}$ cancels its inverse:

Figure 8 shows why: once $\mathbf{H}$ has landed a vector in the column space, projecting again cannot move it.

Trace $p$. The trace of a matrix is the sum of its diagonal entries. Use the cyclic property of the trace, $\operatorname{tr}(\mathbf{B}\mathbf{C}) = \operatorname{tr}(\mathbf{C}\mathbf{B})$ (a standard identity: both sides equal $\sum_{i,j} B_{ij} C_{ji}$), to slide $\mathbf{X}$ around:

Consequently $\mathbf{I} - \mathbf{H}$ is also symmetric and idempotent, with $\operatorname{tr}(\mathbf{I} - \mathbf{H}) = n - p$. That trace is the degrees of freedom for error, and it is the reason the divisor in the variance estimate is $n - p$, proved in 7.5 Estimating the error variance. $\blacksquare$

R and Python¶

The claim that $\mathbf{H}$ depends only on the predictors sounds like a technicality until you try to move it with the response, which is what the next simulation asks you to do.

7.4 The sampling distribution of b¶

Intuition¶

The vector $\mathbf{b}$ came from one sample of 21 cities. A different 21 cities would give a slightly different $\mathbf{b}$, so $\mathbf{b}$ is random, with a mean and a covariance across the samples we could have drawn. Two questions repeat from Chapter 2, now for a whole vector: is $\mathbf{b}$ right on average, and how much does it vary? The answers are clean because $\mathbf{b}$ is a fixed matrix times $\mathbf{Y}$, and 6.7 Random vectors, expectation, and covariance matrices tells us exactly how means and covariances pass through a linear map.

Formula¶

The least-squares vector is unbiased with a covariance matrix built from $(\mathbf{X}'\mathbf{X})^{-1}$:

In words: on average $\mathbf{b}$ hits the true coefficient vector, and its covariance matrix is $\sigma^2$ times the inverse cross-product matrix. The diagonal entries are the variances of the individual coefficients; the off-diagonal entries are the covariances between them, which is why two coefficients can be correlated (the tilted cloud in Figure 10). Replacing $\sigma^2$ by its estimate $\mathrm{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 software prints.

Derivation (mean and covariance of b)¶

Proof. Write $\mathbf{b} = \mathbf{A}\mathbf{Y}$ with the fixed matrix $\mathbf{A} = (\mathbf{X}'\mathbf{X})^{-1}\mathbf{X}'$. For a fixed matrix times a random vector, 6.7 Random vectors, expectation, and covariance matrices gives $E\{\mathbf{A}\mathbf{Y}\} = \mathbf{A}\,E\{\mathbf{Y}\}$ and $\operatorname{Cov}\{\mathbf{A}\mathbf{Y}\} = \mathbf{A}\operatorname{Cov}\{\mathbf{Y}\}\mathbf{A}'$.

Mean. Since $E\{\mathbf{Y}\} = \mathbf{X}\boldsymbol{\beta}$,

Covariance. Since $\operatorname{Cov}\{\mathbf{Y}\} = \sigma^2\mathbf{I}$,

where the middle $\mathbf{X}'\mathbf{X}$ cancels one inverse and the symmetry of $(\mathbf{X}'\mathbf{X})^{-1}$ handles the transpose. $\blacksquare$

If we add the normality assumption $\boldsymbol{\varepsilon} \sim N(\mathbf{0}, \sigma^2\mathbf{I})$, then $\mathbf{Y}$ is multivariate normal, and a fixed linear map of a multivariate normal is again multivariate normal (6.8 The multivariate normal, in brief). So

In words: under normal errors each coefficient is normally distributed around its true value, and this is what powers the $t$ tests and confidence intervals of Chapter 8. Figure 10 confirms the covariance by simulation: refitting many datasets scatters the two Dwaine slopes into a tilted ellipse that matches $\sigma^2 (\mathbf{X}'\mathbf{X})^{-1}$.

R and Python¶

The covariance formula is a prediction about data you have not collected, so the simulation below collects it for you a thousand times and checks.

7.5 Estimating the error variance¶

Intuition¶

The coefficients describe the plane; the model still has one more unknown, the error variance $\sigma^2$, which sets how far the data scatter off the plane. As in Chapter 2 we estimate it from the residuals, but now we must find the right divisor for a model with $p$ parameters instead of 2. The answer is $n - p$, and the reason is the trace we computed in 7.3 The hat matrix: fitting $p$ coefficients uses up $p$ directions of the data, leaving $n - p$ directions for the residual to live in. Figure 12 keeps the bookkeeping in view: the $n$ directions of the data split into a fit part and a residual part, and we average the squared residual over its own $n - p$ directions, not all $n$.

Formula¶

The error sum of squares and the estimator of $\sigma^2$ are

• $\mathrm{SSE}$ is the squared length of the residual vector, the variation the model did not explain.

• $\mathrm{MSE}$ divides $\mathrm{SSE}$ by its degrees of freedom $n - p$, and estimates $\sigma^2$.

• $s$ is the estimated error standard deviation, in the units of $Y$.

Derivation (why the divisor is n - p)¶

We need one fact about the expected value of a quadratic form. It is provable at our level, so we prove it, then apply it.

Lemma (expectation of a quadratic form). For any random vector $\mathbf{Y}$ with mean $\boldsymbol{\mu}$ and covariance $\boldsymbol{\Sigma}$, and any fixed symmetric matrix $\mathbf{A}$,

Proof. The scalar $\mathbf{Y}'\mathbf{A}\mathbf{Y}$ equals its own trace, and the trace is cyclic, so $\mathbf{Y}'\mathbf{A}\mathbf{Y} = \operatorname{tr}(\mathbf{A}\mathbf{Y}\mathbf{Y}')$. Taking expectations and using $E\{\mathbf{Y}\mathbf{Y}'\} = \boldsymbol{\Sigma} + \boldsymbol{\mu}\boldsymbol{\mu}'$ (6.7 Random vectors, expectation, and covariance matrices),

since $\operatorname{tr}(\mathbf{A}\boldsymbol{\mu}\boldsymbol{\mu}') = \boldsymbol{\mu}'\mathbf{A}\boldsymbol{\mu}$ is a scalar. $\blacksquare$

Proof of Theorem 7.9. Apply the lemma with $\mathbf{A} = \mathbf{I} - \mathbf{H}$, $\boldsymbol{\mu} = \mathbf{X}\boldsymbol{\beta}$, and $\boldsymbol{\Sigma} = \sigma^2\mathbf{I}$. The two pieces are:

Trace piece. $\operatorname{tr}\!\big((\mathbf{I} - \mathbf{H})\sigma^2\mathbf{I}\big) = \sigma^2 \operatorname{tr}(\mathbf{I} - \mathbf{H}) = \sigma^2 (n - p)$, using the trace from 7.3 The hat matrix.

Mean piece. The mean term vanishes, because $(\mathbf{I} - \mathbf{H})\mathbf{X} = \mathbf{X} - \mathbf{H}\mathbf{X} = \mathbf{X} - \mathbf{X} = \mathbf{0}$ (the hat matrix leaves the columns of $\mathbf{X}$ unchanged, since they already lie in the column space). Therefore

Combining, $E\{\mathrm{SSE}\} = \sigma^2(n - p) + 0 = (n - p)\sigma^2$, so

The divisor $n - p$, and no other, makes $\mathrm{MSE}$ unbiased. $\blacksquare$

This is the Chapter 2 result $E\{\mathrm{MSE}\} = \sigma^2$ generalized: there $p = 2$ gave the divisor $n - 2$; here any model divides by $n$ minus the number of parameters. Figure 13 shows what goes wrong with the naive divisor $n$: dividing by $n$ instead of $n - p$ biases the estimate low, and the gap is largest when $p$ is a big fraction of $n$.

R and Python¶

We can also read off how well the plane fits overall. The total variation $\mathrm{SSTO} = \sum (Y_i - \bar{Y})^2$ splits into the part the model explains and the part it does not, and $R^2 = 1 - \mathrm{SSE}/\mathrm{SSTO}$ reports the explained fraction, exactly as in 3.6 The analysis of variance approach but now for many predictors.

SSTO <- sum((Y - mean(Y))^2)
SSR <- SSTO - SSE
c(SSTO = round(SSTO, 2), SSR = round(SSR, 2), SSE = round(SSE, 2),
  R2 = round(1 - SSE / SSTO, 4))

      SSTO        SSR        SSE         R2
26196.2100 24015.2800  2180.9300     0.9167

SSTO = np.sum((Y - Y.mean()) ** 2)
SSR = SSTO - SSE
print(round(SSTO, 2), round(SSR, 2), round(SSE, 2), round(1 - SSE / SSTO, 4))

26196.21 24015.28 2180.93 0.9167

The two predictors together explain about $91.7\%$ of the variation in Dwaine sales. Chapter 8 unpacks this number, warns about how it can only rise when predictors are added, and builds the general linear test on top of the $\mathrm{SSE}$ we just computed (8.3 Extra sums of squares).

7.6 The Gauss-Markov theorem¶

Intuition¶

We chose least squares because squared error is a natural penalty, and Chapter 2 rewarded that choice: in simple regression, no other linear unbiased estimator beats $b_1$. The same reward holds for the full vector, and the matrix proof is barely longer than the scalar one. The claim is that any competitor built as a fixed linear combination of the responses, and required to be right on average, must have a covariance matrix at least as large as $\mathbf{b}$'s. No normality is needed; the mean-zero, constant-variance, uncorrelated error assumptions are enough. Figure 14 shows the payoff on Toluca: least squares and a competing unbiased estimator are both centered on the truth, but least squares scatters far less.

Formula¶

Derivation (b is BLUE)¶

Proof. Let $\mathbf{b}^{\ast} = \mathbf{C}\mathbf{Y}$ be any linear estimator with a fixed $p \times n$ matrix $\mathbf{C}$. Its mean is $E\{\mathbf{C}\mathbf{Y}\} = \mathbf{C}\mathbf{X}\boldsymbol{\beta}$. For $\mathbf{b}^{\ast}$ to be unbiased for $\boldsymbol{\beta}$ for every possible $\boldsymbol{\beta}$, we need

Write $\mathbf{C} = (\mathbf{X}'\mathbf{X})^{-1}\mathbf{X}' + \mathbf{D}$, defining $\mathbf{D}$ as the difference from the least-squares matrix. The unbiasedness condition $\mathbf{C}\mathbf{X} = \mathbf{I}$ becomes

Now compute the covariance. Since $\operatorname{Cov}\{\mathbf{Y}\} = \sigma^2\mathbf{I}$, $\operatorname{Cov}\{\mathbf{b}^{\ast}\} = \sigma^2\mathbf{C}\mathbf{C}'$. Expand $\mathbf{C}\mathbf{C}'$ with $\mathbf{C} = (\mathbf{X}'\mathbf{X})^{-1}\mathbf{X}' + \mathbf{D}$. The two cross terms vanish because $\mathbf{D}\mathbf{X} = \mathbf{0}$ (and its transpose $\mathbf{X}'\mathbf{D}' = \mathbf{0}$):

What survives is

Therefore

For any fixed $\mathbf{c}$, the variance of the scalar $\mathbf{c}'\mathbf{b}^{\ast}$ is

The added term $\sigma^2\|\mathbf{D}'\mathbf{c}\|^2 \ge 0$ is a squared length, so every linear unbiased competitor has variance at least that of $\mathbf{c}'\mathbf{b}$, with equality only when $\mathbf{D}'\mathbf{c} = \mathbf{0}$; taking $\mathbf{c}$ to range over the coordinate vectors forces $\mathbf{D} = \mathbf{0}$, that is, the competitor is $\mathbf{b}$ itself. Least squares is BLUE. $\blacksquare$

Setting $\mathbf{X} = [\mathbf{1}, \mathbf{x}]$ and $\mathbf{c} = (0, 1)'$ recovers the Chapter 2 statement that $b_1$ is BLUE (2.5 Sampling behavior and the Gauss-Markov theorem); the matrix proof did every coefficient, and every linear combination of coefficients, at once.

7.7 Simple regression as a special case¶

Intuition¶

Everything in Chapters 2 and 3 should drop out of the matrix formulas when there is one predictor, and it does. Seeing it happen once turns the general formulas from abstractions into trusted tools: the same machine that fit Dwaine’s plane fits Toluca’s line, and returns the exact numbers you already know. Figure 15 is the whole idea in one picture: one formula, two jobs.

Derivation (matrix formulas reduce to Chapter 2)¶

Derivation (the $p = 2$ case). With one predictor, the design matrix is $\mathbf{X} = [\mathbf{1}, \mathbf{x}]$, so

where $S_{xx} = \sum (X_i - \bar{X})^2$ and $S_{xy} = \sum (X_i - \bar{X})(Y_i - \bar{Y})$ are the Chapter 2 sums of squares and cross-products, and $\bar{X}, \bar{Y}$ are the sample means. The inverse of a $2 \times 2$ matrix divides by its determinant $\det(\mathbf{X}'\mathbf{X}) = n\sum X_i^2 - (\sum X_i)^2 = n\,S_{xx}$ (6.4 The inverse and the normal equations),

Multiplying by $\mathbf{X}'\mathbf{Y}$ and simplifying (each entry collapses using $\sum X_i = n\bar{X}$ and the definitions of $S_{xx}, S_{xy}$) gives

which are exactly the Chapter 2 estimators $b_0 = \bar{Y} - b_1\bar{X}$ and $b_1 = S_{xy}/S_{xx}$. The covariance matrix reduces the same way: the bottom-right entry of $\sigma^2(\mathbf{X}'\mathbf{X})^{-1}$ is $\sigma^2 n / (n S_{xx}) = \sigma^2/S_{xx}$, matching $\operatorname{Var}\{b_1\}$, and the top-left entry is $\sigma^2 \sum X_i^2/(n S_{xx}) = \sigma^2(1/n + \bar{X}^2/S_{xx})$, matching $\operatorname{Var}\{b_0\}$ from 2.5 Sampling behavior and the Gauss-Markov theorem. $\blacksquare$

Example 7.5 already confirmed this numerically: the matrix pipeline returned Toluca’s $b_0 = 62.3659$, $b_1 = 3.5702$, and standard errors 26.1774 and 0.3470, the same figures Chapter 2 produced from the scalar formulas. The general linear model does not replace simple regression; it contains it.

7.8 Chapter summary¶

This chapter turned Chapter 6’s matrix arithmetic into a full multiple regression. One equation, $\mathbf{Y} = \mathbf{X}\boldsymbol{\beta} + \boldsymbol{\varepsilon}$, carries any number of predictors. Least squares projects $\mathbf{Y}$ onto the column space of $\mathbf{X}$, giving $\mathbf{b} = (\mathbf{X}'\mathbf{X})^{-1}\mathbf{X}'\mathbf{Y}$ by calculus and by geometry alike. The hat matrix $\mathbf{H}$ packages that projection, and its two properties, being a projection (symmetric and idempotent) with trace $p$, drive the sampling distribution of $\mathbf{b}$, the $n - p$ divisor in $\mathrm{MSE}$, and the Gauss-Markov optimality. Every Chapter 2 formula returns as the $p = 2$ special case.

Key results at a glance.

Key terms. general linear model, design matrix, normal equations, least-squares estimator, column space, hat matrix, idempotent matrix, full column rank, error degrees of freedom, best linear unbiased estimator (BLUE).

You should now be able to.

• Write any linear regression as $\mathbf{Y} = \mathbf{X}\boldsymbol{\beta} + \boldsymbol{\varepsilon}$ and identify the design matrix, coefficient vector, and error vector.

• Derive $\mathbf{b} = (\mathbf{X}'\mathbf{X})^{-1}\mathbf{X}'\mathbf{Y}$ two ways: by calculus and by orthogonal projection.

• Interpret least squares geometrically as projecting $\mathbf{Y}$ onto the column space of $\mathbf{X}$, with the residual perpendicular to it.

• Prove that the hat matrix $\mathbf{H}$ is symmetric and idempotent with trace $p$, and use it to write fitted values and residuals.

• Derive the mean and covariance of $\mathbf{b}$ and state its sampling distribution under normal errors.

• Prove that $s^2 = \mathrm{MSE} = \mathrm{SSE}/(n-p)$ is unbiased for $\sigma^2$ using the trace of $\mathbf{I}-\mathbf{H}$.

• State and prove the Gauss-Markov theorem in matrix form: $\mathbf{b}$ is the best linear unbiased estimator.

• Recover every simple-regression formula from Chapter 2 as the special case $p = 2$.

Where this fits. This chapter is the general engine for the FIT stage of the course workflow (The modeling workflow): it estimates the coefficients and the error variance for a model of any size, replacing the predictor-by-predictor algebra of Chapters 2 and 3 with one formula that a computer can run for ten predictors as easily as one. It also arms the CHECK stage, because the hat matrix $\mathbf{H}$ is the object Chapter 9 (9. Model diagnostics) turns into the leverage values and influence measures, and the USE stage, because $\operatorname{Cov}\{\mathbf{b}\}$ is what Chapter 8 turns into tests and intervals. Chapter 8 takes the Dwaine fit built here, interprets each coefficient with care, and develops the extra-sums-of-squares and general-linear-test machinery (8.3 Extra sums of squares) that decides which predictors a model needs.

7.9 Frequently asked questions¶

Q1. Why does the design matrix need a column of ones? The column of ones lets the intercept $\beta_0$ be handled like any other coefficient: its “predictor value” is 1 for every observation, so $\beta_0 \cdot 1$ is the same baseline for every row. Without it the fitted plane would be forced through the origin. The first normal equation, coming from that column, is exactly $\sum e_i = 0$.

Q2. What does full column rank mean, and what breaks without it? Full column rank means no predictor column is a linear combination of the others (including the ones column). If one is, say a predictor recorded in both meters and feet, then $\mathbf{X}'\mathbf{X}$ is singular and $(\mathbf{X}'\mathbf{X})^{-1}$ does not exist, so there is no unique $\mathbf{b}$: infinitely many coefficient vectors fit equally well. Software then drops a column or returns an error. Near-violations (highly correlated predictors) are multicollinearity, the subject of Chapter 12.

Q3. Is the hat matrix worth computing for real data? As a full $n \times n$ matrix, usually not; software gets fitted values and residuals more cheaply. But its diagonal entries $h_{ii}$, the leverage values, are computed and used constantly in diagnostics, and its algebra (symmetric, idempotent, trace $p$) is what makes the theory work. Think of $\mathbf{H}$ as a proof tool and a source of the leverages, not a thing you print.

Q4. Why is the divisor $n - p$ and not $n - 1$ or $n$? Each parameter you estimate uses one degree of freedom before the residuals are measured. A sample variance estimates one quantity (the mean) and divides by $n - 1$; a regression with $p$ parameters divides by $n - p$. The trace derivation in 7.5 Estimating the error variance shows $E\{\mathrm{SSE}\} = (n - p)\sigma^2$ exactly, so $n - p$ is the unique divisor making $\mathrm{MSE}$ unbiased. Dividing by $n$ (the maximum-likelihood choice) is biased low.

Q5. Do I need normal errors for any of this? No, for most of it. The estimator $\mathbf{b}$, its unbiasedness, its covariance $\sigma^2(\mathbf{X}'\mathbf{X})^{-1}$, the unbiasedness of $\mathrm{MSE}$, and the Gauss-Markov optimality all hold with only mean-zero, constant-variance, uncorrelated errors. Normality is an extra assumption, used to get the exact multivariate normal sampling distribution of $\mathbf{b}$ and the $t$ and $F$ tests of Chapter 8.

Q6. The Dwaine coefficients from the matrix formula and from lm matched exactly. Is that always guaranteed? Yes, up to tiny floating-point differences. Both compute the same $(\mathbf{X}'\mathbf{X})^{-1}\mathbf{X}'\mathbf{Y}$; lm and statsmodels use a more numerically stable routine (a QR or SVD decomposition) instead of forming the inverse explicitly, which matters only for badly conditioned designs. For well-behaved data like these, every printed digit agrees.

Q7. Why are the two Dwaine slope estimates negatively correlated in Figure 10? Because the two predictors are themselves positively correlated across cities. When a sample happens to make the population slope come out high, the fit compensates by making the income slope come out low to avoid double-counting the shared signal. That trade-off is exactly the off-diagonal entry of $\sigma^2(\mathbf{X}'\mathbf{X})^{-1}$, and it is why you cannot read one coefficient’s uncertainty without the whole covariance matrix.

7.10 Practice problems¶

1. (A) Write the Dwaine model in the form $\mathbf{Y} = \mathbf{X}\boldsymbol{\beta} + \boldsymbol{\varepsilon}$, giving the dimensions of each of the four objects and describing the first column of $\mathbf{X}$.

2. (A) Explain in one sentence each what $\mathbf{X}'\mathbf{X}$ and $\mathbf{X}'\mathbf{Y}$ represent, and why the normal equations use them.

3. (A) State the two error assumptions $E\{\boldsymbol{\varepsilon}\} = \mathbf{0}$ and $\operatorname{Cov}\{\boldsymbol{\varepsilon}\} = \sigma^2\mathbf{I}$ in plain words, saying what the diagonal and the off-diagonal of $\sigma^2\mathbf{I}$ each encode.

4. (A) In the projection picture, what is the column space of $\mathbf{X}$, where does $\hat{\mathbf{Y}}$ sit, and in what direction does the residual point? Answer without formulas.

5. (A) Give the three properties of the hat matrix proved in this chapter and say what each one means geometrically.

6. (A) Why is the divisor in $\mathrm{MSE}$ equal to $n - p$? Answer by naming what the number $p$ counts and what $n - p$ counts.

7. (A) A colleague says “with normality the least-squares estimator is BLUE; without it, it is not.” Correct the statement, naming exactly which assumptions the Gauss-Markov theorem uses.

8. (A) Explain why $\mathbf{X}\mathbf{X}'$ (an $n \times n$ matrix) cannot replace $\mathbf{X}'\mathbf{X}$ in the normal equations.

9. (B) Starting from $Q(\mathbf{b}) = (\mathbf{Y} - \mathbf{X}\mathbf{b})'(\mathbf{Y} - \mathbf{X}\mathbf{b})$, expand it, take the gradient, and derive the normal equations and $\mathbf{b} = (\mathbf{X}'\mathbf{X})^{-1}\mathbf{X}'\mathbf{Y}$ (Theorem 7.3). Justify why the stationary point is a minimum.

10. (B) Prove the projection identity $\|\mathbf{Y} - \mathbf{X}\mathbf{b}^{\ast}\|^2 = \|\mathbf{e}\|^2 + \|\mathbf{X}(\mathbf{b} - \mathbf{b}^{\ast})\|^2$ for any $\mathbf{b}^{\ast}$, showing the cross term vanishes, and conclude that $\mathbf{b}$ uniquely minimizes SSE.

11. (B) Prove that $\mathbf{H} = \mathbf{X}(\mathbf{X}'\mathbf{X})^{-1}\mathbf{X}'$ is symmetric and idempotent (Theorem 7.7). Then prove $\mathbf{I} - \mathbf{H}$ is symmetric and idempotent as well.

12. (B) Prove $\operatorname{tr}(\mathbf{H}) = p$ using the cyclic property of the trace, and deduce $\operatorname{tr}(\mathbf{I} - \mathbf{H}) = n - p$.

13. (B) Show $\mathbf{H}\mathbf{X} = \mathbf{X}$ and use it to prove $(\mathbf{I} - \mathbf{H})\mathbf{X} = \mathbf{0}$. Explain in words why the hat matrix leaves the columns of $\mathbf{X}$ unchanged.

14. (B) Derive $E\{\mathbf{b}\} = \boldsymbol{\beta}$ and $\operatorname{Cov}\{\mathbf{b}\} = \sigma^2(\mathbf{X}'\mathbf{X})^{-1}$ (Theorem 7.8), naming the random-vector facts you use at each step. Then, adding the assumption $\boldsymbol{\varepsilon} \sim N(\mathbf{0}, \sigma^2\mathbf{I})$, state the full sampling distribution of $\mathbf{b}$ and explain why the mean and covariance did not need normality but this distributional statement does.

15. (B) Prove the quadratic-form lemma $E\{\mathbf{Y}'\mathbf{A}\mathbf{Y}\} = \operatorname{tr}(\mathbf{A}\boldsymbol{\Sigma}) + \boldsymbol{\mu}'\mathbf{A}\boldsymbol{\mu}$, then use it to show $E\{\mathrm{SSE}\} = (n - p)\sigma^2$ (Theorem 7.9).

16. (B) State and prove the Gauss-Markov theorem in matrix form (Theorem 7.10). Where in the proof does the assumption $\operatorname{Cov}\{\boldsymbol{\varepsilon}\} = \sigma^2\mathbf{I}$ enter, and where does full column rank enter?

17. (B) For the simple-regression design $\mathbf{X} = [\mathbf{1}, \mathbf{x}]$, compute $\mathbf{X}'\mathbf{X}$, its determinant, and $(\mathbf{X}'\mathbf{X})^{-1}$, and derive $b_1 = S_{xy}/S_{xx}$ and $b_0 = \bar{Y} - b_1\bar{X}$ from the matrix formula.

18. (B) Using the reduction in problem 17, derive $\operatorname{Var}\{b_1\} = \sigma^2/S_{xx}$ and $\operatorname{Var}\{b_0\} = \sigma^2(1/n + \bar{X}^2/S_{xx})$ from $\sigma^2(\mathbf{X}'\mathbf{X})^{-1}$.

19. (B) Show that $\operatorname{Cov}\{\mathbf{e}, \hat{\mathbf{Y}}\} = \mathbf{0}$ under the model, using $\hat{\mathbf{Y}} = \mathbf{H}\mathbf{Y}$, $\mathbf{e} = (\mathbf{I} - \mathbf{H})\mathbf{Y}$, and $\mathbf{H}(\mathbf{I} - \mathbf{H}) = \mathbf{0}$. Why does this matter for diagnostics?

20. (B) Prove that adding a predictor (a new column to $\mathbf{X}$) cannot increase $\mathrm{SSE}$, using the projection argument that enlarging the column space can only move $\hat{\mathbf{Y}}$ closer to $\mathbf{Y}$.

21. (C) Read dwaine.csv, build $\mathbf{X}$ and $\mathbf{Y}$, and compute $\mathbf{b}$ from $(\mathbf{X}'\mathbf{X})^{-1}\mathbf{X}'\mathbf{Y}$ in R or Python. Confirm your answer matches lm / smf.ols.

22. (C) For the Dwaine fit, compute the residual vector and verify numerically that $\mathbf{X}'\mathbf{e} = \mathbf{0}$ (all three entries) and that $\sum \hat{Y}_i = \sum Y_i$.

23. (C) Build the Dwaine hat matrix and verify it is symmetric and idempotent and has trace 3. Then report the five largest diagonal entries $h_{ii}$ and say what a large $h_{ii}$ indicates about that city’s predictors.

24. (C) Compute $\mathrm{SSE}$, $\mathrm{MSE}$, and $s$ for Dwaine two ways: from the residual vector, and from summary / .summary2(). Then compute $\mathrm{SSE}/n$ and report the percentage by which it understates $\mathrm{MSE}$.

25. (C) Compute the estimated covariance matrix $\mathrm{MSE}\,(\mathbf{X}'\mathbf{X})^{-1}$ for Dwaine, report the three standard errors, and confirm they match the built-in summary. Report the correlation between the two slope estimates from the off-diagonal entry.

26. (C) Refit Toluca through the matrix formulas and confirm $b_0 = 62.3659$, $b_1 = 3.5702$, $\mathrm{MSE} = 2383.72$, and the two standard errors, matching Chapter 2 to every digit.

27. (C) Write a simulation (seed 4210) that draws 2000 datasets from the fitted Dwaine model, refits each, and reports the mean of the 2000 estimates of each coefficient (checking unbiasedness) and the standard deviation of the target-population slope (checking $\operatorname{Cov}\{\mathbf{b}\}$).

28. (C) Regress Dwaine sales on targtpop alone, then on targtpop and dispoinc together, and confirm numerically that adding dispoinc lowers $\mathrm{SSE}$. Relate the drop to the projection argument in problem 20.

7.11 Exam practice¶

EP 7.1 (interpret output in context). A colleague fits the Dwaine model in R and reads off the coefficient table.

            Estimate Std. Error t value Pr(>|t|)
(Intercept) -68.8571    60.0170 -1.1473   0.2663
targtpop      1.4546     0.2118  6.8682   0.0000
dispoinc      9.3655     4.0640  2.3045   0.0333

Interpret the targtpop coefficient and its $t$ value in the context of studio sales, being explicit about what is held fixed. Then evaluate the colleague’s conclusion: “the intercept has $p = 0.27$, so it is not significantly different from zero, which means the model is misspecified and we should drop the intercept.” Say whether the conclusion is right and why.

EP 7.2 (explain why, and what would change if). Explain in full sentences why the divisor in $\mathrm{MSE} = \mathrm{SSE}/(n - p)$ is $n - p$ and not $n$. Then answer a “what would change if”: suppose an analyst kept adding predictors to the Dwaine data until the model had $p = 21$ parameters, one for every observation. Describe what happens to $\mathrm{SSE}$, to the divisor, and to $\mathrm{MSE}$, and explain why a fit with $\mathrm{SSE} = 0$ is a warning rather than a success.

EP 7.3 (a student claims X, evaluate). A student writes: “Adding a predictor can never increase $\mathrm{SSE}$, and it always raises $R^2$. So the best strategy is to keep adding predictors until $R^2$ is as close to 1 as possible.” Using the projection picture of least squares, evaluate the two factual claims and the strategy they support.

EP 7.4 (interpret computed output in context). The leverage values $h_{ii}$ are the diagonal entries of the hat matrix $\mathbf{H} = \mathbf{X}(\mathbf{X}'\mathbf{X})^{-1}\mathbf{X}'$ for the Dwaine fit. The six largest, and a check on their sum, are below.

row 20  targtpop=  82.7  dispoinc= 19.1  h_ii=0.2788
row 13  targtpop=  88.4  dispoinc= 17.4  h_ii=0.2390
row 15  targtpop=  52.5  dispoinc= 17.8  h_ii=0.2095
row  3  targtpop=  91.3  dispoinc= 18.2  h_ii=0.1737
row  5  targtpop=  46.9  dispoinc= 17.3  h_ii=0.1620
row 19  targtpop=  89.6  dispoinc= 18.1  h_ii=0.1611
sum(h_ii) = 3.0000     mean(targtpop)=62.02  mean(dispoinc)=17.14

Explain what a large $h_{ii}$ says about a city, and interpret why row 20 has the largest leverage value even though its target population (82.7) is not the largest in the data. Then explain why the leverage values sum to exactly 3, and use the rule-of-thumb cutoff $2p/n = 0.286$ to say whether any Dwaine city should be called a high-leverage point.

EP 7.5 (what would change if). An analyst re-expresses target population in hundreds of people instead of thousands, multiplying every targtpop value by 10, and refits. The two coefficient tables are below.

original                              rescaled
             Estimate Std. Error       Estimate Std. Error
(Intercept)  -68.8571    60.0170       -68.8571    60.0170
targtpop       1.4546     0.2118        0.1455     0.0212
dispoinc       9.3655     4.0640        9.3655     4.0640
original  MSE=121.16  R2=0.9167     rescaled  MSE=121.16  R2=0.9167

Explain, in full sentences and with reference to the model algebra, why the target- population coefficient and its standard error both divide by 10 while its $t$ value, the fitted values, $\mathrm{MSE}$, and $R^2$ are all unchanged. What does this tell you about reading the size of a coefficient as a measure of a predictor’s importance?