Skip to contents

GLLVM means generalized linear latent-variable model. A GLLVM analyses several responses together and uses a small number of unobserved shared patterns, called latent variables, to describe how those responses vary together. The word generalized means the model can use response distributions suited to different data, such as continuous measurements, counts, or presence and absence. See Response families for the supported choices and their limits.

gllvmTMB is the R package used in these guides. This glossary translates terms from multivariate statistics, ordination, factor analysis, and mixed models, and points to the page where each idea is demonstrated.

Use Choose your next gllvmTMB step when you know the task you need to complete but not the function or guide that starts it.

These definitions describe the current public interface. They do not imply that every combination of response family, covariance source, and inferential method has equal evidence.

Data and formula terms

Trait and response matrix

A trait is one response variable in the multivariate model. Depending on the application, traits may be species occurrences, morphological measurements, questionnaire items, or study outcomes. Placing units in rows and traits in columns gives the response matrix.

Unit and grouping key

A unit is the entity carrying the multivariate response: for example, a site, individual, species, or study. A grouping key is a data column that identifies which observations share a random-effect or covariance component. In gllvmTMB(), unit identifies the between-unit level, unit_obs identifies within-unit replicate cells, cluster supplies a further ordinary or structured grouping, and optional cluster2 supplies an additional plain-diagonal grouping. These are modelling roles rather than fixed biological meanings: the columns that fill them depend on the study design. Deciding which column of your own data should play the unit role — a checklist, a station, a grid cell — is the subject of What is a site? The unit of analysis.

Long and wide data

Stacked long data have one row per (unit, trait) observation, with a response column such as value, a trait factor, and the grouping columns. This is the package’s internal representation and the most explicit formula form.

Wide data have one row per unit and one column per trait. The traits(...) left-hand side converts that shape to the same stacked representation. These two formulas express the same two-LV model. Long calls name both the trait and unit columns explicitly. Wide calls omit trait = because traits(...) supplies the response columns, but still name unit. Name the unit column explicitly: the temporary implicit "site" fallback is deprecated, and omission otherwise stops with an error.

# Long
gllvmTMB(
  value ~ 0 + trait +
    latent(0 + trait | site, d = 2),
  data = dat_long,
  trait = "trait",
  unit = "site"
)

# Wide
gllvmTMB(
  traits(t1, t2, t3) ~ 1 +
    latent(1 | site, d = 2),
  data = dat_wide,
  unit = "site"
)

The family is the probability distribution used for a response, such as Gaussian, binomial, or Poisson. The link function maps the response mean to the linear-predictor scale: identity for the usual Gaussian model, logit for a common binomial model, and log for a common Poisson model. See Response families before combining families in one fit.

Latent-variable terms

Latent variable and latent score

A latent variable (LV) is an unobserved dimension used to represent shared variation among traits. d = 2 asks for two LVs, conventionally labelled LV1 and LV2.

A latent score is a particular unit’s fitted position on those LVs. In the ordinary latent model, unit ii has score vector 𝐮i\mathbf u_i with the identifying convention Var(𝐮i)=𝐈d\operatorname{Var}(\mathbf u_i)=\mathbf I_d before conditioning on its responses. Fitted scores are conditional estimates from the likelihood, responses, predictors, and covariance structure; they are not directly observed measurements. Use extract_ordination() to obtain scores and loadings together.

Loading and loading matrix

A loading Λtk\Lambda_{tk} describes how strongly trait tt changes with LV kk. The T×dT \times dloading matrix 𝚲\boldsymbol\Lambda contains one row per trait and one column per LV. Large absolute values indicate a strong relationship with that fitted direction on the linear-predictor scale. Their magnitudes are not directly comparable across response families or links. Individual loading signs and axes also depend on orientation, so interpret rotation-invariant summaries first.

Shared low-rank covariance

The product

𝚲𝚲𝖳 \boldsymbol\Lambda\boldsymbol\Lambda^{\mathsf T}

is the shared low-rank covariance generated by the LVs. It is unchanged by an orthogonal rotation of the loading matrix and is therefore safer to compare across refits than raw loading columns.

Psi: tier-specific diagonal variance

𝚿=diag(ψ1,,ψT)\boldsymbol\Psi = \mathrm{diag}(\psi_1, \ldots, \psi_T) is a diagonal covariance component at a named model tier. It captures trait-specific variance at that tier that is not shared through the LVs. Ordinary latent() includes this diagonal companion by default; latent(..., unique = FALSE) requests the loadings-only subset.

Psi is not a generic synonym for every residual variance. In particular, it is distinct from Gaussian observation residual variance, a diagonal component at another grouping tier, and a distribution-specific link-scale residual used by some extractors.

Sigma: model-implied trait covariance

For the ordinary latent decomposition at one tier,

𝚺=𝚲𝚲𝖳+𝚿. \boldsymbol\Sigma = \boldsymbol\Lambda\boldsymbol\Lambda^{\mathsf T} + \boldsymbol\Psi.

𝚺\boldsymbol\Sigma is the model-implied trait covariance at that tier. Use extract_Sigma() for matrices and extract_Sigma_table() for report-ready rows. See Covariance and correlation for the worked interpretation.

Communality

Communality is the fraction of a trait’s variance at a chosen tier and scale that is attributed to the shared LV component:

ct2=(𝚲𝚲𝖳)ttΣtt. c_t^2 = \frac{(\boldsymbol\Lambda\boldsymbol\Lambda^{\mathsf T})_{tt}} {\Sigma_{tt}}.

It is not automatically the fraction of total phenotypic variance across the whole model. If Psi is excluded by construction, communality can equal one mechanically; that is not evidence that the traits are perfectly integrated.

Model-based ordination

Model-based ordination displays fitted latent scores and trait loadings from a probability model. Unlike PCA or NMDS, its mean structure and residual behavior depend on the chosen response family and link. The displayed axes are still orientation-dependent, so conclusions should be checked against rotation-invariant covariance summaries.

Covariance-structure terms

Covariance tier

A covariance tier is one variance component in the model, attached to a grouping key and interpreted at that level. Between-unit and within-unit components are different tiers even when they use the same covariance mode. Do not use tier, grouping key, correlation source, and LV axis as synonyms.

Correlation source

The correlation source says how levels of a grouping factor are related: ordinary independent levels, pedigree/animal relatedness, phylogenetic relatedness, spatial dependence, or a supplied dense kernel. Phylogenetic and spatial labels describe structured covariance sources or variance-share categories; they are not peer data columns alongside unit and unit_obs.

Covariance mode

The covariance mode describes the trait-by-trait covariance within a tier:

  • indep: separate trait variances with zero cross-trait covariance;
  • dep: a full unstructured trait covariance;
  • latent: reduced-rank shared covariance, with the diagonal Psi convention determined by the keyword family and its unique argument.

The canonical source-by-mode interface crosses ordinary, animal, phylogenetic, and spatial sources with three modes: independent, dependent, and latent. common = TRUE on *_indep() is the one-shared-variance modifier; unique = TRUE on *_latent() adds the diagonal Psi companion. A separate dense-kernel interface provides kernel_indep(), kernel_dep(), and kernel_latent(). See the API keyword grid for exact syntax.

Scale and identification terms

The covariance fitted for a random-effect tier is not always the same estimand as a marginal or liability-scale covariance. With link_residual = "none", an extractor reports the fitted model covariance. With link_residual = "auto", supported extractors add a family-specific diagonal convention before standardising or partitioning variance. This is a post-fit scale choice, not an extra likelihood term, and some family-specific additions are approximations.

Rotation equivalence and loading constraints

If 𝐐\mathbf Q is orthogonal, 𝚲\boldsymbol\Lambda and 𝚲𝐐\boldsymbol\Lambda\mathbf Q imply the same shared covariance. The identified object is therefore an equivalence class—or its covariance—not a unique raw loading orientation. A loading constraint or post-hoc rotation chooses an orientation; it does not create information that the data lack.

Weak identification

A parameter is weakly identified when the likelihood changes little over a range of substantively different values. Large standard errors, unstable refits, flat profiles, or collapsing extra LVs can be warning signs. Repeated starts, refits, and likelihood profiles can reveal instability under a fitted scenario, but they do not by themselves prove that the original data identify a particular rank. Testing a rank-selection procedure requires simulations with a predeclared generating rank and a denominator that retains failed or unusable fits.

Boundary

A boundary is the edge of the parameter space, such as a variance equal to zero or a correlation near 1-1 or 11. Ordinary quadratic approximations can be unreliable there. Profile intervals can display asymmetric or one-sided likelihood information, but their reference distribution and calibration are target-specific. See Profile likelihood: routes and fallbacks.

Fit and uncertainty terms

Convergence, gradient, and Hessian

Optimizer convergence means the numerical optimizer stopped according to its criterion. The gradient measures the remaining local slope of the objective. The conservative package verdict requires an unscaled maximum gradient below 0.01; scaled_gradient is descriptive and cannot rescue a failed optimiser or a large raw gradient. The Hessian measures local curvature. A positive-definite Hessian supports local standard-error calculations, but does not prove identifiability, model adequacy, or interval calibration. Use Fit diagnostics before interpreting estimates.

Wald, profile, bootstrap, and Fisher-z intervals

A Wald interval uses local Hessian/delta-method curvature on a chosen working scale. A profile interval inverts a likelihood profile. A parametric bootstrap interval simulates and refits from the fitted model. A Fisher-z interval transforms a correlation with atanh(ρ)\operatorname{atanh}(\rho), uses the approximate standard error 1/neff31 / \sqrt{n_{\mathrm{eff}} - 3}, and back-transforms the limits to [1,1][-1, 1]. In mixed and non-Gaussian latent models, treat this as a sensitivity interval based on an explicit or tier-derived effective sample size, not as Hessian-based inference or a coverage guarantee.

None is universally best. The implemented route, fallback, cutoff, and evidence depend on the target. Inspect method or status fields where the extractor provides them, together with its messages and warnings, rather than assuming the requested method was used.

Missing-data terms

Response mask

A response mask retains the identity of a missing response cell while giving that cell zero contribution to the response likelihood. With miss_control(response = "include"), predict_missing() later returns a plug-in fitted value for each masked cell; it does not currently return a prediction interval.

Modelled missing predictor and conditional mode

mi(x) declares a modelled missing predictor and impute specifies its predictor model. The predictor may live at the unit level, a coarser mi_group() level, or a phylogenetic-tip level. For a continuous Gaussian predictor, imputed() returns a fitted conditional mode and, when available, a local Gaussian/Laplace conditional standard error from the joint Hessian. It is not a conditional mean, a multiple-imputation draw, a Bayesian posterior summary, or a calibrated interval. See Handling missing data.

Species distribution and integrated-data terms

These terms appear across the species-distribution articles (see Start here for the map of those pages).

Intensity

The intensity is the expected number of individuals or records per unit of area and effort — the “surface of expected abundance” those articles are built around. The models there return it up to an unknown overall constant (relative intensity), which is enough for maps and habitat effects but not for absolute abundance.

Arm

An arm is one data source inside an integrated model. A model fed eBird checklists and a structured survey is a two-arm model: one shared surface, each arm keeping its own effort, distortion, and observation type.

Offset

An offset is a predictor whose coefficient is fixed at 1 rather than estimated, used for known scalings: offset(log_effort) tells a count model that doubling the effort doubles the expected count, instead of asking it to discover that.

Point process and quadrature points

A point process models the locations of records in continuous space rather than counts in predefined units — the natural shape for presence-only data. Fitting one in practice lays quadrature points (also called background or dummy points) over the region so the likelihood can be computed as a weighted regression.

cloglog

The cloglog (complementary log-log) link makes a detection probability consistent with an underlying intensity: if individuals are found at Poisson rate μ\mu, the chance of seeing at least one is 1exp(μ)1 - \exp(-\mu), which is exactly what a binomial model with a cloglog link fits. It is the hinge that lets detection/non-detection data share a surface with count data.

Mesh

A mesh is a triangulation of the study area on which a smooth spatial field is represented: the field gets a value at each mesh node and is interpolated between them. Built once with make_mesh() from projected coordinates; the articles’ rule of thumb sets the cutoff near one tenth of the domain’s shorter side.

Where to go next