In healthcare, the “best” metric depends on the cost of being wrong. A model that misses a disease can harm patients; a model that flags too many healthy patients overwhelms clinicians and causes unnecessary testing. Accuracy alone is rarely enough.
This article uses one running example: predicting whether a patient is at high risk of 30-day hospital readmission. The positive class is readmitted = 1, and it is the minority class. Every number quoted comes from running the code shown.
Setup: shared healthcare dataset
The setup creates a synthetic dataset with about 10% readmissions — the features and labels are artificial and only demonstrate the calculations; they are not a clinically valid readmission dataset. Save it as readmission_data.py; every later section imports from it.
Two choices here differ from many tutorials, and both matter for every metric that follows. The model is not class-weighted: class_weight="balanced" inflates predicted probabilities (in this example the mean predicted risk rose from 0.11 to 0.36 for a 10% event), worsened the Brier score from 0.071 to 0.164, and lowered average precision from 0.487 to 0.425. And predictions are not made with predict(), which silently applies a 0.5 cut-off. Instead the operating threshold is chosen explicitly from a clinical requirement — find at least 70% of readmissions — using out-of-fold predictions on the training data, so the test set is used only to report results. The file also defines a bootstrap helper for confidence intervals.
# plot_chance_level on the display classes needs scikit-learn 1.4+
pip install "scikit-learn>=1.4" pandas numpy matplotlib
# Save as readmission_data.py: every later section imports from it.
import numpy as np
import pandas as pd
from sklearn.datasets import make_classification
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import precision_recall_curve
from sklearn.model_selection import (StratifiedKFold, cross_val_predict,
train_test_split)
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler
SEED = 42
TARGET_RECALL = 0.70 # clinical requirement: find at least 70% of cases
X, y = make_classification(
n_samples=5_000, n_features=12, n_informative=7, n_redundant=2,
weights=[0.90, 0.10], class_sep=0.9, random_state=SEED)
feature_names = [
"age", "prior_admissions_12m", "length_of_stay", "num_medications",
"num_comorbidities", "emergency_admission", "icu_stay",
"lab_abnormalities", "followup_scheduled", "social_support_score",
"prior_ed_visits", "discharge_delay_hours"]
patients = pd.DataFrame(X, columns=feature_names)
patients["readmitted"] = y
X_train, X_test, y_train, y_test = train_test_split(
patients[feature_names], patients["readmitted"], test_size=0.25,
stratify=patients["readmitted"], random_state=SEED)
# Unweighted, so predict_proba stays close to real-world prevalence
model = make_pipeline(StandardScaler(), LogisticRegression(max_iter=1000))
def threshold_for_recall(y_true, scores, target):
"""Highest threshold whose recall still meets the target."""
_, recall, thresholds = precision_recall_curve(y_true, scores)
meets = recall[:-1] >= target # recall[i] belongs to thresholds[i]
if not meets.any():
raise ValueError(f"no threshold reaches recall {target}")
return float(thresholds[meets].max())
# Choose the operating threshold on out-of-fold TRAINING predictions,
# so the test set is used only to report results.
cv = StratifiedKFold(n_splits=5, shuffle=True, random_state=SEED)
oof = cross_val_predict(model, X_train, y_train, cv=cv,
method="predict_proba")[:, 1]
THRESHOLD = threshold_for_recall(y_train, oof, TARGET_RECALL)
model.fit(X_train, y_train)
y_score = model.predict_proba(X_test)[:, 1]
y_pred = (y_score >= THRESHOLD).astype(int) # not predict(), which uses 0.5
def bootstrap_ci(metric, y_true, values, n_boot=1_000, alpha=0.05,
seed=SEED):
"""Percentile bootstrap CI for any metric(y_true, values)."""
rng = np.random.default_rng(seed)
y_true, values = np.asarray(y_true), np.asarray(values)
stats = []
for _ in range(n_boot):
idx = rng.integers(0, len(y_true), len(y_true))
if y_true[idx].min() == y_true[idx].max():
continue # a resample with one class has no AUC
stats.append(metric(y_true[idx], values[idx]))
return tuple(np.quantile(stats, [alpha / 2, 1 - alpha / 2]))
if __name__ == "__main__":
print(f"Prevalence {y_test.mean():.3f} threshold {THRESHOLD:.3f}")The chosen threshold is 0.149 — far from 0.5, as is normal for a rare outcome. Many clinically important events (readmission, sepsis, cancer recurrence, adverse drug events) are uncommon, which is why the threshold-dependent metrics below are always quoted with their threshold.
1. Accuracy
Problem statement
A hospital wants a simple first-pass measure of how often its readmission-risk model makes the correct decision.
Concept
Accuracy is the share of all predictions that are correct:
Accuracy = (TP + TN) / (TP + TN + FP + FN)It is easy to communicate and misleading when one class dominates. The code compares the model with a rule that predicts “no readmission” for everyone.
Code
from sklearn.metrics import accuracy_score, confusion_matrix
from readmission_data import y_pred, y_test
accuracy = accuracy_score(y_test, y_pred)
always_no = accuracy_score(y_test, [0] * len(y_test))
tn, fp, fn, tp = confusion_matrix(y_test, y_pred).ravel()
print(f"Model accuracy: {accuracy:.3f}")
print(f"'Never readmitted' rule: {always_no:.3f}")
print(f"TP={tp} FP={fp} FN={fn} TN={tn}")The useful model scored 0.828 accuracy; the useless rule scored 0.896 while finding no readmissions at all. A model tuned to catch 70% of readmissions must accept false positives (176 here), and accuracy punishes it for exactly the trade-off the hospital wants. Show accuracy only with the confusion matrix, the class distribution, and the threshold; do not compare it across datasets with different prevalence; and monitor it by subgroup, site, and period only alongside class-specific metrics.
2. Precision
Problem statement
The care-management team can follow up only with patients the model flags. It wants to know how many flagged patients are actually readmitted.
Concept
Precision answers: of all patients predicted positive, how many were truly positive? High precision means fewer false alarms.
Precision = TP / (TP + FP)Code
from sklearn.metrics import precision_score
from readmission_data import y_pred, y_test
precision = precision_score(y_test, y_pred, zero_division=0)
print(f"Precision: {precision:.3f} "
f"(patients flagged for follow-up: {y_pred.sum()})")At the chosen threshold the model flagged 267 patients with precision 0.341: about one in three flagged patients was readmitted, against one in ten in the population. Precision means little without volume — 95% precision on five patients may be less useful than 75% on 500 — so report it with the number flagged and the exact production threshold, track weekly alert volume and clinician overrides, and evaluate it for important subgroups. Optimizing precision alone makes a model conservative and lets true cases slip through.
3. Recall
Problem statement
The hospital wants to identify as many future readmissions as possible, because a missed high-risk patient may not receive preventive follow-up.
Concept
Recall, also called sensitivity, answers: of all truly positive patients, how many did the model find? High recall means fewer missed cases, which is often critical for screening and early-warning models.
Recall = TP / (TP + FN)Code
from sklearn.metrics import recall_score
from readmission_data import y_pred, y_test
recall = recall_score(y_test, y_pred, zero_division=0)
missed = int(((y_test == 1) & (y_pred == 0)).sum())
print(f"Recall: {recall:.3f} (readmissions missed: {missed})")Recall on the test set was 0.700, matching the requirement set on training data, with 39 readmissions missed. Maximizing recall without limit produces alert fatigue, so choose the threshold from clinical cost and workflow capacity, evaluate recall at the deployed threshold, track it prospectively, audit missed cases through clinical review, and monitor it by age, sex, site, insurance status, and — where appropriate and permitted — ethnicity, since averages can hide worse recall for under-represented groups.
4. F1 score
Problem statement
The hospital wants a single threshold-dependent number that balances false alarms and missed readmissions.
Concept
F1 is the harmonic mean of precision and recall, giving them equal weight. When missing a case is worse than a false alarm, the F-beta score generalizes it: F2 weights recall twice as heavily.
F1 = 2 · Precision · Recall / (Precision + Recall)Code
from sklearn.metrics import f1_score, fbeta_score
from readmission_data import THRESHOLD, y_pred, y_score, y_test
print(f"F1 at threshold {THRESHOLD:.3f}: "
f"{f1_score(y_test, y_pred, zero_division=0):.3f}")
# F2 weights recall twice as heavily as precision
print(f"F2 at threshold {THRESHOLD:.3f}: "
f"{fbeta_score(y_test, y_pred, beta=2, zero_division=0):.3f}")
# F1 is a property of a threshold, not of the model
for t in (0.1, 0.2, 0.3, 0.5):
print(f" F1 @ {t:.1f}: "
f"{f1_score(y_test, y_score >= t, zero_division=0):.3f}")F1 was 0.458 and F2 0.578 at the deployed threshold, and F1 ranged from 0.393 to 0.488 across thresholds for the same model — F1 describes an operating point, not a model. State the threshold whenever you report it, read it next to precision and recall, use macro-averaged or per-class F1 when both classes matter, and do not treat it as a substitute for calibration, decision-curve analysis, or an explicit clinical utility function.
5. ROC-AUC
Problem statement
The analytics team wants a threshold-independent view of how well the model separates patients who will be readmitted from those who will not.
Concept
The ROC curve plots true-positive rate against false-positive rate across all thresholds, and ROC-AUC is the area beneath it: 0.5 is random ranking and 1.0 is perfect. It equals the probability that a randomly chosen positive patient scores higher than a randomly chosen negative one, with ties counted as half. It does not depend on a cut-off and does not measure calibration.
Code
import matplotlib.pyplot as plt
from sklearn.metrics import RocCurveDisplay, roc_auc_score
from readmission_data import bootstrap_ci, y_score, y_test
roc_auc = roc_auc_score(y_test, y_score)
low, high = bootstrap_ci(roc_auc_score, y_test, y_score)
print(f"ROC-AUC: {roc_auc:.3f} (95% CI {low:.3f}-{high:.3f})")
display = RocCurveDisplay.from_predictions(
y_test, y_score, name="Readmission model", plot_chance_level=True)
display.ax_.set_title("ROC curve: 30-day readmission risk")
display.figure_.savefig("roc_curve.png", dpi=150)
plt.close(display.figure_)ROC-AUC was 0.810 with a 95% bootstrap interval of 0.760–0.854 on 1,250 test patients. The interval is the point: two models whose AUCs differ by 0.02 on a test set this size are not distinguishable without a paired comparison. A higher AUC indicates better ranking across all thresholds, which is a reasonable first filter when comparing, say, logistic regression with gradient boosting before clinical validation. It is not a readmission probability, and for rare events it can look strong while precision at useful thresholds is poor — which is why it is paired with PR-AUC and calibration.
6. PR-AUC
Problem statement
With only about 10% of patients readmitted, the hospital wants a metric focused on the minority class and not inflated by the many correctly classified non-readmissions.
Concept
The precision-recall curve plots precision against recall across thresholds; it never uses true negatives. Its usual summary is average precision, which scikit-learn computes as a step-wise sum rather than a trapezoidal area:
AP = Σₙ (Rₙ − Rₙ₋₁) · PₙA random model’s AP equals the prevalence, so the meaningful comparison is with 0.104 here, not with 0.5.
Code
import matplotlib.pyplot as plt
from sklearn.metrics import (PrecisionRecallDisplay,
average_precision_score, precision_score)
from readmission_data import THRESHOLD, bootstrap_ci, y_pred, y_score, y_test
ap = average_precision_score(y_test, y_score)
low, high = bootstrap_ci(average_precision_score, y_test, y_score)
print(f"Average precision: {ap:.3f} (95% CI {low:.3f}-{high:.3f})")
print(f"Chance level = prevalence: {y_test.mean():.3f}")
print(f"Precision at deployed threshold {THRESHOLD:.3f}: "
f"{precision_score(y_test, y_pred, zero_division=0):.3f}")
# Step-wise curve; straight-line interpolation would overstate precision
display = PrecisionRecallDisplay.from_predictions(
y_test, y_score, name="Readmission model", plot_chance_level=True)
display.ax_.set_title("Precision-recall curve: 30-day readmission risk")
display.figure_.savefig("precision_recall_curve.png", dpi=150)
plt.close(display.figure_)AP was 0.487 (95% CI 0.398–0.573) against a chance level of 0.104, and precision at the deployed threshold was 0.341. The curve is drawn as steps: connecting the points with straight lines, as a plain plt.plot does, overstates precision between them. Report AP with prevalence and with precision at the production threshold, never compare it across populations with different prevalence, do not treat it as a probability or calibration measure, and reassess it when patient mix, prevalence, or care processes change.
7. Calibration and clinical utility
Problem statement
Clinicians see a risk of 0.30 for a patient. Does that mean three in ten such patients are readmitted, and does acting on the model do more good than a blanket policy?
Concept
Ranking metrics cannot answer either question. Calibration compares predicted probabilities with observed frequencies; the Brier score (mean squared error of the probabilities) summarizes it and should beat the score of always predicting the prevalence. Decision-curve analysis asks whether acting on the model at a given risk threshold has more net benefit than treating everyone or no one, weighting false positives by the odds of the threshold.
Code
import matplotlib.pyplot as plt
import numpy as np
from sklearn.calibration import CalibrationDisplay
from sklearn.metrics import brier_score_loss
from readmission_data import y_score, y_test
print(f"Brier score: {brier_score_loss(y_test, y_score):.4f} "
f"(always predicting prevalence: "
f"{brier_score_loss(y_test, np.full(len(y_test), y_test.mean())):.4f})")
print(f"Mean predicted risk {y_score.mean():.3f} vs observed "
f"{y_test.mean():.3f}")
display = CalibrationDisplay.from_predictions(
y_test, y_score, n_bins=10, strategy="quantile", name="Model")
display.figure_.savefig("calibration_curve.png", dpi=150)
plt.close(display.figure_)
def net_benefit(y_true, scores, threshold):
"""Decision-curve net benefit of treating patients above threshold."""
y_true, flag = np.asarray(y_true), np.asarray(scores) >= threshold
n = len(y_true)
tp = (flag & (y_true == 1)).sum() / n
fp = (flag & (y_true == 0)).sum() / n
return tp - fp * threshold / (1 - threshold)
for t in (0.05, 0.10, 0.20):
treat_all = net_benefit(y_test, np.ones(len(y_test)), t)
model_nb = net_benefit(y_test, y_score, t)
print(f"risk threshold {t:.2f}: model {model_nb:+.4f} "
f"treat-all {treat_all:+.4f} treat-none +0.0000")The Brier score was 0.0705 against 0.0932 for always predicting the prevalence, and mean predicted risk was 0.109 against an observed 0.104. Net benefit beat both treat-all and treat-none at every threshold tested, and the margin over treat-all grew as the threshold rose (at 0.20, treat-all is clearly harmful). If you do use class weights or resampling, recalibrate before showing probabilities to clinicians, and evaluate calibration separately from discrimination.
Choosing the right metric
| Metric | Question it answers | Threshold? | Best used when |
|---|---|---|---|
| Accuracy | How often are all decisions correct? | Yes | Classes are balanced and errors cost about the same |
| Precision | Of flagged patients, how many are positive? | Yes | False positives are costly or capacity is limited |
| Recall | Of positive patients, how many were found? | Yes | Missed cases are costly |
| F1 / F-beta | How well are precision and recall balanced? | Yes | You need one operating-point number; F2 when recall matters more |
| ROC-AUC | How well are positives ranked above negatives? | No | Comparing ranking quality, with a confidence interval |
| PR-AUC | How well is the positive class ranked? | No | The positive class is rare; compare with prevalence |
| Brier / calibration | Do predicted risks match observed rates? | No | Probabilities are shown to clinicians or drive decisions |
Complete evaluation report
A practical reporting package for the readmission model includes prevalence, the threshold and how it was chosen, precision, recall, and F1 at that threshold, ROC-AUC, PR-AUC, the Brier score, the confusion matrix, and subgroup results — each metric with a confidence interval. Threshold-dependent metrics are computed from the scores at the deployed threshold, so their intervals reflect the same decision rule.
import pandas as pd
from sklearn.metrics import (accuracy_score, average_precision_score,
brier_score_loss, classification_report,
confusion_matrix, f1_score, precision_score,
recall_score, roc_auc_score)
from readmission_data import (THRESHOLD, bootstrap_ci, y_pred, y_score,
y_test)
def at_threshold(fn):
return lambda y_true, scores: fn(y_true, scores >= THRESHOLD,
zero_division=0)
metrics = { # name: (metric(y_true, scores), interpretation)
"accuracy": (lambda t, s: accuracy_score(t, s >= THRESHOLD),
"Overall correct decisions"),
"precision": (at_threshold(precision_score),
"Flagged patients who were readmitted"),
"recall": (at_threshold(recall_score), "Readmissions the model flagged"),
"f1": (at_threshold(f1_score), "Harmonic mean of the two above"),
"roc_auc": (roc_auc_score, "Ranking quality, all thresholds"),
"pr_auc": (average_precision_score, "Positive-class ranking quality"),
"brier": (brier_score_loss, "Probability accuracy (lower = better)"),
}
rows = []
for name, (fn, meaning) in metrics.items():
low, high = bootstrap_ci(fn, y_test, y_score)
rows.append({"metric": name, "value": fn(y_test.to_numpy(), y_score),
"ci_low": low, "ci_high": high, "meaning": meaning})
report = pd.DataFrame(rows).set_index("metric").round(3)
print(f"Prevalence {y_test.mean():.3f} threshold {THRESHOLD:.3f} "
f"n={len(y_test)}")
print(report.to_string())
print(confusion_matrix(y_test, y_pred))
print(classification_report(y_test, y_pred, zero_division=0,
target_names=["No readmission", "Readmission"]))| Metric | Value | 95% CI |
|---|---|---|
| Accuracy | 0.828 | 0.806–0.849 |
| Precision | 0.341 | 0.282–0.397 |
| Recall | 0.700 | 0.617–0.776 |
| F1 | 0.458 | 0.394–0.516 |
| ROC-AUC | 0.810 | 0.760–0.854 |
| PR-AUC | 0.487 | 0.398–0.573 |
| Brier | 0.071 | 0.060–0.082 |
Production checklist
Lock the evaluation dataset, label definition, prediction time, and exclusion criteria before testing.
Use a temporally held-out test set for evolving clinical populations. If the model must generalize to unseen hospitals, hold out whole sites; outcome stratification alone does not test that.
Choose the operating threshold from clinical cost, alert capacity, and patient safety, on data other than the test set; never rely on predict()’s 0.5 default.
Report every metric with a confidence interval, the confusion matrix, prevalence, and calibration; check subgroup performance for inequities.
Monitor data drift, label delay, score distribution, alert volume, and clinician overrides, and revalidate after changes to EHR workflows, coding practices, patient mix, or treatment protocols.
Takeaway: For rare clinical events, report precision, recall, and F1 at a clinically chosen threshold alongside PR-AUC, ROC-AUC, and calibration — each with its uncertainty — rather than accuracy alone.