Thirty-eight cases, not an odds ratio

Turn the model back into probabilities

A logistic coefficient lives on the log-odds scale, which nobody thinks in. The fix is one line and it is the same line every time: predict, and average.

import pandas as pd
import numpy as np
import statsmodels.formula.api as smf

model = smf.logit("completed ~ disability + case_category + age_band + sex"
                  " + service_requested + admin1", data=d).fit(disp=0)

everyone_with = d.assign(disability=1)
everyone_without = d.assign(disability=0)

p1 = model.predict(everyone_with).mean()
p0 = model.predict(everyone_without).mean()
print(f"{p1:.4f} vs {p0:.4f}   average marginal effect {p1 - p0:+.4f}")

library(marginaleffects)

avg_comparisons(model, variables = "disability")

26.11% against 45.52%, a difference of −19.41 points.

That is the average marginal effect: take every case in the data, ask the model what it predicts if that case reported a disability, ask again if it did not, and average the difference. It is the model’s answer in the unit the question was asked in.

Compare it with the crude difference from the last lesson: −19.05 points. The adjusted odds ratio moved from 0.430 to 0.388 and the adjusted risk difference barely moved at all. The average marginal effect is the number that behaves the way a reader expects an adjusted estimate to behave.

Why one odds ratio is several risk differences

This is the property that makes the marginal effect necessary rather than merely convenient.

profile = dict(case_category="child-protection", age_band="0-11", sex="f",
               admin1="Artibonite")
rows = []
for service in ["health", "psychosocial", "legal", "livelihood-support"]:
    with_d = pd.DataFrame([{**profile, "service_requested": service, "disability": 1}])
    without = pd.DataFrame([{**profile, "service_requested": service, "disability": 0}])
    rows.append((service, model.predict(without)[0], model.predict(with_d)[0]))

print(pd.DataFrame(rows, columns=["service", "no disability", "disability"]).round(3))

predictions(model, newdata = datagrid(service_requested = unique(d$service_requested),
                                      disability = 0:1))

Service requested	No disability	Disability reported	Gap
Health	69.2%	46.5%	−22.7 pts
Psychosocial	64.5%	41.3%	−23.2 pts
Legal	41.0%	21.2%	−19.8 pts
Livelihood support	28.7%	13.5%	−15.2 pts

The odds ratio is 0.39 on every one of those rows. The model was fitted with a single disability coefficient, so by construction the odds ratio does not vary — and the risk difference varies from 15.2 points to 23.2 points anyway.

A constant odds ratio is not a constant effect. Where an outcome is already common, the same odds ratio moves more percentage points; where it is rare, fewer. So “the effect of disability” has no single answer in probability terms unless you say at what baseline — which is exactly what the average marginal effect does, by averaging over the baselines the caseload actually has.

An interval on the marginal effect

The marginal effect is a function of every coefficient, so its interval is not in the coefficient table. Two ways to get one, and the second is the one to reach for.

# statsmodels: delta method
margins = model.get_margeff(at="overall")
print(margins.summary())

avg_comparisons(model, variables = "disability")   # delta method, with CI

# bootstrap, when the delta method is awkward or the estimator is custom
rng = np.random.default_rng(20260729)
draws = []
for _ in range(400):
    sample = d.sample(len(d), replace=True, random_state=int(rng.integers(1e9)))
    m = smf.logit(model.model.formula, data=sample).fit(disp=0)
    draws.append(m.predict(sample.assign(disability=1)).mean()
                 - m.predict(sample.assign(disability=0)).mean())
print(np.percentile(draws, [2.5, 97.5]).round(4))

# 400 resamples is enough for a reportable interval; 2,000 for a published one.

−19.41 points, 95% CI −26.1 to −13.2 over 400 bootstrap resamples.

Set the seed and say how many resamples. A bootstrap interval that cannot be reproduced is not an interval, and this platform’s rule about generated numbers applies to the ones a model produces as much as to the ones a generator does.

Multiply it by the caseload

This is the sentence a programme manager reads, and the last course established why: an effect has to arrive in the unit the decision is made in.

n_disability = int((d["disability"] == 1).sum())
comparator = d[d["disability"] == 0]["completed"].mean()
observed = d[d["disability"] == 1]["completed"].sum()
print(f"{n_disability} cases reporting a disability")
print(f"expected at the comparator rate: {comparator * n_disability:.0f}")
print(f"observed: {int(observed)}  shortfall: {comparator * n_disability - observed:.0f}")

# Three lines, and it is the only line of the analysis a manager will quote.

197 cases reporting a disability. 90 would have completed at the comparator rate; 52 did. Thirty-eight cases short.

Thirty-eight is a number a protection team can act on: it is roughly three cases a month, it names a caseload rather than a percentage, and it can be compared against what a fix would cost. “OR 0.39” cannot do any of those things.

State the arithmetic, not just the result. A reader who can see 197 × 45.5% − 52 can check it; a reader given only “38 cases” has to trust it.

Four ways this goes wrong

Predicting at the mean instead of averaging the predictions. Setting every covariate to its mean and predicting once gives the effect for a case that does not exist — a case that is 0.6 female and 0.3 legal. Average the predictions over the real cases instead; that is what at="overall" and avg_comparisons do.

Reporting a marginal effect from a model with an interaction, without saying where. If the model lets the disability effect vary by service, the average is an average over a real difference and the table of profiles is the honest output.

Treating the shortfall as an effect. Thirty-eight cases is what the observed gap corresponds to, not what closing it would deliver. The comparison remains observational.

Extrapolating past the data. The model will happily predict a completion probability for an 80-year-old legal case in a department where none was recorded. Predict on profiles the data contains, and say which ones.

Report it whole

Referral completion by disability status, adjusted

  Logistic model, 1,581 consenting cases with complete covariates.
  Adjusted for case category, age band, sex, service requested, department.

  Average marginal effect   -19.4 points   95% CI -26.1 to -13.2
                                           (400 bootstrap resamples, seed 20260729)
  Predicted completion      26.1% with a disability reported
                            45.5% without

  Applied to the 197 cases reporting a disability, the gap is 38 cases that
  did not reach a service and would have at the comparator rate.

  The gap is not constant: it is 22.7 points for health referrals, where
  completion is common, and 15.2 points for livelihood support, where it is
  not. The odds ratio (0.39) is the same on both.

  Observational. Cases were not randomised, and the model adjusts only for
  what the register records.

The second-to-last paragraph is the one that stops the average being over-applied. A single number is what gets quoted; the sentence that says where it does and does not hold is what keeps the quoting honest.

What comes next

Every covariate so far has been chosen because it looked relevant. The next lesson adds one that is relevant, defensible and wrong — a variable on the causal path that removes 42% of the gap and improves every fit statistic while doing it.