Choose a model for hospital readmissions

What you will be able to do

Take a real prediction problem from the first look at the data to a model that a hospital team can use. You will compare four model families under one fair protocol, and decide with evidence whether the best score is worth its complexity.

Case studyPortfolio piece

Dataset

Diabetes 130-US hospitals, 1999–2008 (UCI)

Licence

CC BY 4.0

Time

About 2 hours

Runs in

Colab, Kaggle or VS Code

Skills used

• Data preparation
• Group split by patient
• Four model families
• Bootstrap confidence intervals
• Capacity threshold
• Fairness slices

Open in ColabOpen in KaggleView on GitHub

Colab and Kaggle open your own copy. To use VS Code, download the notebook from GitHub. Stuck? See the solution notebook.

The notebook runs every code block of this page, with five tasks and checks. Its setup cells download the data (3 MB) and check it, so it needs no local setup. The model cells take about a minute on a laptop with two CPU cores.

Before you startWhat to know and to have readyAbout 120 minutes

What you need to know

This case study uses what you learned in Modules 1 to 8. The ideas that it uses most:

• Train, validation and test sets, and cross-validation, from Train, validation and test sets.
• Average precision and precision at a capacity (precision@k), from The precision–recall trade-off.
• A pipeline with one-hot encoding and scaling, from Preprocess with a pipeline.
• The experiment log, from Machine learning is iterative: The experiment log.
• Decision trees, random forests and boosting, and what complexity costs, from Decision trees and Ensembles: The cost of complexity.
• Grid search and the one-standard-error rule, from Tune hyperparameters: The one-standard-error rule.
• Slices, and when a simple model is the better choice, from Slice analysis and When a simple model is better.
• A bootstrap interval, from Estimate uncertainty: The bootstrap in Essential Mathematics and Statistics for AI.

What you need on your computer (only for the local route; the notebook needs none of it)

• A Python environment with pandas 3.0.6 and scikit-learn 1.9.1, such as the .venv of your project folder from Define a baseline: Set up your project. This case study does not use the Larkfield files.
• A new, empty folder for this case study, for example ~/projects/readmissions, with the notebook in it. The notebook’s setup cell downloads the data into readmissions/data/. To start again, delete that folder and run the setup cell again.

Tested on macOS with Python 3.14.6, pandas 3.0.6 and scikit-learn 1.9.1, in Jupyter. Windows, Colab and Kaggle are not tested yet.

The situation

Nadia Haddad leads the care-transitions team of a hospital group. Her nurses phone patients in the first days after they leave hospital. On the call, a nurse checks the patient’s medicines, books a follow-up visit and asks who helps at home. Calls like these can stop a problem before the patient must come back.

The team cannot call everyone. With four nurses, Nadia can call about one discharged patient in ten. Today, the nurses choose with a rule of thumb: call the patients with the most hospital stays in the past year. Nadia asks you a simple question: “Can a model choose better whom we call?”

Felix Brandt, the quality director, adds his own question: “And if it can, which model? I hear that the newest methods win every competition. I also hear that we must be able to explain what we use. Show me the evidence.”

Each kind of error has a cost:

• A call to a patient who would not come back costs about 20 minutes of a nurse’s time. It does no harm to the patient.
• A patient who comes back and was not called is the expensive error. A readmission is bad for the patient, and it costs the hospital a bed and money.

Because the number of calls is fixed, the question is not “is this patient at risk?”. It is “who are the 10% most at risk?”. That is a ranking problem, like Larkfield’s capacity of 20%.

What you will build

You will build a model that ranks the patients of a discharge list by their risk of coming back to hospital within 30 days. The team calls the top 10%. On the way, you compare four model families that you know: logistic regression, a decision tree, a random forest and gradient boosting. You compare them under one protocol: the same split, the same folds, the same metric. You log every attempt. Then you decide, with a confidence interval and the one-standard-error rule, whether the best model is worth its complexity. You finish with slice checks, one final test, recommendations for Nadia and Felix, and a repository that you can show to an employer.

The data is real: 101,766 hospital stays of patients with diabetes, from 130 US hospitals between 1999 and 2008. Nadia’s hospital group and its team are made up for this case study. A real team would use its own, recent data, but every step would be the same.

The data

The data is the Diabetes 130-US hospitals for years 1999–2008 dataset of the UCI Machine Learning Repository (doi:10.24432/C5230J). Its authors are John Clore, Krzysztof Cios, Jon DeShazo and Beata Strack. It is published under the Creative Commons Attribution 4.0 licence, which allows reuse with credit and a note of the changes. The study that describes it is Strack et al. (2014), “Impact of HbA1c measurement on hospital readmission rates”, BioMed Research International (doi:10.1155/2014/781670).

Nextia Learning hosts a copy in the labs repository, with a dataset card that lists every change. There is one change: in two columns, the text None (“the test was not done”) is an empty cell in our copy. The source de-identified the data: it has no names and no dates. It is still health data about real people, so we never try to find out who anyone is.

The columns you will use

One row is one encounter: one hospital stay of one patient, from 1 to 14 days. The file has 50 columns. You will use these:

Column	What it means
encounter_id, patient_nbr	The number of the stay and of the patient. A patient can have several stays.
readmitted	<30: back in hospital within 30 days; >30: back later; NO: no record of a return.
age	Age band of 10 years, such as [70-80).
admission_type_id, admission_source_id, discharge_disposition_id	Codes for how the patient came in and where they went after the stay. A separate file translates them.
time_in_hospital	Days in hospital.
num_lab_procedures, num_procedures, num_medications, number_diagnoses	Counts during the stay.
number_outpatient, number_emergency, number_inpatient	Visits in the year before the stay: outpatient visits, emergency visits, hospital stays.
medical_specialty	The specialty of the admitting doctor.
diag_1	The main diagnosis, as an ICD-9 code (the international list of diseases that US hospitals used then).
A1Cresult, max_glu_serum	Results of two blood sugar tests, or empty if not measured.
insulin, change, diabetesMed	Changes to insulin, any change to diabetes medicines, and whether any diabetes medicine was given.
race, gender, payer_code	Sensitive columns: used only to check the model, never as inputs (see Frame the problem).

First checks

Load the file. The source writes a missing value as ?, so tell pandas to read ? as missing:

raw = pd.read_csv("data/diabetes_130_hospitals.csv.gz", na_values="?", low_memory=False)
print(raw.shape)
print("encounters:", raw["encounter_id"].nunique(), "| patients:", raw["patient_nbr"].nunique())
print("duplicate rows:", raw.duplicated().sum())

The output:

(101766, 50)
encounters: 101766 | patients: 71518
duplicate rows: 0

What this tells us. Every row is a different stay: 101,766 rows and 101,766 encounter numbers, and no duplicate rows. But there are only 71,518 patients, so about 30,000 rows are second, third or later stays of the same people. Keep this in mind. It decides how we split the data.

Next, the share of missing values in each column that has any:

missing = raw.isna().mean()
missing[missing > 0].sort_values(ascending=False).round(3)

The output:

weight               0.969
max_glu_serum        0.947
A1Cresult            0.833
medical_specialty    0.491
payer_code           0.396
race                 0.022
diag_3               0.014
diag_2               0.004
diag_1               0.000
dtype: float64

What this tells us. Three kinds of “missing” hide in this list, and each one needs a different decision:

• weight is missing in 97% of the stays. A column that is almost always empty cannot help, so we drop it.
• max_glu_serum and A1Cresult are not really missing. Empty means “the test was not done”, and that is information: Strack et al. studied exactly whether doing the HbA1c test changes readmissions. We make “not measured” a value of its own.
• medical_specialty is missing in half of the stays, probably because nobody filled it in. We keep “missing” as a value, so the model can learn whether it matters. diag_1 is missing in only 21 rows.

Now the target:

raw["readmitted"].value_counts()

The output:

readmitted
NO     54864
>30    35545
<30    11357
Name: count, dtype: int64

What this tells us. Only 11,357 stays (11%) end with a return within 30 days. Nadia’s calls happen in the first days after the stay, so a return after two months is not what she can prevent. We predict <30 and treat >30 and NO as “did not come back within 30 days”. The positive class is rare, so accuracy would mislead. A model that says “nobody comes back” is 89% right and useless, as Classification metrics showed.

Translate the codes

Three columns hold numbers, such as discharge_disposition_id = 11. Their meanings are in IDS_mapping.csv, which holds three small tables, one after the other, with an empty line between them. The next block reads them into one dictionary per column:

import csv

codes, column = {}, None
with open("data/IDS_mapping.csv", newline="") as file:
    for row in csv.reader(file):
        if not row or not row[0]:
            continue                          # the empty line between two tables
        if not row[0].isdigit():
            column = row[0]                   # a header line starts a new table
            codes[column] = {}
        else:
            codes[column][int(row[0])] = row[1].strip()
print({name: len(table) for name, table in codes.items()})
print(codes["admission_type_id"])

The output:

{'admission_type_id': 8, 'discharge_disposition_id': 30, 'admission_source_id': 25}
{1: 'Emergency', 2: 'Urgent', 3: 'Elective', 4: 'Newborn', 5: 'Not Available', 6: 'NULL', 7: 'Trauma Center', 8: 'Not Mapped'}

What this tells us. The codes have three tables of 8, 30 and 25 values. Look at the admission types: three of the eight codes (Not Available, NULL, Not Mapped) all mean “we do not know”. A model would treat them as three different things. We will join them into one value, unknown.

Remove the stays that cannot come back

Some stays end with the patient’s death, or with a move to a hospice for care at the end of life. These patients cannot be readmitted, and Nadia’s team would not call them. If we keep them, the model learns that “died in hospital” predicts “did not come back”, which is true and useless. Strack et al. removed them too.

ENDED = [11, 13, 14, 19, 20, 21]   # the patient died or went to a hospice
for code in ENDED:
    print(code, codes["discharge_disposition_id"][code])
data = raw[~raw["discharge_disposition_id"].isin(ENDED)].copy()
print(len(raw) - len(data), "encounters removed,", len(data), "left")

The output:

11 Expired
13 Hospice / home
14 Hospice / medical facility
19 Expired at home. Medicaid only, hospice.
20 Expired in a medical facility. Medicaid only, hospice.
21 Expired, place unknown. Medicaid only, hospice.
2423 encounters removed, 99343 left

What this tells us. 2,423 stays ended in death or hospice care. We keep 99,343. We print the six codes from the code table instead of trusting a list from memory, so anybody can check this step.

Now make the target, a 1 for a return within 30 days:

TARGET = "readmitted_30"
data[TARGET] = (data["readmitted"] == "<30").astype(int)
print(data[TARGET].sum(), "encounters were followed by a readmission within 30 days:",
      f"{data[TARGET].mean():.1%}")

The output:

11314 encounters were followed by a readmission within 30 days: 11.4%

What this tells us. The rate is 11.4%. If the team called 10% of patients at random, about 11 calls in 100 would reach a patient who comes back. Every model must do better than that.

Frame the problem

Before any model, write down what exactly we predict, and when. These choices protect the project from the most common mistakes.

• Unit: one encounter, at the moment the patient leaves hospital. A patient with three stays is three rows, and each stay gets its own score.
• Target: readmitted_30, a return as an inpatient within 30 days.
• Prediction time: the day of discharge. An input is allowed only if the team knows it on that day. Everything about the stay is known: its length, the tests, the medicines, and where the patient goes next. Visits in the year before the stay are known too.
• Not available at prediction time (leaks): readmitted itself, and anything that happens after discharge. A tempting example is “how many stays this patient has in the file”. The file also holds stays that come after this one, so that count uses the future. The failure case shows what it does.
• Not inputs at all: encounter_id and patient_nbr are numbers that identify a stay or a person, not facts about the patient’s health.

Sensitive columns

Three columns describe who a patient is, not how ill they are: race, gender and payer_code (the type of insurance, which also shows income). We keep them out of the inputs. A model that ranks patients by race or by insurance would give fewer calls to some groups for reasons that the team cannot defend. There is a real example of this kind of harm. In 2019, Obermeyer and colleagues published a study in Science about a widely used US health algorithm. It gave Black patients lower risk scores than white patients who were equally ill. The reason: it predicted health costs, and less money had been spent on Black patients.

Leaving these columns out is not enough on its own, because other inputs can carry the same information. So we keep them for slice checks: after training, we check whether the model works equally well, and calls a fair share, for each group. Age is different. It is a medical risk factor, and nurses consider it, so it is an input. We still check the model for each age band.

Prepare the features

The main diagnosis has hundreds of ICD-9 codes. Strack et al. grouped them into nine groups, such as circulatory (codes 390 to 459, and 785) and diabetes (250.xx). We use the same groups:

GROUPS = [("circulatory", 390, 460, "785"), ("respiratory", 460, 520, "786"),
          ("digestive", 520, 580, "787"), ("injury", 800, 1000, None),
          ("musculoskeletal", 710, 740, None), ("genitourinary", 580, 630, "788"),
          ("neoplasms", 140, 240, None)]


def diagnosis_group(code):
    """The group of the main diagnosis (an ICD-9 code), as in Strack et al. (2014)."""
    if pd.isna(code):
        return "missing"
    if code.startswith(("V", "E")):
        return "other"                        # visits for other reasons, and external causes
    if code.startswith("250"):
        return "diabetes"
    number = float(code)
    for name, low, high, extra in GROUPS:
        if low <= number < high or (extra and code.startswith(extra)):
            return name
    return "other"


print([diagnosis_group(code) for code in ["428", "250.83", "786", "V57", "996"]])

The output:

['circulatory', 'diabetes', 'respiratory', 'other', 'injury']

What this tells us. The function gives the expected group for five test codes. 428 (heart failure) is circulatory, 250.83 is diabetes, V57 (a rehabilitation visit) is other, and 996 (a complication of a device) is injury. Testing a function on a few known values takes one line and catches many mistakes.

Now translate the codes, join the “unknown” codes, and make the other features:

UNKNOWN = ["NULL", "Not Available", "Not Mapped", "Unknown/Invalid"]
for column, name in [("admission_type_id", "admission_type"), ("admission_source_id", "admission_source"),
                     ("discharge_disposition_id", "discharge_to")]:
    data[name] = data[column].map(codes[column]).replace(UNKNOWN, "unknown").fillna("unknown")
data["diagnosis_group"] = data["diag_1"].map(diagnosis_group)
data["age_years"] = data["age"].str.extract(r"(\d+)")[0].astype(int) + 5   # "[70-80)" becomes 75
data["specialty"] = data["medical_specialty"].fillna("missing")
tests = ["A1Cresult", "max_glu_serum"]
data[tests] = data[tests].fillna("not measured")
data[["admission_type", "discharge_to", "diagnosis_group", "age_years", "A1Cresult"]].head()

The output:

  admission_type        discharge_to diagnosis_group  age_years     A1Cresult
0        unknown             unknown        diabetes          5  not measured
1      Emergency  Discharged to home           other         15  not measured
2      Emergency  Discharged to home           other         25  not measured
3      Emergency  Discharged to home           other         35  not measured
4      Emergency  Discharged to home       neoplasms         45  not measured

What this tells us. The codes are now words that a nurse can read. The age band became a number (the middle of the band), so “older” means “larger”. The .fillna("unknown") catches codes that are not in the table at all.

Finally, write down the feature lists. The sensitive columns go in their own list, and the code checks that none of them is an input:

CATEGORIES = ["admission_type", "admission_source", "discharge_to", "specialty", "diagnosis_group",
              "A1Cresult", "max_glu_serum", "insulin", "change", "diabetesMed"]
NUMBERS = ["age_years", "time_in_hospital", "num_lab_procedures", "num_procedures", "num_medications",
           "number_outpatient", "number_emergency", "number_inpatient", "number_diagnoses"]
FEATURES = CATEGORIES + NUMBERS
SENSITIVE = ["race", "gender", "payer_code"]   # for slice checks only, never inputs
print(len(FEATURES), "features; sensitive columns among them:", set(SENSITIVE) & set(FEATURES))

The output:

19 features; sensitive columns among them: set()

What this tells us. The model will see 19 features, 10 categories and 9 numbers. set() is an empty set: no sensitive column is an input. To keep the case study small, we leave out the second and third diagnoses and 22 single medicine columns. The Go further list in Make it a portfolio piece suggests adding them.

Explore what the decision needs

Nadia’s rule of thumb uses one fact: hospital stays in the past year. Before we build anything, check how strong that fact is. A rule that already works well sets a high bar for a model.

prior_stays = data["number_inpatient"].clip(upper=5)   # 5 means "5 or more"
data.groupby(prior_stays)[TARGET].agg(["size", "mean"]).round(3)

The output:

                   size   mean
number_inpatient
0                 66245  0.086
1                 18984  0.133
2                  7300  0.179
3                  3271  0.210
4                  1574  0.241
5                  1969  0.371

The share of stays followed by a readmission within 30 days, by the number of hospital stays in the year before. The rate rises with every earlier stay.

Text description of the chart

A column chart. The horizontal axis is the number of hospital stays in the year before: 0, 1, 2, 3, 4, and 5 or more. The vertical axis is the share readmitted within 30 days: 8.6% (66,245 stays), 13.3% (18,984), 17.9% (7,300), 21.0% (3,271), 24.1% (1,574) and 37.1% (1,969).

What this tells us. The rule of thumb is a good one: the rate climbs from 8.6% with no earlier stay to 37.1% with five or more. But look at the sizes. Two stays in three (66,245) had no earlier stay at all, and their rate is 8.6%, not zero. Those patients make up a large share of all returns, and the rule cannot tell them apart. A model can only beat the rule if other facts help there.

Where the patient goes after the stay is the second thing a nurse knows at discharge:

by_discharge = data.groupby("discharge_to")[TARGET].agg(["size", "mean"])
by_discharge[by_discharge["size"] >= 1000].sort_values("mean", ascending=False).round(3)

The output:

                                                     size   mean
discharge_to
Discharged/transferred to another rehab fac inc...   1993  0.277
Discharged/transferred to another type of inpat...   1184  0.209
Discharged/transferred to another short term ho...   2128  0.161
Discharged/transferred to SNF                       13954  0.147
Discharged/transferred to home with home health...  12902  0.127
unknown                                              4680  0.118
Discharged to home                                  60234  0.093

What this tells us. Patients who go to a rehabilitation unit come back at 27.7%, three times as often as those who go home (9.3%). Patients who go to a skilled nursing facility (SNF) or home with home health care are in between. The rule does not use this fact, so a model that combines both should rank better than the rule. Whether it does, and by how much, is what the rest of the case study measures.

Split by patient

Why not split by row

In this data, 16,341 of the 69,990 patients have more than one stay:

encounters = data["patient_nbr"].value_counts()
print(encounters.size, "patients;", (encounters > 1).sum(), "of them have more than one encounter")
print("their encounters:", encounters[encounters > 1].sum(), "of", len(data), "rows")
print("most encounters for one patient:", encounters.max())

The output:

69990 patients; 16341 of them have more than one encounter
their encounters: 45694 of 99343 rows
most encounters for one patient: 40

What this tells us. Almost half of the rows (45,694) belong to patients who appear more than once, and one patient has 40 stays. The rows are not independent: two stays of the same person share that person’s illness, home and habits.

Think of an exam again. If the exam questions are the practice questions with new numbers, a student who memorised the answers looks better than they are. A row split does the same: it can put a patient’s March stay in train and their May stay in test. The model has then “seen” that patient, and the test score answers an easier question than Nadia’s. Her model will meet many patients for the first time. So we split by patient: every patient goes, with all their stays, into exactly one part.

The split

The split gives each patient, not each row, a part: 60% train, 20% valid, 20% test. A random generator with a fixed seed makes it repeatable.

patients = data["patient_nbr"].unique()
rng = np.random.default_rng(0)
part_of = dict(zip(patients, rng.choice(["train", "valid", "test"], size=len(patients), p=[0.6, 0.2, 0.2])))
data["part"] = data["patient_nbr"].map(part_of)
train, valid, test = (data[data["part"] == name].copy() for name in ["train", "valid", "test"])
for name, part in [("train", train), ("valid", valid), ("test", test)]:
    print(f"{name}: {len(part)} encounters, {part['patient_nbr'].nunique()} patients,"
          f" readmitted {part[TARGET].mean():.3f}")
print("patients in more than one part:", (data.groupby("patient_nbr")["part"].nunique() > 1).sum())

The output:

train: 59556 encounters, 41999 patients, readmitted 0.114
valid: 20179 encounters, 14139 patients, readmitted 0.113
test: 19608 encounters, 13852 patients, readmitted 0.116
patients in more than one part: 0

What this tells us. The three parts have about 60%, 20% and 20% of the rows, and almost the same readmission rate (11.3% to 11.6%). No patient is in two parts. From now on, the test set is sealed: we do not look at it until the final test.

Why not a split by time, as at Larkfield? A real deployment always predicts the future, so a time split would be best. But this file has no dates. The encounter numbers probably grow with time, but the source does not say so, and we do not build on a guess. We name this as a limit.

What a row split would do

To see the problem, split the train set by row, as train_test_split does by default (the test set stays sealed):

from sklearn.model_selection import train_test_split

fit_rows, check_rows = train_test_split(train, test_size=0.2, random_state=0)
seen = check_rows["patient_nbr"].isin(fit_rows["patient_nbr"])
print("patients in both parts:", check_rows.loc[seen, "patient_nbr"].nunique())
print(f"check rows of a patient that the model has seen: {seen.sum()} of {len(check_rows)} ({seen.mean():.0%})")

The output:

patients in both parts: 4046
check rows of a patient that the model has seen: 4915 of 11912 (41%)

What this tells us. With a row split, 41% of the rows that we “check” belong to 4,046 patients that the model trained on. Any score on this check set mixes two questions: “how well does the model know these patients?” and “how well does it predict new ones?”. Only the second one matters to Nadia. Later, you will measure how much this changes the score on this data.

Baselines and the protocol

One protocol for every model

“Compare fairly” means: the same data, the same folds, the same metric, the same features. We fix all four now, before any model, so that no model gets an advantage:

• Data: the train set only. Valid confirms the choice; test stays sealed.
• Folds: 5 folds, by patient. scikit-learn’s StratifiedGroupKFold does this. “Group” means that all rows of a patient go to the same fold, and “stratified” means that each fold gets about the same readmission rate.
• Metric: average precision (AP), as in The precision–recall trade-off. It summarises how well the model ranks the rare positives at the top, over all thresholds. It is the right metric to choose a model. For the decision, we add precision in the top 10% later.
• Features and preparation: the 19 features, in one pipeline for every model.

The protocol of this case study. Every model is chosen on the same 5 patient folds of train, confirmed on valid, and only the chosen model meets the test set, once.

Text description of the diagram

The diagram is a flow from left to right (top to bottom on a phone). The 99,343 encounters of 69,990 patients go to Split by patient. The split makes three parts. Train (59,556 encounters) goes to 5 folds by patient. There, we tune each model family and compare the families with the one-standard-error rule. Valid (20,179 encounters) and the result of the folds go to Confirm on valid. There, we compute a bootstrap interval for the gap between models, and the threshold for the 10% capacity. Test (19,608 encounters) and the confirmed model go to Test once.

The next block makes the folds and the pipeline. OneHotEncoder turns each category into 0/1 columns, as in Preprocess with a pipeline. Two settings are new. min_frequency=100 joins categories with fewer than 100 rows into one “infrequent” column, because a rare specialty with 12 stays only teaches noise. handle_unknown="infrequent_if_exist" sends a category that appears only in new data to that same column, instead of failing.

from sklearn.compose import ColumnTransformer
from sklearn.model_selection import StratifiedGroupKFold
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import OneHotEncoder, StandardScaler

folds = list(StratifiedGroupKFold(5, shuffle=True, random_state=0)
             .split(train, train[TARGET], groups=train["patient_nbr"]))


def build_pipeline(model, numbers=NUMBERS):
    """The same preparation for every model: one-hot categories (rare ones grouped) and scaled numbers."""
    prepare = ColumnTransformer([
        ("categories", OneHotEncoder(handle_unknown="infrequent_if_exist", min_frequency=100,
                                     sparse_output=False), CATEGORIES),
        ("numbers", StandardScaler(), numbers)])
    return Pipeline([("prepare", prepare), ("model", model)])

Check that the folds keep patients apart:

for number, (fit, check) in enumerate(folds, 1):
    shared = set(train["patient_nbr"].iloc[fit]) & set(train["patient_nbr"].iloc[check])
    print(f"fold {number}: {len(check)} check rows, readmitted {train[TARGET].iloc[check].mean():.3f},"
          f" patients also in the fit rows: {len(shared)}")

The output:

fold 1: 11912 check rows, readmitted 0.114, patients also in the fit rows: 0
fold 2: 11910 check rows, readmitted 0.114, patients also in the fit rows: 0
fold 3: 11912 check rows, readmitted 0.113, patients also in the fit rows: 0
fold 4: 11911 check rows, readmitted 0.114, patients also in the fit rows: 0
fold 5: 11911 check rows, readmitted 0.114, patients also in the fit rows: 0

What this tells us. Each fold checks about 11,900 rows with the same rate (11.3% to 11.4%), and no patient is on both sides. Every model below uses exactly these folds. The list is fixed, so all models are scored on the same rows.

Baselines

The experiment log is the table from Machine learning is iterative. It has this course’s columns: id, date, change, CV AP and its standard deviation, valid AP, decision and reason. The notebook’s helper log_experiment() writes a row to experiment_log.csv. It also keeps the five fold scores in FOLD_SCORES, so you can compare models fold by fold.

Two baselines come first. The floor gives every patient the same score: it ranks nobody, so its AP equals the readmission rate. The prior-stays rule is Nadia’s rule of thumb, written as a score. The score is the number of stays in the past year; emergency visits break ties. A rule needs no training, so we only measure it on each check fold.

from sklearn.metrics import average_precision_score


def stays_rule(df):
    """The rule of thumb: most hospital stays in the past year first; emergency visits break ties."""
    return df["number_inpatient"] + df["number_emergency"].clip(upper=99) / 100


def rule_folds(score):
    """The AP of a fixed score (nothing to train) on each of the 5 check folds."""
    return np.array([average_precision_score(train[TARGET].iloc[check], score.iloc[check])
                     for _, check in folds])


log_experiment("E01", "Baseline: the same score for every patient", rule_folds(train["number_inpatient"] * 0.0),
               decision="baseline", reason="The floor: its AP is the readmission rate.")
log_experiment("E02", "Baseline: the prior-stays rule", rule_folds(stays_rule(train)),
               decision="baseline", reason="The rule of thumb that a model must beat.")
show_log()

The output:

 id                                     change  cv_ap  cv_std valid_ap decision
E01 Baseline: the same score for every patient  0.114   0.000          baseline
E02             Baseline: the prior-stays rule  0.182   0.009          baseline

What this tells us. The floor is 0.114, the readmission rate. The rule reaches 0.182, which is 60% above the floor. That is the bar: a model that cannot beat 0.182 is not worth building.

Four model families under one protocol

Nobody can say in advance which family fits this data best, as Machine learning is iterative explained. Each family assumes a different shape of pattern, so we try four, from the simplest to the most complex:

• Logistic regression assumes that each fact adds or removes some risk on its own. It is fast and small, and each score splits into exact contributions.
• A decision tree asks a few questions in order. It can find combinations (“many stays and sent to rehabilitation”), and a small tree can be printed for the nurses. But one tree is coarse and unstable.
• A random forest averages many deep trees, each trained on a random sample of rows and features. Averaging lowers the variance of single trees.
• Gradient boosting adds small trees one after another; each one corrects the errors of the ones before. It often ranks best on table data.

Tune each family modestly

Each family has one or two settings that matter most. We tune one setting per family with a small grid, on the same folds, as in Tune hyperparameters. We keep the grids small on purpose. Every extra setting that we try on the same folds adds a little chance to fit the folds’ noise. The helper tune() applies the one-standard-error rule inside each family. Among the settings whose mean AP is within one standard error of the best, it keeps the simplest. Each grid lists its simplest setting first.

from sklearn.model_selection import GridSearchCV


def tune(pipeline, grid):
    """Score each setting of `grid` on the 5 patient folds. Keep the simplest setting whose mean is
    within one standard error of the best (each grid lists the simplest setting first)."""
    search = GridSearchCV(pipeline, grid, cv=folds, scoring="average_precision", refit=False)
    results = pd.DataFrame(search.fit(train[FEATURES], train[TARGET]).cv_results_)
    results["se"] = results["std_test_score"] / np.sqrt(len(folds))
    best = results["mean_test_score"].idxmax()
    bar = results.loc[best, "mean_test_score"] - results.loc[best, "se"]
    chosen = results.index[results["mean_test_score"] >= bar][0]
    print(results[["params", "mean_test_score", "std_test_score"]].round(4).to_string())
    print("one-SE bar:", round(bar, 4), "| chosen:", results.loc[chosen, "params"])
    pipeline.set_params(**results.loc[chosen, "params"])
    return np.array([results.loc[chosen, f"split{i}_test_score"] for i in range(len(folds))])

Logistic regression

The setting is C, the strength of the penalty on large coefficients: a small C means a strong penalty and a simpler model.

from sklearn.linear_model import LogisticRegression

logistic = build_pipeline(LogisticRegression(max_iter=1000))
scores = tune(logistic, {"model__C": [0.01, 0.1, 1.0]})
log_experiment("E03", "Logistic regression; C tuned over 0.01, 0.1, 1", scores, decision="keep",
               reason="Well above the rule; the strongest penalty is within one SE of the best.")

The output:

               params  mean_test_score  std_test_score
0  {'model__C': 0.01}           0.2100          0.0101
1   {'model__C': 0.1}           0.2116          0.0101
2   {'model__C': 1.0}           0.2114          0.0098
one-SE bar: 0.2071 | chosen: {'model__C': 0.01}

What this tells us. All three settings score between 0.2100 and 0.2116, a difference of 0.0016, far less than one standard error (about 0.0045). So the rule keeps the simplest, C = 0.01. With 60,000 rows and only 19 features, the penalty hardly matters. Logistic regression already beats the rule clearly: 0.210 against 0.182.

Decision tree

The setting is the maximum depth, with at least 100 stays in each leaf, so that every leaf’s rate rests on enough patients.

from sklearn.tree import DecisionTreeClassifier

tree = build_pipeline(DecisionTreeClassifier(min_samples_leaf=100, random_state=0))
scores = tune(tree, {"model__max_depth": [4, 6, 8, 10]})
log_experiment("E04", "Decision tree, min leaf 100; depth tuned over 4, 6, 8, 10", scores, decision="keep",
               reason="Above the rule, below logistic regression.")

The output:

                     params  mean_test_score  std_test_score
0   {'model__max_depth': 4}           0.1946          0.0071
1   {'model__max_depth': 6}           0.2020          0.0069
2   {'model__max_depth': 8}           0.2020          0.0061
3  {'model__max_depth': 10}           0.2002          0.0066
one-SE bar: 0.1989 | chosen: {'model__max_depth': 6}

What this tells us. Depth 6 and depth 8 tie at 0.2020, and depth 10 already gets worse: deeper trees start to fit noise. The rule keeps depth 6. Even the best tree (0.202) is below logistic regression (0.210). A tree of depth 6 has up to 64 leaves, already too many for a nurse to read. It also cuts every number into a few steps.

Random forest

We limit the forest to 100 trees so that the notebook runs in a few minutes; more trees change the score very little. The setting is the minimum number of stays in a leaf: larger leaves make simpler trees.

from sklearn.ensemble import RandomForestClassifier

forest = build_pipeline(RandomForestClassifier(n_estimators=100, n_jobs=-1, random_state=0))
scores = tune(forest, {"model__min_samples_leaf": [50, 20, 10]})
log_experiment("E05", "Random forest, 100 trees; min leaf tuned over 50, 20, 10", scores, decision="keep",
               reason="Above logistic regression on the folds.")

The output:

                            params  mean_test_score  std_test_score
0  {'model__min_samples_leaf': 50}           0.2163          0.0080
1  {'model__min_samples_leaf': 20}           0.2206          0.0099
2  {'model__min_samples_leaf': 10}           0.2237          0.0108
one-SE bar: 0.2188 | chosen: {'model__min_samples_leaf': 20}

What this tells us. Smaller leaves score higher: 0.2163, 0.2206, 0.2237. The best is a leaf of 10, but a leaf of 20 is within one standard error, so the rule keeps it. At 0.2206, the forest is the first model above logistic regression.

Gradient boosting

We use scikit-learn’s HistGradientBoostingClassifier with 100 trees. early_stopping=False makes it always build all 100, so every fit does the same amount of work. The setting is the depth of each small tree.

from sklearn.ensemble import HistGradientBoostingClassifier

boosting = build_pipeline(HistGradientBoostingClassifier(early_stopping=False, random_state=0))
scores = tune(boosting, {"model__max_depth": [2, 3, None]})
log_experiment("E06", "Gradient boosting, 100 trees; depth tuned over 2, 3, no limit", scores, decision="keep",
               reason="The best family on the folds, with the forest.")

The output:

                       params  mean_test_score  std_test_score
0     {'model__max_depth': 2}           0.2188          0.0085
1     {'model__max_depth': 3}           0.2226          0.0087
2  {'model__max_depth': None}           0.2212          0.0119
one-SE bar: 0.2188 | chosen: {'model__max_depth': 2}

What this tells us. Depth 3 is the best (0.2226), and depth 2 just reaches the bar (0.2188), so the rule keeps trees of depth 2. Each of these trees asks only two questions. Boosting is also above logistic regression.

Compare the families

Now put the chosen setting of every family side by side, with the standard error of each mean:

compare = pd.DataFrame({id: FOLD_SCORES[id] for id in ["E02", "E03", "E04", "E05", "E06"]},
                       index=[f"fold {i}" for i in range(1, 6)])
summary = pd.DataFrame({"mean": compare.mean(), "se": compare.std(ddof=0) / np.sqrt(len(folds))})
best = summary["mean"].idxmax()
bar = summary.loc[best, "mean"] - summary.loc[best, "se"]
print(summary.round(4).T)
print("best:", best, "| one-SE bar:", round(bar, 4), "| within one SE:", summary.index[summary["mean"] >= bar].tolist())
print("boosting minus logistic, per fold:", (compare["E06"] - compare["E03"]).round(4).tolist())

The output:

         E02     E03     E04     E05     E06
mean  0.1824  0.2100  0.2020  0.2206  0.2188
se    0.0040  0.0045  0.0031  0.0044  0.0038
best: E05 | one-SE bar: 0.2162 | within one SE: ['E05', 'E06']
boosting minus logistic, per fold: [0.0127, 0.007, 0.0136, 0.0065, 0.004]

CV average precision of each experiment on the same 5 patient folds of train, with one standard error on each side. The two ensembles (orange) reach the one-SE bar; logistic regression and the tree are below it.

Text description of the chart

An interval plot. For each experiment, a dot shows the mean CV average precision and a bar shows one standard error on each side. E01, the same score for all: 0.114. E02, the prior-stays rule: 0.182. E03, logistic regression: 0.210. E04, the decision tree: 0.202. E05, the random forest: 0.221. E06, gradient boosting: 0.219. A vertical line marks the one-SE bar (“bar”) at 0.216: only the forest and boosting reach it.

What this tells us. Read it with the one-standard-error rule: keep the simplest model whose mean is within one standard error of the best.

1. The best mean is the forest (0.2206). Its standard error is 0.0044, so the bar is 0.2162. Only the forest and boosting reach it. They tie.
2. Between the two ensembles, keep the simpler one: boosting. It has 100 trees of depth 2, and its file is about 190 times smaller (we measure this below). It is the challenger.
3. Logistic regression is below the bar. Its mean is 0.2100, a gap of 0.0088 to boosting, about two standard errors. And the gap is not the luck of one fold: boosting is higher on all five folds, by 0.004 to 0.014. Here, unlike at Larkfield in When a simple model is better, the complex model wins by more than the noise.

So the rule does not let us call it a tie. Boosting ranks better. One question remains, from Ensembles: The cost of complexity. Is the gain large enough, on the decision that matters, to pay for the complexity?

Two experiments that you must reject

Before the validation set, two more rows go in the log. Both look like improvements, and both are wrong.

E07: folds by row. How much would a row split change the score? Score the same boosting pipeline with ordinary 5-fold CV by row:

from sklearn.model_selection import StratifiedKFold, cross_val_score

row_folds = StratifiedKFold(5, shuffle=True, random_state=0)
row_scores = cross_val_score(boosting, train[FEATURES], train[TARGET], cv=row_folds, scoring="average_precision")
print("boosting, folds by row:    ", row_scores.round(4), round(row_scores.mean(), 4))
print("boosting, folds by patient:", FOLD_SCORES["E06"].round(4), round(FOLD_SCORES["E06"].mean(), 4))
log_experiment("E07", "E06 scored with folds by row (patients in both parts)", row_scores, decision="reject",
               reason="Not a valid protocol: it scores patients that the model has seen.")

The output:

boosting, folds by row:     [0.2087 0.2179 0.2315 0.2171 0.2211] 0.2193
boosting, folds by patient: [0.2191 0.2111 0.2186 0.2109 0.2343] 0.2188

What this tells us. On this data, the row split changes the mean by only 0.0005, much less than the noise. That is an honest result, and it has a reason. Our inputs describe the stay: its length, its tests, the visits before it. They change from one stay to the next, so the model cannot recognise a patient. Suppose an input is fixed for a person, or the data has a rich history of each patient. Then the same mistake gives a large, false gain. The split by patient is right because of what it measures, not because of the size of the change here. A result from a row split answers a question that Nadia did not ask.

E08 adds a feature that looks very useful. It is the failure case at the end of this page, and its row in the log says “reject: a leak”.

Metric and operating point

From average precision to calls

AP chose the families, because it summarises the whole ranking. Nadia’s decision uses only the top of the ranking: the 10% that the team calls. So the metric for the decision is precision@10%: of the patients called, the share who would come back. You met it as precision@k in The precision–recall trade-off. At a fixed capacity, the number of calls is fixed, so precision and recall rise and fall together. Precision is easier for Nadia to use: “1 call in 4 reaches a patient who would come back”.

Now, and only now, fit each tuned pipeline on the whole train set, and score the validation set:

from sklearn.metrics import roc_auc_score


def precision_at(y, score, share=0.10):
    """The precision among the `share` of patients with the highest scores: the patients the team calls."""
    y, score = np.asarray(y), np.asarray(score)
    top = np.argsort(-score, kind="stable")[: int(share * len(y))]
    return y[top].mean()


models = {"E03 logistic": logistic, "E04 tree": tree, "E05 forest": forest, "E06 boosting": boosting}
valid_scores = {"E02 rule": stays_rule(valid).to_numpy()}
for name, pipeline in models.items():
    valid_scores[name] = pipeline.fit(train[FEATURES], train[TARGET]).predict_proba(valid[FEATURES])[:, 1]
table = pd.DataFrame({name: {"AP": average_precision_score(valid[TARGET], score),
                             "ROC AUC": roc_auc_score(valid[TARGET], score),
                             "precision@10%": precision_at(valid[TARGET], score)}
                      for name, score in valid_scores.items()}).T
table.round(3)

The output:

                 AP  ROC AUC  precision@10%
E02 rule      0.175    0.614          0.240
E03 logistic  0.206    0.671          0.261
E04 tree      0.195    0.655          0.251
E05 forest    0.219    0.678          0.258
E06 boosting  0.214    0.673          0.268

Precision@10% on the validation set: the share of the called patients who come back within 30 days. Random calls reach 0.113, the rule 0.240, and the four models 0.251 to 0.268.

Text description of the chart

A column chart of precision@10% on the 20,179 validation stays. Random calls: 0.113. The prior-stays rule: 0.240. Logistic regression: 0.261. Decision tree: 0.251. Random forest: 0.258. Gradient boosting: 0.268.

What this tells us. The order of the families is the same as in CV. The ensembles rank best (AP 0.219 and 0.214), then logistic regression (0.206), the tree (0.195) and the rule (0.175). The ROC AUC tells the same story (0.614 to 0.678). But at the top 10%, the differences shrink. All four models reach 0.251 to 0.268, and even the rule reaches 0.240. The forest has the best AP and only the third-best precision@10%. A difference in AP does not always reach the decision.

Turn the precision into patients. The validation set stands for a period with 20,179 discharges, so the team makes 2,017 calls:

calls = int(0.10 * len(valid))
for name in ["E02 rule", "E03 logistic", "E06 boosting"]:
    reached = round(precision_at(valid[TARGET], valid_scores[name]) * calls)
    print(f"{name:13s} {calls} calls reach {reached} of the {valid[TARGET].sum()} patients who return")

The output:

E02 rule      2017 calls reach 484 of the 2273 patients who return
E03 logistic  2017 calls reach 526 of the 2273 patients who return
E06 boosting  2017 calls reach 540 of the 2273 patients who return

What this tells us. For the same 2,017 calls, logistic regression reaches 42 more patients who would come back than the rule (526 against 484). Boosting reaches 14 more than logistic regression (540 against 526). Are these differences real, or the luck of this validation set? That is the next question.

Is the gain real? A bootstrap interval

The validation set is one sample of patients. Another sample would give slightly different numbers. A bootstrap estimates how much. Draw a new validation set of the same size from the old one, with replacement, and measure the gap again. Repeat this 1,000 times. The middle 95% of the gaps is a 95% interval, as in Compare alternatives: An interval for the difference in the mathematics course. One detail matters here: we draw patients, not rows. A patient’s stays are not independent, so they stay together, for the same reason that we split by patient.

def bootstrap_gap(metric, a, b, repeats=1000):
    """metric(b) - metric(a) on 1,000 resampled validation sets. It resamples patients, not rows,
    because the encounters of one patient are not independent."""
    rows_of = list(valid.groupby("patient_nbr").indices.values())
    y = valid[TARGET].to_numpy()
    rng = np.random.default_rng(0)
    gaps = []
    for _ in range(repeats):
        rows = np.concatenate([rows_of[i] for i in rng.integers(0, len(rows_of), len(rows_of))])
        gaps.append(metric(y[rows], b[rows]) - metric(y[rows], a[rows]))
    return np.array(gaps)


def show_gap(name, metric, a, b):
    gaps = bootstrap_gap(metric, a, b)
    low, high = np.percentile(gaps, [2.5, 97.5])
    observed = metric(valid[TARGET], b) - metric(valid[TARGET], a)
    print(f"{name}: gap {observed:+.4f}, 95% interval {low:+.4f} to {high:+.4f}")

Run it for three gaps. It takes about 10 seconds:

logistic_score, boosting_score = valid_scores["E03 logistic"], valid_scores["E06 boosting"]
show_gap("boosting - logistic, AP           ", average_precision_score, logistic_score, boosting_score)
show_gap("boosting - logistic, precision@10%", precision_at, logistic_score, boosting_score)
show_gap("logistic - rule, precision@10%    ", precision_at, valid_scores["E02 rule"], logistic_score)

The output:

boosting - logistic, AP           : gap +0.0085, 95% interval +0.0039 to +0.0132
boosting - logistic, precision@10%: gap +0.0069, 95% interval -0.0055 to +0.0168
logistic - rule, precision@10%    : gap +0.0208, 95% interval +0.0064 to +0.0360

Bootstrap intervals (1,000 resamples of patients) for the gaps on the validation set. Boosting ranks better than logistic regression overall, but at the 10% capacity the interval includes zero. Logistic regression beats the rule at the capacity.

Text description of the chart

Two panels. The first, ranking by average precision: boosting minus logistic regression is +0.0085, with a 95% interval from +0.0039 to +0.0132, entirely above zero. The second, precision of the top 10%: boosting minus logistic regression is +0.0069, with an interval from −0.0055 to +0.0168, which includes zero. Logistic regression minus the rule is +0.0208, with an interval from +0.0064 to +0.0360, entirely above zero.

What this tells us. Three answers, each with its evidence:

1. Boosting ranks better overall. Its AP gap of +0.0085 has an interval from +0.004 to +0.013, entirely above zero. This confirms the folds.
2. At Nadia’s capacity, we cannot tell them apart. The gap in precision@10% is +0.007 (14 patients in 2,017 calls), but the interval runs from −0.006 to +0.017. On another sample of patients, logistic regression could just as well be ahead.
3. Logistic regression clearly beats the rule. Its gain of +0.021 (42 patients) has an interval from +0.006 to +0.036. A model is worth using.

Is the gain worth the complexity?

Use the checklist of When a simple model is better, question by question.

• Must a person get reasons, or must a reviewer check the model? We know of no law that asks for reasons for a follow-up call. The call is an offer of care, not a refusal (this is not legal advice). But the nurses need a reason to start each call (“you have had four hospital stays this year”). And the hospital’s clinical governance must approve a new tool. Hospitals often use simple point scores for readmission risk for this reason. Two well-known examples are the LACE index (van Walraven et al., 2010) and the HOSPITAL score (Donzé et al., 2013).
• Are the inputs images, text or sound? No, a table of 19 columns.
• Does the complex model win by more than one standard error? In ranking, yes: by about two standard errors in CV, and with a bootstrap interval above zero on valid.
• Is the gain worth the cost? At the decision that matters, the gain is +0.007 in precision, and its interval includes zero. Now look at the cost:

import io
import time

import joblib


def size_and_speed(pipeline):
    """The size of the saved pipeline (KB) and the time to score the validation set (seconds)."""
    buffer = io.BytesIO()
    joblib.dump(pipeline, buffer)
    start = time.perf_counter()
    pipeline.predict_proba(valid[FEATURES])
    return round(buffer.tell() / 1024), round(time.perf_counter() - start, 3)


for name in ["E03 logistic", "E05 forest", "E06 boosting"]:
    print(name, "size KB, seconds to score valid:", size_and_speed(models[name]))

The output (the times depend on your computer):

E03 logistic size KB, seconds to score valid: (12, 0.02)
E05 forest size KB, seconds to score valid: (14936, 0.048)
E06 boosting size KB, seconds to score valid: (79, 0.035)

What this tells us. All three are fast: less than a tenth of a second for 20,000 patients. Size and speed are not the real cost here. The forest is 15 MB, against 79 KB for boosting and 12 KB for logistic regression, which confirms that boosting was the right challenger. The real cost of boosting is elsewhere. Its score comes from 100 trees, so it cannot give each call an exact list of reasons. After a change at the hospital, it is harder to see what changed in 100 trees than in one weighted sum. And clinical governance must review it.

The decision. Felix asked for evidence, and here it is in three sentences. Boosting ranks patients better than logistic regression, by more than the noise. At the 10% capacity, the gain is 14 patients in 2,017 calls, and the interval includes zero. For that uncertain gain, the team would give up exact reasons for each call. So we choose logistic regression (E03) for a pilot, and keep boosting (E06) as a documented challenger. The challenger gets a second chance when the pilot has new data (see Recommendations).

This is the one-standard-error rule and the checklist of When a simple model is better working together. The rule told us that the ensembles really are better at ranking, so we do not pretend that it is a tie. The checklist told us that “better at ranking” is not yet “better for Nadia”. Update the log with the decisions and the valid AP of the two finalists:

log_experiment("E03", "Logistic regression; C tuned over 0.01, 0.1, 1", FOLD_SCORES["E03"],
               valid_ap=table.loc["E03 logistic", "AP"], decision="keep",
               reason="Chosen: at 10% capacity it ties boosting, and each call gets reasons.")
log_experiment("E06", "Gradient boosting, 100 trees; depth tuned over 2, 3, no limit", FOLD_SCORES["E06"],
               valid_ap=table.loc["E06 boosting", "AP"], decision="reject",
               reason="Better AP beyond one SE, but no clear gain at 10%; kept as the challenger.")
log_experiment("E04", "Decision tree, min leaf 100; depth tuned over 4, 6, 8, 10", FOLD_SCORES["E04"],
               decision="reject", reason="Below logistic regression on every fold.")
log_experiment("E05", "Random forest, 100 trees; min leaf tuned over 50, 20, 10", FOLD_SCORES["E05"],
               decision="reject", reason="Ties boosting on the folds; bigger and slower than boosting.")
show_log()

The output:

 id                                                        change  cv_ap  cv_std valid_ap decision
E01                    Baseline: the same score for every patient  0.114   0.000          baseline
E02                                Baseline: the prior-stays rule  0.182   0.009          baseline
E03                Logistic regression; C tuned over 0.01, 0.1, 1  0.210   0.010    0.206     keep
E04     Decision tree, min leaf 100; depth tuned over 4, 6, 8, 10  0.202   0.007            reject
E05      Random forest, 100 trees; min leaf tuned over 50, 20, 10  0.221   0.010            reject
E06 Gradient boosting, 100 trees; depth tuned over 2, 3, no limit  0.219   0.009    0.214   reject
E07         E06 scored with folds by row (patients in both parts)  0.219   0.007            reject
E08     E03 plus encounters_in_file (counts later encounters too)  0.244   0.008            reject

What this tells us. The log tells the whole story in eight rows, including the two traps. Anyone who reads it can see what was tried, how it scored and why it was kept or rejected. The reason column, in the file, says that E06 is the challenger.

The threshold for 10%

The team needs a threshold: call every patient whose score is at or above it. Set it on the validation set, as in Thresholds and calibration: the score of the 2,017th patient.

score = valid_scores["E03 logistic"]
threshold = np.sort(score)[::-1][calls - 1]
valid["score"] = score
valid["called"] = valid["score"] >= threshold
print("threshold:", round(threshold, 3), "| called:", valid["called"].sum(), "of", len(valid))

The output:

threshold: 0.18 | called: 2017 of 20179

What this tells us. A patient with a score of 0.18 or more gets a call. Exactly 2,017 patients in the validation set reach it. On new data, the share will move a little; the final test shows by how much.

A reason for each call

With logistic regression, each score splits into contributions, as in Explain predictions. Take a typical called patient, the one in the middle of the call list:

def contributions(pipeline, patient):
    """Coefficient times prepared value, for each column of one patient (as in Explain predictions)."""
    prepare, model = pipeline.named_steps["prepare"], pipeline.named_steps["model"]
    values = np.asarray(prepare.transform(patient[FEATURES]))[0]
    names = [name.split("__", 1)[1][:40] for name in prepare.get_feature_names_out()]
    table = pd.DataFrame({"value": values, "coefficient": model.coef_[0]}, index=names)
    table["contribution"] = table["value"] * table["coefficient"]
    return table


called = valid[valid["called"]].sort_values("score")
patient = called.iloc[[len(called) // 2]]            # a typical called patient: the middle of the list
print("score:", round(patient["score"].iloc[0], 3), "| readmitted:", patient[TARGET].iloc[0])
print(patient[["age", "number_inpatient", "number_emergency", "num_medications"]].to_string(index=False))
reasons = contributions(logistic, patient)
print(reasons[reasons["value"] != 0].sort_values("contribution", ascending=False).head(5).round(2))

The output:

score: 0.226 | readmitted: 1
    age  number_inpatient  number_emergency  num_medications
[80-90)                 4                 0               20
                                     value  coefficient  contribution
number_inpatient                      2.69         0.33          0.89
diabetesMed_Yes                       1.00         0.07          0.07
age_years                             1.21         0.06          0.07
admission_source_Physician Referral   1.00         0.07          0.07
number_diagnoses                      0.82         0.07          0.06

What this tells us. This patient, in their 80s, had four hospital stays in the past year. That one fact adds +0.89 to the model’s sum, more than ten times any other contribution. The nurse can open the call with it. Remember the warning of Explain predictions: a contribution describes the model, compared with an average patient. It is not the cause of a readmission. This patient did come back within 30 days.

Your turn: a capacity of 20%

Change one condition. Nadia’s budget may grow, so that the team can call 20% of patients. In the notebook, compute precision_at(..., share=0.20) on valid for the rule, logistic regression and boosting. Before you run it, predict: does the gap between the models grow or shrink?

Show the answer

The rule reaches 0.201, logistic regression 0.212, and boosting 0.219. All three fall, because the second 10% of the list has fewer patients at risk. The gap between boosting and logistic regression is 0.007, the same as at 10%. With twice the calls, that is about 28 more patients reached per 4,035 calls. Before you decide, run the bootstrap again at 20%. A gain of the same size is still small against the noise. The point: a decision that depends on the capacity must be checked again when the capacity changes.

Errors and slices

A score for the whole validation set can hide groups for whom the model works badly. Slice analysis showed how to look. Every slice table here has six columns. They are the number of stays, the readmission rate, the share called, the precision and recall of the calls, and the ROC AUC in the group.

By age

def slice_table(df, column):
    """Per group: encounters, readmission rate, share called, precision and recall of the calls, AUC."""
    groups = df.groupby(column, observed=True)
    return pd.DataFrame({
        "encounters": groups.size(),
        "readmitted": groups[TARGET].mean(),
        "called": groups["called"].mean(),
        "precision": groups.apply(lambda g: g.loc[g["called"], TARGET].mean()),
        "recall": groups.apply(lambda g: g.loc[g[TARGET] == 1, "called"].mean()),
        "AUC": groups.apply(lambda g: roc_auc_score(g[TARGET], g["score"]) if g[TARGET].nunique() == 2 else np.nan),
    }).round(3)


valid["age_band"] = valid["age"].replace({"[0-10)": "[0-40)", "[10-20)": "[0-40)", "[20-30)": "[0-40)",
                                          "[30-40)": "[0-40)"})
slice_table(valid, "age_band")

The output:

          encounters  readmitted  called  precision  recall    AUC
age_band
[0-40)          1201       0.097   0.081      0.381   0.319  0.762
[40-50)         1985       0.102   0.078      0.277   0.213  0.713
[50-60)         3427       0.092   0.074      0.249   0.199  0.680
[60-70)         4354       0.110   0.092      0.252   0.211  0.665
[70-80)         5339       0.129   0.113      0.255   0.222  0.647
[80-90)         3343       0.125   0.126      0.276   0.278  0.643
[90-100)         530       0.106   0.172      0.143   0.232  0.559

How well the model ranks patients within each age band of the validation set. The ranking gets weaker with age; for patients aged 90 or more it is close to guessing (0.5).

Text description of the chart

A column chart of the ROC AUC of logistic regression within each age band of the validation set: 0 to 40, 0.76 (1,201 stays); 40 to 50, 0.71 (1,985); 50 to 60, 0.68 (3,427); 60 to 70, 0.67 (4,354); 70 to 80, 0.65 (5,339); 80 to 90, 0.64 (3,343); 90 to 100, 0.56 (530).

What this tells us. The model ranks younger patients best (AUC 0.76 under 40). It gets weaker with every decade, down to 0.56 for patients aged 90 or more, close to guessing. Two things happen at once for the oldest group. The model calls them most often (17.2%), because age and earlier stays raise the score. But its calls are the least precise there. Only 14.3% of the calls reach a patient who comes back, close to the group’s own rate of 10.6%. Those 530 stays are a small slice, so the numbers move by chance. But the trend across all seven bands is clear. For the oldest patients, the nurses should not trust the score much; their own judgement matters more.

By the sensitive columns

Race, gender and payer code are not inputs. Do they still change who gets a call, or how good the calls are?

valid["race_group"] = valid["race"].fillna("missing")
valid["payer_group"] = valid["payer_code"].fillna("missing").where(
    valid["payer_code"].isin(["MC", "HM", "SP", "BC", "MD"]) | valid["payer_code"].isna(), "other")
for column in ["race_group", "gender", "payer_group"]:
    print(slice_table(valid, column).to_string(), end="\n\n")

The output:

                 encounters  readmitted  called  precision  recall    AUC
race_group
AfricanAmerican        3849       0.113   0.101      0.237   0.211  0.665
Asian                   121       0.083   0.083      0.100   0.100  0.724
Caucasian             15058       0.115   0.103      0.269   0.241  0.673
Hispanic                419       0.088   0.062      0.385   0.270  0.674
Other                   293       0.092   0.051      0.200   0.111  0.551
missing                 439       0.084   0.068      0.133   0.108  0.659

                 encounters  readmitted  called  precision  recall    AUC
gender
Female                10914       0.113   0.106      0.247   0.233  0.669
Male                   9264       0.112   0.092      0.280   0.229  0.673
Unknown/Invalid           1       0.000   0.000        NaN     NaN    NaN

             encounters  readmitted  called  precision  recall    AUC
payer_group
BC                  900       0.080   0.050      0.200   0.125  0.685
HM                 1323       0.121   0.106      0.264   0.231  0.650
MC                 6365       0.120   0.134      0.249   0.278  0.649
MD                  659       0.096   0.117      0.273   0.333  0.675
SP                  989       0.100   0.080      0.253   0.202  0.706
missing            8038       0.116   0.089      0.282   0.217  0.688
other              1905       0.096   0.056      0.226   0.131  0.636

What this tells us. Read each table with the counts in mind, as Slice analysis taught.

• Race. The two large groups come back at the same rate (11.3% and 11.5%). The model calls them at the same rate too (10.1% and 10.3%). But the calls are less precise for African American patients (0.237 against 0.269). They also catch fewer of those who return (recall 0.211 against 0.241). With 3,849 stays, that gap is worth watching. The model’s inputs seem to describe the risk of white patients a little better. The smaller groups (121 to 439 stays) are too small to judge.
• Gender. Women and men come back at the same rate (11.3% and 11.2%). The recall is almost the same (0.233 and 0.229). The model calls women a little more often (10.6% against 9.2%), so their calls are a little less precise. One stay has the gender Unknown/Invalid.
• Payer. The file does not say what the codes mean. MC is probably Medicare, the US insurance for people aged 65 or more. MD is probably Medicaid, for people with low income. The model calls the MC group most (13.4%), because it is older. It calls BC and “other” least (5.0% and 5.6%), with low recall. This partly follows their lower readmission rates (8.0% and 9.6%).

None of these gaps comes from a sensitive input, because there is none. They come through other inputs, such as age and earlier stays, and through what the data does not record. The right response is not to hide the table. Report it, watch the race gap in the pilot, and let the nurses add patients whom the list misses.

Who does the model miss?

The calls reach 23% of the patients who come back. Who are the others?

returned = valid[valid[TARGET] == 1]
first = returned["number_inpatient"] == 0
print(f"patients who return: {len(returned)}; with no stay in the past year: {first.sum()} ({first.mean():.0%})")
print(f"recall among them: {returned.loc[first, 'called'].mean():.3f};"
      f" among the others: {returned.loc[~first, 'called'].mean():.3f}")

The output:

patients who return: 2273; with no stay in the past year: 1116 (49%)
recall among them: 0.058; among the others: 0.398

What this tells us. This is the most important error of the model. Half of the patients who come back (1,116 of 2,273) had no hospital stay in the past year. The model calls only 5.8% of them, against 39.8% of the patients with earlier stays. The model, like the rule, mostly finds patients who come back often. For a patient on a first stay, the file has little that predicts a return. That is a limit of the data, not of the model family: all four families depend on the same columns.

The final test, once

Everything is decided: the features, the model (logistic regression, C = 0.01, fitted on train), and the threshold (0.18). Now open the sealed test set, once, as in Thresholds and calibration: Test once. Whatever it shows, we do not go back and change the model.

test_score = logistic.predict_proba(test[FEATURES])[:, 1]
called = test_score >= threshold
returned = test[TARGET].to_numpy() == 1
print(f"test: {len(test)} encounters, {returned.sum()} returned within 30 days")
print(f"called with the valid threshold: {called.sum()} ({called.mean():.1%})")
print(f"precision {returned[called].mean():.3f}, recall {called[returned].mean():.3f}")
print(f"AP {average_precision_score(returned, test_score):.3f}, ROC AUC {roc_auc_score(returned, test_score):.3f},"
      f" precision@10% {precision_at(returned, test_score):.3f}")
print(f"the rule on test: precision@10% {precision_at(returned, stays_rule(test)):.3f}")

The output:

test: 19608 encounters, 2278 returned within 30 days
called with the valid threshold: 1870 (9.5%)
precision 0.270, recall 0.221
AP 0.217, ROC AUC 0.663, precision@10% 0.264
the rule on test: precision@10% 0.247

What this tells us. On 19,608 stays of patients that no step has seen, the model behaves as on valid:

• The threshold from valid calls 9.5% of the test patients, close to the planned 10%. In practice, the share moves from week to week, so the team should count the calls and adjust.
• About 27 calls in 100 reach a patient who comes back (precision 0.270). At random, about 11 in 100 would; with the rule, 25. The calls reach 22% of all returns.
• AP 0.217 and ROC AUC 0.663 are close to the validation numbers (0.206 and 0.671). The model did not overfit the choices that we made on valid.

This is the number to report: precision 0.27 at a 10% capacity, against 0.25 for the rule of thumb and 0.11 at random.

Recommendations

For Nadia and Felix, in their words:

1. Start a pilot with the logistic regression model. Call the patients with a score of 0.18 or more, about 1 in 10. Expect about 1 call in 4 to reach a patient who would come back, against 1 in 9 at random. Compared with today’s rule, that is about 40 more patients reached per 2,000 calls.
2. Give each call its reasons. Show the nurse the two or three largest contributions with the patient’s own values, for example “4 hospital stays this year”. These are the model’s reasons, not a diagnosis.
3. Let the nurses add patients. The model misses most patients who come back without earlier stays. A nurse who knows of a reason (no help at home, a new insulin plan) should be able to add a patient to the list.
4. Measure the pilot properly. Record each called patient’s outcome. The best evidence is a comparison. For a few weeks, choose some calls by the model and some by the rule, and compare the readmission rates. Then decide whether to keep the model.
5. Give the challenger a second chance. Boosting ranked better in every comparison of the whole ranking, but the gain at 10% was within the noise. When the pilot has several thousand new outcomes, run the comparison again on the new data. Switch only if the gain at the capacity is clearly above zero and worth the loss of exact reasons.
6. Watch the slices. Report precision and recall by age, race, gender and payer every month. If the gap for African American patients grows, or the oldest patients are called often with poor precision, review the model.
7. Retrain on your own, recent data. This model learned from US hospitals in 1999–2008. Before any real use, a hospital must train and test it on its own patients.

What not to use the model for. Do not use it to decide who may leave hospital, who gets less care, or what a patient pays. It does not diagnose anything. Do not use it for patients without diabetes, for children, or in a hospital whose data it has not been tested on. A low score does not mean “safe”: half of the patients who come back have no earlier stay, and most of them get low scores.

Limits and next steps

Be honest about the ceiling. This data predicts readmission only weakly. The best model in this case study reaches a ROC AUC of about 0.67, where 0.5 is guessing and 1.0 is perfect. Three facts show that the limit is the data, not the choice of model:

• Four very different families, each tuned, end within 0.02 AP of each other (0.202 to 0.221). If a strong hidden pattern existed in these columns, the flexible models would have found much more of it.
• The rule of thumb, with one column, gets most of the way: 0.240 against 0.261 at the capacity.
• Half of the readmissions come from patients without earlier stays, and the model finds almost none of them.

Readmission depends on things that this file does not record. Examples are who helps at home, whether the patient can get the medicines, and whether a follow-up visit happens. The file also says little about how ill the patient really is: it has no lab values except two diabetes tests. Better data would help more than a better algorithm.

Other limits:

• No dates. We could not split by time, so we do not know how the model ages. A real project must test on the most recent months.
• Old data. US hospitals, 1999–2008, ICD-9 codes. Care, coding and insurance have changed.
• One label. <30 records a return to one of the hospitals in the data. A return to another hospital may look like “no return”.
• Correlated stays. We kept all stays of a patient. A model for first stays only, as in Strack et al., would answer a slightly different question.

Next steps, in order of expected value:

1. Add the second and third diagnoses and a count of chronic conditions (a comorbidity score, as in LACE).
2. Add the medicine columns.
3. Calibrate the score, so that 0.25 means “1 in 4”, as in Thresholds and calibration.
4. Plan the pilot as an experiment.

Failure case: a count that looks into the future

Common mistake

A feature that uses later stays

Symptom: Dev, an engineer on another team, looks at your notebook. “Patients with many stays come back more often. Why not count all their stays?” He adds a feature, encounters_in_file: the number of rows that the patient has in the whole file. The CV AP jumps from 0.210 to 0.244, a bigger gain than any model family gave:

in_file = data["patient_nbr"].map(data["patient_nbr"].value_counts())   # all encounters, also later ones
train["encounters_in_file"] = in_file.loc[train.index]
leaky = build_pipeline(LogisticRegression(C=0.01, max_iter=1000), numbers=NUMBERS + ["encounters_in_file"])
leaky_scores = cross_val_score(leaky, train[FEATURES + ["encounters_in_file"]], train[TARGET], cv=folds,
                               scoring="average_precision")
print("with encounters_in_file:", leaky_scores.round(4), round(leaky_scores.mean(), 4))
print("without (E03):          ", FOLD_SCORES["E03"].round(4), round(FOLD_SCORES["E03"].mean(), 4))
log_experiment("E08", "E03 plus encounters_in_file (counts later encounters too)", leaky_scores,
               decision="reject", reason="A leak: at discharge, nobody knows the later encounters.")

with encounters_in_file: [0.2543 0.2307 0.2389 0.2463 0.2495] 0.2439
without (E03):           [0.2064 0.204  0.205  0.2044 0.2302] 0.21

Cause: the count includes stays that happen after the one we predict. A patient who comes back within 30 days has, by definition, at least one more stay in the file. So the feature partly contains the answer. The folds by patient do not protect against this, because the leak is in time, not across patients. At discharge, the nurses cannot know how many future stays a patient will have.

Fix: ask, for every feature, “would we know this on the day of discharge?” Use only counts of the past, such as number_inpatient, which the source computed for the year before the stay. A gain that is much larger than every honest change in the log is a warning sign: look for a leak before you celebrate. Log the experiment as rejected, with the reason, so that nobody tries it again.

Consequence in AI work: in the pilot, a new patient’s future stays do not exist yet, so the feature would be 1 for most of them. The scores would fall, and the calls would go to the wrong patients. The team would trust a number (0.244) that was never real. Leaks through time are among the most common reasons why a model that looked good in a notebook fails in use.

Make it a portfolio piece

This case study is a complete project that an employer can read in ten minutes. Put it in a public repository.

What to put in the repository:

readmission-model/
  README.md               the summary below
  MODEL_CARD.md           intended use, not for, data, evaluation, slices, limits
  DATA.md                 the dataset card: source, DOI, licence, attribution, your changes
  notebook.ipynb          the notebook, run from top to bottom, with its outputs
  experiment_log.csv      every experiment, E01 to E08
  requirements.txt        pandas==3.0.6, scikit-learn==1.9.1, matplotlib

You do not need to commit the data file itself. CC BY 4.0 allows it, with credit and a note of your changes, but a link and a checksum are enough.

A README template:

Show the README template

# Who should the care team call? A 30-day readmission model

**Problem.** A care-transitions team can call 10% of discharged patients with diabetes.
Which patients are most likely to come back to hospital within 30 days?

**Data.** Diabetes 130-US hospitals 1999–2008 (UCI, doi:10.24432/C5230J, CC BY 4.0;
Strack et al. 2014). 99,343 stays of 69,990 patients after removing deaths and hospice
discharges. Changes: see DATA.md. Race, gender and payer are used for checks only.

**Approach.** Split by patient (60/20/20). Baselines: the same score for all, and the
rule of thumb "most hospital stays in the past year". Four model families in one
pipeline, tuned on the same 5 patient folds (average precision), one-standard-error rule.
Bootstrap intervals (by patient) for the gaps on validation. Test set used once.

**Results (test set, 19,608 stays).** At a 10% capacity, precision 0.27 (random 0.11,
rule of thumb 0.25); ROC AUC 0.66. Gradient boosting ranked better than logistic
regression (AP +0.009), but at the capacity the gain was within the noise, so the
readable model was chosen.

**Limits.** Weak signal (half of the readmissions have no earlier stay and are rarely
found); no dates, so no time split; US data from 1999–2008. Not for decisions about
discharge or care.

The model card: write MODEL_CARD.md with the sections from Document use boundaries:

• Intended use: rank discharged patients with diabetes for a follow-up call, at a 10% capacity, with a nurse who reviews the list.
• Not for: discharge decisions, less care, prices, other patients.
• Data: the dataset card and your changes.
• Why this model: the decision of this case study, with the numbers.
• Evaluation: the test numbers and the baselines.
• Slices: the age table and the three sensitive tables.
• Limits: the ceiling, no dates, old data.

How to describe it in an interview, in two sentences: “I built a 30-day readmission model for a care team that can call 10% of patients, on 100,000 real hospital stays. I compared four model families under one patient-level protocol and chose readable logistic regression, because boosting’s small gain at the team’s capacity was within the noise.”

Expect the follow-up question: “Why not the best model?” Your answer is the evidence in Is the gain worth the complexity?.

Go further:

• Add diag_2, diag_3 and the 22 medicine columns. Does any family gain more than one standard error?
• Build a comorbidity score from the three diagnoses, and test a LACE-like point score as a third baseline.
• Calibrate the logistic regression and draw a reliability chart.
• Train on first stays only, as Strack et al. did, and compare.
• Use encounter_id as an approximate time order, label it clearly as an assumption, and test on the last 20% of stays.
• Write a one-page pilot plan: how many weeks, how to compare model and rule, and what result would make you stop.

Knowledge check

Check your understanding

Boosting has a higher cross-validation AP than logistic regression by about two standard errors, and is higher on all five folds. Felix says: then the one-standard-error rule chooses boosting, so we must use it. What is right and what is missing in his argument?

Show the answer

He is right that the rule does not allow a tie here. Logistic regression (0.2100) is below the bar of 0.2162, and the bootstrap interval for the AP gap on valid (+0.004 to +0.013) is above zero. Boosting really ranks better overall. What is missing is the decision. The rule compares the summary metric; the team uses only the top 10%. There, the gain is +0.007 in precision (14 patients in 2,017 calls), and its interval (−0.006 to +0.017) includes zero. The checklist of When a simple model is better then asks whether the gain is worth the cost. For an uncertain gain, boosting would give up exact reasons for each call and make monitoring and review harder. So the evidence supports logistic regression for the pilot, with boosting as the challenger.

Check your understanding

A colleague splits the stays by row and reports a test AP. Here, folds by row changed the boosting CV AP by only 0.0005. Is a split by row acceptable for this project? Explain.

Show the answer

No. With a row split, 41% of the checked rows belong to patients that the model trained on. The score then answers “how well does the model know these patients?” instead of “how well does it predict new ones?”. On this data the change is small because the inputs describe the stay, not the person. But that is a property of these columns, and nobody can rely on it in advance. A fixed personal input or a richer history would make the leak large. The protocol must measure the question that the team asks, so split by patient (and by time, when there are dates).

Check your understanding

Nadia asks: the model reaches only about 27 calls in 100 that matter. Should we wait for a better algorithm? Use the evidence of this case study.

Show the answer

A better algorithm is unlikely to help much. Four very different families end within 0.02 AP of each other, and the flexible ones gain little at the capacity. The rule of thumb with one column reaches 0.247 on test, close to the model’s 0.270. Half of the patients who return had no earlier stay, and the model finds only 5.8% of them. The file does not record what predicts their return: support at home, medicines, follow-up care. Better data would help more. Meanwhile, the model is clearly better than random (0.11) and better than the rule (interval +0.006 to +0.036 at the capacity). So a pilot is worth it, with nurses able to add patients.

Sources

• Clore, J., Cios, K., DeShazo, J., & Strack, B. (2014). Diabetes 130-US Hospitals for Years 1999-2008 [Dataset]. UCI Machine Learning Repository. doi:10.24432/C5230J. CC BY 4.0.
• Strack, B., DeShazo, J. P., Gennings, C., Olmo, J. L., Ventura, S., Cios, K. J., & Clore, J. N. (2014). Impact of HbA1c measurement on hospital readmission rates: analysis of 70,000 clinical database patient records. BioMed Research International, 781670. doi:10.1155/2014/781670.
• Obermeyer, Z., Powers, B., Vogeli, C., & Mullainathan, S. (2019). Dissecting racial bias in an algorithm used to manage the health of populations. Science, 366(6464), 447–453. doi:10.1126/science.aax2342.
• van Walraven, C., et al. (2010). Derivation and validation of an index to predict early death or unplanned readmission after discharge from hospital to the community (the LACE index). CMAJ, 182(6), 551–557. doi:10.1503/cmaj.091117.
• Donzé, J., Aujesky, D., Williams, D., & Schnipper, J. L. (2013). Potentially avoidable 30-day hospital readmissions in medical patients: derivation and validation of a prediction model (the HOSPITAL score). JAMA Internal Medicine, 173(8), 632–638. doi:10.1001/jamainternmed.2013.3023.

Recap

• Prepare real data with a reason for each step. ? is missing, an empty test is “not measured”, codes become words, and stays that end in death or hospice go.
• Frame the prediction at discharge: inputs known on that day, no future counts, and sensitive columns kept for slice checks only.
• Split by patient, because 46% of the rows belong to repeat patients. Here a row split changed the score little, but it measures the wrong question.
• One protocol for four families: the same 5 patient folds, the same metric (AP), the same pipeline, every attempt in the log. Boosting and the forest beat logistic regression by about two standard errors.
• Decide on the decision. At the 10% capacity, boosting’s gain (+0.007) has a bootstrap interval that includes zero. So the readable model won. On test: precision 0.27, against 0.25 for the rule and 0.11 at random.
• Be honest about the ceiling: this data predicts readmission only weakly, and half of the returns come from patients that the model rarely finds.

Next, a case study without a target. You find groups of online-shop customers by how recently, how often and how much they buy. Then you check the groups against what the customers do next.

Next: Segment online-shop customers