Add percentile-bootstrap confidence intervals for aggregated metrics

MetricsReloaded/processes/overall_process.py

@@ -184,6 +184,7 @@
 from MetricsReloaded.processes.mixed_measures_processes import MultiLabelLocMeasures, MultiLabelPairwiseMeasures, MultiLabelLocSegPairwiseMeasure
 import warnings
 from MetricsReloaded.utility.utils import combine_df, merge_list_df
+from MetricsReloaded.utility.uncertainty import stats_with_ci
 import pandas as pd
 import numpy as np
 
@@ -682,9 +683,13 @@ def label_aggregation(self, option='average',dict_args={}):
 
     def get_stats_res(self):
         """
-        Create summary statistics overall and per label available in self.stats_lab and self.stats_all
+        Create summary statistics overall and per label available in self.stats_lab and self.stats_all.
+        The overall statistics include percentile-bootstrap confidence
+        intervals (rows ci95_low / ci95_high) for each metric's mean, so that
+        aggregated results are reported with their uncertainty rather than as
+        bare point estimates.
         """
-        df_stats_all = self.grouped_lab.describe()
+        df_stats_all = stats_with_ci(self.grouped_lab, exclude=("index", "case"))
         self.stats_all = df_stats_all
         print(self.resdet, self.resseg)
         if len(self.resdet.index)==0 and len(self.resseg.index)==0:

MetricsReloaded/utility/uncertainty.py

@@ -0,0 +1,110 @@
+"""
+Uncertainty quantification - :mod:`MetricsReloaded.utility.uncertainty`
+=======================================================================
+
+This module provides functions for reporting the uncertainty of aggregated
+metric values.
+
+A mean metric over N cases is a point estimate: a mean DSC of 0.85 over 20
+cases and over 2,000 cases support very different conclusions, and validation
+studies routinely compare methods whose intervals overlap entirely. The
+Metrics Reloaded recommendations call for reporting variability alongside
+aggregates; this module provides the standard non-parametric tool for it.
+
+.. currentmodule:: MetricsReloaded.utility.uncertainty
+
+.. autosummary::
+    :nosignatures:
+
+    percentile_bootstrap_ci
+    stats_with_ci
+
+"""
+
+import numpy as np
+import pandas as pd
+
+__all__ = [
+    "percentile_bootstrap_ci",
+    "stats_with_ci",
+]
+
+
+def percentile_bootstrap_ci(values, n_boot=2000, alpha=0.05, seed=42):
+    """Percentile-bootstrap confidence interval for the mean of per-case
+    metric values.
+
+    NaN values are ignored, consistent with the masked aggregations used
+    elsewhere in this package. The resampling is seeded by default so that
+    reported intervals are reproducible across runs of the same evaluation.
+
+    :param values: iterable of per-case metric values (may contain NaN)
+    :param n_boot: number of bootstrap resamples; must be at least 2, and in
+        practice should be far larger. Both bounds are quantiles of the
+        resample distribution, so while n_boot < 2 / alpha they are decided by
+        its most extreme draws alone (fewer than 40 resamples at the default
+        alpha=0.05) and the interval is too narrow rather than merely noisy.
+    :param alpha: 1 - confidence level (0.05 -> 95% interval); must be
+        strictly between 0 and 1
+    :param seed: seed for the resampling generator; pass None for
+        non-deterministic resampling
+    :return: (lower, upper) bounds of the interval; (nan, nan) when fewer
+        than two non-NaN values are available
+    """
+    # n_boot=1 passes any "is it positive?" check and then returns a
+    # zero-width interval, because both quantiles of a single resample are
+    # that resample: a claim of exact precision produced by the least
+    # evidence the function accepts.
+    if n_boot < 2:
+        raise ValueError("n_boot must be at least 2, got %r" % (n_boot,))
+    if not 0.0 < alpha < 1.0:
+        raise ValueError("alpha must be strictly between 0 and 1, got %r" % (alpha,))
+    arr = np.asarray(values, dtype=float)
+    arr = arr[~np.isnan(arr)]
+    if arr.size < 2:
+        return (np.nan, np.nan)
+    rng = np.random.default_rng(seed)
+    indices = rng.integers(0, arr.size, size=(n_boot, arr.size))
+    boot_means = arr[indices].mean(axis=1)
+    return (
+        float(np.quantile(boot_means, alpha / 2)),
+        float(np.quantile(boot_means, 1 - alpha / 2)),
+    )
+
+
+def stats_with_ci(df, n_boot=2000, alpha=0.05, seed=42, exclude=()):
+    """Summary statistics of a per-case results dataframe with
+    confidence-interval rows for each numeric metric column.
+
+    Returns the ``df.describe()`` table augmented with two rows,
+    ``ci95_low`` and ``ci95_high`` (names follow the requested ``alpha``;
+    fractional levels are preserved, e.g. ``alpha=0.025`` -> ``ci97.5_low``),
+    holding the percentile-bootstrap interval for each column's mean.
+
+    :param df: dataframe of per-case metric values (columns = metrics)
+    :param n_boot: number of bootstrap resamples per column
+    :param alpha: 1 - confidence level
+    :param seed: seed for reproducible intervals
+    :param exclude: column names to keep in the describe() table but skip
+        when bootstrapping (identifier columns such as ``case``); their CI
+        cells are left empty
+    :return: describe() dataframe with appended CI rows
+    """
+    described = df.describe()
+    level = (1 - alpha) * 100
+    low_name = "ci%g_low" % level
+    high_name = "ci%g_high" % level
+    lows = {}
+    highs = {}
+    for col in described.columns:
+        if col in exclude:
+            continue
+        series = pd.to_numeric(df[col], errors="coerce")
+        low, high = percentile_bootstrap_ci(
+            series.to_numpy(), n_boot=n_boot, alpha=alpha, seed=seed
+        )
+        lows[col] = low
+        highs[col] = high
+    described.loc[low_name] = pd.Series(lows)
+    described.loc[high_name] = pd.Series(highs)
+    return described

test/test_utility/test_uncertainty.py

@@ -0,0 +1,107 @@
+import numpy as np
+import pandas as pd
+import pytest
+from numpy.testing import assert_allclose
+
+from MetricsReloaded.utility.uncertainty import (
+    percentile_bootstrap_ci,
+    stats_with_ci,
+)
+
+
+def test_interval_brackets_the_mean_and_is_deterministic():
+    values = np.array([0.7, 0.8, 0.85, 0.9, 0.75, 0.82, 0.88, 0.79, 0.81, 0.86])
+    low, high = percentile_bootstrap_ci(values)
+    assert low < values.mean() < high
+    # seeded resampling: identical call gives identical interval
+    assert (low, high) == percentile_bootstrap_ci(values)
+
+
+def test_interval_narrows_with_more_cases():
+    rng = np.random.default_rng(0)
+    small = rng.normal(0.8, 0.1, size=20)
+    large = np.concatenate([small] * 50)
+    low_s, high_s = percentile_bootstrap_ci(small)
+    low_l, high_l = percentile_bootstrap_ci(large)
+    assert (high_l - low_l) < (high_s - low_s)
+
+
+def test_approximates_analytic_interval_for_gaussian_data():
+    rng = np.random.default_rng(1)
+    values = rng.normal(0.5, 0.2, size=400)
+    low, high = percentile_bootstrap_ci(values, n_boot=5000)
+    sem = values.std(ddof=1) / np.sqrt(values.size)
+    assert_allclose(low, values.mean() - 1.96 * sem, atol=3 * sem / 10)
+    assert_allclose(high, values.mean() + 1.96 * sem, atol=3 * sem / 10)
+
+
+def test_nan_values_are_ignored():
+    values = [0.8, np.nan, 0.9, 0.85, np.nan, 0.82]
+    low, high = percentile_bootstrap_ci(values)
+    clean_low, clean_high = percentile_bootstrap_ci([0.8, 0.9, 0.85, 0.82])
+    assert (low, high) == (clean_low, clean_high)
+
+
+def test_degenerate_inputs_return_nan_interval():
+    assert np.isnan(percentile_bootstrap_ci([])[0])
+    assert np.isnan(percentile_bootstrap_ci([0.5])[1])
+    assert np.isnan(percentile_bootstrap_ci([np.nan, np.nan])[0])
+
+
+def test_stats_with_ci_appends_interval_rows():
+    df = pd.DataFrame(
+        {
+            "dsc": [0.7, 0.8, 0.85, 0.9, 0.75],
+            "nsd": [0.6, 0.65, 0.7, 0.72, 0.68],
+        }
+    )
+    stats = stats_with_ci(df)
+    assert "ci95_low" in stats.index
+    assert "ci95_high" in stats.index
+    for col in ("dsc", "nsd"):
+        assert stats.loc["ci95_low", col] < df[col].mean()
+        assert stats.loc["ci95_high", col] > df[col].mean()
+    # describe() content preserved
+    assert_allclose(stats.loc["mean", "dsc"], df["dsc"].mean())
+
+
+def test_invalid_options_are_rejected():
+    values = [0.7, 0.8, 0.9]
+    with pytest.raises(ValueError):
+        percentile_bootstrap_ci(values, n_boot=0)
+    with pytest.raises(ValueError):
+        percentile_bootstrap_ci(values, alpha=1.1)
+    with pytest.raises(ValueError):
+        percentile_bootstrap_ci(values, alpha=0.0)
+
+
+def test_single_resample_is_rejected_rather_than_reported_as_exact():
+    # Both bounds are quantiles of the resample distribution, so with one
+    # resample they are the same number and the interval has zero width --
+    # the strongest possible claim from the weakest possible evidence.
+    with pytest.raises(ValueError):
+        percentile_bootstrap_ci([0.1, 0.5, 0.9, 0.3, 0.7], n_boot=1)
+
+
+def test_fractional_confidence_level_is_labelled_exactly():
+    df = pd.DataFrame({"dsc": [0.7, 0.8, 0.85, 0.9, 0.75]})
+    stats = stats_with_ci(df, alpha=0.025)
+    # a 97.5% interval must not be rounded to "ci98"
+    assert "ci97.5_low" in stats.index
+    assert "ci97.5_high" in stats.index
+
+
+def test_identifier_columns_are_excluded_from_bootstrapping():
+    df = pd.DataFrame(
+        {
+            "case": [1, 2, 3, 4, 5],
+            "dsc": [0.7, 0.8, 0.85, 0.9, 0.75],
+        }
+    )
+    stats = stats_with_ci(df, exclude=("case",))
+    # the identifier keeps its describe() rows but gets no interval
+    assert_allclose(stats.loc["mean", "case"], 3.0)
+    assert np.isnan(stats.loc["ci95_low", "case"])
+    assert np.isnan(stats.loc["ci95_high", "case"])
+    # metric columns still get one
+    assert stats.loc["ci95_low", "dsc"] < df["dsc"].mean()