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,