Validation oracles: what an external package can and cannot check

Advanced evidence note. This page describes the package’s validation architecture; it is not a beginner modelling guide or a certificate for any particular fit.

It explains what a validation oracle is, why gllvmTMB relies on more than its own simulations, and which parts of the covariance grid currently have an independent check and which do not. Read it alongside Fit diagnostics, which asks whether a particular fit is numerically trustworthy; this page instead asks whether the model class itself has been checked against outside evidence.

What a validation oracle is

A validation oracle is an independent implementation of the same statistical model — written by someone else, from their own reading of the same theory, in a different codebase. If gllvmTMB and an external package independently arrive at the same answer on a well-defined quantity (a log-likelihood, a set of loadings up to rotation, a residual correlation), that agreement is real evidence: it is unlikely that two differently written implementations share the same coding mistake.

Why simulation recovery is not enough on its own

Most of gllvmTMB’s own validation evidence is simulation recovery: generate data from a known model, fit it, and check that the estimates land near the truth. That is necessary, but it has a structural blind spot. A simulator and a fitter written by the same person, from the same understanding of a model, can be wrong in the same direction and still agree with each other perfectly — the recovery check would then confirm only that the code reproduces its own assumptions, not that those assumptions are correct.

An externally written package does not share that specific misunderstanding. It was written by different people, from a different codebase, often with a different parameterisation. Agreement with it is evidence that both implementations encode the same right model — not just that gllvmTMB’s fitter is self-consistent with gllvmTMB’s own simulator. This is why the map below matters: it is not a completeness score, it is a record of where that second, independent kind of evidence exists.

Three traps to avoid

Three comparisons look like a straightforward external check and are not. All three apply throughout the map below.

The estimator trap: variational approximation is not Laplace

gllvmTMB fits by the TMB Laplace approximation by default. Its opt-in scalar VA route is a separately fenced research estimator and must be compared to the same target rather than treated as Laplace. gllvm defaults to a variational approximation (method = "VA", falling back to "EVA" for families without a VA implementation) and only switches to Laplace (method = "LA") on request, for the subset of families that support it. A VA fit optimises a lower bound on the log-likelihood, not the log-likelihood itself, and that bound is biased downward. A raw log-likelihood comparison between a gllvmTMB fit and a default gllvm fit is therefore not comparing the same quantity — the two numbers can differ for a reason that has nothing to do with either implementation being wrong.

Two ways around this: request gllvm(..., method = "LA") where that family supports it, so both sides use the same estimator; or, where only VA/EVA is available, loosen the comparison from “do the numbers match” to “did both optimisers reach a comparably good fit,” using a generous relative tolerance rather than an exact match.

The posterior-mean trap: a Bayesian mean is not an MLE

Hmsc and MCMCglmm are both Bayesian: they return a posterior summary (typically a posterior mean), not a maximum-likelihood estimate. A posterior mean and gllvmTMB’s MLE are different estimands, and they carry different priors, defaults, and shrinkage behaviour into the comparison. If the two disagree, that is not automatically evidence of an engine bug — the prior may simply pull the answer elsewhere. If they agree, that is not automatically proof that gllvmTMB’s MLE is correct — only evidence that the two estimates sit in the same neighbourhood. A Bayesian comparator can corroborate that a covariance structure is roughly the right shape; it cannot, by itself, settle a likelihood-level question.

The tmbstan trap: the same objective function is not independence

tmbstan looks like a Stan-based oracle and is not one. It does not re-derive the likelihood; it reuses the same TMB objective function and gradient that gllvmTMB already builds, and hands them to Stan’s sampler. Agreement with tmbstan therefore shows only that TMB’s own objective is internally consistent under a different integration scheme — it cannot catch an error in the objective itself, because it is the same objective.

tmbstan answers a different and genuinely useful question instead: is the Laplace approximation adequate? Its default (laplace = FALSE) samples the fixed and random effects jointly with the No-U-Turn sampler, so it does not use the Laplace approximation at all — which is precisely what makes it a check on Laplace: compare that posterior against the Laplace-based estimates gllvmTMB returns. The laplace = TRUE option does the opposite, integrating the random effects by Laplace and sampling only the fixed effects; the package documents it as generally not recommended, usually slower and less accurate, and it yields no posterior samples for the random effects. The option whose name suggests it is the way to test the Laplace approximation is the one that does not.

tmbstan is available on CRAN (Monnahan & Kristensen 2018, PLoS ONE 13(5): e0197954). Reach for it when the question is whether the Laplace approximation holds — not as an external check on the likelihood.

The oracle map

The grid below is the same source-by-mode grid used throughout gllvmTMB: five correlation sources (no prefix, phylo_, animal_, spatial_, kernel_) crossed with three covariance modes (indep, dep, latent). See the formula keyword grid for the full syntax. The one-shared-variance case (common = TRUE on indep(), which the soft-deprecated scalar() spelling still names) is called out separately inside each indep cell, because several packages that cannot check plain indep or dep can still check this narrower case.

Status legend

• Tested. An automated comparator already exists in gllvmTMB’s own test suite, comparing a gllvmTMB fit against an external package’s fit on simulated data with a known answer, and it passes.
• Confirmed, not yet asserted. The external package’s mechanism was checked directly — by fitting it and reading what it actually estimates — and it matches the cell’s definition. gllvmTMB does not yet run an automated agreement test against it.
• Documented, not yet fit. The external package’s own documentation describes a mechanism that should match the cell’s definition, but neither side has actually been fit and compared.
• Uncertain. Whether any external reference exists for this cell has not been established either way; it needs further scouting before it can be built or ruled out.
• None. No third-party package implements this cell’s model. Simulation recovery is the only evidence currently in hand — a hand-written Stan reference is a separate, unbuilt route, discussed below rather than shown in the table, because it is available in principle for every cell here, not something specific to a “None” verdict.

Cells marked † rest on a Bayesian comparator — MCMCglmm or Hmsc — and so carry the posterior-mean caveat above, rather than repeating it in every cell.

The table

Source × mode: what can check each cell

Source	indep()	dep()	latent()
none	Shared-variance case: confirmed, not yet asserted (a gllvm fit with no similarity matrix supplied was checked directly and returns exactly one shared variance across traits, the same shape as common = TRUE); documented, not yet built (glmmTMB/lme4’s shared-variance random effect). Plain indep: documented, not yet built into the suite (glmmTMB/lme4’s diag(0 + trait | g)).	Documented, not yet built (glmmTMB/lme4’s us(0 + trait | g)).	Tested — glmmTMB::rr()+diag(), and gllvm’s own unconstrained ordination (Poisson and binary responses; loadings compared after a Procrustes rotation).
phylo_	Shared-variance case: tested against glmmTMB::propto(); separately, gllvm’s colMat mechanism is documented, not yet fit — the one-shared-variance shape was verified by fitting without colMat, and the colMat extension follows from gllvm’s documented behaviour rather than from a fit. Plain indep: documented, not yet fit (MCMCglmm’s idh(trait):animal with a tree-derived inverse relationship matrix)†.	Documented, not yet fit (MCMCglmm’s us(trait):animal variant)†.	None.
animal_	Shared-variance case: documented, not yet fit (gllvm’s colMat, with the pedigree supplied as the similarity matrix — the same mechanism as the phylo_ row, and established the same way, from documentation rather than a fit); documented (MCMCglmm’s idv() with a pedigree, its native use case†; glmmTMB::propto() with a pedigree-derived relationship matrix). Plain indep: documented, not yet fit (MCMCglmm’s idh() + pedigree)†.	Documented, not yet fit (MCMCglmm’s us() + pedigree)†.	None.
spatial_	Uncertain. gllvm’s own spatial mechanisms do not match this cell’s definition (see below); whether another package’s mechanism does was not established.	Uncertain, for the same reason.	Documented, not yet fit (gllvm’s lvCor(corExp) / lvCor(corMatern) — and note lvCor is a dense exact Gaussian process while gllvmTMB uses an SPDE approximation, so even a successful comparison would agree only approximately, not exactly; Hmsc’s spatial latent factors†).
kernel_	Shared-variance case: documented, not yet fit (gllvm’s colMat accepting a bare similarity matrix, per its documentation); documented (MCMCglmm’s idv() with a supplied relationship matrix)†. Plain indep: documented, not yet fit (MCMCglmm’s idh() with a supplied relationship matrix)†.	Documented, not yet fit (MCMCglmm’s us() with a supplied relationship matrix)†.	None.

unique = TRUE (the diagonal companion added to a *_latent() term) is not a separate column here. It inherits the same answer as its row’s indep cell — including where the loadings-only part of that row has a partial reference. spatial_latent(..., unique = TRUE), for example, would still have no external check on its diagonal component even though the loadings-only part has a documented, unfit candidate.

Why gllvm cannot check indep() or dep(), for any source

This is worth stating plainly because it is easy to assume otherwise: gllvm’s colMat mechanism — its one route to a source-structured cross-trait covariance — fits exactly one shared variance across all traits, blended by a single relatedness signal parameter. It has no per-trait-distinct variance and no free cross-trait covariance parameter, with or without a similarity matrix supplied. The shape was confirmed by fitting: eight traits simulated with genuinely different slope variances all came back with one identical fitted variance. To be precise about what that fit does and does not show — it was run without colMat; the conclusion that supplying colMat cannot change the shape rests on there being no per-species variance or cross-species covariance parameter anywhere in the fitted object to begin with, not on a second fit. gllvm can therefore corroborate the shared-variance (common = TRUE) case of any source, but it structurally cannot serve as a reference for that source’s plain indep() or dep() cell — no combination of arguments changes that.

Why the spatial_ indep/dep/shared-variance cells are uncertain, not “none”

gllvm has two spatial mechanisms, and neither matches these cells cleanly. Its row.eff term adds the identical realised value to every trait’s linear predictor — every trait shares the literal same spatial field, which is a stricter model than “each trait draws its own value from a shared-variance spatially correlated field.” Its lvCor mechanism instead structures the ordination axis (the latent() column, where it does plausibly match). Whether another package’s mechanism fills the spatial_ indep/dep/ shared-variance gap was not checked. That is an honest gap in current scouting, not a checked absence — unlike the *_latent() row below, which was checked from both directions.

What cannot be externally validated

Read “cannot be externally validated” below as “no third-party package implements this model.” That is narrower than it may sound, and the distinction matters: the next section describes a Stan route that is, in principle, available for every item on this list, unbuilt as it currently is.

Stated plainly, not hedged:

• phylo_latent(), animal_latent(), and kernel_latent() — at any rank, with any response family. No package combines a source-structured relatedness matrix with reduced-rank ordination loadings the way these keywords do. This was checked from both directions: gllvm’s colMat cannot touch the ordination loadings at all (a separate, positive finding, not just an absence), and the one other package with a phylogenetic mechanism structures its trait-regression coefficients, not its ordination loadings. Simulation recovery is the only evidence in hand for these three keywords; see below for the unbuilt Stan route.
• Phylogenetic multinomial (a phylo_latent() term on multinomial category contrasts). This compounds two independent no-third-party-package cells — no package models phylo_latent(), and separately no package models gllvmTMB’s multinomial family. This is not a permanent gap: a hand-written Stan reference can encode the exact combined likelihood, so it is a build target rather than a case with no possible external reference (no such Stan model has been written). It is already partially implemented in gllvmTMB.
• meta_V(), the known-sampling-covariance meta-analysis keyword. None of the packages considered here model a known, externally supplied sampling covariance matrix as a variance component; the closest peers are general-purpose meta-analysis packages, which were out of scope for this comparison and were not evaluated.
• Any delta_* compound family, gengamma(), gamma_mix(), lognormal_mix(), student(), or censored_poisson()/ truncated_poisson(), combined with any source-structured mode. No package considered here has a peer for these families, regardless of the covariance structure chosen alongside them.

A stronger check that has not been built: the Stan route

Stan is a general-purpose probabilistic programming language, not a fixed catalogue of models. It is reachable from R through either rstan or cmdstanr, and a Stan program can hand-encode the exact likelihood of any cell in this document, including every “None” above — nothing about “no third-party package implements this” limits what can be written from scratch.

The check that matters here is not a fitted-estimate comparison, which is what every ✓ and ✓✓ cell above actually does. It is evaluating the log-likelihood at a fixed, shared parameter vector in both gllvmTMB’s TMB implementation and a hand-written Stan model of the same likelihood, and requiring agreement to machine precision. That is a stronger check than comparing fitted point estimates under a tolerance: there is nothing to estimate, so the estimator trap and the posterior-mean trap above cannot arise, and there is no tolerance band for a shared error to hide inside — the two numbers either match to machine precision or they do not. It is also the one route available for every cell in this document at once, independent of whether any third-party package exists at all.

That strength comes with a limitation that this whole page is otherwise warning against, and it deserves the same weight here. A Stan model that gllvmTMB’s own contributors write encodes their reading of the mathematics. It catches implementation errors — a template bug, a wrong parameterisation, a missing Jacobian, an indexing mistake — but not a shared conceptual misunderstanding: if the TMB template and a self-written Stan program both encode the same wrong reading of the model, they still agree to machine precision, and the check reports success on a shared error. A third-party package embodies someone else’s independent reading of the theory; a self-written reference does not. The mitigation is to derive the Stan model from the model’s published definition rather than from gllvmTMB’s own code, and preferably to have it written by someone who did not write the TMB template.

No Stan model has been written for any cell in this document. This section describes a route that exists, not evidence that has been gathered.

Current scope boundary

IN. Two cells have a shipped, passing comparator inside gllvmTMB’s own test suite: ordinary (no-prefix) latent(), checked against both glmmTMB::rr()+diag() and gllvm’s unconstrained ordination; and the phylogenetic shared-variance case, checked against glmmTMB::propto().

PARTIAL. Most of the grid falls here, at varying strength, and the distinction that matters is between a mechanism we have fitted and one we have only read about. Exactly one cell rests on a direct fit: a gllvm run with no similarity matrix supplied was checked and returns one shared variance across traits, matching the no-prefix shared-variance case. Everything else in this tier is documented rather than fitted — the gllvm-colMat shared-variance cases for phylo_, animal_ and kernel_ follow from gllvm’s documented behaviour (the diagnostic fit above was run without colMat), and spatial_latent()’s match to gllvm’s lvCor has not been fitted at all. Most of the phylo_/animal_/kernel_ indep() and dep() cells have only a documented, unfit MCMCglmm mechanism, carrying the posterior-mean caveat above. The spatial_ indep/dep/shared-variance cells are uncertain rather than checked.

PLANNED / NONE. phylo_latent(), animal_latent(), kernel_latent(), meta_V(), and every family named in the previous section have no third-party package reference and are not expected to gain one; nothing here is a near-term build target. Phylogenetic multinomial is the one exception: it has no third-party package peer either, but a hand-written Stan reference is a genuine build target for it (see above), so it is not treated the same as the permanent gaps in this paragraph. Building a gllvm-colMat-based comparator and labelling it a phylogenetic- or animal-latent() check would misrepresent what colMat estimates — that combination is deliberately not planned.

Where to go next

• Formula keyword grid: the full source-by-mode syntax this map is organised around.
• gllvmTMB vocabulary: definitions for loadings, Sigma, Psi, and the other terms used above.
• Fit diagnostics: whether a specific fitted model is numerically trustworthy, as distinct from whether the model class has external evidence behind it.