Linear regression from the statistics view

Before this lesson: lesson 04 (covariance of a random vector, $\operatorname{Cov}(AX) = A\operatorname{Cov}(X)A^\top$), lesson 05 (the multivariate Gaussian), lesson 08 (maximum likelihood), lesson 10 (priors and posteriors) and lesson 11 ($t$-tests and $p$-values). From linear algebra: least squares with QR and projections. Time: about 80 minutes, plus the notebook.

---

Why this matters #

• Is the effect real? Every regression table in a paper or an experiment report shows coefficients, standard errors and $p$-values. A/B tests are often analysed by regression, with covariates to reduce variance (CUPED in lesson 12).
• Uncertainty in predictions. A forecast needs an interval, and the right one is the prediction interval, not the confidence interval. Confusing them gives intervals that are far too narrow.
• Fits in ML research. Scaling laws are linear fits of log loss against log compute; the uncertainty of the fitted exponent decides how far you can extrapolate. Linear probes on network features are regressions too.
• Bayesian linear regression is the simplest Bayesian model with a closed-form posterior. It explains ridge and weight decay as a Gaussian prior, and it is the weight-space view of a Gaussian process (ml-for-science 07).

1. Predict first #

2. The picture #

The running example. Five students report hours of study $x$ and an exam score $y$ (out of 10):

$x$ (hours)	1	2	3	4	5
$y$ (score)	2	4	5	4	5

Least squares gives the line $\hat y = 2.2 + 0.6x$ (section 4 does every step). So each extra hour is worth 0.6 points. But five students are few. Another class of five would give another line. How much would the slope change?

The model. We assume that the data come from a true line plus noise:

where the errors $\varepsilon_1, \dots, \varepsilon_n$ are i.i.d. $\mathcal N(0, \sigma^2)$. The $x_i$ are fixed (we condition on them). Only the noise is random. The fitted coefficients $\hat\beta_0, \hat\beta_1$ are functions of the noisy $y_i$, so they are random variables. Their distribution over repeated data sets is the sampling distribution. You cannot see it from one data set, but you can simulate it:

Interactive figure: loading…

The fan is narrow in the middle of the data and wide at the ends. The fitted lines all pass close to the point $(\bar x, \bar y)$ and rotate around it. A small change of slope moves the line a lot far from $\bar x$ and very little near it. This one fact explains the shape of every band below.

Two questions, two bands. Take a new $x_0$.

• Where is the mean response $\mathbb{E}[y \mid x_0] = \beta_0 + \beta_1 x_0$? The only uncertainty is in the fitted line. This gives the confidence band.
• Where will one new observation $y_0$ fall? Now there are two sources: the uncertainty of the line, plus the new noise $\varepsilon_0$. This gives the wider prediction band.

Interactive figure: loading…

3. Definitions and main results #

The model and the least-squares estimator #

Write the data in matrix form. $y \in \mathbb{R}^n$ is the response, and $X \in \mathbb{R}^{n \times p}$ is the design matrix: row $i$ holds the features of observation $i$, and the first column is all ones (the intercept). In the exam example $p = 2$. The linear model is

We use these assumptions; each result says which ones it needs.

	Assumption
A1	linear mean: $\mathbb{E}[y] = X\beta$
A2	$X$ is fixed and has full column rank $p$ (no column is a combination of the others)
A3	errors have mean 0 and covariance $\sigma^2 I$: uncorrelated, with constant variance (homoscedastic)
A4	errors are Gaussian: $\varepsilon \sim \mathcal N(0, \sigma^2 I)$

The ordinary least squares (OLS) estimator minimises the residual sum of squares $\text{RSS}(\beta) = \|y - X\beta\|^2$. Setting the gradient $-2X^\top(y - X\beta)$ to zero gives the normal equations $X^\top X\hat\beta = X^\top y$, so

In code, never form the inverse: factor $X = QR$ and solve $R\hat\beta = Q^\top y$ (linear algebra 12). The fitted values are $\hat y = X\hat\beta = Hy$ with the hat matrix $H = X(X^\top X)^{-1}X^\top$. $H$ is the orthogonal projection onto the column space of $X$. The residuals are $e = y - \hat y = (I - H)y$. They are orthogonal to every column of $X$: $X^\top e = 0$. With an intercept column this says $\sum_i e_i = 0$.

OLS is the Gaussian maximum likelihood estimator. Under A4 the log-likelihood is $\ell(\beta, \sigma^2) = -\frac n2 \log(2\pi\sigma^2) - \frac{1}{2\sigma^2}\|y - X\beta\|^2$. For any $\sigma^2$, maximising over $\beta$ means minimising $\|y - X\beta\|^2$. So the MLE of $\beta$ is $\hat\beta$ (lesson 08). The MLE of $\sigma^2$ is $\text{RSS}/n$, which is biased; see below.

$\hat\beta$ is a random vector #

Theorem (sampling distribution of OLS). Under A1–A3, $\mathbb{E}[\hat\beta] = \beta$ and $\operatorname{Cov}(\hat\beta) = \sigma^2(X^\top X)^{-1}$. Under A4 in addition,

Proof.

1. Put $y = X\beta + \varepsilon$ into the formula. With $A = (X^\top X)^{-1}X^\top$ (a fixed $p \times n$ matrix): $\hat\beta = A X\beta + A\varepsilon = \beta + A\varepsilon$, because $AX = I$.
2. Mean: $\mathbb{E}[\hat\beta] = \beta + A\,\mathbb{E}[\varepsilon] = \beta$.
3. Covariance (lesson 04, $\operatorname{Cov}(A\varepsilon) = A\operatorname{Cov}(\varepsilon)A^\top$): $\operatorname{Cov}(\hat\beta) = A(\sigma^2 I)A^\top = \sigma^2 (X^\top X)^{-1}X^\top X(X^\top X)^{-1} = \sigma^2(X^\top X)^{-1}$.
4. Under A4, $\hat\beta$ is an affine function of the Gaussian vector $\varepsilon$, so it is Gaussian (lesson 05). $\blacksquare$

Simple regression. For one feature, $X^\top X = \begin{pmatrix} n & \sum x_i \\ \sum x_i & \sum x_i^2\end{pmatrix}$. Inverting it gives

Three levers make the slope precise: less noise (small $\sigma$), more points, and points spread far from $\bar x$ (large $S_{xx}$). If you repeat the same design $k$ times, $S_{xx}$ grows by $k$ and the standard error shrinks by $1/\sqrt k$.

Estimating $\sigma^2$: why $n - p$ #

$\sigma^2$ is unknown. The natural estimate uses the residuals:

Theorem. Under A1–A3, $\mathbb{E}[s^2] = \sigma^2$.

Why $n - p$ and not $n$ (intuition). The fitted line is chosen to be close to the data, so the residuals are smaller than the true errors. The residual vector is orthogonal to the $p$ columns of $X$, so it lives in a space of dimension $n - p$: fitting has used up $p$ directions. With $n = p$ points the fit is exact and all residuals are 0. The $n - p$ is called the residual degrees of freedom.

Proof: the trace of the hat matrix

1. $(I - H)X = X - X = 0$, so $e = (I - H)y = (I - H)(X\beta + \varepsilon) = (I - H)\varepsilon$.
2. $I - H$ is symmetric and idempotent ($(I - H)^2 = I - H$), so $\text{RSS} = e^\top e = \varepsilon^\top (I - H)\varepsilon$.
3. For a random vector with mean 0 and covariance $C$, $\mathbb{E}[\varepsilon^\top M\varepsilon] = \operatorname{tr}(MC)$ (write $\varepsilon^\top M\varepsilon = \operatorname{tr}(M\varepsilon\varepsilon^\top)$ and take the expectation). Here $C = \sigma^2 I$, so $\mathbb{E}[\text{RSS}] = \sigma^2 \operatorname{tr}(I - H)$.
4. $\operatorname{tr}(H) = \operatorname{tr}\big(X(X^\top X)^{-1}X^\top\big) = \operatorname{tr}\big((X^\top X)^{-1}X^\top X\big) = \operatorname{tr}(I_p) = p$ (the trace is invariant under cyclic shifts).
5. So $\mathbb{E}[\text{RSS}] = \sigma^2(n - p)$ and $\mathbb{E}[s^2] = \sigma^2$. $\blacksquare$

Under A4 one can say more: $\text{RSS}/\sigma^2 \sim \chi^2_{n-p}$, and $s^2$ is independent of $\hat\beta$. Independence follows from lesson 05: $\hat\beta - \beta = A\varepsilon$ and $e = (I - H)\varepsilon$ are jointly Gaussian, and $\operatorname{Cov}(A\varepsilon, (I - H)\varepsilon) = \sigma^2 A(I - H) = 0$ because $X^\top(I - H) = 0$.

Standard errors, $t$-tests and confidence intervals #

Replace $\sigma$ by $s$ in the covariance. The standard error of coefficient $j$ is

Theorem. Under A1–A4,

Proof.

1. $Z = (\hat\beta_j - \beta_j)/\big(\sigma\sqrt{[(X^\top X)^{-1}]_{jj}}\big) \sim \mathcal N(0, 1)$ by the theorem above.
2. $V = (n - p)s^2/\sigma^2 \sim \chi^2_{n-p}$, independent of $Z$ (end of the proof block above).
3. $T_j = Z/\sqrt{V/(n - p)}$ (the unknown $\sigma$ cancels). This ratio is the definition of Student's $t_{n-p}$ distribution (lesson 06). $\blacksquare$

This gives the usual tools (lesson 11):

• $t$-test of $H_0: \beta_j = 0$: the $t$-statistic is $t = \hat\beta_j/\operatorname{SE}_j$, and the two-sided $p$-value is $2P(T_{n-p} > |t|)$.
• Confidence interval at level 95%: $\hat\beta_j \pm t_{n-p,\,0.975}\operatorname{SE}_j$, where $t_{n-p,\,0.975}$ is the 0.975 quantile. For large $n - p$ it is close to 1.96; for $n - p = 3$ it is 3.18.

Confidence interval for the mean, prediction interval for a new $y$ #

Let $x_0 \in \mathbb{R}^p$ be a new feature vector (with the leading 1). The prediction is $\hat y_0 = x_0^\top\hat\beta$. Let $h_0 = x_0^\top (X^\top X)^{-1}x_0$.

• Mean response $\mu_0 = x_0^\top\beta$. $\operatorname{Var}(\hat y_0) = \sigma^2 h_0$, so the 95% confidence interval is $\hat y_0 \pm t_{n-p,\,0.975}\, s\sqrt{h_0}$.
• New observation $y_0 = x_0^\top\beta + \varepsilon_0$, with new noise $\varepsilon_0$ independent of the data. The error $y_0 - \hat y_0 = \varepsilon_0 - x_0^\top(\hat\beta - \beta)$ is a sum of two independent terms, so $\operatorname{Var}(y_0 - \hat y_0) = \sigma^2(1 + h_0)$. The 95% prediction interval is $\hat y_0 \pm t_{n-p,\,0.975}\, s\sqrt{1 + h_0}$.

For simple regression $h_0 = \frac1n + \frac{(x_0 - \bar x)^2}{S_{xx}}$: smallest at $x_0 = \bar x$ and growing quadratically away from it. That is the curved shape of the bands. As $n \to \infty$, $h_0 \to 0$: the confidence interval shrinks to a point, while the prediction interval tends to $x_0^\top\beta \pm 1.96\sigma$.

$R^2$, adjusted $R^2$ and the $F$-test #

With an intercept, the total sum of squares $\text{TSS} = \sum_i (y_i - \bar y)^2$ splits into the explained part $\sum_i(\hat y_i - \bar y)^2$ plus RSS (the cross term vanishes because $e$ is orthogonal to $\hat y - \bar y\mathbf 1$). The coefficient of determination is

the fraction of the variance of $y$ that the model explains. In simple regression $R^2 = r^2$, the squared sample correlation of $x$ and $y$. Adding a feature never increases RSS, so $R^2$ never goes down, even for a useless feature. The adjusted $R^2 = 1 - \frac{\text{RSS}/(n - p)}{\text{TSS}/(n - 1)}$ compares variance estimates instead, and can go down.

$F$-test for nested models. A small model with $p_0$ columns is inside a big model with $p_1$ columns (the small model sets $p_1 - p_0$ coefficients to 0). Under $H_0$ (the extra coefficients are 0) and A1–A4,

It asks: does the drop in RSS beat what noise alone would give? For one extra coefficient, $F = t^2$ and the two tests agree. The "overall $F$-test" in regression output compares the model with the intercept-only model.

Significance is not size. In simple regression $t^2 = (n - 2)\frac{R^2}{1 - R^2}$. With $n = 1000$ and $R^2 = 0.02$: $t = \sqrt{998 \cdot 0.02/0.98} = 4.51$ and $p \approx 7 \cdot 10^{-6}$. The slope is surely not 0, but $x$ explains only 2% of the variance of $y$. A small $p$-value says "there is an effect", not "the effect is large" or "the model predicts well".

The Gauss–Markov theorem #

Why least squares, and not another method? Call an estimator linear if it has the form $\tilde\beta = Cy$ for a fixed matrix $C$, and unbiased if $\mathbb{E}[\tilde\beta] = \beta$ for every $\beta$.

Theorem (Gauss–Markov). Under A1–A3 (no Gaussian assumption), OLS is the best linear unbiased estimator (BLUE): for every linear unbiased $\tilde\beta$, the matrix $\operatorname{Cov}(\tilde\beta) - \operatorname{Cov}(\hat\beta)$ is positive semidefinite. In particular, every coefficient and every combination $c^\top\beta$ has the smallest variance with OLS.

Proof

1. Write $C = A + D$ with $A = (X^\top X)^{-1}X^\top$. Unbiasedness for every $\beta$ means $CX\beta = \beta$ for all $\beta$, so $CX = I$. Since $AX = I$, this gives $DX = 0$.
2. $\operatorname{Cov}(\tilde\beta) = \sigma^2 CC^\top = \sigma^2(AA^\top + AD^\top + DA^\top + DD^\top)$.
3. $AD^\top = (X^\top X)^{-1}(DX)^\top = 0$, and $DA^\top = (AD^\top)^\top = 0$.
4. So $\operatorname{Cov}(\tilde\beta) = \sigma^2(X^\top X)^{-1} + \sigma^2 DD^\top$, and $DD^\top$ is positive semidefinite. $\blacksquare$

The theorem has limits. It compares only unbiased estimators: ridge regression is biased and can have a smaller mean squared error. It compares only linear ones: with heavy-tailed noise, a robust (nonlinear) fit can be much better. And it needs A3: with unequal variances, weighted least squares beats OLS.

Multicollinearity #

When two features are strongly correlated, the data cannot tell their effects apart. In the formula, $X^\top X$ is close to singular and $(X^\top X)^{-1}$ is large. A clean way to say it: with an intercept,

where $S_{x_jx_j} = \sum_i (x_{ij} - \bar x_j)^2$ and $R_j^2$ is the $R^2$ of a regression of feature $j$ on all the other features. The factor $\text{VIF}_j = 1/(1 - R_j^2)$ is the variance inflation factor. For two features with correlation $\rho$, $R_j^2 = \rho^2$: $\rho = 0.9$ gives VIF 5.3 (standard errors 2.3 times larger), and $\rho = 0.99$ gives VIF 50 (standard errors 7.1 times larger).

Interactive figure: loading…

Multicollinearity does not bias $\hat\beta$, and it does not hurt predictions at $x$ values like the training data. It makes single coefficients unstable, so do not read much into their signs or sizes. It is the statistical face of a large condition number (linear algebra 13).

Diagnostics: when the assumptions fail #

All the formulas above trust A1–A4. Residual plots check them: plot $e_i$ against the fitted values $\hat y_i$ (or against each feature, or against time). Under the model, this plot is a structureless band around 0.

Assumption fails	What the residual plot shows	What goes wrong	Common fix
A1: the truth is curved	a U shape or a wave	$\hat\beta$ is biased; predictions are wrong	add features ($x^2$, splines), transform $x$ or $y$
A3: variance not constant (heteroscedastic)	a fan: spread grows with $\hat y$	$\hat\beta$ is still unbiased, but the classical SEs are wrong (often too small)	robust standard errors; weighted least squares; $\log y$
A3: errors correlated (time series, groups)	waves in $e_i$ against time; clusters	SEs too small: the data carry less information than $n$ suggests	cluster or HAC standard errors; the block bootstrap (lesson 12)
A4: errors not Gaussian	skew or heavy tails (check a Q–Q plot)	exact $t$ and $F$ tests fail for small $n$; for large $n$ the CLT rescues them	larger $n$, the bootstrap, a robust fit
one point dominates	may show nothing (see below)	the fit follows one point	check leverage and influence

Leverage. The diagonal entries $h_{ii}$ of the hat matrix are the leverages: $\hat y_i = \sum_k h_{ik} y_k$, so $h_{ii}$ measures how much $y_i$ pulls its own fitted value. They lie in $[0, 1]$ and add up to $p$, so the average is $p/n$; a common flag is $h_{ii} > 2p/n$. In simple regression $h_{ii} = \frac1n + \frac{(x_i - \bar x)^2}{S_{xx}}$: points with an extreme $x$ have high leverage. A high-leverage point pulls the line toward itself, so its own residual is small, and the residual plot can hide it. Measures of influence such as Cook's distance combine leverage and the residual.

Interactive figure: loading…

Heteroscedasticity-robust (sandwich) standard errors

If $\operatorname{Cov}(\varepsilon) = \Omega$ is diagonal but not $\sigma^2 I$, step 3 of the proof above gives $\operatorname{Cov}(\hat\beta) = A\Omega A^\top = (X^\top X)^{-1}X^\top\Omega X(X^\top X)^{-1}$. The two outer factors are the "bread" and the middle is the "meat". White (1980) showed that you can replace $\Omega$ by $\operatorname{diag}(e_1^2, \dots, e_n^2)$ (the HC0 estimator) and get a consistent estimate of the covariance, even though each $e_i^2$ is a poor estimate of its own variance. The average over $n$ terms does the work. In small samples HC0 is too small; the HC3 variant uses $e_i^2/(1 - h_{ii})^2$ and is the usual default. Cluster-robust and HAC (Newey–West) estimators extend the same sandwich to correlated errors.

Correlation is not causation. A coefficient $\beta_j$ measures how $\mathbb{E}[y]$ changes with $x_j$ when the other included features are held fixed. If a variable that affects both $x_j$ and $y$ is left out (a confounder), $\hat\beta_j$ absorbs its effect. If the true model is $y = \beta_1 x_1 + \beta_2 x_2 + \varepsilon$ and you fit $y$ on $x_1$ alone, then $\mathbb{E}[\hat\beta_1] = \beta_1 + \beta_2\delta$, where $\delta$ is the slope of $x_2$ on $x_1$ (the omitted-variable bias). Small standard errors do not fix this; see causal inference 01.

Bayesian linear regression #

The Bayesian view (lesson 10) treats $\beta$ itself as random. Keep the Gaussian likelihood with a known noise variance $\sigma^2$, and put a Gaussian prior on the coefficients:

where $\alpha > 0$ is the prior precision (a large $\alpha$ says "the coefficients are small").

Theorem (Bayesian linear regression). The posterior is $\beta \mid y \sim \mathcal N(m_N, S_N)$ with

and the posterior predictive distribution of a new $y_0$ at $x_0$ is $\mathcal N\big(x_0^\top m_N,\ \sigma^2 + x_0^\top S_N x_0\big)$.

Proof.

1. Bayes: $\log p(\beta \mid y) = \log p(y \mid \beta) + \log p(\beta) + \text{const} = -\frac{1}{2\sigma^2}\|y - X\beta\|^2 - \frac\alpha2\|\beta\|^2 + \text{const}$.
2. Expand $\|y - X\beta\|^2 = y^\top y - 2\beta^\top X^\top y + \beta^\top X^\top X\beta$ and collect the terms in $\beta$: $\log p(\beta \mid y) = -\frac12\beta^\top\big(\alpha I + \tfrac{1}{\sigma^2}X^\top X\big)\beta + \beta^\top\tfrac{1}{\sigma^2}X^\top y + \text{const}$.
3. The log-density of $\mathcal N(m, S)$ is $-\frac12\beta^\top S^{-1}\beta + \beta^\top S^{-1}m + \text{const}$. Matching the quadratic terms gives $S_N^{-1}$; matching the linear terms gives $S_N^{-1}m_N = X^\top y/\sigma^2$.
4. Predictive: $y_0 = x_0^\top\beta + \varepsilon_0$, with $\beta \mid y \sim \mathcal N(m_N, S_N)$ and independent $\varepsilon_0 \sim \mathcal N(0, \sigma^2)$. The sum is Gaussian with mean $x_0^\top m_N$ and variance $x_0^\top S_N x_0 + \sigma^2$. $\blacksquare$

The posterior mean is ridge regression. Multiply inside and outside the inverse by $\sigma^2$:

the minimiser of $\|y - X\beta\|^2 + \lambda\|\beta\|^2$ with $\lambda = \alpha\sigma^2$. The posterior is Gaussian, so $m_N$ is also the MAP estimate. Weight decay in a neural network is the same idea: a Gaussian prior on the weights.

Two limits. With a flat prior ($\alpha \to 0$), $m_N \to \hat\beta$ (OLS) and $S_N \to \sigma^2(X^\top X)^{-1}$: the same matrix as the frequentist $\operatorname{Cov}(\hat\beta)$. With known $\sigma$ the 95% credible interval and the 95% confidence interval are then numerically equal, although they answer different questions. With much data, $S_N$ shrinks like $1/n$ and the predictive variance tends to the noise $\sigma^2$, just like the prediction interval.

Interactive figure: loading…

The function-space view: a Gaussian process with a linear kernel

Under the prior, $(\beta, y)$ is jointly Gaussian with $\operatorname{Cov}(y) = \alpha^{-1}XX^\top + \sigma^2 I$ and $\operatorname{Cov}(\beta, y) = \alpha^{-1}X^\top$. Gaussian conditioning (lesson 05) gives $m_N = \alpha^{-1}X^\top(\alpha^{-1}XX^\top + \sigma^2 I)^{-1}y$. The push-through identity $(\alpha I + X^\top X/\sigma^2)^{-1}X^\top/\sigma^2 = \alpha^{-1}X^\top(\alpha^{-1}XX^\top + \sigma^2 I)^{-1}$ shows that the two forms agree. The second form only needs inner products $x_i^\top x_k/\alpha$: this is a Gaussian process with the linear kernel $k(x, x') = x^\top x'/\alpha$. Replace the features by $\phi(x)$, or the kernel by any positive definite kernel, and you get GP regression (ml-for-science 07). The weight-space form solves a $p \times p$ system; the function-space form an $n \times n$ one.

If $\sigma^2$ is unknown too, a normal–inverse-gamma prior on $(\beta, \sigma^2)$ is conjugate. The predictive distribution is then a Student $t$, the Bayesian counterpart of the $t_{n-p}$ above.

Answers to section 1. (1) $4n$ points with the same design make $S_{xx}$ four times larger, so the standard error is divided by $\sqrt 4 = 2$. (2) Close to a fixed width: about $\pm 1.96\sigma$ (the noise of the new point does not shrink). (3) Yes. With $n = 1000$, $R^2 = 0.02$ gives $t = 4.51$ and $p \approx 7 \cdot 10^{-6}$; with larger $n$ the $p$-value is even smaller.

4. Worked example #

Back to the five students. Every number below is checked in section 5.

Step 1: the matrices. With the intercept column,

Step 2: the coefficients. $\det(X^\top X) = 5 \cdot 55 - 15^2 = 50$, so

Step 3: residuals. Fitted values $\hat y = 2.2 + 0.6x = (2.8, 3.4, 4.0, 4.6, 5.2)$, so $e = (-0.8, 0.6, 1.0, -0.6, -0.2)$. Check: $\sum e_i = 0$ and $\sum x_i e_i = -0.8 + 1.2 + 3.0 - 2.4 - 1.0 = 0$ (the residuals are orthogonal to both columns). $\text{RSS} = 0.64 + 0.36 + 1 + 0.36 + 0.04 = 2.4$.

Step 4: noise variance and standard errors. $n - p = 3$, so $s^2 = 2.4/3 = 0.8$. The estimated covariance is

so $\operatorname{SE}(\hat\beta_0) = \sqrt{0.88} = 0.938$ and $\operatorname{SE}(\hat\beta_1) = \sqrt{0.08} = 0.283$. (Check with the simple formula: $S_{xx} = 4 + 1 + 0 + 1 + 4 = 10$ and $\sqrt{0.8/10} = 0.283$.)

Step 5: test and interval for the slope. $t = 0.6/0.283 = 2.12$ with 3 degrees of freedom; the two-sided $p$-value is $0.124$. With $t_{3,\,0.975} = 3.182$, the 95% confidence interval is $0.6 \pm 3.182 \cdot 0.283 = 0.6 \pm 0.90 = (-0.30,\ 1.50)$. It contains 0: five students are not enough to show that studying helps, even though the estimate is positive.

Step 6: $R^2$ and $F$. $\bar y = 4$, so $\text{TSS} = 4 + 0 + 1 + 0 + 1 = 6$ and $R^2 = 1 - 2.4/6 = 0.6$. Adjusted $R^2 = 1 - (2.4/3)/(6/4) = 0.467$. The overall $F = \frac{(6 - 2.4)/1}{2.4/3} = 4.5 = t^2$, with the same $p = 0.124$.

Step 7: the two intervals at a new $x_0$. $h_0 = \frac15 + \frac{(x_0 - 3)^2}{10}$ and $s = \sqrt{0.8} = 0.894$.

$x_0$	$\hat y_0$	$h_0$	SE of the mean $s\sqrt{h_0}$	95% CI for the mean	SE for a new $y$: $s\sqrt{1 + h_0}$	95% prediction interval
3	4.0	0.2	0.40	$4.0 \pm 1.27 = (2.73, 5.27)$	0.98	$4.0 \pm 3.12 = (0.88, 7.12)$
6	5.8	1.1	0.94	$5.8 \pm 2.99 = (2.81, 8.79)$	1.30	$5.8 \pm 4.12 = (1.68, 9.92)$

At $x_0 = 6$ we extrapolate outside the data, and the confidence interval is more than twice as wide as at $\bar x = 3$. The prediction interval is wide in both places, because one student's score is noisy.

Step 8: the Bayesian version. Take one feature without intercept, prior precision $\alpha = 1$, noise variance $\sigma^2 = 1$, data $x = (1, 2)$, $y = (1, 3)$. Then $X^\top X = 1 + 4 = 5$ and $X^\top y = 1 + 6 = 7$.

• $S_N = 1/(\alpha + 5/\sigma^2) = 1/6$ and $m_N = S_N \cdot 7/\sigma^2 = 7/6 \approx 1.167$.
• Ridge with $\lambda = \alpha\sigma^2 = 1$: $7/(5 + 1) = 7/6$. Least squares gives $7/5 = 1.4$: the prior shrinks the slope toward 0.
• Predictive at $x_0$: mean $\frac76 x_0$, variance $1 + x_0^2/6$. At $x_0 = 1.5$: $1.375$, of which only $0.375$ is uncertainty about the slope. At $x_0 = 10$: $17.67$, of which 94% is uncertainty about the slope.

5. In code #

NumPy and SciPy check the worked example. scipy.stats.linregress is the 1-D cross-check.

import numpy as np
from scipy import stats

x = np.array([1., 2, 3, 4, 5])                 # hours of study
y = np.array([2., 4, 5, 4, 5])                 # exam score
X = np.column_stack([np.ones_like(x), x])      # intercept column + x
n, p = X.shape

Q, R = np.linalg.qr(X)                         # least squares through QR, no inverse of X^T X
beta = np.linalg.solve(R, Q.T @ y)
e = y - X @ beta
s2 = e @ e / (n - p)                           # unbiased estimate of sigma^2
XtX_inv = np.linalg.inv(R.T @ R)               # small p: fine for the standard errors
se = np.sqrt(s2 * np.diag(XtX_inv))
t = beta / se
pval = 2 * stats.t.sf(np.abs(t), df=n - p)
tq = stats.t.ppf(0.975, df=n - p)
print(beta.round(3), round(e @ e, 3), round(s2, 3))          # [2.2 0.6] 2.4 0.8
print(se.round(3), t.round(3), pval.round(3))                # [0.938 0.283] [2.345 2.121] [0.101 0.124]
print((beta[1] - tq * se[1]).round(2), (beta[1] + tq * se[1]).round(2))   # -0.3 1.5

R2 = 1 - e @ e / np.sum((y - y.mean())**2)
lr = stats.linregress(x, y)                    # the 1-D check
print(round(R2, 3), round(lr.slope, 3), round(lr.stderr, 3), round(lr.pvalue, 3))   # 0.6 0.6 0.283 0.124

for x0 in (3.0, 6.0):                          # half-widths: CI for the mean, prediction interval
    v = np.array([1, x0]) @ XtX_inv @ np.array([1, x0])
    print(x0, round(beta[0] + beta[1] * x0, 2), (tq * np.sqrt(s2 * v)).round(2), (tq * np.sqrt(s2 * (1 + v))).round(2))
# 3.0 4.0 1.27 3.12
# 6.0 5.8 2.99 4.12

# Bayesian linear regression = ridge with lambda = alpha * sigma^2
alpha, sigma2 = 1.0, 0.8
S_N = np.linalg.inv(alpha * np.eye(p) + X.T @ X / sigma2)
m_N = S_N @ X.T @ y / sigma2
ridge = np.linalg.solve(X.T @ X + alpha * sigma2 * np.eye(p), X.T @ y)
print(m_N.round(3), np.allclose(m_N, ridge))   # [1.277 0.839] True

The last line has a lesson of its own: this prior shrinks the intercept toward 0 too, and the slope grows to compensate. In practice one usually centres $y$ and the features, or leaves the intercept without a prior.

6. Where it shows up #

• Experiments and A/B tests. The treatment effect is the coefficient of a 0/1 treatment feature; adding pre-experiment covariates is regression adjustment (CUPED), which shrinks the standard error by about $\sqrt{1 - \rho^2}$ (lesson 12). A two-sample $t$-test is exactly a regression on a group indicator.
• Scaling laws and empirical fits in ML. A power law $L = aC^{-b}$ is a line in log–log coordinates. The standard error of the slope $b$ and the prediction interval tell you how far the fit can be trusted, and extrapolation is the case $x_0$ far from $\bar x$, where both bands are widest.
• Linear probes and attribution. Regressing a property on network activations is linear regression; with thousands of correlated features, multicollinearity makes single coefficients meaningless, so people use ridge.
• Bayesian last layers and uncertainty. Bayesian linear regression on learned features gives cheap predictive uncertainty; with random Fourier features it approximates a Gaussian process. Thompson sampling with a linear model samples $\beta$ from $\mathcal N(m_N, S_N)$ (reinforcement learning 03).
• Trends in science. A warming trend in temperature data is a regression slope. The residuals of climate and weather series are autocorrelated, so the classical standard error is too small and the trend looks more significant than it is; corrected (HAC) standard errors or an effective sample size are needed.
• Calibration of instruments. Fitting a sensor's response against known standards, then predicting the unknown quantity with an interval, is a textbook use of the prediction interval.

7. Common mistakes #

• Using the confidence interval for a single prediction. It covers the mean response only. A new observation needs the prediction interval, which includes $\sigma^2$.
• Dividing RSS by $n$. That is the MLE and it is biased low; use $n - p$, especially when $p$ is not small compared with $n$.
• Reading a small $p$-value as a large effect. With large $n$, tiny effects are significant. Report the estimate and its confidence interval, and look at $R^2$ or the size of the effect.
• Trusting classical standard errors without looking at the residuals. Curvature biases the coefficients; heteroscedasticity and correlated errors make the standard errors wrong. Plot the residuals first.
• Interpreting single coefficients under strong collinearity. Their signs can flip between similar data sets. Look at the VIF, or at combinations that are well determined.
• Reading a coefficient as a causal effect. It is an association given the included features. Omitted confounders bias it, whatever the standard error.
• Comparing $R^2$ across models with different numbers of features. $R^2$ always rises when features are added; use adjusted $R^2$, an $F$-test or held-out error.
• Forming $(X^\top X)^{-1}$ to fit the model. Solve with QR (or np.linalg.lstsq); the explicit inverse squares the condition number.

8. Check yourself #

1. Under A1–A3, show that the fitted values satisfy $\operatorname{Cov}(\hat y) = \sigma^2 H$ and the residuals $\operatorname{Cov}(e) = \sigma^2(I - H)$. What does this say about the variance of a residual with high leverage?
   Answer

   $\hat y = Hy = X\beta + H\varepsilon$, so $\operatorname{Cov}(\hat y) = H\sigma^2 IH^\top = \sigma^2 H$ ($H$ is symmetric and idempotent). $e = (I - H)\varepsilon$, so $\operatorname{Cov}(e) = \sigma^2(I - H)$ and $\operatorname{Var}(e_i) = \sigma^2(1 - h_{ii})$. A point with leverage close to 1 has a residual with variance close to 0: the line goes through it. This is why studentised residuals divide by $s\sqrt{1 - h_{ii}}$.
2. A regression on $n = 30$ points with $p = 4$ columns has $\text{RSS} = 52$. Estimate $\sigma$. A coefficient has $\hat\beta_j = 1.1$ and $[(X^\top X)^{-1}]_{jj} = 0.25$. Is it significant at the 5% level?
   Answer

   $s^2 = 52/26 = 2$, so $s = 1.414$. $\operatorname{SE}_j = 1.414 \cdot \sqrt{0.25} = 0.707$, $t = 1.1/0.707 = 1.56$ with 26 degrees of freedom. $t_{26,\,0.975} = 2.056 > 1.56$, so no ($p \approx 0.13$).
3. Why is the prediction band of simple regression curved (a hyperbola), and not two parallel lines?
   Answer

   The variance of the fitted line at $x_0$ is $\sigma^2\big(\frac1n + \frac{(x_0 - \bar x)^2}{S_{xx}}\big)$: the uncertain slope moves the line more the farther you are from $\bar x$. The half-width is the square root of a quadratic in $x_0$, which is a hyperbola. Only the noise term $\sigma^2$ is the same everywhere.
4. You fit $y$ on $x_1$ and $x_2$ with correlation $0.95$. Both $t$-statistics are about 1, but the overall $F$-test has $p < 10^{-6}$. Explain.
   Answer

   Each $t$-test asks whether one coefficient is 0 *given* the other feature. With correlated features either one alone explains $y$, so neither adds much given the other, and their SEs are inflated (VIF $\approx 10$). The $F$-test asks whether both are 0, and together they clearly explain $y$. The sum of the coefficients is well determined; the split between them is not.
5. In Bayesian linear regression, what happens to $m_N$ and $S_N$ when $\alpha \to \infty$? And when the data set grows?
   Answer

   $\alpha \to \infty$: $S_N \approx \alpha^{-1}I \to 0$ and $m_N \to 0$; the prior wins and the data are ignored. As $n$ grows (with well-spread $x$), $X^\top X$ grows like $n$, so $S_N \to 0$ like $1/n$ and $m_N \to$ the OLS estimate: the data win. The predictive variance tends to $\sigma^2$.
6. Why does $\hat\beta$ stay unbiased under heteroscedastic errors, while its standard error does not?
   Answer

   Unbiasedness uses only $\mathbb{E}[\varepsilon] = 0$: $\mathbb{E}[\hat\beta] = \beta + A\,\mathbb{E}[\varepsilon] = \beta$. The covariance uses $\operatorname{Cov}(\varepsilon) = \sigma^2 I$; with $\operatorname{Cov}(\varepsilon) = \Omega$ it becomes $A\Omega A^\top$, which the formula $s^2(X^\top X)^{-1}$ does not estimate. The sandwich estimator does.

Recall cardsAnswer in your head first, then reveal. Add the cards to your review deck to see them again just before you would forget.

1. Write the OLS estimator and its covariance under $\operatorname{Cov}(\varepsilon) = \sigma^2 I$.

   Show answer

2. What is the variance of the slope in simple regression?

   Show answer

3. Why does the unbiased noise estimate divide RSS by $n - p$?

   Show answer

4. What is the $t$-statistic for $H_0: \beta_j = 0$, and its null distribution?

   Show answer

5. How do the confidence interval for the mean and the prediction interval for a new $y$ at $x_0$ differ?

   Show answer

6. Define $R^2$.

   Show answer

7. State the Gauss–Markov theorem in one sentence.

   Show answer

8. Which ridge penalty gives the posterior mean of Bayesian linear regression with prior $\mathcal N(0, \alpha^{-1}I)$?

   Show answer

9. What is the predictive variance of Bayesian linear regression at $x_0$?

   Show answer

10. Two features have correlation 0.99. By how much are the standard errors of their coefficients inflated?

    Show answer

Add these cards to my reviewOpen the review →

9. Interview angle #

• "What are the assumptions of linear regression, and which ones matter for what?" — Linear mean (else biased coefficients); uncorrelated, constant-variance errors (else wrong standard errors, though $\hat\beta$ stays unbiased); Gaussian errors (only for exact small-sample $t$ and $F$ tests); full-rank $X$ (else $\hat\beta$ is not unique). Check them with residual plots.
• "Derive the variance of $\hat\beta$." — $\hat\beta = \beta + (X^\top X)^{-1}X^\top\varepsilon$, and $\operatorname{Cov}(A\varepsilon) = \sigma^2 AA^\top = \sigma^2(X^\top X)^{-1}$.
• "Confidence interval vs prediction interval?" — Mean response vs one new observation; the second adds $\sigma^2$ and does not shrink to zero.
• "Your coefficient has $p = 0.001$ but $R^2 = 0.01$. What do you tell the product team?" — There is very likely an association, but $x$ explains about 1% of the variation of $y$; report the effect size with its interval, and do not claim the model predicts well. Check for confounders before claiming a causal effect.
• "What is multicollinearity and how do you detect it?" — Correlated features inflate the variance of single coefficients ($\text{VIF} = 1/(1 - R_j^2)$); predictions are fine. Detect with VIFs or the condition number of $X$; fix with ridge, by combining features, or by dropping one.
• "Derive ridge regression from a Bayesian model." — Gaussian likelihood and prior $\mathcal N(0, \alpha^{-1}I)$; the negative log-posterior is $\frac{1}{2\sigma^2}\|y - X\beta\|^2 + \frac\alpha2\|\beta\|^2$, whose minimiser is ridge with $\lambda = \alpha\sigma^2$.
• "Why is OLS called BLUE?" — Gauss–Markov; also say what it does not claim (biased estimators like ridge can have lower MSE).

10. Practice #

The notebook for this lesson has five exercises, each with a check cell and a folded solution:

1. ols_table: coefficients, standard errors, $t$-values and $p$-values from scratch, checked against scipy.stats.linregress and against a simulation of the coefficient spread.
2. Coverage of the 95% confidence interval for the slope by simulation, with Gaussian and with heteroscedastic errors (classical and robust standard errors).
3. Coverage of the confidence interval for the mean and of the prediction interval for a new observation.
4. blr_posterior and blr_predictive, with the check that $m_N$ equals sklearn.linear_model.Ridge(alpha=α·σ², fit_intercept=False).
5. The diabetes data set (sklearn.datasets.load_diabetes): a coefficient table, VIFs and a residual plot.

11. Watch and read #

13. RegressionMIT OpenCourseWare · Philippe Rigollet, 18.650 Statistics for Applications · 1 h 16 min

The least-squares estimator as a projection, and its distribution under Gaussian noise.

14. Regression (cont.)MIT OpenCourseWare · Philippe Rigollet, 18.650 Statistics for Applications · 1 h 14 min

Tests and confidence intervals for the coefficients, with the $t$ and $\chi^2$ distributions.

Linear Regression, Clearly Explained!!!StatQuest with Josh Starmer · 27 min

A visual tour of least squares, $R^2$ and the $F$-test with small examples.

• Gareth James, Daniela Witten, Trevor Hastie and Robert Tibshirani, An Introduction to Statistical Learning, 2nd edition (free PDF): chapter 3 (standard errors, $t$-tests, $R^2$, the $F$-test, residual plots, outliers, leverage and collinearity).
• Christopher Bishop, Pattern Recognition and Machine Learning (free PDF): §3.1 (least squares as maximum likelihood) and §3.3 (Bayesian linear regression and the predictive distribution).
• Kevin Murphy, Probabilistic Machine Learning: An Introduction (free PDF): §11.2 (least squares and its sampling distribution) and §11.7 (Bayesian linear regression).
• Cosma Shalizi, The Truth About Linear Regression (free lecture notes): what each assumption buys, and what to do when it fails.