sdmTMB_cv: Cross validation with sdmTMB models

View source: R/cross-val.R

sdmTMB_cvR Documentation

Cross validation with sdmTMB models

Description

Facilitates cross validation with sdmTMB models. Returns the log likelihood of left-out data, which is similar in spirit to the ELPD (expected log pointwise predictive density). The function has an option for leave-future-out cross validation. By default, the function creates folds randomly but folds can be manually assigned via the fold_ids argument.

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,
  future_globals = NULL,
  spde = deprecated(),
  ...
)

Arguments

formula

Model formula.

data

A data frame.

mesh_args

Arguments for make_mesh(). If supplied, the mesh will be reconstructed for each fold.

mesh

Output from make_mesh(). If supplied, the mesh will be constant across folds.

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 containing user fold IDs. Can also be a single string, e.g. "fold_id" representing the name of the variable in data. Ignored if lfo is TRUE

lfo

Whether to implement leave-future-out (LFO) cross validation where data are used to predict future folds. time argument in sdmTMB() must be specified. See Details section below.

lfo_forecast

If lfo = TRUE, number of time steps to forecast. Time steps 1, ..., T are used to predict T + lfo_forecast and the last forecasted time step is used for validation. See Details section below.

lfo_validations

If lfo = TRUE, number of times to step through the LFOCV process. 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.

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 argument future.globals in future.apply::future_lapply(). Useful if global objects are used to specify arguments like priors, families, etc.

spde

Depreciated. Use mesh instead.

...

All other arguments required to run sdmTMB() model with the exception of weights, which are used to define the folds.

Details

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.

See example below.

Value

A list:

  • data: Original data plus columns for fold ID, CV predicted value, and CV log likelihood.

  • models: A list of models; one per fold.

  • fold_loglik: Sum of left-out log likelihoods per fold. More positive values are better.

  • sum_loglik: Sum of fold_loglik across all left-out data. More positive values are better.

  • pdHess: Logical vector: Hessian was invertible each fold?

  • converged: Logical: all pdHess TRUE?

  • max_gradients: Max gradient per fold.

Prior to sdmTMB version '0.3.0.9002', elpd was incorrectly returned as the log average likelihood, which is another metric you could compare models with, but not ELPD. For maximum likelihood, ELPD is equivalent in spirit to the sum of the log likelihoods.

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,
  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,
  mesh = mesh,
  lfo = TRUE,
  lfo_forecast = 2,
  lfo_validations = 3,
  family = binomial(),
  spatiotemporal = "off",
  time = "year" # must be specified
)

# See how the LFOCV folds were assigned:
example_data <- m_lfocv$models[[1]]$data
table(example_data$cv_fold, example_data$year)


pbs-assess/sdmTMB documentation built on March 25, 2024, 11:12 p.m.