Formula keyword grid

This page is the syntax lookup for choosing a covariance keyword. The formula grammar has two decisions: what relates the grouping levels, and what covariance among traits is required at that grouping level.

Use Choose your next gllvmTMB step when you know the modelling task but need help locating this keyword reference or the next post-fit function.

The keywords below are live API, but a valid formula does not guarantee that a complex covariance is estimable from a particular data set. Fit health, replication, trait count, and the scientific covariance question still matter.

Most readers only need the two choices below; the optional rho setting on a structured term is covered later, in Advanced: estimating the source strength.

Choose in two steps

The first five sources describe relationships among stable units: for example, individuals, species, locations, or rows of a supplied relatedness matrix. Temporal is the additional source for ordered repeated states within a series. It uses the same three trait-covariance modes, but its series and time coordinates have a different scientific meaning.

No prefix: levels are independent, as in an ordinary random effect.
animal_: individuals are related through a pedigree, A, or Ainv.
phylo_: species are related through a tree or phylogenetic covariance.
spatial_: locations share a mesh-based spatial field, represented with a sparse Gaussian Markov random field (GMRF) approximation so nearby locations can be correlated.
kernel_: levels are related through a supplied dense covariance matrix K — the generic dense-matrix source (the fifth grid row).
temporal_: ordered states within each series share an AR1 or OU kernel. time = names the occasion or elapsed-time column; AR1 preserves integer gaps and OU uses elapsed-time distances.

These sources are alternatives based on study design, not a progression from fine to broad scale. Add a source prefix only when that structured relationship is part of the estimand.

2. What trait covariance is needed?

There are three fundamental cross-trait structures, from simplest to richest:

indep: a separate marginal variance for each trait and zero cross-trait covariance — a diagonal model. Its parsimony modifier common = TRUE ties all trait variances to one shared value (the scalar special case).
dep: a full unstructured trait covariance with (T(T+1)/2) covariance parameters. Use only when trait count and grouping-level replication support a saturated covariance.
latent: a parsimonious reduced-rank cross-trait covariance $\boldsymbol\Lambda\boldsymbol\Lambda^{\mathsf T}$ represented by (d) latent variables. Its modifier unique = TRUE adds the trait-diagonal $\Psi$ companion.

scalar and unique are modifiers, not separate modes. scalar is indep with the variances tied to one shared value; unique = TRUE adds the diagonal- $\Psi$ companion to a latent term. The standalone deprecated unique() keyword remains compatibility syntax for the plain diagonal case. The named scalar(), phylo_scalar(), animal_scalar(), spatial_scalar(), and kernel_scalar() keywords fit exactly this model but are soft-deprecated compatibility syntax (they emit a one-time warning); write indep(..., common = TRUE) in new code.

All fitted keyword covariances are random-effect covariances on the model’s linear-predictor or latent scale. They are not automatically response-scale covariances for non-Gaussian traits.

The 5 × 3 grid, plus temporal

The five stable-unit sources form the familiar 5 × 3 grid. Temporal adds a sixth, within-series row with the same three trait-covariance modes. Thus the formal grammar is a 6 × 3 source-by-mode grid, while this presentation keeps the important distinction visible. Every cell is a live keyword. The scalar and unique sub-cases are modifiers, listed under the table.

Stable-unit sources and the additional temporal source, each paired with a trait-covariance mode
Source	Independent	Dependent	Latent
none	`indep()`	`dep()`	`latent()`
animal	`animal_indep()`	`animal_dep()`	`animal_latent()`
phylogenetic	`phylo_indep()`	`phylo_dep()`	`phylo_latent()`
spatial	`spatial_indep()`	`spatial_dep()`	`spatial_latent()`
kernel	`kernel_indep()`	`kernel_dep()`	`kernel_latent()`
temporal (within series)	`temporal_indep()`	`temporal_dep()`	`temporal_latent()`

Scalar (one shared variance). Add common = TRUE to any indep term — e.g. indep(0 + trait | g, common = TRUE). The soft-deprecated scalar() / *_scalar() keywords fit the same model but warn once and are kept only for compatibility. Omitting a multi-trait covariance term entirely means no such term is fitted.

Unique (diagonal- $\Psi$ companion). Add unique = TRUE to a latent term (see the latent default below). The standalone unique() keyword is soft-deprecated in favour of indep() (for a plain diagonal) and latent(unique = TRUE) (for the paired companion).

The temporal addition

temporal_indep(0 + trait | series, time = occasion), temporal_dep(...), and temporal_latent(..., d = 1) form the temporal row. Here series identifies one independent trajectory, and time = occasion uses the bare occasion column from the data: it is neither a formula term nor a character string. A replicate = measurement column identifies repeated measurements at the same (series, occasion) state. The provider creates that state index privately; it does not change unit or unit_obs, which retain their ordinary stable-unit and within-unit meanings.

AR1 accepts ordered integer occasions and preserves unobserved occasions as gaps; OU accepts elapsed numeric time without automatic rescaling. With temporal_latent(unique = TRUE), Psi is itself correlated through the temporal kernel, not IID occasion noise. The temporal article shows matching long and traits(...) wide calls. In scope: focused Gaussian local fits, source parsing, covariance oracles, simulation, and ordinary unit/unit-observation composition. Partial: recovery is not yet a general claim. The bounded temporal-only forecast, profile, bootstrap, and supplied-candidate comparison helpers have their own contracts; generic forecasts, new-data prediction, generic intervals/profiles, and automatic selection remain unavailable. Four replicated-AR1 additive temporal_indep() pairs with one diagonal kernel, phylogenetic, animal, or spatial source are separately implemented, with pair-specific local evidence only: the kernel combination met its prespecified recovery criteria in one local simulated example, while the phylogenetic, animal, and spatial recovery studies retain failed cases. These pairs require complete Gaussian identity-link data, native TMB/Laplace ML, at least two complete measurements per series–occasion, and a source fixed at rho = 1. Supply the source inside its keyword (K =, tree =/vcv =, pedigree =/A =/Ainv =, or mesh =). The spatial route needs time and space contrasts that are not proportional. OU source combinations, temporal dep/latent combinations, temporal slopes, other families, higher rank, source-by-time interactions, and other structured-source combinations are unavailable. See Temporal covariance for the standalone workflow and the adjacent source-combination limits.

What the modes estimate

For the intercept-only grid, let $K_{gg'}$ be an entry of the relationship or model-implied covariance operator among grouping levels. For the ordinary no-prefix row, $\mathbf K = \mathbf I$ . Animal, phylogenetic, and dense-kernel inputs supply fixed covariance operators and may have non-unit diagonals. Spatial terms instead imply $\mathbf K(\boldsymbol\theta)$ through a mesh-based precision whose range and scale parameters are estimated.

The clearest ordering-free definition is

$\operatorname{Cov}(b_{gt}, b_{g't'}) = \Sigma_{T,tt'} K_{gg'},$

where $\boldsymbol\Sigma_T$ is the trait covariance implied by the selected mode: $\mathrm{diag}(\psi_1,\ldots,\psi_T)$ for indep (collapsing to the shared $\sigma^2\mathbf I_T$ under the common = TRUE / scalar modifier), a free $\boldsymbol\Sigma_T$ for dep, and $\boldsymbol\Lambda\boldsymbol\Lambda^{\mathsf T}$ with an optional $\boldsymbol\Psi$ companion for latent. This elementwise form does not depend on whether software stacks traits or grouping levels first.

The diagonal in animal_indep() or phylo_indep() is diagonal across traits, but observations remain correlated within each trait through the pedigree or phylogeny. It is not independent residual noise.

The important latent default

Here, unique = is an argument of a latent term that toggles its diagonal Psi companion; it is not a separate covariance mode or a call to the deprecated covariance-function syntax.

The diagonal Psi companion is controlled differently across sources:

ordinary latent() includes $\boldsymbol\Lambda\boldsymbol\Lambda^{\mathsf T}+\boldsymbol\Psi$ by default; use unique = FALSE for loadings only;
animal_latent(), phylo_latent(), spatial_latent(), and kernel_latent() are loadings-only by default;
use unique = TRUE on a source-specific or single-kernel latent term when that same source should also carry a trait-diagonal Psi component.

For animal, phylogenetic, and kernel models, source-specific Psi is diagonal across traits but still correlated across individuals, species, or other levels through the supplied relationship or covariance operator. Spatial Psi is also diagonal across traits, while each trait-specific field remains correlated across locations through the mesh-based spatial precision.

These diagonal statements describe intercept-only trait covariance. With an augmented intercept-and-slope term, *_indep(1 + x | group) estimates a per-trait intercept–slope correlation (a double bar, 1 + x || group, drops it). The random-slope section below gives the corresponding blocks.

Ordinary long and wide syntax

The long and wide forms below express the same covariance choices. The Morphometrics article executes and compares both data shapes.

# Decomposed covariance: Lambda Lambda^T + Psi
fit_long <- gllvmTMB(
  value ~ 0 + trait +
    latent(0 + trait | individual, d = 2),
  data = df_long,
  trait = "trait",
  unit = "individual"
)

fit_wide <- gllvmTMB(
  traits(length, mass, wing, tarsus, bill) ~ 1 +
    latent(1 | individual, d = 2),
  data = df_wide,
  unit = "individual"
)

Replace latent(...) with indep(...) for a diagonal baseline or dep(...) for a full covariance. In wide formulas, 1 inside a bar-style term expands to the long-form 0 + trait representation.

Exact source-specific syntax

The latent animal/phylogenetic helpers take the grouping column directly. The indep and dep counterparts use a bar formula, and the one-shared-variance (scalar) case is the common = TRUE modifier on indep. Spatial helpers all use the bar-form orientation.

Animal

animal_indep(0 + trait | id, pedigree = ped, common = TRUE)  # one shared variance
animal_indep(0 + trait | id, pedigree = ped)
animal_dep(0 + trait | id, pedigree = ped)
animal_latent(id, d = 2, pedigree = ped)                 # loadings only
animal_latent(id, d = 2, pedigree = ped, unique = TRUE)  # + animal Psi

Phylogenetic

phylo_indep(0 + trait | species, tree = tree, common = TRUE)  # one shared variance
phylo_indep(0 + trait | species, tree = tree)
phylo_dep(0 + trait | species, tree = tree)
phylo_latent(species, d = 2, tree = tree)                 # loadings only
phylo_latent(species, d = 2, tree = tree, unique = TRUE)  # + phylo Psi

Use tree = for the usual phylogenetic workflow. vcv =, A =, and Ainv = support advanced precomputed covariance/precision inputs; their row and column names must match the grouping levels.

Spatial

spatial_indep(0 + trait | site, mesh = mesh, common = TRUE)  # one shared variance
spatial_indep(0 + trait | site, mesh = mesh)
spatial_dep(0 + trait | site, mesh = mesh)
spatial_latent(0 + trait | site, d = 2, mesh = mesh)
spatial_latent(0 + trait | site, d = 2, mesh = mesh, unique = TRUE)

Spatial models require a mesh built from the model’s coordinate columns with make_mesh(). Pass that mesh either inside the spatial keyword or through the top-level mesh = argument to gllvmTMB(). A source-valid formula still needs enough locations and spatial spread to identify range and covariance parameters.

Dense kernels (the fifth grid row)

A supplied dense covariance matrix K is the generic dense-matrix source. Its keywords follow the same grid as the other sources, with the tied-variance case as the common = TRUE modifier on kernel_indep().

kernel_indep(unit_id, K = K, name = "known", common = TRUE)  # one shared variance x K
kernel_indep(unit_id, K = K, name = "known")
kernel_dep(unit_id, K = K, name = "known")
kernel_latent(unit_id, K = K, d = 2, name = "known")
kernel_latent(
  unit_id, K = K, d = 2, name = "known", unique = TRUE
)

For a single named kernel, unique = TRUE adds that kernel’s diagonal Psi companion. Multiple named-kernel fits currently use the latent-only route; do not assume that each component receives its own Psi. When kernels overlap strongly, the fitted sum may be more stable than attribution to individual kernel components.

Named kernels should encode a relationship fixed before fitting. The model does not turn an analyst-supplied association matrix into evidence of a causal or reciprocal process. When more than one named kernel is fitted, use diagnose_kernel_separability() before interpreting component-specific output; collapse indistinguishable kernels rather than assigning them separate biological stories.

What sits outside the grid

Slope variation by response column

The 6 × 3 grid describes trait covariance for random effects. The slope helpers answer a separate question: how do fixed-predictor coefficients vary across the response columns? They are not a new grid row, column, or mode.

Column relationship	Use
no supplied relationship	`slope(x \| trait)`
a phylogenetic tree	`phylo_slope(x \| trait, tree = tree)`
a pedigree or relationship matrix	`animal_slope(x \| trait, A = A)`
a labelled covariance matrix	`kernel_slope(x \| trait, K = K)`
one coordinate pair per column	`spatial_slope(x \| trait, mesh = column_mesh)`

For this helper family, the right side must be the variable that labels the response columns, supplied through trait = and usually named trait. The helper adds no random intercept. With two or more predictors, | estimates their full covariance and || keeps their predictor variances but fixes their covariance to zero. Retrieve that predictor-by-predictor matrix with extract_Sigma(fit, level = "column_slope").

value ~ 0 + trait + moisture + canopy +
  phylo_slope(moisture + canopy | trait, tree = tree)

This is Gaussian long-format point estimation only. Wide column-slope grammar, non-Gaussian responses, latent predictor covariance, several helpers in one model, and intervals are deliberately deferred. Group-indexed phylo_slope(x | species) and animal_slope(x | id) remain current and warning-free for models whose coefficients vary across sampled units. They answer a different grouping-axis question from the response-column slope helpers described here. In short: use phylo_slope(x | trait, tree = tree) when the slope of x should vary across response columns (e.g. across traits or species treated as responses) and use phylo_slope(x | species) when the slope of x should vary across the sampled grouping units (e.g. across individuals or sites) that the tree relates.

Intercepts and slopes by response column

column_coef(), phylo_coef(), animal_coef(), kernel_coef(), and spatial_coef() extend the coefficient basis to an explicit random intercept, a slope, or both. They are also outside the 6 × 3 grid: their covariance is across coefficient-basis columns, not a new trait-covariance mode.

# IID response-column intercepts and slopes
value ~ 1 + column_coef(1 + moisture | trait)

# Tree-correlated response-column intercepts and slopes; estimate rho
value ~ 1 + phylo_coef(1 + moisture | trait, tree = tree, rho = NULL)

# Pedigree-correlated response-column intercepts and slopes; fix rho
value ~ 1 + animal_coef(1 + moisture | trait, pedigree = pedigree, rho = 1)

# C3/C4 grand intercepts and slopes, plus kernel-correlated
# response-column intercept and slope deviations; estimate rho
value ~ 0 + pathway + moisture:pathway +
  kernel_coef(1 + moisture | trait, K = K, name = "environment", rho = NULL)

# The same hierarchical model from a wide response data frame;
# pass the keyed pathway table through column_data
traits(sp1, sp2, sp3, sp4) ~ 0 + pathway + moisture:pathway +
  kernel_coef(1 + moisture | trait, K = K, name = "environment", rho = NULL)

# Species locations as the response-column relationship; rho is fixed at 1
value ~ 0 + pathway + moisture:pathway +
  spatial_coef(1 + moisture | trait, mesh = column_mesh)

Use | for a full coefficient covariance and || for a diagonal one. Retrieve the ordered basis covariance and source metadata with extract_Sigma(fit, level = "column_coef"). All five helpers support Gaussian point models through long or traits(...) wide input. Numeric rho fixes the phylogenetic, animal, or dense-kernel mixture; phylo_coef(rho = NULL) and kernel_coef(rho = NULL) estimate one interior value. animal_coef() accepts exactly one pedigree, relationship covariance A, or precision Ainv; kernel_coef() accepts one labelled dense positive-definite covariance K and reports its name and supplied scale. spatial_coef() accepts one mesh built with exactly one labelled coordinate pair per response column and currently requires rho = 1; its extractor reports the aligned coordinates, fitted range, and projected unit-diagonal source covariance. Estimated animal rho, spatial IID mixtures, non-Gaussian responses, and intervals remain unavailable. The existing slope-only helpers above remain current, warning-free, and non-deprecated.

For exact compatibility, a no-intercept dense-vcv phylo_coef(..., rho = 1) call uses the released phylo_slope() conditioning K + 1e-8 I. Other fixed/estimated coefficient routes use the raw covariance-scale mixture. At the corresponding no-intercept rho = 1 endpoint, animal_coef() uses the released animal_slope() route exactly. Its dense-A route therefore inherits the released A + 1e-8 I conditioning seam; pedigree and sparse-Ainv endpoints retain their released precision. At the no-intercept rho = 1 kernel endpoint, kernel_coef() uses kernel_slope() exactly and therefore retains raw K without a ridge. At the no-intercept spatial endpoint, spatial_coef() uses spatial_slope() exactly and retains its projected-SPDE normalization and range parameter.

For a worked comparison of a tree among sampled species and a tree among response columns, read Where does the phylogeny belong?.

Known sampling covariance

meta_V() is not a grid row. It adds a known observation/sampling covariance to a meta-analytic model rather than estimating a random-effect trait covariance. The exact-additive route requires the same row-aligned matrix in the formula marker and the top-level argument:

value ~ 0 + trait + latent(0 + trait | study, d = 1) + meta_V(V = V)
# In gllvmTMB(...), also pass known_V = V.

Use A/Ainv for relatedness and V for known sampling covariance; they are different scientific objects. This syntax lookup does not provide the applied estimand, diagnostic, and interpretation workflow needed for a full meta-analytic example. meta_known_V() remains a soft-deprecated alias; new code uses meta_V().

What not to combine

At one grouping level and correlation source, specify one trait-covariance story:

do not combine dep with indep or latent at that source;
do not combine standalone indep with latent at that source;
the same rules apply to ordinary, animal, phylogenetic, spatial, and single-kernel terms.

Different grouping levels can carry different covariance terms—for example, a between-individual latent() term and a within-individual indep() term—when the design identifies both.

Advanced: estimating the source strength

Canonical phylo_*, animal_* and kernel_* indep/dep/latent terms also accept trailing rho = 1. Fixed values between zero and one attenuate between-level source covariance while preserving the resolved source diagonal. This does not add a grid row or mode. With unique = TRUE, the same strength applies to both the loadings and Psi; ordinary intercept effects stay separate. Omitted rho and explicit one preserve existing models, including multi-kernel endpoints. New attenuation currently supports one structured trait-intercept block per fit.

rho = NULL estimates one strength only in the documented complete, replicated Gaussian design with retained observation noise and no competing ordinary covariance. Estimated latent terms require rank one and at least four traits. Spatial terms use the same attenuation weight while range stays separately estimated. Spatial recovery remains unvalidated. Augmented slopes and source-strength intervals are not included. A converged optimizer does not establish identification or recovery. See Strength of a known source for runnable long and wide examples, extraction and interpretation. These defaults do not change the separate response-column coefficient helpers.

If the first choice is too complex

Full dep() is unstable or imprecise: try a lower-rank latent() model, or indep() if cross-trait covariance is not essential.
An added LV is weak or unstable across refits: reduce d and compare fit health and information criteria.
A diagonal variance is at zero: treat it as a boundary estimate; compare the simpler model rather than interpreting a tiny component literally.
A relationship input fails: check grouping-level names, matrix dimensions, positive-semidefiniteness, tree tips, or mesh/data alignment.
You are unsure which source prefix to use: return to the biological sampling design. A more elaborate source is not automatically a better model.

See Fit diagnostics before interpreting the selected covariance and Covariance and correlation for the corresponding estimands.