Get started with gllvmTMB • gllvmTMB

GLLVM means generalized linear latent-variable model. It models several responses together and uses a small number of unobserved variables to describe how those responses vary together after accounting for measured predictors. gllvmTMB fits these models in R. Your responses might be species counts at survey sites, morphological traits of individuals, or environmental measurements at sampling locations. Each site or individual is a sampling unit; each species, trait, or measurement is a response.

Which responses vary together, and how much variation is shared rather than response-specific?

This page works through one small, reproducible simulated example with continuous morphological traits. It loads a prepared morphometrics object, fits the model from a wide trait table first, checks the equivalent long-data call, and extracts the covariance, correlation, communality, and ordination summaries a reader normally wants first. The data are included with the package, and the example shows both wide data (one row per individual) and long data (one row per individual and trait). As they are simulated, the object also contains the known generating values used later to check what the fit recovered.

gllvmTMB is experimental: the formula grammar, defaults, and reported summaries and intervals may still change as the API matures.

Use this example to learn how to fit a model and read its covariance summaries. Successful fitting does not establish that the estimates are accurate for your data or that an interval has its stated coverage. The Gaussian latent model taught here has not passed all of its parameter-recovery checks. For continuous traits, simpler models with independent trait variation (indep()) or a full trait covariance (dep()) have stronger point-estimation evidence in the tested three-trait, complete-data examples.

Before adapting the example to counts, presence/absence, missing observations, or a spatial or phylogenetic model, read Current limitations and boundaries. It explains which results have been checked, where estimates can fail even when the fit appears healthy, and why uncertainty needs a separate check.

First fit: wide individual-by-trait morphology

We have 150 individuals, each measured once on five morphological traits. The biological question: do the traits covary on a small number of latent axes (for example a body-size axis and a shape axis)? With one observation per individual and trait, ordinary latent() estimates both parts: the shared part carried by a few latent axes, and the leftover variance specific to each trait (Psi in the equations later on this page). No spatial or phylogenetic structure here.

library(gllvmTMB)

For a first fit, load the prepared morphometrics example object. It contains simulated long and wide data, the known generating values (truth), and the two formulas used below. The generating code is available in the source repository for reproducibility, but beginners do not need to read it before fitting the model.

morph <- readRDS(system.file(
  "extdata", "examples", "morphometrics-example.rds",
  package = "gllvmTMB"
))
df_wide <- morph$data_wide

morph$story$question
#> [1] "Do five body measurements share two latent axes, size and shape, while retaining trait-specific variation?"
head(df_wide)
#>   individual      length       mass       wing     tarsus        bill
#> 1          1 -0.04447039 -1.2303106  0.2349369  0.6489769 -0.55819209
#> 2          2  0.06037967  0.4328307  0.4848265  0.3284936 -0.02838954
#> 3          3  2.26418385  3.2667958  0.8109460  1.3900777  2.06870978
#> 4          4  2.30693568  2.5671014  1.9187735  1.2025903  0.95885641
#> 5          5  1.29077809  1.6774179  0.8869494  1.2821702  0.89538354
#> 6          6 -0.70460367 -0.5803515 -0.2838886 -0.1574360  0.18059406

The example was generated with a small number of shared axes plus variation specific to each trait. The fit below recovers the individual-level covariance via ordinary latent(); the equations are in “What the fitted model means”, after the covariance table.

In wide data, the formula term is latent(1 | individual, d = 2). In long data, the equivalent term is latent(0 + trait | individual, d = 2). The 1 in the wide form is the compact shorthand for the same trait-specific columns used in the long form.

Fit the model with two latent axes (d = 2) using ordinary latent() at the individual level. That term estimates the shared structure the axes carry plus the trait-specific variance, which ordinary latent() includes by default. In an interactive session (this page hides messages when knitted), the first ordinary latent() fit prints a one-time warning that latent() includes per-trait Psi by default, and the fit may note that sigma_eps was fixed because an internal indep(0 + trait | individual) term already absorbs observation-level residual variance on the stacked rows:

morph$formula_wide
#> traits(length, mass, wing, tarsus, bill) ~ 1 + latent(1 | individual, 
#>     d = 2)

fit <- gllvmTMB(
  morph$formula_wide,
  data = df_wide,
  unit = morph$fit_args$unit,
  family = morph$fit_args$family
)

# A finite likelihood confirms that the model fitted; diagnostics come next.
as.numeric(logLik(fit))
#> [1] -739.6107

That single number says the fit returned a finite log-likelihood, so there is a fitted model to inspect. It does not say the optimiser converged or that the model is right; run the checks under “Check the fitted response distribution” before interpreting anything.

Same model, long data

If your workflow already stores one row per (individual, trait) observation, use the same entry point with the long formula and explicit trait = argument. Internally the wide formula pivots to this stacked form, so the two calls fit the same model — the log-likelihoods agree:

morph$formula_long
#> value ~ 0 + trait + latent(0 + trait | individual, d = 2)

df <- morph$data_long
head(df, 6)
#>   individual  trait       value
#> 1          1 length -0.04447039
#> 2          1   mass -1.23031062
#> 3          1   wing  0.23493689
#> 4          1 tarsus  0.64897692
#> 5          1   bill -0.55819209
#> 6          2 length  0.06037967

fit_long <- gllvmTMB(
  morph$formula_long,
  data = df,
  trait = morph$fit_args$trait,
  unit = morph$fit_args$unit,
  family = morph$fit_args$family
)

# Same fit:
all.equal(as.numeric(logLik(fit)), as.numeric(logLik(fit_long)))
#> [1] TRUE

What if some response cells are missing?

Wide trait tables do not need to be perfectly rectangular. If an individual is missing one trait measurement, leave that response cell as NA; gllvmTMB() can treat only that unit-trait cell as missing and keep the other observed traits for the individual. The same rule applies to long data with NA in the response column. Missing predictors and missing grouping variables are a separate problem: by default they fail loudly because the model cannot construct that row.

Check the fitted response distribution

Before interpreting Sigma (the fitted trait covariance), separate two questions. First, ask whether the numerical fit looks healthy. Then ask whether the fitted response distribution is visibly mismatched to the observed stacked traits. For this Gaussian example, use check_gllvmTMB() first; if a row is not PASS, read its action before interpreting covariance. The residual and predictive-check plots are diagnostic displays here. They do not prove the latent rank or calibrate uncertainty intervals.

health <- check_gllvmTMB(fit)
health[, c("component", "status", "message", "action")]
#>                       component status
#> 1         optimizer_convergence   PASS
#> 2                  max_gradient   PASS
#> 3                      sdreport   PASS
#> 4                    pd_hessian   PASS
#> 5                  hessian_rank   PASS
#> 6                  max_fixed_se   PASS
#> 7               restart_history   PASS
#> 8              selected_restart   PASS
#> 9                boundary_flags   PASS
#> 10     rotation_convention_unit   WARN
#> 11               weak_axis_unit   PASS
#> 12 cross_loading_structure_unit   PASS
#> 13           near_zero_psi_unit   PASS
#> 14           boundary_sigma_eps   PASS
#>                                                        message
#> 1                               optimizer reported convergence
#> 2  largest absolute gradient component at the selected optimum
#> 3                                           sdreport available
#> 4      positive-definite Hessian for curvature-based inference
#> 5  rank of the fixed-parameter covariance matrix from sdreport
#> 6                          largest fixed-effect standard error
#> 7               number of optimizer starts recorded on the fit
#> 8                        restart selected by minimum objective
#> 9                            no simple boundary flags detected
#> 10       Lambda_B is identified up to rotation/sign convention
#> 11              Lambda_B column share of shared loading energy
#> 12      median trait share carried by its dominant latent axis
#> 13        sd_B minimum fitted per-trait psi standard deviation
#> 14     sigma_eps is mapped off by the fitted model/family path
#>                                                                                                            action
#> 1                                    try multiple starts, stronger starts, rescaling, or an alternative optimizer
#> 2                                            tighten optimization, rescale predictors, or inspect weak components
#> 3                                           use point summaries cautiously and prefer profile/bootstrap intervals
#> 4                                check gradients, boundary variances, rank, starts, and profile/bootstrap targets
#> 5                                                            treat rank loss as a Hessian/identifiability warning
#> 6                                                 check collinearity, scaling, or weakly identified fixed effects
#> 7                                                            refit with current gllvmTMB if provenance is missing
#> 8                                                         inspect restart_history for competing likelihood basins
#> 9                                             still inspect profile/bootstrap output for target-specific weakness
#> 10 use Sigma/correlations/communality for invariant summaries; rotate or constrain loadings before comparing axes
#> 11                              compare lower ranks, inspect fit stability, and avoid over-interpreting weak axes
#> 12                              use varimax/promax rotation for interpretation if loadings are spread across axes
#> 13          check whether the trait-specific component is intentionally mapped off, boundary-pinned, or redundant
#> 14                         if estimated near zero, check row-level unique terms or residual-scale identifiability

The WARN on rotation_convention_unit is expected for every latent fit: Lambda_B names the unit-tier loadings, which are identified only up to rotation and sign until you rotate or constrain them. Its action says to read Sigma, correlations, and communality rather than raw loadings, which is the order this page follows.

Each individual-trait cell has one observation, and the fitted diagonal Psi term is indexed at that same cell resolution. Exact conditional Gaussian residuals would therefore be almost zero and would make a misleadingly perfect Q-Q plot. Instead, use simulation-rank residuals from new fitted-model draws with random effects regenerated (condition_on_RE = FALSE).

rq <- residuals(
  fit,
  type = "simulation_rank",
  nsim = 199,
  seed = 1,
  condition_on_RE = FALSE
)
head(rq[, c("trait", "family", "observed", "residual", "status")], 8)
#>    trait   family    observed    residual status
#> 1 length gaussian -0.04447039  0.06502866     ok
#> 2   mass gaussian -1.23031062 -0.92105112     ok
#> 3   wing gaussian  0.23493689  0.29988256     ok
#> 4 tarsus gaussian  0.64897692  0.80934021     ok
#> 5   bill gaussian -0.55819209 -0.69120234     ok
#> 6 length gaussian  0.06037967  0.07455944     ok
#> 7   mass gaussian  0.43283075  0.45555722     ok
#> 8   wing gaussian  0.48482650  0.56560032     ok
stopifnot(diff(range(rq$residual, na.rm = TRUE)) > 1)

Use the plot as a quick fitted-response check before turning to covariance interpretation:

predictive_check(
  fit,
  type = "rq_qq",
  residual_type = "simulation_rank",
  nsim = 199,
  seed = 1,
  condition_on_RE = FALSE
)

Marginal simulation-rank residual Q-Q check for the fitted Gaussian morphometrics model. New random effects are generated for each fitted-model draw. The plot is a diagnostic display, not interval calibration or a posterior predictive check.

Read the covariance summaries first

Start with the fitted Sigma rows. These rows describe the covariance surface implied by the model, not a separate model:

sigma_rows <- extract_Sigma_table(fit, level = "unit")
sigma_rows[, c("trait_i", "trait_j", "estimate", "diagonal")]
#>    trait_i trait_j  estimate diagonal
#> 1   length  length 1.1813835     TRUE
#> 2   length    mass 1.1327882    FALSE
#> 3     mass    mass 1.3287415     TRUE
#> 4   length    wing 0.8660217    FALSE
#> 5     mass    wing 0.8621357    FALSE
#> 6     wing    wing 1.1448261     TRUE
#> 7   length  tarsus 0.7652105    FALSE
#> 8     mass  tarsus 0.7383751    FALSE
#> 9     wing  tarsus 0.8852361    FALSE
#> 10  tarsus  tarsus 0.9995147     TRUE
#> 11  length    bill 0.6693165    FALSE
#> 12    mass    bill 0.8040513    FALSE
#> 13    wing    bill 0.2984215    FALSE
#> 14  tarsus    bill 0.1790603    FALSE
#> 15    bill    bill 0.9340540     TRUE

Each row is one entry of the fitted trait covariance at the individual level; the diagonal column marks the rows that are variances. The diagonal entries combine shared and trait-specific variance; the off-diagonal entries are the covariance the shared axes create between two traits. The full table (sigma_rows) also carries interval columns, which stay empty until you ask for a specific interval route.

What the fitted model means

The fitted Gaussian model starts with observations, not with the covariance matrix. For each individual and trait, the teaching model is:

$\begin{aligned} y_{it} &= \alpha_t + \boldsymbol{\lambda}_t^{\mathsf T}\mathbf{u}_i + e_{it},\\ \mathbf{u}_i &\sim \mathcal{N}(\mathbf{0}, \mathbf{I}_d),\\ e_{it} &\sim \mathcal{N}(0, \psi_t). \end{aligned}$

Here $y_{it}$ is the observed value for individual $i$ and trait $t$ . The intercept $\alpha_t$ is the average for trait $t$ . The vector $\mathbf{u}_i$ is the individual’s latent score: where that individual sits on the unobserved trait axes. The row vector $\boldsymbol{\lambda}_t$ says how strongly trait $t$ responds to those axes. The leftover term $e_{it}$ is the part of trait $t$ not shared with the other traits.

In generalized models, $\eta_{it}$ is the linear predictor and the link function $g()$ connects it to the mean response: $g(\mu_{it}) = \eta_{it}$ . This first example uses a Gaussian family with the identity link, so the linear predictor is the mean itself.

Written in full, the call you ran above is:

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

After fitting, the covariance among traits is the first report-ready summary implied by the latent part of the model:

$\begin{aligned} \operatorname{Cov}(\mathbf{y}_i) &= \boldsymbol{\Sigma},\\ \boldsymbol{\Sigma} &= \boldsymbol{\Lambda}\boldsymbol{\Lambda}^{\mathsf T} + \boldsymbol{\Psi},\\ \boldsymbol{\Psi} &= \operatorname{diag}(\psi_1,\ldots,\psi_T). \end{aligned}$

In words: total trait covariance = shared multivariate structure + response-specific variation. Lambda collects the trait loadings: one row per trait and one column per latent variable. The product Lambda Lambda^T is the shared covariance created by those axes. Psi is the diagonal matrix of trait-specific variance left over after the shared axes. Ordinary latent() estimates both pieces by default.

For a first fit, read Sigma, correlations, and communality before reading individual Lambda entries. The latent variables can be flipped or turned without changing the covariance they imply, so raw loading columns need a rotation or confirmatory constraint before they become a biological story.

The full covariance target is stored in the example object, so this teaching example can show what the fitted model is trying to recover:

truth <- morph$truth
round(truth$Sigma, 2)
#>        length mass wing tarsus bill
#> length   0.97 0.92 0.68   0.58 0.60
#> mass     0.92 1.14 0.72   0.60 0.72
#> wing     0.68 0.72 1.00   0.76 0.24
#> tarsus   0.58 0.60 0.76   0.86 0.12
#> bill     0.60 0.72 0.24   0.12 0.90

Ordination of individuals on the two latent axes. These axes come from the fitted latent model ( $\boldsymbol\Sigma = \boldsymbol\Lambda\boldsymbol\Lambda^\top + \boldsymbol\Psi$ ), not an eigen-decomposition of the raw data — this is model-based ordination, not PCA or NMDS:

ord <- extract_ordination(fit, level = "unit")
plot(ord$scores[, 1], ord$scores[, 2], pch = 19, col = "steelblue",
     xlab = "LV 1", ylab = "LV 2",
     main = "Individuals on latent trait axes")
abline(h = 0, v = 0, lty = 2, col = "grey70")

Individual scores on the two latent trait axes.

That’s the whole model — latent() for shared trait covariance and the diagonal trait-specific $\boldsymbol\Psi$ component, plus trait-specific intercepts.

Per-trait communality $c_t^2$ — the proportion of each trait’s variance shared with the others:

extract_communality(fit, "unit")
#>    length      mass      wing    tarsus      bill 
#> 0.8911051 0.9348032 0.8055330 0.8684273 0.7801615

Pairwise trait correlations are point estimates by default. This is the safest general summary because interval calibration depends on the fitted target and model. method = "fisher-z" adds fast heuristic bounds, while profile and bootstrap routes are slower and remain target-specific rather than universally calibrated.

corr_rows <- extract_correlations(fit, tier = "unit")
corr_rows
#>    tier trait_i trait_j correlation lower upper method interval_status
#> 1  unit  length    mass   0.9041345    NA    NA   none            none
#> 2  unit  length    wing   0.7446695    NA    NA   none            none
#> 3  unit    mass    wing   0.6990136    NA    NA   none            none
#> 4  unit  length  tarsus   0.7041917    NA    NA   none            none
#> 5  unit    mass  tarsus   0.6407110    NA    NA   none            none
#> 6  unit    wing  tarsus   0.8275510    NA    NA   none            none
#> 7  unit  length    bill   0.6371623    NA    NA   none            none
#> 8  unit    mass    bill   0.7217346    NA    NA   none            none
#> 9  unit    wing    bill   0.2885853    NA    NA   none            none
#> 10 unit  tarsus    bill   0.1853184    NA    NA   none            none

For reports, plot the same tidy rows instead of hand-indexing a correlation matrix:

plot_correlations(corr_rows, sort = "magnitude")

Pairwise point estimates of trait correlations from the first fitted model. No interval bounds are requested in this first-fit view.

The matrix form is another display of the same point estimates. This first-fit object does not contain interval bounds. If you explicitly request Fisher-z bounds, treat them as a fast heuristic rather than a generally calibrated interval; the plot itself does not compute uncertainty.

plot_correlations(
  corr_rows,
  style = "heatmap",
  matrix_layout = "by_level",
  label_type = "estimate",
  title = "First-model trait correlations"
)

Point estimates of trait correlations from the first fitted model. The matrix is descriptive; interval routes are target-specific and opt-in.

Then inspect loadings and scores

The recovered loadings are useful, but they are not the first interpretation target. Only $\boldsymbol{\Lambda}\boldsymbol{\Lambda}^{\mathsf T}$ is identified directly from the data; raw loading columns can change sign or rotate.

n_traits <- length(truth$trait_names)
Lambda_hat <- extract_loadings(fit, level = "unit")
rownames(Lambda_hat) <- paste0("trait_", seq_len(n_traits))
colnames(Lambda_hat) <- paste0("LV", seq_len(ncol(Lambda_hat)))
round(Lambda_hat, 3)
#>           LV1    LV2
#> trait_1 1.026  0.000
#> trait_2 1.104 -0.152
#> trait_3 0.844  0.458
#> trait_4 0.746  0.558
#> trait_5 0.652 -0.551

Do not over-interpret the signs of these unrotated loadings. Use them as a quick look at the latent variables, then use Sigma, correlations, and communality for report-ready interpretation unless you have chosen a rotation or confirmatory constraint.

The shared covariance from the loadings can still be compared with the known truth in this simulated teaching object. The off-diagonal pattern is what should agree; the full Sigma diagonal also includes Psi.

Lambda_true <- truth$Lambda
list(true_LL  = round(Lambda_true %*% t(Lambda_true), 2),
     fitted_LL = round(Lambda_hat  %*% t(Lambda_hat),  2))
#> $true_LL
#>        length mass wing tarsus bill
#> length   0.82 0.92 0.68   0.58 0.60
#> mass     0.92 1.04 0.72   0.60 0.72
#> wing     0.68 0.72 0.80   0.76 0.24
#> tarsus   0.58 0.60 0.76   0.74 0.12
#> bill     0.60 0.72 0.24   0.12 0.72
#> 
#> $fitted_LL
#>         trait_1 trait_2 trait_3 trait_4 trait_5
#> trait_1    1.05    1.13    0.87    0.77    0.67
#> trait_2    1.13    1.24    0.86    0.74    0.80
#> trait_3    0.87    0.86    0.92    0.89    0.30
#> trait_4    0.77    0.74    0.89    0.87    0.18
#> trait_5    0.67    0.80    0.30    0.18    0.73

That is the whole first model: trait-specific intercepts, individual latent scores, loadings from those scores to traits, and the diagonal Psi carried by ordinary latent().

Choose your next guide

The next step depends on the question, not on adding every available model component:

Interpret covariance: continue with Individual morphometrics for a fuller Gaussian worked example.
Check the fit: use Can I trust this fit? before relying on coefficients, covariance, or uncertainty.
Choose a specialised model: use the Articles chooser to select repeated behaviours, reaction norms, binary occurrence, multivariate spatial structure, phylogenetic covariance, missing data, rank selection, or a technical reference. The spatial-model guide shows mesh construction and equivalent long- and wide-format fits.

Family, grouping structure, covariance source, covariance mode, missingness, and latent rank answer different scientific questions. Add only the pieces the sampling design and estimand require.

Acknowledgements

gllvmTMB is a standalone TMB-based package. Its design acknowledges several earlier projects; incorporated code is recorded precisely in inst/COPYRIGHTS:

sdmTMB (Anderson, Ward, English, Barnett, Thorson, et al.): the original spatial interface and helper layer were derived from its GPL-3 implementation. The current gllvmTMB helpers were substantially rewritten against the public fmesher API; sdmTMB remains a credited source rather than a runtime dependency. The spatial-model guide explains this provenance boundary (the helper implementation itself is covered by focused tests; broader spatial-family evidence remains partial).
glmmTMB (Brooks, Bolker, Kristensen, Magnusson, Skaug, Nielsen, Berg, van Bentham, et al.): rr / propto / equalto / diag covariance dispatch. The reduced-rank rr() machinery in particular is the contribution of McGillycuddy, Popovic, Bolker & Warton (2025) J. Stat. Softw. 112(1).
GALAMM (Sørensen): mixed-response (gfam) and confirmatory fixed-loading (lambda_constraint) patterns.
MCMCglmm (Hadfield): the sparse $A^{-1}$ phylogeny representation that phylo_latent(tree = ) follows. gllvmTMB constructs it natively (ape + Matrix only) and does not call MCMCglmm at fit time; it is a test-only Suggests, used to cross-check the pedigree precision builder.
TMB (Kristensen et al.) and fmesher (Lindgren): autodiff and finite-element mesh.