sdmTMB_cv: Cross validation with sdmTMB models

In sdmTMB: Spatial and Spatiotemporal SPDE-Based GLMMs with 'TMB'

Cross validation with sdmTMB models

Description

Performs k-fold or leave-future-out cross validation with sdmTMB models. Returns the sum of log likelihoods of held-out data (log predictive density), which can be used to compare models—higher values indicate better out-of-sample prediction. To judge whether differences between models are meaningful, pass the output to loo::elpd() and compare with loo::loo_compare(), which reports the difference and its standard error (see elpd.sdmTMB_cv()). Models must use the same folds (fold_ids) and the same predictive argument. By default, creates folds randomly and stratified by time (set a seed for reproducibility), but folds can be manually assigned via fold_ids. See Ward and Anderson (2025) in the References and the cross-validation vignette.

Usage

sdmTMB_cv(
  formula,
  data,
  mesh_args,
  mesh = NULL,
  time = NULL,
  k_folds = 8,
  fold_ids = NULL,
  lfo = FALSE,
  lfo_forecast = 1,
  lfo_validations = 5,
  parallel = TRUE,
  use_initial_fit = FALSE,
  save_models = TRUE,
  future_globals = NULL,
  predictive = c("mle-mvn", "mle-eb", "joint"),
  nsim = 100L,
  ...
)

Arguments

formula	Model formula.
data	A data frame. Rows sdmTMB() omits for missing values are left out of the validation results.
mesh_args	Arguments for make_mesh(). If supplied, the mesh will be reconstructed for each fold.
mesh	Output from make_mesh(). If supplied, the same mesh will be used for all folds. This is faster and usually what you want.
time	The name of the time column. Leave as NULL if this is only spatial data.
k_folds	Number of folds.
fold_ids	Optional vector of user-supplied fold IDs. Can also be a single string, e.g. "fold_id", naming a variable in data. Ignored if lfo = TRUE.
lfo	Logical. Use leave-future-out (LFO) cross validation? If TRUE, data from earlier time steps are used to predict future time steps. The time argument must be specified. See Details section below.
lfo_forecast	If lfo = TRUE, number of time steps ahead to forecast. For example, lfo_forecast = 1 means fitting to time steps 1 to T and validating on T + 1. See Details section below.
lfo_validations	If lfo = TRUE, number of times to step through the LFO process (i.e., number of validation folds). Defaults to 5. See Details section below.
parallel	If TRUE and a future::plan() is supplied, will be run in parallel.
use_initial_fit	Fit the first fold and use those parameter values as starting values for subsequent folds? Can be faster with many folds.
save_models	Logical. If TRUE (default), the fitted model object for each fold is stored in the output. If FALSE, models are not saved, which can substantially reduce memory usage for large datasets or many folds. When FALSE, functions that require access to the fitted models (e.g., tidy(), cv_to_waywiser()) will not work. To save space, the fold models are fit without the joint precision matrix unless control is supplied (see get_joint_precision in sdmTMBcontrol()); predict(..., nsim = ...) recomputes it when needed.
future_globals	A character vector of global variables used within arguments if an error is returned that future.apply can't find an object. This vector is appended to TRUE and passed to the future.globals in future.apply::future_lapply(). Useful if global objects are used to specify arguments such as priors or families.
predictive	Which predictive density to score held-out data with. "mle-mvn" (default) integrates over the random effects' approximate (Laplace) posterior with fixed effects at their maximum likelihood estimates (MLEs). "mle-eb" (the previous default) evaluates the likelihood with fixed effects at their MLEs and random effects at their empirical Bayes (EB) estimates. These match type in simulate.sdmTMB() and residuals.sdmTMB(). "joint" integrates over the joint approximate posterior of all parameters, as with nsim in predict.sdmTMB(). "mle-eb" ignores uncertainty in the random effects and can be overconfident when random fields are flexible (e.g., fine meshes) or forecast (LFOCV). For models without random effects, "mle-mvn" is identical to "mle-eb". See Details.
nsim	Number of posterior draws per fold if predictive is "mle-mvn" or "joint". The Monte Carlo log-mean-exp estimate is biased low by O(1/\code{nsim}); check that results are stable to increasing nsim.
...	All other arguments required to run the sdmTMB() model. The weights argument is supported and will be combined with the internal fold-assignment mechanism (held-out data are assigned weight 0).

Details

Plug-in vs. integrated predictive densities

With predictive = "mle-eb", held-out observations are scored at the conditional mode of the random effects, \log p(y \mid \hat{u}, \hat{\theta}), which treats the random fields as known (a plug-in predictive density). With "mle-mvn", the score is \log E[p(y \mid u, \hat{\theta})] over u \sim N(\hat{u}, H^{-1}), where H is the inner (Laplace) Hessian of the training fit. With "joint", fixed effects and variance parameters (on their unconstrained scale) are also drawn from the joint precision matrix. Each is averaged over nsim draws per fold separately for each held-out observation, so fold_loglik is \sum_i \log E[p(y_i \mid u)]. If reml = TRUE, fixed effects are treated as random effects and are therefore also drawn under "mle-mvn". For Gaussian observations, "mle-mvn" is exact given \hat{\theta}; for other families it relies on the Laplace approximation's Gaussian approximation to the random effects' posterior.

Parallel processing

Parallel processing can be used by setting a future::plan().

For example:

library(future)
plan(multisession)
# now use sdmTMB_cv() ...

Leave-future-out cross validation (LFOCV)

An example of LFOCV with 9 time steps, lfo_forecast = 1, and lfo_validations = 2:

• Fit data to time steps 1 to 7, predict and validate step 8.

• Fit data to time steps 1 to 8, predict and validate step 9.

An example of LFOCV with 9 time steps, lfo_forecast = 2, and lfo_validations = 3:

• Fit data to time steps 1 to 5, predict and validate step 7.

• Fit data to time steps 1 to 6, predict and validate step 8.

• Fit data to time steps 1 to 7, predict and validate step 9.

Note that these are time steps in the order they appear in the data. For example, in the pcod data example below, steps between observations are not always one year, so lfo_forecast = 2 forecasts two observed time steps ahead, not necessarily two calendar years.

See example below.

Value

A list:

• data: Original data plus columns for fold ID (cv_fold), CV predicted value (cv_predicted), CV log likelihood (cv_loglik), "mle-eb" CV log likelihood (cv_loglik_mle_eb; identical to cv_loglik if predictive = "mle-eb"), and CV deviance residuals (cv_deviance_resid).

• models: A list of fitted models, one per fold. NULL if save_models = FALSE.

• fold_loglik: Sum of cv_loglik per fold (sum of pointwise log predictive densities of the held-out data). More positive values indicate better out-of-sample prediction.

• sum_loglik: Sum of fold_loglik across all folds (total pointwise log predictive density). Use this to compare models; larger values are better.

• sum_loglik_mle_eb: Sum of cv_loglik_mle_eb.

• predictive: The predictive argument used for cv_loglik.

• pdHess: Logical vector: was the Hessian positive definite for each fold?

• converged: Logical: did all folds converge (all pdHess TRUE)?

• max_gradients: Maximum absolute gradient for each fold.

See Also

elpd.sdmTMB_cv(), loo::loo_compare(), loo::elpd()

Examples

mesh <- make_mesh(pcod, c("X", "Y"), cutoff = 25)

# Set parallel processing first if desired with the future package.
# See the Details section above.

m_cv <- sdmTMB_cv(
  density ~ 0 + depth_scaled + depth_scaled2,
  data = pcod, mesh = mesh, spatial = "off",
  family = tweedie(link = "log"), k_folds = 2
)

m_cv$fold_loglik
m_cv$sum_loglik

head(m_cv$data)
m_cv$models[[1]]
m_cv$max_gradients



# Create mesh each fold:
m_cv2 <- sdmTMB_cv(
  density ~ 0 + depth_scaled + depth_scaled2,
  data = pcod, mesh_args = list(xy_cols = c("X", "Y"), cutoff = 20),
  family = tweedie(link = "log"), k_folds = 2
)

# Use fold_ids:
m_cv3 <- sdmTMB_cv(
  density ~ 0 + depth_scaled + depth_scaled2,
  data = pcod, mesh = mesh,
  family = tweedie(link = "log"),
  fold_ids = rep(seq(1, 3), nrow(pcod))[seq(1, nrow(pcod))]
)

# LFOCV:
m_lfocv <- sdmTMB_cv(
  present ~ s(year, k = 4),
  data = pcod,
  lfo = TRUE,
  lfo_forecast = 2,
  lfo_validations = 3,
  family = binomial(),
  mesh = mesh,
  spatial = "off", # fast example
  spatiotemporal = "off", # fast example
  time = "year" # must be specified
)

# See how the LFOCV folds were assigned:
fold_table <- table(m_lfocv$data$cv_fold, m_lfocv$data$year)
fold_table

# Compare models fit to the same folds with the loo package:
m_cv3_null <- sdmTMB_cv(
  density ~ 1,
  data = pcod, mesh = mesh,
  family = tweedie(link = "log"),
  fold_ids = m_cv3$data$cv_fold
)
if (requireNamespace("loo", quietly = TRUE)) {
  loo::loo_compare(list(
    depth = loo::elpd(m_cv3),
    null = loo::elpd(m_cv3_null)
  ))
}

sdmTMB documentation built on Oct. 17, 2026, 1:07 a.m.