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 / kernel / temporal) and three modes (independent / dependent / latent). The common = TRUE modifier on *_indep() gives the one-shared-variance special case:

Usage

gllvmTMB(
  formula,
  data,
  trait = "trait",
  unit = NULL,
  unit_obs = NULL,
  cluster = NULL,
  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,
  estimator = c("ml", "mspl"),
  column_data = 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 three-mode grid covstructs above are supported (plus slope(), phylo_slope(), animal_slope(), kernel_slope(), spatial_slope(), column_coef(), phylo_coef(), animal_coef(), kernel_coef(), spatial_coef(), and meta_V()).

An offset() term is supported for count responses only — poisson(), nbinom1(), nbinom2(), truncated_poisson(), and truncated_nbinom2() — where it is the usual exposure or effort adjustment. Supply it already on the link scale, e.g. offset(log(trap_nights)). It must be its own additive term; it is not interacted with trait.

In long format the offset column carries one value per (unit, trait) row, so it can differ by trait. In wide format, offset(w) names one unit-level column and applies it to every trait, while offset(e1, e2, ...) names one column per trait, in traits() order.

A nonzero offset on a non-count trait is an error naming that trait. Zero is always allowed and does nothing (on the log scale it is a multiplier of one), which is how a mixed-family model gives an offset to its count traits and not the rest.

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. Required — every dataset has its own natural sampling-unit column, and gllvmTMB() aborts naming unit = ... if it is omitted and data has no "site" column to fall back to. For one release, omitting unit when data DOES have a "site" column still works via a deprecated implicit fallback (a one-time warning); pass unit = "site" explicitly to silence it. The implicit fallback is removed in 0.8.0.

unit_obs

Optional. 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 NULL, which resolves to "site_species" (the conventional name in joint species distribution modelling; safe for site × species data). If your data has no such column, omit unit_obs; the engine synthesises it from unit × cluster. For other domains pass e.g. unit_obs = "obs" for behavioural syndromes. If supplied but no covariance keyword uses that column, gllvmTMB() warns; omit the argument or use the column in the intended covariance keyword.

cluster

Optional. 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, and the historical non-trait-RHS phylo_slope route). The response-column form phylo_slope(x | trait, tree = tree) instead takes its column axis from the parsed RHS and does not use cluster. 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 NULL, which resolves to "species". For datasets with no third grouping the column may be absent from data entirely (with or without passing cluster); the engine will synthesise a placeholder. If supplied but no covariance keyword uses that column, gllvmTMB() warns; omit the argument or use the column in the intended covariance keyword.

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() and ordinal_logit() (the gllvmTMB-native ordinal threshold families with sigma_d fixed exactly – 1 for probit, pi^2/3 for logit; 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. Non-unit likelihood weights define a weighted estimating objective, not an ordinary maximum-likelihood fit. The package retains point estimates but refuses likelihood-based model-selection methods; Hessian/Wald uncertainty is not sandwich-corrected and must not be interpreted as calibrated inference for the weighted objective.
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, 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. Either missing-response policy may be used: masked cells (miss_control(response = "include")) contribute nothing to the likelihood, so the restricted likelihood is formed over the observed rows and matches the corresponding "drop" fit. The default FALSE keeps the historical ML fit.

REML is validated for all-Gaussian fits only. For non-Gaussian families an opt-in escape hatch exists (gllvmTMBcontrol(allow_nongaussian_reml = TRUE)) that realises the Cox-Reid adjusted profile likelihood rather than an exact restricted likelihood – this is unvalidated: a hypothesis under test, not an inherited result. Users estimating variance components or heritability on non-Gaussian data should be aware that plain Laplace can carry substantial downward bias on latent standard deviations. An opt-in AGHQ engine exists (gllvmTMBcontrol(aghq = k)) with its own, separately-scoped evidence; AGHQ and REML = TRUE cannot be combined.

mesh

Optional mesh object from make_mesh(). Required for any spatial_*() or spatial() term unless that term supplies its own mesh = argument. A supplied mesh with no spatial term is unused: the fit warns, and a raw fmesher/INLA mesh is rejected, rather than being silently ignored.

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 GLLVModels.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") omits each missing response cell: a unit keeps every trait it does have, so this is not case-wise deletion. miss_control(response = "include") instead keeps those cells and masks them out of the likelihood, preserving original-row accounting in fit$missing_data; it reaches the same optimum as the default and differs in which cells predict_missing() can return.

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

estimator

Estimation criterion. "ml" (the default) preserves the ordinary maximum Laplace-likelihood route. "mspl" opts into the experimental maximum softly penalised Laplace-likelihood point estimator for the currently admitted complete-Bernoulli surface. This means exactly one ordinary latent() block with no free Psi and d = 1 or 2, standalone spatial_indep(), or standalone spatial_latent() with d = 1 or 2; one common logit, probit, or complementary-log-log link; complete single-trial responses; a full-rank resolved fixed-effect design; all-zero offsets; and the native Laplace engine. Unsupported structures fail before optimisation. MSPL does not activate automatically after a separation warning and never changes the default ML fit.

column_data

Optional keyed response-column metadata. Its key column must have the same name as trait and exactly match the response columns. Joined fields are available to fixed effects but cannot be coefficient bases, grouping columns, or covariance sources. In wide input, key this table by the synthetic column name trait and omit the trait argument.

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). S3 methods are registered on the subclasses, not on bare gllvmTMB: for a multi-trait fit (gllvmTMB_multi), tidy(), predict(), summary(), confint() and logLik() are available, and AIC() / BIC() dispatch through logLik() only when the fitted point is an unpenalised native-Laplace ML maximum. AGHQ uses a different integration objective and a loading ridge returns a penalised MAP point, so ordinary likelihood comparisons do not apply. An MSPL result has leading class "gllvmTMB_mspl"; likelihood comparison, standard errors, intervals, profiles, and hypothesis-test methods fail closed because this experimental surface is point estimation only. Trait-level extractors such as extract_ICC_site() and extract_communality() are multi-trait only.

vcov() and coef() are available too — see gllvmTMB_multi-vcov. coef() works on any fit; vcov() needs the sdreport() and raises the same typed errors confint() does when it is missing or non-finite.

History, 2026-08-04: this block used to say vcov() dispatched on gllvmTMB, which was never true of any release — it and coef() existed only for gllvmTMB_va, where they refuse. The wording was corrected first, then the two methods were added, which is why the promise now holds.

Details

source \ mode	independent	dependent	latent
none	`indep()`	`dep()`	`latent()`
animal	`animal_indep()`	`animal_dep()`	`animal_latent()`
phylo	`phylo_indep()`	`phylo_dep()`	`phylo_latent()`
spatial	`spatial_indep()`	`spatial_dep()`	`spatial_latent()`
kernel	`kernel_indep()`	`kernel_dep()`	`kernel_latent()`
temporal	`temporal_indep()`	`temporal_dep()`	`temporal_latent()`

The temporal row is currently a Gaussian identity-link ML/Laplace route: it accepts one ordered series–time provider, AR1 integer occasions or OU elapsed time, and may coexist with ordinary unit and unit_obs covariance. Other structured sources, new-data prediction, intervals, and rank above one remain unavailable for temporal fits.

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.

The five Gaussian long-format response-column slope helpers are slope(), phylo_slope(), animal_slope(), kernel_slope(), and spatial_slope(). Their RHS is the resolved response-column factor, so they are not extra grid modes. Plus 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.

Response-column coefficient helpers are also outside the grid. column_coef() fits IID random intercept/slope bases across response columns; phylo_coef() fits the same bases with fixed or estimated phylogenetic correlation strength; animal_coef() fits them from a pedigree, relationship covariance, or relationship precision with fixed correlation strength; kernel_coef() uses one labelled dense kernel with fixed or estimated strength; and spatial_coef() uses one labelled response-column mesh with rho = 1. All five support Gaussian point models in long and traits(...) wide form. Non-Gaussian coefficient models, spatial IID mixtures, and interval inference remain unavailable.

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(...). The modelled-mi() route is native Laplace only; the variational route refuses mi() predictors. Dense missing-response masks have their own narrower VA fence and do not imply predictor-model support. 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 within-family 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 each have a residual scale parameter, sigma_eps, shared within family. Pure fits retain one slot; a joint Gaussian-lognormal fit has separate raw- and log-scale slots. 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 slot in a pure fit; separate Gaussian raw-scale and lognormal log-scale slots when both coexist. Each slot is shared within family, 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 family 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.

Dispersion parameters in the report slot

The fitted object's $report carries per-family dispersion quantities on the natural scale, each a vector of length n_traits indexed by trait (except sigma_eps, which has one slot normally and two slots only when Gaussian and lognormal responses coexist). Their names follow the engine's internal parameterisation, not standard R distribution arguments — several are easy to misread. What each one IS, and how it maps to standard R arguments (mu = exp(eta) unless noted):

Family (report name)	What it is	Standard-R conversion
gaussian, lognormal (`sigma_eps`)	residual SD / sdlog; one within-family shared slot in pure fits, or separate Gaussian raw-scale and lognormal log-scale slots when both coexist	`pnorm(y, mean = eta, sd = sigma_eps[1])`; `plnorm(y, meanlog = eta, sdlog = sigma_eps[2])` in a joint fit
Gamma (`phi_gamma`)	the shape, not a dispersion	`pgamma(y, shape = phi, scale = mu / phi)`; `CV(y) = 1 / sqrt(phi)`
nbinom2 (`phi_nbinom2`)	the NB `size`	`pnbinom(y, size = phi, mu = mu)`; `Var(y) = mu + mu^2 / phi`
nbinom1 (`phi_nbinom1`)	linear overdispersion, `Var(y) = mu * (1 + phi)`	`pnbinom(y, size = mu / phi, mu = mu)` — the size is mean-dependent
tweedie (`phi_tweedie`, `p_tweedie`)	dispersion and power `p` in (1, 2)	no base-R CDF; `tweedie::ptweedie(y, mu = mu, phi = phi, power = p)`
Beta (`phi_beta`)	the precision	`pbeta(y, shape1 = mu * phi, shape2 = (1 - mu) * phi)` with `mu = plogis(eta)`
betabinomial (`phi_betabinom`)	precision of the Beta mixing	`a = mu * phi`, `b = (1 - mu) * phi` with `mu = plogis(eta)`; no base-R CDF
student (`sigma_student`, `df_student`)	the scale, not the SD	`pt((y - eta) / sigma, df)`; `SD(y) = sigma * sqrt(df / (df - 2))`, undefined for `df <= 2`
truncated_nbinom2 (`phi_truncnb2`)	the NB `size`; a separate vector from `phi_nbinom2`	`pnbinom(y, size = phi, mu = mu)` renormalised by `1 - pnbinom(0, ...)`
delta_lognormal (`sigma_lognormal_delta`)	sdlog of the positive part	`plnorm(y, meanlog = eta, sdlog = sigma)` for `y > 0`; `P(y > 0) = plogis(eta)`
delta_gamma (`phi_gamma_delta`)	the CV of the positive part, not a shape	`pgamma(y, shape = 1 / phi^2, scale = mu * phi^2)` for `y > 0`; `P(y > 0) = plogis(eta)`

Objective provenance

joint_nll report fields vs objective_components: A fitted object carries two different surfaces that both look like "the objective", and they answer different questions.

report$joint_nll_unpenalized and report$joint_nll_penalized are conditional joint quantities computed on the TMB tape (the negative log-likelihood of the data and random effects at their current values, before and after TMB-side MSPL penalty terms are added). They see only penalties applied inside the C++ template – the LA-MSPL penalties documented above – and never see anything computed at the R level.

fit$objective_components is the marginal surface: the value TMB's Laplace approximation integrates the random effects out of, decomposed into likelihood_nll (the unpenalised marginal negative log-likelihood), ridge_penalty (the R-level loading ridge from aghq_ridge, added outside the TMB objective – see .gllvmTMB_penalised_gradient()), optimization_nll (their sum, the objective the optimiser actually minimised), and optimizer_reported (the value nlminb() / optim() returned). This is the only place the R-level loading ridge appears; report$joint_nll_* cannot see it because it never reaches the TMB tape.

The two surfaces are not expected to match, on a ridged fit or otherwise, and a difference between them is not itself an error: they are conditional-joint and marginal quantities from different penalty regimes.

References

Sterzinger, P. and Kosmidis, I. (2023). Maximum softly-penalized likelihood for mixed effects logistic regression. Statistics and Computing, 33, 53. doi:10.1007/s11222-023-10217-3 . This paper is the logistic mixed-effects antecedent for the fixed-effect Jeffreys component. The probit, complementary-log-log, reduced-rank GLLVM, and spatial penalties implemented here are gllvmTMB-specific extensions; that paper does not validate those extensions.

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)
# Long format: one response value per site-species-trait row.
fit_long <- gllvmTMB(
  value ~ 0 + trait + (0 + trait):env_1 + (0 + trait):env_2,
  data   = sim$data,
  family = gaussian(),
  trait  = "trait",
  unit   = "site_species"
)

# Wide format: one row per site-species, with the same three traits.
# `traits()` expands `1 + env_1 + env_2` to the explicit long terms above.
wide_values <- stats::reshape(
  sim$data[c("site_species", "trait", "value")],
  idvar = "site_species", timevar = "trait", direction = "wide"
)
wide <- merge(
  unique(sim$data[c("site_species", "env_1", "env_2")]), wide_values,
  by = "site_species", sort = FALSE
)
names(wide) <- sub("^value[.]", "", names(wide))
fit_wide <- gllvmTMB(
  traits(trait_1, trait_2, trait_3) ~ 1 + env_1 + env_2,
  data   = wide,
  family = gaussian(),
  unit   = "site_species"
)
summary(fit_long)
summary(fit_wide)
} # }