add CAIC and EDF calcs

NAMESPACE

@@ -1,5 +1,6 @@
 # Generated by roxygen2: do not edit by hand
 
+S3method(cAIC,sdmTMB)
 S3method(coef,sdmTMB)
 S3method(confint,sdmTMB)
 S3method(df.residual,sdmTMB)
@@ -25,6 +26,7 @@ S3method(vcov,sdmTMB)
 export(Beta)
 export(add_barrier_mesh)
 export(add_utm_columns)
+export(cAIC)
 export(censored_poisson)
 export(delta_beta)
 export(delta_gamma)

R/caic.R

@@ -0,0 +1,129 @@
+#' Calculate conditional AIC
+#'
+#' Calculates the conditional Akaike Information criterion (cAIC).
+#'
+#' @param object Output from [sdmTMB()].
+#' @param what Whether to return the cAIC or the effective degrees of freedom
+#'   (EDF) for each group of random effects.
+#' @param ... Other arguments for specific methods. Not used.
+#'
+#' @details cAIC is designed to optimize the expected out-of-sample predictive
+#' performance for new data that share the same random effects as the in-sample
+#' (fitted) data, e.g., spatial interpolation.  In this sense, it should be a
+#' fast approximation to optimizing the model structure based on k-fold
+#' cross-validation.
+#'
+#' By contrast, [AIC()] calculates the marginal Akaike Information Criterion,
+#' which is designed to optimize expected predictive performance for new data
+#' that have new random effects, e.g., extrapolation, or inference about
+#' generative parameters.
+#'
+#' cAIC also calculates the effective degrees of freedom (EDF) as a byproduct.
+#' This is the number of fixed effects that would have an equivalent impact on
+#' model flexibility as a given random effect.
+#'
+#' Both cAIC and EDF are calculated using Eq. 6 of Zheng, Cadigan, and Thorson
+#' (2024).
+#'
+#' For models that include profiled fixed effects, these profiles are turned
+#' off.
+#'
+#' @return
+#' Either the cAIC or the effective degrees of freedom (EDF) by group
+#' of random effects depending on the argument `what`.
+#'
+#' @references
+#' **Deriving the general approximation to cAIC used here:**
+#'
+#' Zheng, N., Cadigan, N., & Thorson, J. T. (2024).
+#' A note on numerical evaluation of conditional Akaike information for
+#' nonlinear mixed-effects models (arXiv:2411.14185). arXiv.
+#' \doi{10.48550/arXiv.2411.14185}
+#'
+#' **The utility of EDF to diagnose hierarchical model behaviour:**
+#'
+#' Thorson, J. T. (2024). Measuring complexity for hierarchical
+#' models using effective degrees of freedom. Ecology,
+#' 105(7), e4327 \doi{10.1002/ecy.4327}
+#'
+#' @examples
+#' mesh <- make_mesh(dogfish, c("X", "Y"), cutoff = 15)
+#' fit <- sdmTMB(catch_weight ~ s(log(depth)),
+#'   time_varying = ~1,
+#'   time_varying_type = "ar1",
+#'   time = "year",
+#'   spatiotemporal = "off",
+#'   mesh = mesh,
+#'   family = tweedie(),
+#'   data = dogfish,
+#'   offset = log(dogfish$area_swept)
+#' )
+#' cAIC(fit)
+#' cAIC(fit, what = "EDF")
+#' AIC(fit)
+#' @export
+cAIC <- function(object, what = c("cAIC", "EDF"), ...) {
+  UseMethod("cAIC", object)
+}
+
+#' @exportS3Method
+cAIC.sdmTMB <- function(object, what = c("cAIC", "EDF")) {
+  what <- match.arg(what)
+  what <- tolower(what)
+  tmb_data <- object$tmb_data
+
+  ## Ensure profile = NULL
+  if (is.null(object$control$profile)) {
+    obj <- object$tmb_obj
+  } else {
+    obj <- TMB::MakeADFun(
+      data = tmb_data,
+      parameters = object$parlist,
+      map = object$tmb_map,
+      random = object$tmb_random,
+      DLL = "sdmTMB",
+      profile = NULL #<
+    )
+  }
+
+  ## Make obj_new
+  tmb_data$weights_i[] <- 0
+  obj_new <- TMB::MakeADFun(
+    data = tmb_data,
+    parameters = object$parlist,
+    map = object$tmb_map,
+    random = object$tmb_random,
+    DLL = "sdmTMB",
+    profile = NULL
+  )
+
+  par <- obj$env$parList()
+  parDataMode <- obj$env$last.par
+  indx <- obj$env$lrandom()
+  q <- length(indx)
+  p <- length(object$model$par)
+
+  ## use '-' for Hess because model returns negative loglikelihood
+  Hess_new <- -Matrix::Matrix(obj_new$env$f(parDataMode, order = 1, type = "ADGrad"), sparse = TRUE)
+  Hess_new <- Hess_new[indx, indx] ## marginal precision matrix of REs
+
+  ## Joint hessian etc
+  Hess <- -Matrix::Matrix(obj$env$f(parDataMode, order = 1, type = "ADGrad"), sparse = TRUE)
+  Hess <- Hess[indx, indx]
+  negEDF <- Matrix::diag(Matrix::solve(Hess, Hess_new, sparse = FALSE))
+
+  if (what == "caic") {
+    jnll <- obj$env$f(parDataMode)
+    cnll <- jnll - obj_new$env$f(parDataMode)
+    cAIC <- 2 * cnll + 2 * (p + q) - 2 * sum(negEDF)
+    return(cAIC)
+  } else if (what == "edf") {
+    ## Figure out group for each random-effect coefficient
+    group <- factor(names(object$last.par.best[obj$env$random]))
+    ## Calculate total EDF by group
+    EDF <- tapply(negEDF, INDEX = group, FUN = length) - tapply(negEDF, INDEX = group, FUN = sum)
+    return(EDF)
+  } else {
+    cli_abort("Option not implemented")
+  }
+}

scratch/caic-demo.R

@@ -0,0 +1,16 @@
+
+library(sdmTMB)
+
+# Build a mesh to implement the SPDE approach:
+mesh <- make_mesh(pcod_2011, c("X", "Y"), cutoff = 20)
+
+# Fit a Tweedie spatial random field GLMM with a smoother for depth:
+fit <- sdmTMB(
+  density ~ s(depth),
+  data = pcod_2011, mesh = mesh,
+  family = tweedie(link = "log"),
+  control = sdmTMBcontrol(profile="b_j")
+)
+
+CAIC.sdmTMB(fit, what="CAIC")
+CAIC.sdmTMB(fit, what="EDF")

tests/testthat/test-cAIC.R

@@ -0,0 +1,21 @@
+test_that("cAIC and EDF work", {
+  skip_on_cran()
+  skip_on_ci()
+  mesh <- make_mesh(dogfish, c("X", "Y"), cutoff = 15)
+  suppressMessages(
+    fit <- sdmTMB(catch_weight ~ s(log(depth)),
+      time_varying = ~1,
+      time_varying_type = "ar1",
+      time = "year",
+      spatiotemporal = "off",
+      mesh = mesh,
+      family = tweedie(),
+      data = dogfish,
+      offset = log(dogfish$area_swept)
+    )
+  )
+  expect_equal(AIC(fit), 12192.9613, tolerance = 1e-1)
+  expect_equal(cAIC(fit), 12089.4289, tolerance = 1e-1)
+  edf <- cAIC(fit, what = "EDF")
+  expect_equal(sum(edf), 54.3870, tolerance = 1e-2)
+})