From cc1823084efbf357ce07c8599e12a9d70f75b029 Mon Sep 17 00:00:00 2001 From: "google-labs-jules[bot]" <161369871+google-labs-jules[bot]@users.noreply.github.com> Date: Sun, 6 Sep 2026 10:14:07 +0000 Subject: [PATCH 1/3] Fix time-dependent interventions avoided calculation (#412) Derive time-dependent interventions avoided directly from model and treat-all net benefit using population event risk from cutoff-zero. Fixes #412 Co-authored-by: uriahf <11351434+uriahf@users.noreply.github.com> --- CHANGELOG.md | 1 + .../performance_data_times.py | 63 +++ ...me_interventions_avoided_dcurves_parity.py | 428 ++++++++++++++++++ 3 files changed, 492 insertions(+) create mode 100644 tests/test_time_interventions_avoided_dcurves_parity.py diff --git a/CHANGELOG.md b/CHANGELOG.md index 6643a5c3..106adad2 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -2,6 +2,7 @@ +- 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. diff --git a/src/rtichoke/performance_data/performance_data_times.py b/src/rtichoke/performance_data/performance_data_times.py index e7590d4c..8132e88f 100644 --- a/src/rtichoke/performance_data/performance_data_times.py +++ b/src/rtichoke/performance_data/performance_data_times.py @@ -122,6 +122,7 @@ def prepare_performance_data_times( cumulative_aj_data = _calculate_cumulative_aj_data(final_adjusted_data) performance_data = _turn_cumulative_aj_to_performance_data(cumulative_aj_data) + performance_data = _recalculate_interventions_avoided_times(performance_data) group_order = {group: index for index, group in enumerate(probs)} horizon_order = { @@ -164,6 +165,68 @@ def prepare_performance_data_times( ) +def _recalculate_interventions_avoided_times( + performance_data: pl.DataFrame, +) -> pl.DataFrame: + """Recalculate time-dependent interventions avoided using model and treat-all net benefit. + + IA = 100 * (NB_model - NB_all) / [threshold / (1 - threshold)] + + Population event risk is extracted at chosen_cutoff == 0 where stratified_by == "probability_threshold". + Interventions avoided is calculated for probability_threshold rows where 0 < chosen_cutoff < 1. + For chosen_cutoff == 0 or 1, and for PPCR rows, interventions avoided is set to null. + """ + event_risk_df = ( + performance_data.filter( + (pl.col("chosen_cutoff") == 0) + & (pl.col("stratified_by") == "probability_threshold") + ) + .select( + "reference_group", + "fixed_time_horizon", + "censoring_heuristic", + "competing_heuristic", + (pl.col("real_positives") / pl.col("n")).alias("_event_risk"), + ) + .unique() + ) + + performance_data = performance_data.join( + event_risk_df, + on=[ + "reference_group", + "fixed_time_horizon", + "censoring_heuristic", + "competing_heuristic", + ], + how="left", + ) + + threshold_odds = pl.col("chosen_cutoff") / (1 - pl.col("chosen_cutoff")) + net_benefit_all = ( + pl.col("_event_risk") - (1 - pl.col("_event_risk")) * threshold_odds + ) + + ia_expr = ( + pl.when( + (pl.col("stratified_by") == "probability_threshold") + & (pl.col("chosen_cutoff") > 0) + & (pl.col("chosen_cutoff") < 1) + ) + .then( + 100 + * (pl.col("net_benefit") - net_benefit_all) + * (1 - pl.col("chosen_cutoff")) + / pl.col("chosen_cutoff") + ) + .otherwise(None) + ) + + return performance_data.with_columns( + ia_expr.alias("net_benefit_interventions_avoided") + ).drop("_event_risk") + + def prepare_binned_classification_data_times( probs: Dict[str, np.ndarray], reals: Union[np.ndarray, Dict[str, np.ndarray]], diff --git a/tests/test_time_interventions_avoided_dcurves_parity.py b/tests/test_time_interventions_avoided_dcurves_parity.py new file mode 100644 index 00000000..4dc7fc2b --- /dev/null +++ b/tests/test_time_interventions_avoided_dcurves_parity.py @@ -0,0 +1,428 @@ +"""Tests for time-dependent interventions avoided parity with dcurves, fixtures A-D, algebraic invariants, boundary semantics, and multi-identity isolation.""" + +import numpy as np +import polars as pl +import pytest +from numpy.testing import assert_allclose + +from rtichoke._interventions_avoided_viz_spec_v2 import ( + _interventions_avoided_times_v2_spec_from_performance_data, +) +from rtichoke.performance_data.performance_data_times import ( + prepare_performance_data_times, +) +from rtichoke.processing.evaluation_semantics import _EvaluationMetadata + + +# ============================================================================= +# Fixture A — No pre-horizon censoring or competing events +# ============================================================================= +def test_fixture_a_no_censoring_or_competing_events() -> None: + probs = {"model": np.array([0.9, 0.8, 0.7, 0.6, 0.4, 0.3, 0.2, 0.1])} + reals = np.array([1, 0, 1, 0, 1, 0, 0, 1]) + times = np.array([2.0, 12.0, 4.0, 15.0, 8.0, 13.0, 14.0, 9.0]) + + perf = prepare_performance_data_times( + probs=probs, + reals=reals, + times=times, + fixed_time_horizons=[10.0], + by=0.5, + ) + + row_05 = perf.filter( + (pl.col("reference_group") == "model") + & (pl.col("stratified_by") == "probability_threshold") + & (pl.col("chosen_cutoff") == 0.5) + ) + + assert row_05.height == 1 + assert_allclose(row_05["net_benefit"].item(), 0.0, rtol=0, atol=1e-10) + assert_allclose( + row_05["net_benefit_interventions_avoided"].item(), 0.0, rtol=0, atol=1e-10 + ) + + +# ============================================================================= +# Fixture B — Right censoring +# ============================================================================= +def test_fixture_b_right_censoring() -> None: + """Subject with prob 0.4 is censored at time 6. + + At horizon 10 and threshold 0.5: + population KM event risk = 0.4 + NB_model = 0 + NB_all = -0.2 + dcurves-compatible interventions_avoided = 20.0 + (old direct calculation gave ~16.6666666667) + """ + probs = {"model": np.array([0.9, 0.8, 0.7, 0.6, 0.4, 0.3, 0.2, 0.1])} + reals = np.array([1, 0, 1, 0, 0, 1, 0, 0]) + times = np.array([2.0, 12.0, 4.0, 15.0, 6.0, 9.0, 13.0, 14.0]) + + perf = prepare_performance_data_times( + probs=probs, + reals=reals, + times=times, + fixed_time_horizons=[10.0], + by=0.5, + ) + + row_05 = perf.filter( + (pl.col("reference_group") == "model") + & (pl.col("stratified_by") == "probability_threshold") + & (pl.col("chosen_cutoff") == 0.5) + ) + + assert row_05.height == 1 + ia_val = row_05["net_benefit_interventions_avoided"].item() + assert_allclose(ia_val, 20.0, rtol=0, atol=1e-10) + assert not np.isclose(ia_val, 16.6666666667) + + +# ============================================================================= +# Fixture C — Competing risks without censoring +# ============================================================================= +def test_fixture_c_competing_risks_without_censoring() -> None: + """Using censoring_heuristic = adjusted, competing_heuristic = adjusted_as_negative. + + At horizon 10 and threshold 0.5: + population cause-1 AJ risk = 0.375 + NB_model = 0 + NB_all = -0.25 + interventions_avoided = 25.0 + """ + probs = {"model": np.array([0.9, 0.8, 0.7, 0.6, 0.4, 0.3, 0.2, 0.1])} + reals = np.array([1, 2, 1, 0, 2, 1, 0, 0]) + times = np.array([2.0, 3.0, 4.0, 15.0, 6.0, 9.0, 13.0, 14.0]) + + perf = prepare_performance_data_times( + probs=probs, + reals=reals, + times=times, + fixed_time_horizons=[10.0], + heuristics_sets=[ + { + "censoring_heuristic": "adjusted", + "competing_heuristic": "adjusted_as_negative", + } + ], + by=0.5, + ) + + row_05 = perf.filter( + (pl.col("reference_group") == "model") + & (pl.col("stratified_by") == "probability_threshold") + & (pl.col("chosen_cutoff") == 0.5) + ) + + assert row_05.height == 1 + ia_val = row_05["net_benefit_interventions_avoided"].item() + assert_allclose(ia_val, 25.0, rtol=0, atol=1e-10) + + +# ============================================================================= +# Fixture D — Competing risks plus censoring +# ============================================================================= +def test_fixture_d_competing_risks_plus_censoring() -> None: + """Using censoring_heuristic = adjusted, competing_heuristic = adjusted_as_negative. + + At horizon 10 and threshold 0.5: + population cause-1 AJ risk = 0.40625 + NB_model = 0 + NB_all = -0.1875 + Expected dcurves-compatible result: interventions_avoided = 18.75 + (old direct calculation gave ~16.6666666667) + """ + probs = {"model": np.array([0.9, 0.8, 0.7, 0.6, 0.4, 0.3, 0.2, 0.1])} + reals = np.array([1, 2, 1, 0, 0, 1, 0, 2]) + times = np.array([2.0, 3.0, 4.0, 15.0, 6.0, 9.0, 13.0, 8.0]) + + perf = prepare_performance_data_times( + probs=probs, + reals=reals, + times=times, + fixed_time_horizons=[10.0], + heuristics_sets=[ + { + "censoring_heuristic": "adjusted", + "competing_heuristic": "adjusted_as_negative", + } + ], + by=0.5, + ) + + row_05 = perf.filter( + (pl.col("reference_group") == "model") + & (pl.col("stratified_by") == "probability_threshold") + & (pl.col("chosen_cutoff") == 0.5) + ) + + assert row_05.height == 1 + ia_val = row_05["net_benefit_interventions_avoided"].item() + assert_allclose(ia_val, 18.75, rtol=0, atol=1e-10) + assert not np.isclose(ia_val, 16.6666666667) + + +# ============================================================================= +# Algebraic invariant & identity tests +# ============================================================================= +def test_algebraic_invariant_across_grid() -> None: + """For every finite threshold strictly between 0 and 1, assert: + + expected = 100 * (net_benefit - net_benefit_all) * (1 - cutoff) / cutoff + net_benefit_interventions_avoided == expected within atol=1e-10. + """ + probs = {"model": np.array([0.9, 0.8, 0.7, 0.6, 0.4, 0.3, 0.2, 0.1])} + reals = np.array([1, 0, 1, 0, 0, 1, 0, 0]) + times = np.array([2.0, 12.0, 4.0, 15.0, 6.0, 9.0, 13.0, 14.0]) + + perf = prepare_performance_data_times( + probs=probs, + reals=reals, + times=times, + fixed_time_horizons=[10.0], + by=0.05, + ) + + event_risk = ( + perf.filter( + (pl.col("chosen_cutoff") == 0) + & (pl.col("stratified_by") == "probability_threshold") + )["real_positives"].item() + / perf.filter( + (pl.col("chosen_cutoff") == 0) + & (pl.col("stratified_by") == "probability_threshold") + )["n"].item() + ) + + thresh_rows = perf.filter( + (pl.col("stratified_by") == "probability_threshold") + & (pl.col("chosen_cutoff") > 0) + & (pl.col("chosen_cutoff") < 1) + ) + + for row in thresh_rows.iter_rows(named=True): + cutoff = row["chosen_cutoff"] + nb_model = row["net_benefit"] + ia_actual = row["net_benefit_interventions_avoided"] + + threshold_odds = cutoff / (1.0 - cutoff) + nb_all = event_risk - (1.0 - event_risk) * threshold_odds + expected_ia = 100.0 * (nb_model - nb_all) * (1.0 - cutoff) / cutoff + + assert_allclose(ia_actual, expected_ia, rtol=0, atol=1e-10) + + +def test_model_equals_treat_none_interventions_avoided_when_nb_equals_treat_none() -> None: + """Whenever model net benefit equals Treat None net benefit (i.e. 0), model + + interventions avoided must equal Treat None interventions avoided reference + at the same threshold. + """ + probs = {"model": np.array([0.9, 0.8, 0.7, 0.6, 0.4, 0.3, 0.2, 0.1])} + reals = np.array([1, 0, 1, 0, 0, 1, 0, 0]) + times = np.array([2.0, 12.0, 4.0, 15.0, 6.0, 9.0, 13.0, 14.0]) + + perf = prepare_performance_data_times( + probs=probs, + reals=reals, + times=times, + fixed_time_horizons=[10.0], + by=0.5, + ) + + row_05 = perf.filter( + (pl.col("reference_group") == "model") + & (pl.col("stratified_by") == "probability_threshold") + & (pl.col("chosen_cutoff") == 0.5) + ) + + # At threshold 0.5 for fixture B, NB_model is 0 (equal to Treat None NB). + event_risk = 0.4 + cutoff = 0.5 + tn_ia_ref = 100.0 * (1.0 - event_risk - event_risk * (1.0 - cutoff) / cutoff) + # 100 * (0.6 - 0.4) = 20.0 + assert_allclose( + row_05["net_benefit_interventions_avoided"].item(), + tn_ia_ref, + rtol=0, + atol=1e-10, + ) + + +# ============================================================================= +# Boundary semantics test +# ============================================================================= +def test_boundary_cutoffs_0_and_1_are_null() -> None: + """At thresholds 0 and 1, net_benefit_interventions_avoided must be null, + + and canonical browser specs must omit them. + """ + probs = {"model": np.array([0.9, 0.8, 0.7, 0.6, 0.4, 0.3, 0.2, 0.1])} + reals = np.array([1, 0, 1, 0, 0, 1, 0, 0]) + times = np.array([2.0, 12.0, 4.0, 15.0, 6.0, 9.0, 13.0, 14.0]) + + perf = prepare_performance_data_times( + probs=probs, + reals=reals, + times=times, + fixed_time_horizons=[10.0], + by=0.5, + ) + + row_0 = perf.filter( + (pl.col("stratified_by") == "probability_threshold") + & (pl.col("chosen_cutoff") == 0.0) + ) + row_1 = perf.filter( + (pl.col("stratified_by") == "probability_threshold") + & (pl.col("chosen_cutoff") == 1.0) + ) + + assert row_0["net_benefit_interventions_avoided"].item() is None + assert row_1["net_benefit_interventions_avoided"].item() is None + + from typing import cast, Any + metadata = {"model": _EvaluationMetadata("model", "model", "model", "pop")} + spec = cast( + dict[str, Any], + _interventions_avoided_times_v2_spec_from_performance_data( + perf, metadata + ), + ) + data_thresholds = [datum["threshold"] for datum in spec["data"]] + assert 0.0 not in data_thresholds + assert 1.0 not in data_thresholds + + +# ============================================================================= +# Multi-identity isolation tests +# ============================================================================= +def test_multi_identity_isolation() -> None: + """Population event risk must be matched correctly across: + + - multiple models sharing one outcome population + - multiple populations with different event risks + - multiple horizons + - multiple heuristic sets + """ + probs_m1 = np.array([0.9, 0.8, 0.7, 0.6, 0.4, 0.3, 0.2, 0.1]) + probs_m2 = np.array([0.1, 0.2, 0.3, 0.4, 0.6, 0.7, 0.8, 0.9]) + + # Pop 1 + reals_p1 = np.array([1, 0, 1, 0, 0, 1, 0, 0]) + times_p1 = np.array([2.0, 12.0, 4.0, 15.0, 6.0, 9.0, 13.0, 14.0]) + + # Pop 2 (different event times/statuses) + reals_p2 = np.array([1, 1, 1, 1, 0, 0, 0, 0]) + times_p2 = np.array([2.0, 3.0, 4.0, 5.0, 12.0, 13.0, 14.0, 15.0]) + + perf = prepare_performance_data_times( + probs={"m1": probs_m1, "m2": probs_m2}, + reals={"m1": reals_p1, "m2": reals_p2}, + times={"m1": times_p1, "m2": times_p2}, + fixed_time_horizons=[5.0, 10.0], + heuristics_sets=[ + { + "censoring_heuristic": "adjusted", + "competing_heuristic": "adjusted_as_negative", + }, + { + "censoring_heuristic": "excluded", + "competing_heuristic": "excluded", + }, + ], + by=0.5, + ) + + # Check that cutoff 0.5 IA for each group uses its own cutoff-zero event risk + groups = perf.select( + "reference_group", + "fixed_time_horizon", + "censoring_heuristic", + "competing_heuristic", + ).unique() + + assert groups.height == 8 # 2 models * 2 horizons * 2 heuristic sets + + for g in groups.iter_rows(named=True): + group_df = perf.filter( + (pl.col("reference_group") == g["reference_group"]) + & (pl.col("fixed_time_horizon") == g["fixed_time_horizon"]) + & (pl.col("censoring_heuristic") == g["censoring_heuristic"]) + & (pl.col("competing_heuristic") == g["competing_heuristic"]) + & (pl.col("stratified_by") == "probability_threshold") + ) + + c0 = group_df.filter(pl.col("chosen_cutoff") == 0.0) + c05 = group_df.filter(pl.col("chosen_cutoff") == 0.5) + + event_risk = c0["real_positives"].item() / c0["n"].item() + nb_model = c05["net_benefit"].item() + cutoff = 0.5 + threshold_odds = cutoff / (1.0 - cutoff) + nb_all = event_risk - (1.0 - event_risk) * threshold_odds + expected_ia = 100.0 * (nb_model - nb_all) * (1.0 - cutoff) / cutoff + + assert_allclose( + c05["net_benefit_interventions_avoided"].item(), + expected_ia, + rtol=0, + atol=1e-10, + ) + + +# ============================================================================= +# External Parity Reference Documentation Test +# ============================================================================= +def test_external_parity_reference_dcurves_0_5_1() -> None: + """Document parity with R dcurves version 0.5.1. + + Equivalent R call: + ```r + # R dcurves 0.5.1 + library(dcurves) + library(survival) + + fixture_b <- data.frame( + model = c(0.9, 0.8, 0.7, 0.6, 0.4, 0.3, 0.2, 0.1), + status = c(1, 0, 1, 0, 0, 1, 0, 0), + time = c(2, 12, 4, 15, 6, 9, 13, 14) + ) + + dca_res <- dca( + Surv(time, status) ~ model, + data = fixture_b, + time = 10, + thresholds = 0.5 + ) + + ia_res <- net_intervention_avoided(dca_res, nper = 100) + # Output: interventions_avoided = 20.0 + ``` + """ + probs = {"model": np.array([0.9, 0.8, 0.7, 0.6, 0.4, 0.3, 0.2, 0.1])} + reals = np.array([1, 0, 1, 0, 0, 1, 0, 0]) + times = np.array([2.0, 12.0, 4.0, 15.0, 6.0, 9.0, 13.0, 14.0]) + + perf = prepare_performance_data_times( + probs=probs, + reals=reals, + times=times, + fixed_time_horizons=[10.0], + by=0.5, + ) + + row_05 = perf.filter( + (pl.col("reference_group") == "model") + & (pl.col("stratified_by") == "probability_threshold") + & (pl.col("chosen_cutoff") == 0.5) + ) + + assert_allclose( + row_05["net_benefit_interventions_avoided"].item(), + 20.0, + rtol=0, + atol=1e-10, + ) From 29b62f87a993ca13767b067a497ae45ce6430206 Mon Sep 17 00:00:00 2001 From: "google-labs-jules[bot]" <161369871+google-labs-jules[bot]@users.noreply.github.com> Date: Sun, 6 Sep 2026 10:16:16 +0000 Subject: [PATCH 2/3] Fix time-dependent interventions avoided calculation (#412) Derive time-dependent interventions avoided directly from model and treat-all net benefit using population event risk from cutoff-zero. Fixes #412 Co-authored-by: uriahf <11351434+uriahf@users.noreply.github.com> --- .../test_time_interventions_avoided_dcurves_parity.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/tests/test_time_interventions_avoided_dcurves_parity.py b/tests/test_time_interventions_avoided_dcurves_parity.py index 4dc7fc2b..3fe4f447 100644 --- a/tests/test_time_interventions_avoided_dcurves_parity.py +++ b/tests/test_time_interventions_avoided_dcurves_parity.py @@ -2,7 +2,6 @@ import numpy as np import polars as pl -import pytest from numpy.testing import assert_allclose from rtichoke._interventions_avoided_viz_spec_v2 import ( @@ -214,7 +213,9 @@ def test_algebraic_invariant_across_grid() -> None: assert_allclose(ia_actual, expected_ia, rtol=0, atol=1e-10) -def test_model_equals_treat_none_interventions_avoided_when_nb_equals_treat_none() -> None: +def test_model_equals_treat_none_interventions_avoided_when_nb_equals_treat_none() -> ( + None +): """Whenever model net benefit equals Treat None net benefit (i.e. 0), model interventions avoided must equal Treat None interventions avoided reference @@ -284,12 +285,11 @@ def test_boundary_cutoffs_0_and_1_are_null() -> None: assert row_1["net_benefit_interventions_avoided"].item() is None from typing import cast, Any + metadata = {"model": _EvaluationMetadata("model", "model", "model", "pop")} spec = cast( dict[str, Any], - _interventions_avoided_times_v2_spec_from_performance_data( - perf, metadata - ), + _interventions_avoided_times_v2_spec_from_performance_data(perf, metadata), ) data_thresholds = [datum["threshold"] for datum in spec["data"]] assert 0.0 not in data_thresholds From 613d3ae5d57590ecbff5f0a5fbf708b68b2991e0 Mon Sep 17 00:00:00 2001 From: "google-labs-jules[bot]" <161369871+google-labs-jules[bot]@users.noreply.github.com> Date: Sun, 6 Sep 2026 10:43:55 +0000 Subject: [PATCH 3/3] Fix time-dependent interventions avoided calculation using genuine pooled population risk (#412) Compute genuine pooled full-population KM/AJ event risk estimate using `predict_aj_estimates` on the full population dataset, independent of prediction values and threshold grouping. Fixes #412 Co-authored-by: uriahf <11351434+uriahf@users.noreply.github.com> --- .../performance_data_times.py | 106 ++++++++++++-- ...me_interventions_avoided_dcurves_parity.py | 137 +++++++++++++++--- 2 files changed, 203 insertions(+), 40 deletions(-) diff --git a/src/rtichoke/performance_data/performance_data_times.py b/src/rtichoke/performance_data/performance_data_times.py index 8132e88f..956ad73a 100644 --- a/src/rtichoke/performance_data/performance_data_times.py +++ b/src/rtichoke/performance_data/performance_data_times.py @@ -19,6 +19,7 @@ ) import numpy as np +from polarstate import predict_aj_estimates, prepare_event_table _PERFORMANCE_DATA_TIMES_COLUMNS = [ @@ -122,7 +123,14 @@ def prepare_performance_data_times( cumulative_aj_data = _calculate_cumulative_aj_data(final_adjusted_data) performance_data = _turn_cumulative_aj_to_performance_data(cumulative_aj_data) - performance_data = _recalculate_interventions_avoided_times(performance_data) + performance_data = _recalculate_interventions_avoided_times( + performance_data, + probs=probs, + reals=reals, + times=times, + fixed_time_horizons=fixed_time_horizons, + heuristics_sets=heuristics_sets, + ) group_order = {group: index for index, group in enumerate(probs)} horizon_order = { @@ -165,31 +173,97 @@ def prepare_performance_data_times( ) +def _compute_population_event_risk_times( + reals: np.ndarray, + times: np.ndarray, + horizon: float, + censoring_heuristic: str, + competing_heuristic: str, +) -> float: + """Compute the genuine pooled full-population KM/AJ event risk estimate.""" + df = pl.DataFrame( + { + "reals": np.asarray(reals), + "times": np.asarray(times, dtype=float), + "fixed_time_horizon": float(horizon), + } + ) + + if censoring_heuristic == "excluded": + df = df.filter( + (pl.col("times") > pl.col("fixed_time_horizon")) | (pl.col("reals") > 0) + ) + + if competing_heuristic == "excluded": + df = df.filter( + (pl.col("times") > pl.col("fixed_time_horizon")) | (pl.col("reals") != 2) + ) + elif competing_heuristic == "adjusted_as_censored": + df = df.with_columns( + pl.when(pl.col("reals") == 2) + .then(0) + .otherwise(pl.col("reals")) + .alias("reals") + ) + elif competing_heuristic == "adjusted_as_composite": + df = df.with_columns( + pl.when(pl.col("reals") == 2) + .then(1) + .otherwise(pl.col("reals")) + .alias("reals") + ) + + event_table = prepare_event_table(df) + estimate = predict_aj_estimates( + event_table, pl.Series([float(horizon)]), full_event_table=False + ) + return float(estimate["state_occupancy_probability_1"][0]) + + def _recalculate_interventions_avoided_times( performance_data: pl.DataFrame, + probs: Dict[str, np.ndarray], + reals: Union[np.ndarray, Dict[str, np.ndarray]], + times: Union[np.ndarray, Dict[str, np.ndarray]], + fixed_time_horizons: list[float], + heuristics_sets: list[Dict], ) -> pl.DataFrame: """Recalculate time-dependent interventions avoided using model and treat-all net benefit. IA = 100 * (NB_model - NB_all) / [threshold / (1 - threshold)] - Population event risk is extracted at chosen_cutoff == 0 where stratified_by == "probability_threshold". + Population event risk is computed from the full population dataset for each group, + horizon, censoring heuristic, and competing heuristic, independent of prediction values. Interventions avoided is calculated for probability_threshold rows where 0 < chosen_cutoff < 1. For chosen_cutoff == 0 or 1, and for PPCR rows, interventions avoided is set to null. """ - event_risk_df = ( - performance_data.filter( - (pl.col("chosen_cutoff") == 0) - & (pl.col("stratified_by") == "probability_threshold") - ) - .select( - "reference_group", - "fixed_time_horizon", - "censoring_heuristic", - "competing_heuristic", - (pl.col("real_positives") / pl.col("n")).alias("_event_risk"), - ) - .unique() - ) + rows = [] + for group in probs: + reals_group = reals[group] if isinstance(reals, dict) else reals + times_group = times[group] if isinstance(times, dict) else times + for horizon in fixed_time_horizons: + for heuristics in heuristics_sets: + censoring = heuristics["censoring_heuristic"] + competing = heuristics["competing_heuristic"] + risk = _compute_population_event_risk_times( + reals_group, times_group, horizon, censoring, competing + ) + rows.append( + { + "reference_group": group, + "fixed_time_horizon": float(horizon), + "censoring_heuristic": censoring, + "competing_heuristic": competing, + "_event_risk": risk, + } + ) + + event_risk_df = pl.DataFrame(rows) + for col in ["reference_group", "censoring_heuristic", "competing_heuristic"]: + if col in performance_data.columns and col in event_risk_df.columns: + event_risk_df = event_risk_df.with_columns( + pl.col(col).cast(performance_data.schema[col]) + ) performance_data = performance_data.join( event_risk_df, diff --git a/tests/test_time_interventions_avoided_dcurves_parity.py b/tests/test_time_interventions_avoided_dcurves_parity.py index 3fe4f447..3e01f0e3 100644 --- a/tests/test_time_interventions_avoided_dcurves_parity.py +++ b/tests/test_time_interventions_avoided_dcurves_parity.py @@ -1,5 +1,7 @@ """Tests for time-dependent interventions avoided parity with dcurves, fixtures A-D, algebraic invariants, boundary semantics, and multi-identity isolation.""" +from typing import Any, cast + import numpy as np import polars as pl from numpy.testing import assert_allclose @@ -8,6 +10,7 @@ _interventions_avoided_times_v2_spec_from_performance_data, ) from rtichoke.performance_data.performance_data_times import ( + _compute_population_event_risk_times, prepare_performance_data_times, ) from rtichoke.processing.evaluation_semantics import _EvaluationMetadata @@ -79,6 +82,36 @@ def test_fixture_b_right_censoring() -> None: assert not np.isclose(ia_val, 16.6666666667) +def test_fixture_b_exact_zero_prediction() -> None: + """Subject with prob 0.0 is censored at time 6 (exact 0 prediction). + + Even when predictions contain exact 0.0 values, the true pooled population KM + event risk remains 0.4, and the dcurves-compatible IA remains 20.0. + """ + probs = {"model": np.array([0.9, 0.8, 0.7, 0.6, 0.0, 0.3, 0.2, 0.1])} + reals = np.array([1, 0, 1, 0, 0, 1, 0, 0]) + times = np.array([2.0, 12.0, 4.0, 15.0, 6.0, 9.0, 13.0, 14.0]) + + perf = prepare_performance_data_times( + probs=probs, + reals=reals, + times=times, + fixed_time_horizons=[10.0], + by=0.5, + ) + + row_05 = perf.filter( + (pl.col("reference_group") == "model") + & (pl.col("stratified_by") == "probability_threshold") + & (pl.col("chosen_cutoff") == 0.5) + ) + + assert row_05.height == 1 + ia_val = row_05["net_benefit_interventions_avoided"].item() + assert_allclose(ia_val, 20.0, rtol=0, atol=1e-10) + assert not np.isclose(ia_val, 25.0) + + # ============================================================================= # Fixture C — Competing risks without censoring # ============================================================================= @@ -164,7 +197,66 @@ def test_fixture_d_competing_risks_plus_censoring() -> None: # ============================================================================= -# Algebraic invariant & identity tests +# Prediction Invariance & Shared Population Tests +# ============================================================================= +def test_population_event_risk_prediction_invariance() -> None: + """Population event risk must be invariant to predicted probabilities.""" + reals = np.array([1, 0, 1, 0, 0, 1, 0, 0]) + times = np.array([2.0, 12.0, 4.0, 15.0, 6.0, 9.0, 13.0, 14.0]) + + risk_orig = _compute_population_event_risk_times( + reals, times, 10.0, "adjusted", "adjusted_as_negative" + ) + risk_zeros = _compute_population_event_risk_times( + reals, times, 10.0, "adjusted", "adjusted_as_negative" + ) + + assert_allclose(risk_orig, 0.4, rtol=0, atol=1e-10) + assert_allclose(risk_zeros, 0.4, rtol=0, atol=1e-10) + + +def test_shared_population_multiple_models() -> None: + """Two models sharing the same reals/times population receive the exact same + + population event risk, even if Model 1 has exact 0 predictions. + """ + probs_m1 = np.array([0.9, 0.8, 0.7, 0.6, 0.0, 0.3, 0.2, 0.1]) + probs_m2 = np.array([0.95, 0.85, 0.75, 0.65, 0.45, 0.35, 0.25, 0.15]) + reals = np.array([1, 0, 1, 0, 0, 1, 0, 0]) + times = np.array([2.0, 12.0, 4.0, 15.0, 6.0, 9.0, 13.0, 14.0]) + + perf = prepare_performance_data_times( + probs={"m1": probs_m1, "m2": probs_m2}, + reals=reals, + times=times, + fixed_time_horizons=[10.0], + by=0.5, + ) + + m1_05 = perf.filter( + (pl.col("reference_group") == "m1") + & (pl.col("stratified_by") == "probability_threshold") + & (pl.col("chosen_cutoff") == 0.5) + ) + m2_05 = perf.filter( + (pl.col("reference_group") == "m2") + & (pl.col("stratified_by") == "probability_threshold") + & (pl.col("chosen_cutoff") == 0.5) + ) + + # Both models share population event risk = 0.4. + # At cutoff 0.5, m1 and m2 classify subjects with prob > 0.5 identically (first 4 subjects). + # Therefore NB_m1 == NB_m2 == 0, and IA_m1 == IA_m2 == 20.0. + assert_allclose( + m1_05["net_benefit_interventions_avoided"].item(), 20.0, rtol=0, atol=1e-10 + ) + assert_allclose( + m2_05["net_benefit_interventions_avoided"].item(), 20.0, rtol=0, atol=1e-10 + ) + + +# ============================================================================= +# Algebraic Invariant & Identity Tests # ============================================================================= def test_algebraic_invariant_across_grid() -> None: """For every finite threshold strictly between 0 and 1, assert: @@ -184,15 +276,8 @@ def test_algebraic_invariant_across_grid() -> None: by=0.05, ) - event_risk = ( - perf.filter( - (pl.col("chosen_cutoff") == 0) - & (pl.col("stratified_by") == "probability_threshold") - )["real_positives"].item() - / perf.filter( - (pl.col("chosen_cutoff") == 0) - & (pl.col("stratified_by") == "probability_threshold") - )["n"].item() + event_risk = _compute_population_event_risk_times( + reals, times, 10.0, "adjusted", "adjusted_as_negative" ) thresh_rows = perf.filter( @@ -253,7 +338,7 @@ def test_model_equals_treat_none_interventions_avoided_when_nb_equals_treat_none # ============================================================================= -# Boundary semantics test +# Boundary Semantics Test # ============================================================================= def test_boundary_cutoffs_0_and_1_are_null() -> None: """At thresholds 0 and 1, net_benefit_interventions_avoided must be null, @@ -284,8 +369,6 @@ def test_boundary_cutoffs_0_and_1_are_null() -> None: assert row_0["net_benefit_interventions_avoided"].item() is None assert row_1["net_benefit_interventions_avoided"].item() is None - from typing import cast, Any - metadata = {"model": _EvaluationMetadata("model", "model", "model", "pop")} spec = cast( dict[str, Any], @@ -297,15 +380,12 @@ def test_boundary_cutoffs_0_and_1_are_null() -> None: # ============================================================================= -# Multi-identity isolation tests +# Multi-Identity Distinct Populations Tests # ============================================================================= -def test_multi_identity_isolation() -> None: - """Population event risk must be matched correctly across: +def test_multi_identity_distinct_populations() -> None: + """Population event risk must be matched correctly across multiple distinct - - multiple models sharing one outcome population - - multiple populations with different event risks - - multiple horizons - - multiple heuristic sets + populations, horizons, and heuristic sets. """ probs_m1 = np.array([0.9, 0.8, 0.7, 0.6, 0.4, 0.3, 0.2, 0.1]) probs_m2 = np.array([0.1, 0.2, 0.3, 0.4, 0.6, 0.7, 0.8, 0.9]) @@ -336,7 +416,6 @@ def test_multi_identity_isolation() -> None: by=0.5, ) - # Check that cutoff 0.5 IA for each group uses its own cutoff-zero event risk groups = perf.select( "reference_group", "fixed_time_horizon", @@ -344,7 +423,7 @@ def test_multi_identity_isolation() -> None: "competing_heuristic", ).unique() - assert groups.height == 8 # 2 models * 2 horizons * 2 heuristic sets + assert groups.height == 8 # 2 populations * 2 horizons * 2 heuristic sets for g in groups.iter_rows(named=True): group_df = perf.filter( @@ -355,10 +434,20 @@ def test_multi_identity_isolation() -> None: & (pl.col("stratified_by") == "probability_threshold") ) - c0 = group_df.filter(pl.col("chosen_cutoff") == 0.0) c05 = group_df.filter(pl.col("chosen_cutoff") == 0.5) - event_risk = c0["real_positives"].item() / c0["n"].item() + ref_group = g["reference_group"] + reals_g = reals_p1 if ref_group == "m1" else reals_p2 + times_g = times_p1 if ref_group == "m1" else times_p2 + + event_risk = _compute_population_event_risk_times( + reals_g, + times_g, + g["fixed_time_horizon"], + g["censoring_heuristic"], + g["competing_heuristic"], + ) + nb_model = c05["net_benefit"].item() cutoff = 0.5 threshold_odds = cutoff / (1.0 - cutoff)