Merge branch 'ggeffects-delta'

NAMESPACE

@@ -67,6 +67,7 @@ export(sdmTMB_simulate)
 export(sdmTMB_stacking)
 export(sdmTMBcontrol)
 export(sdmTMBpriors)
+export(set_delta_model)
 export(spread_sims)
 export(student)
 export(tidy)

R/methods.R

@@ -145,6 +145,11 @@ extractAIC.sdmTMB <- function(fit, scale, k = 2, ...) {
 #' @importFrom stats family
 #' @export
 family.sdmTMB <- function (object, ...) {
+  if (.has_delta_attr(object)) {
+    which_model <- attr(object, "delta_model_predict")
+    if (is.na(which_model)) which_model <- 2L # combined; for link
+    return(object$family[[which_model]])
+  }
   if ("visreg_model" %in% names(object)) {
     return(object$family[[object$visreg_model]])
   } else {
@@ -184,8 +189,20 @@ df.residual.sdmTMB <- function(object, ...) {
   nobs(object) - length(object$model$par)
 }
 
+.has_delta_attr <- function(x) {
+  "delta_model_predict" %in% names(attributes(x))
+}
+
 #' @export
 formula.sdmTMB <- function (x, ...) {
+  if (.has_delta_attr(x)) {
+    which_model <- attr(x, "delta_model_predict")
+    if (!identical(x$formula[[1]], x$formula[[2]]) && is.na(which_model)) {
+      cli_abort("Delta component formulas are not the same but ggeffects::ggpredict() is trying to predict on the combined model. For now, predict on one or the other component, or keep the formulas the same, or write your own prediction and plot code.")
+    }
+    if (is.na(which_model)) which_model <- 1L # combined take 1!?
+    return(x$formula[[which_model]])
+  }
   if (length(x$formula) > 1L) {
     if ("visreg_model" %in% names(x)) {
       return(x$formula[[x$visreg_model]])
@@ -233,6 +250,12 @@ Effect.sdmTMB <- function(focal.predictors, mod, ...) {
     cli_abort("Please install the effects package")
   }
 
+  if (is_delta(mod)) {
+    msg <- paste0("Effect() and ggeffects::ggeffect() do not yet work with ",
+      "sdmTMB delta/hurdle models. Please use ggeffects::ggpredict() instead.")
+    cli_abort(msg)
+  }
+
   vc <- vcov(mod)
   b <- tidy(mod, silent = TRUE)
 

R/predict.R

@@ -294,6 +294,9 @@ predict.sdmTMB <- function(object, newdata = NULL,
 
   assert_that(model[[1]] %in% c(NA, 1, 2),
     msg = "`model` argument not valid; should be one of NA, 1, 2")
+  if (missing(model)) {
+    if (.has_delta_attr(object)) model <- attr(object, "delta_model_predict") # for ggpredict
+  }
   model <- model[[1]]
   type <- match.arg(type)
   # FIXME parallel setup here?

R/utils.R

@@ -520,3 +520,68 @@ get_censored_upper <- function(
     upper_bound[prop_removed >= pstar]
   round(high)
 }
+
+#' Set delta model for [ggeffects::ggpredict()]
+#'
+#' Set a delta model component to predict from with [ggeffects::ggpredict()].
+#'
+#' @param x An [sdmTMB::sdmTMB()] model fit with a delta family such as
+#'   [sdmTMB::delta_gamma()].
+#' @param model Which delta/hurdle model component to predict/plot with.
+#'   `NA` does the combined prediction, `1` does the binomial part, and `2`
+#'   does the positive part.
+#'
+#' @details
+#' A complete version of the examples below would be:
+#'
+#' ```
+#' fit <- sdmTMB(density ~ poly(depth_scaled, 2), data = pcod_2011,
+#'   spatial = "off", family = delta_gamma())
+#'
+#' # binomial part:
+#' set_delta_model(fit, model = 1) |>
+#'   ggeffects::ggpredict("depth_scaled [all]")
+#'
+#' # gamma part:
+#' set_delta_model(fit, model = 2) |>
+#'   ggeffects::ggpredict("depth_scaled [all]")
+#'
+#' # combined:
+#' set_delta_model(fit, model = NA) |>
+#'   ggeffects::ggpredict("depth_scaled [all]")
+#' ```
+#'
+#' But cannot be run on CRAN until a version of \pkg{ggeffects} > 1.3.2
+#' is on CRAN. For now, you can install the GitHub version of \pkg{ggeffects}.
+#' <https://github.com/strengejacke/ggeffects>.
+#'
+#' @returns
+#' The fitted model with a new attribute named `delta_model_predict`.
+#' We suggest you use `set_delta_model()` in a pipe (as in the examples)
+#' so that this attribute does not persist. Otherwise, [predict.sdmTMB()]
+#' will choose this model component by default. You can also remove the
+#' attribute yourself after:
+#'
+#' ```
+#' attr(fit, "delta_model_predict") <- NULL
+#' ```
+#'
+#' @examplesIf require("ggeffects", quietly = TRUE)
+#' fit <- sdmTMB(density ~ poly(depth_scaled, 2), data = pcod_2011,
+#'   spatial = "off", family = delta_gamma())
+#'
+#' # binomial part:
+#' set_delta_model(fit, model = 1)
+#'
+#' # gamma part:
+#' set_delta_model(fit, model = 2)
+#'
+#' # combined:
+#' set_delta_model(fit, model = NA)
+#' @export
+set_delta_model <- function(x, model = c(NA, 1, 2)) {
+  assertthat::assert_that(model[[1]] %in% c(NA, 1, 2),
+    msg = "`model` argument not valid; should be one of NA, 1, 2")
+  attr(x, "delta_model_predict") <- model[[1]]
+  x
+}

_pkgdown.yml

@@ -37,6 +37,7 @@ reference:
   - run_extra_optimization
   - residuals.sdmTMB
   - replicate_df
+  - set_delta_model
 
 - title: 'Families'
   desc: |