Run an advanced Laplace-consistency check

gllvmTMB_check_consistency() is an advanced validation helper for the Laplace approximation, not a first-use diagnostic. Start with check_gllvmTMB() or gllvmTMB_diagnose(). Use this function when a fit looks numerically acceptable but you need a direct check of whether the approximate marginal score is centred.

Usage

gllvmTMB_check_consistency(fit, n_sim = 100L, seed = NULL, estimate = FALSE)

Arguments

fit: A fit returned by gllvmTMB().
n_sim: Integer number of simulate-evaluate replicates. Default 100L. Cost is roughly n_sim joint likelihood evaluations; budget 5-30 seconds on a Tier-1 fixture and proportionally longer for larger fits.
seed: Optional integer seed for the simulation RNG. NULL (default) draws a random seed.
estimate: Logical. When TRUE, also refits the model on each simulated dataset and reports the joint-score p_value + bias. Default FALSE (much faster; marginal-score check only).

Value

An object of class gllvmTMB_check_consistency with components:

$marginal_p_value: Joint chi-squared p-value across all marginal-score components. NA if TMB's information matrix could not be inverted (often the case on tiny fixtures).
$marginal_bias: Named numeric vector: per-parameter bias (mean of the simulated score, normalised by its SE). Large absolute values flag specific parameters where the Laplace approximation is biased.
$joint_p_value: Joint-score p-value. NA unless estimate = TRUE. Slower because each replicate refits.
$flagged_parameters: Character names of parameters whose marginal bias exceeds |bias| > 0.5 (a heuristic; finite-n_sim noise can produce small spurious biases, but |0.5| is generally beyond noise for n_sim = 100).
$n_sim: The n_sim actually run.
$threshold: The marginal-bias threshold used.
$diagnostics: Character vector of one or more of: "no_centring_warning" (the finite simulation did not detect non-centring; this is not an accuracy certificate), "marginal_score_non_centred", "joint_score_non_centred", "information_matrix_singular" (TMB couldn't invert it; interpret with caution), "marginal_p_value_unavailable" (TMB returned NA without a captured warning; usually happens on tiny / weakly-identified fixtures).
$raw: The full TMB::checkConsistency() return value in case the user wants the per-replicate gradient matrix.
$call: match.call() of the invocation.

Details

It simulates n_sim datasets from the fitted model under the joint parameter vector at the MLE, then tests whether the approximate marginal score function is centred at zero across the simulations. A non-centred score is a sign that the Laplace approximation is unreliable for this fit – typically because the random-effects posterior is far from Gaussian (saturated binomial / Beta, sparse counts, weakly identified random effects) or the data don't constrain the random effects well.

Scope: local Laplace-consistency diagnostics for a fitted model. This is a diagnostic signal rather than a substitute for bootstrap, profile, or external Bayesian calibration. Broader calibration evidence is not yet available.

This is a complementary signal to sanity_multi() (which checks Hessian definiteness, gradient magnitude, convergence flags) and to a target-explicit known-DGP study of rank-selection behaviour. Slower than the first-line checks – n_sim likelihood evaluations – but the only diagnostic that targets the Laplace approximation itself rather than the optimisation outcome or the parameter identification.

A large marginal-score p-value is an absence of evidence against centring in this finite simulation; it is not a certificate that the Laplace approximation is accurate. If the result is small (p_value <= 0.05) or large flagged-parameter list), consider:

Refitting with a richer random-effects structure (so the conditional posterior is more Gaussian).
Switching to a parametric bootstrap CI (which does not depend on the Laplace approximation being a good fit).
Validating against tmbstan::tmbstan(fit$tmb_obj) (the audit's recommended Bayesian-comparison path).

Examples

set.seed(1)
# Counts of four species at 100 sites, driven by one shared latent axis.
# The check matters most for non-Gaussian responses like these, where
# the Laplace approximation is not exact.
n_sites <- 100
z <- rnorm(n_sites)
loadings <- c(0.8, 0.6, 0.5, 0.3)
dat <- expand.grid(site = factor(seq_len(n_sites)),
                   trait = factor(paste0("sp", 1:4)))
mu <- exp(0.5 + loadings[as.integer(dat$trait)] * z[as.integer(dat$site)])
dat$count <- rpois(nrow(dat), mu)
fit <- gllvmTMB(
  count ~ 0 + trait + latent(0 + trait | site, d = 1, unique = FALSE),
  data = dat, family = poisson(), trait = "trait", unit = "site"
)

# Run the check with 50 replicates. n_sim must be larger than the
# number of fitted parameters (8 here); more replicates give a
# steadier test.
cc <- gllvmTMB_check_consistency(fit, n_sim = 50, seed = 1)
cc
#> 
#> ── gllvmTMB Laplace-consistency check ──────────────────────────────────────────
#> • n_sim: 50
#> • Marginal-score p-value: NA (information matrix singular)
#> • Joint-score p-value: NA (estimate = FALSE; pass estimate = TRUE to compute)
#> ! Diagnostics: "information_matrix_singular"
#> The information matrix could not be inverted. This is common on tiny /
#> weakly-identified fixtures; interpret the marginal p-value with caution and
#> consider increasing `n_sim` or fitting on a richer dataset.
#> 
#> ── Marginal bias per parameter ──
#> 
#>      b_fix      b_fix      b_fix      b_fix theta_rr_B theta_rr_B theta_rr_B 
#>         NA         NA         NA         NA         NA         NA         NA 
#> theta_rr_B 
#>         NA 

# The overall verdict and any parameters whose estimated bias exceeds
# the threshold. "no_centring_warning" means no sign of trouble was
# found (not a guarantee); "information_matrix_singular" means the
# check could not reach a verdict either way.
cc$diagnostics
#> [1] "information_matrix_singular"
cc$flagged_parameters
#> character(0)

Usage

Arguments

Value

Details

See also

Examples