Simple & Multiple Linear Regression

Given paired data $(x_1, y_1), \ldots, (x_n, y_n)$ with the $x_i$ treated as fixed (non-random), the simple linear regression model asserts

where $\beta_0, \beta_1 \in \mathbb{R}$ are unknown parameters (the intercept and slope) and the errors $\epsilon_1, \ldots, \epsilon_n$ are random variables with $E[\epsilon_i] = 0$ and $\text{Var}(\epsilon_i) = \sigma^2$ for all $i$. The errors are assumed uncorrelated: $\text{Cov}(\epsilon_i, \epsilon_j) = 0$ for $i \neq j$.

Galton’s 1886 study of hereditary stature plotted adult-child height against mid-parent height and found the regression line’s slope was less than 1 — tall parents tended to have tall children, but on average less tall than the parents. He called this “regression toward mediocrity” (later softened to “regression toward the mean”). Figure 1 shows the Galton-style scatter; the fit line has slope $\approx 0.65$, cutting through the 45° line at the sample mean. The vertical segments are residuals — the part of each child’s height unexplained by the parental-height linear model.

The model $y_i = \beta_0 + \beta_1 x_i^2 + \beta_2 \log x_i + \epsilon_i$ is a linear regression: the response depends linearly on the unknown parameters $(\beta_0, \beta_1, \beta_2)$, even though $x_i^2$ and $\log x_i$ are nonlinear transformations of the predictor. The matrix form $\mathbf{y} = \mathbf{X}\boldsymbol\beta + \boldsymbol\epsilon$ of §21.4 accommodates any such design matrix. What this topic does not cover is nonlinear regression — models like $y_i = \beta_0 e^{-\beta_1 x_i} + \epsilon_i$ where $\boldsymbol\beta$ enters nonlinearly. Those require iterative fitting (Gauss–Newton, Levenberg–Marquardt) and belong to Topic 22 or later.

The word “regression” comes from Galton’s 1886 observation. A century earlier, Legendre (1805) and Gauss (1809, GAU1823) had introduced least squares for astronomical-orbit estimation. The two traditions — Galton’s statistical “regression toward the mean” and Gauss’s computational least-squares — merged into the modern discipline in the early 20th century, with Fisher’s 1922 foundational paper on estimation theory providing the bridge.

Given data $(x_1, y_1), \ldots, (x_n, y_n)$, the sum-of-squared-residuals loss is

The OLS estimator $(\hat\beta_0, \hat\beta_1)$ is the minimizer of $L$ over $\mathbb{R}^2$.

In the simple-regression setup, let $\hat{\mathbf{y}} = \hat\beta_0 \mathbf{1} + \hat\beta_1 \mathbf{x}$ denote the fitted values (vectors in $\mathbb{R}^n$). The OLS solution $(\hat\beta_0, \hat\beta_1)$ is the unique pair for which the residual vector $\mathbf{e} = \mathbf{y} - \hat{\mathbf{y}}$ satisfies

Equivalently, $\mathbf{e}$ is orthogonal to the column space $\text{col}(\mathbf{1}, \mathbf{x}) \subseteq \mathbb{R}^n$, and $\hat{\mathbf{y}}$ is the orthogonal projection of $\mathbf{y}$ onto that subspace. The OLS solution is unique whenever $\mathbf{x}$ is not constant.

Take the data $\mathbf{x} = (1, 2, 3, 4, 5)$, $\mathbf{y} = (1.2, 1.9, 3.3, 4.1, 5.0)$. Then $\bar x = 3$, $\bar y = 3.1$, $\sum (x_i - \bar x)^2 = 10$, $\sum (x_i - \bar x)(y_i - \bar y) = 9.6$, so

The fitted line $\hat y = 0.22 + 0.96 x$ produces residuals $(0.02, -0.24, 0.20, -0.04, 0.02)$, which sum to $-0.04$ (≈ 0, the small numerical residual from finite arithmetic) and have the signature of an OLS fit: one-half positive, one-half negative, approximately centered on zero.

Three reasons, in order of historical importance. (1) Analytic tractability: squared loss is differentiable and convex, so the minimizer has a closed form via the normal equations. Absolute-deviation loss is neither — LAD requires iterative linear programming or specialized algorithms. (2) Normal-errors connection: under Gaussian $\epsilon_i$, OLS coincides with the MLE (see §21.5 Rem 9). (3) Optimality certificate: Gauss–Markov (§21.6) shows OLS is BLUE under second-order assumptions — no competitor in the linear-unbiased class beats it on variance. The forward-pointer cost is sensitivity to outliers: squared loss magnifies large residuals quadratically, so one bad observation can swing the fit. Topic 15’s M-estimators (Huber, Tukey) trade some of Gauss–Markov’s optimality for robustness; the forward pointer in §21.9 Rem 20 states this explicitly.

If we add the stronger assumption $\epsilon_i \overset{\text{iid}}{\sim} \mathcal{N}(0, \sigma^2)$, the log-likelihood is

Maximizing over $(\beta_0, \beta_1)$ for any fixed $\sigma^2$ is identical to minimizing the sum of squared residuals — OLS and MLE pick the same $(\hat\beta_0, \hat\beta_1)$. Full treatment in §21.5 Rem 9 via Topic 14’s MLE machinery.

Given OLS estimates $(\hat\beta_0, \hat\beta_1)$, define

The coefficient of determination is $R^2 = 1 - \text{SSE}/\text{SST} = \text{SSR}/\text{SST}$, interpreted as the fraction of total variation in $y$ explained by the linear fit.

Under Definition 1 (fixed $x_i$, zero-mean uncorrelated errors with common variance $\sigma^2$),

Continuing Example 2 with $\sigma$ assumed known for illustration at $\sigma = 0.3$: $\sum (x_i - \bar x)^2 = 10$, so $\text{SE}(\hat\beta_1) = 0.3 / \sqrt{10} \approx 0.095$. A 95% Wald CI for $\beta_1$ is $0.96 \pm 1.96 \cdot 0.095 = (0.77, 1.15)$. In practice $\sigma$ is unknown and we plug in $\hat\sigma$ with a t-distribution critical value — see §21.8. The test–CI duality of Topic 19 tells us this CI is exactly the set of $\beta_1$ values a two-sided level-0.05 t-test would not reject.

$R^2$ is the fraction of $y$-variance captured by the linear fit. $R^2 = 1$ means perfect fit; $R^2 = 0$ means the fit is no better than the constant $\bar y$. Two caveats. First, $R^2$ always increases when predictors are added (the enlarged column space cannot shrink the projection’s accuracy), so it is useless as a model-comparison tool across models of different dimension. The adjusted $R^2$, $1 - (1 - R^2)(n - 1)/(n - p - 1)$, penalizes for predictor count but is still ad hoc; principled model selection lives in Topic 24 AIC/BIC/CV. Second, high $R^2$ does not imply the model is correct — the Anscombe quartet of Example 10 exhibits $R^2 \approx 0.67$ for four datasets with wildly different residual structures, only one of which is actually well-modeled by a line. Residual diagnostics (§21.9) are non-negotiable.

The sample Pearson correlation $r = \sum (x_i - \bar x)(y_i - \bar y) / \sqrt{\sum (x_i - \bar x)^2 \sum (y_i - \bar y)^2}$ is related to the slope by $\hat\beta_1 = r \cdot s_y / s_x$ and to $R^2$ by $r^2 = R^2$ (in simple regression). But correlation is symmetric ($r(x, y) = r(y, x)$) while regression is directional — the slope of $y$ on $x$ is not the reciprocal of the slope of $x$ on $y$ unless $r = \pm 1$. Regression asks “how does $E[Y \mid x]$ depend on $x$?”, a conditional-expectation question; correlation asks “how much do $x$ and $y$ co-vary?”, a joint-distribution question.

Given $n$ observations with $p$ non-intercept predictors, the design matrix $\mathbf{X}$ is $n \times (p + 1)$ with a leading column of ones and $i$-th row $(1, x_{i1}, \ldots, x_{ip})$. The parameter vector $\boldsymbol\beta = (\beta_0, \beta_1, \ldots, \beta_p)^\top \in \mathbb{R}^{p+1}$ collects intercept and slopes. The error vector $\boldsymbol\epsilon = (\epsilon_1, \ldots, \epsilon_n)^\top \in \mathbb{R}^n$ has $E[\boldsymbol\epsilon] = \mathbf{0}$ and $\text{Cov}(\boldsymbol\epsilon) = \sigma^2 \mathbf{I}_n$. The linear regression model is

Assuming $\mathbf{X}$ has full column rank $p + 1$, the OLS estimator minimizing $L(\boldsymbol\beta) = \|\mathbf{y} - \mathbf{X} \boldsymbol\beta\|^2$ is the unique solution to the normal equations

explicitly given by $\hat{\boldsymbol\beta} = (\mathbf{X}^\top \mathbf{X})^{-1} \mathbf{X}^\top \mathbf{y}$.

For simple regression ($p = 1$, so $\mathbf{X}$ has columns $(\mathbf{1}, \mathbf{x})$):

Solving $\mathbf{X}^\top \mathbf{X} \hat{\boldsymbol\beta} = \mathbf{X}^\top \mathbf{y}$ by hand (apply the $2 \times 2$ matrix inverse) reproduces

matching §21.3. The sufficient pair $(\mathbf{X}^\top \mathbf{X}, \mathbf{X}^\top \mathbf{y})$ is all the data tells OLS — no other aspect of $(\mathbf{x}, \mathbf{y})$ matters for the point estimate. This discharges the bullet in Topic 16 §16.12: OLS is a function of a Normal-linear-model sufficient statistic.

$(\mathbf{X}^\top \mathbf{X})^{-1}$ exists if and only if $\mathbf{X}$ has full column rank $p + 1$ — equivalently, no predictor is a linear combination of the others. When this fails (perfect collinearity), the normal equations have infinitely many solutions and OLS is undefined. In practice even near-collinearity (multicollinearity) destabilizes $\hat{\boldsymbol\beta}$: variances inflate, standard errors explode, individual coefficients become uninterpretable. The classical diagnostic is the variance inflation factor (VIF); the classical remedy is ridge regression (Topic 23). For Topic 21 we assume full column rank throughout.

When $(Y, \mathbf{X})$ are jointly multivariate Normal with mean $(\mu_Y, \boldsymbol\mu_X)$ and covariance $\begin{pmatrix} \sigma_Y^2 & \boldsymbol\Sigma_{YX} \\ \boldsymbol\Sigma_{XY} & \boldsymbol\Sigma_{XX} \end{pmatrix}$, Topic 8 §8.4 shows the conditional expectation is exactly linear:

So for jointly Normal data, multiple regression is not an approximation — it recovers the population regression line exactly. For non-Normal data, OLS is fitting the best linear approximation to $E[Y \mid \mathbf{X}]$ in mean-squared-error, which may or may not be the true conditional mean. The Gauss–Markov theorem (§21.6) handles the “best” part; Rem 11 under Normality promotes it to “unique UMVUE.”

The Normal linear model extends Definition 4 with the distributional assumption

i.e., the errors are iid $\mathcal{N}(0, \sigma^2)$. The design matrix $\mathbf{X}$ remains fixed (non-random).

Under Definition 5 (Normal linear model, full column rank),

Under the Normal linear model, let $\hat\sigma^2 = \text{SSE} / (n - p - 1)$. Then

The Normal log-likelihood is $\ell(\boldsymbol\beta, \sigma^2) = -\frac{n}{2} \log(2\pi\sigma^2) - \frac{1}{2\sigma^2} \|\mathbf{y} - \mathbf{X}\boldsymbol\beta\|^2$. For any $\sigma^2 > 0$, the $\boldsymbol\beta$ that maximizes $\ell$ is the one that minimizes $\|\mathbf{y} - \mathbf{X}\boldsymbol\beta\|^2$ — exactly the OLS problem. So $\hat{\boldsymbol\beta}_{\text{MLE}} = \hat{\boldsymbol\beta}_{\text{OLS}}$. The MLE of $\sigma^2$ is $\text{SSE}/n$ (biased downward by a factor $(n-p-1)/n$); the unbiased $\hat\sigma^2$ of Theorem 5 divides by $n - p - 1$ instead. See Topic 14 for the general MLE framework.

Combining Theorems 4 and 5: for each coefficient $j$,

The numerator is $\mathcal{N}(0, \sigma^2 [(\mathbf{X}^\top \mathbf{X})^{-1}]_{jj})$; the denominator is $\sqrt{\sigma^2 \chi^2_{n-p-1} / (n-p-1)}$, independent of the numerator. Their ratio is Student’s $t_{n-p-1}$. This pivotal quantity is the foundation of §21.8’s confidence intervals and hypothesis tests — the direct inheritor of Topic 17 §17.7’s one-sample t-test.

An estimator $\tilde{\boldsymbol\beta} = \mathbf{C} \mathbf{y}$ is a linear estimator of $\boldsymbol\beta$. It is unbiased if $E[\tilde{\boldsymbol\beta}] = \boldsymbol\beta$ for all $\boldsymbol\beta \in \mathbb{R}^{p+1}$ — equivalently, $\mathbf{C} \mathbf{X} = \mathbf{I}_{p+1}$. Among all linear unbiased estimators, the Best Linear Unbiased Estimator (BLUE) is the one with the smallest covariance matrix in the PSD sense: $\text{Cov}(\tilde{\boldsymbol\beta}) - \text{Cov}(\hat{\boldsymbol\beta}) \succeq \mathbf{0}$ for every competing $\tilde{\boldsymbol\beta}$.

Under Definition 4 (the linear regression model with $E[\boldsymbol\epsilon] = \mathbf{0}$, $\text{Cov}(\boldsymbol\epsilon) = \sigma^2 \mathbf{I}_n$, full column rank $\mathbf{X}$) — without any Normality assumption — the OLS estimator $\hat{\boldsymbol\beta} = (\mathbf{X}^\top \mathbf{X})^{-1} \mathbf{X}^\top \mathbf{y}$ is BLUE.

Take $n = 4$, $\mathbf{x} = (1, 2, 3, 4)$, $\boldsymbol\beta = (0, 1)^\top$, $\sigma = 1$. OLS gives $\text{Var}(\hat\beta_1) = 1/\sum(x_i - \bar x)^2 = 1/5 = 0.2$. A naive competitor — “use only the first three observations to fit a simple regression” — is also linear and unbiased, with $\text{Var}(\tilde\beta_1) = 1/\sum_{i=1}^3 (x_i - \bar x_3)^2 = 1/2 = 0.5$. Gauss–Markov guarantees OLS wins; the gap here is a factor of 2.5. Figure 6 shows the full MC comparison.

Gauss–Markov says OLS beats every linear unbiased competitor. Could a nonlinear unbiased estimator do better? Under Normal errors, Topic 16 Thm 4 (Lehmann–Scheffé) answers no: $(\mathbf{X}^\top \mathbf{X})^{-1} \mathbf{X}^\top \mathbf{y}$ is the unique UMVUE — minimum variance among all unbiased estimators, linear or not — because it is a function of a complete sufficient statistic for the Normal-linear-model exponential family. Without Normality, BLUE is the best we can claim; with Normality, BLUE = UMVUE. This is the §16.12 Track-6 promise discharged.

Gauss–Markov requires only: zero-mean errors, uncorrelated errors, common variance $\sigma^2$, full column rank $\mathbf{X}$. It does not require Normality, independence (only uncorrelatedness), or even identical distribution — errors can follow different distributions as long as they share variance. The one structural assumption that is load-bearing is homoscedasticity: all errors have the same variance. When this fails — heteroscedastic errors, $\text{Cov}(\boldsymbol\epsilon) = \sigma^2 \text{diag}(w_1, \ldots, w_n)$ with unequal $w_i$ — OLS is no longer BLUE. The weighted least squares (WLS) estimator reweights observations by $1/w_i$ and is BLUE in the heteroscedastic setting; see §21.9 Rem 21 for the one-sentence forward pointer to Topic 22.

Under full column rank, the hat matrix is

It “puts the hat on $\mathbf{y}$”: $\hat{\mathbf{y}} = \mathbf{H} \mathbf{y}$. The diagonal entries $h_{ii}$ are the leverage values of observation $i$ — measuring how much the $i$-th observation’s own response contributes to its own fit.

The hat matrix $\mathbf{H}$ is:

1. Symmetric: $\mathbf{H}^\top = \mathbf{H}$.
2. Idempotent: $\mathbf{H} \mathbf{H} = \mathbf{H}$.
3. Rank and trace: $\text{rank}(\mathbf{H}) = \text{tr}(\mathbf{H}) = p + 1$.
4. Complementary projector: $\mathbf{I} - \mathbf{H}$ is also symmetric idempotent, with $\text{rank}(\mathbf{I} - \mathbf{H}) = n - p - 1$.

Take $\mathbf{x} = (0, 1, 2)^\top$, so $\mathbf{X} = \begin{pmatrix} 1 & 0 \\ 1 & 1 \\ 1 & 2 \end{pmatrix}$. Then $\mathbf{X}^\top \mathbf{X} = \begin{pmatrix} 3 & 3 \\ 3 & 5 \end{pmatrix}$, $(\mathbf{X}^\top \mathbf{X})^{-1} = \tfrac{1}{6} \begin{pmatrix} 5 & -3 \\ -3 & 3 \end{pmatrix}$, and

Verify: row sums equal 1 (as they must when the first column of $\mathbf{X}$ is $\mathbf{1}$), trace $= (5 + 2 + 5)/6 = 2 = p + 1$. The diagonal entries $h_{11} = h_{33} = 5/6 \approx 0.83$ are high — the endpoint observations have the most leverage. The middle observation has $h_{22} = 2/6 \approx 0.33$, the lowest in this 3-point design.

Leverage $h_{ii} \in [0, 1]$ measures how “far” observation $i$‘s predictor vector $\mathbf{x}_i$ is from the center of the predictor cloud. Since $\sum_i h_{ii} = p + 1$, the average leverage is $(p + 1)/n$. The common rule of thumb flags observations with $h_{ii} > 2 (p + 1) / n$ as high-leverage. High leverage is not itself pathological — a well-designed experiment may deliberately include extreme $\mathbf{x}_i$ values to reduce $\text{Var}(\hat{\boldsymbol\beta})$ — but it is a prerequisite for being influential: an observation with low residual but high leverage contributes little to visible fit quality and much to the estimated coefficients. The diagnostic machinery of §21.9 combines leverage with studentized residuals via Cook’s distance.

The ANOVA identity $\text{SST} = \text{SSR} + \text{SSE}$ is the sample analog of Topic 4’s Eve’s law: $\text{Var}(Y) = E[\text{Var}(Y \mid \mathbf{X})] + \text{Var}(E[Y \mid \mathbf{X}])$. The fitted values $\hat y_i$ play the role of the conditional expectation $E[Y \mid \mathbf{X} = \mathbf{x}_i]$; the residuals $e_i$ play the role of the conditional deviation. The decomposition is exact in-sample (Pythagoras: $\|\mathbf{y} - \bar y \mathbf{1}\|^2 = \|\hat{\mathbf{y}} - \bar y \mathbf{1}\|^2 + \|\mathbf{y} - \hat{\mathbf{y}}\|^2$ when the intercept is included, because $\hat{\mathbf{y}} - \bar y \mathbf{1} \in \text{col}(\mathbf{X}) \cap \{\mathbf{1}\}^\perp$) and is the geometric basis of the one-way ANOVA example in §21.8.

For the $j$-th coefficient under the Normal linear model, the t-statistic for testing $H_0: \beta_j = \beta_j^0$ against the two-sided alternative is

Under the Normal linear model and $H_0: \beta_j = \beta_j^0$, $T_j \sim t_{n - p - 1}$. Hence a level-$\alpha$ two-sided test rejects when $|T_j| > t_{n - p - 1, 1 - \alpha/2}$, and the corresponding Wald-t confidence interval at level $1 - \alpha$ is

In a simple regression at $n = 30$, suppose $\hat\beta_1 = 0.45$ with $\widehat{\text{SE}}(\hat\beta_1) = 0.12$. Then $T_1 = 0.45 / 0.12 = 3.75$ on $t_{28}$. The two-sided p-value is $2 \cdot P(T_{28} > 3.75) \approx 8 \times 10^{-4}$, strongly rejecting the null at conventional levels. The 95% Wald-t CI is $0.45 \pm t_{28, 0.975} \cdot 0.12 = 0.45 \pm 2.048 \cdot 0.12 = (0.20, 0.70)$ — the interval excludes zero, consistent with the rejection.

Let $\mathcal{M}_1$ be the full Normal linear model with $p + 1$ parameters and $\mathcal{M}_0 \subset \mathcal{M}_1$ be the reduced model obtained by imposing $k$ linear restrictions $\mathbf{R} \boldsymbol\beta = \mathbf{r}$ (so $\mathcal{M}_0$ has $p + 1 - k$ free parameters). Let $\text{SSE}_1, \text{SSE}_0$ be the residual sums of squares under each. Under $H_0: \mathbf{R} \boldsymbol\beta = \mathbf{r}$ and Normal errors,

exactly (not just asymptotically).

Fit a four-predictor model $y_i = \beta_0 + \beta_1 x_{i1} + \beta_2 x_{i2} + \beta_3 x_{i3} + \beta_4 x_{i4} + \epsilon_i$ at $n = 50$. To test $H_0: \beta_2 = \beta_3 = 0$ (predictors 2 and 3 together add nothing), fit both the full and the reduced (two-predictor) models. With $k = 2$, $n - p - 1 = 45$, suppose $\text{SSE}_0 = 180$, $\text{SSE}_1 = 144$. Then $F = (36/2)/(144/45) = 18/3.2 = 5.625$ on $F_{2, 45}$. Since $F_{2, 45, 0.95} \approx 3.20$, we reject at $\alpha = 0.05$; the p-value is $\approx 0.0065$. Predictors 2 and 3 jointly contribute even though, as individual t-tests might show, neither alone is significant at 5%.

Three groups with $n_1 = n_2 = n_3 = 10$ observations each ($n = 30$). The full model is “group-specific means” — $y_{ij} = \mu_j + \epsilon_{ij}$ for $j \in \{1, 2, 3\}$ — fit by an all-dummy design matrix. The reduced model is “common mean” — $y_{ij} = \mu + \epsilon_{ij}$ — fit by a column of ones. The nested F-test has $k = 2$ (two group-mean differences), $n - p - 1 = 27$. Its $F$ statistic is exactly the classical ANOVA F, and its exact finite-sample distribution under $H_0: \mu_1 = \mu_2 = \mu_3$ is $F_{2, 27}$. Topic 4’s Eve’s law framing is the decomposition: