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..5486c8338 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,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)
 
diff --git a/R/predict.R b/R/predict.R
index 03c64d46e..c77d9613b 100644
--- a/R/predict.R
+++ b/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?
diff --git a/R/utils.R b/R/utils.R
index 076a35cb6..3896d3bc0 100644
--- a/R/utils.R
+++ b/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
+}
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
new file mode 100644
index 000000000..68b5c89a6
--- /dev/null
+++ b/man/set_delta_model.Rd
@@ -0,0 +1,67 @@
+% 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 \code{\link[ggeffects:ggpredict]{ggeffects::ggpredict()}}}
+\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.
+\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{<div class="sourceCode">}}\preformatted{attr(fit, "delta_model_predict") <- NULL
+}\if{html}{\out{</div>}}
+}
+\description{
+Set a delta model component to predict from with \code{\link[ggeffects:ggpredict]{ggeffects::ggpredict()}}.
+}
+\details{
+A complete version of the examples below would be:
+
+\if{html}{\out{<div class="sourceCode">}}\preformatted{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]")
+}\if{html}{\out{</div>}}
+
+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}
+}