Fix Exact-Horizon Censoring Semantics

CHANGELOG.md

@@ -2,6 +2,7 @@
 
 <!--next-version-placeholder-->
 
+- Fixed exact-horizon censoring semantics under `censoring_heuristic="excluded"`: subjects censored exactly at the evaluation horizon (`time == horizon`) are now retained as event-free rather than excluded.
 - Time-dependent interventions avoided is now derived from model-versus-treat-all net benefit, ensuring finite-sample consistency with decision-curve references and parity with dcurves.
 - Added opt-in canonical browser rendering for static Interventions Avoided using the verified `rtichoke_viz v0.7.0` release while preserving Plotly as the default.
 - Fixed Interventions Avoided to apply the per-100 scaling to the full model expression, including the false-negative penalty term.

src/rtichoke/calibration/calibration.py

@@ -156,7 +156,9 @@ def create_calibration_curve_times(
         List of prediction horizons (times) at which to evaluate calibration.
     heuristics_sets : List[Dict[str, str]]
         List of heuristic dictionaries defining censoring and competing risk
-        adjustments.
+        adjustments. Note: under ``censoring_heuristic="excluded"``, observations
+        censored exactly at the evaluation horizon (``time == horizon``) are
+        retained as event-free.
     calibration_type : str, optional
         Type of calibration plot, either ``"discrete"`` (binned) or ``"smooth"``.
         Defaults to ``"discrete"``.
@@ -1145,7 +1147,7 @@ def _apply_heuristics_and_censoring(
 
     # Heuristics for events before or at horizon
     if censoring_heuristic == "excluded":
-        df_adj = df_adj.filter(~((pl.col("real") == 0) & (pl.col("time") <= horizon)))
+        df_adj = df_adj.filter(~((pl.col("real") == 0) & (pl.col("time") < horizon)))
 
     if competing_heuristic == "excluded":
         df_adj = df_adj.filter(~((pl.col("real") == 2) & (pl.col("time") <= horizon)))

src/rtichoke/performance_data/performance_data_times.py

@@ -94,7 +94,9 @@ def prepare_performance_data_times(
         performance. Integer inputs are accepted and normalized to floats.
     heuristics_sets : list[Dict], optional
         A list of dictionaries, each specifying how to handle censored data
-        and competing events. The default is
+        and competing events. Under ``censoring_heuristic="excluded"``,
+        observations censored exactly at the evaluation horizon (``time == horizon``)
+        are retained as event-free. The default is
         ``[{"censoring_heuristic": "adjusted",
         "competing_heuristic": "adjusted_as_negative"}]``.
     stratified_by : Sequence[str], optional
@@ -191,7 +193,7 @@ def _compute_population_event_risk_times(
 
     if censoring_heuristic == "excluded":
         df = df.filter(
-            (pl.col("times") > pl.col("fixed_time_horizon")) | (pl.col("reals") > 0)
+            (pl.col("times") >= pl.col("fixed_time_horizon")) | (pl.col("reals") > 0)
         )
 
     if competing_heuristic == "excluded":

src/rtichoke/processing/adjustments.py

@@ -98,32 +98,27 @@ def _extract_excluded_events(
     censoring_heuristic: str,
     competing_heuristic: str,
 ) -> pl.DataFrame:
-    horizons_df = pl.DataFrame({"times": fixed_time_horizons}).sort("times")
-
-    excluded_events = horizons_df.join_asof(
-        event_table.with_columns(
-            pl.col("count_0").cum_sum().cast(pl.Float64).alias("real_censored_est"),
-            pl.col("count_2").cum_sum().cast(pl.Float64).alias("real_competing_est"),
-        ).select(
-            pl.col("times"),
-            pl.col("real_censored_est"),
-            pl.col("real_competing_est"),
-        ),
-        left_on="times",
-        right_on="times",
-    ).with_columns([pl.col("times").alias("fixed_time_horizon")])
-
-    if censoring_heuristic != "excluded":
-        excluded_events = excluded_events.with_columns(
-            pl.lit(0.0).alias("real_censored_est")
+    rows = []
+    for h in fixed_time_horizons:
+        c_count = (
+            float(event_table.filter(pl.col("times") < float(h))["count_0"].sum())
+            if censoring_heuristic == "excluded"
+            else 0.0
         )
-
-    if competing_heuristic != "excluded":
-        excluded_events = excluded_events.with_columns(
-            pl.lit(0.0).alias("real_competing_est")
+        comp_count = (
+            float(event_table.filter(pl.col("times") <= float(h))["count_2"].sum())
+            if competing_heuristic == "excluded"
+            else 0.0
+        )
+        rows.append(
+            {
+                "fixed_time_horizon": float(h),
+                "real_censored_est": c_count,
+                "real_competing_est": comp_count,
+            }
         )
 
-    return excluded_events
+    return pl.DataFrame(rows)
 
 
 def extract_crude_estimate_polars(data: pl.DataFrame) -> pl.DataFrame:
@@ -498,7 +493,7 @@ def create_adjusted_data(
 def _censored_count(df: pl.DataFrame) -> pl.DataFrame:
     return (
         df.with_columns(
-            ((pl.col("times") <= pl.col("fixed_time_horizon")) & (pl.col("reals") == 0))
+            ((pl.col("times") < pl.col("fixed_time_horizon")) & (pl.col("reals") == 0))
             .cast(pl.Float64)
             .alias("is_censored")
         )
@@ -610,7 +605,16 @@ def _aj_adjusted_events(
     # Special-case: both excluded (faster branch in original)
     if censoring == "excluded" and competing == "excluded":
         non_censored_non_competing = exploded.filter(
-            (pl.col("times") > pl.col("fixed_time_horizon")) | (pl.col("reals") == 1)
+            ~(
+                (
+                    (pl.col("reals") == 0)
+                    & (pl.col("times") < pl.col("fixed_time_horizon"))
+                )
+                | (
+                    (pl.col("reals") == 2)
+                    & (pl.col("times") <= pl.col("fixed_time_horizon"))
+                )
+            )
         )
 
         adjusted = _aj_estimates_per_horizon(
@@ -669,7 +673,7 @@ def _aj_adjusted_events(
     # - "excluded": remove administratively censored observations (use exploded with filter)
     base_df = (
         exploded.filter(
-            (pl.col("times") > pl.col("fixed_time_horizon")) | (pl.col("reals") > 0)
+            (pl.col("times") >= pl.col("fixed_time_horizon")) | (pl.col("reals") > 0)
         )
         if censoring == "excluded"
         else reference_group_data

tests/test_exact_horizon_censoring.py

@@ -0,0 +1,190 @@
+import numpy as np
+import pytest
+import polars as pl
+from polarstate import prepare_event_table
+
+from rtichoke.processing.adjustments import _extract_excluded_events
+from rtichoke.performance_data.performance_data_times import (
+    prepare_performance_data_times,
+    prepare_binned_classification_data_times,
+    _compute_population_event_risk_times,
+)
+
+
+@pytest.fixture
+def fixture_a_data():
+    """Exact eight-subject fixture (Fixture A).
+
+    ID | times | reals
+    A  | 2     | 1
+    B  | 4     | 2
+    C  | 5     | 0
+    D  | 7     | 0
+    E  | 10    | 2
+    F  | 10    | 1
+    G  | 12    | 0
+    H  | 14    | 1
+
+    0 = censored, 1 = primary event, 2 = competing event.
+    """
+    reals = np.array([1, 2, 0, 0, 2, 1, 0, 1])
+    times = np.array([2.0, 4.0, 5.0, 7.0, 10.0, 10.0, 12.0, 14.0])
+    probs = {"model_1": np.array([0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8])}
+    return probs, reals, times
+
+
+def test_extract_excluded_events_unit_and_edge_cases():
+    """Directly test _extract_excluded_events for correct horizon boundary predicates, unsorted horizons, early horizons, and tied event times."""
+    df = pl.DataFrame(
+        {
+            "reals": [1, 2, 0, 0, 2, 1],
+            "times": [2.0, 4.0, 5.0, 7.0, 10.0, 10.0],
+        }
+    )
+    event_table = prepare_event_table(df)
+
+    # Unsorted horizons including early horizon (1.0 < min(times))
+    horizons = [10.0, 1.0, 5.0]
+
+    res = _extract_excluded_events(event_table, horizons, "excluded", "excluded")
+
+    h10 = res.filter(pl.col("fixed_time_horizon") == 10.0).to_dicts()[0]
+    h1 = res.filter(pl.col("fixed_time_horizon") == 1.0).to_dicts()[0]
+    h5 = res.filter(pl.col("fixed_time_horizon") == 5.0).to_dicts()[0]
+
+    # At horizon 1.0: 0 censored (time < 1), 0 competing (time <= 1)
+    assert h1["real_censored_est"] == 0.0
+    assert h1["real_competing_est"] == 0.0
+
+    # At horizon 5.0:
+    # - Censored at t < 5.0: none (t=5.0 is NOT < 5.0) -> 0.0
+    # - Competing at t <= 5.0: t=4.0 -> 1.0
+    assert h5["real_censored_est"] == 0.0
+    assert h5["real_competing_est"] == 1.0
+
+    # At horizon 10.0:
+    # - Censored at t < 10.0: t=5.0 and t=7.0 -> 2.0
+    # - Competing at t <= 10.0: t=4.0 and t=10.0 -> 2.0
+    assert h10["real_censored_est"] == 2.0
+    assert h10["real_competing_est"] == 2.0
+
+
+def test_unconditional_population_event_risk_oracle(fixture_a_data):
+    """Directly assert unconditional full-population KM/AJ event risk estimates across all combinations and horizons."""
+    _, reals, times = fixture_a_data
+
+    # Expected unconditional KM/AJ event risk estimates P(T <= h, event = 1)
+    expected_unconditional_risk = {
+        # Under censoring="excluded", "excluded":
+        # At h=5: 7 retained subjects (A,B,C,D,E,F,G,H except none censored < 5, B excluded <= 5 -> N=7). Risk = 1/7 = 0.142857
+        # At h=10: 4 retained subjects (A,F,G,H). Positives at t<=10 are A(t=2), F(t=10) -> 2 events out of 4 = 2/4 = 0.50
+        ("excluded", "excluded", 5.0): 1 / 7,  # 0.142857
+        ("excluded", "excluded", 10.0): 2 / 4,  # 0.500000
+        ("excluded", "adjusted_as_negative", 5.0): 0.125,
+        ("excluded", "adjusted_as_negative", 10.0): 2 / 6,  # 0.333333
+        ("excluded", "adjusted_as_composite", 5.0): 0.25,
+        ("excluded", "adjusted_as_composite", 10.0): 4 / 6,  # 0.666667
+        ("excluded", "adjusted_as_censored", 5.0): 0.125,
+        ("excluded", "adjusted_as_censored", 10.0): 3
+        / 8,  # 0.375000 (KM with competing as censoring)
+        ("adjusted", "excluded", 5.0): 1 / 7,  # 0.142857
+        ("adjusted", "excluded", 10.0): 4 / 9,  # 0.444444
+        ("adjusted", "adjusted_as_negative", 5.0): 0.125,
+        ("adjusted", "adjusted_as_negative", 10.0): 5 / 16,  # 0.31250
+        ("adjusted", "adjusted_as_composite", 5.0): 0.25,
+        ("adjusted", "adjusted_as_composite", 10.0): 5 / 8,  # 0.625
+        ("adjusted", "adjusted_as_censored", 5.0): 0.125,
+        ("adjusted", "adjusted_as_censored", 10.0): 11 / 32,  # 0.34375
+    }
+
+    for (c, comp, h), exp_risk in expected_unconditional_risk.items():
+        risk = _compute_population_event_risk_times(
+            reals, times, horizon=h, censoring_heuristic=c, competing_heuristic=comp
+        )
+        assert pytest.approx(risk, abs=1e-5) == exp_risk, (
+            f"Unconditional risk mismatch for c={c}, comp={comp}, h={h}"
+        )
+
+
+def test_fixture_a_retained_subjects_and_performance_producer_outputs(fixture_a_data):
+    """Assert exact retained subject identities/counts and performance producer outputs for Fixture A."""
+    probs, reals, times = fixture_a_data
+
+    horizons = [5.0, 10.0]
+    censoring_options = ["excluded", "adjusted"]
+    competing_options = [
+        "excluded",
+        "adjusted_as_negative",
+        "adjusted_as_composite",
+        "adjusted_as_censored",
+    ]
+
+    heuristics_sets = [
+        {"censoring_heuristic": c, "competing_heuristic": comp}
+        for c in censoring_options
+        for comp in competing_options
+    ]
+
+    perf_df = prepare_performance_data_times(
+        probs=probs,
+        reals=reals,
+        times=times,
+        fixed_time_horizons=horizons,
+        heuristics_sets=heuristics_sets,
+        by=0.1,
+    )
+
+    binned_df = prepare_binned_classification_data_times(
+        probs=probs,
+        reals=reals,
+        times=times,
+        fixed_time_horizons=horizons,
+        heuristics_sets=heuristics_sets,
+        by=0.1,
+    )
+    assert binned_df.height > 0
+
+    # 1. Assert exact retained raw cohort size N under censoring="excluded"
+    # At h=5.0: Subject C (censored at t=5.0) is retained because time == 5.0 is NOT < 5.0. N = 8.
+    h5_exc_neg = perf_df.filter(
+        (pl.col("fixed_time_horizon") == 5.0)
+        & (pl.col("censoring_heuristic") == "excluded")
+        & (pl.col("competing_heuristic") == "adjusted_as_negative")
+        & (pl.col("chosen_cutoff") == 0.0)
+    )
+    assert h5_exc_neg["n"][0] == 8.0
+
+    # At h=10.0: Subjects C (t=5.0) and D (t=7.0) are censored at t < 10.0 and excluded. Retained N = 6.
+    h10_exc_neg = perf_df.filter(
+        (pl.col("fixed_time_horizon") == 10.0)
+        & (pl.col("censoring_heuristic") == "excluded")
+        & (pl.col("competing_heuristic") == "adjusted_as_negative")
+        & (pl.col("chosen_cutoff") == 0.0)
+    )
+    assert h10_exc_neg["n"][0] == 6.0
+
+    # At h=10.0 with both censoring="excluded" and competing="excluded":
+    # Excluded: C (t=5, c), D (t=7, c), B (t=4, comp), E (t=10, comp) -> 4 excluded.
+    # Retained: A (t=2, pos), F (t=10, pos), G (t=12, neg), H (t=14, neg) -> 4 retained.
+    h10_both_exc = perf_df.filter(
+        (pl.col("fixed_time_horizon") == 10.0)
+        & (pl.col("censoring_heuristic") == "excluded")
+        & (pl.col("competing_heuristic") == "excluded")
+        & (pl.col("chosen_cutoff") == 0.0)
+    )
+    assert h10_both_exc["n"][0] == 4.0
+    assert h10_both_exc["real_positives"][0] == 2.0  # A and F
+    assert h10_both_exc["real_negatives"][0] == 2.0  # G and H
+
+    # 2. Performance producer cutoff-aggregated expectations
+    # Under censoring_heuristic="adjusted", the performance producer pools within-stratum AJ state risk estimates
+    # over probability cutoff bins, yielding total cutoff-binned mass (e.g. n=13 across 10 threshold bins for single cohort).
+    h10_adj_excl = perf_df.filter(
+        (pl.col("fixed_time_horizon") == 10.0)
+        & (pl.col("censoring_heuristic") == "adjusted")
+        & (pl.col("competing_heuristic") == "excluded")
+        & (pl.col("chosen_cutoff") == 0.0)
+    )
+    assert h10_adj_excl["n"][0] == 13.0
+    assert pytest.approx(h10_adj_excl["real_positives"][0], abs=1e-4) == 5.142857
+    assert pytest.approx(h10_adj_excl["real_negatives"][0], abs=1e-4) == 7.857143

tests/test_heuristics.py

@@ -47,11 +47,13 @@ def test_censoring_excluded(sample_data):
     result = _apply_heuristics_and_censoring(
         sample_data, 10, "excluded", "adjusted_as_negative"
     )
-    # Excludes censored at times 2, 10. Admin censors time > 10. Competing at 3,9 -> 0.
+    # Excludes censored at time 2 (time < 10).
+    # Censored at time 10 (time == 10) is retained as event-free.
+    # Admin censors time > 10. Competing at 3, 9 -> 0.
     expected = pl.DataFrame(
         {
-            "real": [1, 0, 1, 0, 0],
-            "time": [1, 3, 8, 9, 12],
+            "real": [1, 0, 1, 0, 0, 0],
+            "time": [1, 3, 8, 9, 10, 12],
         }
     )
     assert_frame_equal(result.sort("time"), expected.sort("time"))