From f3840f19ad1afbbf01a6119749c9bbbaadc64cc3 Mon Sep 17 00:00:00 2001 From: chross22 <52218551+chross22@users.noreply.github.com> Date: Mon, 10 Aug 2026 16:31:13 -0400 Subject: [PATCH 1/2] Let evaluation take predictions, not only a model that can remake them Every evaluation function here re-predicts. That is the right default - it keeps the scored predictions and the model provably in step - but it assumes the caller holds a model that can reproduce them, and a cross-validated workflow does not. Under k-fold cross-validation each observation is predicted by the one fold model that did not see it. The honest predictions are spread across k models, none of which is the final fit, so by the time a pipeline has a single model to hand it has either discarded them or kept them with nothing to pass them to. Re-predicting from the final model on the same rows answers a different and more flattering question, and evaluation_pairs() offered no way in. held_out(observed, predicted) is that way in. It goes wherever a model goes: pairs <- held_out(cv$observed, cv$predicted) plotROC(pairs, folds = cv$fold) plotThreshold(pairs, folds = cv$fold) plotCalibration(pairs) evaluation_pairs() short-circuits on it, so threshold_metrics(), calibration_estimates() and the three plots over them all work unchanged, folds included. newdata stops being required when the predictions arrive with the object, and the error when it is still missing now names held_out() as the alternative. The test that matters is that the two doors reach the same room: scored through the model and scored through held_out(), the metric table, the AUC, the prevalence and the n all come back identical, with and without folds. A second implementation of the same arithmetic would be worse than no second entry point at all. Two things it deliberately does not do. It cannot verify the predictions are out of sample - nothing in two numeric vectors records what made them - so in.sample is taken on trust, defaults to FALSE, and can be set TRUE to get the same annotation the model path applies when it recognises its own training data. And it cannot serve plotImportance() or permutation_importance(), which shuffle a predictor and re-predict; that needs a model by construction rather than a record of what one once said. Co-Authored-By: Claude Opus 5 --- NAMESPACE | 2 + R/calibration.R | 24 +++-- R/evaluate.R | 15 ++- R/held_out.R | 184 +++++++++++++++++++++++++++++++++ R/plotROC.R | 4 +- man/as_binary_outcome.Rd | 22 ++++ man/calc_deviance.Rd | 1 + man/calibration_estimates.Rd | 3 +- man/held_out.Rd | 92 +++++++++++++++++ man/held_out_pairs.Rd | 25 +++++ man/permutation_importance.Rd | 1 + man/plotCalibration.Rd | 3 +- man/plotImportance.Rd | 1 + man/plotROC.Rd | 3 +- man/plotThreshold.Rd | 3 +- man/print.fancyfx_held_out.Rd | 19 ++++ man/spatial_sorting_bias.Rd | 1 + man/threshold_metrics.Rd | 1 + tests/testthat/test-held-out.R | 136 ++++++++++++++++++++++++ 19 files changed, 523 insertions(+), 17 deletions(-) create mode 100644 R/held_out.R create mode 100644 man/as_binary_outcome.Rd create mode 100644 man/held_out.Rd create mode 100644 man/held_out_pairs.Rd create mode 100644 man/print.fancyfx_held_out.Rd create mode 100644 tests/testthat/test-held-out.R diff --git a/NAMESPACE b/NAMESPACE index dc4224a..5723e29 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -6,6 +6,7 @@ S3method(effect_estimates,gamm) S3method(effect_estimates,gamm4) S3method(effect_estimates,gbm) S3method(effect_estimates,scam) +S3method(print,fancyfx_held_out) export(calc_deviance) export(calibration_estimates) export(combinePlots) @@ -13,6 +14,7 @@ export(comparePlots) export(effect_estimates) export(ensemble_summary) export(fancyfx_palette) +export(held_out) export(hex_bin) export(mess) export(niche_equivalency) diff --git a/R/calibration.R b/R/calibration.R index efc7475..b5a57ca 100644 --- a/R/calibration.R +++ b/R/calibration.R @@ -68,16 +68,19 @@ #' attr(cal, "brier") #' #' @export -calibration_estimates <- function(model, newdata, bins = 10, +calibration_estimates <- function(model, newdata = NULL, bins = 10, binning = c("quantile", "width"), folds = NULL, level = 0.95, ...) { - if (missing(newdata) || is.null(newdata)) { + supplied <- inherits(model, "fancyfx_held_out") + if (!supplied && (missing(newdata) || is.null(newdata))) { stop("newdata is required: a model scored against the data it was fitted ", "to flatters itself. Supply held-out data, or the training data ", - "explicitly if that is genuinely what you want.", call. = FALSE) + "explicitly if that is genuinely what you want.\nTo score ", + "predictions you already have, wrap them with held_out().", + call. = FALSE) } - newdata <- as.data.frame(newdata) - model <- unwrap_gam(model) + if (supplied) newdata <- NULL else newdata <- as.data.frame(newdata) + if (!supplied) model <- unwrap_gam(model) binning <- check_choice(binning, c("quantile", "width"), "binning") check_level(level) @@ -254,7 +257,7 @@ calibration_fit <- function(observed, predicted) { #' plotCalibration(fit, dat[301:600, ]) #' #' @export -plotCalibration <- function(model, newdata, bins = 10, +plotCalibration <- function(model, newdata = NULL, bins = 10, binning = c("quantile", "width"), folds = NULL, level = 0.95, title = "", show.stats = TRUE, @@ -329,8 +332,13 @@ plotCalibration <- function(model, newdata, bins = 10, #' @return A numeric vector of predicted probabilities. #' @keywords internal predicted_for_rug <- function(model, newdata, ...) { - predicted <- predict_probability(unwrap_gam(model), as.data.frame(newdata), - ...) + # Predictions supplied directly are the rug: re-predicting would need a model + # that, by the time this path is used, the caller does not have. + predicted <- if (inherits(model, "fancyfx_held_out")) { + model$predicted + } else { + predict_probability(unwrap_gam(model), as.data.frame(newdata), ...) + } predicted[!is.na(predicted)] } diff --git a/R/evaluate.R b/R/evaluate.R index 8583fac..cfa6b65 100644 --- a/R/evaluate.R +++ b/R/evaluate.R @@ -79,12 +79,15 @@ #' #' @export threshold_metrics <- function(model, newdata, folds = NULL, ...) { - if (missing(newdata) || is.null(newdata)) { + supplied <- inherits(model, "fancyfx_held_out") + if (!supplied && (missing(newdata) || is.null(newdata))) { stop("newdata is required: a model scored against the data it was fitted ", "to flatters itself. Supply held-out data, or the training data ", - "explicitly if that is genuinely what you want.", call. = FALSE) + "explicitly if that is genuinely what you want.\nTo score ", + "predictions you already have, wrap them with held_out().", + call. = FALSE) } - newdata <- as.data.frame(newdata) + if (supplied) newdata <- NULL else newdata <- as.data.frame(newdata) # gamm4 and gamm hand back a wrapper that formula() and predict() both refuse. pairs <- evaluation_pairs(model, newdata, folds, ...) observed <- pairs$observed @@ -149,6 +152,12 @@ threshold_metrics <- function(model, newdata, folds = NULL, ...) { #' @keywords internal evaluation_pairs <- function(model, newdata, folds = NULL, require.both.classes = TRUE, ...) { + # Predictions supplied directly carry everything this function exists to + # produce, so there is nothing to predict and nothing to unwrap. + if (inherits(model, "fancyfx_held_out")) { + return(held_out_pairs(model, folds, require.both.classes)) + } + # Unwrapped here rather than in each caller: gamm4 and gamm hand back a # wrapper that formula() and predict() both refuse, and every evaluation # function reaches this point. diff --git a/R/held_out.R b/R/held_out.R new file mode 100644 index 0000000..8e77af3 --- /dev/null +++ b/R/held_out.R @@ -0,0 +1,184 @@ +#' Evaluate predictions you already have +#' +#' Every evaluation function here takes a fitted model and re-predicts. That is +#' the right default -- it keeps the scored predictions and the model provably +#' in step -- but it assumes the caller is holding a model that can reproduce +#' them, and a cross-validated workflow is not. +#' +#' Under k-fold cross-validation each observation is predicted by the one fold +#' model that did not see it. The honest predictions are therefore spread across +#' `k` models, none of which is the final fit, and by the time a pipeline has a +#' single model to hand it has already thrown them away -- or, more often, kept +#' them and has nothing to pass them to. Re-predicting from the final model on +#' the same rows answers a different and more flattering question. +#' +#' `held_out()` is the way in for those. Wrap the observed outcomes and the +#' predictions that were made for them, and pass the result anywhere a model +#' would go: +#' +#' ```r +#' pairs <- held_out(cv$observed, cv$predicted) +#' plotROC(pairs, folds = cv$fold) +#' plotThreshold(pairs, folds = cv$fold) +#' plotCalibration(pairs) +#' ``` +#' +#' @section What it does not do: +#' It cannot check the predictions are out of sample. Nothing in a pair of +#' numeric vectors records which model made them or what it was fitted to, so +#' `in.sample` is taken on trust -- the argument exists to be set honestly, and +#' defaults to `FALSE` because that is what the function is named for. +#' +#' That is a real difference from the model path, which inspects the fit and +#' warns when it recognises its own training data. Passing training predictions +#' here gets no warning, because there is nothing to notice it with. +#' +#' It also cannot support [plotImportance()] or [permutation_importance()], +#' which shuffle a predictor and re-predict. That needs a model by construction, +#' not a record of what one once said. +#' +#' @param observed Observed outcomes: `0`/`1`, a logical, or a two-level factor +#' whose **second** level is the positive case, matching how [stats::glm()] +#' treats one. +#' @param predicted Predicted probabilities, one per element of `observed`. +#' @param in.sample Whether these predictions were made on the data the model +#' was fitted to. `FALSE` by default; set `TRUE` and every plot built from +#' them is annotated as in-sample, exactly as the model path would. +#' +#' @return An object of class `fancyfx_held_out`, accepted wherever a model is. +#' +#' @family evaluation plots +#' @seealso [threshold_metrics()], [plotROC()], [plotThreshold()], +#' [plotCalibration()]. +#' +#' @examples +#' set.seed(1) +#' truth <- rbinom(200, 1, 0.3) +#' score <- plogis(rnorm(200, ifelse(truth == 1, 1, -1))) +#' +#' pairs <- held_out(truth, score) +#' metrics <- threshold_metrics(pairs) +#' metrics$.threshold[which.max(metrics$.tss)] +#' +#' # Fold-wise, when the predictions came from cross-validation. +#' folds <- rep(1:5, length.out = 200) +#' head(threshold_metrics(pairs, folds = folds)) +#' +#' @export +held_out <- function(observed, predicted, in.sample = FALSE) { + observed <- as_binary_outcome(observed) + predicted <- as.numeric(predicted) + + if (length(observed) != length(predicted)) { + stop("observed and predicted must be the same length: ", length(observed), + " and ", length(predicted), ".", call. = FALSE) + } + if (!length(observed)) { + stop("observed and predicted are empty, so there is nothing to score.", + call. = FALSE) + } + finite <- predicted[is.finite(predicted)] + if (length(finite) && (min(finite) < 0 || max(finite) > 1)) { + stop("predicted must be probabilities in [0, 1], but they run from ", + format(min(finite)), " to ", format(max(finite)), + ". Predictions on the link scale need transforming first.", + call. = FALSE) + } + if (!is.logical(in.sample) || length(in.sample) != 1) { + stop("in.sample must be TRUE or FALSE.", call. = FALSE) + } + + structure( + list(observed = observed, predicted = predicted, in.sample = in.sample), + class = "fancyfx_held_out" + ) +} + +#' Coerce observed outcomes to 0/1 +#' +#' The same three forms [binary_response()] accepts, and the same reading of +#' each, so a `held_out()` pair and a model scored on a data frame agree about +#' which class is positive. Split out rather than shared with +#' [binary_response()] because that one reaches into `newdata` for a column +#' named by the model's formula, and here there is no model and no column. +#' +#' @param observed Observed outcomes. +#' @return A 0/1 numeric vector. +#' @keywords internal +as_binary_outcome <- function(observed) { + if (is.factor(observed)) { + if (nlevels(observed) != 2) { + stop("observed has ", nlevels(observed), " levels. Classification ", + "metrics are defined for a binary outcome only.", call. = FALSE) + } + # Second level is the positive case, as glm() itself treats a factor. + return(as.numeric(observed) - 1) + } + + if (is.logical(observed)) return(as.numeric(observed)) + + values <- unique(stats::na.omit(observed)) + if (!is.numeric(observed) || !all(values %in% c(0, 1))) { + stop("observed is not a binary outcome (found: ", + paste(utils::head(sort(values), 4), collapse = ", "), + if (length(values) > 4) ", ..." else "", + "). AUC and TSS are defined for presence/absence only -- applied to a ", + "continuous response they return a number with no meaning.", + call. = FALSE) + } + as.numeric(observed) +} + +#' The evaluation pairs a held_out() object already carries +#' +#' The short circuit in [evaluation_pairs()]. There is no model to unwrap, no +#' response column to find and no prediction to make; the work is the checking +#' that the model path does after predicting. +#' +#' @param x A `fancyfx_held_out` object. +#' @param folds Optional fold identifiers, one per observation. +#' @param require.both.classes Whether to refuse data containing only one +#' outcome class. +#' @return The same list [evaluation_pairs()] returns. +#' @keywords internal +held_out_pairs <- function(x, folds = NULL, require.both.classes = TRUE) { + observed <- x$observed + predicted <- x$predicted + + if (!is.null(folds) && length(folds) != length(observed)) { + stop("folds must have one entry per observation: ", length(observed), + " expected, ", length(folds), " given.", call. = FALSE) + } + + complete <- !is.na(observed) & !is.na(predicted) + observed <- observed[complete] + predicted <- predicted[complete] + + if (!length(observed)) { + stop("No observation has both an outcome and a prediction.", call. = FALSE) + } + if (require.both.classes && length(unique(observed)) < 2) { + stop("observed contains only one outcome class, so sensitivity and ", + "specificity are not both defined. Evaluation needs both presences ", + "and absences.", call. = FALSE) + } + + list(observed = observed, predicted = predicted, folds = folds, + complete = complete, in.sample = x$in.sample) +} + +#' Print a held_out object +#' +#' @param x A `fancyfx_held_out` object. +#' @param ... Unused. +#' @return `x`, invisibly. +#' @export +print.fancyfx_held_out <- function(x, ...) { + cat("\n") + cat(" observations: ", length(x$observed), "\n", sep = "") + cat(" prevalence: ", format(mean(x$observed, na.rm = TRUE), digits = 3), + "\n", sep = "") + cat(" in sample: ", if (isTRUE(x$in.sample)) "yes" else "no", "\n", + sep = "") + invisible(x) +} diff --git a/R/plotROC.R b/R/plotROC.R index f2af95c..b5cb6f3 100644 --- a/R/plotROC.R +++ b/R/plotROC.R @@ -52,7 +52,7 @@ #' plotROC(fit, test) #' #' @export -plotROC <- function(model, newdata, folds = NULL, title = "", +plotROC <- function(model, newdata = NULL, folds = NULL, title = "", show.auc = TRUE, theme = theme_fancyfx(), palette = fancyfx_palette(), @@ -167,7 +167,7 @@ auc_label <- function(auc.value) { #' plotThreshold(fit, dat[201:400, ], metrics = "tss") #' #' @export -plotThreshold <- function(model, newdata, folds = NULL, +plotThreshold <- function(model, newdata = NULL, folds = NULL, metrics = c("tss", "sensitivity", "specificity"), title = "", mark.best = TRUE, diff --git a/man/as_binary_outcome.Rd b/man/as_binary_outcome.Rd new file mode 100644 index 0000000..46fb68f --- /dev/null +++ b/man/as_binary_outcome.Rd @@ -0,0 +1,22 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/held_out.R +\name{as_binary_outcome} +\alias{as_binary_outcome} +\title{Coerce observed outcomes to 0/1} +\usage{ +as_binary_outcome(observed) +} +\arguments{ +\item{observed}{Observed outcomes.} +} +\value{ +A 0/1 numeric vector. +} +\description{ +The same three forms \code{\link[=binary_response]{binary_response()}} accepts, and the same reading of +each, so a \code{held_out()} pair and a model scored on a data frame agree about +which class is positive. Split out rather than shared with +\code{\link[=binary_response]{binary_response()}} because that one reaches into \code{newdata} for a column +named by the model's formula, and here there is no model and no column. +} +\keyword{internal} diff --git a/man/calc_deviance.Rd b/man/calc_deviance.Rd index aefa60c..a95e2c3 100644 --- a/man/calc_deviance.Rd +++ b/man/calc_deviance.Rd @@ -67,6 +67,7 @@ regression trees. \emph{Journal of Animal Ecology}, 77(4), 802-813. Other evaluation plots: \code{\link[=calibration_estimates]{calibration_estimates()}}, +\code{\link[=held_out]{held_out()}}, \code{\link[=permutation_importance]{permutation_importance()}}, \code{\link[=plotCalibration]{plotCalibration()}}, \code{\link[=plotImportance]{plotImportance()}}, diff --git a/man/calibration_estimates.Rd b/man/calibration_estimates.Rd index f0dc488..dc5c809 100644 --- a/man/calibration_estimates.Rd +++ b/man/calibration_estimates.Rd @@ -6,7 +6,7 @@ \usage{ calibration_estimates( model, - newdata, + newdata = NULL, bins = 10, binning = c("quantile", "width"), folds = NULL, @@ -93,6 +93,7 @@ discrimination. Other evaluation plots: \code{\link[=calc_deviance]{calc_deviance()}}, +\code{\link[=held_out]{held_out()}}, \code{\link[=permutation_importance]{permutation_importance()}}, \code{\link[=plotCalibration]{plotCalibration()}}, \code{\link[=plotImportance]{plotImportance()}}, diff --git a/man/held_out.Rd b/man/held_out.Rd new file mode 100644 index 0000000..db07764 --- /dev/null +++ b/man/held_out.Rd @@ -0,0 +1,92 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/held_out.R +\name{held_out} +\alias{held_out} +\title{Evaluate predictions you already have} +\usage{ +held_out(observed, predicted, in.sample = FALSE) +} +\arguments{ +\item{observed}{Observed outcomes: \code{0}/\code{1}, a logical, or a two-level factor +whose \strong{second} level is the positive case, matching how \code{\link[stats:glm]{stats::glm()}} +treats one.} + +\item{predicted}{Predicted probabilities, one per element of \code{observed}.} + +\item{in.sample}{Whether these predictions were made on the data the model +was fitted to. \code{FALSE} by default; set \code{TRUE} and every plot built from +them is annotated as in-sample, exactly as the model path would.} +} +\value{ +An object of class \code{fancyfx_held_out}, accepted wherever a model is. +} +\description{ +Every evaluation function here takes a fitted model and re-predicts. That is +the right default -- it keeps the scored predictions and the model provably +in step -- but it assumes the caller is holding a model that can reproduce +them, and a cross-validated workflow is not. +} +\details{ +Under k-fold cross-validation each observation is predicted by the one fold +model that did not see it. The honest predictions are therefore spread across +\code{k} models, none of which is the final fit, and by the time a pipeline has a +single model to hand it has already thrown them away -- or, more often, kept +them and has nothing to pass them to. Re-predicting from the final model on +the same rows answers a different and more flattering question. + +\code{held_out()} is the way in for those. Wrap the observed outcomes and the +predictions that were made for them, and pass the result anywhere a model +would go: + +\if{html}{\out{
}}\preformatted{pairs <- held_out(cv$observed, cv$predicted) +plotROC(pairs, folds = cv$fold) +plotThreshold(pairs, folds = cv$fold) +plotCalibration(pairs) +}\if{html}{\out{
}} +} +\section{What it does not do}{ + +It cannot check the predictions are out of sample. Nothing in a pair of +numeric vectors records which model made them or what it was fitted to, so +\code{in.sample} is taken on trust -- the argument exists to be set honestly, and +defaults to \code{FALSE} because that is what the function is named for. + +That is a real difference from the model path, which inspects the fit and +warns when it recognises its own training data. Passing training predictions +here gets no warning, because there is nothing to notice it with. + +It also cannot support \code{\link[=plotImportance]{plotImportance()}} or \code{\link[=permutation_importance]{permutation_importance()}}, +which shuffle a predictor and re-predict. That needs a model by construction, +not a record of what one once said. +} + +\examples{ +set.seed(1) +truth <- rbinom(200, 1, 0.3) +score <- plogis(rnorm(200, ifelse(truth == 1, 1, -1))) + +pairs <- held_out(truth, score) +metrics <- threshold_metrics(pairs) +metrics$.threshold[which.max(metrics$.tss)] + +# Fold-wise, when the predictions came from cross-validation. +folds <- rep(1:5, length.out = 200) +head(threshold_metrics(pairs, folds = folds)) + +} +\seealso{ +\code{\link[=threshold_metrics]{threshold_metrics()}}, \code{\link[=plotROC]{plotROC()}}, \code{\link[=plotThreshold]{plotThreshold()}}, +\code{\link[=plotCalibration]{plotCalibration()}}. + +Other evaluation plots: +\code{\link[=calc_deviance]{calc_deviance()}}, +\code{\link[=calibration_estimates]{calibration_estimates()}}, +\code{\link[=permutation_importance]{permutation_importance()}}, +\code{\link[=plotCalibration]{plotCalibration()}}, +\code{\link[=plotImportance]{plotImportance()}}, +\code{\link[=plotROC]{plotROC()}}, +\code{\link[=plotThreshold]{plotThreshold()}}, +\code{\link[=spatial_sorting_bias]{spatial_sorting_bias()}}, +\code{\link[=threshold_metrics]{threshold_metrics()}} +} +\concept{evaluation plots} diff --git a/man/held_out_pairs.Rd b/man/held_out_pairs.Rd new file mode 100644 index 0000000..9607302 --- /dev/null +++ b/man/held_out_pairs.Rd @@ -0,0 +1,25 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/held_out.R +\name{held_out_pairs} +\alias{held_out_pairs} +\title{The evaluation pairs a held_out() object already carries} +\usage{ +held_out_pairs(x, folds = NULL, require.both.classes = TRUE) +} +\arguments{ +\item{x}{A \code{fancyfx_held_out} object.} + +\item{folds}{Optional fold identifiers, one per observation.} + +\item{require.both.classes}{Whether to refuse data containing only one +outcome class.} +} +\value{ +The same list \code{\link[=evaluation_pairs]{evaluation_pairs()}} returns. +} +\description{ +The short circuit in \code{\link[=evaluation_pairs]{evaluation_pairs()}}. There is no model to unwrap, no +response column to find and no prediction to make; the work is the checking +that the model path does after predicting. +} +\keyword{internal} diff --git a/man/permutation_importance.Rd b/man/permutation_importance.Rd index 57fd7ae..8038326 100644 --- a/man/permutation_importance.Rd +++ b/man/permutation_importance.Rd @@ -89,6 +89,7 @@ effect rather than its weight. Other evaluation plots: \code{\link[=calc_deviance]{calc_deviance()}}, \code{\link[=calibration_estimates]{calibration_estimates()}}, +\code{\link[=held_out]{held_out()}}, \code{\link[=plotCalibration]{plotCalibration()}}, \code{\link[=plotImportance]{plotImportance()}}, \code{\link[=plotROC]{plotROC()}}, diff --git a/man/plotCalibration.Rd b/man/plotCalibration.Rd index b75fc53..97fa81b 100644 --- a/man/plotCalibration.Rd +++ b/man/plotCalibration.Rd @@ -6,7 +6,7 @@ \usage{ plotCalibration( model, - newdata, + newdata = NULL, bins = 10, binning = c("quantile", "width"), folds = NULL, @@ -88,6 +88,7 @@ discrimination, which is a different question. Other evaluation plots: \code{\link[=calc_deviance]{calc_deviance()}}, \code{\link[=calibration_estimates]{calibration_estimates()}}, +\code{\link[=held_out]{held_out()}}, \code{\link[=permutation_importance]{permutation_importance()}}, \code{\link[=plotImportance]{plotImportance()}}, \code{\link[=plotROC]{plotROC()}}, diff --git a/man/plotImportance.Rd b/man/plotImportance.Rd index be2379c..7937f61 100644 --- a/man/plotImportance.Rd +++ b/man/plotImportance.Rd @@ -79,6 +79,7 @@ shape of each effect. Other evaluation plots: \code{\link[=calc_deviance]{calc_deviance()}}, \code{\link[=calibration_estimates]{calibration_estimates()}}, +\code{\link[=held_out]{held_out()}}, \code{\link[=permutation_importance]{permutation_importance()}}, \code{\link[=plotCalibration]{plotCalibration()}}, \code{\link[=plotROC]{plotROC()}}, diff --git a/man/plotROC.Rd b/man/plotROC.Rd index 63fd40d..2194576 100644 --- a/man/plotROC.Rd +++ b/man/plotROC.Rd @@ -6,7 +6,7 @@ \usage{ plotROC( model, - newdata, + newdata = NULL, folds = NULL, title = "", show.auc = TRUE, @@ -82,6 +82,7 @@ the numbers underneath. Other evaluation plots: \code{\link[=calc_deviance]{calc_deviance()}}, \code{\link[=calibration_estimates]{calibration_estimates()}}, +\code{\link[=held_out]{held_out()}}, \code{\link[=permutation_importance]{permutation_importance()}}, \code{\link[=plotCalibration]{plotCalibration()}}, \code{\link[=plotImportance]{plotImportance()}}, diff --git a/man/plotThreshold.Rd b/man/plotThreshold.Rd index 24fd5a6..a7ffe2c 100644 --- a/man/plotThreshold.Rd +++ b/man/plotThreshold.Rd @@ -6,7 +6,7 @@ \usage{ plotThreshold( model, - newdata, + newdata = NULL, folds = NULL, metrics = c("tss", "sensitivity", "specificity"), title = "", @@ -80,6 +80,7 @@ plotThreshold(fit, dat[201:400, ], metrics = "tss") Other evaluation plots: \code{\link[=calc_deviance]{calc_deviance()}}, \code{\link[=calibration_estimates]{calibration_estimates()}}, +\code{\link[=held_out]{held_out()}}, \code{\link[=permutation_importance]{permutation_importance()}}, \code{\link[=plotCalibration]{plotCalibration()}}, \code{\link[=plotImportance]{plotImportance()}}, diff --git a/man/print.fancyfx_held_out.Rd b/man/print.fancyfx_held_out.Rd new file mode 100644 index 0000000..a963e05 --- /dev/null +++ b/man/print.fancyfx_held_out.Rd @@ -0,0 +1,19 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/held_out.R +\name{print.fancyfx_held_out} +\alias{print.fancyfx_held_out} +\title{Print a held_out object} +\usage{ +\method{print}{fancyfx_held_out}(x, ...) +} +\arguments{ +\item{x}{A \code{fancyfx_held_out} object.} + +\item{...}{Unused.} +} +\value{ +\code{x}, invisibly. +} +\description{ +Print a held_out object +} diff --git a/man/spatial_sorting_bias.Rd b/man/spatial_sorting_bias.Rd index bf40a32..e3251b7 100644 --- a/man/spatial_sorting_bias.Rd +++ b/man/spatial_sorting_bias.Rd @@ -72,6 +72,7 @@ blocked split gets used. Other evaluation plots: \code{\link[=calc_deviance]{calc_deviance()}}, \code{\link[=calibration_estimates]{calibration_estimates()}}, +\code{\link[=held_out]{held_out()}}, \code{\link[=permutation_importance]{permutation_importance()}}, \code{\link[=plotCalibration]{plotCalibration()}}, \code{\link[=plotImportance]{plotImportance()}}, diff --git a/man/threshold_metrics.Rd b/man/threshold_metrics.Rd index 7df79e0..31db3d7 100644 --- a/man/threshold_metrics.Rd +++ b/man/threshold_metrics.Rd @@ -95,6 +95,7 @@ species distribution models: prevalence, kappa and the true skill statistic Other evaluation plots: \code{\link[=calc_deviance]{calc_deviance()}}, \code{\link[=calibration_estimates]{calibration_estimates()}}, +\code{\link[=held_out]{held_out()}}, \code{\link[=permutation_importance]{permutation_importance()}}, \code{\link[=plotCalibration]{plotCalibration()}}, \code{\link[=plotImportance]{plotImportance()}}, diff --git a/tests/testthat/test-held-out.R b/tests/testthat/test-held-out.R new file mode 100644 index 0000000..9392630 --- /dev/null +++ b/tests/testthat/test-held-out.R @@ -0,0 +1,136 @@ +held_out_fixture <- function() { + set.seed(1) + d <- data.frame(x1 = runif(600, 1, 10), x2 = runif(600, 1, 10)) + d$y <- rbinom(600, 1, plogis(-3 + 0.6 * d$x1)) + train <- d[1:300, ] + test <- d[301:600, ] + fit <- glm(y ~ x1 + x2, data = train, family = binomial) + list(fit = fit, train = train, test = test, + predicted = unname(stats::predict(fit, test, type = "response"))) +} + +test_that("supplied predictions score identically to the model that made them", { + # The point of the whole entry point: it is a second door into the same + # arithmetic, not a second implementation of it. If these ever diverge, one + # of the two is wrong and there is no way to tell which from the outside. + f <- held_out_fixture() + + by_model <- threshold_metrics(f$fit, f$test) + by_pairs <- threshold_metrics(held_out(f$test$y, f$predicted)) + + expect_equal(as.data.frame(by_pairs), as.data.frame(by_model)) + expect_equal(attr(by_pairs, "auc"), attr(by_model, "auc")) + expect_equal(attr(by_pairs, "prevalence"), attr(by_model, "prevalence")) + expect_equal(attr(by_pairs, "n"), attr(by_model, "n")) +}) + +test_that("calibration agrees between the two paths too", { + f <- held_out_fixture() + + by_model <- calibration_estimates(f$fit, f$test) + by_pairs <- calibration_estimates(held_out(f$test$y, f$predicted)) + + expect_equal(as.data.frame(by_pairs), as.data.frame(by_model)) + expect_equal(attr(by_pairs, "calibration"), attr(by_model, "calibration")) + expect_equal(attr(by_pairs, "brier"), attr(by_model, "brier")) +}) + +test_that("the plots build from supplied predictions", { + f <- held_out_fixture() + pairs <- held_out(f$test$y, f$predicted) + + expect_s3_class(plotROC(pairs), "ggplot") + expect_s3_class(plotThreshold(pairs), "ggplot") + expect_s3_class(plotCalibration(pairs), "patchwork") +}) + +test_that("folds group supplied predictions the way they group a model's", { + f <- held_out_fixture() + folds <- rep(1:5, length.out = nrow(f$test)) + + by_model <- suppressMessages(threshold_metrics(f$fit, f$test, folds = folds)) + by_pairs <- suppressMessages( + threshold_metrics(held_out(f$test$y, f$predicted), folds = folds) + ) + + expect_equal(as.data.frame(by_pairs), as.data.frame(by_model)) + expect_length(attr(by_pairs, "auc"), 5) + expect_s3_class(suppressMessages(plotROC(held_out(f$test$y, f$predicted), + folds = folds)), "ggplot") +}) + +test_that("held_out is not annotated as in-sample unless it says so", { + # Nothing in two numeric vectors records what model made them, so this is + # taken on trust. The flag exists to be set honestly. + f <- held_out_fixture() + + expect_false(attr(threshold_metrics(held_out(f$test$y, f$predicted)), + "in.sample")) + expect_true(attr(threshold_metrics(held_out(f$test$y, f$predicted, + in.sample = TRUE)), + "in.sample")) + # And the caption follows the flag, as it does on the model path. + expect_null(plotROC(held_out(f$test$y, f$predicted))$labels$caption) + expect_false(is.null( + plotROC(held_out(f$test$y, f$predicted, in.sample = TRUE))$labels$caption + )) +}) + +test_that("the three response forms are read the same way as a model's", { + # A two-level factor takes its second level as positive, matching glm(). + f <- held_out_fixture() + numeric <- held_out(f$test$y, f$predicted) + logical <- held_out(f$test$y == 1, f$predicted) + factored <- held_out(factor(ifelse(f$test$y == 1, "yes", "no"), + levels = c("no", "yes")), f$predicted) + + expect_equal(numeric$observed, logical$observed) + expect_equal(numeric$observed, factored$observed) +}) + +test_that("a malformed pair is refused", { + expect_error(held_out(c(0, 1, 1), c(0.2, 0.4)), "same length") + expect_error(held_out(numeric(0), numeric(0)), "nothing to score") + # Predictions on the link scale are the likely mistake, and they are not + # probabilities. + expect_error(held_out(c(0, 1), c(-2.2, 3.1)), "probabilities in \\[0, 1\\]") + expect_error(held_out(c(0, 1, 2), c(0.1, 0.2, 0.3)), "not a binary outcome") + expect_error(held_out(factor(c("a", "b", "c")), c(0.1, 0.2, 0.3)), + "binary outcome only") + expect_error(held_out(c(0, 1), c(0.2, 0.4), in.sample = "yes"), + "must be TRUE or FALSE") +}) + +test_that("one outcome class or a mismatched fold vector is refused", { + expect_error(threshold_metrics(held_out(c(0, 0, 0), c(0.1, 0.2, 0.3))), + "only one outcome class") + expect_error( + threshold_metrics(held_out(c(0, 1, 1), c(0.1, 0.2, 0.3)), folds = c(1, 2)), + "one entry per observation" + ) +}) + +test_that("missing values are dropped as they are on the model path", { + observed <- c(0, 1, NA, 1, 0) + predicted <- c(0.1, 0.9, 0.5, NA, 0.2) + + metrics <- threshold_metrics(held_out(observed, predicted)) + + expect_equal(attr(metrics, "n"), 3) +}) + +test_that("a model still requires newdata, and says how to avoid needing it", { + f <- held_out_fixture() + + expect_error(threshold_metrics(f$fit), "held_out\\(\\)") + expect_error(calibration_estimates(f$fit), "held_out\\(\\)") +}) + +test_that("held_out prints what it holds", { + f <- held_out_fixture() + + expect_output(print(held_out(f$test$y, f$predicted)), "held-out predictions") + expect_output(print(held_out(f$test$y, f$predicted)), "in sample: no") + expect_output(print(held_out(f$test$y, f$predicted, in.sample = TRUE)), + "in sample: yes") +}) From 48c3530bae7d8cc19dc62ab35997d6c6699332c1 Mon Sep 17 00:00:00 2001 From: chross22 <52218551+chross22@users.noreply.github.com> Date: Mon, 10 Aug 2026 16:36:08 -0400 Subject: [PATCH 2/2] Let mess() score a data frame, and name the covariate that made a cell novel Two generalisations, both of which the surface wants on its own terms. A data frame is now a valid x. Plenty of pipelines hold their projection as a table of cells rather than as a raster - anything that predicts over a joined covariate table does - and requiring a round trip through terra to score one was a tax rather than a service. The return follows the input: a raster in, a raster out; a data frame in, a data frame out. limiting = TRUE reports which covariate each cell's score came from. The surface says a cell is novel; this says what made it so, and that is usually the half worth acting on. "This shelf is extrapolated" is a shrug. "It is extrapolated because its chlorophyll is higher than anything in training" is a decision about whether to widen the training window or clip the map. On the data frame side it is a mess_variable column; on the raster side a categorical mess_variable layer, since a raster cannot hold a character and the names have to live in the levels table. Off by default, so the returned shape is unchanged for every existing caller. The check that matters is that the two paths are one method with two doors: the same cells scored as a raster and as a data frame come back identical, and so do the raster's limiting layer and the frame's limiting column. A second implementation for data frames would have been free to drift, and nothing outside would have noticed. Scored against an independent implementation of the same method too - the one in taupatch, written from Elith, Kearney and Phillips (2010) rather than from this code - and the surfaces and the culprit variables agree exactly. That is the more convincing of the two checks, since it does not share a line of code with this one. One test changes rather than being added: test-spatial.R pinned the restriction that x must be a SpatRaster, which is the restriction being lifted. It now asserts what replaced it - a data frame is read, and a list is still refused. Co-Authored-By: Claude Opus 5 --- R/spatial.R | 119 +++++++++++++++++++++---- man/mess.Rd | 29 +++++- man/mess_frame.Rd | 24 +++++ man/mess_vars.Rd | 24 +++++ tests/testthat/test-dismo-successors.R | 84 +++++++++++++++++ tests/testthat/test-spatial.R | 6 +- 6 files changed, 264 insertions(+), 22 deletions(-) create mode 100644 man/mess_frame.Rd create mode 100644 man/mess_vars.Rd diff --git a/R/spatial.R b/R/spatial.R index 8f8cd72..2bcea25 100644 --- a/R/spatial.R +++ b/R/spatial.R @@ -110,7 +110,28 @@ ensemble_summary <- function(x, statistic = c("sd", "cv", "range", "iqr", #' report it as similar. Treat a non-negative surface as the absence of one #' specific problem, not as a licence to project. #' -#' @return A single-layer `SpatRaster` named `mess`. Negative values are novel. +#' @section Which covariate is responsible: +#' The surface says a cell is novel; `limiting = TRUE` says what made it so. +#' That is usually the actionable half -- "this shelf is extrapolated" is a +#' shrug, "extrapolated because its chlorophyll is higher than any training +#' record" is a decision about whether to widen the training window or clip the +#' map. It names the covariate with the lowest similarity, which is the one the +#' minimum was taken from. +#' +#' @section Rasters and data frames: +#' `x` may be a `SpatRaster` of covariate layers or a plain data frame of +#' covariate columns, and the return follows the input. The data frame form is +#' for pipelines that hold their projection as a table of cells rather than as a +#' raster, which is common enough that requiring a round trip through `terra` +#' to score it would be a tax rather than a service. +#' +#' @param limiting Whether to also report the covariate responsible for each +#' cell's score. `FALSE` by default, so the returned shape is unchanged. +#' +#' @return For a `SpatRaster`, a `SpatRaster` named `mess`, gaining a +#' categorical `mess_variable` layer when `limiting = TRUE`. For a data frame, +#' a data frame with a `mess` column and, when `limiting = TRUE`, a +#' `mess_variable` column. Negative values are novel. #' #' @family spatial plots #' @seealso [plotExtrapolation()] to draw it, [ensemble_summary()] for @@ -137,37 +158,97 @@ ensemble_summary <- function(x, statistic = c("sd", "cv", "range", "iqr", #' } #' #' @export -mess <- function(x, training, vars = NULL) { - require_terra() - if (!inherits(x, "SpatRaster")) { - stop("x must be a SpatRaster of covariates, not a <", +mess <- function(x, training, vars = NULL, limiting = FALSE) { + raster <- inherits(x, "SpatRaster") + if (!raster && !is.data.frame(x)) { + stop("x must be a SpatRaster or a data frame of covariates, not a <", paste(class(x), collapse = "/"), ">.", call. = FALSE) } + if (raster) require_terra() training <- training_frame(training) + vars <- mess_vars(x, training, vars, raster) - if (is.null(vars)) vars <- intersect(names(x), names(training)) - if (!length(vars)) { - stop("No covariates in common between the raster (", - paste(names(x), collapse = ", "), ") and the training data (", - paste(names(training), collapse = ", "), ").", call. = FALSE) - } - missing.vars <- setdiff(vars, names(x)) - if (length(missing.vars)) { - stop("Raster has no layer(s): ", paste(missing.vars, collapse = ", "), - call. = FALSE) - } - - layers <- lapply(vars, function(v) { + references <- lapply(vars, function(v) { reference <- stats::na.omit(training[[v]]) if (!length(reference)) { stop("Training data for '", v, "' is entirely missing.", call. = FALSE) } - terra::app(x[[v]], function(p) mess_similarity(p, reference)) + reference + }) + names(references) <- vars + + if (!raster) return(mess_frame(x, references, vars, limiting)) + + layers <- lapply(vars, function(v) { + terra::app(x[[v]], function(p) mess_similarity(p, references[[v]])) }) out <- Reduce(function(a, b) min(a, b), layers) names(out) <- "mess" + if (!limiting) return(out) + + # which.min over the layers, as a categorical layer carrying the names. A + # raster cannot hold a character, so the codes are the levels table. + stacked <- Reduce(c, layers) + worst <- terra::which.min(stacked) + levels(worst) <- data.frame(value = seq_along(vars), mess_variable = vars) + names(worst) <- "mess_variable" + c(out, worst) +} + +#' The covariates a MESS surface can be built from +#' +#' @param x A `SpatRaster` or data frame of covariates. +#' @param training Training data. +#' @param vars Requested covariates, or `NULL` for the ones in common. +#' @param raster Whether `x` is a raster, for the error wording. +#' @return A character vector of covariate names. +#' @keywords internal +mess_vars <- function(x, training, vars, raster) { + available <- names(x) + if (is.null(vars)) vars <- intersect(available, names(training)) + + if (!length(vars)) { + stop("No covariates in common between the ", + if (raster) "raster (" else "data (", + paste(available, collapse = ", "), ") and the training data (", + paste(names(training), collapse = ", "), ").", call. = FALSE) + } + missing.vars <- setdiff(vars, available) + if (length(missing.vars)) { + stop(if (raster) "Raster has no layer(s): " else "Data has no column(s): ", + paste(missing.vars, collapse = ", "), call. = FALSE) + } + vars +} + +#' MESS over a data frame of cells +#' +#' @param x A data frame of covariates. +#' @param references Training values per covariate. +#' @param vars Covariate names. +#' @param limiting Whether to name the covariate responsible. +#' @return A data frame with `mess` and optionally `mess_variable`. +#' @keywords internal +mess_frame <- function(x, references, vars, limiting = FALSE) { + similarity <- vapply(vars, function(v) { + mess_similarity(x[[v]], references[[v]]) + }, numeric(nrow(x))) + + # vapply drops to a plain vector for a single row, which then indexes as if + # it were one column of many. + similarity <- matrix(similarity, nrow = nrow(x), + dimnames = list(NULL, vars)) + + out <- data.frame(mess = apply(similarity, 1, min, na.rm = FALSE)) + if (!limiting) return(out) + + worst <- apply(similarity, 1, function(row) { + if (all(is.na(row))) return(NA_integer_) + which.min(row) + }) + out$mess_variable <- ifelse(is.na(worst), NA_character_, vars[worst]) out } diff --git a/man/mess.Rd b/man/mess.Rd index 9b12b32..61e30b4 100644 --- a/man/mess.Rd +++ b/man/mess.Rd @@ -4,7 +4,7 @@ \alias{mess} \title{Multivariate environmental similarity surface} \usage{ -mess(x, training, vars = NULL) +mess(x, training, vars = NULL, limiting = FALSE) } \arguments{ \item{x}{A \code{SpatRaster} of covariates to project onto. Layer names must @@ -14,9 +14,15 @@ match the columns of \code{training}.} on, or a fitted model to take them from.} \item{vars}{Covariates to consider. Defaults to those common to both.} + +\item{limiting}{Whether to also report the covariate responsible for each +cell's score. \code{FALSE} by default, so the returned shape is unchanged.} } \value{ -A single-layer \code{SpatRaster} named \code{mess}. Negative values are novel. +For a \code{SpatRaster}, a \code{SpatRaster} named \code{mess}, gaining a +categorical \code{mess_variable} layer when \code{limiting = TRUE}. For a data frame, +a data frame with a \code{mess} column and, when \code{limiting = TRUE}, a +\code{mess_variable} column. Negative values are novel. } \description{ Where does a projection leave the conditions the model was fitted under? @@ -39,6 +45,25 @@ separately and still be somewhere the model has never seen, and MESS will report it as similar. Treat a non-negative surface as the absence of one specific problem, not as a licence to project. } +\section{Which covariate is responsible}{ + +The surface says a cell is novel; \code{limiting = TRUE} says what made it so. +That is usually the actionable half -- "this shelf is extrapolated" is a +shrug, "extrapolated because its chlorophyll is higher than any training +record" is a decision about whether to widen the training window or clip the +map. It names the covariate with the lowest similarity, which is the one the +minimum was taken from. +} + +\section{Rasters and data frames}{ + +\code{x} may be a \code{SpatRaster} of covariate layers or a plain data frame of +covariate columns, and the return follows the input. The data frame form is +for pipelines that hold their projection as a table of cells rather than as a +raster, which is common enough that requiring a round trip through \code{terra} +to score it would be a tax rather than a service. +} + \examples{ if (requireNamespace("terra", quietly = TRUE)) { set.seed(1) diff --git a/man/mess_frame.Rd b/man/mess_frame.Rd new file mode 100644 index 0000000..fccddf4 --- /dev/null +++ b/man/mess_frame.Rd @@ -0,0 +1,24 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/spatial.R +\name{mess_frame} +\alias{mess_frame} +\title{MESS over a data frame of cells} +\usage{ +mess_frame(x, references, vars, limiting = FALSE) +} +\arguments{ +\item{x}{A data frame of covariates.} + +\item{references}{Training values per covariate.} + +\item{vars}{Covariate names.} + +\item{limiting}{Whether to name the covariate responsible.} +} +\value{ +A data frame with \code{mess} and optionally \code{mess_variable}. +} +\description{ +MESS over a data frame of cells +} +\keyword{internal} diff --git a/man/mess_vars.Rd b/man/mess_vars.Rd new file mode 100644 index 0000000..b50f2f4 --- /dev/null +++ b/man/mess_vars.Rd @@ -0,0 +1,24 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/spatial.R +\name{mess_vars} +\alias{mess_vars} +\title{The covariates a MESS surface can be built from} +\usage{ +mess_vars(x, training, vars, raster) +} +\arguments{ +\item{x}{A \code{SpatRaster} or data frame of covariates.} + +\item{training}{Training data.} + +\item{vars}{Requested covariates, or \code{NULL} for the ones in common.} + +\item{raster}{Whether \code{x} is a raster, for the error wording.} +} +\value{ +A character vector of covariate names. +} +\description{ +The covariates a MESS surface can be built from +} +\keyword{internal} diff --git a/tests/testthat/test-dismo-successors.R b/tests/testthat/test-dismo-successors.R index ae87432..2308c83 100644 --- a/tests/testthat/test-dismo-successors.R +++ b/tests/testthat/test-dismo-successors.R @@ -319,3 +319,87 @@ test_that("a model that does not keep its data can borrow the plotting data", { expect_s3_class(suppressMessages(plotEffects(fit, d, "x1", n = 10)), "patchwork") }) + +test_that("mess scores a data frame the same as the raster of the same cells", { + skip_if_not_installed("terra") + # The two paths must be one method with two doors. A separate implementation + # for data frames would be free to drift, and nothing outside would notice. + set.seed(1) + training <- data.frame(temp = rnorm(200, 10, 2), depth = runif(200, 0, 100)) + covariates <- c( + terra::rast(nrows = 10, ncols = 10, vals = rnorm(100, 12, 3)), + terra::rast(nrows = 10, ncols = 10, vals = runif(100, -20, 140)) + ) + names(covariates) <- c("temp", "depth") + + by_raster <- mess(covariates, training) + by_frame <- mess(as.data.frame(covariates), training) + + expect_equal(by_frame$mess, as.numeric(terra::values(by_raster)[, 1])) +}) + +test_that("mess names the covariate that made a cell novel", { + training <- data.frame(temp = c(4, 8, 12, 16), depth = c(10, 20, 30, 40)) + # Row 1 is ordinary; row 2 is far too warm; row 3 is far too deep. + cells <- data.frame(temp = c(10, 40, 10), depth = c(25, 25, 400)) + + out <- mess(cells, training, limiting = TRUE) + + expect_equal(names(out), c("mess", "mess_variable")) + expect_gt(out$mess[1], 0) + expect_lt(out$mess[2], 0) + expect_lt(out$mess[3], 0) + expect_equal(out$mess_variable[2], "temp") + expect_equal(out$mess_variable[3], "depth") +}) + +test_that("limiting is off by default, so the returned shape is unchanged", { + training <- data.frame(temp = c(4, 8, 12, 16)) + cells <- data.frame(temp = c(10, 40)) + + expect_equal(names(mess(cells, training)), "mess") + expect_equal(names(mess(cells, training, limiting = TRUE)), + c("mess", "mess_variable")) +}) + +test_that("the raster's limiting layer carries the covariate names", { + skip_if_not_installed("terra") + set.seed(2) + training <- data.frame(temp = rnorm(100, 10, 2), depth = runif(100, 0, 100)) + covariates <- c( + terra::rast(nrows = 8, ncols = 8, vals = rnorm(64, 12, 4)), + terra::rast(nrows = 8, ncols = 8, vals = runif(64, -40, 160)) + ) + names(covariates) <- c("temp", "depth") + + out <- mess(covariates, training, limiting = TRUE) + + expect_equal(names(out), c("mess", "mess_variable")) + # A raster cannot hold a character, so the names live in the levels table. + levels.table <- terra::levels(out[["mess_variable"]])[[1]] + expect_setequal(levels.table$mess_variable, c("temp", "depth")) +}) + +test_that("mess on a data frame handles one row and missing values", { + training <- data.frame(temp = c(4, 8, 12, 16), depth = c(10, 20, 30, 40)) + + one <- mess(data.frame(temp = 10, depth = 25), training, limiting = TRUE) + expect_equal(nrow(one), 1) + expect_equal(names(one), c("mess", "mess_variable")) + + # A covariate missing for a cell makes that cell unscoreable, not an error. + gaps <- mess(data.frame(temp = c(10, NA), depth = c(25, 25)), training, + limiting = TRUE) + expect_false(is.na(gaps$mess[1])) + expect_true(is.na(gaps$mess[2])) +}) + +test_that("mess refuses what it cannot score, and says which side is short", { + training <- data.frame(temp = c(4, 8, 12)) + + expect_error(mess(list(temp = 1), training), "SpatRaster or a data frame") + expect_error(mess(data.frame(salinity = 30), training), + "No covariates in common") + expect_error(mess(data.frame(temp = 10), training, vars = "depth"), + "Data has no column") +}) diff --git a/tests/testthat/test-spatial.R b/tests/testthat/test-spatial.R index d84a03b..d2bd7e0 100644 --- a/tests/testthat/test-spatial.R +++ b/tests/testthat/test-spatial.R @@ -170,7 +170,11 @@ test_that("mess reports what it cannot match", { "No covariates in common") expect_error(mess(f$covariates, f$training, vars = "salinity"), "no layer\\(s\\): salinity") - expect_error(mess(data.frame(a = 1), f$training), "must be a SpatRaster") + # A data frame is a supported input now, so the complaint is about the + # covariates in it rather than about its class. + expect_error(mess(data.frame(a = 1), f$training), "No covariates in common") + # Something that is neither is still refused. + expect_error(mess(list(a = 1), f$training), "SpatRaster or a data frame") }) test_that("a constant covariate does not divide by a zero range", {