sdmTMB.Rd

% Generated by roxygen2: do not edit by hand
% Please edit documentation in R/fit.R
\name{sdmTMB}
\alias{sdmTMB}
\title{Fit a spatial or spatiotemporal GLMM with TMB}
\usage{
sdmTMB(
  formula,
  data,
  mesh,
  time = NULL,
  family = gaussian(link = "identity"),
  spatial = c("on", "off"),
  spatiotemporal = c("iid", "ar1", "rw", "off"),
  share_range = TRUE,
  time_varying = NULL,
  time_varying_type = c("rw", "rw0", "ar1"),
  spatial_varying = NULL,
  weights = NULL,
  offset = NULL,
  extra_time = NULL,
  reml = FALSE,
  silent = TRUE,
  anisotropy = FALSE,
  control = sdmTMBcontrol(),
  priors = sdmTMBpriors(),
  knots = NULL,
  bayesian = FALSE,
  previous_fit = NULL,
  do_fit = TRUE,
  do_index = FALSE,
  predict_args = NULL,
  index_args = NULL,
  experimental = NULL
)
}
\arguments{
\item{formula}{Model formula. IID random intercepts and slopes are possible using
\pkg{lme4} syntax, e.g., \code{+ (1 | g)} or \code{+ (0 + depth | g)} or \code{+ (1 + depth | g)} where \code{g} is a column of class
character or factor representing groups. Penalized splines are possible via
\pkg{mgcv} with \code{s()}. Optionally a list for delta (hurdle) models.  See
examples and details below.}

\item{data}{A data frame.}

\item{mesh}{An object from \code{\link[=make_mesh]{make_mesh()}}.}

\item{time}{An optional time column name (as character). Can be left as
\code{NULL} for a model with only spatial random fields; however, if the data
are actually spatiotemporal and you wish to use \code{\link[=get_index]{get_index()}} or \code{\link[=get_cog]{get_cog()}}
downstream, supply the time argument.}

\item{family}{The family and link. Supports \code{\link[=gaussian]{gaussian()}}, \code{\link[=Gamma]{Gamma()}},
\code{\link[=binomial]{binomial()}}, \code{\link[=poisson]{poisson()}}, \code{\link[sdmTMB:families]{Beta()}},
\code{\link[sdmTMB:families]{nbinom2()}},
\code{\link[sdmTMB:families]{truncated_nbinom2()}},
\code{\link[sdmTMB:families]{nbinom1()}},
\code{\link[sdmTMB:families]{truncated_nbinom1()}},
\code{\link[sdmTMB:families]{censored_poisson()}},
\code{\link[sdmTMB:families]{gamma_mix()}},
\code{\link[sdmTMB:families]{lognormal_mix()}},
\code{\link[sdmTMB:families]{student()}},
\code{\link[sdmTMB:families]{tweedie()}}, and
\code{\link[sdmTMB:families]{gengamma()}}.
Supports the delta/hurdle models:
\code{\link[sdmTMB:families]{delta_beta()}},
\code{\link[sdmTMB:families]{delta_gamma()}},
\code{\link[sdmTMB:families]{delta_gamma_mix()}},
\code{\link[sdmTMB:families]{delta_lognormal_mix()}},
\code{\link[sdmTMB:families]{delta_lognormal()}}, and
\code{\link[sdmTMB:families]{delta_truncated_nbinom2()}},
For binomial family options, see 'Binomial families' in the Details
section below.}

\item{spatial}{Estimate spatial random fields? Options are \code{'on'} / \code{'off'}
or \code{TRUE} / \code{FALSE}. Optionally, a list for delta models, e.g. \code{list('on', 'off')}.}

\item{spatiotemporal}{Estimate the spatiotemporal random fields as \code{'iid'}
(independent and identically distributed; default), stationary \code{'ar1'}
(first-order autoregressive), a random walk (\code{'rw'}), or fixed at 0
\code{'off'}. Will be set to \code{'off'} if \code{time = NULL}. If a delta model, can be
a list. E.g., \code{list('off', 'ar1')}. Note that the spatiotemporal standard
deviation represents the marginal steady-state standard deviation of the
process in the case of the AR1. I.e., it is scaled according to the
correlation. See the \href{https://kaskr.github.io/adcomp/classdensity_1_1AR1__t.html}{TMB documentation}.
If the AR1 correlation coefficient (rho) is estimated close to 1,
say > 0.99, then you may wish to switch to the random walk \code{'rw'}.
Capitalization is ignored. \code{TRUE} gets converted to \code{'iid'} and \code{FALSE}
gets converted to \code{'off'}.}

\item{share_range}{Logical: estimate a shared spatial and spatiotemporal
range parameter (\code{TRUE}, default) or independent range parameters
(\code{FALSE}). If a delta model, can be a list. E.g., \code{list(TRUE, FALSE)}.}

\item{time_varying}{An optional one-sided formula describing covariates
that should be modelled as a time-varying process. Set the type of
process with \code{time_varying_type}. See the help for \code{time_varying_type}
for warnings about modelling the first time step. Structure shared in
delta models.}

\item{time_varying_type}{Type of time-varying process to apply to
\code{time_varying} formula. \code{'rw'} indicates a random walk with the first
time step estimated independently (included for legacy reasons), \code{'rw0'}
indicates a random walk with the first time step estimated with
a mean-zero normal prior, \code{'ar1'} indicates a \href{https://kaskr.github.io/adcomp/classdensity_1_1AR1__t.html}{stationary first-order autoregressive process}
with the first time step estimated with a mean-zero prior. In the case of
\code{'rw'}, be careful not to include covariates (including the intercept) in
both the main and time-varying formula since the first time step is
estimated independently. I.e., in this case, at least one should have \code{~ 0} or \code{~ -1}. Structure shared in delta models.}

\item{spatial_varying}{An optional one-sided formula of coefficients that
should vary in space as random fields. Note that you likely want to include
a fixed effect for the same variable to improve interpretability since the
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.
\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
an intercept is in the case of factor predictors. In this case, if
\code{spatial_varying} excludes the intercept (\code{~ 0} or \code{~ -1}), you should set
\code{spatial = 'off'} to match.  Structure must be shared in delta models.}

\item{weights}{A numeric vector representing optional likelihood weights for
the conditional model. Implemented as in \pkg{glmmTMB}: weights do not have
to sum to one and are not internally modified. Can also be used for trials
with the binomial family; the \code{weights} argument needs to be a vector and not
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.}

\item{extra_time}{Optional extra time slices (e.g., years) to include for
interpolation or forecasting with the predict function. See the Details
section below.}

\item{reml}{Logical: use REML (restricted maximum likelihood) estimation
rather than maximum likelihood? Internally, this adds the fixed effects to
the list of random effects to integrate over.}

\item{silent}{Silent or include optimization details? Helpful to set to
\code{FALSE} for models that take a while to fit.}

\item{anisotropy}{Logical: allow for anisotropy (spatial correlation that is
directionally dependent)? See \code{\link[=plot_anisotropy]{plot_anisotropy()}}.
Must be shared across delta models.}

\item{control}{Optimization control options via \code{\link[=sdmTMBcontrol]{sdmTMBcontrol()}}.}

\item{priors}{Optional penalties/priors via \code{\link[=sdmTMBpriors]{sdmTMBpriors()}}. Must currently
be shared across delta models.}

\item{knots}{Optional named list containing knot values to be used for basis
construction of smoothing terms. See \code{\link[mgcv:gam]{mgcv::gam()}} and \code{\link[mgcv:gamm]{mgcv::gamm()}}.
E.g., \verb{s(x, bs = 'cc', k = 4), knots = list(x = c(1, 2, 3, 4))}}

\item{bayesian}{Logical indicating if the model will be passed to
\pkg{tmbstan}. If \code{TRUE}, Jacobian adjustments are applied to account for
parameter transformations when priors are applied.}

\item{previous_fit}{A previously fitted sdmTMB model to initialize the
optimization with. Can greatly speed up fitting. Note that the model must
be set up \emph{exactly} the same way. However, the data and \code{weights} arguments
can change, which can be useful for cross-validation.}

\item{do_fit}{Fit the model (\code{TRUE}) or return the processed data without
fitting (\code{FALSE})?}

\item{do_index}{Do index standardization calculations while fitting? Saves
memory and time when working with large datasets or projection grids since
the TMB object doesn't have to be rebuilt with \code{\link[=predict.sdmTMB]{predict.sdmTMB()}} and
\code{\link[=get_index]{get_index()}}. If \code{TRUE}, then \code{predict_args} must have a \code{newdata} element
supplied and \code{area} can be supplied to \code{index_args}.
Most users can ignore this option. The fitted object can be passed directly
to \code{\link[=get_index]{get_index()}}.}

\item{predict_args}{A list of arguments to pass to \code{\link[=predict.sdmTMB]{predict.sdmTMB()}} \strong{if}
\code{do_index = TRUE}. Most users can ignore this option.}

\item{index_args}{A list of arguments to pass to \code{\link[=get_index]{get_index()}} \strong{if}
\code{do_index = TRUE}. Currently, only \code{area} is supported. Bias correction
can be done when calling \code{\link[=get_index]{get_index()}} on the resulting fitted object.
Most users can ignore this option.}

\item{experimental}{A named list for esoteric or in-development options. Here
be dragons.}
}
\value{
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}: 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
\item \code{family}: the family object, which includes the inverse link function as \code{family$linkinv()}
\item \code{tmb_params}: The parameters list passed to \code{\link[TMB:MakeADFun]{TMB::MakeADFun()}}
\item \code{tmb_map}: The 'map' list passed to \code{\link[TMB:MakeADFun]{TMB::MakeADFun()}}
\item \code{tmb_data}: The data list passed to \code{\link[TMB:MakeADFun]{TMB::MakeADFun()}}
\item \code{tmb_obj}: The TMB object created by \code{\link[TMB:MakeADFun]{TMB::MakeADFun()}}
}
}
\description{
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}

See the \href{https://pbs-assess.github.io/sdmTMB/articles/model-description.html}{model description}
vignette or the relevant appendix of the preprint on sdmTMB:
\doi{10.1101/2022.03.24.485545}

\strong{Binomial families}

Following the structure of \code{\link[stats:glm]{stats::glm()}} and \pkg{glmmTMB}, a binomial
family can be specified in one of 4 ways: (1) the response may be a factor
(and the model classifies the first level versus all others), (2) the
response may be binomial (0/1), (3) the response can be a matrix of form
\code{cbind(success, failure)}, and (4) the response may be the observed
proportions, and the 'weights' argument is used to specify the Binomial size
(N) parameter (\verb{prob ~ ..., weights = N}).

\strong{Smooth terms}

Smooth terms can be included following GAMs (generalized additive models)
using \code{+ s(x)}, which implements a smooth from \code{\link[mgcv:s]{mgcv::s()}}. \pkg{sdmTMB} uses
penalized smooths, constructed via \code{\link[mgcv:smooth2random]{mgcv::smooth2random()}}. This is a similar
approach implemented in \pkg{gamm4} and \pkg{brms}, among other packages.
Within these smooths, the same syntax commonly used in \code{\link[mgcv:s]{mgcv::s()}} or
\code{\link[mgcv:t2]{mgcv::t2()}} can be applied, e.g. 2-dimensional smooths may be constructed
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 (perhaps with the \code{knots} argument also
be supplied).

\strong{Threshold models}

A linear break-point relationship for a covariate can be included via
\code{+ breakpt(variable)} in the formula, where \code{variable} is a single covariate
corresponding to a column in \code{data}. In this case, the relationship is linear
up to a point and then constant (hockey-stick shaped).

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}.

Note that only a single threshold covariate can be included and the same covariate
is included in both components for the delta families.

\strong{Extra time: forecasting or interpolating}

Extra time slices (e.g., years) can be included for interpolation or
forecasting with the predict function via the \code{extra_time} argument. The
predict function requires all time slices to be defined when fitting the
model to ensure the various time indices are set up correctly. Be careful if
including extra time slices that the model remains identifiable. For example,
including \code{+ as.factor(year)} in \code{formula} will render a model with no data
to inform the expected value in a missing year. \code{\link[=sdmTMB]{sdmTMB()}} makes no attempt
to determine if the model makes sense for forecasting or interpolation. The
options \code{time_varying}, \code{spatiotemporal = "rw"}, \code{spatiotemporal = "ar1"},
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 the gaps between time steps are uneven.

\code{extra_time} can include only extra time steps or all time steps including
those found in the fitted data. This latter option may be simpler.

\strong{Regularization and priors}

You can achieve regularization via penalties (priors) on the fixed effect
parameters. See \code{\link[=sdmTMBpriors]{sdmTMBpriors()}}. You can fit the model once without
penalties and look at the output of \code{print(your_model)} or \code{tidy(your_model)}
or fit the model with \code{do_fit = FALSE} and inspect
\code{head(your_model$tmb_data$X_ij[[1]])} if you want to see how the formula is
translated to the fixed effect model matrix. Also see the
\href{https://pbs-assess.github.io/sdmTMB/articles/bayesian.html}{Bayesian vignette}.

\strong{Delta/hurdle models}

Delta models (also known as hurdle models) can be fit as two separate models
or at the same time by using an appropriate delta family. E.g.:
\code{\link[sdmTMB:families]{delta_gamma()}},
\code{\link[sdmTMB:families]{delta_beta()}},
\code{\link[sdmTMB:families]{delta_lognormal()}}, and
\code{\link[sdmTMB:families]{delta_truncated_nbinom2()}}.
If fit with a delta family, by default the formula, spatial, and spatiotemporal
components are shared. Some elements can be specified independently for the two models
using a list format. These include \code{formula}, \code{spatial}, \code{spatiotemporal},
and \code{share_range}. The first element of the list is for the binomial component
and the second element is for the positive component (e.g., Gamma).
Other elements must be shared for now (e.g., spatially varying coefficients,
time-varying coefficients). Furthermore, there are currently limitations if
specifying two formulas as a list: the two formulas cannot have smoothers or
threshold effects. For now, these must be specified
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, 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) && require("ggeffects", quietly = TRUE)) (if (getRversion() >= "3.4") withAutoprint else force)(\{ # examplesIf}
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
#   x and y coordinates
# - '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)

# Fit a Tweedie spatial random field GLMM with a smoother for depth:
fit <- sdmTMB(
  density ~ s(depth),
  data = pcod_2011, mesh = mesh,
  family = tweedie(link = "log")
)
fit

# Extract coefficients:
tidy(fit, conf.int = TRUE)
tidy(fit, effects = "ran_par", conf.int = TRUE)

# Perform several 'sanity' checks:
sanity(fit)

# Predict on the fitted data; see ?predict.sdmTMB
p <- predict(fit)

# Predict on new data:
p <- predict(fit, newdata = qcs_grid)
head(p)

\donttest{
# Visualize the depth effect with ggeffects:
ggeffects::ggpredict(fit, "depth [all]") |> plot()

# Visualize depth effect with visreg: (see ?visreg_delta)
visreg::visreg(fit, xvar = "depth") # link space; randomized quantile residuals
visreg::visreg(fit, xvar = "depth", scale = "response")
visreg::visreg(fit, xvar = "depth", scale = "response", gg = TRUE, rug = FALSE)

# Add spatiotemporal random fields:
fit <- sdmTMB(
  density ~ 0 + as.factor(year),
  time = "year", #<
  data = pcod_2011, mesh = mesh,
  family = tweedie(link = "log")
)
fit

# Make the fields AR1:
fit <- sdmTMB(
  density ~ s(depth),
  time = "year",
  spatial = "off",
  spatiotemporal = "ar1", #<
  data = pcod_2011, mesh = mesh,
  family = tweedie(link = "log")
)
fit

# Make the fields a random walk:
fit <- sdmTMB(
  density ~ s(depth),
  time = "year",
  spatial = "off",
  spatiotemporal = "rw", #<
  data = pcod_2011, mesh = mesh,
  family = tweedie(link = "log")
)
fit

# Depth smoothers by year:
fit <- sdmTMB(
  density ~ s(depth, by = as.factor(year)), #<
  time = "year",
  spatial = "off",
  spatiotemporal = "rw",
  data = pcod_2011, mesh = mesh,
  family = tweedie(link = "log")
)
fit

# 2D depth-year smoother:
fit <- sdmTMB(
  density ~ s(depth, year), #<
  spatial = "off",
  data = pcod_2011, mesh = mesh,
  family = tweedie(link = "log")
)
fit

# Turn off spatial random fields:
fit <- sdmTMB(
  present ~ poly(log(depth)),
  spatial = "off", #<
  data = pcod_2011, mesh = mesh,
  family = binomial()
)
fit

# Which, matches glm():
fit_glm <- glm(
  present ~ poly(log(depth)),
  data = pcod_2011,
  family = binomial()
)
summary(fit_glm)
AIC(fit, fit_glm)

# Delta/hurdle binomial-Gamma model:
fit_dg <- sdmTMB(
  density ~ poly(log(depth), 2),
  data = pcod_2011, mesh = mesh,
  spatial = "off",
  family = delta_gamma() #<
)
fit_dg

# Delta model with different formulas and spatial structure:
fit_dg <- sdmTMB(
  list(density ~ depth_scaled, density ~ poly(depth_scaled, 2)), #<
  data = pcod_2011, mesh = mesh,
  spatial = list("off", "on"), #<
  family = delta_gamma()
)
fit_dg

# Delta/hurdle truncated NB2:
pcod_2011$count <- round(pcod_2011$density)
fit_nb2 <- sdmTMB(
  count ~ s(depth),
  data = pcod_2011, mesh = mesh,
  spatial = "off",
  family = delta_truncated_nbinom2() #<
)
fit_nb2

# Regular NB2:
fit_nb2 <- sdmTMB(
  count ~ s(depth),
  data = pcod_2011, mesh = mesh,
  spatial = "off",
  family = nbinom2() #<
)
fit_nb2

# IID random intercepts by year:
pcod_2011$fyear <- as.factor(pcod_2011$year)
fit <- sdmTMB(
  density ~ s(depth) + (1 | fyear), #<
  data = pcod_2011, mesh = mesh,
  family = tweedie(link = "log")
)
fit

# IID random slopes and intercepts (implicit) by year:
fit <- sdmTMB(
  density ~ (depth | fyear), #<
  data = pcod_2011, mesh = mesh,
  family = tweedie(link = "log")
)

# Spatially varying coefficient of year:
pcod_2011$year_scaled <- as.numeric(scale(pcod_2011$year))
fit <- sdmTMB(
  density ~ year_scaled,
  spatial_varying = ~ 0 + year_scaled, #<
  data = pcod_2011, mesh = mesh, family = tweedie(), time = "year"
)
fit

# Time-varying effects of depth and depth squared:
fit <- sdmTMB(
  density ~ 0 + as.factor(year),
  time_varying = ~ 0 + depth_scaled + depth_scaled2, #<
  data = pcod_2011, time = "year", mesh = mesh,
  family = tweedie()
)
print(fit)
# Extract values:
est <- as.list(fit$sd_report, "Estimate")
se <- as.list(fit$sd_report, "Std. Error")
est$b_rw_t[, , 1]
se$b_rw_t[, , 1]

# Linear break-point effect of depth:
fit <- sdmTMB(
  density ~ breakpt(depth_scaled), #<
  data = pcod_2011,
  mesh = mesh,
  family = tweedie()
)
fit
}
\dontshow{\}) # examplesIf}
}
\references{
\strong{Main reference introducing the package to cite when using sdmTMB:}

Anderson, S.C., E.J. Ward, P.A. English, L.A.K. Barnett. 2022. sdmTMB: an R
package for fast, flexible, and user-friendly generalized linear mixed effects
models with spatial and spatiotemporal random fields.
bioRxiv 2022.03.24.485545; \doi{10.1101/2022.03.24.485545}.

\emph{Reference for local trends:}

Barnett, L.A.K., E.J. Ward, S.C. Anderson. 2021. Improving estimates of species
distribution change by incorporating local trends. Ecography. 44(3):427-439.
\doi{10.1111/ecog.05176}.

\emph{Further explanation of the model and application to calculating climate
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. Fish and Fisheries. 23(1) 239-255.
\doi{10.1111/faf.12613}.

\emph{Discussion of and illustration of some decision points when fitting these
models:}

Commander, C.J.C., L.A.K. Barnett, E.J. Ward, S.C. Anderson, T.E.
Essington. 2022. The shadow model: how and why small choices in
spatially explicit species distribution models affect predictions. PeerJ 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. 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:}

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{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
Spatio-Temporal (VAST) package in stock, ecosystem, habitat and climate
assessments. Fish. Res. 210:143–161.
\doi{10.1016/j.fishres.2018.10.013}.

\emph{Code for the \code{family} R-to-TMB implementation, selected parameterizations of
the observation likelihoods, general package structure inspiration, and the
idea behind the TMB prediction approach were adapted from the glmmTMB R
package:}

Brooks, M.E., K. Kristensen, K.J. van Benthem, A. Magnusson,
C.W. Berg, A. Nielsen, H.J. Skaug, M. Maechler, B.M. Bolker.
2017. glmmTMB Balances Speed and Flexibility Among Packages for
Zero-inflated Generalized Linear Mixed Modeling. The R Journal, 9(2):378-400.
\doi{10.32614/rj-2017-066}.

\emph{Implementation of geometric anisotropy with the SPDE and use of random
field GLMMs for index standardization}:

Thorson, J.T., A.O. Shelton, E.J. Ward, H.J. Skaug. 2015.
Geostatistical delta-generalized linear mixed models improve precision for
estimated abundance indices for West Coast groundfishes. ICES J. Mar. Sci.
72(5): 1297–1310. \doi{10.1093/icesjms/fsu243}.
}