Add calibration metrics and patient-clustered bootstrap CIs

docs/api/metrics.rst

@@ -18,6 +18,7 @@ fidelity metrics.
     metrics/pyhealth.metrics.multilabel
     metrics/pyhealth.metrics.binary
     metrics/pyhealth.metrics.calibration
+    metrics/pyhealth.metrics.bootstrap
     metrics/pyhealth.metrics.prediction_set
     metrics/pyhealth.metrics.fairness
     metrics/pyhealth.metrics.interpretability

docs/api/metrics/pyhealth.metrics.bootstrap.rst

@@ -0,0 +1,31 @@
+pyhealth.metrics.bootstrap
+===================================
+
+Confidence intervals for evaluation metrics, as clinical reporting guidelines
+such as TRIPOD+AI expect. Samples from one patient are correlated, so pass
+``groups=patient_ids`` to resample whole patients (a cluster bootstrap) rather
+than individual samples. To compare two models, use ``paired_bootstrap_diff``:
+both are scored on identical resamples, so its interval is for their difference.
+
+.. code-block:: python
+
+    from pyhealth.metrics import bootstrap_ci, paired_bootstrap_diff
+
+    ci = bootstrap_ci(y_true, y_prob, "pr_auc", groups=patient_ids, n_boot=1000)
+    # {'estimate': 0.41, 'lower': 0.35, 'upper': 0.47, 'n_boot': 1000, 'n_skipped': 0}
+
+    diff = paired_bootstrap_diff(y_true, y_prob_new, y_prob_baseline, "pr_auc",
+                                 groups=patient_ids)
+    # an interval that excludes 0 indicates a difference at the 5% level
+
+``metric`` is any name accepted by ``binary_metrics_fn``, including the
+calibration metrics ``brier``, ``oe_ratio``, ``calibration_slope`` and
+``calibration_intercept``, or a callable ``metric(y_true, y_prob) -> float``.
+Resamples that contain a single class are skipped and counted in ``n_skipped``.
+See ``examples/calibration_and_bootstrap_ci.py``.
+
+.. currentmodule:: pyhealth.metrics.bootstrap
+
+.. autofunction:: bootstrap_ci
+
+.. autofunction:: paired_bootstrap_diff

examples/calibration_and_bootstrap_ci.py

@@ -0,0 +1,53 @@
+"""Calibration metrics and patient-clustered bootstrap confidence intervals.
+
+Reports what a clinical prediction paper needs for each model: discrimination
+(PR-AUC, ROC-AUC) and calibration (Brier score, observed/expected ratio,
+calibration slope and intercept), each with a 95% interval from a bootstrap
+that resamples whole patients, plus the paired difference between two models.
+
+Uses synthetic predictions (several samples per patient); no model training
+or download needed. Runs in a few seconds.
+
+Usage:
+    python examples/calibration_and_bootstrap_ci.py
+"""
+
+import numpy as np
+
+from pyhealth.metrics import binary_metrics_fn, bootstrap_ci, paired_bootstrap_diff
+
+METRICS = ["pr_auc", "roc_auc", "brier", "oe_ratio", "calibration_slope", "calibration_intercept"]
+
+
+def synthetic_predictions(n_patients=300, seed=0):
+    rng = np.random.default_rng(seed)
+    visits = rng.integers(1, 6, size=n_patients)  # 1-5 samples per patient
+    patient_ids = np.repeat(np.arange(n_patients), visits)
+    risk = np.repeat(rng.beta(1, 6, size=n_patients), visits)  # rare outcome
+    y_true = (rng.uniform(size=risk.size) < risk).astype(int)
+    calibrated = np.clip(risk + rng.normal(0, 0.03, risk.size), 0.001, 0.999)
+    logit = np.log(calibrated / (1 - calibrated))
+    overconfident = 1 / (1 + np.exp(-1.8 * logit))  # same ranking, too extreme
+    return y_true, patient_ids, calibrated, overconfident
+
+
+def main():
+    y_true, patient_ids, model_a, model_b = synthetic_predictions()
+    print(f"{y_true.size} samples from {np.unique(patient_ids).size} patients, "
+          f"{y_true.mean():.1%} positive\n")
+
+    for name, y_prob in (("model A", model_a), ("model B (A made overconfident)", model_b)):
+        print(name)
+        point = binary_metrics_fn(y_true, y_prob, metrics=METRICS)
+        for metric in METRICS:
+            ci = bootstrap_ci(y_true, y_prob, metric, groups=patient_ids, n_boot=500)
+            print(f"  {metric:22s} {point[metric]:6.3f}  [{ci['lower']:.3f}, {ci['upper']:.3f}]")
+        print()
+
+    for metric in ("roc_auc", "brier"):
+        diff = paired_bootstrap_diff(y_true, model_a, model_b, metric, groups=patient_ids, n_boot=500)
+        print(f"A - B {metric:8s} {diff['estimate']:+.3f}  [{diff['lower']:+.3f}, {diff['upper']:+.3f}]")
+
+
+if __name__ == "__main__":
+    main()

pyhealth/metrics/__init__.py

@@ -1,43 +1,46 @@
-from .binary import binary_metrics_fn
-from .drug_recommendation import ddi_rate_score
-from .generative import (
-    calc_membership_inference,
-    calc_nnaar,
-    compute_discriminator_privacy,
-    compute_mle,
-    compute_prevalence_metrics,
-    evaluate_synthetic_ehr,
-)
-from .interpretability import (
-    ComprehensivenessMetric,
-    Evaluator,
-    RemovalBasedMetric,
-    SufficiencyMetric,
-    evaluate_attribution,
-)
-from .fairness import fairness_metrics_fn
-from .multiclass import multiclass_metrics_fn
-from .multilabel import multilabel_metrics_fn
-from .ranking import ranking_metrics_fn
-from .regression import regression_metrics_fn
-
-__all__ = [
-    "binary_metrics_fn",
-    "ddi_rate_score",
-    "calc_nnaar",
-    "calc_membership_inference",
-    "compute_discriminator_privacy",
-    "compute_mle",
-    "compute_prevalence_metrics",
-    "evaluate_synthetic_ehr",
-    "ComprehensivenessMetric",
-    "SufficiencyMetric",
-    "RemovalBasedMetric",
-    "Evaluator",
-    "evaluate_attribution",
-    "fairness_metrics_fn",
-    "multiclass_metrics_fn",
-    "multilabel_metrics_fn",
-    "ranking_metrics_fn",
-    "regression_metrics_fn",
-]
+from .binary import binary_metrics_fn
+from .bootstrap import bootstrap_ci, paired_bootstrap_diff
+from .drug_recommendation import ddi_rate_score
+from .fairness import fairness_metrics_fn
+from .generative import (
+    calc_membership_inference,
+    calc_nnaar,
+    compute_discriminator_privacy,
+    compute_mle,
+    compute_prevalence_metrics,
+    evaluate_synthetic_ehr,
+)
+from .interpretability import (
+    ComprehensivenessMetric,
+    Evaluator,
+    RemovalBasedMetric,
+    SufficiencyMetric,
+    evaluate_attribution,
+)
+from .multiclass import multiclass_metrics_fn
+from .multilabel import multilabel_metrics_fn
+from .ranking import ranking_metrics_fn
+from .regression import regression_metrics_fn
+
+__all__ = [
+    "ComprehensivenessMetric",
+    "Evaluator",
+    "RemovalBasedMetric",
+    "SufficiencyMetric",
+    "binary_metrics_fn",
+    "bootstrap_ci",
+    "calc_membership_inference",
+    "calc_nnaar",
+    "compute_discriminator_privacy",
+    "compute_mle",
+    "compute_prevalence_metrics",
+    "ddi_rate_score",
+    "evaluate_attribution",
+    "evaluate_synthetic_ehr",
+    "fairness_metrics_fn",
+    "multiclass_metrics_fn",
+    "multilabel_metrics_fn",
+    "paired_bootstrap_diff",
+    "ranking_metrics_fn",
+    "regression_metrics_fn",
+]

pyhealth/metrics/binary.py

@@ -6,6 +6,42 @@
 import pyhealth.metrics.calibration as calib
 
 
+def _logit(y_prob: np.ndarray) -> np.ndarray:
+    p = np.clip(np.asarray(y_prob, dtype=float), 1e-12, 1 - 1e-12)
+    return np.log(p / (1 - p))
+
+
+def _fit_logistic(x: np.ndarray, y: np.ndarray, offset: np.ndarray, iters: int = 100):
+    """Unpenalised logistic regression of y on columns x, plus an offset (Newton)."""
+    beta = np.zeros(x.shape[1])
+    for _ in range(iters):
+        mu = 1 / (1 + np.exp(-(x @ beta + offset)))
+        grad = x.T @ (y - mu)
+        hess = (x * (mu * (1 - mu))[:, None]).T @ x
+        try:
+            step = np.linalg.solve(hess, grad)
+        except np.linalg.LinAlgError:
+            return np.full(x.shape[1], np.nan)
+        beta = beta + step
+        if np.max(np.abs(step)) < 1e-10:
+            break
+    return beta
+
+
+def _calibration_slope(y_true: np.ndarray, y_prob: np.ndarray) -> float:
+    """Slope of logistic recalibration: y ~ a + b * logit(p). Ideal: 1."""
+    lp = _logit(y_prob)
+    x = np.column_stack([np.ones_like(lp), lp])
+    return float(_fit_logistic(x, np.asarray(y_true, dtype=float), np.zeros_like(lp))[1])
+
+
+def _calibration_intercept(y_true: np.ndarray, y_prob: np.ndarray) -> float:
+    """Calibration-in-the-large: a in y ~ a + offset(logit(p)). Ideal: 0."""
+    lp = _logit(y_prob)
+    x = np.ones((lp.size, 1))
+    return float(_fit_logistic(x, np.asarray(y_true, dtype=float), lp)[0])
+
+
 def binary_metrics_fn(
     y_true: np.ndarray,
     y_prob: np.ndarray,
@@ -28,6 +64,13 @@ def binary_metrics_fn(
         - jaccard: Jaccard similarity coefficient score
         - ECE: Expected Calibration Error (with 20 equal-width bins). Check :func:`pyhealth.metrics.calibration.ece_confidence_binary`.
         - ECE_adapt: adaptive ECE (with 20 equal-size bins). Check :func:`pyhealth.metrics.calibration.ece_confidence_binary`.
+        - brier: Brier score (mean squared error of the probabilities)
+        - oe_ratio: observed / expected events, sum(y_true) / sum(y_prob); ideal 1
+        - calibration_slope: slope b of the logistic recalibration
+          y ~ a + b * logit(y_prob); ideal 1, below 1 means overconfident
+        - calibration_intercept: calibration-in-the-large, a in
+          y ~ a + offset(logit(y_prob)); ideal 0, above 0 means risks are
+          underestimated
     If no metrics are specified, pr_auc, roc_auc and f1 are computed by default.
 
     This function calls sklearn.metrics functions to compute the metrics. For
@@ -87,6 +130,14 @@ def binary_metrics_fn(
         elif metric == "jaccard":
             jaccard = sklearn_metrics.jaccard_score(y_true, y_pred)
             output["jaccard"] = jaccard
+        elif metric == "brier":
+            output["brier"] = sklearn_metrics.brier_score_loss(y_true, y_prob)
+        elif metric == "oe_ratio":
+            output["oe_ratio"] = float(np.sum(y_true) / np.sum(y_prob))
+        elif metric == "calibration_slope":
+            output["calibration_slope"] = _calibration_slope(y_true, y_prob)
+        elif metric == "calibration_intercept":
+            output["calibration_intercept"] = _calibration_intercept(y_true, y_prob)
         elif metric in {"ECE", "ECE_adapt"}:
             output[metric] = calib.ece_confidence_binary(
                 y_prob, y_true, bins=20, adaptive=metric.endswith("_adapt")