ml_meta_perf.validate

  1"""Validation protocols, baselines, and the term-count sweep.
  2
  3Four protocols are reported, and they answer different questions.
  4
  5* In-sample (IS) fits and scores the equation on all observed rows.
  6* Leave-one-dataset-out (LODO) holds out all rows of one dataset. It answers "what MCC will
  7  these models reach on a dataset nobody has run yet", which is the meta-learning use
  8  case, and it is the harder number because dataset features are constant within a
  9  fold and so the held-out dataset is genuinely unseen.
 10* Leave-one-model-out (LOMO) holds out all rows of one model. It answers "what will a new
 11  model reach on datasets we know".
 12* Doubly held out (DHO) removes every row sharing either the test dataset or the test
 13  model before fitting the coefficients for an observed pair.
 14
 15A random k-fold split is deliberately *not* offered as a headline protocol. With
 16dataset features constant across a dataset's rows, a random split puts the same dataset
 17on both sides of the fold and can therefore leak group identity. ``random_kfold_groups``
 18exists only to measure that diagnostic on the current corpus, not as a reported protocol.
 19
 20Study chapter: [5. Evaluation][study-chapter] -- the rationale, in
 21prose, with the figures.
 22
 23[study-chapter]: https://github.com/mariolpantunes/ml-meta-perf/blob/main/assets/docs/05-evaluation.md
 24"""
 25
 26from __future__ import annotations
 27
 28from collections.abc import Iterator
 29from dataclasses import dataclass, field
 30from math import comb
 31
 32import numpy as np
 33import polars as pl
 34
 35from ml_meta_perf.fit import Standardizer
 36from ml_meta_perf.model import MCC_LOWER, MCC_UPPER, Equation
 37from ml_meta_perf.search import Selector, guided_screen
 38from ml_meta_perf.stats import mae, pearson, r2_score, rmse, smape, spearman
 39from ml_meta_perf.terms import Library
 40
 41#: A fold scoring a negative R2 while its MAE stays under this is a near-constant target, not
 42#: a bad prediction. Set at a third of the corpus's own MCC standard deviation: below that,
 43#: an error is small in absolute terms and any negative R2 is the denominator, not the fit.
 44LOW_VARIANCE_MAE = 0.10
 45
 46
 47def leave_one_group_out(groups: np.ndarray) -> Iterator[tuple[str, np.ndarray, np.ndarray]]:
 48    """Yield ``(held-out label, train mask, test mask)`` for each distinct group."""
 49    for label in np.unique(groups):
 50        test = groups == label
 51        yield str(label), ~test, test
 52
 53
 54def random_kfold_groups(size: int, folds: int = 10, seed: int = 42) -> np.ndarray:
 55    """Fold labels for a random split -- used only to demonstrate the leakage it causes."""
 56    return np.random.default_rng(seed).integers(0, folds, size=size).astype(str)
 57
 58
 59@dataclass(frozen=True)
 60class Scores:
 61    """The metric set reported for every equation and every baseline."""
 62
 63    r2: float
 64    mae: float
 65    rmse: float
 66    smape: float
 67    spearman: float
 68    n: int
 69
 70    def as_dict(self) -> dict[str, float | int]:
 71        return {
 72            "r2": self.r2,
 73            "mae": self.mae,
 74            "rmse": self.rmse,
 75            "smape": self.smape,
 76            "spearman": self.spearman,
 77            "n": self.n,
 78        }
 79
 80
 81def score(truth: np.ndarray, prediction: np.ndarray) -> Scores:
 82    return Scores(
 83        r2=r2_score(truth, prediction),
 84        mae=mae(truth, prediction),
 85        rmse=rmse(truth, prediction),
 86        smape=smape(truth, prediction),
 87        spearman=spearman(truth, prediction),
 88        n=int(truth.shape[0]),
 89    )
 90
 91
 92@dataclass
 93class CrossValidation:
 94    """Out-of-fold predictions plus the terms each fold chose."""
 95
 96    predictions: np.ndarray
 97    per_fold: dict[str, Scores] = field(default_factory=dict)
 98    selected: list[list[str]] = field(default_factory=list)
 99    #: The equation each fold fitted, keyed by the label it held out. Kept so that a
100    #: correction fitted on a fold's training rows -- `ml_meta_perf.identity` fits one -- can
101    #: reuse the search this path already paid for instead of repeating it.
102    equations: dict[str, Equation] = field(default_factory=dict)
103    #: How many out-of-fold predictions `_clip_to_training` moved. Not a diagnostic detail:
104    #: the count is reported for each equation and protocol, since clipping can affect
105    #: the reported predictions. A study that clips has to say how often the clip fired.
106    clipped: int = 0
107
108    def scores(self, truth: np.ndarray) -> Scores:
109        return score(truth, self.predictions)
110
111    def dispersion(self) -> dict[str, float]:
112        """How the per-fold error is spread, next to the pooled number.
113
114        `scores` pools every out-of-fold prediction and scores it against the *global* mean.
115        On this corpus most of the variance is between datasets, which the dataset features
116        capture almost for free, so pooled R² flatters a LODO split: a set
117        reaching a pooled 0.53 was explaining 0.29 of the variance *within* the average
118        dataset. Pooling also hides that a single fold can dominate.
119
120        **Read the MAE columns, not ``worst_fold_r2``.** A per-fold R2 divides by the variance
121        of that fold's own target, and three datasets here have almost none: every one of
122        `5G_Slicing`'s 25 models scores about 0.986, a standard deviation of 0.052, and `NSR`
123        and `DeepSlice` are the same. Dividing an ordinary error by a variance that small
124        returns a large negative number for a nearly perfect prediction -- `5G_Slicing` reads
125        -10.6 while its predictions are all within 0.05 of the truth. Those three are the only
126        strongly negative folds, and the number is an artefact of the metric rather than a
127        failure of the equation.
128
129        MAE has no such denominator, so ``worst_fold_mae`` is the column that actually catches
130        a bad fold. ``worst_fold_r2`` is kept because it is what a reader expects to see, and
131        ``low_variance_folds`` counts the groups on which it cannot be trusted.
132        """
133        if not self.per_fold:
134            return {
135                "median_fold_r2": float("nan"),
136                "worst_fold_r2": float("nan"),
137                "median_fold_mae": float("nan"),
138                "worst_fold_mae": float("nan"),
139                "low_variance_folds": 0.0,
140                "clipped_predictions": float(self.clipped),
141                "folds": 0.0,
142            }
143        r2_values = [fold.r2 for fold in self.per_fold.values()]
144        mae_values = [fold.mae for fold in self.per_fold.values()]
145        return {
146            "median_fold_r2": float(np.median(r2_values)),
147            "worst_fold_r2": float(min(r2_values)),
148            "median_fold_mae": float(np.median(mae_values)),
149            "worst_fold_mae": float(max(mae_values)),
150            "low_variance_folds": float(self._low_variance_folds()),
151            "clipped_predictions": float(self.clipped),
152            "folds": float(len(r2_values)),
153        }
154
155    def _low_variance_folds(self) -> int:
156        """Folds whose per-fold R2 is not worth reading, inferred from the scores themselves.
157
158        A fold with a large negative R2 and a small MAE has a near-constant target, not a bad
159        prediction: R2 is dividing by nothing. The pair of conditions identifies that without
160        the class needing a copy of the target.
161        """
162        return sum(1 for fold in self.per_fold.values() if fold.r2 < 0.0 and fold.mae < LOW_VARIANCE_MAE)
163
164    def stability(self) -> pl.DataFrame:
165        """How often each term was selected across folds.
166
167        A term chosen in nearly every fold has stronger support than one chosen only a few
168        times. Reporting the final all-data equation without this column would present the
169        two identically.
170        """
171        counts: dict[str, int] = {}
172        for names in self.selected:
173            for name in names:
174                counts[name] = counts.get(name, 0) + 1
175        total = max(len(self.selected), 1)
176        return (
177            pl.DataFrame({"term": list(counts), "folds": list(counts.values())})
178            .with_columns((pl.col("folds") / total).alias("frequency"))
179            .sort("folds", descending=True)
180        )
181
182
183def fold_selections(
184    library: Library,
185    target: np.ndarray,
186    groups: np.ndarray,
187    *,
188    n_terms: int,
189    penalty: float,
190    pool_size: int = 250,
191    beam_width: int = 6,
192) -> list[list[str]]:
193    """Which terms selection picks when it is re-run inside each fold.
194
195    **This returns term names and nothing else, deliberately.** Re-selecting inside the folds
196    and scoring the result answers a question about the *discovery procedure*, not about the
197    equation, and the study does not report it; `cross_validate_fixed_form` produces every
198    cross-validated number. Returning no predictions means the re-selecting protocol cannot
199    produce a reported score even by accident.
200
201    What it is for is `CrossValidation.stability`. Fixing an equation's form is a claim that
202    the form describes the phenomenon rather than these 484 rows, and the way to check that is
203    to remove one dataset and see whether the same terms come back. A form whose terms
204    churn fold to fold has not earned the fixed-form protocol.
205    """
206    selections: list[list[str]] = []
207    for _, train, _ in leave_one_group_out(groups):
208        train_matrix = library.matrix[train]
209        standardizer = Standardizer.fit(train_matrix)
210        design = standardizer.apply(train_matrix)
211        pool = guided_screen(library_rows(library, train), target[train], keep=pool_size)
212        selector = Selector(design, target[train], penalty, library.feature_groups)
213        subsets = selector.search(pool, n_terms, beam_width=beam_width)
214        if n_terms in subsets:
215            selections.append([library.terms[index].name for index in subsets[n_terms].indices])
216    return selections
217
218
219def term_stability(selections: list[list[str]], required_terms: tuple[str, ...] = ()) -> pl.DataFrame:
220    """How often each term was selected across folds.
221
222    A term chosen in nearly every fold has stronger support than one chosen only a few times.
223    Reporting the final all-data equation without this column would present the two
224    identically.
225    """
226    counts: dict[str, int] = dict.fromkeys(required_terms, 0)
227    for names in selections:
228        for name in names:
229            counts[name] = counts.get(name, 0) + 1
230    total = max(len(selections), 1)
231    return (
232        pl.DataFrame({"term": list(counts), "folds": list(counts.values())})
233        .with_columns((pl.col("folds") / total).alias("frequency"))
234        .sort("folds", descending=True)
235    )
236
237
238def cross_validate_fixed_form(
239    library: Library,
240    columns: dict[str, np.ndarray],
241    target: np.ndarray,
242    groups: np.ndarray,
243    equations: dict[int, Equation],
244    *,
245    penalty: float,
246) -> dict[int, CrossValidation]:
247    """Cross-validate each equation with its **form fixed**: only the weights are refit.
248
249    **This is the study's cross-validation protocol.** Every reported leave-one-group-out
250    number comes from here. The choice is a position about what the artefact is rather than a
251    statistical convenience. The form of an equation is its
252    conceptual claim -- a statement about which quantities govern how well a learner does on a
253    dataset -- and for a simpler problem one would write that form down from domain expertise
254    and never search for it. What cross-validation then tests is whether the claim survives
255    data it has not seen, with its constants recalibrated.
256
257    `fold_selections` re-runs selection inside every fold and so answers a different
258    question -- whether the *discovery procedure* generalises. **It cannot be used for a reported score** -- it
259    returns no predictions. Twenty folds fit twenty different equations, so its pooled R2 is an
260    average over twenty models: it moved by up to 0.3 when the requested length changed by
261    two, while the fixed form varies by 0.014 over the same range. Its one remaining job is
262    `CrossValidation.stability`.
263
264    The obligation this creates is that the form must not be an artefact of the sample, and
265    `term_stability` over a `fold_selections` run is how that is checked --
266    separately, and not as part of the reported score.
267
268    Predictions are bounded by each training fold's own target range, as in
269    `fold_selections`; see `_clip_to_training`.
270    """
271    index = {name: position for position, name in enumerate(library.names)}
272    results: dict[int, CrossValidation] = {}
273    for size, equation in equations.items():
274        positions = [index[term.name] for term in equation.terms]
275        if not positions:
276            continue
277        matrix = library.matrix[:, positions]
278        outcome = CrossValidation(predictions=np.zeros_like(target))
279        for label, train, test in leave_one_group_out(groups):
280            standardizer = Standardizer.fit(matrix[train])
281            design = standardizer.apply(matrix[train])
282            gram = design.T @ design + penalty * np.eye(len(positions))
283            offset = float(target[train].mean())
284            weights = np.linalg.solve(gram, design.T @ (target[train] - offset))
285            held = standardizer.apply(matrix[test]) @ weights + offset
286            bounded = _clip_to_training(held, target[train])
287            outcome.clipped += int(np.count_nonzero(bounded != held))
288            outcome.predictions[test] = bounded
289            outcome.per_fold[label] = score(target[test], outcome.predictions[test])
290            outcome.selected.append([term.name for term in equation.terms])
291            # The same terms every fold, carrying that fold's weights, with the standardisation
292            # folded back in exactly as `fit.to_equation` does. Recorded so a downstream
293            # correction -- `identity.correct_out_of_fold` fits one -- can reuse the fold's own
294            # equation instead of the all-data one.
295            raw = weights / standardizer.scale
296            outcome.equations[label] = Equation(
297                intercept=offset - float(raw @ standardizer.mean),
298                terms=tuple(equation.terms),
299                weights=tuple(float(value) for value in raw),
300                standardized_weights=tuple(float(value) for value in weights),
301                name=f"fold_{label}",
302            )
303        results[size] = outcome
304    return results
305
306
307def cross_validate_doubly_held_out(
308    library: Library,
309    columns: dict[str, np.ndarray],
310    target: np.ndarray,
311    first: np.ndarray,
312    second: np.ndarray,
313    equations: dict[int, Equation],
314    *,
315    penalty: float,
316) -> dict[int, CrossValidation]:
317    """Refit each equation with **both** groups of a cell removed, and predict that cell.
318
319    This is the protocol for the question the study is actually for: *what MCC will this
320    model reach on this dataset*, when neither has been run. Neither single-group protocol
321    answers it. LODO holds the dataset out and keeps the model -- it has
322    seen that learner on nineteen other problems. LOMO does the reverse. For a
323    cell to be genuinely unseen, the whole dataset row-block *and* the whole model
324    column-block have to leave the training set, which is what this does.
325
326    **It is also the only protocol under which the comparison with the trivial baselines is
327    fair, and under it they do not exist.** A per-model mean or median answers "how well does
328    this model usually do", and a model held out of every fold has no rows to average. So a
329    baseline that beats the equation on ranking under LODO is a baseline
330    being handed the model identity the equation is denied; here it cannot be computed at
331    all, while the equation still predicts from features.
332
333    Costlier than the single-group protocols -- one solve per observed cell rather than per
334    group -- but each solve is a small ridge on a fixed design, so it is seconds, not minutes.
335    """
336    index = {name: position for position, name in enumerate(library.names)}
337    results: dict[int, CrossValidation] = {}
338    for size, equation in equations.items():
339        positions = [index[term.name] for term in equation.terms]
340        if not positions:
341            continue
342        matrix = library.matrix[:, positions]
343        outcome = CrossValidation(predictions=np.zeros_like(target))
344        for row_label in np.unique(first):
345            for column_label in np.unique(second):
346                test = (first == row_label) & (second == column_label)
347                if not test.any():
348                    continue
349                train = (first != row_label) & (second != column_label)
350                if not train.any():
351                    continue
352                standardizer = Standardizer.fit(matrix[train])
353                design = standardizer.apply(matrix[train])
354                gram = design.T @ design + penalty * np.eye(len(positions))
355                offset = float(target[train].mean())
356                weights = np.linalg.solve(gram, design.T @ (target[train] - offset))
357                held = standardizer.apply(matrix[test]) @ weights + offset
358                bounded = _clip_to_training(held, target[train])
359                outcome.clipped += int(np.count_nonzero(bounded != held))
360                outcome.predictions[test] = bounded
361        # Scored per held-out *dataset* rather than per cell: a cell is one to a few rows and
362        # a per-cell R2 has no variance to divide by. The dataset grouping is what the
363        # dispersion columns and the paired tests use everywhere else.
364        for label, _, test in leave_one_group_out(first):
365            outcome.per_fold[label] = score(target[test], outcome.predictions[test])
366            outcome.selected.append([term.name for term in equation.terms])
367        results[size] = outcome
368    return results
369
370
371def _clip_to_training(prediction: np.ndarray, training_target: np.ndarray) -> np.ndarray:
372    """Bound a fold's predictions by the target range that fold was trained on.
373
374    `model.predict` already clips to ``[MCC_LOWER, MCC_UPPER]``, the range MCC can take in
375    principle. That is the right bound for the equation as published and too loose here: this
376    corpus has exactly one negative row (-0.29), so a floor at -1.0 leaves a held-out fold free
377    to be predicted at a value nothing in training ever took.
378
379    The training fold's own observed range uses no held-out information -- it is exactly what
380    the fitted equation was shown. This bounds extrapolation during coefficient refitting;
381    it does not remove optimism from global form selection.
382    """
383    return np.clip(prediction, float(training_target.min()), float(training_target.max()))
384
385
386def library_rows(library: Library, mask: np.ndarray) -> Library:
387    """A library restricted to a row subset, without re-evaluating any term."""
388    clone = object.__new__(Library)
389    clone.terms = library.terms
390    clone.matrix = library.matrix[mask]
391    return clone
392
393
394def baseline_group_centre(
395    target: np.ndarray,
396    outer: np.ndarray,
397    inner: np.ndarray | None = None,
398    *,
399    centre: str = "mean",
400) -> np.ndarray:
401    """Predict a training-fold centre, optionally conditioned on a second grouping.
402
403    With ``inner`` set to the model column under a LODO split, this is
404    "what does this model usually score", the baseline any meta-model has to beat to be
405    worth writing down.
406
407    ``centre`` picks the summary. **Which one is the fair comparison depends on the metric
408    being reported, and using the mean for all of them understates the baseline.** The mean
409    minimises squared error, so it is the right opponent for R2 and RMSE. The *median*
410    minimises absolute error, so an MAE quoted against a mean baseline is quoted against a
411    baseline that is not even trying -- and SMAPE, being an absolute-error ratio, behaves the
412    same way. Both are reported so each metric is read against the baseline that is hardest
413    to beat on it.
414
415    The leave-one-group-out loop is not decoration: with 17-25 rows per group, including the
416    row being predicted inflates the per-model mean's R2 by 0.082. A naive group centre
417    computed over all rows is not this function.
418    """
419    summarise = np.median if centre == "median" else np.mean
420    predictions = np.zeros_like(target)
421    for _, train, test in leave_one_group_out(outer):
422        fallback = float(summarise(target[train]))
423        if inner is None:
424            predictions[test] = fallback
425            continue
426        for label in np.unique(inner[test]):
427            selected = test & (inner == label)
428            source = train & (inner == label)
429            predictions[selected] = float(summarise(target[source])) if source.any() else fallback
430    return np.clip(predictions, MCC_LOWER, MCC_UPPER)
431
432
433def baseline_group_mean(
434    target: np.ndarray,
435    outer: np.ndarray,
436    inner: np.ndarray | None = None,
437) -> np.ndarray:
438    """`baseline_group_centre` at the mean. Kept as the name the rest of the study uses."""
439    return baseline_group_centre(target, outer, inner, centre="mean")
440
441
442def additive_mean_reference(target: np.ndarray, first: np.ndarray, second: np.ndarray) -> np.ndarray:
443    """Dataset and model target means combined into a descriptive IS reference.
444
445    The calculation adds each observed dataset mean and model mean, subtracts the global
446    mean, and clips the result to the MCC range. It uses target identities and is therefore
447    not a deployable predictor. With an incomplete dataset-by-model grid it is also not the
448    jointly fitted least-squares optimum or a ceiling for additive equations.
449    """
450    grand = float(target.mean())
451    prediction = np.full_like(target, grand)
452    for group in (first, second):
453        for label in np.unique(group):
454            mask = group == label
455            prediction[mask] += float(target[mask].mean()) - grand
456    return np.clip(prediction, MCC_LOWER, MCC_UPPER)
457
458
459def interaction_oracle(
460    target: np.ndarray,
461    first: np.ndarray,
462    second: np.ndarray,
463    rank: int,
464) -> np.ndarray:
465    """The additive mean-based reference plus the best rank-``rank`` approximation of what it misses.
466
467    ``additive_mean_reference`` describes how much MCC follows dataset and model means.
468    The obvious next question is what the leftover looks like, and the leftover is not
469    noise: it is a (dataset x model) matrix of interactions, and interaction matrices are
470    usually dominated by a few components.
471
472    So the residual is decomposed by SVD and its leading components added back. This is
473    the AMMI model -- additive main effects, multiplicative interaction -- long used for
474    genotype-by-environment trials, which is structurally the same problem: a grid of
475    subjects crossed with conditions where some pairings suit each other.
476
477    The result is a *ladder* rather than a single ceiling. ``rank=0`` is the additive
478    mean-based reference; each further component is one more pattern of "this kind of model suits this
479    kind of dataset". At full rank it reproduces every observed cell and R2 is 1, which is
480    why the interesting question is how fast the ladder climbs, not where it ends.
481
482    Like ``additive_mean_reference`` this is fitted with the true values and predicts nothing.
483    Unobserved cells (24 of the 500 here) contribute zero residual, so they neither
484    distort the decomposition nor are counted in any score.
485    """
486    rows = np.unique(first)
487    columns = np.unique(second)
488    row_index = {label: position for position, label in enumerate(rows)}
489    column_index = {label: position for position, label in enumerate(columns)}
490
491    grid = np.full((rows.shape[0], columns.shape[0]), np.nan)
492    for value, row, column in zip(target, first, second, strict=True):
493        grid[row_index[row], column_index[column]] = value
494    observed = ~np.isnan(grid)
495
496    grand = float(np.nanmean(grid))
497    row_effect = np.nanmean(grid, axis=1) - grand
498    column_effect = np.nanmean(grid, axis=0) - grand
499    additive = grand + row_effect[:, None] + column_effect[None, :]
500
501    if rank > 0:
502        residual = np.where(observed, grid - additive, 0.0)
503        left, values, right = np.linalg.svd(residual, full_matrices=False)
504        additive = additive + (left[:, :rank] * values[:rank]) @ right[:rank]
505
506    prediction = np.array(
507        [additive[row_index[row], column_index[column]] for row, column in zip(first, second, strict=True)]
508    )
509    return np.clip(prediction, MCC_LOWER, MCC_UPPER)
510
511
512def oracle_ladder(
513    target: np.ndarray,
514    first: np.ndarray,
515    second: np.ndarray,
516    ranks: tuple[int, ...] = (0, 1, 2, 3, 4, 6, 8),
517) -> pl.DataFrame:
518    """How much each additional interaction component would be worth."""
519    rows: list[dict[str, object]] = []
520    previous: float | None = None
521    for rank in ranks:
522        value = r2_score(target, interaction_oracle(target, first, second, rank))
523        rows.append(
524            {
525                "interaction_rank": rank,
526                "r2": value,
527                "gain": value - previous if previous is not None else float("nan"),
528            }
529        )
530        previous = value
531    return pl.DataFrame(rows)
532
533
534#: A model is a right answer for the ranking if its MCC is within this of its dataset's best.
535#: See `ranking_report`: a fixed top-k relevance set misscores the datasets whose best models
536#: are genuinely tied.
537RELEVANCE_TOLERANCE = 0.01
538
539
540def average_precision(labels: np.ndarray, scores: np.ndarray) -> float:
541    """Area under the precision-recall curve, by the step-wise sum.
542
543    ``sum (R_n - R_{n-1}) * P_n`` over the ranking, the definition that does not interpolate
544    and so cannot flatter a short list. Returns NaN when no label is positive, which is a
545    real state here: at a threshold of 0.9 some datasets have no model above it.
546    """
547    if not labels.any():
548        return float("nan")
549    order = np.argsort(-scores, kind="stable")
550    hits = np.cumsum(labels[order])
551    precision = hits / np.arange(1, labels.size + 1)
552    return float((precision * labels[order]).sum() / labels.sum())
553
554
555def decision_report(
556    target: np.ndarray,
557    prediction: np.ndarray,
558    groups: np.ndarray | None = None,
559    thresholds: tuple[float, ...] = (0.5, 0.6, 0.7, 0.8, 0.9),
560) -> pl.DataFrame:
561    """Quality of the go/no-go decision the equation supports.
562
563    An R2 of 0.47 sounds too weak to act on, and read as "how precisely can I state the
564    MCC" it is. But the question a practitioner actually asks is coarser -- *will this
565    model work on this data* -- and a regression too imprecise for the first question can
566    be accurate on the second, because thresholding discards exactly the precision it
567    lacks.
568
569    Each row thresholds both the truth and the prediction at the same value and scores the
570    resulting binary decision. ``majority`` is the accuracy of always answering with the
571    larger class, which is the bar any such rule has to clear to be worth running.
572
573    ``f1`` is reported beside accuracy because the classes are far from balanced at the outer
574    thresholds -- at 0.9 a constant predictor reaches 0.49 accuracy and an F1 of exactly zero,
575    and only the pair distinguishes a rule that works from one that has guessed the majority.
576
577    ``map`` needs ``groups`` and is the mean over datasets of the average precision of the
578    ranking the prediction induces. It is the threshold-free companion: accuracy grades one
579    cut of an ordering, average precision grades the whole ordering. The mean is taken over
580    datasets rather than pooled over all rows because the question is always "which model for
581    *this* data" -- pooling would let a correct ordering between datasets hide a wrong one
582    within a dataset, and nobody ever chooses between datasets.
583    """
584    rows: list[dict[str, object]] = []
585    for threshold in thresholds:
586        actual = target >= threshold
587        predicted = prediction >= threshold
588        hits = int(np.sum(actual & predicted))
589        correct_rejections = int(np.sum(~actual & ~predicted))
590        false_alarms = int(np.sum(~actual & predicted))
591        misses = int(np.sum(actual & ~predicted))
592
593        positives = hits + false_alarms
594        actual_positives = hits + misses
595        denominator = float(
596            np.sqrt(
597                float(positives)
598                * float(actual_positives)
599                * float(correct_rejections + false_alarms)
600                * float(correct_rejections + misses)
601            )
602        )
603        share = float(np.mean(actual))
604        rows.append(
605            {
606                "threshold": threshold,
607                "accuracy": (hits + correct_rejections) / target.shape[0],
608                "majority": max(share, 1.0 - share),
609                "precision": hits / positives if positives else float("nan"),
610                "recall": hits / actual_positives if actual_positives else float("nan"),
611                "mcc": (
612                    (hits * correct_rejections - false_alarms * misses) / denominator
613                    if denominator > 0.0
614                    else float("nan")
615                ),
616                "f1": (
617                    2 * hits / (2 * hits + false_alarms + misses)
618                    if (2 * hits + false_alarms + misses)
619                    else float("nan")
620                ),
621                "map": _grouped_average_precision(target, prediction, groups, threshold),
622                "n_positive": actual_positives,
623            }
624        )
625    return pl.DataFrame(rows)
626
627
628def _grouped_average_precision(
629    target: np.ndarray, prediction: np.ndarray, groups: np.ndarray | None, threshold: float
630) -> float:
631    if groups is None:
632        return float("nan")
633    scores = [
634        average_precision(target[groups == label] >= threshold, prediction[groups == label])
635        for label in np.unique(groups)
636    ]
637    usable = [value for value in scores if not np.isnan(value)]
638    return float(np.mean(usable)) if usable else float("nan")
639
640
641def ranking_report(
642    target: np.ndarray,
643    prediction: np.ndarray,
644    groups: np.ndarray,
645) -> pl.DataFrame:
646    """Per-group ranking quality, scored the way a search engine is scored.
647
648    Only the head of the list is ever used: a practitioner tries the top few models and never
649    sees the tail, so a metric rewarding a correct rank twenty is measuring something nobody
650    reads. Hence ``mrr``, ``hit@1`` and ``regret`` -- all head-weighted -- with ``ap`` grading
651    the whole ordering and ``spearman`` kept for continuity.
652
653    **Read ``ap``, ``mrr``, ``hit@1`` and ``regret``; treat ``spearman`` as weak evidence.**
654    Measured on this corpus it sits between 0.63 and 0.73 for every predictor *and* every
655    baseline, including a constant, so it cannot separate the things this study compares.
656
657    A model counts as a right answer if its MCC is within ``RELEVANCE_TOLERANCE`` of the best
658    on its dataset, rather than by a fixed top-k cut. That is not a convenience: multiple observed
659    rows can share MCC 1.0, so datasets may have several tied best models, and
660    a top-3 rule would score a correct answer as a miss.
661
662    Regret is the practical question in the target's own units: if you pick the model this
663    equation ranks first, how much MCC do you give up against the best you could have picked?
664    """
665    rows: list[dict[str, object]] = []
666    for label in np.unique(groups):
667        mask = groups == label
668        truth, predicted = target[mask], prediction[mask]
669        relevant = truth >= truth.max() - RELEVANCE_TOLERANCE
670        order = np.argsort(-predicted, kind="stable")
671        ranks = [position for position, index in enumerate(order, 1) if relevant[index]]
672        rows.append(
673            {
674                "group": str(label),
675                "ap": average_precision(relevant, predicted),
676                "mrr": 1.0 / ranks[0] if ranks else 0.0,
677                "hit_at_1": float(relevant[order[0]]),
678                "regret": float(truth.max() - truth[int(order[0])]),
679                "spearman": spearman(predicted, truth),
680            }
681        )
682    return pl.DataFrame(rows).sort("group")
683
684
685#: Resamples for the paired bootstrap. Twenty groups is a small enough sample that the
686#: interval is the point of the exercise, and 20k resamples costs microseconds.
687BOOTSTRAP_ROUNDS = 20_000
688BOOTSTRAP_SEED = 42
689
690
691@dataclass(frozen=True)
692class PairedResult:
693    """Whether a per-group difference is real, or inside the spread between folds.
694
695    Every headline in this study is a mean over twenty held-out datasets, and a mean over
696    twenty groups is not a measurement until something says how much of it one group could
697    have produced. A higher mean can come from a few large wins while the other predictor
698    wins on more datasets. Hit@1 is coarser still: over twenty datasets it moves in steps of
699    0.05, so an apparently visible change can be only one or two datasets changing their top
700    pick.
701
702    So: a **sign test** on how many groups improved, which no single group can swing, and
703    a **paired bootstrap** interval on the mean, which shows what the mean is worth. The
704    two answer different questions and disagreeing is informative -- a significant sign
705    test with an interval spanning zero means a consistent but tiny effect.
706
707    ``mean`` and the interval are in the metric's own units, positive meaning the first
708    argument is better. Errors must therefore be negated by the caller; `paired_comparison`
709    does it for you when told the metric is an error.
710    """
711
712    mean: float
713    wins: int
714    losses: int
715    n: int
716    p_value: float
717    low: float
718    high: float
719
720    @property
721    def significant(self) -> bool:
722        """Whether the bootstrap interval excludes zero."""
723        return self.low > 0.0 or self.high < 0.0
724
725    def as_dict(self) -> dict[str, float | int]:
726        return {
727            "mean": self.mean,
728            "wins": self.wins,
729            "losses": self.losses,
730            "n_differing": self.n,
731            "p_value": self.p_value,
732            "ci_low": self.low,
733            "ci_high": self.high,
734        }
735
736
737def sign_test(differences: np.ndarray) -> tuple[int, int, float]:
738    """Wins, comparisons that differ, and the exact two-sided binomial p-value.
739
740    Exact rather than normal-approximate because the sample is twenty groups, where the
741    approximation is poor exactly when the answer matters. Summing binomial coefficients
742    directly is also why this needs no scipy: ``math.comb`` is in the standard library and
743    twenty choose ten is not a large number.
744
745    Ties are dropped rather than split, which is the conservative convention: a fold where
746    two predictors score identically is evidence for neither.
747    """
748    wins = int((differences > 0.0).sum())
749    n = int((differences != 0.0).sum())
750    if n == 0:
751        return 0, 0, 1.0
752    tail = sum(comb(n, k) for k in range(min(wins, n - wins) + 1))
753    return wins, n, min(1.0, 2.0 * tail / 2**n)
754
755
756def paired_comparison(
757    first: np.ndarray,
758    second: np.ndarray,
759    *,
760    lower_is_better: bool = False,
761    rounds: int = BOOTSTRAP_ROUNDS,
762    seed: int = BOOTSTRAP_SEED,
763) -> PairedResult:
764    """Compare two predictors group by group, on the same folds.
765
766    ``first`` and ``second`` are per-group scores in the same group order -- a column of
767    `ranking_report`, or per-fold errors. Pass ``lower_is_better`` for an error metric and
768    the sign is handled here, so a positive ``mean`` always means ``first`` won.
769
770    The bootstrap resamples *groups*, not rows. Rows within a dataset are not exchangeable
771    with rows in another one, and the question is whether the result survives a different
772    draw of datasets, which is the sampling this study's twenty-dataset corpus is a draw
773    from.
774    """
775    if first.shape != second.shape:
776        raise ValueError("paired comparison needs one score per group on both sides")
777    differences = (second - first) if lower_is_better else (first - second)
778    wins, n, p_value = sign_test(differences)
779    losses = int((differences < 0.0).sum())
780    if differences.size == 0:
781        return PairedResult(float("nan"), 0, 0, 0, 1.0, float("nan"), float("nan"))
782    generator = np.random.default_rng(seed)
783    draws = generator.integers(0, differences.size, size=(rounds, differences.size))
784    means = differences[draws].mean(axis=1)
785    return PairedResult(
786        mean=float(differences.mean()),
787        wins=wins,
788        losses=losses,
789        n=n,
790        p_value=p_value,
791        low=float(np.percentile(means, 2.5)),
792        high=float(np.percentile(means, 97.5)),
793    )
794
795
796def interaction_capture(
797    target: np.ndarray,
798    prediction: np.ndarray,
799    first: np.ndarray,
800    second: np.ndarray,
801    rank: int = 1,
802) -> dict[str, float]:
803    """How much of the leading interaction pattern the equation actually reaches.
804
805    `interaction_oracle` measures what a rank-``k`` interaction would be *worth* if
806    someone could predict it. It says nothing about whether the fitted equation gets any
807    of it, and chapter 5 originally answered that by comparing two R2 values -- which
808    cannot distinguish an equation that misses the pattern from one that finds it and is
809    inaccurate elsewhere.
810
811    This compares the two interaction *structures* directly. Both the truth and the
812    prediction are laid on the (dataset x model) grid and stripped of their own additive
813    part, leaving each side's interaction residual; the truth's is then reduced to its
814    leading ``rank`` components. ``alignment`` is the squared correlation between the two
815    over observed cells: 1.0 means the equation's interactions lie exactly along the
816    pattern the oracle found, 0.0 means they are unrelated to it.
817
818    Only observed cells count. The 24 absent ones contribute to neither side.
819    """
820    rows, columns = np.unique(first), np.unique(second)
821    row_index = {label: position for position, label in enumerate(rows)}
822    column_index = {label: position for position, label in enumerate(columns)}
823
824    def grid_of(values: np.ndarray) -> np.ndarray:
825        grid = np.full((rows.shape[0], columns.shape[0]), np.nan)
826        for value, row, column in zip(values, first, second, strict=True):
827            grid[row_index[row], column_index[column]] = value
828        return grid
829
830    def interaction(grid: np.ndarray) -> np.ndarray:
831        grand = float(np.nanmean(grid))
832        row_effect = np.nanmean(grid, axis=1) - grand
833        column_effect = np.nanmean(grid, axis=0) - grand
834        return grid - (grand + row_effect[:, None] + column_effect[None, :])
835
836    truth_grid = grid_of(target)
837    observed = ~np.isnan(truth_grid)
838    truth_interaction = np.where(observed, interaction(truth_grid), 0.0)
839    model_interaction = np.where(observed, interaction(grid_of(prediction)), 0.0)
840
841    left, values, right = np.linalg.svd(truth_interaction, full_matrices=False)
842    leading = (left[:, :rank] * values[:rank]) @ right[:rank]
843
844    wanted, reached = leading[observed], model_interaction[observed]
845    total = float((truth_interaction[observed] ** 2).sum())
846    return {
847        "alignment": float(pearson(wanted, reached) ** 2),
848        "leading_share": float((wanted**2).sum() / total) if total > 0.0 else float("nan"),
849        "interaction_share": total / float(((truth_grid[observed] - np.nanmean(truth_grid)) ** 2).sum())
850        if total > 0.0
851        else float("nan"),
852    }