Chapter 1 — Linear Regression: Drawing the Best Line

By the end of this chapter you can fit a line by hand and in code, prove why it is the best line, test every assumption it rests on, correct for biased data, and answer any regression question an interviewer asks. One idea carries all of it: choose the line that makes the squared vertical misses as small as possible.

How to read it. The chapter is a chain of 13 questions in three acts. Each section opens with its question and a 30-second version, and ends with the question it leaves open — so you always know why you're reading what you're reading.

Two passes. On the first pass, read only the question, the 30-second version, and the figures of each section. You'll follow the whole arc in under an hour. On the second pass, read everything, doing the Remember box at the end of each section from memory before you look. Practice lives in its own tabs: Derivations (what to prove on a whiteboard) and Interview Bank (every question by level, with the intuition behind each answer).

1. The Problem and the Intuition

The question: five houses sold nearby. What is mine worth — and which line through the dots should we trust?

The 30-second version: draw the line whose vertical misses, squared and added up, are as small as possible. Squaring treats over- and under-shooting the same, punishes big misses, and keeps the math smooth enough to solve. First read: all of it.

The problem

You want to sell your 1,500 sq ft house. Five neighbors sold recently:

Bigger houses sell for more, but not perfectly: house C is the same size as yours and sold for less than the trend suggests. You need a rule that uses all five sales, not just the nearest one. The question — "how much does price rise per square foot, and what does that say about my house?" — is a regression question. We'll use these five houses for the entire chapter, so every formula can be checked by hand.

The intuition

Plot the houses as dots and lay a ruler across them. Each dot sits a little above or below the ruler. That vertical gap is the residual: actual price minus the line's price. A good ruler makes the residuals small.

We judge "small" by squaring each residual and adding them up. Three reasons, remembered as S-B-S:

• Sign: a miss of −20 is as bad as +20; squaring makes both positive so they can't cancel.
• Big misses hurt more: a 30k miss costs 900, three 10k misses cost only 300. Squaring prefers many small errors to one large one.
• Smooth: a squared function is a smooth bowl, so calculus can find its lowest point exactly.

The best line runs through the middle of the cloud. Notice the blue sticks above and below it cancel out exactly — that's not luck, and Section 3 shows why.

Remember: Regression = the line that makes the sum of squared vertical misses as small as possible. Everything else in this chapter is a consequence of that one sentence.

Next question → We have a rule for "best line" in words. Before computing it, we need to name every part precisely and say what we're assuming about the world. (Section 2)

2. The Model and the Loss

The question: what exactly is the model, and what is the loss measuring?

The 30-second version: y = Xw + b + \text{noise}. "Linear" means linear in the weights, so curves like x^2 are allowed. The squared loss makes the line predict the average y at each x; an absolute loss would predict the median. First read: the vocabulary table, the model, what "linear" means. Skim the loss-family table.

Vocabulary you'll be asked about

The same thing has three names depending on whether you learned it from statistics, ML, or economics. Interviewers switch between them freely.

Two traps in that table. First, ML calls b the "bias", which has nothing to do with statistical bias (Section 5). Second, errors and residuals are different things. Memory hook: Errors come from the truth and are never seen; Residuals are what Remains after our fit. Residuals are our only window onto the errors, which is why every diagnostic in Section 6 plots residuals.

The model

For one feature (simple linear regression):

y_i

the actual sale price of house i

x_i

the size of house i (the feature)

i

which house (row) we're talking about: 1, 2, …, n

\beta_0

the true intercept: the price the real-world line gives a 0 sq ft house. Unknown

\beta_1

the true slope: how much price really rises per extra sq ft. Unknown

\varepsilon_i

the error: how far house i's price lands from the true line (luck, haggling, a nice kitchen). Never observed

\hat{y}_i

our line's predicted price for house i (the hat means "estimated")

b

our intercept: the estimate of \beta_0

w

our weight (slope): the estimate of \beta_1, in $k per sq ft

For many features (multiple linear regression), with p features:

p

the number of features

x_{i1}

feature 1 of row i (first subscript = row, second = which feature)

x_{ip}

the last (p-th) feature of row i

w_1

the weight on feature 1: change in \hat y per unit of x_1, holding the other features fixed

w_p

the weight on the last feature

x_{i2}

feature 2 of row i

Multivariate regression is different again: it predicts several y's at once. Interviewers like this distinction.

What "linear" really means

Linear means linear in the weights, not in the features. All of these are linear regression, because each is a weighted sum where the weights are plain multipliers:

• \hat y = b + w_1 x + w_2 x^2 (a curve — still linear regression)
• \hat y = b + w_1 \log(x)
• \hat y = b + w_1 x_1 + w_2 x_1 x_2 (an interaction)

This is not: \hat y = b \cdot e^{wx}, because w sits inside the exponent. (Taking logs of both sides turns it into one that is.)

The test: "Can I compute every column from the data first, then just multiply each by a weight and add?" If yes, it's linear regression.

The loss, and its family

The loss turns "how wrong is this line?" into one number. Ours is the mean squared error:

\text{MSE}(w, b)

mean squared error: one number scoring how wrong the line with slope w and intercept b is

n

the number of rows (houses)

\sum_{i=1}^{n}

add up the term after it for every row, i = 1 to n

(y_i - \hat{y}_i)^2

house i's miss (residual), squared so misses can't cancel and big ones count extra

RSS (residual sum of squares, also called SSE) is the same sum without the 1/n. RMSE is its square root, which puts it back in dollars. All three have the same minimizer, so they give the same line.

Squared error isn't the only choice. Each loss quietly decides what kind of average the model predicts:

The two robust losses in the table, written out (e = y - \hat y, threshold \delta, quantile q):

\ell_\delta(e)

the Huber penalty for a miss of size e: squared when small, straight-line when large

e

the residual y - \hat y: how far the line misses a house's price, in $k

\delta

the threshold (in $k) where Huber switches from squaring the miss to counting it linearly

\ell_q(e)

the pinball penalty: charges under-predictions by q and over-predictions by 1 - q

q

the target quantile, e.g. 0.9 to predict the 90th-percentile price

Why squared gives the mean: the number c that minimizes \sum (y_i - c)^2 is exactly the average of the y_i (take the derivative, set it to zero). Regression just lets that average depend on x.

Remember: Square → mean, Absolute → median. Linear = linear in the weights. Residuals are what Remain; Errors are the unseen truth.

Next question → The loss is written down. Can we find the weights that minimize it with a formula? (Section 3)

3. Solving It Exactly: the Closed-Form Solution

The question: can we find the best weights with one formula?

The 30-second version: yes. The loss is a bowl; set its slope to zero and solve X^\top X w = X^\top y. Geometrically, the fit drops a perpendicular from the actual prices onto the space of predictions a line can make, so the residuals sit at right angles to every feature. First read: what "closed form" means, the by-hand example, the normal equation. Second read: the geometry and when the formula breaks.

What "closed form" means

A closed-form solution is an answer you can write as one finite formula — adds, multiplies, a square root, a matrix inverse — and evaluate in a single shot. No guessing, no repeating, no learning rate, no "stop when it's close enough".

The school example: the quadratic formula solves ax^2 + bx + c = 0 in closed form. The alternative is guess-and-check — try x = 1, too big, try 0.5, closer… That's an iterative (numerical) method, and it's what Section 4 does.

Why least squares gets one. The loss is quadratic in w → its gradient is linear in w → setting a linear thing to zero is a system of linear equations → linear equations can be solved exactly. Memory chain: quadratic loss → linear gradient → one solve.

The simple case, by hand

Setting both partial derivatives of the MSE to zero:

L

the loss: here the MSE, viewed as a function of w and b

\frac{\partial L}{\partial b}

"how fast the loss changes as you nudge b, holding w fixed": a slope, not a fraction to divide. Zero at the bottom of the bowl

\sum_i

add up over every house i

y_i - w x_i - b

house i's residual: actual price minus the line's prediction

\bar{y}

the average price over all n houses (bar = mean)

\bar{x}

the average size

w

the slope we're solving for

b

the intercept we're solving for

n

the number of houses

\frac{\partial L}{\partial w}

how fast the loss changes as you nudge the slope w, holding b fixed

x_i - \bar{x}

how far house i's size is from the average size

S_{xy}

the co-movement sum: \sum (x_i - \bar x)(y_i - \bar y). Positive when big houses cost more

S_{xx}

the spread of sizes: \sum (x_i - \bar x)^2

x_i

house i's size in sq ft

y_i - \bar{y}

how far house i's price is from the average price

With our five houses, \bar x = 1{,}520 sq ft and \bar y = $416k:

Plugging in:

\hat y(1{,}500)

the line's predicted price, in $k, for a 1,500 sq ft house like yours

Three ways to read the slope, all equal:

\text{Cov}(x, y)

the covariance of size and price: how they move together (positive = bigger costs more)

\text{Var}(x)

the variance of the sizes: how much they spread around their average

r

the correlation between size and price, from −1 to 1 and unit-free

s_y

the standard deviation of the prices, in $k

s_x

the standard deviation of the sizes, in sq ft

R^2

the share of the spread in prices that the line explains (Section 9)

In words: how x and y move together, over how x moves alone — or the correlation rescaled into y-units per x-unit. The slope always has the same sign as the correlation.

The general case: the normal equation

Stack the data into a design matrix X: one row per house, one column per feature, plus a first column of 1s so the intercept becomes just another weight.

L(\mathbf{w})

the loss (sum of squared misses) as a function of the whole weight vector

\mathbf{w}

all the weights stacked in one column: intercept first, then slope (bold = vector)

\mathbf{y}

all n sale prices stacked into one column vector

X

the design matrix: one row per house, a column of 1s, then the size column

\lVert \mathbf{y} - X\mathbf{w} \rVert^2

the squared length of the residual vector, i.e. the sum of squared misses

X^\top

X transposed: rows become columns, so each row is now one feature

X^\top X

a small table of how the feature columns overlap with each other (one row and column per weight)

\nabla_{\mathbf{w}} L

the gradient: the vector of slopes of the loss, one per weight; points uphill

\hat{\mathbf{w}}

the best-fit weights the formula returns (hat = estimated from data)

(X^\top X)^{-1}

the matrix inverse of X^\top X: undoes it, the matrix version of dividing by it

The middle equation, X^\top X w = X^\top y, is the set of normal equations. Because the loss is a bowl (X^\top X is positive semi-definite, so the loss is convex), the point where the gradient is zero is the global minimum — there are no other valleys to get stuck in.

The geometry: why it's called "normal"

Think of the n prices as one point y in n-dimensional space. Every possible line produces a prediction vector Xw, and all of those live on a flat plane (the column space of X). Least squares picks the point on that plane closest to y — and the closest point is always straight below, at a right angle.

The residual is perpendicular (normal) to the plane, so X^\top e = 0. The matrix that does the dropping is the hat matrix, H = X(X^\top X)^{-1}X^\top, because \hat y = Hy — it "puts the hat on y". Its diagonal entries are each point's leverage (Section 6).

H

the hat matrix: turns the actual prices into the line's predicted prices

\hat{\mathbf y}

all n predicted prices stacked in a vector

\text{tr}(H)

the trace of H: the sum of its diagonal entries (the leverages)

p

the number of weights fitted, here 2 (intercept and slope)

X^\top e = 0 hands you three free facts about every OLS fit with an intercept:

1. The residuals sum to zero. Ours: 1.8 + 10.4 - 31.9 + 25.9 - 6.3 = 0 (up to rounding).
2. The residuals are uncorrelated with every feature in the training data.
3. The line passes through the point of means (\bar x, \bar y) = (1{,}520,\ 416).

Fact 2 hides an interview trap: you cannot check whether a feature is correlated with the error by correlating it with the residuals. OLS forces that correlation to zero by construction, whether the assumption holds or not. (Section 6, assumption R.)

When the closed form breaks

The formula needs (X^\top X)^{-1} to exist. It doesn't when:

• Perfect multicollinearity: one column is an exact combination of others — size in sq ft and size in sq m, or one-hot encoding every category plus an intercept (the dummy variable trap).
• More features than rows (p > n): infinitely many lines fit perfectly.

Fixes: drop the redundant column; use the pseudo-inverse, which picks the smallest-norm solution among the infinitely many; or use ridge, since X^\top X + \lambda I is always invertible for \lambda > 0.

In practice, never compute the inverse. Explicit inversion is slow and numerically fragile — forming X^\top X squares the condition number, so small rounding errors grow. Libraries solve with a QR decomposition or SVD instead (np.linalg.lstsq). Cost: about n \cdot p^2 to form or factor, plus p^3 to solve. That p^3 is why the closed form stops being practical past tens of thousands of features.

Remember: Quadratic loss → linear gradient → one solve. The answer is a projection: X^\top e = 0, so residuals sum to zero and the line passes through the means. Never invert — use lstsq.

Next question → The formula costs about p^3 and doesn't exist for most models. What do we do then? (Section 4)

4. Solving It Iteratively: Gradient Descent

The question: how do we find the weights when the formula is too expensive or doesn't exist?

The 30-second version: start anywhere, feel the slope, step downhill: w \leftarrow w - \eta \nabla L. Scale the features so the bowl is round. On this bowl it lands exactly where the formula does — and it's the same procedure that trains every neural network. First read: the idea, Figure 1.3, the learning-rate table. Second read: scaling and the three flavors.

Why iterate when a formula exists?

Three reasons, and the third is the one that matters for the rest of the book:

1. Too many features. The closed form costs about p^3; at a million features that's 10^{18} operations.
2. Too much data. The closed form needs all of X at once. Gradient descent can stream it in small batches.
3. No formula exists. Logistic regression, Lasso, and every neural network have no closed form. Gradient descent is the engine that trains all of them, including LLMs. Learning it on a problem where we know the right answer is how we learn to trust it.

The idea

You're on a foggy hillside and want the valley floor. You can't see it, but you can feel which way the ground slopes under your feet. Take a step downhill. Feel again. Repeat. The gradient is the direction of steepest uphill, so we step the opposite way:

\nabla_{\mathbf{w}} L

the gradient: the vector of slopes of the loss, one per weight; points uphill

n

the number of houses (rows)

X

the design matrix: one row per house, a column of 1s, then the size column

X^\top

X transposed: rows become columns, so each feature column meets every row's miss

\mathbf{y} - X\mathbf{w}

how wrong the current line is on each house: the residual vector

\mathbf{y}

the vector of all n actual prices

\mathbf{w}

the current weights (intercept and slope) as one vector (bold = vector)

\eta

the learning rate: how big a step to take, a small number you choose

\leftarrow

"replace with": compute the right side, then overwrite \mathbf w with it

Read the gradient in words: each feature column, weighted by how wrong we are on each row. Features that line up with big errors get big corrections.

The step size is \eta \times \text{slope}, so the steps shrink on their own as the ground flattens — you don't need to slow down by hand. Because the least-squares loss is a single bowl (convex), gradient descent always reaches the same bottom as the closed form.

The learning rate η

For least squares there's an exact limit, set by the steepest direction of the bowl — steeper bowl, smaller safe step:

\lambda_{\max}

the largest eigenvalue: how sharply the bowl curves in its steepest direction

\tfrac{2}{n}X^\top X

the bowl's curvature matrix (the Hessian of the MSE)

Why feature scaling matters

If one feature is in thousands (sq ft) and another in single digits (bedrooms), the bowl is a long, narrow ravine. A step size small enough for the steep direction is far too small for the shallow one, so gradient descent zigzags down the ravine walls. Standardizing each feature (subtract its mean, divide by its standard deviation) makes the bowl round, and the steps point straight at the bottom. The ratio of steepest to shallowest curvature is the condition number; scaling shrinks it.

Scaling doesn't change the closed-form fit's predictions, but it changes the weights' units — remember to convert them back.

Three flavors

Mini-batch is what everyone actually uses, including for LLMs. Two words interviewers check: an iteration is one step (one batch); an epoch is one full pass through the data.

When to stop: the gradient's size is near zero, the loss stops improving by more than a tiny tolerance, or — the ML answer — the validation loss stops improving (early stopping).

Preview of Part IV: momentum keeps a running average of past steps so you roll through ravines instead of zigzagging; Adam also adapts the step size per weight. Both are this same update with memory added.

Remember: \text{Step} = -\eta \times \text{slope}. Scale features first. Batch / mini-batch / stochastic = all / some / one row per step. Convex bowl → same answer as the formula.

Next question → We can compute the line two ways. But it came from just five sales. Would different sales give a different line — and is ours aimed at the truth? (Section 5)

5. Is Our Line Any Good? Estimators, Bias, and Unbiasedness

The question: is our estimate aimed at the truth, and how much would it wobble with different data?

The 30-second version: the fitted slope is one random draw. OLS is unbiased if the model is linear, X has full rank and the errors are exogenous. Its wobble shrinks with less noise, more data, and more spread in x. Estimate the noise by dividing by n - p, not n. Under the Gauss–Markov conditions it's the best linear unbiased estimator. First read: the big shift, Figure 1.4, the unbiasedness proof and what it used. Second read: n - p and Gauss–Markov.

The big shift: the line itself is random

If five different neighbors had sold, you'd get a different line. So 0.2074 isn't "the" slope — it's one draw from a distribution of slopes you could have gotten. The true slope \beta_1 is a fixed, unknown number; our \hat w is a random variable that tries to land near it.

Bias and variance of an estimator

\text{Bias}(\hat\theta)

how far the estimator's average lands from the truth: its aim error

\hat\theta

any estimator, e.g. our fitted slope \hat w (hat = estimated from data)

\theta

the true value it tries to hit, e.g. the real slope \beta_1

\mathbb{E}[\hat\theta]

the estimator's average over many repeated samples of house sales

\mathbb{E}

the average over repeated samples (a fresh set of house sales each time)

\text{Var}(\hat\theta)

how much the estimate wobbles from sample to sample

An estimator is unbiased when its bias is zero: averaged over infinitely many repeated samples, it lands exactly on the truth. It does not mean any single estimate is right — just that it isn't systematically too high or too low.

The two combine into one number, the estimator's mean squared error:

\text{MSE}(\hat\theta)

the estimator's mean squared error: its average squared distance from the truth

This is the bias–variance decomposition. It's why a little bias can be worth it: if accepting some bias cuts the variance a lot, total error falls. That trade is the entire idea behind regularization (Section 10).

Three more properties interviewers name-drop:

• Consistent: the estimate converges to the truth as n \to \infty. Different from unbiased: an estimator can be biased but consistent (dividing by n instead of n - 1), or unbiased but inconsistent (using only the first data point).
• Efficient: lowest variance among unbiased estimators.
• Asymptotically normal: for large n, the sampling distribution looks Gaussian — which is what makes confidence intervals work.

Proof that OLS is unbiased

Substitute the true model y = X\beta + \varepsilon into the OLS formula:

\hat{\mathbf{w}}

the OLS weights: intercept and slope fitted from the data (hat = estimated)

X

the design matrix: a column of 1s plus the size column, one row per house

X^\top

X transposed: rows become columns

(X^\top X)^{-1}

the matrix inverse of X^\top X: undoes it

\boldsymbol\beta

the true weights (intercept and slope) we're trying to recover (bold = vector)

\boldsymbol\varepsilon

the true errors of all the houses stacked in a vector: luck, haggling, kitchens

\mathbb{E}[\hat{\mathbf{w}} \mid X]

the fitted weights averaged over repeated samples, given (holding fixed) the house sizes in X

\mathbb{E}[\boldsymbol\varepsilon \mid X]

the average error at each size; zero means nothing hidden tilts prices (exogeneity)

Look at what the proof used: the model is linear (y = X\beta + \varepsilon), X^\top X is invertible (no perfect multicollinearity), and the errors average zero at every x (exogeneity). Look at what it didn't use: constant variance, independence, normality. Those affect how precise and how testable the estimate is — never whether it's centered on the truth. This is the single most useful fact for answering assumption questions.

How precise is it? The variance of OLS

If the errors also have constant variance \sigma^2 and are uncorrelated:

\text{Var}(\hat{\mathbf{w}} \mid X)

how much the fitted weights wobble across samples, given the house sizes (a matrix)

\sigma^2

the noise variance: the typical squared size of a house's error, in $k²

\text{Var}(\hat w)

how much the fitted slope alone wobbles from sample to sample

S_{xx}

the spread of sizes: \sum (x_i - \bar x)^2

n

the number of houses

\text{Var}(x)

the variance of the house sizes: how spread out they are

Read it as a recipe for a precise slope: less noise (\sigma^2 small), more data (n big), more spread in x (\text{Var}(x) big). That last one is a design lesson — to learn how price depends on size, sample houses of very different sizes. It also shows why multicollinearity hurts: near-duplicate columns make X^\top X nearly singular, so its inverse — and the variance — blows up.

Estimating the noise: why divide by n − p

We don't know \sigma^2, so we estimate it from the residuals. The unbiased version is:

s^2

our estimate of the noise variance \sigma^2, built from the residuals

p

the number of fitted parameters: 2 here (w and b)

\sum_{i=1}^{n}

add up over every house, i = 1 to n

e_i^2

house i's residual, squared

Why not divide by n? The fitted line chases the data, so residuals are systematically smaller than the true errors. Each of the p estimated parameters imposes one constraint on the residuals (that's X^\top e = 0 again), leaving only n - p free — the degrees of freedom. One can show E[\text{RSS}] = (n - p)\sigma^2, so dividing by n - p makes s^2 unbiased. It's the same reason the sample variance divides by n - 1 (Bessel's correction): there, one parameter (the mean) was estimated.

Our houses: \text{RSS} = 1{,}837, n - p = 5 - 2 = 3, so s^2 \approx 612 and s \approx $24.7k — a typical miss. Dividing by n would give 367 and understate the noise by 40%.

Dividing by n is what maximum likelihood gives — so maximum likelihood estimates are not always unbiased. A favorite interview question.

Gauss–Markov: OLS is BLUE

Gauss–Markov theorem. If the model is linear, X has full rank, errors have mean zero given X, constant variance, and are uncorrelated, then OLS has the smallest variance of any linear unbiased estimator. It is BLUE: the Best Linear Unbiased Estimator.

Three things to say when asked about it:

1. Normality is not required. Add it, and OLS becomes the best of all unbiased estimators, linear or not, and exact t- and F-tests become valid.
2. "Best" only among unbiased. A biased estimator such as ridge can beat OLS on total error. Gauss–Markov is a reason to use OLS when you need unbiasedness, not proof that it's the best predictor.
3. Break constant variance and OLS stays unbiased but stops being best. Weighted or generalized least squares (WLS/GLS) becomes BLUE, and the usual standard errors become wrong.

Remember: Bias = where you aim, variance = how steady your hand. OLS is unbiased if the model is linear, X has full rank, and errors are exogenous. Add independent, equal-variance errors and it's BLUE. Add normality for exact tests. Divide RSS by n - p, not n.

Next question → Every guarantee in this section rested on assumptions. Which ones exactly, how do we check them, and what do we do when they fail? (Section 6)

6. The Assumptions: Remember, See, Test, Fix

The question: which assumptions do the guarantees rest on, and how do we check each one?

The 30-second version: the model's name is the checklist, L-I-N-E-A-R. L, A and R keep the estimate on target; I and E keep the error bars honest; N makes small-sample tests exact. For each: look at a residual plot first, run the formal test second, apply the fix third. First read: the LINEAR table, the ladder, Figure 1.5. Second read: the test table and influence.

The mnemonic: the model's own name, L-I-N-E-A-R

Now attach each letter to what it buys you. This ladder is the answer to almost every "what happens if…" question:

See it: every failure has a fingerprint

You always look before you test. Plot the residuals — they're your window onto the unseen errors — and each violation draws its own shape:

The four plots statsmodels and R draw by default map straight onto the letters: Residuals vs Fitted (L), Normal Q-Q (N), Scale-Location, which plots \sqrt{\lvert \text{standardized residual} \rvert} against fitted values (E), and Residuals vs Leverage (influential points). For I, plot residuals in time order or their autocorrelation (ACF).

Test it: the formal checks

Every test below has the null hypothesis "the assumption holds". So: if p is low, the assumption must go.

The formulas behind the main checks:

\text{VIF}_j

the variance inflation factor: how many times feature j's slope variance is blown up by overlap

R_j^2

the R^2 from regressing feature j on all the other features

DW

the Durbin–Watson statistic: near 2 means no lag-1 autocorrelation

\sum_t

add up over the rows in time order, t = 1, 2, \ldots

e_t

the residual of the row at time t

e_{t-1}

the residual of the row just before it

\hat\rho_1

the estimated correlation between each residual and the previous one (lag 1)

e_i^2

house i's squared residual, a stand-in for its own noise variance

\mathbf x_i

house i's features as a vector (bold = vector)

LM

the Breusch–Pagan test statistic (a Lagrange multiplier statistic)

n

the number of rows (houses)

R^2_{\text{aux}}

the R^2 of that auxiliary regression of e_i^2 on the features

\sim

"is distributed as" (when the null hypothesis holds)

\chi^2_{p-1}

a chi-squared distribution with p - 1 degrees of freedom

h_{ii}

house i's leverage: the i-th diagonal entry of the hat matrix

\mathbf x_i^\top

house i's feature vector, transposed (rows become columns) into a row

X

the design matrix: one row per house, a column of 1s, then the features

(X^\top X)^{-1}

the matrix inverse of X^\top X: undoes it

p

the number of fitted parameters, intercept included

D_i

Cook's distance for house i: how far the whole line moves if you delete it

s^2

the estimated noise variance, \text{RSS} / (n - p)

Two cautions about the tests. With a huge n they flag tiny, harmless deviations; with a small n they miss real ones. Plots first, tests second.

Outliers, leverage, and influence

Not an assumption, but asked in every interview. Memory hook: outlier = weird y; leverage = weird x; influence = weird and it matters.

A point needs both leverage and a large residual to be influential: a weird x that sits right on the line pulls nothing. What to do: never delete automatically. Check for a data-entry error first, then report results with and without it, or switch to a loss that shrugs off outliers (Huber, or RANSAC, which fits on random subsets and keeps the consensus).

Run every check in one go

import statsmodels.api as sm
from statsmodels.stats.diagnostic import het_breuschpagan, linear_reset
from statsmodels.stats.stattools import durbin_watson
from statsmodels.stats.outliers_influence import variance_inflation_factor, OLSInfluence
from scipy.stats import shapiro

m = sm.OLS(y, sm.add_constant(X)).fit()
r = m.resid
print("L  RESET p      :", linear_reset(m, power=2, use_f=True).pvalue)
print("I  Durbin-Watson:", durbin_watson(r))
print("N  Shapiro p    :", shapiro(r).pvalue)
print("E  Breusch-Pagan p:", het_breuschpagan(r, m.model.exog)[1])
print("A  VIFs         :", [variance_inflation_factor(m.model.exog, i) for i in range(1, m.model.exog.shape[1])])
print("   max Cook's D :", OLSInfluence(m).cooks_distance[0].max())
# R: no test from residuals — reason about confounders

Remember: The model's name is the checklist: L-I-N-E-A-R. LAR to aim right, IE to know how sure, N for exact tests. Plot residuals first; every failure has a shape. If p is low, the assumption must go. Outlier = weird y, leverage = weird x, influence = both.

Next question → Assumptions are about the model. But what if the rows themselves are unrepresentative, incomplete, or mismeasured? (Section 7)

7. Bias in the Data, and How to Correct It

The question: what if the data itself is biased — the wrong rows, a missing variable, bad measurements?

The 30-second version: ask Who, What, How, When. A missing confounder leaks its effect into the features you kept (bias = \beta_2 \delta). Noisy x flattens slopes; noisy y doesn't. Selecting rows on y biases; selecting on x doesn't. Missing data and drift each have a named fix. First read: the four questions, omitted variables, measurement error, selection. Second read: missing data, drift, fairness, the toolbox.

Section 5 was about bias in the estimator: does the recipe aim at the truth when the data are fine? This section is about bias in the data: the rows you have don't represent the world you'll predict for, or they were measured wrong. No estimator can fix bad data by itself — OLS on a biased sample is a perfectly unbiased estimate of the wrong thing.

The mnemonic: four questions before you trust a dataset

Who, What, How, When.

1. WHO is in the data? Selection, survivorship, non-response, coverage, imbalance.
2. WHAT is missing? Omitted variables (confounders) and missing values.
3. HOW was it measured? Measurement error, proxies, biased labels.
4. WHEN was it collected? Drift over time, and leakage from the future.

WHAT is missing: omitted-variable bias

The most-asked one. Suppose the truth is y = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + \varepsilon, but you leave x_2 out. Then:

\mathbb{E}[\hat w_1]

the fitted slope on x_1, averaged over many repeated samples

\beta_1

the true effect of x_1 (say, size) on price

\beta_2

the true effect of the left-out feature x_2 on price

\delta

how much x_2 rises per unit of x_1: the slope from regressing x_2 on x_1

x_2

the omitted feature, e.g. how good the neighborhood is

x_1

the feature you kept, e.g. size

The left-out variable's effect leaks into whichever included feature it's correlated with. Bias needs both ingredients: x_2 must affect y (\beta_2 \ne 0) and be correlated with x_1 (\delta \ne 0). That's exactly what a confounder is.

Direction of the bias = sign of \beta_2 \times sign of the correlation:

The extreme version is Simpson's paradox: the slope flips sign once you control for a group variable. Example: across all cities, bigger homes look cheaper per sq ft, because big homes are mostly in cheap cities; within each city, the slope is positive.

And the flip side: adding a control variable is not always good. Don't control for something caused by x (a mediator, such as "number of bathrooms" if you want size's total effect), or for a common effect of x and y (a collider) — both create bias rather than remove it.

HOW was it measured: measurement error

Why noisy x flattens the slope: if we observe x = x^* + u (true value plus noise), then

\text{plim}\;\hat w

what the fitted slope settles to with unlimited data (its probability limit)

\beta

the true slope: the real $k per sq ft

\text{Var}(x^*)

the spread of the true sizes x^* (star = the true, unseen value)

\text{Var}(u)

the spread of the measurement noise u added to each size

\lambda

the reliability: the share of the measured spread in sizes that is real signal

The noise spreads the points sideways, and a line through a sideways-smeared cloud is flatter. Fixes: measure better or average repeated measurements; divide by a known reliability \lambda; errors-in-variables models (Deming or total least squares, which measure distance perpendicular to the line rather than vertically); or an instrumental variable.

WHO is in the data: selection

The surprising rule: selecting rows on x is safe; selecting on y is not.

• Keep only houses between 1,000 and 2,000 sq ft (selection on x): the slope inside that range is still unbiased. You just can't extrapolate.
• Keep only houses that sold above $300k (selection on y): low-priced small houses vanish, the left end of the cloud gets lifted, and the slope flattens. This is truncation.
• Keep only houses that sold at all: over-priced listings that never sold are missing, so the data look like the market is more efficient than it is. Survivorship bias is selection on the outcome by another name.

Missing values

Memory hook: Completely random → drop; At random → impute; Not at random → model it.

Mean-imputation is the common wrong answer: it shrinks variance and biases correlations toward zero.

WHEN was it collected: drift and leakage

• Drift. A model fit on 2019 sales underprices 2026 sales. Covariate shift = the x's change (bigger homes being built); concept drift = the relationship itself changes ($ per sq ft rises). Correct with time-based validation (train on the past, test on the most recent period), time features, regular retraining, and monitoring residuals over time.
• Leakage. A feature contains information you won't have at prediction time — "days on market" or "final appraisal" when predicting list price. Symptom: suspiciously great test scores. Correct with a feature audit ("would I know this at the moment of prediction?") and splits that respect time and groups.

Biased labels and fairness

If historical prices or appraisals were systematically lower in certain neighborhoods for reasons unrelated to the house, the model learns that pattern and repeats it. Features like zip code act as proxies even when sensitive attributes are excluded. Correct by auditing residuals and errors by group, questioning proxy features, and, where needed, fairness constraints. Part VIII returns to this for LLMs.

The correction toolbox, ranked

1. Fix the data: collect the missing variable, measure better, sample the missing group. Always the best fix.
2. Design around it: randomized experiments, instrumental variables, difference-in-differences, regression discontinuity, fixed effects.
3. Control for it: add the confounders as features (not mediators or colliders).
4. Reweight: post-stratification, inverse probability weights, importance weights.
5. Model the bias: Heckman, errors-in-variables, explicit missingness models.
6. Bound it and say so: sensitivity analysis — "how strong would a hidden confounder need to be to erase this effect?"

Remember: Ask Who, What, How, When before trusting any dataset. Omitted-variable bias = \beta_2 \times \delta, sign = sign × sign. Noisy x flattens; noisy y blurs. Select on x, safe; select on y, biased. Completely random → drop, At random → impute, Not at random → model.

Next question → Once we trust the estimate, how do we say how sure we are about it — and about a single prediction? (Section 8)

8. Inference: How Sure Are We?

The question: how sure are we about the slope, and about one prediction?

The 30-second version: the standard error measures the wobble from Section 5. t = \text{estimate} \div \text{SE}. The p-value says how surprising the data would be if the true effect were zero. A confidence interval covers the average price at a size; a prediction interval covers one new house and is always wider, because it adds that house's own noise. First read: the t-test and CI vs PI. Second read: the F-test and the summary walkthrough.

Prediction asks "what's the price?" Inference asks "how sure are we about the slope, and could it really be zero?" Inference is where assumptions I, E, and N earn their keep.

Standard error and the t-test

The slope's standard error estimates how much \hat w would wobble across samples (Section 5), with s in place of the unknown \sigma:

\text{SE}(\hat w)

the standard error of the slope: how much \hat w would wobble across samples

s

the estimated noise standard deviation, $24.75k, standing in for the unknown \sigma

S_{xx}

the spread of sizes, \sum (x_i - \bar x)^2 = 788{,}000

To test "size has no effect" (H_0: \beta_1 = 0), divide the estimate by its wobble:

t

the t-statistic: how many standard errors the slope sits away from zero

\hat w

our fitted slope, 0.2074 (hat = estimated)

n

the number of houses, 5

p

the number of fitted parameters, 2 (w and b)

A t of 7.44 means the slope is 7.44 standard errors away from zero. With 3 degrees of freedom, p \approx 0.005.

What a p-value is: the probability of seeing a t at least this extreme if H_0 were true. Interviewers listen for what it is not:

• Not the probability that H_0 is true.
• Not the probability the result is due to chance.
• Not a measure of effect size — with a million rows, a $0.01-per-sq-ft effect is "significant" and useless. Statistical significance \ne practical significance.
• Test 20 useless features at 0.05 and expect one to look significant. Correct for multiple testing (Bonferroni: divide 0.05 by the number of tests; or control the false discovery rate with Benjamini–Hochberg).

Confidence interval for a coefficient

t_{0.975,\,\text{df}}

the t-table cutoff leaving 2.5% in each tail with df degrees of freedom: 3.182 here

Correct reading: if we repeated this whole procedure on many samples, 95% of the intervals built this way would contain the true slope. The interval excludes 0, which matches p < 0.05 — the two always agree.

Confidence interval vs prediction interval

The most common inference question in interviews. For your 1,500 sq ft house:

\hat y_0

the line's predicted price for the new house ($412k for yours)

x_0

the new house's size, here 1,500 sq ft

t

the t-table cutoff, 3.182 for 95% with 3 degrees of freedom

\bar x

the average size of the five houses, 1,520 sq ft (bar = mean)

\mathbf{1}

the extra 1: your house's own noise, on top of the uncertainty about the line

Memory hook: the prediction interval adds the 1 — one extra \sigma^2 for your house's own luck, on top of uncertainty about the line. More data shrinks the CI toward zero width, but the PI never shrinks below the noise. Both widen as x_0 moves away from \bar x, which is the math behind "don't extrapolate".

The F-test

Tests whether all slopes are zero at once (H_0: the model is no better than predicting \bar y):

F

the F-statistic: variance explained per slope, divided by leftover variance per degree of freedom

Ours: F = 55.3, p = 0.005. With one feature, F = t^2 (7.44^2 \approx 55.3). With many features, a significant F but no significant individual t's is a classic sign of multicollinearity: together they matter, but the model can't tell which one.

Reading a regression summary, line by line

This is our five-house model in statsmodels:

                            OLS Regression Results
==============================================================================
Dep. Variable:                      y   R-squared:                       0.949
Model:                            OLS   Adj. R-squared:                  0.931
Method:                 Least Squares   F-statistic:                     55.32
No. Observations:                   5   Prob (F-statistic):            0.00503
Df Residuals:                       3   Log-Likelihood:                -21.861
Df Model:                           1   AIC:                             47.72
Covariance Type:            nonrobust   BIC:                             46.94
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const        100.8122     43.797      2.302      0.105     -38.568     240.192
x1             0.2074      0.028      7.438      0.005       0.119       0.296
==============================================================================
Omnibus:                          nan   Durbin-Watson:                   3.392
Prob(Omnibus):                    nan   Jarque-Bera (JB):                0.260
Skew:                          -0.397   Prob(JB):                        0.878
Kurtosis:                       2.215   Cond. No.                     6.22e+03
==============================================================================

Remember: t = estimate ÷ its wobble. p = how surprising the data would be if nothing were going on — not the chance H_0 is true. The prediction interval adds the 1. With one feature, F = t^2. Lower AIC/BIC wins.

Next question → Being sure about a coefficient isn't the same as predicting well. How good is the model overall, and compared with what? (Section 9)

9. Evaluating the Model

The question: how good is the model, and compared with what?

The 30-second version: RMSE punishes big misses, MAE reports the typical miss, R^2 compares you with just predicting the average. Numbers computed on the training data flatter you; score on held-out data and against a simple baseline. First read: all of it.

Every metric, on our five houses

The same metrics as formulas:

\frac1n

divide by n, the number of houses (5): turns the sum into an average

\sum_i

add up over every house i

e_i

house i's residual: actual price minus predicted, in $k

n

the number of houses, 5

y_i

house i's actual sale price

R^2

the share of the spread in prices the line explains, compared with just predicting the average

\bar y

the average price, $416k (bar = mean)

R^2_{\text{adj}}

adjusted R^2: charges each extra feature rent, so useless ones lower it

p

the number of fitted parameters, 2 here

\hat y_i

the line's predicted price for house i (hat = estimated)

\text{RMSE} \ge \text{MAE} always; the bigger the gap, the more your errors include a few large misses.

What R² really says, and what it hides

R^2 = the fraction of the variance in y that the model explains, compared with just predicting the mean. \text{TSS} = 35{,}720 is the total spread of prices around $416k; the line leaves \text{RSS} = 1{,}837 unexplained, so R^2 = 1 - 1{,}837 / 35{,}720 = 0.949.

Six traps, all asked:

1. R^2 never goes down when you add a feature — even pure noise. That's why adjusted R^2 exists: it charges each feature rent and falls if the feature doesn't pay it.
2. R^2 can be negative on test data, meaning the model is worse than predicting the training mean.
3. High R^2 doesn't mean the model is right. Anscombe's quartet: four datasets with the same line and R^2 \approx 0.67 — one linear, one curved, one with an outlier, one with a single leverage point. Always plot.
4. Low R^2 isn't always bad. In noisy fields (behavior, finance), R^2 = 0.1 can be valuable if the coefficient is precise and meaningful.
5. R^2 depends on the spread of x in your sample. Restrict to houses between 1,400 and 1,600 sq ft and R^2 collapses, though the slope is unchanged.
6. R^2 says nothing about causation or about whether the assumptions hold.

In-sample vs out-of-sample

Every number above was computed on the same five houses we fit on, so it's optimistic. The honest measure is error on data the model never saw: a held-out test set, or k-fold cross-validation (split into k parts, train on k - 1, test on the last, rotate). Chapter 5 covers this properly. For time-ordered data, always test on the most recent period.

And always compare against a baseline: predicting the mean (R^2 = 0) or last year's price. A model is only as good as the margin by which it beats the dumbest reasonable alternative.

Remember: RMSE punishes, MAE forgives, MAPE breaks at zero. R^2 compares you to the mean; adjusted R^2 charges rent for each feature. In-sample numbers flatter — score on held-out data, against a baseline, and plot.

Next question → What if a straight line on the raw features isn't good enough? (Section 10)

10. Extensions: Bending the Straight Line

The question: how do we bend the straight line without leaving linear regression?

The 30-second version: transform features (logs change what a coefficient means: percent instead of units), encode categories as k - 1 dummies, add interactions and splines. When features are many or collinear, regularize: L2 shrinks every weight, L1 sets some to zero. First read: transforms, categoricals, regularization with Figure 1.7. Second read: interactions, splines, relatives.

Transforms, and how to read the coefficient

Logs fix curved relationships, funnel-shaped residuals, and skewed targets in one move. But they change what w means:

y

the sale price ($k); logged on the left in the log– rows

x

the house size (sq ft); logged on the right in the –log rows

b

the intercept

w

the slope; what it means depends on which sides are logged

\Delta y

the change in price

\Delta x

the change in size

\%\Delta y

the percent change in price

\%\Delta x

the percent change in size

e^{w}

e \approx 2.718 raised to the power w; e^{w} - 1 is the exact fractional change per unit of x

Memory hook: log on the left → percent in y; log on the right → percent in x; both → percent for percent.

One trap: if you fit \ln y and predict e^{\text{predicted } \ln y}, you predict the median, which underestimates the mean. Multiply by e^{s^2/2} (if errors are normal) or use Duan's smearing estimate. Box–Cox generalizes the choice of transform by searching over powers of y.

To compare the importance of features in different units, standardize them first; the coefficients then mean "change in y per one standard deviation of x".

Categorical features

A zip code isn't a number. One-hot (dummy) encode it: k categories become k - 1 columns of 0s and 1s, and the left-out category is the baseline absorbed by the intercept. Each coefficient is then "the price difference versus the baseline, all else equal."

• Include all k dummies and an intercept and the columns sum to the intercept column — perfect multicollinearity, the dummy variable trap. Drop one (or drop the intercept).
• Don't encode categories as 1, 2, 3 unless they really are evenly spaced; the model will assume they are.
• Hundreds of categories: regularize, group rare ones, or use target encoding computed out-of-fold (otherwise it leaks the answer).

Interactions

If a pool adds more value to a big house than to a small one, the slope of size depends on the pool:

\hat y

the predicted price, in $k

w_1

the extra price per sq ft for a house without a pool

w_2

the pool's price premium (read at size 0, which is why we center size)

w_3

the extra price per sq ft that a pool adds: the interaction

The size slope is w_1 without a pool and w_1 + w_3 with one. Keep the main effects whenever you keep the interaction (the hierarchy principle), and center variables before multiplying them — otherwise x and x \cdot z, or x and x^2, are highly correlated and VIFs explode.

Polynomials and splines

Adding x^2, x^3… bends the line, but high degrees wiggle wildly between points and explode at the edges. Splines (piecewise polynomials joined smoothly at knots) and GAMs (a smooth curve per feature, added up) bend gently and stay linear in their weights.

Regularization: buying lower variance with a little bias

Section 5 showed \text{MSE} = \text{bias}^2 + \text{variance}, and that unbiased isn't always best. Regularization adds a penalty on the size of the weights:

\min_{\mathbf w}

find the weight vector \mathbf w that makes what follows as small as possible

\mathbf w

all the weights as one vector (bold = vector)

\mathbf y

all the sale prices as one vector

X

the design matrix: one row per house, one column per feature

\lVert \mathbf y - X\mathbf w \rVert^2

the sum of squared misses (RSS), written as a squared vector length

\lambda

the penalty strength: 0 gives plain OLS, bigger shrinks harder; tuned by cross-validation

\lVert \mathbf w \rVert_2^2

the squared L2 length of the weights: the sum of their squares

\hat{\mathbf w}

the ridge weights (hat = estimated from data)

X^\top

X transposed: rows become columns

(X^\top X + \lambda I)^{-1}

the matrix inverse of X^\top X + \lambda I (I = identity matrix); always exists when \lambda > 0

\lVert \mathbf w \rVert_1

the L1 length of the weights: the sum of their absolute values

The same problem can be written as "minimize the error, subject to the weights staying inside a region": a circle for ridge, a diamond for lasso. The solution is where the smallest error contour first touches the region.

Practical rules: standardize features first (the penalty is unfair to features in small units); don't penalize the intercept; choose \lambda by cross-validation. \lambda = 0 gives OLS; \lambda \to \infty pushes every weight to 0. In neural networks, ridge is called weight decay, and it's a standard setting when training LLMs.

Memory hook: L1 kills, L2 shrinks. Diamond corners sit on the axes; circles have no corners.

Other relatives, one line each

• Weighted least squares (WLS): each row weighted by 1 / its variance; the fix for a funnel.
• Generalized least squares (GLS): handles correlated errors too.
• Quantile regression: predicts the 10th or 90th percentile instead of the mean; great for price ranges.
• Robust regression: Huber loss, RANSAC, Theil–Sen — resistant to outliers.
• Bayesian linear regression: gives a full distribution over weights instead of one estimate; ridge is its special case.
• Principal component regression / partial least squares: regress on a few combined features; for many collinear features.
• Generalized linear models (GLMs): same linear core, different output link — logistic regression (next chapter) for yes/no, Poisson for counts.

Remember: Log left → % in y, log right → % in x, both → elasticity. k categories → k - 1 dummies. Center before multiplying. L1 kills, L2 shrinks — standardize first, tune \lambda by cross-validation.

Next question → We have the whole theory. Can we build it ourselves and match a library to the last digit? (Section 11)

11. Build It From Scratch

The question: can we rebuild everything from scratch and prove we understand it?

The 30-second version: add a column of 1s, solve with lstsq, compute residuals, estimate noise with n - p, then standard errors, t, p and intervals. About 50 lines of NumPy reproduce the statsmodels summary exactly. First read: run the code once.

One class, NumPy and SciPy only, that does everything this chapter derived: the closed form, gradient descent, and the full inference table. Every line maps to a section above.

import numpy as np
from scipy import stats

class LinearRegressionScratch:
    """Ordinary least squares, two ways, plus the inference statsmodels gives you."""

    def fit_closed_form(self, X, y):
        X1 = np.column_stack([np.ones(len(X)), X])          # add the intercept column
        self.w, *_ = np.linalg.lstsq(X1, y, rcond=None)      # QR/SVD under the hood, never inv()
        return self._diagnostics(X1, y)

    def fit_gradient_descent(self, X, y, lr=0.1, steps=2000):
        mu, sd = X.mean(0), X.std(0)
        Z = np.column_stack([np.ones(len(X)), (X - mu) / sd])  # scale, then add intercept
        w = np.zeros(Z.shape[1])
        for _ in range(steps):
            grad = -2 / len(y) * Z.T @ (y - Z @ w)              # the gradient from Section 4
            w -= lr * grad
        slopes = w[1:] / sd                                     # undo the scaling
        self.w = np.concatenate([[w[0] - slopes @ mu], slopes])
        return self._diagnostics(np.column_stack([np.ones(len(X)), X]), y)

    def _diagnostics(self, X1, y):
        n, p = X1.shape
        e = y - X1 @ self.w                                     # residuals
        rss, tss = e @ e, ((y - y.mean()) ** 2).sum()
        s2 = rss / (n - p)                                      # unbiased noise variance (Section 5)
        se = np.sqrt(np.diag(s2 * np.linalg.inv(X1.T @ X1)))    # small p, so inv is fine here
        t = self.w / se
        tcrit = stats.t.ppf(0.975, n - p)
        return {
            "w": self.w.round(4), "se": se.round(4), "t": t.round(2),
            "p": (2 * stats.t.sf(np.abs(t), n - p)).round(4),
            "ci95": np.column_stack([self.w - tcrit * se, self.w + tcrit * se]).round(3),
            "R2": round(1 - rss / tss, 4), "adjR2": round(1 - (rss / (n - p)) / (tss / (n - 1)), 4),
            "residuals_sum": round(e.sum(), 8), "X_dot_e": (X1.T @ e).round(8),
        }

X = np.array([[1000], [1200], [1500], [1800], [2100]], dtype=float)
y = np.array([310, 360, 380, 500, 530], dtype=float)
m = LinearRegressionScratch()
print(m.fit_closed_form(X, y))
print(m.fit_gradient_descent(X, y)["w"])

What it prints, and what to notice:

• Both methods give w = [100.8122, 0.2074] — two completely different routes, one answer.
• se, t, p, and ci95 match the statsmodels summary in Section 8 to every printed digit. You've rebuilt the summary table yourself.
• residuals_sum = 0 and X_dot_e = [0, 0]: the three free facts from Section 3, confirmed numerically.

Common mistake: running the gradient-descent version on raw square footage with lr = 0.1. The gradient is in the hundreds of thousands, the first step overshoots by orders of magnitude, and the weights hit infinity in a handful of steps. Delete the scaling line and watch it happen once — you'll never forget it.

Try it: add a second column that's size in square meters (\text{size} \times 0.0929). lstsq still returns an answer (the minimum-norm one), but np.linalg.inv(X1.T @ X1) in the diagnostics fails or returns garbage. That's perfect multicollinearity (Section 3) in four lines.

Remember: If you can write this class from memory, you understand linear regression: add a 1s column, solve with lstsq, residuals → s^2 with n - p → SE → t → p and CI.

Next question → How does this look in real code at work? (Section 12)

12. Use It Like a Pro

The question: how do practitioners actually use it?

The 30-second version: scikit-learn to predict (pipelines, cross-validation, regularization); statsmodels to explain (robust standard errors, full summaries). Then a checklist before anything ships. First read: the two-jobs table and the checklist.

Two jobs, two libraries

Prediction with scikit-learn

import numpy as np
from sklearn.datasets import fetch_california_housing
from sklearn.model_selection import train_test_split, cross_val_score
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import LinearRegression, RidgeCV, LassoCV
from sklearn.metrics import root_mean_squared_error, mean_absolute_error, r2_score

df = fetch_california_housing(as_frame=True).frame          # target: MedHouseVal, in $100k
X, y = df.drop(columns="MedHouseVal"), df["MedHouseVal"]
X_tr, X_te, y_tr, y_te = train_test_split(X, y, test_size=0.2, random_state=0)

models = {
    "OLS":   make_pipeline(StandardScaler(), LinearRegression()),
    "Ridge": make_pipeline(StandardScaler(), RidgeCV(alphas=np.logspace(-3, 3, 13))),
    "Lasso": make_pipeline(StandardScaler(), LassoCV(cv=5, random_state=0)),
}
for name, model in models.items():
    cv = cross_val_score(model, X_tr, y_tr, cv=5, scoring="neg_root_mean_squared_error")
    model.fit(X_tr, y_tr)
    pred = model.predict(X_te)
    print(f"{name:5s} CV RMSE {-cv.mean():.3f} | test RMSE {root_mean_squared_error(y_te, pred):.3f}"
          f"  MAE {mean_absolute_error(y_te, pred):.3f}  R2 {r2_score(y_te, pred):.3f}")

Three pro habits in those lines:

1. The scaler lives inside the pipeline. Cross-validation then re-fits it on each training fold, so test-fold statistics never leak into training.
2. Model selection uses cross-validation on the training set only. The test set is touched once, at the end.
3. Compare against simpler and regularized versions. If ridge and OLS score the same, you don't have a variance problem; if lasso zeros a feature, ask why.

Inference with statsmodels

import statsmodels.formula.api as smf

train = X_tr.assign(MedHouseVal=y_tr)
fit = smf.ols("np.log(MedHouseVal) ~ MedInc + HouseAge + AveRooms + AveBedrms + Latitude + Longitude",
              data=train).fit(cov_type="HC3")                 # heteroscedasticity-robust SEs
print(fit.summary())

Notes for reading it: with cov_type="HC3" the table shows z instead of t (robust SEs rely on large-sample theory). Log of the target means each coefficient reads as a percent change (Section 10). And in this dataset AveRooms and AveBedrms are strongly correlated, so check their VIFs before believing either coefficient's sign.

Before you ship a regression model

• [ ] Plotted residuals vs fitted, Q-Q, scale-location, leverage (Section 6)
• [ ] Asked Who, What, How, When of the data (Section 7)
• [ ] Scored on held-out data, against a baseline (Section 9)
• [ ] Validated on the most recent time period if time matters
• [ ] Checked residuals by important subgroup (region, price band) for systematic misses
• [ ] Used robust SEs if you report p-values or intervals
• [ ] Wrote down the valid input range — no extrapolation in production
• [ ] Set up monitoring for drift in inputs and residuals

Remember: sklearn predicts, statsmodels explains. Scaler inside the pipeline. Cross-validate on train, touch test once. Robust SEs by default when you report inference.

Next question → Where does this show up in the world, and how does it grow into an LLM? (Section 13)

13. In the Real World, and How It Connects to LLMs

The question: where is linear regression used, and how does it become an LLM?

The 30-second version: it's the industry baseline and the standard language for measuring effects. Its model, loss, optimizer, regularizer and evaluation are exactly the pieces an LLM uses, at a much larger scale. First read: the connections table.

Where linear regression earns its keep

Where it loses: on large tabular datasets with interactions and thresholds, gradient-boosted trees (Chapter 10) usually win on accuracy. On images, audio, and text, neural networks win. Linear regression stays in the toolbox because it's fast, explainable, a strong baseline, and the honest way to estimate an effect.

How this chapter becomes an LLM

Nothing in Part VII is new except the scale. Every piece of this chapter reappears:

Remember: Linear regression is the baseline in industry and the blueprint for deep learning. The model, the loss, the optimizer, the regularizer, and the evaluation in an LLM are this chapter, scaled up.

Next → The whole chapter on one page (Section 14). Then practice: the Derivations tab, then the Interview Bank.

14. The Chapter on One Page

One map of the whole chapter and every memory hook in one table. Use it before an interview, or to check a second read: if any row feels unfamiliar, go back to its section.

The memory map

Every hook in one place

Remember: if you can redraw the memory map on a blank page and say the hook for every row of the table, you have the chapter. Now practice.

15. Practice

The chapter ends here. Practice lives in two tabs, so this page stays a story you can read straight through.

• Ch 1 · Derivations — 13 derivation cards, each with the key move that unlocks it and the common slip; 18 formulas to memorize; and a mastery tracker.
• Ch 1 · Interview Bank — 70 questions in four levels, each answered in four beats (intuition, answer, why, what they'll ask next); two worked research-level problems with figures; 12 drills; and rapid-fire flashcards for the day before.

Suggested order: finish the chapter's second read → Derivations cards 1–4 → Interview Bank Levels 1–2 → cards 5–13 → Levels 3–4.