sdmTMB model structure
An sdmTMB model combines relationships with measured predictors,
variation across space and time, and an observation distribution
describing variation around the expected response. For example, a model
of species density might use depth as a predictor, a spatial field for
persistent differences among locations, and a spatiotemporal field for
spatial patterns that change from year to year. With a log link, a
simple version is
\[
\log(\mu_{\mathbf{s},t}) = \beta_0 +
\beta_{\mathrm{depth}}\,\mathrm{depth}_{\mathbf{s},t} +
\omega_{\mathbf{s}} + \epsilon_{\mathbf{s},t}.
\]
Here, \(\mu_{\mathbf{s},t}\) is the
expected density conditional on the random fields. The spatial field
\(\omega_{\mathbf{s}}\) describes
spatial variation remaining after accounting for depth that persists
through time, while \(\epsilon_{\mathbf{s},t}\) describes
additional spatial deviations that change through time. All terms on the
right are on the log scale: a field value of \(\log(2)\) doubles the expected density,
holding the other terms constant. The observation distribution describes
variation around this conditional mean. The fields describe spatial and
spatiotemporal patterns, not their causes.
The general model structure is
\[
\begin{aligned}
y_{\mathbf{s},t} &\sim \mathcal{D} \left( \mu_{\mathbf{s},t}, \phi
\right),\\
\mu_{\mathbf{s},t} &= f^{-1} \left( \eta_{\mathbf{s},t} \right),\\
\eta_{\mathbf{s},t} &=
\underbrace{\mathbf{X}^{\mathrm{main}}_{\mathbf{s},t}
\boldsymbol{\beta}}_{\substack{\text{main}\\\text{effects}}} +
\underbrace{O_{\mathbf{s},t}}_{\text{offset}} +
\underbrace{\sum_{k} \mathbf{z}_{k,\mathbf{s},t}^\top
\mathbf{b}_{k,g}}_{\substack{\text{random intercepts}\\\text{and
slopes}}} +
\underbrace{\mathbf{X}^{\mathrm{tvc}}_{\mathbf{s},t}
\boldsymbol{\gamma}_{t}}_{\substack{\text{time-varying}\\\text{coefficients}}}
+
\underbrace{\mathbf{X}^{\mathrm{svc}}_{\mathbf{s},t}
\boldsymbol{\zeta}_{\mathbf{s}}}_{\substack{\text{spatially
varying}\\\text{coefficients}}} +
\underbrace{\omega_{\mathbf{s}}}_{\substack{\text{spatial}\\\text{field}}}
+
\underbrace{\epsilon_{\mathbf{s},t}}_{\substack{\text{spatiotemporal}\\\text{field}}},
\end{aligned}
\]
where
- \(y_{\mathbf{s},t}\) represents the
response data at point \(\mathbf{s}\)
and time \(t\);
- \(\mathcal{D}\) represents the
observation distribution (family) with parameter \(\mu_{\mathbf{s},t}\), dispersion parameter
\(\phi\) where applicable, and, for
some families, additional parameters (e.g., the Tweedie power
parameter);
- \(\mu = f^{-1}(\eta)\) represents
the family’s main parameter, which for most families is the conditional
mean of the response (i.e., \(\mathbb{E}[y_{\mathbf{s},t} \mid \text{random
effects}] = \mu_{\mathbf{s},t}\));
- \(f\) represents a link function
(e.g., log or logit) and \(f^{-1}\)
represents its inverse;
- \(\eta\) represents the linear
predictor (i.e., in link space, before applying \(f^{-1}\));
- \(\mathbf{X}^{\mathrm{main}}\),
\(\mathbf{X}^{\mathrm{tvc}}\), and
\(\mathbf{X}^{\mathrm{svc}}\) represent
design matrices (the superscript identifiers ‘main’ = main effects,
‘tvc’ = time-varying coefficients, and ‘svc’ = spatially varying
coefficients);
- \(\boldsymbol{\beta}\) represents a
vector of fixed-effect coefficients;
- \(O_{\mathbf{s},t}\) represents an
offset: a covariate with a coefficient fixed at one that enters on the
linear predictor scale (e.g., log effort with a log link);
- \(\mathbf{b}_{k,g}\) represents the
vector of random intercepts and/or slopes for level \(g\) of grouping factor \(k\) (with \(g\) being the level of factor \(k\) for this observation) and \(\mathbf{z}_{k,\mathbf{s},t}\) the
corresponding covariate values (a row of the random-effect design matrix
\(\mathbf{Z}\));
- \(\boldsymbol{\gamma}_{t}\)
represents a vector of time-varying coefficients, where each coefficient
\(\gamma_{p,t}\) can follow a random
walk or AR(1) process;
- \(\boldsymbol{\zeta}_{\mathbf{s}}\)
represents the values of the spatially varying coefficients at location
\(\mathbf{s}\), where each coefficient
field \(\boldsymbol{\zeta}_{l} \sim
\mathrm{MVNormal}(\boldsymbol{0},\boldsymbol{\Sigma}_{\zeta,l})\)
is a Gaussian random field;
- \(\omega_{\mathbf{s}}\) represents
a spatial random field, \(\boldsymbol{\omega}\sim
\mathrm{MVNormal}(\boldsymbol{0},\boldsymbol{\Sigma}_\omega)\),
with \(\boldsymbol{\Sigma}_\omega\)
approximating a Matérn covariance via the SPDE approach (described below); and
- \(\epsilon_{\mathbf{s},t}\)
represents a spatiotemporal random field with IID, AR(1), or random walk
dynamics through time; its spatial covariance is based on the Matérn
model, and its marginal variance depends on the temporal process.
For binomial and beta-binomial counts, \(\mu\) is a success probability and the
count mean is \(N\mu\). For truncated
negative binomials, ordered beta, and mixture families, \(\mu\) describes an underlying distribution
or component rather than the full response mean; the family sections
explain these distinctions. In the predictor equations, scalar field
terms such as \(\omega_{\mathbf{s}}\)
mean the field evaluated at the observation location. We use stars
explicitly when distinguishing mesh-vertex values from their
interpolated values (e.g., \(\boldsymbol{\omega}^* =
\mathbf{A}\boldsymbol{\omega}\)).
Penalized smoothers are included in \(\mathbf{X}^{\mathrm{main}}_{\mathbf{s},t}
\boldsymbol{\beta}\) for brevity (their penalized coefficients
are random effects; see Smoothers), and nonlocal (distributed lag)
terms are omitted from the overview. Delta (hurdle) families have
two linear predictors, each of which can contain the above components
(see Delta models).
A single sdmTMB model will rarely, if ever, contain all of the above
components.
Linear predictor components
This section describes each component of the linear predictor in more
detail, using ‘\(\ldots\)’ to represent
the other optional components.
Main effects
\[
\begin{aligned}
\mu_{\mathbf{s},t} &= f^{-1} \left(
\mathbf{X}^{\mathrm{main}}_{\mathbf{s},t} \boldsymbol{\beta}+ \ldots
\right)
\end{aligned}
\]
Within sdmTMB(), \(\mathbf{X}^{\mathrm{main}}_{\mathbf{s},t}
\boldsymbol{\beta}\) is defined by the formula
argument and represents the main-effect model matrix and a corresponding
vector of coefficients. This main effect formula can contain optional
penalized smoothers or non-linear functions as defined below.
Smoothers
Smoothers in sdmTMB are implemented with the same formula syntax
familiar to mgcv (Wood 2017) users fitting
GAMs (generalized additive models). Smooths are implemented in the
formula using + s(x), which implements a smooth from
mgcv::s(). Within these smooths, the same syntax commonly
used in mgcv::s() can be applied, e.g. 2-dimensional
smooths may be constructed with + s(x, y); smooths can be
specific to various factor levels, + s(x, by = group);
smooths can vary according to a continuous variable,
+ s(x, by = x2); the basis function dimensions may be
specified, e.g. + s(x, k = 4) (see
?mgcv::choose.k); and various types of splines may be
constructed such as cyclic splines to model seasonality,
e.g. + s(month, bs = "cc", k = 12).
sdmTMB supports penalized smooths, which balance fit to the data
against excessive wiggliness. The function
mgcv::smooth2random() represents each supported smooth
using unpenalized coefficients and penalized coefficients treated as
random effects, with associated design matrices. This allows the smooth
to be estimated in a mixed-effects modelling framework. This is the same
approach as is implemented in the R packages gamm4 (Wood & Scheipl 2020) and brms (Bürkner 2017).
Linear break-point threshold models
The linear break-point or “hockey stick” model can be used to
describe threshold or asymptotic responses. This function consists of
two pieces, so that for \(x < c\),
\(s(x) = x \cdot \lambda\), and for
\(x \ge c\), \(s(x) = c \cdot \lambda\). Here, \(\lambda\) represents the slope of the
function up to the threshold (break point) \(c\), and the product \(c \cdot \lambda\) represents the value at
the asymptote. No constraints are placed on \(\lambda\) or \(c\).
These models can be fit by including + breakpt(x) in the
model formula, where x is a covariate. The formula can
contain a single break-point covariate.
Logistic threshold models
Models with logistic threshold relationships between a predictor and
the response can be fit with the form
\[
s(x)=\tau + \psi\ { \left[ 1+{ e }^{ -\ln \left(19\right) \cdot \left(
x-s50 \right)
/ \left(s95 - s50 \right) } \right] }^{-1},
\]
where \(s(x)\) denotes the threshold
function, \(\psi\) is a scaling
parameter (controlling the height and sign of the response curve;
unconstrained), \(\tau\) is an
intercept, \(s50\) is the value of
\(x\) at which the function has reached
50% of its total change (\(s(x) = \tau +
0.5\psi\)), and \(s95\) is the
value of \(x\) at which it has reached
95% of its total change (\(s(x) = \tau +
0.95\psi\)). The parameter \(s50\) is unconstrained but \(s95\) is constrained to be larger than
\(s50\).
These models can be fit by including + logistic(x) in
the model formula, where x is a covariate. The formula can
contain a single logistic covariate.
Offset terms
Offset terms can be included through the offset argument
in sdmTMB(). These are included in the linear predictor
as
\[
\begin{aligned}
\mu_{\mathbf{s},t} &= f^{-1} \left( \ldots + O_{\mathbf{s},t} +
\ldots \right),
\end{aligned}
\]
where \(O_{\mathbf{s},t}\) is an
offset term—a variable on the linear predictor scale with its
coefficient fixed at one. With a log link, the offset is usually a
log-transformed measure of effort (e.g.,
log(area_swept)), so that \(\mu_{\mathbf{s},t}\) is proportional to
effort.
When predicting, predict.sdmTMB() without
newdata uses the offset from the fitted model. With
newdata, the offset is 0 unless an offset
vector is supplied to predict.sdmTMB(). Therefore, if the
offset represents log effort, predictions on newdata are
for one unit of effort (log(1) = 0) by default.
Random intercepts and slopes
Multilevel (hierarchical) intercepts and slopes follow the familiar
lme4/glmmTMB formulation. Each grouping factor \(k\) can carry an intercept-only term or a
vector of random effects (intercept plus slopes) with a full covariance
matrix:
\[
\begin{aligned}
\mu_{\mathbf{s},t} &= f^{-1} \left( \ldots + \sum_{k}
\mathbf{z}_{k,\mathbf{s},t}^\top \mathbf{b}_{k,g} + \ldots \right),\\
\mathbf{b}_{k,g} &\sim \mathrm{MVNormal} \left( \boldsymbol{0},
\boldsymbol{\Sigma}_{k} \right),
\end{aligned}
\]
where \(g\) is the level of grouping
factor \(k\) for the observation at
\((\mathbf{s}, t)\), \(\mathbf{b}_{k,g}\) is the vector of random
effects for that level, and \(\mathbf{z}_{k,\mathbf{s},t}\) contains the
corresponding covariate values (1 for an intercept and the covariate
value for each slope). Random effects are independent across levels and
grouping factors. The scalar special case (1 | g) is a
random intercept \(b_{g} \sim
\mathrm{Normal}(0, \sigma_{b}^2)\).
Use lme4-style syntax in formula,
e.g. (1 | g) for random intercepts,
(1 + x | g) for correlated intercepts and slopes, or
(1 | g1) + (1 + x | g2) to give each grouping factor its
own (co)variance.
Internally, standard deviations are estimated on the log scale, and
correlations are represented by unconstrained parameters of a Cholesky
factor of the correlation matrix (TMB’s UNSTRUCTURED_CORR),
which guarantees a positive-definite \(\boldsymbol{\Sigma}_{k}\) for each grouping
factor.
Time-varying regression parameters
Parameters can be modelled as time-varying according to a random walk
or first-order autoregressive, AR(1), process. The time-series model is
defined by time_varying_type. For all types:
\[
\begin{aligned}
\mu_{\mathbf{s},t} &= f^{-1} \left( \ldots
+ \mathbf{X}^{\mathrm{tvc}}_{\mathbf{s},t} \boldsymbol{\gamma}_{t} +
\ldots \right),
\end{aligned}
\] where \(\boldsymbol{\gamma}_{t}\) is an optional
vector of time-varying regression parameters and \(\mathbf{X}^{\mathrm{tvc}}_{\mathbf{s},t}\)
is the corresponding model matrix with covariate values. This is defined
via the time_varying argument, assuming that the
time argument is also supplied a column name.
time_varying takes a one-sided formula.
~ 1 implies a time-varying intercept.
The following equations describe the time-series process for a single
time-varying coefficient \(\gamma_{p,t}\), where \(p\) indexes the time-varying coefficient.
When multiple time-varying coefficients are specified, each follows its
time-series process independently with its own variance parameter \(\sigma_{\gamma,p}^2\) and (for AR1)
correlation parameter \(\rho_{\gamma,p}\).
For time_varying_type = 'rw', the initial value is
unpenalized (equivalent to an improper flat prior); only changes between
successive time steps receive a Gaussian penalty:
\[
\begin{aligned}
\gamma_{p,t>1} &\sim \mathrm{Normal} \left(\gamma_{p,t-1},
\sigma_{\gamma,p}^2 \right).
\end{aligned}
\]
Because the initial value is estimated as part of the time-varying
formula, the same variable should not appear in the fixed-effects
formula. The formula time_varying = ~ 1 implicitly
represents a time-varying intercept (assuming the time
argument has been supplied) and, in this case, the intercept should be
omitted from the main-effects formula (e.g.,
formula = y ~ 0 + ... or
formula = y ~ -1 + ...).
For time_varying_type = 'rw0', the first time step is
estimated from a mean-zero prior:
\[
\begin{aligned}
\gamma_{p,t=1} &\sim \mathrm{Normal} \left(0, \sigma_{\gamma,p}^2
\right),\\
\gamma_{p,t>1} &\sim \mathrm{Normal} \left(\gamma_{p,t-1},
\sigma_{\gamma,p}^2 \right).
\end{aligned}
\] In this case, the time-varying variable (including the
intercept) should be included in the main effects.
For time_varying_type = 'ar1':
\[
\begin{aligned}
\gamma_{p,t=1} &\sim \mathrm{Normal} \left(0, \sigma_{\gamma,p}^2
\right),\\
\gamma_{p,t>1} &\sim \mathrm{Normal}
\left(\rho_{\gamma,p}\gamma_{p,t-1}, (1 - \rho_{\gamma,p}^2)
\sigma_{\gamma,p}^2 \right),
\end{aligned}
\] where \(\rho_{\gamma,p}\) is
the correlation between subsequent time steps. Here, \(\sigma_{\gamma,p}\) is the stationary
marginal standard deviation; the innovation standard deviation is \(\sigma_{\gamma,p}\sqrt{1 -
\rho_{\gamma,p}^2}\). As with rw0, include the
corresponding main effect to allow a non-zero average coefficient.
Spatial random fields
Spatial random fields, \(\omega_{\mathbf{s}}\), are included if
spatial = 'on' (or TRUE) and omitted if
spatial = 'off' (or FALSE).
\[
\begin{aligned}
\mu_{\mathbf{s},t} &= f^{-1} \left( \ldots + \omega_{\mathbf{s}} +
\ldots \right),\\
\boldsymbol{\omega}&\sim \mathrm{MVNormal} \left( \boldsymbol{0},
\boldsymbol{\Sigma}_\omega \right),
\end{aligned}
\]
where \(\boldsymbol{\omega}\) is the
vector of field values at the mesh vertices. The marginal standard
deviation of \(\boldsymbol{\omega}\),
\(\sigma_\omega\), is indicated by
Spatial SD in the printed model output or as
sigma_O in the output of
sdmTMB::tidy(fit, "ran_pars"). The ‘O’ is for ‘omega’
(\(\omega\)).
Internally, the random field is a Gaussian Markov random field (GMRF)
defined by a sparse precision (inverse covariance) matrix \(\mathbf{Q}_\omega\):
\[
\boldsymbol{\omega}\sim \mathrm{MVNormal}\left(\boldsymbol{0},
\mathbf{Q}^{-1}_\omega\right).
\]
\(\mathbf{Q}_\omega\) depends on the
mesh and on the parameters \(\kappa_\omega\) and \(\tau_\omega\), which determine the spatial
range and \(\sigma_\omega\) (see Matérn parameterization).
Spatiotemporal random fields
Spatiotemporal random fields are included by default if there are
multiple time elements (time argument is not
NULL) and can be set to IID (independent and identically
distributed, 'iid'; default), AR(1) ('ar1'),
random walk ('rw'), or off ('off') via the
spatiotemporal argument. These text values are case
insensitive.
Spatiotemporal random fields are represented by \(\boldsymbol{\epsilon}_t\) within sdmTMB.
This has been chosen to match the representation in VAST (Thorson 2019). The standard deviation parameter
\(\sigma_\epsilon\) is indicated by
Spatiotemporal SD in the printed model output or as
sigma_E in the output of
sdmTMB::tidy(fit, "ran_pars"). The ‘E’ is for ‘epsilon’
(\(\epsilon\)). For IID and AR(1)
fields, this is the marginal standard deviation; for a random walk, it
is the innovation standard deviation, as described below.
IID spatiotemporal random fields
IID spatiotemporal random fields
(spatiotemporal = 'iid') can be represented as
\[
\begin{aligned}
\mu_{\mathbf{s},t} &= f^{-1} \left( \ldots + \epsilon_{\mathbf{s},t}
+ \ldots \right),\\
\boldsymbol{\epsilon}_{t} &\sim \mathrm{MVNormal} \left(
\boldsymbol{0}, \boldsymbol{\Sigma}_{\epsilon} \right).
\end{aligned}
\]
where \(\epsilon_{\mathbf{s},t}\)
represents random field deviations at point \(\mathbf{s}\) and time \(t\). The random fields are assumed
independent across time steps.
Similarly to the spatial random fields, these IID spatiotemporal
random fields are parameterized internally with a sparse precision
matrix \(\mathbf{Q}_\epsilon\) (with
parameters \(\kappa_\epsilon\) and
\(\tau_\epsilon\))
\[
\boldsymbol{\epsilon}_{t} \sim \mathrm{MVNormal}\left(\boldsymbol{0},
\mathbf{Q}^{-1}_\epsilon\right).
\]
AR(1) spatiotemporal random fields
First-order autoregressive, AR(1), spatiotemporal random fields
(spatiotemporal = 'ar1') add a parameter defining the
correlation between random field deviations from one time step to the
next. They are defined as
\[
\begin{aligned}
\mu_{\mathbf{s},t} &= f^{-1} \left( \ldots + \epsilon_{\mathbf{s},t}
+ \ldots \right),\\
\boldsymbol{\epsilon}_{t=1} &\sim \mathrm{MVNormal} (\boldsymbol{0},
\boldsymbol{\Sigma}_{\epsilon}),\\
\boldsymbol{\epsilon}_{t>1} &= \rho \boldsymbol{\epsilon}_{t-1} +
\sqrt{1 - \rho^2} \boldsymbol{\delta}_{t}, \:
\boldsymbol{\delta}_{t} \sim \mathrm{MVNormal} \left(\boldsymbol{0},
\boldsymbol{\Sigma}_{\epsilon} \right),
\end{aligned}
\] where \(\rho\) is the
correlation between subsequent spatiotemporal random fields and \(\boldsymbol{\delta}_t\) are IID spatial
innovations. Because the first time step is drawn from the stationary
distribution and the \(\sqrt{1 -
\rho^2}\) factor scales the innovations, every time step has the
same marginal variance, \(\sigma_\epsilon^2\); the process is
stationary. The correlation \(\rho\)
allows for mean-reverting spatiotemporal fields, and is constrained to
be \(-1 < \rho < 1\). Internally,
the parameter is estimated as ar1_phi, which is
unconstrained. The parameter ar1_phi is transformed to
\(\rho\) with \(\rho = 2 \left(
\mathrm{logit}^{-1}(\texttt{ar1\_phi}) \right) - 1\), mapping
ar1_phi to the range \((-1,
1)\).
Random walk spatiotemporal random fields (RW)
Random walk spatiotemporal random fields
(spatiotemporal = 'rw') represent a model where the
difference in spatiotemporal deviations from one time step to the next
are IID. They are defined as
\[
\begin{aligned}
\mu_{\mathbf{s},t} &= f^{-1} \left( \ldots + \epsilon_{\mathbf{s},t}
+ \ldots \right),\\
\boldsymbol{\epsilon}_{t=1} &\sim \mathrm{MVNormal} (\boldsymbol{0},
\boldsymbol{\Sigma}_{\epsilon}),\\
\boldsymbol{\epsilon}_{t>1} &= \boldsymbol{\epsilon}_{t-1}
+ \boldsymbol{\delta}_{t}, \:
\boldsymbol{\delta}_{t} \sim \mathrm{MVNormal} \left(\boldsymbol{0},
\boldsymbol{\Sigma}_{\epsilon} \right),
\end{aligned}
\]
where the distribution of the spatiotemporal field in the initial
time step is the same as for the AR(1) model, but the absence of the
\(\rho\) parameter allows the
spatiotemporal field to be non-stationary in time. A random walk has no
stationary variance: under this initialization, \(\mathrm{Cov}(\boldsymbol{\epsilon}_t) =
t\boldsymbol{\Sigma}_\epsilon\), so the marginal variance grows
with each time step. Therefore, \(\sigma_\epsilon\) is the standard deviation
of the innovations \(\boldsymbol{\delta}_t\) (and of the field
in the first time step), not the marginal standard deviation of the
field in later time steps.
Spatially varying coefficients (SVC)
Spatially varying coefficient models are defined as
\[
\begin{aligned}
\mu_{\mathbf{s},t} &= f^{-1} \left( \ldots +
\mathbf{X}^{\mathrm{svc}}_{\mathbf{s}, t}
\boldsymbol{\zeta}_{\mathbf{s}} + \ldots \right),\\
\boldsymbol{\zeta}_{l} &\sim \mathrm{MVNormal} \left(
\boldsymbol{0}, \boldsymbol{\Sigma}_{\zeta,l} \right),
\end{aligned}
\]
where \(\zeta_{l,\mathbf{s}}\)
represents the value of the \(l\)-th
spatially varying coefficient field \(\boldsymbol{\zeta}_l\) at location \(\mathbf{s}\). When multiple spatially
varying coefficients are specified, each is estimated as an independent
spatial random field with its own marginal variance \(\sigma_{\zeta,l}^2\). \(\mathbf{X}^{\mathrm{svc}}_{\mathbf{s}, t}\)
represents the corresponding design matrix for the spatially varying
coefficient terms, as defined by a one-sided formula supplied to
spatial_varying. For example
spatial_varying = ~ 0 + x, where 0 omits the
intercept. Because \(\zeta_{l,\mathbf{s}}\) is a mean-zero
deviation, the main effect of x should usually also be
included in formula so that the coefficient at location
\(\mathbf{s}\) is \(\beta_x + \zeta_{l,\mathbf{s}}\).
The random fields are parameterized internally with a sparse
precision matrix \(\mathbf{Q}_{\zeta,l}\):
\[
\boldsymbol{\zeta}_{l} \sim \mathrm{MVNormal}\left(\boldsymbol{0},
\mathbf{Q}^{-1}_{\zeta,l}\right).
\]
Each \(\mathbf{Q}_{\zeta,l}\) has
its own \(\tau_{\zeta,l}\) (and
therefore \(\sigma_{\zeta,l}\)) but
shares \(\kappa\) with the spatial
random field \(\boldsymbol{\omega}\).
Nonlocal formula (covariate diffusion)
Nonlocal (distributed lag) terms let the effect of a covariate be
spread over nearby locations and/or previous time steps, rather than
acting only where and when the covariate was measured. This can be
useful when ecological responses are delayed, transported, or
accumulated through space and time. They are added with the
nonlocal_formula argument to sdmTMB() as a
one-sided formula using diffusion() (space) and/or
time_lag() (time) wrappers, e.g.,
~ diffusion(x) + time_lag(x). time_lag() terms
require a time argument. See (and cite) Lindmark et al. (2025) for the spatial
diffusion model and Thorson et al.
(2026) for the extension to time; those papers give the full
formulation and guidance on interpretation.
Each covariate \(x_r\) in
nonlocal_formula is transformed into a smoothed covariate
\(\tilde{x}_r\), which enters the
linear predictor with an estimated coefficient:
\[
\mu_{\mathbf{s},t} = f^{-1} \left(\ldots + \sum_{r}
\beta^{\mathrm{nl}}_{r}\,\tilde{x}_{r,\mathbf{s},t} + \ldots \right).
\]
If the same covariate appears in both wrappers, as in
~ diffusion(x) + time_lag(x), both apply to one transformed
covariate with one coefficient. If the covariates differ, as in
~ diffusion(x1) + time_lag(x2), the model has two
transformed covariates and two coefficients.
For diffusion(x), the covariate at time \(t\), \(\mathbf{x}_t\) (evaluated at the mesh
vertices), is smoothed in space by solving
\[
\left(\mathbf{C} +
\kappa_{S}^{-2}\mathbf{G}\right)\tilde{\mathbf{x}}_{t} =
\mathbf{C}\mathbf{x}_{t},
\]
where \(\mathbf{C}\) and \(\mathbf{G}\) are the SPDE mass and
stiffness matrices (see Matérn
parameterization) and \(\kappa_S^{-1}\) controls the spatial scale
of the diffusion. The output reports this scale as RMSDK
(the root-mean-squared displacement of the diffusion kernel, \(2/\kappa_S\), in the units of the mesh
coordinates).
For time_lag(x), the covariate at each vertex is an
exponentially weighted average of current and past values:
\[
\tilde{x}_{\mathbf{s},t} = (1 - \rho_T)\,x_{\mathbf{s},t} +
\rho_T\,\tilde{x}_{\mathbf{s},t-1}, \qquad \rho_T = \frac{\kappa_T}{1 +
\kappa_T} \in [0, 1),
\]
By default, time_lag(x) uses
start = "stationary", assuming that the covariate held its
first observed value before the series began: \(\tilde{x}_{\mathbf{s},0} =
x_{\mathbf{s},1}\). Use time_lag(x, start = "zero")
for \(\tilde{x}_{\mathbf{s},0} = 0\).
The choice affects the transformed covariate near the start of the
series. Larger \(\rho_T\) means a
longer memory of past covariate values.
When both wrappers are applied to the same covariate, spatial
smoothing and temporal propagation act on the same evolving field (this
corresponds to kappaST = 0 in Thorson et al. (2026)), and the
reported MSDK and RMSDK account for both. With
both wrappers, the stationary initial value is the spatially smoothed
first covariate slice. A separate \(\kappa_S\) is estimated for each covariate
with spatial diffusion, and a separate \(\kappa_T\) for each covariate with a
temporal lag.
Gaussian random fields
This section describes how the spatial (\(\boldsymbol{\omega}\)), spatiotemporal
(\(\boldsymbol{\epsilon}_t\)), and
spatially varying coefficient (\(\boldsymbol{\zeta}_l\)) random fields are
constructed. How each field enters the model is described under Linear predictor components.
Interpreting range and standard deviation
A spatial random field is a collection of random effects whose values
are correlated across locations. For the usual isotropic Matérn model,
correlation depends on distance, with nearby locations more strongly
correlated than distant ones. The range describes the
distance over which this correlation declines: at one range, correlation
is approximately 0.13, not zero. The marginal standard
deviation describes the typical size of field deviations at a
location, on the link scale. A larger range produces broader spatial
patterns; a larger standard deviation allows larger departures from the
predictor effects. For spatiotemporal fields, the reported standard
deviation describes the field itself under IID or AR(1) dynamics, but
the changes between time steps under a random walk.
In a model with both spatial and spatiotemporal fields,
share_range = TRUE (the default) estimates a common range
with separate standard deviations for the two fields. Sharing a range
can be useful when the data contain limited information about separate
ranges. Use share_range = FALSE to estimate separate
ranges. Spatially varying coefficient fields share the spatial
range.
The following construction details are optional for understanding the
model components and interpreting these parameters.
Matérn parameterization
For continuous-space models without barriers, sdmTMB approximates
Gaussian random fields with Matérn spatial covariance. The Matérn
covariance between the field values at locations \(\mathbf{s}_j\) and \(\mathbf{s}_k\), separated by distance \(d_{jk} = \lVert \mathbf{s}_j - \mathbf{s}_k
\rVert\), is
\[
\mathrm{Cov}\left( \omega_{\mathbf{s}_j}, \omega_{\mathbf{s}_k} \right)
= \frac{\sigma^2}{2^{\nu - 1}\Gamma(\nu)} (\kappa d_{jk})^\nu K_\nu
\left( \kappa d_{jk} \right),
\]
where \(\sigma^2\) is the marginal
variance, \(\nu\) controls the
smoothness of the field, \(\Gamma\) is
the gamma function, \(K_\nu\) is the
modified Bessel function of the second kind, and \(\kappa\) is a scaling parameter (larger
\(\kappa\) means correlation decays
faster with distance).
Computing with a dense Matérn covariance matrix is slow for many
locations. Instead, sdmTMB uses the stochastic partial differential
equation (SPDE) approach of Lindgren et
al. (2011). A Gaussian field with Matérn covariance is the
solution to the SPDE \((\kappa^2 -
\Delta)^{\alpha/2}(\tau\,\omega(\mathbf{s})) =
\mathcal{W}(\mathbf{s})\), where \(\Delta\) is the Laplacian, \(\mathcal{W}\) is Gaussian white noise, and
\(\alpha = \nu + d/2\) for \(d\) spatial dimensions. sdmTMB uses \(d = 2\) and \(\alpha = 2\), so \(\nu = 1\). Approximating the solution on a
triangulated mesh with piecewise-linear basis functions (the finite
element method) gives field values at the mesh vertices with a sparse
precision matrix
\[
\mathbf{Q}= \tau^2 \left( \kappa^4 \mathbf{C} + 2 \kappa^2 \mathbf{G} +
\mathbf{G} \mathbf{C}^{-1} \mathbf{G} \right),
\]
where \(\mathbf{C}\) is a diagonal
approximation to the finite-element mass matrix (called mass lumping),
and \(\mathbf{G}\) is a sparse
stiffness matrix. Both matrices depend only on the mesh. The sparsity of
\(\mathbf{Q}\) (most entries are zero)
is what makes the computations fast.
The parameters \(\kappa\) and \(\tau\) determine the reported range and
marginal standard deviation:
\[
\textrm{range} = \frac{\sqrt{8}}{\kappa}, \qquad
\sigma = \frac{1}{\sqrt{4 \pi}\, \tau \kappa}.
\]
These relationships describe the stationary Matérn field; the finite
mesh and its boundaries introduce approximation error. For a fixed \(\tau\), increasing \(\kappa\) decreases \(\sigma\). Sharing a range means sharing
\(\kappa\), while separate \(\tau\) parameters allow different field
standard deviations. Throughout, \(\mathbf{Q}_\omega\), \(\mathbf{Q}_\epsilon\), and \(\mathbf{Q}_{\zeta,l}\) denote full
precision matrices, including the \(\tau^2\) scaling (the code applies this
scaling separately when evaluating the field densities).
Projection \(\mathbf{A}\)
matrix
In the SPDE approach, the random field at any location is defined by
the piecewise-linear basis expansion \(\omega(\mathbf{s}) = \sum_k \psi_k(\mathbf{s})\,
\omega_k\), where \(\omega_k\)
is the field value at mesh vertex \(k\)
and \(\psi_k\) is a “tent” function
that is 1 at vertex \(k\) and decreases
linearly to 0 at the neighbouring vertices. Evaluating the field at
observed or prediction locations is therefore a multiplication by a
sparse projection matrix \(\mathbf{A}\)
with elements \(A_{ik} =
\psi_k(\mathbf{s}_i)\) (Lindgren & Rue
2015):
\[
\boldsymbol{\omega}^* = \mathbf{A}\boldsymbol{\omega},
\] where \(\boldsymbol{\omega}^*\) represents the
values of the spatial random fields at the observed locations or
predicted data locations. The matrix \(\mathbf{A}\) has a row for each data point
or prediction point and a column for each mesh vertex. Each row has at
most three non-zero elements: the barycentric weights of the vertices of
the triangle containing that location. This interpolation is part of the
model definition rather than a separate approximation step. The same
interpolation happens for any spatiotemporal random fields
\[
\boldsymbol{\epsilon}_t^* = \mathbf{A}\boldsymbol{\epsilon}_t.
\]
Anisotropy
TMB allows for anisotropy, where spatial covariance may depend on
direction (full
details). Anisotropy can be turned on or off with the logical
anisotropy argument to sdmTMB(). There are a
number of ways to implement anisotropic covariance (Fuglstad et al. 2015), and we adopt
geometric anisotropy defined by a 2-parameter matrix \(\mathbf{H}\), which replaces Euclidean
distance with \(d_{jk} = \sqrt{(\mathbf{s}_j -
\mathbf{s}_k)^\top \mathbf{H}^{-1} (\mathbf{s}_j -
\mathbf{s}_k)}\). The elements of \(\mathbf{H}\) are defined by the parameters
\(h_1 > 0\) and \(h_2\) so that \(H_{1,1} = h_{1}\), \(H_{1,2} = H_{2,1} = h_{2}\) and \(H_{2,2} = (1 + h_{2}^2) / h_{1}\). \(\mathbf{H}\) is symmetric and positive
definite with determinant 1, so it stretches distances along one axis
and shrinks them along the perpendicular axis without changing area. The
result is a different range in each direction, with the orientation of
the major axis estimated.
Once a model is fitted with sdmTMB(), the anisotropy
relationships may be plotted using the plot_anisotropy()
function, which takes the fitted object as an argument. If a barrier
mesh is used, anisotropy is disabled.
Incorporating physical barriers into the SPDE
In some cases the spatial domain of interest may be complex and
bounded by some barrier such as by land or water (e.g., coastlines,
islands, lakes). SPDE models allow for physical barriers to be
incorporated into the modelling (Bakka et
al. 2019). With sdmTMB() models, the mesh
construction occurs in two steps: the user (1) constructs a mesh with a
call to sdmTMB::make_mesh(), and (2) passes the mesh to
sdmTMBextra::add_barrier_mesh(). The barriers must be
constructed as sf objects (Pebesma
2018) with polygons defining the barriers. See
?sdmTMBextra::add_barrier_mesh for an example.
The barrier implementation requires the user to select a fraction
value (range_fraction argument) that defines the fraction
of the usual spatial range when crossing the barrier (Bakka et al. 2019). For example, if
the range was estimated at 10 km, range_fraction = 0.2
would assign a 2 km range to the triangles marked as barrier. That makes
correlation decay much faster through the barrier than around it, which
effectively discourages spatial smoothing from “jumping” across
landmasses. From experimentation, values around 0.1 or 0.2 seem to work
well but values much lower than 0.1 can result in convergence issues.
The range_fraction is fixed, not estimated.
Following Bakka et al. (2019),
the barrier model lets the range vary in space through the SPDE
\[
\omega(\mathbf{s}) - \nabla \cdot \left( \frac{r(\mathbf{s})^2}{8}
\nabla \omega(\mathbf{s}) \right) = r(\mathbf{s}) \sqrt{\frac{\pi}{2}}\,
\sigma\, \mathcal{W}(\mathbf{s}),
\]
where \(\nabla\) is the gradient,
\(r(\mathbf{s}) = r\) in the normal
part of the domain and \(r(\mathbf{s}) = r
\times\) range_fraction in barrier triangles. With a
constant range \(r = \sqrt{8}/\kappa\),
this is the Matérn SPDE given above, with marginal standard deviation
\(\sigma\). With barriers, the field is
no longer a stationary Matérn field, so its precision matrix is not the
\(\mathbf{Q}\) given above. sdmTMB
builds it with the formulation used in the INLAspacetime
package, parameterized directly by the range and marginal standard
deviation. The reported range is still \(\sqrt{8}/\kappa\), and \(\sigma\) is still computed from \(\kappa\) and \(\tau\) as above. Both describe the field
away from barriers; near barriers, the variance and correlation are
modified.
This
website by Francesco Serafini and Haakon Bakka provides an
illustration with INLA. The original TMB implementation in sdmTMB was
adapted from code written by Olav Nikolai Breivik and Hans Skaug at the
TMB Case Studies
GitHub site.
Areal autoregressive models (SAR/CAR)
In addition to continuous-space SPDE random fields, sdmTMB can fit
areal spatial autoregressive models with
spatial_model = "sar" or
spatial_model = "car". These require an areal domain
supplied to mesh with the help of
make_areal_domain().
For both areal options, the latent spatial and spatiotemporal effects
are Gaussian Markov random fields with precision defined by the areal
weight matrix \(\mathbf{W}\). For the
spatial field and IID or AR(1) spatiotemporal fields,
\[
\begin{aligned}
\boldsymbol{\omega}&\sim \mathrm{MVNormal}\left(\boldsymbol{0},
\sigma_\omega^2 \mathbf{Q}^{-1}\right),\\
\boldsymbol{\epsilon}_t &\sim \mathrm{MVNormal}\left(\boldsymbol{0},
\sigma_\epsilon^2 \mathbf{Q}^{-1}\right).
\end{aligned}
\]
For a random walk, \(\sigma_\epsilon^2
\mathbf{Q}^{-1}\) instead describes the covariance of the initial
field and subsequent innovations. Unlike the SPDE case, \(\sigma_\omega\) and \(\sigma_\epsilon\) here are scale
parameters, not marginal standard deviations: the diagonal entries of
\(\sigma^2\mathbf{Q}^{-1}\) depend on
\(\mathbf{W}\) and the autocorrelation
parameter and generally differ among areas.
For SAR (spatial_model = "sar"), the precision matrix
is
\[
\mathbf{Q}_{\mathrm{SAR}} = \left(\mathbf{I} -
\rho_{\mathrm{SAR}}\mathbf{W}\right)^\top
\left(\mathbf{I} - \rho_{\mathrm{SAR}}\mathbf{W}\right),
\]
where \(\rho_{\mathrm{SAR}} \in (-1,
1)\) is reported as rho_sar. By default, SAR uses a
row-normalized \(\mathbf{W}\)
(sar_weight_style = "row"), with an option to use raw
weights. With sar_weight_style = "raw", SAR uses the
unnormalized adjacency/weight matrix directly, so neighbour influence
depends on the original edge weights (and, if unweighted, on the number
of neighbours) rather than being scaled to sum to 1 within each row.
For CAR (spatial_model = "car"), the precision matrix
is
\[
\mathbf{Q}_{\mathrm{CAR}} = \mathbf{D} -
\alpha_{\mathrm{CAR}}\mathbf{W},
\]
where \(\mathbf{D}\) is diagonal
with elements \(D_{ii} = \sum_j
W_{ij}\) (the number of neighbours for unweighted adjacency; set
to 1 for areas with no neighbours) and \(\alpha_{\mathrm{CAR}} \in (0, 1)\) is
reported as alpha_car. CAR requires a symmetric adjacency
matrix.
Delta models
sdmTMB allows for several different kinds of delta (also known as
hurdle) models. These families are implemented by specifying the family
as a delta distribution. For example:
sdmTMB(
...,
family = delta_gamma()
)
The list of supported families is included in the documentation on additional
families, delta
models, and Poisson-link
delta models. By default, the delta_* families don’t
use the Poisson link, but this structure can be specified with
delta_gamma(type = "poisson-link"), which follows the
formulation in Thorson (2018).
In the “standard” delta model implementation, sdmTMB constructs two
internal models (formally, two “linear predictors”), with the first
model representing presence-absence and the second model representing
the positive component (such as catch rates in fisheries applications).
The positive component can be continuous (e.g.,
delta_gamma(), delta_lognormal(),
delta_gengamma()), a proportion in \((0, 1)\) (delta_beta()), or a
positive count from a zero-truncated negative binomial
(delta_truncated_nbinom2(),
delta_truncated_nbinom1()). The expected value is the
product of the probability of a non-zero observation and the expected
value of the positive component. Default links associated with each
family can be inspected with delta_lognormal() and
equivalent functions. For the standard delta models, the default first
linear predictor link is logit. For the Poisson-link type, the first
linear predictor link is log.
The formula, spatial,
spatiotemporal, and share_range arguments of
sdmTMB() can be specified independently as a 2-element
list. For example, the spatial random field might be estimated for only
the first linear predictor with:
sdmTMB(
...,
family = delta_gamma(),
spatial = list("on", "off")
)
Or we may want separate main-effects formulas. For example:
sdmTMB(
formula = list(
y ~ depth + I(depth^2),
y ~ 1),
family = delta_gamma()
)
All other arguments are shared between the linear predictors.
Currently if delta models contain smoothers, both components must
have the same main-effects formula.