Cross validation with sdmTMB models

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).

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.

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.

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
)
#> Running fits with `future.apply()`.
#> Set a parallel `future::plan()` to use parallel processing.

m_cv$fold_loglik
#> [1] -4462.402 -4176.627
m_cv$sum_loglik
#> [1] -8639.029

head(m_cv$data)
#> # A tibble: 6 × 17
#>    year     X     Y depth density present   lat   lon depth_mean depth_sd
#>   <int> <dbl> <dbl> <dbl>   <dbl>   <dbl> <dbl> <dbl>      <dbl>    <dbl>
#> 1  2003  446. 5793.   201   113.        1  52.3 -130.       5.16    0.445
#> 2  2003  446. 5800.   212    41.7       1  52.3 -130.       5.16    0.445
#> 3  2003  449. 5802.   220     0         0  52.4 -130.       5.16    0.445
#> 4  2003  437. 5802.   197    15.7       1  52.4 -130.       5.16    0.445
#> 5  2003  421. 5771.   256     0         0  52.1 -130.       5.16    0.445
#> 6  2003  418. 5772.   293     0         0  52.1 -130.       5.16    0.445
#> # ℹ 7 more variables: depth_scaled <dbl>, depth_scaled2 <dbl>, cv_fold <int>,
#> #   cv_predicted <dbl>, cv_deviance_resid <dbl>, cv_loglik_mle_eb <dbl>,
#> #   cv_loglik <dbl>
m_cv$models[[1]]
#> Model fit by ML ['sdmTMB']
#> Formula: density ~ 0 + depth_scaled + depth_scaled2
#> Mesh: sdmTMBmesh (isotropic covariance)
#> Data: tbl_df
#> Family: tweedie(link = 'log')
#>  
#> Conditional model:
#>               coef.est coef.se
#> depth_scaled     -0.95    0.14
#> depth_scaled2     1.55    0.17
#> 
#> Dispersion parameter: 61.99
#> Tweedie p: 1.85
#> ML criterion at convergence: 4131.555
#> 
#> See ?tidy.sdmTMB to extract these values as a data frame.
m_cv$max_gradients
#> [1] 7.323093e-12 1.741236e-10


# \donttest{
# 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
)
#> Running fits with `future.apply()`.
#> Set a parallel `future::plan()` to use parallel processing.

# 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))]
)
#> Running fits with `future.apply()`.
#> Set a parallel `future::plan()` to use parallel processing.

# 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
)
#> Running fits with `future.apply()`.
#> Set a parallel `future::plan()` to use parallel processing.

# See how the LFOCV folds were assigned:
fold_table <- table(m_lfocv$data$cv_fold, m_lfocv$data$year)
fold_table
#>    
#>     2013 2015 2017
#>   3  240    0    0
#>   4    0  238    0
#>   5    0    0  240

# 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
)
#> Running fits with `future.apply()`.
#> Set a parallel `future::plan()` to use parallel processing.
if (requireNamespace("loo", quietly = TRUE)) {
  loo::loo_compare(list(
    depth = loo::elpd(m_cv3),
    null = loo::elpd(m_cv3_null)
  ))
}
#>  model elpd_diff se_diff p_worse diag_diff diag_elpd
#>  depth       0.0     0.0      NA                    
#>   null    -199.1    33.2    1.00                    
# }

Usage

Arguments

Value

Details

See also

Examples