diff --git a/CHANGELOG.md b/CHANGELOG.md index 106adad2..b288a893 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -2,6 +2,7 @@ +- 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. diff --git a/src/rtichoke/calibration/calibration.py b/src/rtichoke/calibration/calibration.py index 0baba7da..581c237a 100644 --- a/src/rtichoke/calibration/calibration.py +++ b/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))) diff --git a/src/rtichoke/performance_data/performance_data_times.py b/src/rtichoke/performance_data/performance_data_times.py index 956ad73a..04ac5f2c 100644 --- a/src/rtichoke/performance_data/performance_data_times.py +++ b/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": diff --git a/src/rtichoke/processing/adjustments.py b/src/rtichoke/processing/adjustments.py index 519b7823..76459e17 100644 --- a/src/rtichoke/processing/adjustments.py +++ b/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 diff --git a/tests/test_exact_horizon_censoring.py b/tests/test_exact_horizon_censoring.py new file mode 100644 index 00000000..38b4b647 --- /dev/null +++ b/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 diff --git a/tests/test_heuristics.py b/tests/test_heuristics.py index 6730c1c9..f59c5523 100644 --- a/tests/test_heuristics.py +++ b/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"))