From 6231927fce1706c229260d8e3cd3fbf16d8479f6 Mon Sep 17 00:00:00 2001 From: Sean Anderson Date: Thu, 19 Oct 2023 11:11:36 -0700 Subject: [PATCH] Copy edit docs/clean code --- R/fit.R | 130 ++++++++++++------------------------------ R/predict.R | 30 +++++----- man/predict.sdmTMB.Rd | 27 +++++---- man/sdmTMB.Rd | 71 ++++++++++++----------- 4 files changed, 103 insertions(+), 155 deletions(-) diff --git a/R/fit.R b/R/fit.R index 85f222855..4db725a76 100644 --- a/R/fit.R +++ b/R/fit.R @@ -3,11 +3,9 @@ NULL #' Fit a spatial or spatiotemporal GLMM with TMB #' -#' Fit a spatial or spatiotemporal Gaussian random field generalized linear -#' mixed effects model (GLMM) with the TMB (Template Model Builder) R package and -#' the SPDE (stochastic partial differential equation) approach. This can be -#' useful for (dynamic) species distribution models and relative abundance index -#' standardization among many other uses. +#' Fit a spatial or spatiotemporal generalized linear mixed effects model (GLMM) +#' with the TMB (Template Model Builder) R package and the SPDE (stochastic +#' partial differential equation) approximation to Gaussian random fields. #' #' @param formula Model formula. IID random intercepts are possible using #' \pkg{lme4} syntax, e.g., `+ (1 | g)` where `g` is a column of class @@ -82,8 +80,7 @@ NULL #' \doi{10.1111/ecog.05176} and the [spatial trends #' vignette](https://pbs-assess.github.io/sdmTMB/articles/spatial-trend-models.html). #' Note this predictor should usually be centered to have mean zero and have a -#' standard deviation of approximately 1 and should likely also be included as -#' a main effect. +#' standard deviation of approximately 1. #' **The spatial intercept is controlled by the `spatial` argument**; therefore, #' include or exclude the spatial intercept by setting `spatial = 'on'` or #' `'off'`. The only time when it matters whether `spatial_varying` excludes @@ -97,8 +94,8 @@ NULL #' a name of the variable in the data frame. See the Details section below. #' @param offset A numeric vector representing the model offset *or* a character #' value representing the column name of the offset. In delta/hurdle models, -#' this applies only to the positive component. Usually a log -#' transformed variable. +#' this applies only to the positive component. Usually a log transformed +#' variable. #' @param extra_time Optional extra time slices (e.g., years) to include for #' interpolation or forecasting with the predict function. See the Details #' section below. @@ -163,7 +160,7 @@ NULL #' An object (list) of class `sdmTMB`. Useful elements include: #' #' * `sd_report`: output from [TMB::sdreport()] -#' * `gradients`: log likelihood gradients with respect to each fixed effect +#' * `gradients`: marginal log likelihood gradients with respect to each fixed effect #' * `model`: output from [stats::nlminb()] #' * `data`: the fitted data #' * `mesh`: the object that was supplied to the `mesh` argument @@ -202,8 +199,8 @@ NULL #' with `+ s(x, y)` or `+ t2(x, y)`; smooths can be specific to various factor #' levels, `+ s(x, by = group)`; the basis function dimensions may be specified, #' e.g. `+ s(x, k = 4)`; and various types of splines may be constructed such as -#' cyclic splines to model seasonality, `+ s(month, bs = "cc", k = 12)` (perhaps -#' with the `knots` argument also be supplied). +#' cyclic splines to model seasonality (perhaps with the `knots` argument also +#' be supplied). #' #' **Threshold models** #' @@ -216,7 +213,8 @@ NULL #' `+ logistic(variable)`. This option models the relationship as a logistic #' function of the 50% and 95% values. This is similar to length- or size-based #' selectivity in fisheries, and is parameterized by the points at which f(x) = -#' 0.5 or 0.95. See the [threshold vignette](https://pbs-assess.github.io/sdmTMB/articles/threshold-models.html). +#' 0.5 or 0.95. See the +#' [threshold vignette](https://pbs-assess.github.io/sdmTMB/articles/threshold-models.html). #' #' Note that only a single threshold covariate can be included and the same covariate #' is included in both components for the delta families. @@ -236,16 +234,7 @@ NULL #' time slices with process error. #' #' `extra_time` can also be used to fill in missing time steps for the purposes -#' of a random walk or AR(1) process if their inclusion makes the gaps between -#' time steps even. -#' -#' **Index standardization** -#' -#' For index standardization, you may wish to include `0 + as.factor(year)` -#' (or whatever the time column is called) in the formula. See a basic -#' example of index standardization in the relevant -#' [package vignette](https://pbs-assess.github.io/sdmTMB/articles/index-standardization.html). -#' You will need to specify the `time` argument. See [get_index()]. +#' of a random walk or AR(1) process if the gaps between time steps are uneven. #' #' **Regularization and priors** #' @@ -279,11 +268,19 @@ NULL #' The main advantage of specifying such models using a delta family (compared #' to fitting two separate models) is (1) coding simplicity and (2) calculation #' of uncertainty on derived quantities such as an index of abundance with -#' [get_index()] using the generalized delta method within TMB. Also, parameters -#' can be shared across the models. +#' [get_index()] using the generalized delta method within TMB. Also, selected +#' parameters can be shared across the models. #' #' See the [delta-model vignette](https://pbs-assess.github.io/sdmTMB/articles/delta-models.html). #' +#' **Index standardization** +#' +#' For index standardization, you may wish to include `0 + as.factor(year)` +#' (or whatever the time column is called) in the formula. See a basic +#' example of index standardization in the relevant +#' [package vignette](https://pbs-assess.github.io/sdmTMB/articles/index-standardization.html). +#' You will need to specify the `time` argument. See [get_index()]. +#' #' @references #' #' **Main reference introducing the package to cite when using sdmTMB:** @@ -305,7 +302,7 @@ NULL #' English, P., E.J. Ward, C.N. Rooper, R.E. Forrest, L.A. Rogers, K.L. Hunter, #' A.M. Edwards, B.M. Connors, S.C. Anderson. 2021. Contrasting climate velocity #' impacts in warm and cool locations show that effects of marine warming are -#' worse in already warmer temperate waters. In press at Fish and Fisheries. +#' worse in already warmer temperate waters. Fish and Fisheries. 23(1) 239-255. #' \doi{10.1111/faf.12613}. #' #' *Discussion of and illustration of some decision points when fitting these @@ -319,9 +316,10 @@ NULL #' *Application and description of threshold/break-point models:* #' #' Essington, T.E. S.C. Anderson, L.A.K. Barnett, H.M. Berger, S.A. Siedlecki, -#' E.J. Ward. Advancing statistical models to reveal the effect of dissolved -#' oxygen on the spatial distribution of marine taxa using thresholds and a -#' physiologically based index. In press at Ecography. \doi{10.1111/ecog.06249}. +#' E.J. Ward. 2022. Advancing statistical models to reveal the effect of +#' dissolved oxygen on the spatial distribution of marine taxa using thresholds +#' and a physiologically based index. Ecography. 2022: e06249 +#' \doi{10.1111/ecog.06249}. #' #' *Application to fish body condition:* #' @@ -329,7 +327,7 @@ NULL #' spatiotemporal individual condition of a bottom-associated marine fish. #' bioRxiv 2022.04.19.488709. \doi{10.1101/2022.04.19.488709}. #' -#' *A number of sections of the original TMB model code were adapted from the +#' *Several sections of the original TMB model code were adapted from the #' VAST R package:* #' #' Thorson, J.T., 2019. Guidance for decisions using the Vector Autoregressive @@ -363,13 +361,14 @@ NULL #' #' # Build a mesh to implement the SPDE approach: #' mesh <- make_mesh(pcod_2011, c("X", "Y"), cutoff = 20) -#' # * this example uses a fairly coarse mesh so these examples run quickly -#' # * `cutoff` is the minimum distance between mesh vertices in units of the +#' +#' # - this example uses a fairly coarse mesh so these examples run quickly +#' # - 'cutoff' is the minimum distance between mesh vertices in units of the #' # x and y coordinates -#' # * `cutoff = 10` or `cutoff = 15` might make more sense in applied situations -#' # for this dataset -#' # * or build any mesh in 'fmesher' and pass it to the `mesh` argument in `make_mesh()` -#' # * not needed if you will be turning off all spatial/spatiotemporal random fields +#' # - 'cutoff = 10' might make more sense in applied situations for this dataset +#' # - or build any mesh in 'fmesher' and pass it to the 'mesh' argument in make_mesh()` +#' # - the mesh is not needed if you will be turning off all +#' # spatial/spatiotemporal random fields #' #' # Quick mesh plot: #' plot(mesh) @@ -607,14 +606,10 @@ sdmTMB <- function( spatiotemporal <- rep("iid", n_m) } - if (is.null(time)) { spatial_only <- rep(TRUE, n_m) } else { spatial_only <- ifelse(spatiotemporal == "off", TRUE, FALSE) - # if(all(spatiotemporal == "off")) { - # cli_abort("Time needs to be null if spatiotemporal fields are not included") - # } } if (is.list(spatial)) { @@ -634,8 +629,6 @@ sdmTMB <- function( } if (!include_spatial && all(spatiotemporal == "off") || !include_spatial && all(spatial_only)) { - # message("Both spatial and spatiotemporal fields are set to 'off'.") - # control$map_rf <- TRUE no_spatial <- TRUE if (missing(mesh)) { mesh <- sdmTMB::pcod_mesh_2011 # internal data; fake! @@ -787,20 +780,6 @@ sdmTMB <- function( # "As of version 0.3.1, sdmTMB turns off the constant spatial field `omega_s` when `spatial_varying` is specified so that the intercept or factor-level means are fully described by the spatially varying random fields `zeta_s`.") cli_inform(paste(msg, collapse = " ")) } - # if (!omit_spatial_intercept & !length(attr(z_i, "contrasts"))) { - # msg <- c("The spatial intercept is now in the first element of the spatially varying random fields `zeta_s` instead of the constant spatial random field `omega_s`. This change in the output format occurred in version 0.3.1.") - # cli_inform(msg) - # } -# .int <- sum(grep("(Intercept)", colnames(z_i)) > 0) -# if (.int && !omit_spatial_intercept) { -# # msg <- c("Detected an intercept in `spatial_varying`.", -# # "Make sure you have `spatial = 'off'` set since this also represents a spatial intercept.") -# # cli_warn(msg) -# # actually, just do it! -# omit_spatial_intercept <- TRUE -# include_spatial <- TRUE -# spatial <- "on" -# } .int <- grep("(Intercept)", colnames(z_i)) if (sum(.int) > 0) z_i <- z_i[,-.int,drop=FALSE] spatial_varying <- colnames(z_i) @@ -1042,7 +1021,6 @@ sdmTMB <- function( # TODO: make this cleaner X_ij_list <- list() - #X_ij_array <- array(data = NA, dim = c(nrow(X_ij[[1]]), ncol(X_ij[[1]]), n_m)) for (i in seq_len(n_m)) X_ij_list[[i]] <- X_ij[[i]] n_t <- length(unique(data[[time]])) @@ -1174,8 +1152,6 @@ sdmTMB <- function( tmb_params$b_j <- stats::coef(temp) } - #if (delta && !is.null(thresh$threshold_parameter)) cli_abort("Thresholds not implemented with delta models yet.") # TODO DELTA - # Map off parameters not needed tmb_map <- map_all_params(tmb_params) tmb_map$b_j <- NULL @@ -1203,14 +1179,6 @@ sdmTMB <- function( tmb_map$ln_phi <- as.factor(tmb_map$ln_phi) if (!is.null(thresh[[1]]$threshold_parameter)) tmb_map$b_threshold <- NULL - # optional models on spatiotemporal sd parameter - # if (est_epsilon_re == 0L) { - # tmb_map <- c(tmb_map, - # list( - # ln_epsilon_re_sigma = factor(rep(NA, n_m)), - # epsilon_re = factor(rep(NA, tmb_data$n_t)) - # )) - # } if (est_epsilon_re == 1L) { tmb_map <- unmap(tmb_map, c("ln_epsilon_re_sigma","epsilon_re")) } @@ -1222,7 +1190,7 @@ sdmTMB <- function( original_tmb_data <- tmb_data - # # much faster on first phase!? + # much faster on first phase!? tmb_data$no_spatial <- 1L # tmb_data$include_spatial <- 0L tmb_data$include_spatial <- rep(0L, length(spatial)) # for 1st phase @@ -1287,11 +1255,6 @@ sdmTMB <- function( if (reml) tmb_random <- c(tmb_random, "b_j") if (reml && delta) tmb_random <- c(tmb_random, "b_j2") -## if (est_epsilon_model >= 2) { -## # model 2 = re model, model 3 = loglinear-re -## tmb_random <- c(tmb_random, "epsilon_rw") -## } - if (sm$has_smooths) { if (reml) tmb_random <- c(tmb_random, "bs") tmb_random <- c(tmb_random, "b_smooth") # smooth random effects @@ -1380,8 +1343,6 @@ sdmTMB <- function( prof <- c("b_j") if (delta) prof <- c(prof, "b_j2") - - out_structure <- structure(list( data = data, spde = spde, @@ -1445,7 +1406,6 @@ sdmTMB <- function( out_structure$do_index <- FALSE } - tmb_obj <- TMB::MakeADFun( data = tmb_data, parameters = tmb_params, map = tmb_map, profile = if (control$profile) prof else NULL, @@ -1546,23 +1506,6 @@ set_limits <- function(tmb_obj, lower, upper, loc = NULL, silent = TRUE) { .upper["ar1_phi"] <- stats::qlogis((0.999 + 1) / 2) } - # if (!is.null(loc) && !"ln_kappa" %in% union(names(lower), names(upper))) { - # .dist <- stats::dist(loc) - # range_limits <- c(min(.dist) * 0.25, max(.dist) * 4) - # .upper[names(.upper) == "ln_kappa"] <- log(sqrt(8) / range_limits[1]) - # .lower[names(.upper) == "ln_kappa"] <- log(sqrt(8) / range_limits[2]) - # message("Setting limits on range to ", - # round(range_limits[1], 1), " and ", - # round(range_limits[2], 1), ":\n", - # "half the minimum and twice the maximum knot distance." - # ) - # if (.upper[names(.upper) == "ln_kappa"][1] <= tmb_obj$par[["ln_kappa"]][1]) { - # .upper[names(.upper) == "ln_kappa"] <- tmb_obj$par[["ln_kappa"]][1] + 0.1 - # } - # if (.lower[names(.lower) == "ln_kappa"][1] >= tmb_obj$par[["ln_kappa"]][1]) { - # .lower[names(.lower) == "ln_kappa"] <- tmb_obj$par[["ln_kappa"]][1] - 0.1 - # } - # } list(lower = .lower, upper = .upper) } @@ -1624,7 +1567,6 @@ check_irregalar_time <- function(data, time, spatiotemporal, time_varying) { find_missing_time <- function(x) { if (!is.factor(x)) { ti <- sort(unique(x)) - # mindiff <- min(diff(ti)) mindiff <- 1L allx <- seq(min(ti), max(ti), by = mindiff) setdiff(allx, ti) diff --git a/R/predict.R b/R/predict.R index 20a9693db..28ff3ac24 100644 --- a/R/predict.R +++ b/R/predict.R @@ -4,8 +4,8 @@ #' Predict from an sdmTMB model #' -#' Make predictions from an sdmTMB model; can predict on the original or new -#' data. +#' Make predictions from an \pkg{sdmTMB} model; can predict on the original or +#' new data. #' #' @param object A model fitted with [sdmTMB()]. #' @param newdata A data frame to make predictions on. This should be a data @@ -29,32 +29,36 @@ #' 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 *does* affects [get_index()]. -#' @param nsim Experimental: If `> 0`, simulate from the joint precision +#' @param nsim If `> 0`, simulate from the joint precision #' matrix with `nsim` draws. Returns a matrix of `nrow(data)` by `nsim` #' 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 -#' characterize uncertainty on predictions in space with sdmTMB. +#' Can be useful for deriving uncertainty on predictions +#' (e.g., `apply(x, 1, sd)`) or propagating uncertainty. This is currently +#' the fastest way to characterize uncertainty on predictions in space with +#' sdmTMB. #' @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"`, +#' Defaults to 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 Deprecated. See `mcmc_samples`. #' @param mcmc_samples See `extract_mcmc()` in the #' \href{https://github.com/pbs-assess/sdmTMBextra}{sdmTMBextra} package for -#' more details and the Bayesian vignette. 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. +#' more details and the +#' \href{https://pbs-assess.github.io/sdmTMB/articles/web_only/bayesian.html}{Bayesian vignette}. +#' 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 *and* `nsim > 0` or #' `mcmc_samples` is supplied: `NA` returns the combined prediction from both #' components on the link scale for the positive component; `1` or `2` return #' the first or second model component only on the link or response scale -#' depending on the argument `type`. +#' depending on the argument `type`. For regular prediction from delta models, +#' both sets of predictions are returned. #' @param offset A numeric vector of optional offset values. If left at default #' `NULL`, the offset is implicitly left at 0. #' @param return_tmb_report Logical: return the output from the TMB -#' report? For regular prediction this is all the reported variables +#' report? For regular prediction, this is all the reported variables #' at the MLE parameter values. For `nsim > 0` or when `mcmc_samples` #' 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. @@ -499,8 +503,6 @@ predict.sdmTMB <- function(object, newdata = NULL, tmb_data$calc_se <- as.integer(se_fit) tmb_data$pop_pred <- as.integer(pop_pred) tmb_data$exclude_RE <- exclude_RE - # tmb_data$calc_index_totals <- as.integer(!se_fit) - # tmb_data$calc_cog <- as.integer(!se_fit) tmb_data$proj_spatial_index <- newdata$sdm_spatial_id tmb_data$proj_Zs <- sm$Zs tmb_data$proj_Xs <- sm$Xs diff --git a/man/predict.sdmTMB.Rd b/man/predict.sdmTMB.Rd index 88d5da1b8..53ea4bfaa 100644 --- a/man/predict.sdmTMB.Rd +++ b/man/predict.sdmTMB.Rd @@ -51,15 +51,17 @@ predictions. \code{~0} or \code{NA} for population-level predictions. No other options (e.g., some but not all random intercepts) are implemented yet. Only affects predictions with \code{newdata}. This \emph{does} affects \code{\link[=get_index]{get_index()}}.} -\item{nsim}{Experimental: If \verb{> 0}, simulate from the joint precision +\item{nsim}{If \verb{> 0}, simulate from the joint precision matrix with \code{nsim} draws. Returns a matrix of \code{nrow(data)} by \code{nsim} representing the estimates of the linear predictor (i.e., in link space). -Can be useful for deriving uncertainty on predictions (e.g., \code{apply(x, 1, sd)}) or propagating uncertainty. This is currently the fastest way to -characterize uncertainty on predictions in space with sdmTMB.} +Can be useful for deriving uncertainty on predictions +(e.g., \code{apply(x, 1, sd)}) or propagating uncertainty. This is currently +the fastest way to characterize uncertainty on predictions in space with +sdmTMB.} \item{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: \code{"omega_s"}, +Defaults to link-space predictions. Options include: \code{"omega_s"}, \code{"zeta_s"}, \code{"epsilon_st"}, and \code{"est_rf"} (as described below). Other options will be passed verbatim.} @@ -67,23 +69,26 @@ Other options will be passed verbatim.} \code{mcmc_samples} is supplied: \code{NA} returns the combined prediction from both components on the link scale for the positive component; \code{1} or \code{2} return the first or second model component only on the link or response scale -depending on the argument \code{type}.} +depending on the argument \code{type}. For regular prediction from delta models, +both sets of predictions are returned.} \item{offset}{A numeric vector of optional offset values. If left at default \code{NULL}, the offset is implicitly left at 0.} \item{mcmc_samples}{See \code{extract_mcmc()} in the \href{https://github.com/pbs-assess/sdmTMBextra}{sdmTMBextra} package for -more details and the Bayesian vignette. If specified, the predict function -will return a matrix of a similar form as if \code{nsim > 0} but representing -Bayesian posterior samples from the Stan model.} +more details and the +\href{https://pbs-assess.github.io/sdmTMB/articles/web_only/bayesian.html}{Bayesian vignette}. +If specified, the predict function will return a matrix of a similar form +as if \code{nsim > 0} but representing Bayesian posterior samples from the Stan +model.} \item{return_tmb_object}{Logical. If \code{TRUE}, will include the TMB object in a list format output. Necessary for the \code{\link[=get_index]{get_index()}} or \code{\link[=get_cog]{get_cog()}} functions.} \item{return_tmb_report}{Logical: return the output from the TMB -report? For regular prediction this is all the reported variables +report? For regular prediction, this is all the reported variables at the MLE parameter values. For \code{nsim > 0} or when \code{mcmc_samples} 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.} @@ -135,8 +140,8 @@ A matrix: } } \description{ -Make predictions from an sdmTMB model; can predict on the original or new -data. +Make predictions from an \pkg{sdmTMB} model; can predict on the original or +new data. } \examples{ \dontshow{if (ggplot2_installed()) (if (getRversion() >= "3.4") withAutoprint else force)(\{ # examplesIf} diff --git a/man/sdmTMB.Rd b/man/sdmTMB.Rd index 1176c5c1e..2bf9be1ba 100644 --- a/man/sdmTMB.Rd +++ b/man/sdmTMB.Rd @@ -113,8 +113,7 @@ random field is assumed to have a mean of 0. If a (scaled) time column is used, it will represent a local-time-trend model. See \doi{10.1111/ecog.05176} and the \href{https://pbs-assess.github.io/sdmTMB/articles/spatial-trend-models.html}{spatial trends vignette}. Note this predictor should usually be centered to have mean zero and have a -standard deviation of approximately 1 and should likely also be included as -a main effect. +standard deviation of approximately 1. \strong{The spatial intercept is controlled by the \code{spatial} argument}; therefore, include or exclude the spatial intercept by setting \code{spatial = 'on'} or \code{'off'}. The only time when it matters whether \code{spatial_varying} excludes @@ -130,8 +129,8 @@ a name of the variable in the data frame. See the Details section below.} \item{offset}{A numeric vector representing the model offset \emph{or} a character value representing the column name of the offset. In delta/hurdle models, -this applies only to the positive component. Usually a log -transformed variable.} +this applies only to the positive component. Usually a log transformed +variable.} \item{extra_time}{Optional extra time slices (e.g., years) to include for interpolation or forecasting with the predict function. See the Details @@ -189,7 +188,7 @@ be dragons.} An object (list) of class \code{sdmTMB}. Useful elements include: \itemize{ \item \code{sd_report}: output from \code{\link[TMB:sdreport]{TMB::sdreport()}} -\item \code{gradients}: log likelihood gradients with respect to each fixed effect +\item \code{gradients}: marginal log likelihood gradients with respect to each fixed effect \item \code{model}: output from \code{\link[stats:nlminb]{stats::nlminb()}} \item \code{data}: the fitted data \item \code{mesh}: the object that was supplied to the \code{mesh} argument @@ -201,11 +200,9 @@ An object (list) of class \code{sdmTMB}. Useful elements include: } } \description{ -Fit a spatial or spatiotemporal Gaussian random field generalized linear -mixed effects model (GLMM) with the TMB (Template Model Builder) R package and -the SPDE (stochastic partial differential equation) approach. This can be -useful for (dynamic) species distribution models and relative abundance index -standardization among many other uses. +Fit a spatial or spatiotemporal generalized linear mixed effects model (GLMM) +with the TMB (Template Model Builder) R package and the SPDE (stochastic +partial differential equation) approximation to Gaussian random fields. } \details{ \strong{Model description} @@ -235,8 +232,8 @@ Within these smooths, the same syntax commonly used in \code{\link[mgcv:s]{mgcv: with \code{+ s(x, y)} or \code{+ t2(x, y)}; smooths can be specific to various factor levels, \code{+ s(x, by = group)}; the basis function dimensions may be specified, e.g. \code{+ s(x, k = 4)}; and various types of splines may be constructed such as -cyclic splines to model seasonality, \code{+ s(month, bs = "cc", k = 12)} (perhaps -with the \code{knots} argument also be supplied). +cyclic splines to model seasonality (perhaps with the \code{knots} argument also +be supplied). \strong{Threshold models} @@ -249,7 +246,8 @@ Similarly, a logistic-function threshold model can be included via \code{+ logistic(variable)}. This option models the relationship as a logistic function of the 50\% and 95\% values. This is similar to length- or size-based selectivity in fisheries, and is parameterized by the points at which f(x) = -0.5 or 0.95. See the \href{https://pbs-assess.github.io/sdmTMB/articles/threshold-models.html}{threshold vignette}. +0.5 or 0.95. See the +\href{https://pbs-assess.github.io/sdmTMB/articles/threshold-models.html}{threshold vignette}. Note that only a single threshold covariate can be included and the same covariate is included in both components for the delta families. @@ -269,16 +267,7 @@ or a smoother on the time column provide mechanisms to predict over missing time slices with process error. \code{extra_time} can also be used to fill in missing time steps for the purposes -of a random walk or AR(1) process if their inclusion makes the gaps between -time steps even. - -\strong{Index standardization} - -For index standardization, you may wish to include \code{0 + as.factor(year)} -(or whatever the time column is called) in the formula. See a basic -example of index standardization in the relevant -\href{https://pbs-assess.github.io/sdmTMB/articles/index-standardization.html}{package vignette}. -You will need to specify the \code{time} argument. See \code{\link[=get_index]{get_index()}}. +of a random walk or AR(1) process if the gaps between time steps are uneven. \strong{Regularization and priors} @@ -312,10 +301,18 @@ through a single formula that is shared across the two models. The main advantage of specifying such models using a delta family (compared to fitting two separate models) is (1) coding simplicity and (2) calculation of uncertainty on derived quantities such as an index of abundance with -\code{\link[=get_index]{get_index()}} using the generalized delta method within TMB. Also, parameters -can be shared across the models. +\code{\link[=get_index]{get_index()}} using the generalized delta method within TMB. Also, selected +parameters can be shared across the models. See the \href{https://pbs-assess.github.io/sdmTMB/articles/delta-models.html}{delta-model vignette}. + +\strong{Index standardization} + +For index standardization, you may wish to include \code{0 + as.factor(year)} +(or whatever the time column is called) in the formula. See a basic +example of index standardization in the relevant +\href{https://pbs-assess.github.io/sdmTMB/articles/index-standardization.html}{package vignette}. +You will need to specify the \code{time} argument. See \code{\link[=get_index]{get_index()}}. } \examples{ \dontshow{if (require("visreg", quietly = TRUE)) (if (getRversion() >= "3.4") withAutoprint else force)(\{ # examplesIf} @@ -323,13 +320,14 @@ library(sdmTMB) # Build a mesh to implement the SPDE approach: mesh <- make_mesh(pcod_2011, c("X", "Y"), cutoff = 20) -# * this example uses a fairly coarse mesh so these examples run quickly -# * `cutoff` is the minimum distance between mesh vertices in units of the + +# - this example uses a fairly coarse mesh so these examples run quickly +# - 'cutoff' is the minimum distance between mesh vertices in units of the # x and y coordinates -# * `cutoff = 10` or `cutoff = 15` might make more sense in applied situations -# for this dataset -# * or build any mesh in 'fmesher' and pass it to the `mesh` argument in `make_mesh()` -# * not needed if you will be turning off all spatial/spatiotemporal random fields +# - 'cutoff = 10' might make more sense in applied situations for this dataset +# - or build any mesh in 'fmesher' and pass it to the 'mesh' argument in make_mesh()` +# - the mesh is not needed if you will be turning off all +# spatial/spatiotemporal random fields # Quick mesh plot: plot(mesh) @@ -531,7 +529,7 @@ velocities:} English, P., E.J. Ward, C.N. Rooper, R.E. Forrest, L.A. Rogers, K.L. Hunter, A.M. Edwards, B.M. Connors, S.C. Anderson. 2021. Contrasting climate velocity impacts in warm and cool locations show that effects of marine warming are -worse in already warmer temperate waters. In press at Fish and Fisheries. +worse in already warmer temperate waters. Fish and Fisheries. 23(1) 239-255. \doi{10.1111/faf.12613}. \emph{Discussion of and illustration of some decision points when fitting these @@ -545,9 +543,10 @@ e12783. \doi{10.7717/peerj.12783}. \emph{Application and description of threshold/break-point models:} Essington, T.E. S.C. Anderson, L.A.K. Barnett, H.M. Berger, S.A. Siedlecki, -E.J. Ward. Advancing statistical models to reveal the effect of dissolved -oxygen on the spatial distribution of marine taxa using thresholds and a -physiologically based index. In press at Ecography. \doi{10.1111/ecog.06249}. +E.J. Ward. 2022. Advancing statistical models to reveal the effect of +dissolved oxygen on the spatial distribution of marine taxa using thresholds +and a physiologically based index. Ecography. 2022: e06249 +\doi{10.1111/ecog.06249}. \emph{Application to fish body condition:} @@ -555,7 +554,7 @@ Lindmark, M., S.C. Anderson, M. Gogina, M. Casini. Evaluating drivers of spatiotemporal individual condition of a bottom-associated marine fish. bioRxiv 2022.04.19.488709. \doi{10.1101/2022.04.19.488709}. -\emph{A number of sections of the original TMB model code were adapted from the +\emph{Several sections of the original TMB model code were adapted from the VAST R package:} Thorson, J.T., 2019. Guidance for decisions using the Vector Autoregressive