Can I trust this fit? • gllvmTMB

Before interpreting a fitted GLLVM, separate two questions. Did the optimizer find a numerically healthy point estimate? Does the fitted response distribution show an obvious mismatch to the three behavioural responses arranged in long format? This article answers those questions for a small count example before moving on to covariance, ordination, repeatability, or interval claims.

Use these diagnostics to decide what to inspect or refit next, not as a certificate that every inferential target is calibrated.

The example is deliberately modest: three behaviours counted on each individual, fitted with one shared latent variable. latent() estimates a low-rank axis of between-individual covariance; this is a diagnostic example, not the full ordinary latent() covariance-decomposition tutorial from Morphometrics.

Setup

library(gllvmTMB)
library(ggplot2)

Fit A Small Model

The data have one row per (individual, trait) observation. The working question is whether the three behaviours show a shared between-individual axis and whether the fitted Poisson response model looks reasonable enough to inspect.

set.seed(20260609)

n_ind <- 48L
trait_names <- c("forage", "alarm", "explore")
n_traits <- length(trait_names)

z <- rnorm(n_ind, sd = 0.45)
eta <- cbind(
  forage = 0.25 + 0.75 * z,
  alarm = -0.15 - 0.45 * z,
  explore = 0.05 + 0.55 * z
)
y <- matrix(
  rpois(n_ind * n_traits, lambda = exp(as.vector(eta))),
  nrow = n_ind,
  ncol = n_traits
)

df_long <- data.frame(
  individual = factor(rep(seq_len(n_ind), each = n_traits)),
  trait = factor(rep(trait_names, n_ind), levels = trait_names),
  count = as.vector(t(y))
)

head(df_long, 8)
#>   individual   trait count
#> 1          1  forage     0
#> 2          1   alarm     0
#> 3          1 explore     1
#> 4          2  forage     3
#> 5          2   alarm     0
#> 6          2 explore     0
#> 7          3  forage     1
#> 8          3   alarm     2

The long-format call is canonical. It names the response, the trait column, and the unit column explicitly.

fit_long <- gllvmTMB(
  count ~ 0 + trait +
    latent(0 + trait | individual, d = 1),
  data = df_long,
  trait = "trait",
  unit = "individual",
  family = poisson(),
  control = gllvmTMBcontrol(n_init = 2, init_jitter = 0.03)
)
#> 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.

fit <- fit_long
fit
#> Stacked-trait gllvmTMB fit
#>   Traits = 3, individual = 48 
#>   Covstructs: latent_unit, indep_unit 
#>   Fixed effects (b_fix): 3
#>   ML log L = -190.257   convergence = 0
#>   Run gllvmTMB_diagnose(fit) for a full health check, or summary(fit) for parameter estimates.

The same model can be fitted from a wide data frame with the traits(...) left-hand side. The wide call is useful when the data start as one row per individual and one column per behaviour.

df_wide <- stats::reshape(
  df_long,
  idvar = "individual",
  timevar = "trait",
  direction = "wide"
)
names(df_wide) <- sub("^count\\.", "", names(df_wide))

fit_wide <- gllvmTMB(
  traits(forage, alarm, explore) ~ 1 +
    latent(1 | individual, d = 1),
  data = df_wide,
  unit = "individual",
  family = poisson(),
  control = gllvmTMBcontrol(n_init = 2, init_jitter = 0.03)
)

signif(
  abs(as.numeric(logLik(fit_long)) - as.numeric(logLik(fit_wide))),
  3
)
#> [1] 2.14e-09

The observed difference is about 2.14e-09, which is negligible at the reported precision: the two calls reach the same stacked-trait likelihood.

Read Fit Health

The conservative convergence verdict is fit$fit_health$converged. It requires an optimiser success code, a finite objective, and an unscaled maximum gradient below 0.01. scaled_gradient is retained as a descriptive diagnostic, but objective scaling cannot turn a failed optimiser or a large raw gradient into a usable optimum. Hessian curvature and sdreport() remain separate inference checks.

data.frame(
  converged = isTRUE(fit$fit_health$converged),
  optimizer_converged = isTRUE(fit$fit_health$optimizer_converged),
  stationary_by_gradient = isTRUE(fit$fit_health$stationary_by_gradient),
  raw_max_gradient = signif(fit$fit_health$max_gradient, 4),
  raw_gradient_threshold = 0.01,
  scaled_gradient_descriptive = signif(fit$fit_health$scaled_gradient, 4)
)
#>   converged optimizer_converged stationary_by_gradient raw_max_gradient
#> 1      TRUE                TRUE                   TRUE         1.72e-06
#>   raw_gradient_threshold scaled_gradient_descriptive
#> 1                   0.01                   8.995e-09

check_gllvmTMB() returns the detailed, machine-readable fit-health table used in simulations, tests, and reports. gllvmTMB_diagnose() instead prints a guided narrative for an interactive first look. Start with either one before interpreting coefficients, latent variables, or covariance summaries.

health <- check_gllvmTMB(fit)
key_components <- c(
  "optimizer_convergence",
  "max_gradient",
  "sdreport",
  "pd_hessian",
  "weak_axis_unit",
  "boundary_sigma_eps"
)
health[
  health$component %in% key_components,
  c("component", "status", "message")
]
#>                component status
#> 1  optimizer_convergence   PASS
#> 2           max_gradient   PASS
#> 3               sdreport   PASS
#> 4             pd_hessian   PASS
#> 11        weak_axis_unit   PASS
#> 13    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
#> 11              Lambda_B column share of shared loading energy
#> 13     sigma_eps is mapped off by the fitted model/family path

These six displayed rows pass for this example. The full table contains more checks, so always inspect every warning rather than treating a selected subset as the whole verdict. Here are all non-passing rows, including their thresholds, meaning, and suggested action.

health[
  health$status != "PASS",
  c("component", "status", "value", "threshold", "message", "action")
]
#>             component status          value threshold
#> 9      boundary_flags   WARN near_zero_sd_B      none
#> 12 near_zero_psi_unit   WARN       1.03e-07     1e-04
#>                                                 message
#> 9  near-boundary loading or variance component detected
#> 12 sd_B minimum fitted per-trait psi standard deviation
#>                                                                                                   action
#> 9                                            consider lower rank, simpler covariance, or stronger starts
#> 12 check whether the trait-specific component is intentionally mapped off, boundary-pinned, or redundant

In this fit, the two warnings describe the same feature: the fitted trait-specific diagonal standard deviation is almost zero. The simulation planted a shared latent component without an additional trait-specific between-individual component, so a boundary estimate is plausible here. It still matters for interpretation: this example supports a shared axis, but it does not demonstrate recovery of a separate diagonal variance. “Mapped off” means a parameter was deliberately fixed or excluded from the fitted model, rather than estimated near zero.

Read the rows literally. Optimizer and gradient rows speak to the point estimate. sdreport and pd_hessian speak to Hessian-based standard errors and Wald-style inference. A weak-axis row says a latent variable explains only a small share of the fitted common covariance and may therefore be hard to interpret; it does not by itself prove that the fitted response distribution is wrong.

boundary_sigma_eps checks the continuous-family residual scale when that parameter is active. A warning means the estimate is below the reported threshold, so the residual scale is close to its boundary; it does not apply when the model or response family maps that parameter off.

The same table records intentionally skipped standard errors. A no-SE fit can still be useful for point estimates or refit-based uncertainty, such as a bootstrap that repeatedly fits the model to new samples, but it should not be reported as if Wald standard errors were computed.

fit_no_se <- gllvmTMB(
  count ~ 0 + trait +
    latent(0 + trait | individual, d = 1),
  data = df_long,
  trait = "trait",
  unit = "individual",
  family = poisson(),
  control = gllvmTMBcontrol(se = FALSE)
)

no_se_health <- check_gllvmTMB(fit_no_se)
no_se_health[
  no_se_health$component %in% c("sdreport", "pd_hessian", "max_fixed_se"),
  c("component", "status", "message")
]
#>      component status
#> 3     sdreport   WARN
#> 4   pd_hessian   WARN
#> 6 max_fixed_se   WARN
#>                                                             message
#> 3 standard-error calculation skipped by gllvmTMBcontrol(se = FALSE)
#> 4           positive-definite Hessian for curvature-based inference
#> 6                               largest fixed-effect standard error

For hard fits, move from this table to the Convergence and start values article before making scientific claims.

Check Response Diagnostics

Numerical convergence is not enough. For Gaussian, Poisson, and NB2 rows, residuals(type = "randomized_quantile") evaluates each fitted family’s cumulative distribution function at the observed response. These are exact CDF calculations at the fitted parameter values and fitted random-effect modes, but they do not propagate parameter, refit, or random-effect uncertainty and are not a formal calibration test. NB1 uses the same implemented exact-CDF route and has a display smoke test, but it has not yet reached the same depth of evidence as Gaussian, Poisson, and NB2. diagnostic_table() turns the attached metadata into ordinary data frames for reports and articles.

rq <- residuals(
  fit,
  type = "randomized_quantile",
  seed = 20260609
)

diagnostic_table(rq, table = "row_status")
#>   status   n
#> 1     ok 144
diagnostic_table(rq, table = "fit_health_status")
#>   status  n
#> 1   PASS 11
#> 2   WARN  2

head(
  diagnostic_table(rq, table = "data")[
    , c("trait", "family", "observed", "residual", "status")
  ],
  8
)
#>     trait  family observed   residual status
#> 1  forage poisson        0 -1.2801242     ok
#> 2   alarm poisson        0 -0.7317520     ok
#> 3 explore poisson        1  0.4856749     ok
#> 4  forage poisson        3  1.4685245     ok
#> 5   alarm poisson        0 -0.8456795     ok
#> 6 explore poisson        0 -2.1057465     ok
#> 7  forage poisson        1 -0.2994617     ok
#> 8   alarm poisson        2  0.6812327     ok

The row_status table tells you whether every model row produced a usable residual. If a row is marked unsupported_family, missing_response, or nonfinite_observed, keep that status in the report rather than silently dropping the row.

The Q-Q points below stay broadly near the reference line through the centre, but the alarm and forage panels depart in the upper tail. The display is useful because it does not turn a mostly straight pattern into a blanket pass: those tails deserve inspection, and the plot is not evidence that the latent rank or uncertainty estimates are calibrated.

predictive_check(
  fit,
  type = "rq_qq",
  seed = 20260609
)

Randomized-quantile residual Q-Q check for the fitted Poisson GLLVM. This is a fitted-model diagnostic display, not interval calibration, latent-rank proof, or a Bayesian posterior predictive check.

For Poisson and NB2 count models, a rootogram is often easier to read than a Q-Q plot. Bars near zero mean the observed and fitted-model simulated count frequencies are similar on the square-root scale. The simulations default to condition_on_RE = TRUE, so they condition on the fitted random-effect modes; they check the fitted response pattern, not the marginal random-effect distribution. The plot below shows noticeable alternating discrepancies, especially for forage counts 3 and 5. They do not all point in one direction, but they are large enough to investigate rather than dismiss. Simulation noise remains, even with the larger draw count used here.

predictive_check(
  fit,
  type = "rootogram",
  ndraws = 200,
  seed = 20260609,
  max_count = 7L
)

Rootogram comparing observed behavioural counts with simulations conditional on the fitted random-effect modes. Persistent departures can flag a response-family, link, covariate, offset, or fitted-linear-predictor mismatch; this display does not test the marginal random-effect distribution.

What To Do Next

Use the diagnostics as routing information, not as a single pass/fail stamp.

What you see	What it means	Next step
The four numerical rows shown above all pass	The numerical fit is a reasonable starting point, provided the remaining rows contain no unexplained warnings	Inspect every non-PASS row, then inspect `Sigma`, correlations, ordination, or fitted effects.
Optimizer or gradient rows warn/fail	The point estimate may be unstable	Try more starts, rescaling, simpler rank, or an alternative optimizer.
`sdreport` or `pd_hessian` warns	Wald uncertainty is unavailable or local curvature may be unreliable	Inspect the scientific target and use a supported profile or refit-based bootstrap route where available.
Residual rows have non-`ok` status	Some rows were not diagnosed	Inspect missing values, non-finite responses, and family support.
Q-Q or rootogram shows systematic structure	The fitted conditional response pattern may be wrong	Revisit family, link, covariates, offsets, or the fitted linear predictor; investigate random effects separately.

If the fit-health and response checks are acceptable but the latent rank is still a modelling choice, continue to How many latent dimensions should I fit?. That article compares one ML candidate set with AIC/BIC and uses REML only after a Gaussian model has been chosen.

What A Simulation Check Can Establish

Simulation is useful only when its question is named. Three common checks answer different questions:

Check	What it can show	What it cannot show
Fitted-model predictive draw	Selected response patterns	Biological truth, latent rank, or interval calibration
Fitted-model simulate–refit	Local pipeline self-consistency	Real-data truth or protection against misspecification
Predeclared known-DGP grid	Bias, recovery, and interval performance in the studied regime	Performance outside the regimes actually studied

One simulated dataset is an illustration, not recovery evidence. A validation study needs repeated independent datasets, an estimand that matches the fitted model, and separate counts for attempted fits, healthy fits, usable intervals, and successful coverage. It should also report Monte Carlo uncertainty rather than turning a small number of replicates into a pass/fail certificate.

This page covers fit health and no-SE status, fitted-model predictive diagnostics and residuals for Gaussian, Poisson, and NB2 rows, and report-ready metadata tables. NB1 uses an implemented exact-CDF route, but its evidence is currently shallower than for those three families.

The page does not calibrate confidence intervals, prove the latent rank, run formal residual tests, test the marginal random-effect distribution, or draw from a Bayesian posterior. Exact residual support for delta, hurdle, truncated, ordinal, and mixture-family rows remains outside its current claims.