Diabetes Readmission

Hospitals in the US are penalised when patients come back within 30 days. Flagging high-risk diabetic patients at discharge would let a team call them, review their medication or schedule an early visit. The UCI Diabetes 130-US Hospitals dataset (101,766 encounters, 1999–2008) is the standard public benchmark for this — and it is a hard one.

The question here is not “which algorithm wins a leaderboard” but how far a well-specified logistic regression goes, what gradient boosting adds on top, and what each metric actually reveals when the target is rare and the signal is weak.

Almost half of the encounters belong to patients seen more than once. A random split puts the same person in training and test, which can let a flexible model recognise people instead of learning risk. So the held-out test set is 13,998 patients never seen in training, hyperparameters are tuned with 5-fold StratifiedGroupKFold, and confidence intervals come from a patient-level bootstrap.

Then I measured what the shortcut would have cost:

An honest null result. Splitting encounters instead of patients inflated CV AUC by at most 0.0016 — even for a forest with leaves of size 1, built to memorise. No feature identifies the patient, so there is little to memorise. The grouped split is still the right one (it is the only one that answers “how does this work on new patients?”), but here the cost of getting it wrong was small. Test AUC sits 0.006–0.009 above CV AUC for every model alike: a property of this particular split, not of any model.

Feature engineering is a row-wise function (age bands → years, ICD-9 codes → chapters, discharge IDs → destinations, medication counts) that sits at the start of every scikit-learn Pipeline. Everything learned from data — scaling, spline bases, rare-level pooling — is fitted inside each CV fold. Every model is tuned on the same grouped folds by log loss, a proper score, not accuracy.

The selected elastic net is pure L1 and keeps 142 of 201 coefficients — sparsity without loss. Note the class-weighted LR: same features, same ranking, twice the log loss. Section 06 explains why.

On held-out patients, LightGBM reaches ROC AUC 0.681 and the main logistic regression 0.670. Measured from chance (0.5), the logistic model keeps 94% of the lift that boosting achieves — with about 200 coefficients that can be read, audited and shipped as a formula. One feature alone (prior inpatient stays) already gets 0.611.

Can the logistic model close the gap? Boosting gets two things for free: high-cardinality codes and interactions. I gave the logistic regression both — the 3-digit ICD-9 codes (~700 levels, rare ones pooled) and interactions suggested by the trees (pairs of log-counts, prior utilisation × discharge destination):

The richest logistic model closes about a third of the gap (+0.004 of 0.011 AUC) and gains most in PR AUC and log loss, but LightGBM stays significantly ahead. The rest of the signal lives in higher-order interactions that a hand-specified linear predictor does not reach: on this dataset, interpretability costs about 0.007 AUC. Whether that is worth paying is a deployment decision, not a modelling one.

Threshold metrics depend on a decision, so the decision is made explicit. Three thresholds are chosen on out-of-fold training predictions — never on the test set — and the fourth is the default 0.5 that most tutorials use. Pick a policy:

Two lessons. First, at a sensible threshold every model except the one-feature baseline lands in a narrow band — F1 0.279–0.297, MCC 0.163–0.187: threshold metrics barely separate models, while proper scores and paired bootstrap differences do. Second, the accuracy paradox in plain sight: the best “accuracy” in the 0.5 tab belongs to models that flag almost nobody.

A risk score is only useful for decisions if “20%” means 20%. The main LR and LightGBM both track the diagonal (ECE 0.007). The class-weighted LR has the same AUC but predicts a mean risk of 47% for an outcome that happens 11% of the time:

Why it happens: weighting the positives by w shifts the logit intercept by about log w and leaves the slopes nearly untouched. Ranking is preserved, probabilities are not. Moving the threshold achieves the same trade-off without destroying the probabilities. A smaller lesson: the splines did not improve discrimination (0.6698 vs 0.6698) but they did improve calibration — slope 1.04 vs 1.08, maximum decile error 1.3 vs 2.7 points, Hosmer–Lemeshow p 0.077 vs < 0.001.

For interpretation I fitted a separate logistic GLM with readable terms (counts per unit, most frequent level as reference, 68 parameters) using IRLS written from scratch in NumPy — Newton–Raphson on the log-likelihood, with model-based, sandwich and cluster-robust (by patient) covariance. Coefficients and clustered standard errors match statsmodels to about 10⁻¹⁴. Terms with |z| > 3.5:

• Prior inpatient stays dominate: OR 1.30 per stay in the previous year (z = 35.6).
• Where the patient goes next is the second signal: rehab OR 3.7, another inpatient facility 2.0, SNF/ICF 1.5 (vs. home).
• Measured HbA1c > 8 is associated with lower risk (OR 0.86) — in line with Strack et al.: measuring HbA1c goes with closer follow-up.
• Clustering matters unevenly: the median clustered/model SE ratio is 1.02, but for prior ER visits it is 1.54 — the patients with many ER visits are exactly the ones who come back.

Associations, not causal effects. The predictive model (section 03) and this inferential model share the data, not the purpose.

A decision curve converts the threshold into a cost ratio: acting at τ means a missed readmission costs (1 − τ)/τ times a needless follow-up. Net benefit beats “treat everyone” from τ ≈ 0.05 upward but shrinks toward zero above ~0.25. The model pays for itself for cheap interventions — a phone call, a pharmacist review — not expensive ones.

And the data-size question: refitting both models on random subsets of patients shows where the extra flexibility starts to pay.

• Code and every number: GitHub repository ↗ — four scripts (download, experiment, close-the-gap, figures) write results/results.json; nothing on this page is typed by hand.
• Model: Hugging Face Hub ↗ — the main pipeline serialised with skops (no pickle), with a model card, metrics and a usage example.
• Tests: 23 pytest cases — ICD-9 grouping, cohort rules, IRLS vs scikit-learn and statsmodels (coefficients, model and clustered SEs), calibration metrics on simulated data, decision-curve formula.
• Limits: 1999–2008 US data with ICD-9 codes; readmissions to other hospitals are invisible, so the label under-counts; encounters ending in death or hospice were excluded and must not be scored. Not for clinical use.

Data: Strack et al. (2014), UCI Machine Learning Repository, doi:10.24432/C5230J, CC BY 4.0 — downloaded at runtime, not redistributed.