From 4b9ed9686296f133eb608e75c1cefa8aa59acf84 Mon Sep 17 00:00:00 2001 From: solarsys Date: Fri, 2 Oct 2026 20:58:13 -0400 Subject: [PATCH] Add calibration metrics and patient-clustered bootstrap CIs binary_metrics_fn had discrimination metrics and ECE, but none of the calibration measures TRIPOD+AI asks for, and PyHealth had no way to put confidence intervals on a metric or on the difference between models. - binary_metrics_fn: "brier" (sklearn), "oe_ratio" (observed/expected events), "calibration_slope" (b in y ~ a + b*logit(p)) and "calibration_intercept" (calibration-in-the-large: a with logit(p) as an offset). The logistic fits are a small unpenalised Newton solve in NumPy; no new dependency. - New pyhealth.metrics.bootstrap: bootstrap_ci and paired_bootstrap_diff. Percentile intervals; `groups=` resamples whole patients (cluster bootstrap); the paired difference scores both models on identical resamples; single-class resamples are skipped and counted; results are deterministic given `seed`. Metric = a binary_metrics_fn name or a callable. Exported from pyhealth.metrics. - tests/core/test_calibration_bootstrap.py: brier vs sklearn, O/E, a calibrated model (slope ~1, intercept ~0, O/E ~1), an overconfident one (slope ~0.5), an under-predicting one (intercept ~+1), slope equal to sklearn's unpenalised logistic regression to 4 places; bootstrap determinism, whole-patient resampling, skipped single-class resamples, identical resamples in the paired difference. - docs: pyhealth.metrics.bootstrap page; metric names in the docstring. - examples/calibration_and_bootstrap_ci.py. Co-Authored-By: Claude Opus 5.5 --- docs/api/metrics.rst | 1 + .../metrics/pyhealth.metrics.bootstrap.rst | 31 ++++ examples/calibration_and_bootstrap_ci.py | 53 ++++++ pyhealth/metrics/__init__.py | 89 +++++----- pyhealth/metrics/binary.py | 51 ++++++ pyhealth/metrics/bootstrap.py | 156 ++++++++++++++++++ tests/core/test_calibration_bootstrap.py | 124 ++++++++++++++ 7 files changed, 462 insertions(+), 43 deletions(-) create mode 100644 docs/api/metrics/pyhealth.metrics.bootstrap.rst create mode 100644 examples/calibration_and_bootstrap_ci.py create mode 100644 pyhealth/metrics/bootstrap.py create mode 100644 tests/core/test_calibration_bootstrap.py diff --git a/docs/api/metrics.rst b/docs/api/metrics.rst index 9e6bc160a..3a036453c 100644 --- a/docs/api/metrics.rst +++ b/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 diff --git a/docs/api/metrics/pyhealth.metrics.bootstrap.rst b/docs/api/metrics/pyhealth.metrics.bootstrap.rst new file mode 100644 index 000000000..8317d0d88 --- /dev/null +++ b/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 diff --git a/examples/calibration_and_bootstrap_ci.py b/examples/calibration_and_bootstrap_ci.py new file mode 100644 index 000000000..d1371bbd5 --- /dev/null +++ b/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() diff --git a/pyhealth/metrics/__init__.py b/pyhealth/metrics/__init__.py index 42d922c72..444bc8212 100644 --- a/pyhealth/metrics/__init__.py +++ b/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", +] diff --git a/pyhealth/metrics/binary.py b/pyhealth/metrics/binary.py index ea7d125f7..ee340dc4c 100644 --- a/pyhealth/metrics/binary.py +++ b/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") diff --git a/pyhealth/metrics/bootstrap.py b/pyhealth/metrics/bootstrap.py new file mode 100644 index 000000000..520a10617 --- /dev/null +++ b/pyhealth/metrics/bootstrap.py @@ -0,0 +1,156 @@ +"""Bootstrap confidence intervals for evaluation metrics. + +Clinical evaluations usually report a metric with a confidence interval and, +when comparing models, the interval of their difference on the same data. +Samples from one patient are correlated, so resample whole patients +(``groups=patient_ids``) rather than individual samples. +""" + +from collections.abc import Callable, Iterator + +import numpy as np + +from .binary import binary_metrics_fn + +Metric = str | Callable[[np.ndarray, np.ndarray], float] + + +def _metric_fn(metric: Metric) -> Callable[[np.ndarray, np.ndarray], float]: + if callable(metric): + return metric + return lambda y_true, y_prob: binary_metrics_fn(y_true, y_prob, metrics=[metric])[ + metric + ] + + +def _resample_indices( + n: int, groups: np.ndarray | None, n_boot: int, seed: int +) -> Iterator[np.ndarray]: + """Yields n_boot index arrays, resampling samples or whole groups.""" + rng = np.random.default_rng(seed) + if groups is None: + for _ in range(n_boot): + yield rng.integers(0, n, size=n) + return + _, inverse, counts = np.unique(groups, return_inverse=True, return_counts=True) + order = np.argsort(inverse, kind="stable") + members = np.split(order, np.cumsum(counts)[:-1]) + for _ in range(n_boot): + picked = rng.integers(0, len(members), size=len(members)) + yield np.concatenate([members[i] for i in picked]) + + +def _summarise(estimate: float, stats: list, n_skipped: int, alpha: float) -> dict: + if stats: + lower, upper = np.quantile(stats, [alpha / 2, 1 - alpha / 2]) + else: + lower = upper = float("nan") + return { + "estimate": float(estimate), + "lower": float(lower), + "upper": float(upper), + "n_boot": len(stats), + "n_skipped": n_skipped, + } + + +def bootstrap_ci( + y_true: np.ndarray, + y_prob: np.ndarray, + metric: Metric, + groups: np.ndarray | None = None, + n_boot: int = 1000, + seed: int = 0, + alpha: float = 0.05, +) -> dict: + """Percentile bootstrap confidence interval for a binary metric. + + Args: + y_true: True binary labels, shape (n,). + y_prob: Predicted probabilities, shape (n,). + metric: A name accepted by ``binary_metrics_fn`` (e.g. ``"roc_auc"``, + ``"brier"``) or a callable ``metric(y_true, y_prob) -> float``. + groups: Optional cluster ids, shape (n,), e.g. patient ids. Whole + groups are resampled with replacement, keeping all their samples. + n_boot: Number of resamples drawn. + seed: Random seed; results are deterministic given it. + alpha: 1 - confidence level (0.05 gives a 95% interval). + + Returns: + dict with ``estimate`` (metric on all data), ``lower`` and ``upper`` + (percentile bounds), ``n_boot`` (resamples used) and ``n_skipped`` + (resamples with a single class, which are skipped). + + Examples: + >>> import numpy as np + >>> from pyhealth.metrics.bootstrap import bootstrap_ci + >>> y_true = np.array([0, 0, 1, 1, 0, 1, 0, 1]) + >>> y_prob = np.array([0.1, 0.3, 0.7, 0.8, 0.4, 0.6, 0.2, 0.9]) + >>> patients = np.array([1, 1, 2, 2, 3, 3, 4, 4]) + >>> ci = bootstrap_ci(y_true, y_prob, "brier", groups=patients, n_boot=200) + >>> sorted(ci) + ['estimate', 'lower', 'n_boot', 'n_skipped', 'upper'] + """ + y_true, y_prob = np.asarray(y_true), np.asarray(y_prob) + fn = _metric_fn(metric) + stats, n_skipped = [], 0 + for idx in _resample_indices(len(y_true), groups, n_boot, seed): + if np.unique(y_true[idx]).size < 2: + n_skipped += 1 + continue + stats.append(fn(y_true[idx], y_prob[idx])) + return _summarise(fn(y_true, y_prob), stats, n_skipped, alpha) + + +def paired_bootstrap_diff( + y_true: np.ndarray, + y_prob_a: np.ndarray, + y_prob_b: np.ndarray, + metric: Metric, + groups: np.ndarray | None = None, + n_boot: int = 1000, + seed: int = 0, + alpha: float = 0.05, +) -> dict: + """Bootstrap interval for metric(model A) - metric(model B) on the same data. + + Both models are scored on identical resamples, so the interval reflects + their difference rather than two independent uncertainties. An interval + that excludes 0 indicates a difference at level ``alpha``. + + Args: + y_true: True binary labels, shape (n,). + y_prob_a: Model A's predicted probabilities, shape (n,). + y_prob_b: Model B's predicted probabilities, shape (n,). + metric: A ``binary_metrics_fn`` name or a callable, as in + :func:`bootstrap_ci`. + groups: Optional cluster ids (e.g. patient ids) for a cluster bootstrap. + n_boot: Number of resamples drawn. + seed: Random seed; results are deterministic given it. + alpha: 1 - confidence level. + + Returns: + dict with ``estimate`` (A - B on all data), ``lower``, ``upper``, + ``n_boot`` and ``n_skipped``, as in :func:`bootstrap_ci`. + + Examples: + >>> import numpy as np + >>> from pyhealth.metrics.bootstrap import paired_bootstrap_diff + >>> y_true = np.array([0, 0, 1, 1, 0, 1, 0, 1]) + >>> model_a = np.array([0.1, 0.3, 0.7, 0.8, 0.4, 0.6, 0.2, 0.9]) + >>> model_b = np.array([0.4, 0.5, 0.5, 0.6, 0.5, 0.4, 0.3, 0.7]) + >>> diff = paired_bootstrap_diff(y_true, model_a, model_b, "roc_auc", n_boot=200) + >>> diff["estimate"] > 0 + True + """ + y_true = np.asarray(y_true) + y_prob_a, y_prob_b = np.asarray(y_prob_a), np.asarray(y_prob_b) + fn = _metric_fn(metric) + stats, n_skipped = [], 0 + for idx in _resample_indices(len(y_true), groups, n_boot, seed): + if np.unique(y_true[idx]).size < 2: + n_skipped += 1 + continue + stats.append(fn(y_true[idx], y_prob_a[idx]) - fn(y_true[idx], y_prob_b[idx])) + estimate = fn(y_true, y_prob_a) - fn(y_true, y_prob_b) + return _summarise(estimate, stats, n_skipped, alpha) diff --git a/tests/core/test_calibration_bootstrap.py b/tests/core/test_calibration_bootstrap.py new file mode 100644 index 000000000..489ba797c --- /dev/null +++ b/tests/core/test_calibration_bootstrap.py @@ -0,0 +1,124 @@ +"""Calibration metrics in binary_metrics_fn, and patient-clustered bootstrap CIs.""" + +import unittest + +import numpy as np +from sklearn.linear_model import LogisticRegression +from sklearn.metrics import brier_score_loss, roc_auc_score + +from pyhealth.metrics import binary_metrics_fn +from pyhealth.metrics.bootstrap import bootstrap_ci, paired_bootstrap_diff + + +def _logit(p): + return np.log(p / (1 - p)) + + +def _sigmoid(x): + return 1 / (1 + np.exp(-x)) + + +def _calibrated(n=50_000, seed=0): + rng = np.random.default_rng(seed) + p = rng.uniform(0.02, 0.98, size=n) + y = (rng.uniform(size=n) < p).astype(int) + return y, p + + +class TestCalibrationMetrics(unittest.TestCase): + def setUp(self): + self.y, self.p = _calibrated() + + def _metrics(self, p, names): + return binary_metrics_fn(self.y, p, metrics=names) + + def test_brier_matches_sklearn(self): + out = self._metrics(self.p, ["brier"]) + self.assertAlmostEqual(out["brier"], brier_score_loss(self.y, self.p)) + + def test_oe_ratio_is_observed_over_expected(self): + out = self._metrics(self.p, ["oe_ratio"]) + self.assertAlmostEqual(out["oe_ratio"], self.y.sum() / self.p.sum()) + + def test_calibrated_model(self): + out = self._metrics( + self.p, ["oe_ratio", "calibration_slope", "calibration_intercept"] + ) + self.assertAlmostEqual(out["oe_ratio"], 1.0, delta=0.02) + self.assertAlmostEqual(out["calibration_slope"], 1.0, delta=0.05) + self.assertAlmostEqual(out["calibration_intercept"], 0.0, delta=0.05) + + def test_overconfident_model_has_slope_below_one(self): + overconfident = _sigmoid(2 * _logit(self.p)) + out = self._metrics(overconfident, ["calibration_slope"]) + self.assertAlmostEqual(out["calibration_slope"], 0.5, delta=0.03) + + def test_underpredicting_model_has_positive_intercept(self): + too_low = _sigmoid(_logit(self.p) - 1.0) + out = self._metrics(too_low, ["calibration_intercept", "oe_ratio"]) + self.assertAlmostEqual(out["calibration_intercept"], 1.0, delta=0.05) + self.assertGreater(out["oe_ratio"], 1.0) + + def test_slope_matches_unpenalised_logistic_regression(self): + y, p = _calibrated(n=3_000, seed=1) + noisy = _sigmoid(1.3 * _logit(p) + 0.2) + lr = LogisticRegression(penalty=None, tol=1e-10, max_iter=1000) + lr.fit(_logit(noisy).reshape(-1, 1), y) + out = binary_metrics_fn(y, noisy, metrics=["calibration_slope"]) + self.assertAlmostEqual(out["calibration_slope"], lr.coef_[0][0], places=4) + + +class TestBootstrap(unittest.TestCase): + def setUp(self): + self.y, self.p = _calibrated(n=600, seed=2) + + def test_deterministic_and_contains_estimate(self): + a = bootstrap_ci(self.y, self.p, "roc_auc", n_boot=200, seed=7) + b = bootstrap_ci(self.y, self.p, "roc_auc", n_boot=200, seed=7) + self.assertEqual(a, b) + self.assertAlmostEqual(a["estimate"], roc_auc_score(self.y, self.p)) + self.assertLessEqual(a["lower"], a["estimate"]) + self.assertLessEqual(a["estimate"], a["upper"]) + self.assertEqual(a["n_boot"], 200) + + def test_callable_metric(self): + out = bootstrap_ci(self.y, self.p, roc_auc_score, n_boot=50, seed=0) + self.assertAlmostEqual(out["estimate"], roc_auc_score(self.y, self.p)) + + def test_groups_resample_whole_patients(self): + # Two samples per patient, each patient with its own probability. + y = np.repeat(self.y[:200], 2) + p = np.repeat(self.p[:200], 2) + groups = np.repeat(np.arange(200), 2) + + def every_patient_whole(y_res, p_res): + _, counts = np.unique(p_res, return_counts=True) + if np.any(counts % 2): + raise AssertionError("a patient's samples were split up") + return roc_auc_score(y_res, p_res) + + out = bootstrap_ci(y, p, every_patient_whole, groups=groups, n_boot=100, seed=0) + self.assertEqual(out["n_boot"], 100) + + def test_single_class_resamples_are_skipped_and_counted(self): + y = np.zeros(40, dtype=int) + y[0] = 1 + p = np.linspace(0.01, 0.99, 40) + out = bootstrap_ci(y, p, "roc_auc", n_boot=300, seed=0) + self.assertGreater(out["n_skipped"], 0) + self.assertEqual(out["n_boot"] + out["n_skipped"], 300) + + def test_paired_difference_uses_identical_resamples(self): + same = paired_bootstrap_diff(self.y, self.p, self.p, "roc_auc", n_boot=100, seed=0) + self.assertEqual((same["estimate"], same["lower"], same["upper"]), (0.0, 0.0, 0.0)) + + worse = np.clip(self.p + np.random.default_rng(3).normal(0, 0.3, self.p.size), 0.01, 0.99) + diff = paired_bootstrap_diff(self.y, self.p, worse, "roc_auc", n_boot=300, seed=0) + self.assertAlmostEqual( + diff["estimate"], roc_auc_score(self.y, self.p) - roc_auc_score(self.y, worse) + ) + self.assertGreater(diff["lower"], 0.0) # the clean model is clearly better + + +if __name__ == "__main__": + unittest.main()