Cross-validated comparison of spatial models

Fits and cross-validates one or more model types, returning a unified comparison table. Unlike compare_models(), this function does perform fitting (inside CV folds), because CV inherently requires repeated fitting.

Usage

compare_models_cv(
  data_sf,
  response_var,
  predictor_vars,
  models = c("GWR", "Bayesian"),
  k = 5,
  seed = 123,
  folds = NULL,
  boundary = NULL,
  pointize = "auto",
  gwr_args = list(),
  bayes_args = list(),
  rf_args = list(),
  summary = c("mean", "median"),
  quiet = FALSE,
  block_size = NULL,
  auto_range = FALSE,
  metrics = NULL
)

Arguments

data_sf: An sf object.
response_var: Response column name.
predictor_vars: Predictor column names.
models: Character vector: any subset of c("GWR", "Bayesian", "RF"), in any order. Each is cross-validated on the same folds. Names outside that set raise a warning and are dropped; if nothing recognised remains, this is an error. There is no silent fallback. A recognised model whose backend package is not installed is dropped with a message so the call still returns the models that could run. But if none of the requested backends is installed, nothing is left to compare and the call errors with "no viable models.". Guard with requireNamespace() when the model set is not known in advance.
k: Number of folds. Default 5.
seed: RNG seed. Default 123.
folds: Optional fold definitions: a make_folds() return value, or a bare list of list(train =, test =) pairs of ..row_id values. Train and test must be disjoint (a fold that trains on its own test rows is not a cross-validation split and is refused with an error), and IDs naming no row in the prepared data are dropped with a logged count (expected when rows were removed for missing values; a sign the folds came from other data when they were not).
boundary: Optional polygon sf/sfc.
pointize: Geometry coercion strategy. It also decides where a polygon or line row falls in the shared blocks, so each row is placed by the point every model is fitted at.
gwr_args: Extra arguments for cv_gwr. Only names that are formal arguments of cv_gwr() are forwarded (it has no ...), so entries meant for fit_gwr_model() alone (e.g. longlat) cannot be passed this way. Anything dropped is named in a warning; call cv_gwr() directly if you need it.
bayes_args: Extra arguments for fit_bayesian_spatial_model(). Forwarded whole as cv_bayes(fit_args = ), so an unrecognised name raises an error from fit_bayesian_spatial_model() and is not dropped silently. compute_loo, boundary and pointize are overridden by the CV internals.
rf_args: Extra arguments for cv_rf, which passes anything it does not recognise on to fit_rf_model and thence to ranger::ranger().
summary: "mean" or "median" for Bayesian predictions.
quiet: Logical; suppress this function's progress message()s. It does not silence R warnings, nor the package's console log echo (see spatialkit_quiet for that). Default FALSE.
block_size: Optional minimum block edge length for the shared spatial CV blocks (projected CRS units), passed to make_folds() when folds is NULL. Default NULL.
auto_range: Logical. If TRUE and folds is NULL, the autocorrelation range of the response (detrended on predictor_vars) is estimated and used as the minimum block size of the shared folds, as in make_folds(). Default FALSE: geometric blocks, as before this argument existed. Either way the fold set is built once and every backend is scored on it. When it cannot be built (a block_size or estimated range that leaves a single block, say) the call is an error, as it is for each backend on its own; no model is scored on a design other than the one asked for.
metrics: Optional scoring function of your own, handed to every backend's cv_*(): a function(y, yhat) returning a named numeric vector, applied per fold and to each backend's pooled predictions, whose names become columns of by_fold and overall beside the built-in ones. See Your own metrics on cv_spatial() for the contract. Because the backends are scored on the same folds and, in overall, on the same rows (see Value), the columns are comparable across rows of overall.

Value

A list with overall, by_fold, and per-model cv_results (gwr_cv, bayes_cv, rf_cv for the models that ran). overall has one row per model with the pooled metrics, n_pred (the rows they are computed on), the coverage and CRPS columns described above when a Bayesian model ran, and model as its last column. Shared folds do not guarantee shared rows: a model that fails on a fold (GWR with a fixed bandwidth across a gap in the data, say) or predicts NA for some rows pools fewer rows, usually without the hardest ones. When the models predicted different rows, the function warns and recomputes every model's pooled metrics, your own metrics included, on the rows all of them predicted, so n_pred is the same on every row that has predictions; a model that predicted nothing stays an NA row. Each model's metrics over all the rows it predicted stay in its *_cv element and in attr(overall, "all_rows"), a table of the same shape. by_fold and the Bayesian coverage and CRPS columns are per fold and are not recomputed. Only the models that actually ran appear, so check which names are present; there is not always one entry per requested model, because a backend whose package is missing is dropped with a message. When no requested backend is available there is nothing to return and the function errors with "no viable models." instead of returning an empty comparison.

Coverage and CRPS in the overall table

A model can predict well on average and still be badly calibrated, so a comparison read from RMSE alone can prefer the model whose uncertainty is wrong. Heaton et al. (2019) found good point prediction routinely alongside poor interval coverage. When "Bayesian" is among the models that ran, overall therefore also carries the columns of cv_bayes()$predictive_coverage: coverage_50, coverage_80, coverage_95 (the share of held-out rows inside the posterior predictive interval at each level, averaged across folds weighted by each fold's n_pred) and mean_CRPS (the continuous ranked probability score, lower is better). They are NA on the GWR and RF rows, which produce point predictions and no draws, and absent when no Bayesian model ran. Read coverage against its nominal level: 0.95 at coverage_95 is calibrated, well below it is overconfident, well above it is wider than it needs to be.

Percentage errors on responses with zeros

MAPE divides by the observed value and SMAPE by $|y| + |\hat{y}|$, so neither is defined where its denominator is zero. Neither returns Inf or NaN. Both are averaged over the rows whose denominator is non-zero, and are NA when no row qualifies. Non-zero is judged at the scale of the data: a denominator no larger than 100 machine epsilons times the largest one counts as zero, so the rule does not depend on the units of the response. The n_MAPE and n_SMAPE columns record how many rows that was; the n column counts finite observation/prediction pairs. Read a percentage error next to its count: when n_MAPE < n, MAPE is an average over a subset of the data, whatever its value.

This bites on any response taking exact zeros: counts, rainfall, abundance, claim amounts. On a zero-inflated response with 62 zeros out of 120, MAPE is an average over the 58 non-zero rows, which n_MAPE = 58 now says. SMAPE fails differently and more subtly: it drops the rows where observation and prediction are both near zero (which on a well-fitted zero-inflated model are the rows it got right), so it averages the harder rows only and reads worse than the fit deserves; n_SMAPE shows how many rows it kept, and the count is only a label, not a repair.

RMSE, MAE and $R^2$ use every finite row and are unaffected; prefer them whenever the response can be zero. For a Bayesian fit, cv_bayes() additionally reports CRPS and interval coverage, which are proper scoring rules and have no such failure mode.

Which metrics survive a non-Gaussian response

RMSE and MAE are defined for any numeric response and are what to read for a count, a rate or a bounded outcome. MAPE and SMAPE assume a response that is rarely zero (see the previous section), and R-squared and adjusted R-squared compare residual variance to total variance, which is the right comparison for a Gaussian response and a loose one for anything whose variance tracks its mean. None of the four is wrong to compute; each is Gaussian-shaped thinking, and on a Poisson or zero-inflated response should be read as a rough summary rather than a score.

For the Bayesian backend, cv_bayes() additionally reports CRPS and interval coverage at 50, 80 and 95 percent. Both are proper scoring rules computed from posterior draws, so they are meaningful for any family that predicts one number per row (a count, a rate, a binary or bounded outcome), and they are the numbers to compare when the response is not Gaussian. A categorical or ordinal family predicts a probability per response category instead, so cv_bayes() refuses one before fitting anything. When every fold fails, the fold_metrics frame cv_bayes() returns carries the CRPS column but not the coverage_* columns, so code that reads those columns must tolerate their absence.

References

Heaton, M. J., Datta, A., Finley, A. O., Furrer, R., Guinness, J., Guhaniyogi, R., Gerber, F., Gramacy, R. B., Hammerling, D., Katzfuss, M., Lindgren, F., Nychka, D. W., Sun, F. and Zammit-Mangion, A. (2019). A case study competition among methods for analyzing large spatial data. Journal of Agricultural, Biological and Environmental Statistics, 24(3), 398–425. doi:10.1007/s13253-018-00348-w

Examples

if (requireNamespace("ranger", quietly = TRUE)) {
  library(sf)
  set.seed(1)
  n <- 120
  dat <- st_as_sf(
    data.frame(x = 5e5 + runif(n, 0, 1000), y = 5e6 + runif(n, 0, 1000),
               elev = rnorm(n)),
    coords = c("x", "y"), crs = 32632
  )
  dat$price <- 10 + 0.01 * (st_coordinates(dat)[, 1] - 5e5) +
    2 * dat$elev + rnorm(n)
  # The random forest, and GWR beside it when GWmodel is installed, on the
  # same folds.  "Bayesian" is left out here because a Stan fit per fold
  # takes minutes; add it to `models` for a run to report.
  models <- c("RF", if (requireNamespace("GWmodel", quietly = TRUE) &&
                       requireNamespace("sp", quietly = TRUE)) "GWR")
  cmp <- compare_models_cv(dat, "price", "elev", models = models, k = 3,
                           rf_args = list(num_trees = 100),
                           gwr_args = list(bandwidth = 30))
  cmp$overall
}
#> compare_models_cv(): running CV for GWR ...
#> compare_models_cv(): running CV for RF ...
#>       RMSE      MAE     MAPE    SMAPE        R2 Adj_R2 n_pred n_MAPE n_SMAPE
#> 1 2.281133 1.853980 12.79336 12.53050 0.6738657     NA    120    120     120
#> 2 3.711068 3.128498 21.71902 20.99972 0.1368361     NA    120    120     120
#>   model
#> 1   GWR
#> 2    RF