From ee80b6140e38565c0823eadf3541a758a04caafb Mon Sep 17 00:00:00 2001 From: chross22 <52218551+chross22@users.noreply.github.com> Date: Mon, 10 Aug 2026 16:23:20 -0400 Subject: [PATCH 1/5] Plan what taupatch should take from fancyfx, and what it should not fancyfx has gone 0.2.0 to 0.8.0 and now exports 26 functions, a good number of which overlap what taupatch built for itself. This is the survey: what to take, what to leave, what becomes redundant, and the two changes fancyfx needs first. The finding that shapes it is that fancyfx's four evaluation plots re-predict. They route through evaluation_pairs(), which calls predict_probability(model, newdata) and offers no way to hand it predictions that already exist - while taupatch's ROC, threshold, calibration and PR curves come from pooled out-of-fold predictions, each station scored by the fold model that did not see it. Handing fancyfx the final model and its own training stations would produce something weaker, and fancyfx says so itself: it warns, and captions the figure "In-sample: optimistic, and not validation". On a recipes workflow it does not even get that far, since the engine's response is the baked ..y rather than patch, and the call errors. So the four are not a swap as things stand. They become one with a small generalisation upstream - an entry point taking observed and predicted rather than a model - which fancyfx wants anyway, since it is what makes the package usable by anything that cross-validates rather than only by callers holding a single fitted model. Same for mess(): it duplicates novelty_surface() exactly, but takes only a SpatRaster and drops the novel_variable column that says which predictor put a cell outside the training range, which is the half worth having. What survives the "streamlining and improvement" test without any upstream work is the projection panels: plotUncertainty() and plotExtrapolation() take rasters, so none of the above applies, and they delete two hand-built ggplots while adding downsampling a real grid needs. Three are deliberately excluded rather than deferred. thin_points() changes which stations get modelled, which is a study's decision and not a package default. calc_deviance() is defined for glm and gam but not rf or brt, and taupatch reports metrics that mean the same thing across all four on purpose. plotHexbin() solves a problem the station map does not currently have. No code yet. The redundancy pass is deliberately scheduled last, against the code as it stands after the integration rather than against this document. Co-Authored-By: Claude Opus 5 --- docs/fancyfx_plan.md | 220 +++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 220 insertions(+) create mode 100644 docs/fancyfx_plan.md diff --git a/docs/fancyfx_plan.md b/docs/fancyfx_plan.md new file mode 100644 index 0000000..61c971f --- /dev/null +++ b/docs/fancyfx_plan.md @@ -0,0 +1,220 @@ +# Taking what fancyfx now offers + +`fancyfx` has gone from 0.2.0 to 0.8.0 and now exports 26 functions, a good +number of which overlap what taupatch built for itself. This is the plan for +adopting it: what to take, what to leave, what becomes redundant, and the one +thing that cannot be taken as it stands. + +Written against `fancyfx 0.8.0` (`d429991`) and taupatch 0.2.0 (`983cf55`). + +## The finding that shapes everything else + +**fancyfx's four evaluation plots cannot consume taupatch's cross-validated +predictions, and adopting them as they stand would replace honest diagnostics +with optimistic ones.** + +This was worth checking rather than assuming, and the answer is not the one the +function signatures suggest. `plotROC()`, `plotThreshold()`, `plotImportance()` +and `plotCalibration()` all route through `fancyfx:::evaluation_pairs()`, which +does this and only this: + +```r +observed <- binary_response(model, newdata) +predicted <- predict_probability(model, newdata, ...) +``` + +There is no path to hand it predictions that already exist. It re-predicts, +always. + +taupatch's ROC, threshold, calibration and PR curves are drawn from **pooled +out-of-fold predictions**: every station is scored by the one fold model that +never saw it, collected in `model$predictions` with `.row`, `id`, `patch` and +`.pred_patch`. That is cross-validated performance. Handing the final fitted +model and its own training stations to fancyfx instead produces something +different — and fancyfx knows it, which is to its credit: + +* `evaluation_pairs()` warns: *"Evaluating on the data the model was fitted to. + These metrics are optimistic and are not validation."* +* `in_sample_caption()` prints it onto the figure: *"In-sample: scored on the + data the model was fitted to. Optimistic, and not validation."* + +The `folds` argument does not rescue this. It splits the *scoring* of one +already-fitted model by fold; it does not make each fold's predictions come from +a model that held that fold out. On taupatch's data `is_training_data()` returns +`TRUE`, so every plot would carry the caption above. + +A second, blunter obstacle sits on top of the first. taupatch fits through a +`recipes` workflow, so the engine's response is the baked `..y`, not `patch`. +Handing `model_data` to `threshold_metrics()` does not produce an optimistic +number — it errors outright: + +``` +newdata has no column '..y', the model's response. +``` + +So these four cannot be adopted **as fancyfx currently stands**. The fix is not +in taupatch, and it is not a workaround: it is a generalisation fancyfx wants +anyway. + +### Upstream change 1: let evaluation take predictions + +An entry point that accepts predictions rather than a model — the shape any +cross-validated workflow already has, taupatch included: + +```r +threshold_metrics(observed = , predicted = , folds = ) +``` + +`evaluation_pairs()` returns early when handed those, with `in.sample = FALSE`, +and every plot built on it works unchanged. This is the whole change. It costs +fancyfx nothing, it removes the `..y` problem along with the in-sample one, and +it makes fancyfx usable by *any* package that cross-validates rather than only +by ones that keep a single fitted model around. + +Only with this in place do the four evaluation plots become a streamlining: +taupatch deletes four hand-built ggplots and keeps its out-of-fold numbers. +Without it they are a downgrade, and would not be worth taking. + +## What to take now + +### New capability — nothing in taupatch does these + +| function | what it adds | why it belongs here | +|---|---|---| +| `spatial_sorting_bias()` | Hijmans (2012). How much of the AUC is an artefact of presences and absences being differently distributed in space | ECOMON stations are clustered along transects. A spatially sorted sample inflates AUC, and nothing in the current evaluation says by how much. This is a genuine gap | +| `niche_overlap()` | Schoener's *D* and Warren's *I* between two projected surfaces | taupatch models three species and can now compare runs. "Do *cfin* and *ctyp* occupy the same habitat?" and "do rf and gam draw the same map?" are the natural next questions, and `compare_runs()` deliberately answers neither — it compares performance, not surfaces | +| `niche_equivalency()` | A permutation test for whether two niches are the same | The inferential half of the above | + +### Better versions of what exists + +| function | replaces | why it is a clean swap | +|---|---|---| +| `plotUncertainty()` | the spread panel of `plot_projection_uncertainty()` | Takes a raster, not a model, so the in-sample problem does not arise. Adds `cv`/`range`/`iqr` beside `sd`, and `max.cells` downsampling, which the current panel needs on a real grid | +| `plotExtrapolation()` | the novelty panel of the same | Same. Adds `novel.only` | +| `comparePlots()` | nothing yet | Overlays several models' effect curves for one covariate — exactly the per-member view the new ensemble wants | +| `mess()` | `novelty_surface()` | The same method implemented twice. Requires upstream change 2 below | + +### Upstream change 2: let mess() take a data frame, and name the culprit + +`fancyfx::mess()` requires a `SpatRaster` and returns the surface alone: + +``` +Error: x must be a SpatRaster of covariates, not a . +``` + +taupatch's projections are data frames at that point in the pipeline, and its +`novelty_surface()` returns a second column — `novel_variable`, the predictor +that put the cell outside the training range. That column is the actionable +half: "this shelf is extrapolated" is a shrug, "extrapolated because its +chlorophyll is higher than any station saw" is a decision. It is used in the run +log, the projection CSV, and the plot subtitle. + +Two small generalisations, then, and both make `mess()` better on its own terms: +accept a data frame as well as a raster, and return which variable was +responsible. With those, `novelty_surface()` deletes and taupatch calls +`mess()`. + +### Deliberately not taken + +Excluded by the "streamlining **and** improvement" test rather than by any fault +of theirs: + +* **`thin_points()`** — spatial thinning changes which stations are modelled. + That is a methodological decision for a study to make, not a default for a + package to acquire, and it streamlines nothing here. +* **`calc_deviance()`** — defined for `glm` and `gam` and not for `rf` or + `brt`. taupatch reports metrics that mean the same thing for all four types + on purpose; a column that is `NA` for half of them works against that. +* **`plotHexbin()`** — the station map is not currently a problem worth a new + code path. + +## Redundancy assessment + +To be done **after** the integration lands, against the code as it then stands, +not against this plan. The candidates, in the order they are worth examining: + +1. **`novelty_surface()` against `fancyfx::mess()`.** These are the same method + — the multivariate environmental similarity surface of Elith et al. (2010) — + implemented twice. One should go. taupatch's returns `novel_variable`, the + predictor responsible, which is the actionable half and must survive the + merge whichever way it goes. +2. **`plot_projection_uncertainty()`** becomes a thin arrangement of two + fancyfx panels rather than two hand-built ggplots. +3. **`ensemble_spread()` against `fancyfx::ensemble_summary()`.** Same + statistics over a stack of member predictions. taupatch's returns a data + frame aligned to the projection rows; fancyfx's takes a raster. Whether they + can be one function depends on which side the alignment lives. +4. **`permutation_importance()`** — taupatch's is deliberately model-agnostic + and computed by prediction so all four model types are comparable. fancyfx's + takes `(model, newdata)` and has the same in-sample problem as the plots. + Blocked with them, and possibly not worth merging even after. +5. **`plot_importance()`, `plot_calibration()`, `plot_roc_curve()`, + `plot_threshold_performance()`** — not candidates. Their fancyfx + counterparts answer a different question on different evidence. + +The rule for the assessment: **a duplicate is only redundant if the replacement +answers the same question on the same evidence.** `mess()` and +`novelty_surface()` do. `plotROC()` and `plot_roc_curve()` currently do not, +because one is in-sample and the other is not, and that is a difference in what +the number means rather than in how it is drawn. + +## The app + +Into the existing **Diagnostics** tab rather than a new one, since that is where +a user already goes for this material: + +* spatial sorting bias beside the evaluation table, with the one-line + explanation of what a high value means for the AUC above it +* explained deviance in the same block +* the ensemble's per-member effect curves, via `comparePlots()`, when the run + fitted an ensemble +* niche overlap between the active species and one other, when more than one has + been run + +The app must not gain a dependency the package does not have: `fancyfx` is a +Suggests, so every panel here is behind `has_fancyfx()` and degrades to the +current view without it. + +## Order of work + +Two of these are changes to fancyfx. They come first, because the taupatch side +of each is small and pointless without them. + +**In fancyfx** + +1. `threshold_metrics(observed=, predicted=, folds=)` — upstream change 1. +2. `mess()` on a data frame, returning the responsible variable — upstream + change 2. + +**In taupatch** + +3. `plotUncertainty()` / `plotExtrapolation()` into + `plot_projection_uncertainty()`. Needs nothing upstream; two hand-built + ggplots delete. +4. `comparePlots()` for ensemble members. +5. `spatial_sorting_bias()` into the evaluation table. +6. `niche_overlap()` / `niche_equivalency()` as a projection comparison, + alongside `compare_runs()`. +7. `novelty_surface()` replaced by `mess()`, once (2) lands. +8. The four evaluation plots replaced, once (1) lands. +9. App: Diagnostics tab. +10. Redundancy pass, against the code as it then stands. +The four evaluation plots are not on this list, and are not expected to join +it. + +## References + +Elith J, Kearney M, Phillips S (2010). The art of modelling range-shifting +species. *Methods in Ecology and Evolution* **1**(4), 330–342. +doi:10.1111/j.2041-210X.2010.00036.x — MESS + +Hijmans RJ (2012). Cross-validation of species distribution models: removing +spatial sorting bias and calibration with a null model. *Ecology* **93**(3), +679–688. doi:10.1890/11-0826.1 — spatial sorting bias + +Schoener TW (1968). The *Anolis* lizards of Bimini: resource partitioning in a +complex fauna. *Ecology* **49**(4), 704–726. doi:10.2307/1935534 — *D* + +Warren DL, Glor RE, Turelli M (2008). Environmental niche equivalency versus +conservatism. *Evolution* **62**(11), 2868–2883. +doi:10.1111/j.1558-5646.2008.00482.x — *I*, and the equivalency test From 10a2e24decc3fdc9a3755d9c58dce04e01b39d9a Mon Sep 17 00:00:00 2001 From: chross22 <52218551+chross22@users.noreply.github.com> Date: Mon, 10 Aug 2026 16:51:21 -0400 Subject: [PATCH 2/5] Draw the uncertainty panels with fancyfx instead of by hand plot_projection_uncertainty() was ninety lines of ggplot2 building two map panels and a diverging colour ramp, plus a helper whose whole job was to pin that ramp's break to the true zero. fancyfx draws both panels now, and all of it deletes: 101 lines out, 85 in, and the 85 are mostly the documentation of why the panels stay separate. It is a better picture as well as less code. plotUncertainty() and plotExtrapolation() downsample above max.cells, which a real Copernicus grid needs and the hand-built version did not do at all - it handed every cell to geom_raster and hoped. The extrapolation ramp diverges about zero in colours chosen to survive the common colour vision deficiencies, which the hand-rolled six-colour ramp was reaching for and did not achieve. Nothing about what is drawn changes. The panels are still separate rather than blended into one "confidence" layer, for the reason they always were: a cell can be stable across every ensemble member and still be extrapolated, and merging them averages that case away. The novelty subtitle still counts the cells outside the training range and names the covariate responsible, which is taupatch's own and not something fancyfx offers. Both panels take a single-layer raster, and the projection is a table of cells at that point, so projection_raster() changes the container. No resampling: project_patch_model() already builds its GeoTIFF from the same frame. fancyfx is a Suggests, so without it there is a message and no file, exactly as plot_gam_smooths() behaves. That is a real change for anyone who does not have it - the panels used to be drawn unconditionally - and it is the trade the deletion asks for. Worth revisiting whether fancyfx should simply be an Imports, since taupatch already takes datamatch from GitHub that way. Checked by drawing it: a mock ensemble run with projection.uncertainty on writes all three panels - resample spread, algorithm disagreement, novelty - stacked under one title, with the subtitle reporting 1 of 221 cells outside the training range and naming SSS. Co-Authored-By: Claude Opus 5 --- R/plotting.R | 186 +++++++++++++---------------- README.md | 2 +- man/novelty_scale_positions.Rd | 23 ---- man/novelty_subtitle.Rd | 18 +++ man/plot_projection_uncertainty.Rd | 32 +++-- man/projection_raster.Rd | 23 ++++ tests/testthat/test-effects.R | 54 +++++++++ 7 files changed, 202 insertions(+), 136 deletions(-) delete mode 100644 man/novelty_scale_positions.Rd create mode 100644 man/novelty_subtitle.Rd create mode 100644 man/projection_raster.Rd diff --git a/R/plotting.R b/R/plotting.R index 0e3cdf9..879c5f6 100644 --- a/R/plotting.R +++ b/R/plotting.R @@ -23,23 +23,31 @@ plot_projection <- function(predicted, year, month, species, path) { invisible(path) } -#' Plot how far a monthly projection can be trusted +#' Draw the layers that say how far to trust a projection #' -#' The companion to [plot_projection()]: the same cells, panelled by what is -#' uncertain about them rather than by what is predicted. Whichever of the -#' uncertainty surfaces the run produced are drawn, and nothing else. +#' Two or three panels above the map, one per question the projection can be +#' doubted on. Only the ones a run actually produced are drawn. #' -#' The two panels are deliberately not merged into a single "confidence" layer. +#' The panels are deliberately not merged into a single "confidence" layer. #' They measure different things and are free to disagree in either direction, #' and the disagreement is the informative part: a cell can be stable across #' every member and still be extrapolated, because agreement between members -#' trained on the same data is not evidence about ground the data never covered. -#' Blending the two would average that case away instead of showing it. -#' -#' Diverging colour on the novelty panel, centred at zero, because zero is the -#' meaningful break: above it the model is interpolating, below it the cell is -#' outside the training range on some predictor and the model has no evidence -#' for what it says there. +#' trained on the same data is not evidence about ground the data never +#' covered. Blending them would average that case away instead of showing it. +#' +#' @section Drawn by fancyfx: +#' The panels themselves come from [fancyfx::plotUncertainty()] and +#' [fancyfx::plotExtrapolation()] rather than being built here. That is worth a +#' sentence because it used to be ninety lines of `ggplot2` in this file, and +#' the replacement is better in two ways this package would otherwise have had +#' to write for itself: cells are downsampled above `max.cells`, which a real +#' Copernicus grid needs and the hand-built version did not do, and the +#' extrapolation ramp diverges about zero in colours that survive the common +#' colour vision deficiencies. +#' +#' `fancyfx` is a Suggests. Without it there is a message and no file, the same +#' way [plot_gam_smooths()] behaves — the projection itself, and every number +#' behind these panels, is written either way. #' #' @param predicted a projection from `predict_grid()` with uncertainty columns #' @param year year being projected @@ -50,83 +58,49 @@ plot_projection <- function(predicted, year, month, species, path) { #' @seealso [novelty_surface()] for how the novelty panel is computed #' @export plot_projection_uncertainty <- function(predicted, year, month, species, path) { + drawable <- intersect(c("suitability_sd", "algorithm_sd", "novelty"), + names(predicted)) + if (length(drawable) == 0) return(invisible(NULL)) + + if (!has_fancyfx()) { + message(" skipping the uncertainty panels: install fancyfx to draw them ", + "(remotes::install_github('chross22/fancyfx'))") + return(invisible(NULL)) + } + panels <- list() - if ("suitability_sd" %in% names(predicted)) { - panels$spread <- ggplot2::ggplot( - predicted, ggplot2::aes(x = .data$lon, y = .data$lat, - fill = .data$suitability_sd)) + - ggplot2::geom_raster() + - ggplot2::scale_fill_viridis_c(option = "magma", na.value = "white", - name = "SD") + - ggplot2::labs(subtitle = "Spread across the ensemble") + if ("suitability_sd" %in% drawable) { + panels$spread <- fancyfx::plotUncertainty( + projection_raster(predicted, "suitability_sd"), legend.lab = "SD" + ) + 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") + # A multi-algorithm ensemble only. Kept separate from the panel above on + # purpose: that one is a single algorithm refitted on resampled stations, + # this one is different algorithms on the same stations, and a cell where the + # forest and the GLM disagree is not the same worry as a cell where the + # forest is unstable. + if ("algorithm_sd" %in% drawable) { + panels$algorithms <- fancyfx::plotUncertainty( + projection_raster(predicted, "algorithm_sd"), legend.lab = "SD", + option = "cividis" + ) + 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. - # So the visual weight goes to the negatives: saturated red below zero, - # against a muted ramp above it. The alarming colour marks the cells to - # distrust rather than, as a plain diverging scale does, the safest ones. - outside <- sum(!is.na(predicted$novelty) & predicted$novelty < 0) - scored <- sum(!is.na(predicted$novelty)) - - subtitle <- if (outside == 0) { - "Every cell is inside the training range" - } else { - culprit <- names(sort(table( - predicted$novel_variable[!is.na(predicted$novelty) & predicted$novelty < 0] - ), decreasing = TRUE))[1] - sprintf("%d of %d cells (%d%%) are outside the training range, mostly %s", - outside, scored, max(1, round(100 * outside / scored)), culprit) - } - - panels$novelty <- ggplot2::ggplot( - predicted, ggplot2::aes(x = .data$lon, y = .data$lat, - fill = .data$novelty)) + - ggplot2::geom_raster() + - ggplot2::scale_fill_gradientn( - colours = c("#67001f", "#d6604d", "#fddbc7", "#e8eef2", "#7fa8c4", - "#3a6b8a"), - # Zero sits where the warm colours end. rescale puts the break at the - # true zero of this month's range rather than at the midpoint of it, - # so the boundary means the same thing on every map in a run. - values = novelty_scale_positions(predicted$novelty), - na.value = "white", name = "Similarity" - ) + - ggplot2::labs(subtitle = subtitle) + if ("novelty" %in% drawable) { + # Already a MESS surface, so it is passed as one rather than recomputed + # from the covariates - which is what `training = NULL` means here. + panels$novelty <- fancyfx::plotExtrapolation( + projection_raster(predicted, "novelty"), legend.lab = "Similarity" + ) + ggplot2::labs(subtitle = novelty_subtitle(predicted)) } - if (length(panels) == 0) return(invisible(NULL)) - - panels <- lapply(panels, function(p) { - p + ggplot2::coord_quickmap() + - ggplot2::labs(x = NULL, y = NULL) + - ggplot2::theme_bw() + - ggplot2::theme(panel.grid = ggplot2::element_blank(), - legend.position = "bottom") - }) - title <- paste0(species, " - ", month.name[month], " ", year, ": how far to trust it") # Stacked with a shared title rather than side by side, so each panel keeps - # the aspect ratio the coastline has and neither is squashed. patchwork is a + # the aspect ratio the coastline has and none is squashed. patchwork is a # Suggests, and its `/` operator is the whole reason it is wanted here - so # without it, fall back to the panel that carries more of the answer rather # than failing. Novelty is that panel: the spread cannot tell you a cell is @@ -142,37 +116,47 @@ plot_projection_uncertainty <- function(predicted, year, month, species, path) { } ggplot2::ggsave(path, plot = combined, width = 7, - height = if (stacked) 11 else 7, dpi = 150) + height = if (stacked) 4 + 3.5 * length(panels) else 7, + dpi = 150) invisible(path) } -#' Where zero falls on a novelty colour ramp +#' One column of a projection, as a raster fancyfx can draw #' -#' `scale_fill_gradientn()` places its colours at positions in `[0, 1]` across -#' the data range, so zero — the only value on a novelty surface that means -#' anything fixed — lands somewhere different on every map unless it is put -#' there deliberately. This returns positions that pin the warm-to-cool break to -#' the true zero, so red means extrapolated on every month of a run rather than -#' meaning "low for this month". +#' The projection is a table of cells at this point in the pipeline, and every +#' fancyfx map takes a `SpatRaster`. The cells are already on a regular grid — +#' `project_patch_model()` builds the GeoTIFF from the same frame — so this is +#' a change of container rather than any resampling. #' -#' @param novelty the novelty values being plotted -#' @return a vector of six positions in `[0, 1]`, for the six ramp colours +#' @param predicted a projection from `predict_grid()` +#' @param column which column to rasterize +#' @return a single-layer `SpatRaster` #' @keywords internal -novelty_scale_positions <- function(novelty) { - finite <- novelty[is.finite(novelty)] - # Nothing to anchor: an all-positive or empty surface gets an even ramp. - if (length(finite) == 0) return(seq(0, 1, length.out = 6)) - - low <- min(finite) - high <- max(finite) - if (high <= low) return(seq(0, 1, length.out = 6)) - - zero <- (0 - low) / (high - low) - # Clamped off the ends so the break stays visible when everything is on one - # side of zero, which is the common case. - zero <- min(max(zero, 0.02), 0.98) +projection_raster <- function(predicted, column) { + out <- terra::rast(as.data.frame(predicted[c("lon", "lat", column)]), + type = "xyz", crs = "EPSG:4326") + names(out) <- column + out +} - c(0, zero * 0.5, zero * 0.95, zero, zero + (1 - zero) * 0.5, 1) +#' The one line that says whether a novelty panel needed looking at +#' +#' @param predicted a projection carrying `novelty` and `novel_variable` +#' @return a subtitle string +#' @keywords internal +novelty_subtitle <- function(predicted) { + outside <- sum(!is.na(predicted$novelty) & predicted$novelty < 0) + scored <- sum(!is.na(predicted$novelty)) + if (outside == 0) return("Every cell is inside the training range") + + culprit <- if ("novel_variable" %in% names(predicted)) { + names(sort(table( + predicted$novel_variable[!is.na(predicted$novelty) & predicted$novelty < 0] + ), decreasing = TRUE))[1] + } + sprintf("%d of %d cells (%d%%) are outside the training range%s", + outside, scored, max(1, round(100 * outside / scored)), + if (is.null(culprit)) "" else paste(", mostly", culprit)) } #' Build a leaflet map of a projection GeoTIFF diff --git a/README.md b/README.md index 7f89191..08e4536 100644 --- a/README.md +++ b/README.md @@ -1263,7 +1263,7 @@ projections/suitability.csv every cell of every month: species, year, mont projections/suitability.grd the same, as one raster with a layer per month (projection.write_grd) projections/__.tif one layer, or one per surface with projection.uncertainty plots/__.png -plots/___uncertainty.png the spread and novelty panels, with projection.uncertainty +plots/___uncertainty.png the spread and novelty panels, with projection.uncertainty (needs fancyfx) covariates/monthly_means.csv study-area mean per covariate, month, and year covariates/_heatmap.png month-by-year heatmap bathymetry/ marmap's cached NOAA download, if used diff --git a/man/novelty_scale_positions.Rd b/man/novelty_scale_positions.Rd deleted file mode 100644 index 966855b..0000000 --- a/man/novelty_scale_positions.Rd +++ /dev/null @@ -1,23 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/plotting.R -\name{novelty_scale_positions} -\alias{novelty_scale_positions} -\title{Where zero falls on a novelty colour ramp} -\usage{ -novelty_scale_positions(novelty) -} -\arguments{ -\item{novelty}{the novelty values being plotted} -} -\value{ -a vector of six positions in \verb{[0, 1]}, for the six ramp colours -} -\description{ -\code{scale_fill_gradientn()} places its colours at positions in \verb{[0, 1]} across -the data range, so zero — the only value on a novelty surface that means -anything fixed — lands somewhere different on every map unless it is put -there deliberately. This returns positions that pin the warm-to-cool break to -the true zero, so red means extrapolated on every month of a run rather than -meaning "low for this month". -} -\keyword{internal} diff --git a/man/novelty_subtitle.Rd b/man/novelty_subtitle.Rd new file mode 100644 index 0000000..5832b16 --- /dev/null +++ b/man/novelty_subtitle.Rd @@ -0,0 +1,18 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/plotting.R +\name{novelty_subtitle} +\alias{novelty_subtitle} +\title{The one line that says whether a novelty panel needed looking at} +\usage{ +novelty_subtitle(predicted) +} +\arguments{ +\item{predicted}{a projection carrying \code{novelty} and \code{novel_variable}} +} +\value{ +a subtitle string +} +\description{ +The one line that says whether a novelty panel needed looking at +} +\keyword{internal} diff --git a/man/plot_projection_uncertainty.Rd b/man/plot_projection_uncertainty.Rd index 86f7717..c131a86 100644 --- a/man/plot_projection_uncertainty.Rd +++ b/man/plot_projection_uncertainty.Rd @@ -2,7 +2,7 @@ % Please edit documentation in R/plotting.R \name{plot_projection_uncertainty} \alias{plot_projection_uncertainty} -\title{Plot how far a monthly projection can be trusted} +\title{Draw the layers that say how far to trust a projection} \usage{ plot_projection_uncertainty(predicted, year, month, species, path) } @@ -21,23 +21,33 @@ plot_projection_uncertainty(predicted, year, month, species, path) \code{path} invisibly, or \code{NULL} when there is nothing to draw } \description{ -The companion to \code{\link[=plot_projection]{plot_projection()}}: the same cells, panelled by what is -uncertain about them rather than by what is predicted. Whichever of the -uncertainty surfaces the run produced are drawn, and nothing else. +Two or three panels above the map, one per question the projection can be +doubted on. Only the ones a run actually produced are drawn. } \details{ -The two panels are deliberately not merged into a single "confidence" layer. +The panels are deliberately not merged into a single "confidence" layer. They measure different things and are free to disagree in either direction, and the disagreement is the informative part: a cell can be stable across every member and still be extrapolated, because agreement between members -trained on the same data is not evidence about ground the data never covered. -Blending the two would average that case away instead of showing it. +trained on the same data is not evidence about ground the data never +covered. Blending them would average that case away instead of showing it. +} +\section{Drawn by fancyfx}{ + +The panels themselves come from \code{\link[fancyfx:plotUncertainty]{fancyfx::plotUncertainty()}} and +\code{\link[fancyfx:plotExtrapolation]{fancyfx::plotExtrapolation()}} rather than being built here. That is worth a +sentence because it used to be ninety lines of \code{ggplot2} in this file, and +the replacement is better in two ways this package would otherwise have had +to write for itself: cells are downsampled above \code{max.cells}, which a real +Copernicus grid needs and the hand-built version did not do, and the +extrapolation ramp diverges about zero in colours that survive the common +colour vision deficiencies. -Diverging colour on the novelty panel, centred at zero, because zero is the -meaningful break: above it the model is interpolating, below it the cell is -outside the training range on some predictor and the model has no evidence -for what it says there. +\code{fancyfx} is a Suggests. Without it there is a message and no file, the same +way \code{\link[=plot_gam_smooths]{plot_gam_smooths()}} behaves — the projection itself, and every number +behind these panels, is written either way. } + \seealso{ \code{\link[=novelty_surface]{novelty_surface()}} for how the novelty panel is computed } diff --git a/man/projection_raster.Rd b/man/projection_raster.Rd new file mode 100644 index 0000000..24d05e2 --- /dev/null +++ b/man/projection_raster.Rd @@ -0,0 +1,23 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/plotting.R +\name{projection_raster} +\alias{projection_raster} +\title{One column of a projection, as a raster fancyfx can draw} +\usage{ +projection_raster(predicted, column) +} +\arguments{ +\item{predicted}{a projection from \code{predict_grid()}} + +\item{column}{which column to rasterize} +} +\value{ +a single-layer \code{SpatRaster} +} +\description{ +The projection is a table of cells at this point in the pipeline, and every +fancyfx map takes a \code{SpatRaster}. The cells are already on a regular grid — +\code{project_patch_model()} builds the GeoTIFF from the same frame — so this is +a change of container rather than any resampling. +} +\keyword{internal} diff --git a/tests/testthat/test-effects.R b/tests/testthat/test-effects.R index 8df1153..ec53c4a 100644 --- a/tests/testthat/test-effects.R +++ b/tests/testthat/test-effects.R @@ -266,3 +266,57 @@ test_that("abundance over the record is one continuous series", { expect_false("GeomLine" %in% geoms) expect_true(all(geoms %in% c("GeomPoint", "GeomJitter"))) }) + +test_that("the uncertainty panels are drawn for whichever surfaces exist", { + skip_on_cran() + skip_if_not_installed("fancyfx") + cells <- expand.grid(lon = seq(-70, -66, by = 0.5), + lat = seq(41, 44, by = 0.5)) + set.seed(1) + cells$suitability_sd <- runif(nrow(cells), 0, 0.1) + cells$novelty <- runif(nrow(cells), -20, 90) + cells$novel_variable <- "SST" + + path <- tempfile(fileext = ".png") + expect_equal(plot_projection_uncertainty(cells, 2018, 6, "cfin", path), path) + expect_true(file.exists(path)) + + # An algorithm ensemble adds a third panel; a run without one does not get it. + cells$algorithm_sd <- runif(nrow(cells), 0, 0.2) + three <- tempfile(fileext = ".png") + plot_projection_uncertainty(cells, 2018, 6, "cfin", three) + expect_gt(file.info(three)$size, file.info(path)$size) +}) + +test_that("a projection with no uncertainty surfaces draws nothing", { + cells <- data.frame(lon = c(-70, -69), lat = c(41, 42), + suitability = c(0.2, 0.8)) + + expect_null(plot_projection_uncertainty(cells, 2018, 6, "cfin", + tempfile(fileext = ".png"))) +}) + +test_that("without fancyfx the panels are skipped, not failed", { + # fancyfx is a Suggests, and the projection itself plus every number behind + # these panels is written either way. + cells <- data.frame(lon = c(-70, -69), lat = c(41, 42), + novelty = c(10, -5), novel_variable = c("SST", "SST")) + path <- tempfile(fileext = ".png") + + local_mocked_bindings(has_fancyfx = function() FALSE) + expect_message(result <- plot_projection_uncertainty(cells, 2018, 6, "cfin", + path), + "install fancyfx") + expect_null(result) + expect_false(file.exists(path)) +}) + +test_that("the novelty subtitle counts the cells and names the culprit", { + clean <- data.frame(novelty = c(10, 40, 90)) + expect_match(novelty_subtitle(clean), "Every cell is inside") + + mixed <- data.frame(novelty = c(-5, 40, -2, 90), + novel_variable = c("CHL", NA, "CHL", NA)) + expect_match(novelty_subtitle(mixed), "2 of 4 cells") + expect_match(novelty_subtitle(mixed), "mostly CHL") +}) From 84d8454616c260ab7dffab0b7758005724cb9ba5 Mon Sep 17 00:00:00 2001 From: chross22 <52218551+chross22@users.noreply.github.com> Date: Tue, 11 Aug 2026 13:39:55 -0400 Subject: [PATCH 3/5] Make fancyfx an Imports, and let it own the MESS arithmetic Two changes that only make sense together, because the first is what makes the second worth doing. fancyfx moves from Suggests to Imports. It was already in Remotes, so it was already being installed; what the Suggests bought was a set of branches for the case where it was not. Those branches were the cost: has_fancyfx() and its three call sites, an error in plot_gam_smooths(), a silent skip in write_effect_plots(), a message in plot_projection_uncertainty(), a reactive in the app, and two tests mocking an absence that can no longer happen. All of it goes. datamatch is already a hard dependency taken from GitHub the same way, so this is the arrangement the package already had for a dependency it cannot work without - and after the previous commit, the uncertainty panels cannot be drawn without this one. novelty_surface() now calls fancyfx::mess() rather than implementing MESS a second time. variable_similarity(), sixty lines of per-predictor percentile arithmetic, deletes. The two implementations agreed to the last decimal on every case tried, including which predictor gets named as the culprit, which is the argument for there being one of them. What does not change is the interface. novelty_surface() keeps its name, its arguments, and above all its column names: `novelty` and `novel_variable` are what the projection CSV and the GeoTIFF layers are written under, and a run's output should not lose them to an internal tidy-up. mess() returns `mess` and `mess_variable`; the renaming happens here, once. The behaviour is checked rather than assumed: an ordinary cell scores positive, a cell too warm or too cold names SST, one too green names CHL, a missing covariate gives NA, a grid with no shared predictors gives all NA, and a single-row grid keeps its shape - the case a vapply collapse used to get wrong. A cell where one predictor is NA and another is not still names the one that is worst among those it could score, exactly as before. The tests for variable_similarity() go with it. They tested arithmetic this package no longer owns, and fancyfx tests it. What replaces them tests what taupatch still asserts to its readers: that the scale reads 100 at the median, falls to zero at the edge of the training range, and goes negative outside it in proportion to how far. Co-Authored-By: Claude Opus 5 --- DESCRIPTION | 2 +- R/plot_effects.R | 48 +++---------------- R/plotting.R | 10 ---- R/uncertainty.R | 77 +++++++----------------------- README.md | 19 ++++---- inst/shiny/app.R | 8 ++-- man/has_fancyfx.Rd | 31 ------------ man/novelty_surface.Rd | 10 ++++ man/plot_gam_smooths.Rd | 3 +- man/plot_projection_uncertainty.Rd | 4 -- man/variable_similarity.Rd | 22 --------- tests/testthat/test-effects.R | 33 ------------- tests/testthat/test-uncertainty.R | 30 +++++------- 13 files changed, 62 insertions(+), 235 deletions(-) delete mode 100644 man/has_fancyfx.Rd delete mode 100644 man/variable_similarity.Rd diff --git a/DESCRIPTION b/DESCRIPTION index 1bd0629..f6b6a93 100644 --- a/DESCRIPTION +++ b/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, diff --git a/R/plot_effects.R b/R/plot_effects.R index 0997a89..e27ba3d 100644 --- a/R/plot_effects.R +++ b/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) diff --git a/R/plotting.R b/R/plotting.R index 879c5f6..a239eac 100644 --- a/R/plotting.R +++ b/R/plotting.R @@ -45,10 +45,6 @@ plot_projection <- function(predicted, year, month, species, path) { #' extrapolation ramp diverges about zero in colours that survive the common #' colour vision deficiencies. #' -#' `fancyfx` is a Suggests. Without it there is a message and no file, the same -#' way [plot_gam_smooths()] behaves — the projection itself, and every number -#' behind these panels, is written either way. -#' #' @param predicted a projection from `predict_grid()` with uncertainty columns #' @param year year being projected #' @param month month being projected @@ -62,12 +58,6 @@ plot_projection_uncertainty <- function(predicted, year, month, species, path) { names(predicted)) if (length(drawable) == 0) return(invisible(NULL)) - if (!has_fancyfx()) { - message(" skipping the uncertainty panels: install fancyfx to draw them ", - "(remotes::install_github('chross22/fancyfx'))") - return(invisible(NULL)) - } - panels <- list() if ("suitability_sd" %in% drawable) { diff --git a/R/uncertainty.R b/R/uncertainty.R index b618789..e75e35a 100644 --- a/R/uncertainty.R +++ b/R/uncertainty.R @@ -231,6 +231,14 @@ ensemble_spread <- function(ensemble, newdata, level = 0.9) { #' because its chlorophyll is higher than anything a station saw" is a decision #' about whether to widen the training window or clip the map. #' +#' @section Where the numbers come from: +#' [fancyfx::mess()] computes both columns. taupatch used to implement MESS +#' itself, and the two implementations agreed to the last decimal on every case +#' tested -- which is the argument for there being one of them rather than two. +#' What survives here is the interface: the `novelty` and `novel_variable` +#' names, which the projection CSV and the GeoTIFF layers are written under and +#' which a run\'s output should not lose to an internal tidy-up. +#' #' @param grid the cells to score, with one column per predictor #' @param model_data the data the model was fitted on #' @param predictors predictor column names @@ -256,67 +264,16 @@ novelty_surface <- function(grid, model_data, predictors) { stringsAsFactors = FALSE)) } - similarity <- vapply(predictors, function(v) { - variable_similarity(grid[[v]], model_data[[v]]) - }, numeric(nrow(grid))) - # vapply drops to a vector when the grid has one row, which then indexes - # wrongly below. - similarity <- matrix(similarity, nrow = nrow(grid), - dimnames = list(NULL, predictors)) - - worst <- apply(similarity, 1, function(row) { - if (all(is.na(row))) return(NA_integer_) - which.min(row) - }) - - data.frame( - novelty = apply(similarity, 1, min, na.rm = FALSE), - novel_variable = ifelse(is.na(worst), NA_character_, predictors[worst]), - stringsAsFactors = FALSE - ) -} - -#' Similarity of values to a training distribution, for one predictor -#' -#' The per-variable half of [novelty_surface()]. Negative below the training -#' minimum and above its maximum, scaled by the training range; inside, twice -#' the distance to the nearer tail in percentile terms, so the median scores 100. -#' -#' @param values the values to score -#' @param train the training values for the same predictor -#' @return a numeric vector the length of `values` -#' @keywords internal -variable_similarity <- function(values, train) { - train <- sort(train[is.finite(train)]) - n <- length(train) - if (n == 0) return(rep(NA_real_, length(values))) - - low <- train[1] - high <- train[n] - span <- high - low - - # A predictor that never varied in training carries no information about - # what is ordinary, so it can only say same-or-not. - if (span == 0) { - return(ifelse(is.na(values), NA_real_, ifelse(values == low, 100, -Inf))) - } - - # Percentile of each value within the training distribution. findInterval - # with left.open counts the training values strictly below, which is the - # count this needs and is a binary search rather than a full comparison - # against every training row - the difference between a grid that scores in - # a moment and one that does not. - below <- findInterval(values, train, left.open = TRUE) - percentile <- 100 * below / n + # The arithmetic is fancyfx's. This function is the name taupatch's output + # columns are written under, and the place the scale is explained; keeping it + # as a thin call means the projection CSV and the GeoTIFF layer keep the + # names they have always had. + out <- fancyfx::mess(as.data.frame(grid[predictors]), + as.data.frame(model_data[predictors]), + vars = predictors, limiting = TRUE) - out <- ifelse( - percentile == 0, 100 * (values - low) / span, - ifelse(percentile <= 50, 2 * percentile, - ifelse(percentile < 100, 2 * (100 - percentile), - 100 * (high - values) / span)) - ) - out[is.na(values)] <- NA_real_ - out + data.frame(novelty = out$mess, novel_variable = out$mess_variable, + stringsAsFactors = FALSE) } #' Summarise a novelty surface for the run log diff --git a/README.md b/README.md index 08e4536..67bb115 100644 --- a/README.md +++ b/README.md @@ -145,6 +145,9 @@ remotes::install_github("chross22/datamatch") `datamatch` no longer depends on `BigelowLab/copernicus` — it calls the Copernicus Marine Toolbox directly — so that package no longer needs installing. +`fancyfx` draws the effect and uncertainty figures and is installed with the +package, from the `Remotes` field, so it needs no separate step. + Derived covariates (`covariates.derivoce`) additionally need: ```r @@ -878,14 +881,10 @@ On top of that, each model contributes what only it can: | `gam` | Effective degrees of freedom per smooth. An `edf` of 1 means the smooth collapsed to a line | | `rf` / `brt` | None. The partial effect curve *is* their answer | -With [`fancyfx`](https://github.com/chross22/fancyfx) installed, a GAM also -gets its **fitted smooths** drawn — each term with its standard error band and a -rug showing where the data actually is. Those carry uncertainty, which a partial -dependence curve cannot: - -```r -remotes::install_github("chross22/fancyfx") -``` +A GAM also gets its **fitted smooths** drawn by +[`fancyfx`](https://github.com/chross22/fancyfx) — each term with its standard +error band and a rug showing where the data actually is. Those carry +uncertainty, which a partial dependence curve cannot. Their x axes read in standard deviations, because the smooths belong to the model and the model was fitted on the recipe's output. Set `covariates.normalize: false` @@ -1256,14 +1255,14 @@ diagnostics/cv_predictions.csv held-out predictions, for any metric not tabu 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 fancyfx: fitted smooths with error bands +diagnostics/gam_smooths.png gam only: fitted smooths with error bands, drawn by fancyfx 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) projections/__.tif one layer, or one per surface with projection.uncertainty plots/__.png -plots/___uncertainty.png the spread and novelty panels, with projection.uncertainty (needs fancyfx) +plots/___uncertainty.png the spread and novelty panels, with projection.uncertainty covariates/monthly_means.csv study-area mean per covariate, month, and year covariates/_heatmap.png month-by-year heatmap bathymetry/ marmap's cached NOAA download, if used diff --git a/inst/shiny/app.R b/inst/shiny/app.R index ac3e084..70956ad 100644 --- a/inst/shiny/app.R +++ b/inst/shiny/app.R @@ -1210,9 +1210,9 @@ server <- function(input, output, session) { # A GAM has its own partial effects, and they are better than the generic # ones: read out of the fitted model rather than reconstructed by prediction, # so they carry the uncertainty a partial dependence curve cannot. + # fancyfx is an Imports now, so this is only ever about the model type. use_fancyfx <- reactive({ - identical(run_result()$model$type, "gam") && - requireNamespace("fancyfx", quietly = TRUE) + identical(run_result()$model$type, "gam") }) output$partial_effects_note <- renderUI({ @@ -1256,8 +1256,8 @@ server <- function(input, output, session) { )) } if (identical(model$type, "gam")) { - # The smooths themselves are the partial effects panel above when - # fancyfx is present, so they are not repeated here. + # The smooths themselves are the partial effects panel above, so they + # are not repeated here. return(tagList( h4("Smooth terms"), helpText("Effective degrees of freedom per smooth. An edf of 1 means", diff --git a/man/has_fancyfx.Rd b/man/has_fancyfx.Rd deleted file mode 100644 index a35d08d..0000000 --- a/man/has_fancyfx.Rd +++ /dev/null @@ -1,31 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/plot_effects.R -\name{has_fancyfx} -\alias{has_fancyfx} -\title{Whether fancyfx is available to draw smooths} -\usage{ -has_fancyfx() -} -\value{ -\code{TRUE} when fancyfx is installed -} -\description{ -Its own function so the optional path can be exercised in tests without -mocking \code{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: \code{chross22/fancygam} still resolves on -GitHub, so \code{Remotes: chross22/fancygam} kept installing — but what it -installs now declares \code{Package: fancyfx}, so \code{requireNamespace("fancygam")} -returned \code{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. -} - -\keyword{internal} diff --git a/man/novelty_surface.Rd b/man/novelty_surface.Rd index 9295ed5..4c076a1 100644 --- a/man/novelty_surface.Rd +++ b/man/novelty_surface.Rd @@ -40,6 +40,16 @@ half: "this coast is extrapolated" is a shrug, and "this coast is extrapolated because its chlorophyll is higher than anything a station saw" is a decision about whether to widen the training window or clip the map. } +\section{Where the numbers come from}{ + +\code{\link[fancyfx:mess]{fancyfx::mess()}} computes both columns. taupatch used to implement MESS +itself, and the two implementations agreed to the last decimal on every case +tested -- which is the argument for there being one of them rather than two. +What survives here is the interface: the \code{novelty} and \code{novel_variable} +names, which the projection CSV and the GeoTIFF layers are written under and +which a run\'s output should not lose to an internal tidy-up. +} + \examples{ train <- data.frame(SST = c(4, 8, 12, 16), CHL = c(0.2, 0.5, 1.0, 2.0)) grid <- data.frame(SST = c(10, 25), CHL = c(0.6, 0.6)) diff --git a/man/plot_gam_smooths.Rd b/man/plot_gam_smooths.Rd index e0831bb..30cd0ee 100644 --- a/man/plot_gam_smooths.Rd +++ b/man/plot_gam_smooths.Rd @@ -24,8 +24,7 @@ rather than reconstructed by prediction, so it carries uncertainty, which a partial dependence curve cannot. } \details{ -Drawn by \href{https://github.com/chross22/fancyfx}{fancyfx}, which is a Suggests -— a run without it still gets the generic partial effect curves. +Drawn by \href{https://github.com/chross22/fancyfx}{fancyfx}. } \section{Why the axes read in standard deviations}{ diff --git a/man/plot_projection_uncertainty.Rd b/man/plot_projection_uncertainty.Rd index c131a86..163c8f4 100644 --- a/man/plot_projection_uncertainty.Rd +++ b/man/plot_projection_uncertainty.Rd @@ -42,10 +42,6 @@ to write for itself: cells are downsampled above \code{max.cells}, which a real Copernicus grid needs and the hand-built version did not do, and the extrapolation ramp diverges about zero in colours that survive the common colour vision deficiencies. - -\code{fancyfx} is a Suggests. Without it there is a message and no file, the same -way \code{\link[=plot_gam_smooths]{plot_gam_smooths()}} behaves — the projection itself, and every number -behind these panels, is written either way. } \seealso{ diff --git a/man/variable_similarity.Rd b/man/variable_similarity.Rd deleted file mode 100644 index 6728270..0000000 --- a/man/variable_similarity.Rd +++ /dev/null @@ -1,22 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/uncertainty.R -\name{variable_similarity} -\alias{variable_similarity} -\title{Similarity of values to a training distribution, for one predictor} -\usage{ -variable_similarity(values, train) -} -\arguments{ -\item{values}{the values to score} - -\item{train}{the training values for the same predictor} -} -\value{ -a numeric vector the length of \code{values} -} -\description{ -The per-variable half of \code{\link[=novelty_surface]{novelty_surface()}}. Negative below the training -minimum and above its maximum, scaled by the training range; inside, twice -the distance to the nearer tail in percentile terms, so the median scores 100. -} -\keyword{internal} diff --git a/tests/testthat/test-effects.R b/tests/testthat/test-effects.R index ec53c4a..ee3342e 100644 --- a/tests/testthat/test-effects.R +++ b/tests/testthat/test-effects.R @@ -176,7 +176,6 @@ test_that("the smoothed variables are read back off the fitted model", { test_that("fancyfx draws the fitted smooths for a GAM", { skip_on_cran() skip_if_not_installed("mgcv") - skip_if_not_installed("fancyfx") model <- fitted_of_type("gam") plot <- suppressMessages(plot_gam_smooths(model)) @@ -194,7 +193,6 @@ test_that("fancyfx draws the fitted smooths for a GAM", { test_that("smooth plots are drawn against the data the model actually saw", { skip_on_cran() skip_if_not_installed("mgcv") - skip_if_not_installed("fancyfx") model <- fitted_of_type("gam") # gratia reports the smooths in the recipe's output units, so the rug has to @@ -207,21 +205,6 @@ test_that("smooth plots are drawn against the data the model actually saw", { expect_gt(abs(mean(model$model_data$SST)), 1) }) -test_that("without fancyfx the run still gets its generic curves", { - skip_on_cran() - skip_if_not_installed("mgcv") - # fancyfx is a Suggests: its absence removes the extra plot, not the - # diagnostics. - local_mocked_bindings(has_fancyfx = function() FALSE) - out <- tempfile("diag"); dir.create(out) - suppressMessages(suppressWarnings(write_effect_plots(fitted_of_type("gam"), out))) - written <- list.files(out) - - expect_false("gam_smooths.png" %in% written) - expect_true("partial_effects.png" %in% written) - expect_true("smooth_terms.csv" %in% written) -}) - test_that("maps draw, not just build", { skip_on_cran() # geom_sf() refuses to work under any coord but coord_sf(), and refuses when @@ -269,7 +252,6 @@ test_that("abundance over the record is one continuous series", { test_that("the uncertainty panels are drawn for whichever surfaces exist", { skip_on_cran() - skip_if_not_installed("fancyfx") cells <- expand.grid(lon = seq(-70, -66, by = 0.5), lat = seq(41, 44, by = 0.5)) set.seed(1) @@ -296,21 +278,6 @@ test_that("a projection with no uncertainty surfaces draws nothing", { tempfile(fileext = ".png"))) }) -test_that("without fancyfx the panels are skipped, not failed", { - # fancyfx is a Suggests, and the projection itself plus every number behind - # these panels is written either way. - cells <- data.frame(lon = c(-70, -69), lat = c(41, 42), - novelty = c(10, -5), novel_variable = c("SST", "SST")) - path <- tempfile(fileext = ".png") - - local_mocked_bindings(has_fancyfx = function() FALSE) - expect_message(result <- plot_projection_uncertainty(cells, 2018, 6, "cfin", - path), - "install fancyfx") - expect_null(result) - expect_false(file.exists(path)) -}) - test_that("the novelty subtitle counts the cells and names the culprit", { clean <- data.frame(novelty = c(10, 40, 90)) expect_match(novelty_subtitle(clean), "Every cell is inside") diff --git a/tests/testthat/test-uncertainty.R b/tests/testthat/test-uncertainty.R index 19cb0f8..bcf11f1 100644 --- a/tests/testthat/test-uncertainty.R +++ b/tests/testthat/test-uncertainty.R @@ -42,26 +42,22 @@ test_that("malformed uncertainty settings are refused at config load", { # ---- novelty ---------------------------------------------------------------- -test_that("similarity is negative outside the training range and scaled by it", { - train <- c(0, 10) # range of 10 +test_that("the MESS scale still reads the way the documentation says", { + # The arithmetic belongs to fancyfx::mess() now, and is tested there. What is + # checked here is the property this package documents and its readers rely + # on: 100 at the median, falling to 0 at the edge of the training range, and + # negative outside it in proportion to how far. + train <- data.frame(x = 0:100) - # Half a range below the minimum, and one range above the maximum. - expect_equal(unname(variable_similarity(-5, train)), -50) - expect_equal(unname(variable_similarity(20, train)), -100) + middle <- novelty_surface(data.frame(x = 50), train, "x")$novelty + edge <- novelty_surface(data.frame(x = 99), train, "x")$novelty + outside <- novelty_surface(data.frame(x = c(-50, 150)), train, "x")$novelty - # Inside, it is positive. The endpoints are the edge of the range, not - # outside it, so they are not negative. - expect_gte(variable_similarity(5, train), 0) -}) - -test_that("similarity peaks at the middle of the training data", { - train <- 1:101 - - middle <- variable_similarity(51, train) - edge <- variable_similarity(95, train) - - expect_gt(middle, edge) expect_lte(middle, 100) + expect_gt(middle, edge) + expect_gte(edge, 0) + # Half a training range below the minimum, and half a range above the max. + expect_equal(outside, c(-50, -50)) }) test_that("a cell is as novel as its worst predictor, and says which", { From 98d0dbde13d8ac48abee079b79453aabd2f15e23 Mon Sep 17 00:00:00 2001 From: chross22 <52218551+chross22@users.noreply.github.com> Date: Tue, 11 Aug 2026 14:40:29 -0400 Subject: [PATCH 4/5] Report the spatial sorting bias, and whether two maps are the same map Two things taupatch could not say before, both from fancyfx. spatial_bias() is the spatial sorting bias of Hijmans (2012): how much closer the held-out patches sit to the training patches than the held-out non-patches do. It matters here more than in most places. ECOMON stations are not scattered at random - they sit on transects, revisited season after season - so a randomly held-out station usually has a near neighbour in the training folds, often the same station in another year. A model can then score well by recognising where it has already been rather than what makes a patch, and nothing in the ROC curve tells the two apart. It goes in evals.csv as an `ssb` row rather than into a file of its own, because it is the caveat on the roc_auc three rows above it and a reader who has to go looking for it will not. The run says it out loud too. The per-fold table is written to spatial_bias.csv, so one odd fold is visible rather than averaged into the answer. On the mock config it comes out at 0.2, which says that run's 0.85 AUC is optimistic - a real finding the evaluation could not previously produce, and an argument for spatially blocked folds rather than for a different model. This needed the fitted model to know where its stations were. model_data deliberately holds predictors only, and coordinates must not become predictors, so they are carried beside it and indexed by the same .row the held-out predictions use. projection_overlap() answers what compare_runs() cannot: not which surface scores better, but whether the two are the same surface. Schoener's D and Warren's I. Two runs can rank stations equally well and disagree completely about where the habitat is, and a comparison of AUCs reports them as equivalent. Cells are matched on coordinates rather than row order, and projections covering different ground are intersected with a warning, the same contract compare_runs() uses. Both statistics run high, and that is documented rather than left to be misread: two surfaces of independent noise over the same grid score about 0.7 on D and 0.9 on I, because both spread their probability over the same cells in similar proportions. Those are the floor for surfaces of this shape, not near-identity. niche_equivalency() is not wired up. It wants occurrence points for two entities and a refitting harness, which does not map onto a config-driven pipeline without inventing an interface for it, so it stays unused until there is a reason. Co-Authored-By: Claude Opus 5 --- NAMESPACE | 2 + R/model.R | 43 +++++++++- R/pipeline.R | 7 ++ R/power.R | 114 +++++++++++++++++++++++++ R/spatial_bias.R | 133 ++++++++++++++++++++++++++++++ README.md | 13 +++ man/evaluation_table.Rd | 2 +- man/overall_ssb.Rd | 18 ++++ man/projection_cells.Rd | 20 +++++ man/projection_overlap.Rd | 73 ++++++++++++++++ man/spatial_bias.Rd | 70 ++++++++++++++++ man/spatial_bias_note.Rd | 21 +++++ man/station_coordinates.Rd | 19 +++++ tests/testthat/test-power.R | 74 +++++++++++++++++ tests/testthat/test-uncertainty.R | 85 +++++++++++++++++++ tools/citations.csv | 3 + 16 files changed, 693 insertions(+), 4 deletions(-) create mode 100644 R/spatial_bias.R create mode 100644 man/overall_ssb.Rd create mode 100644 man/projection_cells.Rd create mode 100644 man/projection_overlap.Rd create mode 100644 man/spatial_bias.Rd create mode 100644 man/spatial_bias_note.Rd create mode 100644 man/station_coordinates.Rd diff --git a/NAMESPACE b/NAMESPACE index c903863..f484643 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -69,11 +69,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) diff --git a/R/model.R b/R/model.R index 49f1860..bc98227 100644 --- a/R/model.R +++ b/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 diff --git a/R/pipeline.R b/R/pipeline.R index 0241f8e..c9b384c 100644 --- a/R/pipeline.R +++ b/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) diff --git a/R/power.R b/R/power.R index 711a250..fae105d 100644 --- a/R/power.R +++ b/R/power.R @@ -485,3 +485,117 @@ power_point <- function(model_data, configs, predictors, types, fraction, stringsAsFactors = FALSE ) } + +#' How much two projected surfaces agree about the habitat +#' +#' [compare_runs()] asks which run scores better. This asks something the +#' scores cannot: whether the two maps are the *same map*. Two runs can rank +#' stations equally well and disagree completely about where the habitat is, +#' and a comparison of AUCs would report them as equivalent. +#' +#' Schoener's *D* and Warren's *I*, both on 0 to 1. **1** is identical +#' surfaces; **0** is no overlap at all. There is no threshold at which two +#' niches become "the same" — the number is a description, and what counts as +#' close depends on what the maps are for. +#' +#' @section Read them as high numbers: +#' Both statistics compare the suitability *distributions* cell by cell after +#' normalising them, so they run high. Two surfaces of independent random noise +#' over the same grid score around 0.7 on *D* and 0.9 on *I* — not because they +#' agree about anything, but because both spread their probability over the +#' same cells in similar proportions. Treat those as the floor for surfaces of +#' this shape rather than reading 0.9 as near-identity, and compare overlaps +#' against each other rather than against 1. +#' +#' @section What it is good for: +#' Two questions, mostly. *Do two species occupy the same habitat?* — run +#' `cfin` and `ctyp` over the same months and compare their surfaces. And *does +#' the choice of algorithm change the map?* — the more useful one when a +#' comparison of scores came back inconclusive, because a small difference in +#' AUC with a low overlap means the two models are doing genuinely different +#' things and the tie in performance is hiding it. +#' +#' @section Matching the cells: +#' The two surfaces must be the same cells in the same order, or the overlap is +#' measured between places that have nothing to do with each other. Cells are +#' matched on their coordinates rather than on row position, and cells present +#' in one projection and not the other are dropped with a warning — a +#' projection that covers less ground usually has a reason, and averaging over +#' the difference would hide it. +#' +#' @param x,y projections to compare: data frames with `lon`, `lat` and +#' `suitability`, as [project_patch_model()] writes to `suitability.csv` +#' @param statistic `"both"`, `"D"`, or `"I"` +#' @param digits how many decimal places the coordinates are matched on. +#' Projections from the same grid agree exactly; this is for two that came +#' off slightly different pipelines +#' @return a named numeric vector, with an `n_cells` attribute +#' @examples +#' \dontrun{ +#' rf <- readr::read_csv("output/rf/projections/suitability.csv") +#' gam <- readr::read_csv("output/gam/projections/suitability.csv") +#' +#' projection_overlap(rf[rf$month == 7, ], gam[gam$month == 7, ]) +#' } +#' @references +#' Schoener TW (1968). The *Anolis* lizards of Bimini: resource partitioning in +#' a complex fauna. *Ecology* **49**(4), 704-726. \doi{10.2307/1935534} — *D* +#' +#' Warren DL, Glor RE, Turelli M (2008). Environmental niche equivalency versus +#' conservatism: quantitative approaches to niche evolution. *Evolution* +#' **62**(11), 2868-2883. \doi{10.1111/j.1558-5646.2008.00482.x} — *I* +#' @seealso [compare_runs()], which asks which surface is *better* rather than +#' whether they are the same +#' @export +projection_overlap <- function(x, y, statistic = c("both", "D", "I"), + digits = 6) { + statistic <- match.arg(statistic) + x <- projection_cells(x, "x") + y <- projection_cells(y, "y") + + key <- function(cells) { + paste(round(cells$lon, digits), round(cells$lat, digits), sep = "\r") + } + shared <- intersect(key(x), key(y)) + dropped <- length(union(key(x), key(y))) - length(shared) + + if (length(shared) < 2) { + stop("The two projections share ", length(shared), " cells, so there is ", + "no common ground to measure overlap over.\nThey were probably ", + "projected onto different grids; `covariates.grid` decides that.", + call. = FALSE) + } + if (dropped > 0) { + warning("The two projections do not cover the same cells: ", dropped, + " of ", length(shared) + dropped, " are in one and not the other, ", + "and the overlap uses the ", length(shared), " they share.", + call. = FALSE) + } + + out <- fancyfx::niche_overlap( + x$suitability[match(shared, key(x))], + y$suitability[match(shared, key(y))], + statistic = statistic + ) + attr(out, "n_cells") <- length(shared) + out +} + +#' The cells of a projection, checked +#' +#' @param projection a projection table +#' @param label which argument it came from, for the error +#' @return a data frame of `lon`, `lat` and `suitability` +#' @keywords internal +projection_cells <- function(projection, label) { + projection <- as.data.frame(projection) + needed <- c("lon", "lat", "suitability") + missing <- setdiff(needed, names(projection)) + if (length(missing) > 0) { + stop("`", label, "` is missing: ", paste(missing, collapse = ", "), + ".\nprojection_overlap() reads the table project_patch_model() ", + "writes to suitability.csv; one month of it at a time, since two ", + "months stacked are two surfaces rather than one.", call. = FALSE) + } + projection[needed] +} diff --git a/R/spatial_bias.R b/R/spatial_bias.R new file mode 100644 index 0000000..b76611c --- /dev/null +++ b/R/spatial_bias.R @@ -0,0 +1,133 @@ +#' How much of the AUC is an artefact of where the stations are +#' +#' A cross-validated AUC is optimistic when the held-out patches sit closer to +#' the training patches than the held-out non-patches do. The model can then +#' score well by recognising *where* it has already been rather than *what* +#' makes a patch, and nothing in the ROC curve distinguishes the two. +#' +#' This is the spatial sorting bias of Hijmans (2012): the mean distance from +#' each held-out patch to the nearest training patch, over the same distance +#' for the held-out non-patches. +#' +#' \deqn{SSB = \frac{\overline{d}(\mathrm{held\ out\ patch},\ \mathrm{training\ patch})}{\overline{d}(\mathrm{held\ out\ non\ patch},\ \mathrm{training\ patch})}} +#' +#' * **Near 1** — held-out patches and non-patches are equally far from the +#' training patches. The split is spatially fair and the AUC means what it +#' appears to. +#' * **Near 0** — held-out patches are much closer to training patches than the +#' non-patches are. **The AUC is inflated**, and by an amount this number +#' does not tell you. +#' +#' It is a property of the *split*, not of the model. Every model type fitted on +#' the same folds gets the same value, which is the point: it says how much to +#' trust the comparison between them, not which of them won. +#' +#' @section Why an ECOMON run should look at it: +#' Survey stations are not scattered at random. They sit on transects, revisited +#' season after season, so a randomly held-out station usually has a near +#' neighbour in the training folds — often the same station in another year. +#' That is the situation this measures, and random `model.cv_folds` cannot avoid +#' it. A low value is an argument for spatially blocked folds, not for a +#' different model. +#' +#' @param model a fitted model from [fit_patch_model()] or +#' [fit_patch_ensemble()], carrying `coordinates` and held-out `predictions` +#' @param geo whether to measure great-circle distances. `TRUE` by default, +#' because the coordinates are longitude and latitude, and a degree of +#' longitude is not a degree of latitude anywhere but the equator +#' @return a data frame with one row per fold and a `fold` of `"overall"` for +#' the mean, carrying `patch_distance`, `non_patch_distance` and `ssb` +#' @examples +#' \dontrun{ +#' model <- fit_patch_model(dat, config) +#' spatial_bias(model) +#' } +#' @references +#' Hijmans RJ (2012). Cross-validation of species distribution models: removing +#' spatial sorting bias and calibration with a null model. *Ecology* **93**(3), +#' 679-688. \doi{10.1890/11-0826.1} +#' @seealso [evaluation_table()], which reports the summary beside the AUC it +#' qualifies +#' @export +spatial_bias <- function(model, geo = TRUE) { + coordinates <- model$coordinates + predictions <- model$predictions + + if (is.null(coordinates) || is.null(predictions) || + !all(c(".row", "id", "patch") %in% names(predictions))) { + return(NULL) + } + + folds <- split(seq_len(nrow(predictions)), predictions$id) + rows <- lapply(names(folds), function(fold) { + at <- folds[[fold]] + held_out <- predictions$.row[at] + is_patch <- predictions$patch[at] == "patch" + + # The reference is the training patches: every patch station except the + # ones held out in this fold. A held-out patch that is close to one of them + # is a patch the model has, in effect, already seen. + all_patches <- predictions$.row[predictions$patch == "patch"] + training <- setdiff(all_patches, held_out[is_patch]) + + if (!any(is_patch) || !any(!is_patch) || length(training) == 0) { + return(NULL) + } + + out <- fancyfx::spatial_sorting_bias( + presence = coordinates[held_out[is_patch], , drop = FALSE], + absence = coordinates[held_out[!is_patch], , drop = FALSE], + reference = coordinates[training, , drop = FALSE], + geo = geo + ) + data.frame(fold = fold, patch_distance = out[["presence"]], + non_patch_distance = out[["absence"]], ssb = out[["ssb"]], + stringsAsFactors = FALSE) + }) + + out <- do.call(rbind, Filter(Negate(is.null), rows)) + if (is.null(out) || nrow(out) == 0) return(NULL) + + # One number to report beside the AUC, and the folds it came from kept so a + # single odd fold is visible rather than averaged into the answer. + overall <- data.frame( + fold = "overall", + patch_distance = mean(out$patch_distance, na.rm = TRUE), + non_patch_distance = mean(out$non_patch_distance, na.rm = TRUE), + ssb = mean(out$ssb, na.rm = TRUE), + stringsAsFactors = FALSE + ) + rbind(out, overall) +} + +#' The one-line reading of a spatial sorting bias +#' +#' Written for the run log, where a bare ratio would be ignored. The thresholds +#' are conventional rather than derived — Hijmans (2012) offers no cutoff, and +#' presenting one as though he did would be worse than a rule of thumb labelled +#' as one. +#' +#' @param ssb the overall spatial sorting bias +#' @return a single string +#' @keywords internal +spatial_bias_note <- function(ssb) { + if (is.na(ssb)) return("spatial sorting bias could not be computed") + if (ssb >= 0.8) { + return(paste0("spatial sorting bias ", signif(ssb, 2), + ": the folds are spatially fair")) + } + paste0("spatial sorting bias ", signif(ssb, 2), + ": held-out patches sit closer to training patches than non-patches ", + "do, so the AUC above is optimistic. Consider spatially blocked folds") +} + +#' The single spatial sorting bias from a per-fold table +#' +#' @param bias the result of [spatial_bias()], or `NULL` +#' @return the overall ratio, or `NA_real_` +#' @keywords internal +overall_ssb <- function(bias) { + if (is.null(bias) || !("fold" %in% names(bias))) return(NA_real_) + value <- bias$ssb[bias$fold == "overall"] + if (length(value) == 1) value else NA_real_ +} diff --git a/README.md b/README.md index 67bb115..f7ec6c2 100644 --- a/README.md +++ b/README.md @@ -1513,6 +1513,19 @@ Earth](https://www.naturalearthdata.com/), public domain, via `rnaturalearth`. 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` and `compare_runs()` +- Hijmans RJ (2012). Cross-validation of species distribution models: removing + spatial sorting bias and calibration with a null model. *Ecology* **93**(3), + 679–688. [doi:10.1890/11-0826.1](https://doi.org/10.1890/11-0826.1) — the + `ssb` row in `evals.csv` +- Schoener TW (1968). The *Anolis* lizards of Bimini: resource partitioning in + a complex fauna. *Ecology* **49**(4), 704–726. + [doi:10.2307/1935534](https://doi.org/10.2307/1935534) — Schoener's *D* in + `projection_overlap()` +- Warren DL, Glor RE, Turelli M (2008). Environmental niche equivalency versus + conservatism: quantitative approaches to niche evolution. *Evolution* + **62**(11), 2868–2883. + [doi:10.1111/j.1558-5646.2008.00482.x](https://doi.org/10.1111/j.1558-5646.2008.00482.x) + — Warren's *I* - Hoenig JM, Heisey DM (2001). The abuse of power: the pervasive fallacy of power calculations for data analysis. *The American Statistician* **55**(1), 19–24. diff --git a/man/evaluation_table.Rd b/man/evaluation_table.Rd index 3aaa506..fd32c5e 100644 --- a/man/evaluation_table.Rd +++ b/man/evaluation_table.Rd @@ -4,7 +4,7 @@ \alias{evaluation_table} \title{Assemble a self-explanatory evaluation table} \usage{ -evaluation_table(predictions, cv_metrics, cutoff, bounds = NULL) +evaluation_table(predictions, cv_metrics, cutoff, bounds = NULL, ssb = NULL) } \arguments{ \item{predictions}{held-out predictions from resampling} diff --git a/man/overall_ssb.Rd b/man/overall_ssb.Rd new file mode 100644 index 0000000..177ed53 --- /dev/null +++ b/man/overall_ssb.Rd @@ -0,0 +1,18 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/spatial_bias.R +\name{overall_ssb} +\alias{overall_ssb} +\title{The single spatial sorting bias from a per-fold table} +\usage{ +overall_ssb(bias) +} +\arguments{ +\item{bias}{the result of \code{\link[=spatial_bias]{spatial_bias()}}, or \code{NULL}} +} +\value{ +the overall ratio, or \code{NA_real_} +} +\description{ +The single spatial sorting bias from a per-fold table +} +\keyword{internal} diff --git a/man/projection_cells.Rd b/man/projection_cells.Rd new file mode 100644 index 0000000..710269a --- /dev/null +++ b/man/projection_cells.Rd @@ -0,0 +1,20 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/power.R +\name{projection_cells} +\alias{projection_cells} +\title{The cells of a projection, checked} +\usage{ +projection_cells(projection, label) +} +\arguments{ +\item{projection}{a projection table} + +\item{label}{which argument it came from, for the error} +} +\value{ +a data frame of \code{lon}, \code{lat} and \code{suitability} +} +\description{ +The cells of a projection, checked +} +\keyword{internal} diff --git a/man/projection_overlap.Rd b/man/projection_overlap.Rd new file mode 100644 index 0000000..1d4f8cd --- /dev/null +++ b/man/projection_overlap.Rd @@ -0,0 +1,73 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/power.R +\name{projection_overlap} +\alias{projection_overlap} +\title{How much two projected surfaces agree about the habitat} +\usage{ +projection_overlap(x, y, statistic = c("both", "D", "I"), digits = 6) +} +\arguments{ +\item{x, y}{projections to compare: data frames with \code{lon}, \code{lat} and +\code{suitability}, as \code{\link[=project_patch_model]{project_patch_model()}} writes to \code{suitability.csv}} + +\item{statistic}{\code{"both"}, \code{"D"}, or \code{"I"}} + +\item{digits}{how many decimal places the coordinates are matched on. +Projections from the same grid agree exactly; this is for two that came +off slightly different pipelines} +} +\value{ +a named numeric vector, with an \code{n_cells} attribute +} +\description{ +\code{\link[=compare_runs]{compare_runs()}} asks which run scores better. This asks something the +scores cannot: whether the two maps are the \emph{same map}. Two runs can rank +stations equally well and disagree completely about where the habitat is, +and a comparison of AUCs would report them as equivalent. +} +\details{ +Schoener's \emph{D} and Warren's \emph{I}, both on 0 to 1. \strong{1} is identical +surfaces; \strong{0} is no overlap at all. There is no threshold at which two +niches become "the same" — the number is a description, and what counts as +close depends on what the maps are for. +} +\section{What it is good for}{ + +Two questions, mostly. \emph{Do two species occupy the same habitat?} — run +\code{cfin} and \code{ctyp} over the same months and compare their surfaces. And \emph{does +the choice of algorithm change the map?} — the more useful one when a +comparison of scores came back inconclusive, because a small difference in +AUC with a low overlap means the two models are doing genuinely different +things and the tie in performance is hiding it. +} + +\section{Matching the cells}{ + +The two surfaces must be the same cells in the same order, or the overlap is +measured between places that have nothing to do with each other. Cells are +matched on their coordinates rather than on row position, and cells present +in one projection and not the other are dropped with a warning — a +projection that covers less ground usually has a reason, and averaging over +the difference would hide it. +} + +\examples{ +\dontrun{ +rf <- readr::read_csv("output/rf/projections/suitability.csv") +gam <- readr::read_csv("output/gam/projections/suitability.csv") + +projection_overlap(rf[rf$month == 7, ], gam[gam$month == 7, ]) +} +} +\references{ +Schoener TW (1968). The \emph{Anolis} lizards of Bimini: resource partitioning in +a complex fauna. \emph{Ecology} \strong{49}(4), 704-726. \doi{10.2307/1935534} — \emph{D} + +Warren DL, Glor RE, Turelli M (2008). Environmental niche equivalency versus +conservatism: quantitative approaches to niche evolution. \emph{Evolution} +\strong{62}(11), 2868-2883. \doi{10.1111/j.1558-5646.2008.00482.x} — \emph{I} +} +\seealso{ +\code{\link[=compare_runs]{compare_runs()}}, which asks which surface is \emph{better} rather than +whether they are the same +} diff --git a/man/spatial_bias.Rd b/man/spatial_bias.Rd new file mode 100644 index 0000000..0c6224c --- /dev/null +++ b/man/spatial_bias.Rd @@ -0,0 +1,70 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/spatial_bias.R +\name{spatial_bias} +\alias{spatial_bias} +\title{How much of the AUC is an artefact of where the stations are} +\usage{ +spatial_bias(model, geo = TRUE) +} +\arguments{ +\item{model}{a fitted model from \code{\link[=fit_patch_model]{fit_patch_model()}} or +\code{\link[=fit_patch_ensemble]{fit_patch_ensemble()}}, carrying \code{coordinates} and held-out \code{predictions}} + +\item{geo}{whether to measure great-circle distances. \code{TRUE} by default, +because the coordinates are longitude and latitude, and a degree of +longitude is not a degree of latitude anywhere but the equator} +} +\value{ +a data frame with one row per fold and a \code{fold} of \code{"overall"} for +the mean, carrying \code{patch_distance}, \code{non_patch_distance} and \code{ssb} +} +\description{ +A cross-validated AUC is optimistic when the held-out patches sit closer to +the training patches than the held-out non-patches do. The model can then +score well by recognising \emph{where} it has already been rather than \emph{what} +makes a patch, and nothing in the ROC curve distinguishes the two. +} +\details{ +This is the spatial sorting bias of Hijmans (2012): the mean distance from +each held-out patch to the nearest training patch, over the same distance +for the held-out non-patches. + +\deqn{SSB = \frac{\overline{d}(\mathrm{held\ out\ patch},\ \mathrm{training\ patch})}{\overline{d}(\mathrm{held\ out\ non\ patch},\ \mathrm{training\ patch})}} +\itemize{ +\item \strong{Near 1} — held-out patches and non-patches are equally far from the +training patches. The split is spatially fair and the AUC means what it +appears to. +\item \strong{Near 0} — held-out patches are much closer to training patches than the +non-patches are. \strong{The AUC is inflated}, and by an amount this number +does not tell you. +} + +It is a property of the \emph{split}, not of the model. Every model type fitted on +the same folds gets the same value, which is the point: it says how much to +trust the comparison between them, not which of them won. +} +\section{Why an ECOMON run should look at it}{ + +Survey stations are not scattered at random. They sit on transects, revisited +season after season, so a randomly held-out station usually has a near +neighbour in the training folds — often the same station in another year. +That is the situation this measures, and random \code{model.cv_folds} cannot avoid +it. A low value is an argument for spatially blocked folds, not for a +different model. +} + +\examples{ +\dontrun{ +model <- fit_patch_model(dat, config) +spatial_bias(model) +} +} +\references{ +Hijmans RJ (2012). Cross-validation of species distribution models: removing +spatial sorting bias and calibration with a null model. \emph{Ecology} \strong{93}(3), +679-688. \doi{10.1890/11-0826.1} +} +\seealso{ +\code{\link[=evaluation_table]{evaluation_table()}}, which reports the summary beside the AUC it +qualifies +} diff --git a/man/spatial_bias_note.Rd b/man/spatial_bias_note.Rd new file mode 100644 index 0000000..bd90bd4 --- /dev/null +++ b/man/spatial_bias_note.Rd @@ -0,0 +1,21 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/spatial_bias.R +\name{spatial_bias_note} +\alias{spatial_bias_note} +\title{The one-line reading of a spatial sorting bias} +\usage{ +spatial_bias_note(ssb) +} +\arguments{ +\item{ssb}{the overall spatial sorting bias} +} +\value{ +a single string +} +\description{ +Written for the run log, where a bare ratio would be ignored. The thresholds +are conventional rather than derived — Hijmans (2012) offers no cutoff, and +presenting one as though he did would be worse than a rule of thumb labelled +as one. +} +\keyword{internal} diff --git a/man/station_coordinates.Rd b/man/station_coordinates.Rd new file mode 100644 index 0000000..7008018 --- /dev/null +++ b/man/station_coordinates.Rd @@ -0,0 +1,19 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/model.R +\name{station_coordinates} +\alias{station_coordinates} +\title{Station coordinates, in the order the modelling data is in} +\usage{ +station_coordinates(dat) +} +\arguments{ +\item{dat}{labeled modeling data} +} +\value{ +a two-column data frame of \code{lon} and \code{lat}, or \code{NULL} +} +\description{ +\code{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. +} +\keyword{internal} diff --git a/tests/testthat/test-power.R b/tests/testthat/test-power.R index 3e32779..c6abd3c 100644 --- a/tests/testthat/test-power.R +++ b/tests/testthat/test-power.R @@ -211,3 +211,77 @@ test_that("a power curve is reproducible", { expect_equal(once$difference, twice$difference) expect_equal(once$std_err, twice$std_err) }) + +# ---- projection overlap ------------------------------------------------------ + +overlap_cells <- function(seed = 1) { + cells <- expand.grid(lon = seq(-70, -66, by = 0.5), + lat = seq(41, 44, by = 0.5)) + set.seed(seed) + cells$suitability <- stats::runif(nrow(cells)) + cells +} + +test_that("a surface overlaps itself completely", { + cells <- overlap_cells() + + out <- projection_overlap(cells, cells) + + expect_equal(unname(out[["D"]]), 1) + expect_equal(unname(out[["I"]]), 1) + expect_equal(attr(out, "n_cells"), nrow(cells)) +}) + +test_that("a surface overlaps its own inverse less than itself", { + cells <- overlap_cells() + flipped <- cells + flipped$suitability <- 1 - cells$suitability + + expect_lt(projection_overlap(cells, flipped)[["D"]], + projection_overlap(cells, cells)[["D"]]) +}) + +test_that("one statistic can be asked for on its own", { + cells <- overlap_cells() + + expect_named(projection_overlap(cells, cells, statistic = "D"), "D") + expect_named(projection_overlap(cells, cells, statistic = "I"), "I") + expect_setequal(names(projection_overlap(cells, cells)), c("D", "I")) +}) + +test_that("cells are matched on coordinates rather than on row order", { + # Two projections written by different runs need not be in the same order, + # and comparing row i of one with row i of the other would be comparing two + # different places. + cells <- overlap_cells() + shuffled <- cells[order(stats::runif(nrow(cells))), ] + + expect_equal(unname(projection_overlap(cells, shuffled)[["D"]]), 1) +}) + +test_that("projections covering different ground are intersected, loudly", { + cells <- overlap_cells() + partial <- cells[cells$lat <= 43, ] + + expect_warning(out <- projection_overlap(cells, partial), + "do not cover the same cells") + expect_equal(attr(out, "n_cells"), nrow(partial)) +}) + +test_that("projections with no common ground are refused", { + cells <- overlap_cells() + elsewhere <- cells + elsewhere$lon <- elsewhere$lon + 100 + + expect_error(suppressWarnings(projection_overlap(cells, elsewhere)), + "share 0 cells") +}) + +test_that("a table that is not a projection says what it is missing", { + cells <- overlap_cells() + + expect_error(projection_overlap(cells[, c("lon", "lat")], cells), + "missing: suitability") + expect_error(projection_overlap(cells, cells[, c("lon", "suitability")]), + "missing: lat") +}) diff --git a/tests/testthat/test-uncertainty.R b/tests/testthat/test-uncertainty.R index bcf11f1..5314c68 100644 --- a/tests/testthat/test-uncertainty.R +++ b/tests/testthat/test-uncertainty.R @@ -169,3 +169,88 @@ test_that("the extra surfaces reach the GeoTIFF as named layers", { expect_true(any(grepl("_uncertainty\\.png$", plots))) expect_true(any(!grepl("_uncertainty\\.png$", plots))) }) + +# ---- spatial sorting bias --------------------------------------------------- + +test_that("spatially fair folds score near 1 and sorted ones score near 0", { + # Built rather than fitted, so the answer is known in advance. Patches and + # non-patches drawn from the same places is the fair case; patches clustered + # away from the non-patches is the sorted one. + set.seed(1) + n <- 200 + fair <- data.frame(lon = runif(n, -70, -66), lat = runif(n, 41, 44)) + # Shuffled rather than alternating: `rep` of a 2-cycle across a 4-fold cycle + # puts every patch in the odd folds, leaving each fold with one class and + # nothing to compare. + labels <- sample(rep(c("patch", "non_patch"), length.out = n)) + predictions <- data.frame( + .row = seq_len(n), + id = rep(paste0("Fold", 1:4), length.out = n), + patch = factor(labels, levels = c("patch", "non_patch")) + ) + + unbiased <- spatial_bias(list(coordinates = fair, predictions = predictions)) + expect_equal(overall_ssb(unbiased), 1, tolerance = 0.35) + + # Now put every patch in one corner and every non-patch in another. + sorted <- fair + is_patch <- predictions$patch == "patch" + sorted$lon[is_patch] <- runif(sum(is_patch), -70, -69.5) + sorted$lon[!is_patch] <- runif(sum(!is_patch), -66.5, -66) + + biased <- spatial_bias(list(coordinates = sorted, predictions = predictions)) + expect_lt(overall_ssb(biased), overall_ssb(unbiased)) + expect_lt(overall_ssb(biased), 0.3) +}) + +test_that("spatial bias reports every fold as well as the overall", { + set.seed(2) + n <- 120 + coordinates <- data.frame(lon = runif(n, -70, -66), lat = runif(n, 41, 44)) + predictions <- data.frame( + .row = seq_len(n), + id = rep(paste0("Fold", 1:3), length.out = n), + patch = factor(sample(rep(c("patch", "non_patch"), length.out = n)), + levels = c("patch", "non_patch")) + ) + + out <- spatial_bias(list(coordinates = coordinates, predictions = predictions)) + + expect_equal(nrow(out), 4) + expect_equal(out$fold, c("Fold1", "Fold2", "Fold3", "overall")) + expect_true(all(out$ssb > 0)) +}) + +test_that("spatial bias declines to answer without coordinates", { + # A model fitted on a hand-built frame has none, and that is not an error. + predictions <- data.frame(.row = 1:4, id = "Fold1", + patch = factor(c("patch", "non_patch"), + levels = c("patch", "non_patch"))) + + expect_null(spatial_bias(list(coordinates = NULL, predictions = predictions))) + expect_null(spatial_bias(list(coordinates = data.frame(lon = 1, lat = 1), + predictions = NULL))) + expect_true(is.na(overall_ssb(NULL))) +}) + +test_that("the spatial bias note says what a low value means for the AUC", { + expect_match(spatial_bias_note(0.95), "spatially fair") + expect_match(spatial_bias_note(0.2), "optimistic") + expect_match(spatial_bias_note(NA_real_), "could not be computed") +}) + +test_that("the fitted model carries its coordinates and its spatial bias", { + skip_on_cran() + config <- mock_config() + config$model$trees <- 50 + config$model$cv_folds <- 5 + dat <- labeled_mock_data(config) + + model <- fit_patch_model(dat, config) + + expect_equal(nrow(model$coordinates), nrow(model$model_data)) + expect_setequal(names(model$coordinates), c("lon", "lat")) + expect_false(is.null(model$spatial_bias)) + # And it reaches the evaluation table, beside the metrics it qualifies. + expect_true("ssb" %in% model$evaluation$metric) +}) diff --git a/tools/citations.csv b/tools/citations.csv index 3ece160..c9c1d7d 100644 --- a/tools/citations.csv +++ b/tools/citations.csv @@ -41,3 +41,6 @@ elith2011,10.1111/j.1472-4642.2010.00725.x,,Elith,2011,A statistical explanation nadeau2003,10.1023/A:1024068626366,,Nadeau,2003,Inference for the Generalization Error,Machine Learning,52,239-281,crossref, dietterich1998,10.1162/089976698300017197,,Dietterich,1998,Approximate Statistical Tests for Comparing Supervised Classification Learning Algorithms,Neural Computation,10,1895-1923,crossref, hoenig2001,10.1198/000313001300339897,,Hoenig,2001,The Abuse of Power,The American Statistician,55,19-24,crossref, +hijmans2012,10.1890/11-0826.1,,Hijmans,2012,Cross‐validation of species distribution models: removing spatial sorting bias and calibration with a null model,Ecology,93,679-688,crossref, +schoener1968,10.2307/1935534,,Schoener,1968,The Anolis Lizards of Bimini: Resource Partitioning in a Complex Fauna,Ecology,49,704-726,crossref, +warren2008,10.1111/j.1558-5646.2008.00482.x,,Warren,2008,ENVIRONMENTAL NICHE EQUIVALENCY VERSUS CONSERVATISM: QUANTITATIVE APPROACHES TO NICHE EVOLUTION,Evolution,62,2868-2883,crossref, From 6857a089e5723b1ee0bdc1be4878621906375385 Mon Sep 17 00:00:00 2001 From: chross22 <52218551+chross22@users.noreply.github.com> Date: Wed, 12 Aug 2026 13:34:16 -0400 Subject: [PATCH 5/5] Offer a derived covariate whether or not its source is in the model Two gaps in the same function, both of which made a reasonable request unaskable through the app. Integrated chlorophyll without chlorophyll itself is an ordinary thing to want. The accumulated bloom is what feeds the animals sampled a month later; the instantaneous value at the moment of sampling is a different quantity and often not the interesting one. It could not be asked for. derivoce_choices() built its per-covariate options by looping over the *selected* covariates, so CHL_int appeared only once CHL had been made a predictor - and there was then no way to take CHL back out. The loop now runs over the whole catalogue, and the options that would cost a download say so in their group rather than their label, so the picker shows the cost without every entry carrying a parenthesis. The app already puts such ingredients in covariates.exclude, so they are fetched and not modelled, which is what that mechanism was built for. The second gap is the gradient of the current speed. UO and VO could each be differentiated and the speed could be computed, but the speed could not then be differentiated - and the gradient of the speed is the quantity a front in the flow actually is. Two components can each be changing steeply while the speed is constant, which is a turn rather than a shear. It is also the original pipeline's uv_grad, so the package could produce it from a hand-written config and the app could not. That needed a notion the choices did not have: a step that depends on another step rather than on a download. `depends` is that, and with_dependencies() follows it, so choosing speed_grad pulls current_speed in. Order needs no sorting - the choices already list a dependency before anything that reads it, and steps are emitted in that order. Every per-covariate step is offered on the speed, not only the gradient, and EKE gets a gradient too. A dependency is an ingredient in the same sense as a download, so derivoce_dependency_columns() names the derived columns a selection produced without anyone choosing them, and the app excludes those too. Picking speed_grad alone computes the speed and does not model it; picking both models both. Without this the speed would have arrived as a predictor nobody asked for, which is the bug the download side already avoided. Co-Authored-By: Claude Opus 5 --- NAMESPACE | 1 + R/derivoce.R | 119 ++++++++++++++++++++++++++++- inst/shiny/app.R | 10 ++- man/derivoce_dependency_columns.Rd | 42 ++++++++++ man/projection_overlap.Rd | 11 +++ man/with_dependencies.Rd | 29 +++++++ tests/testthat/test-derivoce.R | 63 +++++++++++++++ 7 files changed, 270 insertions(+), 5 deletions(-) create mode 100644 man/derivoce_dependency_columns.Rd create mode 100644 man/with_dependencies.Rd diff --git a/NAMESPACE b/NAMESPACE index f484643..7637f4c 100644 --- a/NAMESPACE +++ b/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) diff --git a/R/derivoce.R b/R/derivoce.R index 8fea5d7..d65e8b5 100644 --- a/R/derivoce.R +++ b/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))) +} diff --git a/inst/shiny/app.R b/inst/shiny/app.R index 70956ad..325184e 100644 --- a/inst/shiny/app.R +++ b/inst/shiny/app.R @@ -920,7 +920,15 @@ server <- function(input, output, session) { input$bathymetry %||% character()) config$covariates$selected <- union(input$covariates %||% character(), ingredients) - config$covariates$exclude <- ingredients + # A derived covariate can also be an ingredient: the gradient of current + # speed computes the speed on the way, and that column is no more a + # predictor than the velocity components behind it are. + config$covariates$exclude <- union( + ingredients, + derivoce_dependency_columns(input$derived %||% character(), + input$covariates %||% character(), + input$bathymetry %||% character()) + ) config$covariates$bathymetry <- input$bathymetry %||% character() config$covariates$climate <- input$climate %||% character() diff --git a/man/derivoce_dependency_columns.Rd b/man/derivoce_dependency_columns.Rd new file mode 100644 index 0000000..d9a6dd3 --- /dev/null +++ b/man/derivoce_dependency_columns.Rd @@ -0,0 +1,42 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/derivoce.R +\name{derivoce_dependency_columns} +\alias{derivoce_dependency_columns} +\title{Derived columns pulled in only as ingredients} +\usage{ +derivoce_dependency_columns(ids, selected, bathymetry = character()) +} +\arguments{ +\item{ids}{chosen derived covariate ids} + +\item{selected}{time-varying covariate names} + +\item{bathymetry}{static seafloor covariate names} +} +\value{ +character vector of derived column names, possibly empty +} +\description{ +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. +} +\details{ +This names the derived columns a selection produced without anyone choosing +them, so they can go into \code{covariates.exclude} alongside the downloads that +\code{\link[=derivoce_required_inputs]{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. +} +\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{ +\code{\link[=derivoce_required_inputs]{derivoce_required_inputs()}}, which does the same for downloads +} diff --git a/man/projection_overlap.Rd b/man/projection_overlap.Rd index 1d4f8cd..5e587bf 100644 --- a/man/projection_overlap.Rd +++ b/man/projection_overlap.Rd @@ -31,6 +31,17 @@ surfaces; \strong{0} is no overlap at all. There is no threshold at which two niches become "the same" — the number is a description, and what counts as close depends on what the maps are for. } +\section{Read them as high numbers}{ + +Both statistics compare the suitability \emph{distributions} cell by cell after +normalising them, so they run high. Two surfaces of independent random noise +over the same grid score around 0.7 on \emph{D} and 0.9 on \emph{I} — not because they +agree about anything, but because both spread their probability over the +same cells in similar proportions. Treat those as the floor for surfaces of +this shape rather than reading 0.9 as near-identity, and compare overlaps +against each other rather than against 1. +} + \section{What it is good for}{ Two questions, mostly. \emph{Do two species occupy the same habitat?} — run diff --git a/man/with_dependencies.Rd b/man/with_dependencies.Rd new file mode 100644 index 0000000..873a98d --- /dev/null +++ b/man/with_dependencies.Rd @@ -0,0 +1,29 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/derivoce.R +\name{with_dependencies} +\alias{with_dependencies} +\title{The chosen derived covariates, plus the ones they are computed from} +\usage{ +with_dependencies(ids, choices) +} +\arguments{ +\item{ids}{chosen derived covariate ids} + +\item{choices}{the \code{\link[=derivoce_choices]{derivoce_choices()}} list} +} +\value{ +\code{ids} with any dependencies added +} +\description{ +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. +} +\details{ +Order is not this function's problem. \code{derivoce_choices()} lists a +dependency before anything that depends on it, and steps are emitted in that +order, so \code{current_speed} runs before the gradient of \code{speed} without +anything having to sort them. +} +\keyword{internal} diff --git a/tests/testthat/test-derivoce.R b/tests/testthat/test-derivoce.R index 183bf2a..6d665d6 100644 --- a/tests/testthat/test-derivoce.R +++ b/tests/testthat/test-derivoce.R @@ -274,3 +274,66 @@ test_that("a run with no derivoce block is untouched", { env <- fetch_covariates(config, years = 2018, months = 6) expect_identical(add_derivoce_covariates(env, config), env) }) + +test_that("a derived covariate is offered whether or not its source is modelled", { + # Integrated chlorophyll without chlorophyll itself is an ordinary thing to + # want - the accumulated bloom is what feeds the animals, the instantaneous + # value is not - and it used to be unaskable: the derived form appeared only + # once its source was a predictor, and then the source could not be removed. + ids <- vapply(derivoce_choices(c("SST", "SSS")), function(x) x$id, + character(1)) + + expect_true("CHL_int" %in% ids) + expect_true("CHL_lag1" %in% ids) + expect_true("SST_int" %in% ids) + # And the source is fetched, so it can be excluded from the predictors. + expect_equal(derivoce_required_inputs("CHL_int", c("SST", "SSS")), "CHL") + # A source that *is* selected needs no extra download. + expect_length(derivoce_required_inputs("SST_int", c("SST", "SSS")), 0) +}) + +test_that("the extra choices say in their group that they cost a download", { + choices <- derivoce_choices(c("SST")) + group_of <- function(id) { + Filter(function(x) x$id == id, choices)[[1]]$group + } + + expect_equal(group_of("SST_int"), "Temporal") + expect_match(group_of("CHL_int"), "downloads CHL") +}) + +test_that("current speed can be differentiated, not just its components", { + # The gradient of the speed is the quantity a front in the flow is. The + # components can each be changing steeply while the speed is constant, which + # is a turn rather than a shear. + ids <- vapply(derivoce_choices(c("SST")), function(x) x$id, character(1)) + expect_true("speed_grad" %in% ids) + expect_true("EKE_grad" %in% ids) + + steps <- derivoce_steps_for("speed_grad", c("SST")) + + # The speed has to be computed before it can be differentiated, and the + # config is read in order. + expect_equal(vapply(steps, function(s) s$type, character(1)), + c("current_speed", "horizontal_gradient")) + expect_equal(steps[[2]]$vars, "speed") + expect_setequal(derivoce_required_inputs("speed_grad", c("SST")), + c("UO", "VO")) +}) + +test_that("a dependency pulled in on the way is not made a predictor", { + expect_equal(derivoce_dependency_columns("speed_grad", c("SST")), "speed") + # Unless it was asked for in its own right. + expect_length(derivoce_dependency_columns(c("speed", "speed_grad"), c("SST")), + 0) + # A download is not a derived column; derivoce_required_inputs covers those. + expect_length(derivoce_dependency_columns("CHL_int", c("SST")), 0) + expect_length(derivoce_dependency_columns(character(), c("SST")), 0) +}) + +test_that("asking for a derived covariate twice over does not duplicate a step", { + steps <- derivoce_steps_for(c("speed", "speed_grad"), c("SST")) + + expect_equal(sum(vapply(steps, function(s) s$type, character(1)) == + "current_speed"), 1) +})