From 01eb4a1db4fcc8c6a4a0b18cb18536dbbdf63d7b Mon Sep 17 00:00:00 2001
From: chross22 <52218551+chross22@users.noreply.github.com>
Date: Mon, 10 Aug 2026 15:34:20 -0400
Subject: [PATCH] Say whether the gap between two runs is real, and what it
 would take to see it

Two runs come back with two numbers - ROC AUC 0.854 against 0.892 - and
nothing in either says whether that is a difference between the models or a
difference between the stations the survey happened to visit. Two functions,
answering that and the question that should follow it.

compare_runs() takes fitted runs and refits nothing. Every run already stores
its held-out predictions with the station index and the fold, so the per-fold
metrics are already there to be paired; all this adds is the comparison. The
test is corrected_paired_test(), the same Nadeau-Bengio correction the
covariate jackknife uses and for the same reason - two cross-validation
training sets share most of their rows, so folds are not independent - but
two-sided, since asking which of two models is better has no direction the way
leaving a covariate out does.

The column worth having is `detectable`: the smallest true difference the
comparison would have found, at 80% power, given the fold-to-fold variance it
actually saw. "Not significant" on its own conflates "these models perform
alike" with "this survey could not tell them apart", and that column is what
separates them. It is deliberately not observed power, which is a
deterministic function of the p-value and tells a reader nothing new
(Hoenig and Heisey 2001).

Runs are paired on the station index rather than on row order, because tune
returns folds in its own order and two runs need not agree on it. Runs that
cover different stations - which is what different covariate sets produce, via
different missingness - are compared on what they share, and warn with the
count rather than intersecting silently.

power_curve() answers the other question by refitting at several subsample
sizes and tracing power against n. Everything is refitted at every size: a
model trained on half the stations is a different model, not the same model
evaluated on fewer, and a curve built the cheap way would measure the wrong
thing. Both runs see the same subsample and the same folds at every point,
which is what keeps the comparison paired the whole way down. It goes through
the jackknife's fold_scores(), so it pays for fitting and scoring and for
nothing else - no bootstrap intervals, no importance, no projection.

Both were checked against a run deliberately starved of its ocean covariates,
which compare_runs() puts 0.20 of AUC below the full model at p = 0.005, and
whose power curve climbs from 0.43 at a quarter of the stations to 1.00 at all
of them while the detectable difference falls from 0.43 to 0.13. Those are the
directions the functions exist to produce, and they are asserted rather than
eyeballed.

Two limits are documented rather than hidden. The curve's target difference
defaults to the one observed on the full data, so a gap that is itself mostly
noise gives a curve about a size that may not be real. And these are
differences in a bounded metric, so the normal-theory interval behind
`detectable` degrades as AUC approaches 1.

corrected_paired_test() gains an `alternative` argument, defaulting to the
one-sided behaviour the jackknife already relied on.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
---
 NAMESPACE                    |   2 +
 R/jackknife.R                |  34 ++-
 R/power.R                    | 487 +++++++++++++++++++++++++++++++++++
 README.md                    |  83 +++++-
 man/achieved_power.Rd        |  25 ++
 man/align_predictions.Rd     |  20 ++
 man/compare_one.Rd           |  38 +++
 man/compare_runs.Rd          | 109 ++++++++
 man/corrected_paired_test.Rd |  13 +-
 man/critical_t.Rd            |  20 ++
 man/minimum_detectable.Rd    |  33 +++
 man/power_curve.Rd           | 109 ++++++++
 man/power_point.Rd           |  41 +++
 man/run_predictions.Rd       |  18 ++
 man/yaml_bool_marker.Rd      |  18 +-
 tests/testthat/test-power.R  | 213 +++++++++++++++
 tools/citations.csv          |   1 +
 17 files changed, 1243 insertions(+), 21 deletions(-)
 create mode 100644 R/power.R
 create mode 100644 man/achieved_power.Rd
 create mode 100644 man/align_predictions.Rd
 create mode 100644 man/compare_one.Rd
 create mode 100644 man/compare_runs.Rd
 create mode 100644 man/critical_t.Rd
 create mode 100644 man/minimum_detectable.Rd
 create mode 100644 man/power_curve.Rd
 create mode 100644 man/power_point.Rd
 create mode 100644 man/run_predictions.Rd
 create mode 100644 tests/testthat/test-power.R

diff --git a/NAMESPACE b/NAMESPACE
index 1d1b454..c903863 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -10,6 +10,7 @@ export(available_stages)
 export(bathymetry_covariates)
 export(bootstrap_evaluation)
 export(climate_index_covariates)
+export(compare_runs)
 export(copernicus_client)
 export(copernicus_covariates)
 export(covariate_grid)
@@ -64,6 +65,7 @@ export(plot_roc_curve)
 export(plot_station_map)
 export(plot_station_series)
 export(plot_threshold_performance)
+export(power_curve)
 export(prejoin_steps)
 export(project_patch_model)
 export(projection_map)
diff --git a/R/jackknife.R b/R/jackknife.R
index ee9dc7e..69c66c4 100644
--- a/R/jackknife.R
+++ b/R/jackknife.R
@@ -469,11 +469,15 @@ subset_workflow <- function(train, vars, config, type) {
 #' test-set to training-set size in k-fold — which is the standard correction
 #' and costs roughly a factor of `sqrt(2)` off the statistic.
 #'
-#' One-sided, because the hypothesis is directional: the question is whether
-#' removing the covariate makes the model *worse*, and a covariate whose removal
-#' improves the model has failed the test rather than passed a different one.
+#' The default is one-sided, because the jackknife's hypothesis is directional:
+#' the question is whether removing the covariate makes the model *worse*, and a
+#' covariate whose removal improves the model has failed the test rather than
+#' passed a different one. Comparing two models is not directional in that way -
+#' either may be the better - so [compare_runs()] asks for `two.sided`.
 #'
 #' @param differences per-fold score of the full model minus the reduced one
+#' @param alternative `"greater"` for a directional hypothesis, `"two.sided"`
+#'   when either sign is a finding
 #' @return a list of `estimate`, `std_err`, `statistic`, `df`, `p_value`, `n`
 #' @references
 #' Nadeau C, Bengio Y (2003). Inference for the generalization error. *Machine
@@ -487,7 +491,9 @@ subset_workflow <- function(train, vars, config, type) {
 #' contains an underscore, and the citation checker's DOI pattern treats one as
 #' a terminator.
 #' @keywords internal
-corrected_paired_test <- function(differences) {
+corrected_paired_test <- function(differences,
+                                  alternative = c("greater", "two.sided")) {
+  alternative <- match.arg(alternative)
   usable <- differences[is.finite(differences)]
   k <- length(usable)
   none <- list(estimate = NA_real_, std_err = NA_real_, statistic = NA_real_,
@@ -500,18 +506,26 @@ corrected_paired_test <- function(differences) {
   std_err <- sqrt(stats::var(usable) * (1 / k + 1 / (k - 1)))
 
   if (!is.finite(std_err) || std_err == 0) {
-    # Identical on every fold. Either the covariate did exactly nothing, or it
-    # did the same thing everywhere - and only the second is evidence.
+    # Identical on every fold. Either the difference was exactly nothing, or it
+    # was the same everywhere - and only the second is evidence.
+    certain <- if (identical(alternative, "two.sided")) {
+      estimate != 0
+    } else {
+      estimate > 0
+    }
     return(list(estimate = estimate, std_err = 0,
                 statistic = if (estimate > 0) Inf else -Inf, df = k - 1,
-                p_value = if (estimate > 0) 0 else 1, n = k))
+                p_value = if (certain) 0 else 1, n = k))
   }
 
   statistic <- estimate / std_err
+  p_value <- if (identical(alternative, "two.sided")) {
+    2 * stats::pt(abs(statistic), df = k - 1, lower.tail = FALSE)
+  } else {
+    stats::pt(statistic, df = k - 1, lower.tail = FALSE)
+  }
   list(estimate = estimate, std_err = std_err, statistic = statistic,
-       df = k - 1,
-       p_value = stats::pt(statistic, df = k - 1, lower.tail = FALSE),
-       n = k)
+       df = k - 1, p_value = p_value, n = k)
 }
 
 #' The likelihood-based test, where the model type has one
diff --git a/R/power.R b/R/power.R
new file mode 100644
index 0000000..711a250
--- /dev/null
+++ b/R/power.R
@@ -0,0 +1,487 @@
+#' Is the gap between two model runs real?
+#'
+#' Two runs come back with two numbers — ROC AUC 0.857 against 0.871 — and
+#' nothing in either says whether the gap is a difference between the models or
+#' a difference between the stations the survey happened to visit. This answers
+#' that, and answers the question that should be asked next when the gap is not
+#' significant: **how large would a difference have had to be before this study
+#' could have seen it?**
+#'
+#' Those two are reported together deliberately. "Not significant" on its own is
+#' the least informative result in modelling — it conflates *these models
+#' perform alike* with *this survey could not have told them apart*, and
+#' `detectable` is what separates the two. A run that cannot detect anything
+#' smaller than 0.09 in AUC has not shown that a 0.014 gap is absent.
+#'
+#' @section How the comparison is paired:
+#' Every run cross-validates, so each carries a metric per fold rather than one
+#' number, and two runs on the same stations can be compared fold by fold. That
+#' pairing is most of the statistical power available: the folds vary a great
+#' deal between themselves and much less between two models scored on the *same*
+#' fold, and an unpaired comparison throws that away.
+#'
+#' Runs are matched on `.row`, the station index, not on position. Two runs with
+#' different covariates drop different stations to missingness, so the
+#' comparison is made on the stations both actually scored and the number
+#' dropped is reported. A run that drops many is telling you something, which is
+#' why this warns rather than silently intersecting.
+#'
+#' @section The test, and why it is not a plain t-test:
+#' The per-fold differences go through [corrected_paired_test()], the same
+#' Nadeau and Bengio (2003) correction the covariate jackknife uses, and for the
+#' same reason: any two cross-validation training sets share most of their rows,
+#' so folds are not independent and an uncorrected paired t-test finds
+#' significance that is not there.
+#'
+#' Two-sided here, unlike the jackknife. Leaving a covariate out has a direction
+#' worth testing against; asking which of two models is better does not.
+#'
+#' @section What "detectable" means:
+#' The smallest true difference this comparison would have found significant at
+#' `level`, with probability `power`, given the fold-to-fold variability it
+#' actually saw:
+#'
+#' \deqn{d_{min} = SE \times (t_{1-\alpha/2, df} + t_{power, df})}
+#'
+#' It is a property of **this** design — this many folds, these stations, this
+#' much variance between folds — not a general statement about the models. It
+#' says nothing about whether a smaller difference exists, only that this study
+#' would probably have missed it.
+#'
+#' @param runs a named list of two or more fitted runs, from
+#'   [fit_patch_model()] or [fit_patch_ensemble()]. The first is the reference
+#'   every other is compared against.
+#' @param metric `"roc_auc"` or `"pr_auc"`; both are threshold-free, which is
+#'   what lets them be compared fold by fold without a cutoff moving underneath
+#' @param level confidence level for the interval and the test
+#' @param power the power `detectable` is computed at
+#' @return a data frame with one row per comparison against the reference:
+#'   `reference`, `comparison`, `metric`, `reference_score`,
+#'   `comparison_score`, `difference` (comparison minus reference), `lower`,
+#'   `upper`, `statistic`, `df`, `p_value`, `detectable`, `n_folds`,
+#'   `n_stations`, and `n_dropped`
+#' @examples
+#' \dontrun{
+#' rf <- fit_patch_model(dat, within_config(config, type = "rf"))
+#' gam <- fit_patch_model(dat, within_config(config, type = "gam"))
+#'
+#' compare_runs(list(rf = rf, gam = gam))
+#' }
+#' @references
+#' Nadeau C, Bengio Y (2003). Inference for the generalization error. *Machine
+#' Learning* **52**(3), 239-281. \doi{10.1023/A:1024068626366} — the variance
+#' correction
+#'
+#' Dietterich TG (1998). Approximate statistical tests for comparing supervised
+#' classification learning algorithms. *Neural Computation* **10**(7),
+#' 1895-1923. \doi{10.1162/089976698300017197} — why comparing learning
+#' algorithms on shared folds needs one
+#'
+#' Hoenig JM, Heisey DM (2001). The abuse of power: the pervasive fallacy of
+#' power calculations for data analysis. *The American Statistician* **55**(1),
+#' 19-24. \doi{10.1198/000313001300339897} — why `detectable` is reported
+#' rather than the observed-power statistic it is often confused with
+#' @seealso [power_curve()] for how the answer changes with more stations
+#' @export
+compare_runs <- function(runs, metric = "roc_auc", level = 0.95, power = 0.8) {
+  if (!is.list(runs) || length(runs) < 2) {
+    stop("compare_runs() needs a list of at least 2 fitted runs; got ",
+         length(runs), ".", call. = FALSE)
+  }
+  if (is.null(names(runs)) || any(!nzchar(names(runs)))) {
+    stop("Name the runs, so the comparison table can say which is which: ",
+         "compare_runs(list(rf = rf_run, gam = gam_run)).", call. = FALSE)
+  }
+  if (!(metric %in% c("roc_auc", "pr_auc"))) {
+    stop("compare_runs() metric must be 'roc_auc' or 'pr_auc', got '", metric,
+         "'.\nBoth are threshold-free, which is what lets them be compared ",
+         "fold by fold.", call. = FALSE)
+  }
+
+  predictions <- lapply(runs, run_predictions)
+  reference <- names(runs)[1]
+
+  rows <- lapply(names(runs)[-1], function(name) {
+    compare_one(predictions[[reference]], predictions[[name]], reference, name,
+                metric = metric, level = level, power = power)
+  })
+
+  out <- do.call(rbind, rows)
+  rownames(out) <- NULL
+  out
+}
+
+#' The held-out predictions a run kept, checked
+#'
+#' @param run a [fit_patch_model()] or [fit_patch_ensemble()] result
+#' @return the predictions data frame
+#' @keywords internal
+run_predictions <- function(run) {
+  predictions <- run$predictions
+  needed <- c(".row", ".pred_patch", "patch", "id")
+  missing <- setdiff(needed, names(predictions))
+  if (is.null(predictions) || length(missing) > 0) {
+    stop("A run has no usable held-out predictions (missing: ",
+         paste(missing, collapse = ", "),
+         ").\ncompare_runs() reads the cross-validated predictions each run ",
+         "stores; a model fitted some other way cannot be compared this way.",
+         call. = FALSE)
+  }
+  predictions
+}
+
+#' Compare one run against the reference
+#'
+#' @param reference_predictions held-out predictions of the reference run
+#' @param other_predictions held-out predictions of the run being compared
+#' @param reference the reference run's name
+#' @param comparison the other run's name
+#' @param metric `"roc_auc"` or `"pr_auc"`
+#' @param level confidence level
+#' @param power the power `detectable` is computed at
+#' @return a one-row data frame
+#' @keywords internal
+compare_one <- function(reference_predictions, other_predictions, reference,
+                        comparison, metric = "roc_auc", level = 0.95,
+                        power = 0.8) {
+  shared <- intersect(reference_predictions$.row, other_predictions$.row)
+  dropped <- length(union(reference_predictions$.row,
+                          other_predictions$.row)) - length(shared)
+  if (length(shared) < 2) {
+    stop("'", reference, "' and '", comparison, "' share ", length(shared),
+         " stations, so there is nothing to compare.\nThey were probably ",
+         "fitted on different data rather than on different models.",
+         call. = FALSE)
+  }
+  if (dropped > 0) {
+    warning("'", reference, "' and '", comparison, "' do not cover the same ",
+            "stations: ", dropped, " of ", length(shared) + dropped,
+            " are in one run and not the other, and the comparison uses the ",
+            length(shared), " they share. Different covariates drop different ",
+            "stations to missingness.", call. = FALSE)
+  }
+
+  a <- align_predictions(reference_predictions, shared)
+  b <- align_predictions(other_predictions, shared)
+
+  # Fold membership is the reference run's. Two runs seeded alike on identical
+  # data fold identically, but a run that dropped rows folded differently, and
+  # then only one of the two labellings can define the pairing.
+  if (!identical(as.character(a$id), as.character(b$id))) {
+    warning("'", reference, "' and '", comparison, "' assigned these stations ",
+            "to different cross-validation folds, so the pairing uses '",
+            reference, "'s. The comparison stays valid - both models are ",
+            "scored on the same stations - but each fold is a held-out set ",
+            "for one model and not necessarily for the other.", call. = FALSE)
+  }
+
+  folds <- split(seq_along(shared), a$id)
+  score <- if (identical(metric, "pr_auc")) {
+    yardstick::pr_auc_vec
+  } else {
+    yardstick::roc_auc_vec
+  }
+  scored <- function(truth, probability) {
+    if (length(unique(truth)) < 2) return(NA_real_)
+    tryCatch(score(truth, probability), error = function(e) NA_real_,
+             warning = function(w) NA_real_)
+  }
+
+  per_fold <- vapply(folds, function(rows) {
+    c(reference = scored(a$patch[rows], a$.pred_patch[rows]),
+      comparison = scored(b$patch[rows], b$.pred_patch[rows]))
+  }, numeric(2))
+
+  differences <- per_fold["comparison", ] - per_fold["reference", ]
+  test <- corrected_paired_test(differences, alternative = "two.sided")
+
+  data.frame(
+    reference = reference,
+    comparison = comparison,
+    metric = metric,
+    reference_score = mean(per_fold["reference", ], na.rm = TRUE),
+    comparison_score = mean(per_fold["comparison", ], na.rm = TRUE),
+    difference = test$estimate,
+    lower = test$estimate - critical_t(level, test$df) * test$std_err,
+    upper = test$estimate + critical_t(level, test$df) * test$std_err,
+    statistic = test$statistic,
+    df = test$df,
+    p_value = test$p_value,
+    detectable = minimum_detectable(test$std_err, test$df, level, power),
+    n_folds = test$n,
+    n_stations = length(shared),
+    n_dropped = dropped,
+    stringsAsFactors = FALSE
+  )
+}
+
+#' One run's predictions, restricted to shared stations and put in their order
+#'
+#' @param predictions a run's held-out predictions
+#' @param shared the station indices to keep
+#' @return the matching rows, in `shared` order
+#' @keywords internal
+align_predictions <- function(predictions, shared) {
+  predictions[match(shared, predictions$.row), , drop = FALSE]
+}
+
+#' The two-sided critical value
+#'
+#' @param level confidence level
+#' @param df degrees of freedom
+#' @return the critical `t`, or `NA_real_`
+#' @keywords internal
+critical_t <- function(level, df) {
+  if (is.na(df) || df < 1) return(NA_real_)
+  stats::qt(1 - (1 - level) / 2, df = df)
+}
+
+#' The smallest difference a design could have detected
+#'
+#' Reported instead of "observed power", which is the statistic this is usually
+#' confused with and which carries no information a p-value does not — it is a
+#' deterministic function of it (Hoenig and Heisey 2001). The minimum detectable
+#' difference is about the *design* rather than about the result, which is what
+#' makes it worth reading beside a null finding.
+#'
+#' @param std_err the corrected standard error of the difference
+#' @param df degrees of freedom
+#' @param level confidence level
+#' @param power the power to solve at
+#' @return the smallest detectable difference, or `NA_real_`
+#' @references
+#' Hoenig JM, Heisey DM (2001). The abuse of power: the pervasive fallacy of
+#' power calculations for data analysis. *The American Statistician* **55**(1),
+#' 19-24. \doi{10.1198/000313001300339897}
+#' @keywords internal
+minimum_detectable <- function(std_err, df, level = 0.95, power = 0.8) {
+  if (is.na(std_err) || is.na(df) || df < 1 || !is.finite(std_err)) {
+    return(NA_real_)
+  }
+  std_err * (critical_t(level, df) + stats::qt(power, df = df))
+}
+
+#' The power a design has against a stated difference
+#'
+#' The complement of [minimum_detectable()]: given how variable the folds were,
+#' how often would a true difference of `difference` be called significant?
+#'
+#' @param difference the true difference to detect
+#' @param std_err the corrected standard error of the difference
+#' @param df degrees of freedom
+#' @param level confidence level
+#' @return a probability, or `NA_real_`
+#' @keywords internal
+achieved_power <- function(difference, std_err, df, level = 0.95) {
+  if (is.na(std_err) || is.na(df) || df < 1 || std_err <= 0) return(NA_real_)
+  critical <- critical_t(level, df)
+  ncp <- abs(difference) / std_err
+  # Both tails, so a large difference in either direction counts as detected.
+  stats::pt(-critical, df = df, ncp = ncp) +
+    stats::pt(critical, df = df, ncp = ncp, lower.tail = FALSE)
+}
+
+#' How the comparison would improve with more stations
+#'
+#' [compare_runs()] answers what this survey could see. This answers what a
+#' larger one would: it refits both runs on subsamples of the stations, at
+#' several sizes, and traces how the power to detect a difference grows with
+#' `n`.
+#'
+#' The curve is the useful artefact rather than any single number on it. Power
+#' against sample size is steeply non-linear, and where a study sits on that
+#' curve decides what the next survey is worth: a comparison at 0.35 power is
+#' one more season away from being decisive, and one at 0.9 will not be improved
+#' by more stations because it is already there.
+#'
+#' @section What it costs, and what it therefore skips:
+#' `fractions × replicates × runs` cross-validations. Everything is refitted at
+#' every size — the point is precisely that a model trained on half the stations
+#' is a different model, not the same model evaluated on fewer — so this is the
+#' expensive function in the package and parallelises over the whole grid.
+#'
+#' It goes through the same fold-scoring path the covariate jackknife uses,
+#' which fits and scores and stops there. No bootstrap intervals, no variable
+#' importance, no projection: none of it enters the curve, and all of it would
+#' be paid for at every point.
+#'
+#' @section Subsampling stations, not folds:
+#' Rows are drawn without replacement, and the folds are then built inside each
+#' subsample. Reusing the full run's folds and thinning them would shrink the
+#' held-out sets while leaving the training sets nearly whole, which measures
+#' something else entirely — the curve has to come from models that were
+#' actually trained on less.
+#'
+#' Both runs see the **same** subsample and the **same** folds at every point,
+#' which is what keeps the comparison paired all the way down the curve.
+#'
+#' @section Reading it honestly:
+#' The target difference defaults to the one observed on the full data, and that
+#' is an estimate, not a truth. If the observed gap is itself mostly noise, the
+#' curve answers "how many stations to reliably detect a difference this size"
+#' for a size that may not be real. It is a projection under an assumption, and
+#' the assumption is the observed effect.
+#'
+#' @param dat labeled modeling data from `label_patch()` with covariates attached
+#' @param configs a named list of two or more configs to compare. The first is
+#'   the reference; each is fitted exactly as its own run would be.
+#' @param fractions the shares of the stations to fit at
+#' @param replicates how many subsamples per fraction; the spread across them is
+#'   what stops one unlucky draw from setting a point on the curve
+#' @param difference the true difference to compute power against; `NULL` uses
+#'   the one observed at the largest fraction
+#' @param metric `"roc_auc"` or `"pr_auc"`
+#' @param level confidence level the test would use
+#' @param workers how many workers; see [resolve_workers()]
+#' @param seed a seed, so a curve is reproducible
+#' @return a data frame with one row per fraction: `fraction`, `n_stations`,
+#'   `replicates`, `difference` (mean observed), `std_err`, `df`, `power`, and
+#'   `detectable`. The target difference is on it as a `difference` attribute.
+#' @examples
+#' \dontrun{
+#' rf <- config; rf$model$type <- "rf"
+#' gam <- config; gam$model$type <- "gam"
+#'
+#' curve <- power_curve(dat, list(rf = rf, gam = gam))
+#' curve[c("n_stations", "power", "detectable")]
+#' }
+#' @references
+#' Nadeau C, Bengio Y (2003). Inference for the generalization error. *Machine
+#' Learning* **52**(3), 239-281. \doi{10.1023/A:1024068626366}
+#' @seealso [compare_runs()], which answers the same question for the data you
+#'   already have
+#' @export
+power_curve <- function(dat, configs, fractions = c(0.25, 0.5, 0.75, 1),
+                        replicates = 5, difference = NULL,
+                        metric = "roc_auc", level = 0.95, workers = NULL,
+                        seed = 42) {
+  if (!is.list(configs) || length(configs) < 2) {
+    stop("power_curve() needs at least 2 configs to compare; got ",
+         length(configs), ".", call. = FALSE)
+  }
+  if (is.null(names(configs)) || any(!nzchar(names(configs)))) {
+    stop("Name the configs, so the curve can say what it compared: ",
+         "power_curve(dat, list(rf = rf_config, gam = gam_config)).",
+         call. = FALSE)
+  }
+  fractions <- sort(unique(fractions[fractions > 0 & fractions <= 1]))
+  if (length(fractions) == 0) {
+    stop("fractions must be shares of the stations, in (0, 1].", call. = FALSE)
+  }
+
+  # Every run is scored on the rows every run can use, so a config whose
+  # covariates are missing somewhere does not quietly get a different sample.
+  predictors <- lapply(configs, function(config) predictor_names(dat, config))
+  shared <- Reduce(union, predictors)
+  model_data <- as.data.frame(dat[c(shared, "patch")])
+  model_data <- model_data[stats::complete.cases(model_data), , drop = FALSE]
+  if (nrow(model_data) < 10) {
+    stop("Only ", nrow(model_data), " stations are complete across every ",
+         "config's covariates, which is too few to build a power curve from.",
+         call. = FALSE)
+  }
+
+  types <- vapply(configs, resolve_model_type, character(1))
+  for (type in unique(types)) check_model_packages(type)
+
+  tasks <- expand.grid(fraction = fractions, replicate = seq_len(replicates),
+                       KEEP.OUT.ATTRS = FALSE)
+  workers <- resolve_workers(workers, nrow(tasks))
+  message("  power curve: ", nrow(tasks), " subsamples x ", length(configs),
+          " runs across ", workers, if (workers == 1) " worker" else " workers")
+
+  drawn <- taupatch_lapply(seq_len(nrow(tasks)), function(i) {
+    power_point(model_data, configs, predictors, types,
+                fraction = tasks$fraction[i], replicate = tasks$replicate[i],
+                metric = metric, seed = seed)
+  }, workers = workers, seed = seed)
+
+  drawn <- do.call(rbind, Filter(Negate(is.null), drawn))
+  if (is.null(drawn) || nrow(drawn) == 0) {
+    stop("No subsample produced a usable comparison. The runs may be failing ",
+         "to fit on a fraction of these stations.", call. = FALSE)
+  }
+
+  # The difference to have power against. Taken from the largest fraction,
+  # which is the best estimate of it available.
+  target <- difference %||% mean(
+    drawn$difference[drawn$fraction == max(drawn$fraction)], na.rm = TRUE
+  )
+
+  out <- do.call(rbind, lapply(fractions, function(fraction) {
+    at <- drawn[drawn$fraction == fraction, , drop = FALSE]
+    usable <- at[is.finite(at$std_err) & is.finite(at$df), , drop = FALSE]
+    if (nrow(usable) == 0) return(NULL)
+
+    # Averaged across replicates rather than pooled: each replicate is its own
+    # complete comparison, and averaging their standard errors is what keeps one
+    # unlucky draw from setting the point.
+    std_err <- mean(usable$std_err)
+    df <- mean(usable$df)
+
+    data.frame(
+      fraction = fraction,
+      n_stations = round(mean(usable$n_stations)),
+      replicates = nrow(usable),
+      difference = mean(usable$difference, na.rm = TRUE),
+      std_err = std_err,
+      df = df,
+      power = achieved_power(target, std_err, df, level),
+      detectable = minimum_detectable(std_err, df, level, power = 0.8),
+      stringsAsFactors = FALSE
+    )
+  }))
+
+  rownames(out) <- NULL
+  attr(out, "difference") <- target
+  out
+}
+
+#' One point on the power curve
+#'
+#' Draws a subsample, folds it, and scores every config on those same folds.
+#'
+#' @param model_data the complete-case modeling data
+#' @param configs the configs being compared
+#' @param predictors each config's predictors
+#' @param types each config's model type
+#' @param fraction the share of stations to draw
+#' @param replicate which draw this is
+#' @param metric `"roc_auc"` or `"pr_auc"`
+#' @param seed the run's seed
+#' @return a one-row data frame, or `NULL` when the draw could not be scored
+#' @keywords internal
+power_point <- function(model_data, configs, predictors, types, fraction,
+                        replicate, metric = "roc_auc", seed = 42) {
+  # Reproducible per point rather than per call, so a curve is the same however
+  # its tasks were spread across workers.
+  set.seed(seed + replicate * 1000L + round(fraction * 100))
+
+  n <- max(10L, round(nrow(model_data) * fraction))
+  if (n > nrow(model_data)) n <- nrow(model_data)
+  drawn <- model_data[sample.int(nrow(model_data), n), , drop = FALSE]
+  if (length(unique(drawn$patch)) < 2) return(NULL)
+
+  folds <- tryCatch(
+    rsample::vfold_cv(drawn, v = configs[[1]]$model$cv_folds %||% 10,
+                      strata = "patch"),
+    error = function(e) NULL
+  )
+  if (is.null(folds)) return(NULL)
+
+  scores <- lapply(names(configs), function(name) {
+    fold_scores(intersect(predictors[[name]], names(drawn)), folds,
+                configs[[name]], types[[name]], metric)
+  })
+  names(scores) <- names(configs)
+
+  reference <- scores[[1]]
+  differences <- scores[[2]] - reference
+  test <- corrected_paired_test(differences, alternative = "two.sided")
+
+  data.frame(
+    fraction = fraction, replicate = replicate, n_stations = n,
+    difference = test$estimate, std_err = test$std_err, df = test$df,
+    stringsAsFactors = FALSE
+  )
+}
diff --git a/README.md b/README.md
index b6868e0..faa0357 100644
--- a/README.md
+++ b/README.md
@@ -24,6 +24,7 @@
   - [Testing which covariates earn their place](#testing-which-covariates-earn-their-place)
   - [Model type](#model-type)
   - [Fitting all of them at once](#fitting-all-of-them-at-once)
+  - [Comparing two runs, and knowing when you cannot](#comparing-two-runs-and-knowing-when-you-cannot)
   - [Training and projection windows](#training-and-projection-windows)
   - [How far to trust a map](#how-far-to-trust-a-map)
   - [Reading the evaluation](#reading-the-evaluation)
@@ -1031,6 +1032,79 @@ those replicates are pooled in proportion to member weight, and the algorithm
 disagreement is reported on top in its own column — so a projection carries both
 without either standing in for the other.
 
+### Comparing two runs, and knowing when you cannot
+
+Two runs come back with two numbers — ROC AUC 0.854 against 0.892 — and nothing
+in either says whether that gap is a difference between the models or a
+difference between the stations the survey happened to visit. `compare_runs()`
+answers that, from the held-out predictions each run already stores, so it
+refits nothing:
+
+```r
+rf   <- fit_patch_model(dat, rf_config)
+gam  <- fit_patch_model(dat, gam_config)
+weak <- fit_patch_model(dat, no_ocean_config)   # jday only, no SST or SSS
+
+compare_runs(list(rf = rf, gam = gam, weak = weak))
+#>   comparison reference_score comparison_score difference   lower   upper
+#> 1        gam           0.854            0.892     0.0384  0.0141  0.0627
+#> 2       weak           0.854            0.654    -0.1997 -0.2964 -0.1030
+#>   p_value detectable n_stations
+#> 1 0.01177     0.0325        684
+#> 2 0.00458     0.1294        684
+```
+
+The first run in the list is the reference; every other is compared against it,
+and `difference` is *comparison minus reference*, so a positive number means the
+comparison won.
+
+**`detectable` is the column that stops a null result being misread.** It is the
+smallest true difference this comparison would have found, at 80% power — a
+property of the design rather than of the models. "Not significant" on its own
+conflates *these two models perform alike* with *this survey could not have told
+them apart*, and this separates them. A comparison that cannot see anything
+below 0.13 has not shown that a 0.03 gap is absent.
+
+The test is paired fold by fold, using the same Nadeau–Bengio correction as the
+[jackknife](#testing-which-covariates-earn-their-place) and for the same reason.
+It is two-sided here: leaving a covariate out has a direction worth testing
+against, but asking which of two models is better does not.
+
+Runs are matched on the station index, not on row order. Two runs with different
+covariates drop different stations to missingness, so the comparison uses the
+stations both actually scored and **warns** with the count — a run that drops
+many is telling you something.
+
+**How many stations would settle it?** `power_curve()` refits at several
+subsample sizes and traces power against `n`:
+
+```r
+power_curve(dat, list(full = config, starved = starved),
+            fractions = c(0.25, 0.5, 1), replicates = 3)
+#>   fraction n_stations replicates difference std_err df power detectable
+#> 1     0.25        171          3     -0.316  0.1169  4 0.428      0.434
+#> 2     0.50        342          3     -0.267  0.0893  4 0.632      0.332
+#> 3     1.00        684          3     -0.272  0.0339  4 1.000      0.126
+```
+
+The curve is the artefact, not any point on it. Power against sample size is
+steeply non-linear, and where a study sits on that curve is what decides whether
+another season of sampling is worth it: a comparison at 0.43 is one survey away
+from being decisive, one at 1.00 will not be improved by more stations.
+
+Everything is refitted at every size — a model trained on half the stations is a
+different model, not the same one evaluated on fewer — so this is the expensive
+function here, and it parallelises over the whole grid. Both runs see the same
+subsample and the same folds at every point, which is what keeps the comparison
+paired all the way down.
+
+Two honest limits. The target difference defaults to the one observed on the
+full data, which is an estimate rather than a truth; if the observed gap is
+itself mostly noise, the curve is answering a question about a size that may not
+be real. And these are differences in a bounded metric, so the normal-theory
+interval behind `detectable` is an approximation that gets worse as AUC
+approaches 1.
+
 ### Training and projection windows
 
 These are separate. Fitting on a long history and projecting a shorter or later
@@ -1210,6 +1284,7 @@ R/model.R               fit_patch_model()
 R/model_types.R         model_types(), permutation_importance()
 R/jackknife.R           jackknife_covariates(), the leave-one-out covariate test
 R/ensemble.R            fit_patch_ensemble(), combining several model types
+R/power.R               compare_runs(), power_curve()
 R/parallel.R            the worker pool both of those run on
 R/plot_effects.R        partial_effects(), glm_coefficients(), gam_smooth_terms()
 R/uncertainty.R         novelty_surface(), the projection interval
@@ -1438,7 +1513,13 @@ Earth](https://www.naturalearthdata.com/), public domain, via `rnaturalearth`.
 - Nadeau C, Bengio Y (2003). Inference for the generalization error. *Machine
   Learning* **52**(3), 239–281.
   [doi:10.1023/A:1024068626366](https://doi.org/10.1023/A:1024068626366) — the
-  variance correction behind the jackknife's `p_value`
+  variance correction behind the jackknife's `p_value` and `compare_runs()`
+- Hoenig JM, Heisey DM (2001). The abuse of power: the pervasive fallacy of
+  power calculations for data analysis. *The American Statistician* **55**(1),
+  19–24.
+  [doi:10.1198/000313001300339897](https://doi.org/10.1198/000313001300339897)
+  — why `compare_runs()` reports a minimum detectable difference rather than
+  the “observed power” it is often confused with
 - Fisher A, Rudin C, Dominici F (2019). All models are wrong, but many are useful:
   learning a variable's importance by studying an entire class of prediction
   models simultaneously. *Journal of Machine Learning Research* **20**(177), 1–81.
diff --git a/man/achieved_power.Rd b/man/achieved_power.Rd
new file mode 100644
index 0000000..76ef464
--- /dev/null
+++ b/man/achieved_power.Rd
@@ -0,0 +1,25 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/power.R
+\name{achieved_power}
+\alias{achieved_power}
+\title{The power a design has against a stated difference}
+\usage{
+achieved_power(difference, std_err, df, level = 0.95)
+}
+\arguments{
+\item{difference}{the true difference to detect}
+
+\item{std_err}{the corrected standard error of the difference}
+
+\item{df}{degrees of freedom}
+
+\item{level}{confidence level}
+}
+\value{
+a probability, or \code{NA_real_}
+}
+\description{
+The complement of \code{\link[=minimum_detectable]{minimum_detectable()}}: given how variable the folds were,
+how often would a true difference of \code{difference} be called significant?
+}
+\keyword{internal}
diff --git a/man/align_predictions.Rd b/man/align_predictions.Rd
new file mode 100644
index 0000000..cbd3d3e
--- /dev/null
+++ b/man/align_predictions.Rd
@@ -0,0 +1,20 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/power.R
+\name{align_predictions}
+\alias{align_predictions}
+\title{One run's predictions, restricted to shared stations and put in their order}
+\usage{
+align_predictions(predictions, shared)
+}
+\arguments{
+\item{predictions}{a run's held-out predictions}
+
+\item{shared}{the station indices to keep}
+}
+\value{
+the matching rows, in \code{shared} order
+}
+\description{
+One run's predictions, restricted to shared stations and put in their order
+}
+\keyword{internal}
diff --git a/man/compare_one.Rd b/man/compare_one.Rd
new file mode 100644
index 0000000..d22a1cf
--- /dev/null
+++ b/man/compare_one.Rd
@@ -0,0 +1,38 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/power.R
+\name{compare_one}
+\alias{compare_one}
+\title{Compare one run against the reference}
+\usage{
+compare_one(
+  reference_predictions,
+  other_predictions,
+  reference,
+  comparison,
+  metric = "roc_auc",
+  level = 0.95,
+  power = 0.8
+)
+}
+\arguments{
+\item{reference_predictions}{held-out predictions of the reference run}
+
+\item{other_predictions}{held-out predictions of the run being compared}
+
+\item{reference}{the reference run's name}
+
+\item{comparison}{the other run's name}
+
+\item{metric}{\code{"roc_auc"} or \code{"pr_auc"}}
+
+\item{level}{confidence level}
+
+\item{power}{the power \code{detectable} is computed at}
+}
+\value{
+a one-row data frame
+}
+\description{
+Compare one run against the reference
+}
+\keyword{internal}
diff --git a/man/compare_runs.Rd b/man/compare_runs.Rd
new file mode 100644
index 0000000..a60b90d
--- /dev/null
+++ b/man/compare_runs.Rd
@@ -0,0 +1,109 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/power.R
+\name{compare_runs}
+\alias{compare_runs}
+\title{Is the gap between two model runs real?}
+\usage{
+compare_runs(runs, metric = "roc_auc", level = 0.95, power = 0.8)
+}
+\arguments{
+\item{runs}{a named list of two or more fitted runs, from
+\code{\link[=fit_patch_model]{fit_patch_model()}} or \code{\link[=fit_patch_ensemble]{fit_patch_ensemble()}}. The first is the reference
+every other is compared against.}
+
+\item{metric}{\code{"roc_auc"} or \code{"pr_auc"}; both are threshold-free, which is
+what lets them be compared fold by fold without a cutoff moving underneath}
+
+\item{level}{confidence level for the interval and the test}
+
+\item{power}{the power \code{detectable} is computed at}
+}
+\value{
+a data frame with one row per comparison against the reference:
+\code{reference}, \code{comparison}, \code{metric}, \code{reference_score},
+\code{comparison_score}, \code{difference} (comparison minus reference), \code{lower},
+\code{upper}, \code{statistic}, \code{df}, \code{p_value}, \code{detectable}, \code{n_folds},
+\code{n_stations}, and \code{n_dropped}
+}
+\description{
+Two runs come back with two numbers — ROC AUC 0.857 against 0.871 — and
+nothing in either says whether the gap is a difference between the models or
+a difference between the stations the survey happened to visit. This answers
+that, and answers the question that should be asked next when the gap is not
+significant: \strong{how large would a difference have had to be before this study
+could have seen it?}
+}
+\details{
+Those two are reported together deliberately. "Not significant" on its own is
+the least informative result in modelling — it conflates \emph{these models
+perform alike} with \emph{this survey could not have told them apart}, and
+\code{detectable} is what separates the two. A run that cannot detect anything
+smaller than 0.09 in AUC has not shown that a 0.014 gap is absent.
+}
+\section{How the comparison is paired}{
+
+Every run cross-validates, so each carries a metric per fold rather than one
+number, and two runs on the same stations can be compared fold by fold. That
+pairing is most of the statistical power available: the folds vary a great
+deal between themselves and much less between two models scored on the \emph{same}
+fold, and an unpaired comparison throws that away.
+
+Runs are matched on \code{.row}, the station index, not on position. Two runs with
+different covariates drop different stations to missingness, so the
+comparison is made on the stations both actually scored and the number
+dropped is reported. A run that drops many is telling you something, which is
+why this warns rather than silently intersecting.
+}
+
+\section{The test, and why it is not a plain t-test}{
+
+The per-fold differences go through \code{\link[=corrected_paired_test]{corrected_paired_test()}}, the same
+Nadeau and Bengio (2003) correction the covariate jackknife uses, and for the
+same reason: any two cross-validation training sets share most of their rows,
+so folds are not independent and an uncorrected paired t-test finds
+significance that is not there.
+
+Two-sided here, unlike the jackknife. Leaving a covariate out has a direction
+worth testing against; asking which of two models is better does not.
+}
+
+\section{What "detectable" means}{
+
+The smallest true difference this comparison would have found significant at
+\code{level}, with probability \code{power}, given the fold-to-fold variability it
+actually saw:
+
+\deqn{d_{min} = SE \times (t_{1-\alpha/2, df} + t_{power, df})}
+
+It is a property of \strong{this} design — this many folds, these stations, this
+much variance between folds — not a general statement about the models. It
+says nothing about whether a smaller difference exists, only that this study
+would probably have missed it.
+}
+
+\examples{
+\dontrun{
+rf <- fit_patch_model(dat, within_config(config, type = "rf"))
+gam <- fit_patch_model(dat, within_config(config, type = "gam"))
+
+compare_runs(list(rf = rf, gam = gam))
+}
+}
+\references{
+Nadeau C, Bengio Y (2003). Inference for the generalization error. \emph{Machine
+Learning} \strong{52}(3), 239-281. \doi{10.1023/A:1024068626366} — the variance
+correction
+
+Dietterich TG (1998). Approximate statistical tests for comparing supervised
+classification learning algorithms. \emph{Neural Computation} \strong{10}(7),
+1895-1923. \doi{10.1162/089976698300017197} — why comparing learning
+algorithms on shared folds needs one
+
+Hoenig JM, Heisey DM (2001). The abuse of power: the pervasive fallacy of
+power calculations for data analysis. \emph{The American Statistician} \strong{55}(1),
+19-24. \doi{10.1198/000313001300339897} — why \code{detectable} is reported
+rather than the observed-power statistic it is often confused with
+}
+\seealso{
+\code{\link[=power_curve]{power_curve()}} for how the answer changes with more stations
+}
diff --git a/man/corrected_paired_test.Rd b/man/corrected_paired_test.Rd
index ca885c1..523d1f2 100644
--- a/man/corrected_paired_test.Rd
+++ b/man/corrected_paired_test.Rd
@@ -4,10 +4,13 @@
 \alias{corrected_paired_test}
 \title{A paired test across folds, with the Nadeau-Bengio variance correction}
 \usage{
-corrected_paired_test(differences)
+corrected_paired_test(differences, alternative = c("greater", "two.sided"))
 }
 \arguments{
 \item{differences}{per-fold score of the full model minus the reduced one}
+
+\item{alternative}{\code{"greater"} for a directional hypothesis, \code{"two.sided"}
+when either sign is a finding}
 }
 \value{
 a list of \code{estimate}, \code{std_err}, \code{statistic}, \code{df}, \code{p_value}, \code{n}
@@ -21,9 +24,11 @@ test-set to training-set size in k-fold — which is the standard correction
 and costs roughly a factor of \code{sqrt(2)} off the statistic.
 }
 \details{
-One-sided, because the hypothesis is directional: the question is whether
-removing the covariate makes the model \emph{worse}, and a covariate whose removal
-improves the model has failed the test rather than passed a different one.
+The default is one-sided, because the jackknife's hypothesis is directional:
+the question is whether removing the covariate makes the model \emph{worse}, and a
+covariate whose removal improves the model has failed the test rather than
+passed a different one. Comparing two models is not directional in that way -
+either may be the better - so \code{\link[=compare_runs]{compare_runs()}} asks for \code{two.sided}.
 }
 \references{
 Nadeau C, Bengio Y (2003). Inference for the generalization error. \emph{Machine
diff --git a/man/critical_t.Rd b/man/critical_t.Rd
new file mode 100644
index 0000000..1b1dddb
--- /dev/null
+++ b/man/critical_t.Rd
@@ -0,0 +1,20 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/power.R
+\name{critical_t}
+\alias{critical_t}
+\title{The two-sided critical value}
+\usage{
+critical_t(level, df)
+}
+\arguments{
+\item{level}{confidence level}
+
+\item{df}{degrees of freedom}
+}
+\value{
+the critical \code{t}, or \code{NA_real_}
+}
+\description{
+The two-sided critical value
+}
+\keyword{internal}
diff --git a/man/minimum_detectable.Rd b/man/minimum_detectable.Rd
new file mode 100644
index 0000000..f54fbb9
--- /dev/null
+++ b/man/minimum_detectable.Rd
@@ -0,0 +1,33 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/power.R
+\name{minimum_detectable}
+\alias{minimum_detectable}
+\title{The smallest difference a design could have detected}
+\usage{
+minimum_detectable(std_err, df, level = 0.95, power = 0.8)
+}
+\arguments{
+\item{std_err}{the corrected standard error of the difference}
+
+\item{df}{degrees of freedom}
+
+\item{level}{confidence level}
+
+\item{power}{the power to solve at}
+}
+\value{
+the smallest detectable difference, or \code{NA_real_}
+}
+\description{
+Reported instead of "observed power", which is the statistic this is usually
+confused with and which carries no information a p-value does not — it is a
+deterministic function of it (Hoenig and Heisey 2001). The minimum detectable
+difference is about the \emph{design} rather than about the result, which is what
+makes it worth reading beside a null finding.
+}
+\references{
+Hoenig JM, Heisey DM (2001). The abuse of power: the pervasive fallacy of
+power calculations for data analysis. \emph{The American Statistician} \strong{55}(1),
+19-24. \doi{10.1198/000313001300339897}
+}
+\keyword{internal}
diff --git a/man/power_curve.Rd b/man/power_curve.Rd
new file mode 100644
index 0000000..82bfdbc
--- /dev/null
+++ b/man/power_curve.Rd
@@ -0,0 +1,109 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/power.R
+\name{power_curve}
+\alias{power_curve}
+\title{How the comparison would improve with more stations}
+\usage{
+power_curve(
+  dat,
+  configs,
+  fractions = c(0.25, 0.5, 0.75, 1),
+  replicates = 5,
+  difference = NULL,
+  metric = "roc_auc",
+  level = 0.95,
+  workers = NULL,
+  seed = 42
+)
+}
+\arguments{
+\item{dat}{labeled modeling data from \code{label_patch()} with covariates attached}
+
+\item{configs}{a named list of two or more configs to compare. The first is
+the reference; each is fitted exactly as its own run would be.}
+
+\item{fractions}{the shares of the stations to fit at}
+
+\item{replicates}{how many subsamples per fraction; the spread across them is
+what stops one unlucky draw from setting a point on the curve}
+
+\item{difference}{the true difference to compute power against; \code{NULL} uses
+the one observed at the largest fraction}
+
+\item{metric}{\code{"roc_auc"} or \code{"pr_auc"}}
+
+\item{level}{confidence level the test would use}
+
+\item{workers}{how many workers; see \code{\link[=resolve_workers]{resolve_workers()}}}
+
+\item{seed}{a seed, so a curve is reproducible}
+}
+\value{
+a data frame with one row per fraction: \code{fraction}, \code{n_stations},
+\code{replicates}, \code{difference} (mean observed), \code{std_err}, \code{df}, \code{power}, and
+\code{detectable}. The target difference is on it as a \code{difference} attribute.
+}
+\description{
+\code{\link[=compare_runs]{compare_runs()}} answers what this survey could see. This answers what a
+larger one would: it refits both runs on subsamples of the stations, at
+several sizes, and traces how the power to detect a difference grows with
+\code{n}.
+}
+\details{
+The curve is the useful artefact rather than any single number on it. Power
+against sample size is steeply non-linear, and where a study sits on that
+curve decides what the next survey is worth: a comparison at 0.35 power is
+one more season away from being decisive, and one at 0.9 will not be improved
+by more stations because it is already there.
+}
+\section{What it costs, and what it therefore skips}{
+
+\verb{fractions × replicates × runs} cross-validations. Everything is refitted at
+every size — the point is precisely that a model trained on half the stations
+is a different model, not the same model evaluated on fewer — so this is the
+expensive function in the package and parallelises over the whole grid.
+
+It goes through the same fold-scoring path the covariate jackknife uses,
+which fits and scores and stops there. No bootstrap intervals, no variable
+importance, no projection: none of it enters the curve, and all of it would
+be paid for at every point.
+}
+
+\section{Subsampling stations, not folds}{
+
+Rows are drawn without replacement, and the folds are then built inside each
+subsample. Reusing the full run's folds and thinning them would shrink the
+held-out sets while leaving the training sets nearly whole, which measures
+something else entirely — the curve has to come from models that were
+actually trained on less.
+
+Both runs see the \strong{same} subsample and the \strong{same} folds at every point,
+which is what keeps the comparison paired all the way down the curve.
+}
+
+\section{Reading it honestly}{
+
+The target difference defaults to the one observed on the full data, and that
+is an estimate, not a truth. If the observed gap is itself mostly noise, the
+curve answers "how many stations to reliably detect a difference this size"
+for a size that may not be real. It is a projection under an assumption, and
+the assumption is the observed effect.
+}
+
+\examples{
+\dontrun{
+rf <- config; rf$model$type <- "rf"
+gam <- config; gam$model$type <- "gam"
+
+curve <- power_curve(dat, list(rf = rf, gam = gam))
+curve[c("n_stations", "power", "detectable")]
+}
+}
+\references{
+Nadeau C, Bengio Y (2003). Inference for the generalization error. \emph{Machine
+Learning} \strong{52}(3), 239-281. \doi{10.1023/A:1024068626366}
+}
+\seealso{
+\code{\link[=compare_runs]{compare_runs()}}, which answers the same question for the data you
+already have
+}
diff --git a/man/power_point.Rd b/man/power_point.Rd
new file mode 100644
index 0000000..7aac54f
--- /dev/null
+++ b/man/power_point.Rd
@@ -0,0 +1,41 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/power.R
+\name{power_point}
+\alias{power_point}
+\title{One point on the power curve}
+\usage{
+power_point(
+  model_data,
+  configs,
+  predictors,
+  types,
+  fraction,
+  replicate,
+  metric = "roc_auc",
+  seed = 42
+)
+}
+\arguments{
+\item{model_data}{the complete-case modeling data}
+
+\item{configs}{the configs being compared}
+
+\item{predictors}{each config's predictors}
+
+\item{types}{each config's model type}
+
+\item{fraction}{the share of stations to draw}
+
+\item{replicate}{which draw this is}
+
+\item{metric}{\code{"roc_auc"} or \code{"pr_auc"}}
+
+\item{seed}{the run's seed}
+}
+\value{
+a one-row data frame, or \code{NULL} when the draw could not be scored
+}
+\description{
+Draws a subsample, folds it, and scores every config on those same folds.
+}
+\keyword{internal}
diff --git a/man/run_predictions.Rd b/man/run_predictions.Rd
new file mode 100644
index 0000000..b40c876
--- /dev/null
+++ b/man/run_predictions.Rd
@@ -0,0 +1,18 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/power.R
+\name{run_predictions}
+\alias{run_predictions}
+\title{The held-out predictions a run kept, checked}
+\usage{
+run_predictions(run)
+}
+\arguments{
+\item{run}{a \code{\link[=fit_patch_model]{fit_patch_model()}} or \code{\link[=fit_patch_ensemble]{fit_patch_ensemble()}} result}
+}
+\value{
+the predictions data frame
+}
+\description{
+The held-out predictions a run kept, checked
+}
+\keyword{internal}
diff --git a/man/yaml_bool_marker.Rd b/man/yaml_bool_marker.Rd
index 620e0a0..e77af26 100644
--- a/man/yaml_bool_marker.Rd
+++ b/man/yaml_bool_marker.Rd
@@ -7,13 +7,19 @@
 yaml_bool_marker
 }
 \description{
-A control character, so it cannot collide with anything a YAML file could
-legitimately contain — including a quoted string that was meant to be the
-text \code{"true"}, which must survive as that text and not become a logical.
-}
-\details{
-A string rather than an attribute or a class because it has to survive
+A prefix rather than an attribute or a class, because it has to survive
 \code{yaml} collapsing a sequence of scalars into an atomic vector, which drops
 attributes. \verb{[true, false]} would otherwise come back as two strings.
 }
+\details{
+Deliberately plain ASCII. The first version used control characters, on the
+reasoning that nothing could collide with them. That was true, and it made
+\code{file(1)} report the whole of \code{config.R} as binary rather than as source, so
+editors and diff viewers presented the file as corrupt.
+
+Spelling it out costs nothing. Only the boolean handlers ever prepend this,
+and they only ever see scalars YAML itself resolved as booleans — a quoted
+\code{"true"} carries no boolean tag and is never marked. The one way left to
+collide is a config value that genuinely begins with this text.
+}
 \keyword{internal}
diff --git a/tests/testthat/test-power.R b/tests/testthat/test-power.R
new file mode 100644
index 0000000..3e32779
--- /dev/null
+++ b/tests/testthat/test-power.R
@@ -0,0 +1,213 @@
+# A run with only what compare_runs() reads: held-out predictions carrying the
+# station index, the fold, the truth and the probability. Building these by hand
+# is what lets the pairing and the arithmetic be tested without fitting
+# anything, and lets a "better" run be better by construction rather than by
+# luck.
+stub_run <- function(skill = 0.3, n = 200, folds = 5, seed = 1, rows = NULL) {
+  set.seed(seed)
+  rows <- rows %||% seq_len(n)
+  n <- length(rows)
+  is_patch <- rep(c(TRUE, FALSE), length.out = n)[order(stats::runif(n))]
+
+  # A larger `skill` separates the classes further, so ROC AUC rises with it.
+  probability <- stats::plogis(ifelse(is_patch, skill, -skill) +
+                                 stats::rnorm(n, sd = 0.5))
+  list(predictions = data.frame(
+    .row = rows,
+    id = paste0("Fold", rep(seq_len(folds), length.out = n)),
+    patch = factor(ifelse(is_patch, "patch", "non_patch"),
+                   levels = c("patch", "non_patch")),
+    .pred_patch = probability,
+    stringsAsFactors = FALSE
+  ))
+}
+
+test_that("compare_runs refuses what it cannot compare", {
+  run <- stub_run()
+
+  expect_error(compare_runs(list(a = run)), "at least 2 fitted runs")
+  expect_error(compare_runs(list(run, run)), "Name the runs")
+  expect_error(compare_runs(list(a = run, b = run), metric = "tss"),
+               "must be 'roc_auc' or 'pr_auc'")
+  # A model fitted outside the package has no held-out predictions to read.
+  expect_error(compare_runs(list(a = run, b = list(predictions = NULL))),
+               "no usable held-out predictions")
+})
+
+test_that("a run compared with itself shows no difference", {
+  run <- stub_run()
+
+  out <- compare_runs(list(a = run, b = run))
+
+  expect_equal(nrow(out), 1)
+  expect_equal(out$difference, 0)
+  expect_equal(out$p_value, 1)
+  expect_equal(out$reference_score, out$comparison_score)
+  expect_equal(out$n_dropped, 0)
+})
+
+test_that("a genuinely better run is found to be better", {
+  # Better by construction: the same stations, the same folds, a wider
+  # separation between the classes.
+  worse <- stub_run(skill = 0.2, seed = 3)
+  better <- stub_run(skill = 2.0, seed = 3)
+
+  out <- compare_runs(list(worse = worse, better = better))
+
+  expect_gt(out$difference, 0)
+  expect_lt(out$p_value, 0.05)
+  expect_gt(out$comparison_score, out$reference_score)
+  # The interval excludes zero, which is the same statement the p-value makes.
+  expect_gt(out$lower, 0)
+})
+
+test_that("runs are paired on the station index, not on row order", {
+  # tune returns folds in its own order, and two runs need not agree on it.
+  # Pairing by position would silently compare station i of one run with a
+  # different station of the other.
+  run <- stub_run(skill = 1.2, seed = 7)
+  shuffled <- run
+  shuffled$predictions <- run$predictions[order(stats::runif(nrow(run$predictions))), ]
+
+  out <- compare_runs(list(a = run, b = shuffled))
+
+  expect_equal(out$difference, 0)
+  expect_equal(out$reference_score, out$comparison_score)
+})
+
+test_that("runs covering different stations are intersected, loudly", {
+  # Different covariates drop different stations to missingness. Comparing on
+  # what they share is right; doing it silently is not.
+  full <- stub_run(n = 200, seed = 11, rows = 1:200)
+  partial <- stub_run(n = 160, seed = 11, rows = 1:160)
+
+  expect_warning(out <- compare_runs(list(full = full, partial = partial)),
+                 "do not cover the same stations")
+  expect_equal(out$n_stations, 160)
+  expect_equal(out$n_dropped, 40)
+})
+
+test_that("two runs sharing almost nothing are refused rather than compared", {
+  a <- stub_run(n = 100, rows = 1:100)
+  b <- stub_run(n = 100, rows = 500:599)
+
+  expect_error(suppressWarnings(compare_runs(list(a = a, b = b))),
+               "share 0 stations")
+})
+
+test_that("every run after the first is compared against the first", {
+  runs <- list(ref = stub_run(seed = 2), one = stub_run(seed = 2),
+               two = stub_run(seed = 2))
+
+  out <- compare_runs(runs)
+
+  expect_equal(nrow(out), 2)
+  expect_equal(out$comparison, c("one", "two"))
+  expect_true(all(out$reference == "ref"))
+})
+
+test_that("the minimum detectable difference behaves like one", {
+  # Bigger when the folds disagree more, smaller with more folds, and larger
+  # for more power - each of which is the direction that makes it usable.
+  expect_gt(minimum_detectable(0.04, df = 4), minimum_detectable(0.02, df = 4))
+  expect_gt(minimum_detectable(0.03, df = 2), minimum_detectable(0.03, df = 20))
+  expect_gt(minimum_detectable(0.03, df = 9, power = 0.9),
+            minimum_detectable(0.03, df = 9, power = 0.8))
+
+  # The formula it claims to be.
+  se <- 0.03; df <- 9
+  expect_equal(minimum_detectable(se, df, level = 0.95, power = 0.8),
+               se * (stats::qt(0.975, df) + stats::qt(0.8, df)))
+
+  expect_true(is.na(minimum_detectable(NA_real_, 4)))
+  expect_true(is.na(minimum_detectable(0.03, NA_real_)))
+})
+
+test_that("power rises with the difference and falls with the noise", {
+  expect_gt(achieved_power(0.10, 0.02, df = 9), achieved_power(0.02, 0.02, df = 9))
+  expect_gt(achieved_power(0.05, 0.01, df = 9), achieved_power(0.05, 0.05, df = 9))
+  expect_true(all(achieved_power(c(0, 0.5), 0.02, df = 9) %in% c(0, 1) |
+                    (achieved_power(c(0, 0.5), 0.02, df = 9) >= 0 &
+                       achieved_power(c(0, 0.5), 0.02, df = 9) <= 1)))
+
+  # At exactly the minimum detectable difference, power is the power it was
+  # solved for. The two functions are inverses and must agree.
+  se <- 0.03; df <- 9
+  mde <- minimum_detectable(se, df, level = 0.95, power = 0.8)
+  expect_equal(achieved_power(mde, se, df, level = 0.95), 0.8, tolerance = 0.02)
+})
+
+test_that("the paired test is two-sided when asked", {
+  differences <- c(0.03, 0.01, 0.04, 0.02, 0.025)
+
+  one <- corrected_paired_test(differences)
+  two <- corrected_paired_test(differences, alternative = "two.sided")
+
+  expect_equal(two$p_value, 2 * one$p_value)
+  expect_equal(two$estimate, one$estimate)
+  # A difference in the other direction is a finding for two-sided and not for
+  # one-sided, which is the whole reason compare_runs() asks for it.
+  flipped <- corrected_paired_test(-differences, alternative = "two.sided")
+  expect_equal(flipped$p_value, two$p_value)
+  expect_gt(corrected_paired_test(-differences)$p_value, 0.5)
+})
+
+test_that("power_curve refuses what it cannot trace", {
+  config <- mock_config()
+
+  expect_error(power_curve(data.frame(), list(a = config)), "at least 2 configs")
+  expect_error(power_curve(data.frame(), list(config, config)),
+               "Name the configs")
+})
+
+test_that("power_curve traces power rising with the number of stations", {
+  skip_on_cran()
+  config <- mock_config()
+  config$model$trees <- 50
+  config$model$cv_folds <- 5
+  dat <- labeled_mock_data(config)
+
+  # A covariate-starved run against the full one, so there is a real difference
+  # for the curve to have power against.
+  starved <- config
+  starved$covariates$exclude <- c("SST", "SSS")
+
+  curve <- suppressMessages(
+    power_curve(dat, list(full = config, starved = starved),
+                fractions = c(0.25, 1), replicates = 2, workers = 1)
+  )
+
+  expect_equal(nrow(curve), 2)
+  expect_true(all(c("fraction", "n_stations", "power", "detectable",
+                    "difference", "std_err") %in% names(curve)))
+  # More stations, more power and a smaller detectable difference. This is the
+  # shape the whole function exists to produce.
+  expect_gt(curve$n_stations[2], curve$n_stations[1])
+  expect_gte(curve$power[2], curve$power[1])
+  expect_lt(curve$detectable[2], curve$detectable[1])
+  expect_true(all(curve$power >= 0 & curve$power <= 1))
+  # The starved run is worse, so the difference is negative throughout.
+  expect_true(all(curve$difference < 0))
+  expect_true(is.numeric(attr(curve, "difference")))
+})
+
+test_that("a power curve is reproducible", {
+  skip_on_cran()
+  config <- mock_config()
+  config$model$trees <- 50
+  config$model$cv_folds <- 5
+  dat <- labeled_mock_data(config)
+  starved <- config
+  starved$covariates$exclude <- "SST"
+
+  runs <- list(full = config, starved = starved)
+  once <- suppressMessages(power_curve(dat, runs, fractions = 0.5,
+                                        replicates = 2, workers = 1))
+  twice <- suppressMessages(power_curve(dat, runs, fractions = 0.5,
+                                         replicates = 2, workers = 1))
+
+  # Seeded per point rather than per call, so the answer does not depend on how
+  # the tasks happened to be spread across workers.
+  expect_equal(once$difference, twice$difference)
+  expect_equal(once$std_err, twice$std_err)
+})
diff --git a/tools/citations.csv b/tools/citations.csv
index db3ef07..3ece160 100644
--- a/tools/citations.csv
+++ b/tools/citations.csv
@@ -40,3 +40,4 @@ marmion2009,10.1111/j.1472-4642.2008.00491.x,,Marmion,2009,Evaluation of consens
 elith2011,10.1111/j.1472-4642.2010.00725.x,,Elith,2011,A statistical explanation of MaxEnt for ecologists,Diversity and Distributions,17,43-57,crossref,
 nadeau2003,10.1023/A:1024068626366,,Nadeau,2003,Inference for the Generalization Error,Machine Learning,52,239-281,crossref,
 dietterich1998,10.1162/089976698300017197,,Dietterich,1998,Approximate Statistical Tests for Comparing Supervised Classification Learning Algorithms,Neural Computation,10,1895-1923,crossref,
+hoenig2001,10.1198/000313001300339897,,Hoenig,2001,The Abuse of Power,The American Statistician,55,19-24,crossref,