cAIC refinements, see notes

- style code
- copy edit docs
- define cAIC generic
- CAIC -> cAIC
- add example
- add unit test
- don't require() Matrix

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)
@@ -23,9 +24,9 @@ S3method(terms,sdmTMB)
 S3method(tidy,sdmTMB)
 S3method(vcov,sdmTMB)
 export(Beta)
-export(CAIC.sdmTMB)
 export(add_barrier_mesh)
 export(add_utm_columns)
+export(cAIC)
 export(censored_poisson)
 export(delta_beta)
 export(delta_gamma)

R/caic.R

@@ -1,110 +1,129 @@
-
-#' @title Calculate conditional AIC
+#' Calculate conditional AIC
 #'
-#' @description
 #' Calculates the conditional Akaike Information criterion (cAIC).
 #'
-#' @param object Output from \code{\link{sdmTMB}}
+#' @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.
+#'   (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.
 #'
-#' @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 crossvalidation.
-#' By contrast, \code{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.
+#' 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 as a byproduct the effective degrees of freedom,
-#' i.e., the number of fixed effects that would have an equivalent impact on
+#' 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 Thorson 2024.
+#' Both cAIC and EDF are calculated using Eq. 6 of Zheng, Cadigan, and Thorson
+#' (2024).
 #'
-#' Note that, for models that include profiled fixed effects, these profiles
-#' are turned off.
+#' 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
+#' 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**
+#' **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 behavior**
+#' **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.sdmTMB <-
-function( object,
-          what = c("CAIC","EDF") ){
+cAIC <- function(object, what = c("cAIC", "EDF"), ...) {
+  UseMethod("cAIC", object)
+}
 
-  what = match.arg(what)
-  require(Matrix)
-  tmb_data = object$tmb_data
+#' @exportS3Method
+cAIC.sdmTMB <- function(object, what = c("cAIC", "EDF")) {
+  what <- match.arg(what)
+  what <- tolower(what)
+  tmb_data <- object$tmb_data
 
-  # Make sure 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 )
+  ## 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 )
+  ## 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()
+  par <- obj$env$parList()
   parDataMode <- obj$env$last.par
-  indx = obj$env$lrandom()
-  q = length(indx)
-  p = length(object$model$par)
+  indx <- obj$env$lrandom()
+  q <- length(indx)
+  p <- length(object$model$par)
 
-  ## use - for Hess because model returns negative loglikelihood;
-  #cov_Psi_inv = -Hess_new[indx,indx]; ## this is the marginal prec mat of REs;
-  Hess_new = -Matrix(obj_new$env$f(parDataMode,order=1,type="ADGrad"),sparse = TRUE)
-  Hess_new = Hess_new[indx,indx]
+  ## 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(obj$env$f(parDataMode,order=1,type="ADGrad"),sparse = TRUE)
-  Hess = Hess[indx,indx]
-  negEDF = diag(solve(Hess, Hess_new))
+  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)
+  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)
-  }
-  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)
+  } 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")
   }
 }

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)
+})