Main effects

\[ \begin{aligned} \mu_{\boldsymbol{s},t} &= f^{-1} \left( \boldsymbol{X}^{\mathrm{main}}_{\boldsymbol{s},t} \boldsymbol{\beta} \ldots \right) \end{aligned} \]

Within sdmTMB(), \(\boldsymbol{X}^{\mathrm{main}}_{\boldsymbol{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.

Spatial random fields

Spatial random fields, \(\omega_{\boldsymbol{s}}\), are included if spatial = 'on' (or TRUE) and omitted if spatial = 'off' (or FALSE).

\[ \begin{aligned} \mu_{\boldsymbol{s},t} &= f^{-1} \left( \ldots + \omega_{\boldsymbol{s}} + \ldots \right),\\ \boldsymbol{\omega} &\sim \operatorname{MVNormal} \left( \boldsymbol{0}, \boldsymbol{\Sigma}_\omega \right),\\ \end{aligned} \] The marginal standard deviation of \(\boldsymbol{\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 fields follow a Gaussian Markov random field (GMRF)

\[ \boldsymbol{\omega} \sim \mathrm{MVNormal}\left(\boldsymbol{0}, \sigma_\omega^2 \boldsymbol{Q}^{-1}_\omega\right), \] where \(\boldsymbol{Q}_\omega\) is a sparse precision matrix and \(\sigma_\omega^2\) is the marginal variance.

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 marginal standard deviation of \(\boldsymbol{\epsilon}_t\) 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\)).

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_{\boldsymbol{s},t} &= f^{-1} \left( \ldots + \boldsymbol{X}^{\mathrm{tvc}}_{\boldsymbol{s},t} \boldsymbol{\gamma_{t}} + \ldots \right), \end{aligned} \] where \(\boldsymbol{\gamma_t}\) is an optional vector of time-varying regression parameters and \(\boldsymbol{X}^{\mathrm{tvc}}_{\boldsymbol{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.

For time_varying_type = 'rw', the first time step is estimated independently:

\[ \begin{aligned} \gamma_{t=1} &\sim \operatorname{Uniform} \left(-\infty, \infty \right),\\ \gamma_{t>1} &\sim \operatorname{Normal} \left(\gamma_{t-1}, \sigma^2_{\gamma} \right). \end{aligned} \]

In this case, the first time-step value is given an implicit uniform prior. I.e., the same variable should not appear in the fixed effect formula since the initial value is estimated as part of the time-varying formula. The formula time_varying = ~ 1 implicitly represents a time-varying intercept (assuming the time argument has been supplied) and, this case, the intercept should be omitted from the main effects (formula ~ + 0 + ... or formula ~ -1 + ...).

For time_varying_type = 'rw0', the first time step is estimated from a mean-zero prior:

\[ \begin{aligned} \gamma_{t=1} &\sim \operatorname{Normal} \left(0, \sigma^2_{\gamma} \right),\\ \gamma_{t>1} &\sim \operatorname{Normal} \left(\gamma_{t-1}, \sigma^2_{\gamma} \right). \end{aligned} \] In this case, the time-varying variable (including the intercept) should be included in the main effects. We suggest using this formulation, but leave the 'rw' option so that legacy code works.

For time_varying_type = 'ar1':

\[ \begin{aligned} \gamma_{t=1} &\sim \operatorname{Normal} \left(0, \sigma^2_{\gamma} \right),\\ \gamma_{t>1} &\sim \operatorname{Normal} \left(\rho_\gamma\gamma_{t-1}, \sqrt{1 - \rho_\gamma^2} \sigma^2_{\gamma} \right), \end{aligned} \] where \(\rho_{\gamma}\) is the correlation between subsequent time steps. The first time step is given a mean-zero prior.

Spatially varying coefficients (SVC)

Spatially varying coefficient models are defined as

\[ \begin{aligned} \mu_{\boldsymbol{s},t} &= f^{-1} \left( \ldots + \boldsymbol{X}^{\mathrm{svc}}_{\boldsymbol{s}, t} \zeta_{\boldsymbol{s}} + \ldots \right),\\ \boldsymbol{\zeta} &\sim \operatorname{MVNormal} \left( \boldsymbol{0}, \boldsymbol{\Sigma}_\zeta \right), \end{aligned} \]

where \(\boldsymbol{\zeta}\) is a random field representing a spatially varying coefficient. Usually, \(\boldsymbol{X}^{\mathrm{svc}}_{\boldsymbol{s}, t}\) would represent a prediction matrix that is constant spatially for a given time \(t\) as defined by a one-sided formula supplied to spatial_varying. For example spatial_varying = ~ 0 + x, where 0 omits the intercept.

The random fields are parameterized internally with a sparse precision matrix (\(\boldsymbol{Q}_\zeta\))

\[ \boldsymbol{\zeta} \sim \mathrm{MVNormal}\left(\boldsymbol{0}, \sigma_\zeta^2 \boldsymbol{Q}^{-1}_\zeta\right). \]

IID random or multi-level intercepts

Multilevel/hierarchical intercepts are defined as

\[ \begin{aligned} \mu_{\boldsymbol{s},t} &= f^{-1} \left( \ldots + \alpha_{g} + \ldots \right),\\ \alpha_g &\sim \operatorname{Normal} \left(0, \sigma_\alpha^2 \right),\\ \end{aligned} \]

where \(\alpha_g\) is an example optional “random” intercept—an intercept with mean zero that varies by level \(g\) and is constrained by \(\sigma_\alpha\). This is defined by the formula argument via the (1 | g) syntax as in lme4 or glmmTMB. There can be multiple random intercepts, despite only showing one above. E.g., (1 | g1) + (1 | g2), in which case they are assumed independent and uncorrelated from each other.

Offset terms

Offset terms can be included through the offset argument in sdmTMB(). These are included in the linear predictor as

\[ \begin{aligned} \mu_{\boldsymbol{s},t} &= f^{-1} \left( \ldots + O_{\boldsymbol{s},t} + \ldots \right), \end{aligned} \]

where \(O_{\boldsymbol{s},t}\) is an offset term—a log transformed variable without a coefficient (assuming a log link). The offset is not included in the prediction. Therefore, if offset represents a measure of effort, for example, the prediction is for one unit of effort (log(1) = 0).

Binomial

\[ \operatorname{Binomial} \left(N, \mu \right) \] where \(N\) is the size or number of trials, and \(\mu\) is the probability of success for each trial. If \(N = 1\), the distribution becomes the Bernoulli distribution. Internally, the distribution is parameterized as the robust version in TMB, which is numerically stable when probabilities approach 0 or 1. Following the structure of stats::glm(), lme4, and 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 cbind(success, failure), or
4. the response may be the observed proportions, and the weights argument is used to specify the Binomial size (\(N\)) parameter (probabilty ~ ..., weights = N).

Code defined within TMB.

Example: family = binomial(link = "logit")

Beta

\[ \operatorname{Beta} \left(\mu \phi, (1 - \mu) \phi \right) \] where \(\mu\) is the mean and \(\phi\) is a precision parameter. This parametrization follows Ferrari & Cribari-Neto (2004) and the betareg R package (Cribari-Neto & Zeileis 2010). The variance is \(\mu (1 - \mu) / (\phi + 1)\).

Code defined within TMB.

Example: family = Beta(link = "logit")

Gamma

\[ \operatorname{Gamma} \left( \phi, \frac{\mu}{\phi} \right) \] where \(\phi\) represents the Gamma shape and \(\mu / \phi\) represents the scale. The mean is \(\mu\) and variance is \(\mu \cdot \phi^2\).

Code defined within TMB.

Example: family = Gamma(link = "log")

Gaussian

\[ \operatorname{Normal} \left( \mu, \phi^2 \right) \] where \(\mu\) is the mean and \(\phi\) is the standard deviation. The variance is \(\phi^2\).

Example: family = Gaussian(link = "identity")

Code defined within TMB.

Lognormal

sdmTMB uses the bias-corrected lognormal distribution where \(\phi\) represents the standard deviation in log-space:

\[ \operatorname{Lognormal} \left( \log \mu - \frac{\phi^2}{2}, \phi^2 \right). \] Because of the bias correction, \(\mathbb{E}[y] = \mu\) and \(\mathrm{Var}[\log y] = \phi^2\).

Code defined within sdmTMB based on the TMB dnorm() normal density.

Example: family = lognormal(link = "log")

Negative Binomial 1 (NB1)

\[ \operatorname{NB1} \left( \mu, \phi \right) \]

where \(\mu\) is the mean and \(\phi\) is the dispersion parameter. The variance scales linearly with the mean \(\mathrm{Var}[y] = \mu + \mu / \phi\) (Hilbe 2011). Internally, the distribution is parameterized as the robust version in TMB.

Code defined within sdmTMB based on NB2 and borrowed from glmmTMB.

Example: family = nbinom1(link = "log")

Negative Binomial 2 (NB2)

\[ \operatorname{NB2} \left( \mu, \phi \right) \]

where \(\mu\) is the mean and \(\phi\) is the dispersion parameter. The variance scales quadratically with the mean \(\mathrm{Var}[y] = \mu + \mu^2 / \phi\) (Hilbe 2011). The NB2 parametrization is more commonly seen in ecology than the NB1. Internally, the distribution is parameterized as the robust version in TMB.

Code defined within TMB.

Example: family = nbinom2(link = "log")

Poisson

\[ \operatorname{Poisson} \left( \mu \right) \] where \(\mu\) represents the mean and \(\mathrm{Var}[y] = \mu\).

Code defined within TMB.

Example: family = poisson(link = "log")

Student-t

\[ \operatorname{Student-t} \left( \mu, \phi, \nu \right) \]

where \(\nu\), the degrees of freedom (df), is a user-supplied fixed parameter. Lower values of \(\nu\) result in heavier tails compared to the Gaussian distribution. Above approximately df = 20, the distribution becomes very similar to the Gaussian. The Student-t distribution with a low degrees of freedom (e.g., \(\nu \le 7\)) can be helpful for modelling data that would otherwise be suitable for Gaussian but needs an approach that is robust to outliers (e.g., Anderson et al. 2017).

Code defined within sdmTMB based on the dt() distribution in TMB.

Example: family = student(link = "log", df = 7)

Tweedie

\[ \operatorname{Tweedie} \left(\mu, p, \phi \right), \: 1 < p < 2 \]

where \(\mu\) is the mean, \(p\) is the power parameter constrained between 1 and 2, and \(\phi\) is the dispersion parameter. The Tweedie distribution can be helpful for modelling data that are positive and continuous but also contain zeros.

Internally, \(p\) is transformed from \(\mathrm{logit}^{-1} (\texttt{thetaf}) + 1\) to constrain it between 1 and 2 and is estimated as an unconstrained variable.

The source code is implemented as in the cplm package (Zhang 2013) and is based on Dunn & Smyth (2005). The TMB version is defined here.

Example: family = tweedie(link = "log")

Gamma mixture

This is a 2 component mixture that extends the Gamma distribution,

\[ (1 - p) \cdot \operatorname{Gamma} \left( \phi, \frac{\mu_{1}}{\phi} \right) + p \cdot \operatorname{Gamma} \left( \phi, \frac{\mu_{2}}{\phi} \right), \] where \(\phi\) represents the Gamma shape, \(\mu_{1} / \phi\) represents the scale for the first (smaller component) of the distribution, \(\mu_{2} / \phi\) represents the scale for the second (larger component) of the distribution, and \(p\) controls the contribution of each component to the mixture (also interpreted as the probability of larger events).

The mean is \((1-p) \cdot \mu_{1} + p \cdot \mu_{2}\) and the variance is \((1-p) ^ 2 \cdot \mu_{1} \cdot \phi^2 + (p) ^ 2 \cdot \mu_{2} \cdot \phi^2\).

Here, and for the other mixture distributions, the probability of the larger mean can be obtained from plogis(fit$model$par[["logit_p_mix"]]) and the ratio of the larger mean to the smaller mean can be obtained from 1 + exp(fit$model$par[["log_ratio_mix"]]). The standard errors are available in the TMB sdreport: fit$sd_report.

If you wish to fix the probability of a large (i.e., extreme) mean, which can be hard to estimate, you can use the map list. E.g.:

sdmTMB(...,
  control = sdmTMBcontrol(
    start = list(logit_p_mix = qlogis(0.05)), # 5% probability of 'mu2'
    map = list(logit_p_mix = factor(NA)) # don't estimate
  )
)

Example: family = gamma_mix(link = "log"). See also family = delta_gamma_mix() for an extension incorporating this distribution with delta models.

Lognormal mixture

This is a 2 component mixture that extends the lognormal distribution,

\[ (1 - p) \cdot \operatorname{Lognormal} \left( \log \mu_{1} - \frac{\phi^2}{2}, \phi^2 \right) + p \cdot \operatorname{Lognormal} \left( \log \mu_{2} - \frac{\phi^2}{2}, \phi^2 \right). \]

Because of the bias correction, \(\mathbb{E}[y] = (1-p) \cdot \mu_{1} + p \cdot \mu_{2}\) and \(\mathrm{Var}[\log y] = (1-p)^2 \cdot \phi^2 + p^2 \cdot \phi^2\).

As with the Gamma mixture, \(p\) controls the contribution of each component to the mixture (also interpreted as the probability of larger events).

Example: family = lognormal_mix(link = "log"). See also family = delta_lognormal_mix() for an extension incorporating this distribution with delta models.

Negative binomial 2 mixture

This is a 2 component mixture that extends the NB2 distribution,

\[ (1 - p) \cdot \operatorname{NB2} \left( \mu_1, \phi \right) + p \cdot \operatorname{NB2} \left( \mu_2, \phi \right) \]

where \(\mu_1\) is the mean of the first (smaller component) of the distribution, \(\mu_2\) is the mean of the larger component, and \(p\) controls the contribution of each component to the mixture.

Example: family = nbinom2_mix(link = "log")

Matérn parameterization

The Matérn defines the covariance \(\Phi \left( s_j, s_k \right)\) between spatial locations \(s_j\) and \(s_k\) as

\[ \Phi\left( s_j,s_k \right) = \tau^2/\Gamma(\nu)2^{\nu - 1} (\kappa d_{jk})^\nu K_\nu \left( \kappa d_{jk} \right), \]

where \(\tau^2\) controls the spatial variance, \(\nu\) controls the smoothness, \(\Gamma\) represents the Gamma function, \(d_{jk}\) represents the distance between locations \(s_j\) and \(s_k\), \(K_\nu\) represents the modified Bessel function of the second kind, and \(\kappa\) represents the decorrelation rate. The parameter \(\nu\) is set to 1 to take advantage of the Stochastic Partial Differential Equation (SPDE) approximation to the GRF to greatly increase computational efficiency (Lindgren et al. 2011). Internally, the parameters \(\kappa\) and \(\tau\) are converted to range and marginal standard deviation \(\sigma\) as \(\textrm{range} = \sqrt{8} / \kappa\) and \(\sigma = 1 / \sqrt{4 \pi \exp \left(2 \log(\tau) + 2 \log(\kappa) \right) }\).

In the case of a spatiotemporal model with both spatial and spatiotemporal fields, if share_range = TRUE in sdmTMB() (the default), then a single \(\kappa\) and range are estimated with separate \(\sigma_\omega\) and \(\sigma_\epsilon\). This often makes sense since data are often only weakly informative about \(\kappa\). If share_range = FALSE, then separate \(\kappa_\omega\) and \(\kappa_\epsilon\) are estimated. The spatially varying coefficient field always shares \(\kappa\) with the spatial random field.

Projection \(\boldsymbol{A}\) matrix

The values of the spatial variables at the knots are multiplied by a projection matrix \(\boldsymbol{A}\) that bilinearly interpolates from the knot locations to the values at the locations of the observed or predicted data (Lindgren & Rue 2015)

\[ \boldsymbol{\omega}^* = \boldsymbol{A} \boldsymbol{\omega}, \] where \(\boldsymbol{\omega}^*\) represents the values of the spatial random fields at the observed locations or predicted data locations. The matrix \(\boldsymbol{A}\) has a row for each data point or prediction point and a column for each knot. Three non-zero elements on each row define the weight of the neighbouring 3 knot locations for location \(\boldsymbol{s}\). The same bilinear interpolation happens for any spatiotemporal random fields

\[ \boldsymbol{\epsilon}_t^* = \boldsymbol{A} \boldsymbol{\epsilon}_t. \]

Anisotropy

TMB allows for anisotropy, where spatial covariance may be asymmetric with respect to latitude and longitude (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 a 2-parameter rotation matrix \(\textbf{H}\). The elements of \(\textbf{H}\) are defined by the parameter vector \(\boldsymbol{x}\) so that \(H_{1,1} = x_{1}\), \(H_{1,2} = H_{2,1} = x_{2}\) and \(H_{2,2} = (1 + x_{2}^2) / x_{1}\).

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 sdmTMB::add_barrier_mesh(). The barriers must be constructed as sf objects (Pebesma 2018) with polygons defining the barriers. See ?sdmTMB::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 assume that the range was 2 km across the barrier. This would let the spatial correlation decay 5 times faster with distance. 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.

This website by Francesco Serafini and Haakon Bakka provides an illustration with INLA. The implementation within TMB was borrowed from code written by Olav Nikolai Breivik and Hans Skaug at the TMB Case Studies Github site.

Optimization details

The sdmTMB model is fit by maximum marginal likelihood. Internally, a TMB (Kristensen et al. 2016) model template calculates the marginal log likelihood and its gradient, and the negative log likelihood is minimized via the non-linear optimization routine stats::nlminb() in R (Gay 1990; R Core Team 2021). Random effects are estimated at values that maximize the log likelihood conditional on the estimated fixed effects and are integrated over via the Laplace approximation (Kristensen et al. 2016).

Like AD Model Builder (Fournier et al. 2012), TMB allows for parameters to be fit in phases and we include the multiphase argument in sdmTMB::sdmTMBcontrol() to allow this. For high-dimensional models (many fixed and random effects), phased estimation may be faster and result in more stable convergence. In sdmTMB, phased estimation proceeds by first estimating all fixed-effect parameters contributing to the likelihood (holding random effects constant at initial values). In the second phase, the random-effect parameters (and their variances) are also estimated. Fixed-effect parameters are also estimated in the second phase and are initialized at their estimates from the first phase.

In some cases, a single call to stats::nlminb() may not be result in convergence (e.g., the maximum gradient of the marginal likelihood with respect to fixed-effect parameters is not small enough yet), and the algorithm may need to be run multiple times. In the sdmTMB::sdmTMBcontrol() function, we include an argument nlminb_loops that will restart the optimization at the previous best values. The number of nlminb_loops should generally be small (e.g., 2 or 3 initially), and defaults to 1. For some sdmTMB models, the Hessian may also be unstable and need to be re-evaluated. We do this optionally with the stats::optimHess() routine after the call to stats::nlminb(). The stats::optimHess() function implements a Newton optimization routine to find the Hessian, and we include the argument newton_loops in sdmTMB::sdmTMBcontrol() to allow for multiple function evaluations (each starting at the previous best value). By default, this is not included and newton_loops is set to 0. If a model is already fit, the function sdmTMB::run_extra_optimization() can run additional optimization loops with either routine to further reduce the maximum gradient.

Assessing convergence

Much of the guidance around diagnostics and glmmTMB also applies to sdmTMB, e.g. the glmmTMB vignette on troubleshooting. Optimization with stats::nlminb() involves specifying the number of iterations and evaluations (eval.max and iter.max) and the tolerances (abs.tol, rel.tol, x.tol, xf.tol)—a greater number of iterations and smaller tolerance thresholds increase the chance that the optimal solution is found, but more evaluations translates into longer computation time. Warnings of non-positive-definite Hessian matrices (accompanied by parameters with NAs for standard errors) often mean models are improperly specified given the data. Standard errors can be observed in the output of print.sdmTMB() or by checking fit$sd_report. The maximum gradient of the marginal likelihood with respect to fixed-effect parameters can be checked by inspecting (fit$gradients). Guidance varies, but the maximum gradient should likely be at least \(< 0.001\) before assuming the fitting routine is consistent with convergence. If maximum gradients are already relatively small, they can often be reduced further with additional optimization calls beginning at the previous best parameter vector as described above with sdmTMB::run_extra_optimization().