diff --git a/DESCRIPTION b/DESCRIPTION index 862d12f7..197b1ac7 100644 --- a/DESCRIPTION +++ b/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"), diff --git a/NEWS.md b/NEWS.md index 9574e32f..c321ec74 100644 --- a/NEWS.md +++ b/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 ; also see . J.T. Thorson wrote the function code. diff --git a/R/caic.R b/R/caic.R index 8279795d..64f9b684 100644 --- a/R/caic.R +++ b/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) diff --git a/R/print.R b/R/print.R index 54b33baf..90b907d7 100644 --- a/R/print.R +++ b/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)