Skip to contents

Run an advanced Laplace-consistency check

Source: R/check-consistency.R
gllvmTMB_check_consistency.Rd

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. Must be larger than the number of fixed parameters of the fit (length(fit$tmb_obj$par)): the covariance of the score vector from n_sim replicates has rank at most n_sim - 1, so a smaller n_sim cannot be inverted. The function stops with an error naming the parameter count otherwise. 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 (passed to TMB::checkConsistency()). The refits are kept in $raw; TMB's deviance test of them is summary(cc$raw)$estimate$deviance.p.value. The wrapper does not report it separately. Default FALSE (much faster).

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: estimated bias of each parameter, from TMB's -solve(var(score)) %*% mean(score). It is in each parameter's own units, not a z-score, so the fixed |bias| > 0.5 flag threshold is a heuristic that depends on the parameter's scale.

$joint_p_value

p-value of TMB's joint-score test. This test checks that the simulation itself is correct; it uses the joint gradient and so needs simulated random effects as well as simulated data. NA when the template does not simulate the random effects. It does not depend on estimate.

$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; any other error from TMB's summary() is listed in $warnings).

$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. This needs the TMB template to simulate new responses (a SIMULATE block). The gllvmTMB template does not have one yet (#1582), so every "simulated" dataset would be the observed data and every score vector identical. The function detects this and stops with a "simulation not available" error rather than report a verdict it cannot compute. 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 the flagged-parameter list is large, 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).

See also

sanity_multi() (structural / convergence checks), check_auto_residual() (mixed-family safeguard), gllvmTMB_diagnose() (holistic fit summary), TMB::checkConsistency() (the underlying TMB call).

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"
)

# n_sim must be larger than the number of fitted parameters (8 here);
# a smaller value stops with an error that names the count.
tryCatch(
  gllvmTMB_check_consistency(fit, n_sim = 5, seed = 1),
  error = function(e) cat(conditionMessage(e), "\n")
)
#> `n_sim` = 5 is too small for this fit, which has 8 fixed parameters.
#> ℹ Use `n_sim` of at least 9; more replicates give a steadier test. 

# The template cannot yet simulate new responses (#1582), so with a
# valid n_sim the check stops with a "simulation not available" error.
# Once simulation is available, `cc$diagnostics` gives the overall
# verdict and `cc$flagged_parameters` the parameters whose estimated
# bias exceeds the threshold.
tryCatch(
  gllvmTMB_check_consistency(fit, n_sim = 50, seed = 1),
  error = function(e) cat(conditionMessage(e), "\n")
)
#> The consistency check cannot run: simulation not available for this
#> template.
#> ✖ `TMB::checkConsistency()` needs the template to simulate new responses, but
#>   the gllvmTMB template has no `SIMULATE` block: every replicate would re-score
#>   the observed data and the score covariance is exactly zero.
#> ℹ This is a limitation of the tool, not a property of your data (#1582).
#> → Use `gllvmTMB_diagnose()` to check convergence and the gradient of this fit
#>   instead.