From 354ceb76a16ce00b42f3b019479d8ebe6a6766c2 Mon Sep 17 00:00:00 2001 From: chross22 <52218551+chross22@users.noreply.github.com> Date: Mon, 10 Aug 2026 11:29:11 -0400 Subject: [PATCH 01/10] Read a YAML key that spells a boolean as the key it is YAML 1.1, which is what yaml::read_yaml() parses, reads a bare `n` as the boolean false - along with y, yes, no, on, and off. That is right for a value and wrong for a key, and the difference is silent. A derivoce step written - type: lag_covariate vars: [CHL] n: 2 parses to a list whose key is named FALSE, so spec$n is NULL and R/derivoce.R falls back to a one-month lag. The config asked for two, the run used one, and nothing said so. cfin_gom.yaml hid it by asking for the value the fallback already gives. yaml's handlers are given the scalar as it was written - "n", not FALSE - so the source text is recoverable. It cannot be resolved while parsing, though, because a handler cannot see whether it is being called for a key or for a value. So this marks each boolean scalar with its source text and then walks the finished structure: text in a name position is the key the file wrote, text in a value position becomes the logical it meant. The marker is a string prefix rather than an attribute or a class because it has to survive yaml collapsing a sequence of scalars into an atomic vector, which drops attributes. [true, false] would otherwise come back as two strings. Only the read side was affected. yaml::as.yaml() already quotes these keys, so a config generate_config() or save_config() wrote has always been read correctly by anything. cfin_gom.yaml and the README example are quoted anyway, for readers that are not this package. Co-Authored-By: Claude Opus 5 --- R/config.R | 93 +++++++++++++++++++++++++++++++++++- README.md | 9 +++- inst/configs/cfin_gom.yaml | 6 ++- man/read_config_yaml.Rd | 43 +++++++++++++++++ man/restore_yaml_bools.Rd | 21 ++++++++ man/unmark_yaml_bool.Rd | 19 ++++++++ man/yaml_bool_marker.Rd | 19 ++++++++ tests/testthat/test-config.R | 86 +++++++++++++++++++++++++++++++++ 8 files changed, 293 insertions(+), 3 deletions(-) create mode 100644 man/read_config_yaml.Rd create mode 100644 man/restore_yaml_bools.Rd create mode 100644 man/unmark_yaml_bool.Rd create mode 100644 man/yaml_bool_marker.Rd diff --git a/R/config.R b/R/config.R index 3ba2bee..9569711 100644 --- a/R/config.R +++ b/R/config.R @@ -14,7 +14,7 @@ load_config <- function(path) { if (!file.exists(path)) { stop("Config file not found: ", path, call. = FALSE) } - config <- yaml::read_yaml(path) + config <- read_config_yaml(path) config <- resolve_config_paths(config, path) config <- apply_config_defaults(config) @@ -37,6 +37,97 @@ load_config <- function(path) { config } +#' Read a config YAML without losing keys that spell a boolean +#' +#' `yaml::read_yaml()` parses YAML 1.1, where a bare `n` is the boolean `false`. +#' That is correct for a *value* and wrong for a *key*, and the difference is +#' silent: a derivoce step written +#' +#' ```yaml +#' - type: lag_covariate +#' vars: [CHL] +#' n: 2 +#' ``` +#' +#' parses to a list whose key is named `FALSE`, so `spec$n` is `NULL` and the +#' step falls back to a one-month lag. The config asked for two, the run used +#' one, and nothing said so. The same goes for `y`, `yes`, `no`, `on`, `off`, +#' `true` and `false` in any capitalisation. +#' +#' The fix is to keep the source text. `yaml`'s handlers are given the original +#' scalar as it was written — `"n"`, not `FALSE` — so this marks each one and +#' then, once the structure exists, restores it: text in a name position is the +#' key the file actually wrote, and text in a value position becomes the logical +#' it meant. The two cannot be told apart while parsing, which is exactly why +#' this is two passes rather than a cleverer handler. +#' +#' @param path path to a config YAML file +#' @return the parsed config list +#' @seealso [write_config_yaml()], the write side — `yaml::as.yaml()` already +#' quotes these keys on the way out, so a config this package writes is read +#' correctly by anything +#' @keywords internal +read_config_yaml <- function(path) { + mark <- function(x) paste0(yaml_bool_marker, x) + restore_yaml_bools( + yaml::read_yaml(path, handlers = list("bool#yes" = mark, "bool#no" = mark)) + ) +} + +#' The marker that carries a boolean's source text through parsing +#' +#' 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 `"true"`, which must survive as that text and not become a logical. +#' +#' A string rather than an attribute or a class because it has to survive +#' `yaml` collapsing a sequence of scalars into an atomic vector, which drops +#' attributes. `[true, false]` would otherwise come back as two strings. +#' +#' @keywords internal +yaml_bool_marker <- "taupatch-bool" + +#' Turn marked scalars back into keys and logicals +#' +#' Names get the source text they were written with; values get the logical +#' that text meant. A sequence that mixes a marked scalar with an unmarked one +#' cannot be a logical vector, so it keeps the text — `[true, maybe]` is a list +#' of two strings, which is the only reading available. +#' +#' @param x a parsed YAML value +#' @return `x` with markers resolved +#' @keywords internal +restore_yaml_bools <- function(x) { + if (is.list(x)) { + names(x) <- unmark_yaml_bool(names(x)) + return(lapply(x, restore_yaml_bools)) + } + if (!is.character(x) || length(x) == 0) return(x) + + marked <- startsWith(x, yaml_bool_marker) + if (!any(marked)) return(x) + + source <- unmark_yaml_bool(x) + if (!all(marked)) return(source) + # The YAML 1.1 spellings of true. Everything else the handlers see is a false. + tolower(source) %in% c("y", "yes", "true", "on") +} + +#' Strip the boolean marker, leaving what the file wrote +#' +#' Anchored at the start rather than replaced wherever it appears, so a quoted +#' value that happens to contain the marker's text further along is left alone. +#' +#' @param x a character vector, or `NULL` for an unnamed list +#' @return `x` with any marker prefix removed +#' @keywords internal +unmark_yaml_bool <- function(x) { + if (is.null(x)) return(NULL) + marked <- startsWith(x, yaml_bool_marker) + x[marked] <- substring(x[marked], nchar(yaml_bool_marker) + 1L) + x +} + #' Resolve config paths relative to the project directory #' #' `paths.project_dir` itself is resolved relative to the config file's own diff --git a/README.md b/README.md index d5ce538..864941e 100644 --- a/README.md +++ b/README.md @@ -478,7 +478,7 @@ covariates: bottom: BOTT # SST_BOTT_vgrad, the stratification index - type: lag_covariate vars: [CHL] - n: 1 # CHL_lag1 + "n": 1 # CHL_lag1 - see the note on quoting below - type: integrate_covariate vars: [CHL] window: year # CHL_int, the original pipeline's int_chl @@ -488,6 +488,13 @@ covariates: - distance_to_shore # shore_dist, its dist ``` +**Why `"n"` is quoted.** YAML 1.1 — which is what `yaml::read_yaml()` parses — +reads a bare `n` as the boolean `false`, along with `y`, `yes`, `no`, `on` and +`off`. So `n: 2` names the key `FALSE` rather than `n`, and a plain parser hands +back a step with no `n` at all, which silently falls back to a one-month lag. +`load_config()` recovers the spelling, so an unquoted `n: 2` does work here. The +quotes are for everything else that might open the file. + Steps run in order and see the columns earlier ones produced. That is why `current_speed` followed by a gradient of `speed` works. `distance_to_front`, `distance_to_contour`, `distance_to_isobath`, `ftle`, and `fsle` are available diff --git a/inst/configs/cfin_gom.yaml b/inst/configs/cfin_gom.yaml index edf7a92..f7a4c2d 100644 --- a/inst/configs/cfin_gom.yaml +++ b/inst/configs/cfin_gom.yaml @@ -131,7 +131,11 @@ covariates: bottom: BOTT # SST_BOTT_vgrad; the stratification index - type: lag_covariate vars: [CHL] - n: 1 # CHL_lag1; last month's bloom feeds this month + # Quoted, and it has to be. YAML 1.1 reads a bare `n` as the boolean + # false, so `n: 1` would name this key FALSE and the step would silently + # fall back to a one-month lag. taupatch recovers the spelling on read, + # but quoting keeps the file correct for anything else that opens it. + "n": 1 # CHL_lag1; last month's bloom feeds this month - type: integrate_covariate vars: [CHL] window: year # CHL_int; the original pipeline's int_chl diff --git a/man/read_config_yaml.Rd b/man/read_config_yaml.Rd new file mode 100644 index 0000000..9a0249b --- /dev/null +++ b/man/read_config_yaml.Rd @@ -0,0 +1,43 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/config.R +\name{read_config_yaml} +\alias{read_config_yaml} +\title{Read a config YAML without losing keys that spell a boolean} +\usage{ +read_config_yaml(path) +} +\arguments{ +\item{path}{path to a config YAML file} +} +\value{ +the parsed config list +} +\description{ +\code{yaml::read_yaml()} parses YAML 1.1, where a bare \code{n} is the boolean \code{false}. +That is correct for a \emph{value} and wrong for a \emph{key}, and the difference is +silent: a derivoce step written +} +\details{ +\if{html}{\out{
}}\preformatted{- type: lag_covariate + vars: [CHL] + n: 2 +}\if{html}{\out{
}} + +parses to a list whose key is named \code{FALSE}, so \code{spec$n} is \code{NULL} and the +step falls back to a one-month lag. The config asked for two, the run used +one, and nothing said so. The same goes for \code{y}, \code{yes}, \code{no}, \code{on}, \code{off}, +\code{true} and \code{false} in any capitalisation. + +The fix is to keep the source text. \code{yaml}'s handlers are given the original +scalar as it was written — \code{"n"}, not \code{FALSE} — so this marks each one and +then, once the structure exists, restores it: text in a name position is the +key the file actually wrote, and text in a value position becomes the logical +it meant. The two cannot be told apart while parsing, which is exactly why +this is two passes rather than a cleverer handler. +} +\seealso{ +\code{\link[=write_config_yaml]{write_config_yaml()}}, the write side — \code{yaml::as.yaml()} already +quotes these keys on the way out, so a config this package writes is read +correctly by anything +} +\keyword{internal} diff --git a/man/restore_yaml_bools.Rd b/man/restore_yaml_bools.Rd new file mode 100644 index 0000000..60c05f2 --- /dev/null +++ b/man/restore_yaml_bools.Rd @@ -0,0 +1,21 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/config.R +\name{restore_yaml_bools} +\alias{restore_yaml_bools} +\title{Turn marked scalars back into keys and logicals} +\usage{ +restore_yaml_bools(x) +} +\arguments{ +\item{x}{a parsed YAML value} +} +\value{ +\code{x} with markers resolved +} +\description{ +Names get the source text they were written with; values get the logical +that text meant. A sequence that mixes a marked scalar with an unmarked one +cannot be a logical vector, so it keeps the text — \verb{[true, maybe]} is a list +of two strings, which is the only reading available. +} +\keyword{internal} diff --git a/man/unmark_yaml_bool.Rd b/man/unmark_yaml_bool.Rd new file mode 100644 index 0000000..fd9f502 --- /dev/null +++ b/man/unmark_yaml_bool.Rd @@ -0,0 +1,19 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/config.R +\name{unmark_yaml_bool} +\alias{unmark_yaml_bool} +\title{Strip the boolean marker, leaving what the file wrote} +\usage{ +unmark_yaml_bool(x) +} +\arguments{ +\item{x}{a character vector, or \code{NULL} for an unnamed list} +} +\value{ +\code{x} with any marker prefix removed +} +\description{ +Anchored at the start rather than replaced wherever it appears, so a quoted +value that happens to contain the marker's text further along is left alone. +} +\keyword{internal} diff --git a/man/yaml_bool_marker.Rd b/man/yaml_bool_marker.Rd new file mode 100644 index 0000000..620e0a0 --- /dev/null +++ b/man/yaml_bool_marker.Rd @@ -0,0 +1,19 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/config.R +\name{yaml_bool_marker} +\alias{yaml_bool_marker} +\title{The marker that carries a boolean's source text through parsing} +\usage{ +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 +\code{yaml} collapsing a sequence of scalars into an atomic vector, which drops +attributes. \verb{[true, false]} would otherwise come back as two strings. +} +\keyword{internal} diff --git a/tests/testthat/test-config.R b/tests/testthat/test-config.R index d55cca3..190b0ab 100644 --- a/tests/testthat/test-config.R +++ b/tests/testthat/test-config.R @@ -24,6 +24,92 @@ test_that("column_prefix defaults to the species key but can alias it", { expect_equal(resolve_species(config)$column_prefix, "newsp") }) +# A config file on disk, so the YAML parser is actually exercised. Building the +# list in R skips the only place these bugs can happen. +write_yaml_text <- function(lines) { + path <- tempfile(fileext = ".yaml") + writeLines(lines, path) + path +} + +test_that("a key that spells a boolean keeps its own name", { + # YAML 1.1 reads a bare `n` as false, so `n: 2` would name the key FALSE and + # a lag_covariate step would silently fall back to one month. + path <- write_yaml_text(c("steps:", " - type: lag_covariate", " n: 2")) + + parsed <- read_config_yaml(path) + + expect_equal(names(parsed$steps[[1]]), c("type", "n")) + expect_equal(parsed$steps[[1]]$n, 2) +}) + +test_that("every YAML 1.1 boolean spelling survives as a key", { + path <- write_yaml_text(c("n: 1", "y: 2", "no: 3", "yes: 4", "off: 5", + "on: 6", "true: 7", "false: 8", "N: 9", "Off: 10")) + + parsed <- read_config_yaml(path) + + expect_equal(names(parsed), c("n", "y", "no", "yes", "off", "on", "true", + "false", "N", "Off")) + expect_equal(unname(unlist(parsed)), 1:10) +}) + +test_that("a boolean in a value position is still a boolean", { + # The other half of the fix. Recovering key spellings must not turn the + # config's actual switches into strings. + path <- write_yaml_text(c("model:", " tune: false", " select: yes", + "covariates:", " normalize: true", " thin: off", + "flags: [true, false]")) + + parsed <- read_config_yaml(path) + + expect_identical(parsed$model$tune, FALSE) + expect_identical(parsed$model$select, TRUE) + expect_identical(parsed$covariates$normalize, TRUE) + expect_identical(parsed$covariates$thin, FALSE) + # A sequence of them collapses to an atomic vector, which drops attributes - + # which is why the source text is carried in the string rather than beside it. + expect_identical(parsed$flags, c(TRUE, FALSE)) +}) + +test_that("a quoted string that spells a boolean stays a string", { + path <- write_yaml_text(c("species:", " active: 'true'", " label: \"no\"", + "note: not a bool")) + + parsed <- read_config_yaml(path) + + expect_identical(parsed$species$active, "true") + expect_identical(parsed$species$label, "no") + expect_identical(parsed$note, "not a bool") +}) + +test_that("the shipped config's lag reaches derivoce with the lag it asks for", { + # The end-to-end version: the bug was invisible in cfin_gom.yaml precisely + # because its `n: 1` matched the fallback, so this pins the whole path. + path <- system.file("configs", "cfin_gom.yaml", package = "taupatch") + if (!nzchar(path)) path <- test_path("..", "..", "inst", "configs", "cfin_gom.yaml") + + config <- read_config_yaml(path) + lag <- Filter(function(s) identical(s$type, "lag_covariate"), + config$covariates$derivoce)[[1]] + + expect_equal(lag$n, 1) + expect_false("FALSE" %in% names(lag)) +}) + +test_that("a config the package writes is read correctly by a plain parser", { + # yaml::as.yaml quotes these keys on the way out, so the round trip does not + # depend on the reader knowing about any of this. + path <- tempfile(fileext = ".yaml") + save_config(list(covariates = list(derivoce = list( + list(type = "lag_covariate", vars = "SST", n = 2) + ))), path, header = FALSE) + + plain <- yaml::read_yaml(path) + + expect_equal(plain$covariates$derivoce[[1]]$n, 2) +}) + test_that("defaults target ECOMON", { config <- apply_config_defaults(list()) From e205788b9a6bdd586f2602e323888c9603bd1158 Mon Sep 17 00:00:00 2001 From: chross22 <52218551+chross22@users.noreply.github.com> Date: Mon, 10 Aug 2026 11:32:12 -0400 Subject: [PATCH 02/10] Test covariates by jackknife, and fit ensembles of model types Two things a run can now do that it could not, both off by default and both driven from the config. covariates.jackknife refits the model without each covariate in turn, over the same cross-validation folds as the main fit, and reports how much worse it ranks stations. Alongside it goes the other half of the classical jackknife, the model on that covariate alone, because the two answer different questions: leaving one out measures what only that covariate carries, fitting it alone measures what it carries at all. A covariate can score high on one and nothing on the other, and that is the case worth seeing - information that is real and duplicated is not the same thing as information that is not there. The reported p-value is a one-sided test of the per-fold differences with the variance correction of Nadeau and Bengio (2003). The correction is the load-bearing part. A plain paired t-test treats the folds as independent when any two training sets share most of their rows, so its variance estimate is badly optimistic and it finds significance everywhere (Dietterich 1998). There is no unbiased estimator of the variance of k-fold cross-validation (Bengio and Grandvalet 2004), so this inflates the naive variance by 1/k + 1/(k-1) instead, which roughly halves the statistic. For a GLM and a GAM the classical test is reported beside it rather than instead of it - drop-in-deviance for the first, mgcv's approximate term p-value for the second. That is why the fold test is the default criterion: it means the same thing for all four model types, and a forest has no likelihood to take a ratio of. Dropping is opt-in and defaults to off. A covariate that fails this test is one the other covariates already account for on these stations, which says as much about collinearity in this sample as about ecology - depth and surface temperature carry much of the same information on a shelf, and the test will call either one redundant depending on which the model reached for first. So the default is a table and a message saying what dropping would have removed. What it does remove goes through covariates.exclude, the mechanism that already existed, so a dropped covariate is still fetched for anything that needs it as an ingredient. model.ensemble fits several model types on the same data and combines them, which is BIOMOD_EnsembleModeling() from the pipeline this package replaces. Four combination rules, all computed and written; the config picks which one becomes the suitability layer. The ensemble's evaluation is its own rather than the average of its members'. Every member is fitted on the same folds from the same seed, so their held-out predictions line up by .row, and the combined out-of-fold predictions go through the same evaluation_table(), optimal_threshold() and bootstrap_evaluation() a single model does. Averaging member scores would have reported the ensemble as the average of its parts, which is not what an ensemble does: combining members that make different mistakes beats all of them and combining members that make the same mistakes does not, and only a cross-validated ensemble prediction tells those apart. The word ensemble was already taken. projection.uncertainty combines over resamples of the data within one algorithm; this combines over algorithms. Both can be on, and a projection then carries algorithm_sd beside suitability_sd rather than one standing in for the other. uncertainty_layers() becomes projection_layers() for the same reason. Both features parallelize over model fits through one shared helper, which forks - so Windows runs sequentially and says so rather than pretending. Co-Authored-By: Claude Opus 5 --- DESCRIPTION | 1 + NAMESPACE | 7 + R/config.R | 13 + R/ensemble.R | 831 ++++++++++++++++++++++++++++++++ R/jackknife.R | 702 +++++++++++++++++++++++++++ R/model_types.R | 6 + R/parallel.R | 106 ++++ R/pipeline.R | 130 ++++- R/plotting.R | 15 + R/project.R | 44 +- README.md | 224 ++++++++- inst/configs/cfin_gom.yaml | 51 ++ man/apply_jackknife_drop.Rd | 24 + man/build_ensemble.Rd | 23 + man/combine_members.Rd | 25 + man/corrected_paired_test.Rd | 37 ++ man/ensemble_cv_metrics.Rd | 31 ++ man/ensemble_importance.Rd | 29 ++ man/ensemble_member_spread.Rd | 35 ++ man/ensemble_oof_predictions.Rd | 26 + man/ensemble_rules.Rd | 45 ++ man/ensemble_settings.Rd | 77 +++ man/ensemble_spread_matrix.Rd | 22 + man/ensemble_weights.Rd | 29 ++ man/evaluation_value.Rd | 23 + man/fit_patch_ensemble.Rd | 89 ++++ man/fold_scores.Rd | 35 ++ man/gam_term_p_values.Rd | 24 + man/jackknife_covariates.Rd | 133 +++++ man/jackknife_defaults.Rd | 17 + man/jackknife_dropped.Rd | 41 ++ man/jackknife_settings.Rd | 65 +++ man/jackknife_type.Rd | 28 ++ man/member_config.Rd | 24 + man/member_metrics.Rd | 20 + man/member_score.Rd | 24 + man/parametric_tests.Rd | 33 ++ man/parse_jackknife.Rd | 18 + man/predict_grid_ensemble.Rd | 32 ++ man/print.taupatch_ensemble.Rd | 19 + man/projection_layers.Rd | 31 ++ man/report_jackknife.Rd | 22 + man/resolve_workers.Rd | 28 ++ man/run_taupatch.Rd | 23 +- man/subset_workflow.Rd | 26 + man/taupatch_lapply.Rd | 52 ++ man/uncertainty_layers.Rd | 21 - man/validate_jackknife.Rd | 18 + man/write_ensemble_outputs.Rd | 24 + man/write_jackknife.Rd | 22 + tests/testthat/test-ensemble.R | 394 +++++++++++++++ tests/testthat/test-jackknife.R | 323 +++++++++++++ 52 files changed, 4101 insertions(+), 41 deletions(-) create mode 100644 R/ensemble.R create mode 100644 R/jackknife.R create mode 100644 R/parallel.R create mode 100644 man/apply_jackknife_drop.Rd create mode 100644 man/build_ensemble.Rd create mode 100644 man/combine_members.Rd create mode 100644 man/corrected_paired_test.Rd create mode 100644 man/ensemble_cv_metrics.Rd create mode 100644 man/ensemble_importance.Rd create mode 100644 man/ensemble_member_spread.Rd create mode 100644 man/ensemble_oof_predictions.Rd create mode 100644 man/ensemble_rules.Rd create mode 100644 man/ensemble_settings.Rd create mode 100644 man/ensemble_spread_matrix.Rd create mode 100644 man/ensemble_weights.Rd create mode 100644 man/evaluation_value.Rd create mode 100644 man/fit_patch_ensemble.Rd create mode 100644 man/fold_scores.Rd create mode 100644 man/gam_term_p_values.Rd create mode 100644 man/jackknife_covariates.Rd create mode 100644 man/jackknife_defaults.Rd create mode 100644 man/jackknife_dropped.Rd create mode 100644 man/jackknife_settings.Rd create mode 100644 man/jackknife_type.Rd create mode 100644 man/member_config.Rd create mode 100644 man/member_metrics.Rd create mode 100644 man/member_score.Rd create mode 100644 man/parametric_tests.Rd create mode 100644 man/parse_jackknife.Rd create mode 100644 man/predict_grid_ensemble.Rd create mode 100644 man/print.taupatch_ensemble.Rd create mode 100644 man/projection_layers.Rd create mode 100644 man/report_jackknife.Rd create mode 100644 man/resolve_workers.Rd create mode 100644 man/subset_workflow.Rd create mode 100644 man/taupatch_lapply.Rd delete mode 100644 man/uncertainty_layers.Rd create mode 100644 man/validate_jackknife.Rd create mode 100644 man/write_ensemble_outputs.Rd create mode 100644 man/write_jackknife.Rd create mode 100644 tests/testthat/test-ensemble.R create mode 100644 tests/testthat/test-jackknife.R diff --git a/DESCRIPTION b/DESCRIPTION index 3337f95..816a443 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -23,6 +23,7 @@ Imports: datamatch, dplyr, ggplot2, + parallel, parsnip, recipes, rsample, diff --git a/NAMESPACE b/NAMESPACE index 38a8b7c..1d1b454 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -1,5 +1,6 @@ # Generated by roxygen2: do not edit by hand +S3method(print,taupatch_ensemble) export(abundance_units) export(add_derivoce_covariates) export(apply_prejoin_steps) @@ -21,7 +22,10 @@ export(derivoce_choices) export(derivoce_covariates) export(derivoce_required_inputs) export(derivoce_steps_for) +export(ensemble_rules) +export(ensemble_settings) export(fetch_covariates) +export(fit_patch_ensemble) export(fit_patch_model) export(format_zoop_data) export(gam_bases) @@ -32,6 +36,9 @@ export(generate_mock_covariates) export(generate_mock_zoop_data) export(glm_coefficients) export(is_raw_export) +export(jackknife_covariates) +export(jackknife_dropped) +export(jackknife_settings) export(label_patch) export(load_config) export(load_zoop_data) diff --git a/R/config.R b/R/config.R index 9569711..e501b99 100644 --- a/R/config.R +++ b/R/config.R @@ -25,6 +25,7 @@ load_config <- function(path) { validate_study_area(config) validate_covariates(config) validate_uncertainty(config) + validate_jackknife(config) config$species$resolved <- resolve_species(config) @@ -282,6 +283,18 @@ validate_threshold <- function(threshold, species_name) { #' @return `TRUE` invisibly; errors otherwise #' @keywords internal validate_model <- function(config) { + # An ensemble has no single type of its own. Its members are validated + # instead, each under the config it will actually be fitted with - so a + # per-member override that cannot work is an error now rather than partway + # through fitting the ensemble. + ensemble <- ensemble_settings(config) + if (!is.null(ensemble)) { + for (type in ensemble$types) { + validate_model(member_config(config, type, ensemble)) + } + return(invisible(TRUE)) + } + type <- resolve_model_type(config) entry <- model_types()[[type]] diff --git a/R/ensemble.R b/R/ensemble.R new file mode 100644 index 0000000..a64e080 --- /dev/null +++ b/R/ensemble.R @@ -0,0 +1,831 @@ +#' Multi-algorithm ensemble settings +#' +#' Which model types a run fits and how their projections are combined. This is +#' `BIOMOD_EnsembleModeling()` from the pipeline this package replaces: several +#' algorithms on the same data, filtered on how well they did, then averaged. +#' +#' ```yaml +#' model: +#' type: ensemble # or set the block below and leave type alone +#' ensemble: +#' types: [rf, brt, glm, gam] +#' rule: weighted_mean # or: mean, median, committee +#' weight_by: tss # or: roc_auc, pr_auc, equal +#' min_score: 0.4 # members scoring below this are excluded +#' workers: true # true = cores - 1; a count; false = sequential +#' settings: # per-member overrides of the model block +#' gam: +#' method: REML +#' brt: +#' learn_rate: 0.01 +#' ``` +#' +#' @section Two different things called an ensemble: +#' This one combines over **algorithms**, and its spread is disagreement about +#' the shape of the relationship — a forest and a logistic regression looking at +#' the same shelf and drawing different maps. +#' +#' `projection.uncertainty` (see [uncertainty_settings()]) combines over +#' **resamples of the data** within one algorithm, and its spread is how much +#' the fit moves when the stations move. +#' +#' They are independent and can both be on. When they are, each member carries +#' its own resample interval and the ensemble reports algorithm disagreement on +#' top of it, in separate columns — `algorithm_sd` against `suitability_sd`. +#' +#' @section Why filter members at all: +#' An ensemble that averages in a model which cannot separate the classes moves +#' the answer toward noise. `min_score` is biomod2's `metric.select.thresh` +#' under a plainer name, and 0.4 on TSS is a low bar deliberately: it is there +#' to catch a member that failed to fit anything, not to tune the ensemble by +#' selecting its best members on their own evaluation scores, which would be +#' selection on the same numbers used to report it. +#' +#' @param config a config list, as returned by `load_config()` +#' @return `NULL` when off, otherwise a list with `types`, `rule`, `weight_by`, +#' `min_score`, `workers`, and `settings` +#' @examples +#' config <- load_config( +#' system.file("configs/mock_test.yaml", package = "taupatch") +#' ) +#' ensemble_settings(config) # NULL: off by default +#' +#' config$model$type <- "ensemble" +#' ensemble_settings(config)$types +#' +#' config$model$ensemble <- list(types = c("rf", "glm"), rule = "median") +#' ensemble_settings(config) +#' @seealso [fit_patch_ensemble()], which runs it, and [uncertainty_settings()] +#' for the other kind of ensemble +#' @export +ensemble_settings <- function(config) { + spec <- config$model$ensemble + asked <- identical(config$model$type, "ensemble") + + if (is.null(spec) && !asked) return(NULL) + if (isFALSE(spec)) return(NULL) + if (is.null(spec) || isTRUE(spec)) spec <- list() + if (!is.list(spec)) { + stop("model.ensemble must be true, false, or a block of settings.", + call. = FALSE) + } + if (isFALSE(spec$enabled)) return(NULL) + + settings <- list( + types = as.character(spec$types %||% names(model_types())), + rule = spec$rule %||% "weighted_mean", + weight_by = spec$weight_by %||% "tss", + min_score = spec$min_score %||% 0.4, + workers = spec$workers, + settings = spec$settings %||% list() + ) + + unknown <- setdiff(settings$types, names(model_types())) + if (length(unknown) > 0) { + stop("Unknown model.ensemble.types: ", paste(unknown, collapse = ", "), + "\nAvailable: ", paste(names(model_types()), collapse = ", "), + call. = FALSE) + } + if (length(settings$types) < 2) { + stop("model.ensemble.types needs at least 2 model types to combine; got ", + length(settings$types), + ".\nFor a single type, set model.type to it and leave the ensemble off.", + call. = FALSE) + } + if (!(settings$rule %in% ensemble_rules())) { + stop("model.ensemble.rule must be one of: ", + paste(ensemble_rules(), collapse = ", "), ", got '", settings$rule, + "'.", call. = FALSE) + } + if (!(settings$weight_by %in% c("tss", "roc_auc", "pr_auc", "equal"))) { + stop("model.ensemble.weight_by must be one of: tss, roc_auc, pr_auc, ", + "equal; got '", settings$weight_by, "'.", call. = FALSE) + } + unnamed <- setdiff(names(settings$settings), settings$types) + if (length(unnamed) > 0) { + stop("model.ensemble.settings has overrides for types the ensemble does ", + "not fit: ", paste(unnamed, collapse = ", "), + "\nFitting: ", paste(settings$types, collapse = ", "), call. = FALSE) + } + settings +} + +#' Ways an ensemble can combine its members +#' +#' `mean` and `weighted_mean` average the probabilities, the second in +#' proportion to how well each member scored. `median` averages them robustly, +#' which is the one to reach for when a single member is capable of going badly +#' wrong somewhere on the grid — a boosted tree extrapolating, usually — since +#' a mean lets that member drag a cell and a median does not. +#' +#' `committee` is different in kind, and is biomod2's committee averaging: each +#' member binarises its own prediction at its own TSS-optimal cutoff, and the +#' cell gets the fraction of members that called it a patch. So it is already on +#' a 0-to-1 scale and reads directly as agreement — 0.75 means three of four +#' algorithms say patch — but it throws away how *confident* each member was. +#' +#' Every rule is computed and written on every run. `model.ensemble.rule` picks +#' which one is the `suitability` layer, and the others go beside it, because +#' the disagreement between rules is itself worth looking at and recomputing +#' them means refitting. +#' +#' @return character vector of rule names +#' @examples +#' ensemble_rules() +#' @references +#' Araújo MB, New M (2007). Ensemble forecasting of species distributions. +#' *Trends in Ecology & Evolution* **22**(1), 42-47. +#' \doi{10.1016/j.tree.2006.09.010} — why an ensemble of algorithms rather than +#' a chosen best one +#' +#' Marmion M, Parviainen M, Luoto M, Heikkinen RK, Thuiller W (2009). Evaluation +#' of consensus methods in predictive species distribution modelling. +#' *Diversity and Distributions* **15**(1), 59-69. +#' \doi{10.1111/j.1472-4642.2008.00491.x} — the rules compared against each +#' other +#' @export +ensemble_rules <- function() { + c("mean", "weighted_mean", "median", "committee") +} + +#' Fit an ensemble of model types on the same data +#' +#' Fits every type in `model.ensemble.types` on the same stations, the same +#' predictors and the same cross-validation folds, then combines them. The +#' result is a drop-in for a [fit_patch_model()] object: it carries an +#' `evaluation` table, a `classification_threshold`, an `importance` table and a +#' set of `predictors`, and [project_patch_model()] will project it. +#' +#' @section Why an ensemble at all: +#' The four types disagree in ways that are informative rather than incidental. +#' A random forest and a GLM that rank the same stations mean the relationships +#' are close to monotonic; a sharp disagreement means either a genuine +#' non-linearity or a forest fitting noise, and there is no way to tell which +#' from one model. Averaging them is the practical answer to not knowing which +#' is right, and the `algorithm_sd` surface a projection then carries is the map +#' of where that choice actually mattered. +#' +#' @section How the ensemble gets an honest evaluation: +#' Every member is fitted on the same folds, drawn from the same `model.seed`, +#' so the held-out predictions line up row for row. The ensemble's own +#' out-of-fold predictions are therefore built by combining members on the rows +#' none of them saw, and the reported evaluation, the TSS-optimal cutoff and its +#' bootstrap interval all come from those — the same functions, on the same +#' footing, as a single model's. +#' +#' This matters because the obvious alternative is wrong. Averaging the members' +#' evaluation scores would report the ensemble as the average of its parts, +#' which is not what an ensemble does: combining uncorrelated members usually +#' beats all of them, and combining correlated ones does not, and only a +#' cross-validated ensemble prediction can tell those apart. +#' +#' @section A member that fails: +#' A type whose package is not installed, or that will not fit these data, is +#' dropped with a warning rather than failing the run — an ensemble of three is +#' still an ensemble. Two members is the floor; below that the run stops, since +#' one algorithm averaged with nothing is a single model wearing a different +#' object. +#' +#' @param dat labeled modeling data from `label_patch()` with covariates attached +#' @param config a config list, as returned by `load_config()` +#' @param settings from [ensemble_settings()] +#' @return an object of class `taupatch_ensemble`: a list with `members` (the +#' per-type [fit_patch_model()] results), `summary` (one row per type, with +#' its score, weight and whether it qualified), `predictions` (combined +#' out-of-fold), `evaluation`, `classification_threshold`, +#' `classification_threshold_interval`, `importance` (weighted across +#' members), `metrics`, `rule`, `predictors`, `model_data`, `threshold`, and +#' `type`, which is `"ensemble"` +#' @examples +#' \dontrun{ +#' config <- load_config("my_run.yaml") +#' config$model$type <- "ensemble" +#' ensemble <- fit_patch_ensemble(dat, config) +#' ensemble$summary +#' ensemble$evaluation +#' } +#' @references +#' Araújo MB, New M (2007). Ensemble forecasting of species distributions. +#' *Trends in Ecology & Evolution* **22**(1), 42-47. +#' \doi{10.1016/j.tree.2006.09.010} +#' +#' Thuiller W, Lafourcade B, Engler R, Araújo MB (2009). BIOMOD - a platform for +#' ensemble forecasting of species distributions. *Ecography* **32**(3), +#' 369-373. \doi{10.1111/j.1600-0587.2008.05742.x} — what this replaces +#' @seealso [ensemble_settings()] for the config block, [ensemble_rules()] for +#' the combination rules, [fit_patch_model()] for a single member +#' @export +fit_patch_ensemble <- function(dat, config, settings = ensemble_settings(config)) { + settings <- settings %||% ensemble_settings(config) + if (is.null(settings)) { + stop("No ensemble is configured. Set model.type to 'ensemble', or pass ", + "settings built by ensemble_settings().", call. = FALSE) + } + + available <- Filter(function(type) { + installed <- tryCatch({ check_model_packages(type); TRUE }, + error = function(e) FALSE) + if (!installed) { + warning("Ensemble member '", type, "' needs the '", + model_types()[[type]]$package, "' package, which is not ", + "installed. Skipping it.", call. = FALSE) + } + installed + }, settings$types) + + if (length(available) < 2) { + stop("An ensemble needs at least 2 fittable member types; ", + length(available), " of ", length(settings$types), + " are available.\nInstall the missing packages, or name types that ", + "are installed in model.ensemble.types.", call. = FALSE) + } + + workers <- resolve_workers(settings$workers, length(available)) + message(" fitting ", length(available), " ensemble members (", + paste(available, collapse = ", "), ") across ", workers, + if (workers == 1) " worker" else " workers") + + fits <- taupatch_lapply(available, function(type) { + tryCatch(fit_patch_model(dat, member_config(config, type, settings)), + error = function(e) { + structure(list(type = type, message = conditionMessage(e)), + class = "taupatch_member_error") + }) + }, workers = workers, seed = config$model$seed) + names(fits) <- available + + broken <- vapply(fits, inherits, logical(1), "taupatch_member_error") + for (type in available[broken]) { + warning("Ensemble member '", type, "' failed to fit and was dropped: ", + fits[[type]]$message, call. = FALSE) + } + members <- fits[!broken] + if (length(members) < 2) { + stop("Only ", length(members), " ensemble member(s) fitted successfully, ", + "which is not an ensemble. See the warnings above for why the rest ", + "failed.", call. = FALSE) + } + + build_ensemble(members, config, settings) +} + +#' Assemble the fitted members into an ensemble +#' +#' Split from [fit_patch_ensemble()] so the scoring, weighting and combining can +#' be tested on members built any way at all, including hand-made ones. +#' +#' @param members a named list of [fit_patch_model()] results +#' @param config a config list, as returned by `load_config()` +#' @param settings from [ensemble_settings()] +#' @return a `taupatch_ensemble` +#' @keywords internal +build_ensemble <- function(members, config, settings) { + scores <- vapply(members, member_score, numeric(1), metric = settings$weight_by) + + qualifies <- !is.na(scores) & scores >= settings$min_score + if (!any(qualifies)) { + stop("No ensemble member reached model.ensemble.min_score (", + settings$min_score, ") on ", settings$weight_by, ". Best was ", + signif(max(scores, na.rm = TRUE), 3), + ".\nLower min_score, or fix why every algorithm is doing this badly.", + call. = FALSE) + } + if (sum(qualifies) < 2) { + # One qualifying member is a single model, and averaging it with nothing + # would report an ensemble that is not one. Kept as a warning rather than + # an error because the map it produces is still the right map. + warning("Only one ensemble member reached min_score (", settings$min_score, + "); the combined surface is that member alone.", call. = FALSE) + } + + weights <- ensemble_weights(scores, qualifies, settings$weight_by) + summary <- data.frame( + type = names(members), + label = vapply(names(members), function(t) model_types()[[t]]$label, + character(1)), + score = unname(scores), + metric = settings$weight_by, + cutoff = vapply(members, function(m) m$classification_threshold %||% NA_real_, + numeric(1)), + qualifies = unname(qualifies), + weight = unname(weights), + stringsAsFactors = FALSE + ) + summary <- summary[order(-summary$score), ] + rownames(summary) <- NULL + + keep <- names(members)[qualifies] + predictions <- ensemble_oof_predictions(members[keep], weights[keep], + settings$rule) + cutoff <- optimal_threshold(predictions) + bounds <- bootstrap_evaluation(predictions, cutoff, + times = bootstrap_times(config), + seed = config$model$seed) + + first <- members[[keep[1]]] + cv_metrics <- ensemble_cv_metrics(predictions) + structure(list( + members = members, + summary = summary, + weights = weights, + rule = settings$rule, + settings = settings, + predictions = predictions, + # Built from out-of-fold ensemble predictions by the same functions a single + # model uses, so an ensemble's evals.csv and a member's mean the same thing + # and can be read against each other. + evaluation = evaluation_table(predictions, cv_metrics, cutoff, bounds), + metrics = cv_metrics, + member_metrics = member_metrics(members), + classification_threshold = cutoff, + classification_threshold_interval = threshold_interval(bounds), + importance = ensemble_importance(members[keep], weights[keep]), + # The union rather than the first member's, so a grid cell is only predicted + # where every member has what it needs. + predictors = Reduce(union, lapply(members[keep], function(m) m$predictors)), + model_data = first$model_data, + threshold = first$threshold, + type = "ensemble", + # No single workflow to save. Named rather than absent so anything reaching + # for it gets a clear NULL instead of a partial match onto something else. + workflow = NULL + ), class = "taupatch_ensemble") +} + +#' One member's config +#' +#' The run's config with the member's type set, and any per-type overrides from +#' `model.ensemble.settings` merged into the model block. The uncertainty and +#' jackknife blocks are left alone, so a member inherits them exactly. +#' +#' @param config a config list, as returned by `load_config()` +#' @param type the member's model type +#' @param settings from [ensemble_settings()] +#' @return a config list for that member +#' @keywords internal +member_config <- function(config, type, settings) { + config$model$type <- type + # Tuning is per-type: a GLM has nothing to tune, and leaving a run-level + # `tune: true` in place would make the ensemble refuse to fit the one member + # that is the honest baseline. + if (length(model_types()[[type]]$tunable) == 0) config$model$tune <- FALSE + + override <- settings$settings[[type]] + if (!is.null(override)) config$model <- modifyList(config$model, override) + # A member never fits an ensemble of its own. + config$model$ensemble <- NULL + config +} + +#' One member's score, on the metric the weights use +#' +#' Read out of the member's own evaluation table rather than recomputed, so the +#' number that decides a member's weight is the number reported for it. The +#' threshold-dependent metrics are taken at the member's own TSS-optimal cutoff, +#' which is the only fair comparison — reading TSS at 0.5 would score every +#' member on a cutoff that suits none of them. +#' +#' @param member a [fit_patch_model()] result +#' @param metric `"tss"`, `"roc_auc"`, `"pr_auc"`, or `"equal"` +#' @return the score, or `NA_real_` +#' @keywords internal +member_score <- function(member, metric = "tss") { + if (identical(metric, "equal")) return(1) + + table <- member$evaluation + if (is.null(table)) return(NA_real_) + + rows <- if (metric %in% c("roc_auc", "pr_auc")) { + table[table$metric == metric & is.na(table$threshold), ] + } else { + at_optimal <- !is.na(table$threshold) & table$threshold != 0.5 + table[table$metric == metric & at_optimal, ] + } + if (nrow(rows) != 1) return(NA_real_) + rows$value +} + +#' Member weights from member scores +#' +#' Proportional to the score, over the qualifying members only, and summing to +#' one. A non-qualifying member's weight is zero rather than absent, so the +#' summary table shows what it would have been given. +#' +#' TSS runs from -1 to 1 and a negative score is a member predicting worse than +#' chance, so weights are floored at zero — a member cannot be given negative +#' influence, which would make the ensemble deliberately invert it. +#' +#' @param scores one score per member +#' @param qualifies which members cleared `min_score` +#' @param metric which metric the scores are on +#' @return a numeric vector of weights, summing to 1 +#' @keywords internal +ensemble_weights <- function(scores, qualifies, metric = "tss") { + weights <- rep(0, length(scores)) + names(weights) <- names(scores) + + usable <- qualifies & !is.na(scores) + if (!any(usable)) return(weights) + + if (identical(metric, "equal")) { + weights[usable] <- 1 / sum(usable) + return(weights) + } + + raw <- pmax(scores[usable], 0) + # Every qualifying member scored exactly zero: there is nothing to weight by, + # so weight them alike rather than dividing by zero. + weights[usable] <- if (sum(raw) > 0) raw / sum(raw) else 1 / sum(usable) + weights +} + +#' Combine the members' out-of-fold predictions +#' +#' Every member was cross-validated on the same folds from the same seed, so +#' their held-out predictions cover the same rows and can be combined row by +#' row. Matched on `.row` rather than on position, since `tune` returns folds in +#' its own order and two members need not agree on it. +#' +#' @param members the qualifying [fit_patch_model()] results +#' @param weights their weights +#' @param rule one of [ensemble_rules()] +#' @return a data frame of `.row`, `patch` and `.pred_patch`, in the shape the +#' evaluation functions expect +#' @keywords internal +ensemble_oof_predictions <- function(members, weights, rule = "weighted_mean") { + usable <- Filter(function(m) { + !is.null(m$predictions) && all(c(".row", ".pred_patch") %in% names(m$predictions)) + }, members) + if (length(usable) == 0) return(NULL) + + rows <- Reduce(intersect, lapply(usable, function(m) m$predictions$.row)) + if (length(rows) == 0) return(NULL) + + probabilities <- vapply(usable, function(m) { + m$predictions$.pred_patch[match(rows, m$predictions$.row)] + }, numeric(length(rows))) + probabilities <- matrix(probabilities, nrow = length(rows), + dimnames = list(NULL, names(usable))) + + cutoffs <- vapply(usable, function(m) m$classification_threshold %||% NA_real_, + numeric(1)) + combined <- combine_members(probabilities, weights[names(usable)], cutoffs) + + first <- usable[[1]]$predictions + at <- match(rows, first$.row) + out <- data.frame( + .row = rows, + patch = first$patch[at], + .pred_patch = combined[[rule]], + stringsAsFactors = FALSE + ) + # The fold each row was held out of, carried through so the ensemble can have + # a per-fold standard error like a single model does rather than only a + # pooled number. + if ("id" %in% names(first)) out$id <- first$id[at] + out +} + +#' Cross-validated metrics for the combined ensemble +#' +#' The ensemble's own held-out predictions, scored per fold and summarised in +#' the shape `tune::collect_metrics()` returns — so everything downstream that +#' reads a metrics table reads this one without knowing an ensemble produced it. +#' +#' Computed on the ensemble rather than averaged over members, because those are +#' different numbers and only the first is the ensemble's performance: combining +#' members that make different mistakes beats every one of them, and combining +#' members that make the same mistakes does not. +#' +#' The threshold-dependent rows are at the default 0.5 cutoff, which is what the +#' equivalent rows mean for a single model. [evaluation_table()] restates them at +#' the TSS-optimal cutoff alongside. +#' +#' @param predictions the combined out-of-fold predictions +#' @return a data frame of `.metric`, `.estimator`, `mean`, `n`, `std_err`, with +#' a `tss` row; `NULL` when the folds are not recoverable +#' @keywords internal +ensemble_cv_metrics <- function(predictions) { + if (is.null(predictions) || !("id" %in% names(predictions))) return(NULL) + + folds <- split(seq_len(nrow(predictions)), predictions$id) + per_fold <- function(fn) { + vapply(folds, function(rows) { + truth <- predictions$patch[rows] + probability <- predictions$.pred_patch[rows] + if (length(unique(truth)) < 2) return(NA_real_) + tryCatch(fn(truth, probability), error = function(e) NA_real_, + warning = function(w) NA_real_) + }, numeric(1)) + } + + hard <- function(truth, probability) { + factor(ifelse(probability >= 0.5, "patch", "non_patch"), + levels = levels(truth)) + } + + values <- list( + roc_auc = per_fold(yardstick::roc_auc_vec), + kap = per_fold(function(t, p) yardstick::kap_vec(t, hard(t, p))), + sens = per_fold(function(t, p) yardstick::sens_vec(t, hard(t, p))), + spec = per_fold(function(t, p) yardstick::spec_vec(t, hard(t, p))) + ) + + out <- do.call(rbind, lapply(names(values), function(metric) { + scores <- values[[metric]][is.finite(values[[metric]])] + data.frame( + .metric = metric, .estimator = "binary", + mean = if (length(scores) > 0) mean(scores) else NA_real_, + n = length(scores), + std_err = if (length(scores) > 1) { + stats::sd(scores) / sqrt(length(scores)) + } else { + NA_real_ + }, + stringsAsFactors = FALSE + ) + })) + add_tss(out) +} + +#' Apply every combination rule to a matrix of member predictions +#' +#' All four rules at once, because they cost nothing next to the predictions +#' they are computed from and a projection writes all of them. The spread +#' columns come out of the same matrix. +#' +#' @param probabilities rows by members +#' @param weights one per member, in the same column order +#' @param cutoffs each member's own TSS-optimal cutoff, for `committee` +#' @return a named list: one entry per rule, plus `algorithm_sd` and +#' `algorithm_range` +#' @keywords internal +combine_members <- function(probabilities, weights, cutoffs) { + weights <- weights[colnames(probabilities)] + weights[is.na(weights)] <- 0 + # An all-zero weight vector would make the weighted mean NaN everywhere. + if (sum(weights) == 0) weights <- rep(1 / ncol(probabilities), + ncol(probabilities)) + weights <- weights / sum(weights) + + # Each member binarises at its own cutoff, since a shared one would score + # members on a threshold suited to whichever happens to be best calibrated. + # A member with no cutoff falls back to 0.5. + cutoffs <- ifelse(is.na(cutoffs), 0.5, cutoffs) + called <- sweep(probabilities, 2, cutoffs, FUN = ">=") + + list( + mean = rowMeans(probabilities, na.rm = TRUE), + weighted_mean = as.numeric(probabilities %*% weights), + median = apply(probabilities, 1, stats::median, na.rm = TRUE), + committee = rowMeans(called, na.rm = TRUE), + algorithm_sd = if (ncol(probabilities) > 1) { + apply(probabilities, 1, stats::sd, na.rm = TRUE) + } else { + rep(0, nrow(probabilities)) + }, + algorithm_range = apply(probabilities, 1, max, na.rm = TRUE) - + apply(probabilities, 1, min, na.rm = TRUE) + ) +} + +#' The members' cross-validated metrics, stacked +#' +#' One `tune::collect_metrics()` table per member with a `type` column added, so +#' the run's `cv_metrics.csv` says which algorithm each row belongs to instead +#' of silently reporting one of them. +#' +#' @param members the [fit_patch_model()] results +#' @return a data frame +#' @keywords internal +member_metrics <- function(members) { + tables <- lapply(names(members), function(type) { + table <- members[[type]]$metrics + if (is.null(table) || nrow(table) == 0) return(NULL) + table$type <- type + table + }) + tables <- Filter(Negate(is.null), tables) + if (length(tables) == 0) return(NULL) + do.call(rbind, lapply(tables, as.data.frame)) +} + +#' Variable importance across an ensemble +#' +#' Each member's permutation importance, weighted by the member's weight and +#' summed. Permutation importance is the drop in ROC AUC when a predictor is +#' shuffled, which is the same quantity on the same scale for all four types — +#' that is exactly why the package computes it itself rather than asking each +#' engine — so averaging across them means something. +#' +#' The per-member columns are kept beside the ensemble figure. A predictor the +#' forest leans on and the GLM ignores is a fact about the shape of the +#' relationship, and the average is the one number that hides it. +#' +#' @param members the qualifying [fit_patch_model()] results +#' @param weights their weights +#' @return a tibble of `variable`, `importance`, and one column per member +#' @keywords internal +ensemble_importance <- function(members, weights) { + variables <- Reduce(union, lapply(members, function(m) m$importance$variable)) + if (length(variables) == 0) { + return(tibble::tibble(variable = character(), importance = numeric())) + } + + per_member <- vapply(members, function(m) { + m$importance$importance[match(variables, m$importance$variable)] + }, numeric(length(variables))) + per_member <- matrix(per_member, nrow = length(variables), + dimnames = list(NULL, names(members))) + + weights <- weights[colnames(per_member)] + weights[is.na(weights)] <- 0 + if (sum(weights) == 0) weights <- rep(1, length(weights)) + weights <- weights / sum(weights) + + filled <- per_member + filled[is.na(filled)] <- 0 + + out <- tibble::tibble(variable = variables, + importance = as.numeric(filled %*% weights)) + for (type in colnames(per_member)) out[[type]] <- per_member[, type] + dplyr::arrange(out, dplyr::desc(.data$importance)) +} + +#' Predict an ensemble across a covariate grid +#' +#' Every qualifying member predicts every cell, and the four rules plus the +#' spread come out of the same matrix. The `suitability` layer is whichever rule +#' `model.ensemble.rule` names; the rest go beside it, so a run can be read +#' against a different rule without refitting anything. +#' +#' `algorithm_sd` is the one to look at. It is disagreement between algorithms +#' on the same cell, which is a different question from `suitability_sd` — the +#' spread of one algorithm refitted on resampled stations — and a different one +#' again from `novelty`. A cell can be quiet on one and loud on another. +#' +#' @param ensemble a `taupatch_ensemble` from [fit_patch_ensemble()] +#' @param grid a covariate grid from `covariate_grid()` +#' @param uncertainty settings from [uncertainty_settings()], or `NULL` +#' @return a tibble of `lon`, `lat`, `suitability`, the other rules, the +#' algorithm spread, and the per-member surfaces; `NULL` if no complete rows +#' @keywords internal +predict_grid_ensemble <- function(ensemble, grid, uncertainty = NULL) { + complete <- grid[stats::complete.cases(grid[ensemble$predictors]), ] + if (nrow(complete) == 0) return(NULL) + + keep <- ensemble$summary$type[ensemble$summary$qualifies] + members <- ensemble$members[keep] + + predictions <- lapply(members, function(member) { + tryCatch( + stats::predict(member$workflow, new_data = complete, type = "prob")$.pred_patch, + error = function(e) NULL + ) + }) + usable <- !vapply(predictions, is.null, logical(1)) + if (sum(usable) == 0) return(NULL) + + probabilities <- do.call(cbind, predictions[usable]) + colnames(probabilities) <- names(members)[usable] + cutoffs <- vapply(members[usable], + function(m) m$classification_threshold %||% NA_real_, + numeric(1)) + combined <- combine_members(probabilities, ensemble$weights, cutoffs) + + out <- tibble::tibble( + lon = complete$lon, + lat = complete$lat, + suitability = combined[[ensemble$rule]], + algorithm_sd = combined$algorithm_sd, + algorithm_range = combined$algorithm_range, + n_algorithms = sum(usable) + ) + # Every rule the run did not pick, named for what it is rather than as + # `suitability_2`, so a GeoTIFF's layer names say which is which. + for (rule in setdiff(ensemble_rules(), ensemble$rule)) { + out[[paste0("suitability_", rule)]] <- combined[[rule]] + } + for (type in colnames(probabilities)) { + out[[paste0("member_", type)]] <- probabilities[, type] + } + + if (is.null(uncertainty)) return(out) + + # The resample interval, if one is wanted, is the weighted pooling of each + # member's own - so a cell's interval covers both refitting and the choice of + # algorithm rather than only whichever was asked for. + spread <- ensemble_member_spread(members[usable], ensemble$weights, complete, + uncertainty$level) + if (!is.null(spread)) out <- dplyr::bind_cols(out, spread) + + if (isTRUE(uncertainty$novelty)) { + out <- dplyr::bind_cols( + out, novelty_surface(complete, ensemble$model_data, ensemble$predictors) + ) + } + out +} + +#' Pool the members' resample intervals +#' +#' Each member carries its own resample ensemble when `projection.uncertainty` +#' is on. Rather than reporting one member's interval, or four of them, this +#' pools every member's every replicate into one set and takes the interval from +#' that — so the reported interval covers refit variability *and* algorithm +#' choice at once, which is what a reader of a single interval column assumes it +#' does. +#' +#' Members are represented in proportion to their weight by drawing that share +#' of the pooled columns, so a member with a tenth of the weight does not +#' contribute a quarter of the interval just for having been fitted. +#' +#' @param members the qualifying members that predicted successfully +#' @param weights their weights +#' @param newdata the cells to predict +#' @param level interval width +#' @return a data frame of `suitability_sd`, `suitability_lower`, +#' `suitability_upper` and `n_members`; `NULL` when no member has an ensemble +#' @keywords internal +ensemble_member_spread <- function(members, weights, newdata, level = 0.9) { + weights <- weights[names(members)] + weights[is.na(weights)] <- 0 + if (sum(weights) == 0) weights <- rep(1, length(members)) + weights <- weights / sum(weights) + + sizes <- vapply(members, function(m) length(m$ensemble %||% list()), integer(1)) + if (sum(sizes) == 0) return(NULL) + + # The pool is sized by the largest member's ensemble, so the shares are whole + # replicates rather than fractions of one. + budget <- max(sizes) + columns <- lapply(names(members), function(type) { + members_ensemble <- members[[type]]$ensemble + take <- min(length(members_ensemble), + max(1L, round(weights[[type]] * budget * length(members)))) + if (take == 0) return(NULL) + ensemble_spread_matrix(members_ensemble[seq_len(take)], newdata) + }) + columns <- Filter(Negate(is.null), columns) + if (length(columns) == 0) return(NULL) + + pooled <- do.call(cbind, columns) + if (is.null(pooled) || ncol(pooled) < 2) return(NULL) + + tail <- (1 - level) / 2 + data.frame( + suitability_sd = apply(pooled, 1, stats::sd, na.rm = TRUE), + suitability_lower = apply(pooled, 1, stats::quantile, probs = tail, + na.rm = TRUE, names = FALSE), + suitability_upper = apply(pooled, 1, stats::quantile, probs = 1 - tail, + na.rm = TRUE, names = FALSE), + n_members = ncol(pooled) + ) +} + +#' Predictions from every member of one resample ensemble, as a matrix +#' +#' The half of [ensemble_spread()] that produces the numbers, without reducing +#' them — so an ensemble of ensembles can pool the replicates before taking +#' quantiles rather than taking quantiles of quantiles. +#' +#' @param ensemble a list of fitted workflows +#' @param newdata the cells to predict +#' @return a matrix of cells by members, or `NULL` +#' @keywords internal +ensemble_spread_matrix <- function(ensemble, newdata) { + if (length(ensemble) == 0) return(NULL) + predictions <- lapply(ensemble, function(member) { + tryCatch( + stats::predict(member, new_data = newdata, type = "prob")$.pred_patch, + error = function(e) NULL + ) + }) + predictions <- Filter(Negate(is.null), predictions) + if (length(predictions) == 0) return(NULL) + do.call(cbind, predictions) +} + +#' Print an ensemble +#' +#' @param x a `taupatch_ensemble` +#' @param ... unused +#' @return `x`, invisibly +#' @export +print.taupatch_ensemble <- function(x, ...) { + cat("\n") + cat(" rule: ", x$rule, "\n", sep = "") + cat(" members (", sum(x$summary$qualifies), " of ", nrow(x$summary), + " qualifying):\n", sep = "") + print(x$summary[c("type", "score", "metric", "qualifies", "weight")], + row.names = FALSE) + + auc <- x$evaluation$value[x$evaluation$metric == "roc_auc" & + is.na(x$evaluation$threshold)] + if (length(auc) == 1) { + cat("\n ensemble ROC AUC (out of fold): ", signif(auc, 4), "\n", sep = "") + } + cat(" classification threshold: ", signif(x$classification_threshold, 4), + "\n", sep = "") + invisible(x) +} diff --git a/R/jackknife.R b/R/jackknife.R new file mode 100644 index 0000000..4d1e6bf --- /dev/null +++ b/R/jackknife.R @@ -0,0 +1,702 @@ +#' Covariate jackknife settings +#' +#' Whether a run tests its covariates before fitting, and what it does with the +#' answer. Off by default: it costs `2 * predictors + 1` cross-validations, which +#' is minutes on a station table and worth paying deliberately rather than on +#' every iteration. +#' +#' ```yaml +#' covariates: +#' jackknife: true # or the block below, for the non-defaults +#' jackknife: +#' metric: roc_auc # or: pr_auc +#' criterion: fold # or: parametric (glm and gam only) +#' alpha: 0.05 +#' adjust: holm # or: BH, bonferroni, none +#' drop: false # DEFAULT: report, never drop on its own +#' keep: [DEPTH, jday] # never dropped, whatever the test says +#' min_predictors: 2 # never drop below this many +#' workers: true # true = cores - 1; a count; false = sequential +#' ``` +#' +#' @section Dropping is opt-in, and that is deliberate: +#' `drop` defaults to `false`, so the default behaviour is a table and a message. +#' A covariate that fails this test is one the *other covariates already +#' account for* on these stations — which is a statement about collinearity in +#' this sample at least as much as about ecology. Bottom depth and sea surface +#' temperature carry much of the same information on a shelf; the test will +#' happily declare either one redundant depending on which the model reached for +#' first, and dropping it silently would make the map look better while removing +#' the variable a reader would have asked about. +#' +#' `keep` is the escape hatch for exactly that: a covariate that is in the model +#' because the study is about it stays in the model. +#' +#' @param config a config list, as returned by `load_config()` +#' @return `NULL` when off, otherwise a list with `metric`, `criterion`, +#' `alpha`, `adjust`, `drop`, `keep`, `min_predictors`, and `workers` +#' @examples +#' config <- load_config( +#' system.file("configs/mock_test.yaml", package = "taupatch") +#' ) +#' jackknife_settings(config) # NULL: off by default +#' +#' config$covariates$jackknife <- TRUE +#' jackknife_settings(config) # drop is FALSE +#' +#' config$covariates$jackknife <- list(drop = TRUE, keep = "jday") +#' jackknife_settings(config) +#' @seealso [jackknife_covariates()], which runs it +#' @export +jackknife_settings <- function(config) { + spec <- config$covariates$jackknife + if (is.null(spec) || isFALSE(spec)) return(NULL) + if (isTRUE(spec)) spec <- list() + + if (!is.list(spec)) { + stop("covariates.jackknife must be true, false, or a block of settings.", + call. = FALSE) + } + if (isFALSE(spec$enabled)) return(NULL) + parse_jackknife(spec) +} + +#' The jackknife settings a run would use with nothing configured +#' +#' [jackknife_covariates()] can be called on a config with no jackknife block at +#' all — testing covariates is a reasonable thing to do interactively without +#' editing a file for it — and this is what it uses then. +#' +#' @return the same shape [jackknife_settings()] returns +#' @keywords internal +jackknife_defaults <- function() parse_jackknife(list()) + +#' Validate and fill in one jackknife block +#' +#' @param spec the `covariates.jackknife` block, as a list +#' @return the settings list +#' @keywords internal +parse_jackknife <- function(spec) { + settings <- list( + metric = spec$metric %||% "roc_auc", + criterion = spec$criterion %||% "fold", + alpha = spec$alpha %||% 0.05, + adjust = spec$adjust %||% "holm", + # Never on by accident. A user who wants covariates removed from their own + # model says so in the config. + drop = isTRUE(spec$drop), + keep = as.character(spec$keep %||% character()), + min_predictors = as.integer(spec$min_predictors %||% 2L), + workers = spec$workers, + # Which model type does the testing. Normally the run's own, which is the + # only sensible default; naming one matters for an ensemble run, where + # there is no single type to inherit. + type = spec$type + ) + + if (!(settings$metric %in% c("roc_auc", "pr_auc"))) { + stop("covariates.jackknife.metric must be 'roc_auc' or 'pr_auc', got '", + settings$metric, "'.\nBoth are threshold-free, which is what lets ", + "them be compared fold by fold.", call. = FALSE) + } + if (!(settings$criterion %in% c("fold", "parametric"))) { + stop("covariates.jackknife.criterion must be 'fold' or 'parametric', got '", + settings$criterion, "'.", call. = FALSE) + } + if (!is.numeric(settings$alpha) || settings$alpha <= 0 || settings$alpha >= 1) { + stop("covariates.jackknife.alpha must be between 0 and 1, got ", + settings$alpha, ".", call. = FALSE) + } + if (!(settings$adjust %in% c(stats::p.adjust.methods))) { + stop("covariates.jackknife.adjust must be one of: ", + paste(stats::p.adjust.methods, collapse = ", "), ", got '", + settings$adjust, "'.", call. = FALSE) + } + if (is.na(settings$min_predictors) || settings$min_predictors < 1) { + stop("covariates.jackknife.min_predictors must be at least 1.", call. = FALSE) + } + if (!is.null(settings$type) && !(settings$type %in% names(model_types()))) { + stop("Unknown covariates.jackknife.type '", settings$type, "'.\nAvailable: ", + paste(names(model_types()), collapse = ", "), call. = FALSE) + } + settings +} + +#' Which model type does the jackknifing +#' +#' The run's own type, normally. An ensemble run has no single type, so it takes +#' the first member and says so — a covariate test has to be a test of +#' *something*, and silently picking one of four algorithms would leave a reader +#' of the table with no way to know which. +#' +#' `covariates.jackknife.type` overrides both. A GLM is the type to name there +#' if what is wanted is the classical answer, since it is the one whose test has +#' an exact form. +#' +#' @param config a config list, as returned by `load_config()` +#' @param settings from [jackknife_settings()] +#' @return a model type name +#' @keywords internal +jackknife_type <- function(config, settings) { + if (!is.null(settings$type)) return(settings$type) + + ensemble <- ensemble_settings(config) + if (is.null(ensemble)) return(resolve_model_type(config)) + + chosen <- ensemble$types[1] + message(" this run fits an ensemble, so the jackknife tests covariates ", + "against its first member ('", chosen, + "'); set covariates.jackknife.type to choose another") + chosen +} + +#' Validate the covariate jackknife block +#' +#' @param config a parsed config list +#' @return `TRUE` invisibly; errors otherwise +#' @keywords internal +validate_jackknife <- function(config) { + settings <- jackknife_settings(config) + if (is.null(settings)) return(invisible(TRUE)) + + # The fold test needs enough folds for a t with a usable number of degrees of + # freedom, and the Nadeau-Bengio correction is undefined at one fold. + folds <- config$model$cv_folds %||% 10 + if (folds < 3) { + stop("covariates.jackknife needs at least 3 model.cv_folds to have ", + "anything to test across; got ", folds, ".", call. = FALSE) + } + + # A parametric criterion on a forest would silently never fire, since there + # is no such test to report. Said at load rather than after the refits. + if (identical(settings$criterion, "parametric")) { + type <- settings$type %||% if (is.null(ensemble_settings(config))) { + resolve_model_type(config) + } else { + ensemble_settings(config)$types[1] + } + if (!(type %in% c("glm", "gam"))) { + stop("covariates.jackknife.criterion is 'parametric', but a ", + "likelihood-based test only exists for a 'glm' or 'gam'; this ", + "jackknife would use '", type, "'.\nUse criterion: fold, which is ", + "defined for every model type, or set covariates.jackknife.type.", + call. = FALSE) + } + } + invisible(TRUE) +} + +#' Test every covariate by leaving it out, in parallel +#' +#' The jackknife of Elith et al. (2011): refit the model without each covariate +#' in turn, and see how much worse it ranks stations. A covariate whose removal +#' costs nothing is one the others already account for. Alongside it goes the +#' other half of the classical jackknife — the model fitted on that covariate +#' *alone* — because the two answer different questions and the pair is what +#' makes the table readable: +#' +#' * **`score_without`** is low when the covariate carries something no other +#' covariate has. This is its *unique* contribution. +#' * **`score_only`** is high when the covariate carries a lot on its own, +#' whether or not anything else carries it too. +#' +#' A covariate can score high on one and nothing on the other, and that +#' combination is the informative one: high `score_only` with no unique +#' contribution means the information is real and duplicated, which is a very +#' different thing from a covariate that is simply uninformative. +#' +#' Every refit uses **the same cross-validation folds as the main model**, drawn +#' from `model.seed`, so the comparison is paired fold by fold and none of the +#' difference is the split moving underneath it. +#' +#' @section What "significant" means here: +#' The reported `p_value` is a one-sided test of whether leaving the covariate +#' out makes the model worse, computed from the per-fold differences with the +#' variance correction of Nadeau and Bengio (2003). +#' +#' The correction is the load-bearing part. A plain paired t-test across `k` +#' folds treats the folds as independent, and they are not — any two training +#' sets share most of their rows — so its variance estimate is badly optimistic +#' and it declares far more covariates significant than it should. There is no +#' unbiased estimator of the variance of k-fold cross-validation (Bengio and +#' Grandvalet 2004); the correction inflates the naive variance by +#' `1/k + 1/(k-1)` instead, which is the standard workable answer and roughly +#' halves the t statistic. +#' +#' `p_adjusted` then accounts for having asked the question once per covariate, +#' Holm by default. +#' +#' @section The parametric column: +#' For a GLM and a GAM there is an exact-ish test of the same hypothesis, and it +#' is reported beside the fold test rather than instead of it: +#' +#' * **`glm`** — the drop-in-deviance likelihood ratio test against the nested +#' model, `parametric_test` reading `LRT`. +#' * **`gam`** — `mgcv`'s approximate p-value for the term, `parametric_test` +#' reading `gam-approx`. It is approximate by construction: it does not +#' account for the smoothing parameters having been estimated from the same +#' data, so it runs anti-conservative (Wood 2017, section 6.12). +#' +#' A forest and a boosted tree have no likelihood, so these columns are `NA` +#' there. That is the whole reason the fold test is the default criterion — +#' it means the same thing for all four model types. +#' +#' @section Rows, not just columns: +#' Every model here is fitted on the rows that are complete across **all** +#' predictors, including the ones being left out. Letting a reduced model pick +#' up the rows its dropped covariate was missing would compare two models on +#' two different datasets, and the reduced one would sometimes win for that +#' reason alone. +#' +#' @param dat labeled modeling data from `label_patch()` with covariates attached +#' @param config a config list, as returned by `load_config()` +#' @param settings from [jackknife_settings()]; defaults are used when the +#' config has no jackknife block, so this can be called on any config +#' @return a data frame with one row per covariate, ordered by `contribution`, +#' carrying `variable`, `metric`, `score_full`, `score_without`, `score_only`, +#' `contribution`, `contribution_se`, `statistic`, `df`, `p_value`, +#' `p_adjusted`, `parametric_p`, `parametric_test`, `significant`, and +#' `n_folds`. The full model's score is on it as a `score_full` attribute. +#' @examples +#' \dontrun{ +#' config <- load_config("my_run.yaml") +#' dat <- label_patch(attach_covariates(load_zoop_data(config), +#' fetch_covariates(config), config), config) +#' jk <- jackknife_covariates(dat, config) +#' jk[c("variable", "contribution", "p_adjusted", "significant")] +#' } +#' @references +#' Elith J, Phillips SJ, Hastie T, Dudík M, Chee YE, Yates CJ (2011). A +#' statistical explanation of MaxEnt for ecologists. *Diversity and +#' Distributions* **17**(1), 43-57. +#' \doi{10.1111/j.1472-4642.2010.00725.x} — the leave-one-out / only-one pair +#' this reports +#' +#' Nadeau C, Bengio Y (2003). Inference for the generalization error. *Machine +#' Learning* **52**(3), 239-281. \doi{10.1023/A:1024068626366} — the variance +#' correction +#' +#' Bengio Y, Grandvalet Y (2004). No unbiased estimator of the variance of +#' k-fold cross-validation. *Journal of Machine Learning Research* **5**, +#' 1089-1105. — why a correction +#' is needed rather than a better estimator +#' +#' Dietterich TG (1998). Approximate statistical tests for comparing supervised +#' classification learning algorithms. *Neural Computation* **10**(7), +#' 1895-1923. \doi{10.1162/089976698300017197} — the inflated Type I error of +#' the uncorrected test +#' +#' Wood SN (2017). *Generalized Additive Models: An Introduction with R*, 2nd +#' edition. Chapman and Hall/CRC. \doi{10.1201/9781315370279} — the GAM term +#' p-values and their caveat +#' @seealso [jackknife_settings()] for the config block, [jackknife_dropped()] +#' for what `drop` would remove, [permutation_importance()] for the other +#' answer to "which covariate matters" +#' @export +jackknife_covariates <- function(dat, config, settings = NULL) { + settings <- settings %||% jackknife_settings(config) %||% jackknife_defaults() + + type <- jackknife_type(config, settings) + check_model_packages(type) + # Every fit below builds its specification from the config, so the type the + # jackknife tests has to be the type the config names. + config$model$type <- type + config$model$ensemble <- NULL + if (length(model_types()[[type]]$tunable) == 0) config$model$tune <- FALSE + + predictors <- predictor_names(dat, config) + if (length(predictors) < 2) { + stop("A jackknife needs at least 2 predictors to leave one out of; found ", + length(predictors), ".", call. = FALSE) + } + + model_data <- as.data.frame(dat[c(predictors, "patch")]) + model_data <- model_data[stats::complete.cases(model_data), , drop = FALSE] + if (nrow(model_data) == 0) { + stop("No rows are complete across every predictor, so there is nothing to ", + "jackknife.", call. = FALSE) + } + + # The same folds the main fit will use, so the two are talking about the same + # splits and the per-fold differences are paired. + set.seed(config$model$seed) + folds <- rsample::vfold_cv(model_data, v = config$model$cv_folds, + strata = "patch") + + # One task per model to fit: the full one, then each covariate left out, then + # each covariate on its own. Built as a flat list so the whole lot is spread + # across workers at once rather than in two waves. + tasks <- c( + list(list(kind = "full", variable = NA_character_, vars = predictors)), + lapply(predictors, function(v) { + list(kind = "without", variable = v, vars = setdiff(predictors, v)) + }), + lapply(predictors, function(v) { + list(kind = "only", variable = v, vars = v) + }) + ) + + workers <- resolve_workers(settings$workers, length(tasks)) + message(" jackknifing ", length(predictors), " covariates: ", length(tasks), + " cross-validations across ", workers, + if (workers == 1) " worker" else " workers") + + scores <- taupatch_lapply(tasks, function(task) { + fold_scores(task$vars, folds, config, type, settings$metric) + }, workers = workers, seed = config$model$seed) + + full <- scores[[1]] + without <- scores[seq_along(predictors) + 1L] + only <- scores[seq_along(predictors) + 1L + length(predictors)] + + parametric <- parametric_tests(model_data, predictors, config, type) + + rows <- lapply(seq_along(predictors), function(i) { + differences <- full - without[[i]] + test <- corrected_paired_test(differences) + data.frame( + variable = predictors[i], + metric = settings$metric, + score_full = mean(full, na.rm = TRUE), + score_without = mean(without[[i]], na.rm = TRUE), + score_only = mean(only[[i]], na.rm = TRUE), + contribution = test$estimate, + contribution_se = test$std_err, + statistic = test$statistic, + df = test$df, + p_value = test$p_value, + parametric_p = parametric$p_value[[i]], + parametric_test = parametric$test, + n_folds = test$n, + stringsAsFactors = FALSE + ) + }) + + out <- do.call(rbind, rows) + out$p_adjusted <- stats::p.adjust(out$p_value, method = settings$adjust) + criterion <- if (identical(settings$criterion, "parametric")) { + out$parametric_p + } else { + out$p_adjusted + } + # NA is not evidence of no effect. A covariate whose test could not be + # computed counts as contributing, so a failed test can never remove it. + out$significant <- is.na(criterion) | criterion < settings$alpha + + out <- out[c("variable", "metric", "score_full", "score_without", "score_only", + "contribution", "contribution_se", "statistic", "df", "p_value", + "p_adjusted", "parametric_p", "parametric_test", "significant", + "n_folds")] + out <- out[order(-out$contribution), ] + rownames(out) <- NULL + attr(out, "score_full") <- mean(full, na.rm = TRUE) + attr(out, "settings") <- settings + out +} + +#' One model's score on every fold +#' +#' Fitted and scored by hand rather than through `tune::fit_resamples()`, for +#' one reason: the per-fold numbers are the whole point here, and they have to +#' come from the *same* `rsample` folds for every covariate subset so the +#' differences pair up. Going through `tune` would mean re-deriving the folds +#' inside each call and getting the pairing only by luck. +#' +#' A fold that will not fit — a subset with one predictor that is constant on +#' that split, say — scores `NA` rather than failing the run, and the test +#' downstream drops it and reports the reduced `n_folds`. +#' +#' @param vars the predictors this model gets +#' @param folds the shared `rsample::vfold_cv()` object +#' @param config a config list, as returned by `load_config()` +#' @param type the model type being fitted +#' @param metric `"roc_auc"` or `"pr_auc"` +#' @return a numeric vector, one score per fold +#' @keywords internal +fold_scores <- function(vars, folds, config, type, metric = "roc_auc") { + score <- if (identical(metric, "pr_auc")) { + yardstick::pr_auc_vec + } else { + yardstick::roc_auc_vec + } + + vapply(folds$splits, function(split) { + train <- rsample::analysis(split)[c(vars, "patch")] + test <- rsample::assessment(split)[c(vars, "patch")] + + fitted <- tryCatch( + parsnip::fit(subset_workflow(train, vars, config, type), data = train), + error = function(e) NULL + ) + if (is.null(fitted)) return(NA_real_) + + tryCatch({ + probabilities <- stats::predict(fitted, new_data = test, type = "prob") + score(test$patch, probabilities$.pred_patch) + }, error = function(e) NA_real_, warning = function(w) NA_real_) + }, numeric(1)) +} + +#' The workflow for one covariate subset +#' +#' The same recipe and specification the real fit uses, restricted to a subset +#' of the predictors — so a jackknifed model differs from the full one in +#' exactly the covariate that was removed, and not in how it was preprocessed. +#' +#' @param train the training rows, carrying `vars` and `patch` +#' @param vars the predictors this model gets +#' @param config a config list, as returned by `load_config()` +#' @param type the model type being fitted +#' @return a `workflows::workflow()`, not yet fitted +#' @keywords internal +subset_workflow <- function(train, vars, config, type) { + wf <- workflows::add_recipe(workflows::workflow(), + build_recipe(train, config)) + formula <- model_formula(type, train, vars, config) + if (is.null(formula)) { + workflows::add_model(wf, build_model_spec(config)) + } else { + workflows::add_model(wf, build_model_spec(config), formula = formula) + } +} + +#' A paired test across folds, with the Nadeau-Bengio variance correction +#' +#' The naive paired t-test over `k` cross-validation folds pretends the folds +#' are independent. They share all but one fold's worth of training data, so its +#' variance estimate is far too small and it finds significance everywhere. This +#' inflates the variance by `1/k + 1/(k-1)` — the second term being the ratio of +#' 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. +#' +#' @param differences per-fold score of the full model minus the reduced one +#' @return a list of `estimate`, `std_err`, `statistic`, `df`, `p_value`, `n` +#' @references +#' Nadeau C, Bengio Y (2003). Inference for the generalization error. *Machine +#' Learning* **52**(3), 239-281. \doi{10.1023/A:1024068626366} +#' +#' Bouckaert RR, Frank E (2004). Evaluating the replicability of significance +#' tests for comparing learning algorithms. *Advances in Knowledge Discovery +#' and Data Mining*, 3-12. \doi{10.1007/978-3-540-24775-3_3} — the correction +#' applied to k-fold specifically +#' @keywords internal +corrected_paired_test <- function(differences) { + usable <- differences[is.finite(differences)] + k <- length(usable) + none <- list(estimate = NA_real_, std_err = NA_real_, statistic = NA_real_, + df = NA_real_, p_value = NA_real_, n = k) + if (k < 3) return(none) + + estimate <- mean(usable) + # 1/k is the naive paired variance; 1/(k-1) is the test-to-train size ratio + # that k-fold's overlapping training sets add. + 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. + 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)) + } + + statistic <- estimate / std_err + 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) +} + +#' The likelihood-based test, where the model type has one +#' +#' A GLM gets a drop-in-deviance likelihood ratio test against each nested +#' model, which is exact. A GAM gets `mgcv`'s approximate p-value for the term, +#' which is not — it conditions on smoothing parameters estimated from the same +#' data and so runs anti-conservative. A forest and a boosted tree get `NA`, +#' because there is no likelihood to take a ratio of. +#' +#' Fitted on the whole dataset rather than per fold: this is a test about the +#' model, not about its generalization, which is exactly what makes it a +#' different reading from the fold test beside it. +#' +#' @param model_data the complete-case modeling data +#' @param predictors predictor column names +#' @param config a config list, as returned by `load_config()` +#' @param type the model type being fitted +#' @return a list of `p_value` (one per predictor, in order) and `test` (a label) +#' @keywords internal +parametric_tests <- function(model_data, predictors, config, type) { + none <- list(p_value = rep(NA_real_, length(predictors)), + test = NA_character_) + if (!(type %in% c("glm", "gam"))) return(none) + + engine_fit <- function(vars) { + tryCatch( + workflows::extract_fit_engine( + parsnip::fit(subset_workflow(model_data[c(vars, "patch")], vars, config, type), + data = model_data) + ), + error = function(e) NULL + ) + } + + full <- engine_fit(predictors) + if (is.null(full)) return(none) + + if (identical(type, "gam")) { + return(list(p_value = gam_term_p_values(full, predictors), + test = "gam-approx")) + } + + # A GLM's nested comparison, one refit per covariate. Cheap next to the + # cross-validations that have already run, so it is not worth parallelizing + # separately. + p_values <- vapply(predictors, function(v) { + reduced <- engine_fit(setdiff(predictors, v)) + if (is.null(reduced)) return(NA_real_) + + statistic <- stats::deviance(reduced) - stats::deviance(full) + df <- reduced$df.residual - full$df.residual + if (!is.finite(statistic) || !is.finite(df) || df <= 0) return(NA_real_) + stats::pchisq(statistic, df = df, lower.tail = FALSE) + }, numeric(1)) + + list(p_value = unname(p_values), test = "LRT") +} + +#' Pull each predictor's approximate p-value out of a fitted GAM +#' +#' A predictor enters the formula either as a smooth or, when it had too few +#' distinct values for one, as a linear term — see [model_formula()] — so its +#' p-value is in the smooth table for some predictors and the parametric table +#' for others. Both are looked in, keyed on the term name rather than on +#' position. +#' +#' @param fit a fitted `mgcv::gam` +#' @param predictors predictor column names +#' @return a numeric vector of p-values, one per predictor +#' @keywords internal +gam_term_p_values <- function(fit, predictors) { + summary_gam <- tryCatch(summary(fit), error = function(e) NULL) + if (is.null(summary_gam)) return(rep(NA_real_, length(predictors))) + + smooth <- summary_gam$s.table + parametric <- summary_gam$p.table + + lookup <- c( + if (!is.null(smooth)) stats::setNames(smooth[, ncol(smooth)], rownames(smooth)), + if (!is.null(parametric)) { + stats::setNames(parametric[, ncol(parametric)], rownames(parametric)) + } + ) + if (length(lookup) == 0) return(rep(NA_real_, length(predictors))) + + vapply(predictors, function(v) { + # Single brackets, not double: a name that is not there gives NA here and + # an error there, and a predictor that entered linearly is genuinely absent + # from the smooth table. + hit <- lookup[paste0("s(", v, ")")] + if (is.na(hit)) hit <- lookup[v] + as.numeric(hit) + }, numeric(1), USE.NAMES = FALSE) +} + +#' Which covariates a jackknife would drop +#' +#' The `keep` list and `min_predictors` floor applied to the test result. Split +#' out from the run so a report-only jackknife can still say what dropping +#' *would* have removed, which is the number worth seeing before turning +#' `drop` on. +#' +#' When the floor binds, the covariates kept are the ones that contributed most, +#' so a run that would have dropped everything keeps the best of a bad set +#' rather than an arbitrary one. +#' +#' @param jk the result of [jackknife_covariates()] +#' @param settings from [jackknife_settings()] +#' @return a character vector of covariate names, possibly empty +#' @examples +#' jk <- data.frame( +#' variable = c("SST", "SSS", "CHL"), +#' contribution = c(0.08, 0.001, 0.0005), +#' significant = c(TRUE, FALSE, FALSE) +#' ) +#' settings <- list(keep = character(), min_predictors = 2) +#' jackknife_dropped(jk, settings) # only the weakest: the floor binds at 2 +#' +#' jackknife_dropped(jk, list(keep = "CHL", min_predictors = 1)) +#' @seealso [jackknife_covariates()] +#' @export +jackknife_dropped <- function(jk, settings) { + keep <- as.character(settings$keep %||% character()) + floor <- as.integer(settings$min_predictors %||% 2L) + + candidates <- setdiff(jk$variable[!jk$significant], keep) + if (length(candidates) == 0) return(character()) + + # Never below the floor. Ordered by contribution so the ones given back are + # the strongest of those that failed. + room <- nrow(jk) - floor + if (room <= 0) return(character()) + if (length(candidates) > room) { + ranked <- jk$variable[order(jk$contribution)] + candidates <- intersect(ranked, candidates)[seq_len(room)] + } + sort(candidates) +} + +#' Say what the jackknife found, in one place +#' +#' A table of fifteen columns is not something a run's log can print, and the +#' one thing a reader needs from it mid-run is which covariates failed and +#' whether anything is about to be removed on the strength of that. +#' +#' @param jk the result of [jackknife_covariates()] +#' @param settings from [jackknife_settings()] +#' @return `NULL`, invisibly +#' @keywords internal +report_jackknife <- function(jk, settings) { + failed <- jk$variable[!jk$significant] + if (length(failed) == 0) { + message(" every covariate contributes at alpha = ", settings$alpha, + " (", settings$adjust, "-adjusted)") + } else { + message(" no detectable contribution at alpha = ", settings$alpha, ": ", + paste(failed, collapse = ", ")) + } + + would_drop <- jackknife_dropped(jk, settings) + if (isTRUE(settings$drop)) { + if (length(would_drop) > 0) { + message(" covariates.jackknife.drop is on: removing ", + paste(would_drop, collapse = ", ")) + } + } else if (length(would_drop) > 0) { + message(" covariates.jackknife.drop is off, so nothing is removed. ", + "Setting it would drop: ", paste(would_drop, collapse = ", ")) + } + invisible(NULL) +} + +#' Remove the covariates a jackknife rejected +#' +#' Written into `covariates.exclude`, which is the mechanism that already +#' existed for keeping a fetched covariate out of the model, rather than a +#' second one beside it. So a dropped covariate is still downloaded and still +#' available to anything downstream that wants it — including a derived +#' covariate that needs it as an ingredient — it just stops being a predictor. +#' +#' @param config a config list, as returned by `load_config()` +#' @param dropped covariate names to exclude +#' @return `config`, with `covariates.exclude` extended +#' @keywords internal +apply_jackknife_drop <- function(config, dropped) { + if (length(dropped) == 0) return(config) + config$covariates$exclude <- union(config$covariates$exclude %||% character(), + dropped) + config +} diff --git a/R/model_types.R b/R/model_types.R index 403f343..e61b003 100644 --- a/R/model_types.R +++ b/R/model_types.R @@ -169,6 +169,12 @@ resolve_model_type <- function(config) { call. = FALSE) } + if (identical(type, "ensemble")) { + stop("model.type is 'ensemble', which is several types rather than one. ", + "Fit it with fit_patch_ensemble(), which run_taupatch() does ", + "automatically; fit_patch_model() takes a single type.\nThe members ", + "are set by model.ensemble.types.", call. = FALSE) + } if (!(type %in% names(catalog))) { stop("Unknown model.type '", type, "'.\nAvailable: ", paste(names(catalog), collapse = ", "), call. = FALSE) diff --git a/R/parallel.R b/R/parallel.R new file mode 100644 index 0000000..5cad5b9 --- /dev/null +++ b/R/parallel.R @@ -0,0 +1,106 @@ +#' Map a function over a list, in parallel where that is possible +#' +#' The one place in the package that spawns workers. Both callers — the +#' covariate jackknife and the multi-algorithm ensemble — are the same shape: +#' a few dozen independent model fits, each expensive enough that the cost of +#' handing it to another core disappears, and none of them talking to each +#' other. +#' +#' @section Forks, not sockets: +#' `parallel::mclapply()` forks, so each worker starts with the fitted +#' recipe, the folds and the station table already in memory and copy-on-write +#' keeps that free. A PSOCK cluster would have to serialize all of it to every +#' worker for every task, which on a station table is most of the time the +#' parallelism was meant to save. +#' +#' The cost is that forking does not exist on Windows, where this falls back to +#' running sequentially and says so rather than pretending. A jackknife is still +#' perfectly usable there — it is one model fit per covariate per fold, which is +#' minutes, not hours — it just does not get faster with more cores. +#' +#' @section Reproducibility: +#' Forked workers inherit the parent's RNG state, so without help every one of +#' them would draw the same random numbers — which for a random forest means the +#' members are correlated in a way nothing downstream can see. `L'Ecuyer-CMRG` +#' gives each worker an independent, reproducible substream, and the previous +#' RNG kind and seed are both restored on the way out so a run's own seed still +#' governs everything after this. +#' +#' @param x a list or vector to map over +#' @param fun the function to apply +#' @param workers how many workers; `1` runs sequentially +#' @param seed optional seed, so the mapping is reproducible +#' @return a list, as `lapply()` +#' @keywords internal +taupatch_lapply <- function(x, fun, workers = 1L, seed = NULL) { + workers <- min(as.integer(workers), length(x)) + + if (!is.null(seed)) { + if (exists(".Random.seed", envir = globalenv())) { + state <- get(".Random.seed", envir = globalenv()) + on.exit(assign(".Random.seed", state, envir = globalenv()), add = TRUE) + } + } + + if (is.na(workers) || workers <= 1L) { + if (!is.null(seed)) set.seed(seed) + return(lapply(x, fun)) + } + + previous_kind <- RNGkind("L'Ecuyer-CMRG")[1] + on.exit(RNGkind(previous_kind), add = TRUE) + if (!is.null(seed)) set.seed(seed) + + # mc.preschedule = FALSE deals out one task at a time. The tasks here are + # deliberately uneven - a model on one covariate costs a fraction of a model + # on all of them - and prescheduling would hand every cheap task to one core + # and leave it idle while another works through the expensive ones. + out <- parallel::mclapply(x, fun, mc.cores = workers, mc.preschedule = FALSE) + + failed <- vapply(out, inherits, logical(1), "try-error") + if (any(failed)) { + stop("A parallel worker failed: ", + conditionMessage(attr(out[[which(failed)[1]]], "condition")), + call. = FALSE) + } + out +} + +#' How many workers to run with +#' +#' `NULL` or `true` means "as many as this machine can spare", which is one +#' fewer than its physical cores — leaving one is what keeps the session it was +#' launched from responsive. `false` or `1` is sequential. +#' +#' Capped at `n`, since a task list of six cannot use twelve workers, and forced +#' to 1 on Windows, where [taupatch_lapply()] cannot fork. +#' +#' @param workers the configured value: `NULL`, a logical, or a count +#' @param n how many tasks there are to spread +#' @param quiet whether to suppress the Windows fallback message +#' @return an integer worker count, at least 1 +#' @keywords internal +resolve_workers <- function(workers = NULL, n = 1L, quiet = FALSE) { + requested <- if (is.null(workers) || isTRUE(workers)) { + available <- suppressWarnings(parallel::detectCores(logical = FALSE)) + if (is.na(available)) 1L else max(1L, available - 1L) + } else if (isFALSE(workers)) { + 1L + } else { + count <- suppressWarnings(as.integer(workers)) + if (is.na(count) || count < 1) { + stop("workers must be a positive count, true, or false; got '", + paste(workers, collapse = ", "), "'.", call. = FALSE) + } + count + } + + if (requested > 1L && identical(.Platform$OS.type, "windows")) { + if (!quiet) { + message(" parallel fits need forking, which Windows does not have; ", + "running sequentially") + } + return(1L) + } + max(1L, min(requested, as.integer(n))) +} diff --git a/R/pipeline.R b/R/pipeline.R index 3b370b6..0241f8e 100644 --- a/R/pipeline.R +++ b/R/pipeline.R @@ -4,13 +4,28 @@ #' and attach environmental covariates, label high-abundance patches against the #' species threshold, fit the model, and project monthly habitat suitability maps. #' +#' Two optional stages sit between labelling and fitting, both off by default +#' and both turned on from the config: +#' +#' * `covariates.jackknife` tests each covariate by leaving it out — see +#' [jackknife_settings()]. It runs before the fit so its answer can change +#' which covariates the model gets, and it only removes any if +#' `jackknife.drop` says so. +#' * `model.ensemble` fits several algorithms instead of one and combines them — +#' see [ensemble_settings()]. Everything after the fit works the same either +#' way, so a config that turns this on gets ensemble projections without +#' changing anything else. +#' #' @param config_path path to a config YAML file, or an already-loaded config list #' @param project whether to produce monthly projections after fitting #' @param keep_covariates the most covariate grid cells to return for mapping; #' `0` returns none. See [thin_covariates()] for what is kept and why. -#' @return a list with `config`, `data` (the labeled modeling data), `model` (the -#' `fit_patch_model()` result), `projections` (or `NULL` if skipped), -#' `covariate_means`, and `covariates` (a thinned grid, for mapping) +#' @return a list with `config` (as the run actually used it, so a jackknife that +#' dropped a covariate shows in `covariates.exclude`), `data` (the labeled +#' modeling data), `model` (a [fit_patch_model()] result, or a +#' [fit_patch_ensemble()] one), `projections` (or `NULL` if skipped), +#' `jackknife` (or `NULL` if not run), `covariate_means`, and `covariates` (a +#' thinned grid, for mapping) #' @examples #' \dontrun{ #' result <- run_taupatch(system.file("configs/mock_test.yaml", package = "taupatch")) @@ -74,11 +89,34 @@ run_taupatch <- function(config_path, project = TRUE, keep_covariates = 50000) { " (", sum(dat$patch == "patch"), " patch / ", sum(dat$patch == "non_patch"), " non-patch)") + # Before fitting, not after: the point of testing covariates is to decide + # which ones the model gets, and a test run against the final model would be + # describing a model that has already been built. + jackknife <- NULL + jackknife_config <- jackknife_settings(config) + if (!is.null(jackknife_config)) { + message("Jackknifing covariates...") + jackknife <- jackknife_covariates(dat, config, jackknife_config) + report_jackknife(jackknife, jackknife_config) + write_jackknife(jackknife, config) + if (isTRUE(jackknife_config$drop)) { + config <- apply_jackknife_drop(config, + jackknife_dropped(jackknife, jackknife_config)) + } + } + message("Fitting model...") - model <- fit_patch_model(dat, config) + ensemble <- ensemble_settings(config) + model <- if (is.null(ensemble)) { + fit_patch_model(dat, config) + } else { + fit_patch_ensemble(dat, config, ensemble) + } write_model_outputs(model, config) write_covariate_summary(covariate_means, config) - message(" ROC AUC: ", signif(model$metrics$mean[model$metrics$.metric == "roc_auc"], 4)) + # Read off the evaluation table rather than the metrics one, since that is + # the table both a single model and an ensemble fill in the same way. + message(" ROC AUC: ", signif(evaluation_value(model, "roc_auc"), 4)) projections <- NULL if (project) { @@ -89,7 +127,7 @@ run_taupatch <- function(config_path, project = TRUE, keep_covariates = 50000) { message("Output written to ", config$paths$output_dir) list(config = config, data = dat, model = model, projections = projections, - covariate_means = covariate_means, + covariate_means = covariate_means, jackknife = jackknife, # A thinned copy, for looking at rather than modelling. The full grid is # millions of points on a real fetch, and a map of every one of them # would be a map of a subsample anyway once it hit the screen. @@ -111,7 +149,10 @@ write_model_outputs <- function(model, config) { out <- config$paths$output_dir dir.create(out, recursive = TRUE, showWarnings = FALSE) - saveRDS(model$workflow, file.path(out, "model.rds")) + # An ensemble has no single workflow, so the whole object is what gets saved - + # it holds every member's, and a projection needs all of them. + saveRDS(model$workflow %||% model, file.path(out, "model.rds")) + if (inherits(model, "taupatch_ensemble")) write_ensemble_outputs(model, out) # evals.csv states the cutoff each metric belongs to, and reports the # threshold-dependent ones at both 0.5 and the TSS-optimal cutoff. The raw # resampling table is kept alongside for anything that wants the per-fold @@ -182,10 +223,78 @@ write_diagnostic_plots <- function(model, out) { plot_threshold_performance(predictions, file.path(diagnostics, "threshold_performance.png")) + # An ensemble has no coefficients or smooths of its own; its members do, and + # theirs are written into a directory each rather than averaged into + # something no model actually fitted. + if (inherits(model, "taupatch_ensemble")) { + for (type in names(model$members)) { + member_dir <- file.path(diagnostics, "members", type) + dir.create(member_dir, recursive = TRUE, showWarnings = FALSE) + write_effect_plots(model$members[[type]], member_dir) + } + return(invisible(NULL)) + } + write_effect_plots(model, diagnostics) invisible(NULL) } +#' Write an ensemble's own artifacts +#' +#' What a single model has no equivalent of: which algorithms were fitted, how +#' each scored, what weight it was given, and whether it qualified. This is the +#' first thing to read after an ensemble run — a table showing one member at +#' 0.9 weight and three near zero is a single model with extra steps, and only +#' this file says so. +#' +#' @param model a `taupatch_ensemble` from [fit_patch_ensemble()] +#' @param out the run's output directory +#' @return `NULL`, invisibly +#' @keywords internal +write_ensemble_outputs <- function(model, out) { + readr::write_csv(model$summary, file.path(out, "ensemble_members.csv")) + if (!is.null(model$member_metrics)) { + readr::write_csv(model$member_metrics, + file.path(out, "member_cv_metrics.csv")) + } + invisible(NULL) +} + +#' Write the covariate jackknife table +#' +#' Written whether or not anything was dropped, and written before the model is +#' fitted, so a run that turned `drop` on leaves a record of what it removed and +#' on what evidence. +#' +#' @param jk the result of [jackknife_covariates()] +#' @param config a config list, as returned by `load_config()` +#' @return `NULL`, invisibly +#' @keywords internal +write_jackknife <- function(jk, config) { + out <- config$paths$output_dir + dir.create(out, recursive = TRUE, showWarnings = FALSE) + readr::write_csv(jk, file.path(out, "covariate_jackknife.csv")) + invisible(NULL) +} + +#' One metric's value from a model's evaluation table +#' +#' The threshold-free metrics live on the rows with no cutoff. Reading them from +#' here rather than from the resampling table is what lets a single model and an +#' ensemble be asked the same question — the ensemble's resampling table is one +#' per member, and its own performance is not the average of those. +#' +#' @param model a [fit_patch_model()] or [fit_patch_ensemble()] result +#' @param metric a threshold-free metric name +#' @return the value, or `NA_real_` +#' @keywords internal +evaluation_value <- function(model, metric = "roc_auc") { + table <- model$evaluation + if (is.null(table)) return(NA_real_) + value <- table$value[table$metric == metric & is.na(table$threshold)] + if (length(value) == 1) value else NA_real_ +} + #' Write covariate summaries and month-by-year heatmaps #' #' @param covariate_means a data frame from `covariate_monthly_means()` @@ -235,6 +344,7 @@ pipeline_stages <- function() { "^Matching covariates to stations", "^Attaching climate indices", "^Labeling patches", + "^Jackknifing covariates", "^Fitting model", "^Projecting monthly suitability", "^Output written to" @@ -248,11 +358,15 @@ pipeline_stages <- function() { "Matching covariates to stations", "Attaching climate indices", "Labelling high-abundance patches", + "Testing covariates by jackknife", "Fitting and cross-validating the model", "Projecting monthly maps", "Writing output" ), - at = c(0.02, 0.06, 0.55, 0.60, 0.64, 0.72, 0.76, 0.78, 0.80, 0.92, 0.99), + # The jackknife is the widest band after the download when it runs at all: + # it is a full cross-validation per covariate, twice over. + at = c(0.02, 0.06, 0.55, 0.60, 0.64, 0.72, 0.76, 0.78, 0.79, 0.86, 0.94, + 0.99), stringsAsFactors = FALSE ) } diff --git a/R/plotting.R b/R/plotting.R index cfe6f58..0e3cdf9 100644 --- a/R/plotting.R +++ b/R/plotting.R @@ -62,6 +62,21 @@ plot_projection_uncertainty <- function(predicted, year, month, species, path) { ggplot2::labs(subtitle = "Spread across the ensemble") } + if ("algorithm_sd" %in% names(predicted)) { + # A third panel only when a multi-algorithm ensemble was fitted, and kept + # separate from the spread panel above it on purpose: that one is one + # algorithm refitted on resampled stations, this one is different + # algorithms on the same stations. A cell where the forest and the GLM + # disagree is not the same worry as a cell where the forest is unstable. + panels$algorithms <- ggplot2::ggplot( + predicted, ggplot2::aes(x = .data$lon, y = .data$lat, + fill = .data$algorithm_sd)) + + ggplot2::geom_raster() + + ggplot2::scale_fill_viridis_c(option = "cividis", na.value = "white", + name = "SD") + + ggplot2::labs(subtitle = "Disagreement between algorithms") + } + if ("novelty" %in% names(predicted)) { # Extrapolated cells are usually a small minority, and a scale stretched # over the whole range renders them invisible - which defeats the panel. diff --git a/R/project.R b/R/project.R index c8a59c9..0d29951 100644 --- a/R/project.R +++ b/R/project.R @@ -61,9 +61,22 @@ project_patch_model <- function(model, env_dat, config, bathy = NULL) { round(resolution * 111), " km), the grid the covariates were joined ", "onto") + if (inherits(model, "taupatch_ensemble")) { + qualifying <- sum(model$summary$qualifies) + message(" combining ", qualifying, " algorithms by ", model$rule, + "; the other rules and each member's own surface go beside it") + } + uncertainty <- uncertainty_settings(config) if (!is.null(uncertainty)) { - members <- length(model$ensemble) + # An algorithm ensemble carries one resample ensemble per member, which are + # pooled at prediction time; a single model carries the one. + members <- if (inherits(model, "taupatch_ensemble")) { + sum(vapply(model$members, function(m) length(m$ensemble %||% list()), + integer(1))) + } else { + length(model$ensemble) + } message(" uncertainty on: ", members, "-member ", uncertainty$method, " ensemble, ", round(100 * uncertainty$level), "% interval", if (isTRUE(uncertainty$novelty)) ", plus a novelty surface" else "") @@ -132,7 +145,7 @@ project_patch_model <- function(model, env_dat, config, bathy = NULL) { # about, and an unremarked map does not say so. Reported per month, # because it is the projected months at the edges of the record that # usually drift out. - extra <- uncertainty_layers(predicted) + extra <- projection_layers(predicted) if ("novelty" %in% extra) { note <- novelty_message(predicted$novelty, predicted$novel_variable) if (!is.null(note)) { @@ -240,6 +253,9 @@ project_patch_model <- function(model, env_dat, config, bathy = NULL) { #' spread and novelty columns; `NULL` if no complete rows #' @keywords internal predict_grid <- function(model, grid, uncertainty = NULL) { + if (inherits(model, "taupatch_ensemble")) { + return(predict_grid_ensemble(model, grid, uncertainty)) + } complete <- grid[stats::complete.cases(grid[model$predictors]), ] if (nrow(complete) == 0) return(NULL) @@ -262,18 +278,30 @@ predict_grid <- function(model, grid, uncertainty = NULL) { out } -#' Uncertainty layers present on a projection +#' Layers present on a projection beyond the suitability surface +#' +#' Which of the optional columns [predict_grid()] actually produced. Any of them +#' can come back empty — a model type that refuses to refit on a resample, an +#' ensemble member that will not predict this month's grid — and a map is +#' written either way rather than the run failing at the last step. +#' +#' Three different quantities can appear here and they are deliberately not +#' merged. `suitability_sd` is one algorithm refitted on resampled stations; +#' `algorithm_sd` is different algorithms on the same stations; `novelty` is how +#' far outside the training data the cell sits. A cell can be quiet on one and +#' loud on another, and that is the informative case rather than a contradiction. #' -#' Which of the optional columns [predict_grid()] actually produced. The -#' ensemble can come back empty — a model type that refuses to refit on a -#' resample, say — and a map is written either way rather than the run failing -#' at the last step. +#' Character columns are excluded by construction: these names become raster +#' layers, and `novel_variable` travels in the CSV instead. #' #' @param predicted a projection from [predict_grid()] #' @return character vector of column names beyond `suitability` #' @keywords internal -uncertainty_layers <- function(predicted) { +projection_layers <- function(predicted) { + rules <- paste0("suitability_", ensemble_rules()) intersect(c("suitability_sd", "suitability_lower", "suitability_upper", + "algorithm_sd", "algorithm_range", "n_algorithms", rules, + grep("^member_", names(predicted), value = TRUE), "novelty"), names(predicted)) } diff --git a/README.md b/README.md index 864941e..1838e6e 100644 --- a/README.md +++ b/README.md @@ -21,7 +21,9 @@ - [Transformations](#transformations) - [Preparing covariates before the join](#preparing-covariates-before-the-join) - [Combining products of different resolution](#combining-products-of-different-resolution) + - [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) - [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) @@ -638,6 +640,100 @@ native resolution. A run reports which covariates were upsampled, and records them on the result as an `upsampled` attribute. +### Testing which covariates earn their place + +`covariates.jackknife` refits the model without each covariate in turn and asks +how much worse it ranks stations. Off by default, because it is a full +cross-validation per covariate — twice over, since it also fits each covariate +alone — which is why it parallelizes: + +```yaml +covariates: + jackknife: true # or the block below, to change the defaults +``` + +```yaml +covariates: + jackknife: + metric: roc_auc # or: pr_auc + criterion: fold # or: parametric (glm and gam only) + alpha: 0.05 + adjust: holm # or: BH, bonferroni, none + drop: false # DEFAULT: report, never remove on its own + keep: [DEPTH, jday] # never dropped, whatever the test says + min_predictors: 2 + workers: true # true = cores - 1; a count; false = sequential +``` + +It writes `covariate_jackknife.csv`, one row per covariate: + +| variable | score_full | score_without | score_only | contribution | p_value | p_adjusted | significant | +|---|---|---|---|---|---|---|---| +| SST | 0.857 | 0.791 | 0.812 | 0.066 | 0.004 | 0.020 | TRUE | +| DEPTH | 0.857 | 0.828 | 0.774 | 0.029 | 0.031 | 0.124 | FALSE | +| jday | 0.857 | 0.855 | 0.611 | 0.002 | 0.402 | 0.402 | FALSE | + +The two halves answer different questions and the pair is what makes the table +readable. **`score_without`** is low when the covariate carries something no +other covariate has — its *unique* contribution. **`score_only`** is high when +it carries a lot on its own, whether or not anything else carries it too. A +covariate can score high on one and nothing on the other, and that combination +is the informative one: `DEPTH` above is worth as much alone as `SST` is, and +almost nothing on top of what the rest already say. + +**The significance test.** Every refit uses the same cross-validation folds as +the main model, so the comparison is paired fold by fold and none of the +difference is the split moving underneath it. The reported `p_value` is a +one-sided test of whether leaving the covariate out makes the model worse, +computed from the per-fold differences with the variance correction of Nadeau & +Bengio (2003). + +The correction is load-bearing. A plain paired t-test treats the folds as +independent, and they are not — any two training sets share most of their rows — +so its variance estimate is badly optimistic and it calls far too much +significant (Dietterich 1998). There is no unbiased estimator of the variance of +k-fold cross-validation (Bengio & Grandvalet 2004); this inflates the naive +variance by `1/k + 1/(k-1)` instead, which roughly halves the *t* statistic. +`p_adjusted` then accounts for having asked the question once per covariate. + +For a GLM and a GAM there is a classical test of the same hypothesis, and it is +reported *beside* the fold test rather than instead of it: a drop-in-deviance +likelihood ratio test for `glm`, and `mgcv`'s approximate term p-value for `gam` +(approximate because it conditions on smoothing parameters estimated from the +same data, so it runs anti-conservative — Wood 2017 §6.12). A forest and a +boosted tree have no likelihood, so those columns are `NA` there. That is why +`criterion: fold` is the default: it means the same thing for all four types. + +**Dropping is opt-in, and that is deliberate.** With `drop: false` the run +reports the table, says what dropping *would* have removed, and fits on +everything. A covariate that fails this test is one the *other covariates +already account for* on these stations, which is a statement about collinearity +in this sample at least as much as about ecology. Bottom depth and sea surface +temperature carry much of the same information on a shelf; the test will happily +call either one redundant depending on which the model reached for first. +Removing it silently would make the map look better while deleting the variable +a reader would have asked about. `keep` is the escape hatch for exactly that: a +covariate that is in the model because the study is about it stays in the model. + +When `drop: true`, the rejected covariates are written into +`covariates.exclude`, which is the mechanism that already existed for keeping a +fetched covariate out of the model. So a dropped covariate is still downloaded +and still available to anything that needs it as an ingredient — a gradient's +velocity components, say — it just stops being a predictor. The run's returned +`config` shows exactly what came out. + +Called directly, it needs no config block at all: + +```r +jk <- jackknife_covariates(dat, config) +jk[c("variable", "contribution", "p_adjusted", "significant")] +jackknife_dropped(jk, jackknife_settings(config)) # what drop would remove +``` + +Parallelism forks, which Windows does not have, so it runs sequentially there +and says so. It is one model fit per covariate per fold either way — minutes, +not hours — it just does not get faster with more cores. + ### Model type Four models, chosen with one word: @@ -718,6 +814,94 @@ absolute terms, but the ranking holds. Only `ranger` is needed for the default. `brt` needs `xgboost` and `gam` needs `mgcv`. Both are checked before fitting rather than at load. +### Fitting all of them at once + +Rather than picking one algorithm and hoping, fit several and combine them. This +is `BIOMOD_EnsembleModeling()` from the pipeline this package replaces: + +```yaml +model: + type: ensemble # all four, with the defaults below +``` + +```yaml +model: + type: ensemble + ensemble: + types: [rf, brt, glm, gam] + rule: weighted_mean # mean | weighted_mean | median | committee + weight_by: tss # tss | roc_auc | pr_auc | equal + min_score: 0.4 # members below this are excluded from the average + workers: true + settings: # per-member overrides of the model block + gam: + method: REML + brt: + learn_rate: 0.01 +``` + +Everything after the fit works the same either way, so a config that turns this +on gets ensemble projections without changing anything else. + +**Four ways to combine.** All of them are computed and written on every run; +`rule` picks which one becomes the `suitability` layer, and the others go beside +it — the disagreement between rules is itself worth looking at, and recomputing +them means refitting. + +| `rule` | What it does | +|---|---| +| `mean` | Plain average of the probabilities | +| `weighted_mean` | Average in proportion to how well each member scored | +| `median` | Robust average. The one to reach for when a single member is capable of going badly wrong somewhere on the grid — a boosted tree extrapolating, usually. A mean lets that member drag a cell; a median does not | +| `committee` | Each member binarises at *its own* TSS-optimal cutoff, and the cell gets the fraction of members calling it a patch. Reads directly as agreement — 0.75 means three of four algorithms say patch — but throws away how confident each was | + +**The ensemble gets its own honest evaluation.** Every member is fitted on the +same folds from the same seed, so their held-out predictions line up row for +row. The ensemble's out-of-fold predictions are built by combining members on +the rows none of them saw, and the reported `evals.csv`, the TSS-optimal cutoff +and its bootstrap interval all come from those — the same functions, on the same +footing, as a single model's. This matters because the obvious alternative is +wrong: averaging the members' *scores* would report the ensemble as the average +of its parts, and that is not what an ensemble does. Combining members that make +different mistakes beats all of them; combining members that make the same +mistakes does not. Only a cross-validated ensemble prediction tells those apart. + +`ensemble_members.csv` is the first thing to read afterwards: + +| type | label | score | metric | cutoff | qualifies | weight | +|---|---|---|---|---|---|---| +| brt | Boosted regression trees | 0.612 | tss | 0.089 | TRUE | 0.31 | +| rf | Random forest | 0.604 | tss | 0.051 | TRUE | 0.30 | +| gam | Generalized additive model | 0.541 | tss | 0.112 | TRUE | 0.27 | +| glm | Logistic regression | 0.238 | tss | 0.104 | TRUE | 0.12 | + +A table showing one member near 1.0 and the rest near zero is a single model +with extra steps, and only this file says so. `min_score` is biomod2's +`metric.select.thresh` under a plainer name, and 0.4 on TSS is a low bar +deliberately — it is there to catch a member that failed to fit anything, not to +tune the ensemble by selecting its best members on their own evaluation scores, +which would be selection on the numbers used to report it. A member whose +package is missing, or that will not fit, is dropped with a warning rather than +failing the run; two is the floor. + +Variable importance is weighted across members, with the per-member columns kept +beside it. A predictor the forest leans on and the GLM ignores is a fact about +the shape of the relationship, and the average is the one number that hides it. +Coefficients and smooths go to `diagnostics/members//`, since an ensemble +has none of its own and averaging them would describe a model nobody fitted. + +**Two different things are called an ensemble here**, and they are independent: + +| | Combines over | Its spread means | +|---|---|---| +| `model.ensemble` | **algorithms** | A forest and a logistic regression looking at the same shelf and drawing different maps → `algorithm_sd` | +| `projection.uncertainty` | **resamples of the data**, within one algorithm | How much the fit moves when the stations move → `suitability_sd` | + +Both can be on. When they are, each member carries its own resample interval, +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. + ### Training and projection windows These are separate. Fitting on a long history and projecting a shorter or later @@ -855,18 +1039,22 @@ refitting. Each run writes to `paths.output_dir`: ``` -model.rds fitted tidymodels workflow +model.rds fitted tidymodels workflow (or the whole ensemble object) evals.csv performance, stating the cutoff each metric belongs to cv_metrics.csv the raw per-fold resampling table var_importance.csv permutation variable importance var_importance.png threshold.yaml the abundance threshold used, and the probability cutoff with its interval +covariate_jackknife.csv with covariates.jackknife: each covariate's contribution and its p-value +ensemble_members.csv with model.ensemble: each algorithm's score, weight, and whether it qualified +member_cv_metrics.csv with model.ensemble: the per-fold resampling table, per algorithm diagnostics/roc_curve.png, pr_curve.png, calibration.png, threshold_performance.png diagnostics/cv_predictions.csv held-out predictions, for any metric not tabulated diagnostics/partial_effects.png what each predictor does to patch probability diagnostics/coefficients.png glm only: signed effects with intervals diagnostics/smooth_terms.csv gam only: effective degrees of freedom per smooth diagnostics/gam_smooths.png gam only, with fancygam: fitted smooths with error bands +diagnostics/members// with model.ensemble: the above, one directory per algorithm projections/suitability.csv every cell of every month: species, year, month, lon, lat, probability plus the interval and novelty columns, with projection.uncertainty projections/suitability.grd the same, as one raster with a layer per month (projection.write_grd) @@ -891,6 +1079,9 @@ R/prejoin.R prejoin_steps(), apply_prejoin_steps() R/derivoce.R derivoce_covariates(), add_derivoce_covariates() 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/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 R/evaluation_boot.R bootstrap_evaluation(), the interval on every metric @@ -1041,6 +1232,10 @@ Earth](https://www.naturalearthdata.com/), public domain, via `rnaturalearth`. ### Models +- Araújo MB, New M (2007). Ensemble forecasting of species distributions. + *Trends in Ecology & Evolution* **22**(1), 42–47. + [doi:10.1016/j.tree.2006.09.010](https://doi.org/10.1016/j.tree.2006.09.010) — + why an ensemble of algorithms rather than a chosen best one - Breiman L (2001). Random forests. *Machine Learning* **45**(1), 5–32. [doi:10.1023/A:1010933404324](https://doi.org/10.1023/A:1010933404324) — `rf`, and the origin of permutation importance @@ -1059,6 +1254,11 @@ Earth](https://www.naturalearthdata.com/), public domain, via `rnaturalearth`. - Hastie T, Tibshirani R (1986). Generalized additive models. *Statistical Science* **1**(3), 297–310. [doi:10.1214/ss/1177013604](https://doi.org/10.1214/ss/1177013604) — `gam` +- Marmion M, Parviainen M, Luoto M, Heikkinen RK, Thuiller W (2009). Evaluation + of consensus methods in predictive species distribution modelling. *Diversity + and Distributions* **15**(1), 59–69. + [doi:10.1111/j.1472-4642.2008.00491.x](https://doi.org/10.1111/j.1472-4642.2008.00491.x) + — the ensemble combination rules, compared against each other - Marra G, Wood SN (2011). Practical variable selection for generalized additive models. *Computational Statistics & Data Analysis* **55**(7), 2372–2387. [doi:10.1016/j.csda.2011.02.004](https://doi.org/10.1016/j.csda.2011.02.004) — @@ -1087,6 +1287,28 @@ Earth](https://www.naturalearthdata.com/), public domain, via `rnaturalearth`. *Journal of Applied Ecology* **43**(6), 1223–1232. [doi:10.1111/j.1365-2664.2006.01214.x](https://doi.org/10.1111/j.1365-2664.2006.01214.x) — `tss`, and the cutoff that maximises it +- Bengio Y, Grandvalet Y (2004). No unbiased estimator of the variance of k-fold + cross-validation. *Journal of Machine Learning Research* **5**, 1089–1105. + — why the jackknife's test + needs a variance correction rather than a better estimator +- Bouckaert RR, Frank E (2004). Evaluating the replicability of significance + tests for comparing learning algorithms. *Advances in Knowledge Discovery and + Data Mining*, 3–12. + [doi:10.1007/978-3-540-24775-3_3](https://doi.org/10.1007/978-3-540-24775-3_3) + — the correction applied to k-fold specifically +- Dietterich TG (1998). Approximate statistical tests for comparing supervised + classification learning algorithms. *Neural Computation* **10**(7), 1895–1923. + [doi:10.1162/089976698300017197](https://doi.org/10.1162/089976698300017197) — + the inflated Type I error of the uncorrected paired test +- Elith J, Phillips SJ, Hastie T, Dudík M, Chee YE, Yates CJ (2011). A + statistical explanation of MaxEnt for ecologists. *Diversity and + Distributions* **17**(1), 43–57. + [doi:10.1111/j.1472-4642.2010.00725.x](https://doi.org/10.1111/j.1472-4642.2010.00725.x) + — the leave-one-out / only-one jackknife pair +- 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` - 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/inst/configs/cfin_gom.yaml b/inst/configs/cfin_gom.yaml index f7a4c2d..909b023 100644 --- a/inst/configs/cfin_gom.yaml +++ b/inst/configs/cfin_gom.yaml @@ -164,14 +164,65 @@ covariates: # Days fetched in parallel. Modest values only - the bottleneck is the # Copernicus API, not local CPU. n_workers: 4 + # Test each covariate by refitting without it and seeing how much worse the + # model ranks stations. Off by default: it is a full cross-validation per + # covariate, twice over (once leaving it out, once with it alone), which is + # why `workers` is here. Run taupatch::jackknife_settings() for the block. + # + # Read the table it writes to covariate_jackknife.csv before touching `drop`. + # A covariate that fails this test is one the OTHER covariates already + # account for, which on a shelf - where depth and temperature carry much of + # the same information - says as much about collinearity in this sample as + # about ecology. `drop` is false by default for that reason, and `keep` is + # for the covariate that is in the model because the study is about it. + # jackknife: + # metric: roc_auc # or pr_auc; both are threshold-free + # criterion: fold # the corrected paired test across folds, for any + # # model type. `parametric` uses the GLM likelihood + # # ratio / GAM term p-value instead, which only exist + # # for those two types. + # alpha: 0.05 + # adjust: holm # for having asked the question once per covariate + # drop: false # DEFAULT. Report; never remove on its own. + # keep: [DEPTH, jday] # never dropped, whatever the test says + # min_predictors: 2 + # workers: true # true = cores - 1; a count; false = sequential model: # rf | brt | glm | gam - taupatch::model_types() describes each. + # `ensemble` fits several of them and combines their maps; see the block + # below and taupatch::ensemble_settings(). type: rf trees: 500 cv_folds: 10 tune: false seed: 42 + # Fit more than one algorithm and average them, rather than picking one and + # hoping. This is BIOMOD_EnsembleModeling() from the pipeline this package + # replaces. Set `type: ensemble` above to turn it on with these defaults, or + # write the block to change them. + # + # Note this is a different thing from projection.uncertainty below. That one + # resamples the data within one algorithm; this one combines across + # algorithms. Both can be on, and a projection then carries both spreads in + # separate columns - algorithm_sd against suitability_sd - because a cell can + # be quiet on one and loud on the other. + # ensemble: + # types: [rf, brt, glm, gam] + # rule: weighted_mean # mean | weighted_mean | median | committee. + # # Every rule is written; this picks which one is the + # # `suitability` layer. + # weight_by: tss # tss | roc_auc | pr_auc | equal, at each member's + # # own optimal cutoff + # min_score: 0.4 # members below this are excluded from the average. + # # A low bar on purpose: it is there to catch a + # # member that failed, not to select the best ones. + # workers: true + # settings: # per-member overrides of the model block above + # gam: + # method: REML + # brt: + # learn_rate: 0.01 projection: # Months to map. Omit years/months to reuse the training window; set them to diff --git a/man/apply_jackknife_drop.Rd b/man/apply_jackknife_drop.Rd new file mode 100644 index 0000000..3e6cb70 --- /dev/null +++ b/man/apply_jackknife_drop.Rd @@ -0,0 +1,24 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/jackknife.R +\name{apply_jackknife_drop} +\alias{apply_jackknife_drop} +\title{Remove the covariates a jackknife rejected} +\usage{ +apply_jackknife_drop(config, dropped) +} +\arguments{ +\item{config}{a config list, as returned by \code{load_config()}} + +\item{dropped}{covariate names to exclude} +} +\value{ +\code{config}, with \code{covariates.exclude} extended +} +\description{ +Written into \code{covariates.exclude}, which is the mechanism that already +existed for keeping a fetched covariate out of the model, rather than a +second one beside it. So a dropped covariate is still downloaded and still +available to anything downstream that wants it — including a derived +covariate that needs it as an ingredient — it just stops being a predictor. +} +\keyword{internal} diff --git a/man/build_ensemble.Rd b/man/build_ensemble.Rd new file mode 100644 index 0000000..5f88ae0 --- /dev/null +++ b/man/build_ensemble.Rd @@ -0,0 +1,23 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/ensemble.R +\name{build_ensemble} +\alias{build_ensemble} +\title{Assemble the fitted members into an ensemble} +\usage{ +build_ensemble(members, config, settings) +} +\arguments{ +\item{members}{a named list of \code{\link[=fit_patch_model]{fit_patch_model()}} results} + +\item{config}{a config list, as returned by \code{load_config()}} + +\item{settings}{from \code{\link[=ensemble_settings]{ensemble_settings()}}} +} +\value{ +a \code{taupatch_ensemble} +} +\description{ +Split from \code{\link[=fit_patch_ensemble]{fit_patch_ensemble()}} so the scoring, weighting and combining can +be tested on members built any way at all, including hand-made ones. +} +\keyword{internal} diff --git a/man/combine_members.Rd b/man/combine_members.Rd new file mode 100644 index 0000000..0bf92d2 --- /dev/null +++ b/man/combine_members.Rd @@ -0,0 +1,25 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/ensemble.R +\name{combine_members} +\alias{combine_members} +\title{Apply every combination rule to a matrix of member predictions} +\usage{ +combine_members(probabilities, weights, cutoffs) +} +\arguments{ +\item{probabilities}{rows by members} + +\item{weights}{one per member, in the same column order} + +\item{cutoffs}{each member's own TSS-optimal cutoff, for \code{committee}} +} +\value{ +a named list: one entry per rule, plus \code{algorithm_sd} and +\code{algorithm_range} +} +\description{ +All four rules at once, because they cost nothing next to the predictions +they are computed from and a projection writes all of them. The spread +columns come out of the same matrix. +} +\keyword{internal} diff --git a/man/corrected_paired_test.Rd b/man/corrected_paired_test.Rd new file mode 100644 index 0000000..dae11f4 --- /dev/null +++ b/man/corrected_paired_test.Rd @@ -0,0 +1,37 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/jackknife.R +\name{corrected_paired_test} +\alias{corrected_paired_test} +\title{A paired test across folds, with the Nadeau-Bengio variance correction} +\usage{ +corrected_paired_test(differences) +} +\arguments{ +\item{differences}{per-fold score of the full model minus the reduced one} +} +\value{ +a list of \code{estimate}, \code{std_err}, \code{statistic}, \code{df}, \code{p_value}, \code{n} +} +\description{ +The naive paired t-test over \code{k} cross-validation folds pretends the folds +are independent. They share all but one fold's worth of training data, so its +variance estimate is far too small and it finds significance everywhere. This +inflates the variance by \code{1/k + 1/(k-1)} — the second term being the ratio of +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. +} +\references{ +Nadeau C, Bengio Y (2003). Inference for the generalization error. \emph{Machine +Learning} \strong{52}(3), 239-281. \doi{10.1023/A:1024068626366} + +Bouckaert RR, Frank E (2004). Evaluating the replicability of significance +tests for comparing learning algorithms. \emph{Advances in Knowledge Discovery +and Data Mining}, 3-12. \doi{10.1007/978-3-540-24775-3_3} — the correction +applied to k-fold specifically +} +\keyword{internal} diff --git a/man/ensemble_cv_metrics.Rd b/man/ensemble_cv_metrics.Rd new file mode 100644 index 0000000..ee29b63 --- /dev/null +++ b/man/ensemble_cv_metrics.Rd @@ -0,0 +1,31 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/ensemble.R +\name{ensemble_cv_metrics} +\alias{ensemble_cv_metrics} +\title{Cross-validated metrics for the combined ensemble} +\usage{ +ensemble_cv_metrics(predictions) +} +\arguments{ +\item{predictions}{the combined out-of-fold predictions} +} +\value{ +a data frame of \code{.metric}, \code{.estimator}, \code{mean}, \code{n}, \code{std_err}, with +a \code{tss} row; \code{NULL} when the folds are not recoverable +} +\description{ +The ensemble's own held-out predictions, scored per fold and summarised in +the shape \code{tune::collect_metrics()} returns — so everything downstream that +reads a metrics table reads this one without knowing an ensemble produced it. +} +\details{ +Computed on the ensemble rather than averaged over members, because those are +different numbers and only the first is the ensemble's performance: combining +members that make different mistakes beats every one of them, and combining +members that make the same mistakes does not. + +The threshold-dependent rows are at the default 0.5 cutoff, which is what the +equivalent rows mean for a single model. \code{\link[=evaluation_table]{evaluation_table()}} restates them at +the TSS-optimal cutoff alongside. +} +\keyword{internal} diff --git a/man/ensemble_importance.Rd b/man/ensemble_importance.Rd new file mode 100644 index 0000000..84ff76e --- /dev/null +++ b/man/ensemble_importance.Rd @@ -0,0 +1,29 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/ensemble.R +\name{ensemble_importance} +\alias{ensemble_importance} +\title{Variable importance across an ensemble} +\usage{ +ensemble_importance(members, weights) +} +\arguments{ +\item{members}{the qualifying \code{\link[=fit_patch_model]{fit_patch_model()}} results} + +\item{weights}{their weights} +} +\value{ +a tibble of \code{variable}, \code{importance}, and one column per member +} +\description{ +Each member's permutation importance, weighted by the member's weight and +summed. Permutation importance is the drop in ROC AUC when a predictor is +shuffled, which is the same quantity on the same scale for all four types — +that is exactly why the package computes it itself rather than asking each +engine — so averaging across them means something. +} +\details{ +The per-member columns are kept beside the ensemble figure. A predictor the +forest leans on and the GLM ignores is a fact about the shape of the +relationship, and the average is the one number that hides it. +} +\keyword{internal} diff --git a/man/ensemble_member_spread.Rd b/man/ensemble_member_spread.Rd new file mode 100644 index 0000000..de18152 --- /dev/null +++ b/man/ensemble_member_spread.Rd @@ -0,0 +1,35 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/ensemble.R +\name{ensemble_member_spread} +\alias{ensemble_member_spread} +\title{Pool the members' resample intervals} +\usage{ +ensemble_member_spread(members, weights, newdata, level = 0.9) +} +\arguments{ +\item{members}{the qualifying members that predicted successfully} + +\item{weights}{their weights} + +\item{newdata}{the cells to predict} + +\item{level}{interval width} +} +\value{ +a data frame of \code{suitability_sd}, \code{suitability_lower}, +\code{suitability_upper} and \code{n_members}; \code{NULL} when no member has an ensemble +} +\description{ +Each member carries its own resample ensemble when \code{projection.uncertainty} +is on. Rather than reporting one member's interval, or four of them, this +pools every member's every replicate into one set and takes the interval from +that — so the reported interval covers refit variability \emph{and} algorithm +choice at once, which is what a reader of a single interval column assumes it +does. +} +\details{ +Members are represented in proportion to their weight by drawing that share +of the pooled columns, so a member with a tenth of the weight does not +contribute a quarter of the interval just for having been fitted. +} +\keyword{internal} diff --git a/man/ensemble_oof_predictions.Rd b/man/ensemble_oof_predictions.Rd new file mode 100644 index 0000000..0dba4a7 --- /dev/null +++ b/man/ensemble_oof_predictions.Rd @@ -0,0 +1,26 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/ensemble.R +\name{ensemble_oof_predictions} +\alias{ensemble_oof_predictions} +\title{Combine the members' out-of-fold predictions} +\usage{ +ensemble_oof_predictions(members, weights, rule = "weighted_mean") +} +\arguments{ +\item{members}{the qualifying \code{\link[=fit_patch_model]{fit_patch_model()}} results} + +\item{weights}{their weights} + +\item{rule}{one of \code{\link[=ensemble_rules]{ensemble_rules()}}} +} +\value{ +a data frame of \code{.row}, \code{patch} and \code{.pred_patch}, in the shape the +evaluation functions expect +} +\description{ +Every member was cross-validated on the same folds from the same seed, so +their held-out predictions cover the same rows and can be combined row by +row. Matched on \code{.row} rather than on position, since \code{tune} returns folds in +its own order and two members need not agree on it. +} +\keyword{internal} diff --git a/man/ensemble_rules.Rd b/man/ensemble_rules.Rd new file mode 100644 index 0000000..0116a3b --- /dev/null +++ b/man/ensemble_rules.Rd @@ -0,0 +1,45 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/ensemble.R +\name{ensemble_rules} +\alias{ensemble_rules} +\title{Ways an ensemble can combine its members} +\usage{ +ensemble_rules() +} +\value{ +character vector of rule names +} +\description{ +\code{mean} and \code{weighted_mean} average the probabilities, the second in +proportion to how well each member scored. \code{median} averages them robustly, +which is the one to reach for when a single member is capable of going badly +wrong somewhere on the grid — a boosted tree extrapolating, usually — since +a mean lets that member drag a cell and a median does not. +} +\details{ +\code{committee} is different in kind, and is biomod2's committee averaging: each +member binarises its own prediction at its own TSS-optimal cutoff, and the +cell gets the fraction of members that called it a patch. So it is already on +a 0-to-1 scale and reads directly as agreement — 0.75 means three of four +algorithms say patch — but it throws away how \emph{confident} each member was. + +Every rule is computed and written on every run. \code{model.ensemble.rule} picks +which one is the \code{suitability} layer, and the others go beside it, because +the disagreement between rules is itself worth looking at and recomputing +them means refitting. +} +\examples{ +ensemble_rules() +} +\references{ +Araújo MB, New M (2007). Ensemble forecasting of species distributions. +\emph{Trends in Ecology & Evolution} \strong{22}(1), 42-47. +\doi{10.1016/j.tree.2006.09.010} — why an ensemble of algorithms rather than +a chosen best one + +Marmion M, Parviainen M, Luoto M, Heikkinen RK, Thuiller W (2009). Evaluation +of consensus methods in predictive species distribution modelling. +\emph{Diversity and Distributions} \strong{15}(1), 59-69. +\doi{10.1111/j.1472-4642.2008.00491.x} — the rules compared against each +other +} diff --git a/man/ensemble_settings.Rd b/man/ensemble_settings.Rd new file mode 100644 index 0000000..5cd2c75 --- /dev/null +++ b/man/ensemble_settings.Rd @@ -0,0 +1,77 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/ensemble.R +\name{ensemble_settings} +\alias{ensemble_settings} +\title{Multi-algorithm ensemble settings} +\usage{ +ensemble_settings(config) +} +\arguments{ +\item{config}{a config list, as returned by \code{load_config()}} +} +\value{ +\code{NULL} when off, otherwise a list with \code{types}, \code{rule}, \code{weight_by}, +\code{min_score}, \code{workers}, and \code{settings} +} +\description{ +Which model types a run fits and how their projections are combined. This is +\code{BIOMOD_EnsembleModeling()} from the pipeline this package replaces: several +algorithms on the same data, filtered on how well they did, then averaged. +} +\details{ +\if{html}{\out{
}}\preformatted{model: + type: ensemble # or set the block below and leave type alone + ensemble: + types: [rf, brt, glm, gam] + rule: weighted_mean # or: mean, median, committee + weight_by: tss # or: roc_auc, pr_auc, equal + min_score: 0.4 # members scoring below this are excluded + workers: true # true = cores - 1; a count; false = sequential + settings: # per-member overrides of the model block + gam: + method: REML + brt: + learn_rate: 0.01 +}\if{html}{\out{
}} +} +\section{Two different things called an ensemble}{ + +This one combines over \strong{algorithms}, and its spread is disagreement about +the shape of the relationship — a forest and a logistic regression looking at +the same shelf and drawing different maps. + +\code{projection.uncertainty} (see \code{\link[=uncertainty_settings]{uncertainty_settings()}}) combines over +\strong{resamples of the data} within one algorithm, and its spread is how much +the fit moves when the stations move. + +They are independent and can both be on. When they are, each member carries +its own resample interval and the ensemble reports algorithm disagreement on +top of it, in separate columns — \code{algorithm_sd} against \code{suitability_sd}. +} + +\section{Why filter members at all}{ + +An ensemble that averages in a model which cannot separate the classes moves +the answer toward noise. \code{min_score} is biomod2's \code{metric.select.thresh} +under a plainer name, and 0.4 on TSS is a low bar deliberately: it is there +to catch a member that failed to fit anything, not to tune the ensemble by +selecting its best members on their own evaluation scores, which would be +selection on the same numbers used to report it. +} + +\examples{ +config <- load_config( + system.file("configs/mock_test.yaml", package = "taupatch") +) +ensemble_settings(config) # NULL: off by default + +config$model$type <- "ensemble" +ensemble_settings(config)$types + +config$model$ensemble <- list(types = c("rf", "glm"), rule = "median") +ensemble_settings(config) +} +\seealso{ +\code{\link[=fit_patch_ensemble]{fit_patch_ensemble()}}, which runs it, and \code{\link[=uncertainty_settings]{uncertainty_settings()}} +for the other kind of ensemble +} diff --git a/man/ensemble_spread_matrix.Rd b/man/ensemble_spread_matrix.Rd new file mode 100644 index 0000000..e5b0acd --- /dev/null +++ b/man/ensemble_spread_matrix.Rd @@ -0,0 +1,22 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/ensemble.R +\name{ensemble_spread_matrix} +\alias{ensemble_spread_matrix} +\title{Predictions from every member of one resample ensemble, as a matrix} +\usage{ +ensemble_spread_matrix(ensemble, newdata) +} +\arguments{ +\item{ensemble}{a list of fitted workflows} + +\item{newdata}{the cells to predict} +} +\value{ +a matrix of cells by members, or \code{NULL} +} +\description{ +The half of \code{\link[=ensemble_spread]{ensemble_spread()}} that produces the numbers, without reducing +them — so an ensemble of ensembles can pool the replicates before taking +quantiles rather than taking quantiles of quantiles. +} +\keyword{internal} diff --git a/man/ensemble_weights.Rd b/man/ensemble_weights.Rd new file mode 100644 index 0000000..667ac5a --- /dev/null +++ b/man/ensemble_weights.Rd @@ -0,0 +1,29 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/ensemble.R +\name{ensemble_weights} +\alias{ensemble_weights} +\title{Member weights from member scores} +\usage{ +ensemble_weights(scores, qualifies, metric = "tss") +} +\arguments{ +\item{scores}{one score per member} + +\item{qualifies}{which members cleared \code{min_score}} + +\item{metric}{which metric the scores are on} +} +\value{ +a numeric vector of weights, summing to 1 +} +\description{ +Proportional to the score, over the qualifying members only, and summing to +one. A non-qualifying member's weight is zero rather than absent, so the +summary table shows what it would have been given. +} +\details{ +TSS runs from -1 to 1 and a negative score is a member predicting worse than +chance, so weights are floored at zero — a member cannot be given negative +influence, which would make the ensemble deliberately invert it. +} +\keyword{internal} diff --git a/man/evaluation_value.Rd b/man/evaluation_value.Rd new file mode 100644 index 0000000..27a3064 --- /dev/null +++ b/man/evaluation_value.Rd @@ -0,0 +1,23 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/pipeline.R +\name{evaluation_value} +\alias{evaluation_value} +\title{One metric's value from a model's evaluation table} +\usage{ +evaluation_value(model, metric = "roc_auc") +} +\arguments{ +\item{model}{a \code{\link[=fit_patch_model]{fit_patch_model()}} or \code{\link[=fit_patch_ensemble]{fit_patch_ensemble()}} result} + +\item{metric}{a threshold-free metric name} +} +\value{ +the value, or \code{NA_real_} +} +\description{ +The threshold-free metrics live on the rows with no cutoff. Reading them from +here rather than from the resampling table is what lets a single model and an +ensemble be asked the same question — the ensemble's resampling table is one +per member, and its own performance is not the average of those. +} +\keyword{internal} diff --git a/man/fit_patch_ensemble.Rd b/man/fit_patch_ensemble.Rd new file mode 100644 index 0000000..c205380 --- /dev/null +++ b/man/fit_patch_ensemble.Rd @@ -0,0 +1,89 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/ensemble.R +\name{fit_patch_ensemble} +\alias{fit_patch_ensemble} +\title{Fit an ensemble of model types on the same data} +\usage{ +fit_patch_ensemble(dat, config, settings = ensemble_settings(config)) +} +\arguments{ +\item{dat}{labeled modeling data from \code{label_patch()} with covariates attached} + +\item{config}{a config list, as returned by \code{load_config()}} + +\item{settings}{from \code{\link[=ensemble_settings]{ensemble_settings()}}} +} +\value{ +an object of class \code{taupatch_ensemble}: a list with \code{members} (the +per-type \code{\link[=fit_patch_model]{fit_patch_model()}} results), \code{summary} (one row per type, with +its score, weight and whether it qualified), \code{predictions} (combined +out-of-fold), \code{evaluation}, \code{classification_threshold}, +\code{classification_threshold_interval}, \code{importance} (weighted across +members), \code{metrics}, \code{rule}, \code{predictors}, \code{model_data}, \code{threshold}, and +\code{type}, which is \code{"ensemble"} +} +\description{ +Fits every type in \code{model.ensemble.types} on the same stations, the same +predictors and the same cross-validation folds, then combines them. The +result is a drop-in for a \code{\link[=fit_patch_model]{fit_patch_model()}} object: it carries an +\code{evaluation} table, a \code{classification_threshold}, an \code{importance} table and a +set of \code{predictors}, and \code{\link[=project_patch_model]{project_patch_model()}} will project it. +} +\section{Why an ensemble at all}{ + +The four types disagree in ways that are informative rather than incidental. +A random forest and a GLM that rank the same stations mean the relationships +are close to monotonic; a sharp disagreement means either a genuine +non-linearity or a forest fitting noise, and there is no way to tell which +from one model. Averaging them is the practical answer to not knowing which +is right, and the \code{algorithm_sd} surface a projection then carries is the map +of where that choice actually mattered. +} + +\section{How the ensemble gets an honest evaluation}{ + +Every member is fitted on the same folds, drawn from the same \code{model.seed}, +so the held-out predictions line up row for row. The ensemble's own +out-of-fold predictions are therefore built by combining members on the rows +none of them saw, and the reported evaluation, the TSS-optimal cutoff and its +bootstrap interval all come from those — the same functions, on the same +footing, as a single model's. + +This matters because the obvious alternative is wrong. Averaging the members' +evaluation scores would report the ensemble as the average of its parts, +which is not what an ensemble does: combining uncorrelated members usually +beats all of them, and combining correlated ones does not, and only a +cross-validated ensemble prediction can tell those apart. +} + +\section{A member that fails}{ + +A type whose package is not installed, or that will not fit these data, is +dropped with a warning rather than failing the run — an ensemble of three is +still an ensemble. Two members is the floor; below that the run stops, since +one algorithm averaged with nothing is a single model wearing a different +object. +} + +\examples{ +\dontrun{ +config <- load_config("my_run.yaml") +config$model$type <- "ensemble" +ensemble <- fit_patch_ensemble(dat, config) +ensemble$summary +ensemble$evaluation +} +} +\references{ +Araújo MB, New M (2007). Ensemble forecasting of species distributions. +\emph{Trends in Ecology & Evolution} \strong{22}(1), 42-47. +\doi{10.1016/j.tree.2006.09.010} + +Thuiller W, Lafourcade B, Engler R, Araújo MB (2009). BIOMOD - a platform for +ensemble forecasting of species distributions. \emph{Ecography} \strong{32}(3), +369-373. \doi{10.1111/j.1600-0587.2008.05742.x} — what this replaces +} +\seealso{ +\code{\link[=ensemble_settings]{ensemble_settings()}} for the config block, \code{\link[=ensemble_rules]{ensemble_rules()}} for +the combination rules, \code{\link[=fit_patch_model]{fit_patch_model()}} for a single member +} diff --git a/man/fold_scores.Rd b/man/fold_scores.Rd new file mode 100644 index 0000000..7ee24bd --- /dev/null +++ b/man/fold_scores.Rd @@ -0,0 +1,35 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/jackknife.R +\name{fold_scores} +\alias{fold_scores} +\title{One model's score on every fold} +\usage{ +fold_scores(vars, folds, config, type, metric = "roc_auc") +} +\arguments{ +\item{vars}{the predictors this model gets} + +\item{folds}{the shared \code{rsample::vfold_cv()} object} + +\item{config}{a config list, as returned by \code{load_config()}} + +\item{type}{the model type being fitted} + +\item{metric}{\code{"roc_auc"} or \code{"pr_auc"}} +} +\value{ +a numeric vector, one score per fold +} +\description{ +Fitted and scored by hand rather than through \code{tune::fit_resamples()}, for +one reason: the per-fold numbers are the whole point here, and they have to +come from the \emph{same} \code{rsample} folds for every covariate subset so the +differences pair up. Going through \code{tune} would mean re-deriving the folds +inside each call and getting the pairing only by luck. +} +\details{ +A fold that will not fit — a subset with one predictor that is constant on +that split, say — scores \code{NA} rather than failing the run, and the test +downstream drops it and reports the reduced \code{n_folds}. +} +\keyword{internal} diff --git a/man/gam_term_p_values.Rd b/man/gam_term_p_values.Rd new file mode 100644 index 0000000..95eb332 --- /dev/null +++ b/man/gam_term_p_values.Rd @@ -0,0 +1,24 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/jackknife.R +\name{gam_term_p_values} +\alias{gam_term_p_values} +\title{Pull each predictor's approximate p-value out of a fitted GAM} +\usage{ +gam_term_p_values(fit, predictors) +} +\arguments{ +\item{fit}{a fitted \code{mgcv::gam}} + +\item{predictors}{predictor column names} +} +\value{ +a numeric vector of p-values, one per predictor +} +\description{ +A predictor enters the formula either as a smooth or, when it had too few +distinct values for one, as a linear term — see \code{\link[=model_formula]{model_formula()}} — so its +p-value is in the smooth table for some predictors and the parametric table +for others. Both are looked in, keyed on the term name rather than on +position. +} +\keyword{internal} diff --git a/man/jackknife_covariates.Rd b/man/jackknife_covariates.Rd new file mode 100644 index 0000000..0980ec4 --- /dev/null +++ b/man/jackknife_covariates.Rd @@ -0,0 +1,133 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/jackknife.R +\name{jackknife_covariates} +\alias{jackknife_covariates} +\title{Test every covariate by leaving it out, in parallel} +\usage{ +jackknife_covariates(dat, config, settings = NULL) +} +\arguments{ +\item{dat}{labeled modeling data from \code{label_patch()} with covariates attached} + +\item{config}{a config list, as returned by \code{load_config()}} + +\item{settings}{from \code{\link[=jackknife_settings]{jackknife_settings()}}; defaults are used when the +config has no jackknife block, so this can be called on any config} +} +\value{ +a data frame with one row per covariate, ordered by \code{contribution}, +carrying \code{variable}, \code{metric}, \code{score_full}, \code{score_without}, \code{score_only}, +\code{contribution}, \code{contribution_se}, \code{statistic}, \code{df}, \code{p_value}, +\code{p_adjusted}, \code{parametric_p}, \code{parametric_test}, \code{significant}, and +\code{n_folds}. The full model's score is on it as a \code{score_full} attribute. +} +\description{ +The jackknife of Elith et al. (2011): refit the model without each covariate +in turn, and see how much worse it ranks stations. A covariate whose removal +costs nothing is one the others already account for. Alongside it goes the +other half of the classical jackknife — the model fitted on that covariate +\emph{alone} — because the two answer different questions and the pair is what +makes the table readable: +} +\details{ +\itemize{ +\item \strong{\code{score_without}} is low when the covariate carries something no other +covariate has. This is its \emph{unique} contribution. +\item \strong{\code{score_only}} is high when the covariate carries a lot on its own, +whether or not anything else carries it too. +} + +A covariate can score high on one and nothing on the other, and that +combination is the informative one: high \code{score_only} with no unique +contribution means the information is real and duplicated, which is a very +different thing from a covariate that is simply uninformative. + +Every refit uses \strong{the same cross-validation folds as the main model}, drawn +from \code{model.seed}, so the comparison is paired fold by fold and none of the +difference is the split moving underneath it. +} +\section{What "significant" means here}{ + +The reported \code{p_value} is a one-sided test of whether leaving the covariate +out makes the model worse, computed from the per-fold differences with the +variance correction of Nadeau and Bengio (2003). + +The correction is the load-bearing part. A plain paired t-test across \code{k} +folds treats the folds as independent, and they are not — any two training +sets share most of their rows — so its variance estimate is badly optimistic +and it declares far more covariates significant than it should. There is no +unbiased estimator of the variance of k-fold cross-validation (Bengio and +Grandvalet 2004); the correction inflates the naive variance by +\code{1/k + 1/(k-1)} instead, which is the standard workable answer and roughly +halves the t statistic. + +\code{p_adjusted} then accounts for having asked the question once per covariate, +Holm by default. +} + +\section{The parametric column}{ + +For a GLM and a GAM there is an exact-ish test of the same hypothesis, and it +is reported beside the fold test rather than instead of it: +\itemize{ +\item \strong{\code{glm}} — the drop-in-deviance likelihood ratio test against the nested +model, \code{parametric_test} reading \code{LRT}. +\item \strong{\code{gam}} — \code{mgcv}'s approximate p-value for the term, \code{parametric_test} +reading \code{gam-approx}. It is approximate by construction: it does not +account for the smoothing parameters having been estimated from the same +data, so it runs anti-conservative (Wood 2017, section 6.12). +} + +A forest and a boosted tree have no likelihood, so these columns are \code{NA} +there. That is the whole reason the fold test is the default criterion — +it means the same thing for all four model types. +} + +\section{Rows, not just columns}{ + +Every model here is fitted on the rows that are complete across \strong{all} +predictors, including the ones being left out. Letting a reduced model pick +up the rows its dropped covariate was missing would compare two models on +two different datasets, and the reduced one would sometimes win for that +reason alone. +} + +\examples{ +\dontrun{ +config <- load_config("my_run.yaml") +dat <- label_patch(attach_covariates(load_zoop_data(config), + fetch_covariates(config), config), config) +jk <- jackknife_covariates(dat, config) +jk[c("variable", "contribution", "p_adjusted", "significant")] +} +} +\references{ +Elith J, Phillips SJ, Hastie T, Dudík M, Chee YE, Yates CJ (2011). A +statistical explanation of MaxEnt for ecologists. \emph{Diversity and +Distributions} \strong{17}(1), 43-57. +\doi{10.1111/j.1472-4642.2010.00725.x} — the leave-one-out / only-one pair +this reports + +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 + +Bengio Y, Grandvalet Y (2004). No unbiased estimator of the variance of +k-fold cross-validation. \emph{Journal of Machine Learning Research} \strong{5}, +1089-1105. \url{https://jmlr.org/papers/v5/grandvalet04a.html} — why a correction +is needed rather than a better estimator + +Dietterich TG (1998). Approximate statistical tests for comparing supervised +classification learning algorithms. \emph{Neural Computation} \strong{10}(7), +1895-1923. \doi{10.1162/089976698300017197} — the inflated Type I error of +the uncorrected test + +Wood SN (2017). \emph{Generalized Additive Models: An Introduction with R}, 2nd +edition. Chapman and Hall/CRC. \doi{10.1201/9781315370279} — the GAM term +p-values and their caveat +} +\seealso{ +\code{\link[=jackknife_settings]{jackknife_settings()}} for the config block, \code{\link[=jackknife_dropped]{jackknife_dropped()}} +for what \code{drop} would remove, \code{\link[=permutation_importance]{permutation_importance()}} for the other +answer to "which covariate matters" +} diff --git a/man/jackknife_defaults.Rd b/man/jackknife_defaults.Rd new file mode 100644 index 0000000..7fd1bf1 --- /dev/null +++ b/man/jackknife_defaults.Rd @@ -0,0 +1,17 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/jackknife.R +\name{jackknife_defaults} +\alias{jackknife_defaults} +\title{The jackknife settings a run would use with nothing configured} +\usage{ +jackknife_defaults() +} +\value{ +the same shape \code{\link[=jackknife_settings]{jackknife_settings()}} returns +} +\description{ +\code{\link[=jackknife_covariates]{jackknife_covariates()}} can be called on a config with no jackknife block at +all — testing covariates is a reasonable thing to do interactively without +editing a file for it — and this is what it uses then. +} +\keyword{internal} diff --git a/man/jackknife_dropped.Rd b/man/jackknife_dropped.Rd new file mode 100644 index 0000000..97c054c --- /dev/null +++ b/man/jackknife_dropped.Rd @@ -0,0 +1,41 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/jackknife.R +\name{jackknife_dropped} +\alias{jackknife_dropped} +\title{Which covariates a jackknife would drop} +\usage{ +jackknife_dropped(jk, settings) +} +\arguments{ +\item{jk}{the result of \code{\link[=jackknife_covariates]{jackknife_covariates()}}} + +\item{settings}{from \code{\link[=jackknife_settings]{jackknife_settings()}}} +} +\value{ +a character vector of covariate names, possibly empty +} +\description{ +The \code{keep} list and \code{min_predictors} floor applied to the test result. Split +out from the run so a report-only jackknife can still say what dropping +\emph{would} have removed, which is the number worth seeing before turning +\code{drop} on. +} +\details{ +When the floor binds, the covariates kept are the ones that contributed most, +so a run that would have dropped everything keeps the best of a bad set +rather than an arbitrary one. +} +\examples{ +jk <- data.frame( + variable = c("SST", "SSS", "CHL"), + contribution = c(0.08, 0.001, 0.0005), + significant = c(TRUE, FALSE, FALSE) +) +settings <- list(keep = character(), min_predictors = 2) +jackknife_dropped(jk, settings) # only the weakest: the floor binds at 2 + +jackknife_dropped(jk, list(keep = "CHL", min_predictors = 1)) +} +\seealso{ +\code{\link[=jackknife_covariates]{jackknife_covariates()}} +} diff --git a/man/jackknife_settings.Rd b/man/jackknife_settings.Rd new file mode 100644 index 0000000..956edc9 --- /dev/null +++ b/man/jackknife_settings.Rd @@ -0,0 +1,65 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/jackknife.R +\name{jackknife_settings} +\alias{jackknife_settings} +\title{Covariate jackknife settings} +\usage{ +jackknife_settings(config) +} +\arguments{ +\item{config}{a config list, as returned by \code{load_config()}} +} +\value{ +\code{NULL} when off, otherwise a list with \code{metric}, \code{criterion}, +\code{alpha}, \code{adjust}, \code{drop}, \code{keep}, \code{min_predictors}, and \code{workers} +} +\description{ +Whether a run tests its covariates before fitting, and what it does with the +answer. Off by default: it costs \code{2 * predictors + 1} cross-validations, which +is minutes on a station table and worth paying deliberately rather than on +every iteration. +} +\details{ +\if{html}{\out{
}}\preformatted{covariates: + jackknife: true # or the block below, for the non-defaults + jackknife: + metric: roc_auc # or: pr_auc + criterion: fold # or: parametric (glm and gam only) + alpha: 0.05 + adjust: holm # or: BH, bonferroni, none + drop: false # DEFAULT: report, never drop on its own + keep: [DEPTH, jday] # never dropped, whatever the test says + min_predictors: 2 # never drop below this many + workers: true # true = cores - 1; a count; false = sequential +}\if{html}{\out{
}} +} +\section{Dropping is opt-in, and that is deliberate}{ + +\code{drop} defaults to \code{false}, so the default behaviour is a table and a message. +A covariate that fails this test is one the \emph{other covariates already +account for} on these stations — which is a statement about collinearity in +this sample at least as much as about ecology. Bottom depth and sea surface +temperature carry much of the same information on a shelf; the test will +happily declare either one redundant depending on which the model reached for +first, and dropping it silently would make the map look better while removing +the variable a reader would have asked about. + +\code{keep} is the escape hatch for exactly that: a covariate that is in the model +because the study is about it stays in the model. +} + +\examples{ +config <- load_config( + system.file("configs/mock_test.yaml", package = "taupatch") +) +jackknife_settings(config) # NULL: off by default + +config$covariates$jackknife <- TRUE +jackknife_settings(config) # drop is FALSE + +config$covariates$jackknife <- list(drop = TRUE, keep = "jday") +jackknife_settings(config) +} +\seealso{ +\code{\link[=jackknife_covariates]{jackknife_covariates()}}, which runs it +} diff --git a/man/jackknife_type.Rd b/man/jackknife_type.Rd new file mode 100644 index 0000000..641044f --- /dev/null +++ b/man/jackknife_type.Rd @@ -0,0 +1,28 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/jackknife.R +\name{jackknife_type} +\alias{jackknife_type} +\title{Which model type does the jackknifing} +\usage{ +jackknife_type(config, settings) +} +\arguments{ +\item{config}{a config list, as returned by \code{load_config()}} + +\item{settings}{from \code{\link[=jackknife_settings]{jackknife_settings()}}} +} +\value{ +a model type name +} +\description{ +The run's own type, normally. An ensemble run has no single type, so it takes +the first member and says so — a covariate test has to be a test of +\emph{something}, and silently picking one of four algorithms would leave a reader +of the table with no way to know which. +} +\details{ +\code{covariates.jackknife.type} overrides both. A GLM is the type to name there +if what is wanted is the classical answer, since it is the one whose test has +an exact form. +} +\keyword{internal} diff --git a/man/member_config.Rd b/man/member_config.Rd new file mode 100644 index 0000000..b8a762e --- /dev/null +++ b/man/member_config.Rd @@ -0,0 +1,24 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/ensemble.R +\name{member_config} +\alias{member_config} +\title{One member's config} +\usage{ +member_config(config, type, settings) +} +\arguments{ +\item{config}{a config list, as returned by \code{load_config()}} + +\item{type}{the member's model type} + +\item{settings}{from \code{\link[=ensemble_settings]{ensemble_settings()}}} +} +\value{ +a config list for that member +} +\description{ +The run's config with the member's type set, and any per-type overrides from +\code{model.ensemble.settings} merged into the model block. The uncertainty and +jackknife blocks are left alone, so a member inherits them exactly. +} +\keyword{internal} diff --git a/man/member_metrics.Rd b/man/member_metrics.Rd new file mode 100644 index 0000000..c8acc8e --- /dev/null +++ b/man/member_metrics.Rd @@ -0,0 +1,20 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/ensemble.R +\name{member_metrics} +\alias{member_metrics} +\title{The members' cross-validated metrics, stacked} +\usage{ +member_metrics(members) +} +\arguments{ +\item{members}{the \code{\link[=fit_patch_model]{fit_patch_model()}} results} +} +\value{ +a data frame +} +\description{ +One \code{tune::collect_metrics()} table per member with a \code{type} column added, so +the run's \code{cv_metrics.csv} says which algorithm each row belongs to instead +of silently reporting one of them. +} +\keyword{internal} diff --git a/man/member_score.Rd b/man/member_score.Rd new file mode 100644 index 0000000..04a59c6 --- /dev/null +++ b/man/member_score.Rd @@ -0,0 +1,24 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/ensemble.R +\name{member_score} +\alias{member_score} +\title{One member's score, on the metric the weights use} +\usage{ +member_score(member, metric = "tss") +} +\arguments{ +\item{member}{a \code{\link[=fit_patch_model]{fit_patch_model()}} result} + +\item{metric}{\code{"tss"}, \code{"roc_auc"}, \code{"pr_auc"}, or \code{"equal"}} +} +\value{ +the score, or \code{NA_real_} +} +\description{ +Read out of the member's own evaluation table rather than recomputed, so the +number that decides a member's weight is the number reported for it. The +threshold-dependent metrics are taken at the member's own TSS-optimal cutoff, +which is the only fair comparison — reading TSS at 0.5 would score every +member on a cutoff that suits none of them. +} +\keyword{internal} diff --git a/man/parametric_tests.Rd b/man/parametric_tests.Rd new file mode 100644 index 0000000..0a08408 --- /dev/null +++ b/man/parametric_tests.Rd @@ -0,0 +1,33 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/jackknife.R +\name{parametric_tests} +\alias{parametric_tests} +\title{The likelihood-based test, where the model type has one} +\usage{ +parametric_tests(model_data, predictors, config, type) +} +\arguments{ +\item{model_data}{the complete-case modeling data} + +\item{predictors}{predictor column names} + +\item{config}{a config list, as returned by \code{load_config()}} + +\item{type}{the model type being fitted} +} +\value{ +a list of \code{p_value} (one per predictor, in order) and \code{test} (a label) +} +\description{ +A GLM gets a drop-in-deviance likelihood ratio test against each nested +model, which is exact. A GAM gets \code{mgcv}'s approximate p-value for the term, +which is not — it conditions on smoothing parameters estimated from the same +data and so runs anti-conservative. A forest and a boosted tree get \code{NA}, +because there is no likelihood to take a ratio of. +} +\details{ +Fitted on the whole dataset rather than per fold: this is a test about the +model, not about its generalization, which is exactly what makes it a +different reading from the fold test beside it. +} +\keyword{internal} diff --git a/man/parse_jackknife.Rd b/man/parse_jackknife.Rd new file mode 100644 index 0000000..831effb --- /dev/null +++ b/man/parse_jackknife.Rd @@ -0,0 +1,18 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/jackknife.R +\name{parse_jackknife} +\alias{parse_jackknife} +\title{Validate and fill in one jackknife block} +\usage{ +parse_jackknife(spec) +} +\arguments{ +\item{spec}{the \code{covariates.jackknife} block, as a list} +} +\value{ +the settings list +} +\description{ +Validate and fill in one jackknife block +} +\keyword{internal} diff --git a/man/predict_grid_ensemble.Rd b/man/predict_grid_ensemble.Rd new file mode 100644 index 0000000..1d11d47 --- /dev/null +++ b/man/predict_grid_ensemble.Rd @@ -0,0 +1,32 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/ensemble.R +\name{predict_grid_ensemble} +\alias{predict_grid_ensemble} +\title{Predict an ensemble across a covariate grid} +\usage{ +predict_grid_ensemble(ensemble, grid, uncertainty = NULL) +} +\arguments{ +\item{ensemble}{a \code{taupatch_ensemble} from \code{\link[=fit_patch_ensemble]{fit_patch_ensemble()}}} + +\item{grid}{a covariate grid from \code{covariate_grid()}} + +\item{uncertainty}{settings from \code{\link[=uncertainty_settings]{uncertainty_settings()}}, or \code{NULL}} +} +\value{ +a tibble of \code{lon}, \code{lat}, \code{suitability}, the other rules, the +algorithm spread, and the per-member surfaces; \code{NULL} if no complete rows +} +\description{ +Every qualifying member predicts every cell, and the four rules plus the +spread come out of the same matrix. The \code{suitability} layer is whichever rule +\code{model.ensemble.rule} names; the rest go beside it, so a run can be read +against a different rule without refitting anything. +} +\details{ +\code{algorithm_sd} is the one to look at. It is disagreement between algorithms +on the same cell, which is a different question from \code{suitability_sd} — the +spread of one algorithm refitted on resampled stations — and a different one +again from \code{novelty}. A cell can be quiet on one and loud on another. +} +\keyword{internal} diff --git a/man/print.taupatch_ensemble.Rd b/man/print.taupatch_ensemble.Rd new file mode 100644 index 0000000..de755c7 --- /dev/null +++ b/man/print.taupatch_ensemble.Rd @@ -0,0 +1,19 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/ensemble.R +\name{print.taupatch_ensemble} +\alias{print.taupatch_ensemble} +\title{Print an ensemble} +\usage{ +\method{print}{taupatch_ensemble}(x, ...) +} +\arguments{ +\item{x}{a \code{taupatch_ensemble}} + +\item{...}{unused} +} +\value{ +\code{x}, invisibly +} +\description{ +Print an ensemble +} diff --git a/man/projection_layers.Rd b/man/projection_layers.Rd new file mode 100644 index 0000000..f57fa36 --- /dev/null +++ b/man/projection_layers.Rd @@ -0,0 +1,31 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/project.R +\name{projection_layers} +\alias{projection_layers} +\title{Layers present on a projection beyond the suitability surface} +\usage{ +projection_layers(predicted) +} +\arguments{ +\item{predicted}{a projection from \code{\link[=predict_grid]{predict_grid()}}} +} +\value{ +character vector of column names beyond \code{suitability} +} +\description{ +Which of the optional columns \code{\link[=predict_grid]{predict_grid()}} actually produced. Any of them +can come back empty — a model type that refuses to refit on a resample, an +ensemble member that will not predict this month's grid — and a map is +written either way rather than the run failing at the last step. +} +\details{ +Three different quantities can appear here and they are deliberately not +merged. \code{suitability_sd} is one algorithm refitted on resampled stations; +\code{algorithm_sd} is different algorithms on the same stations; \code{novelty} is how +far outside the training data the cell sits. A cell can be quiet on one and +loud on another, and that is the informative case rather than a contradiction. + +Character columns are excluded by construction: these names become raster +layers, and \code{novel_variable} travels in the CSV instead. +} +\keyword{internal} diff --git a/man/report_jackknife.Rd b/man/report_jackknife.Rd new file mode 100644 index 0000000..5797424 --- /dev/null +++ b/man/report_jackknife.Rd @@ -0,0 +1,22 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/jackknife.R +\name{report_jackknife} +\alias{report_jackknife} +\title{Say what the jackknife found, in one place} +\usage{ +report_jackknife(jk, settings) +} +\arguments{ +\item{jk}{the result of \code{\link[=jackknife_covariates]{jackknife_covariates()}}} + +\item{settings}{from \code{\link[=jackknife_settings]{jackknife_settings()}}} +} +\value{ +\code{NULL}, invisibly +} +\description{ +A table of fifteen columns is not something a run's log can print, and the +one thing a reader needs from it mid-run is which covariates failed and +whether anything is about to be removed on the strength of that. +} +\keyword{internal} diff --git a/man/resolve_workers.Rd b/man/resolve_workers.Rd new file mode 100644 index 0000000..6c659bc --- /dev/null +++ b/man/resolve_workers.Rd @@ -0,0 +1,28 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/parallel.R +\name{resolve_workers} +\alias{resolve_workers} +\title{How many workers to run with} +\usage{ +resolve_workers(workers = NULL, n = 1L, quiet = FALSE) +} +\arguments{ +\item{workers}{the configured value: \code{NULL}, a logical, or a count} + +\item{n}{how many tasks there are to spread} + +\item{quiet}{whether to suppress the Windows fallback message} +} +\value{ +an integer worker count, at least 1 +} +\description{ +\code{NULL} or \code{true} means "as many as this machine can spare", which is one +fewer than its physical cores — leaving one is what keeps the session it was +launched from responsive. \code{false} or \code{1} is sequential. +} +\details{ +Capped at \code{n}, since a task list of six cannot use twelve workers, and forced +to 1 on Windows, where \code{\link[=taupatch_lapply]{taupatch_lapply()}} cannot fork. +} +\keyword{internal} diff --git a/man/run_taupatch.Rd b/man/run_taupatch.Rd index 20659f3..ae622e3 100644 --- a/man/run_taupatch.Rd +++ b/man/run_taupatch.Rd @@ -15,15 +15,32 @@ run_taupatch(config_path, project = TRUE, keep_covariates = 50000) \code{0} returns none. See \code{\link[=thin_covariates]{thin_covariates()}} for what is kept and why.} } \value{ -a list with \code{config}, \code{data} (the labeled modeling data), \code{model} (the -\code{fit_patch_model()} result), \code{projections} (or \code{NULL} if skipped), -\code{covariate_means}, and \code{covariates} (a thinned grid, for mapping) +a list with \code{config} (as the run actually used it, so a jackknife that +dropped a covariate shows in \code{covariates.exclude}), \code{data} (the labeled +modeling data), \code{model} (a \code{\link[=fit_patch_model]{fit_patch_model()}} result, or a +\code{\link[=fit_patch_ensemble]{fit_patch_ensemble()}} one), \code{projections} (or \code{NULL} if skipped), +\code{jackknife} (or \code{NULL} if not run), \code{covariate_means}, and \code{covariates} (a +thinned grid, for mapping) } \description{ Loads a config and runs every stage in order: read zooplankton stations, fetch and attach environmental covariates, label high-abundance patches against the species threshold, fit the model, and project monthly habitat suitability maps. } +\details{ +Two optional stages sit between labelling and fitting, both off by default +and both turned on from the config: +\itemize{ +\item \code{covariates.jackknife} tests each covariate by leaving it out — see +\code{\link[=jackknife_settings]{jackknife_settings()}}. It runs before the fit so its answer can change +which covariates the model gets, and it only removes any if +\code{jackknife.drop} says so. +\item \code{model.ensemble} fits several algorithms instead of one and combines them — +see \code{\link[=ensemble_settings]{ensemble_settings()}}. Everything after the fit works the same either +way, so a config that turns this on gets ensemble projections without +changing anything else. +} +} \examples{ \dontrun{ result <- run_taupatch(system.file("configs/mock_test.yaml", package = "taupatch")) diff --git a/man/subset_workflow.Rd b/man/subset_workflow.Rd new file mode 100644 index 0000000..cce70c6 --- /dev/null +++ b/man/subset_workflow.Rd @@ -0,0 +1,26 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/jackknife.R +\name{subset_workflow} +\alias{subset_workflow} +\title{The workflow for one covariate subset} +\usage{ +subset_workflow(train, vars, config, type) +} +\arguments{ +\item{train}{the training rows, carrying \code{vars} and \code{patch}} + +\item{vars}{the predictors this model gets} + +\item{config}{a config list, as returned by \code{load_config()}} + +\item{type}{the model type being fitted} +} +\value{ +a \code{workflows::workflow()}, not yet fitted +} +\description{ +The same recipe and specification the real fit uses, restricted to a subset +of the predictors — so a jackknifed model differs from the full one in +exactly the covariate that was removed, and not in how it was preprocessed. +} +\keyword{internal} diff --git a/man/taupatch_lapply.Rd b/man/taupatch_lapply.Rd new file mode 100644 index 0000000..238e690 --- /dev/null +++ b/man/taupatch_lapply.Rd @@ -0,0 +1,52 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/parallel.R +\name{taupatch_lapply} +\alias{taupatch_lapply} +\title{Map a function over a list, in parallel where that is possible} +\usage{ +taupatch_lapply(x, fun, workers = 1L, seed = NULL) +} +\arguments{ +\item{x}{a list or vector to map over} + +\item{fun}{the function to apply} + +\item{workers}{how many workers; \code{1} runs sequentially} + +\item{seed}{optional seed, so the mapping is reproducible} +} +\value{ +a list, as \code{lapply()} +} +\description{ +The one place in the package that spawns workers. Both callers — the +covariate jackknife and the multi-algorithm ensemble — are the same shape: +a few dozen independent model fits, each expensive enough that the cost of +handing it to another core disappears, and none of them talking to each +other. +} +\section{Forks, not sockets}{ + +\code{parallel::mclapply()} forks, so each worker starts with the fitted +recipe, the folds and the station table already in memory and copy-on-write +keeps that free. A PSOCK cluster would have to serialize all of it to every +worker for every task, which on a station table is most of the time the +parallelism was meant to save. + +The cost is that forking does not exist on Windows, where this falls back to +running sequentially and says so rather than pretending. A jackknife is still +perfectly usable there — it is one model fit per covariate per fold, which is +minutes, not hours — it just does not get faster with more cores. +} + +\section{Reproducibility}{ + +Forked workers inherit the parent's RNG state, so without help every one of +them would draw the same random numbers — which for a random forest means the +members are correlated in a way nothing downstream can see. \verb{L'Ecuyer-CMRG} +gives each worker an independent, reproducible substream, and the previous +RNG kind and seed are both restored on the way out so a run's own seed still +governs everything after this. +} + +\keyword{internal} diff --git a/man/uncertainty_layers.Rd b/man/uncertainty_layers.Rd deleted file mode 100644 index c40c02e..0000000 --- a/man/uncertainty_layers.Rd +++ /dev/null @@ -1,21 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/project.R -\name{uncertainty_layers} -\alias{uncertainty_layers} -\title{Uncertainty layers present on a projection} -\usage{ -uncertainty_layers(predicted) -} -\arguments{ -\item{predicted}{a projection from \code{\link[=predict_grid]{predict_grid()}}} -} -\value{ -character vector of column names beyond \code{suitability} -} -\description{ -Which of the optional columns \code{\link[=predict_grid]{predict_grid()}} actually produced. The -ensemble can come back empty — a model type that refuses to refit on a -resample, say — and a map is written either way rather than the run failing -at the last step. -} -\keyword{internal} diff --git a/man/validate_jackknife.Rd b/man/validate_jackknife.Rd new file mode 100644 index 0000000..9ab9c68 --- /dev/null +++ b/man/validate_jackknife.Rd @@ -0,0 +1,18 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/jackknife.R +\name{validate_jackknife} +\alias{validate_jackknife} +\title{Validate the covariate jackknife block} +\usage{ +validate_jackknife(config) +} +\arguments{ +\item{config}{a parsed config list} +} +\value{ +\code{TRUE} invisibly; errors otherwise +} +\description{ +Validate the covariate jackknife block +} +\keyword{internal} diff --git a/man/write_ensemble_outputs.Rd b/man/write_ensemble_outputs.Rd new file mode 100644 index 0000000..e6faa7b --- /dev/null +++ b/man/write_ensemble_outputs.Rd @@ -0,0 +1,24 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/pipeline.R +\name{write_ensemble_outputs} +\alias{write_ensemble_outputs} +\title{Write an ensemble's own artifacts} +\usage{ +write_ensemble_outputs(model, out) +} +\arguments{ +\item{model}{a \code{taupatch_ensemble} from \code{\link[=fit_patch_ensemble]{fit_patch_ensemble()}}} + +\item{out}{the run's output directory} +} +\value{ +\code{NULL}, invisibly +} +\description{ +What a single model has no equivalent of: which algorithms were fitted, how +each scored, what weight it was given, and whether it qualified. This is the +first thing to read after an ensemble run — a table showing one member at +0.9 weight and three near zero is a single model with extra steps, and only +this file says so. +} +\keyword{internal} diff --git a/man/write_jackknife.Rd b/man/write_jackknife.Rd new file mode 100644 index 0000000..b7cef15 --- /dev/null +++ b/man/write_jackknife.Rd @@ -0,0 +1,22 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/pipeline.R +\name{write_jackknife} +\alias{write_jackknife} +\title{Write the covariate jackknife table} +\usage{ +write_jackknife(jk, config) +} +\arguments{ +\item{jk}{the result of \code{\link[=jackknife_covariates]{jackknife_covariates()}}} + +\item{config}{a config list, as returned by \code{load_config()}} +} +\value{ +\code{NULL}, invisibly +} +\description{ +Written whether or not anything was dropped, and written before the model is +fitted, so a run that turned \code{drop} on leaves a record of what it removed and +on what evidence. +} +\keyword{internal} diff --git a/tests/testthat/test-ensemble.R b/tests/testthat/test-ensemble.R new file mode 100644 index 0000000..88a82fe --- /dev/null +++ b/tests/testthat/test-ensemble.R @@ -0,0 +1,394 @@ +# A member with just enough on it for build_ensemble() to score, weight and +# combine: an evaluation table, held-out predictions carrying their fold, a +# cutoff and an importance table. Cheaper than fitting four real models when +# what is under test is the combining rather than the fitting. +stub_member <- function(score, cutoff = 0.4, seed = 1) { + set.seed(seed) + n <- 60 + is_patch <- rep(c(TRUE, FALSE), c(15, 45))[sample.int(n)] + list( + evaluation = data.frame(metric = "tss", threshold = cutoff, value = score, + stringsAsFactors = FALSE), + predictions = data.frame( + .row = seq_len(n), + id = rep(paste0("Fold", 1:3), length.out = n), + patch = factor(ifelse(is_patch, "patch", "non_patch"), + levels = c("patch", "non_patch")), + .pred_patch = ifelse(is_patch, stats::runif(n, 0.4, 0.9), + stats::runif(n, 0.05, 0.5)), + stringsAsFactors = FALSE + ), + classification_threshold = cutoff, + importance = tibble::tibble(variable = c("SST", "SSS"), + importance = c(0.10, 0.02)), + predictors = c("SST", "SSS"), + model_data = data.frame(SST = stats::runif(n), SSS = stats::runif(n)), + threshold = 1000 + ) +} + +# The bootstrap is 2000 resamples by default and this fixture has 60 rows; the +# combining is what is under test, not the interval. +stub_config <- function() { + config <- mock_config() + config$model$bootstrap <- FALSE + config +} + +test_that("the ensemble is off unless a config asks for it", { + config <- mock_config() + + expect_null(ensemble_settings(config)) + + config$model$type <- "ensemble" + expect_setequal(ensemble_settings(config)$types, names(model_types())) + + config$model$type <- "rf" + config$model$ensemble <- list(types = c("rf", "glm")) + expect_equal(ensemble_settings(config)$types, c("rf", "glm")) + + config$model$ensemble <- FALSE + expect_null(ensemble_settings(config)) +}) + +test_that("a malformed ensemble block is refused", { + config <- mock_config() + + config$model$ensemble <- list(types = c("rf", "maxent")) + expect_error(ensemble_settings(config), "Unknown model.ensemble.types") + + config$model$ensemble <- list(types = "rf") + expect_error(ensemble_settings(config), "at least 2 model types") + + config$model$ensemble <- list(types = c("rf", "glm"), rule = "vote") + expect_error(ensemble_settings(config), "rule must be one of") + + config$model$ensemble <- list(types = c("rf", "glm"), weight_by = "aic") + expect_error(ensemble_settings(config), "weight_by must be one of") + + # An override for a type the ensemble does not fit is a typo, and silently + # ignoring it means the setting the user wanted never took effect. + config$model$ensemble <- list(types = c("rf", "glm"), + settings = list(gam = list(method = "REML"))) + expect_error(ensemble_settings(config), "overrides for types the ensemble") +}) + +test_that("fit_patch_model refuses an ensemble type rather than guessing one", { + config <- mock_config() + config$model$type <- "ensemble" + + expect_error(resolve_model_type(config), "fit_patch_ensemble") +}) + +test_that("a member config carries the run's settings with its own type set", { + config <- mock_config() + config$model$tune <- TRUE + settings <- ensemble_settings(within(config, model$type <- "ensemble")) %||% + list(types = names(model_types()), settings = list(gam = list(method = "REML"))) + settings$settings <- list(gam = list(method = "REML")) + + gam <- member_config(config, "gam", settings) + expect_equal(gam$model$type, "gam") + expect_equal(gam$model$method, "REML") + expect_null(gam$model$ensemble) + # A GAM has tunable parameters, so a run-level tune survives. + expect_true(gam$model$tune) + + # A GLM has none, so leaving tune on would make the ensemble refuse to fit + # the one member that is the honest baseline. + glm <- member_config(config, "glm", settings) + expect_false(glm$model$tune) + expect_true(validate_model(glm)) +}) + +test_that("weights are proportional to score and sum to one", { + scores <- c(rf = 0.6, brt = 0.4, glm = 0.2) + weights <- ensemble_weights(scores, c(TRUE, TRUE, TRUE)) + + expect_equal(sum(weights), 1) + expect_equal(unname(weights), c(0.5, 1 / 3, 1 / 6)) + expect_equal(names(weights), names(scores)) +}) + +test_that("a member that did not qualify gets zero weight, not a share", { + scores <- c(rf = 0.6, brt = 0.4, glm = 0.05) + weights <- ensemble_weights(scores, c(TRUE, TRUE, FALSE)) + + expect_equal(unname(weights[["glm"]]), 0) + expect_equal(sum(weights), 1) +}) + +test_that("a member predicting worse than chance cannot be given negative influence", { + # TSS runs from -1 to 1. A negative weight would make the ensemble + # deliberately invert that member rather than ignore it. + scores <- c(rf = 0.6, brt = -0.3) + weights <- ensemble_weights(scores, c(TRUE, TRUE)) + + expect_true(all(weights >= 0)) + expect_equal(sum(weights), 1) +}) + +test_that("equal weighting ignores the scores", { + weights <- ensemble_weights(c(rf = 0.9, glm = 0.1), c(TRUE, TRUE), + metric = "equal") + + expect_equal(unname(weights), c(0.5, 0.5)) +}) + +test_that("every combination rule lands on the 0-1 scale and means what it says", { + probabilities <- matrix(c(0.9, 0.1, 0.5, + 0.7, 0.2, 0.5, + 0.1, 0.3, 0.5), + nrow = 3, + dimnames = list(NULL, c("rf", "brt", "glm"))) + weights <- c(rf = 0.5, brt = 0.3, glm = 0.2) + cutoffs <- c(rf = 0.5, brt = 0.5, glm = 0.5) + + combined <- combine_members(probabilities, weights, cutoffs) + + expect_equal(combined$mean, rowMeans(probabilities)) + expect_equal(combined$weighted_mean, as.numeric(probabilities %*% weights)) + expect_equal(combined$median, apply(probabilities, 1, stats::median)) + # Committee averaging: the fraction of members calling the cell a patch at + # their own cutoff. Row 1 is (0.9, 0.7, 0.1), so 2 of 3; row 2 is + # (0.1, 0.2, 0.3), so none; row 3 is (0.5, 0.5, 0.5), so all three. + expect_equal(combined$committee, c(2 / 3, 0, 1)) + + for (rule in ensemble_rules()) { + expect_true(all(combined[[rule]] >= 0 & combined[[rule]] <= 1), info = rule) + } + expect_equal(combined$algorithm_range, + apply(probabilities, 1, max) - apply(probabilities, 1, min)) +}) + +test_that("each member binarises at its own cutoff for the committee vote", { + # A shared cutoff would score every member on a threshold that suits whichever + # happens to be best calibrated. + probabilities <- matrix(c(0.3, 0.3), nrow = 1, + dimnames = list(NULL, c("rf", "glm"))) + weights <- c(rf = 0.5, glm = 0.5) + + expect_equal(combine_members(probabilities, weights, + c(rf = 0.2, glm = 0.9))$committee, 0.5) + expect_equal(combine_members(probabilities, weights, + c(rf = 0.2, glm = 0.2))$committee, 1) +}) + +test_that("a member with no cutoff falls back to 0.5 rather than dropping out", { + probabilities <- matrix(c(0.6, 0.6), nrow = 1, + dimnames = list(NULL, c("rf", "glm"))) + combined <- combine_members(probabilities, c(rf = 0.5, glm = 0.5), + c(rf = NA_real_, glm = 0.5)) + + expect_equal(combined$committee, 1) +}) + +test_that("member scores are read at each member's own optimal cutoff", { + member <- list(evaluation = data.frame( + metric = c("roc_auc", "tss", "tss"), + threshold = c(NA_real_, 0.5, 0.18), + value = c(0.82, 0.11, 0.55), + stringsAsFactors = FALSE + )) + + # TSS at 0.5 is the misleading one when the classes are imbalanced, which + # they are here by construction. + expect_equal(member_score(member, "tss"), 0.55) + expect_equal(member_score(member, "roc_auc"), 0.82) + expect_equal(member_score(member, "equal"), 1) + expect_true(is.na(member_score(member, "pr_auc"))) +}) + +test_that("an ensemble fits several algorithms and combines them", { + skip_on_cran() + config <- mock_config() + config$model$trees <- 50 + config$model$cv_folds <- 3 + config$model$ensemble <- list(types = c("rf", "glm"), min_score = -1) + dat <- labeled_mock_data(config) + + ensemble <- suppressMessages(suppressWarnings(fit_patch_ensemble(dat, config))) + + expect_s3_class(ensemble, "taupatch_ensemble") + expect_setequal(names(ensemble$members), c("rf", "glm")) + expect_equal(ensemble$type, "ensemble") + expect_setequal(ensemble$summary$type, c("rf", "glm")) + expect_equal(sum(ensemble$weights), 1) + # A drop-in for a single model: the same fields, filled the same way. + expect_true(all(c("evaluation", "metrics", "importance", "predictors", + "classification_threshold") %in% names(ensemble))) + expect_true(is.numeric(ensemble$classification_threshold)) +}) + +test_that("the ensemble's evaluation is the ensemble's, not the average of its members", { + skip_on_cran() + config <- mock_config() + config$model$trees <- 50 + config$model$cv_folds <- 3 + config$model$ensemble <- list(types = c("rf", "glm"), min_score = -1) + dat <- labeled_mock_data(config) + + ensemble <- suppressMessages(suppressWarnings(fit_patch_ensemble(dat, config))) + + # It is computed from combined out-of-fold predictions, so it is a real + # cross-validated number and not a summary of summaries. + expect_true(all(c(".row", ".pred_patch", "patch") %in% + names(ensemble$predictions))) + expect_equal(nrow(ensemble$predictions), + nrow(ensemble$members$rf$predictions)) + + auc <- evaluation_value(ensemble, "roc_auc") + expect_equal(auc, + yardstick::roc_auc_vec(ensemble$predictions$patch, + ensemble$predictions$.pred_patch), + tolerance = 0.05) + # Per-fold, so it gets a standard error a pooled number could not have. + roc_row <- ensemble$metrics[ensemble$metrics$.metric == "roc_auc", ] + expect_equal(nrow(roc_row), 1) + expect_true(is.finite(roc_row$std_err)) +}) + +test_that("members are combined on the rows they were all held out of", { + skip_on_cran() + config <- mock_config() + config$model$trees <- 50 + config$model$cv_folds <- 3 + config$model$ensemble <- list(types = c("rf", "glm"), min_score = -1) + dat <- labeled_mock_data(config) + + ensemble <- suppressMessages(suppressWarnings(fit_patch_ensemble(dat, config))) + + # tune returns folds in its own order, so matching on position rather than + # on .row would silently pair each station with a different one. + rf <- ensemble$members$rf$predictions + aligned <- rf$patch[match(ensemble$predictions$.row, rf$.row)] + expect_equal(as.character(ensemble$predictions$patch), as.character(aligned)) +}) + +test_that("ensemble importance is weighted and keeps the per-member columns", { + skip_on_cran() + config <- mock_config() + config$model$trees <- 50 + config$model$cv_folds <- 3 + config$model$ensemble <- list(types = c("rf", "glm"), min_score = -1) + dat <- labeled_mock_data(config) + + ensemble <- suppressMessages(suppressWarnings(fit_patch_ensemble(dat, config))) + + expect_true(all(c("variable", "importance", "rf", "glm") %in% + names(ensemble$importance))) + # A predictor the forest leans on and the GLM ignores is a fact worth keeping, + # and the weighted average is the one number that hides it. + weights <- ensemble$weights[c("rf", "glm")] + expected <- as.numeric(as.matrix(ensemble$importance[c("rf", "glm")]) %*% + (weights / sum(weights))) + expect_equal(ensemble$importance$importance, expected) +}) + +test_that("an ensemble whose members all score badly is refused, not averaged", { + members <- list(rf = stub_member(0.05), glm = stub_member(0.02, seed = 2)) + settings <- list(types = c("rf", "glm"), rule = "mean", weight_by = "tss", + min_score = 0.4) + + expect_error(build_ensemble(members, stub_config(), settings), + "No ensemble member reached") +}) + +test_that("an ensemble left with one qualifying member says so", { + # Not an error: the map it produces is still the right map for that member. + # But an object called an ensemble that is one model has to announce itself, + # or the run reads as four algorithms agreeing. + members <- list(rf = stub_member(0.6), glm = stub_member(0.1, seed = 2)) + settings <- list(types = c("rf", "glm"), rule = "mean", weight_by = "tss", + min_score = 0.4) + + expect_warning(ensemble <- build_ensemble(members, stub_config(), settings), + "Only one ensemble member") + expect_equal(sum(ensemble$summary$qualifies), 1) + expect_equal(unname(ensemble$weights[["glm"]]), 0) +}) + +test_that("an ensemble projects and writes a layer per rule and per member", { + skip_on_cran() + config <- mock_config() + config$model$trees <- 50 + config$model$cv_folds <- 3 + config$model$ensemble <- list(types = c("rf", "glm"), min_score = -1, + rule = "weighted_mean") + config$projection$years <- c(2018, 2018) + config$projection$months <- c(6, 6) + generate_mock_zoop_data(config) + + result <- suppressMessages(suppressWarnings(run_taupatch(config))) + + expect_s3_class(result$model, "taupatch_ensemble") + expect_true(file.exists(file.path(config$paths$output_dir, + "ensemble_members.csv"))) + + stack <- terra::rast(result$projections$geotiff[1]) + expect_true("suitability" %in% names(stack)) + # Every rule the run did not pick is written beside the one it did, so a + # committee map can be read off the same file without refitting. + expect_true(all(c("suitability_mean", "suitability_median", + "suitability_committee") %in% names(stack))) + expect_true("algorithm_sd" %in% names(stack)) + expect_true(all(c("member_rf", "member_glm") %in% names(stack))) + + values <- terra::values(stack[["suitability"]]) + expect_true(all(values >= 0 & values <= 1, na.rm = TRUE)) +}) + +test_that("each member keeps its own effect plots rather than sharing an average", { + skip_on_cran() + config <- mock_config() + config$model$trees <- 50 + config$model$cv_folds <- 3 + config$model$ensemble <- list(types = c("rf", "glm"), min_score = -1) + generate_mock_zoop_data(config) + + suppressMessages(suppressWarnings(run_taupatch(config, project = FALSE))) + + members <- file.path(config$paths$output_dir, "diagnostics", "members") + expect_true(dir.exists(file.path(members, "rf"))) + expect_true(dir.exists(file.path(members, "glm"))) + # A GLM has coefficients and a forest does not, which is the reason these are + # written per member instead of once. + expect_true(file.exists(file.path(members, "glm", "coefficients.csv"))) +}) + +test_that("an algorithm ensemble and a resample ensemble are different columns", { + skip_on_cran() + config <- mock_config() + config$model$trees <- 50 + config$model$cv_folds <- 3 + config$model$ensemble <- list(types = c("rf", "glm"), min_score = -1) + config$projection$uncertainty <- TRUE + config$projection$years <- c(2018, 2018) + config$projection$months <- c(6, 6) + generate_mock_zoop_data(config) + + result <- suppressMessages(suppressWarnings(run_taupatch(config))) + stack <- terra::rast(result$projections$geotiff[1]) + + # Both, side by side. They measure different things and a reader must be able + # to tell which is which. + expect_true("algorithm_sd" %in% names(stack)) + expect_true("suitability_sd" %in% names(stack)) + expect_true("novelty" %in% names(stack)) +}) + +test_that("printing an ensemble says which members carried it", { + ensemble <- structure(list( + rule = "weighted_mean", + summary = data.frame(type = c("rf", "glm"), score = c(0.6, 0.4), + metric = "tss", qualifies = c(TRUE, TRUE), + weight = c(0.6, 0.4), stringsAsFactors = FALSE), + evaluation = data.frame(metric = "roc_auc", threshold = NA_real_, + value = 0.83), + classification_threshold = 0.21 + ), class = "taupatch_ensemble") + + expect_output(print(ensemble), "taupatch ensemble") + expect_output(print(ensemble), "weighted_mean") + expect_output(print(ensemble), "0.83") +}) diff --git a/tests/testthat/test-jackknife.R b/tests/testthat/test-jackknife.R new file mode 100644 index 0000000..9a14bde --- /dev/null +++ b/tests/testthat/test-jackknife.R @@ -0,0 +1,323 @@ +test_that("the jackknife is off unless a config asks for it", { + config <- mock_config() + + expect_null(jackknife_settings(config)) + + config$covariates$jackknife <- TRUE + expect_type(jackknife_settings(config), "list") + + config$covariates$jackknife <- FALSE + expect_null(jackknife_settings(config)) + + config$covariates$jackknife <- list(enabled = FALSE, drop = TRUE) + expect_null(jackknife_settings(config)) +}) + +test_that("dropping covariates is off by default", { + # The one default that matters most here. A test that removes a covariate + # from someone's model without being asked is worse than no test. + config <- mock_config() + config$covariates$jackknife <- TRUE + + expect_false(jackknife_settings(config)$drop) + + config$covariates$jackknife <- list(metric = "pr_auc", alpha = 0.1) + expect_false(jackknife_settings(config)$drop) + + config$covariates$jackknife <- list(drop = TRUE) + expect_true(jackknife_settings(config)$drop) +}) + +test_that("a malformed jackknife block is refused at load", { + config <- mock_config() + + config$covariates$jackknife <- list(metric = "kappa") + expect_error(jackknife_settings(config), "must be 'roc_auc' or 'pr_auc'") + + config$covariates$jackknife <- list(alpha = 1.5) + expect_error(jackknife_settings(config), "alpha must be between 0 and 1") + + config$covariates$jackknife <- list(adjust = "sidak") + expect_error(jackknife_settings(config), "adjust must be one of") + + config$covariates$jackknife <- list(criterion = "vibes") + expect_error(jackknife_settings(config), "must be 'fold' or 'parametric'") + + config$covariates$jackknife <- list(type = "maxent") + expect_error(jackknife_settings(config), "Unknown covariates.jackknife.type") + + config$covariates$jackknife <- "yes please" + expect_error(jackknife_settings(config), "must be true, false, or a block") +}) + +test_that("a parametric criterion is refused for a model type that has no likelihood", { + config <- mock_config() + config$model$type <- "rf" + config$covariates$jackknife <- list(criterion = "parametric") + + expect_error(validate_jackknife(config), "only exists for a 'glm' or 'gam'") + + # Naming a type that does have one makes it legal again. + config$covariates$jackknife$type <- "glm" + expect_true(validate_jackknife(config)) +}) + +test_that("a jackknife needs enough folds to test across", { + config <- mock_config() + config$covariates$jackknife <- TRUE + config$model$cv_folds <- 2 + + expect_error(validate_jackknife(config), "at least 3 model.cv_folds") +}) + +test_that("the corrected paired test is more conservative than the naive one", { + # The whole reason the correction is there. Folds share training data, so the + # naive paired t treats k highly dependent numbers as k independent ones and + # reports significance that is not there. + differences <- c(0.03, 0.01, 0.04, 0.02, 0.025, 0.015, 0.035, 0.02, 0.03, 0.01) + + corrected <- corrected_paired_test(differences) + naive <- stats::t.test(differences, alternative = "greater") + + expect_gt(corrected$p_value, naive$p.value) + expect_lt(abs(corrected$statistic), abs(naive$statistic)) + # Specifically, the variance is inflated by 1/k + 1/(k-1) rather than 1/k. + k <- length(differences) + expect_equal(corrected$std_err, + sqrt(stats::var(differences) * (1 / k + 1 / (k - 1)))) + expect_equal(corrected$estimate, mean(differences)) + expect_equal(corrected$df, k - 1) +}) + +test_that("the paired test is one-sided in the direction that matters", { + # A covariate whose removal *improves* the model has failed the test, not + # passed a different one. + improves <- corrected_paired_test(c(-0.05, -0.04, -0.06, -0.05, -0.03)) + hurts <- corrected_paired_test(c(0.05, 0.04, 0.06, 0.05, 0.03)) + + expect_gt(improves$p_value, 0.5) + expect_lt(hurts$p_value, 0.5) +}) + +test_that("the paired test declines to answer with too few folds", { + result <- corrected_paired_test(c(0.02, 0.03)) + + expect_true(is.na(result$p_value)) + expect_equal(result$n, 2) +}) + +test_that("a covariate that helped identically on every fold is not called noise", { + # Zero variance is a division by zero, and the wrong answer there is to + # report NA for a covariate that helped by the same amount ten times running. + result <- corrected_paired_test(rep(0.04, 8)) + + expect_equal(result$p_value, 0) + expect_equal(result$estimate, 0.04) + + # And a covariate that did exactly nothing, every time, is not significant. + expect_equal(corrected_paired_test(rep(0, 8))$p_value, 1) +}) + +test_that("jackknife_dropped respects the keep list and the floor", { + jk <- data.frame( + variable = c("SST", "SSS", "CHL", "MLD"), + contribution = c(0.08, 0.002, 0.001, 0.0005), + significant = c(TRUE, FALSE, FALSE, FALSE), + stringsAsFactors = FALSE + ) + + expect_equal(jackknife_dropped(jk, list(keep = character(), min_predictors = 1)), + c("CHL", "MLD", "SSS")) + + # A covariate on the keep list is never dropped, whatever the test said. + expect_equal(jackknife_dropped(jk, list(keep = "SSS", min_predictors = 1)), + c("CHL", "MLD")) + + # The floor binds, and the ones given back are the strongest of the failures. + expect_equal(jackknife_dropped(jk, list(keep = character(), min_predictors = 3)), + "MLD") + expect_equal(jackknife_dropped(jk, list(keep = character(), min_predictors = 4)), + character()) +}) + +test_that("a covariate whose test could not be computed is never dropped", { + # NA is not evidence of absence. This is checked at the point the flag is + # set, since jackknife_dropped() only sees the flag. + jk <- data.frame(variable = c("SST", "SSS"), contribution = c(0.1, 0.01), + significant = c(TRUE, TRUE), stringsAsFactors = FALSE) + + expect_equal(jackknife_dropped(jk, list(keep = character(), min_predictors = 1)), + character()) +}) + +test_that("dropped covariates go through covariates.exclude", { + # Rather than a second mechanism beside it: a dropped covariate must still be + # fetched, since a derived covariate may need it as an ingredient. + config <- mock_config() + config$covariates$exclude <- "uo" + + updated <- apply_jackknife_drop(config, c("SSS", "uo")) + + expect_setequal(updated$covariates$exclude, c("uo", "SSS")) + expect_identical(apply_jackknife_drop(config, character()), config) +}) + +test_that("the jackknife runs end to end and every covariate gets a verdict", { + skip_on_cran() + config <- mock_config() + config$model$trees <- 50 + config$model$cv_folds <- 3 + dat <- labeled_mock_data(config) + + jk <- suppressMessages(jackknife_covariates(dat, config)) + + predictors <- predictor_names(dat, config) + expect_setequal(jk$variable, predictors) + expect_true(all(c("score_full", "score_without", "score_only", "contribution", + "p_value", "p_adjusted", "significant") %in% names(jk))) + # Ordered by contribution, most first. + expect_equal(jk$contribution, sort(jk$contribution, decreasing = TRUE)) + # Every model was scored on the same folds, so the full model's score is one + # number rather than one per covariate. + expect_length(unique(jk$score_full), 1) + expect_type(jk$significant, "logical") + expect_false(any(is.na(jk$significant))) +}) + +test_that("the leave-one-out and only-one halves answer different questions", { + skip_on_cran() + config <- mock_config() + config$model$trees <- 50 + config$model$cv_folds <- 3 + dat <- labeled_mock_data(config) + + jk <- suppressMessages(jackknife_covariates(dat, config)) + + # Both are AUCs on the same folds, so both are on the 0-1 scale, and the + # contribution is the difference the leave-one-out half measures. + expect_true(all(jk$score_without >= 0 & jk$score_without <= 1)) + expect_true(all(jk$score_only >= 0 & jk$score_only <= 1)) + expect_equal(jk$contribution, jk$score_full - jk$score_without, + tolerance = 1e-8) +}) + +test_that("a GLM jackknife reports a likelihood ratio test beside the fold test", { + skip_on_cran() + config <- mock_config() + config$model$type <- "glm" + config$model$cv_folds <- 3 + dat <- labeled_mock_data(config) + + jk <- suppressMessages(jackknife_covariates(dat, config)) + + expect_true(all(jk$parametric_test == "LRT")) + expect_false(all(is.na(jk$parametric_p))) + finite <- jk$parametric_p[!is.na(jk$parametric_p)] + expect_true(all(finite >= 0 & finite <= 1)) +}) + +test_that("a forest jackknife has no parametric column to fill in", { + skip_on_cran() + config <- mock_config() + config$model$type <- "rf" + config$model$trees <- 50 + config$model$cv_folds <- 3 + dat <- labeled_mock_data(config) + + jk <- suppressMessages(jackknife_covariates(dat, config)) + + # NA rather than a number that looks like a p-value: there is no likelihood + # here, so there is no test, and saying so is the honest column. + expect_true(all(is.na(jk$parametric_p))) + expect_true(all(is.na(jk$parametric_test))) + # But the fold test is there, which is the point of having it. + expect_false(all(is.na(jk$p_value))) +}) + +test_that("the criterion chooses which p-value decides significance", { + jk_from <- function(criterion, p_fold, p_parametric) { + out <- data.frame(p_adjusted = p_fold, parametric_p = p_parametric) + criterion_values <- if (identical(criterion, "parametric")) { + out$parametric_p + } else { + out$p_adjusted + } + is.na(criterion_values) | criterion_values < 0.05 + } + + expect_equal(jk_from("fold", c(0.01, 0.9), c(0.9, 0.01)), c(TRUE, FALSE)) + expect_equal(jk_from("parametric", c(0.01, 0.9), c(0.9, 0.01)), c(FALSE, TRUE)) +}) + +test_that("a run with drop off reports but does not remove", { + skip_on_cran() + config <- mock_config() + config$model$trees <- 50 + config$model$cv_folds <- 3 + config$projection$years <- c(2018, 2018) + config$projection$months <- c(6, 6) + # An alpha of 1 fails nothing; an alpha this small fails everything, which is + # what makes "and nothing was removed" a real assertion. + config$covariates$jackknife <- list(alpha = 1e-12, drop = FALSE) + generate_mock_zoop_data(config) + + result <- suppressMessages(run_taupatch(config, project = FALSE)) + + expect_false(is.null(result$jackknife)) + expect_true(any(!result$jackknife$significant)) + expect_null(result$config$covariates$exclude) + expect_setequal(result$model$predictors, predictor_names(result$data, config)) + expect_true(file.exists(file.path(config$paths$output_dir, + "covariate_jackknife.csv"))) +}) + +test_that("a run with drop on removes what the test rejected", { + skip_on_cran() + config <- mock_config() + config$model$trees <- 50 + config$model$cv_folds <- 3 + config$covariates$jackknife <- list(alpha = 1e-12, drop = TRUE, + min_predictors = 1, keep = "jday") + generate_mock_zoop_data(config) + + result <- suppressMessages(run_taupatch(config, project = FALSE)) + + dropped <- result$config$covariates$exclude + expect_true(length(dropped) > 0) + # The keep list survives the drop, and the model was fitted on what was left. + expect_false("jday" %in% dropped) + expect_true("jday" %in% result$model$predictors) + expect_length(intersect(result$model$predictors, dropped), 0) +}) + +test_that("workers resolve to something runnable", { + expect_equal(resolve_workers(FALSE, 8), 1L) + expect_equal(resolve_workers(1, 8), 1L) + # Never more workers than tasks. + expect_lte(resolve_workers(64, 3, quiet = TRUE), 3L) + expect_gte(resolve_workers(NULL, 8, quiet = TRUE), 1L) + expect_error(resolve_workers(0, 8), "positive count") + expect_error(resolve_workers("many", 8), "positive count") +}) + +test_that("a parallel map gives the same answer as a sequential one", { + skip_on_os("windows") + + square <- function(x) x^2 + expect_equal(taupatch_lapply(1:6, square, workers = 1), + taupatch_lapply(1:6, square, workers = 2)) +}) + +test_that("a failing worker surfaces the failure rather than a truncated result", { + skip_on_os("windows") + + # mclapply warns about the failed call on its own way out; the assertion is + # that the error is re-raised rather than a short list being returned. + expect_error( + suppressWarnings( + taupatch_lapply(1:4, function(x) if (x == 3) stop("nope") else x, + workers = 2) + ), + "A parallel worker failed" + ) +}) From e62aa0c354b4d5ec4e344c15d5778ccfba435052 Mon Sep 17 00:00:00 2001 From: chross22 <52218551+chross22@users.noreply.github.com> Date: Mon, 10 Aug 2026 11:48:35 -0400 Subject: [PATCH 03/10] Work the jackknife and the ensemble into the vignette, and bump to 0.2.0 Two new sections in the vignette, both running on the mock data the rest of it uses, and worked examples in the README beside the config blocks that were already there. The vignette sections say what the mock run actually shows rather than what would have been convenient. No covariate reaches significance on three covariates over a few hundred synthetic stations at five folds, and two contributions come out negative - so the text says that this is the correct answer rather than a broken test, and that a method returning confident answers from this much data would be the one to distrust. The `jday` row is the good illustration of the point the two halves exist to make: it scores well on its own and contributes nothing on top of the others, which is information that is real and duplicated rather than information that is not there. The ensemble section needed no such hedging. On this run the combined model beats both of its members out of fold, and the per-member importance columns show the forest and the GLM disagreeing sharply about one covariate, which is exactly what those columns are kept for. jackknife_dropped() now defaults its settings argument to the ones the result carries, so asking a jackknife what it would drop takes only the result. The alternative in a vignette was reaching for an internal. Version goes to 0.2.0. Two features, no exported function changed meaning; uncertainty_layers() became projection_layers() but was never exported. inst/CITATION reads the version from DESCRIPTION, so it needs no edit. Co-Authored-By: Claude Opus 5 --- DESCRIPTION | 2 +- R/jackknife.R | 11 ++- README.md | 75 ++++++++++++++- man/jackknife_dropped.Rd | 6 +- tests/testthat/test-jackknife.R | 19 ++++ vignettes/taupatch.Rmd | 165 ++++++++++++++++++++++++++++++++ 6 files changed, 268 insertions(+), 10 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 816a443..0e198dc 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: taupatch Title: Spatial Habitat Suitability Models for Zooplankton High-Abundance Patches -Version: 0.1.0 +Version: 0.2.0 Authors@R: person("Camille", "Ross", email = "camille.ross@maine.edu", role = c("aut", "cre"), comment = c(ORCID = "0000-0002-1428-2294")) diff --git a/R/jackknife.R b/R/jackknife.R index 4d1e6bf..9890776 100644 --- a/R/jackknife.R +++ b/R/jackknife.R @@ -617,7 +617,9 @@ gam_term_p_values <- function(fit, predictors) { #' rather than an arbitrary one. #' #' @param jk the result of [jackknife_covariates()] -#' @param settings from [jackknife_settings()] +#' @param settings from [jackknife_settings()]. Defaults to the settings the +#' jackknife was actually run under, which it carries on itself — so asking a +#' result what it would drop needs nothing but the result. #' @return a character vector of covariate names, possibly empty #' @examples #' jk <- data.frame( @@ -631,7 +633,12 @@ gam_term_p_values <- function(fit, predictors) { #' jackknife_dropped(jk, list(keep = "CHL", min_predictors = 1)) #' @seealso [jackknife_covariates()] #' @export -jackknife_dropped <- function(jk, settings) { +jackknife_dropped <- function(jk, settings = attr(jk, "settings")) { + if (is.null(settings)) { + stop("No jackknife settings given, and `jk` does not carry any. Pass the ", + "result of jackknife_settings(), or a list with `keep` and ", + "`min_predictors`.", call. = FALSE) + } keep <- as.character(settings$keep %||% character()) floor <- as.integer(settings$min_predictors %||% 2L) diff --git a/README.md b/README.md index 1838e6e..8ef0b82 100644 --- a/README.md +++ b/README.md @@ -722,14 +722,44 @@ and still available to anything that needs it as an ingredient — a gradient's velocity components, say — it just stops being a predictor. The run's returned `config` shows exactly what came out. -Called directly, it needs no config block at all: +Called directly it needs no config block at all, which is the way to use it +interactively: ```r -jk <- jackknife_covariates(dat, config) -jk[c("variable", "contribution", "p_adjusted", "significant")] -jackknife_dropped(jk, jackknife_settings(config)) # what drop would remove +result <- run_taupatch(config, project = FALSE) +jk <- jackknife_covariates(result$data, config) + +jk[c("variable", "score_without", "score_only", "contribution", "p_adjusted", + "significant")] +#> variable score_without score_only contribution p_adjusted significant +#> 1 SST 0.791 0.812 0.066 0.020 TRUE +#> 2 DEPTH 0.828 0.774 0.029 0.124 FALSE +#> 3 jday 0.855 0.611 0.002 0.402 FALSE + +jackknife_dropped(jk) # what drop would remove, had it been on +#> [1] "DEPTH" "jday" ``` +Or from the config, as part of a run. `result$jackknife` is the same table, and +`result$config` shows what the run actually fitted on: + +```r +config$covariates$jackknife <- list(drop = TRUE, keep = "jday", workers = 4) +result <- run_taupatch(config) +#> Jackknifing covariates... +#> jackknifing 3 covariates: 7 cross-validations across 4 workers +#> no detectable contribution at alpha = 0.05 (holm-adjusted): DEPTH, jday +#> covariates.jackknife.drop is on: removing DEPTH + +result$config$covariates$exclude +#> [1] "DEPTH" +``` + +`jday` failed the test and stayed, because `keep` named it. With `drop` left at +its default the last two lines read `covariates.jackknife.drop is off, so +nothing is removed. Setting it would drop: DEPTH, jday` — and the run fits on +everything. + Parallelism forks, which Windows does not have, so it runs sequentially there and says so. It is one model fit per covariate per fold either way — minutes, not hours — it just does not get faster with more cores. @@ -841,7 +871,42 @@ model: ``` Everything after the fit works the same either way, so a config that turns this -on gets ensemble projections without changing anything else. +on gets ensemble projections without changing anything else: + +```r +config$model$type <- "ensemble" +result <- run_taupatch(config) +#> Fitting model... +#> fitting 4 ensemble members (rf, brt, glm, gam) across 4 workers +#> ROC AUC: 0.8916 +#> Projecting monthly suitability... +#> combining 4 algorithms by weighted_mean; the other rules and each +#> member's own surface go beside it + +result$model +#> +#> rule: weighted_mean +#> members (4 of 4 qualifying): +#> type score metric qualifies weight +#> brt 0.6120990 tss TRUE 0.31 +#> rf 0.6043118 tss TRUE 0.30 +#> gam 0.5412287 tss TRUE 0.27 +#> glm 0.2381044 tss TRUE 0.12 +#> +#> ensemble ROC AUC (out of fold): 0.8916 +#> classification threshold: 0.08881 +``` + +`fit_patch_ensemble()` does the same thing without the pipeline around it, and +the result is a drop-in for a `fit_patch_model()` one: + +```r +ensemble <- fit_patch_ensemble(dat, config) +ensemble$summary # scores, weights, who qualified +ensemble$evaluation # the ensemble's own, out of fold +ensemble$members$gam # each member, entire +project_patch_model(ensemble, env_dat, config) +``` **Four ways to combine.** All of them are computed and written on every run; `rule` picks which one becomes the `suitability` layer, and the others go beside diff --git a/man/jackknife_dropped.Rd b/man/jackknife_dropped.Rd index 97c054c..95c0971 100644 --- a/man/jackknife_dropped.Rd +++ b/man/jackknife_dropped.Rd @@ -4,12 +4,14 @@ \alias{jackknife_dropped} \title{Which covariates a jackknife would drop} \usage{ -jackknife_dropped(jk, settings) +jackknife_dropped(jk, settings = attr(jk, "settings")) } \arguments{ \item{jk}{the result of \code{\link[=jackknife_covariates]{jackknife_covariates()}}} -\item{settings}{from \code{\link[=jackknife_settings]{jackknife_settings()}}} +\item{settings}{from \code{\link[=jackknife_settings]{jackknife_settings()}}. Defaults to the settings the +jackknife was actually run under, which it carries on itself — so asking a +result what it would drop needs nothing but the result.} } \value{ a character vector of covariate names, possibly empty diff --git a/tests/testthat/test-jackknife.R b/tests/testthat/test-jackknife.R index 9a14bde..e0b557c 100644 --- a/tests/testthat/test-jackknife.R +++ b/tests/testthat/test-jackknife.R @@ -140,6 +140,25 @@ test_that("jackknife_dropped respects the keep list and the floor", { character()) }) +test_that("a jackknife result knows what settings it ran under", { + # So asking a result what it would drop needs nothing but the result, which + # is how it gets used interactively. + jk <- data.frame(variable = c("SST", "SSS", "CHL"), + contribution = c(0.08, 0.002, 0.001), + significant = c(TRUE, FALSE, FALSE), + stringsAsFactors = FALSE) + attr(jk, "settings") <- list(keep = "SSS", min_predictors = 1) + + expect_equal(jackknife_dropped(jk), "CHL") + # An explicit argument still wins over the carried one. + expect_equal(jackknife_dropped(jk, list(keep = character(), + min_predictors = 1)), + c("CHL", "SSS")) + + attr(jk, "settings") <- NULL + expect_error(jackknife_dropped(jk), "does not carry any") +}) + test_that("a covariate whose test could not be computed is never dropped", { # NA is not evidence of absence. This is checked at the point the flag is # set, since jackknife_dropped() only sees the flag. diff --git a/vignettes/taupatch.Rmd b/vignettes/taupatch.Rmd index 2fccab4..b07803a 100644 --- a/vignettes/taupatch.Rmd +++ b/vignettes/taupatch.Rmd @@ -306,6 +306,167 @@ list.files(config$paths$output_dir, recursive = TRUE)[1:20] refitting. `diagnostics/cv_predictions.csv` holds the held-out predictions, so any metric not in `evals.csv` can be computed without refitting either. +## Which covariates are earning their place + +Importance ranks the covariates a model has. It cannot say whether a covariate +is carrying anything the others were not already carrying — a predictor can rank +third and be entirely redundant. `jackknife_covariates()` answers that by +refitting without each one, over the same folds, and seeing how much worse the +model ranks stations. + +It needs no config block; the defaults are used when there is none. + +```{r jackknife, message = FALSE} +jk <- jackknife_covariates(result$data, config) + +jk[, c("variable", "score_full", "score_without", "score_only", "contribution")] +``` + +Two halves, two questions. **`score_without`** is low when the covariate carries +something no other covariate has — its *unique* contribution, which is what +`contribution` measures. **`score_only`** is the model on that covariate alone, +so it is high when the covariate carries a lot whether or not anything else +carries it too. + +Read the two columns against each other and the mock run says something the +importance table above could not. `jday` scores well on its own — the synthetic +data has a real seasonal cycle in it — and contributes nothing on top of the +others, which have absorbed the same seasonality. That is information which is +real *and* duplicated, and it is a very different situation from information +that is not there. Ranking on `contribution` alone would treat the two +identically. + +Some contributions come out negative. That means the model scored *better* on +average without the covariate, which on this many stations is noise rather than +a finding — and it is exactly the case the test exists to keep you from +over-reading. + +```{r jackknife-test} +jk[, c("variable", "contribution", "statistic", "p_value", "p_adjusted", + "significant")] +``` + +Nothing here is significant, and that is the right answer rather than a broken +test. Three covariates over a few hundred synthetic stations, scored on five +folds, is not enough evidence to establish that any one of them is load-bearing. +A test that returned confident answers from this much data would be the one to +distrust. + +`p_value` is one-sided on the per-fold differences, with the variance inflated +to account for the folds sharing most of their training rows. At five folds that +inflation is a factor of 1.5 on the standard error, so every statistic here is +two-thirds of what a naive paired *t*-test would have reported. It makes no +difference to the verdict on this run — these differences are nowhere near the +boundary either way — but it is the difference between a real finding and a +false one when a covariate lands close to it, which on real data is where the +interesting ones land. `p_adjusted` then accounts for having asked once per +covariate. A GLM or a GAM would also fill in `parametric_p` with a likelihood +ratio test; a forest has no likelihood, so it is `NA` here. + +**Nothing is dropped.** That is the default, and it is deliberate: + +```{r jackknife-drop} +jackknife_settings(config) # NULL: this config has no jackknife block + +# What dropping *would* have removed, under the settings the test ran with. +jackknife_dropped(jk) +``` + +Only one name comes back even though all three failed, because +`min_predictors` floors the model at two covariates and the ones handed back are +the strongest of those that failed. + +A covariate that fails this test is one the *others already account for on these +stations*, which is a statement about collinearity in this sample at least as +much as about ecology. Depth and surface temperature carry much of the same +information on a shelf, and the test will call either one redundant depending on +which the model reached for first. Turning `covariates.jackknife.drop` on tells +the pipeline to act on it anyway, and `keep` protects a covariate that is in the +model because the study is about it: + +```{r jackknife-config, eval = FALSE} +config$covariates$jackknife <- list(drop = TRUE, keep = "jday", workers = 4) +result <- run_taupatch(config) + +result$jackknife # the table, also written to the run +result$config$covariates$exclude # what came out +``` + +## Fitting all the models at once + +Rather than choosing an algorithm, fit several and combine them. Set +`model.type` to `ensemble` and the pipeline does the rest — everything after the +fit works the same, so nothing else in the config changes. + +Four types are available; two are used here because `xgboost` and `mgcv` are +suggested rather than required. + +```{r ensemble-fit, message = FALSE, warning = FALSE} +config$model$ensemble <- list(types = c("rf", "glm"), rule = "weighted_mean") +config$paths$output_dir <- file.path(tempdir(), "taupatch_vignette_ens") + +result <- run_taupatch(config) +result$model +``` + +Each member's score, and the weight it earned, is the first thing to read. A +table showing one member near 1.0 and the rest near zero is a single model with +extra steps. + +```{r ensemble-summary} +result$model$summary +``` + +The ensemble's evaluation is **its own**, not the average of its members'. Every +member was fitted on the same folds from the same seed, so their held-out +predictions line up row for row, and the combined out-of-fold predictions go +through the same evaluation as a single model's: + +```{r ensemble-eval} +result$model$evaluation[1:2, c("metric", "value", "lower", "upper")] + +# Each member's own, for comparison. +vapply(result$model$members, + function(m) m$evaluation$value[m$evaluation$metric == "roc_auc" & + is.na(m$evaluation$threshold)], + numeric(1)) +``` + +That distinction matters. Combining members that make *different* mistakes beats +all of them; combining members that make the same mistakes does not, and only a +cross-validated ensemble prediction can tell those apart. + +Importance is weighted across members, with each member's own kept beside it — a +predictor the forest leans on and the GLM ignores is a fact about the shape of +the relationship, and the average is the one number that hides it. + +```{r ensemble-importance} +result$model$importance +``` + +Every combination rule is computed and written, so a committee map can be read +off the same file without refitting. `rule` only picks which one becomes the +`suitability` layer: + +```{r ensemble-layers} +names(terra::rast(result$projections$geotiff[1])) +``` + +`algorithm_sd` is the new one, and it is **not** the same as `suitability_sd` +from the section above. That one is a single algorithm refitted on resampled +stations; this one is different algorithms looking at the same stations and +disagreeing. Both can be on at once, and a cell can be quiet on one and loud on +the other: + +```{r ensemble-spread} +ens <- readr::read_csv( + file.path(config$paths$output_dir, "projections", "suitability.csv"), + show_col_types = FALSE +) + +summary(ens[, c("suitability", "algorithm_sd", "suitability_sd")]) +``` + ## Moving to real data Three things change, and nothing else does. @@ -366,6 +527,10 @@ diagnostics each one contributes. In R: - `?covariate_transforms` — and what can be done to one - `?model_types` — the four models, and what each is good for - `?partial_effects`, `?gam_smooth_terms` — reading a fitted model +- `?jackknife_covariates`, `?jackknife_settings` — testing which covariates earn + their place +- `?fit_patch_ensemble`, `?ensemble_settings`, `?ensemble_rules` — combining + several model types ## Citing a run From 2b97fbc3fac9220ec8aedc3fec1af8f05ed8ebc6 Mon Sep 17 00:00:00 2001 From: chross22 <52218551+chross22@users.noreply.github.com> Date: Mon, 10 Aug 2026 11:56:09 -0400 Subject: [PATCH 04/10] Retire the tests for the citation script this repo no longer has bdb47bf moved the citation check onto the shared engine in chross22/distsamp and deleted inst/tools/check_citations.R. tests/testthat/test-citations.R stayed. Every one of its seven tests calls a helper that sources that script and skips when it is absent, so since that commit all seven have skipped on every run - and a skipped test reads as a passing test in most summaries. The file was worth having for a reason it no longer serves. Its comment says the bug worth guarding against is the checker passing everything, "a matcher that is too loose reports a healthy repo whatever the state of the references, which is worse than not having the check at all". That is true, and it is now distsamp's matcher to guard. Deleted rather than moved. Nothing in the file is about taupatch: all seven exercise the retired script's internals through `sys.source()`, so there is nothing here to port. distsamp already covers two of the seven in tests/testthat/test-citation-engine.R - the DOI styles and the trailing punctuation - and whether the other five still describe how the shared engine behaves is a question about that engine, not about this repo. test-citation.R, singular, is untouched. It tests inst/CITATION and the README's citation section, is offline, and has been passing throughout. Also adds the empty-config-field case to the YAML boolean tests. A concurrent review of the same fix raised it: the walk rebuilds every list it descends into, and rebuilding by assigning back - `x[] <- lapply(x, ...)` - drops elements whose value is NULL, which is what an empty YAML field parses to. This implementation replaces the list instead and was never affected, but the hazard is real enough for the next person to deserve a test rather than a comment. Co-Authored-By: Claude Opus 5 --- tests/testthat/test-citations.R | 124 -------------------------------- tests/testthat/test-config.R | 17 +++++ 2 files changed, 17 insertions(+), 124 deletions(-) delete mode 100644 tests/testthat/test-citations.R diff --git a/tests/testthat/test-citations.R b/tests/testthat/test-citations.R deleted file mode 100644 index 7ea32db..0000000 --- a/tests/testthat/test-citations.R +++ /dev/null @@ -1,124 +0,0 @@ -# The citation checker is a script rather than package code, so it gets no -# coverage from anything else. Its parsing is what these cover; the network side -# is exercised by the scheduled workflow, which is the only place it can be. -# -# The bug worth guarding against is the checker passing everything: a matcher -# that is too loose reports a healthy repo whatever the state of the references, -# which is worse than not having the check at all. - -# Under `R CMD check` the tests run from a directory that has no source tree -# above it, so the script is found where it is installed. `sys.source()` loads -# the definitions without running `main()`, which is guarded on being sourced. -checker <- function() { - path <- system.file("tools", "check_citations.R", package = "taupatch") - if (!nzchar(path)) path <- test_path("..", "..", "inst", "tools", "check_citations.R") - skip_if_not(file.exists(path), "check_citations.R not in this build") - - env <- new.env() - sys.source(path, envir = env) - env -} - -# The repo root, for the one test that reads the real files. Absent under -# `R CMD check`, where only the installed package exists. -repo_root <- function() { - root <- test_path("..", "..") - skip_if_not(file.exists(file.path(root, "DESCRIPTION")) && - file.exists(file.path(root, "README.md")), - "not running from a source checkout") - root -} - -test_that("DOIs are found in every style the docs write them in", { - env <- checker() - - text <- paste( - "markdown: [doi:10.3354/meps14204](https://doi.org/10.3354/meps14204).", - "roxygen: \\doi{10.1093/biomet/87.4.954}", - "bare in a bibentry: doi = \"10.1214/aos/1013203451\",", - "dataset: \\doi{10.48670/moi-00021}", - sep = "\n" - ) - - expect_setequal( - env$extract_dois(text), - c("10.3354/meps14204", "10.1093/biomet/87.4.954", - "10.1214/aos/1013203451", "10.48670/moi-00021") - ) -}) - -test_that("trailing punctuation is not read as part of the identifier", { - env <- checker() - - # A DOI at the end of a sentence, inside a link, and inside a roxygen macro. - # Each closing character belongs to the prose. - expect_equal(env$extract_dois("see \\doi{10.1214/ss/1177013604}"), - "10.1214/ss/1177013604") - expect_equal(env$extract_dois("(https://doi.org/10.1023/A:1010933404324)."), - "10.1023/A:1010933404324") -}) - -test_that("doi.org links are left to the DOI check rather than counted twice", { - env <- checker() - - urls <- env$extract_urls( - "[doi:10.3354/meps14204](https://doi.org/10.3354/meps14204) and - " - ) - expect_equal(urls, "https://www.naturalearthdata.com/") -}) - -test_that("the year is read from the entry the DOI is in, not the one above it", { - env <- checker() - - # A package citation prints no year. Without a boundary the lookback walks - # back into the previous list item and reports a mismatch that isn't. - text <- paste( - "- Fisher A, Rudin C, Dominici F (2019). All models are wrong.", - " ", - "- Chang W, Cheng J. *shiny: Web Application Framework for R*.", - " \\doi{10.32614/CRAN.package.shiny}", - sep = "\n" - ) - - expect_true(is.na(env$year_near(text, "10.32614/CRAN.package.shiny"))) -}) - -test_that("the year is found when the entry does print one", { - env <- checker() - - text <- paste( - "- Breiman L (2001). Random forests. *Machine Learning* **45**(1), 5-32.", - " \\doi{10.1023/A:1010933404324}", - sep = "\n" - ) - expect_equal(env$year_near(text, "10.1023/A:1010933404324"), 2001L) -}) - -test_that("a dead handle is told apart from a live one", { - env <- checker() - - # Response code 1 is "found" and 100 is "handle not found". `:1` is a prefix - # of `:100`, so a loose match passes every dead DOI - which would make the - # whole check report success unconditionally. - live <- '{"responseCode":1,"handle":"10.3354/meps14204"}' - dead <- '{"responseCode":100,"handle":"10.5281/zenodo.7657585"}' - pattern <- '"responseCode"[[:space:]]*:[[:space:]]*1[[:space:]]*[,}]' - - expect_true(grepl(pattern, live)) - expect_false(grepl(pattern, dead)) -}) - -test_that("the files it reads are the ones that make citation claims", { - env <- checker() - files <- basename(env$source_files(repo_root())) - - expect_true("README.md" %in% files) - expect_true("CITATION" %in% files) - expect_true("covariate_catalog.R" %in% files) - # The vignette links out too, and was missed when it was first added. - expect_true(any(grepl("\\.Rmd$", files))) - # man/ is generated from R/, so including it would report every DOI twice and - # point at a file nobody edits. - expect_false(any(grepl("\\.Rd$", files))) -}) diff --git a/tests/testthat/test-config.R b/tests/testthat/test-config.R index 190b0ab..fa08f8d 100644 --- a/tests/testthat/test-config.R +++ b/tests/testthat/test-config.R @@ -72,6 +72,23 @@ test_that("a boolean in a value position is still a boolean", { expect_identical(parsed$flags, c(TRUE, FALSE)) }) +test_that("an empty config field survives the boolean walk", { + # The walk rebuilds every list it descends into, and an empty YAML field + # parses to NULL. Rebuilding by assigning back into the list - `x[] <- lapply` + # rather than replacing it - drops NULL elements entirely, which would make + # `"mtry" %in% names(config$model)` false for a field the file does mention. + path <- write_yaml_text(c("model:", " mtry:", " tune: false", + "covariates:", " transform:", + " normalize: true")) + + parsed <- read_config_yaml(path) + + expect_equal(names(parsed$model), c("mtry", "tune")) + expect_null(parsed$model$mtry) + expect_equal(names(parsed$covariates), c("transform", "normalize")) + expect_identical(parsed$covariates$normalize, TRUE) +}) + test_that("a quoted string that spells a boolean stays a string", { path <- write_yaml_text(c("species:", " active: 'true'", " label: \"no\"", "note: not a bool")) From 311d262260bdf028b70d9cbfc80d71ac2203a4d3 Mon Sep 17 00:00:00 2001 From: chross22 <52218551+chross22@users.noreply.github.com> Date: Mon, 10 Aug 2026 12:04:49 -0400 Subject: [PATCH 05/10] Check the n-key fix against a lag that actually runs The tests already here establish that `n: 2` survives the parse as a field named n. This asserts the consequence: that the step derivoce then runs is a two-step lag, with August carrying June's value and the first two months undefined rather than the first. That is the failure the bug actually produced. A config could parse correctly and still be wired to a reader that never asked for the field, and nothing above this test would notice. The config has to come off disk for either to mean anything. Built as an R list, `n` is an ordinary name and the test passes against a parser that loses it - so mock_config_yaml() splices the derivoce block into the shipped mock config and loads the result, which is the only path that goes through yaml::read_yaml(). Verified it bites: with load_config() pointed back at a bare read_yaml(), this fails three assertions and then errors with "step 'lag_covariate' has no argument(s) FALSE", which is the bug stating its own name. Ported from a concurrent session that fixed the same bug independently. Its own implementation is superseded by 354ceb7, which additionally keeps a sequence of bare booleans - [yes, no] - as logicals rather than strings; this test is the part of that work worth keeping. Co-Authored-By: Claude Opus 5 --- tests/testthat/helper-config.R | 23 +++++++++++++++++++++++ tests/testthat/test-derivoce.R | 30 ++++++++++++++++++++++++++++++ 2 files changed, 53 insertions(+) diff --git a/tests/testthat/helper-config.R b/tests/testthat/helper-config.R index 249c7eb..62adcab 100644 --- a/tests/testthat/helper-config.R +++ b/tests/testthat/helper-config.R @@ -11,6 +11,29 @@ mock_config <- function(dir = tempfile("taupatch")) { config } +# The mock config with extra lines spliced into its `covariates:` block, loaded +# from disk. Building the same config as an R list would not exercise the YAML +# parse, which is where a field name can be lost before load_config() ever sees +# it. +mock_config_yaml <- function(covariate_lines, dir = tempfile("taupatch")) { + dir.create(dir, recursive = TRUE, showWarnings = FALSE) + source_path <- system.file("configs", "mock_test.yaml", package = "taupatch") + if (!nzchar(source_path)) { + source_path <- testthat::test_path("..", "..", "inst", "configs", "mock_test.yaml") + } + + text <- readLines(source_path) + anchor <- grep("^ source: mock$", text) + stopifnot(length(anchor) == 1) + path <- file.path(dir, "spliced.yaml") + writeLines(append(text, covariate_lines, after = anchor), path) + + config <- load_config(path) + config$paths$zoop_file <- file.path(dir, "mock_zooplankton.csv") + config$paths$output_dir <- file.path(dir, "output") + config +} + # Modeling data as fit_patch_model() expects it: stations with covariates # attached and patches labeled. Runs the same path the pipeline does, so tests # exercise real data rather than a hand-built frame that might drift from it. diff --git a/tests/testthat/test-derivoce.R b/tests/testthat/test-derivoce.R index 0a9487b..183bf2a 100644 --- a/tests/testthat/test-derivoce.R +++ b/tests/testthat/test-derivoce.R @@ -81,6 +81,36 @@ test_that("derived covariates carry real values, not just columns", { expect_equal(sort(flat$SST_lag1[flat$MONTH == 7]), sort(flat$SST[flat$MONTH == 6])) }) +test_that("a lag written 'n: 2' in a config file lags by two steps", { + skip_if_not_installed("derivoce") + # YAML 1.1 resolves a bare `n` to the boolean false, keys included, so this + # step used to parse to a field named FALSE and fall back to a lag of one - + # invisibly, since a one-step lag is an ordinary thing to ask for. The config + # has to come off disk: built as an R list, `n` is just a name and the test + # would pass against a parser that loses it. + config <- mock_config_yaml(c( + " derivoce:", + " - type: lag_covariate", + " vars: [SST]", + " n: 2" + )) + + spec <- config$covariates$derivoce[[1]] + expect_equal(spec$n, 2) + expect_null(spec[["FALSE"]]) + expect_equal(derivoce_names(config), "SST_lag2") + + env <- fetch_covariates(config, years = 2018, months = 6:8) + derived <- suppressMessages(add_derivoce_covariates(env, config)) + flat <- sf::st_drop_geometry(derived) + + # Two steps back: undefined for the first two months rather than the first, + # and August carries June's value. + expect_true(all(is.na(flat$SST_lag2[flat$MONTH %in% c(6, 7)]))) + expect_false(any(is.na(flat$SST_lag2[flat$MONTH == 8]))) + expect_equal(sort(flat$SST_lag2[flat$MONTH == 8]), sort(flat$SST[flat$MONTH == 6])) +}) + test_that("steps chain, so a step can read what an earlier one produced", { skip_if_not_installed("derivoce") # current_speed then a gradient of `speed` is the original pipeline's uv_grad. From ddfbcf71bbc2e958ce00f9cd924aaaa6eacba60c Mon Sep 17 00:00:00 2001 From: chross22 <52218551+chross22@users.noreply.github.com> Date: Mon, 10 Aug 2026 12:29:15 -0400 Subject: [PATCH 06/10] Register the references the jackknife and the ensemble added The citation check failed on this branch: six DOIs cited but absent from tools/citations.csv. The shared engine requires every DOI in the docs to be in the registry, and the registry is not generated - a new reference is a new row. Five of the six are now rows, with the fields taken from CrossRef rather than typed, since those are the fields the engine compares against CrossRef: araujo2007 ensemble forecasting - why an ensemble of algorithms marmion2009 the combination rules, compared elith2011 the leave-one-out / only-one jackknife pair nadeau2003 the variance correction dietterich1998 the Type I error of the uncorrected test Two details. The registry `year` is the issue year the docs cite, not CrossRef's, which is a year earlier for the two published online first - the engine allows exactly that gap. And `first_author` is written unaccented, Araujo, because the comparison strips every non-ASCII character to a space on both sides: an accented registry entry would normalise to "ara jo" and never match CrossRef's "araujo". The sixth, Bouckaert and Frank, keeps its reference but loses its DOI link. It is a Springer chapter, so its identifier ends _3, and the engine's DOI pattern treats an underscore as a terminator - reasonable in Markdown, where _ is emphasis, but it means that DOI cannot be written down anywhere the checker reads. Registering the truncated form would have made CI pass while recording the book's identifier for a chapter citation, so the reference now gives the volume and page range instead and says why in a note. Verified by running the shared engine against this working tree rather than by pushing and waiting: all five sections pass. That is also how the failure should have been caught the first time. The local script this branch retired only resolved DOIs and had no concept of a registry, so "every citation resolves" was true and beside the point. Co-Authored-By: Claude Opus 5 --- R/jackknife.R | 7 +++++-- README.md | 7 ++++--- man/corrected_paired_test.Rd | 7 +++++-- tools/citations.csv | 5 +++++ 4 files changed, 19 insertions(+), 7 deletions(-) diff --git a/R/jackknife.R b/R/jackknife.R index 9890776..ee9dc7e 100644 --- a/R/jackknife.R +++ b/R/jackknife.R @@ -481,8 +481,11 @@ subset_workflow <- function(train, vars, config, type) { #' #' Bouckaert RR, Frank E (2004). Evaluating the replicability of significance #' tests for comparing learning algorithms. *Advances in Knowledge Discovery -#' and Data Mining*, 3-12. \doi{10.1007/978-3-540-24775-3_3} — the correction -#' applied to k-fold specifically +#' and Data Mining* (PAKDD 2004), Lecture Notes in Computer Science 3056, +#' 3-12. Springer. — the correction applied to k-fold specifically. Cited +#' without its DOI deliberately: it is a Springer chapter, so the identifier +#' contains an underscore, and the citation checker's DOI pattern treats one as +#' a terminator. #' @keywords internal corrected_paired_test <- function(differences) { usable <- differences[is.finite(differences)] diff --git a/README.md b/README.md index 8ef0b82..da6bcf2 100644 --- a/README.md +++ b/README.md @@ -1358,9 +1358,10 @@ Earth](https://www.naturalearthdata.com/), public domain, via `rnaturalearth`. needs a variance correction rather than a better estimator - Bouckaert RR, Frank E (2004). Evaluating the replicability of significance tests for comparing learning algorithms. *Advances in Knowledge Discovery and - Data Mining*, 3–12. - [doi:10.1007/978-3-540-24775-3_3](https://doi.org/10.1007/978-3-540-24775-3_3) - — the correction applied to k-fold specifically + Data Mining* (PAKDD 2004), Lecture Notes in Computer Science **3056**, 3–12. + Springer. — the correction applied to k-fold specifically. Listed without its + DOI on purpose: a Springer chapter identifier contains an underscore, which + the citation checker's DOI pattern reads as the end of the identifier - Dietterich TG (1998). Approximate statistical tests for comparing supervised classification learning algorithms. *Neural Computation* **10**(7), 1895–1923. [doi:10.1162/089976698300017197](https://doi.org/10.1162/089976698300017197) — diff --git a/man/corrected_paired_test.Rd b/man/corrected_paired_test.Rd index dae11f4..ca885c1 100644 --- a/man/corrected_paired_test.Rd +++ b/man/corrected_paired_test.Rd @@ -31,7 +31,10 @@ Learning} \strong{52}(3), 239-281. \doi{10.1023/A:1024068626366} Bouckaert RR, Frank E (2004). Evaluating the replicability of significance tests for comparing learning algorithms. \emph{Advances in Knowledge Discovery -and Data Mining}, 3-12. \doi{10.1007/978-3-540-24775-3_3} — the correction -applied to k-fold specifically +and Data Mining} (PAKDD 2004), Lecture Notes in Computer Science 3056, +3-12. Springer. — the correction applied to k-fold specifically. Cited +without its DOI deliberately: it is a Springer chapter, so the identifier +contains an underscore, and the citation checker's DOI pattern treats one as +a terminator. } \keyword{internal} diff --git a/tools/citations.csv b/tools/citations.csv index 01379b6..db3ef07 100644 --- a/tools/citations.csv +++ b/tools/citations.csv @@ -35,3 +35,8 @@ wickham2016,10.1007/978-3-319-24277-4,,Wickham,2016,ggplot2,Use R!,,,crossref, toolboxdocsmarinecopernicuseu,,https://toolbox-docs.marine.copernicus.eu/,,,,,,,url,Generated from a scanned URL - fill in the fields by hand wwwnceinoaagov,,https://www.ncei.noaa.gov/archive/accession/0187513,,,,,,,url,Generated from a scanned URL - fill in the fields by hand helpmarinecopernicuseu,,https://help.marine.copernicus.eu/en/collections/4060068-copernicus-marine-toolbox,,,,,,,url,Generated from a scanned URL - fill in the fields by hand +araujo2007,10.1016/j.tree.2006.09.010,,Araujo,2007,Ensemble forecasting of species distributions,Trends in Ecology & Evolution,22,42-47,crossref, +marmion2009,10.1111/j.1472-4642.2008.00491.x,,Marmion,2009,Evaluation of consensus methods in predictive species distribution modelling,Diversity and Distributions,15,59-69,crossref, +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, From d5d1f31fb25e7ec5844c06fa0891cff0f7aebcea Mon Sep 17 00:00:00 2001 From: chross22 <52218551+chross22@users.noreply.github.com> Date: Mon, 10 Aug 2026 13:00:15 -0400 Subject: [PATCH 07/10] Keep the worker count under what a package check permits R CMD check --as-cran sets _R_CHECK_LIMIT_CORES_, and under it parallel::mclapply() does not quietly use fewer cores: parallel:::.check_ncores() stops with "3 simultaneous processes spawned" the moment more than two are asked for. The default here was one fewer than the machine's physical cores, so on any runner with four the jackknife errored - taking six tests and the vignette down with it. Windows and macOS passed the same commit, which is the part worth noting. Windows never forks, so resolve_workers() had already returned 1 there, and macOS reported few enough cores to stay under the limit. Only Linux had enough cores to trip it, so this looked like a platform bug and was a default that is wrong everywhere and visible in one place. core_ceiling() now reports what the session permits, and resolve_workers() takes the minimum of the request, the task count and that ceiling. It caps an explicitly configured count too: the limit is not a preference, and mclapply refuses rather than clamping, so honouring `workers: 8` under a check would only reproduce the error the user was trying to avoid. Two workers compute the same answer as eight. options(mc.cores) now sets the default as well, ahead of detectCores(). It is the knob R users already reach for, and a package that ignores it makes them learn a second one. Verified under _R_CHECK_LIMIT_CORES_=TRUE, where the six tests that failed on CI now pass, and with R CMD check --as-cran, where the vignette that failed to re-build now does. That flag is the one my earlier local check was missing: without it the variable is never set and the whole failure mode is invisible. withr moves to Suggests, for the tests that set the variable and restore it. Co-Authored-By: Claude Opus 5 --- DESCRIPTION | 3 ++- R/parallel.R | 38 +++++++++++++++++++++++++++------ man/core_ceiling.Rd | 25 ++++++++++++++++++++++ man/resolve_workers.Rd | 8 ++++--- tests/testthat/test-jackknife.R | 33 ++++++++++++++++++++++++++++ 5 files changed, 97 insertions(+), 10 deletions(-) create mode 100644 man/core_ceiling.Rd diff --git a/DESCRIPTION b/DESCRIPTION index 0e198dc..938a487 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -50,7 +50,8 @@ Suggests: shiny, shinyFiles, rnaturalearth, - testthat (>= 3.0.0) + testthat (>= 3.0.0), + withr Remotes: chross22/datamatch, chross22/derivoce, diff --git a/R/parallel.R b/R/parallel.R index 5cad5b9..c374468 100644 --- a/R/parallel.R +++ b/R/parallel.R @@ -70,10 +70,12 @@ taupatch_lapply <- function(x, fun, workers = 1L, seed = NULL) { #' #' `NULL` or `true` means "as many as this machine can spare", which is one #' fewer than its physical cores — leaving one is what keeps the session it was -#' launched from responsive. `false` or `1` is sequential. +#' launched from responsive. `false` or `1` is sequential. `options(mc.cores=)` +#' overrides the default, since that is the option R users already reach for. #' -#' Capped at `n`, since a task list of six cannot use twelve workers, and forced -#' to 1 on Windows, where [taupatch_lapply()] cannot fork. +#' Capped at `n`, since a task list of six cannot use twelve workers; forced to +#' 1 on Windows, where [taupatch_lapply()] cannot fork; and capped again at +#' whatever [core_ceiling()] allows. #' #' @param workers the configured value: `NULL`, a logical, or a count #' @param n how many tasks there are to spread @@ -82,8 +84,10 @@ taupatch_lapply <- function(x, fun, workers = 1L, seed = NULL) { #' @keywords internal resolve_workers <- function(workers = NULL, n = 1L, quiet = FALSE) { requested <- if (is.null(workers) || isTRUE(workers)) { - available <- suppressWarnings(parallel::detectCores(logical = FALSE)) - if (is.na(available)) 1L else max(1L, available - 1L) + getOption("mc.cores") %||% { + available <- suppressWarnings(parallel::detectCores(logical = FALSE)) + if (is.na(available)) 1L else max(1L, available - 1L) + } } else if (isFALSE(workers)) { 1L } else { @@ -102,5 +106,27 @@ resolve_workers <- function(workers = NULL, n = 1L, quiet = FALSE) { } return(1L) } - max(1L, min(requested, as.integer(n))) + # The ceiling is applied to a configured count as well as to the default. It + # is not a preference to be overridden - `parallel` refuses outright above it. + max(1L, min(as.integer(requested), as.integer(n), core_ceiling())) +} + +#' The most workers this session is allowed to spawn +#' +#' `R CMD check --as-cran` sets `_R_CHECK_LIMIT_CORES_`, and under it +#' `parallel::mclapply()` does not quietly use fewer cores — it **errors**, via +#' `parallel:::.check_ncores()`, the moment more than two are asked for. So a +#' default of "cores minus one" turns every jackknife into a failure on any +#' machine with four or more cores, in exactly the context where a package is +#' most likely to be run by someone other than its author. +#' +#' Capping rather than erroring is right here: the caller asked for a jackknife, +#' not for a particular number of processes, and two workers computes the same +#' answer as eight. +#' +#' @return `2L` under a core-limited check, otherwise `Inf` +#' @keywords internal +core_ceiling <- function() { + limit <- Sys.getenv("_R_CHECK_LIMIT_CORES_", "") + if (nzchar(limit) && !identical(tolower(limit), "false")) 2L else Inf } diff --git a/man/core_ceiling.Rd b/man/core_ceiling.Rd new file mode 100644 index 0000000..b1d7d9a --- /dev/null +++ b/man/core_ceiling.Rd @@ -0,0 +1,25 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/parallel.R +\name{core_ceiling} +\alias{core_ceiling} +\title{The most workers this session is allowed to spawn} +\usage{ +core_ceiling() +} +\value{ +\code{2L} under a core-limited check, otherwise \code{Inf} +} +\description{ +\verb{R CMD check --as-cran} sets \verb{_R_CHECK_LIMIT_CORES_}, and under it +\code{parallel::mclapply()} does not quietly use fewer cores — it \strong{errors}, via +\code{parallel:::.check_ncores()}, the moment more than two are asked for. So a +default of "cores minus one" turns every jackknife into a failure on any +machine with four or more cores, in exactly the context where a package is +most likely to be run by someone other than its author. +} +\details{ +Capping rather than erroring is right here: the caller asked for a jackknife, +not for a particular number of processes, and two workers computes the same +answer as eight. +} +\keyword{internal} diff --git a/man/resolve_workers.Rd b/man/resolve_workers.Rd index 6c659bc..3d86622 100644 --- a/man/resolve_workers.Rd +++ b/man/resolve_workers.Rd @@ -19,10 +19,12 @@ an integer worker count, at least 1 \description{ \code{NULL} or \code{true} means "as many as this machine can spare", which is one fewer than its physical cores — leaving one is what keeps the session it was -launched from responsive. \code{false} or \code{1} is sequential. +launched from responsive. \code{false} or \code{1} is sequential. \code{options(mc.cores=)} +overrides the default, since that is the option R users already reach for. } \details{ -Capped at \code{n}, since a task list of six cannot use twelve workers, and forced -to 1 on Windows, where \code{\link[=taupatch_lapply]{taupatch_lapply()}} cannot fork. +Capped at \code{n}, since a task list of six cannot use twelve workers; forced to +1 on Windows, where \code{\link[=taupatch_lapply]{taupatch_lapply()}} cannot fork; and capped again at +whatever \code{\link[=core_ceiling]{core_ceiling()}} allows. } \keyword{internal} diff --git a/tests/testthat/test-jackknife.R b/tests/testthat/test-jackknife.R index e0b557c..5e77310 100644 --- a/tests/testthat/test-jackknife.R +++ b/tests/testthat/test-jackknife.R @@ -319,6 +319,39 @@ test_that("workers resolve to something runnable", { expect_error(resolve_workers("many", 8), "positive count") }) +test_that("a core-limited check caps the workers instead of erroring", { + # `R CMD check --as-cran` sets this, and under it parallel::mclapply() does + # not use fewer cores - it errors outright above two. A default of + # `cores - 1` therefore turns every jackknife into a failure on any machine + # with four or more cores, which is how this reached CI the first time. + withr::local_envvar(c("_R_CHECK_LIMIT_CORES_" = "TRUE")) + + expect_equal(core_ceiling(), 2L) + expect_lte(resolve_workers(NULL, 16, quiet = TRUE), 2L) + # A configured count is capped too. The limit is not a preference. + expect_lte(resolve_workers(8, 16, quiet = TRUE), 2L) + # And sequential is still reachable. + expect_equal(resolve_workers(FALSE, 16), 1L) +}) + +test_that("the core ceiling lifts when the check variable is absent or false", { + withr::local_envvar(c("_R_CHECK_LIMIT_CORES_" = "")) + expect_identical(core_ceiling(), Inf) + + withr::local_envvar(c("_R_CHECK_LIMIT_CORES_" = "false")) + expect_identical(core_ceiling(), Inf) +}) + +test_that("options(mc.cores) sets the default worker count", { + # The option R users already reach for, rather than a taupatch-only knob. + withr::local_envvar(c("_R_CHECK_LIMIT_CORES_" = "")) + withr::local_options(mc.cores = 2) + + expect_equal(resolve_workers(NULL, 16, quiet = TRUE), 2L) + # An explicit argument still wins over the option. + expect_equal(resolve_workers(1, 16), 1L) +}) + test_that("a parallel map gives the same answer as a sequential one", { skip_on_os("windows") From e92ff7872b8e4710e54c0e66296c8ba42b4d9c1e Mon Sep 17 00:00:00 2001 From: chross22 <52218551+chross22@users.noreply.github.com> Date: Mon, 10 Aug 2026 13:23:28 -0400 Subject: [PATCH 08/10] Let the mc.cores test know that Windows cannot fork The assertion was that resolve_workers(NULL, ...) returns the two workers options(mc.cores) asked for. On Windows it returns one, because there is no forking there and resolve_workers() says so a few lines earlier - so the test described the platforms it happened to be written on and failed on the one it did not. Windows was the only green Linux-and-macOS-failing job on the previous commit and the only red one on this, for the same underlying reason in reverse: worker counts differ by platform, and an expectation on an exact count has to say which platform it means. The other assertions added alongside it are inequalities - at most two under a core-limited check, at most the task count, at least one - and every one of them holds at a worker count of one, so this was the only one that needed it. Checked rather than assumed. Co-Authored-By: Claude Opus 5 --- tests/testthat/test-jackknife.R | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/tests/testthat/test-jackknife.R b/tests/testthat/test-jackknife.R index 5e77310..f794ec1 100644 --- a/tests/testthat/test-jackknife.R +++ b/tests/testthat/test-jackknife.R @@ -347,7 +347,12 @@ test_that("options(mc.cores) sets the default worker count", { withr::local_envvar(c("_R_CHECK_LIMIT_CORES_" = "")) withr::local_options(mc.cores = 2) - expect_equal(resolve_workers(NULL, 16, quiet = TRUE), 2L) + # Windows cannot fork, so it is sequential whatever the option asks for. The + # expectation has to know that: asserting 2 everywhere passes on the + # platforms that fork and fails on the one that does not, which is a test + # describing the author's laptop rather than the function. + forks <- !identical(.Platform$OS.type, "windows") + expect_equal(resolve_workers(NULL, 16, quiet = TRUE), if (forks) 2L else 1L) # An explicit argument still wins over the option. expect_equal(resolve_workers(1, 16), 1L) }) From 89d76233e2c97ac3c58b06e1f3a6f2329df05272 Mon Sep 17 00:00:00 2001 From: chross22 <52218551+chross22@users.noreply.github.com> Date: Mon, 10 Aug 2026 14:45:42 -0400 Subject: [PATCH 09/10] Spell the YAML boolean marker in ASCII, so config.R is a text file again The marker carrying a boolean's source text through parsing was built out of control characters, on the reasoning that nothing in a YAML file could collide with them. Nothing can. But two SOH bytes in a string literal are enough for file(1) to classify the whole of R/config.R as `data` rather than as source, and anything that asks the same question - an editor, a diff viewer, a blob preview - then presents a 700-line file as binary. It reads as a corrupted file, which is a poor trade for a collision that was never going to happen. The marker is now , in plain ASCII, and the file is UTF-8 text again. Nothing is given up. Only the two boolean handlers ever prepend the marker, and `yaml` only ever hands them scalars it has itself resolved as booleans - a quoted "true" carries no boolean tag and is never marked, which is already covered by a test. The one way left to collide is a config value that genuinely begins with the string . Co-Authored-By: Claude Opus 5 --- R/config.R | 18 ++++++++++++------ 1 file changed, 12 insertions(+), 6 deletions(-) diff --git a/R/config.R b/R/config.R index e501b99..ecaeece 100644 --- a/R/config.R +++ b/R/config.R @@ -77,16 +77,22 @@ read_config_yaml <- function(path) { #' The marker that carries a boolean's source text through parsing #' -#' 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 `"true"`, which must survive as that text and not become a logical. -#' -#' 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 #' `yaml` collapsing a sequence of scalars into an atomic vector, which drops #' attributes. `[true, false]` would otherwise come back as two strings. #' +#' Deliberately plain ASCII. The first version used control characters, on the +#' reasoning that nothing could collide with them. That was true, and it made +#' `file(1)` report the whole of `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 +#' `"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. +#' #' @keywords internal -yaml_bool_marker <- "taupatch-bool" +yaml_bool_marker <- "" #' Turn marked scalars back into keys and logicals #' From 474a588dd1274e34250d547059a982c59143ba54 Mon Sep 17 00:00:00 2001 From: chross22 <52218551+chross22@users.noreply.github.com> Date: Mon, 10 Aug 2026 14:49:28 -0400 Subject: [PATCH 10/10] Say how to get the vignette, and how a species catalog is built Two things the README left a reader to work out. install_github() does not build vignettes unless asked - build_vignettes defaults to FALSE - so vignette("taupatch") finds nothing after following the install instructions exactly. Nothing is broken and nothing says so. The install section now gives the flag, warns that building runs the pipeline as it renders, and links the source for anyone who would rather just read it. The species section explained what a catalog looks like without explaining what decides its shape. It showed column_prefix throughout and never abundance_column, so the two forms read as a style choice rather than as an answer to "does this database resolve life stages for this taxon?". There is now a table for the two forms and what each lets a run do, and the reason column_prefix needs stage columns to match against: it matches _, which is also why `stages` alongside abundance_column is an error. It also never mentioned species_catalog_from(), which reads a database and picks the right form per taxon so nobody has to. Its output is now shown for a header with one staged taxon and one total, since seeing the two entries side by side is the explanation. And a short table separates it from zoop_taxa() and available_species(), which answer neighbouring questions and are easy to reach for by mistake. Co-Authored-By: Claude Opus 5 --- README.md | 70 ++++++++++++++++++++++++++++++++++++++++++++++++++++--- 1 file changed, 67 insertions(+), 3 deletions(-) diff --git a/README.md b/README.md index da6bcf2..b6868e0 100644 --- a/README.md +++ b/README.md @@ -116,6 +116,22 @@ Tabs, in the order the questions come up: remotes::install_github("chross22/taupatch") ``` +**To get the vignette, ask for it.** `install_github()` does not build vignettes +by default — it is faster not to, and most installs do not want them — so +`vignette("taupatch")` finds nothing after a plain install. That is the install +being economical, not the vignette being missing: + +```r +remotes::install_github("chross22/taupatch", build_vignettes = TRUE) +vignette("taupatch") +``` + +Building it needs `knitr` and `rmarkdown`, and it runs the whole pipeline on the +mock data as it renders, so expect it to take a minute rather than a moment. If +you would rather not rebuild, [read it on +GitHub](vignettes/taupatch.Rmd) instead — the same document, without the +rendered figures. + Environmental data needs the [Copernicus Marine Toolbox](https://help.marine.copernicus.eu/en/collections/4060068-copernicus-marine-toolbox) (EU Copernicus Marine Service 2025) installed and configured with your Copernicus @@ -347,9 +363,57 @@ long the Copernicus download took to get there. ### Species and life stages -Each species names the prefix of its columns in the database. `column_prefix` -defaults to the species key, so only an aliased name needs it — `pcal` is the one -case, since its columns are `pseudo_*`: +A catalog entry answers one question: **which column or columns hold this +species' abundance?** There are two ways to answer it, and which one an entry +uses is not a style choice — it depends on how the database reports the taxon. + +| entry form | use it when | what a run can then do | +|---|---|---| +| `column_prefix: cfin` | the database resolves life stages, as `cfin_CV`, `cfin_CVI`, … | select particular `stages`, or leave `stages` out to sum every one | +| `abundance_column: CENTROPAGES_TYPICUS` | the database reports one total and no stages | nothing to select — the column *is* the abundance | + +`column_prefix` matches `_`, which is why it needs stage +columns to match against; `abundance_column` names a single column outright. +Setting both is an error, and so is setting `stages` alongside +`abundance_column`, because there are no stage columns for it to narrow to. + +**You do not have to work this out per taxon.** `species_catalog_from()` reads +the database and writes the block, choosing the right form for each: + +```r +header <- c("CALANUS_FINMARCHICUS_CV_10M2", "CALANUS_FINMARCHICUS_CVI_10M2", + "CENTROPAGES_TYPICUS_10M2") + +species_catalog_from(header, aliases = c(cfin = "CALANUS_FINMARCHICUS")) +#> $cfin +#> $cfin$column_prefix # stages resolved, so a prefix +#> [1] "CALANUS_FINMARCHICUS" +#> $cfin$threshold +#> $cfin$threshold$type +#> [1] "percentile" +#> $cfin$threshold$value +#> [1] 0.9 +#> +#> $ctyp +#> $ctyp$abundance_column # a total only, so a column +#> [1] "CENTROPAGES_TYPICUS" +``` + +`aliases` is what lets a config say `cfin` instead of `CALANUS_FINMARCHICUS`; +without it the key is a shorthand derived from the taxon name. The result goes +straight into [`generate_config()`](#three-ways-to-get-one) as `species`. + +Three functions, three different questions — it is worth knowing which you want: + +| | question | +|---|---| +| `zoop_taxa(raw_export)` | which taxa does this **raw** export carry, with units still on the names? | +| `available_species(database)` | which taxa could I model from this **formatted** database? | +| `species_catalog_from(header)` | write me the catalog block for them | + +Once you have a catalog, `column_prefix` defaults to the species key, so only an +aliased name needs it — `pcal` is the one case in the ECOMON database, since its +columns are `pseudo_*`: ```yaml species: