5.10.2 Regression in matrix form

In this lesson

The model $\mathbf{y} = X\boldsymbol{\beta} + \boldsymbol{\varepsilon}$ , with a design matrix $X$ : a column of ones and one column per regressor
The estimate $\hat{\boldsymbol{\beta}} = (X^\top X)^{-1}X^\top\mathbf{y}$ from the normal equations of lesson 3.5.3
What a coefficient means when there are several regressors, and why leaving one out changes the others
Curves and groups: "linear" means linear in $\boldsymbol{\beta}$
The hat matrix $H$ , which turns $\mathbf{y}$ into the fitted values: the projection of lesson 3.5.2, and the leverage of each point
$R^2$ with several regressors, and why it never falls when a column is added
Checking every step in Python

You need simple regression from lesson 5.10.1, and least squares with the normal equations from lesson 3.5.3. The hat matrix is lesson 3.5.2's projection matrix.

The model in matrix form

Lesson 5.10.1 explained $y$ with one regressor $x$ . Most questions involve several: $x_1, x_2, \ldots$ , each measured on every observation. Write one equation per observation and stack them, as lesson 3.5.3 stacked the equations of a line.

Rule: Linear regression

\begin{gathered} \mathbf{y} = X\boldsymbol{\beta} + \boldsymbol{\varepsilon} \\[4pt] y_i = \beta_0 + \beta_1x_{i1} + \cdots \\ + \beta_px_{ip} + \varepsilon_i \end{gathered}

$\mathbf{y}$ holds the $n$ responses and $\boldsymbol{\varepsilon}$ the $n$ errors, each with mean 0 and variance $\sigma^2$ . The design matrix $X$ is $n \times k$ : a column of ones, then one column for each of the $p$ regressors, so $k = p + 1$ . $\boldsymbol{\beta} = (\beta_0, \ldots, \beta_p)$ holds the $k$ coefficients, in the order of the columns.

Simple regression is the case $p = 1$ : $X$ has the columns $\mathbf{1}$ and $\mathbf{x}$ , as in lesson 3.5.3. Here is an example with two regressors. The second, $x_2$ , is a dummy variable: 1 for the four points of a second group and 0 for the first.

$i$	1	2	3	4	5	6	7	8
$x_{i1}$	1	2	3	4	3	4	5	6
$x_{i2}$	0	0	0	0	1	1	1	1
$y_i$	3.5	6	5.5	7	1.5	1.5	2.5	4.5

Row $i$ of $X$ is $(1, x_{i1}, x_{i2})$ : the first row is $(1, 1, 0)$ and the last $(1, 6, 1)$ . So $X$ is $8 \times 3$ and $\boldsymbol{\beta} = (\beta_0, \beta_1, \beta_2)$ .

The normal equations

Least squares chooses $\hat{\boldsymbol{\beta}}$ to make $\lVert \mathbf{y} - X\mathbf{b} \rVert^2$ as small as possible over all $\mathbf{b}$ , the sum of squared residuals. That is lesson 3.5.3's problem with $A = X$ , so its answer carries over.

Rule: Least-squares estimate

\begin{gathered} X^\top X\,\hat{\boldsymbol{\beta}} = X^\top\mathbf{y} \\[4pt] \hat{\boldsymbol{\beta}} = (X^\top X)^{-1}X^\top\mathbf{y} \end{gathered}

The inverse exists when the columns of $X$ are independent (lesson 3.5.1). The residuals $\mathbf{e} = \mathbf{y} - X\hat{\boldsymbol{\beta}}$ are orthogonal to every column: $X^\top\mathbf{e} = \mathbf{0}$ . Compute $\hat{\boldsymbol{\beta}}$ with np.linalg.lstsq, not by inverting $X^\top X$ (lesson 3.5.3).

Each entry of $X^\top X$ is a dot product of two columns: the count of points, the regressors' sums, and their sums of squares and products. $X^\top\mathbf{y}$ holds the sum of the $y_i$ and the sums of each regressor times $y$ . For the eight points:

\begin{gathered} X^\top X = \begin{pmatrix} 8 & 28 & 4 \\ 28 & 116 & 18 \\ 4 & 18 & 4 \end{pmatrix} \\[6pt] X^\top\mathbf{y} = \begin{pmatrix} 32 \\ 110 \\ 10 \end{pmatrix} \end{gathered}

The 8 is the number of points, 28 the sum of $x_{i1}$ and 4 the number in the second group. The 18 is the sum of $x_{i1}$ over that group, and the last 4 is the sum of $x_{i2}^2$ , again the group's size.

Solving gives $\hat{\boldsymbol{\beta}} = (3, 1, -5)$ . The fitted values are $4, 5, 6, 7$ for the first group and $1, 2, 3, 4$ for the second. The residuals are $-0.5, 1, -0.5, 0, 0.5, -0.5, -0.5, 0.5$ , and their sum of squares is 2.5.

With $p = 1$ , the two normal equations are those of lesson 3.5.3's line, and their solution is lesson 5.10.1's $\hat{\beta}_1 = \frac{s_{xy}}{s_x^2}$ and $\hat{\beta}_0 = \bar{y} - \hat{\beta}_1\bar{x}$ .

X

import numpy as np x1 = np.array([1, 2, 3, 4, 3, 4, 5, 6.0]) x2 = np.array([0, 0, 0, 0, 1, 1, 1, 1.0]) y = np.array([3.5, 6, 5.5, 7, 1.5, 1.5, 2.5, 4.5]) X = np.column_stack([np.ones(8), x1, x2]) print(X.T @ X, X.T @ y) print(np.round(np.linalg.solve(X.T @ X, X.T @ y), 4)) beta, ss, rank, _ = np.linalg.lstsq(X, y, rcond=None) e = y - X @ beta print(np.round(beta, 4), rank) print(np.round(e, 4), np.round(X.T @ e, 10) + 0.0, round(e @ e, 4))

⌘/Ctrl + Enter runs the cell

What a coefficient means

With several regressors, each coefficient is the change in the mean of $y$ per unit of its own regressor, with the others held fixed.

Here $\hat{\beta}_1 = 1$ : within either group, one more unit of $x_1$ goes with one more unit of $y$ . And $\hat{\beta}_2 = -5$ : at the same $x_1$ , the second group sits 5 lower. The fit is two parallel lines, $y = 3 + x_1$ for the first group and $y = -2 + x_1$ for the second.

Inside each group $y$ rises one for one with $x_1$ . The second group has larger $x_1$ and lower $y$ , so the dashed line, fitted to $x_1$ alone, mixes the two effects and slopes down.

Leaving a regressor out

Fit $y$ on $x_1$ alone and the slope is $-\frac{1}{9} \approx -0.1111$ , with $R^2 = 0.0073$ . On its own, $x_1$ seems to do nothing, or even to lower $y$ . With $x_2$ in the model, $R^2$ is 0.918.

The reason is that $x_1$ and $x_2$ move together: their correlation is 0.6667, as the second group's $x_1$ values average 4.5 against 2.5. The slope of $x_2$ on $x_1$ is $\hat{\delta} = \frac{2}{9} \approx 0.2222$ ( $\delta$ is delta). Leaving $x_2$ out credits $x_1$ with part of $x_2$ 's effect.

Rule: Leaving a regressor out

\begin{gathered} \text{slope of } y \text{ on } x_1 \text{ alone} \\ = \hat{\beta}_1 + \hat{\beta}_2\,\hat{\delta} \end{gathered}

$\hat{\beta}_1$ and $\hat{\beta}_2$ come from the regression on both, and $\hat{\delta}$ is the slope of $x_2$ on $x_1$ . It holds exactly for least-squares fits. The slope alone equals $\hat{\beta}_1$ only when $\hat{\delta} = 0$ (the regressors are uncorrelated) or $\hat{\beta}_2 = 0$ .

Here $1 + (-5) \times \frac{2}{9} = -\frac{1}{9}$ . The gap is called omitted-variable bias. It is why a coefficient means nothing until you know which other regressors are in the model.

y

import numpy as np x1 = np.array([1, 2, 3, 4, 3, 4, 5, 6.0]) x2 = np.array([0, 0, 0, 0, 1, 1, 1, 1.0]) y = np.array([3.5, 6, 5.5, 7, 1.5, 1.5, 2.5, 4.5]) X = np.column_stack([np.ones(8), x1, x2]) A = np.column_stack([np.ones(8), x1]) def r2(M): e = y - M @ np.linalg.lstsq(M, y, rcond=None)[0] return 1 - (e @ e) / ((y - y.mean()) ** 2).sum() b = np.linalg.lstsq(X, y, rcond=None)[0] alone = np.linalg.lstsq(A, y, rcond=None)[0] print(np.round(alone, 4), round(r2(A), 4), round(r2(X), 4), ((y - y.mean()) ** 2).sum()) delta = np.linalg.lstsq(A, x2, rcond=None)[0][1] print(round(delta, 4), round(np.corrcoef(x1, x2)[0, 1], 4), x1[x2 == 0].mean(), x1[x2 == 1].mean(), round(b[1] + b[2] * delta, 4))

⌘/Ctrl + Enter runs the cell

Linear in the coefficients

"Linear" regression means that $\mathbf{y}$ is linear in $\boldsymbol{\beta}$ , not that it is a straight line in $x$ . A column may hold $x^2$ , $\ln x$ or a dummy, and least squares works the same.

A curve. The points $(-2, 5)$ , $(-1, 2)$ , $(0, 1)$ , $(1, 1)$ and $(2, 4)$ bend upwards. A straight line fits them badly: $R^2 = 0.0682$ . Add the column $x^2$ and the fit is $\hat{y} = 0.7429 - 0.3x + 0.9286x^2$ , with $R^2 = 0.9827$ .

np.polyfit(x, y, 2) gives the same coefficients, highest power first.

Groups. A dummy shifts the line for its group, as $x_2$ did. With $g$ groups and a column of ones, use $g - 1$ dummies. A dummy for every group adds up to the column of ones, so the columns are dependent (lesson 3.5.1).

For the eight points, the columns $\mathbf{1}$ , $x_2$ and $1 - x_2$ have rank 2. lstsq returns one of infinitely many $\hat{\boldsymbol{\beta}}$ , as lesson 3.5.3 found for dependent columns. Its fitted values, the two group means 5.5 and 2.5, are still unique.

x^2

import numpy as np x = np.array([-2, -1, 0, 1, 2.0]) y = np.array([5, 2, 1, 1, 4.0]) tss = ((y - y.mean()) ** 2).sum() for M in (np.column_stack([np.ones(5), x]), np.column_stack([np.ones(5), x, x**2])): b = np.linalg.lstsq(M, y, rcond=None)[0] e = y - M @ b print(np.round(b, 4), round(1 - (e @ e) / tss, 4)) print(np.round(np.polyfit(x, y, 2), 4)) x2 = np.array([0, 0, 0, 0, 1, 1, 1, 1.0]) y8 = np.array([3.5, 6, 5.5, 7, 1.5, 1.5, 2.5, 4.5]) D = np.column_stack([np.ones(8), x2, 1 - x2]) b, ss, rank, _ = np.linalg.lstsq(D, y8, rcond=None) print(rank, np.round(b, 4), np.round(D @ b, 4))

⌘/Ctrl + Enter runs the cell

The hat matrix

The fitted values are $\hat{\mathbf{y}} = X\hat{\boldsymbol{\beta}} = X(X^\top X)^{-1}X^\top\mathbf{y}$ . The matrix in front of $\mathbf{y}$ is the projection matrix of lesson 3.5.2 onto the column space of $X$ . Regression calls it the hat matrix $H$ , because it puts the hat on $\mathbf{y}$ .

Rule: The hat matrix

\begin{gathered} H = X(X^\top X)^{-1}X^\top \\[4pt] \hat{\mathbf{y}} = H\mathbf{y}, \quad \mathbf{e} = (I - H)\mathbf{y} \\[4pt] H^2 = H, \quad H^\top = H, \quad \operatorname{trace} H = k \end{gathered}

$H$ is $n \times n$ . The fitted values are the projection of $\mathbf{y}$ onto the column space, and the residuals are what is left, orthogonal to it. The trace, the sum of the diagonal entries, equals the number of columns $k$ .

For lesson 3.5.3's five points, $(0, 1)$ to $(4, 6)$ , with $X = [\mathbf{1}\ \mathbf{x}]$ :

H = \frac{1}{10}\begin{pmatrix} 6 & 4 & 2 & 0 & -2 \\ 4 & 3 & 2 & 1 & 0 \\ 2 & 2 & 2 & 2 & 2 \\ 0 & 1 & 2 & 3 & 4 \\ -2 & 0 & 2 & 4 & 6 \end{pmatrix}

$H\mathbf{y} = (1, 2.2, 3.4, 4.6, 5.8)$ , lesson 3.5.3's fitted values. Each fitted value is a weighted sum of all the $y_j$ , with row $i$ of $H$ as the weights. The middle row gives every point the weight 0.2: the fitted value at $\bar{x} = 2$ is the mean $\bar{y}$ .

Leverage

The diagonal entry $h_{ii}$ is the leverage of point $i$ : the weight of $y_i$ in its own fitted value. Each lies between 0 and 1, and they add up to $k$ , so they average $\frac{k}{n}$ .

For the five points the leverages are $0.6, 0.3, 0.2, 0.3, 0.6$ , adding up to 2. In simple regression $h_{ii} = \frac{1}{n} + \frac{(x_i - \bar{x})^2}{\sum_j (x_j - \bar{x})^2}$ : points far from $\bar{x}$ pull hardest.

Move the last point from $x = 4$ to $x = 12$ and its leverage becomes 0.9571. The line then has to pass close to it, whatever its $y$ . So check the leverages when one $x$ sits far from the rest: a large error there moves the whole fit.

H

import numpy as np x = np.array([0, 1, 2, 3, 4.0]) y = np.array([1, 2, 4, 4, 6.0]) X = np.column_stack([np.ones(5), x]) H = X @ np.linalg.solve(X.T @ X, X.T) print(np.round(10 * H, 10) + 0.0) print(np.allclose(H @ H, H), np.allclose(H, H.T), round(np.trace(H), 10)) print(np.round(H @ y, 4), np.round(y - H @ y, 4)) print(np.round(np.diag(H), 4), np.round(1 / 5 + (x - 2) ** 2 / 10, 4)) far = np.column_stack([np.ones(5), [0, 1, 2, 3, 12.0]]) print(np.round(np.diag(far @ np.linalg.solve(far.T @ far, far.T)), 4))

⌘/Ctrl + Enter runs the cell

R² with several regressors

$R^2 = 1 - \frac{\operatorname{RSS}}{\operatorname{TSS}}$ as in lesson 5.10.1, with $X$ including the column of ones. For the eight points, $\operatorname{RSS} = 2.5$ and $\operatorname{TSS} = 30.5$ , so $R^2 = 0.918$ .

Adding a column can never raise RSS. The column space only grows, so the nearest point in it can only come closer to $\mathbf{y}$ (lesson 3.5.2). So $R^2$ never falls, even when the new column is pure noise.

Lesson 3.5.3's 20 points have $R^2 = 0.8522$ . Add a column of 20 random numbers, unrelated to anything, and $R^2$ rises to 0.853. The adjusted $R^2$ charges for each column, and here it falls, from 0.844 to 0.8358.

Rule: Adjusted R squared

\bar{R}^2 = 1 - \frac{\operatorname{RSS}/(n - k)}{\operatorname{TSS}/(n - 1)}

Each sum is divided by its degrees of freedom: $n - k$ for the residuals and $n - 1$ for $y$ around its mean. Lesson 5.10.1's $s^2$ used $n - 2$ , the case $k = 2$ . A column that adds little lowers $\bar{R}^2$ ; lesson 5.10.3 tests whether regressors help.

R^2

import numpy as np rng = np.random.default_rng(353) x = np.arange(1, 21, dtype=float) y = 2 + 0.5 * x + rng.normal(0, 1, 20) noise = np.random.default_rng(5102).standard_normal(20) def fit(X): b = np.linalg.lstsq(X, y, rcond=None)[0] e = y - X @ b tss = ((y - y.mean()) ** 2).sum() n, k = X.shape return np.round(b, 4), round(float(1 - (e @ e) / tss), 4), round(float(1 - (e @ e / (n - k)) / (tss / (n - 1))), 4) print(*fit(np.column_stack([np.ones(20), x]))) print(*fit(np.column_stack([np.ones(20), x, noise])))

⌘/Ctrl + Enter runs the cell

Lesson 5.10.3 adds the statistics of $\hat{\boldsymbol{\beta}}$ : its standard errors, and tests of each coefficient and of all of them together. Case study 5.10.5 fits a stock's returns on two factors at once.

Common mistakes

Keep the column of ones. Without it the fit is forced through the origin, and $R^2$ as defined here no longer makes sense.
A coefficient depends on the other columns. $x_1$ 's slope is 1 with $x_2$ in the model and $-0.1111$ without it.
Use one dummy fewer than the number of groups. A dummy for every group, plus the ones, makes the columns dependent.
The order of $\hat{\boldsymbol{\beta}}$ follows the columns of $X$ . np.polyfit lists the highest power first.
Do not invert $X^\top X$ . Use lstsq and read its rank (lesson 3.5.3).
A higher $R^2$ is not proof that a new regressor matters. Even noise raises it; compare $\bar{R}^2$ or test it (lesson 5.10.3).
Look for high leverage. A point far from the other $x$ values can move the whole line.

Regression in matrix form

0 of 12 correct

1.
A regression has 50 observations, an intercept and 3 regressors. How many columns does $X$ have?
2.
$X$ has a column of ones and the column $\mathbf{x} = (1, 2, 3)$ . Find the entry in row 2, column 2 of $X^\top X$
3.
The normal equations are $X^\top X = \begin{pmatrix} 2 & 0 \\ 0 & 8 \end{pmatrix}$ and $X^\top\mathbf{y} = (6, 4)$ . Find $\hat{\beta}_1$ , the second entry of $\hat{\boldsymbol{\beta}}$
4.
A fit has $\hat{\boldsymbol{\beta}} = (3, 1, -5)$ for the columns $\mathbf{1}$ , $x_1$ , $x_2$ . Find the fitted value at $x_1 = 4$ , $x_2 = 1$
5.
With both regressors, $\hat{\beta}_1 = 2$ and $\hat{\beta}_2 = 3$ . The slope of $x_2$ on $x_1$ is 0.5. Find the slope of $y$ on $x_1$ alone
6.
With the same $\hat{\beta}_1 = 2$ and $\hat{\beta}_2 = 3$ , $x_2$ is uncorrelated with $x_1$ . Find the slope of $y$ on $x_1$ alone
7.
A regression with an intercept compares four groups. How many dummy columns should it have?
8.
$X$ is $30 \times 4$ with independent columns. Find the trace of the hat matrix $H$
9.
In simple regression with $n = 10$ , a point has $x_i - \bar{x} = 3$ and $\sum_j (x_j - \bar{x})^2 = 90$ . Find its leverage
10.
A regression has $n = 20$ and $k = 2$ . What is the average leverage of its points?
11.
A fit has $n = 25$ , $k = 5$ , $\operatorname{RSS} = 20$ and $\operatorname{TSS} = 100$ . Find the adjusted $R^2$
12.
A column is added to a regression, and $R^2$ was 0.6. Which is impossible for the new $R^2$ : 0.58, 0.6 or 0.65? Type the number

How to type answers

Powers with ^: x^7, 2^(1/3), a^-2.
Multiply with * or by writing side by side: 3x^2, 2(x+1).
ln(x) is the natural log, log(x, 2) is log base 2, e is Euler's number, sqrt(x) is a square root.
sin(x), cos(x) and tan(x) take radians; write pi for π and sin(x)^2 for sin²x.
Your first answer to each question is saved: it counts towards your score and brings the question back for review.
Exact answers like ln(20)/ln(3) are best. Decimals need about 4 significant figures.

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

Regression in matrix form

The model in matrix form

Rule: Linear regression

The normal equations

Rule: Least-squares estimate

What a coefficient means

Leaving a regressor out

Rule: Leaving a regressor out

Linear in the coefficients

The hat matrix

Rule: The hat matrix

Leverage

R² with several regressors

Rule: Adjusted R squared

Common mistakes

Regression in matrix form

Summary card

Free further reading