Common pitfalls and how to avoid them

Six checks for when a multivariate fit looks wrong. First ask whether the data coding, target summary, formula, or covariance structure is mismatched before treating the result as an optimizer bug. Each pitfall starts with a symptom, gives one concrete check, and ends with a rule of thumb.

Use this page to diagnose data and model mismatches after confirming that the fitted point surface is stable; it is not a substitute for convergence or uncertainty checks.

The examples use the long-format API because it makes each diagnostic target explicit. The same checks apply when the model is written in wide data-frame form with traits(...). See Get started for the side-by-side long and wide calls.

Before diagnosing a pitfall. Check fit$fit_health$converged, fit$fit_health$optimizer_converged, and the raw fit$fit_health$max_gradient, then run check_gllvmTMB(fit). The objective-scaled fit$fit_health$scaled_gradient is a secondary descriptive signal; it cannot override a failed optimiser code or a large raw gradient. A failed Hessian or standard-error calculation warns against Wald inference; it does not by itself prove that the point estimate failed. Conversely, profiles or bootstraps cannot rescue an unstable point surface. See Fit diagnostics and Convergence and start values for the full decision path.

library(gllvmTMB)

1. Keep trait labels and matrix rows in the same order

Symptom. A fitted trait covariance, correlation, or loading matrix has the right labels but the rows do not line up with the reference matrix.

Diagnosis. factor(rep(traits, n)) sorts levels alphabetically by default. If your declared order is boldness, exploration, activity, aggression, sociability, the model sees a different order unless you pin the levels. The fit can be correct while the comparison is wrong.

traits <- covex$truth$trait_names

# MISMATCHED: alphabetical levels
levels(factor(rep(traits, 3)))
#> [1] "activity"    "aggression"  "boldness"    "exploration" "sociability"

# MATCHED: pin levels to the declared order
levels(factor(rep(traits, 3), levels = traits))
#> [1] "boldness"    "exploration" "activity"    "aggression"  "sociability"

Rule of thumb. Whenever you compare model output with a known matrix or a manually ordered display, pin the trait factor with levels = <your trait vector>. Otherwise the comparison can be meaningless even when the fit is right.

2. Compare the estimate with the estimand you asked for

Symptom. A variance or covariance summary looks too small or too large by a familiar amount, such as a residual variance or a diagonal trait-specific component.

Diagnosis. First check what the extractor returns. Here the formula asks for the loadings-only subset with latent(..., unique = FALSE), so extract_Sigma(level = "unit") returns the between-unit covariance $\boldsymbol\Lambda\boldsymbol\Lambda^{\!\top}$ . It does not return the marginal observation covariance

$\operatorname{Var}(\mathbf y_i) = \boldsymbol\Lambda\boldsymbol\Lambda^{\!\top} + \sigma_\varepsilon^2\mathbf I.$

The apparently missing amount can therefore be real observation noise rather than bias in the between-unit covariance.

set.seed(1)
sim <- simulate_site_trait(
  n_sites = 60, n_species = 1, n_traits = 3,
  mean_species_per_site = 1,
  Lambda_B = matrix(c(1.0, 0.4, -0.3,
                      0.0, 0.8,  0.5), 3, 2),
  psi_B = NULL, beta = matrix(0, 3, 2), seed = 1
)
fit <- gllvmTMB(
  value ~ 0 + trait + latent(0 + trait | site, d = 2, unique = FALSE),
  data = sim$data,
  trait = "trait",
  unit = "site"
)

S_hat    <- extract_Sigma(fit, level = "unit")$Sigma
LLt_true <- tcrossprod(sim$truth$Lambda_B)
round(rbind(diag_LLt_true   = diag(LLt_true),
            diag_Sigma_hat = diag(S_hat)), 2)
#>                trait_1 trait_2 trait_3
#> diag_LLt_true     1.00    0.80    0.34
#> diag_Sigma_hat    1.03    0.66    0.58

Rule of thumb. Name the estimand before judging the estimate. If the scientific target includes trait-specific diagonal $\boldsymbol\Psi$ , use the default latent() fit, or add indep(0 + trait | site) to a loadings-only formula, and inspect the total between-unit covariance. Observation residual variance is a separate component; do not silently mix it into that target. Small finite-sample differences are expected even when the estimand matches.

3. Make the data-generating story match the formula

Symptom. A simulation-recovery check shows extra variation in a component that should have been clean.

Diagnosis. simulate_site_trait() generates nonzero environmental coefficients when beta = NULL, whereas psi_B = NULL means that no between-unit diagonal component is simulated. Observation residual variance is generated separately through sigma2_eps. If the simulator contains an environmental or diagonal component that the fitted formula omits, that variation has to be absorbed elsewhere. The optimizer can behave correctly while the fitted model answers a different question from the simulation.

# A clean latent-only data-generating process
sim_clean <- simulate_site_trait(
  n_sites = 50, n_species = 1, n_traits = 3,
  mean_species_per_site = 1,
  Lambda_B = matrix(c(1.0, 0.4, -0.3,
                      0.0, 0.8,  0.5), 3, 2),
  beta = matrix(0, 3, 2),
  psi_B = NULL,
  seed = 1
)

Rule of thumb. For a simulation check, either remove a component from the data-generating process or include its matching term in the formula. For this latent-only example, set beta = matrix(0, n_traits, n_predictors) and leave psi_B = NULL. If the simulation includes environmental effects, include the corresponding trait-specific slopes in the fit.

4. Interpret rotation-invariant summaries, not arbitrary coordinates

Symptom. Two fits on identical data return different loading coordinates.

Diagnosis. Reduced-rank factor models have rotational ambiguity: $\boldsymbol\Lambda$ and $\boldsymbol\Lambda\mathbf Q$ imply the same shared covariance for any orthogonal $\mathbf Q$ . The covariance is rotation-invariant; the loading coordinates are not. Rotation invariance does not guarantee precise estimation, so check fit health and uncertainty separately.

# Sigma is rotation-invariant
round(extract_Sigma(fit, level = "unit")$Sigma, 2)
#>         trait_1 trait_2 trait_3
#> trait_1    1.03    0.19   -0.40
#> trait_2    0.19    0.66    0.44
#> trait_3   -0.40    0.44    0.58

# Varimax provides one reproducible presentation convention
rot <- rotate_loadings(fit, level = "unit", method = "varimax")
round(rot$Lambda, 2)
#>           LV1  LV2
#> trait_1  1.01 0.02
#> trait_2  0.17 0.80
#> trait_3 -0.41 0.64

# Or obtain an identification scaffold for a constrained refit
sug <- suggest_lambda_constraint(
  value ~ 0 + trait + latent(0 + trait | site, d = 2),
  data = sim$data, convention = "lower_triangular"
)
#> Warning: ! Ordinary `latent()` now includes a per-trait Psi by default (Sigma = Lambda
#>   Lambda^T + Psi).
#> ℹ This changed in gllvmTMB 0.2.0; earlier `latent()` was loadings-only (Lambda
#>   Lambda^T).
#> → Pass `latent(..., unique = FALSE)` for the old rotation-invariant
#>   loadings-only fit.
sug$constraint
#>         f1 f2
#> trait_1 NA  0
#> trait_2 NA NA
#> trait_3 NA NA

Rule of thumb. Use extract_Sigma() for rotation-invariant covariance questions. Use rotate_loadings(..., method = "varimax") for one ordered, sign-anchored presentation of the axes, not a uniquely true biological solution. suggest_lambda_constraint() supplies an identification convention; confirmatory claims require constraints specified from the scientific hypothesis rather than chosen from the fitted pattern.

5. Declare the observation-grouping column

Symptom. Error: Column site not found in data.

Diagnosis. Every stacked-trait model needs a column identifying the observational unit. gllvmTMB() defaults to unit = "site", but the unit may instead be an individual, species, or study. A grouping term in the formula does not replace unit =: the argument tells the data-stacking and prediction machinery which rows belong to the same observational unit.

df <- sim$data
df$individual <- df$site
df$site <- NULL

# MISMATCHED (uncomment to see the error):
# gllvmTMB(value ~ 0 + trait + latent(0 + trait | individual, d = 2),
#          data = df, trait = "trait")

# MATCHED: name the observational unit explicitly
fit_ind <- gllvmTMB(
  value ~ 0 + trait + latent(0 + trait | individual, d = 2),
  data = df, trait = "trait", unit = "individual"
)
fit_ind$fit_health$converged
#> [1] TRUE

Rule of thumb. If the observational-unit column is not literally named site, pass unit = "<your_column>". The morphometrics article shows the canonical individual × trait setup.

6. Separate relatedness covariance from sampling-error covariance

Symptom. You have a known matrix but are unsure whether it belongs in a relatedness term or a measurement-error term.

Diagnosis. First name what the matrix indexes and what variation it represents. A, Ainv, and V are not interchangeable.

Symbol	Meaning	Model route	Required alignment
A	Known relatedness covariance among individuals	`A = A` inside an `animal_()` or `phylo_()` term	Symmetric covariance aligned to individual labels
Ainv	Compatible sparse relatedness precision; it may include unobserved pedigree ancestors	`Ainv = Ainv` inside an `animal_()` or `phylo_()` term	Labelled precision aligned to all indexed pedigree nodes
V	Known sampling-error covariance among observation rows	`meta_V(V = V)` plus top-level `known_V = V`	$n_{obs} \times n_{obs}$ covariance aligned to data rows

For independent sampling errors, first assemble the row-aligned covariance matrix, for example V <- diag(sampling_variance). Supply that matrix in the formula as meta_V(V = V) and to gllvmTMB() as known_V = V. A pedigree covariance instead belongs in A = (or its compatible precision in Ainv =) inside an animal_*() term.

Rule of thumb. Use A or Ainv for relatedness and V for sampling-error covariance. Name what each matrix indexes before choosing a keyword. See the Vocabulary glossary for plain-language definitions.