Current limitations and boundaries

Can I use gllvmTMB for this analysis?

gllvmTMB models several responses measured on the same sampling units, such as species counts at sites or traits measured on individuals. A generalized linear latent-variable model (GLLVM) uses a few unobserved variables to describe how responses vary together after accounting for measured predictors. It can help ask which species respond similarly to an environmental gradient, or how much variation is shared among traits. Shared variation alone does not identify a cause or demonstrate an interaction between species.

For a first runnable example, follow Get started: fit continuous morphological traits. It loads data included with the package, fits a model, checks the fit, and reads trait covariances. Both wide and long data are shown. This is a teaching example; running it successfully does not establish reliability for a new dataset.

gllvmTMB is experimental. Its formulas and summaries may change, and support for a model does not mean that all its estimates and intervals have been shown to be reliable. Before reporting a result:

• Match the response model to your data. Continuous measurements, counts, and presence/absence need different choices; see the response-family guide.
• Check the fitted model. A fit can appear numerically healthy and still give poor estimates, particularly for some negative-binomial count models and models with missing responses.
• Read the evidence for your model in the table below. Results for a simple model do not automatically apply after adding a phylogeny, spatial structure, or another source of variation.
• Treat uncertainty separately from point estimates. An available confidence interval does not guarantee that it contains the true value as often as its stated confidence level implies. See Current uncertainty limits.

For continuous traits, the strongest starting evidence here is for Gaussian models with independent trait variation (indep()) or a full trait covariance (dep()), in the tested three-trait complete-data examples. The Gaussian latent model used in the introductory tutorial has not passed all its recovery checks. For counts and binary responses, use the matching row below; the Gaussian results do not establish reliability for those data.

Choose a model within the evidence

This page describes the 0.7.1 source version. It separates three practical decisions:

• A safer starting point within the tested conditions: point-estimate checks passed for the model and data design stated. This does not establish reliability for other designs or for confidence intervals.
• Exploratory only: you can fit and describe the stated model, but its accuracy is incomplete or a relevant check failed.
• Do not use it for this purpose: the requested conclusion is unsupported, or known failures make that interpretation unreliable.

Match the whole model to the evidence: its response family, number of responses, sample size, covariance structure, and missing-data pattern all matter. Adding a phylogeny, spatial term, or other component does not inherit the evidence for the simpler model. The technical record behind these checks is available in Validation oracles.

Choose a response model

For the ordinary latent route, latent(..., unique = TRUE) is the default and represents $\Sigma = \Lambda\Lambda^\mathsf{T} + \Psi$, where $\Psi$ is diagonal. Setting unique = FALSE requests the narrower loadings-only subset.

Your model or question	What you can use	Important limit
Gaussian indep()	Safer starting point: trait-specific variance point estimates after fit-health checks, for complete three-trait data in the tested n = 60 and n = 240 designs	No calibrated uncertainty, broad robustness claim, or guarantee at other sample sizes or boundary conditions
Gaussian dep()	Safer starting point: a full trait-covariance point estimate after fit-health checks, for complete three-trait data in the tested n = 60 and n = 240 designs	A positive-definite Hessian does not show that every covariance is accurate, and this evidence does not transfer to structured models
Gaussian latent()	Exploratory: rotation-invariant covariance (Sigma), correlations, communality, and an ordination after diagnostics	The n = 60 central design was not ready for release and correlation-stress designs failed accuracy checks. Do not attach biological meaning to raw loading signs or orientation, or assume reliability near weak-variance or boundary cases
Poisson-log rank-1 latent()	Exploratory point estimates: rotation-invariant covariance summaries for clean three-trait data in the tested n = 100 and n = 300 designs	A rare-level design had excessive fixed-effect bias. This does not establish general robustness, interval reliability, or dependable failure detection
NB2-log rank-1 latent()	Exploratory: report point results with their exact sample size and dispersion diagnostics; the n = 300 central design passed	Do not make a dependable point-estimation claim: n = 100 and missing-response designs had severe dispersion errors that diagnostics often missed
NB2-log latent model with missing responses	Do not use for dependable results. The tested model is useful only for studying the failure mechanism	Standard optimizer, Hessian, and boundary diagnostics missed 30 of 41 severe truth errors in the tested missing-response design
Binomial(10)-logit rank-1 latent()	Exploratory point estimates: rotation-invariant covariance summaries for clean three-trait data in the tested n = 100 and n = 300 designs	No broad robustness claim, Bernoulli equivalence, separation protection, or reliability in the small-sample designs where failures could be silent
Other response families	Exploratory: only within the conditions in the response-family and model-specific guides	Constructing a family, or obtaining one successful fit, does not establish broad recovery
Ordinal models	Do not use for a validated conclusion. Exploratory point fits need model-specific validation	Nominal convergence does not rule out a degenerate fit, and binomial diagnostic thresholds do not apply to ordinal fits
Mixed response families across traits	Exploratory: documented combinations with one family and link per trait	Do not directly compare correlations that live on incompatible latent or response scales
Family or link varying between rows of one trait	Do not use. The fit now stops before optimisation	An accepted per-row family identifier does not make a shared loading coherent across Poisson-log and binomial-logit/probit. The planned Poisson-log/binomial-cloglog integrated route is not yet validated

Check additional model components

More complex structure: covariance, time, groups, and missing data

Your model or question	Availability	What you can use	Important limit
Non-unit likelihood weights	Do not use for likelihood-based inference	Weighted-objective point estimates only, with the weighting rule reported	Ordinary Hessian/Wald uncertainty, logLik(), AIC, BIC, or likelihood-ratio interpretation; no sandwich route is certified
Phylogenetic, animal, spatial, or kernel covariance	Exploratory	Point estimates for the exact documented route after structure-specific diagnostics	That ordinary-model evidence transfers to a structured source, that spatial tau or field SD has validated absolute-scale recovery, or that structured intervals are calibrated
Gaussian temporal covariance	Exploratory	Use the temporal guide only for repeated Gaussian measurements with complete panels, at least three traits, and three occasions per series	General recovery, calibrated intervals, temporal slopes, higher rank, or non-Gaussian support
Temporal plus another source	Exploratory	A narrow Gaussian AR1 model plus one kernel, phylogenetic, animal, or spatial source; see the temporal guide for the required complete-panel design	Broad recovery or coverage, OU combinations, multiple structured sources, source-by-time interactions, and spatial forecasting
Random slopes or several grouping tiers	Exploratory	Only the exact family, mode, slope-count, and grouping combinations with documented recovery evidence, with the covariance block named explicitly	A grid-wide claim across families, modes, slope counts, or grouping levels
Missing response cells	Exploratory	Point estimates after confirming the intended cells were omitted from the likelihood; the strongest worked evidence is Gaussian	That every family has equal recovery depth or that missingness is ignorable by default
One modelled missing predictor with mi()	Exploratory	Native-Laplace results for one explicitly modelled predictor in the documented predictor-family routes	VA support, several mi() predictors, general MNAR handling, or missing grouping variables

Alternative fitting methods and the optional Julia connection

Your model or question	Availability	What you can use	Important limit
Variational approximation	Exploratory	Opt-in exploratory point estimates for the documented families and ranks, using the current automatic GH evaluator with H = 7, and report the method metadata and exact family/rank	General point recovery outside those families and ranks; its ELBO is not a log-likelihood, so ordinary logLik(), AIC, BIC, and likelihood-ratio comparisons are undefined; do not infer calibrated covariance, Wald intervals, or evidence equivalent to Laplace
Adaptive Gauss-Hermite quadrature	Exploratory	Opt-in sensitivity analysis in its documented scalar-random-effect regime, stating the ridge setting	Ordinary ML interpretation: AGHQ defaults to a tau = 2 penalised MAP fit, so ordinary AIC, BIC, and likelihood-ratio interpretations do not apply; the evaluated node/ridge rule is not a generally calibrated estimator
EVA	Do not use for this claim	No public estimator is available	That an implementation detail implies an available estimator
Julia engine	Exploratory	Optional-bridge results for the documented model types; start with Julia bridge setup	Full R/TMB feature parity or a required dependency for the ordinary workflow

Specialised response-column, integrated, and multinomial routes

Your model or question	Availability	What you can use	Important limit
Response-column coefficient models (column_coef(), phylo_coef(), animal_coef(), kernel_coef(), spatial_coef())	Exploratory	Gaussian point estimates for a random intercept, slope, or both, varying across response columns, fitted through long or traits(...) wide data. phylo_coef() and kernel_coef() can fix a numeric mixture strength or estimate one interior value; animal_coef() fixes it; spatial_coef() fixes it at 1 and estimates the spatial range	Calibrated intervals, non-Gaussian responses, an estimated mixture strength for the animal or spatial routes, or recovery outside the tested Gaussian cells
Prediction at new locations for integrated (iSDM) fits (predict() on isdm_sources() fits)	Exploratory	Point predictions, including maps, at new in-hull coordinates for the ordinary intercept-only spatial tier; training-row predictions reproduce the fitted values exactly	Held-out map accuracy, which has not been demonstrated in the evaluated design; a standard error or confidence interval on newdata predictions (se.fit is refused there); or predictions for spatial slope tiers, which remain unsupported
Predictor-informed latent scores (latent(..., lv = ~ x))	Exploratory	Point estimates of the rotation-invariant effect of a predictor acting through the shared latent axes, combining traits from more than one response family in one model, in the specific tested construction and evidence cells	Broad recovery or calibration across arbitrary family combinations; support for missing latent-score predictors, combined fixed-and-latent-score predictors, REML, or a structured (phylogenetic/spatial/kernel) source on this term
Multinomial structured routes	Exploratory	Fixed-effect point estimates and probability predictions for one unordered multinomial response are well tested. Several structured extensions (a phylogenetic route, a shared-latent cross-family route, pedigree- or kernel-sourced routes, and ordinary group random intercepts) exist and each has its own narrower tested evidence	Fitting more than one multinomial response in the same model; augmented (intercept-and-slope) structured forms; or calibrated intervals on any structured multinomial route
Estimated source strength (rho = NULL on phylo_*/animal_*/kernel_*/spatial_* terms)	Exploratory	Bounded point-recovery evidence for a native Gaussian model with complete replicated trait vectors, one known relatedness source, retained observation-level residual variance, and no competing ordinary covariance term	General recovery across other source/covariance-mode combinations — most combinations tested so far did not pass; any confidence interval on the estimated strength, which does not yet exist; or recovery for augmented (intercept-and-slope) forms

What to interpret first

For latent models, start from quantities that do not change when the latent axes are rotated:

1. total covariance, Sigma;
2. correlations derived from Sigma;
3. shared covariance, Lambda Lambda^T;
4. the diagonal unique component, Psi;
5. communality.

Raw loading columns can change sign, order, or orientation without changing the fitted covariance. Use them for an ordination only after applying a stated rotation convention; do not treat an arbitrary raw orientation as biology.

Current uncertainty limits

Wald, bootstrap, and profile methods are available for named targets, but their availability does not establish repeated-sampling coverage. Report the target, method, and returned status; calibrated uncertainty is not recommended as a general claim. The narrow loading-interval exception is described below.

profile_ci_total_variance() is a penalty-profile approximation, not a certified exact likelihood-ratio profile. Historical ordinary Gaussian n_units = 150, d in {1, 2} cells and the newer n_units = 400, d = 2 cell met their prespecified numerical coverage criteria, but the retained endpoints do not prove constrained-refit convergence and exact target attainment. Every computed row remains an approximate calculation rather than a certified exact profile. The earlier 0.94 coverage result does not justify exact-profile language until the mechanism is repaired and recalibrated.

There is one exact loading-interval exception. Standardized symmetric joint-delta Wald intervals are certified for the structurally free strict-lower targets in native pinned unrotated ordinary Gaussian three-trait cells (n_units=150,d=2), (n_units=400,d=1), and (n_units=400,d=2), at 95%. This is one frozen DGP: trait intercepts (-0.20, 0.10, 0.25), unique standard deviations (0.70, 0.80, 0.90), and loading vector (0.80, 0.45, -0.35) for d=1, with second column (0, 0.70, 0.40) for d=2. Coverage is conditional on eligible fits (optimizer convergence, converged fit health, available sdreport(), and a positive-definite Hessian); availability was 98.82%, 93.38%, and 96.18% in the three certified cells. The (n_units=150,d=1) cell failed. No other loading, unique-variance, or intercept regime inherits the result. Pinned diagnostic rows, Fisher-z Wald, arbitrary constraints, rotated fits, and every neighbouring family, tier, rank, sample size, or confidence level remain uncalibrated.

For spatial and phylo_diag models, current simulation and bootstrap fallbacks may condition on fitted random effects. The returned intervals can therefore be too narrow. Read the warning and do not report those intervals as calibrated unconditional uncertainty.

VA interval methods remain labelled calibrated = FALSE. AGHQ’s node and ridge choices were evaluated for numerical failure and runaway avoidance in a narrow regime, not for general interval coverage.

Check every fit before interpretation

At minimum, run:

check_gllvmTMB(fit)
diagnostic_table(fit)

Inspect optimizer status, raw and scaled gradients, Hessian status, boundary flags, weak latent axes, non-finite estimates, restart history, and the diagnostics specific to the response family. Procedural convergence and a positive-definite Hessian are useful signals, but neither proves that the scientific estimand was recovered.

In the tested NB2 examples, the observable terminal-status detector had sensitivity 0.2375 and specificity 0.9627: it missed 61 of 80 catastrophic truth errors. Treat diagnostics as triage signals, not as certificates that a scientifically plausible estimate is accurate.

See Can I trust this fit? for the first diagnostic workflow and Convergence and start values when the fit is weak, unstable, or boundary-pinned.

When another package is a better choice

Your main task	Start with
A single-response GLMM	glmmTMB
A single-response spatial model	sdmTMB
One- or two-response distributional regression	drmTMB
An established matrix-interface GLLVM workflow using VA, EVA, or Laplace	gllvm
Bayesian multivariate posterior inference	MCMCglmm or brms
A stacked-trait multivariate likelihood with explicit covariance modes	gllvmTMB, within the boundary above

Next, fit the simplest Gaussian GLLVM, use the function map to find an operation, or consult the response-family guide before choosing a non-Gaussian likelihood.