Logistic regression as a GLM, and how to compare models

A generalised linear model is a linear regression wearing two adapters: a link function that maps the linear predictor onto the scale of the mean, and a likelihood that matches the kind of data you have. For a yes/no outcome:

$$ y_i \sim \text{Bernoulli}(p_i), \qquad \text{logit}(p_i) = \log\frac{p_i}{1-p_i} = \alpha + x_i^\top \beta $$

Writing that down in PyMC takes four lines. Everything that is hard about it - and everything this notebook is about - comes from the link: coefficients live on the log-odds scale, where neither priors nor posteriors mean what your intuition says they mean.

#1 · Question and data

Many wells in Bangladesh are contaminated with natural arsenic. A research team measured the wells in one area, told every household its level, and encouraged those above the safety limit (0.5 in the units used here, hundreds of micrograms per litre) to switch to a safe well nearby. A few years later they came back and recorded who had switched.

What drives the decision to switch - and by how much, in units a public-health planner can use? The data are from Gelman & Hill (2007), who use them for exactly this lesson.

Only households with an unsafe well are in the data, which is why arsenic starts at 0.51. distance is the distance in metres to the nearest known safe well, education the years of schooling of the head of household, association whether anyone in the household takes part in community organisations.

A scatter plot of a 0/1 outcome is two stripes of dots, so we plot the switching rate in bins instead. Further to walk: less switching. More arsenic: more switching, but the curve flattens - remember that for later.

#2 · Scale the predictors before you think about priors

A prior is a statement about a coefficient, and a coefficient is "change in log-odds per unit of x". In metres, the distance coefficient would be around -0.006: what prior is "weakly informative" for that? Nobody knows without doing arithmetic. So change the units until a coefficient of 1 is a big but conceivable effect:

• dist100: distance in units of 100 m (most households are within 100 m of a safe well).
• educ4: education in units of 4 years.
• arsenic is already of order 1.

We also centre every predictor. Then $\alpha$ is the log-odds of switching for an average household - something we can reason about - instead of for a household with zero arsenic living on top of a safe well, which does not exist. Centring also makes main effects readable when we add interactions.

A model in this notebook is a list of terms; "a:b" means the product of two centred columns.

#3 · The model, and what priors look like through a logit link

One function builds every model we fit today. The two prior scales are arguments because we are about to compare two choices.

"When in doubt, make the prior wide" is a fine instinct on the scale of the outcome. On the logit scale it backfires. Let us simulate from a "vague" Normal(0, 10) on both parameters and from the defaults above, and look at what each claims about probabilities.

The "vague" prior is anything but. A log-odds of ±10 is a probability of 0.00005 or 0.99995, so Normal(0, 10) says: almost every household either certainly switches or certainly does not, and moving a few metres flips the village from one to the other. That is a strong, bizarre belief - the prior on the probability scale is a bathtub, with two thirds of its mass below 1% or above 99%.

The weakly-informative version spreads over the whole probability scale, allows steep curves in both directions, and only rules out cliffs. With 3020 observations either prior will be overwhelmed for the main effects. With few observations, rare outcomes, or perfectly separating predictors it will not, and the bathtub prior actively hurts. Always check GLM priors on the outcome scale.

#4 · Fit, diagnose, interpret: distance only

The helper below samples, optionally adds posterior predictive draws (for the calibration checks in section 7), and computes PSIS-LOO, which section 8 explains. LOO needs the pointwise log-likelihood, a 4000 x 3020 array of about 100 MB per model. We fit five models, and a notebook kernel never frees anything you still hold a name for, so the helper keeps the small LOO result and deletes the big array straight away. On a laptop this is the difference between a notebook that peaks above 3 GB and one that stays well under 2.

Diagnostics are clean (r_hat of 1.00, ESS in the thousands, no divergences), so we may interpret.

The intercept. Predictors are centred, so expit(alpha) is the switching probability for a household at the average distance of 48 m.

The slope, and the divide-by-4 rule. The logistic curve is steepest at $p = 0.5$, where its slope is $\beta/4$. So $\beta / 4$ is an upper bound on the change in probability per unit of $x$ - and a good approximation whenever probabilities are between roughly 0.25 and 0.75, as they are here.

#5 · Add arsenic, then ask a better question than "what is beta?"

Both effects have the expected sign, and the distance coefficient became more negative once arsenic entered (-0.62 before, -0.89 now). That happens when predictors are correlated in the "masking" direction: households far from a safe well tend to have somewhat higher arsenic (correlation 0.18, printed below), so in the distance-only model the two effects partly cancelled.

Which predictor matters more? Coefficients cannot say: they are per unit, and "1 unit of arsenic" and "100 m" are not comparable. They are also on the wrong scale. The planner's question is about probabilities, and because of the non-linear link the effect of distance on probability is different for every household.

The honest summary is an average predictive difference: take every household as it is, set the predictor to a low value, then to a high value, leave everything else alone, and average the change in predicted probability over households. Do it once per posterior draw and you get a posterior distribution for the answer.

Note the counterfactuals are set on the raw columns (distance in metres, arsenic) and pushed through prepare and design again, so derived columns such as the logarithm or an interaction stay consistent. Compare the first line with the divide-by-4 bound from the coefficient in the table above (-0.89 / 4 = -0.22): close, and a little smaller, as it must be.

#6 · An interaction

Does distance put people off less when the arsenic level is frightening? That is an interaction: the slope of one predictor depends on the value of another. With centred inputs the main effects keep their meaning (the effect of distance at average arsenic).

The interaction is probably negative, but the 94% interval is wide and just includes zero. Hold on to that thought: whether this term earns its place is a model comparison question, and we will get a quantitative answer in section 8.

#7 · Posterior predictive checks that work for binary data

The standard density overlay (az.plot_ppc_dist) is useless here. The data are zeros and ones; any model with an intercept reproduces the overall share of ones, so the check always passes. Two checks that can fail:

#Calibration

When the model says 70%, do 70% of those households switch? az.plot_ppc_pava estimates the conditional event probability as a function of the predicted probability (with isotonic regression, the pool-adjacent-violators algorithm - so no arbitrary bins) and compares it with the diagonal. The band shows where the curve could fall if the model were perfectly calibrated.

The curve tracks the diagonal through the middle, but at both ends it hugs the edge of the band. Look at the right: above predicted probabilities of about 0.85 it runs along the lower edge - the households the model is most sure about switch a little less often than promised. Suggestive rather than damning; we need a sharper tool.

#Binned residuals

Calibration plots show that something is off, not where. Raw residuals $y_i - \hat p_i$ of a binary model are unreadable (two stripes again). Average them within bins of a predictor and they become informative: if the model is right, about 95% of the bin averages should fall inside $\pm 2\sqrt{\bar p(1-\bar p)/n}$.

Against distance the residuals look like noise: four bins poke just outside the bounds where chance alone would give about two of 40, but their signs alternate and there is no trend. Against arsenic there is a pattern: the three lowest-arsenic bins are negative, the first of them by 27 percentage points (the model predicts more switching than happens among households barely over the limit), the middle is mostly positive, and five of the six highest bins are negative again. A straight line in arsenic rises too slowly at first and too fast later - the signature of diminishing returns, which is what the very first plot hinted at. The standard fix is to put arsenic on the log scale: going from 0.5 to 1 then matters as much as going from 2 to 4.

Better, not perfect. The wave in the residuals has gone and the calibration curve now stays inside its band at the top end, but the very first bin - households just above the safety limit - still switches far less than the model expects. (Gelman & Hill suspect a psychological threshold: "barely unsafe" does not feel unsafe. A smooth function of arsenic cannot capture that.)

One more step before comparing everything: education and community association, and whether education changes how people respond to distance.

#8 · Model comparison with PSIS-LOO

We have five models. Binned residuals told us how some of them fail; LOO tells us how much that matters for predicting households we have not seen.

The quantity being estimated is the expected log pointwise predictive density, elpd: the sum over observations of the log probability the model assigns to $y_i$ when $y_i$ was not used for fitting. Refitting 3020 times is out of the question, so PSIS-LOO re-weights the draws of the one posterior you have (Pareto-smoothed importance sampling). That is still the slow step of this notebook (a 4000 x 3020 matrix of weights to smooth per model), so the fit helper computed it once per model, as az.loo(idata, pointwise=True), and kept only the result; az.compare accepts those results in place of the DataTrees. pointwise=True keeps the per-household values, which we will want later. Start by reading a single model:

How to read this:

• elpd_loo is only meaningful relative to another model on the same data. Higher (less negative) is better. For scale: a coin flip would score $3020 \times \log 0.5 \approx -2093$.
• p_loo is the effective number of parameters. Here it is close to the actual count (six coefficients), as it should be for a well-specified model with weak priors. A p_loo far above the parameter count is a warning sign of misspecification.
• Pareto $k$ is a per-observation reliability diagnostic for the importance sampling itself. $k < 0.7$: fine. Larger: that observation is so influential that removing it changes the posterior more than re-weighting can mimic, and the elpd estimate cannot be trusted. With 3020 observations and six parameters no single household is influential, so every $k$ is tiny. (You will meet bad ones in the challenges - they usually point at outliers or at an over-flexible model.)

Reading the table:

• Rows are sorted best first. elpd_diff is each model's elpd minus the best model's; p is p_loo.
• dse is the standard error of that difference. It is not derived from the two se columns: the pointwise elpd values of two models on the same data are highly correlated (a household that is hard to predict is hard for both), so the difference is known far more precisely than either total. Always compare elpd_diff with dse, never with se.
• Rules of thumb: a difference under about 4 is too small to care about whatever its dse; beyond that, take it seriously when it is several dse from zero. p_worse turns the pair into a probability under a normal approximation - treat it as a rough guide.
• dse is always relative to the top row. To compare two models further down, call az.compare on just those two, as we do below.
• weight is a stacking weight: how to mix the models' predictions to predict best. It is not a posterior probability that the model is true. Near-duplicates share or lose weight arbitrarily.

What all this says:

1. Adding arsenic to distance is a huge improvement (about 70 elpd).
2. The log transformation beats raw arsenic by about 16 elpd, almost four times its dse, with the same number of parameters. LOO agrees with the binned residuals.
3. The distance x arsenic interaction buys nothing: half a unit of elpd, a fraction of its dse, for one extra parameter. A coefficient with a 96% posterior probability of being negative can still be irrelevant for prediction. Drop it, unless the interaction itself is the scientific question.
4. The full model is best: 15 elpd ahead of the runner-up, about 2.5 dse. Convincing, if not overwhelming.

LOO also works pointwise, which makes it a diagnostic rather than only a scoreboard. Where does the log model gain over the linear one?

Most of the gain - roughly 13 of the 16 units - comes from the two end bins: the barely-unsafe wells and the most contaminated ones, exactly where the binned residuals of the linear model were negative. In between, each bin moves by about one unit or less in either direction: there the two models predict almost equally well.

#9 · Answer the planner

Back to probabilities, with the best model. The interaction with education makes "the effect of distance" depend on who you ask, which is precisely the situation average predictive differences are built for.

For a planner:

• Having to walk 100 m instead of next door lowers the switching probability by about 22 percentage points on average. The deterrent is strongest among the least educated; with 12 years of schooling the curve is much flatter (and, with few such households far from a safe well, much less certain). New safe wells are worth most where the nearest alternative is far away and education is low.
• The choice of scale for arsenic was not a technicality. The linear model put the effect of going from 0.5 to 1.0 at under 6 percentage points; the log model says 15 - the same as going from 1.0 to 2.0. The better-fitting model tells a different story about who responds to information.
• The surprise is community association: if anything its members switch less, though the interval reaches zero.

A caution the model cannot give you: this is an observational survey. "Setting distance to 0" in a computation is not the same as drilling a well next door; calling a predictive difference causal takes assumptions that no amount of LOO can check.

#Try it yourself

1. Refit the distance-only model with the Normal(0, 10) priors on all 3020 households and then on a random subsample of 30. How much do the two priors disagree in each case? Plot the posterior curves for both.
2. The binned residuals still show too few switchers just above the safety limit. Add a term that lets the model treat "barely unsafe" wells differently (for example an indicator for arsenic < 0.6, or a spline in log arsenic) and compare with az.compare. Is the improvement worth a parameter nobody can interpret?
3. Swap the link: replace logit_p=eta with p=pm.math.invprobit(eta). Coefficients shrink by a factor of about 1.6 - do the average predictive differences change? Does LOO care?

Next: challenge C03 moves from binary outcomes to counts, and gives LOO's Pareto-$k$ diagnostic real work to do.