tinyVAST (J. T. Thorson et al. 2025) is a bivariate extension of a generalized linear mixed model (see Tables 1 and 2 for notation), which includes two linear predictors that each include four additive components:
Spatial interactions among variables: The user can specify interactions among variables at a given site (or for spatially correlated latent variables) using arrow notation derived from path analysis, based on the interface from R package sem;
Temporal interaction among variables: The user can specify simultaneous and lagged interactions among variables over time using an expanded arrow-and-lag notation that is derived from R package dsem, where these interactions the annual intercept for a given variable therefore apply uniformly for all locations.
Spatio-temporal interactions among variables: The user can specify simultaneous and lagged interactions among variables over time, where these interactions occur on a site-by-site basis.
Generalized additive model: The user specifies a formula for a generalized additive model (GAM) that is interpreted by R package mgcv. If other model components are missing, tinyVAST estimates parameters that are similar to mgcv, with small differences resulting from different methods for parameter estimation;
These four components are assembled into two linear predictors:
\[ \begin{aligned} p_{\mathrm 1,i} &= \underbrace{\mathbf X_{\mathrm 1,i} \mathbf\alpha_{\mathrm 1} + \mathbf Z_{\mathrm 1,i} \mathbf\gamma_{\mathrm 1}}_\text{formula} + \underbrace{\mathbf A_{i} \mathbf\Omega_{\mathrm 1,c[i]}}_\text{space_term} + \underbrace{\mathbf D_{\mathrm 1,c[i],t[i]}}_\text{time_term} + \underbrace{\mathbf A_i \mathbf E_{\mathrm 1,c[i],t[i]}}_\text{spacetime_term} + \underbrace{\mathbf W_{\mathrm 1,i} (\mathbf A_{i} \mathbf\Xi_{\mathrm 1})^{\mathrm T} }_\text{spatially_varying} \\ p_{\mathrm 2,i} &= \underbrace{\mathbf X_{\mathrm 2,i} \mathbf\alpha_{\mathrm 2} + \mathbf Z_{\mathrm 2,i} \mathbf\gamma_{\mathrm 2}}_\text{formula} + \underbrace{\mathbf A_{i} \mathbf\Omega_{\mathrm 2,c[i]}}_\text{space_term} + \underbrace{\mathbf D_{\mathrm 2,c[i],t[i]}}_\text{time_term} + \underbrace{\mathbf A_i \mathbf E_{\mathrm 2,c[i],t[i]}}_\text{spacetime_term} + \underbrace{\mathbf W_{\mathrm 2,i} (\mathbf A_{i} \mathbf\Xi_{\mathrm 2})^{\mathrm T} }_\text{spatially_varying} \end{aligned} \]
where
mgcv;and terms are defined similarly for the second linear predictor \(p_{\mathrm 2,i}\) except using subscript-2. The linear predictors are then passed through a bivariate inverse-link function to specify the distribution for errors:
\[ y_i \sim f_{e[i]}( g_{e[i]}^{-1} (p_{\mathrm 1,i}, p_{\mathrm 2,i}), \theta_{e[i]} ) \]
where
In the simple case, the distribution only requires a single linear predictor such that \(\mathbf p_{\mathrm 2} = \mathbf 0\) by construction and drops out of the model. In this case, the model collapses to a generalized linear mixed model. For example we might have a log-link and Poisson distribution, such that it collapses to a log-linked Poisson GLMM, \(y_i \sim \mathrm{Poisson}(e^{p_{\mathrm 1,i}})\).
However, tinyVAST can also handle a delta-model using either logit-log or Poisson-linked link functions:
\[ g_{e[i]}^{-1} (p_{\mathrm 1,i}, p_{\mathrm 2,i}) = ( \mu_{\mathrm 1,i}, \mu_{\mathrm 2,i} ) \]
where \(\mu_{\mathrm 1,i}\) and \(\mu_{\mathrm 2,i}\) are the two linear predictors, such that \(\mathbb{E}[y_i] = \mu_{\mathrm 1,i} \mu_{\mathrm 2,i}\);
For the conventional logit-log bivariate link function we use a logit-link for encounter probabilities and a log-link for positive catch rates:
\[ \begin{aligned} \mu_{\mathrm 1,i} &= \frac{e^{p_{\mathrm 1,i}}}{1+e^{p_{\mathrm 1,i}}} \\ \mu_{\mathrm 2,i} &= e^{p_{\mathrm 2,i}} \end{aligned} \]
while for the Poisson-linked link function (J. T. Thorson 2018) we specify a complemetary log-log link for encounter probabilities and define the second link function such that \(\mu_{\mathrm 1,i} \mu_{\mathrm 2,i} = e^{p_{\mathrm 1,i}} e^{p_{\mathrm 2,i}}\):
\[ \begin{aligned} \mu_{\mathrm 1,i} &= 1 - e^{-e^{p_{\mathrm 1,i}}} \\ \mu_{\mathrm 2,i} &= \frac{e^{p_{\mathrm 1,i}}}{\mu_{\mathrm 1,i}} e^{p_{\mathrm 2,i}} \end{aligned} \] where \(e^{p_{\mathrm 1,i}}\) is the density for an underlying point process, and \(e^{p_{\mathrm 2,i}}\) is biomass per point (i.e., animal or group).
In either the conventional or Poisson-linked delta model, \(\mu_{\mathrm 1,i} = \mathrm{Pr}(Y>0)\) is the encounter probability (and is fitted with a Bernoulli distribution), while \(\mu_{\mathrm 2,i}\) is the central tendancy for positive values, i.e., \(\mu_{\mathrm 2,i} = \mathbb{E}(Y | Y>0)\).
Usefully, the interface allows analysts to specify a different distribution \(e_i\) for different partitions of the data. This allows the model to jointly analyze different data types (Grüss and Thorson 2019), or data with different magnitudes of sampling error.
Linear predictors include spatially autocorrelated latent variables. These variables are treated as Gaussian Markov random fields (GMRFs), and evaluating the probability density of GMRFs involves calculating the precision matrix \(\mathbf Q_{\mathrm{domain}} = \mathbf\Sigma^{-1}\) as the inverse of the spatial covariance matrix (J. Thorson and Kristensen 2024). tinyVAST involves three options for specifying this spatial precision:
The analyst can approximate spatial variation over a continuous
two-dimensional surface by constructing finite element mesh (FEM),
treating the value at vertices as a GMRF, and then using bilinear
interpolation (i.e., a piecewise linear approximation) to interpolate
from vertices to the spatial domain. This method was developed by (F. Lindgren, Rue, and Lindström 2011), as
popularized in the software R-INLA (Lindgren
2012), first implemented in TMB by (J. T.
Thorson et al. 2014), and using elements constructed by R-package
fmesher (F. Lindgren 2023).
In this case, the precision is constructed as:
\[ \mathbf Q_{\mathrm{domain}} = \tau^2 ( \kappa^4 \mathbf M_{\mathrm 0} + 2\kappa^2 \mathbf M_{\mathrm 1} + \mathbf M_{\mathrm 2} ) \] where every row \(\mathbf A_i\) of the interpolation matrix \(\mathbf A\) is nonzero for only the three vertices of the triangle that contains sample \(i\)
Alternatively, the analyst can apply a simultaneous autoregressive process (Ver Hoef, Hanks, and Hooten 2018), specifying an areal model that represents the value within spatial strata:
\[ \mathbf Q_{\mathrm{domain}} = \tau^2 (\mathbf I - \kappa \mathbf A^*)^2 \] where \(\mathbf A^*\) is the adjacency matrix of the graph specified by the analyst, \(\kappa\) is the estimated partial correlation between adjacent areas, and each row \(\mathbf A_i\) of interpolation matrix \(\mathbf A\) is nonzero only for the single spatial stratum containing sample \(i\) (and noting that adjacency matrix \(\mathbf A^*\) is different from interpolation matrix \(\mathbf A\), but we use the same notation for both due to a collision in notational standards).
Finally, the analyst can specify that sites are partially correlated if they are adjacent along a stream network (while ignoring flow direction). This results in an Onstein-Uhlenbeck process along the stream network (Charsley et al. 2023), or an exponential tail-down model:
\[ \mathbf{Q}_{\mathrm{domain}} = \mathbf{(I+D-P)}^T \mathbf{I+D}^{-1} \mathbf{(I+D-P)} \] where \(\mathbf D\) is a sparse diagonal matrix with diagonal elements
\[ d_{i,j} =\mathrm{exp}(-2 \theta |\mathbf s_i, \mathbf s_j|) / (1 - \mathrm{exp}(-2 \theta |\mathbf s_i, \mathbf s_j|)) \] where \(|\mathbf s_i, \mathbf s_j|\) is the distance from downstream (parent) node \(j\) to upstream (child) node \(i\), and \(\theta\) is the O-U parameter governing the decorrelation rate. Similarly, \(\mathbf P\) is a sparse matrix containing values \(\rho_{i,j}\), where: \[ \rho_{i,j} = \mathrm{exp}(-\theta |\mathbf s_i, \mathbf s_j|) / (1 - \mathrm{exp}(-2 \theta |\mathbf s_i, \mathbf s_j|)) \] such that \(\rho_{i,j}\) the regression slope predicting upstream node \(i\) from downstream node \(j\). The spatial interpolation matrix \(\mathbf A\) again has row \(\mathbf A_i\) for each sampling or predictive location, and \(\mathbf A_i\) is nonzero for the two nodes immediately up and downstream from a given location, with values such that a given location is predicted as a weighted average based upon the distance from those two nodes.
For the SPDE and SAR spatial domains, the term \(\tau\) is defined such that precision \(\mathbf Q_{\mathrm{domain}}\) has unit variance. This is done because the spatial domain always arises in combination with other parameters, which are used to define the variance of the associated spatial variable.
tinyVAST also involves specifying a structural equation model (SEM). This SEM be viewed either:
To specify a SEM, the user uses arrow notation derived from
package sem (Fox 2006), as
described in TMB by (J. Thorson and Kristensen
2024). For example, to specify a linear model this involves:
This then estimates a single slope parameter (represented with a one-headed arrow), as well as the variance of \(W_1\) and \(W_2\) (specified with two-headed arrows). In a more complicated case, \(W_1\) might cause \(W_2\), which in turn causes \(W_3\). This is then represented as:
# Path, parameter_name, start value
w1 -> w2, b_12, 0
w2 -> w3, b_23, 0
w1 <-> w1, s_1, 1
w2 <-> w2, s_2, 1
w3 <-> w3, s_3, 1SEM interactions can be as complicated or simple as desired, and can include:
parameter_name is provided as NA and the
starting value that follows is the fixed value;parameter_name is provided for multiple rows of the text
file.In this preceding example, path coefficients for one-headed arrows then define path matrix \(\mathbf P\):
\[ \mathbf P = \begin{pmatrix} 0 & 0 & 0 \\ b_{12} & 0 & 0 \\ 0 & b_{23} & 0 \end{pmatrix} \]
and coefficents for two-headed arrows define the Cholesky \(\mathbf G\) of the exnogenous covariance matrix \(\mathbf G^T \mathbf G\):
\[ \mathbf G = \begin{pmatrix} s_{1} & 0 & 0 \\ 0 & s_{2} & 0 \\ 0 & 0 & s_{3} \end{pmatrix} \]
These matrices are define a simultaneous equation model:
\[ \mathbf{ w = P w + \epsilon} \\ \mathbf\epsilon \sim \mathrm{MVN}( \mathbf 0, \mathbf G^T \mathbf G ) \] where the variance \(\mathrm{Var}(\mathbf w) = (\mathbf{I - P})^{-1} \mathbf G^2 (\mathbf{I - P}^T)^{-1}\). This then results in a sparse precision matrix:
\[ \mathbf Q = (\mathbf{I - P}^T) \mathbf G^{-1} \mathbf G^{-T} (\mathbf{I - P}) \]
where this precision matrix \(\mathbf Q\) is then used as a modular component in the larger tinyVAST model.
Similarly, tinyVAST involves specifying dynamic structural equation models (DSEM) (J. T. Thorson et al. 2024). To specify a DSEM, the user uses arrow-and-lag notation. For example, to specify a univariate first-order autoregressive process:
If there were four time-intervals (\(T=4\)) this would then result in the path matrix:
\[ \mathbf P = \begin{pmatrix} 0 & 0 & 0 & 0 \\ \rho & 0 & 0 & 0 \\ 0 & \rho & 0 & 0 \\ 0 & 0 & \rho & 0 \end{pmatrix} \]
and when the DSEM involves multiple times and variables, the sparse precision is formed by summing across the Kronecker product of time-lag and interaction matrices. This DSEM defines a GMRF over a nonseparable interaction of time and variables, represented by a matrix \(\mathbf Q_{\mathrm{time\_term}}\) with dimension \(CT \times CT\). The user can specify a separate arrow-and-lag notation to define the precision matrix for the time-variable interaction \(\mathbf Q_{\mathrm{time\_term}}\) and for the space-time-variable interaction \(\mathbf Q_{\mathrm{spacetime\_term}}\).
The precision matrix \(\mathbf Q_{\mathrm{time\_term}}\) for the time term \(\mathbf D\) is used to define the time-varying intercept for each variable:
\[ \mathrm{vec}(\mathbf D) \sim \mathrm{MVN}(\mathbf 0, \mathbf Q_{\mathrm{time\_term}}) \] Meanwhile, the space-time term \(\mathbf Q_{\mathrm{spacetime\_term}}\) is combined with the spatial precision \(\mathbf Q_{\mathrm{space\_term}}\) as we explain next.
tinyVAST uses the SEM and DSEM notation to construct the joint precision for the space-variable interaction \(\mathbf\Omega\) with dimension \(S \times C\), and the space-time-variable interaction \(\mathbf E\) with dimension \(S \times C \times T\). To do so, it constructs the separable precision for each process:
\[ \mathrm{vec}(\mathbf E) \sim \mathrm{MVN}(\mathbf 0, \mathbf Q_{\mathrm{spacetime\_term}} \otimes \mathbf Q_{\mathrm{domain}}) \]
where the precision matrix \(\mathbf Q_{\mathrm{spacetime\_term}} \otimes \mathbf Q_{\mathrm{domain}}\) has dimension \(STC \times STC\) to match the length of \(\mathrm{vec}(\mathbf E)\), and \[ \mathrm{vec}(\mathbf\Omega) \sim \mathrm{MVN}(\mathbf 0, \mathbf Q_{\mathrm{space\_term}} \otimes \mathbf Q_{\mathrm{domain}}) \]
where the precision matrix \(\mathbf Q_{\mathrm{space\_term}} \otimes \mathbf Q_{\mathrm{domain}}\) has dimension \(SC \times SC\), to match the length of \(\mathrm{vec}(\mathbf\Omega)\). This specification generalizes spatial factor analysis (Thorson et al. 2015) and spatial dynamic factor analysis (J. T. Thorson et al. 2016).
The analyst can specify a generalized additive model using syntax
from package mgcv (Wood
2006), and parsed into TMB using an interface developed by
sdmTMB (Anderson et al.
n.d.). For example this might involve:
If year and species are factors and
depth and log_area are continuous, then this
would specify a fixed effect for each level year, a spline
smoother for depth, using log_area as an
offset, and estimating a random intercept for each level of
species. This formula is parsed internally to assemble
fixed effects in design matrix \(\mathbf
X\) and the basis functions for spline smoothers and random
effects in design matrix \(\mathbf Z\).
The coefficients \(\mathbf\gamma\)
associated with smoothers and random effects are then specified to
follow a GMRF:
\[ \mathbf\gamma \sim \mathrm{GMRF}( \mathbf 0, \mathbf Q_{\mathrm{gam}}) \] where \(\mathbf Q_{\mathrm{gam}}\) is a blockwise diagonal matrix, assembled from estimated variance parameters and matrices constructed by mgcv.
Finally, the analyst can account for spatial variation in the relationship between a specified covariate and the estimated response. This is done by specifying a spatially varying coefficient (SVC) (J. T. Thorson et al. 2023), and this can be used e.g. to account for spatially varying differences in catchability/detectability (Grüss et al. 2023) or to account for nonlocal impacts of regional oceanographic indices (J. T. Thorson 2019). This involves estimating each SVC \(l \in \{1,2,...,L\}\) as a zero-centered Gaussian Markov random field while estimating its corresponding variance \(\sigma_{\mathrm{\xi},l}^2\):
\[ \Xi_{\mathrm 1,l} \sim \mathrm{GMRF}( \mathbf 0, \sigma_{\mathrm{\xi},l}^{-2} \mathbf Q_{\mathrm{domain}} ) \]
Any covariate specified as an SVC it is also typically specified in
the formula to estimate a non-zero mean for the SVC. If the
estimated variance of the SVC approaches zero, it then suggests that the
covariate slope does not vary spatially.
| Symbol | Description |
|---|---|
| \(i\) | Index for each sample, \(i\) in \((1,2,...,I)\) |
| \(s[i]\) | spatial coordinate for sample \(i\), \(s\) in \((1,2,...,S)\) |
| \(t[i]\) | time-interval for sample \(i\), \(t\) in \((1,2,...,T)\) |
| \(c[i]\) | category for sample \(i\), \(c\) in \((1,2,...,C)\) |
| \(e[i]\) | error distribution and link function for sample \(i\) |
| Symbol | Code | Description |
|---|---|---|
| \(y\) | y_i |
Observed response data |
| \(p_1\) | p_i |
first linear predictor |
| \(p_2\) | p2_i |
second linear predictor |