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

While mgcv can fit unpenalized (e.g., B-splines) or penalized splines (P-splines), sdmTMB only implements penalized splines. The penalized splines are constructed in sdmTMB using the function mgcv::smooth2random(), which transforms splines into random effects (and associated design matrices) that are estimable 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 < b_{1}\), \(s(x) = x \cdot b_{0}\), and for \(x > b_{1}\), \(s(x) = b_{1} \cdot b_{0}\). In both cases, \(b_{0}\) represents the slope of the function up to a threshold, and the product \(b_{1} \cdot b_{0}\) represents the value at the asymptote. No constraints are placed on parameters \(b_{0}\) or \(b_{1}\).

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\) represents the logistic function, \(\psi\) is a scaling parameter (controlling the height of the y-axis for the response; unconstrained), \(\tau\) is an intercept, \(s50\) is a parameter controlling the point at which the function reaches 50% of the maximum (\(\psi\)), and \(s95\) is a parameter controlling the point at which the function reaches 95% of the maximum. 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.

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}\) represent 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 spatiotemporal random fields (including all versions described below) are parameterized internally with a sparse precision matrix (\(\mathbf{Q}_\epsilon\))

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

AR(1) spatiotemporal random fields

First-order auto regressive, 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 + \delta_{\mathbf{s},t} + \ldots \right),\\ \boldsymbol{\delta}_{t=1} &\sim \mathrm{MVNormal} (\boldsymbol{0}, \boldsymbol{\Sigma}_{\epsilon}),\\ \boldsymbol{\delta}_{t>1} &= \rho \boldsymbol{\delta}_{t-1} + \sqrt{1 - \rho^2} \boldsymbol{\epsilon}_{t}, \: \boldsymbol{\epsilon}_{t} \sim \mathrm{MVNormal} \left(\boldsymbol{0}, \boldsymbol{\Sigma}_{\epsilon} \right), \end{aligned} \] where \(\rho\) is the correlation between subsequent spatiotemporal random fields. The \(\rho \boldsymbol{\delta}_{t-1} + \sqrt{1 - \rho^2}\) term scales the spatiotemporal variance by the correlation such that it represents the steady-state marginal variance. 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 + \delta_{\mathbf{s},t} + \ldots \right),\\ \boldsymbol{\delta}_{t=1} &\sim \mathrm{MVNormal} (\boldsymbol{0}, \boldsymbol{\Sigma}_{\epsilon}),\\ \boldsymbol{\delta}_{t>1} &= \boldsymbol{\delta}_{t-1} + \boldsymbol{\epsilon}_{t}, \: \boldsymbol{\epsilon}_{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. Note that, in contrast to the AR(1) parametrization, the variance is no longer the steady-state marginal variance.

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 (probability ~ ..., weights = N).

Code defined within TMB.

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

Beta-binomial

\[ \operatorname{BetaBinomial} \left(N, \mu, \phi \right) \] where \(N\) is the number of trials, \(\mu\) is the mean success probability, and \(\phi\) is a precision parameter that controls overdispersion relative to a Binomial. The implied Beta parameters are \(\alpha = \mu \phi\) and \(\beta = (1 - \mu)\phi\), and the variance is $ [y] = N (1 - ), , $ which exceeds the Binomial variance unless \(\phi \to \infty\). Available links are logit and cloglog. This family is useful for overdispersed counts of successes/failures (e.g., aggregated Bernoulli data, proportions with extra-binomial variation).

Code defined within sdmTMB and implemented in the TMB likelihood at src/sdmTMB.cpp (e.g., lines 974–982).

Example: family = betabinomial(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")

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")

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")

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 mixture

This is a 2 component mixture that extends the NB2 distribution, following mixture-distribution approaches for shoaling/aggregation in survey data (Thorson et al. 2011).

\[ (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")

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^2 / \phi\).

Code defined within TMB.

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

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")

Generalized gamma

\[ \operatorname{GenGamma} \left( \mu, \phi, Q^{\mathrm{gg}}\right) \]

sdmTMB implements the Prentice (1974) parameterization introduced for spatiotemporal models and index standardization by Dunic et al. (2025), with parameters mean \(\mu\), scale \(\phi\) (standard deviation on the log scale), and a shape parameter \(Q^{\mathrm{gg}}\) (reported as Generalized gamma Q).

Here, \(\mu\) is the mean on the data scale, \(\phi\) is the log-scale standard deviation, and \(Q^{\mathrm{gg}}\) controls the shape (\(Q^{\mathrm{gg}}\to 0\) yields the lognormal; \(Q^{\mathrm{gg}}= \phi\) yields the gamma). See Dunic et al. (2025) for the full PDF in the Prentice formulation as implemented. Links available: identity, log, inverse. This flexibility is useful for right-skewed positive responses with tails heavier or lighter than gamma/lognormal, and is often paired in a hurdle/mixture as delta_gengamma() for zero-inflated biomass or catch data (see Dunic et al. (2025)).

Code defined within sdmTMB.

Example: family = gengamma(link = "log"); delta/hurdle: family = delta_gengamma(link1 = "logit", link2 = "log").

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, motivated by mixture-distribution treatments of aggregation in survey data (Thorson et al. 2011),

\[ (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}\). The variance follows the usual mixture formula: \((1-p)\left(\mu_1^2 / \phi + \mu_1^2\right) + p\left(\mu_2^2 / \phi + \mu_2^2\right) - \left[(1-p)\mu_1 + p\mu_2\right]^2\).

Here, and for the other mixture distributions, the probability of the larger mean can be obtained from plogis(fit$model$par[["logit_p_extreme"]]) 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 fix this value and pass this to the family:

sdmTMB(...,
  family = gamma_mix(link = "log", p_extreme = 0.01)
)

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, again in the spirit of mixture approaches for aggregating/ shoaling data (Thorson et al. 2011),

\[ (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}\). The log-scale variance of the mixture is not simply \(\phi^2\); it can be obtained with the standard mixture-variance formula using component log-means \(\log \mu_i - \phi^2/2\) and log-variance \(\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. Like with the gamma mixture, fixed probabilities of extreme events (\(p\) in notation above) can be passed in, e.g.

sdmTMB(...,
  family = delta_lognormal_mix(p_extreme = 0.01)
)

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.

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)

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). 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 for 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.