Get started with gllvmTMB

gllvmTMB fits multivariate models for data where each site, individual, species, or study has several responses. The first question is usually not mathematical:

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

This page does one small, reproducible example. 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.

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

The online guide follows the development branch and may describe features absent from this 0.7.1 release. Use this installed vignette and package help for the matching release syntax. The linked online Current limitations and boundaries page is supplementary and may include changes made after 0.7.1. For the scope and restrictions of 0.7.1, use this release’s NEWS and the help and vignettes installed with the package. The Gaussian latent model below is retained as a teaching example for the covariance decomposition and rotation-invariant summaries. Use this fit to learn the decomposition and interpret covariance point estimates. For a first point estimate in the narrow tested Gaussian setting, start with indep() or dep() and check fit health. The example does not certify general interval coverage.

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.

The R formula for that model is compact:

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.

First fit: wide individual-by-trait morphology

We have n individuals, each measured once on T morphological traits. The biological question: do the traits covary on a small number of latent variables (e.g. a body-size axis, a shape axis)? With one observation per individual-trait pair, ordinary latent() estimates both the shared low-rank component and the trait-specific diagonal Psi. No spatial or phylogenetic structure here.

library(gllvmTMB)

For a first fit, load the prepared morphometrics example object. It contains long data, wide data, and the truth table. The reproducible generator is data-raw/examples/make-morphometrics-example.R, but beginners do not need to read the generator before fitting the model.

example_path <- system.file(
  "extdata", "examples", "morphometrics-example.rds",
  package = "gllvmTMB"
)
if (!nzchar(example_path)) {
  stop("Install gllvmTMB to use its morphometrics example data.")
}
morph <- readRDS(example_path)
df <- morph$data_long
df_wide <- morph$data_wide
truth <- morph$truth
trait_names <- truth$trait_names
n_traits <- length(trait_names)
Lambda_true <- truth$Lambda
morph$story$question
#> [1] "Do five body measurements share two latent axes, size and shape, while retaining trait-specific variation?"
head(df_wide, 3)
#>   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
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

The example has shared low-rank covariance \(\boldsymbol\Lambda\boldsymbol\Lambda^\top\) plus per-trait diagonal variance \(\boldsymbol\Psi\). The fit below recovers the individual-level covariance via ordinary latent().

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 a d = 2 reduced-rank GLLVM with ordinary latent() at the individual level: the shared low-rank structure plus trait-specific unique variance, which ordinary latent() carries by default:

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

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

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:

fit_long <- gllvmTMB(
  value ~ 0 + trait + latent(0 + trait | individual, d = 2),
  data = df,
  trait = "trait",
  unit = "individual",
  family = gaussian()
)

# 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, 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

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.

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
#>                     estimand trait_i trait_j i j level component matrix
#> 1  Sigma_unit[length,length]  length  length 1 1  unit     total  Sigma
#> 2    Sigma_unit[length,mass]  length    mass 1 2  unit     total  Sigma
#> 3      Sigma_unit[mass,mass]    mass    mass 2 2  unit     total  Sigma
#> 4    Sigma_unit[length,wing]  length    wing 1 3  unit     total  Sigma
#> 5      Sigma_unit[mass,wing]    mass    wing 2 3  unit     total  Sigma
#> 6      Sigma_unit[wing,wing]    wing    wing 3 3  unit     total  Sigma
#> 7  Sigma_unit[length,tarsus]  length  tarsus 1 4  unit     total  Sigma
#> 8    Sigma_unit[mass,tarsus]    mass  tarsus 2 4  unit     total  Sigma
#> 9    Sigma_unit[wing,tarsus]    wing  tarsus 3 4  unit     total  Sigma
#> 10 Sigma_unit[tarsus,tarsus]  tarsus  tarsus 4 4  unit     total  Sigma
#> 11   Sigma_unit[length,bill]  length    bill 1 5  unit     total  Sigma
#> 12     Sigma_unit[mass,bill]    mass    bill 2 5  unit     total  Sigma
#> 13     Sigma_unit[wing,bill]    wing    bill 3 5  unit     total  Sigma
#> 14   Sigma_unit[tarsus,bill]  tarsus    bill 4 5  unit     total  Sigma
#> 15     Sigma_unit[bill,bill]    bill    bill 5 5  unit     total  Sigma
#>     estimate lower upper interval_method interval_status  scale diagonal
#> 1  1.1813814    NA    NA            none            none latent     TRUE
#> 2  1.1327864    NA    NA            none            none latent    FALSE
#> 3  1.3287406    NA    NA            none            none latent     TRUE
#> 4  0.8660207    NA    NA            none            none latent    FALSE
#> 5  0.8621340    NA    NA            none            none latent    FALSE
#> 6  1.1448231    NA    NA            none            none latent     TRUE
#> 7  0.7652080    NA    NA            none            none latent    FALSE
#> 8  0.7383726    NA    NA            none            none latent    FALSE
#> 9  0.8852321    NA    NA            none            none latent    FALSE
#> 10 0.9995113    NA    NA            none            none latent     TRUE
#> 11 0.6693162    NA    NA            none            none latent    FALSE
#> 12 0.8040512    NA    NA            none            none latent    FALSE
#> 13 0.2984193    NA    NA            none            none latent    FALSE
#> 14 0.1790593    NA    NA            none            none latent    FALSE
#> 15 0.9340536    NA    NA            none            none latent     TRUE
#>    triangle
#> 1  diagonal
#> 2     upper
#> 3  diagonal
#> 4     upper
#> 5     upper
#> 6  diagonal
#> 7     upper
#> 8     upper
#> 9     upper
#> 10 diagonal
#> 11    upper
#> 12    upper
#> 13    upper
#> 14    upper
#> 15 diagonal

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

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.

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\) is the proportion of each trait’s variance shared with the others:

extract_communality(fit, "unit")
#>    length      mass      wing    tarsus      bill 
#> 0.8911057 0.9348019 0.8055344 0.8684222 0.7801664

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, and method = "bootstrap" is a slower target-specific option. Nonlinear profile intervals for correlations are unavailable. Neither supported interval route is universally calibrated.

corr_rows <- extract_correlations(fit, tier = "unit")
corr_rows
#>    tier trait_i trait_j correlation lower upper method interval_status
#> 1     B  length    mass   0.9041341    NA    NA   none            none
#> 2     B  length    wing   0.7446703    NA    NA   none            none
#> 3     B    mass    wing   0.6990134    NA    NA   none            none
#> 4     B  length  tarsus   0.7041912    NA    NA   none            none
#> 5     B    mass  tarsus   0.6407102    NA    NA   none            none
#> 6     B    wing  tarsus   0.8275497    NA    NA   none            none
#> 7     B  length    bill   0.6371627    NA    NA   none            none
#> 8     B    mass    bill   0.7217349    NA    NA   none            none
#> 9     B    wing    bill   0.2885837    NA    NA   none            none
#> 10    B  tarsus    bill   0.1853178    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.

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

The matrix form displays the same point estimates by trait pair. No interval bounds were requested. The plot does not compute or calibrate uncertainty itself.

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.

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

The correlation results above are point estimates only. Profile intervals for correlations are unavailable. The separate total-variance route targets the diagonal of Sigma, not correlation. This example profiles one trait, chosen from the simulated truth table, and prints the interval status with its bounds. With real data, choose the trait to profile from the scientific question and pass its index with trait_idx; the simulated truth is used here only to make this example’s choice reproducible.

profile_trait <- which.max(diag(truth$Sigma))
total_variance_ci <- profile_ci_total_variance(
  fit,
  tier = "unit",
  trait_idx = profile_trait
)
total_variance_ci[, c(
  "trait", "estimate", "lower", "upper", "interval_status"
)]
#>   trait estimate    lower   upper interval_status
#> 1  mass 1.328741 1.079394 1.67447  certified-0.94

certified-0.94 means the fit falls inside one simulated Gaussian regime (unpenalised native-Laplace, n_units >= 150, d <= 2, and a two-sided 95% input). Its measured 0.9467 marginal coverage is a 0.94 floor, not nominal 95% coverage or a guarantee for an individual interval. route-only means the bounds were computed outside that measured regime. The installed help ?profile_ci_total_variance lists the exact target and status rules. These total-variance intervals do not apply to the correlations shown above.

Choose a covariance source and level

unit is the between-unit covariance grouping, called B in older output. unit_obs is the within-unit observation grouping, called W in older output; each unit_obs level must be nested within one unit. These groupings are different from covariance source families:

Component Keyword family What identifies it
Between-unit or within-unit covariance indep(), dep(), latent() Use unit or unit_obs as the grouping factor in the formula.
Animal relatedness animal_indep(), animal_dep(), animal_latent() A pedigree, A, or Ainv supplies relatedness for the animal IDs.
Phylogenetic relatedness phylo_indep(), phylo_dep(), phylo_latent() A tree or vcv matrix supplies relatedness on the cluster axis, usually species.
Spatial fields spatial_indep(), spatial_dep(), spatial_latent() A mesh built from the fitted locations supplies the spatial structure.
Known kernel kernel_indep(), kernel_dep(), kernel_latent() A supplied matrix K supplies relatedness or similarity.

The mode describes covariance among traits: indep() gives separate trait variances with no cross-trait covariance; dep() estimates an unstructured trait covariance; and latent(..., d = K) estimates a rank-K shared component. Ordinary latent() includes diagonal Psi by default. A source-specific *_latent() term is loadings-only by default; use unique = TRUE to add its source-specific diagonal Psi. The older ordinary pair latent(..., unique = FALSE) + unique(...) remains accepted for compatibility. Check the installed help for the selected keyword before combining covariance sources.

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.

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.

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:

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: