Multiple Linear Regression in Matrix Form: Orthogonal Projections, Hat Matrix & Partial Correlations

In observational econometrics, agronomy, clinical biostatistics, and industrial experimentation, a response variable $Y$ is almost always governed by several simultaneous explanatory variables $X_1, X_2, \dots, X_k$. Omitting relevant covariates that are correlated with the included predictors induces omitted-variable bias, whereas incorporating them explicitly allows the analyst to isolate the partial effect of each regressor while holding all other covariates constant. Suppose we observe $n$ experimental units, where unit $i$ yields a continuous response $Y_i$ along with $k$ predictor measurements $(x_{i1}, x_{i2}, \dots, x_{ik})$. The Multiple Linear Regression Model with an intercept $\beta_0$ and $k$ slope parameters $(\beta_1, \dots, \beta_k)$ is written for each unit $i = 1, 2, \dots, n$ as:

Setting $p = k + 1$ to denote the total number of regression parameters (including the intercept $\beta_0$, which is associated with a constant dummy predictor $x_{i0} \equiv 1$), we assemble all $n$ equations compactly into a single matrix equation:

where the constituent vectors and matrices have dimensions:

Here, $\mathbf{Y}$ is the $n \times 1$ observable response vector, $\mathbf{X} = [\mathbf{1}_n, \mathbf{x}_{(1)}, \dots, \mathbf{x}_{(k)}]$ is the $n \times p$ model (or design) matrix whose rows are $\mathbf{x}_i^T = (1, x_{i1}, \dots, x_{ik})$, $\boldsymbol{\beta}$ is the $p \times 1$ vector of unknown population regression parameters, and $\boldsymbol{\varepsilon}$ is the $n \times 1$ unobservable random disturbance vector.

The classical matrix Gauss–Markov assumptions state that:

1. Expectation Specification: $\mathbb{E}[\boldsymbol{\varepsilon} \mid \mathbf{X}] = \mathbf{0}_n$, which implies $\mathbb{E}[\mathbf{Y} \mid \mathbf{X}] = \mathbf{X}\boldsymbol{\beta}$.

2. Spherical Error Covariance: $\text{Cov}(\boldsymbol{\varepsilon} \mid \mathbf{X}) = \mathbb{E}[\boldsymbol{\varepsilon}\boldsymbol{\varepsilon}^T \mid \mathbf{X}] = \sigma^2 \mathbf{I}_n$, where $\mathbf{I}_n$ is the $n \times n$ identity matrix and $\sigma^2 > 0$ is the common error variance. Consequently, $\text{Cov}(\mathbf{Y} \mid \mathbf{X}) = \sigma^2 \mathbf{I}_n$.

3. Full Column Rank Condition: The design matrix $\mathbf{X}$ has full column rank:

The full column rank condition $\text{rank}(\mathbf{X}) = p$ has a clear algebraic and empirical meaning: the $p$ column vectors of $\mathbf{X}$ are linearly independent in $\mathbb{R}^n$, meaning no predictor column can be written as an exact linear combination of the remaining columns (no exact multicollinearity), and the sample size $n$ is at least as large as the number of parameters $p$. Under $\text{rank}(\mathbf{X}) = p$, the $p \times p$ symmetric Gram (cross-product) matrix $\mathbf{X}^T\mathbf{X}$ is strictly positive definite, because for any non-zero vector $\mathbf{v} \in \mathbb{R}^p \setminus \{\mathbf{0}\}$, we have $\mathbf{X}\mathbf{v} \neq \mathbf{0}$ and therefore:

Since every positive definite matrix is nonsingular, the inverse matrix $(\mathbf{X}^T\mathbf{X})^{-1}$ exists uniquely and is itself symmetric and positive definite.

To derive the Ordinary Least Squares estimator $\hat{\boldsymbol{\beta}}$, we minimize the sum of squared residuals $Q(\mathbf{b})$ over all candidate parameter vectors $\mathbf{b} \in \mathbb{R}^p$:

where we used the scalar transpose identity $\mathbf{Y}^T \mathbf{X} \mathbf{b} = (\mathbf{Y}^T \mathbf{X} \mathbf{b})^T = \mathbf{b}^T \mathbf{X}^T \mathbf{Y}$. Using standard rules of vector calculus ($\nabla_{\mathbf{b}}(\mathbf{b}^T \mathbf{a}) = \mathbf{a}$ and $\nabla_{\mathbf{b}}(\mathbf{b}^T \mathbf{A} \mathbf{b}) = 2 \mathbf{A}\mathbf{b}$ for symmetric $\mathbf{A}$), we differentiate $Q(\mathbf{b})$ with respect to $\mathbf{b}$ and set the gradient vector equal to $\mathbf{0}_p$:

Replacing $\mathbf{b}$ by $\hat{\boldsymbol{\beta}}$ yields the Matrix Normal Equations:

Premultiplying both sides by $(\mathbf{X}^T\mathbf{X})^{-1}$ produces the unique closed-form Ordinary Least Squares Estimator:

Furthermore, the Hessian matrix of second derivatives is $\nabla_{\mathbf{b}}^2 Q(\mathbf{b}) = 2 \mathbf{X}^T\mathbf{X}$, which is strictly positive definite, confirming that $\hat{\boldsymbol{\beta}}$ is the unique global minimizer of $Q(\mathbf{b})$. Alternatively, without using calculus, one can complete the square algebraically for any $\mathbf{b} \in \mathbb{R}^p$:

because the cross-term $2(\hat{\boldsymbol{\beta}} - \mathbf{b})^T \mathbf{X}^T (\mathbf{Y} - \mathbf{X}\hat{\boldsymbol{\beta}}) = 2(\hat{\boldsymbol{\beta}} - \mathbf{b})^T (\mathbf{X}^T\mathbf{Y} - \mathbf{X}^T\mathbf{X}\hat{\boldsymbol{\beta}}) = 0$ vanishes by the normal equations. Since $\mathbf{X}^T\mathbf{X}$ is positive definite, the second term is strictly positive whenever $\mathbf{b} \neq \hat{\boldsymbol{\beta}}$ and equals zero uniquely at $\mathbf{b} = \hat{\boldsymbol{\beta}}$.

Numerical Linear Algebra of Normal Equations: Cholesky vs QR Factorization

Although the closed-form matrix expression writes the Ordinary Least Squares estimator in terms of the explicit inverse of the cross-product Gram matrix, computational statisticians never invert the Gram matrix directly in finite-precision floating-point arithmetic. Forming the cross-product matrix squares the spectral condition number of the design matrix, which doubles the loss of significant digits when predictors are moderately correlated. Instead, statistical computing environments such as R and Python's SciPy library solve the least-squares problem either by computing the Cholesky factorization of the positive-definite Gram matrix or, far more stably, by applying an orthogonal-triangular QR factorization directly to the design matrix itself. Expressing the design matrix as the product of an column-orthonormal matrix and an upper-triangular matrix reduces the normal equations to an upper-triangular linear system that can be solved rapidly and accurately by back-substitution without ever forming the cross-product matrix.

While matrix calculus provides a rapid derivation of $\hat{\boldsymbol{\beta}}$, the deepest insights into multiple linear regression emerge from its **Euclidean geometry in $\mathbb{R}^n$**. Rather than viewing the data as $n$ points scattered in a $(k+1)$-dimensional variable space, the geometric perspective views the response vector $\mathbf{Y} \in \mathbb{R}^n$ and the $p$ columns of the design matrix $\mathbf{X} = [\mathbf{x}_{(0)}, \mathbf{x}_{(1)}, \dots, \mathbf{x}_{(k)}]$ as vectors residing in the $n$-dimensional observation space $\mathbb{R}^n$. As the parameter vector $\mathbf{b} \in \mathbb{R}^p$ varies over all possible values, the candidate mean vector $\boldsymbol{\mu}(\mathbf{b}) = \mathbf{X}\mathbf{b} = \sum_{j=0}^k b_j \mathbf{x}_{(j)}$ sweeps out a $p$-dimensional linear subspace of $\mathbb{R}^n$ known as the **column space** (or estimation space) of $\mathbf{X}$: $$\mathcal{C}(\mathbf{X}) = \text{Col}(\mathbf{X}) = \left\{ \mathbf{X}\mathbf{b} : \mathbf{b} \in \mathbb{R}^p \right\} \subset \mathbb{R}^n$$ Because of random disturbances $\boldsymbol{\varepsilon}$, the observed response vector $\mathbf{Y}$ almost never lies inside the $p$-dimensional subspace $\mathcal{C}(\mathbf{X})$. The Ordinary Least Squares optimization problem—minimizing $\|\mathbf{Y} - \mathbf{X}\mathbf{b}\|_2^2$—asks for the vector $\hat{\mathbf{Y}} = \mathbf{X}\hat{\boldsymbol{\beta}} \in \mathcal{C}(\mathbf{X})$ whose Euclidean distance to $\mathbf{Y}$ is minimal. By the classical **Orthogonal Projection Theorem** of Hilbert spaces, the unique point $\hat{\mathbf{Y}} \in \mathcal{C}(\mathbf{X})$ closest to $\mathbf{Y}$ is the orthogonal projection of $\mathbf{Y}$ onto $\mathcal{C}(\mathbf{X})$, characterized by the condition that the residual vector $\mathbf{e} = \mathbf{Y} - \hat{\mathbf{Y}}$ must be orthogonal to every vector in $\mathcal{C}(\mathbf{X})$ (and hence orthogonal to every column of $\mathbf{X}$): $$\mathbf{X}^T \mathbf{e} = \mathbf{0}_p \iff \mathbf{X}^T (\mathbf{Y} - \mathbf{X}\hat{\boldsymbol{\beta}}) = \mathbf{0}_p \iff (\mathbf{X}^T\mathbf{X})\hat{\boldsymbol{\beta}} = \mathbf{X}^T\mathbf{Y}$$ Thus, the normal equations are simply the algebraic statement that the residual vector $\mathbf{e}$ lies in the $(n - p)$-dimensional **orthogonal complement** (or left null space) $\mathcal{C}(\mathbf{X})^\perp = \mathcal{N}(\mathbf{X}^T)$! Substituting $\hat{\boldsymbol{\beta}} = (\mathbf{X}^T\mathbf{X})^{-1}\mathbf{X}^T\mathbf{Y}$ into the fitted response vector $\hat{\mathbf{Y}} = \mathbf{X}\hat{\boldsymbol{\beta}}$ yields: $$\hat{\mathbf{Y}} = \mathbf{X} (\mathbf{X}^T \mathbf{X})^{-1} \mathbf{X}^T \mathbf{Y} = \mathbf{H} \mathbf{Y}$$ where the $n \times n$ matrix $\mathbf{H}$ is defined as the **Hat Matrix** (so named by John Tukey because it maps $\mathbf{Y}$ into "$\mathbf{Y}$-hat"): $$\mathbf{H} = \mathbf{X} (\mathbf{X}^T \mathbf{X})^{-1} \mathbf{X}^T \in \mathbb{R}^{n \times n}$$ Correspondingly, the OLS residual vector $\mathbf{e} = \mathbf{Y} - \hat{\mathbf{Y}}$ is obtained via the **residual-maker (annihilator) matrix** $\mathbf{M} = \mathbf{I}_n - \mathbf{H}$: $$\mathbf{e} = \mathbf{Y} - \mathbf{H}\mathbf{Y} = (\mathbf{I}_n - \mathbf{H})\mathbf{Y} = \mathbf{M}\mathbf{Y}$$ The matrices $\mathbf{H}$ and $\mathbf{M} = \mathbf{I}_n - \mathbf{H}$ are the orthogonal projection operators onto $\mathcal{C}(\mathbf{X})$ and $\mathcal{C}(\mathbf{X})^\perp$, respectively, and satisfy six fundamental algebraic properties: Furthermore, because $\mathbf{M}\mathbf{X} = \mathbf{0}_{n \times p}$, substituting $\mathbf{Y} = \mathbf{X}\boldsymbol{\beta} + \boldsymbol{\varepsilon}$ into $\mathbf{e} = \mathbf{M}\mathbf{Y}$ proves that the observable residual vector is an exact linear transformation of the unobservable error vector $\boldsymbol{\varepsilon}$ alone, completely free of $\boldsymbol{\beta}$: $$\mathbf{e} = \mathbf{M}(\mathbf{X}\boldsymbol{\beta} + \boldsymbol{\varepsilon}) = \mathbf{M}\mathbf{X}\boldsymbol{\beta} + \mathbf{M}\boldsymbol{\varepsilon} = \mathbf{M}\boldsymbol{\varepsilon} = (\mathbf{I}_n - \mathbf{H})\boldsymbol{\varepsilon}$$ Consequently, $\mathbb{E}[\mathbf{e}] = \mathbf{0}_n$ and $\text{Cov}(\mathbf{e}) = \mathbf{M} \text{Cov}(\boldsymbol{\varepsilon}) \mathbf{M}^T = \sigma^2 \mathbf{M} \mathbf{I}_n \mathbf{M} = \sigma^2 \mathbf{M} = \sigma^2 (\mathbf{I}_n - \mathbf{H})$. Notably, even though the true errors $\varepsilon_i$ are homoscedastic and uncorrelated ($\text{Cov}(\boldsymbol{\varepsilon}) = \sigma^2 \mathbf{I}_n$), the OLS residuals $e_i$ are **heteroscedastic and correlated**, with variance $\text{Var}(e_i) = \sigma^2(1 - h_{ii})$ and covariance $\text{Cov}(e_i, e_j) = -\sigma^2 h_{ij}$. ### Invariance of the Orthogonal Projection Operator Under Nonsingular Reparameterizations A profound geometric property of the Hat matrix and the residual-maker matrix is their complete invariance under any nonsingular linear reparameterization of the design matrix. Suppose an analyst rescales predictors from metric to imperial units, centers covariates around their sample means, or replaces raw variables by orthogonal linear combinations, which corresponds algebraically to post-multiplying the design matrix by an invertible square transformation matrix. While the individual parameter estimates and their covariance matrix change according to the transformation, the column space spanned by the transformed design matrix remains identical to the original column space. Consequently, the orthogonal projection matrix onto that subspace—the Hat matrix—as well as the fitted value vector, the residual vector, the leverages, and the residual sum of squares remain strictly unchanged. This geometric invariance proves that linear rescaling of predictors never alters the predictive fit or the degrees of freedom of a multiple linear regression model.

Using the matrix representation $\hat{\boldsymbol{\beta}} = (\mathbf{X}^T\mathbf{X})^{-1}\mathbf{X}^T\mathbf{Y}$, we now establish the exact finite-sample expectation vector, variance-covariance matrix, and Gauss–Markov optimality of $\hat{\boldsymbol{\beta}}$ in the general multiple regression model. Substituting $\mathbf{Y} = \mathbf{X}\boldsymbol{\beta} + \boldsymbol{\varepsilon}$ into $\hat{\boldsymbol{\beta}}$ gives the fundamental decomposition:

Taking expectations conditionally on $\mathbf{X}$ under $\mathbb{E}[\boldsymbol{\varepsilon}] = \mathbf{0}_n$ establishes exact unbiasedness:

Next, let $\mathbf{A} = (\mathbf{X}^T\mathbf{X})^{-1}\mathbf{X}^T \in \mathbb{R}^{p \times n}$ so that $\hat{\boldsymbol{\beta}} - \boldsymbol{\beta} = \mathbf{A}\boldsymbol{\varepsilon}$. Under the spherical error covariance assumption $\text{Cov}(\boldsymbol{\varepsilon}) = \mathbb{E}[\boldsymbol{\varepsilon}\boldsymbol{\varepsilon}^T] = \sigma^2 \mathbf{I}_n$, the $p \times p$ variance-covariance matrix of $\hat{\boldsymbol{\beta}}$ is:

The $j$-th diagonal element of $\sigma^2(\mathbf{X}^T\mathbf{X})^{-1}$ is the sampling variance $\text{Var}(\hat{\beta}_j) = \sigma^2 [(\mathbf{X}^T\mathbf{X})^{-1}]_{jj}$, and the $(j, m)$-th off-diagonal element is the covariance $\text{Cov}(\hat{\beta}_j, \hat{\beta}_m) = \sigma^2 [(\mathbf{X}^T\mathbf{X})^{-1}]_{jm}$.

To compare two $p \times p$ covariance matrices in multiple dimensions, mathematical statistics employs the Loewner Partial Ordering of symmetric matrices: for two $p \times p$ symmetric matrices $\mathbf{V}_1$ and $\mathbf{V}_2$, we write $\mathbf{V}_1 \succeq \mathbf{V}_2$ if and only if the difference matrix $\boldsymbol{\Delta} = \mathbf{V}_1 - \mathbf{V}_2$ is positive semi-definite (PSD), meaning $\mathbf{c}^T (\mathbf{V}_1 - \mathbf{V}_2) \mathbf{c} \ge 0$ for all vectors $\mathbf{c} \in \mathbb{R}^p$.

Theorem 2.1 (The Matrix Gauss–Markov Theorem)

Under the linear model $\mathbf{Y} = \mathbf{X}\boldsymbol{\beta} + \boldsymbol{\varepsilon}$ with $\mathbb{E}[\boldsymbol{\varepsilon}] = \mathbf{0}_n$, $\text{Cov}(\boldsymbol{\varepsilon}) = \sigma^2 \mathbf{I}_n$, and $\text{rank}(\mathbf{X}) = p$, let $\tilde{\boldsymbol{\beta}} = \mathbf{C}\mathbf{Y}$ be any competing linear unbiased estimator of $\boldsymbol{\beta}$, where $\mathbf{C} \in \mathbb{R}^{p \times n}$ is a deterministic matrix. Then the covariance matrix of $\tilde{\boldsymbol{\beta}}$ exceeds the covariance matrix of the OLS estimator $\hat{\boldsymbol{\beta}}$ by a positive semi-definite matrix:

Equivalently, for every vector $\mathbf{c} \in \mathbb{R}^p$, the scalar linear combination $\mathbf{c}^T \hat{\boldsymbol{\beta}}$ is the unique Best Linear Unbiased Estimator (BLUE) of $\mathbf{c}^T \boldsymbol{\beta}$:

with equality for all $\mathbf{c} \in \mathbb{R}^p$ if and only if $\mathbf{C} = (\mathbf{X}^T\mathbf{X})^{-1}\mathbf{X}^T$, i.e., $\tilde{\boldsymbol{\beta}} = \hat{\boldsymbol{\beta}}$.

Complete Matrix Proof:

Let $\tilde{\boldsymbol{\beta}} = \mathbf{C}\mathbf{Y}$. For $\tilde{\boldsymbol{\beta}}$ to be unbiased for $\boldsymbol{\beta}$ across all $\boldsymbol{\beta} \in \mathbb{R}^p$, we require:

Define the discrepancy matrix $\mathbf{D} = \mathbf{C} - (\mathbf{X}^T\mathbf{X})^{-1}\mathbf{X}^T \in \mathbb{R}^{p \times n}$, so that $\mathbf{C} = (\mathbf{X}^T\mathbf{X})^{-1}\mathbf{X}^T + \mathbf{D}$. Multiplying $\mathbf{D}$ on the right by $\mathbf{X}$ and applying the unbiasedness condition $\mathbf{C}\mathbf{X} = \mathbf{I}_p$ reveals the matrix orthogonality relation:

Now compute the covariance matrix of $\tilde{\boldsymbol{\beta}} = \mathbf{C}\mathbf{Y}$:

Expanding the four matrix products and observing that the two cross-terms $(\mathbf{X}^T\mathbf{X})^{-1}\mathbf{X}^T\mathbf{D}^T = (\mathbf{D}\mathbf{X}(\mathbf{X}^T\mathbf{X})^{-1})^T = \mathbf{0}_{p \times p}$ and $\mathbf{D}\mathbf{X}(\mathbf{X}^T\mathbf{X})^{-1} = \mathbf{0}_{p \times p}$ vanish identically because $\mathbf{D}\mathbf{X} = \mathbf{0}_{p \times p}$, we obtain:

For any arbitrary vector $\mathbf{c} \in \mathbb{R}^p$, the variance of $\mathbf{c}^T \tilde{\boldsymbol{\beta}}$ is:

Because $\|\mathbf{D}^T\mathbf{c}\|_2^2 \ge 0$, the matrix $\sigma^2 \mathbf{D}\mathbf{D}^T$ is positive semi-definite, and vanishes for all $\mathbf{c}$ if and only if $\mathbf{D} = \mathbf{0}_{p \times n}$. Consequently, $\hat{\boldsymbol{\beta}}$ simultaneously minimizes the individual variance $\text{Var}(\hat{\beta}_j)$ of every coordinate, the trace $\text{tr}(\text{Cov}(\tilde{\boldsymbol{\beta}}))$, and the generalized variance $\det(\text{Cov}(\tilde{\boldsymbol{\beta}}))$.

Geometric and Decision-Theoretic Meaning of the Loewner Positive Semi-Definite Ordering

The Matrix Gauss–Markov Theorem elevates optimality from a scalar variance inequality to a matrix ordering in the space of symmetric positive semi-definite matrices. Because the difference between the covariance matrix of any competing linear unbiased estimator and the covariance matrix of the Ordinary Least Squares estimator is positive semi-definite, three powerful consequences follow simultaneously. First, every diagonal element of the difference matrix is non-negative, which guarantees that each individual partial regression coefficient is estimated with minimum variance. Second, for any linear combination of the parameters—such as a predicted mean response at an arbitrary covariate vector or a contrast between two treatment slopes—the plug-in least-squares estimator achieves the smallest possible sampling variance. Third, under multivariate normal disturbances, the concentration ellipsoid of the Ordinary Least Squares estimator is strictly contained inside the concentration ellipsoid of any other linear unbiased estimator for every probability level.

Moreover, the positive semi-definite matrix ordering established by the Matrix Gauss–Markov Theorem directly implies that the generalized variance (the determinant of the variance-covariance matrix) and the total variance (the trace of the variance-covariance matrix) are simultaneously minimized by the Ordinary Least Squares estimator across the entire class of linear unbiased estimators. Consequently, whether an investigator evaluates estimation precision via D-optimality (minimizing the determinant of the covariance matrix) or A-optimality (minimizing the trace of the covariance matrix), no linear unbiased competitor can outperform Ordinary Least Squares when the disturbances are homoscedastic and mutually uncorrelated.

To derive the exact expectation and sampling distribution of the Residual Sum of Squares $\text{SSE} = \mathbf{e}^T\mathbf{e}$ in multiple linear regression, we express sums of squares as quadratic forms $\mathbf{Y}^T \mathbf{A} \mathbf{Y}$ where $\mathbf{A} \in \mathbb{R}^{n \times n}$ is a symmetric matrix.

First, let us establish the distribution-free Trace Formula for the Expectation of a Quadratic Form: if $\mathbf{Y}$ is any random vector with mean $\mathbb{E}[\mathbf{Y}] = \boldsymbol{\mu}$ and covariance matrix $\text{Cov}(\mathbf{Y}) = \boldsymbol{\Sigma}$, and $\mathbf{A}$ is any symmetric $n \times n$ matrix, then because $\mathbf{Y}^T \mathbf{A} \mathbf{Y}$ is a $1 \times 1$ scalar equal to its own trace and the trace operator is cyclic ($\text{tr}(\mathbf{U}\mathbf{V}) = \text{tr}(\mathbf{V}\mathbf{U})$):

Let us apply this powerful identity directly to the Residual Sum of Squares $\text{SSE} = \mathbf{e}^T\mathbf{e}$. Since $\mathbf{e} = (\mathbf{I}_n - \mathbf{H})\mathbf{Y} = \mathbf{M}\mathbf{Y}$ and $\mathbf{M}$ is symmetric and idempotent ($\mathbf{M}^T\mathbf{M} = \mathbf{M}^2 = \mathbf{M}$):

where the last equality holds because $(\mathbf{I}_n - \mathbf{H})\mathbf{X} = \mathbf{0}_{n \times p}$. Applying the trace formula to $\text{SSE} = \boldsymbol{\varepsilon}^T (\mathbf{I}_n - \mathbf{H}) \boldsymbol{\varepsilon}$ with $\mathbb{E}[\boldsymbol{\varepsilon}] = \mathbf{0}_n$ and $\boldsymbol{\Sigma} = \sigma^2 \mathbf{I}_n$:

Using the cyclic trace property on $\mathbf{H} = \mathbf{X}(\mathbf{X}^T\mathbf{X})^{-1}\mathbf{X}^T$, we have:

Therefore, $\mathbb{E}[\text{SSE}] = (n - p)\sigma^2$, which proves that dividing $\text{SSE}$ by $n - p = n - k - 1$ degrees of freedom yields an unbiased estimator of $\sigma^2$:

Now suppose we invoke the normality assumption $\boldsymbol{\varepsilon} \sim \mathcal{N}_n(\mathbf{0}_n, \sigma^2 \mathbf{I}_n)$, so $\mathbf{Y} \sim \mathcal{N}_n(\mathbf{X}\boldsymbol{\beta}, \sigma^2 \mathbf{I}_n)$. Three foundational distribution theorems for normal vectors and symmetric idempotent matrices govern exact inference:

1. Idempotent Quadratic Form Chi-Square Theorem: Because $\mathbf{M} = \mathbf{I}_n - \mathbf{H}$ is symmetric and idempotent of rank $n - p$, its spectral decomposition is $\mathbf{M} = \mathbf{Q} \boldsymbol{\Lambda} \mathbf{Q}^T$, where $\mathbf{Q} \in \mathbb{R}^{n \times n}$ is orthogonal ($\mathbf{Q}^T\mathbf{Q} = \mathbf{I}_n$) and $\boldsymbol{\Lambda} = \text{diag}(1, \dots, 1, 0, \dots, 0)$ has $n - p$ ones and $p$ zeros on the diagonal. Define the rotated standard normal vector $\mathbf{Z} = \frac{1}{\sigma} \mathbf{Q}^T \boldsymbol{\varepsilon} \sim \mathcal{N}_n(\mathbf{0}_n, \mathbf{I}_n)$. Then:

Consequently, $\text{Var}(s^2) = \frac{\sigma^4}{(n-p)^2} \text{Var}(\chi_{n-p}^2) = \frac{2\sigma^4}{n - p}$.

2. Statistical Independence of $\hat{\boldsymbol{\beta}}$ and $\text{SSE}$ (Craig–Sakamoto Lemma): Both $\hat{\boldsymbol{\beta}} = (\mathbf{X}^T\mathbf{X})^{-1}\mathbf{X}^T\mathbf{Y}$ and the residual vector $\mathbf{e} = (\mathbf{I}_n - \mathbf{H})\mathbf{Y}$ are linear transformations of the multivariate normal vector $\mathbf{Y}$. Their $p \times n$ cross-covariance matrix is:

Because jointly normal vectors with zero cross-covariance are statistically independent, $\hat{\boldsymbol{\beta}}$ is completely independent of $\mathbf{e}$, and hence independent of $s^2 = \mathbf{e}^T\mathbf{e} / (n - p)$!

3. Cochran's Theorem on Partitioning Quadratic Forms: Let $\mathbf{I}_n - \frac{1}{n}\mathbf{J}_n$ be the centering matrix (where $\mathbf{J}_n = \mathbf{1}_n \mathbf{1}_n^T$ is the $n \times n$ matrix of ones). Then the ANOVA identity $\text{SST} = \text{SSR} + \text{SSE}$ corresponds to the algebraic matrix decomposition:

Because each of the three matrices on the right is symmetric and idempotent, their product $\left(\mathbf{H} - \frac{1}{n}\mathbf{J}_n\right)(\mathbf{I}_n - \mathbf{H}) = \mathbf{0}_{n \times n}$ vanishes, and their ranks sum to $(p - 1) + (n - p) = n - 1$, Cochran's Theorem guarantees that $\text{SSR} = \mathbf{Y}^T(\mathbf{H} - \frac{1}{n}\mathbf{J}_n)\mathbf{Y}$ and $\text{SSE} = \mathbf{Y}^T(\mathbf{I}_n - \mathbf{H})\mathbf{Y}$ are statistically independent, and under the global null hypothesis $H_0: \beta_1 = \dots = \beta_k = 0$, $\text{SSR}/\sigma^2 \sim \chi_k^2$.

Role of Cochran's Theorem in Partitioning Degrees of Freedom Across Subspaces

Cochran's Theorem serves as the rigorous mathematical bridge between linear algebra and chi-square distribution theory in linear models. Whenever the identity matrix (or a symmetric idempotent projection matrix) is decomposed into a sum of symmetric matrices whose ranks add up to the rank of the sum, those constituent matrices are automatically idempotent and pairwise mutually orthogonal. In the context of multiple linear regression, the total corrected sum of squares projects the response vector onto the $(n-1)$-dimensional subspace orthogonal to the constant vector of ones. This subspace splits orthogonally into the $k$-dimensional regression subspace spanned by the centered predictors and the $(n-p)$-dimensional residual subspace orthogonal to all columns of the design matrix. Because the corresponding projection matrices multiply to the zero matrix, the regression sum of squares and the residual sum of squares are statistically independent quadratic forms under Gaussian errors, justifying the exact central $F$-distribution of the analysis-of-variance test statistic.

Furthermore, the trace identity for the expectation of a quadratic form demonstrates that the unbiasedness of the residual mean square does not depend on the normality assumption. Even when the disturbances follow a skewed or heavy-tailed distribution, so long as they have zero conditional mean and spherical covariance matrix, the expected value of the residual sum of squares equals the trace of the residual-maker matrix multiplied by the error variance, which is strictly $(n-p)\sigma^2$. Gaussianity is invoked only when we move from first-moment unbiasedness to the exact chi-square distribution and exact statistical independence between the least-squares coefficient vector and the residual sum of squares.

What is the exact algebraic meaning of a single partial regression coefficient $\hat{\beta}_j$ inside a multiple linear regression model when the predictors are mutually correlated? The definitive answer is provided by the Frisch–Waugh–Lovell (FWL) Theorem (originally established by Ragnar Frisch and Frederick Waugh in 1933 for time-series detrending and generalized by Michael Lovell in 1963 to arbitrary partitioned regressions).

Partition the $n \times p$ full-rank design matrix $\mathbf{X}$ and the $p \times 1$ parameter vector $\boldsymbol{\beta}$ into two blocks of columns, $\mathbf{X} = [\mathbf{X}_1, \mathbf{X}_2]$ and $\boldsymbol{\beta} = (\boldsymbol{\beta}_1^T, \boldsymbol{\beta}_2^T)^T$, where $\mathbf{X}_1 \in \mathbb{R}^{n \times p_1}$ represents $p_1$ baseline covariates (such as an intercept and control variables) and $\mathbf{X}_2 \in \mathbb{R}^{n \times p_2}$ represents $p_2$ predictors of primary interest ($p_1 + p_2 = p$). The partitioned linear model and its partitioned OLS normal equations $(\mathbf{X}^T\mathbf{X})\hat{\boldsymbol{\beta}} = \mathbf{X}^T\mathbf{Y}$ are:

Let $\mathbf{H}_1 = \mathbf{X}_1(\mathbf{X}_1^T\mathbf{X}_1)^{-1}\mathbf{X}_1^T$ denote the orthogonal projection matrix onto $\mathcal{C}(\mathbf{X}_1)$, and let $\mathbf{M}_1 = \mathbf{I}_n - \mathbf{H}_1$ denote the residual-maker matrix that sweeps out (partials out) the linear influence of $\mathbf{X}_1$.

Theorem 2.2 (The Frisch–Waugh–Lovell Theorem)

In the partitioned OLS regression $\mathbf{Y} = \mathbf{X}_1 \hat{\boldsymbol{\beta}}_1 + \mathbf{X}_2 \hat{\boldsymbol{\beta}}_2 + \mathbf{e}$:

1. The subvector estimator $\hat{\boldsymbol{\beta}}_2$ obtained from the full multiple regression is algebraically identical to the estimator obtained by first regressing both $\mathbf{Y}$ and $\mathbf{X}_2$ on $\mathbf{X}_1$, extracting their residualized counterparts $\tilde{\mathbf{Y}} = \mathbf{M}_1 \mathbf{Y}$ and $\tilde{\mathbf{X}}_2 = \mathbf{M}_1 \mathbf{X}_2$, and then regressing $\tilde{\mathbf{Y}}$ on $\tilde{\mathbf{X}}_2$:

1. Furthermore, the residual vector $\tilde{\mathbf{e}}$ from the residualized regression of $\tilde{\mathbf{Y}}$ on $\tilde{\mathbf{X}}_2$ is identically equal to the residual vector $\mathbf{e} = (\mathbf{I}_n - \mathbf{H})\mathbf{Y}$ of the full multiple regression:

Proof via Block Gaussian Elimination:

From the first block of normal equations, $\mathbf{X}_1^T \mathbf{X}_1 \hat{\boldsymbol{\beta}}_1 + \mathbf{X}_1^T \mathbf{X}_2 \hat{\boldsymbol{\beta}}_2 = \mathbf{X}_1^T \mathbf{Y}$. Premultiplying by $(\mathbf{X}_1^T\mathbf{X}_1)^{-1}$ expresses $\hat{\boldsymbol{\beta}}_1$ in terms of $\hat{\boldsymbol{\beta}}_2$:

Substituting this expression for $\hat{\boldsymbol{\beta}}_1$ into the second block of normal equations $\mathbf{X}_2^T \mathbf{X}_1 \hat{\boldsymbol{\beta}}_1 + \mathbf{X}_2^T \mathbf{X}_2 \hat{\boldsymbol{\beta}}_2 = \mathbf{X}_2^T \mathbf{Y}$ yields:

Recognizing $\mathbf{H}_1 = \mathbf{X}_1(\mathbf{X}_1^T\mathbf{X}_1)^{-1}\mathbf{X}_1^T$ and collecting terms in $\hat{\boldsymbol{\beta}}_2$:

Because $\mathbf{M}_1$ is symmetric and idempotent ($\mathbf{M}_1 = \mathbf{M}_1^T \mathbf{M}_1$), we have $\mathbf{X}_2^T \mathbf{M}_1 \mathbf{X}_2 = (\mathbf{M}_1\mathbf{X}_2)^T(\mathbf{M}_1\mathbf{X}_2) = \tilde{\mathbf{X}}_2^T \tilde{\mathbf{X}}_2$ and $\mathbf{X}_2^T \mathbf{M}_1 \mathbf{Y} = \tilde{\mathbf{X}}_2^T \tilde{\mathbf{Y}}$. Moreover, substituting $\hat{\boldsymbol{\beta}}_1$ into the full residual vector $\mathbf{e} = \mathbf{Y} - \mathbf{X}_1\hat{\boldsymbol{\beta}}_1 - \mathbf{X}_2\hat{\boldsymbol{\beta}}_2$ gives:

Three Immediate Corollaries of the FWL Theorem:

1. Centering the Design Matrix: Take $\mathbf{X}_1 = \mathbf{1}_n$ (the intercept column). Then $\mathbf{M}_1 = \mathbf{I}_n - \frac{1}{n}\mathbf{1}_n\mathbf{1}_n^T$ centers every variable by subtracting its sample mean! Thus, the $k \times 1$ slope vector $\hat{\boldsymbol{\beta}}_{\text{slopes}}$ is given by $\hat{\boldsymbol{\beta}}_{\text{slopes}} = (\mathbf{X}_c^T \mathbf{X}_c)^{-1} \mathbf{X}_c^T \mathbf{y}_c$, where $\mathbf{X}_c$ is the $n \times k$ mean-centered predictor matrix.

2. Interpretation of Individual Partial Slopes $\hat{\beta}_j$: Taking $\mathbf{X}_2 = \mathbf{x}_{(j)}$ (a single predictor column) and $\mathbf{X}_1 = \mathbf{X}_{(-j)}$ (all other $p - 1$ columns), the scalar slope $\hat{\beta}_j$ is simply the bivariate OLS slope of $\mathbf{Y}$ (or $\tilde{\mathbf{Y}}_{(-j)}$) on the residual vector $\mathbf{e}_{j \mid (-j)} = \mathbf{M}_{(-j)} \mathbf{x}_{(j)}$:

That is, $\hat{\beta}_j$ measures the linear effect of the unique variation in $X_j$ that is orthogonal to all other predictors!

3. Orthogonal Predictors: If $\mathbf{X}_1^T \mathbf{X}_2 = \mathbf{0}_{p_1 \times p_2}$ (so the two predictor blocks are mutually orthogonal, as in balanced factorial designs), then $\mathbf{M}_1 \mathbf{X}_2 = \mathbf{X}_2$, and $\hat{\boldsymbol{\beta}}_2 = (\mathbf{X}_2^T\mathbf{X}_2)^{-1}\mathbf{X}_2^T\mathbf{Y}$ is completely unaffected by whether $\mathbf{X}_1$ is included or omitted.

Methodological Significance of the Frisch-Waugh-Lovell Theorem in Econometrics and Data Science

The Frisch–Waugh–Lovell Theorem is one of the most celebrated results in statistical methodology because it demystifies what "controlling for" or "holding constant" a set of covariates actually accomplishes algebraically. When an econometrician includes seasonal indicator variables or a linear time trend alongside economic regressors in a single multiple regression, the Frisch–Waugh–Lovell Theorem guarantees that the estimated slopes on the economic regressors are identical to those obtained by first deseasonalizing and detrending both the response series and each regressor series via Ordinary Least Squares and then regressing the detrended response on the detrended predictors. Similarly, in high-dimensional fixed-effects panel data regressions containing thousands of individual or firm-level dummy variables, inverting the full Gram matrix directly would be computationally prohibitive; applying the Frisch–Waugh–Lovell annihilator matrix sweeps out the group means in linear time ("within-group demeaning"), allowing exact estimation of the structural slope parameters.

In simple linear regression, there is only one explanatory variable, so the Regression Sum of Squares $\text{SSR}$ unambiguously belongs to $X_1$. In multiple linear regression with correlated predictors ($X_1, \dots, X_k$), the total explained sum of squares $\text{SSR}(X_1, \dots, X_k) = \hat{\boldsymbol{\beta}}^T \mathbf{X}^T \mathbf{Y} - n\bar{Y}^2$ cannot be decomposed into unique, order-independent additive components for each individual variable unless the predictors are mutually orthogonal. Instead, statistical theory defines Extra Sums of Squares that measure the marginal reduction in $\text{SSE}$ (or equivalent increase in $\text{SSR}$) when a subset of predictors is added to an existing baseline model:

By the Frisch–Waugh–Lovell Theorem and the orthogonal decomposition of projections $\mathcal{C}(\mathbf{X}_1, \mathbf{X}_2) = \mathcal{C}(\mathbf{X}_1) \oplus \mathcal{C}(\mathbf{M}_1 \mathbf{X}_2)$, the projection matrix of the full model decomposes additively as $\mathbf{H} = \mathbf{H}_1 + \mathbf{P}_{\tilde{\mathbf{X}}_2}$, where $\mathbf{P}_{\tilde{\mathbf{X}}_2} = \mathbf{M}_1\mathbf{X}_2(\mathbf{X}_2^T\mathbf{M}_1\mathbf{X}_2)^{-1}\mathbf{X}_2^T\mathbf{M}_1$, and the extra sum of squares is the quadratic form:

This distinction generates two primary sum-of-squares conventions reported in statistical software (such as R's anova() vs car::Anova() and SAS):

1. Sequential (Type I) Sums of Squares: Predictors are entered in a specified sequential order $X_1, X_2, \dots, X_k$. The sequential sums of squares partition the total model $\text{SSR}$ additively:

Type I sums of squares always add up to $\text{SSR}(X_1, \dots, X_k)$, but their values depend critically on the order in which the variables are listed!

2. Partial / Marginal (Type III) Sums of Squares: For each variable $X_j$, the Type III sum of squares $\text{SSR}(X_j \mid X_1, \dots, X_{j-1}, X_{j+1}, \dots, X_k)$ measures the marginal reduction in $\text{SSE}$ when $X_j$ is entered last into the model containing all other $k - 1$ predictors. Type III sums of squares are invariant to variable ordering and correspond directly to the individual partial $t$-tests ($F_j = \text{SSR}(X_j \mid \mathbf{X}_{(-j)}) / s^2 = t_j^2$), but they do not sum to $\text{SSR}(X_1, \dots, X_k)$ when predictors are correlated.

Coefficient of Partial Determination and Adjusted $R_a^2$

The Coefficient of Partial Determination $R_{Y j \cdot (-j)}^2$ measures the proportional reduction in the unexplained sum of squares $\text{SSE}(\mathbf{X}_{(-j)})$ achieved by adding predictor $X_j$ to a model that already contains all other covariates $\mathbf{X}_{(-j)}$:

Taking the signed square root (matching the sign of $\hat{\beta}_j$) gives the partial correlation coefficient $r_{Y j \cdot (-j)}$, which is identically the simple correlation between the FWL residuals $\mathbf{e}_{Y \mid (-j)}$ and $\mathbf{e}_{j \mid (-j)}$.

Finally, in multiple regression, the ordinary Coefficient of Multiple Determination $R^2 = 1 - \text{SSE}/\text{SST}$ is a monotonically non-decreasing function of the number of predictors $k$: adding any arbitrary noise variable $X_{k+1}$ cannot increase the minimized residual sum of squares $\text{SSE}$, so $R^2$ never decreases even when useless predictors are added. To penalize superfluous parameters, Henri Theil introduced the Adjusted Coefficient of Determination, $R_a^2$ (or $\bar{R}^2$), which replaces $\text{SSE}$ and $\text{SST}$ by their unbiased variance estimators $\text{MSE} = \text{SSE}/(n - p)$ and $\text{MST} = \text{SST}/(n - 1)$:

Algebraically, adding a single new predictor $X_j$ increases $R_a^2$ if and only if its partial $F$-statistic exceeds $1$ ($F_j = t_j^2 > 1$, or $|t_j| > 1$), because only when $F_j > 1$ does the reduction in $\text{SSE}$ more than compensate for the loss of one error degree of freedom.

Reconciling Sequential (Type I) and Marginal (Type III) Sums of Squares in Applied Modeling

Understanding why Sequential (Type I) and Partial (Type III) sums of squares differ when predictors are correlated prevents serious misinterpretations in applied data analysis. When two explanatory variables share substantial overlapping variation, the variable entered first in a Sequential (Type I) decomposition is credited with the entire shared variation in addition to its own unique variation, whereas the variable entered second is credited only with its incremental contribution beyond the first variable. Consequently, Type I tests are ideal when variables possess a natural hierarchical or causal ordering—such as polynomial powers where the linear term must precede the quadratic term, or nuisance baseline covariates that must be adjusted for prior to testing experimental treatments. Conversely, Type III partial sums of squares treat all predictors symmetrically by evaluating each variable as if it were entered last, making Type III tests the standard choice for assessing whether any single predictor can be dropped from a non-hierarchical multiple regression model.

In addition, comparing the unadjusted coefficient of multiple determination against Theil's adjusted coefficient of determination provides an essential diagnostic against overparameterization in small samples. When the number of columns in the design matrix approaches the sample size, the $p$-dimensional column space of the predictors spans nearly the entire $n$-dimensional observation space, forcing the unadjusted coefficient of determination toward unity even if every predictor is pure random noise independent of the response. Because the adjusted coefficient of determination divides both the residual sum of squares and the total sum of squares by their respective degrees of freedom, it penalizes the loss of residual degrees of freedom and can even become negative when a set of explanatory variables explains less variation than expected by chance alone.