A probability cannot follow a straight line, because it must stay between 0 and 1. Its log-odds can. Logistic regression fits a straight line on the log-odds scale and bends it back into an S on the probability scale.
Code: data and style used in this chapter
import matplotlib.pyplot as pltimport numpy as npimport pandas as pdimport statsmodels.api as smimport statsmodels.formula.api as smffrom matplotlib.axes import Axesfrom scipy import statsfrom sklearn import datasetsfrom sklearn.linear_model import LogisticRegressionfrom statsmodels.base.model import LikelihoodModelResultsfrom statly import stylestyle.use()rng = np.random.default_rng(1)anes = sm.datasets.anes96.load_pandas().data # 944 voters, US presidential election 1996anes["dole"] = anes.vote # 1 = expected to vote for Dole, 0 = for Clintoniris = datasets.load_iris(as_frame=True).frame # flowers of three iris speciesSCALE = np.linspace(0, 8, 200)logistic = smf.logit("dole ~ selfLR", anes).fit(disp=0)shares = anes.groupby("selfLR").dole.agg(["mean", "count"])def inverse_logit(log_odds: np.ndarray) -> np.ndarray:"""Turn log-odds back into probabilities: 1 / (1 + e^(−log-odds))."""return1/ (1+ np.exp(-log_odds))
4.4.1 A line for a yes/no outcome
Before the 1996 US presidential election, 944 voters placed themselves on a scale from 1 (extremely liberal) to 7 (extremely conservative) and said whom they would vote for: Bill Clinton or Bob Dole. How does the chance of a Dole vote change along the scale?
Code the answer as 1 for Dole and 0 for Clinton. A straight line through these zeros and ones, the linear probability model, predicts probabilities below 0 and above 1 at the ends. Its noise is also far from normal and far from constant: each answer is 0 or 1, never something in between.
Code
straight = smf.ols("dole ~ selfLR", anes).fit()def plot_curves(ax: Axes) ->None:"""The share voting Dole at each step, with a line and a logistic curve.""" ax.scatter(shares.index, shares["mean"], s=shares["count"], color=style.BLUE, alpha=0.6) ax.plot( SCALE, straight.predict(pd.DataFrame({"selfLR": SCALE})), color=style.INK, ls="--", lw=1.5, label="straight line", ) ax.plot( SCALE, logistic.predict(pd.DataFrame({"selfLR": SCALE})), color=style.CORAL, lw=2.2, label="logistic", ) ax.axhspan(-0.4, 0, color=style.GREY, alpha=0.12, lw=0) ax.axhspan(1, 1.4, color=style.GREY, alpha=0.12, lw=0) ax.yaxis.set_major_formatter(style.PERCENT) ax.set_xticks(range(1, 8)) ax.legend(fontsize=8, loc="upper left") ax.set( xlabel="self-placement (1 = extremely liberal, 7 = extremely conservative)", ylabel="chance of a Dole vote", ylim=(-0.25, 1.25), xlim=(0, 8), )fig, ax = plt.subplots(figsize=(8, 3.6))plot_curves(ax)plt.show()
Figure 4.4.1: Share of voters expecting to vote for Dole at each point of the scale (dots, sized by the number of voters), with a straight line (dashed), which leaves the range 0 to 100%, and a logistic curve (red), which does not.
4.4.2 Into logit space, step by step
The way out is to change the scale of the probability until a straight line fits. Two steps do it. First, turn the probability p into odds, p / (1 − p): a 75% chance is odds of 3 to 1, or 3; a 20% chance is odds of 0.25. Odds run from 0 to infinity. Second, take the logarithm: the log-odds or logit run from minus infinity to plus infinity, and 50% sits at 0.
Now a straight line has room. Here is the whole journey with the voters:
Code
log_odds_of_shares = np.log(shares["mean"] / (1- shares["mean"]))line = logistic.params.Intercept + logistic.params.selfLR * SCALEdef step_answers(ax: Axes) ->None:"""Each answer as a dot at 0 or 1.""" x = anes.selfLR + rng.uniform(-0.25, 0.25, len(anes)) y = anes.dole + rng.uniform(-0.06, 0.06, len(anes)) ax.scatter(x, y, s=3, color=style.BLUE, alpha=0.25) ax.set( title="1. Answers", xlabel="self-placement", yticks=[0, 1], yticklabels=["Clinton", "Dole"], )def step_shares(ax: Axes) ->None:"""Shares voting Dole per step.""" ax.scatter(shares.index, shares["mean"], s=shares["count"] /3, color=style.BLUE, alpha=0.7) ax.yaxis.set_major_formatter(style.PERCENT) ax.set( title="2. Shares", xlabel="self-placement", ylabel="share for Dole", ylim=(-0.05, 1.05) )def step_logits(ax: Axes) ->None:"""Shares on the log-odds scale, with the fitted straight line.""" ax.scatter( shares.index, log_odds_of_shares, s=shares["count"] /3, color=style.GREEN, alpha=0.8 ) ax.plot(SCALE, line, color=style.CORAL) ax.axhline(0, color=style.INK, lw=0.8) ax.set( title="3. Log-odds: a line", xlabel="self-placement", ylabel="log-odds of Dole", xlim=(0.5, 7.5), )def step_back(ax: Axes) ->None:"""The fitted line on the probability scale.""" ax.scatter(shares.index, shares["mean"], s=shares["count"] /3, color=style.BLUE, alpha=0.7) ax.plot(SCALE, inverse_logit(line), color=style.CORAL) ax.yaxis.set_major_formatter(style.PERCENT) ax.set(title="4. Back: an S", xlabel="self-placement", ylim=(-0.05, 1.05), xlim=(0.5, 7.5))fig, axes = style.subplots(4, cols=4, height=2.8, width=11.5)for ax, step inzip(axes, [step_answers, step_shares, step_logits, step_back]): step(ax)fig.tight_layout()plt.show()
Figure 4.4.2: The road into logit space. 1: each voter’s answer, 0 or 1, against self-placement (jittered so the dots do not hide each other). 2: the share voting Dole at each step, between 0 and 1. 3: each share turned into log-odds, log(p / (1 − p)); on this scale the seven points lie close to a straight line, which the model fits. 4: the fitted line turned back into probabilities is the S-curve.
So logistic regression is a straight line on the log-odds scale:
The right-hand form is the logistic function, the S. The logit is called the link: it links the probability to the straight line (Chapter 4.6 meets other links). The model is not fitted to the seven shares, which only make the picture; it is fitted to all 944 answers at once (Section 4.4.4).
4.4.3 Odds ratios
On the log-odds scale every step of the predictor adds the same amount, . Adding to log-odds is multiplying the odds, so every step multiplies the odds by the same factor, : the odds ratio.
Code
odds_ratio = np.exp(logistic.params.selfLR)def plot_log_odds(ax: Axes) ->None:"""The fitted straight line on the log-odds scale.""" ax.plot(SCALE, line, color=style.CORAL) ax.axhline(0, color=style.INK, lw=0.8) ax.set_xticks(range(1, 8)) ax.set( title="Log-odds: a line", xlabel="self-placement", ylabel="log-odds of a Dole vote", xlim=(0, 8), ) ax.text(0.3, 0.3, "50% chance", fontsize=9, color=style.INK)def plot_probability(ax: Axes) ->None:"""The same line turned into probabilities.""" ax.plot(SCALE, inverse_logit(line), color=style.CORAL) ax.axhline(0.5, color=style.INK, lw=0.8) ax.yaxis.set_major_formatter(style.PERCENT) ax.set_xticks(range(1, 8)) style.note(ax, f"each step: odds × {odds_ratio:.1f}", x=0.03, color=style.INK) ax.set( title="Probability: an S", xlabel="self-placement", ylabel="chance of a Dole vote", xlim=(0, 8), )fig, axes = style.subplots(2, height=3.2, width=9.5)plot_log_odds(axes[0])plot_probability(axes[1])fig.tight_layout()plt.show()
Figure 4.4.3: The same model on two scales. Left: on the log-odds scale it is a straight line, rising by the same amount per step. Right: turned back into probabilities, the line becomes the S-curve. Each step multiplies the odds of a Dole vote by the same factor.
One step to the right on the scale multiplies the odds of a Dole vote by 3.3. On the probability scale the same step does different things: it moves a voter from 1% to 4% at the liberal end, but from 28% to 56% in the middle. The odds ratio is constant; the change in probability is not. That is why logistic results are reported as odds ratios, and why “the effect on the probability” needs a starting point. A useful rule of thumb: near 50%, a change of 1 in log-odds is a change of about 25 percentage points in probability (the slope of the S at its middle is ¼).
4.4.4 How the model is fitted
There is no sum of squares to minimise; the answers are 0 or 1. Instead, maximum likelihood (Chapter 2.3): for each voter, the model gives a probability of the answer Dole. The likelihood of all the answers together is the product of for Dole voters and for Clinton voters, the Bernoulli likelihood:
Each pair of coefficients gives a different S and a different log-likelihood. The fit is the top of that landscape:
Code
intercepts = np.linspace(-7.5, -3.5, 120)slopes = np.linspace(0.7, 1.7, 120)x_voters, y_voters = anes.selfLR.to_numpy(), anes.dole.to_numpy()def log_likelihood(b0: float, b1: float) ->float:"""Bernoulli log-likelihood of all answers for one intercept and slope.""" p = inverse_logit(b0 + b1 * x_voters)returnfloat(np.sum(y_voters * np.log(p) + (1- y_voters) * np.log(1- p)))surface = np.array([[log_likelihood(b0, b1) for b0 in intercepts] for b1 in slopes])OTHERS = [(-4.5, 0.9, style.GOLD), (-7, 1.6, style.GREEN)]def plot_surface(ax: Axes) ->None:"""Contours of the log-likelihood with the maximum and two other points.""" top = surface.max() ax.contour( intercepts, slopes, surface, levels=top - np.array([200, 100, 50, 20, 8, 2]), colors=style.INK, linewidths=0.8, linestyles="solid", ) ax.plot( logistic.params.Intercept, logistic.params.selfLR,"x", ms=10, mew=2.5, color=style.CORAL, )for b0, b1, color in OTHERS: ax.plot(b0, b1, "o", ms=7, color=color) ax.set(title="The likelihood hill", xlabel="intercept β₀", ylabel="slope β₁")def plot_candidates(ax: Axes) ->None:"""The S-curves of the fitted and the two other points.""" ax.scatter(shares.index, shares["mean"], s=shares["count"] /3, color=style.BLUE, alpha=0.6)for b0, b1, color in OTHERS + [ (logistic.params.Intercept, logistic.params.selfLR, style.CORAL) ]: ax.plot( SCALE, inverse_logit(b0 + b1 * SCALE), color=color, lw=2if color == style.CORAL else1.4, ) ax.yaxis.set_major_formatter(style.PERCENT) ax.set( title="Three S-curves", xlabel="self-placement", ylabel="chance of a Dole vote", xlim=(0.5, 7.5), )fig, axes = style.subplots(2, height=3.3, width=9.5)plot_surface(axes[0])plot_candidates(axes[1])fig.tight_layout()plt.show()
Figure 4.4.4: Left: the log-likelihood of the voters’ answers for every pair of intercept and slope; the contours are like a map of a hill, and the cross marks the top, the fitted model. Right: S-curves for three points on the map; the fitted one (red) follows the shares best.
There is no formula for the top, so the computer climbs: start somewhere, look at the slope of the hill, step uphill, repeat. (The usual method, iteratively reweighted least squares, solves a weighted regression at every step.) The standard errors come from how sharply the hill curves at the top, as in Chapter 2.3.
4.4.5 Tests and deviance
Two tests answer “is this coefficient 0?”:
The Wald test divides the estimate by its standard error, , and compares it with the normal curve. It is what the regression table shows.
The likelihood ratio test fits the model with and without the coefficient and compares the log-likelihoods (Chapter 2.3). It is more reliable, especially with small samples or very large effects.
The likelihood ratio test is usually written with the deviance, , the logistic version of the residual sum of squares: smaller is better. The null deviance belongs to the model with only an intercept (everyone has the same chance); the drop from it to the model’s deviance is the likelihood ratio statistic, compared with a chi-square with one degree of freedom per coefficient added.
Code
full = smf.logit("dole ~ selfLR + age + educ + income", anes).fit(disp=0)table = pd.DataFrame( {"log-odds": full.params,"standard error": full.bse,"z": full.tvalues,"odds ratio": np.exp(full.params),"95% from": np.exp(full.conf_int()[0]),"95% to": np.exp(full.conf_int()[1]), }).round(3)table = table.astype(object)table["p-value"] = ["< 0.001"if p <0.001elsef"{p:.3f}"for p in full.pvalues]table.index = ["intercept", "self-placement", "age", "education", "income"]print(table.to_markdown(colalign=["left"] + ["right"] *7))
Table 4.4.1: Logistic regression of a Dole vote on self-placement, age, education and income: coefficients on the log-odds scale, odds ratios with 95% intervals, and Wald p-values.
log-odds
standard error
z
odds ratio
95% from
95% to
p-value
intercept
-8.182
0.618
-13.242
0
0
0.001
< 0.001
self-placement
1.221
0.079
15.418
3.392
2.904
3.962
< 0.001
age
0.006
0.005
1.194
1.006
0.996
1.017
0.233
education
0.167
0.058
2.859
1.181
1.054
1.324
0.004
income
0.077
0.016
4.685
1.08
1.046
1.115
< 0.001
The two tests look at the same log-likelihood hill from two sides. Take the coefficient of age in the model below. The likelihood ratio test measures how far the hill drops between its top and the point where the age coefficient is held at 0 (with all other coefficients refitted). The Wald test instead fits a parabola to the top of the hill, using its curvature (the standard error), and reads the drop from the parabola:
Code
age_fit = smf.logit("dole ~ selfLR + age + educ + income", anes).fit(disp=0)age_hat, age_se = age_fit.params["age"], age_fit.bse["age"]others ="dole ~ selfLR + educ + income"def profile(value: float) ->float:"""Log-likelihood with the age coefficient fixed at `value`, the others refitted.""" model = smf.logit(others, anes, offset=value * anes.age)returnfloat(model.fit(disp=0).llf)grid = np.linspace(age_hat -4* age_se, age_hat +4* age_se, 41)curve = np.array([profile(v) for v in grid])top = age_fit.llflr_statistic =2* (top - profile(0.0))fig, axes = style.subplots(2, height=3.2, width=10.5)axes[0].plot(grid, curve, color=style.BLUE, lw=2, label="log-likelihood")axes[0].plot( grid, top -0.5* ((grid - age_hat) / age_se) **2, color=style.INK, ls="--", lw=1.2, label="Wald's parabola",)axes[0].axvline(0, color=style.CORAL, lw=1)axes[0].legend(fontsize=8)style.note( axes[0],f"Wald z² = {(age_hat / age_se) **2:.2f}\nLR = {lr_statistic:.2f}", x=0.03, y=0.97, color=style.INK,)axes[0].set(xlabel="age coefficient (log-odds per year)", ylabel="log-likelihood")style.null_panel( axes[1], lr_statistic, dist=stats.chi2(1), tail="right", xlabel="likelihood ratio statistic")axes[1].set_xlim(0, max(8, lr_statistic *1.3))fig.tight_layout()plt.show()
Figure 4.4.5: Testing the age coefficient in the model of Table 1.4.4.1. Left: the log-likelihood with the age coefficient fixed at each value and the other coefficients refitted (blue), and the parabola the Wald test assumes (dashed). Here the parabola fits the hill well, so the two tests agree. The drop from the top to the value at 0, doubled, is the likelihood ratio statistic. Right: that statistic against the chi-square distribution with 1 degree of freedom, its yardstick when the true coefficient is 0; the red tail is the p-value.
When the hill is lopsided, as with small samples or coefficients near separation (Section 4.4.7), the parabola misfits and the likelihood ratio test is the one to trust.
For this model the null deviance is 1,282 and the model’s deviance 853: a drop of 429 for four coefficients, far beyond a chi-square with 4 degrees of freedom. McFadden’s pseudo-R², , is 0.33. It is not a share of variance explained, and values of 0.2 to 0.4 already indicate a good fit.
4.4.6 Several predictors
As in Chapter 4.2, each coefficient is the effect of its predictor with the others held fixed, here on the log-odds scale. Odds ratios are easiest to compare in one picture, and predicted probabilities show what they mean for real people:
Code
LABELS = {"selfLR": "one step right on the scale","educ": "one step more education","income": "one step higher income","age": "one year older",}def plot_odds_ratios(ax: Axes, model: LikelihoodModelResults) ->None:"""Odds ratios with their intervals, on a log scale, against 1.""" intervals = np.exp(model.conf_int())for i, (name, label) inenumerate(LABELS.items()): low, high = intervals.loc[name] ax.plot([low, high], [i, i], color=style.CORAL, lw=2.5) ax.plot(np.exp(model.params[name]), i, "o", ms=7, color=style.CORAL) ax.axvline(1, color=style.INK, lw=1) ax.set_xscale("log") ticks = [0.8, 1, 1.5, 2, 3, 4] ax.set_xticks(ticks, [str(t) for t in ticks]) ax.minorticks_off() ax.set_yticks(range(len(LABELS)), list(LABELS.values())) ax.invert_yaxis() ax.set(xlabel="odds ratio (log scale)")def plot_profiles(ax: Axes) ->None:"""Predicted probability along the scale for two education levels."""for level, color, label in ( (1, style.GOLD, "1 to 8 grades"), (6, style.PURPLE, "degree or more"), ): grid = pd.DataFrame( {"selfLR": SCALE,"age": anes.age.mean(),"educ": level,"income": anes.income.mean(), } ) ax.plot(SCALE, full.predict(grid), color=color, lw=2, label=label) ax.yaxis.set_major_formatter(style.PERCENT) ax.legend(fontsize=8, title="education", title_fontsize=8) ax.set(xlabel="self-placement", ylabel="chance of a Dole vote", xlim=(1, 7))fig = plt.figure(figsize=(10.5, 3.2))grid = fig.add_gridspec(1, 2, width_ratios=[1.2, 1])plot_odds_ratios(fig.add_subplot(grid[0]), full)plot_profiles(fig.add_subplot(grid[1]))fig.tight_layout()plt.show()
Figure 4.4.6: Left: odds ratios for a Dole vote from the model with four predictors, with 95% intervals, on a log scale; 1 means no effect. Right: predicted chance of a Dole vote along the scale for voters of average age and income with little (1 to 8 grades) and much (a degree or more) education. The same odds ratio for education moves the probability most in the middle of the scale.
4.4.7 Checking a logistic model
Residuals of 0/1 answers are not much use one at a time: each is either or . The checks therefore group cases with similar predicted probabilities.
Figure 4.4.7: Left: calibration. The voters are grouped into ten bins by predicted chance of a Dole vote; each dot is a bin’s average prediction against the share who actually said Dole, with a 95% interval. A well-calibrated model lies on the diagonal. Right: binned residuals, the average residual per bin against the average prediction, with ±2 standard error bands (dashed); most bins should fall inside.
The model is well calibrated: when it says 70%, about 70% of those voters said Dole. Calibration plots also check classifiers (Chapter 6.1); the Hosmer-Lemeshow test turns the same binning into a chi-square test.
One failure is special to logistic regression. If a predictor separates the two outcomes perfectly, every voter above some value saying Dole and everyone below saying Clinton, the S can always be made steeper and fit better, so the slope runs off to infinity. This is separation. The software warns about enormous standard errors; the fix is a penalised fit (Chapter 6.4) or Firth’s correction.
4.4.8 More than two outcomes
With three or more unordered categories, multinomial logistic regression gives each category but one (the reference) its own line on the log-odds scale against the reference, and turns them into probabilities that add up to 1. For iris species from petal length:
Figure 4.4.8: Multinomial logistic regression of iris species on petal length. The three probabilities add up to 100% at every length; the dots at the top and bottom show the flowers of each species. Short petals mean setosa, long petals virginica, and versicolor wins in between. Petal length separates setosa perfectly, so the plain fit runs into separation; this one uses a mild penalty (Chapter 6.4).
When the categories are ordered (disagree, neutral, agree), an ordinal logistic model is simpler: one line on the log-odds scale and a set of cut points along it.
4.4.9 In one picture
A straight line cannot model a probability; its log-odds can: .
Back on the probability scale the line is an S: .
is the odds ratio: constant per step, while the change in probability depends on where you start.
The model is fitted by maximum likelihood on the Bernoulli likelihood.
Test with Wald z (the table) or the likelihood ratio (deviance), which is more reliable.
Check calibration and binned residuals; beware of separation.
More than two categories: multinomial (unordered) or ordinal logistic regression.