Chapter 3 — Logistic Regression and Classification: Predicting Yes or No

Chapter 1 drew a line to predict a number, a house price. Chapter 2 showed how to find the bottom of a loss that has no formula, and previewed the logistic loss. Now the answer is yes or no: will this applicant default on the loan? A line predicts a number. What predicts a probability? By the end of this chapter you can fit a logistic regression by hand and in code, read its coefficients correctly to a credit committee, test every assumption it rests on, correct for data that were sampled on the outcome or selected by an earlier decision, check that its probabilities mean what they say, turn them into approve/decline decisions, and extend the model to many classes and counts. One idea carries all of it: keep the line, but let it predict the log-odds, and choose the line that makes the yeses and noes you actually saw as probable as possible.

How to read it. The chapter is a chain of 15 questions in three acts: Build it (Sections 1–5), Trust it (Sections 6–9) and Use it (Sections 10–15). Each section opens with its question and a 30-second version, and ends with the question it leaves open.

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. Why Not Draw a Line Through the 0s and 1s?

The question: ten past applicants, four of whom defaulted. Can Chapter 1's straight line tell us Applicant K's chance of default?

The 30-second version: fit a line to 0/1 labels (the linear probability model, LPM) and it predicts 0.91 for Applicant A and −0.11 for Applicant J: a "probability" below zero. The line also assumes the noise is the same size everywhere, but a 0/1 outcome's variance is p(1 - p), largest at 50% and tiny near 0 or 1. And the true relationship is S-shaped: risk can't keep falling once it is near zero. Three failures, remembered as B-S-S: Bounds, Spread, Shape. We need a function that bends a score into a probability. First read: all of it.

The problem

You work on the credit-risk team of a fintech lender that makes unsecured personal loans. The outcome is default: y = 1 if the loan was charged off within 12 months of booking, y = 0 if it was repaid. The credit committee keeps a review sample of ten past files:

DTI is debt-to-income: monthly debt payments as a percent of monthly income. The feature we use by hand is x, credit score in "hundreds of points above 680". A new applicant, Applicant K, has score 680 (so x = 0), DTI 36% and income $70k. What is K's chance of default?

Lower scores default more, but not perfectly: C (620) repaid and F (690) defaulted. We'll use these ten files for every hand computation in the chapter. For code and diagnostics we add a simulated book of 20,000 applicants with a known true model (Section 6 introduces it), so we can check every method against the truth.

What the line does

Treat y as a number and fit Chapter 1's least squares (Ch 1, Section 3). The fit is \hat y = 0.382 - 0.444x, so K (x = 0) gets 0.382. The other predictions:

Three failures: B-S-S

Chapter 1's reasons for squaring the errors were S-B-S. The line fails yes/no labels for three reasons, B-S-S:

• Bounds. A line has no floor and no ceiling. J gets −0.106, a negative probability; extend the line to a 530 score and it passes 1. You cannot price a loan off a probability of −11%.
• Spread. A 0/1 outcome with probability p has variance p(1 - p): 0.25 at p = 0.5 but only 0.09 at p = 0.1. The noise is built to be unequal, so Chapter 1's assumption E (equal variance) fails by construction, and the usual standard errors are wrong.
• Shape. Real risk is S-shaped. Going from a 560 score to a 600 score matters a lot; going from 780 to 820 barely matters, because risk is already near zero and can't fall much further. A line forces the same change everywhere.

When people still use the line. Economists often fit the LPM on purpose, for average effects with a randomized treatment. The slope is then a clean average change in probability, it is easy to explain, and it works well when fitted values stay between about 0.2 and 0.8. Two conditions: report heteroscedasticity-robust standard errors (Ch 1, Section 8), and don't use it to produce individual probabilities. Ch 1 Q8 met the same −0.1 "probability"; this chapter is the fix.

Remember: a line has no floor or ceiling; probabilities do. The line fails 0/1 labels on B-S-S: Bounds, Spread, Shape.

Next question → How do we turn an unbounded score into a probability without giving up the line? (Section 2)

2. From a Score to a Probability: Sigmoid, Odds and Log-Odds

The question: how do we turn the lender's linear risk score into a number between 0 and 1?

The 30-second version: climb the P-O-L ladder: Probability → Odds = p/(1 - p) → Log-odds = \ln(\text{odds}). Log-odds run from -\infty to +\infty, so a line can predict them. Walk back down with the sigmoid, \sigma(z) = 1/(1 + e^{-z}). The boundary where p = 0.5 is where the score is 0: a straight line (a hyperplane), so logistic regression is a linear classifier with a curved probability. A second picture: a hidden "risk" equal to the score plus logistic noise crosses zero → default; Gaussian noise instead gives probit. First read: the ladder table, Figure 3.2, the decision boundary.

The ladder

A probability lives in [0, 1]. A line lives on (-\infty, +\infty). Two moves stretch one onto the other:

K's numbers come from the fitted model in Section 4. Odds of 1 mean 50/50; odds below 1 mean "less likely than not". Log-odds are 0 at 50/50, negative below it, and symmetric: p = 0.1 and p = 0.9 have log-odds −2.197 and +2.197.

The sigmoid: walking back down

The inverse of the logit is the sigmoid (also called the logistic function):

\sigma

the sigmoid (logistic function): squashes any score into a probability between 0 and 1

z

any score on the log-odds scale, from -\infty to +\infty (K's is -0.824)

e

Euler's number, about 2.718, the base of the natural log

\text{logit}(p)

the log-odds of p: turn it into odds p/(1-p), then take the natural log

p

a probability of default, between 0 and 1

\sigma^{-1}

the inverse of the sigmoid: takes a probability back down the ladder to log-odds

Three facts used again and again (Card 1); the first makes every gradient in this chapter clean, the third is Section 5's divide-by-4 rule:

\sigma'

the sigmoid's slope: how fast the probability rises as you nudge the score z

The model. Logistic regression predicts the log-odds with a line, then reads off the probability:

p_i

applicant i's probability of default

i

which applicant (row): A, B, …, J in the ten files

\eta_i

applicant i's score (linear predictor): the log-odds the line predicts

w_0

the intercept: the log-odds of default when every feature is 0 (a 680 score)

w_1

the weight on feature 1: log-odds added per unit of that feature

x_{i1}

applicant i's first feature, e.g. credit score in hundreds of points above 680

w_k

the weight on the last feature

x_{ik}

applicant i's last feature; k is the number of features

\mathbf x_i

applicant i's features stacked into one vector, with a 1 first for the intercept (bold = vector)

\mathbf x_i^\top \mathbf w

applicant i's score as a weighted sum (^\top = transpose: the column of features laid flat as a row)

\mathbf w

all the weights w_0, w_1, \dots, w_k stacked into one vector

P(y_i = 1 \mid \mathbf x_i)

the probability that applicant i defaults (y_i = 1), given their features

\eta ("eta") is the linear predictor, or score. It is the same weighted sum as Chapter 1; only what it predicts has changed.

The decision boundary is a straight line. Where is p = 0.5? Where \eta = 0. So the set of applicants the model calls a coin flip is w^\top x = 0: a line in two features, a plane in three, a hyperplane in general. Logistic regression is a linear classifier: the boundary is straight even though the probability bends. Every threshold t gives another parallel boundary, \eta = \operatorname{logit}(t). Curved boundaries need curved features (x^2, splines, interactions), exactly as in Ch 1, Section 10.

The latent-variable view: logit and probit

A second way to get the same model. Imagine each applicant has a hidden financial-stress number y^* = x^\top w + \varepsilon, and they default when y^* > 0. If \varepsilon has a logistic distribution (bell-shaped, with heavier tails than a normal and variance \pi^2/3 \approx 3.29), then P(y = 1 \mid x) = \sigma(x^\top w): logistic regression. If \varepsilon is standard normal, P(y = 1 \mid x) = \Phi(x^\top w): probit regression.

The two curves are nearly the same shape at different scales, so their coefficients differ by a near-constant factor. Theory gives three versions of that factor: matching the slopes at 0 gives \sqrt{8/\pi} = 1.60, the best overall fit gives 1.70, and matching standard deviations gives \pi/\sqrt{3} = 1.81. On our simulated book the logit/probit ratios are 1.79 to 1.92, and they grow as the base rate moves away from 50%, because then the fit relies more on the tails, where the two curves differ. The latent view matters again in Section 8: the noise sets the scale of every coefficient.

Vocabulary you'll be asked about

Remember: climb the P-O-L ladder with p/(1 − p), then ln; come down with the sigmoid. Linear in log-odds, S-shaped in probability, straight boundary.

Next question → Many lines give valid probabilities. Which one is best: what is the loss? (Section 3)

3. The Loss: Maximum Likelihood Becomes Cross-Entropy

The question: which weights make the four defaults and six repayments we saw most probable, and why not just minimize squared error?

The 30-second version: each applicant contributes the probability the model gave to what actually happened: p if they defaulted, 1 − p if they repaid. Multiply, take logs, flip the sign: that is the log loss, also called binary cross-entropy. It's −log of the probability you gave the truth, so a confident wrong answer costs without limit. Squared error through a sigmoid is non-convex and its gradient vanishes exactly when the model is confidently wrong. Why cross-entropy: the 3 Cs, Convex, Confident mistakes punished, Calibrated. Same move as Ch 1: OLS was maximum likelihood under Gaussian noise; this is maximum likelihood under a coin flip. First read: the likelihood in words, the per-point loss table, the 3 Cs.

Likelihood in words

Treat each applicant as a biased coin with its own probability of landing "default". For any weights w, ask: how probable would these exact ten outcomes be? That number is the likelihood, and maximum likelihood picks the w that makes it largest.

Start at w = 0. Every p is 0.5, so the probability of the observed ten outcomes is 0.5^{10} = 0.000977. Its log is 10 \times \ln 0.5 = -6.93, a mean loss of \ln 2 = 0.693 per applicant. Any useful model must beat that. The best intercept-only model (p = 0.4 for everyone) scores −6.73, and the fitted model in Section 4 reaches −4.36.

From the coin flip to the loss

A Bernoulli outcome y \in \{0, 1\} with probability p has probability p^y(1 - p)^{1-y}: that is p when y = 1 and 1 - p when y = 0. Multiply over independent rows and take logs (Card 2):

L(\mathbf w)

the likelihood: how probable the ten observed outcomes are under weights \mathbf w

\mathbf w

the weights (intercept and slopes) as one vector; bold = vector

\prod_{i=1}^n

multiply the term after it over every applicant, i = 1 to n

i

which applicant (row)

n

the number of applicants (10 in the toy)

p_i

the model's probability that applicant i defaults

y_i

applicant i's outcome: 1 if they defaulted, 0 if they repaid

\ell(\mathbf w)

the log-likelihood: the log of L, a sum instead of a product (−4.36 at the fit)

\sum_{i=1}^n

add up the term after it over every applicant, i = 1 to n

Maximizing \ell is the same as minimizing its negative average, the log loss or binary cross-entropy:

\frac1n

divide by n, the number of applicants, to get an average

\text{BCE}(\mathbf w)

binary cross-entropy (log loss): the average −log of the probability given to what happened

\sigma

the sigmoid: turns a score into a probability between 0 and 1

\mathbf x_i

applicant i's features as a vector, with a 1 for the intercept

\mathbf x_i^\top \mathbf w

applicant i's score: features times weights, summed (^\top = transpose)

In one line: log loss = −log of the probability you gave to what happened. Per point:

Squared error can never exceed 1. Log loss goes to infinity as a confident prediction turns out wrong.

Why not squared error?

Push the sigmoid through squared error and two things break (Card 7):

1. The gradient goes deaf. The derivative of (y - \sigma(\eta))^2 with respect to \eta is -2(y - \sigma)\sigma(1 - \sigma). The \sigma(1 - \sigma) factor goes to zero when the model is confident. Take a repayer (y = 0) the model scores at p = 0.99: cross-entropy's gradient on \eta is p - y = 0.99; squared error's is 2 \times 0.99 \times 0.99 \times 0.01 = 0.0196, fifty times weaker. The worst mistakes get the smallest corrections.
2. The bowl gets bumps. Squared error through a sigmoid is not convex in w. On our ten files, with the intercept at its best value, the squared-error curve in w_1 is concave over the whole stretch from −15 to −4.3 and again from 1.1 to 5, where it flattens into a plateau. At w_1 = -15 its slope is −0.0048, against −0.087 for cross-entropy: almost no signal pointing home. Cross-entropy is a bowl everywhere.

Why cross-entropy: the 3 Cs

• Convex. Its Hessian is X^\top WX \succeq 0 (Section 4), so there is one bottom and no local minima.
• Confident mistakes punished. -\ln p grows without bound; the gradient p - y is largest exactly when the model is most wrong.
• Calibrated. It is the likelihood, so its minimizer is the maximum-likelihood fit, and its score equations force the predicted probabilities to add up to the observed defaults (Section 9). Cross-entropy is a proper scoring rule: in expectation it is minimized by reporting the true probability (Ch 5).

Entropy, cross-entropy and KL in one paragraph. The entropy of a coin with probability \pi is H(\pi) = -\pi \ln \pi - (1 - \pi) \ln(1 - \pi): the loss of the best constant prediction. On the simulated book's training data that is H(0.121) = 0.369 nats; the fitted model gets 0.287 on the test split, and the true probabilities themselves get 0.2867. Cross-entropy H(q, p) = -\sum q \ln p is the average loss when the truth is q and you predict p. It splits as H(q, p) = H(q) + \operatorname{KL}(q \Vert p), and H(q) doesn't depend on the model, so minimizing cross-entropy = minimizing the KL divergence from the data to the model. That sentence is also the definition of LLM pretraining (Section 15).

The same loss with labels in \{-1, +1\}. ML papers often write the labels as \tilde y = 2y - 1 \in \{-1, +1\} and the loss as the logistic loss of the margin \tilde y\eta:

y

one applicant's label: 1 for default, 0 for repaid

\eta

that applicant's score (log-odds), \mathbf x^\top \mathbf w

e

Euler's number, about 2.718

\tilde y

the same label recoded as +1 (default) or -1 (repaid); tilde = recoded version

Same loss, different bookkeeping. A large positive margin (right side, confident) costs almost nothing; a negative margin costs about \lvert\eta\rvert. The hinge loss of SVMs, \max(0, 1 - \tilde y\eta), is its piecewise-linear cousin (Ch 11).

Remember: log loss = −log of the probability you gave to what happened. It's maximum likelihood for coin flips. The 3 Cs: Convex, Confident mistakes punished, Calibrated. Squared error through a sigmoid is non-convex and deaf to confident mistakes.

Next question → There's no formula for the minimum. How do we find it, and how fast? (Section 4)

4. Fitting It: Gradient, Hessian, and Newton = IRLS

The question: there's no closed form. How do we find the weights for our ten applicants, and how many steps does it take?

The 30-second version: the gradient is X^\top(p - y): each feature column weighted by how wrong we are on each row, exactly the shape of OLS's X^\top(\hat y - y). Setting it to zero gives equations that are nonlinear in w, so there is no closed form. The Hessian X^\top WX with W = \operatorname{diag}(p(1 - p)) is positive semi-definite, so the loss is a bowl with one bottom (if it exists: Section 6). Newton's method on this bowl is a weighted least-squares fit to a "working response", repeated: IRLS, weigh, regress, repeat. On the ten files, 5 Newton steps land within 1.2\times 10^{-8} of the answer; gradient descent on standardized x needs 51 steps at \eta = 2 and 191 at \eta = 0.5. At scale, use L-BFGS or SGD (Ch 2). First read: the gradient in words, Figure 3.4, the IRLS recipe.

The gradient

Differentiate the total log loss -\ell(w). The chain rule and \sigma' = \sigma(1 - \sigma) collapse everything (Card 3):

\frac{\partial(-\ell)}{\partial \eta_i}

how fast the total log loss changes as you nudge applicant i's score, holding the rest fixed: a slope, not a fraction to divide

\ell

the log-likelihood, so -\ell is the total log loss over the ten files

\eta_i

applicant i's score (log-odds), \mathbf x_i^\top \mathbf w

i

which applicant (row)

p_i

the model's probability that applicant i defaults

y_i

applicant i's outcome: 1 = defaulted, 0 = repaid

\nabla_{\mathbf w}(-\ell)

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

\mathbf w

the weights (intercept w_0 and slope w_1) as one vector

\sum_i

add up over every applicant i

\mathbf x_i

applicant i's features as a vector: a 1 for the intercept, then x_i

X

the feature matrix: one row per applicant, first column all 1s, second the score x

X^\top

the transpose of X: rows become columns, so each row of it is one feature across all applicants

\mathbf p

all ten predicted probabilities stacked into one vector (bold = vector)

\mathbf y

all ten 0/1 outcomes stacked into one vector

Read it as in Chapter 1: features × misses. Ch 1's least-squares gradient was X^\top(\hat y - y) up to the factor 2/n; here \hat y is replaced by p. Rows the model gets badly wrong pull hardest, and the pull is in the direction of their features.

By hand at w = 0, every p = 0.5, so p - y is −0.5 for the four defaulters and +0.5 for the six repayers:

x_i

applicant i's credit score in hundreds of points above 680

using \sum x = -0.4 and \sum xy = -2.3. Both entries are positive, so the loss falls by moving w_0 and w_1 down: fewer defaults overall, and fewer at high scores.

Why there is no closed form. Set the gradient to zero: X^\top(\sigma(Xw) - y) = 0. In OLS the unknown w appeared linearly (X^\top Xw = X^\top y), so one linear solve finished the job. Here w sits inside \sigma, so these are nonlinear equations in w. Ch 1's chain was quadratic loss → linear gradient → one solve; the log loss isn't quadratic, so the chain breaks at the first link. We iterate.

The Hessian, and why the loss is a bowl

Differentiate once more (Card 4):

\nabla^2_{\mathbf w}(-\ell)

the Hessian: how each slope of the log loss changes as you nudge each weight; the loss's curvature

\mathbf x_i \mathbf x_i^\top

applicant i's feature vector times itself: a small square matrix (outer product)

X^\top W X

the Hessian as a matrix product: the features, each row weighted by its p_i(1 - p_i)

W

a diagonal matrix holding each applicant's weight p_i(1 - p_i), the Bernoulli variance

p_1

the first applicant's (A's) predicted probability

p_n

the last applicant's predicted probability

n

the number of applicants (rows), 10 in the toy

For any direction v, v^\top X^\top WXv = \sum p_i(1 - p_i)(x_i^\top v)^2 \ge 0, so the Hessian is positive semi-definite and the log loss is convex (Ch 2, Section 3). If X has full column rank it is strictly convex, so there is at most one minimum. At most: Section 6 shows data where the bottom sits at infinity. At w = 0 on our toy, W = 0.25\,I and

\sum

add up over all ten applicants

Newton's step is a weighted regression: IRLS

Newton's method (Ch 2, Section 9) replaces the loss with its local quadratic and jumps to that quadratic's bottom: w \leftarrow w - H^{-1}g. For logistic regression the step rearranges into a weighted least-squares problem (Card 5):

\mathbf w_{\text{new}}

the weights after one Newton (IRLS) step

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

the matrix inverse of the Hessian: undoes it, like dividing by the curvature

\mathbf z

the working response: each applicant's current score nudged by how wrong it is

W^{-1}

the inverse of W: divides each applicant's miss by its p_i(1 - p_i)

The working response z_i = \eta_i + (y_i - p_i)/(p_i(1 - p_i)) is the current score nudged by how wrong the row is, rescaled to the log-odds scale. Each Newton step is a weighted regression of z on X with weights p_i(1 - p_i), and the weights change every step. Hence the name, iteratively reweighted least squares, and the hook: weigh, regress, repeat.

Rows near p = 0.5 get the most weight; rows the model is already sure about get almost none. That is the Bernoulli's variance again (B-S-S's Spread), now used on purpose.

The first step by hand. At w = 0 every weight is 0.25 and every working response is (y - 0.5)/0.25 = \pm 2: +2 for defaulters, −2 for repayers. Solving with the Hessian and gradient above:

\mathbf w_1

the weights after the first Newton step

\mathbf 0

the starting weights: all zeros, so every applicant gets p = 0.5

Repeat and the path runs (−0.471, −1.774) → (−0.722, −2.498) → (−0.815, −2.743) → (−0.824, −2.766). The distance to the answer after each step is 0.991, 0.268, 0.023, 0.00019, then 1.2\times 10^{-8} after step 5: the number of correct digits roughly doubles each step (quadratic convergence). A rule that stops when every coordinate of the step is below 10^{-8} fires on step 7. The answer:

\hat w_0

the fitted intercept: K's log-odds of default at a 680 score (hat = estimated from data)

\hat w_1

the fitted slope: change in log-odds per 100 points of credit score

\ell(\hat{\mathbf w})

the log-likelihood at the fitted weights: the best achievable on the ten files

Gradient descent vs IRLS

Gradient descent (Ch 1, Section 4) takes w \leftarrow w - \eta\nabla(\text{mean loss}). Because p(1 - p) \le 1/4, the mean-loss Hessian is never larger than X^\top X/(4n), which gives a safe step size of 2/\lambda_{\max} of that matrix (Ch 2, Section 4). On standardized x that bound is \eta = 8. Stopping when the log-likelihood changes by less than 10^{-8}, gradient descent takes 191 steps at \eta = 0.5, 51 at \eta = 2 and 13 at \eta = 8; at \eta = 16 it oscillates forever without diverging. IRLS needs 5.

On the 20,000-applicant book IRLS converges in 8 iterations, and sklearn's L-BFGS in 8 on standardized features but 49 on raw score, DTI and income (Section 14).

Free facts from the score equations

At the MLE the gradient is zero, X^\top(y - \hat p) = 0: one equation per column, and like Ch 1's X^\top e = 0 it hands out free facts (Card 6). The intercept's equation gives \sum\hat p_i = \sum y_i: the fitted probabilities in Section 1's table add up to exactly 4, the number of defaults. Each feature's equation gives \sum\hat p_ix_{ij} = \sum y_ix_{ij}: on the toy, \sum\hat p_ix_i = -2.3 = \sum y_ix_i. So in training data the average prediction equals the observed default rate, overall and along every feature: calibration in the large, built in (Section 9).

Remember: gradient = X^\top(p - y), the OLS shape with p in place of \hat y: features × misses. Hessian X^\top WX \succeq 0, so convex. Newton = IRLS: weigh, regress, repeat, 5 steps on the toy.

Next question → The fit says −2.77 per 100 points of credit score. What does that number mean to the credit committee? (Section 5)

5. Reading the Coefficients: Odds Ratios and Marginal Effects

The question: \hat w_1 = -2.77 per 100 points of credit score. What do we tell the credit committee it means?

The 30-second version: on the log-odds scale it adds: +100 points lowers log-odds by 2.77. On the odds scale it multiplies: odds × e^{-2.77} = 0.063 per 100 points (× 0.758 per 10 points). On the probability scale it depends on where you start: the change is w\cdot p(1 - p), largest at p = 0.5, where it is w/4 (the divide-by-4 rule). Hook: log-odds add, odds multiply, probabilities depend. Trap: a coefficient is not a change in probability, and an odds ratio is not a risk ratio unless the outcome is rare. First read: the three-scales table and the traps.

Three scales, one coefficient

Move Applicant K from a 680 score to 780 (x from 0 to 1), holding everything else fixed:

Memory hook: log-odds add, odds multiply, probabilities depend. The hook is the P-O-L ladder read sideways: a constant step on the bottom rung is a constant factor on the middle rung and a varying distance on the top one.

The odds ratio, and its confidence interval

The odds ratio per unit of x is e^w. Per 100 points it is 0.063: each extra 100 points cuts the odds of default to 6.3% of what they were. Per 10 points it is e^{-0.2766} = 0.758, a 24% cut in the odds. Rescaling x divides the coefficient; it never changes the model.

For a confidence interval, build it on the coefficient's scale and exponentiate the endpoints. Odds ratios are skewed, so never add \pm to the odds ratio itself. On the ten files the Wald interval for \hat w_1 is [−6.12, 0.585], so the OR per 100 points runs from 0.0022 to 1.80. It includes 1 because the Wald test fails to reject here (Section 7). The profile-likelihood interval, [0.00049, 0.80], does not.

Marginal effects: the probability scale

The derivative of the probability with respect to a feature (Card 8):

\frac{\partial p}{\partial x_j}

how fast the default probability changes as you nudge feature j, holding the rest fixed: a slope, not a fraction to divide

p

the applicant's predicted probability of default

x_j

feature j, e.g. credit score in hundreds of points above 680

j

which feature

w_j

feature j's coefficient (log-odds per unit), e.g. \hat w_1 = -2.766

\text{AME}_j

average marginal effect of feature j: its probability slope averaged over all applicants

n

the number of applicants

\sum_{i=1}^n

add up over every applicant, i = 1 to n

i

which applicant

\hat p_i

applicant i's fitted probability of default (hat = estimated)

At K, \hat w_1p(1 - p) = -2.766 \times 0.305 \times 0.695 = -0.586 per 100 points: about −5.9 percentage points per 10 points, for an applicant near K. The divide-by-4 rule bounds it: no applicant's probability moves faster than \hat w_1/4 = -0.69 per 100 points, the rate at p = 0.5. Quick mental math for any logistic coefficient.

Two ways to summarize one number for the whole population:

• Marginal effect at the mean (MEM): plug the average applicant into p(1 - p). Toy: −0.610.
• Average marginal effect (AME): compute the effect for every applicant, then average. Toy: −0.389.

Report the AME to humans. The average applicant may not exist (an average of 0/1 features is no one), and because p(1 - p) is curved, the effect at the average isn't the average effect. On the book the gap is large: the probability at the average features is 6.8%, while the average probability is 12.0%.

On the 20,000-applicant book (Section 6), with z = (\text{score} - 680)/100, d = (\text{DTI} - 36)/10 and \ell = \ln(\text{income}/75):

Doubling income multiplies the odds by e^{-0.544 \times \ln 2} = 0.686. On the toy the score coefficient is −2.77 against the book's −1.51: the ten files use score alone, so score also carries some of DTI's and income's effect (Ch 1's omitted-variable bias), and with an SE of 1.71 the toy estimate is very noisy anyway.

Odds ratio is not a risk ratio. The risk ratio (relative risk) is p_1/p_0; the odds ratio is [p_1/(1 - p_1)]/[p_0/(1 - p_0)]. They agree only when both probabilities are small, because then 1 - p \approx 1. For K at +100 points, the odds fall to 0.063 of their value but the probability falls from 30.5% to 2.7%, a risk ratio near 0.09. An OR of 1.84 per 10 DTI points does not mean "84% more likely to default". When the outcome is common, the OR always sits further from 1 than the risk ratio.

Other readings interviewers check

• The intercept is the log-odds when every feature is 0. On the toy, x = 0 is a 680 score, so \hat w_0 = -0.824 is exactly K's log-odds. With raw score as the feature, the intercept would be the log-odds at a score of 0, which no one has. Center features you want to read.
• Standardized coefficients (per one SD of each feature) compare features in different units. On the book: score −0.90, DTI +0.57, log income −0.26 per SD. Score matters most.
• Categorical features work exactly as in Ch 1, Section 10: k - 1 dummies, each an odds ratio against the baseline category; the dummy-variable trap is the same.
• Interactions do not carry over between scales. An interaction coefficient on the log-odds scale is not the interaction on the probability scale. A model with no interaction term still has probability effects that depend on the other features, through p(1 - p), and the interaction's effect on probability can even have the opposite sign to its coefficient (Ai and Norton, 2003). Compute the effect on the scale the question is asked on.

Two traps in all this. First, "the coefficient is −2.77, so default falls by 2.77 percentage points" mixes up the bottom and top rungs of the ladder. Second, "an OR of 0.063 means a 6.3% chance" turns a ratio of odds into a probability (Q3).

Remember: log-odds add, odds multiply, probabilities depend. The fastest a probability can move is w/4. Exponentiate CI endpoints, not \pm on the OR. Report AMEs to humans.

Next question → Every reading above assumes the model is right. What must be true, and how do we check? (Section 6)

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

The question: what must be true for these coefficients to mean anything, and how do we check on the 20,000-applicant book?

The 30-second version: the hook is L-I-N-E-S: Linearity in the log-odds, Independent rows, No extreme multicollinearity, Enough events per feature (the large-sample MLE needs them), and no perfect Separation. Two of Ch 1's L-I-N-E-A-R letters disappear, Normality and Equal variance, because the Bernoulli fixes the variance at p(1 - p). Exogeneity (R) carries over unchanged for causal claims. Residuals are binned, never a raw scatter of 0/1. Separation is the new one: if a line splits defaulters from repayers perfectly, the likelihood keeps rising forever and the coefficient runs to infinity. In our toy, DTI alone does it: every applicant above 39% DTI defaulted. Fix it with Firth's correction or a penalty. First read: the LINES table, Figure 3.5, separation. Second read: the test grid.

The simulated book

To test diagnostics you need data where you know the truth. The chapter's book is 20,000 simulated applicants (seed 303), drawn in this order:

The base rate is 12.0% (2,406 defaults). The MLE recovers (−2.62, −1.51, 0.61, −0.54), each within about one SE of the truth. Score and DTI are correlated at −0.51. Variants of this book, each with its own seed, test one failure at a time.

The mnemonic: L-I-N-E-S

What disappeared from Chapter 1: Normality and Equal variance. There is no error term to be normal, and the variance is p(1 - p) by definition, so it can't be "constant" or "wrong". What carries over: Ch 1's R, exogeneity. Every causal reading of a logistic coefficient needs E[\text{unobserved risk} \mid x] not to depend on x, and as in Ch 1 no residual plot can show it fails.

The ladder, what each letter buys:

Residuals for a 0/1 outcome

A raw residual y_i - \hat p_i is either 1 - \hat p_i or -\hat p_i, so plotting it against anything draws two curved bands, whatever the truth. Three fixes:

• Pearson residual: (y_i - \hat p_i)/\sqrt{\hat p_i(1 - \hat p_i)}, the miss in standard deviations.
• Deviance residual: \operatorname{sign}(y_i - \hat p_i)\cdot\sqrt{-2[y_i \ln \hat p_i + (1 - y_i) \ln(1 - \hat p_i)]}. Their squares add up to the deviance (Section 7).
• Binned residuals (the one to plot): sort rows by \hat p or by a feature, cut into about 20 bins, and plot each bin's mean of y - \hat p with a \pm 2 SE band of 2\sqrt{\sum\hat p(1 - \hat p)}/n_{\text{bin}}. About 95% of bins should sit inside the band.

See it, test it, fix it

As in Ch 1: plots first, tests second, and if p is low, the assumption must go.

Goodness of fit: Hosmer–Lemeshow. Sort rows into 10 groups by \hat p and compare observed with expected defaults:

\text{HL}

the Hosmer–Lemeshow statistic: total mismatch between observed and expected defaults across groups

\sum_{g=1}^{10}

add up over the 10 groups (deciles of \hat p)

g

which group: rows sorted by \hat p and cut into 10

O_g

the observed number of defaults in group g

E_g

the expected number of defaults in group g: the sum of its predicted probabilities

n_g

the number of applicants in group g

\sim

is distributed as

\chi^2_{8}

a chi-squared distribution with 8 degrees of freedom (10 groups minus 2)

\sum_{i \in g}

add up over every applicant i in group g

i

which applicant

\hat p_i

applicant i's fitted probability of default (hat = estimated)

On the book, \text{HL} = 10.4 with p = 0.24: no evidence of misfit. The test shares the weakness of every Ch 1 test. At 2 million rows, a model that misses by at most 0.6 percentage points in any decile gets p \approx 10^{-30}; at 20,000 rows the same misspecification gets p = 0.067. HL answers "can I detect misfit with this n?", not "is the misfit big enough to matter?". Plots first.

Linearity failing, live

The NL variant (seed 304) adds a hockey stick: risk accelerates above 43% DTI, \eta_{\text{NL}} = \eta + 2.5\cdot\max(0, (\text{DTI} - 43)/10). We fit the linear model anyway.

• See: binned residuals by DTI sit above zero below about 35% DTI, dip below it from the high 30s to the mid 40s, then shoot up above 48%: a hockey stick the straight line can't follow; 75% of the 20 bins fall outside the ±2 SE band (5% on the main book).
• Test: Box–Tidwell z = 21.3, p \approx 10^{-100}. A 3-df restricted-cubic-spline term in DTI improves the log-likelihood by \text{LR} = 455 (p \approx 10^{-98}). HL gives p \approx 10^{-27}.
• Damage: the linear fit compensates by steepening the DTI slope to 1.56 (truth 0.6, 2.6 times as steep) and is still wrong at the extremes: at 60% DTI it predicts log-odds 1.17 where the truth is 3.09.
• Fix: with the spline, HL rises to p = 0.75 and only 5% of bins fall outside the band.

Separation: when the MLE does not exist

Look at DTI in the ten files. The repayers have DTI 18, 22, 25, 30, 33 and 38; the defaulters 40, 42, 45 and 48. Every applicant above 39% defaulted, every one below repaid. That is complete separation: some line puts every defaulter on one side and every repayer on the other. (Quasi-complete separation is the same with some ties exactly on the line.)

Why it breaks the MLE: take the separating direction w^* and scale it up by c. Every defaulter's p rises toward 1 and every repayer's falls toward 0, so the likelihood keeps rising toward its ceiling of 1 (\log L \to 0) and never peaks. Along that ray on the toy, \log L goes −6.93 (c = 0), −3.90 (c = 1), −1.27 (c = 5), −0.256 (c = 20), −0.013 (c = 50). The "best" coefficient is infinite (Card 11).

What IRLS does: on DTI alone, \lVert w\rVert grows 1.6–1.8× per iteration for the first five, then by a constant \approx 10.8 per iteration, while \log L creeps to 0, the weights p(1 - p) underflow to 0, and X^\top WX becomes numerically singular after 41 iterations (on score + DTI: +8.2 per iteration, singular after 40). statsmodels on score + DTI returns coefficients (−58.5, 32.6, 151) with SEs of about 44,000, 54,000 and 106,000, plus a ConvergenceWarning and a PerfectSeparationWarning. The Wald z is tiny, because the SE grows faster than the coefficient: the test calls the strongest predictor in the data useless (Section 7).

Fixes:

• Firth's correction (Firth, 1993): maximize \ell(w) + \tfrac12 \ln\lvert I(w)\rvert, the log-likelihood plus the log of a Jeffreys prior. The penalty grows as the fit approaches the separating edge, so the estimate is always finite, and it removes the O(1/n) small-sample bias (Section 7). On the toy with d = (\text{DTI} - 35)/10, Firth gives (−0.827, 2.83); on score + DTI it gives (−0.563, −0.104, 2.18).
• L2 penalty (Section 12): also finite. sklearn's default C = 1 on standardized score and DTI gives (−0.695, −0.535, 1.22).
• Report a bound ("the odds ratio is at least …" from a profile likelihood), or merge the sparse categories that cause quasi-separation.

Influence. Ch 1's leverage carries over with weights: the logistic hat matrix is H = W^{1/2}X(X^\top WX)^{-1}X^\top W^{1/2}, and its diagonal is each row's leverage. Pregibon's \Delta\beta (the logistic Cook's distance) measures how far the coefficients move when a row is dropped. A rare defaulter with an excellent score has high influence: a single surprising default can tilt a small-sample fit.

Run every check in one go

import numpy as np
import statsmodels.api as sm
from scipy.stats import chi2
from scipy.optimize import linprog
from statsmodels.stats.outliers_influence import variance_inflation_factor

def make_book(seed=303, n=20000, hinge=0.0):
    """The chapter's simulated book (seed 303); seed 304 with hinge=2.5 gives the NL variant."""
    rng = np.random.default_rng(seed)
    score = np.clip(rng.normal(680, 60, n), 300, 850); z = (score - 680) / 100
    dti = np.clip(36 - 8 * z + rng.normal(0, 8, n), 0, 80); d = (dti - 36) / 10
    l = 0.25 * z + rng.normal(0, 0.45, n)
    rng.normal(0, 0.3, n)                                   # self-reported income noise (unused here)
    eta = -2.6 - 1.5 * z + 0.6 * d - 0.5 * l + hinge * np.maximum(0, (dti - 43) / 10)
    y = (rng.random(n) < 1 / (1 + np.exp(-eta))).astype(float)
    return z, d, l, dti, y

def hosmer_lemeshow(y, p, g=10):
    groups = np.array_split(np.argsort(p), g)
    stat = sum((y[i].sum() - p[i].sum())**2 / (p[i].sum() * (1 - p[i].mean())) for i in groups)
    return stat, chi2.sf(stat, g - 2)

def lines_checks(z, d, l, dti, y):
    X = sm.add_constant(np.column_stack([z, d, l]))
    m = sm.Logit(y, X).fit(disp=0); p = m.predict(X)
    # L: binned residuals by DTI (20 bins) and the Box-Tidwell term x ln x, x = DTI/10
    bins = np.array_split(np.argsort(dti, kind="stable"), 20)
    outside = np.mean([abs((y[b] - p[b]).mean()) > 2 * np.sqrt((p[b] * (1 - p[b])).sum()) / len(b) for b in bins])
    x = dti / 10; xlx = np.where(x > 0, x * np.log(np.where(x > 0, x, 1)), 0)
    bt = sm.Logit(y, np.column_stack([X, xlx])).fit(disp=0)
    # N: VIFs on the design; E: events per parameter; fit: Hosmer-Lemeshow
    vif = [variance_inflation_factor(X, j) for j in (1, 2, 3)]
    print(f"bins outside 2SE {outside:.2f} | Box-Tidwell p {bt.pvalues[-1]:.2g} | "
          f"VIF {np.round(vif, 2)} | EPV {y.sum() / 3:.0f} | HL p {hosmer_lemeshow(y, p)[1]:.3g}")

def separated(X, y):
    """S: is there a w with (2y - 1) * x'w >= 1 for every row? (a linear-programming feasibility check)"""
    s = 2 * y - 1
    res = linprog(np.zeros(X.shape[1]), A_ub=-(s[:, None] * X), b_ub=-np.ones(len(y)), bounds=(None, None))
    return res.status == 0

lines_checks(*make_book(303))                 # main book: 5% of bins outside, VIFs 1.47 / 1.36 / 1.11, EPV 802, HL p 0.237
lines_checks(*make_book(304, hinge=2.5))      # NL book: 75% of bins outside, Box-Tidwell p ~ 1e-100, HL p ~ 1e-27
toy_x = np.array([-1.2, -0.8, -0.6, -0.4, -0.2, 0.1, 0.3, 0.5, 0.8, 1.1])
toy_dti = np.array([48, 42, 30, 45, 38, 40, 25, 33, 22, 18.])
toy_y = np.array([1, 1, 0, 1, 0, 1, 0, 0, 0, 0.])
print("toy separated on score:", separated(sm.add_constant(toy_x), toy_y))       # False: C and F overlap
print("toy separated on DTI:  ", separated(sm.add_constant(toy_dti), toy_y))     # True: the 39% line

The EPV on the book is 802 (2,406 defaults over 3 slopes), so E holds with room to spare. The I check needs a group variable this book doesn't have; on real loan data, cluster by branch or by origination month.

Remember: L-I-N-E-S: Linearity in log-odds, Independent rows, No extreme multicollinearity, Enough events, no Separation. No N, no equal variance: the Bernoulli sets the spread. Bin your residuals. Separation sends \hat w to infinity: Firth or a penalty.

Next question → The assumptions hold. How sure are we about each coefficient, and why can two tests on the same data disagree? (Section 7)

7. Inference: Wald, Likelihood-Ratio and Score Tests

The question: is credit score's effect real in our ten files, and why does one test say yes and another say no?

The 30-second version: the MLE is approximately normal with variance (X^\top WX)^{-1}, the inverse of the Fisher information, the curvature of the log-likelihood. A Bernoulli carries information p(1 - p) per row on the log-odds scale, so rows near 50% teach the most. Three tests ask the same question three ways: Wald measures Width (how far \hat w is from 0 in SEs), LR measures Height (how much log-likelihood is lost by forcing w = 0), Score measures Slope (how steep the log-likelihood is at w = 0). Hook: W-H-S. They agree in big samples. On our ten files they don't: Wald p = 0.106, LR p = 0.029. Trust LR in small samples and whenever effects are large; Wald fails outright under separation. First read: Figure 3.6, the toy table. Second read: deviance, MLE properties, the summary.

Fisher information and standard errors

The Fisher information measures how sharply the log-likelihood curves around its peak: a sharp peak pins the parameter down, a flat one doesn't. For logistic regression it is exactly the Hessian of Section 4 (Card 9):

I(\mathbf w)

the Fisher information: how sharply the log-likelihood curves at weights \mathbf w

\mathbf w

the weights (intercept and slopes) as one vector

X

the feature matrix: one row per applicant, first column all 1s

X^\top

the transpose of X: rows become columns

W

a diagonal matrix of each applicant's p_i(1 - p_i): the information each row carries

\widehat{\text{Var}}(\hat{\mathbf w})

the estimated variance (and covariances) of the fitted weights across repeated samples

\hat W

W computed at the fitted probabilities (hat = estimated)

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

the matrix inverse of the information: undoes it, so sharp curvature means small variance

\text{SE}(\hat w_j)

the standard error of coefficient j: its typical wobble from sample to sample

\hat w_j

fitted coefficient j (hat = estimated from data)

j

which coefficient: 0 for the intercept, 1 for the score slope

\big[(X^\top \hat W X)^{-1}\big]_{jj}

the j-th diagonal entry of that matrix: coefficient j's variance

Compare Ch 1's \sigma^2(X^\top X)^{-1}: the weights p(1 - p) replace 1/\sigma^2. Each row's information on the log-odds scale is p(1 - p): 0.25 at 50%, 0.09 at 10%, 0.0099 at 1%. A book of 1% default rates needs about 25 times as many loans as a 50/50 book to pin down a coefficient equally well. For the intercept-only model, \operatorname{Var}(\operatorname{logit} \hat p) \approx 1/(n p(1 - p)): at a 12% base rate, pinning the log-odds to \pm 0.05 takes 14,551 loans (Drill 3).

On the ten files: \operatorname{SE}(\hat w_0) = 0.931 and \operatorname{SE}(\hat w_1) = 1.710.

The three tests: W-H-S

All three test H_0: w_1 = 0, and all are compared with a \chi^2 distribution with 1 degree of freedom (Card 10):

W = \left(\frac{\hat w_1}{\text{SE}(\hat w_1)}\right)^2

the Wald statistic: squared distance of the slope from 0, in standard errors

\hat w_1

the fitted score slope, −2.766 per 100 points

\text{SE}(\hat w_1)

its standard error, 1.710 on the toy

LR

the likelihood-ratio statistic: twice the log-likelihood lost by forcing w_1 = 0

\ell(\hat{\mathbf w})

the log-likelihood at the full fit (−4.361)

\ell(\hat{\mathbf w}_0)

the log-likelihood with w_1 forced to 0, intercept only (−6.730)

\hat{\mathbf w}_0

the restricted fit: intercept only, w_1 = 0

S

the score (Rao) statistic: squared slope at the restricted fit, scaled by the information there

U(\hat{\mathbf w}_0)

the score: the slope of the log-likelihood in w_1 at the restricted fit (-2.14)

I(\hat{\mathbf w}_0)

the Fisher information (curvature) at the restricted fit

Here \hat w_0 is the fit with w_1 forced to 0 (intercept only), U is the slope of the log-likelihood in w_1 at that point, and I is the information there.

Why they disagree. All three treat the log-likelihood as a parabola. Wald fits the parabola at the peak, score at zero, and LR doesn't need the parabola at all, it reads the actual height. On ten rows the log-likelihood is lopsided: it falls steeply toward w_1 = 0 but flattens toward -\infty. Wald's parabola is too wide, its SE too big, and its p too large.

Which to trust. LR, in small samples and with large effects. Wald has a known pathology, the Hauck–Donner effect (Hauck and Donner, 1977): as the true effect grows, the SE grows faster than the coefficient, so the Wald z falls. Separation is its endpoint: Wald's p goes to 1 for a perfect predictor (Drill 19). Wald's advantage is cost: one fit, and it comes free in every summary table. In large samples all three agree. On the book, testing DTI's coefficient gives Wald 393.6, LR 413.1 and score 406.6, all with p below 10^{-86}.

Deviance. The deviance is -2 \log L, measured from the saturated model, which fits every row perfectly. For 0/1 data the saturated \log L is 0, so the deviance is just -2\ell:

D

the deviance: twice how far the model's log-likelihood falls short of a perfect fit

\ell_{\text{sat}}

the log-likelihood of the saturated model that fits every row perfectly: 0 for 0/1 data

D_{\text{null}}

the deviance of the intercept-only model (13.46 on the toy)

D_{\text{model}}

the deviance of the fitted model (8.72 on the toy)

Toy: null deviance 13.46, residual deviance 8.72, drop 4.74, the LR statistic. Deviance plays the role of Ch 1's RSS: nested models are compared by the drop in deviance, which is \chi^2 with as many df as dropped parameters. For ungrouped 0/1 data, the deviance itself is not \chi^2(n - k) and says nothing about fit by itself; use HL or binned residuals.

Confidence intervals: Wald vs profile likelihood. A Wald interval is \hat w \pm 1.96 \operatorname{SE}: on the toy, [−6.12, 0.585]. A profile-likelihood interval collects every w_1 whose LR test against \hat w_1 isn't rejected, that is, where \log L is within 1.92 (half of \chi^2_1's 3.84) of its peak: on the toy, [−7.63, −0.223]. It is lopsided like the likelihood, excludes 0 (matching the LR test), and is what R's confint() reports for a GLM. Use it in small samples.

Properties of the MLE

Under L-I-N-E-S, as n grows the MLE is:

• Consistent: \hat w \to w.
• Asymptotically normal: \hat w \approx N(w, I(w)^{-1}). That is where every SE and Wald test comes from.
• Asymptotically efficient: no well-behaved (regular) consistent estimator has a smaller asymptotic variance (the Cramér–Rao bound; Ch 6 formalizes it).
• But biased away from zero in small samples. On a one-feature version of the book (true slope −2.03 once DTI and income are left out), 2,000 repeated samples give:

At n = 50 the average sample has only 5.9 defaults. The bias is the O(1/n) bias of every MLE, not separation: only 0.3% of the n = 50 samples separate. Firth's penalty removes the O(1/n) term (the n = 50 Firth mean is −2.01, 1% from the truth). This is Ch 1's lesson again: maximum likelihood estimates are not always unbiased.

Reading a Logit summary, line by line

The book's fit in statsmodels:

                           Logit Regression Results
==============================================================================
Dep. Variable:                      y   No. Observations:                20000
Model:                          Logit   Df Residuals:                    19996
Method:                           MLE   Df Model:                            3
Date:                Sun, 18 Oct 2026   Pseudo R-squ.:                  0.2065
Time:                        09:53:51   Log-Likelihood:                -5832.2
converged:                       True   LL-Null:                       -7350.4
Covariance Type:            nonrobust   LLR p-value:                     0.000
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
const         -2.6227      0.033    -78.607      0.000      -2.688      -2.557
z             -1.5117      0.052    -28.885      0.000      -1.614      -1.409
d              0.6123      0.031     19.839      0.000       0.552       0.673
l             -0.5441      0.054    -10.100      0.000      -0.650      -0.438
==============================================================================

\text{AIC} = 2k - 2\ell = 11{,}672.5 and BIC = 11,704.1 (Ch 1, Section 8); lower is better. On the toy McFadden R^2 is 0.35 and AIC 12.7.

Robust and clustered SEs work exactly as in Ch 1, Section 8 (fit(cov_type="cluster", cov_kwds={"groups": g})). For grouped binomial data with more spread than p(1 - p) allows (overdispersion), fit a quasi-binomial model, which scales the SEs by \sqrt{\text{Pearson} \chi^2/\text{df}}.

Remember: \operatorname{Var}(\hat w) \approx (X^\top WX)^{-1}, the inverse Fisher information; each row adds p(1 - p). W-H-S: Wald = Width, LR = Height, Score = Slope. Small n or big effects: trust LR and profile CIs. Small-sample MLE is biased away from zero; Firth fixes it.

Next question → All of this assumed the rows represent the applicants we'll score. Ours don't. (Section 8)

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

The question: our ten files over-represent defaults, our 20,000-loan book contains only approved applicants, some defaults were never recorded, and recent loans haven't had time to default. What does each do to the model?

The 30-second version: reuse Ch 1's Who, What, How, When. WHO (sampling on y): keep all defaults and a fraction of repayers and only the intercept moves, by \ln(s_1/s_0); subtract it back. Our toy's naive 30.5% for Applicant K becomes 5.6%. Class weights do the same thing to probabilities. WHO (selection): training only on approved loans hides how rejected applicants behave: the approved-only model predicts 10.9% default for rejects who truly default at 27.7% (reject inference). WHAT (omitted): leaving out a risk factor shrinks every coefficient even if it is independent of the others (non-collapsibility). HOW (labels): misrecorded defaults cap probabilities and flatten slopes; noisy self-reported income flattens its own coefficient (Ch 1). WHEN: a loan booked last month can't have defaulted yet. Hook: sample on y, shift the intercept. First read: case-control with the toy, reject inference, non-collapsibility.

WHO-1 · Sampling on the outcome: case-control and undersampling

The ten files aren't a random sample. They came from a batch of 48 applicants: the committee pulled all 4 defaults and 6 of the 44 repayers. The batch's default rate was 4/48 = 8.3%; the files show 40%. In Ch 1, selecting rows on y biased the slope (Ch 1, Section 7). In logistic regression it doesn't.

Why only the intercept moves. Bayes' rule on odds: posterior odds = prior odds × likelihood ratio. Sampling defaulters at rate s_1 and repayers at rate s_0 multiplies everyone's odds by the same factor s_1/s_0, whatever their x (Card 12):

P(y = 1 \mid \mathbf x, \text{sampled})

the chance an applicant with features \mathbf x is a default, among the files we kept

P(y = 0 \mid \mathbf x, \text{sampled})

the chance they repaid, among the files we kept

P(y = 1 \mid \mathbf x)

the chance an applicant with features \mathbf x defaults, in the full population

P(y = 0 \mid \mathbf x)

the chance they repay, in the full population

y

the outcome: 1 = default, 0 = repaid

\mathbf x

an applicant's features as a vector (bold = vector)

s_1

the share of defaulters kept in the sample (all of them: 1)

s_0

the share of repayers kept (6 of 44)

w_0 + \ln\frac{s_1}{s_0}

the new intercept: the population intercept w_0 plus the sampling offset

\mathbf x^\top \mathbf w

the slope part of the score: features times weights, summed (^\top = transpose)

A constant factor on the odds is a constant shift on the log-odds: the intercept absorbs it and every slope survives. Hook: sample on y, shift the intercept.

The toy, corrected. Here s_1 = 1 and s_0 = 6/44, so the offset is \ln(44/6) = 1.992. Subtract it:

\hat w_0^{\text{corrected}}

the fitted intercept with the sampling offset \ln(s_1/s_0) = 1.992 subtracted

p_K

Applicant K's corrected probability of default in the population

\sigma

the sigmoid: turns log-odds into a probability

Applicant K's 30.5% was an artifact of an enriched sample; in the population K's chance is 5.6%. Section 10 shows that this flips the decision.

The book, simulated. Keep every default and 20% of repayers: 5,896 rows with a 40.8% default rate. The fitted intercept is −0.985 and the slopes (−1.49, 0.574, −0.479) match the full-data fit. The intercept sits 1.638 above the full-data fit's, against \ln(1/0.2) = \ln 5 = 1.609 in theory. Over 200 repeated samples the shift averages 1.608 (SD 0.039) and the slopes average (−1.501, 0.604, −0.505), the truth. Uncorrected, the model's mean predicted default rate is 32.5%; after subtracting \ln 5 it is 11.9%, against a base rate of 12.0%.

In production, either fit with an offset of \ln(s_1/s_0) (statsmodels offset=), or correct after fitting (Card 13):

p_{\text{pop}}

the corrected probability for the real population of applicants

p_s

the probability output by the model fitted on the sample

\eta_{\text{pop}}

the corrected score (log-odds) for the population

\eta_s

the score from the model fitted on the sample

If you know the population base rate \pi but not the sampling rates, use s_1/s_0 = [\pi_s/(1 - \pi_s)]/[\pi/(1 - \pi)] (King and Zeng, 2001; Elkan, 2001). Ch 1 Drill 27 did this for one number; WP1 in the bank proves it.

What doesn't survive: a wrong model. The shift is exact only when the logistic form is right. Apply the same sampling to the NL book (fitted linearly) and the shift is 1.709, not 1.609, and the DTI slope moves by −0.25. Sampling changes which rows the misfit is averaged over.

WHO-2 · Class imbalance and class weights

class_weight="balanced" weights each defaulter by n/(2n_1) and each repayer by n/(2n_0): on the book, 4.16 and 0.568. A weight is a fractional copy of a row, so weighting is sampling in expectation, with s_1/s_0 = w_1/w_0 = (1 - \pi)/\pi. The intercept therefore shifts by \ln((1 - \pi)/\pi) = \ln 7.31 = 1.99 (simulated: 2.01), and the slopes stay put.

On the book, the mean predicted default rate becomes 38.9% against a base rate of 12.0%, so a finance forecast of expected loss built on these probabilities is inflated 3.2×. Recall at the 0.5 cutoff rises from 17% to 74%, not because the model got better but because every probability moved up. The ranking is untouched: AUC 0.8118 with weights, 0.8118 without, since a monotone shift can't reorder anyone.

Use weights when you only care about ranking, or when a rare class starves SGD of signal. Never use them for a probability-of-default model without recalibrating afterwards (Section 9). Moving the threshold (Section 10) does the same job without touching the probabilities.

WHO-3 · Selection by an earlier decision: reject inference

The model is trained on loans the lender approved, because only approved loans have outcomes. But approval depended on things the model can't see. In the RI variant (seed 307) an underwriter sees a signal u \sim N(0, 1) the model doesn't (verified employment, a phone call), the true log-odds are \eta - 0.8u, and an applicant is approved if z + 0.8u > -0.3.

• 61.7% are approved. They default at 5.1%; the rejects would default at 27.7%; through the door it's 13.8%.
• Among the approved, score and u are correlated at −0.38 (−0.01 in the full population): a low-score applicant was approved only if their u was strong. This is Ch 1's selection-on-a-collider (Berkson) pattern. So low-score approvals look safer than low-score applicants really are, and the score coefficient flattens from −1.35 (fit on everyone) to −0.92 (approved only).
• The approved-only model predicts 10.9% default for the rejects. Lower the cutoff on that basis and losses arrive at more than twice the forecast.

Fixes, ranked:

1. An explore arm. Approve a small random share of would-be rejects and weight them by 1/P(\text{approved}). Here 5% of rejects (394 loans) were approved at random. Their observed default rate, 27.2%, measures the rejects' risk directly. The reweighted model's score coefficient is −1.42, close to the through-the-door −1.35, but it predicts only 20.3% for rejects (bootstrap SE 1.4 points).
2. Outside data: bureau records of how rejected applicants performed on loans from other lenders.
3. Augmentation or parceling: score the rejects with the approved model and add them as fractional defaulters. This is valid only if selection was on observed features, and here it wasn't. Worse, it adds nothing: rows labeled with the model's own \hat p contribute exactly zero to the score equation, so the refit reproduces the approved-only model to the last digit.
4. Bounds: report the reject default rate under best- and worst-case assumptions.

The hard lesson: no model without u can reach 27.7% for the rejects, not even one fit on every applicant, which predicts 21.0%. The rejects were chosen by u. Only randomized outcomes on rejects measure them (Q28).

WHAT is missing · Non-collapsibility: the omitted variable that is independent

In Ch 1, an omitted variable biased a slope only if it was correlated with the included features (\text{bias} = \beta_2\delta). In logistic regression, leaving out any risk factor shrinks every coefficient toward zero, even one independent of all the others.

Why: the logistic model's scale is set by its unexplained noise (Section 2's latent variable). Omit a factor Z and its variation joins the noise. The latent noise is now bigger than the standard logistic, and the model, which must keep the noise at the standard logistic's variance, divides every coefficient by the extra spread. The probit approximation \sigma(t) \approx \Phi(t\sqrt{\pi/8}) turns that into a formula (Card 14):

\mathbb E_Z

the average over the omitted factor Z's values: the population-average probability

a

the score from everything kept in the model, e.g. -2.6 - 1.5z

\gamma

the omitted factor's true coefficient (1.0 in the NC variant)

Z

the omitted risk factor, independent of score; Z \sim N(0, 1.5^2)

\pi

the circle constant, 3.14159 (not a base rate here)

\sigma_Z^2

the variance of Z (2.25); here \sigma is a standard deviation, not the sigmoid

\beta_{\text{omit}}

a coefficient as fitted with Z left out

\beta

the same coefficient with Z in the model (e.g. −1.507 for score)

The NC variant (seed 309): Z \sim N(0, 1.5^2) independent of z, true logit -2.6 - 1.5z + 1.0Z, base rate 14.5%.

Both formulas are approximations, and neither is uniformly closer. The true ratio depends on the base rate: 0.84 at a 5.5% base rate, 0.77 at 15%, 0.73 at 34%, 0.72 at 50%. Here, at 14.5%, the latent-variance match (Ch 1 Drill 26) wins; near 50%, \pi/8 wins. The average marginal effect barely moves: risk differences are collapsible, odds ratios are not.

Two consequences interviewers probe. First, you can't compare odds ratios across models with different controls, even in a randomized trial: adding a prognostic covariate that is independent of treatment makes the treatment's log-OR larger in magnitude (Gail, Wieand and Piantadosi, 1984; WP2 and Q30). Neither is "biased": they answer different questions, conditional vs population-average. Second, report AMEs or risk differences when a stakeholder needs one number.

HOW was it measured · Misrecorded labels and noisy features

Label noise. Suppose a true default is recorded as "repaid" with probability \alpha_1 (collections sold the debt, the record never updated), and a true repayment is recorded as "default" with probability \alpha_0. Then the observed label follows (Card 15):

P(\tilde y = 1 \mid \mathbf x)

the chance the recorded label says default, for an applicant with features \mathbf x

\tilde y

the recorded label, which may be wrong (tilde = noisy version of the true y)

\alpha_0

the chance a true repayment is recorded as a default (0.02 in the LN variant)

\alpha_1

the chance a true default is recorded as repaid (0.15 in the LN variant)

\mathbf x^\top\mathbf w

the applicant's score: features times weights, summed

That curve lives in [\alpha_0, 1 - \alpha_1], not [0, 1], and it isn't logistic. In the LN variant (\alpha_1 = 0.15, \alpha_0 = 0.02) the observed probabilities are trapped between 2% and 85%. The naive fit flattens every slope (ratios 0.79, 0.80 and 0.86 against the clean fit) and still predicts up to 91.6% where the labels can't exceed 85%: the logistic shape can't represent a ceiling. If an audit gives you \alpha_0 and \alpha_1, maximize the likelihood of the formula above (Hausman, Abrevaya and Scott-Morton, 1998); here that recovers (−2.61, −1.53, 0.621, −0.489) against the truth (−2.6, −1.5, 0.6, −0.5). Otherwise, clean the labels on an audited subsample.

Noisy features. Self-reported income has reporting noise: in the book, \operatorname{Var}(\text{true log income})/\operatorname{Var}(\text{reported}) = 0.715. Fit on the reported value and its coefficient attenuates from −0.544 to −0.385 (ratio 0.71), Ch 1's noisy x flattens. The correlated score coefficient moves too (−1.51 → −1.55) as it picks up part of income's lost signal. Fix as in Ch 1: verified income, repeated measures, or a reliability correction.

WHEN was it collected · Immature labels and drift

The label "charged off within 12 months" isn't known until 12 months after booking. On the book, loans observed for only 3 months show 3.1% defaults and those observed for 6 months 6.2%, against 12% at full term. Train on them as if final and the model learns that recent applicants are safe. Fixes: a fixed performance window (train only on loans at least 12 months old), vintage analysis (default curves by booking month), or survival models that use partial follow-up properly (discrete-time logistic or cloglog, Section 11; Ch 14). Drift in the base rate is the calibration problem of Section 9.

The correction table

Remember: Who, What, How, When. Sample on y, shift the intercept; weights are sampling. Approved-only data understates rejects' risk, and only an explore arm measures it. Omitting even an independent factor shrinks logistic coefficients, not AMEs.

Next question → Even an unbiased model has to output probabilities we can use at face value. Does it? (Section 9)

9. Calibration: When the Model Says 10%, Is It 10%?

The question: among applicants the model scores at 10%, do 10% default, and what breaks that?

The 30-second version: a model is calibrated when P(\text{default} \mid \hat p = q) = q for every q. Check it with a reliability diagram: bin by \hat p, plot each bin's default rate against its average \hat p. An unpenalized logistic regression with an intercept is calibrated on average by construction, because its intercept's score equation forces \sum\hat p = \sum y, and along every feature (\sum\hat px = \sum yx). Three things break it, Shrink, Sample, Shift: regularization pulls probabilities toward the base rate, resampling or class weights inflate them, and a change in base rate moves them. Fix on held-out data: Platt is a line, isotonic is a staircase; for prior shift, correct the intercept. Calibration is about the probabilities, not the ranking: any monotone recalibration leaves AUC unchanged. First read: the definition, Figure 3.8, Shrink/Sample/Shift.

What calibrated means

Take every applicant the model scores between 8% and 12%. If about 10% of them default, the model is calibrated there. Formally:

P

the actual share (probability) of what's in the brackets

y

the outcome: 1 = default, 0 = repaid

\hat p(\mathbf x)

the model's predicted default probability for an applicant with features \mathbf x

q

any probability level, e.g. 10%

Calibration and discrimination are different properties. A model that predicts 12% for everyone is perfectly calibrated in the large and useless at ranking. A model that ranks perfectly but outputs \hat p^2 is a great ranker with wrong numbers. Lending needs both, because \text{expected loss} = \text{PD} \times \text{LGD} \times \text{EAD} (probability of default × loss given default × exposure at default) uses the PD as a number, not as a rank.

The reliability diagram. Sort the test rows by \hat p, cut into 10 equal-size bins, and plot each bin's observed default rate against its mean \hat p. A calibrated model sits on the diagonal. One number summarizes it: the expected calibration error, \text{ECE} = \sum (n_b/n)\cdot\lvert \text{observed}_b - \text{mean } \hat p_b\rvert. Ch 5 owns ECE and its pitfalls; here it is a ruler.

The calibration slope is the coefficient from regressing y on \operatorname{logit}(\hat p) in held-out data. A slope of 1 is ideal; above 1 means the predictions are too timid (squashed toward the mean), below 1 too extreme (overfit).

Why logistic regression is calibrated by construction

Section 4's score equations say that at the unpenalized MLE, \sum\hat p_i = \sum y_i and \sum\hat p_ix_{ij} = \sum y_ix_{ij} (Card 6). So on the training data the average prediction equals the default rate, overall and along every feature direction. That is calibration in the large, and it is why the log loss's third C is Calibrated. What it doesn't promise: calibration inside every bin if the model is misspecified (the NL model's top HL decile predicted 73% for a group that defaulted at 79.5%), or anything on data from a different population.

What breaks it: Shrink, Sample, Shift

Shrink: regularization. An L2 penalty pulls slopes toward zero, so predictions bunch toward the base rate. On the calibration variant (seed 312: 14,000 training rows, 3,000 calibration rows, 6,000 test rows; test default rate 12.2%):

The mean doesn't move, because sklearn's lbfgs does not penalize the intercept, so the intercept's score equation still holds: on training data, mean \hat p = \bar y = 0.119 at every C. The spread collapses and the slope climbs above 1. Two traps: at n = 14{,}000 the default C = 1 is practically no penalty, because sklearn weighs the penalty \tfrac12\lVert w\rVert^2 against C × the summed loss, so what matters is C\cdot n. And solver="liblinear" penalizes the intercept too: at C = 0.001 its mean training prediction is 0.217 against \bar y = 0.119, and even calibration in the large is gone.

Sample: resampling and class weights. Section 8's balanced weights moved the mean \hat p to 0.391 and the ECE to 0.269. Every bin sits far below the diagonal.

Shift: the base rate changes. The SHIFT variant (seed 313) is next year's book: the true intercept falls by 0.8 and the default rate from 11.9% to 6.6%. The old model still predicts 11.9% on average (ECE 0.053). Its ranking is untouched (AUC 0.8115), but every probability is too high. The prior-shift correction adds the change in log prior odds to every score (Card 13):

\eta_{\text{new}}

an applicant's corrected score (log-odds) for next year's book

\eta_{\text{old}}

the old model's score for that applicant

\pi'

the new base rate: 6.57% default next year

\pi

the old base rate the model was trained on: 11.88%

ECE falls from 0.053 to 0.006, and the AUC is still 0.8115. One caveat for Drill 17: the formula assumes label shift (P(x \mid y) fixed, only the base rate changes), while this variant changed P(y \mid x) with P(x) fixed. The true intercept change was −0.8, so a 20% prediction should become 10.1%, while the correction gives 11.5%. Close, not exact. Misspecification breaks calibration too; that's the fourth way.

The fixes

All of them are fitted on data the model was not trained on: a held-out calibration split, or cross-validation folds (sklearn's CalibratedClassifierCV(cv=5)).

• Platt scaling (Platt, 1999): fit a one-feature logistic regression of y on the model's logit, q = \sigma(a\cdot\eta + b). It is this chapter's model again. On the balanced-weights model it learns a = 0.981 and b = -2.035, almost exactly undoing the weights' +2.0 intercept shift. ECE: 0.269 → 0.014.
• Isotonic regression (Zadrozny and Elkan, 2002): fit the best non-decreasing step function from \hat p to y by pool-adjacent-violators. More flexible: it fixes any monotone distortion, not just a logit-linear one. Here it uses 23 steps and gets ECE 0.011. But it needs data: with only 300 calibration rows it set 45% of test predictions to exactly 0, and its log loss was 0.421 against Platt's 0.301. Hook: Platt is a line, isotonic is a staircase. Below about 1,000 calibration rows, prefer the line.
• Temperature scaling (Guo et al., 2017): divide every logit by one number T. Platt without the intercept, and the standard fix for multiclass and neural-network outputs (Section 15).
• Intercept correction for a known base-rate change: the formula above.

Recalibration with any monotone map leaves AUC unchanged: Platt 0.805, isotonic 0.805, balanced 0.805.

Remember: calibrated = "10% means 10%". Unpenalized logistic MLE is calibrated in the large by its own score equation. Shrink, Sample, Shift break it. Platt is a line, isotonic is a staircase: fit them on held-out data.

Next question → With honest probabilities, how do we turn 5.6% into approve or decline? (Section 10)

10. From Probability to Decision: Thresholds and Costs

The question: Applicant K has a 5.6% predicted chance of default. Approve or decline?

The 30-second version: a classifier is a probability plus a threshold. Decline when the expected loss from approving beats the expected profit forgone by declining: decline if p > t^* = c_{\text{FP}}/(c_{\text{FP}} + c_{\text{FN}}). Declining a good borrower (a false positive) costs $600 of lost profit; approving a defaulter (a false negative) costs $6,000 of lost principal. So t^* = 600/6{,}600 = 9.1\%, and K (5.6%) is approved. At the default 0.5 cutoff almost no one would be declined. The threshold needs calibrated probabilities, and moving it never requires refitting. First read: all of it.

The decision rule from expected cost

Call "predict default → decline" a positive. For an applicant with default probability p:

\mathbb E[\text{cost} \mid \text{approve}]

the average cost of approving this applicant, over many applicants like them

\mathbb E[\text{cost} \mid \text{decline}]

the average cost of declining them

p

the applicant's probability of default (K: 5.6%)

c_{FN}

cost of a false negative: approving a defaulter, $6,000 of lost principal

c_{FP}

cost of a false positive: declining a good borrower, $600 of lost profit

t^*

the cost-minimizing threshold, 9.1% here (star = the optimum)

Hook: cost of a false alarm over total cost (Card 16). With c_{\text{FP}} = $600 and c_{\text{FN}} = $6,000, t^* = 0.091. For K the expected costs are:

The sampling correction of Section 8 is the difference between declining and approving the same person. That's why calibration comes before thresholds.

The confusion matrix, at 0.5 and at t∗

The four outcomes on the book's 6,000-applicant test split:

The cost-optimal cutoff halves the cost, by declining 1,901 good borrowers to avoid 478 defaults. Scanning thresholds on the test split, the cheapest is 0.10 at $304.2: t^* from the formula is within a hair of it, as it should be when the probabilities are calibrated. Using the true probabilities at t^* would cost $308.4, so the fitted model is nearly as good as knowing the truth.

0.5 is a convention, not a rule. It is the cost-optimal cutoff only when the two errors cost the same and the probabilities are calibrated. With 10:1 costs, it approves 97% of applicants. Accuracy tells the same lie: 89.1% at 0.5 against 88.1% for approving everyone (Q5). Ch 5 owns precision, recall, ROC and PR curves; all are built from this matrix.

Three actions: approve, refer, decline. Add a manual review that costs $40 and catches 80% of the defaulters it sees. Each action's expected cost is now a line in p: approve 6,000p, refer 40 + 1,200p, decline 600(1 − p). The cheapest action switches at two thresholds: approve below 0.83%, refer between 0.83% and 31.1%, decline above 31.1% (Drill 16). With the refer band, K (5.6%) is referred, not approved.

When the base rate shifts (Section 9), either correct the probabilities and keep t^*, or keep the old probabilities and move the threshold by the same log-odds amount. Do one or the other, never both.

Remember: cut at cost of a false alarm over total cost: t^* = c_{\text{FP}}/(c_{\text{FP}} + c_{\text{FN}}). Calibrate first, then threshold. Moving a threshold never needs a refit.

Next question → What if there are more than two outcomes, or the outcome is a count? (Section 11)

11. Beyond Yes or No: Softmax, Ordinal Models and GLMs

The question: collections wants to predict prepaid / on time / defaulted, risk wants current / 30 days late / defaulted, and finance wants the number of missed payments. Does the same machinery work?

The 30-second version: for K unordered classes, give each class its own score line and turn scores into probabilities with the softmax, p_k \propto e^{\eta_k}. The loss is still −log of the probability given to what happened, and the gradient is still X^\top(P - Y). One class's weights are redundant, so one class is the reference. One-vs-rest fits K separate logistic models instead; their probabilities don't sum to 1. Ordered classes need only one line plus K - 1 cut-points (proportional odds). All of these are GLMs: a Distribution, a Line, and a Link. Poisson regression = Poisson + line + log, for counts, with \log(\text{exposure}) as an offset. Multiclass hook: one-vs-rest votes, softmax shares, ordinal climbs. First read: softmax, the GLM table.

Softmax and its gradient

Each class k gets a score \eta_k = x^\top w_k. The softmax turns K scores into K probabilities that are positive and sum to 1:

p_k

the predicted probability of class k (e.g. prepaid, on time, defaulted)

k

which class

e^{\eta_k}

e (about 2.718) raised to class k's score: always positive

\eta_k

class k's score from its own line, \mathbf x^\top \mathbf w_k

\sum_{j=1}^K

add up over all K classes

j

a class index that runs over every class

K

the number of classes (3 here)

\eta_j

class j's score

e^{\eta_j}

e raised to class j's score

p_y

the probability given to the true class

y

the observed outcome: the applicant's true class

\eta_y

the true class's score

\sum_j

add up over every class j

\frac{\partial(-\ln p_y)}{\partial \eta_k}

how fast the loss changes as you nudge class k's score, holding the rest fixed: a slope, not a fraction to divide

\mathbb 1[k = y]

1 if k is the true class, else 0

By hand, with scores (2, 1, 0): probabilities (0.665, 0.245, 0.090). If the true class is the second, the loss is -\ln 0.245 = 1.408 and the gradient on the scores is p - e_y = (0.665, -0.755, 0.090): push the true class up, push the others down in proportion to their probability. Stacked over rows, the gradient for the weights is X^\top(P - Y), features × misses again (Card 17).

Identifiability. Adding the same constant to every score changes nothing (scores (7, 6, 5) give the same probabilities), so one class's weights can be fixed at 0: the reference class. With K = 2 the softmax is exactly the sigmoid of the score difference: scores (1.3, 0.2) give 0.750 either way. Binary logistic regression is two-class softmax.

On the MC variant (seed 314: on time 63.2%, prepaid 27.5%, defaulted 9.2%; "on time" is the reference), the multinomial fit recovers the prepay line (−1.00, 0.835, −0.349) and the default line (−2.70, −1.64, 0.552), against truths (−1.0, 0.8, −0.3) and (−2.6, −1.5, 0.6): 6 parameters.

Multinomial vs one-vs-rest. One-vs-rest (OvR) fits K separate binary models ("prepaid vs not", …), each with its own normalization. On the MC data, OvR's three probabilities sum to anywhere from 0.92 to 1.73 for the same applicant (mean miss 0.069). Normalized, its log loss is 0.7828 against the multinomial's 0.7821 (truth 0.7811, base rates only 0.8708): nearly as good on log loss, but not a coherent probability model. In scikit-learn 1.9, lbfgs fits the true multinomial for more than two classes, and the old multi_class argument is gone (Section 14).

Ordered outcomes: proportional odds

Current < 30 days late < default has an order. The proportional-odds (ordinal logistic) model keeps one line and climbs it with cut-points: a latent severity s^* = x^\top w + logistic noise, and the class is decided by which cut-points s^* passes:

P(Y > k \mid \mathbf x)

the chance the outcome lands above class k (worse than it), given features \mathbf x

Y

the ordered outcome: current < 30 days late < default

\mathbf x

an applicant's features as a vector

\sigma

the sigmoid: turns a score into a probability

\mathbf x^\top \mathbf w

the one shared score line: features times weights, summed (^\top = transpose)

\theta_k

the cut-point the severity must pass to land above class k

\theta_1

the first cut-point (fitted 1.77)

\theta_2

the second cut-point (fitted 3.00)

\theta_{K-1}

the last cut-point: K classes need K - 1 of them

On the ORD variant (seed 315: cut-points 1.8 and 3.0, true slopes (−1.5, 0.6), shares 79.8% / 11.9% / 8.4%), the fit gives slopes (−1.52, 0.586) and cut-points (1.77, 3.00) with 4 parameters, against 6 for multinomial; AIC 21,689.5 vs 21,715.1. The assumption is that the slopes are the same at every cut: separate binary fits of Y \ge 1 and Y \ge 2 give (−1.520, 0.582) and (−1.536, 0.595), consistent with it. The Brant test formalizes that check. Hook: one-vs-rest votes, softmax shares, ordinal climbs.

The family: generalized linear models

Every model in this chapter, and OLS too, is one recipe: a Distribution, a Line, and a Link.

When the link is the distribution's canonical link (identity, logit, log for these three), the score equations are always X^\top(y - \mu) = 0. So every canonical GLM with an intercept is "calibrated on totals", \sum\hat{\mu}_i = \sum y_i, and its gradient is features × misses (Card 17).

Poisson regression for counts, with exposure

Missed payments depend on how long the loan has been on the books. A loan observed 24 months has four times the chances of one observed 6 months. Model the rate and put \log(\text{exposure}) in as an offset, a feature whose coefficient is fixed at 1:

\mathbb E[\text{count}_i]

borrower i's expected number of missed payments (the average over repeated samples)

i

which borrower

\text{months}_i

how many months borrower i's loan has been observed: the exposure

\mathbf x_i^\top \mathbf w

borrower i's score: the log of missed payments per month (^\top = transpose)

e^{\mathbf x_i^\top\mathbf w}

e (about 2.718) raised to the score: borrower i's expected missed payments per month

On the POIS variant (seed 316: months uniform on 6–24, true log-rate -3.0 - 0.9z + 0.4d), the fit recovers (−3.00, −0.890, 0.385); e^\beta gives rate ratios: 0.41 per 100 score points, 1.47 per 10 DTI points. \sum\hat{\mu} = \sum y holds to 10^{-11}. Two borrowers with identical features, observed 6 and 24 months, expect 0.30 and 1.19 missed payments (Drill 20).

Overdispersion. A Poisson's variance equals its mean. Check with Pearson \chi^2/\text{df}: 1.00 on Poisson data, 1.52 on negative-binomial data. Ignoring it, the naive Poisson SEs are too small: on the overdispersed data the true sampling SD of the score coefficient is 1.61× the naive SE. Fix with quasi-Poisson (SEs × \sqrt{\text{dispersion}}, here ×1.23), robust SEs (×1.48), or a negative binomial model (×1.37).

Remember: one-vs-rest votes, softmax shares, ordinal climbs. Softmax gradient = p - e_y. Every GLM = a Distribution, a Line, and a Link; a canonical link gives the gradient X^\top(y - \mu) and \sum\hat{\mu} = \sum y. Counts over time need a log-exposure offset.

Next question → With hundreds of candidate features, how do we stop overfitting, and is modeling the boundary the only way to classify? (Section 12)

12. Regularized Logistic Regression, and Generative vs Discriminative

The question: risk has 300 candidate features and the regulator wants a short, stable scorecard. And would modeling what defaulters look like do better than modeling the boundary?

The 30-second version: add Ch 1's penalty to the log loss. L1 kills, L2 shrinks (Ch 1): L1 gives a sparse scorecard; L2 keeps every estimate finite even under separation. In scikit-learn the knob is C = 1/\lambda: big C, little penalty, and it's on by default. The penalty is a prior: Gaussian → L2, Laplace → L1. It costs calibration (Section 9). Logistic regression is discriminative: it models P(y \mid x) directly. Generative classifiers (LDA, Naive Bayes) model P(x \mid y) and the base rate, then apply Bayes. If both classes are Gaussian with a shared covariance, LDA's posterior is a logistic regression; LDA is then a little more efficient, and logistic regression is safer when the Gaussian assumption fails. Hook: generative draws the classes, discriminative draws the line. First read: C, the L1 scorecard, the generative/discriminative table.

The penalized log loss

\min_{w_0,\,\mathbf w}

find the intercept and weights that make what follows smallest

w_0

the intercept (not penalized)

\mathbf w

the feature weights as one vector

C

sklearn's inverse penalty strength, 1/\lambda: big C, little penalty

\sum_{i=1}^n

add up over every applicant, i = 1 to n

i

which applicant

n

the number of training applicants

y_i

applicant i's outcome: 1 = default, 0 = repaid

p_i

the model's probability that applicant i defaults

\lVert\mathbf w\rVert_2^2

the squared length of the weight vector: the sum of squared weights (the L2 penalty)

\lVert\mathbf w\rVert_1

the sum of the weights' absolute sizes (the L1 penalty)

C multiplies the summed data loss, so it is 1/\lambda: big C, little penalty. Three rules from Ch 1, Section 10 still hold: standardize first (the penalty is unfair to features in small units), don't penalize the intercept, and choose C by cross-validation (Ch 4 owns the details). As a prior: the L2 fit is the MAP estimate under a Gaussian prior on w, the L1 fit under a Laplace prior (Ch 6).

An L1 scorecard

The L1 variant (seed 318) adds 20 pure-noise features and three near-copies of the score (correlation 0.986) to the book: 26 candidates. Along the L1 path from C = 0.001 to 10, the number of nonzero coefficients climbs from 2 to 26.

• CV's choice (C = 0.024, best 5-fold log loss) keeps 9: score, DTI, income, one score copy, and 5 noise features.
• A sparser rule (the smallest C whose CV log loss is within 0.001 of the best; C = 0.0092) keeps 4: score, DTI, income and one score copy, at a cost of 0.0005 in log loss.
• Stability, over 100 bootstrap refits at the sparser C: score, DTI and income are kept every time; noise features 4.6% of the time on average; the score copies 44%, 62% and 47% of the time, interchangeably. The model size ranges from 3 to 9 features.

L1 picks among near-duplicates almost arbitrarily (Ch 1's "picks one"), so selection is unstable even when predictions aren't. For a regulator-facing scorecard, merge or pre-select near-duplicates first, and report selection frequencies. Elastic net (L1 + L2) keeps correlated groups together.

Scorecard points. Credit scorecards convert log-odds into points with a fixed scaling, such as 600 points at good:bad odds of 50:1 and 20 points to double the odds (PDO):

Applicant K's corrected good:bad odds are 1/0.0598 = 16.7, so K scores 487.1 + 28.85 \times \ln 16.7 = 568 points (511 on the naive toy odds). Each feature's contribution is its coefficient × value × factor, which is how adverse-action reasons are generated (Section 15). Real scorecards also bin each feature into monotone weight-of-evidence (WOE) bins, so every bin's points move in the expected direction.

Generative vs discriminative

When they are the same model. If x \mid y = k \sim N(\mu_k, \Sigma) with a shared \Sigma, Bayes' rule gives log-odds that are linear in x (Card 18):

P(y = 1 \mid \mathbf x)

the chance an applicant with features \mathbf x defaults

P(y = 0 \mid \mathbf x)

the chance they repay

\mathbf x^\top

the applicant's features laid flat as a row (transpose: rows become columns)

\Sigma

the covariance matrix of the features, shared by both classes

\Sigma^{-1}

the matrix inverse of \Sigma: undoes it

\boldsymbol\mu_1

the average feature vector of defaulters (bold = vector)

\boldsymbol\mu_0

the average feature vector of repayers

\boldsymbol\mu_1^\top

\boldsymbol\mu_1 transposed into a row

\boldsymbol\mu_0^\top

\boldsymbol\mu_0 transposed into a row

\pi_1

the base rate of defaulters: the prior chance of default

\pi_0

the share of repayers, 1 - \pi_1

In the LDA variant (seed 317: one feature, x \mid y = 0 \sim N(0, 1), x \mid y = 1 \sim N(1, 1), \pi_1 = 0.2), theory says slope 1 and intercept -0.5 + \ln 0.25 = -1.886. On 50,000 rows logistic regression gets (−1.886, 1.004) and LDA (−1.885, 1.001). The cost of not assuming the Gaussian is small: at n = 200 the logistic slope's variance is 1.12× LDA's.

When they differ. Give the defaulters a wider spread, x \mid y = 1 \sim N(1, 2^2), and the log-odds become quadratic: the x^2 coefficient is \tfrac12(1/\sigma_0^2 - 1/\sigma_1^2) = 0.375 (QDA). Logistic regression with an x^2 feature recovers 0.382; without it, the linear fit is misspecified. Ch 8 covers Naive Bayes, LDA and QDA in depth.

Remember: big C, little penalty, and it's on by default. L1 kills (scorecards), L2 shrinks (stability, separation). Selection among near-duplicates is unstable. Generative draws the classes, discriminative draws the line; shared-covariance Gaussians make them the same model.

Next question → Can we build it ourselves and match the library to the last digit? (Section 13)

13. Build It From Scratch

The question: can we rebuild logistic regression and its inference in about 50 lines of NumPy?

The 30-second version: add a 1s column; loop: p = \sigma(Xw), W = p(1 - p), solve the weighted least-squares step; stop when the step is tiny. Then SEs from (X^\top WX)^{-1}, Wald z, and LR against the intercept-only model. Same answer as statsmodels on the toy; gradient descent gets there too, more slowly. First read: run the code once. Then check what it prints against the three bullets below it.

import numpy as np
from scipy import stats

class LogisticRegressionScratch:
    """Logistic regression two ways (Newton/IRLS and gradient descent), plus the inference statsmodels prints."""

    def fit_irls(self, X, y, tol=1e-10, max_iter=50):
        X1 = np.column_stack([np.ones(len(X)), X])                # add the intercept column
        w, norms = np.zeros(X1.shape[1]), [0.0]
        for it in range(1, max_iter + 1):
            p = 1 / (1 + np.exp(-X1 @ w))
            W = p * (1 - p)                                       # weigh
            try:                                                  # regress: the WLS fit of the working
                step = np.linalg.solve(X1.T @ (W[:, None] * X1), X1.T @ (y - p))  # response, as a step
            except np.linalg.LinAlgError:
                print(f"X'WX singular after {it - 1} iterations: separation, no MLE (Section 6)")
                return None
            w = w + step                                          # repeat
            norms.append(np.linalg.norm(w))
            if it == 1:
                print("first IRLS step:", w.round(3))
            if it == 10 and norms[-1] > norms[-2] + 1:            # still growing fast after 10 steps
                print("warning: |w| keeps growing -> separation? (Section 6)")
            if np.max(np.abs(step)) < tol:
                break
        self.w, self.n_iter = w, it
        return self._inference(X1, y)

    def fit_gd(self, X, y, lr=2.0, steps=5000):
        mu, sd = X.mean(0), X.std(0)
        Z = np.column_stack([np.ones(len(X)), (X - mu) / sd])     # standardize, then add the 1s
        w = np.zeros(Z.shape[1])
        for _ in range(steps):
            w -= lr * Z.T @ (1 / (1 + np.exp(-Z @ w)) - y) / len(y)   # gradient of the mean log loss
        slopes = w[1:] / sd                                       # undo the scaling
        return np.concatenate([[w[0] - slopes @ mu], slopes])

    def _inference(self, X1, y):
        loglik = lambda eta: np.sum(y * eta - np.logaddexp(0, eta))
        p = 1 / (1 + np.exp(-X1 @ self.w))
        se = np.sqrt(np.diag(np.linalg.inv(X1.T @ ((p * (1 - p))[:, None] * X1))))  # inverse Fisher info
        z = self.w / se
        ll, ll0 = loglik(X1 @ self.w), loglik(np.full(len(y), np.log(y.mean() / (1 - y.mean()))))
        lr = 2 * (ll - ll0)                                       # LR test against intercept-only
        return {"w": self.w.round(4), "se": se.round(4), "z": z.round(3),
                "wald_p": (2 * stats.norm.sf(np.abs(z))).round(4), "loglik": round(ll, 4),
                "loglik_null": round(ll0, 4), "LR": round(lr, 3), "LR_p": round(stats.chi2.sf(lr, X1.shape[1] - 1), 4),
                "mcfadden_R2": round(1 - ll / ll0, 4), "iterations": self.n_iter,
                "sum_p_minus_sum_y": round(p.sum() - y.sum(), 10),     # score equation, intercept
                "XT_y_minus_p": (X1.T @ (y - p)).round(10)}            # score equations, all columns

x = np.array([-1.2, -0.8, -0.6, -0.4, -0.2, 0.1, 0.3, 0.5, 0.8, 1.1])   # (score - 680) / 100
y = np.array([1, 1, 0, 1, 0, 1, 0, 0, 0, 0], dtype=float)               # 1 = defaulted
m = LogisticRegressionScratch()
print(m.fit_irls(x[:, None], y))
print("gradient descent:", m.fit_gd(x[:, None], y).round(4))

• first IRLS step: [−0.471, −1.774] is Section 4's hand computation; iterations: 7 is the step-below-tolerance rule firing on step 7, although step 5 was already within 1.2\times 10^{-8}.
• w = [−0.8241, −2.7658], se = [0.9311, 1.7098], Wald p = 0.1057, LR = 4.738 with p = 0.0295, McFadden R^2 = 0.352 match the statsmodels summary of the toy digit for digit. Gradient descent lands on the same w: two routes, one answer.
• sum_p_minus_sum_y = 0 and XT_y_minus_p = [0, 0]: the score equations of Section 4, the calibration in the large of Section 9.

Common mistake: run fit_irls with DTI as the only feature, m.fit_irls(((dti - 35) / 10)[:, None], y) with dti = np.array([48, 42, 30, 45, 38, 40, 25, 33, 22, 18.]). The first step looks harmless, (−0.251, 1.659). Then \lVert w\rVert grows 1.6–1.8× per iteration for five iterations and by a steady \approx 10.8 per iteration after that, \log L creeps toward 0, the weights p(1 - p) underflow, the growth warning fires at iteration 10, and X^\top WX becomes singular after 41 iterations. That's separation (Section 6), live. There is no MLE to make inference about, so the method returns None.

Try it: the ten files are a case-control sample (Section 8). Fit them, then subtract \ln(44/6) = 1.992 from m.w[0] and recompute K's probability: 0.305 becomes 0.0564.

Remember: add 1s, then weigh, regress, repeat; SEs from (X^\top WX)^{-1}; LR from two log-likelihoods; the score equations are your unit test.

Next question → How does this look in production code? (Section 14)

14. Use It Like a Pro

The question: how do practitioners fit, explain and ship logistic regression?

The 30-second version: scikit-learn to predict, statsmodels to explain (Ch 1's rule). Three defaults in sklearn's LogisticRegression bite: the 3 Ps: Penalized by default (L2, C = 1.0), Predict() cuts at 0.5, Plain lbfgs needs scaled features and enough iterations. First read: the 3 Ps and the checklist.

Predict with scikit-learn, explain with statsmodels

import warnings
warnings.filterwarnings("ignore", category=FutureWarning)   # LogisticRegressionCV 1.9 transition notices
import numpy as np
import statsmodels.api as sm
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import LogisticRegression, LogisticRegressionCV
from sklearn.calibration import CalibratedClassifierCV
from sklearn.metrics import log_loss, roc_auc_score

rng = np.random.default_rng(303)                       # the chapter's 20,000-applicant book
n = 20000
score = np.clip(rng.normal(680, 60, n), 300, 850); z = (score - 680) / 100
dti = np.clip(36 - 8 * z + rng.normal(0, 8, n), 0, 80); d = (dti - 36) / 10
income = 75 * np.exp(0.25 * z + rng.normal(0, 0.45, n)); rng.normal(0, 0.3, n)
y = (rng.random(n) < 1 / (1 + np.exp(-(-2.6 - 1.5 * z + 0.6 * d - 0.5 * np.log(income / 75))))).astype(int)
X = np.column_stack([score, dti, np.log(income)])      # raw units: points, %, log $k
tr, te = slice(0, 14000), slice(14000, None)           # first 14,000 train, last 6,000 test

models = {
    "MLE (C=inf)":   make_pipeline(StandardScaler(), LogisticRegression(C=np.inf)),
    "default C=1":   make_pipeline(StandardScaler(), LogisticRegression()),
    "CV over C":     make_pipeline(StandardScaler(), LogisticRegressionCV(Cs=10, cv=5, scoring="neg_log_loss")),
    "balanced+Platt": CalibratedClassifierCV(make_pipeline(StandardScaler(),
                      LogisticRegression(C=np.inf, class_weight="balanced")), method="sigmoid", cv=5),
}
for name, mdl in models.items():
    p = mdl.fit(X[tr], y[tr]).predict_proba(X[te])[:, 1]          # column 1 = classes_[1] = default
    print(f"{name:15s} log loss {log_loss(y[te], p):.4f}  AUC {roc_auc_score(y[te], p):.3f}  "
          f"mean p {p.mean():.3f}  approve rate at .predict() {np.mean(mdl.predict(X[te]) == 0):.3f}")
t_star = 600 / (600 + 6000)                                       # Section 10, not 0.5
p = models["MLE (C=inf)"].predict_proba(X[te])[:, 1]
print("approve rate at t* = 0.091:", np.mean(p <= t_star).round(3))

# statsmodels to explain: same MLE, with SEs, z, LR test and average marginal effects
Xs = sm.add_constant(np.column_stack([z, d, np.log(income / 75)]))
fit = sm.Logit(y, Xs).fit(disp=0)
print(fit.summary(xname=["const", "z", "d", "l"]))
print(fit.get_margeff(at="overall").summary_frame().round(4))

What to notice. With 14,000 rows, the default C = 1 gives the same log loss as the MLE (0.2868), because the penalty is weighed against C × the summed loss (Section 9). What matters is C\cdot n: on 400 rows a copied C = 0.003 holds the median top prediction to 22% (Q23). Every model approves 97.2% at .predict(), the 0.5 cutoff; the cost-optimal cutoff approves 58.5% (Section 10). Platt on top of balanced weights gives back calibrated probabilities (mean \hat p 0.121). And statsmodels reproduces Section 7's summary and Section 5's AMEs (−0.1306, 0.0529, −0.0470).

Gotchas (scikit-learn 1.9.1, statsmodels 0.15.0)

Before you ship a probability-of-default model

• [ ] L-I-N-E-S checked: binned residuals per feature, VIFs, EPV, separation (Section 6)
• [ ] Who, What, How, When asked of the training data: sampling, approval selection, label quality, label maturity (Section 8)
• [ ] Calibrated on held-out, recent data; reliability diagram in the model document (Section 9)
• [ ] Threshold set from costs, not 0.5; refer band if manual review exists (Section 10)
• [ ] Reject-inference strategy documented, ideally with an explore arm
• [ ] Adverse-action reason codes generated and checked for stability (Section 15)
• [ ] Fairness review for proxies and disparate impact
• [ ] Monitoring of the base rate, calibration and input drift after launch

Remember: sklearn predicts, statsmodels explains. The 3 Ps: Penalized by default, Predict cuts at 0.5, Plain lbfgs needs scaling. Scale, set the penalty on purpose, and threshold by cost.

Next question → Where does logistic regression earn its keep, and how does it grow into an LLM? (Section 15)

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

The question: where is logistic regression used, and why is an LLM's last layer the same model?

The 30-second version: it's the workhorse of credit scoring, click prediction and clinical risk scores because it is fast, calibrated and explainable. An LLM predicts the next token with a softmax over its vocabulary on top of a learned representation: multinomial logistic regression on learned features, trained with this chapter's cross-entropy. Reward models for RLHF are logistic regressions on the difference of two rewards. First read: both tables.

Where logistic regression earns its keep

At a 0.5% click rate, keeping 1% of non-clicks makes the training sample 33.4% clicks; a training-sample prediction of 0.5 maps back to q = 0.0099, a log-odds shift of \ln 0.01 = -4.61 (Q34).

Where it loses. With strong interactions and thresholds in tabular data, gradient-boosted trees (Ch 10) usually win on accuracy. Not always: on this chapter's hockey-stick NL book, a default HistGradientBoosting scored AUC 0.880 against the linear logistic model's 0.883. A common hybrid is a GBM for ranking with monotone constraints, or logistic regression on GBM-derived bins, to keep explanations and calibration (Q32).

How this chapter becomes an LLM

Remember: logistic regression is the scorecard, the click model, and the last layer of every LLM. Cross-entropy is the LLM loss; softmax is the LLM's output; a reward model is a logistic regression.

Next → The whole chapter on one page (Section 16). Then practice. Our classifier fits the ten files and the 20,000-applicant book. Will it work on next year's applicants? That's Chapter 4: Generalization.

16. 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.

17. Practice

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

• Ch 3 · Derivations: 18 derivation cards, from the sigmoid's derivative to LDA's logistic posterior, each with the key move that unlocks it (toolkit hook L-O-G-I-T: Likelihood, Odds, Gradient, Information, Taylor) and the common slip; about 22 formulas to memorize; and a mastery tracker.
• Ch 3 · Interview Bank: 39 questions in three levels, each answered in four beats (intuition, answer, why, what they'll ask next); two worked research-level problems with figures (case-control sampling, and the independent covariate in a randomized test); 20 drills; and 36 rapid-fire flashcards for the day before.

Suggested order: finish the chapter's second read → Derivations Cards 1–6 → Interview Bank Levels 1–2 → Cards 7–18 → Levels 3–4.