Multiple regression in matrix form · Probability & Statistics

One predictor is rarely enough. Exam scores depend on study hours and prior grades and sleep. Writing the model with p predictors in scalar notation gets unmanageable fast — in matrix form it stays as simple as the one-variable case.

The model

$Y_i = \beta_0 + \beta_1x_{i1} + \beta_2x_{i2} + \cdots + \beta_px_{ip} + \varepsilon_i$

Stack the observations:

$\mathbf{y} = \mathbf{X}\boldsymbol\beta + \boldsymbol\varepsilon$

$\underbrace{\begin{pmatrix}y_1\\y_2\\\vdots\\y_n\end{pmatrix}}_{n \times 1} = \underbrace{\begin{pmatrix}1 & x_{11} & \cdots & x_{1p}\\ 1 & x_{21} & \cdots & x_{2p}\\ \vdots & \vdots & & \vdots \\ 1 & x_{n1} & \cdots & x_{np}\end{pmatrix}}_{n \times (p+1)}\underbrace{\begin{pmatrix}\beta_0\\\beta_1\\\vdots\\\beta_p\end{pmatrix}}_{(p+1) \times 1} + \boldsymbol\varepsilon$

\mathbf{X} is the design matrix. Its first column of 1s is what produces the intercept — the intercept is just a coefficient on a constant predictor.

The least squares solution

Minimise the squared residuals, now a vector norm:

$S(\boldsymbol\beta) = \|\mathbf{y} - \mathbf{X}\boldsymbol\beta\|^2 = (\mathbf{y} - \mathbf{X}\boldsymbol\beta)^\top(\mathbf{y} - \mathbf{X}\boldsymbol\beta)$

Differentiate with respect to \boldsymbol\beta and set to zero:

$\frac{\partial S}{\partial\boldsymbol\beta} = -2\mathbf{X}^\top(\mathbf{y} - \mathbf{X}\boldsymbol\beta) = \mathbf{0}$

giving the normal equations $\mathbf{X}^\top\mathbf{X}\boldsymbol\beta = \mathbf{X}^\top\mathbf{y}$, and

$\boxed{\;\hat{\boldsymbol\beta} = (\mathbf{X}^\top\mathbf{X})^{-1}\mathbf{X}^\top\mathbf{y}\;}$

One formula, any number of predictors. With p = 1 it reduces exactly to S_{xy}/S_{xx} from the last lesson.

The geometry

This is the clearest way to understand least squares.

The fitted values are

$\hat{\mathbf{y}} = \mathbf{X}\hat{\boldsymbol\beta} = \underbrace{\mathbf{X}(\mathbf{X}^\top\mathbf{X})^{-1}\mathbf{X}^\top}_{\mathbf{H}}\mathbf{y}$

\mathbf{H} is the hat matrix (it puts the hat on \mathbf{y}), and it is a projection matrix: \mathbf{H}^2 = \mathbf{H} and $\mathbf{H}^\top = \mathbf{H}$.

So least squares is orthogonal projection of \mathbf{y} onto the column space of \mathbf{X} — the closest point in the space of everything the predictors can produce.

The residuals $\mathbf{e} = \mathbf{y} - \hat{\mathbf{y}} = (\mathbf{I} - \mathbf{H})\mathbf{y}$ are orthogonal to that space:

$\mathbf{X}^\top\mathbf{e} = \mathbf{0}$

Which says: the residuals are uncorrelated with every predictor, by construction. If they weren't, there'd be signal left to extract. That's the whole content of "least squares" in one equation.

Properties

Under the assumptions E[\boldsymbol\varepsilon] = \mathbf 0 and \operatorname{Cov}(\boldsymbol\varepsilon) = \sigma^2\mathbf{I}:

$E[\hat{\boldsymbol\beta}] = \boldsymbol\beta \qquad \text{(unbiased)}$

$\operatorname{Cov}(\hat{\boldsymbol\beta}) = \sigma^2(\mathbf{X}^\top\mathbf{X})^{-1}$

$\hat\sigma^2 = \frac{\mathbf{e}^\top\mathbf{e}}{n - p - 1}$

The n - p - 1 is n minus the number of estimated coefficients (including the intercept).

Gauss–Markov theorem: among all linear unbiased estimators, \hat{\boldsymbol\beta} has the smallest variance. It is BLUE — Best Linear Unbiased Estimator. Note the qualifiers: linear and unbiased. Dropping unbiasedness lets you do better in MSE, which is exactly what ridge regression does two lessons from now.

Standard errors come off the diagonal:

$\operatorname{SE}(\hat\beta_j) = \hat\sigma\sqrt{\big[(\mathbf{X}^\top\mathbf{X})^{-1}\big]_{jj}}, \qquad t_j = \frac{\hat\beta_j}{\operatorname{SE}(\hat\beta_j)} \sim t_{n-p-1}$

Interpreting coefficients

\beta_j is the expected change in y per unit increase in x_j, holding all other predictors fixed.

That last clause changes everything, and it is the single biggest difference from simple regression.

A predictor's coefficient in a multiple regression can differ in size — or even in sign — from its coefficient alone. Ice cream sales and drownings are positively related; add temperature as a predictor and the ice cream coefficient collapses toward zero, because temperature explains both.

Corollary: "holding others fixed" must be possible for the interpretation to be meaningful. If two predictors always move together in your data, no observation ever holds one fixed while the other varies, and the coefficients are estimated from almost nothing. That's multicollinearity — next lesson.

When \mathbf{X}^\top\mathbf{X} can't be inverted

The formula needs \mathbf{X} to have full column rank. It fails when:

p + 1 > n — more predictors than observations. There are infinitely many perfect fits.
Exact collinearity — one predictor is a linear combination of others. Including both "height in cm" and "height in inches" does it; so does the dummy variable trap, where you include an indicator for every category and an intercept, since the dummies sum to the constant column.

Near-collinearity doesn't break the inverse but makes it numerically unstable and the standard errors enormous.

In practice, never compute the inverse. numpy.linalg.lstsq and scipy.linalg.lstsq solve the system by QR or SVD, which is faster and far better conditioned than forming (\mathbf{X}^\top\mathbf{X})^{-1} explicitly.

Worked example

Predict exam score from hours studied (x_1) and prior GPA (x_2):

x_1	2	3	5	7	8
x_2	3.0	3.2	3.5	3.8	3.9
y	65	70	75	85	90

$\mathbf{X} = \begin{pmatrix}1 & 2 & 3.0\\ 1 & 3 & 3.2\\ 1 & 5 & 3.5\\ 1 & 7 & 3.8\\ 1 & 8 & 3.9\end{pmatrix}, \qquad \mathbf{y} = \begin{pmatrix}65\\70\\75\\85\\90\end{pmatrix}$

Solving $\hat{\boldsymbol\beta} = (\mathbf{X}^\top\mathbf{X})^{-1} \mathbf{X}^\top\mathbf{y}$ gives

$\hat y = 102.3 + 6.54x_1 - 16.67x_2$

Look at that GPA coefficient: it's negative. Taken at face value, the model says higher prior GPA lowers your exam score by 16.7 points per GPA point — which is nonsense, and the opposite of what a simple regression on GPA alone would say.

Now look at the standard errors:

term	estimate	std err	t	p
intercept	102.31	79.83	1.28	0.33
hours	6.54	4.40	1.49	0.28
gpa	−16.67	29.24	−0.57	0.63

They're enormous, and not one coefficient is significant — even though the model has R^2 = 0.988.

The reason: x_1 and x_2 have a correlation of 0.997. They carry essentially the same information, so the data cannot tell whether the effect belongs to study hours or to GPA. Any trade-off between them fits about equally well, including ones with absurd signs. The combination predicts beautifully; the individual coefficients are nearly unidentifiable.

This is the classic multicollinearity signature — a model that fits superbly with no significant predictors, and coefficients that flip sign — and it's why the next lesson exists.

Doing it in Python

The normal equations, three ways:

import numpy as np

x1 = np.array([2, 3, 5, 7, 8], dtype=float)
x2 = np.array([3.0, 3.2, 3.5, 3.8, 3.9])
y = np.array([65, 70, 75, 85, 90], dtype=float)

X = np.column_stack([np.ones(len(x1)), x1, x2])     # the design matrix

# 1. The textbook formula (never use this in production)
beta_inv = np.linalg.inv(X.T @ X) @ X.T @ y

# 2. Solve the normal equations without inverting
beta_solve = np.linalg.solve(X.T @ X, X.T @ y)

# 3. What you should actually use: least squares via SVD
beta_lstsq, *_ = np.linalg.lstsq(X, y, rcond=None)

print("inv formula :", beta_inv.round(4))
print("solve       :", beta_solve.round(4))
print("lstsq       :", beta_lstsq.round(4))
print("all agree   :", np.allclose(beta_inv, beta_lstsq))

The projection geometry, made concrete:

import numpy as np

x1 = np.array([2, 3, 5, 7, 8], dtype=float)
x2 = np.array([3.0, 3.2, 3.5, 3.8, 3.9])
y = np.array([65, 70, 75, 85, 90], dtype=float)
X = np.column_stack([np.ones(len(x1)), x1, x2])

H = X @ np.linalg.inv(X.T @ X) @ X.T          # the hat matrix
y_hat = H @ y
resid = y - y_hat

print("H is idempotent (H@H == H):", np.allclose(H @ H, H))
print("H is symmetric            :", np.allclose(H, H.T))
print("trace(H) =", round(np.trace(H), 4), "= number of parameters")
print()
print("residuals orthogonal to every predictor:")
print("  X.T @ e =", (X.T @ resid).round(10))
print("  so residuals are uncorrelated with the fitted values by construction")

Full inference — coefficients, standard errors and t statistics from scratch:

import numpy as np
from scipy.stats import t as tdist

x1 = np.array([2, 3, 5, 7, 8], dtype=float)
x2 = np.array([3.0, 3.2, 3.5, 3.8, 3.9])
y = np.array([65, 70, 75, 85, 90], dtype=float)
X = np.column_stack([np.ones(len(x1)), x1, x2])

n, k = X.shape
beta, *_ = np.linalg.lstsq(X, y, rcond=None)
resid = y - X @ beta
sigma2 = (resid @ resid) / (n - k)
cov = sigma2 * np.linalg.inv(X.T @ X)
se = np.sqrt(np.diag(cov))
tvals = beta / se
pvals = 2 * (1 - tdist.cdf(np.abs(tvals), df=n - k))

print(f"{'term':>10} {'estimate':>10} {'std err':>10} {'t':>8} {'p':>8}")
for name, b, s, tv, pv in zip(["intercept", "hours", "gpa"], beta, se, tvals, pvals):
    print(f"{name:>10} {b:>10.4f} {s:>10.4f} {tv:>8.3f} {pv:>8.4f}")

r2 = 1 - (resid @ resid) / ((y - y.mean()) ** 2).sum()
print(f"\nR^2 = {r2:.4f}")
print(f"correlation between the two predictors: {np.corrcoef(x1, x2)[0,1]:.4f}")
print("\nGreat fit, no significant predictors -- the multicollinearity signature.")

And the sign-flip that makes "holding others fixed" concrete:

import numpy as np

rng = np.random.default_rng(0)
n = 2_000
temperature = rng.normal(25, 8, n)
ice_cream = 2 * temperature + rng.normal(0, 5, n)
drownings = 0.5 * temperature + rng.normal(0, 3, n)      # ice cream has NO effect

# Simple regression: ice cream looks like it causes drownings
simple = np.polyfit(ice_cream, drownings, 1)[0]

# Multiple regression, controlling for temperature
X = np.column_stack([np.ones(n), ice_cream, temperature])
beta, *_ = np.linalg.lstsq(X, drownings, rcond=None)

print(f"ice cream coefficient, alone            : {simple:.4f}")
print(f"ice cream coefficient, given temperature: {beta[1]:.4f}")
print(f"temperature coefficient                 : {beta[2]:.4f}  (true value 0.5)")
print("\nThe apparent effect vanishes once the confounder is included.")

Your turn

1. n = 100 observations, 4 predictors plus an intercept. What are the dimensions of \mathbf{X}, \hat{\boldsymbol\beta}, and the degrees of freedom for \hat\sigma^2?

2. Why must \mathbf{X} have full column rank?

3. A predictor is significant alone but not in the multiple regression. What does that suggest?

Solutions

1. With n = 100 and p = 4 predictors plus an intercept:

\mathbf{X} is 100 \times 5 (the extra column is the 1s).
\hat{\boldsymbol\beta} is 5 \times 1.
df = n - p - 1 = 100 - 5 = 95.

Equivalently, df = n - \operatorname{rank}(\mathbf X), which is the version that stays correct when columns are collinear.

2. Because otherwise \mathbf{X}^\top\mathbf{X} is singular and cannot be inverted — there is no unique solution.

The deeper reason is identifiability. If column j is a linear combination of the others, then infinitely many \boldsymbol\beta vectors give exactly the same fitted values. The data cannot distinguish them, so "the coefficient on x_j" isn't a well-defined quantity.

Concretely, with height in both cm and inches, moving 1 unit of weight from one coefficient to the other (suitably scaled) changes nothing about the predictions. There's no fact of the matter about which one "the effect" belongs to.

Note the predictions remain unique even when the coefficients aren't — the projection onto the column space is still well defined. It's only the decomposition into individual contributions that fails.

3. Most likely multicollinearity — the predictor is strongly correlated with others in the model, so once they're included it has little unique information left to contribute. Its coefficient is estimated from the small part of its variation that's independent of the rest, which inflates the standard error.

Two other possibilities worth distinguishing:

Confounding, correctly handled. The simple relationship may have been spurious, and the multiple regression is right to kill it — the ice cream and drownings case. Here the change is a finding, not a problem.
Mediation. If x_j affects y through another included predictor, controlling for the mediator removes the path you were trying to measure. Statistically identical to the confounding case; causally the opposite conclusion.

The statistics alone cannot tell these apart. Check the correlation matrix and the VIFs (next lesson) for multicollinearity, and reason about the causal structure for the other two.

Check yourself in code

Fit the two-predictor model, confirm the hat matrix is a projection, and verify the residuals are orthogonal to the predictors.

Print exactly this:

beta [102.3077   6.5385 -16.6667]
H idempotent: True
residuals orthogonal: True
R squared 0.9881

Round the coefficients to 4 decimal places and R^2 to 4. Use numpy.linalg.lstsq. Treat orthogonality as satisfied if every entry of X.T @ resid is below 10^{-8} in absolute value.

import numpy as np

x1 = np.array([2, 3, 5, 7, 8], dtype=float)
x2 = np.array([3.0, 3.2, 3.5, 3.8, 3.9])
y = np.array([65, 70, 75, 85, 90], dtype=float)
X = np.column_stack([np.ones(len(x1)), x1, x2])

beta, *_ = np.linalg.lstsq(X, y, rcond=None)
print("beta", beta.round(4))

# Build the hat matrix, check it is idempotent, check the residuals are
# orthogonal to X, and report R squared.

import numpy as np

x1 = np.array([2, 3, 5, 7, 8], dtype=float)
x2 = np.array([3.0, 3.2, 3.5, 3.8, 3.9])
y = np.array([65, 70, 75, 85, 90], dtype=float)
X = np.column_stack([np.ones(len(x1)), x1, x2])

beta, *_ = np.linalg.lstsq(X, y, rcond=None)
print("beta", beta.round(4))

H = X @ np.linalg.inv(X.T @ X) @ X.T
print("H idempotent:", bool(np.allclose(H @ H, H)))

resid = y - X @ beta
print("residuals orthogonal:", bool(np.all(np.abs(X.T @ resid) < 1e-8)))

r2 = 1 - (resid @ resid) / ((y - y.mean()) ** 2).sum()
print("R squared", round(r2, 4))

In matrix form, multiple regression is one equation: $\hat{\boldsymbol\beta} = (\mathbf{X}^\top\mathbf{X})^{-1}\mathbf{X}^\top \mathbf{y}$, which is orthogonal projection of \mathbf y onto the column space of \mathbf X. Coefficients mean "per unit, holding the others fixed" — a clause that only makes sense when the predictors can actually vary independently.

Next: how to check whether the model's assumptions actually hold.

41. Multiple regression in matrix form