diff --git a/DESCRIPTION b/DESCRIPTION index 3337f95..938a487 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")) @@ -23,6 +23,7 @@ Imports: datamatch, dplyr, ggplot2, + parallel, parsnip, recipes, rsample, @@ -49,7 +50,8 @@ Suggests: shiny, shinyFiles, rnaturalearth, - testthat (>= 3.0.0) + testthat (>= 3.0.0), + withr Remotes: chross22/datamatch, chross22/derivoce, 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 3ba2bee..ecaeece 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) @@ -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) @@ -37,6 +38,103 @@ 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 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 <- "" + +#' 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 @@ -191,6 +289,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..ee9dc7e --- /dev/null +++ b/R/jackknife.R @@ -0,0 +1,712 @@ +#' 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* (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)] + 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()]. 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( +#' 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 = 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) + + 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..c374468 --- /dev/null +++ b/R/parallel.R @@ -0,0 +1,132 @@ +#' 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. `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; 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 +#' @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)) { + 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 { + 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) + } + # 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/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 d5ce538..b6868e0 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) @@ -114,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 @@ -345,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: @@ -478,7 +544,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 +554,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 @@ -631,6 +704,130 @@ 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, which is the way to use it +interactively: + +```r +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. + ### Model type Four models, chosen with one word: @@ -711,6 +908,129 @@ 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: + +```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 +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 @@ -848,18 +1168,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) @@ -884,6 +1208,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 @@ -1034,6 +1361,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 @@ -1052,6 +1383,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) — @@ -1080,6 +1416,29 @@ 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* (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) — + 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 edf7a92..909b023 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 @@ -160,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/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/corrected_paired_test.Rd b/man/corrected_paired_test.Rd new file mode 100644 index 0000000..ca885c1 --- /dev/null +++ b/man/corrected_paired_test.Rd @@ -0,0 +1,40 @@ +% 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} (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/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..95c0971 --- /dev/null +++ b/man/jackknife_dropped.Rd @@ -0,0 +1,43 @@ +% 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 = attr(jk, "settings")) +} +\arguments{ +\item{jk}{the result of \code{\link[=jackknife_covariates]{jackknife_covariates()}}} + +\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 +} +\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/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/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..3d86622 --- /dev/null +++ b/man/resolve_workers.Rd @@ -0,0 +1,30 @@ +% 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. \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; 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/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/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/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/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/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/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-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 d55cca3..fa08f8d 100644 --- a/tests/testthat/test-config.R +++ b/tests/testthat/test-config.R @@ -24,6 +24,109 @@ 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("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")) + + 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()) 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. 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..f794ec1 --- /dev/null +++ b/tests/testthat/test-jackknife.R @@ -0,0 +1,380 @@ +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 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. + 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 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) + + # 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) +}) + +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" + ) +}) 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, 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