Convergence and hard fits • gllvmTMB

Some multivariate latent models are hard because the likelihood has several plausible basins—a basin is a local solution reached from one set of starting values—weakly identified variance splits, or a flat curvature direction near a boundary. The practical question is therefore not just “did the optimizer print convergence code 0?” but “which part of the fit is healthy, which part is weak, and what should I try next?”

This article covers check_gllvmTMB(), restart history, protected or skipped TMB::sdreport() status, and gllvmTMBcontrol(se = FALSE). Bootstrap and profile uncertainty are follow-up paths when Wald standard errors are unsafe, but they have narrower evidence: bootstrap summaries depend on successful refits, and profile intervals apply only to supported targets with stable profiles. The start strategies below are options to try, not established defaults or universal convergence recommendations.

The most important distinction is simple. fit$fit_health$converged is the package’s conservative point-stationarity conjunction: optimiser success, a finite objective, and a small raw maximum gradient. By contrast, pd_hessian reports TMB’s curvature result. FALSE means ordinary Wald uncertainty is not supported; NA after se = FALSE means the Hessian was not computed. Keep the point-fit and inference checks separate, and do not use objective scaling to override a failed optimiser code.

library(gllvmTMB)
set.seed(20260519)

Fit a small model

The example below uses a two-trait Gaussian model. The scientific target is the trait covariance Sigma = Lambda Lambda^T + Psi, where Lambda contains the shared latent loadings and Psi is the diagonal trait-specific variance. In current ordinary latent() models, that diagonal Psi companion is included by default; latent(..., unique = FALSE) requests the older no-Psi subset.

Use this alignment while reading the diagnostics:

Symbol	R syntax or object	Plain meaning
`Lambda`	`latent(0 + trait \| site, d = 1)`; `Lambda_B` in the simulator	The shared loading axis that creates cross-trait covariance.
`Psi`	default companion inside `latent(0 + trait \| site, d = 1)`; `psi_B` in the simulator	Trait-specific variance that is not shared with the other trait.
`Sigma`	`extract_Sigma_table(fit, level = "unit", part = "total")`	The full trait covariance, `Lambda Lambda^T + Psi`, used for interpretation.

sim <- simulate_site_trait(
  n_sites = 24, n_species = 1, n_traits = 2,
  mean_species_per_site = 1,
  Lambda_B = matrix(c(0.7, 0.4), nrow = 2, ncol = 1),
  psi_B = c(0.3, 0.3),
  seed = 20260519
)
df_long <- sim$data

The long-format call is the canonical interface: one row per (unit, trait) observation, with the trait and unit columns named explicitly.

fit <- gllvmTMB(
  value ~ 0 + trait +
    latent(0 + trait | site, d = 1),
  data = df_long,
  trait = "trait",
  unit = "site",
  control = gllvmTMBcontrol(n_init = 2, init_jitter = 0.05)
)
#> 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.

The same model can be written with the wide data-frame formula interface. This route is useful when your data start as one row per unit and one column per trait.

df_wide <- stats::reshape(
  df_long[c("site", "trait", "value")],
  idvar = "site", timevar = "trait", direction = "wide"
)
names(df_wide) <- sub("^value\\.", "", names(df_wide))

fit_wide <- gllvmTMB(
  traits(trait_1, trait_2) ~ 1 +
    latent(1 | site, d = 1),
  data = df_wide,
  unit = "site",
  control = gllvmTMBcontrol(n_init = 2, init_jitter = 0.05)
)

signif(abs(as.numeric(logLik(fit)) - as.numeric(logLik(fit_wide))), 3)
#> [1] 0

The printed difference is exactly zero at the shown precision, so the long and wide calls reached the same likelihood for this example. That agreement does not replace the health checks below.

Read fit health

Start with the conservative point-stationarity verdict. It requires an optimiser success code, a finite objective, and an unscaled maximum gradient below 0.01. The objective-scaled gradient is retained as a descriptive signal, not a substitute for those 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),
  objective = signif(fit$fit_health$objective, 6),
  optimizer_code = fit$fit_health$convergence
)
#>   converged optimizer_converged stationary_by_gradient raw_max_gradient
#> 1      TRUE                TRUE                   TRUE        1.037e-06
#>   raw_gradient_threshold scaled_gradient_descriptive objective optimizer_code
#> 1                   0.01                   1.608e-08   63.4643              0

health <- check_gllvmTMB(fit)
triage_components <- c(
  "optimizer_convergence", "max_gradient", "sdreport", "pd_hessian",
  "hessian_rank", "restart_history", "boundary_flags", "weak_axis_unit",
  "near_zero_psi_unit"
)
health[
  health$component %in% triage_components,
  c("component", "status", "value")
]
#>                component status           value
#> 1  optimizer_convergence   PASS               0
#> 2           max_gradient   PASS       1.037e-06
#> 3               sdreport   PASS            TRUE
#> 4             pd_hessian   WARN           FALSE
#> 5           hessian_rank   PASS             6/6
#> 7        restart_history   PASS               2
#> 9         boundary_flags   PASS            none
#> 11        weak_axis_unit   PASS min=1; shares=1
#> 12    near_zero_psi_unit   PASS          0.7266

Read the components in three lanes:

Lane	Components	What it answers	Typical next action
Convergence	optimizer code, raw maximum gradient, restart history; scaled gradient is descriptive	Is the selected point sufficiently stationary, and did starts find competing objective values?	Try more starts, rescale predictors, or compare optimizers.
Hessian-based inference	`sdreport`, `pd_hessian`, Hessian rank, fixed-effect SE	Is local curvature usable for Wald uncertainty?	If unavailable or weak, inspect the target and use a supported profile or bootstrap route.
Identifiability and interpretation	weak axes, boundary flags, near-zero `Psi`, rotation convention	Is a covariance component or latent axis weak, redundant, or convention-dependent?	Compare lower rank, simpler covariance, and rotation-invariant `Sigma`; use the advanced identifiability check when appropriate.

In this baseline fit, the conservative convergence conjunction is TRUE, both restarts reach the same displayed objective, sdreport is available, and the reported Hessian rank is full. However, pd_hessian is FALSE. That combination says the selected point is stationary but local curvature still warns against a blanket Wald-uncertainty claim. It does not by itself prove either model failure or a benign rotation ridge; inspect the target and its profile or refit behaviour.

This separation matters. A fit can be stationary while Wald uncertainty is unavailable, or it can have a positive-definite Hessian while an extra latent axis remains scientifically weak. Do not flatten every warning into model failure.

After the numerical triage, inspect the fitted response distribution before leaning on covariance interpretation. residuals(fit, type = "randomized_quantile") and predictive_check(fit, type = "rq_qq") are available for the scoped Gaussian, Poisson, and NB2 diagnostic paths. They serve the same role as the fit-health table: useful for surfacing mismatch, not proof of interval calibration, latent-rank adequacy, or Bayesian posterior prediction. The worked diagnostic sequence is shown in Can I trust this fit?.

Convergence is not identifiability

Convergence asks whether optimization reached a sufficiently stationary point. Identifiability asks whether the data can distinguish the parameters or latent structure being interpreted. Raw Lambda loadings also depend on a rotation convention, whereas Sigma is rotation invariant and is therefore the safer cross-start comparison.

Simulation can assess a predeclared rank-selection rule only when the data- generating rank is fixed independently of the fitted model. Simulating from a rank-two fit and successfully refitting rank two tests stability under that fitted scenario; it cannot establish that rank two was required by the original data. For a study-level validation, predeclare candidate ranks, simulate from known rank-one and rank-two models over realistic sample sizes, and report how often the full selection workflow recovers each generating rank.

Choose a start strategy

Start values are opt-in tools. The actual default is one nlminb fit using historical starts with se = TRUE. The strategies below are experimental: their recovery benefit across families and ranks is not yet established by simulation, so treat the table as a menu of things to try on a hard fit, not a proven policy.

What you observe	Strategy	Control
No warning or competing objective yet	Use the default first	`gllvmTMBcontrol()`
Starts reach different objective values	Add restarts and modest jitter	`gllvmTMBcontrol(n_init = 5, init_jitter = 0.1)`
A non-Gaussian reduced-rank fit repeatedly stops poorly	Seed loadings and scores from residual structure	`gllvmTMBcontrol(start_method = list(method = "res", jitter.sd = 0.2))`
A Gaussian two-level latent fit is sensitive to its initial covariance	Warm-start from the matching independent model	`gllvmTMBcontrol(start_method = list(method = "indep"))`
A scientifically sensible simpler fit is already available	Copy compatible starting parameters	`gllvmTMBcontrol(start_from = simpler_fit)`
A fitted `phi` parameter is unstable, especially for NB1/NB2	Warm up dispersion trait by trait	`gllvmTMBcontrol(init_strategy = "single_trait_warmup")`
`nlminb` does not reach a stationary point	Compare `optim` with BFGS	`gllvmTMBcontrol(optimizer = "optim", optArgs = list(method = "BFGS"))`
Point estimates are stable but Wald calculation is the fragile step	Deliberately skip standard errors	`gllvmTMBcontrol(se = FALSE)`

For a difficult model, start with more than one restart and inspect fit$restart_history. This is currently a low-level fitted-object diagnostic field, so its columns may evolve during the experimental lifecycle.

fit$restart_history[, c("restart", "start_method", "objective",
                        "convergence", "selected")]
#>   restart start_method objective convergence selected
#> 1       1      default  63.46429           0     TRUE
#> 2       2      default  63.46429           0    FALSE

Lower objective values are preferred. Here the two displayed objectives are equal at the shown precision and restart 1 is selected. That supports, but does not prove, that the starts reached the same solution. Do not privilege a start simply because it converged: compare the selected fit’s optimiser code and raw gradient and rotation-invariant targets such as Sigma. Compare raw Lambda only after alignment or rotation handling.

Use no-SE fits honestly

gllvmTMBcontrol(se = FALSE) intentionally skips TMB::sdreport(). Use it only after an ordinary fit has shown that point estimates are stable but standard-error calculation is the fragile step. The skipped status is NA, not evidence that a computed Hessian was non-positive-definite.

fit_no_se <- gllvmTMB(
  value ~ 0 + trait +
    latent(0 + trait | site, d = 1),
  data = df_long,
  trait = "trait",
  unit = "site",
  control = gllvmTMBcontrol(se = FALSE, n_init = 2, init_jitter = 0.05)
)

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

This option is not a way to hide convergence problems. It is a way to avoid making Hessian-based uncertainty claims when the Hessian path is the fragile part of the workflow.

Bootstrap after a one-core smoke test

Bootstrap refits the model many times. It is one practical uncertainty path when pdHess is false, when sdreport() fails, or when the target is nonlinear, but it is not a cure for an unstable fitted surface. bootstrap_Sigma() exposes n_cores so independent refits can run in parallel.

The Gaussian bootstrap_Sigma() route has direct simulate-refit contract and plumbing tests, but this article does not claim uniform interval coverage across all targets. Non-Gaussian and mixed-family refit plumbing is shallower; treat those intervals as an audit trail until the relevant family and target have coverage evidence.

boot_smoke <- bootstrap_Sigma(
  fit,
  n_boot = 25,
  level = "unit",
  what = c("Sigma", "R", "communality"),
  seed = 1,
  n_cores = 1
)

boot_smoke$n_failed
extract_Sigma_table(boot_smoke, level = "unit", entries = "upper")

The 25-replicate call is only a plumbing smoke test. Once the model refits cleanly, increase to at least 200 replicates and optionally use multiple cores. Even 200 is illustrative for percentile endpoints because only a handful of draws determine each tail; increase n_boot further for final reporting and check that endpoints are stable. Always record the failure count. A bootstrap CI with many failed refits is a diagnostic result, not just an inconvenience; report the failure count and simplify the model before treating the interval as a stable uncertainty summary.

Profile when the target is scalar

Profile likelihood is a targeted route for scalar quantities such as repeatability, communality, selected correlations, or a single variance share. Use it when the target has a stable profile and the model is not bouncing between likelihood basins. It is not a general substitute for bootstrap across every derived table entry.

For full latent covariance models, some derived targets still fall back to Wald or bootstrap while constrained-profile machinery is expanded. The output method column should be read literally: if it says "wald", the interval is Hessian-based even if the requested method was "profile".

The generic confint(method = "profile") and confint(method = "bootstrap") routes are contract-tested, but route availability is not the same as empirical interval calibration. Target-specific profile helpers are narrower: check profile_targets(fit) before treating a derived quantity as profile-ready, and read the returned method and status literally. For the interval workflow, continue to Profile-likelihood confidence intervals; for failed or irregular profiles, use its troubleshooting decision table.

Troubleshooting ladder

Use this order when a hard model does not give clean inference:

Read check_gllvmTMB(fit) and fit$restart_history.
Check whether the optimizer, gradient, pdHess, or sdreport component is the first weak link.
Use residuals() or predictive_check() to check fitted-response mismatch before interpreting covariance summaries.
Inspect boundary flags for near-zero variance, redundant rank, or unstable dispersion.
Try more starts and modest jitter.
Try residual starts for non-Gaussian reduced-rank models.
Try independent or simpler-model starts for Gaussian two-level models.
Try optimizer = "optim" with optArgs = list(method = "BFGS"), or a longer nlminb budget such as gllvmTMBcontrol(optArgs = list(control = list(eval.max = 5000, iter.max = 4000))).
If point estimates are stable but pdHess remains false, use bootstrap or profile uncertainty.
If the target still changes across starts, simplify rank or covariance structure before making a scientific claim.

The bottom line

The carry-forward principle is that hard fits need honest triage. A weak Hessian should change the uncertainty workflow, not automatically discard every point estimate the model produced.