Does cross validation (sdmTMB_cv) work with delta models?

[ericward-noaa]

Question via email: does cross validation (sdmTMB_cv) work with delta models?

[ericward-noaa]

The delta model features are relatively new, and compatibility with sdmTMB_cv() has not yet been added. That said, if you have a large dataset, you can manually do a test/train split or cross validation. For example, with the pcod dataset, we can holdout 10% of the observations as a test set and fit a delta-Gamma model to cod densities

data(pcod)
w <- ifelse(runif(nrow(pcod)) < 0.9, 1, 0) # weights, 10% in test set
pcod_spde <- make_mesh(pcod, c("X", "Y"), cutoff = 15)
fit <- sdmTMB(density ~ 1,
              data = pcod, mesh = pcod_spde, spatial = "on",
              spatiotemporal = list("off", "off"), time = "year",
              family = delta_gamma(),
              weights = w # notice weights argument
)

Next, we have to make predictions for the binomial component -- and calculate the likelihood

pred <- predict(fit, newdata = pcod[which(w==1),])
bin_like <- dbinom(pred$present, size=1, prob=plogis(pred$est1), log=TRUE)

The gamma likelihood is a little trickier,

s1 = exp(fit$sd_report$par.fixed[["ln_phi"]])# shape
s2 = pred$est2 / s1 # scale
gamma_like <- dgamma(pred$density, shape = s1, scale = s2, log = TRUE)

The total likelihood is just the sum of the 2 components,
tot_like <- sum(bin_like) + sum(gamma_like[which(is.finite(gamma_like))])

The 'tot_like' is the total out of sample predictive density, and it's something you want to maximize via model selection. Another issue to be worked out is how to best select the training data. The above approach is totally random, which might be fine. But there's some recent work (see blockCV or spatialsample R packages) that can do spatial blocking -- to reduce autocorrelation between test / training data. This is a very active area of current research.

[Raquel-RuizDiaz]

Good morning,

I was wondering if doing cross-validation with sdmTMB_cv is still not compatible with delta-gamma model. if so, could I manually do a test/train "leaving the future out" (i.e. leaving the last 3 years of my dataset out and seen how good the model predicts these years) and then using the script you shared here to calculate likelihood of the predictions?

Thanks for your help.
Raquel Ruiz.

[ericward-noaa]

Hi -- this is still not implemented for delta models. I think maybe an even easier approach compared to the above script would just be to use the sdmTMB_cv() function for each of the 2 model components separately. We have added some arguments that you'll want to pay attention to:
lfo - you'll want to make this TRUE
lfo_forecast - how many time steps to forecast (defaults to 1)
lfo_validations - how many validation years to use, in your case 3

[seananderson]

Cross validation is now working with delta models:

library(sdmTMB)
future::plan(future::multisession)
out_dg <- sdmTMB_cv(
  density ~ depth_scaled,
  data = pcod_2011, 
  mesh = make_mesh(pcod_2011, c("X", "Y"), cutoff = 10),
  family = delta_gamma(), 
  k_folds = 5
)
#> Running fits with `future.apply()`.
#> Set a parallel `future::plan()` to use parallel processing.
out_dg$fold_loglik
#> [1] -531.3995 -702.3048 -591.7539 -608.1418 -611.2388
out_dg$sum_loglik
#> [1] -3044.839