diff --git a/src/rtichoke/performance_data/performance_data.py b/src/rtichoke/performance_data/performance_data.py index 1e37d053..d8902248 100644 --- a/src/rtichoke/performance_data/performance_data.py +++ b/src/rtichoke/performance_data/performance_data.py @@ -23,6 +23,7 @@ "reference_group", "stratified_by", "chosen_cutoff", + "probability_threshold", "true_positives", "true_negatives", "false_positives", @@ -244,9 +245,11 @@ def prepare_performance_data( Returns ------- pl.DataFrame - A Polars DataFrame where each row corresponds to a probability cutoff - for a given model/dataset. Columns include the cutoff value and a rich - set of performance metrics (e.g., `tpr`, `fpr`, `ppv`, `net_benefit`). + A Polars DataFrame where each row corresponds to an operating point for a + given model/dataset. Columns include `chosen_cutoff` (requested grid value), + `probability_threshold` (effective score boundary corresponding to the grid + point), `ppcr`, and a rich set of performance metrics (e.g., `sensitivity`, + `specificity`, `ppv`, `net_benefit`). Examples -------- @@ -274,6 +277,44 @@ def prepare_performance_data( performance_data = _turn_cumulative_aj_to_performance_data(cumulative_aj_data) + prob_threshold_rows = [] + for group_name, probs_array in probs.items(): + probs_g = np.asarray(probs_array, dtype=float) + unique_cutoffs = performance_data.filter( + pl.col("reference_group") == group_name + )["chosen_cutoff"].unique() + for c in unique_cutoffs: + eff_cutoff = float(np.quantile(probs_g, 1.0 - float(c), method="linear")) + prob_threshold_rows.append( + { + "reference_group": group_name, + "stratified_by": "ppcr", + "chosen_cutoff": float(c), + "probability_threshold": eff_cutoff, + } + ) + + ppcr_threshold_df = pl.DataFrame( + prob_threshold_rows, + schema={ + "reference_group": performance_data.schema["reference_group"], + "stratified_by": performance_data.schema["stratified_by"], + "chosen_cutoff": pl.Float64, + "probability_threshold": pl.Float64, + }, + ) + + performance_data = performance_data.join( + ppcr_threshold_df, + on=["reference_group", "stratified_by", "chosen_cutoff"], + how="left", + ).with_columns( + pl.when(pl.col("stratified_by") == "probability_threshold") + .then(pl.col("chosen_cutoff")) + .otherwise(pl.col("probability_threshold")) + .alias("probability_threshold") + ) + return performance_data.select(_PERFORMANCE_DATA_COLUMNS).sort( ["reference_group", "stratified_by", "chosen_cutoff"] ) diff --git a/src/rtichoke/performance_data/probs_distribution.py b/src/rtichoke/performance_data/probs_distribution.py index 46500b27..697497f9 100644 --- a/src/rtichoke/performance_data/probs_distribution.py +++ b/src/rtichoke/performance_data/probs_distribution.py @@ -202,7 +202,11 @@ def _prepare_probs_distribution_data( Returns ------- _PredictionDistributionData - TypedDict containing ``bins`` and ``operating_points`` Polars DataFrames. + TypedDict containing ``bins``, ``operating_points``, and ``rank_bins`` + Polars DataFrames. In ``operating_points``, ``value`` is the requested grid + point (requested PPCR or probability threshold), ``cutoff`` is the effective + predicted-probability boundary (from ``probability_threshold``), and + ``realized_ppcr`` is the empirical proportion classified positive. """ if not isinstance(stratified_by, (list, tuple)) or len(stratified_by) != 1: raise ValueError( @@ -270,7 +274,7 @@ def _prepare_probs_distribution_data( requested_value = float( row["ppcr"] if stratification_type == "ppcr" else row["chosen_cutoff"] ) - effective_cutoff = float(row["chosen_cutoff"]) + effective_cutoff = float(row["probability_threshold"]) n_observations = int(row["n"]) predicted_positives = int(row["predicted_positives"]) realized_ppcr = ( @@ -292,7 +296,11 @@ def _prepare_probs_distribution_data( ) # Build interval boundaries from effective cutoffs - cutoffs = evaluation_performance_data["chosen_cutoff"].to_numpy().astype(float) + cutoffs = ( + evaluation_performance_data["probability_threshold"] + .to_numpy() + .astype(float) + ) interval_boundaries = np.unique(np.concatenate(([0.0, 1.0], cutoffs))) interval_boundaries.sort() diff --git a/src/rtichoke/processing/transforms.py b/src/rtichoke/processing/transforms.py index 594ac996..0a0356a1 100644 --- a/src/rtichoke/processing/transforms.py +++ b/src/rtichoke/processing/transforms.py @@ -396,7 +396,10 @@ def _cast_and_join_adjusted_data_binary( & (pl.col("stratified_by") == "probability_threshold") ) | ( - ((1 - pl.col("chosen_cutoff")) >= pl.col("mid_point")) + ( + (1 - pl.col("chosen_cutoff")).round(10) + >= pl.col("mid_point").round(10) + ) & (pl.col("stratified_by") == "ppcr") ) ) diff --git a/tests/test_performance_schema_order.py b/tests/test_performance_schema_order.py index 5a21bc9d..de57cc5b 100644 --- a/tests/test_performance_schema_order.py +++ b/tests/test_performance_schema_order.py @@ -13,6 +13,7 @@ "reference_group", "stratified_by", "chosen_cutoff", + "probability_threshold", "true_positives", "true_negatives", "false_positives", diff --git a/tests/test_probs_distribution.py b/tests/test_probs_distribution.py index 4d139b79..8cb0de79 100644 --- a/tests/test_probs_distribution.py +++ b/tests/test_probs_distribution.py @@ -512,3 +512,239 @@ def test_primary_golden_fixture_order_invariance(): res_perm = _prepare_probs_distribution_data({"m1": probs_perm}, reals_perm, by=0.20) assert res_orig["rank_bins"].equals(res_perm["rank_bins"]) + + +def test_frozen_r_oracle_ppcr_tied_fixture(): + probs = {"m1": np.array([0.00, 0.15, 0.30, 0.50, 0.50, 0.50, 0.65, 0.80, 1.00])} + reals = np.array([0, 1, 0, 1, 0, 1, 1, 0, 1]) + + perf_df = prepare_performance_data(probs, reals, stratified_by=("ppcr",), by=0.20) + dist = _prepare_probs_distribution_data( + probs, reals, stratified_by=("ppcr",), by=0.20 + ) + ops = dist["operating_points"] + + # Independently hardcoded frozen R oracle expected values + expected = [ + # req_ppcr, eff_cutoff, realized_ppcr, tp, fp, tn, fn, sens, spec, ppv, npv, lift + (0.00, 1.00, 0.00 / 9.0, 0, 0, 4, 5, 0.0, 1.0, np.nan, 4 / 9, np.nan), + (0.20, 0.71, 2.00 / 9.0, 1, 1, 3, 4, 0.2, 0.75, 0.5, 3 / 7, 0.9), + (0.40, 0.50, 3.00 / 9.0, 2, 1, 3, 3, 0.4, 0.75, 2 / 3, 0.5, 1.2), + (0.60, 0.50, 3.00 / 9.0, 2, 1, 3, 3, 0.4, 0.75, 2 / 3, 0.5, 1.2), + (0.80, 0.24, 7.00 / 9.0, 4, 3, 1, 1, 0.8, 0.25, 4 / 7, 0.5, 36 / 35), + (1.00, 0.00, 9.00 / 9.0, 5, 4, 0, 0, 1.0, 0.0, 5 / 9, np.nan, 1.0), + ] + + for ( + req_ppcr, + eff_cutoff, + real_ppcr, + tp, + fp, + tn, + fn, + sens, + spec, + ppv, + npv, + lift, + ) in expected: + p_row = perf_df.filter(pl.col("chosen_cutoff") == req_ppcr).row(0, named=True) + op_row = ops.filter(pl.col("value") == req_ppcr).row(0, named=True) + + assert p_row["ppcr"] == pytest.approx(req_ppcr) + assert p_row["chosen_cutoff"] == pytest.approx(req_ppcr) + assert p_row["probability_threshold"] == pytest.approx(eff_cutoff, abs=1e-4) + + assert op_row["type"] == "ppcr" + assert op_row["value"] == pytest.approx(req_ppcr) + assert op_row["cutoff"] == pytest.approx(eff_cutoff, abs=1e-4) + assert op_row["realized_ppcr"] == pytest.approx(real_ppcr) + + assert p_row["true_positives"] == tp + assert p_row["false_positives"] == fp + assert p_row["true_negatives"] == tn + assert p_row["false_negatives"] == fn + + assert p_row["sensitivity"] == pytest.approx(sens) + assert p_row["specificity"] == pytest.approx(spec) + + if np.isnan(ppv): + assert np.isnan(p_row["ppv"]) + else: + assert p_row["ppv"] == pytest.approx(ppv) + + if np.isnan(npv): + assert np.isnan(p_row["npv"]) + else: + assert p_row["npv"] == pytest.approx(npv) + + if np.isnan(lift): + assert np.isnan(p_row["lift"]) + else: + assert p_row["lift"] == pytest.approx(lift) + + # Concrete assertion that operating_point.cutoff != operating_point.value for tied PPCR point + op_02 = ops.filter(pl.col("value") == 0.20).row(0, named=True) + assert op_02["cutoff"] != op_02["value"] + assert op_02["cutoff"] == pytest.approx(0.71) + assert op_02["value"] == pytest.approx(0.20) + + +def test_ppcr_repeated_effective_cutoffs_under_ties(): + probs = {"m1": np.array([0.00, 0.15, 0.30, 0.50, 0.50, 0.50, 0.65, 0.80, 1.00])} + reals = np.array([0, 1, 0, 1, 0, 1, 1, 0, 1]) + + perf_df = prepare_performance_data(probs, reals, stratified_by=("ppcr",), by=0.20) + dist = _prepare_probs_distribution_data( + probs, reals, stratified_by=("ppcr",), by=0.20 + ) + ops = dist["operating_points"] + + op_04 = ops.filter(pl.col("value") == 0.40).row(0, named=True) + op_06 = ops.filter(pl.col("value") == 0.60).row(0, named=True) + + perf_04 = perf_df.filter(pl.col("chosen_cutoff") == 0.40).row(0, named=True) + perf_06 = perf_df.filter(pl.col("chosen_cutoff") == 0.60).row(0, named=True) + + # Different requested value + assert op_04["value"] == pytest.approx(0.40) + assert op_06["value"] == pytest.approx(0.60) + assert op_04["value"] != op_06["value"] + + # Same effective cutoff, realized PPCR, and confusion matrix + assert op_04["cutoff"] == pytest.approx(op_06["cutoff"]) == pytest.approx(0.50) + assert ( + op_04["realized_ppcr"] + == pytest.approx(op_06["realized_ppcr"]) + == pytest.approx(3 / 9) + ) + + # Explicitly compare TP/FP/TN/FN + assert perf_04["true_positives"] == perf_06["true_positives"] == 2 + assert perf_04["false_positives"] == perf_06["false_positives"] == 1 + assert perf_04["true_negatives"] == perf_06["true_negatives"] == 3 + assert perf_04["false_negatives"] == perf_06["false_negatives"] == 3 + + +def test_ppcr_n_less_than_q(): + probs = {"m1": np.array([0.10, 0.50, 0.90])} + reals = np.array([0, 1, 1]) + + res = _prepare_probs_distribution_data( + probs, reals, stratified_by=("ppcr",), by=0.20 + ) + rb = res["rank_bins"] + ops = res["operating_points"] + + # Explicit empty rank strata preserved (q = 5) + assert len(rb) == 5 + assert rb["n_positive"].to_list() == [0, 0, 1, 0, 1] + assert rb["n_negative"].to_list() == [1, 0, 0, 0, 0] + + # Verify operating point cutoffs for N < q + cutoffs = ops["cutoff"].to_list() + assert len(cutoffs) == 6 + assert np.allclose(cutoffs, [0.9, 0.74, 0.58, 0.42, 0.26, 0.1]) + + +def test_ppcr_distinct_predictions_type7_boundaries(): + probs = {"m1": np.array([0.10, 0.30, 0.70, 0.90])} + reals = np.array([0, 1, 0, 1]) + + res = _prepare_probs_distribution_data( + probs, reals, stratified_by=("ppcr",), by=0.25 + ) + ops = res["operating_points"] + + # Expected Type-7 linear quantiles at 1.0, 0.75, 0.50, 0.25, 0.0 + expected_cutoffs = [0.90, 0.75, 0.50, 0.25, 0.10] + for row, exp_c in zip(ops.iter_rows(named=True), expected_cutoffs): + assert row["cutoff"] == pytest.approx(exp_c) + + +def test_ppcr_exact_scores_zero_and_one_endpoints(): + probs = {"m1": np.array([0.0, 0.0, 0.5, 1.0, 1.0])} + reals = np.array([0, 1, 0, 1, 0]) + + res = _prepare_probs_distribution_data( + probs, reals, stratified_by=("ppcr",), by=0.20 + ) + ops = res["operating_points"] + + op_0 = ops.filter(pl.col("value") == 0.0).row(0, named=True) + assert op_0["cutoff"] == pytest.approx(1.0) + assert op_0["realized_ppcr"] == pytest.approx(0.0) + + op_1 = ops.filter(pl.col("value") == 1.0).row(0, named=True) + assert op_1["cutoff"] == pytest.approx(0.0) + assert op_1["realized_ppcr"] == pytest.approx(1.0) + + +def test_ppcr_multiple_models_shared_outcomes(): + probs = { + "m1": np.array([0.1, 0.4, 0.6, 0.9]), + "m2": np.array([0.2, 0.5, 0.8, 1.0]), + } + reals = np.array([0, 1, 0, 1]) + + res = _prepare_probs_distribution_data( + probs, reals, stratified_by=("ppcr",), by=0.5 + ) + ops = res["operating_points"] + + ops_m1 = ops.filter(pl.col("model") == "m1") + ops_m2 = ops.filter(pl.col("model") == "m2") + + assert len(ops_m1) == 3 + assert len(ops_m2) == 3 + + # Per-model effective cutoffs computed independently + m1_c05 = ops_m1.filter(pl.col("value") == 0.5).row(0, named=True)["cutoff"] + m2_c05 = ops_m2.filter(pl.col("value") == 0.5).row(0, named=True)["cutoff"] + + assert m1_c05 == pytest.approx(0.5) + assert m2_c05 == pytest.approx(0.65) + + +def test_ppcr_multiple_populations_unequal_sample_sizes(): + probs = { + "pop1": np.array([0.1, 0.5, 0.9]), + "pop2": np.array([0.2, 0.4, 0.6, 0.8, 1.0]), + } + reals = { + "pop1": np.array([0, 1, 1]), + "pop2": np.array([0, 1, 0, 1, 0]), + } + + res = _prepare_probs_distribution_data( + probs, reals, stratified_by=("ppcr",), by=0.5 + ) + ops = res["operating_points"] + + p1_ops = ops.filter(pl.col("population") == "pop1") + p2_ops = ops.filter(pl.col("population") == "pop2") + + assert len(p1_ops) == 3 + assert len(p2_ops) == 3 + + +def test_ppcr_pre_post_non_regression(): + probs = {"m1": np.array([0.00, 0.15, 0.30, 0.50, 0.50, 0.50, 0.65, 0.80, 1.00])} + reals = np.array([0, 1, 0, 1, 0, 1, 1, 0, 1]) + + perf_df = prepare_performance_data(probs, reals, stratified_by=("ppcr",), by=0.20) + dist = _prepare_probs_distribution_data( + probs, reals, stratified_by=("ppcr",), by=0.20 + ) + + # Rank bins schema and masses frozen + rb = dist["rank_bins"] + assert rb["n_positive"].to_list() == [1, 2, 0, 1, 1] + assert rb["n_negative"].to_list() == [1, 2, 0, 0, 1] + + # Metrics frozen matching R oracle + assert perf_df["true_positives"].to_list() == [0, 1, 2, 2, 4, 5] + assert perf_df["false_positives"].to_list() == [0, 1, 1, 1, 3, 4] + assert perf_df["true_negatives"].to_list() == [4, 3, 3, 3, 1, 0] + assert perf_df["false_negatives"].to_list() == [5, 4, 3, 3, 1, 0]