Take what fancyfx does better, and drop what it duplicates

DESCRIPTION

@@ -22,6 +22,7 @@ Depends:
 Imports:
     datamatch,
     dplyr,
+    fancyfx,
     ggplot2,
     parallel,
     parsnip,
@@ -38,7 +39,6 @@ Imports:
     yardstick
 Suggests:
     derivoce,
-    fancyfx,
     knitr,
     leaflet,
     ranger,

NAMESPACE

@@ -21,6 +21,7 @@ export(default_copernicus_datasets)
 export(default_species_catalog)
 export(derivoce_choices)
 export(derivoce_covariates)
+export(derivoce_dependency_columns)
 export(derivoce_required_inputs)
 export(derivoce_steps_for)
 export(ensemble_rules)
@@ -69,11 +70,13 @@ export(power_curve)
 export(prejoin_steps)
 export(project_patch_model)
 export(projection_map)
+export(projection_overlap)
 export(raw_abundance_suffix)
 export(run_taupatch)
 export(run_taupatch_app)
 export(save_config)
 export(single_stages)
+export(spatial_bias)
 export(species_catalog_from)
 export(split_dates)
 export(stage_suffix_pattern)

R/derivoce.R

@@ -542,9 +542,10 @@ derivoce_choices <- function(selected, bathymetry = character(),
   available <- union(selected, as.character(fetchable))
 
   candidate <- function(id, label, group, step, expensive = FALSE,
-                        requires = character()) {
+                        requires = character(), depends = character()) {
     list(id = id, label = label, group = group, expensive = expensive,
-         step = step, requires = setdiff(requires, selected))
+         step = step, requires = setdiff(requires, selected),
+         depends = depends)
   }
 
   # One per fetched covariate. These are the cheap ones, and the ones a habitat
@@ -575,6 +576,16 @@ derivoce_choices <- function(selected, bathymetry = character(),
          step = function(v) list(type = "distance_to_front", var = v))
   )
 
+  # Selected covariates first, then the rest of the catalogue. Deriving from a
+  # covariate you are not modelling is an ordinary thing to want - integrated
+  # chlorophyll without chlorophyll itself, because the accumulated bloom is
+  # what feeds the animals and the instantaneous value is not - and there was
+  # no way to ask for it: the derived covariate was only offered once its
+  # source had been made a predictor, and then there was no way to take the
+  # source back out.
+  #
+  # The extra ones say so in their group rather than their label, so the cost
+  # is visible in the picker without every label carrying a parenthesis.
   out <- list()
   for (entry in per_covariate) {
     for (v in selected) {
@@ -583,6 +594,13 @@ derivoce_choices <- function(selected, bathymetry = character(),
         entry$step(v), entry$expensive
       )
     }
+    for (v in setdiff(available, selected)) {
+      out[[length(out) + 1]] <- candidate(
+        paste0(v, entry$suffix), sprintf(entry$label, v),
+        paste0(entry$group, " (downloads ", v, ")"),
+        entry$step(v), entry$expensive, requires = v
+      )
+    }
   }
 
   # Steps that read particular covariates, offered only when those were fetched.
@@ -602,6 +620,28 @@ derivoce_choices <- function(selected, bathymetry = character(),
       "EKE", "Eddy kinetic energy", "Flow",
       list(type = "eke"), requires = c("UO", "VO")
     )
+    # Steps see what earlier steps produced, so a gradient can be taken of the
+    # speed rather than of the two components it came from. That is the
+    # quantity a front in the flow actually is - the components can each be
+    # changing steeply while the speed is constant, which is a turn and not a
+    # shear - and it is the original pipeline's uv_grad.
+    #
+    # `depends` rather than `requires` because what is needed is another step,
+    # not another download: picking this pulls current_speed in whether or not
+    # the speed itself was asked for as a predictor.
+    for (entry in per_covariate) {
+      if (entry$type == "distance_to_front") next
+      out[[length(out) + 1]] <- candidate(
+        paste0("speed", entry$suffix), sprintf(entry$label, "current speed"),
+        "Flow", entry$step("speed"), entry$expensive,
+        requires = c("UO", "VO"), depends = "speed"
+      )
+    }
+    out[[length(out) + 1]] <- candidate(
+      "EKE_grad", "Spatial gradient of eddy kinetic energy (per km)", "Flow",
+      list(type = "horizontal_gradient", vars = "EKE"),
+      requires = c("UO", "VO"), depends = "EKE"
+    )
     # Backward rather than forward, because backward finds the attracting
     # structures where water converges and plankton accumulate, which is the
     # question a habitat model is asking. Forward finds transport barriers, and
@@ -647,10 +687,40 @@ derivoce_steps_for <- function(ids, selected, bathymetry = character()) {
   if (length(ids) == 0) return(list())
 
   choices <- derivoce_choices(selected, bathymetry)
-  chosen <- Filter(function(x) x$id %in% ids, choices)
+  chosen <- Filter(function(x) x$id %in% with_dependencies(ids, choices),
+                   choices)
   lapply(chosen, function(x) x$step)
 }
 
+#' The chosen derived covariates, plus the ones they are computed from
+#'
+#' A gradient of current speed needs the speed, and the speed is itself a
+#' derived covariate rather than a download. Choosing the gradient therefore
+#' has to pull in the step that produces what it reads — otherwise the config
+#' asks derivoce for a gradient of a column that was never computed.
+#'
+#' Order is not this function's problem. `derivoce_choices()` lists a
+#' dependency before anything that depends on it, and steps are emitted in that
+#' order, so `current_speed` runs before the gradient of `speed` without
+#' anything having to sort them.
+#'
+#' @param ids chosen derived covariate ids
+#' @param choices the [derivoce_choices()] list
+#' @return `ids` with any dependencies added
+#' @keywords internal
+with_dependencies <- function(ids, choices) {
+  needed <- ids
+  repeat {
+    depends <- unlist(lapply(
+      Filter(function(x) x$id %in% needed, choices),
+      function(x) x$depends %||% character()
+    ))
+    grown <- union(needed, depends %||% character())
+    if (length(grown) == length(needed)) return(grown)
+    needed <- grown
+  }
+}
+
 #' Covariates a set of derived choices needs fetching
 #'
 #' A derived covariate is computed from others, and those have to be downloaded
@@ -670,6 +740,47 @@ derivoce_required_inputs <- function(ids, selected, bathymetry = character()) {
   if (length(ids) == 0) return(character())
 
   choices <- derivoce_choices(selected, bathymetry)
-  chosen <- Filter(function(x) x$id %in% ids, choices)
+  chosen <- Filter(function(x) x$id %in% with_dependencies(ids, choices),
+                   choices)
   unique(unlist(lapply(chosen, function(x) x$requires))) %||% character()
 }
+
+#' Derived columns pulled in only as ingredients
+#'
+#' Choosing the gradient of current speed computes the speed on the way, and
+#' that column then sits in the modelling data looking exactly like one that was
+#' asked for. It was not: the ingredient of a derived covariate is no more a
+#' predictor than the velocity components behind an FSLE are.
+#'
+#' This names the derived columns a selection produced without anyone choosing
+#' them, so they can go into `covariates.exclude` alongside the downloads that
+#' [derivoce_required_inputs()] finds. A column named here is still computed and
+#' still available to anything later that reads it; it just does not become a
+#' predictor.
+#'
+#' @param ids chosen derived covariate ids
+#' @param selected time-varying covariate names
+#' @param bathymetry static seafloor covariate names
+#' @return character vector of derived column names, possibly empty
+#' @examples
+#' # Asking for the gradient of current speed computes the speed too, and that
+#' # is an ingredient rather than a request.
+#' derivoce_dependency_columns("speed_grad", c("SST", "SSS"))
+#'
+#' # Asking for both makes the speed a request, so it is not excluded.
+#' derivoce_dependency_columns(c("speed", "speed_grad"), c("SST", "SSS"))
+#' @seealso [derivoce_required_inputs()], which does the same for downloads
+#' @export
+derivoce_dependency_columns <- function(ids, selected,
+                                        bathymetry = character()) {
+  if (length(ids) == 0) return(character())
+
+  choices <- derivoce_choices(selected, bathymetry)
+  pulled <- setdiff(with_dependencies(ids, choices), ids)
+  if (length(pulled) == 0) return(character())
+
+  # The id of a per-covariate candidate is the column it produces, which is
+  # what `exclude` has to name.
+  chosen <- Filter(function(x) x$id %in% pulled, choices)
+  unique(vapply(chosen, function(x) x$id, character(1)))
+}

R/model.R

@@ -35,6 +35,10 @@ fit_patch_model <- function(dat, config) {
   }
 
   model_data <- dat[c(predictors, "patch")] |> as.data.frame()
+  # Kept beside the modelling data rather than in it: the coordinates are not
+  # predictors and must not become any, but spatial_bias() needs to know where
+  # each station was, and `.row` in the held-out predictions indexes this.
+  coordinates <- station_coordinates(dat)
 
   rec <- build_recipe(model_data, config)
   spec <- build_model_spec(config)
@@ -92,6 +96,12 @@ fit_patch_model <- function(dat, config) {
                                  times = bootstrap_times(config),
                                  seed = config$model$seed)
 
+  # Computed before the evaluation table so it can go in it. It is a property
+  # of the folds rather than of the fit, so every model type on these folds
+  # gets the same number - which is what makes it a caveat on a comparison
+  # rather than a score in one.
+  bias <- spatial_bias(list(coordinates = coordinates, predictions = predictions))
+
   # The point estimate stays the model fitted on everything. The ensemble only
   # ever adds columns beside it, so turning uncertainty on never moves a map.
   ensemble <- if (is.null(uncertainty)) {
@@ -103,7 +113,9 @@ fit_patch_model <- function(dat, config) {
   list(
     workflow = fitted,
     metrics = cv_metrics,
-    evaluation = evaluation_table(predictions, cv_metrics, cutoff, bounds),
+    evaluation = evaluation_table(predictions, cv_metrics, cutoff, bounds,
+                                  ssb = overall_ssb(bias)),
+    spatial_bias = bias,
     predictions = predictions,
     classification_threshold = cutoff,
     # The cutoff is estimated, and how much it moves decides whether a binarised
@@ -117,11 +129,25 @@ fit_patch_model <- function(dat, config) {
     # can be produced without refitting. It is the station table, so hundreds to
     # a few thousand rows.
     model_data = model_data,
+    coordinates = coordinates,
     predictors = predictors,
     threshold = attr(dat, "threshold")
   )
 }
 
+#' Station coordinates, in the order the modelling data is in
+#'
+#' `NULL` when the station table has no coordinates, which is the case for a
+#' hand-built frame in a test rather than for anything a run produces.
+#'
+#' @param dat labeled modeling data
+#' @return a two-column data frame of `lon` and `lat`, or `NULL`
+#' @keywords internal
+station_coordinates <- function(dat) {
+  if (!all(c("lon", "lat") %in% names(dat))) return(NULL)
+  as.data.frame(dat[c("lon", "lat")])
+}
+
 #' Assemble a self-explanatory evaluation table
 #'
 #' The cross-validated metrics table reports sensitivity, specificity and kappa
@@ -171,7 +197,8 @@ fit_patch_model <- function(dat, config) {
 #' \doi{10.1111/2041-210X.13140} — the same argument for rare events in species
 #' distribution models, which is what a patch is
 #' @keywords internal
-evaluation_table <- function(predictions, cv_metrics, cutoff, bounds = NULL) {
+evaluation_table <- function(predictions, cv_metrics, cutoff, bounds = NULL,
+                             ssb = NULL) {
   threshold_free <- c("roc_auc", "pr_auc")
 
   ranking <- data.frame(
@@ -195,7 +222,17 @@ evaluation_table <- function(predictions, cv_metrics, cutoff, bounds = NULL) {
   at_best$std_err <- NA_real_
   at_best$note <- "TSS-optimal cutoff; use this one to binarise a projection"
 
-  out <- rbind(ranking, at_default, at_best)
+  # Reported beside the ranking metrics rather than in a file of its own,
+  # because it is the number that says how much to believe them. A reader who
+  # sees roc_auc 0.86 and has to go looking for the caveat will not.
+  bias <- if (is.null(ssb) || is.na(ssb)) NULL else data.frame(
+    metric = "ssb", threshold = NA_real_, value = ssb, std_err = NA_real_,
+    note = paste("spatial sorting bias of the folds, not a model score;",
+                 "1 is fair, near 0 means the metrics above are optimistic"),
+    stringsAsFactors = FALSE
+  )
+
+  out <- rbind(ranking, at_default, at_best, bias)
 
   # Every row gets an interval, including the ones that never had a standard
   # error - which is the point. Column order puts the two uncertainty measures

R/pipeline.R

@@ -117,6 +117,10 @@ run_taupatch <- function(config_path, project = TRUE, keep_covariates = 50000) {
   # 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))
+  # Said out loud rather than left in a column, because it is the caveat on the
+  # line above it and a reader who has to go looking for it will not.
+  ssb <- overall_ssb(model$spatial_bias)
+  if (!is.na(ssb)) message("  ", spatial_bias_note(ssb))
 
   projections <- NULL
   if (project) {
@@ -160,6 +164,9 @@ write_model_outputs <- function(model, config) {
   readr::write_csv(model$evaluation, file.path(out, "evals.csv"))
   readr::write_csv(model$metrics, file.path(out, "cv_metrics.csv"))
   readr::write_csv(model$importance, file.path(out, "var_importance.csv"))
+  if (!is.null(model$spatial_bias)) {
+    readr::write_csv(model$spatial_bias, file.path(out, "spatial_bias.csv"))
+  }
   plot_importance(model$importance, file.path(out, "var_importance.png"))
   write_diagnostic_plots(model, out)
 

R/plot_effects.R

@@ -239,30 +239,6 @@ gam_smooth_terms <- function(fitted) {
   out[order(-out$edf), ]
 }
 
-#' Whether fancyfx is available to draw smooths
-#'
-#' Its own function so the optional path can be exercised in tests without
-#' mocking `requireNamespace()` itself, which every package that loads a
-#' graphics device also goes through.
-#'
-#' @section It used to be called fancygam:
-#' The package was renamed when it grew past GAMs. The rename is why this
-#' matters more than a find-and-replace: `chross22/fancygam` still resolves on
-#' GitHub, so `Remotes: chross22/fancygam` kept installing — but what it
-#' installs now declares `Package: fancyfx`, so `requireNamespace("fancygam")`
-#' returned `FALSE` on every fresh install and the smooth plots were skipped in
-#' silence. Anyone with the old package still sitting in their library saw
-#' nothing wrong.
-#'
-#' That is the failure mode to watch for here: this function gates a diagnostic
-#' rather than the run, so a wrong answer costs a plot and no error.
-#'
-#' @return `TRUE` when fancyfx is installed
-#' @keywords internal
-has_fancyfx <- function() {
-  requireNamespace("fancyfx", quietly = TRUE)
-}
-
 #' Variables a fitted GAM gave a smooth to
 #'
 #' Not every predictor gets one — [model_formula()] gives a linear term to any
@@ -287,8 +263,7 @@ gam_smoothed_variables <- function(fitted) {
 #' rather than reconstructed by prediction, so it carries uncertainty, which a
 #' partial dependence curve cannot.
 #'
-#' Drawn by [fancyfx](https://github.com/chross22/fancyfx), which is a Suggests
-#' — a run without it still gets the generic partial effect curves.
+#' Drawn by [fancyfx](https://github.com/chross22/fancyfx).
 #'
 #' @section Why the axes read in standard deviations:
 #' The smooths belong to the model, and the model was fitted on the recipe's
@@ -309,11 +284,6 @@ gam_smoothed_variables <- function(fitted) {
 #' @seealso [gam_smooth_terms()] for the numbers behind these
 #' @export
 plot_gam_smooths <- function(model, vars = NULL, path = NULL) {
-  if (!has_fancyfx()) {
-    stop("The 'fancyfx' package is required to plot GAM smooths. ",
-         "Install it with remotes::install_github('chross22/fancyfx').",
-         call. = FALSE)
-  }
   workflow <- if (inherits(model, "workflow")) model else model$workflow
 
   vars <- vars %||% gam_smoothed_variables(workflow)
@@ -382,16 +352,12 @@ write_effect_plots <- function(model, out) {
       written <- c(written, "smooth_terms.csv")
     }
 
-    # The fitted smooths themselves, with their uncertainty. Skipped without a
-    # word when fancyfx is absent: it is a Suggests, and the generic partial
-    # effect curves above already cover the question.
-    if (has_fancyfx()) {
-      smooths <- try_diagnostic(
-        plot_gam_smooths(model, path = file.path(out, "gam_smooths.png")),
-        "GAM smooth plots"
-      )
-      if (!is.null(smooths)) written <- c(written, "gam_smooths.png")
-    }
+    # The fitted smooths themselves, with their uncertainty.
+    smooths <- try_diagnostic(
+      plot_gam_smooths(model, path = file.path(out, "gam_smooths.png")),
+      "GAM smooth plots"
+    )
+    if (!is.null(smooths)) written <- c(written, "gam_smooths.png")
   }
 
   invisible(written)