Control parameters for gllvmTMB() — gllvmTMBcontrol • gllvmTMB

Control parameters for gllvmTMB()

Usage

gllvmTMBcontrol(
  d_B = NULL,
  d_W = NULL,
  spde_mode = c("per_trait", "shared"),
  n_init = 1L,
  optimizer = c("nlminb", "optim"),
  optArgs = list(),
  init_jitter = 0.3,
  init_strategy = c("default", "single_trait_warmup"),
  start_method = list(method = NULL, jitter.sd = 0),
  start_from = NULL,
  se = TRUE,
  verbose = FALSE,
  integration = c("laplace", "va"),
  va_H = 7L,
  va_eval_method = c("auto", "jj", "gh"),
  aghq = FALSE,
  aghq_iter_cap = 1L,
  aghq_n_adapt = 400L,
  aghq_ridge = 2,
  aghq_continuation = TRUE,
  aghq_shift_tol = 1e-04,
  aghq_grad_tol = 1e-04,
  aghq_grad_tol_rel = 1e-06,
  aghq_multistart = TRUE,
  aghq_f_tol = 1e-09,
  aghq_escalate_patience = 3L,
  aghq_rho_min = 1/64,
  warn_runaway = TRUE,
  allow_nongaussian_reml = FALSE,
  loading_ridge = NULL,
  optimizer_passes = 1L,
  svd_start = TRUE,
  va_polish = TRUE,
  ...
)

Arguments

d_B, d_W

Not read by the fitter. Retained in the control list for backward compatibility only. Set latent rank in the formula with latent(..., d = K) (or the d argument on the keyword helper), not here.

spde_mode

Not read by the fitter. Retained for backward compatibility only. Spatial field sharing is determined by the spatial_*() formula keyword, not by this control argument.

n_init

Number of random-start replicates. Reduced-rank GLLVMs are often multimodal because the latent-factor likelihood has many equivalent local maxima; running several restarts with different starting values and keeping the fit with the lowest -logLik reduces the risk of settling on a suboptimal solution. Default 1 (the default start, plus the svd_start restart when that applies). Increase to 5–10 for two-level rr models, and to check any fit whose log-likelihood will be compared: a single start can stop at a lower local optimum while reporting convergence = 0, a positive-definite Hessian and a small gradient. The jittered restarts use R's random number generator, so call set.seed() first for a reproducible fit.

optimizer

Optimiser for the native Laplace and AGHQ routes. For the experimental VA route, omission keeps the family registry's automatic primary solver; explicitly set "nlminb" to pin that solver or "optim" to select "BFGS" or "L-BFGS-B" through optArgs$method.

optArgs

Native optimiser arguments. Existing Laplace and AGHQ behavior is unchanged. For VA, nlminb accepts control, lower, upper, and scale; optim accepts a nested control list and method = "BFGS" or "L-BFGS-B" (the latter also accepts bounds). Put native stopping controls such as maxit inside control; unknown or method-incompatible values are rejected before the VA objective is built.

init_jitter

Standard deviation of N(0, sigma) jitter applied to the starting parameter vector across the n_init restarts. Default 0.3.

init_strategy

One of "default" (current behaviour) or "single_trait_warmup". The warmup option fits an intercept-only univariate GLM per trait — with that trait's family — and seeds the matching log_phi_* entries before MakeADFun(). Recommended for count families (especially nbinom2) where the default initialisation can leave the optimiser walking the $(\psi_t, \phi_t)$ trade-off ridge. No effect for traits whose family doesn't carry a phi parameter (Gaussian, Poisson, binomial). Default "default".

start_method

Optional reduced-rank starting-value method, modelled after glmmTMB::glmmTMBControl(start_method = ...) but extended for two-level gllvmTMB models. Use list(method = "res", jitter.sd = 0.2) to seed latent() loadings and latent scores from a reduced-rank decomposition of fixed-effect residuals. Use list(method = "indep") to first fit the matching independent diagonal GLMM/GLLVM and copy its estimated fixed effects, per-trait variance starts, and random effects into the full latent covariance fit. jitter.sd adds Normal jitter to residual-start latent scores. Default list(method = NULL, jitter.sd = 0) keeps the historical starts.

start_from

Optional fitted gllvmTMB object, usually a simpler model such as one latent() tier or an independent diagonal model (indep()). Any estimated TMB parameters with shapes matching the current model are copied into the starting parameter list before optimisation. This implements the "fit simpler, then use it as starting values" workflow recommended for complex reduced-rank models.

se

Logical; if TRUE, compute TMB::sdreport() after optimisation so Wald standard errors and Hessian diagnostics are available. If FALSE, skip standard-error calculation and return the point-estimate fit with sd_report = NULL and a diagnostic note. This is useful for hard models where the point estimate is needed and uncertainty will be obtained by bootstrap or profile methods. Default TRUE.

verbose

If TRUE, prints a one-line summary per restart so the user can see which seed led to the winning fit. Default FALSE.

integration

Which method evaluates the latent-variable integral. "laplace" (default) is the Laplace approximation. For an unpenalised native-Laplace fit, logLik(), AIC(), BIC() and likelihood-ratio tests have their ordinary maximum-likelihood interpretation. AGHQ-refined fits also expose these methods, but use a different integration objective and must not be compared across engines; with a loading ridge, the fitted point is penalised MAP and ordinary likelihood-comparison interpretations do not apply.

"va" selects an opt-in research route whose objective is an evidence lower bound, not a marginal likelihood. Fixed-effect vcov() and confint() use profiled-Schur VA-Wald information and are labelled calibrated = FALSE; getLV() can return a variational posterior SD, which is not a frequentist random-effect standard error. No nominal coverage claim is made. A bound must not be compared across ranks or models, so it cannot be used for model or rank selection.

It is currently admitted only for these model conditions: latent(..., unique = FALSE), the 18 scalar family/link cells in the response-family registry, d up to 2, up to 80 responses, at least 100 units, and the native TMB engine. Requesting it outside these conditions is an error, not a warning. The d limit is retained from the historical explicit-JJ preregistered recovery gate: d = 4 was measured and refused because, with few responses, the planted axes collapsed more often than that gate allowed. This is an implementation admission fence, not recovery evidence for the current automatic GH-H7 route. In the current 36-cell family-by-rank campaign only the Poisson-log q = 5 cell passed the overall point route, and that rank lies outside the public d <= 2 fence; VA therefore remains experimental and calibrated = FALSE. It cannot be combined with aghq, which is an alternative evaluation of the same integral rather than an additional layer.

Offering this value advertises nothing about its accuracy.

The ordination surface is available for a "va" fit: extract_ordination(), getLV(), getLoadings() and extract_loadings() return latent scores and loadings as POINT ESTIMATES. Fixed-effect confint() and vcov() expose uncalibrated VA-Wald uncertainty; getLV(se = TRUE) returns a variational posterior SD, not a frequentist standard error, and only when eval_method resolves to "gh". Loading/covariance intervals remain unavailable.

Historical explicit-JJ evidence. Recovery of the ordination below was measured with the earlier JJ evaluator against planted truth at the admitted cells (d = 2, 8 responses, n = 150 and 400, 50 seeds per cell), with the Laplace route run on the same simulated data as a control. This campaign predates the current automatic Gauss-Hermite default and must not be used as evidence for the accuracy of the current GH-default route:

Gaussian: the loading-implied variance is recovered at 0.98–1.02 of truth, and the mean paired difference in latent-score correlation against Laplace is below 2e-07 at both n.
Poisson: also 0.98–1.02 of truth, with a paired latent-score difference of ~4e-04 — negligible in magnitude but statistically detectable, and the loading scale sits ~1% below Laplace.
Binomial: latent scores recover at about r = 0.59. That ceiling is a limit of binary data at these cells rather than of the historical JJ route — the Laplace default reaches about r = 0.56–0.59 on the same data, and a third estimator (Gauss-Hermite) lands in the same place. The loading SCALE, however, is biased low (0.58–0.67 of truth) and this is a real bias relative to Laplace, whose median brackets 1. Read a binomial ordination's axes and relative positions, not the absolute magnitude of its loadings.

The comparison is not effort-matched: the historical JJ arm ran four starts behind an agreement-and-gradient health gate, the Laplace arm a single start admitted on a positive-definite Hessian.

"va" returns an object of class "gllvmTMB_va" (see gllvmTMB_va-methods) rather than an ordinary fit, so that every method which would treat its objective as a likelihood fails loudly instead of returning a number. Because the engine runs its own multi-start and optimiser policy, optimizer and optArgs select the VA primary solver, while va_polish controls the existing secondary-solver sequence. If optimizer is omitted, the family registry still selects the primary solver. For a fitted VA object, per-start diagnostics are in fit$engine_result$starts[[i]]: primary_optimizer names the selected primary method and final_optimizer names the method whose parameter vector and objective were retained. polish_optimizer summarizes the polish methods actually attempted ("disabled", "not_attempted", "nlminb_attempted", or "nlminb_then_L-BFGS-B_attempted"). The ordered optimizer_attempts records include an accepted flag. The primary result starts accepted; a polish record is accepted only when it replaces the retained result. Acceptance records replacement, not convergence or a healthy fit; inspect convergence and healthy separately. Thus the last accepted attempt identifies final_optimizer; a failed primary has no retained solver and records NA. The settings n_init, start_from, init_* and se have no effect on this route. Any model structure the route cannot represent (a latent term away from the unit grouping, a constrained ordination, an offset, weights, REML, lambda_constraint, Xcoef_fixed, or a further random effect) is an error, never a silent omission.

There is no "eva" value. The EVA engine exists in this package and is reachable as a research route, but it is not wired to gllvmTMB(), and an argument value that could only ever raise an error would advertise a capability the package does not have. No general coefficient, covariance, interval, or calibration claim is made for EVA. Existing measurements are research diagnostics only and do not establish a public estimator route.

va_H

Gauss-Hermite quadrature order for the variational route, a single odd integer >= 3 (default 7). Ignored unless integration = "va".

Quadrature dominates the cost of a variational fit and the work is linear in H, so this is the main speed control on that route. On a binomial probit cell with 20 traits, H = 7 was indistinguishable from H = 61 in both the recovered loading scale and the recovered linear predictor, at two and at five latent dimensions, while running 3.4-6.7 times faster. H = 5 was measurably worse at five latent dimensions. These results support H = 7 as the default; H = 61 remains available as a diagnostic.

va_eval_method

How the variational route evaluates the per-observation expectation. "auto" (default) selects Gauss-Hermite quadrature for every admitted scalar family/link cell and is what you want unless you are comparing methods.

"gh" uses Gauss-Hermite quadrature. "jj" uses the Jaakkola-Jordan bound, which is defined only for pure binomial-logit fits and is substantially cheaper. The two differ in more than speed: on binary data the bound's recovered loading magnitudes are biased downward and do not improve as the sample grows, whereas quadrature's converge. Latent-variable scores, and therefore ordinations, are effectively unaffected by the choice. Ignored unless integration = "va".

aghq

Adaptive Gauss-Hermite quadrature for the between-unit latent block. FALSE (default) fits by Laplace approximation. A positive integer requests that many quadrature nodes. "auto" lets the package decide, and it declines to Laplace whenever the model is ineligible or the expected gain does not justify the cost. Opt-in and experimental: no capability claim is made for quadrature-fitted models, and eligibility is narrow — the random part must be a single ordinary latent() block on the unit tier.

aghq_iter_cap

Optimiser iteration budget for the first adaptation pass. Default 1L. Later passes escalate when continuation is enabled.

aghq_n_adapt

Maximum number of adaptation passes. Default 400L.

aghq_ridge

Ridge penalty on the loadings, as the scale tau; Inf disables it. Default 2. Naming a numeric value explicitly also makes the Laplace + ridge control reachable without quadrature. The default does not penalise an ordinary Laplace fit unless the caller explicitly names aghq_ridge. AGHQ itself is opt-in, but once enabled it uses the default tau = 2 penalised MAP route unless Inf is requested. "auto" is an opt-in experimental route for pure single-trial Bernoulli models with one ordinary unit-tier latent() block. It uses an unpenalised 9-node multi-start AGHQ pilot, then tau = min(6, max(1, ||Lambda_pilot||_F / sqrt(pq))). If the pilot or scale-aware final fit is unusable, the returned fit transparently falls back to the shipped tau = 2 route; inspect fit$aghq$ridge_auto. Evidence supports failure/runaway avoidance in this scope, not a broad loading-accuracy improvement. The measured grid used a logit link with p = 6, q = 2, and n = 100, 400, or 1600; other links and dimensions are extrapolations. The "auto" route always uses 9-node multi-start AGHQ for both pilot and final fits; a conflicting node/start control is replaced with a warning, while an explicit aghq = FALSE is incompatible.

aghq_continuation

If TRUE (default), the adaptation loop may raise aghq_iter_cap across passes. FALSE pins the cap and disables escalation.

aghq_shift_tol, aghq_grad_tol, aghq_f_tol

Convergence tolerances for the adaptation loop: the quadrature-point shift, the gradient at the adaptation point, and the change in objective. Defaults 1e-4, 1e-4, 1e-9.

aghq_grad_tol_rel

Relative gradient tolerance, tested as max|grad| / max(1, |objective|). Default 1e-6. The gradient of a likelihood summed over many observations grows with the sample, so the absolute aghq_grad_tol on its own becomes unreachable at scale. The two are combined with OR, so this can only ever admit more fits as converged, never fewer; set it to 0 to test the absolute tolerance only. Read fit$aghq$converged for the verdict — not fit$opt$convergence, which on the quadrature path records the optimiser's per-pass iteration cap and so reports a limit even on a healthy fit.

aghq_multistart

If TRUE (default), the quadrature is run from two starting points — the Laplace optimum and a data-driven alternative that uses no knowledge of the truth — and the fit with the better final objective is kept. FALSE uses the Laplace warm start only, reproducing the behaviour of earlier versions. The extra start costs one additional adaptation run, and exists because the Laplace optimum is sometimes itself a runaway that the quadrature would otherwise inherit.

aghq_escalate_patience

Number of accepted passes before the iteration cap is escalated. Default 3L.

aghq_rho_min

Minimum trust-region-style acceptance ratio for a pass. Default 1/64.

warn_runaway

If TRUE (default), warn when a binomial latent-variable fit triggers the package's existing runaway-loading diagnostic. The warning is shown once per session for each distinct dataset and model (response, family and link, and loading-matrix shape), so a second, different runaway fit in the same session still warns while a re-fit of the same data and model does not. Set FALSE to silence the fit-time warning; the diagnostic remains available through gllvmTMB_diagnose(). The same switch also governs the multinomial contrast-degeneracy warning (collapsed contrast variance, rail-correlated contrasts, or a collapsed spatial range), which uses its own once-per-session slot so neither family's warning can suppress the other's. Ordinal fits emit no fit-time warning: the ordinal_probit() loading screen reports through check_gllvmTMB() only.

allow_nongaussian_reml

Permit REML = TRUE for non-Gaussian families. Default FALSE; the restricted likelihood is defined by a Gaussian marginalisation, so this is an escape hatch, not an endorsement.

loading_ridge

Integration-neutral alias for aghq_ridge. Use this spelling to request the same loading MAP penalty on a Laplace fit without suggesting that the penalty belongs to AGHQ. The default NULL preserves historical unpenalised Laplace fits. Supply at most one of loading_ridge and aghq_ridge; neither may be combined with estimator = "mspl". A finite positive value is the standard deviation tau of independent zero-mean normal priors on free ordinary between-unit loadings: loading_ridge = 2 means variance 2^2, smaller positive values apply stronger shrinkage, and Inf disables the penalty. For unstable binomial fits with ordinary between-unit latent loadings, an explicit loading_ridge = 2 fit is a starting sensitivity analysis; this does not change the ordinary Laplace default. Estimates are penalised MAP points. Retained penalised-Hessian LV uncertainty is a local posterior approximation conditional on tau, with delta-method propagation for trait effects, not calibrated repeated-sampling coverage. Older fits without the required metadata withhold LV uncertainty.

optimizer_passes

Number of exact-gradient optimisation passes from the preceding estimate. The default 1 retains the historical single pass. A later pass is retained only when it converges and does not increase the objective. It is available for native Laplace fits with aghq = FALSE. The setting is saved in the public call and replayed by update().

svd_start

Logical. If TRUE (default), a fit with an ordinary latent() block is also optimised from a second, data-informed start whose loadings and latent scores come from one singular value decomposition of the residual matrix, and the start with the lower objective is kept (restart_history labels it "svd"). The default start alone sometimes stopped at a lower local optimum with clean diagnostics. This adds one optimisation per fit and uses no random numbers. It is skipped when start_from or start_method is supplied and on the estimator = "mspl" route. Set FALSE for the single-start behaviour of gllvmTMB 0.8.0 and earlier.

va_polish

For VA fits, whether to run the existing secondary solver polishing sequence after each primary start. Defaults to TRUE to retain prior behavior; FALSE disables every secondary solver and is useful for controlled optimizer comparisons. The four primary starts and their existing agreement rule are unchanged. BFGS uses a dense inverse-Hessian approximation, so its memory use grows quadratically with the parameter count; use it as an opt-in method for modest models.

...

Reserved for future use. Currently ignored with a warning.

Value

A gllvmTMBcontrol list.

Details

Recommended workflow for two-level latent models. Reduced-rank likelihoods are notoriously multimodal in two-level fits like value ~ 0 + trait + latent(0+trait|ID, d_B) + latent(0+trait|obs_ID, d_W). The recommended workflow is:

Run several restarts (n_init = 5 to 10) and keep the one with the lowest -logLik.
If that fails to converge cleanly, switch the optimiser to optim with BFGS: gllvmTMBcontrol(optimizer = "optim", optArgs = list(method = "BFGS")).
start_method = list(method = "res", jitter.sd = 0.2) is soft-deprecated as of 0.6.0 and warns once per session. It still fits: it decomposes the fixed-effect residual matrix into starting values for $\Lambda$ and the latent scores. Do not reach for it.

It was retired on measurement. Across 89 simulated fits — Gaussian, Poisson and negative-binomial, d = 1 to 3, three and five traits — it was never materially better than the default start (its three best margins were 0.07, 0.29 and 0.66 log-likelihood units) and was materially worse eight times, once by 14.6 units. At d >= 2 it made no difference at all. Every failure reported convergence == 0 and a positive-definite Hessian on both sides, so the worse fit was silent, and restarts did not detect it: in one such fit all five n_init = 5 restarts returned the same worse optimum.

If you use it anyway, compare fit$opt$objective against a default-start fit and keep the lower value.
For Gaussian two-level models, prefer start_method = list(method = "indep") or manually fit a simpler model and pass it through start_from = simpler_fit. This is a GLMM warm start rather than a fixed-effect-only GLM residual start.

Examples

# gllvmTMBcontrol() is a pure-R constructor: it builds and validates a
# control list without fitting anything.
gllvmTMBcontrol()                                  # historical defaults
#> $d_B
#> NULL
#> 
#> $d_W
#> NULL
#> 
#> $spde_mode
#> [1] "per_trait"
#> 
#> $integration
#> [1] "laplace"
#> 
#> $va_H
#> [1] 7
#> 
#> $va_eval_method
#> [1] "auto"
#> 
#> $n_init
#> [1] 1
#> 
#> $n_init_explicit
#> [1] FALSE
#> 
#> $svd_start
#> [1] TRUE
#> 
#> $svd_start_explicit
#> [1] FALSE
#> 
#> $optimizer
#> [1] "nlminb"
#> 
#> $optimizer_explicit
#> [1] FALSE
#> 
#> $optArgs
#> list()
#> 
#> $va_polish
#> [1] TRUE
#> 
#> $optimizer_passes
#> [1] 1
#> 
#> $init_jitter
#> [1] 0.3
#> 
#> $init_strategy
#> [1] "default"
#> 
#> $start_method
#> $start_method$method
#> NULL
#> 
#> $start_method$jitter.sd
#> [1] 0
#> 
#> 
#> $start_from
#> NULL
#> 
#> $se
#> [1] TRUE
#> 
#> $verbose
#> [1] FALSE
#> 
#> $aghq
#> [1] FALSE
#> 
#> $aghq_iter_cap
#> [1] 1
#> 
#> $aghq_n_adapt
#> [1] 400
#> 
#> $aghq_ridge
#> [1] 2
#> 
#> $aghq_ridge_explicit
#> [1] FALSE
#> 
#> $loading_ridge
#> NULL
#> 
#> $loading_ridge_explicit
#> [1] FALSE
#> 
#> $aghq_continuation
#> [1] TRUE
#> 
#> $aghq_shift_tol
#> [1] 1e-04
#> 
#> $aghq_grad_tol
#> [1] 1e-04
#> 
#> $aghq_grad_tol_rel
#> [1] 1e-06
#> 
#> $aghq_multistart
#> [1] TRUE
#> 
#> $aghq_f_tol
#> [1] 1e-09
#> 
#> $aghq_escalate_patience
#> [1] 3
#> 
#> $aghq_rho_min
#> [1] 0.015625
#> 
#> $warn_runaway
#> [1] TRUE
#> 
#> $allow_nongaussian_reml
#> [1] FALSE
#> 

# Multi-start to guard against multimodal reduced-rank likelihoods, with
# jitter on the starting parameter vector across restarts.
gllvmTMBcontrol(n_init = 3, init_jitter = 0.2)
#> $d_B
#> NULL
#> 
#> $d_W
#> NULL
#> 
#> $spde_mode
#> [1] "per_trait"
#> 
#> $integration
#> [1] "laplace"
#> 
#> $va_H
#> [1] 7
#> 
#> $va_eval_method
#> [1] "auto"
#> 
#> $n_init
#> [1] 3
#> 
#> $n_init_explicit
#> [1] TRUE
#> 
#> $svd_start
#> [1] TRUE
#> 
#> $svd_start_explicit
#> [1] FALSE
#> 
#> $optimizer
#> [1] "nlminb"
#> 
#> $optimizer_explicit
#> [1] FALSE
#> 
#> $optArgs
#> list()
#> 
#> $va_polish
#> [1] TRUE
#> 
#> $optimizer_passes
#> [1] 1
#> 
#> $init_jitter
#> [1] 0.2
#> 
#> $init_strategy
#> [1] "default"
#> 
#> $start_method
#> $start_method$method
#> NULL
#> 
#> $start_method$jitter.sd
#> [1] 0
#> 
#> 
#> $start_from
#> NULL
#> 
#> $se
#> [1] TRUE
#> 
#> $verbose
#> [1] FALSE
#> 
#> $aghq
#> [1] FALSE
#> 
#> $aghq_iter_cap
#> [1] 1
#> 
#> $aghq_n_adapt
#> [1] 400
#> 
#> $aghq_ridge
#> [1] 2
#> 
#> $aghq_ridge_explicit
#> [1] FALSE
#> 
#> $loading_ridge
#> NULL
#> 
#> $loading_ridge_explicit
#> [1] FALSE
#> 
#> $aghq_continuation
#> [1] TRUE
#> 
#> $aghq_shift_tol
#> [1] 1e-04
#> 
#> $aghq_grad_tol
#> [1] 1e-04
#> 
#> $aghq_grad_tol_rel
#> [1] 1e-06
#> 
#> $aghq_multistart
#> [1] TRUE
#> 
#> $aghq_f_tol
#> [1] 1e-09
#> 
#> $aghq_escalate_patience
#> [1] 3
#> 
#> $aghq_rho_min
#> [1] 0.015625
#> 
#> $warn_runaway
#> [1] TRUE
#> 
#> $allow_nongaussian_reml
#> [1] FALSE
#> 

# Switch the optimiser to optim + BFGS for finicky two-level rr fits.
gllvmTMBcontrol(optimizer = "optim", optArgs = list(method = "BFGS"))
#> $d_B
#> NULL
#> 
#> $d_W
#> NULL
#> 
#> $spde_mode
#> [1] "per_trait"
#> 
#> $integration
#> [1] "laplace"
#> 
#> $va_H
#> [1] 7
#> 
#> $va_eval_method
#> [1] "auto"
#> 
#> $n_init
#> [1] 1
#> 
#> $n_init_explicit
#> [1] FALSE
#> 
#> $svd_start
#> [1] TRUE
#> 
#> $svd_start_explicit
#> [1] FALSE
#> 
#> $optimizer
#> [1] "optim"
#> 
#> $optimizer_explicit
#> [1] TRUE
#> 
#> $optArgs
#> $optArgs$method
#> [1] "BFGS"
#> 
#> 
#> $va_polish
#> [1] TRUE
#> 
#> $optimizer_passes
#> [1] 1
#> 
#> $init_jitter
#> [1] 0.3
#> 
#> $init_strategy
#> [1] "default"
#> 
#> $start_method
#> $start_method$method
#> NULL
#> 
#> $start_method$jitter.sd
#> [1] 0
#> 
#> 
#> $start_from
#> NULL
#> 
#> $se
#> [1] TRUE
#> 
#> $verbose
#> [1] FALSE
#> 
#> $aghq
#> [1] FALSE
#> 
#> $aghq_iter_cap
#> [1] 1
#> 
#> $aghq_n_adapt
#> [1] 400
#> 
#> $aghq_ridge
#> [1] 2
#> 
#> $aghq_ridge_explicit
#> [1] FALSE
#> 
#> $loading_ridge
#> NULL
#> 
#> $loading_ridge_explicit
#> [1] FALSE
#> 
#> $aghq_continuation
#> [1] TRUE
#> 
#> $aghq_shift_tol
#> [1] 1e-04
#> 
#> $aghq_grad_tol
#> [1] 1e-04
#> 
#> $aghq_grad_tol_rel
#> [1] 1e-06
#> 
#> $aghq_multistart
#> [1] TRUE
#> 
#> $aghq_f_tol
#> [1] 1e-09
#> 
#> $aghq_escalate_patience
#> [1] 3
#> 
#> $aghq_rho_min
#> [1] 0.015625
#> 
#> $warn_runaway
#> [1] TRUE
#> 
#> $allow_nongaussian_reml
#> [1] FALSE
#> 

# Compare VA primary solvers on the same model and four seeded starts.
# This example is deliberately not run automatically because each method
# performs four variational fits. Reset the seed for each method and disable
# polish so the comparison isolates the selected primary solver.
if (FALSE) { # \dontrun{
set.seed(20260802L)
n <- 120L
p <- 4L
va_data <- data.frame(
  y = rbinom(n * p, 1L, 0.5),
  trait = factor(rep(seq_len(p), times = n)),
  site = factor(rep(seq_len(n), each = p))
)
va_formula <- y ~ 0 + trait + latent(0 + trait | site, d = 1,
                                     unique = FALSE)
va_methods <- list(
  nlminb = list(optimizer = "nlminb",
                optArgs = list(control = list(iter.max = 2000L,
                                               eval.max = 2000L))),
  BFGS = list(optimizer = "optim",
              optArgs = list(method = "BFGS",
                              control = list(maxit = 2000L))),
  `L-BFGS-B` = list(optimizer = "optim",
                    optArgs = list(method = "L-BFGS-B",
                                    control = list(maxit = 2000L)))
)
va_fits <- lapply(va_methods, function(method) {
  set.seed(1595L)
  gllvmTMB(
    va_formula, data = va_data,
    family = stats::binomial(link = "probit"), unit = "site",
    control = gllvmTMBcontrol(
      integration = "va", va_H = 7L, va_eval_method = "gh",
      optimizer = method$optimizer, optArgs = method$optArgs,
      va_polish = FALSE
    )
  )
})
lapply(va_fits, function(fit) {
  starts <- fit$engine_result$starts
  data.frame(
    start = vapply(starts, `[[`, integer(1), "start"),
    primary = vapply(starts, `[[`, character(1), "primary_optimizer"),
    final = vapply(starts, `[[`, character(1), "final_optimizer"),
    polish = vapply(starts, `[[`, character(1), "polish_optimizer"),
    objective = vapply(starts, `[[`, numeric(1), "objective"),
    convergence = vapply(starts, `[[`, integer(1), "convergence")
  )
})
} # }

# Experimental, opt-in scale-aware ridge for the evaluated Bernoulli AGHQ
# scope. The route uses a fixed 9-node multi-start rule; this is not an
# interval-calibration certificate.
gllvmTMBcontrol(aghq_ridge = "auto")
#> $d_B
#> NULL
#> 
#> $d_W
#> NULL
#> 
#> $spde_mode
#> [1] "per_trait"
#> 
#> $integration
#> [1] "laplace"
#> 
#> $va_H
#> [1] 7
#> 
#> $va_eval_method
#> [1] "auto"
#> 
#> $n_init
#> [1] 1
#> 
#> $n_init_explicit
#> [1] FALSE
#> 
#> $svd_start
#> [1] TRUE
#> 
#> $svd_start_explicit
#> [1] FALSE
#> 
#> $optimizer
#> [1] "nlminb"
#> 
#> $optimizer_explicit
#> [1] FALSE
#> 
#> $optArgs
#> list()
#> 
#> $va_polish
#> [1] TRUE
#> 
#> $optimizer_passes
#> [1] 1
#> 
#> $init_jitter
#> [1] 0.3
#> 
#> $init_strategy
#> [1] "default"
#> 
#> $start_method
#> $start_method$method
#> NULL
#> 
#> $start_method$jitter.sd
#> [1] 0
#> 
#> 
#> $start_from
#> NULL
#> 
#> $se
#> [1] TRUE
#> 
#> $verbose
#> [1] FALSE
#> 
#> $aghq
#> [1] 9
#> 
#> $aghq_iter_cap
#> [1] 1
#> 
#> $aghq_n_adapt
#> [1] 400
#> 
#> $aghq_ridge
#> [1] "auto"
#> 
#> $aghq_ridge_explicit
#> [1] TRUE
#> 
#> $loading_ridge
#> NULL
#> 
#> $loading_ridge_explicit
#> [1] FALSE
#> 
#> $aghq_continuation
#> [1] TRUE
#> 
#> $aghq_shift_tol
#> [1] 1e-04
#> 
#> $aghq_grad_tol
#> [1] 1e-04
#> 
#> $aghq_grad_tol_rel
#> [1] 1e-06
#> 
#> $aghq_multistart
#> [1] TRUE
#> 
#> $aghq_f_tol
#> [1] 1e-09
#> 
#> $aghq_escalate_patience
#> [1] 3
#> 
#> $aghq_rho_min
#> [1] 0.015625
#> 
#> $warn_runaway
#> [1] TRUE
#> 
#> $allow_nongaussian_reml
#> [1] FALSE
#> 
# The same fixed loading MAP penalty, named independently of integration.
gllvmTMBcontrol(loading_ridge = 2)
#> $d_B
#> NULL
#> 
#> $d_W
#> NULL
#> 
#> $spde_mode
#> [1] "per_trait"
#> 
#> $integration
#> [1] "laplace"
#> 
#> $va_H
#> [1] 7
#> 
#> $va_eval_method
#> [1] "auto"
#> 
#> $n_init
#> [1] 1
#> 
#> $n_init_explicit
#> [1] FALSE
#> 
#> $svd_start
#> [1] TRUE
#> 
#> $svd_start_explicit
#> [1] FALSE
#> 
#> $optimizer
#> [1] "nlminb"
#> 
#> $optimizer_explicit
#> [1] FALSE
#> 
#> $optArgs
#> list()
#> 
#> $va_polish
#> [1] TRUE
#> 
#> $optimizer_passes
#> [1] 1
#> 
#> $init_jitter
#> [1] 0.3
#> 
#> $init_strategy
#> [1] "default"
#> 
#> $start_method
#> $start_method$method
#> NULL
#> 
#> $start_method$jitter.sd
#> [1] 0
#> 
#> 
#> $start_from
#> NULL
#> 
#> $se
#> [1] TRUE
#> 
#> $verbose
#> [1] FALSE
#> 
#> $aghq
#> [1] FALSE
#> 
#> $aghq_iter_cap
#> [1] 1
#> 
#> $aghq_n_adapt
#> [1] 400
#> 
#> $aghq_ridge
#> [1] 2
#> 
#> $aghq_ridge_explicit
#> [1] TRUE
#> 
#> $loading_ridge
#> [1] 2
#> 
#> $loading_ridge_explicit
#> [1] TRUE
#> 
#> $aghq_continuation
#> [1] TRUE
#> 
#> $aghq_shift_tol
#> [1] 1e-04
#> 
#> $aghq_grad_tol
#> [1] 1e-04
#> 
#> $aghq_grad_tol_rel
#> [1] 1e-06
#> 
#> $aghq_multistart
#> [1] TRUE
#> 
#> $aghq_f_tol
#> [1] 1e-09
#> 
#> $aghq_escalate_patience
#> [1] 3
#> 
#> $aghq_rho_min
#> [1] 0.015625
#> 
#> $warn_runaway
#> [1] TRUE
#> 
#> $allow_nongaussian_reml
#> [1] FALSE
#> 

# Experimental scalar variational route: public auto uses GH with H = 7.
gllvmTMBcontrol(integration = "va")
#> $d_B
#> NULL
#> 
#> $d_W
#> NULL
#> 
#> $spde_mode
#> [1] "per_trait"
#> 
#> $integration
#> [1] "va"
#> 
#> $va_H
#> [1] 7
#> 
#> $va_eval_method
#> [1] "auto"
#> 
#> $n_init
#> [1] 1
#> 
#> $n_init_explicit
#> [1] FALSE
#> 
#> $svd_start
#> [1] TRUE
#> 
#> $svd_start_explicit
#> [1] FALSE
#> 
#> $optimizer
#> [1] "nlminb"
#> 
#> $optimizer_explicit
#> [1] FALSE
#> 
#> $optArgs
#> list()
#> 
#> $va_polish
#> [1] TRUE
#> 
#> $optimizer_passes
#> [1] 1
#> 
#> $init_jitter
#> [1] 0.3
#> 
#> $init_strategy
#> [1] "default"
#> 
#> $start_method
#> $start_method$method
#> NULL
#> 
#> $start_method$jitter.sd
#> [1] 0
#> 
#> 
#> $start_from
#> NULL
#> 
#> $se
#> [1] TRUE
#> 
#> $verbose
#> [1] FALSE
#> 
#> $aghq
#> [1] FALSE
#> 
#> $aghq_iter_cap
#> [1] 1
#> 
#> $aghq_n_adapt
#> [1] 400
#> 
#> $aghq_ridge
#> [1] 2
#> 
#> $aghq_ridge_explicit
#> [1] FALSE
#> 
#> $loading_ridge
#> NULL
#> 
#> $loading_ridge_explicit
#> [1] FALSE
#> 
#> $aghq_continuation
#> [1] TRUE
#> 
#> $aghq_shift_tol
#> [1] 1e-04
#> 
#> $aghq_grad_tol
#> [1] 1e-04
#> 
#> $aghq_grad_tol_rel
#> [1] 1e-06
#> 
#> $aghq_multistart
#> [1] TRUE
#> 
#> $aghq_f_tol
#> [1] 1e-09
#> 
#> $aghq_escalate_patience
#> [1] 3
#> 
#> $aghq_rho_min
#> [1] 0.015625
#> 
#> $warn_runaway
#> [1] TRUE
#> 
#> $allow_nongaussian_reml
#> [1] FALSE
#> 
# H = 61 remains a diagnostic; JJ is explicit binomial-logit only.
gllvmTMBcontrol(integration = "va", va_H = 61, va_eval_method = "gh")
#> $d_B
#> NULL
#> 
#> $d_W
#> NULL
#> 
#> $spde_mode
#> [1] "per_trait"
#> 
#> $integration
#> [1] "va"
#> 
#> $va_H
#> [1] 61
#> 
#> $va_eval_method
#> [1] "gh"
#> 
#> $n_init
#> [1] 1
#> 
#> $n_init_explicit
#> [1] FALSE
#> 
#> $svd_start
#> [1] TRUE
#> 
#> $svd_start_explicit
#> [1] FALSE
#> 
#> $optimizer
#> [1] "nlminb"
#> 
#> $optimizer_explicit
#> [1] FALSE
#> 
#> $optArgs
#> list()
#> 
#> $va_polish
#> [1] TRUE
#> 
#> $optimizer_passes
#> [1] 1
#> 
#> $init_jitter
#> [1] 0.3
#> 
#> $init_strategy
#> [1] "default"
#> 
#> $start_method
#> $start_method$method
#> NULL
#> 
#> $start_method$jitter.sd
#> [1] 0
#> 
#> 
#> $start_from
#> NULL
#> 
#> $se
#> [1] TRUE
#> 
#> $verbose
#> [1] FALSE
#> 
#> $aghq
#> [1] FALSE
#> 
#> $aghq_iter_cap
#> [1] 1
#> 
#> $aghq_n_adapt
#> [1] 400
#> 
#> $aghq_ridge
#> [1] 2
#> 
#> $aghq_ridge_explicit
#> [1] FALSE
#> 
#> $loading_ridge
#> NULL
#> 
#> $loading_ridge_explicit
#> [1] FALSE
#> 
#> $aghq_continuation
#> [1] TRUE
#> 
#> $aghq_shift_tol
#> [1] 1e-04
#> 
#> $aghq_grad_tol
#> [1] 1e-04
#> 
#> $aghq_grad_tol_rel
#> [1] 1e-06
#> 
#> $aghq_multistart
#> [1] TRUE
#> 
#> $aghq_f_tol
#> [1] 1e-09
#> 
#> $aghq_escalate_patience
#> [1] 3
#> 
#> $aghq_rho_min
#> [1] 0.015625
#> 
#> $warn_runaway
#> [1] TRUE
#> 
#> $allow_nongaussian_reml
#> [1] FALSE
#>