Compare 76310c4 ... +249 ... 39a1e97

#' Extract MCMC samples from a model fit with [tmbstan::tmbstan()].

#'

#' @description

#' `r lifecycle::badge("experimental")`

#'

#' Returns a matrix of parameter samples. Rows correspond to the order

#' of `your_model$tmb_obj$env$last.par.best`. Columns correspond to

#' posterior samples. Is used internally by [predict.sdmTMB()] to make

#' fully Bayesian predictions. See the `tmbstan_model` argument

#' in [predict.sdmTMB()].

#'

#' @param object Output from [tmbstan::tmbstan()] run on the `tmb_obj`

#'   element of an [sdmTMB()] model. E.g., `tmbstan(your_model$tmb_obj)`.

#' @examples

#'

#' \dontrun{

#' pcod_spde <- make_mesh(pcod, c("X", "Y"), cutoff = 30)

#' plot(pcod_spde)

#'

#' # here we will fix the random field parameters at their approximate

#' # MLEs (maximum likelihood estimates) from a previous fit

#' # to improve speed of convergence:

#' m_tmb <- sdmTMB(density ~ 0 + as.factor(year),

#'   data = pcod, mesh = pcod_spde, family = tweedie(link = "log"), time = "year",

#'   control = sdmTMBcontrol(start = list(ln_kappa = rep(-1.58, 2),

#'     ln_tau_E = -0.15, ln_tau_O = -0.65),

#'     map = list(ln_kappa = rep(factor(NA), 2),

#'       ln_tau_E = factor(NA), ln_tau_O = factor(NA))))

#' m_tmb

#'

#' # will take 3-5 minutes:

#' library(tmbstan)

#' m_stan <- tmbstan(m_tmb$tmb_obj, iter = 200, chains = 1)

#' print(m_stan, pars = c("b_j", "thetaf", "ln_phi", "omega_s[1]", "epsilon_st[1]"))

#'

#' post <- extract_mcmc(m_stan)

#' dim(post)

#'

#' p <- predict(m_tmb, newdata = qcs_grid, tmbstan_model = m_stan)

#' p_last <- p[qcs_grid$year == max(qcs_grid$year), ] # just plot last year

#' pred <- qcs_grid[qcs_grid$year == max(qcs_grid$year), ]

#' pred$est <- apply(exp(p_last), 1, median)

#' pred$lwr <- apply(exp(p_last), 1, quantile, probs = 0.1)

#' pred$upr <- apply(exp(p_last), 1, quantile, probs = 0.9)

#' pred$cv <- apply(exp(p_last), 1, function(x) sd(x) / mean(x))

#'

#' library(ggplot2)

#' ggplot(pred, aes(X, Y, fill = est)) + geom_raster() +

#'   scale_fill_viridis_c(trans = "log")

#' ggplot(pred, aes(X, Y, fill = cv)) + geom_raster() +

#'   scale_fill_viridis_c(trans = "log")

#'

#' index_quantiles <- get_index_sims(p)

#' ggplot(index_quantiles, aes(year, est, ymin = lwr, ymax = upr)) +

#'   geom_line() + geom_ribbon(alpha = 0.5)

#'

#' index_samples <- get_index_sims(p, return_sims = TRUE)

#' ggplot(index_samples, aes(as.factor(year), .value)) +

#'   geom_violin()

#' }

#' @export

extract_mcmc <- function(object) {

  }

      } else {

      }

    }

  })

}

#' Plot sdmTMB models with the \pkg{visreg} package

#'

#' sdmTMB models fit with regular (non-delta) families can be passed to

#' [visreg::visreg()] or [visreg::visreg2d()] directly. Examples are shown

#' below. Delta models can use the helper functions `visreg_delta()` or

#' `visreg2d_delta()` described here.

#'

#' @param object Fit from [sdmTMB()]

#' @param model 1st or 2nd delta model

#' @param ... Any arguments passed to [visreg::visreg()] or

#'   [visreg::visreg2d()]

#'

#' @details

#' Note the residuals are currently randomized quantile residuals,

#' *not* deviance residuals as is usual for GLMs with \pkg{visreg}.

#'

#' @return

#' A plot from the visreg package. Optionally, the data plotted invisibly if

#' `plot = FALSE`. This is useful if you want to make your own plot after.

#'

#' @export

#' @rdname visreg_delta

#'

#' @examples

#' if (inla_installed() &&

#'   require("ggplot2", quietly = TRUE) &&

#'   require("visreg", quietly = TRUE)) {

#'

#'   pcod_2011$fyear <- as.factor(pcod_2011$year)

#'   fit <- sdmTMB(

#'     density ~ s(depth_scaled) + fyear,

#'     data = pcod_2011, mesh = pcod_mesh_2011,

#'     spatial = "off",

#'     family = tweedie()

#'   )

#'   visreg::visreg(fit, xvar = "depth_scaled")

#'   visreg::visreg(fit, xvar = "fyear")

#'   visreg::visreg(fit, xvar = "depth_scaled", scale = "response")

#'   visreg::visreg2d(fit, xvar = "fyear", yvar = "depth_scaled")

#'   v <- visreg::visreg(fit, xvar = "depth_scaled")

#'   head(v$fit)

#'   # now use ggplot2 etc. if desired

#'

#'   # Delta model example:

#'   fit_dg <- sdmTMB(

#'     density ~ s(depth_scaled, year, k = 8),

#'     data = pcod_2011, mesh = pcod_mesh_2011,

#'     spatial = "off",

#'     family = delta_gamma()

#'   )

#'   visreg_delta(fit_dg, xvar = "depth_scaled", model = 1, gg = TRUE)

#'   visreg_delta(fit_dg, xvar = "depth_scaled", model = 2, gg = TRUE)

#'   visreg_delta(fit_dg,

#'     xvar = "depth_scaled", model = 1,

#'     scale = "response", gg = TRUE

#'   )

#'   visreg_delta(fit_dg,

#'     xvar = "depth_scaled", model = 2,

#'     scale = "response"

#'   )

#'   visreg_delta(fit_dg,

#'     xvar = "depth_scaled", model = 2,

#'     scale = "response", gg = TRUE, rug = FALSE

#'   )

#'   visreg2d_delta(fit_dg,

#'     xvar = "depth_scaled", yvar = "year",

#'     model = 2, scale = "response"

#'   )

#'   visreg2d_delta(fit_dg,

#'     xvar = "depth_scaled", yvar = "year",

#'     model = 1, scale = "response", plot.type = "persp"

#'   )

#'   visreg2d_delta(fit_dg,

#'     xvar = "depth_scaled", yvar = "year",

#'     model = 2, scale = "response", plot.type = "gg"

#'   )

#' }

visreg_delta <- function(object, ..., model = c(1, 2)) {

}

#' @export

#' @rdname visreg_delta

visreg2d_delta <- function(object, ..., model = c(1, 2)) {

}

check_model_arg <- function(model) {

}

#' @param model Which model if a delta model.

#'

#' @return A plot of eigenvectors illustrating the estimated anisotropy. A list

#'   of the plotted data is invisibly returned.

#' @references Code adapted from VAST R package

#' if (inla_installed()) {

#'   d <- pcod

#'   m <- sdmTMB(

#'     data = d,

#'     formula = density ~ 0 + as.factor(year),

#'     time = "year", mesh = make_mesh(d, c("X", "Y"), n_knots = 80, type = "kmeans"),

#'     family = tweedie(link = "log"), anisotropy = TRUE,

#'     spatial = "off"

#'   )

#'   plot_anisotropy(m)

#' }

plot_anisotropy <- function(object, model = 1) {

  )

#' Plot a smooth term from an sdmTMB model

#'

#' @param object An [sdmTMB()] model.

#' @param select The smoother term to plot.

#' @param n The number of equally spaced points to evaluate the smoother along.

#' @param level The confidence level.

#' @param ggplot Logical: use the \pkg{ggplot2} package?

#' @param rug Logical: add rug lines along the lower axis?

#' @param return_data Logical: return the predicted data instead of making a plot?

#' @export

#'

#' @details

#' Note:

#' * Any numeric predictor is set to its mean

#' * Any factor predictor is set to its first-level value

#' * The time element (if present) is set to its minimum value

#' * The x and y coordinates are set to their mean values

#'

#' @examples

#' if (inla_installed()) {

#'   d <- subset(pcod, year >= 2000 & density > 0)

#'   pcod_spde <- make_mesh(d, c("X", "Y"), cutoff = 30)

#'   m <- sdmTMB(

#'     data = d,

#'     formula = log(density) ~ s(depth_scaled) + s(year, k = 5),

#'     mesh = pcod_spde

#'   )

#'   plot_smooth(m)

#' }

plot_smooth <- function(object, select = 1, n = 100, level = 0.95,

                        ggplot = FALSE, rug = TRUE, return_data = FALSE) {

  )

    }

  }

  }

      }

    }

  }

      } else {

      }

    }

  }

  }

    } else {

    }

    )

    }

  } else {

    )) +

        )

    }

  }

}

#' Make predictions from an sdmTMB model; can predict on the original or new

#' data.

#' @param newdata A data frame to make predictions on. This should be a data

#'   frame with the same predictor columns as in the fitted data and a time

#'   column (if this is a spatiotemporal model) with the same name as in the

#'   fitted data. There should be predictor data for each year in the original

#'   data set.

#' @param type Should the `est` column be in link (default) or response space?

#' @param se_fit Should standard errors on predictions at the new locations

#'   given by `newdata` be calculated? Warning: the current implementation can

#'   be very slow for large data sets or high-resolution projections. A *much*

#'   faster option is to use the `nsim` argument below and calculate uncertainty

#'   on the simulations from the joint precision matrix.

#' @param return_tmb_object Logical. If `TRUE`, will include the TMB object in a

#'   list format output. Necessary for the [get_index()] or [get_cog()]

#'   functions.

#' @param re_form `NULL` to specify including all spatial/spatiotemporal random

#'   effects in predictions. `~0` or `NA` for population-level predictions. Note

#'   that unlike lme4 or glmmTMB, this only affects what the standard errors are

#'   calculated on if `se_fit = TRUE`. This does not affect [get_index()]

#'   calculations.

#' @param re_form_iid `NULL` to specify including all random intercepts in the

#'   predictions. `~0` or `NA` for population-level predictions. No other

#'   options (e.g., some but not all random intercepts) are implemented yet.

#'   Only affects predictions with `newdata`. This also affects [get_index()].

#' @param nsim **Experimental.** If > 0, simulate from the joint precision matrix with `sims`

#'   draws Returns a matrix of `nrow(data)` by `sim` representing the estimates

#'   of the linear predictor (i.e., in link space). Can be useful for deriving

#'   uncertainty on predictions (e.g., `apply(x, 1, sd)`) or propagating

#'   uncertainty. This is currently the fastest way to generate estimates of

#'   uncertainty on predictions in space with sdmTMB.

#' @param area **Deprecated**. Please use `area` in [get_index()].

#' @param sims **Deprecated**. Please use `nsim` instead.

#' @param sims_var **Experimental.** Which TMB reported variable from the model

#'   should be extracted from the joint precision matrix simulation draws?

#'   Defaults to the link-space predictions. Options include: `"omega_s"`,

#'   `"zeta_s"`, `"epsilon_st"`, and `"est_rf"` (as described below).

#'   Other options will be passed verbatim.

#' @param tmbstan_model A model fit with [tmbstan::tmbstan()]. See

#'   [extract_mcmc()] for more details and an example. If specified, the

#'   predict function will return a matrix of a similar form as if `nsim > 0`

#'   but representing Bayesian posterior samples from the Stan model.

#' @param model Type of prediction if a delta/hurdle model:

#'   `NA` returns the combined prediction from both components on

#'   the response scale; `1` or `2` return the first or second model

#'   component only on the link or response scale depending on the argument

#'   `type`.

#' @param return_tmb_report Logical: return the output from the TMB

#'   report? For regular prediction this is all the reported variables

#'   at the MLE parameter values. For `nsim > 0` or when `tmbstan_model`

#'   is supplied, this is a list where each element is a sample and the

#'   contents of each element is the output of the report for that sample.

#' If `return_tmb_object = FALSE` (and `nsim = 0` and `tmbstan_model = NULL`):

#'

#' * `epsilon_st`: Spatiotemporal (intercept) random fields, could be

#'    off (zero), IID, AR1, or random walk

#'

#' If `return_tmb_object = TRUE` (and `nsim = 0` and `tmbstan_model = NULL`):

#' * `obj`: The TMB object returned from the prediction run

#' * `fit_obj`: The original TMB model object

#'

#' In this case, you likely only need the `data` element as an end user.

#' The other elements are included for other functions.

#'

#' If `nsim > 0` or `tmbstan_model` is not `NULL`:

#' A matrix:

#'

#' * Columns represent samples

#' * Rows represent predictions with one row per row of `newdata`

#' if (require("ggplot2", quietly = TRUE) && inla_installed()) {

#' d <- pcod_2011

#' mesh <- make_mesh(d, c("X", "Y"), cutoff = 30) # a coarse mesh for example speed

#'  time = "year", mesh = mesh, family = tweedie(link = "log")

#'

#' qcs_grid_2011 <- subset(qcs_grid, year >= min(pcod_2011$year))

#' predictions <- predict(m, newdata = qcs_grid_2011)

#' # Because this is a spatiotemporal model, you'll need at least one time

#' # element. If time isn't also a fixed effect then it doesn't matter what you pick:

#' nd$year <- 2011L # L: integer to match original data

#'  data = d, formula = density ~ 0 + as.factor(year) + s(depth_scaled, k = 5),

#'  time = "year", mesh = mesh, family = tweedie(link = "log")

#' plot_smooth(m_gam)

#'

#' # or manually:

#' nd$year <- 2011L

#' mesh <- make_mesh(d, c("X", "Y"), cutoff = 15)

#' unique(d$year)

#'   data = d, formula = density ~ 1,

#'   spatiotemporal = "AR1", # using an AR1 to have something to forecast with

#'   extra_time = 2019L, # `L` for integer to match our data

#'   spatial = "off",

#'   time = "year", mesh = mesh, family = tweedie(link = "log")

#' grid2019 <- qcs_grid_2011[qcs_grid_2011$year == max(qcs_grid_2011$year), ]

#' grid2019$year <- 2019L # `L` because `year` is an integer in the data

#' qcsgrid_forecast <- rbind(qcs_grid_2011, grid2019)

#' d <- pcod

#' d$year_scaled <- as.numeric(scale(d$year))

#' mesh <- make_mesh(pcod, c("X", "Y"), cutoff = 25)

#' m <- sdmTMB(data = d, formula = density ~ depth_scaled + depth_scaled2,

#'   mesh = mesh, family = tweedie(link = "log"),

#'   spatial_varying = ~ 0 + year_scaled, time = "year", spatiotemporal = "off")

#' nd <- qcs_grid

#' nd$year_scaled <- (nd$year - mean(d$year)) / sd(d$year)

#' p <- predict(m, newdata = nd)

predict.sdmTMB <- function(object, newdata = object$data,

  type = c("link", "response"),

  se_fit = FALSE,

  return_tmb_object = FALSE,

  area = deprecated(), re_form = NULL, re_form_iid = NULL, nsim = 0,

  sims = deprecated(),

  tmbstan_model = NULL,

  sims_var = "est",

  model = c(NA, 1, 2),

  return_tmb_report = FALSE,

  ...) {

  } else {

  }

  } else {

  }

  } else {

  }

  } else {

  }

  }

  }

  } else {

  }

      )

    # newdata$sdm_orig_id <- seq(1, nrow(newdata))

    # fake_newdata <- unique(newdata[,xy_cols])

    # fake_newdata[["sdm_spatial_id"]] <- seq(1, nrow(fake_newdata)) - 1L

    #

    # newdata <- base::merge(newdata, fake_newdata, by = xy_cols,

    #   all.x = TRUE, all.y = FALSE)

    # newdata <- newdata[order(newdata$sdm_orig_id),, drop = FALSE]

      # this formula has breakpt() etc. in it:

    } else {

    }

    # deal with prediction IID random intercepts:

    ## not checking so that not all factors need to be in prediction:

    # fct_check <- vapply(RE_names, function(x) check_valid_factor_levels(data[[x]], .name = x), TRUE)

    # TODO DELTA hardcoded to 1:

    }

    # tmb_data$calc_index_totals <- as.integer(!se_fit)

    # tmb_data$calc_cog <- as.integer(!se_fit)

    # SVC:

    } else {

    }

      # covariate vector dimensioned by number of time steps

      }

    }

      } else {

      }

      } else {

      }

    }

        } else {

        }

      }

    }

        }

        }

        } else {

        }

      } else { # not a delta model: