diff --git a/MetricsReloaded/processes/overall_process.py b/MetricsReloaded/processes/overall_process.py index 9e945a0..ac2cdb8 100644 --- a/MetricsReloaded/processes/overall_process.py +++ b/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: diff --git a/MetricsReloaded/utility/uncertainty.py b/MetricsReloaded/utility/uncertainty.py new file mode 100644 index 0000000..ed56618 --- /dev/null +++ b/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 diff --git a/test/test_utility/test_uncertainty.py b/test/test_utility/test_uncertainty.py new file mode 100644 index 0000000..c2f1ade --- /dev/null +++ b/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()