Study Notes

Linear Models

The organizing claim: regression is projection, and the whole course is the hat matrix H=X(XTX)−1XTH = X(X^TX)^{-1}X^T viewed from different angles — geometry (track 2 §2 cashed in), inference (track 3’s pivots on its diagonal), computation (QR/SVD/Cholesky), diagnostics (HiiH_{ii}), model selection (tr⁡H\operatorname{tr} H = df), and shrinkage (modify HH‘s eigenvalues). The second load-bearing idea: every multiple-regression coefficient is a simple regression on an adjusted variable — which explains both what coefficients mean and why correlated predictors ruin them.

1. The frame

  • The target is the regression function f(x)=E(y∣x)f(x) = E(y \mid x); the linear model is its first-order Taylor approximation — almost never the truth, almost always the right amount of model for the data you have. Watch n/pn/p; basis expansions buy flexibility while staying linear in β\beta.
  • The original “regression”: for bivariate normal (mother, daughter) heights with correlation ρ\rho,

E(HD∣Hm)=(1−ρ)μ+ρHmE(H_D \mid H_m) = (1-\rho)\mu + \rho H_m

— the prediction is shrunk toward the mean by 1−ρ1 - \rho. Regression to the mean is a property of conditioning, not a causal force.

  • Geometric view of the humble t-test: write z=zˉ1+(z−zˉ1)z = \bar{z}\mathbf{1} + (z - \bar{z}\mathbf{1}) = projection onto 1\mathbf{1} + residual. The t statistic is monotone in cos⁡θ\cos\theta between zz and 1\mathbf{1}:

{t>c}  ⟺  {n zˉ∥z∥>c0}  ⟺  {θ<θ0}\{t > c\} \iff \left\{\frac{\sqrt{n}\,\bar{z}}{\|z\|} > c_0\right\} \iff \{\theta < \theta_0\}

Under H0H_0, z/∥z∥z/\|z\| is uniform on the sphere — tests are statements about angles. This picture scales to everything below.

Core competency set

  • State f(x)=E(y∣x)f(x) = E(y|x) and the Taylor justification for linearity; the regression-to-the-mean formula and its non-causal reading.
  • Re-derive the t-test-as-angle picture.

2. The assumptions ladder

Least squares makes sense at three levels; each rung buys more:

RungAssumptionWhat you get
0nonepure approximation: β^\hat{\beta} minimizes ∥y−Xβ∥2\|y - X\beta\|^2, fit = projection
1Eε=0E\varepsilon = 0, Cov⁡(ε)=σ2I\operatorname{Cov}(\varepsilon) = \sigma^2I, XX fixedunbiasedness, Cov⁡(β^)=σ2(XTX)−1\operatorname{Cov}(\hat{\beta}) = \sigma^2(X^TX)^{-1}, Gauss–Markov
2ε∼N(0,σ2I)\varepsilon \sim N(0, \sigma^2I)exact distributions: β^\hat{\beta} Gaussian, t/F inference (track 3’s pivots)
  • Rung 0 is track 2 §2 verbatim: normal equations XTXβ^=XTyX^TX\hat{\beta} = X^Ty ⟺ residual ⊥ columns.
  • Rung 1’s variance, derived from the definition of covariance. The estimator is the normal-equations solution from rung 0, β^=(XTX)−1XTy\hat\beta = (X^TX)^{-1}X^Ty; substitute the rung-1 model y=Xβ+εy = X\beta + \varepsilon (so a true β\beta exists) and distribute:

β^=(XTX)−1XTy=(XTX)−1XT(Xβ+ε)=(XTX)−1XTX β+(XTX)−1XTε\hat\beta = (X^TX)^{-1}X^Ty = (X^TX)^{-1}X^T(X\beta + \varepsilon) = (X^TX)^{-1}X^TX\,\beta + (X^TX)^{-1}X^T\varepsilon

The first term collapses since (XTX)−1XTX=I(X^TX)^{-1}X^TX = I:

β^=β+(XTX)−1XTε\hat\beta = \beta + (X^TX)^{-1}X^T\varepsilon

This is a rung-1 identity, not a definition — β^\hat\beta splits into truth plus a linear function of the noise (at rung 0 there is no true β\beta and the split is meaningless). All the sampling variability lives in that second piece. So the estimation error is a fixed linear image of the noise, β^−β=Aε\hat\beta - \beta = A\varepsilon with A=(XTX)−1XTA = (X^TX)^{-1}X^T (XX is fixed, so AA is constant). The covariance of a linear image of a random vector is Cov⁡(Aε)=ACov⁡(ε)AT\operatorname{Cov}(A\varepsilon) = A\operatorname{Cov}(\varepsilon)A^T:

Cov⁡(β^)=A Cov⁡(ε) AT=(XTX)−1XT (σ2I) X(XTX)−1\operatorname{Cov}(\hat\beta) = A\,\operatorname{Cov}(\varepsilon)\,A^T = (X^TX)^{-1}X^T\,\big(\sigma^2 I\big)\,X(X^TX)^{-1}

Pull out the scalar σ2\sigma^2 and collapse the middle, using (XTX)−1(X^TX)^{-1} symmetric and XTX(XTX)−1=IX^TX(X^TX)^{-1} = I:

Cov⁡(β^)=σ2 (XTX)−1XTX(XTX)−1⏟I=σ2(XTX)−1\operatorname{Cov}(\hat\beta) = \sigma^2\,(X^TX)^{-1}\underbrace{X^TX(X^TX)^{-1}}_{I} = \sigma^2(X^TX)^{-1}

Simple-regression special case, to see what the matrix hides. The slope is a linear combination of the responses, β^1=∑iciyi\hat\beta_1 = \sum_i c_i y_i with ci=xi−xˉ∑j(xj−xˉ)2c_i = \frac{x_i - \bar x}{\sum_j(x_j - \bar x)^2}. The yiy_i are independent with variance σ2\sigma^2, so variance pulls through the squared weights:

Var⁡(β^1)=Var⁡(∑iciyi)=∑ici2 Var⁡(yi)=σ2∑ici2=σ2∥c∥2\operatorname{Var}(\hat\beta_1) = \operatorname{Var}\Big(\sum_i c_i y_i\Big) = \sum_i c_i^2\,\operatorname{Var}(y_i) = \sigma^2\sum_i c_i^2 = \sigma^2\|c\|^2

Evaluate ∥c∥2\|c\|^2 — the denominator is common to every cic_i, so it squares out and one power cancels against the numerator sum:

∥c∥2=∑i(xi−xˉ)2[∑j(xj−xˉ)2]2=∑i(xi−xˉ)2[∑j(xj−xˉ)2]2=1∑i(xi−xˉ)2\|c\|^2 = \sum_i \frac{(x_i - \bar x)^2}{\big[\sum_j(x_j - \bar x)^2\big]^2} = \frac{\sum_i(x_i - \bar x)^2}{\big[\sum_j(x_j - \bar x)^2\big]^2} = \frac{1}{\sum_i(x_i - \bar x)^2}

giving Var⁡(β^1)=σ2∑i(xi−xˉ)2\operatorname{Var}(\hat\beta_1) = \dfrac{\sigma^2}{\sum_i(x_i - \bar x)^2}. Spread-out xx‘s are leverage: the denominator is the design’s information about the slope. (Fisher information arithmetic from track 3, in design language.)

  • Rung 2: β^j∼N(βj,σ2vj)\hat{\beta}_j \sim N(\beta_j, \sigma^2v_j) with vj=[(XTX)−1]jjv_j = [(X^TX)^{-1}]_{jj}, s2=RSSn−ps^2 = \frac{RSS}{n-p} unbiased, and β^j−βjSE^∼tn−p\frac{\hat{\beta}_j - \beta_j}{\hat{SE}} \sim t_{n-p} — the track-1 pivot architecture assembled for coefficients.

Core competency set

  • Recite the ladder and what each rung buys; derive σ2(XTX)−1\sigma^2(X^TX)^{-1} and the simple-slope variance.
  • Say why spread in xx is the design’s information.

3. Estimability and Gauss–Markov

When XX is rank-deficient (r<pr < p), y^=Hy\hat{y} = Hy is still unique but β^\hat{\beta} is not (β^+c\hat{\beta} + c for any c∈N(X)c \in N(X) also solves).

  • A linear function cTβc^T\beta is estimable ⟺ cc lies in the row space of XX ⟺ c=aTXc = a^TX for some aa, and then aTya^Ty is a linear unbiased estimate. (One-way ANOVA: μ\mu and the αi\alpha_i are not estimable; the contrast α1−α2\alpha_1 - \alpha_2 is — differences of fits are unique even when parameters aren’t.)
  • Fixes, roughly equivalent: restrict to estimable functions; impose side conditions (∑αi=0\sum\alpha_i = 0 — must be non-estimable constraints, enough to reach rank pp without shrinking the model); reparametrize to full rank.
  • Identifiability is the underlying concept: β\beta identifiable ⟺ Xβ1=Xβ2⇒β1=β2X\beta_1 = X\beta_2 \Rightarrow \beta_1 = \beta_2 ⟺ XX full column rank. Don’t estimate what isn’t identified — you don’t know what you’re estimating. (n<pn < p: never identifiable without further structure.)

Gauss–Markov theorem. Under rung 1, among all linear unbiased estimates of an estimable cTβc^T\beta, the least-squares one cTβ^c^T\hat\beta has the smallest variance (BLUE). The proof is the orthogonal-split move (track 2 §6’s min-norm twin): write any competitor as the LS estimate plus an orthogonal piece, and show the piece only adds variance.

A linear estimate is aTya^Ty. Unbiasedness for every β\beta pins down a condition on aa:

E[aTy]=aTXβ=cTβ  for all β⇒aTX=cTE[a^Ty] = a^TX\beta = c^T\beta \ \text{ for all }\beta \quad\Rightarrow\quad a^TX = c^T

Its variance, since Cov⁡(y)=σ2I\operatorname{Cov}(y) = \sigma^2 I:

Var⁡(aTy)=aTCov⁡(y) a=σ2∥a∥2\operatorname{Var}(a^Ty) = a^T\operatorname{Cov}(y)\,a = \sigma^2\|a\|^2

Split aa orthogonally into its projection onto the column space C(X)C(X) and the remainder: a=a∗+(a−a∗)a = a^* + (a - a^*) with a∗∈C(X)a^*\in C(X) and (a−a∗)⊥C(X)(a - a^*)\perp C(X). Two facts make the split pay. First, the remainder is invisible to the constraint: (a−a∗)⊥C(X)(a - a^*)\perp C(X) means (a−a∗)TX=0(a - a^*)^TX = 0, so a∗TX=aTX=cTa^{*T}X = a^TX = c^T — the projection alone already satisfies unbiasedness, and a∗Ty=cTβ^a^{*T}y = c^T\hat\beta is exactly the least-squares estimate. Second, the remainder is orthogonal to the projection, so Pythagoras splits the norm:

∥a∥2=∥a∗∥2+∥a−a∗∥2\|a\|^2 = \|a^*\|^2 + \|a - a^*\|^2

Multiply by σ2\sigma^2 and drop the non-negative remainder term:

Var⁡(aTy)=σ2∥a∥2=σ2∥a∗∥2+σ2∥a−a∗∥2  ≥  σ2∥a∗∥2=Var⁡(cTβ^)\operatorname{Var}(a^Ty) = \sigma^2\|a\|^2 = \sigma^2\|a^*\|^2 + \sigma^2\|a - a^*\|^2 \;\geq\; \sigma^2\|a^*\|^2 = \operatorname{Var}(c^T\hat\beta)

with equality iff a=a∗a = a^*. The off-column-space part of any competitor is pure variance and buys no unbiasedness — same structure as the min-norm A+bA^+b argument (track 2 §6): split orthogonally, one part does the work, the other is pure cost.

Implications

  • BLUE is a linear, unbiased crown — biased estimators (ridge, §8) can and do beat it in MSE; this is track 3’s Cramér–Rao loophole, now concrete.
  • Estimability via SVD: Var⁡(vjTβ^)=σ2/dj2\operatorname{Var}(v_j^T\hat{\beta}) = \sigma^2/d_j^2 — contrasts along small singular directions are barely estimable; at dj=0d_j = 0, not at all. Estimability is quantitative, not binary.

Core competency set

  • Define estimable/identifiable; produce the ANOVA contrast example.
  • Run the Gauss–Markov projection proof and name its track-2 twin.

4. Computation: QR, and what MGS teaches about coefficients

  • Solve LS by orthogonalizing, never by inverting: ∥y−Xβ∥2=∥QTy−Rβ∥2\|y - X\beta\|^2 = \|Q^Ty - R\beta\|^2; with RR triangular, back-substitute R~β^=y1∗\tilde{R}\hat{\beta} = y_1^* and RSS=∥y2∗∥2RSS = \|y_2^*\|^2. Methods: Householder reflections (Q=I−2uuTQ = I - 2uu^T; map each column to ∥x∥e1\|x\|e_1), Givens rotations (zero one entry at a time — wins on sparse XX), modified Gram–Schmidt (orthogonalize against each new qq immediately — numerically safer than classical GS).
  • Never form (XTX)−1(X^TX)^{-1}: squaring XX squares its condition number κ=d1/dp\kappa = d_1/d_p (track 2 §6) — normal-equation methods need double precision to match QR. QR also hands you the Cholesky of XTXX^TX for free (XTX=R~TR~X^TX = \tilde{R}^T\tilde{R}).
  • Rank deficiency in practice: pivot columns by size; declare pseudo-rank when ∣Rjj∣/∣R11∣<τ|R_{jj}|/|R_{11}| < \tau (machine-epsilon scaled); standardize columns first.

The payoff theorem — what Gram–Schmidt means statistically:

β^j=⟨y, xj⋅rest⟩∥xj⋅rest∥2,Var⁡(β^j)=σ2∥xj⋅rest∥2\hat{\beta}_j = \frac{\langle y,\ x_{j\cdot\text{rest}}\rangle}{\|x_{j\cdot\text{rest}}\|^2}, \qquad \operatorname{Var}(\hat{\beta}_j) = \frac{\sigma^2}{\|x_{j\cdot\text{rest}}\|^2}

where xj⋅restx_{j\cdot\text{rest}} is xjx_j with the other predictors regressed out.

  • Each multiple-regression coefficient is a simple regression on the adjusted variable — “the effect of xjx_j after accounting for the others” is literal, not a figure of speech.
  • And the variance formula is the whole collinearity story: correlated predictors shrink ∥xj⋅rest∥\|x_{j\cdot\text{rest}}\|, blowing up Var⁡(β^j)\operatorname{Var}(\hat{\beta}_j). Orthogonal predictors ⟹ multiple = univariate.
  • Numbers — two standardized predictors with correlation ρ\rho:

(XTX)−1=11−ρ2[1−ρ−ρ1](X^TX)^{-1} = \frac{1}{1-\rho^2}\begin{bmatrix} 1 & -\rho \\ -\rho & 1 \end{bmatrix}

At ρ=0.9\rho = 0.9: variances inflated by 11−0.81≈5.3×\frac{1}{1 - 0.81} \approx 5.3\times, and positively correlated predictors get negatively correlated coefficients — the seesaw you see in real fits. (VIF is exactly this diagonal, per variable.)

Implications — what falls out of adjust-then-regress:

  • Omitted-variable bias has a formula: dropping a correlated predictor changes what the kept coefficients mean — they absorb the omitted variable’s effect through the correlation. “Effect of xjx_j” is always relative to what else was adjusted for.
  • Forward stepwise is just this move repeated: at each step, regress yy on each adjusted candidate and take the best — one MGS sweep per step (§7).
  • The standing counterexample — the suppressor: a predictor with coefficient ≈ 0 univariately can be strong and significant multivariately. Picture: two clusters offset by a binary x2x_2, each with a clean x1x_1–yy slope; the marginal regression must span both clumps and flattens to nothing, but adjusting for x2x_2 recovers the within-cluster slope. Univariate screening can discard jointly essential variables — and the reverse seesaw means jointly insignificant t-stats can hide a strong pair (drop both and the model collapses).

Core competency set

  • Name the three QR routes and when each wins; state why normal equations are numerically forbidden.
  • Recite the adjusted-variable theorem for β^j\hat{\beta}_j and its variance; run the ρ=0.9\rho = 0.9 numbers and the coefficient-seesaw implication.

5. The hat matrix: fits, residuals, leverage, influence

H=X(XTX)−1XTH = X(X^TX)^{-1}X^T; y^=Hy\hat{y} = Hy, r=(I−H)yr = (I-H)y. Symmetric idempotent (track 2’s projection test), tr⁡(H)=p\operatorname{tr}(H) = p (cyclic trace), eigenvalues = pp ones and n−pn-p zeros.

  • Covariances under rung 1: Cov⁡(y^)=σ2H\operatorname{Cov}(\hat{y}) = \sigma^2H, Cov⁡(r)=σ2(I−H)\operatorname{Cov}(r) = \sigma^2(I - H), Cov⁡(y^,r)=0\operatorname{Cov}(\hat{y}, r) = 0. Variance is conserved and split: a point’s variance goes to its fit or its residual — HiiH_{ii} near 1 means the fit chases the point and the residual goes blind to gross errors there.
  • Leverage is Mahalanobis distance. Split the design into the intercept column and the centered rest, X=[1 X~]X = [\mathbf 1\ \tilde X], where centering means 1TX~=0\mathbf 1^T\tilde X = 0 (each non-intercept column has mean zero). The cross-block vanishes, so XTXX^TX is block diagonal and inverts block-by-block:

XTX=[1T11TX~X~T1X~TX~]=[n00X~TX~]  ⇒  (XTX)−1=[1/n00(X~TX~)−1]X^TX = \begin{bmatrix} \mathbf 1^T\mathbf 1 & \mathbf 1^T\tilde X \\ \tilde X^T\mathbf 1 & \tilde X^T\tilde X \end{bmatrix} = \begin{bmatrix} n & 0 \\ 0 & \tilde X^T\tilde X \end{bmatrix} \;\Rightarrow\; (X^TX)^{-1} = \begin{bmatrix} 1/n & 0 \\ 0 & (\tilde X^T\tilde X)^{-1} \end{bmatrix}

Write observation ii‘s row as x⃗i=(1, x~i)\vec x_i = (1,\ \tilde x_i) and form the leverage; the block-diagonal inverse splits the quadratic form into the intercept piece and the centered piece:

Hii=x⃗iT(XTX)−1x⃗i=1⋅1n⋅1⏟intercept+x~iT(X~TX~)−1x~i=1n+x~iT(X~TX~)−1x~iH_{ii} = \vec x_i^T(X^TX)^{-1}\vec x_i = \underbrace{1\cdot\tfrac1n\cdot 1}_{\text{intercept}} + \tilde x_i^T(\tilde X^T\tilde X)^{-1}\tilde x_i = \frac1n + \tilde x_i^T(\tilde X^T\tilde X)^{-1}\tilde x_i

Factor out nn via the sample covariance Σ^=1nX~TX~\hat\Sigma = \frac1n\tilde X^T\tilde X, so that (X~TX~)−1=1nΣ^−1(\tilde X^T\tilde X)^{-1} = \frac1n\hat\Sigma^{-1}:

Hii=1n+1n x~iTΣ^−1x~i=1n(1+x~iTΣ^−1x~i)H_{ii} = \frac1n + \frac1n\,\tilde x_i^T\hat\Sigma^{-1}\tilde x_i = \frac{1}{n}\Big(1 + \tilde x_i^T\hat\Sigma^{-1}\tilde x_i\Big)

— the quadratic form x~iTΣ^−1x~i\tilde x_i^T\hat\Sigma^{-1}\tilde x_i is the Mahalanobis distance to the center of the xx-cloud.

— leverage is distance from the center of the xx-cloud in the metric of its shape. The picture to hold: an elongated cloud of xx-points with concentric covariance ellipses — a point at the tip of the long axis can sit far from xˉ\bar{x} with modest leverage, while a point the same Euclidean distance off the short axis has huge leverage. Rules of thumb: average is p/np/n; flag Hii>2p/nH_{ii} > 2p/n; dread >0.5> 0.5.

  • Leave-one-out is free. Removing observation ii is a rank-one downdate of XTXX^TX; Sherman–Morrison–Woodbury inverts it from the full inverse you already have:

(X(i)TX(i))−1=(XTX)−1+(XTX)−1xixiT(XTX)−11−Hii\big(X_{(i)}^TX_{(i)}\big)^{-1} = (X^TX)^{-1} + \frac{(X^TX)^{-1}x_ix_i^T(X^TX)^{-1}}{1 - H_{ii}}

Pushing this through the fit gives the two identities to keep:

y^i=(1−Hii)y^−i+Hiiyi⇒yi−y^−i=ri1−Hii\hat{y}_i = (1 - H_{ii})\hat{y}_{-i} + H_{ii}y_i \qquad\Rightarrow\qquad y_i - \hat{y}_{-i} = \frac{r_i}{1 - H_{ii}}

HiiH_{ii} is literally the weight the fit puts on a point’s own yy — self-influence. The jackknife residual exposes what the raw residual hides at high leverage.

  • Standardized residuals ri∗=riσ^1−Hiir_i^* = \frac{r_i}{\hat{\sigma}\sqrt{1 - H_{ii}}}; Cook’s distance packages leverage × outlyingness:

Di=∥y^−y^(i)∥2p σ^2=1p ri∗2 Hii1−HiiD_i = \frac{\|\hat{y} - \hat{y}_{(i)}\|^2}{p\,\hat\sigma^2} = \frac{1}{p}\,r_i^{*2}\,\frac{H_{ii}}{1 - H_{ii}}

Numbers: ri∗=2r_i^* = 2, Hii=0.4H_{ii} = 0.4, p=5p = 5: Di=15(4)0.40.6≈0.53D_i = \frac{1}{5}(4)\frac{0.4}{0.6} \approx 0.53 — a moderately outlying point at moderately high leverage already moves the whole fit by half a standard-error ball. High leverage is good when the point agrees with the model (free precision) and poison when it doesn’t — Cook’s distance is the arbiter.

  • Residual plots: structure in residuals vs fits/covariates = missing mean structure or non-constant variance; fix with transformations, WLS (§9), or basis expansion; test lack-of-fit by F against the saturated group-means model.

Implications — what falls out of the hat matrix:

  • LOOCV is closed-form for linear smoothers — the jackknife residual makes “refit nn times” a single fit; the same identity underlies GCV and fast influence diagnostics.
  • High-leverage gross errors are invisible in raw residuals: the fit chases the point (Var⁡(ri)=σ2(1−Hii)→0\operatorname{Var}(r_i) = \sigma^2(1 - H_{ii}) \to 0), so always read standardized or jackknife residuals, never raw ones, at high leverage.
  • Good leverage is free sample size: a high-HiiH_{ii} point consistent with the model tightens every coefficient — efficiency from design, the §2 spread-is-leverage fact at the level of single points.

Core competency set

  • Recite HH‘s properties and the covariance split; narrate “variance goes to fit or residual.”
  • Derive leverage-as-Mahalanobis (block split) and the LOO identity’s meaning; compute Cook’s distance from ri∗r_i^* and HiiH_{ii}.

6. The SVD view: PCA, collinearity, and where variance hides

X=UDVTX = UDV^T (track 2 §6 working at full tilt):

  • β^=VD−1UTy\hat{\beta} = VD^{-1}U^Ty — regression is trivial in the singular basis: project yy on each uju_j, divide by djd_j. Small djd_j ⟹ that coordinate’s noise is amplified by 1/dj1/d_j; rotating back by VV smears the variance across all coefficients. This is collinearity, exactly.
  • Var⁡(vjTβ^)=σ2/dj2\operatorname{Var}(v_j^T\hat{\beta}) = \sigma^2/d_j^2: the well-estimated contrasts live along big singular directions.
  • PCA = SVD of the centered XX: Σ^=1nX~TX~=1nVD2VT\hat{\Sigma} = \frac{1}{n}\tilde{X}^T\tilde{X} = \frac{1}{n}VD^2V^T — columns of VV = principal directions, variances dj2/nd_j^2/n. (Uncentered SVD finds directions through the origin, not directions of variance — centering is not optional.) Equivalent characterizations, same Lagrange/eigen computation: max-variance directions of Σ^\hat\Sigma ⟺ best-fitting affine hyperplane (closest rank-rr manifold, α0=xˉ\alpha_0 = \bar{x}).
  • Collinearity detection: scale columns to equal norm, then read the small singular values and their vjv_j‘s — large entries name the involved variables.
  • Total least squares (errors in XX too): stack [y:X][y : X], SVD, kill the last singular value — a “skew projection” for attribution problems; for pure prediction keep OLS (you’ll only have noisy XX at test time anyway).

Core competency set

  • Write β^=VD−1UTy\hat{\beta} = VD^{-1}U^Ty and tell the variance-amplification story.
  • State PCA = SVD-of-centered-X with both characterizations; explain why centering matters.

7. Model selection and the bias–variance ledger

  • The decomposition — add-and-subtract twice, the cross-terms dying for two different reasons. Setup: a new point (x,y)(x,y) with y=f(x)+εy = f(x) + \varepsilon, independent of the training set τ\tau that built f^\hat f; abbreviate f=f(x)f = f(x), f^=f^(x)\hat f = \hat f(x).

    First, add-and-subtract ff inside the prediction error and expand the square into three terms:

Err=EτEx,y[(y−f)+(f−f^)]2=E[(y−f)2]+2E[(y−f)(f−f^)]+E[(f−f^)2]\mathrm{Err} = E_\tau E_{x,y}\big[(y - f) + (f - \hat f)\big]^2 = E\big[(y - f)^2\big] + 2E\big[(y - f)(f - \hat f)\big] + E\big[(f - \hat f)^2\big]

The first term is E[ε2]=σ2E[\varepsilon^2] = \sigma^2. The cross-term dies because the new point’s noise ε=y−f\varepsilon = y - f is independent of f^\hat f (which saw only τ\tau) and has mean zero, so E[(y−f)(f−f^)]=E[ε] E[f−f^]=0E[(y-f)(f-\hat f)] = E[\varepsilon]\,E[f-\hat f] = 0. Left with:

Err=σ2+ExEτ(f−f^)2\mathrm{Err} = \sigma^2 + E_x E_\tau(f - \hat f)^2

Second, add-and-subtract the average fit Eτf^E_\tau\hat f inside the remaining term and expand:

Eτ(f−f^)2=Eτ[(f−Eτf^)+(Eτf^−f^)]2=(f−Eτf^)2+Eτ(f^−Eτf^)2E_\tau(f - \hat f)^2 = E_\tau\big[(f - E_\tau\hat f) + (E_\tau\hat f - \hat f)\big]^2 = (f - E_\tau\hat f)^2 + E_\tau(\hat f - E_\tau\hat f)^2

(this cross-term dies by definition of the mean: the factor f−Eτf^f - E_\tau\hat f is constant in τ\tau, and Eτ[Eτf^−f^]=0E_\tau[E_\tau\hat f - \hat f] = 0). Assembled:

Err=σ2⏟irreducible+Ex[(f−Eτf^)2⏟bias2+Eτ(f^−Eτf^)2⏟variance]\mathrm{Err} = \underbrace{\sigma^2}_{\text{irreducible}} + E_x\big[\underbrace{(f - E_\tau\hat f)^2}_{\text{bias}^2} + \underbrace{E_\tau(\hat f - E_\tau\hat f)^2}_{\text{variance}}\big]

  • Training error lies: RSS/nRSS/n is biased down for Err and monotone in complexity. The correction is priced by degrees of freedom:

df=1σ2∑iCov⁡(y^i,yi)df = \frac{1}{\sigma^2}\sum_i\operatorname{Cov}(\hat{y}_i, y_i)

— how much the fit eavesdrops on its own noise. For linear smoothers df=tr⁡(Hλ)df = \operatorname{tr}(H_\lambda) (OLS: pp; interpolation: nn). Cp=RSS/n+2 df σ2/nC_p = RSS/n + 2\,df\,\sigma^2/n restores unbiasedness; AIC is its likelihood twin.

  • Search procedures: best subset (combinatorial; aggressive — pays hidden df for the hunt, overfits noisy problems), forward stepwise (one MGS sweep finds the next variable: regress yy on each adjusted candidate), lasso (§8 — searches and shrinks at once). Use CV to pick the size λ\lambda, then refit the chosen size on all data; CV’s standard errors are optimistic.
  • Multiple-testing guard rails when reading coefficient t-stats: Bonferroni α/k\alpha/k; and remember the §4 seesaw — correlated predictors can all look insignificant while jointly essential.

Core competency set

  • Run the bias–variance derivation and name the move; define df by the covariance formula and compute it for OLS.
  • Explain why RSS can’t select models and what CpC_p/CV each fix.

8. Shrinkage: ridge and lasso

Ridge (min⁡∥y−Xβ∥2+λ∥β∥2\min \|y - X\beta\|^2 + \lambda\|\beta\|^2 — scale the columns first; not equivariant):

  • Closed form β^λ=(XTX+λI)−1XTy\hat{\beta}_\lambda = (X^TX + \lambda I)^{-1}X^Ty; substitute the SVD X=UDVTX = UDV^T and watch it diagonalize. Write each piece in the singular basis, using VTV=VVT=IV^TV = VV^T = I so that λI=λVVT\lambda I = \lambda VV^T:

XTX+λI=VD2VT+λVVT=V(D2+λI)VTX^TX + \lambda I = VD^2V^T + \lambda VV^T = V(D^2 + \lambda I)V^T

Invert (it’s V⋅diagonal⋅VTV\cdot\text{diagonal}\cdot V^T, so the inverse just inverts the diagonal block) and multiply by XT=VDUTX^T = VDU^T, collapsing VTV=IV^TV = I:

(XTX+λI)−1XT=V(D2+λI)−1VT VDUT=V(D2+λI)−1DUT(X^TX + \lambda I)^{-1}X^T = V(D^2 + \lambda I)^{-1}V^T\,VDU^T = V(D^2 + \lambda I)^{-1}DU^T

The middle factor (D2+λI)−1D(D^2 + \lambda I)^{-1}D is diagonal with entries djdj2+λ\frac{d_j}{d_j^2 + \lambda}, so applying to yy:

β^λ=V(D2+λI)−1DUTy=∑jvj djdj2+λ ujTy\hat{\beta}_\lambda = V(D^2 + \lambda I)^{-1}DU^Ty = \sum_j v_j\,\frac{d_j}{d_j^2 + \lambda}\,u_j^Ty

Ridge shrinks each principal coordinate smoothly by dj2dj2+λ\frac{d_j^2}{d_j^2 + \lambda}:

— big-dd directions barely touched, small-dd (high-variance) directions crushed. Ridge is targeted variance repair for exactly the directions §6 diagnosed. One SVD gives the whole λ\lambda-path. df(λ)=tr⁡(Hλ)=∑jdj2dj2+λdf(\lambda) = \operatorname{tr}(H_\lambda) = \sum_j\frac{d_j^2}{d_j^2 + \lambda}.

  • β^\hat\beta always lies in the row space of XX (compare min-norm, track 2).

Lasso (λ∥β∥1\lambda\|\beta\|_1): selection + shrinkage.

  • The stationarity conditions are the working object: on the active set AA, every active variable’s correlation with the residual is tied at λ\lambda, inactive ones below:

∣⟨xj,r⟩∣=λ (j∈A),∣⟨xj,r⟩∣≤λ (j∉A)|\langle x_j, r\rangle| = \lambda \ (j \in A), \qquad |\langle x_j, r\rangle| \leq \lambda \ (j \notin A)

Consequences: the path is piecewise linear in λ\lambda (LARS exploits this); the first variable to enter is arg⁡max⁡j∣⟨xj,y⟩∣\arg\max_j|\langle x_j, y\rangle| and λmax=max⁡j∣⟨xj,y⟩∣\lambda_{max} = \max_j|\langle x_j, y\rangle|.

  • Coordinate descent = soft thresholding. Take one standardized predictor (∥x∥=1\|x\| = 1, so its OLS coefficient is β^=xTy\hat\beta = x^Ty); the lasso objective in that one coordinate reduces to g(β)=12(β−β^)2+λ∣β∣g(\beta) = \frac{1}{2}(\beta - \hat\beta)^2 + \lambda|\beta|. Differentiate away from the kink, where ddβ∣β∣=sign⁡(β)\frac{d}{d\beta}|\beta| = \operatorname{sign}(\beta):

g′(β)=(β−β^)+λ sign⁡(β)=0g'(\beta) = (\beta - \hat\beta) + \lambda\,\operatorname{sign}(\beta) = 0

Solve in each sign branch and check the sign is consistent:

β>0: β=β^−λ  (valid iff β^>λ);β<0: β=β^+λ  (valid iff β^<−λ)\beta > 0:\ \beta = \hat\beta - \lambda\ \text{ (valid iff } \hat\beta > \lambda); \qquad \beta < 0:\ \beta = \hat\beta + \lambda\ \text{ (valid iff } \hat\beta < -\lambda)

When ∣β^∣≤λ|\hat\beta| \leq \lambda neither branch is consistent, so the minimum sits at the non-differentiable kink β=0\beta = 0. Assembled into one formula:

β^(λ)=S(β^,λ)=sign(β^) (∣β^∣−λ)+\hat{\beta}(\lambda) = S(\hat{\beta}, \lambda) = \mathrm{sign}(\hat{\beta})\,(|\hat{\beta}| - \lambda)_+

Numbers: β^=3,λ=1⇒2\hat{\beta} = 3, \lambda = 1 \Rightarrow 2; β^=0.8,λ=1⇒0\hat{\beta} = 0.8, \lambda = 1 \Rightarrow 0 — shrink the survivors, kill the small.

— slide toward zero by λ\lambda, clip at zero. The picture to hold: the soft-threshold function — a 45° line with a flat dead zone of width 2λ2\lambda at the origin — vs ridge’s uniformly tilted line (shrink everything, zero nothing) and best subset’s hard threshold (keep or kill). Multiple predictors: cycle, each time soft-thresholding the simple coefficient on the partial residual. Warm starts down a λ\lambda-grid + strong rules ⟹ the fast path algorithm everyone actually uses.

  • df(lasso) ≈ ∣A(λ)∣|A(\lambda)| — remarkably, the post-search shrinkage pays back the df spent hunting (best subset doesn’t; that’s why aggressive selection loses on noisy data). Relaxed lasso: refit (or blend) OLS on the active set to undo over-shrinkage.

PCR vs PLS vs ridge: PCR hard-thresholds principal components, ridge shrinks them smoothly, PLS lets yy pick directions (maximizes Var⁡(Xγ)Cor⁡2(Xγ,y)\operatorname{Var}(X\gamma)\operatorname{Cor}^2(X\gamma, y)). Ridge is usually the well-behaved default.

Implications — what falls out of shrinkage:

  • Every penalty is a prior (track 7 §3 makes it exact): ridge ↔ Gaussian, lasso ↔ Laplace. Choosing a regularizer is choosing what you believe about coefficients.
  • The df ledger explains the selection-method ranking on noisy data: best subset spends df hunting and keeps unshrunk coefficients (high variance); lasso hunts but then shrinks, paying the hunt back — which is why df(lasso) ≈ ∣A∣|A| and why it beats subset selection away from the low-noise regime.
  • Ridge never zeroes, lasso zeroes exactly (the dead zone) — use ridge for correlated groups you want averaged, lasso when you need a short list; relaxed lasso splits the difference by refitting the survivors.

Core competency set

  • Write ridge’s SVD shrinkage factor and its df; explain “targeted variance repair.”
  • State the lasso stationarity conditions and read off path-linearity and λmax\lambda_{max}; derive soft thresholding.
  • Compare df accounting: subset vs lasso vs ridge.

9. Weighted least squares and robustness

Correlated/heteroscedastic errors (Cov⁡(ε)=C\operatorname{Cov}(\varepsilon) = C): whiten, then OLS. Factor C=BBTC = BB^T (Cholesky, track 2 §5) and set A=B−1A = B^{-1}, so the transformed errors are white:

Cov⁡(Aε)=ACAT=B−1(BBT)B−T=I\operatorname{Cov}(A\varepsilon) = ACA^T = B^{-1}(BB^T)B^{-T} = I

Apply AA to the model and run OLS on (Ay, AX)(Ay,\ AX). The transformed objective is a C−1C^{-1}-weighted version of the original, via ATA=B−TB−1=(BBT)−1=C−1A^TA = B^{-T}B^{-1} = (BB^T)^{-1} = C^{-1}:

∥Ay−AXβ∥2=(y−Xβ)TATA(y−Xβ)=(y−Xβ)TC−1(y−Xβ)\|Ay - AX\beta\|^2 = (y - X\beta)^TA^TA(y - X\beta) = (y - X\beta)^TC^{-1}(y - X\beta)

Minimizing is OLS on the transformed data, so its normal equations (AX)T(AX)β^=(AX)TAy(AX)^T(AX)\hat\beta = (AX)^TAy become, using ATA=C−1A^TA = C^{-1} again:

XTC−1X β^=XTC−1y  ⇒  β^=(XTC−1X)−1XTC−1yX^TC^{-1}X\,\hat\beta = X^TC^{-1}y \;\Rightarrow\; \hat{\beta} = (X^TC^{-1}X)^{-1}X^TC^{-1}y

  • = MLE under Gaussian errors; scaling of CC cancels. Diagonal CC = weights wi=1/σi2w_i = 1/\sigma_i^2 (e.g. yiy_i an average of nin_i raw points ⟹ wi=niw_i = n_i, identical to OLS on the raw data). Mixed-effects/longitudinal structure ⟹ block-correlation CC.

Outliers in yy (quadratic loss amplifies them):

  • Huber/M-estimation: quadratic near 0, linear beyond kk; fit by iteratively reweighted LS; k=1.345k = 1.345 keeps 95% Gaussian efficiency. The fitted weights are themselves the diagnostic.
  • L1 / quantile (pinball) regression: LαL_\alpha population minimizer satisfies P(Y≥γ)=1−αP(Y \geq \gamma) = 1 - \alpha — α=12\alpha = \frac{1}{2} gives the conditional median, general α\alpha gives conditional quantiles.
  • LMS/LTS: trim the largest residuals; breakdown point ~50% (half the data can be arbitrarily corrupted); LTS preferred for efficiency.
  • Interpretation guard rails: a strong coefficient can be a proxy for an unmeasured cause (Tukey’s warning — undetectable from inside the model); with observational data, coefficients are comparisons across units, not effects of intervening. That gap is track 9’s whole subject.

Core competency set

  • Run the whitening derivation; map heteroscedasticity to weights and averages to wi=niw_i = n_i.
  • Match each robust method to its failure mode; derive the pinball-quantile connection.

10. Memorize cold

The instant-recall layer — these should cost nothing to produce; the fluency flashcards drill exactly this list.

  • The ladder: rung 0 projection / rung 1 Gauss–Markov / rung 2 t & F.
  • Cov⁡(β^)=σ2(XTX)−1\operatorname{Cov}(\hat{\beta}) = \sigma^2(X^TX)^{-1}; Var⁡(β^1)=σ2∑(xi−xˉ)2\operatorname{Var}(\hat{\beta}_1) = \frac{\sigma^2}{\sum(x_i - \bar{x})^2}; s2=RSSn−ps^2 = \frac{RSS}{n-p}; β^j−βjSE^∼tn−p\frac{\hat\beta_j - \beta_j}{\hat{SE}} \sim t_{n-p}.
  • Estimable ⟺ cc ∈ row space of XX; Gauss–Markov: cTβ^c^T\hat\beta is BLUE (projection proof).
  • Adjusted-variable theorem: β^j=⟨y,xj⋅rest⟩∥xj⋅rest∥2\hat{\beta}_j = \frac{\langle y, x_{j\cdot\text{rest}}\rangle}{\|x_{j\cdot\text{rest}}\|^2}, Var⁡=σ2∥xj⋅rest∥2\operatorname{Var} = \frac{\sigma^2}{\|x_{j\cdot\text{rest}}\|^2}; ρ=0.9\rho = 0.9 ⟹ ~5.3× variance inflation, coefficient seesaw.
  • HH: symmetric idempotent, tr⁡H=p\operatorname{tr} H = p; Cov⁡(y^)=σ2H\operatorname{Cov}(\hat y) = \sigma^2H, Cov⁡(r)=σ2(I−H)\operatorname{Cov}(r) = \sigma^2(I-H); Hii=1n(1+Mahalanobis2)H_{ii} = \frac{1}{n}(1 + \text{Mahalanobis}^2); flag Hii>2p/nH_{ii} > 2p/n.
  • LOO: yi−y^−i=ri1−Hiiy_i - \hat{y}_{-i} = \frac{r_i}{1 - H_{ii}}; Cook’s Di=1pri∗2Hii1−HiiD_i = \frac{1}{p}r_i^{*2}\frac{H_{ii}}{1-H_{ii}}; SMW for rank-one updates.
  • SVD regression: β^=VD−1UTy\hat{\beta} = VD^{-1}U^Ty; Var⁡(vjTβ^)=σ2/dj2\operatorname{Var}(v_j^T\hat\beta) = \sigma^2/d_j^2; PCA = SVD of centered XX.
  • Bias–variance: Err=σ2+bias2+var\mathrm{Err} = \sigma^2 + \text{bias}^2 + \text{var}; df=1σ2∑Cov⁡(y^i,yi)=tr⁡(Hλ)df = \frac{1}{\sigma^2}\sum\operatorname{Cov}(\hat y_i, y_i) = \operatorname{tr}(H_\lambda); Cp=RSSn+2 df σ2nC_p = \frac{RSS}{n} + \frac{2\,df\,\sigma^2}{n}.
  • Ridge: (XTX+λI)−1XTy(X^TX + \lambda I)^{-1}X^Ty; SVD factor djdj2+λ\frac{d_j}{d_j^2 + \lambda}; df=∑dj2dj2+λdf = \sum\frac{d_j^2}{d_j^2+\lambda}.
  • Lasso: ∣⟨xj,r⟩∣=λ|\langle x_j, r\rangle| = \lambda on the active set; λmax=max⁡∣⟨xj,y⟩∣\lambda_{max} = \max|\langle x_j, y\rangle|; soft threshold S(β^,λ)=sign(β^)(∣β^∣−λ)+S(\hat\beta, \lambda) = \mathrm{sign}(\hat\beta)(|\hat\beta| - \lambda)_+; df ≈ ∣A∣|A|.
  • WLS: whiten via C=BBTC = BB^T; β^=(XTC−1X)−1XTC−1y\hat\beta = (X^TC^{-1}X)^{-1}X^TC^{-1}y; Huber k=1.345k = 1.345; pinball α\alpha ⟹ conditional α\alpha-quantile.
  • Named moves: adjust-then-regress (MGS, stepwise, partial coefficients); split-orthogonally-and-drop (Gauss–Markov, min-norm); whiten-then-OLS; rank-one-update (SMW/LOO); add-and-subtract (bias–variance — fourth appearance); eavesdropping-df (covariance definition).