Add EDF printing to smoothers #383 #387

Author: seananderson

DESCRIPTION

@@ -1,7 +1,7 @@
 Type: Package
 Package: sdmTMB
 Title: Spatial and Spatiotemporal SPDE-Based GLMMs with 'TMB'
-Version: 0.6.0.9017
+Version: 0.6.0.9018
 Authors@R: c(
     person(c("Sean", "C."), "Anderson", , "sean@seananderson.ca",
            role = c("aut", "cre"),

NEWS.md

@@ -1,5 +1,9 @@
 # sdmTMB (development version)
 
+* Add EDF (effective degrees of freedom) printing to smoothers with 
+  `print.sdmTMB()` and `summary.sdmTMB()` if argument `edf = TRUE` 
+  is included. E.g. `print.sdmTMB(fit, edf = TRUE)`. #383 #387
+
 * Add `cAIC()` for calculating *conditional* AIC. Theory based on
   <https://arxiv.org/abs/2411.14185>; also see
   <https://doi.org/10.1002/ecy.4327>. J.T. Thorson wrote the function code.

R/caic.R

@@ -68,7 +68,9 @@ cAIC <- function(object, what = c("cAIC", "EDF"), ...) {
 
 #' @exportS3Method
 cAIC.sdmTMB <- function(object, what = c("cAIC", "EDF"), ...) {
-  what <- match.arg(what)
+
+  what <- tolower(what)
+  what <- match.arg(what, choices = c("caic", "edf"))
   what <- tolower(what)
   tmb_data <- object$tmb_data
 
@@ -104,6 +106,10 @@ cAIC.sdmTMB <- function(object, what = c("cAIC", "EDF"), ...) {
   p <- length(object$model$par)
 
   ## use '-' for Hess because model returns negative loglikelihood
+  if (is.null(object$tmb_random)) {
+    cli_inform(c("This model has no random effects.", "cAIC and EDF only apply to models with random effects."))
+    return(invisible(NULL))
+  }
   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
 
@@ -119,7 +125,21 @@ cAIC.sdmTMB <- function(object, what = c("cAIC", "EDF"), ...) {
     return(cAIC_out)
   } else if (what == "edf") {
     ## Figure out group for each random-effect coefficient
-    group <- factor(names(object$last.par.best[obj$env$random]))
+    group <- names(object$last.par.best[obj$env$random])
+
+    convert_bsmooth2names <- function(object, model = 1) {
+      sn <- row.names(print_smooth_effects(object, m = model, silent = TRUE)$smooth_sds)
+      sn <- gsub("^sd", "", sn)
+      dms <- object$smoothers$sm_dims
+      unlist(lapply(seq_along(dms), \(i) rep(sn[i], dms[i])))
+
+    }
+    s_groups <- convert_bsmooth2names(object)
+    # smoothers always shared in delta models
+    if (is_delta(object)) s_groups <- c(paste0("1LP-", s_groups), paste0("2LP-", s_groups))
+    group[group == "b_smooth"] <- s_groups
+    group <- factor(group)
+
     ## Calculate total EDF by group
     EDF <- tapply(negEDF, INDEX = group, FUN = length) - tapply(negEDF, INDEX = group, FUN = sum)
     return(EDF)

R/print.R

@@ -88,7 +88,7 @@ print_main_effects <- function(x, m = 1) {
   mm
 }
 
-print_smooth_effects <- function(x, m = 1) {
+print_smooth_effects <- function(x, m = 1, edf = NULL, silent = FALSE) {
   sr <- x$sd_report
   sr_se <- as.list(sr, "Std. Error")
   sr_est <- as.list(sr, "Estimate")
@@ -122,7 +122,7 @@ print_smooth_effects <- function(x, m = 1) {
           "This does not affect model fitting.",
           "We'll use generic covariate names ('scovariate') here intead."
         )
-        cli_warn(msg)
+        if (!silent) cli_warn(msg)
         row.names(mm_sm) <- paste0("scovariate-", seq_len(nrow(mm_sm)))
       } else {
         mm_sm <- NULL
@@ -136,6 +136,16 @@ print_smooth_effects <- function(x, m = 1) {
     re_sm_mat[, 1] <- smooth_sds
     rownames(re_sm_mat) <- sm_names_sds
     colnames(re_sm_mat) <- "Std. Dev."
+
+    if (!is.null(edf)) {
+      if (is_delta(x)) {
+        lp_regex <- paste0("^", m, "LP-")
+        edf <- edf[grepl(lp_regex, names(edf))]
+      }
+      edf <- round(edf, 2)
+      re_sm_mat <- cbind(re_sm_mat, matrix(edf, ncol = 1))
+      colnames(re_sm_mat)[2] <- "EDF"
+    }
   } else {
     re_sm_mat <- NULL
     mm_sm <- NULL
@@ -312,10 +322,15 @@ print_header <- function(x) {
   cat(info$overall_family)
 }
 
-print_one_model <- function(x, m = 1) {
+print_one_model <- function(x, m = 1, edf = FALSE, silent = FALSE) {
+  if (edf) {
+    .edf <- suppressMessages(cAIC(x, what = "EDF"))
+  } else {
+    .edf <- NULL
+  }
   info <- print_model_info(x)
   main <- print_main_effects(x, m = m)
-  smooth <- print_smooth_effects(x, m = m)
+  smooth <- print_smooth_effects(x, m = m, edf = .edf, silent = silent)
   iid_re <- print_iid_re(x, m = m)
   tv <- print_time_varying(x, m = m)
   range <- print_range(x, m = m)
@@ -381,10 +396,10 @@ print.sdmTMB <- function(x, ...) {
   delta <- isTRUE(x$family$delta)
   print_header(x)
   if (delta) cat("\nDelta/hurdle model 1: -----------------------------------\n")
-  print_one_model(x, 1)
+  print_one_model(x, 1, ...)
   if (delta) {
     cat("\nDelta/hurdle model 2: -----------------------------------\n")
-    print_one_model(x, 2)
+    print_one_model(x, 2, ...)
   }
   if (delta) cat("\n")
   print_footer(x)