From 3e2f82baebaab36772c3a0a4fbefdc4227e8a7d8 Mon Sep 17 00:00:00 2001 From: Sean Anderson Date: Tue, 24 Oct 2023 11:19:52 -0700 Subject: [PATCH 1/4] Allow delta models to work with ggpredict --- NAMESPACE | 1 + R/methods.R | 12 ++++++++++++ R/predict.R | 2 ++ R/utils.R | 32 ++++++++++++++++++++++++++++++++ man/set_delta_model.Rd | 35 +++++++++++++++++++++++++++++++++++ 5 files changed, 82 insertions(+) create mode 100644 man/set_delta_model.Rd diff --git a/NAMESPACE b/NAMESPACE index 690bfa480..b97077bdc 100644 --- a/NAMESPACE +++ b/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) diff --git a/R/methods.R b/R/methods.R index 49a57c84d..095875645 100644 --- a/R/methods.R +++ b/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,15 @@ 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)) { + return(x$formula[[attr(fit, "delta_model_predict")]]) + } if (length(x$formula) > 1L) { if ("visreg_model" %in% names(x)) { return(x$formula[[x$visreg_model]]) diff --git a/R/predict.R b/R/predict.R index f81352ad4..6fd66da5e 100644 --- a/R/predict.R +++ b/R/predict.R @@ -292,9 +292,11 @@ predict.sdmTMB <- function(object, newdata = NULL, sims <- nsim } + # TODO check if missing and only use attribute then!! assert_that(model[[1]] %in% c(NA, 1, 2), msg = "`model` argument not valid; should be one of NA, 1, 2") model <- model[[1]] + if (.has_delta_attr(object)) model <- attr(object, "delta_model_predict") # for ggpredict type <- match.arg(type) # FIXME parallel setup here? diff --git a/R/utils.R b/R/utils.R index 076a35cb6..77852f783 100644 --- a/R/utils.R +++ b/R/utils.R @@ -520,3 +520,35 @@ get_censored_upper <- function( upper_bound[prop_removed >= pstar] round(high) } + +#' Set delta model for prediction/\pkg{ggeffects} +#' +#' @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. +#' +#' @export +#' @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) |> +#' 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]") + +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 +} diff --git a/man/set_delta_model.Rd b/man/set_delta_model.Rd new file mode 100644 index 000000000..7cf0e11fc --- /dev/null +++ b/man/set_delta_model.Rd @@ -0,0 +1,35 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/utils.R +\name{set_delta_model} +\alias{set_delta_model} +\title{Set delta model for prediction/\pkg{ggeffects}} +\usage{ +set_delta_model(x, model = c(NA, 1, 2)) +} +\arguments{ +\item{x}{An \code{\link[=sdmTMB]{sdmTMB()}} model fit with a delta family such as +\code{\link[=delta_gamma]{delta_gamma()}}.} + +\item{model}{Which delta/hurdle model component to predict/plot with.} +} +\description{ +Set delta model for prediction/\pkg{ggeffects} +} +\examples{ +\dontshow{if (require("ggeffects", quietly = TRUE)) (if (getRversion() >= "3.4") withAutoprint else force)(\{ # examplesIf} +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]") +\dontshow{\}) # examplesIf} +} From b10ce3d884538440a1b5922a2f9a1d7f9e21cba9 Mon Sep 17 00:00:00 2001 From: Sean Anderson Date: Tue, 24 Oct 2023 12:03:19 -0700 Subject: [PATCH 2/4] Check formula for delta NA --- R/methods.R | 13 ++++++++++++- 1 file changed, 12 insertions(+), 1 deletion(-) diff --git a/R/methods.R b/R/methods.R index 095875645..5486c8338 100644 --- a/R/methods.R +++ b/R/methods.R @@ -196,7 +196,12 @@ df.residual.sdmTMB <- function(object, ...) { #' @export formula.sdmTMB <- function (x, ...) { if (.has_delta_attr(x)) { - return(x$formula[[attr(fit, "delta_model_predict")]]) + 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)) { @@ -245,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) From 858e60f3baf4af096172e59aeb134c9e2902f9ab Mon Sep 17 00:00:00 2001 From: Sean Anderson Date: Wed, 1 Nov 2023 09:16:12 -0700 Subject: [PATCH 3/4] Conditionally check attribute --- R/predict.R | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/R/predict.R b/R/predict.R index 6fd66da5e..d57984896 100644 --- a/R/predict.R +++ b/R/predict.R @@ -292,11 +292,12 @@ predict.sdmTMB <- function(object, newdata = NULL, sims <- nsim } - # TODO check if missing and only use attribute then!! 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]] - if (.has_delta_attr(object)) model <- attr(object, "delta_model_predict") # for ggpredict type <- match.arg(type) # FIXME parallel setup here? From 03b95aa329e4fadb4fa679f62d55869577e7eb9b Mon Sep 17 00:00:00 2001 From: Sean Anderson Date: Wed, 1 Nov 2023 10:25:38 -0700 Subject: [PATCH 4/4] Docs --- R/utils.R | 41 +++++++++++++++++++++++++++++++++++---- _pkgdown.yml | 1 + man/set_delta_model.Rd | 44 ++++++++++++++++++++++++++++++++++++------ 3 files changed, 76 insertions(+), 10 deletions(-) diff --git a/R/utils.R b/R/utils.R index 77852f783..3896d3bc0 100644 --- a/R/utils.R +++ b/R/utils.R @@ -521,7 +521,9 @@ get_censored_upper <- function( round(high) } -#' Set delta model for prediction/\pkg{ggeffects} +#' 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()]. @@ -529,8 +531,10 @@ get_censored_upper <- function( #' `NA` does the combined prediction, `1` does the binomial part, and `2` #' does the positive part. #' -#' @export -#' @examplesIf require("ggeffects", quietly = TRUE) +#' @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()) #' @@ -545,7 +549,36 @@ get_censored_upper <- function( #' # 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}. +#' . +#' +#' @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") diff --git a/_pkgdown.yml b/_pkgdown.yml index b92e27840..947bb1b9b 100644 --- a/_pkgdown.yml +++ b/_pkgdown.yml @@ -37,6 +37,7 @@ reference: - run_extra_optimization - residuals.sdmTMB - replicate_df + - set_delta_model - title: 'Families' desc: | diff --git a/man/set_delta_model.Rd b/man/set_delta_model.Rd index 7cf0e11fc..68b5c89a6 100644 --- a/man/set_delta_model.Rd +++ b/man/set_delta_model.Rd @@ -2,7 +2,7 @@ % Please edit documentation in R/utils.R \name{set_delta_model} \alias{set_delta_model} -\title{Set delta model for prediction/\pkg{ggeffects}} +\title{Set delta model for \code{\link[ggeffects:ggpredict]{ggeffects::ggpredict()}}} \usage{ set_delta_model(x, model = c(NA, 1, 2)) } @@ -10,14 +10,27 @@ set_delta_model(x, model = c(NA, 1, 2)) \item{x}{An \code{\link[=sdmTMB]{sdmTMB()}} model fit with a delta family such as \code{\link[=delta_gamma]{delta_gamma()}}.} -\item{model}{Which delta/hurdle model component to predict/plot with.} +\item{model}{Which delta/hurdle model component to predict/plot with. +\code{NA} does the combined prediction, \code{1} does the binomial part, and \code{2} +does the positive part.} +} +\value{ +The fitted model with a new attribute named \code{delta_model_predict}. +We suggest you use \code{set_delta_model()} in a pipe (as in the examples) +so that this attribute does not persist. Otherwise, \code{\link[=predict.sdmTMB]{predict.sdmTMB()}} +will choose this model component by default. You can also remove the +attribute yourself after: + +\if{html}{\out{
}}\preformatted{attr(fit, "delta_model_predict") <- NULL +}\if{html}{\out{
}} } \description{ -Set delta model for prediction/\pkg{ggeffects} +Set a delta model component to predict from with \code{\link[ggeffects:ggpredict]{ggeffects::ggpredict()}}. } -\examples{ -\dontshow{if (require("ggeffects", quietly = TRUE)) (if (getRversion() >= "3.4") withAutoprint else force)(\{ # examplesIf} -fit <- sdmTMB(density ~ poly(depth_scaled, 2), data = pcod_2011, +\details{ +A complete version of the examples below would be: + +\if{html}{\out{
}}\preformatted{fit <- sdmTMB(density ~ poly(depth_scaled, 2), data = pcod_2011, spatial = "off", family = delta_gamma()) # binomial part: @@ -31,5 +44,24 @@ set_delta_model(fit, model = 2) |> # combined: set_delta_model(fit, model = NA) |> ggeffects::ggpredict("depth_scaled [all]") +}\if{html}{\out{
}} + +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}. +\url{https://github.com/strengejacke/ggeffects}. +} +\examples{ +\dontshow{if (require("ggeffects", quietly = TRUE)) (if (getRversion() >= "3.4") withAutoprint else force)(\{ # examplesIf} +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) \dontshow{\}) # examplesIf} }