Fit multivariate response models from wide or long data

Top-level entry point for models with several responses per site, individual, species, or study. Start from a wide data frame with traits() on the formula left-hand side, or from already-stacked long data with one row per (unit, trait) observation. Both routes estimate the same trait covariance, pairwise correlations, shared latent axes, and trait-specific variance. The formula syntax also supports fixed effects plus covariance-structure keywords organised by correlation source (none / animal / phylo / spatial) and mode (scalar / independent / dependent / latent):

Usage

gllvmTMB(
  formula,
  data,
  trait = "trait",
  unit = "site",
  unit_obs = "site_species",
  cluster = "species",
  cluster2 = NULL,
  family = gaussian(),
  weights = NULL,
  REML = FALSE,
  mesh = NULL,
  phylo_vcv = NULL,
  phylo_tree = NULL,
  known_V = NULL,
  lambda_constraint = NULL,
  Xcoef_fixed = NULL,
  control = gllvmTMBcontrol(),
  missing = miss_control(),
  impute = NULL,
  silent = TRUE,
  engine = c("tmb", "julia"),
  ci_method = c("none", "wald", "profile", "bootstrap"),
  ci_level = 0.95,
  ci_nboot = 200L,
  ci_seed = 0L,
  site = NULL,
  species = NULL
)

Arguments

formula

A glmmTMB-style formula, e.g. value ~ 0 + trait + (0 + trait):env_temp + (0 + trait):env_precip. Fixed effects and any of the four-mode grid covstructs above are supported (plus phylo_slope(), animal_slope(), and meta_V()).

data

A data frame. With an ordinary response LHS such as value ~ ..., data is long: one row per (unit, trait) observation, with the trait column named by trait. With a traits() LHS, data is wide: one row per unit and one column per response named inside traits(...). Missing response values are allowed: by default they are dropped before fitting, while missing = miss_control(response = "include") masks them out of the likelihood and keeps original-row accounting for predict_missing(). Ordinary missing predictors, grouping variables, and design-matrix values still error; explicitly modelled missing predictors use mi(x) with missing = miss_control(predictor = "model") and impute = list(...).

trait

Name of the column holding the trait factor (the "trait" dimension of the unit × trait response matrix). Default "trait".

unit

Name of the column holding the between-unit grouping factor (the "unit" dimension of the unit × trait response matrix). Examples: "site" for site × species data, "individual" for behavioural-syndrome data, "species" for PGLLVM, "paper" for systematic mapping. Default "site".

unit_obs

Name of the column holding the within-unit grouping factor — one level per (unit, replicate) cell — used by latent(0 + trait | unit_obs, ...) and indep(0 + trait | unit_obs) for the W-tier covariance. Default "site_species" (the conventional name in joint species distribution modelling; safe for site × species data). For other domains pass e.g. unit_obs = "obs" for behavioural syndromes.

cluster

Name of the column holding the third grouping factor (the "cluster" slot). When the column name matches the row/column names of phylo_vcv (or the tip labels of phylo_tree), this slot also drives the phylogenetic random effects (phylo_latent, phylo_scalar, phylo_indep, phylo_slope). When the column does not match a phylogenetic correlation, the slot still functions as a regular crossed/nested third grouping (e.g. cluster = "population" for 3-level personality data, cluster = "study" for multi-study meta-analysis). Default "species". For datasets with no third grouping the column may be absent from data; the engine will synthesise a placeholder.

The engine does not enforce nesting between unit_obs, unit, and cluster — crossed and nested designs both fit. Two canonical patterns:

Functional biogeography (crossed): unit = "site", unit_obs = "site_species", cluster = "species". Sites and species are crossed; site_species indexes the intersection cell.
3-level personality (strictly nested): unit = "individual", unit_obs = "session_id", cluster = "population".

cluster2

Optional name of the column holding a second independent grouping factor (the "cluster2" slot), for fitting two crossed (or nested) plain diagonal per-trait variance components at once. Default NULL (slot inactive). It is a plain crossed/nested diagonal grouping only – it carries no phylogenetic or spatial correlation (those stay bound to cluster / coords). A indep(0 + trait | <cluster2 col>) term then fits a per-trait variance at this second grouping, exactly as cluster does at the third slot. Example: cluster = "site", cluster2 = "year" to fit a site variance and a year variance simultaneously. As with the other slots, nesting is not enforced (crossed and nested both fit). The cluster2 column must be disjoint from the unit / unit_obs / cluster columns (a diagonal term routes to whichever slot its grouping column matches).

family

A family object. The multivariate engine (formulas with latent(), indep(), etc.) supports gaussian(), binomial() (logit / probit / cloglog), poisson() (log link), ordinal_probit() (the gllvmTMB-native ordinal threshold family with sigma_d = 1 fixed exactly; no delta-method approximation), lognormal() (log link), Gamma(link = "log"), nbinom2() (NB2 negative binomial; log link), nbinom1() (NB1 negative binomial; log link), tweedie() (compound Poisson-Gamma; log link), Beta() (logit link, mean-precision parameterisation; y in (0, 1)), betabinomial() (logit link, mean-precision parameterisation; k-of-n data), student() (heavy-tailed continuous; identity link), truncated_poisson() (zero-truncated count; log link), truncated_nbinom2() (zero-truncated NB2; log link), and the hurdle / two-part families delta_lognormal() / delta_gamma() (Bernoulli for presence + Lognormal/Gamma for the positive component; share one linear predictor under the current implementation – see Details). For multi-trial binomial or beta-binomial (k-of-n) data, write the LHS of formula as cbind(successes, failures) (the canonical R / glm() / glmmTMB convention); the engine then uses the per-row trial count successes + failures as the binomial size. A weights = numeric vector of trial counts is also accepted as an alternative API (see Details). The negative-binomial families fit one log-dispersion parameter per trait: nbinom1() uses Var = mu * (1 + phi), while nbinom2() uses Var = mu + mu^2 / phi. The tweedie() family fits one log-phi and one logit-(p-1) per trait, with p in (1, 2). The Beta() and betabinomial() families each fit one log-precision per trait (a = mu*phi, b = (1-mu)*phi; Smithson & Verkuilen 2006). The student() family fits one log-sigma and one log(df-1) per trait (df > 1; pass student(df = 3) to fix df at a known value). The truncated_*() families require strictly positive integer responses (y >= 1). The delta families fit one log-dispersion of the positive component per trait (sigma_lognormal or phi_gamma).

weights

Optional numeric vector of length nrow(data). Family- conditional semantics (lme4 / glmmTMB convention):

For non-binomial families, weights is interpreted as per- observation likelihood multipliers — each row's log-likelihood contribution is multiplied by weights[i]. This matches the weights argument of lme4::lmer(), glmmTMB::glmmTMB() and stats::glm(). Use this for relative-abundance weighting (sqrt(p_si)), inverse-variance weighting (1 / var_i for heteroscedastic Gaussian fits), or to down-weight outliers. Must be non-negative and finite; weights[i] = 0 zeroes that row's contribution to the joint NLL.
For binomial families fit without a cbind(succ, fail) LHS, weights continues to be interpreted as the per-row trial count (binomial size; alternative API to cbind(succ, fail)). This pre-existing semantics is preserved on a per-row basis: in mixed- family fits, binomial rows use weights[i] as their trial count and non-binomial rows use it as the likelihood multiplier. Default NULL is equivalent to a length-nrow(data) vector of ones (unweighted Bernoulli on binomial rows). Wide matrix calls through gllvmTMB_wide() and wide data-frame calls through traits() normalise their accepted weight shapes to this same stacked response vector before fitting.

REML

Logical; use restricted maximum likelihood for Gaussian-only fits. The current REML pilot is deliberately narrow: all response rows must be Gaussian, observation weights are not supported, missing responses must use the default miss_control(response = "drop"), and mi() predictor models, predictor-informed latent(..., lv = ~ x) scores, and Xcoef_fixed maps are not supported. The observed fixed design must be full rank with positive residual degrees of freedom. The default FALSE keeps the historical ML fit.

mesh

Optional mesh object from make_mesh(). Required for any spatial_*() or spatial() term unless that term supplies its own mesh = argument; ignored only when the model has no spatial term.

phylo_vcv

(legacy global) Optional tip-only n_species × n_species phylogenetic correlation matrix. The canonical syntax is vcv = inside each phylo_*() keyword. — prefer tree =. The dense path inverts via Matrix::solve(), giving the same MLE as the sparse path but at O(n^2) memory and O(n^3) Cholesky cost. Use only when you have a Cphy in hand and no Newick tree (e.g. comparing against nlme::corPagel).

phylo_tree

(legacy global) Optional ape::phylo tree. The current syntax is to pass tree = inside each phylo_*() keyword (e.g. phylo_latent(species, d = K, tree = tree)); this argument is the older outer-level fallback and is slated for soft deprecation. When supplied (in either form) the phylo_*() terms build a sparse $\mathbf{A}^{-1}$ over tips + internal nodes natively, using only ape and Matrix – no MCMCglmm dependency. The construction is the deterministic Hadfield & Nakagawa (2010) sparse phylogenetic inverse (appendix eqs. 26-29), adopted from that method and ported from the sister package drmTMB (see inst/COPYRIGHTS). The result is ~5n non-zeros for an n-tip tree, vs n^2 for the dense path. This is the recommended path at any n_species; the speedup grows to ~24× at n_species = 1000.

known_V

Optional list of block-diagonal sampling-error matrices V_t. Used by the equalto() two-stage workflow when sampling-error matrices are available from a prior stage of estimation.

lambda_constraint

Optional list with elements unit and / or unit_obs (plus phy, spde for structural levels), each an n_traits × d matrix of confirmatory loading constraints (galamm-style). NA entries are estimated; numerical entries are pinned. Upper-triangle entries are silently ignored — the engine's lower-triangular parameterisation already fixes those at zero. Default NULL uses the engine's exploratory lower-triangular convention. See confirmatory_lambda() to build the matrix from functional-group membership, or suggest_lambda_constraint() for a minimum statistical identification scaffold when you have no biological hypothesis. Deprecated legacy element names B and W are still accepted (with a one-shot soft deprecation message) and map to unit and unit_obs respectively.

Xcoef_fixed

Optional named numeric vector of fixed-effect coefficient constraints. Names must match the expanded fixed-effect design columns (fit$X_fix_names); values must currently be 0, pinning those coefficients exactly at structural zero. Use this when a predictor is meaningful for some responses but should be fixed at zero for others. Native TMB fits use a parameter map; admitted engine = "julia" fixed-effect-X rows pass the same zero mask to GLLVM.jl. This is ML-only: REML = TRUE stops loudly. Native fitted objects still report all fixed-effect rows; pinned rows have estimate = 0, std.error = NA, and status = "fixed" in tidy(fit, "fixed").

control

Output of gllvmTMBcontrol().

missing

Output of miss_control() configuring missing-data handling. The default miss_control() (response = "drop", predictor = "fail") is the historical complete-case behaviour. miss_control(response = "include") keeps rows with a missing response and masks them out of the likelihood, preserving original-row accounting in fit$missing_data.

impute

Optional specification of the covariate model for a predictor declared missing with mi(x) in formula, used only when missing = miss_control(predictor = "model"). Supply a two-sided predictor-model formula (for Gaussian sugar, e.g. x ~ z) or an impute_model() object for an explicit predictor family. The default NULL is appropriate when no mi() term is present.

silent

Logical; suppress TMB and gllvmTMB chatter. Default TRUE.

engine

Character; "tmb" (default) fits with the native TMB engine, "julia" routes the fit through the experimental GLLVM.jl bridge fitting path via JuliaCall (see R/julia-bridge.R). The Julia path currently maps the unconstrained-ordination core (a single latent() block + per-trait intercepts) and errors on structures it does not yet support.

ci_method

Confidence-interval route requested at fit time for admitted engine = "julia" no-X rows. One of "none" (default), "wald", "profile", or "bootstrap". Native engine = "tmb" fits use confint() after fitting; non-default ci_* arguments therefore error unless engine = "julia". Grouped-dispersion rows, per-trait ordinal rows, response masks, mixed-family vectors, and fixed-effect-X rows remain gated.

ci_level

Nominal confidence level when ci_method != "none" on the Julia bridge.

ci_nboot

Number of parametric bootstrap replicates when ci_method = "bootstrap" on the Julia bridge.

ci_seed

Seed passed to the Julia bootstrap CI route.

site

(deprecated) alias for unit. Kept for backward compatibility. Use unit = ... in new code.

species

(deprecated) alias for cluster. Kept for backward compatibility. Use cluster = ... in new code.

Value

A gllvmTMB object. With no covariance-structure terms in the formula the result has class "gllvmTMB" (single-response engine); with latent() or other covstruct sugar it has class c("gllvmTMB_multi", "gllvmTMB") (multi-trait engine). Either way, S3 methods such as tidy(), predict(), vcov(), logLik() etc. dispatch on gllvmTMB; gllvmTMB_multi- specific methods (e.g. trait-level extract_ICC_site(), extract_communality()) are available for multi-trait fits.

Details

source \ mode	scalar	independent	dependent	latent
none	(omit)	`indep()`	`dep()`	`latent()`
animal	`animal_scalar()`	`animal_indep()`	`animal_dep()`	`animal_latent()`
phylo	`phylo_scalar()`	`phylo_indep()`	`phylo_dep()`	`phylo_latent()`
spatial	`spatial_scalar()`	`spatial_indep()`	`spatial_dep()`	`spatial_latent()`

The three covariance modes (indep / dep / latent) encode covstruct intent across traits:

latent — the decomposition mode $\boldsymbol\Sigma = \boldsymbol\Lambda \boldsymbol\Lambda^\top + \boldsymbol\Psi$: a low-rank shared part plus a diagonal $\boldsymbol\Psi$ companion. Ordinary latent() carries $\boldsymbol\Psi$ by default; the unique = argument controls it (latent(..., unique = FALSE) for the loadings-only subset, *_latent(..., unique = TRUE) to fold the $\boldsymbol\Psi$ companion into a source-specific term).
indep — the marginal-only mode: each trait gets its own variance, no cross-trait covariance.
dep — the full unstructured mode: $\boldsymbol\Sigma$ is free with $T(T+1)/2$ parameters via a Cholesky factor.

Plus the supporting phylo_slope() / animal_slope() (random slopes), meta_V() (known-V meta-analytic; meta_known_V() is a deprecated alias), and the engine-internal propto() / equalto() covstructs (used by the canonical keywords above; not typically called directly). See the formula keyword grid article for a one-paragraph tour of each cell.

gllvmTMB() parses the glmmTMB-style formula, converts wide traits() input to the same internal stacked-trait representation as explicit long data, and dispatches to the underlying TMB template. Covariance-structure terms (latent(), indep(), propto(), equalto(), spatial()) are processed by extending the formula parser and the TMB template.

Per the manuscript, when stacking traits one should set dispformula = ~ 0 so that no implicit residual variance competes with structured indep(0 + trait | …) terms. gllvmTMB() enforces this internally.

Multi-trial binomial. The TMB engine evaluates dbinom(y, n_trials, p) so binomial fits are not restricted to Bernoulli (size = 1). To pass a per-row trial count, use either:

cbind(successes, failures) ~ ... on the LHS of the formula (canonical R / glmmTMB convention; recommended), or
a flat response y with weights = n_trials (the alternative glmmTMB API).

Both interfaces produce identical fits. For Bernoulli data, omit weights and use a 0/1 response; the engine sets n_trials = 1 and the likelihood is identical to the previous Bernoulli-only behaviour.

Missing responses and predictors. NA response cells are treated as unobserved unit-trait cells. Under the default miss_control(response = "drop"), they are dropped before the TMB likelihood is built; with miss_control(response = "include"), they are kept in the model data, masked out of the likelihood, and available through predict_missing(). For cbind(successes, failures) binomial responses, a row is treated as missing when either response component is missing. Observation weights, when supplied, are subset to retained likelihood rows before validation. Ordinary missing predictors, grouping variables, or fixed-effect design values still error because the model cannot construct a design row; explicitly modelled missing predictors use mi(x), miss_control(predictor = "model"), and impute = list(...). These contracts are covered by the package's missing-data validation tests, both for responses and for predictors.

Delta (hurdle) families. delta_lognormal() and delta_gamma() use a single linear predictor for both components: presence is $\Pr(y > 0) = \mathrm{invlogit}(\eta)$ and the positive-component mean is $E[y \mid y > 0] = \exp(\eta)$. This matches sdmTMB's standard delta default (one fixed-effects matrix, one set of trait random effects). The shared-predictor scheme is the simplest hurdle formulation; a future release may decouple the two predictors so presence and abundance can have independent fixed and random effects. Only type = "standard" (the default) is currently wired in the multivariate engine; type = "poisson-link" is not implemented. Each delta family carries one per-trait dispersion of the positive component (no extra Bernoulli dispersion). The response must be non-negative.

Per-trait residual variance: when does it activate?

"Residual" in a mixed-effects model is scale-relative: what counts as residual variance shifts as you add levels to the model. In a Gaussian fit without a per-row diagonal term, the residual is row-level noise captured by a single shared sigma_eps. Once you add a per-row indep(0 + trait | obs) term, the row-level residual is now T per-trait random-effect variances and sigma_eps is auto-suppressed to avoid double-counting. If you also add a site-level diagonal term on top, the row-level term remains the residual and the site-level term is now an additional, higher-level random effect, not a residual at all. The dispatch table below records this explicitly for the configurations the engine supports today.

Gaussian and lognormal responses have one residual scale parameter, sigma_eps. Ordinary Gamma responses instead carry a per-trait shape phi_gamma (CV = 1 / sqrt(phi_gamma)). Per-trait residual variances only appear if you explicitly add a per-row indep(...) term. The dispatch is automatic:

No per-row indep, Gaussian/lognormal present: One shared sigma_eps across Gaussian/lognormal rows. Not per-trait.
No per-row indep, no Gaussian/lognormal rows: sigma_eps is mapped off; the family's intrinsic dispersion handles the residual (for ordinary Gamma, phi_gamma).
indep(0 + trait | g) where g has fewer levels than rows (e.g. g = "site"): sigma_eps is still estimated as the row-level residual; the diagonal term adds a per-trait random effect at level g on top.
indep(0 + trait | obs) at the per-row level (one level per row), Gaussian/lognormal rows fitted: sigma_eps is auto-suppressed (mapped off, fixed at a tiny stabiliser); the T per-trait diagonal random effects are the residual. A one-shot cli::cli_inform fires at fit time announcing the auto-suppression.
indep(0 + trait | obs) at the per-row level, non-Gaussian or mixed-family fit: Treated as observation-level random effects (OLRE). For Bernoulli traits the OLRE is statistically unidentifiable and is mapped off; for hurdle / delta families a warning is emitted (see "Per-family-aware OLRE selection" below).

Mnemonic: Gaussian/lognormal sigma_eps is the default; ordinary Gamma uses phi_gamma; per-row indep replaces the Gaussian/lognormal scalar residual with T per-trait residuals; non-per-row indep adds a higher-level random effect on top of the family residual; non-continuous families never carry sigma_eps regardless.

Per-family-aware OLRE selection. When indep(0 + trait | <unit_obs>) is at per-row resolution, i.e. one row per (trait, unit_obs) cell, the resulting per-trait random effects on the linear predictor are an observation-level random effect (OLRE). The engine now decides per trait what to do with the OLRE variance based on the trait's response family:

single-trial Bernoulli (binomial(), all rows have n_trials == 1): theta_diag_W[t] and the corresponding s_W column are mapped off. The OLRE is statistically unidentifiable here (no within-cell replication) and the MLE is the trivial boundary $\sigma_W \to 0$; pinning the parameter at $\sigma_W \approx 10^{-6}$ removes the spurious free parameter and emits a one-shot informational message. Multi-trial binomial (via cbind(succ, fail) or weights = n_trials) is identifiable and fit normally.
delta / hurdle families (delta_lognormal(), delta_gamma()): OLRE is fit as before, but a one-shot warning is emitted because the OLRE acts on the shared linear predictor of the hurdle and mixes presence and positive-component noise.
all other families: OLRE is fit as before.

In a mixed-family fit (per-row family list) the per-trait decision applies only to the affected trait. References: Nakagawa & Schielzeth (2010) Biol. Rev. 85: 935-956; Nakagawa, Johnson & Schielzeth (2017) J. R. Soc. Interface 14: 20170213.

Examples

if (FALSE) { # \dontrun{
set.seed(1)
sim <- simulate_site_trait(n_sites = 50, n_species = 8, n_traits = 3,
                           mean_species_per_site = 5)
fit <- gllvmTMB(
  value ~ 0 + trait + (0 + trait):env_1 + (0 + trait):env_2,
  data   = sim$data,
  family = gaussian(),
  trait  = "trait",
  unit   = "site"
)
summary(fit)
} # }