feat: add binary classification and proper-scoring alert metrics

chap_core/assessment/metrics/__init__.py

@@ -19,6 +19,7 @@
     DEFAULT_OUTPUT_DIMENSIONS,
     AggregationOp,
     DeterministicMetric,
+    GlobalOnlyMetric,
     Metric,
     MetricSpec,
     OptimizationDirection,
@@ -86,7 +87,9 @@ def _discover_metrics():
         example_metric,
         mae,
         mape,
+        outbreak_classification,
         outbreak_detection,
+        outbreak_probability,
         percentile_coverage,
         rmse,
         test_metrics,
@@ -104,11 +107,22 @@ def _discover_metrics():
 from chap_core.assessment.metrics.example_metric import ExampleMetric
 from chap_core.assessment.metrics.mae import MAEMetric
 from chap_core.assessment.metrics.mape import MAPEMetric
+from chap_core.assessment.metrics.outbreak_classification import (
+    FalseAlarmRateMetric,
+    MatthewsCorrelationMetric,
+    OutbreakF1Metric,
+    OutbreakPrecisionMetric,
+)
 from chap_core.assessment.metrics.outbreak_detection import (
     OutbreakAccuracyMetric,
     SensitivityMetric,
     SpecificityMetric,
 )
+from chap_core.assessment.metrics.outbreak_probability import (
+    BrierScoreMetric,
+    BrierSkillScoreMetric,
+    OutbreakLogScoreMetric,
+)
 from chap_core.assessment.metrics.peak_diff import PeakPeriodLagMetric, PeakValueDiffMetric
 from chap_core.assessment.metrics.percentile_coverage import (
     Coverage10_90Metric,
@@ -132,6 +146,8 @@ def _discover_metrics():
 __all__ = [
     "DEFAULT_OUTPUT_DIMENSIONS",
     "AggregationOp",
+    "BrierScoreMetric",
+    "BrierSkillScoreMetric",
     "CRPSLog1pMetric",
     "CRPSMetric",
     "CRPSNormMetric",
@@ -140,11 +156,17 @@ def _discover_metrics():
     "DataDimension",
     "DeterministicMetric",
     "ExampleMetric",
+    "FalseAlarmRateMetric",
+    "GlobalOnlyMetric",
     "MAEMetric",
     "MAPEMetric",
+    "MatthewsCorrelationMetric",
     "Metric",
     "MetricSpec",
     "OutbreakAccuracyMetric",
+    "OutbreakF1Metric",
+    "OutbreakLogScoreMetric",
+    "OutbreakPrecisionMetric",
     "PeakPeriodLagMetric",
     "PeakValueDiffMetric",
     "PercentileCoverageMetric",
@@ -226,6 +248,9 @@ def compute_all_detailed_metrics(evaluation: Evaluation) -> pd.DataFrame:
         metric = metric_factory(historical_observations=historical_df)
         if not metric.is_applicable(flat_data.observations):
             continue
+        if not metric.spec.output_dimensions:
+            # Global-only metric (F1, skill scores): no per-cell value to export.
+            continue
         try:
             detailed = metric.get_detailed_metric(flat_data.observations, flat_data.forecasts)
         except Exception:

chap_core/assessment/metrics/base.py

@@ -235,6 +235,44 @@ def get_description(self) -> str:
         return self.spec.description
 
 
+class GlobalOnlyMetric(Metric):
+    """Base for metrics that cannot be computed per cell and then averaged.
+
+    F1, Matthews correlation and skill scores are ratios of aggregates: there is
+    no per-cell value whose mean recovers the metric. Such a metric is defined
+    only over a whole set of scored cells, so it reports one global value and
+    refuses to be broken down by dimension rather than returning a number that
+    looks per-cell but is not.
+
+    Subclasses implement :meth:`compute_global` and declare
+    ``output_dimensions=()`` in their spec.
+    """
+
+    def get_metric(
+        self,
+        observations: pa.typing.DataFrame[FlatObserved],
+        forecasts: pa.typing.DataFrame[FlatForecasts],
+        dimensions: tuple[DataDimension, ...] = (),
+    ) -> pd.DataFrame:
+        if dimensions:
+            named = ", ".join(d.value for d in dimensions)
+            raise ValueError(
+                f"{type(self).__name__} is defined over a whole set of scored cells, so it cannot be broken down by {named}."
+            )
+        null_mask = observations.disease_cases.isnull()  # type: ignore[attr-defined]
+        observations = observations[~null_mask]  # type: ignore[index]
+        value = self.compute_global(observations, forecasts)  # type: ignore[arg-type]
+        return self._validate_output(pd.DataFrame({"metric": [value]}), ())
+
+    @abstractmethod
+    def compute_global(self, observations: pd.DataFrame, forecasts: pd.DataFrame) -> float:
+        """Compute the metric over every scored cell at once."""
+        raise NotImplementedError
+
+    def compute_detailed(self, observations: pd.DataFrame, forecasts: pd.DataFrame) -> pd.DataFrame:
+        raise NotImplementedError(f"{type(self).__name__} has no per-cell value.")
+
+
 class DeterministicMetric(Metric):
     """
     Base class for deterministic metrics that operate on the median of samples.

chap_core/assessment/metrics/outbreak_classification.py

@@ -0,0 +1,143 @@
+"""Binary classification metrics for outbreak alerts.
+
+These score the alert decision, after the exceedance probability has been cut at
+:data:`~chap_core.assessment.metrics.outbreak_detection.ALERT_SAMPLE_FRACTION`.
+
+Only F1 and Matthews correlation carry an ``optimization_direction``. Precision
+and false-alarm rate, like sensitivity and specificity before them, are each
+maximised by a degenerate model -- precision by alerting once and only when
+certain, false-alarm rate by never alerting at all -- so neither is valid as a
+standalone objective even though both are informative alongside the others.
+"""
+
+from __future__ import annotations
+
+from typing import TYPE_CHECKING
+
+import numpy as np
+
+from chap_core.assessment.metrics import metric
+from chap_core.assessment.metrics.base import (
+    AggregationOp,
+    GlobalOnlyMetric,
+    Metric,
+    MetricSpec,
+    OptimizationDirection,
+)
+from chap_core.assessment.metrics.outbreak_detection import OutbreakScoredMixin, _as_metric
+
+if TYPE_CHECKING:
+    import pandas as pd
+
+
+def _confusion(frame: pd.DataFrame) -> tuple[int, int, int, int]:
+    """Counts of (true positive, false positive, true negative, false negative)."""
+    alert = frame["alert"] == 1.0
+    outbreak = frame["outbreak"] == 1.0
+    return (
+        int((alert & outbreak).sum()),
+        int((alert & ~outbreak).sum()),
+        int((~alert & ~outbreak).sum()),
+        int((~alert & outbreak).sum()),
+    )
+
+
+@metric()
+class OutbreakPrecisionMetric(OutbreakScoredMixin, Metric):
+    """Precision: the share of raised alerts that turned out to be real outbreaks.
+
+    The counterpart to sensitivity. A pipeline that alerts constantly scores well
+    on sensitivity and badly here; one that alerts almost never does the reverse.
+    """
+
+    spec = MetricSpec(
+        metric_id="outbreak_precision",
+        metric_name="Outbreak Precision",
+        aggregation_op=AggregationOp.MEAN,
+        description="Share of raised alerts that were real outbreaks",
+    )
+
+    def compute_detailed(self, observations: pd.DataFrame, forecasts: pd.DataFrame) -> pd.DataFrame:
+        frame = self.outbreak_frame(observations, forecasts)
+        alerted = frame[frame["alert"] == 1.0]
+        return _as_metric(alerted, alerted["outbreak"])
+
+
+@metric()
+class FalseAlarmRateMetric(OutbreakScoredMixin, Metric):
+    """False-alarm rate: the share of quiet periods that were alerted anyway.
+
+    One minus specificity, reported directly because it is the number an
+    operational team feels -- how often a warning turns out to be nothing.
+    """
+
+    spec = MetricSpec(
+        metric_id="false_alarm_rate",
+        metric_name="False Alarm Rate",
+        aggregation_op=AggregationOp.MEAN,
+        description="Share of non-outbreak periods that raised an alert",
+    )
+
+    def compute_detailed(self, observations: pd.DataFrame, forecasts: pd.DataFrame) -> pd.DataFrame:
+        frame = self.outbreak_frame(observations, forecasts)
+        quiet = frame[frame["outbreak"] == 0.0]
+        return _as_metric(quiet, quiet["alert"])
+
+
+@metric()
+class OutbreakF1Metric(OutbreakScoredMixin, GlobalOnlyMetric):
+    """F1: the harmonic mean of precision and sensitivity.
+
+    A ratio of aggregates, so it has no per-cell value. Unlike its two parts it
+    prices both error types, which makes it usable as an objective.
+    """
+
+    spec = MetricSpec(
+        metric_id="outbreak_f1",
+        metric_name="Outbreak F1",
+        output_dimensions=(),
+        aggregation_op=AggregationOp.MEAN,
+        description="Harmonic mean of outbreak precision and sensitivity",
+        optimization_direction=OptimizationDirection.MAXIMIZE,
+    )
+
+    def compute_global(self, observations: pd.DataFrame, forecasts: pd.DataFrame) -> float:
+        true_positive, false_positive, _, false_negative = _confusion(self.outbreak_frame(observations, forecasts))
+        denominator = 2 * true_positive + false_positive + false_negative
+        if denominator == 0:
+            return float("nan")
+        return 2 * true_positive / denominator
+
+
+@metric()
+class MatthewsCorrelationMetric(OutbreakScoredMixin, GlobalOnlyMetric):
+    """Matthews correlation between alerts and outbreaks, in [-1, 1].
+
+    Uses all four cells of the confusion matrix, so unlike F1 it does not ignore
+    true negatives -- worth having when outbreaks are rare and quiet periods
+    dominate. Zero means no better than chance.
+    """
+
+    spec = MetricSpec(
+        metric_id="outbreak_mcc",
+        metric_name="Outbreak Matthews Correlation",
+        output_dimensions=(),
+        aggregation_op=AggregationOp.MEAN,
+        description="Correlation between alert and outbreak status across the full confusion matrix",
+        optimization_direction=OptimizationDirection.MAXIMIZE,
+    )
+
+    def compute_global(self, observations: pd.DataFrame, forecasts: pd.DataFrame) -> float:
+        true_positive, false_positive, true_negative, false_negative = _confusion(
+            self.outbreak_frame(observations, forecasts)
+        )
+        denominator = np.sqrt(
+            float(true_positive + false_positive)
+            * float(true_positive + false_negative)
+            * float(true_negative + false_positive)
+            * float(true_negative + false_negative)
+        )
+        if denominator == 0:
+            # A row or column of the matrix is empty, so correlation is undefined.
+            return float("nan")
+        return float((true_positive * true_negative - false_positive * false_negative) / denominator)

chap_core/assessment/metrics/outbreak_detection.py

@@ -35,7 +35,7 @@
 #: Fraction of forecast samples that must exceed the threshold for an alert to be raised.
 ALERT_SAMPLE_FRACTION = 0.5
 
-_OUTBREAK_COLUMNS = ["location", "time_period", "horizon_distance", "outbreak", "alert"]
+_OUTBREAK_COLUMNS = ["location", "time_period", "horizon_distance", "outbreak", "alert", "probability"]
 _METRIC_DIMENSIONS = ["location", "time_period", "horizon_distance"]
 
 
@@ -75,9 +75,10 @@ def outbreak_and_alert(
     Returns:
         One row per ``(location, time_period, horizon_distance)`` that has both a
         computable threshold and a forecast, with columns
-        ``[location, time_period, horizon_distance, outbreak, alert]``. ``outbreak``
-        is 1.0 where observed cases exceed the channel; ``alert`` is 1.0 where more
-        than :data:`ALERT_SAMPLE_FRACTION` of the samples exceed it. Cells whose
+        ``[location, time_period, horizon_distance, outbreak, alert, probability]``.
+        ``outbreak`` is 1.0 where observed cases exceed the channel; ``probability``
+        is the share of samples above it and ``alert`` is 1.0 where that share
+        exceeds :data:`ALERT_SAMPLE_FRACTION`. Cells whose
         threshold cannot be computed (a single historical value gives an undefined
         standard deviation) are dropped.
     """
@@ -108,7 +109,7 @@ def outbreak_and_alert(
     labelled["outbreak"] = (labelled["disease_cases"] > labelled["threshold"]).astype(float)
 
     merged = labelled[["location", "time_period", "outbreak"]].merge(
-        probabilities[[*_METRIC_DIMENSIONS, "alert"]],
+        probabilities[[*_METRIC_DIMENSIONS, "alert", "probability"]],
         on=["location", "time_period"],
         how="inner",
     )
@@ -122,18 +123,25 @@ def _as_metric(rows: pd.DataFrame, values: pd.Series) -> pd.DataFrame:
     return result
 
 
-class _OutbreakMetric(Metric):
-    """Shared applicability rule for metrics scored against a seasonal threshold."""
+class OutbreakScoredMixin:
+    """Applicability rule and frame access shared by every metric scored against the channel.
+
+    Mixed in ahead of :class:`Metric` or :class:`GlobalOnlyMetric` depending on
+    whether the metric has a per-cell value.
+    """
+
+    historical_observations: pd.DataFrame | None
 
     def is_applicable(self, observations: pa.typing.DataFrame[FlatObserved]) -> bool:
         return self.historical_observations is not None and has_season_buckets(observations)
 
-    def _frame(self, observations: pd.DataFrame, forecasts: pd.DataFrame) -> pd.DataFrame:
+    def outbreak_frame(self, observations: pd.DataFrame, forecasts: pd.DataFrame) -> pd.DataFrame:
+        """Labelled and scored cells: see :func:`outbreak_and_alert`."""
         return outbreak_and_alert(self.historical_observations, observations, forecasts)
 
 
 @metric()
-class SensitivityMetric(_OutbreakMetric):
+class SensitivityMetric(OutbreakScoredMixin, Metric):
     """Sensitivity (true positive rate) for outbreak detection.
 
     Measures the proportion of actual outbreaks that were correctly
@@ -151,13 +159,13 @@ class SensitivityMetric(_OutbreakMetric):
     )
 
     def compute_detailed(self, observations: pd.DataFrame, forecasts: pd.DataFrame) -> pd.DataFrame:
-        frame = self._frame(observations, forecasts)
+        frame = self.outbreak_frame(observations, forecasts)
         outbreaks = frame[frame["outbreak"] == 1.0]
         return _as_metric(outbreaks, outbreaks["alert"])
 
 
 @metric()
-class SpecificityMetric(_OutbreakMetric):
+class SpecificityMetric(OutbreakScoredMixin, Metric):
     """Specificity (true negative rate) for outbreak detection.
 
     Measures the proportion of non-outbreak periods that were correctly
@@ -175,13 +183,13 @@ class SpecificityMetric(_OutbreakMetric):
     )
 
     def compute_detailed(self, observations: pd.DataFrame, forecasts: pd.DataFrame) -> pd.DataFrame:
-        frame = self._frame(observations, forecasts)
+        frame = self.outbreak_frame(observations, forecasts)
         quiet = frame[frame["outbreak"] == 0.0]
         return _as_metric(quiet, 1.0 - quiet["alert"])
 
 
 @metric()
-class OutbreakAccuracyMetric(_OutbreakMetric):
+class OutbreakAccuracyMetric(OutbreakScoredMixin, Metric):
     """Accuracy for outbreak detection.
 
     Measures the proportion of all periods where the alert status
@@ -200,5 +208,5 @@ class OutbreakAccuracyMetric(_OutbreakMetric):
     )
 
     def compute_detailed(self, observations: pd.DataFrame, forecasts: pd.DataFrame) -> pd.DataFrame:
-        frame = self._frame(observations, forecasts)
+        frame = self.outbreak_frame(observations, forecasts)
         return _as_metric(frame, (frame["alert"] == frame["outbreak"]).astype(float))