Validation evidence • hsquared

This article is the honest answer to one question: what does it mean when hsquared says a model is validated, and what is the actual evidence? It is written test-first. Every claim below points at a fixture, a test, or a study that lives in the repository, and nothing is advertised as working beyond what those checks demonstrate. One practical fact up front: the default univariate fit needs a local Julia engine (HSquared.jl via JuliaCall); the pure-R reference checks below run without it, and the “What runs where” section gives the full split. When in doubt about evidence atoms, the machine-readable table is validation_status(). What you may report is listed on Can I fit and report this? — this table is not the complete user list. Its columns are the capability, its status (one of planned, partial, or covered, defined just below), and the claim boundary attached to each.

Three-rung validation ladder from planned to partial to covered. A side note marks internal self-agreement as insufficient for covered status.

library(hsquared)
#> hsquared is an experimental 0.9.0 release (not production; CRAN availability is tracked separately). Default fitting needs Julia + HSquared.jl; preview with hs_control(engine = "validate").
#> What you may report is listed on Can I fit and report this? -- not in validation_status(). Intervals are not coverage-calibrated.
#> See https://itchyshin.github.io/hsquared/articles/current-limits.html
validation_status()
#> <hs_validation_status>
#>   validation: status table only; checks are run by tests and CI
#>   public claims: developer evidence table, not the user-facing list
#>   of what you may report; see Can I fit and report this?
#>   this table is not the full list of covered routes -- opt-in
#>   random_regression (k = 2) and direct_maternal are covered with no
#>   row yet; see the current-limits article
#>   shown below: `capability_label` (current wording); look rows up by
#>   the stable `capability` id that dated evidence records cite
#>                                                                                                                              capability_label
#>                                                                                                               tiny deterministic Ainv fixture
#>                                                                                                               Mrode9 pedigree Ainv comparator
#>                                                                                                       supplied-variance Henderson MME fixture
#>                                                                                                               sparse REML likelihood identity
#>                                                                                                         Mrode-style supplied-variance outputs
#>                                                                                                   experimental sparse REML estimator (opt-in)
#>                                                                                                 experimental repeatability estimator (opt-in)
#>  two-effect / arbitrary-N independent-effect estimator (opt-in; covered: common-env + (1|g) iid / A2=I; experimental: maternal / A2=pedigree)
#>                                                                   experimental supplied-relationship estimator (opt-in: genomic, single-step)
#>                                                       experimental SNP-BLUP marker-effect model (opt-in; supplied-variance or REML-estimated)
#>                                                                                      experimental multivariate REML estimator (default route)
#>                                                                                  univariate Gaussian animal-model fit (default path, AI-REML)
#>                                                                                       external published-REML recovery (gryphon, R reference)
#>                                                                                     known-truth DGP variance-component recovery (R reference)
#>                                                                                                             Mrode fitted animal-model outputs
#>                                                                                                                      ASReml comparison policy
#>                                                                                                          BLUPF90/DMU/WOMBAT comparison policy
#>                                                                                                                         XSim simulation truth
#>                                                                                                               genomic and QTL/eQTL validation
#>                                                                                                           GLLVM-style multivariate validation
#>                                                                                                                    CPU/GPU backend comparison
#>     phase  status
#>   Phase 1 partial
#>   Phase 1 partial
#>   Phase 1 partial
#>   Phase 1 partial
#>   Phase 1 partial
#>   Phase 1 partial
#>   Phase 2 partial
#>   Phase 2 covered
#>   Phase 5 covered
#>   Phase 5 partial
#>   Phase 3 covered
#>   Phase 1 covered
#>   Phase 1 covered
#>   Phase 1 covered
#>   Phase 1 planned
#>   Phase 1 planned
#>   Phase 1 planned
#>  Phase 5+ planned
#>  Phase 5+ planned
#>   Phase 6 planned
#>  Phase 7+ planned

What “validated” means here

hsquared treats validation as a product surface, not a footnote. A public capability is only advertised as working once evidence for it exists. The status ladder has three rungs, and the meanings are deliberately narrow:

planned — named in the roadmap, no implementation and no evidence. The parser may even reserve the syntax, but it errors as not implemented.
partial — implemented, but not yet through every requirement for a covered claim. A partial route may have limited evidence (for example, tests of shape or internal consistency, or a narrow evidence atom), but that evidence does not establish its full declared scope. It remains experimental; its routing controls and limitations are stated row by row.
covered — has declared scope evidence: documented known-truth recovery and/or a same-estimand external anchor at the stated tolerances. Evidence scope and routing are separate: a covered route can be default-routed or opt-in. Coverage is not a production guarantee or a nominal-interval guarantee; use Model status and route scope for each route’s controls and boundaries.

The gate is real. The v0.1 univariate Gaussian animal model could not move from partial to covered until a written, binding promotion rule (the predicate) was met and audited (see the V0.1 contract). That predicate explicitly rejects the easy mistakes: internal self-agreement (an R reference solve matching the engine at supplied variances), start-independence, and dense-vs-sparse-vs-pure-R agreement all prove the optimizer is reproducible on the same objective — they are not known-truth evidence and never satisfied the gate on their own.

A second discipline runs underneath: a test of tests. A fixture that asserts “the reference matches this pinned number” can be vacuous if the pinned number came from the same code, or if the tolerance is wide enough to accept anything. The negative-control suite feeds each real comparison a deliberately wrong value and proves the comparison rejects it. Validation that cannot fail is not validation.

The v0.1 evidence

The promoted capability is the univariate Gaussian animal model, y ~ fixed + animal(1 | id, pedigree = ped), Gaussian, REML, one record per animal. Here is the evidence behind that promotion, weakest-to-strongest, with the file each piece lives in.

Published-anchor recovery (gryphon birth weight)

The signed-off external anchor is the gryphon birth-weight univariate animal model from Wilson et al. (2010, J. Anim. Ecol. 79:13–26), with published REML estimates VA = 3.3954, VE = 3.8286, h² = 0.470 (data shipped in the CRAN package enhancer). Two independent recoveries are checked:

hsquared’s own pure-R REML reference recovers the published VA, VE, and h² to within the signed-off band (VA/VE within ~2%, h² within ~0.02).
the Julia engine (fit_sparse_reml and fit_ai_reml) recovers the same published values within that band, when supplied the published relationship matrix A_gryphon. The open Wilson et al. 2010 gryphon pedigree has no ancestral loops; the signed-off anchor still uses the supplied A. Engine-vs-pure-R agreement is to machine precision.

This is recovery within a maintainer-signed-off band, not a bit-exact match, and the gryphon population is teaching/simulated data. Tests: test-validation-fixtures.R (“hsquared’s R REML reference recovers the published gryphon estimates”; “the Julia engine recovers the published gryphon estimates via supplied A”).

Known-truth recovery (simulated DGP)

The strongest correctness evidence is a known-truth recovery study following the ADEMP framework (Morris, White & Crowther 2019) — a standard checklist for simulation studies: Aims, Data-generating mechanism, Estimands, Methods, Performance measures. Data are simulated from a univariate Gaussian animal model with known variance components over a clean, non-selfing generational pedigree (n = 420 animals; founders unrelated). True breeding values are drawn as u = sqrt(sigma_a2) * U' z with A = U'U, so Cov(u) = sigma_a2 * A. The estimator must recover the generating values and produce EBVs that track the true breeding values.

Recorded result (engine ai_reml, 120 replicates, fixed master seed 20240613, n = 420, all converged):

σ²ₐ, σ²ₑ, and h² are near-unbiased — 0 lies inside bias ± 2·MCSE for all three (e.g. h² bias −0.0049, MCSE 0.0073);
EBVs track the true breeding values, mean cor(EBV, true) ≈ 0.74 at h² = 0.4;
the engine matches the independent pure-R reference to machine precision.

Recovery holds across an h² grid (0.2 / 0.4 / 0.6, 100 reps per cell): all near-unbiased, with EBV accuracy rising with h² (≈ 0.60 / 0.74 / 0.83). The near-boundary cell (h² = 0.1) is reported honestly, not hidden: mild upward bias, 94% convergence, and 5% boundary pinning of σ²ₐ (the estimate hitting its lower bound of zero). Recovery also holds for a model with a fixed effect (y ~ x + animal): h² near-unbiased and the slope recovered (b_x ≈ 0.99 vs 1.0). The full study is reproducible evidence in data-raw/dgp-recovery-study.R (it is .Rbuildignored, not part of the build). A faster pure-R regression test guards a small-N case in CI (test-validation-fixtures.R, “REML recovers known variance components from a simulated DGP”).

External-package agreement

Two external R packages are used as comparators:

sommer is the signed-off V1-COMPARATORS agreement check. On the gryphon anchor it agrees with the pure-R reference within the signed-off band (variance components ~1–2% relative, h² ~0.01–0.02 absolute). The check is robust to sommer API churn — it skips, rather than fails, if the API differs.
pedigreemm provides a one-sided floor only: on a deterministic replicated dataset, hsquared’s REML solution is verified to be at least as good by REML log-likelihood (it reaches the true optimum; pedigreemm’s optimizer lands slightly off on these pedigree models, and cannot fit the saturated one-record-per-animal design at all — hence the replicated design). This is a lower bound, not two-sided parity.

Tests: test-validation-fixtures.R (“hsquared’s REML solution is at least as good as the pedigreemm comparator”).

For the executed side-by-side numbers (published vs sommer vs hsquared, plus the pedigreemm log-likelihood floor) with reproducing code, see the Benchmark: hsquared vs sommer and pedigreemm article.

Supplied-variance Henderson / Mrode-style fixtures

At supplied variance components, the engine and an independent R reference are pinned against textbook-style mixed-model-equation outputs:

a five-animal Henderson MME fixture (fixed effects, EBVs, fitted values, h²);
a twelve-animal Mrode-style fixture (Ainv, fixed effects, EBVs, fitted values, PEV, reliability, h², ML log-likelihood, and dense/sparse REML log-likelihood);
a tiny three-founder likelihood-identity fixture with closed-form ML and REML targets (-0.5(2·log2π + 3·log2 + log1.5 + 1) for REML), checked against Julia’s dense and sparse REML evaluators.

These solve the BLUP/MME at given variances; they do not estimate variance components. Tests: test-validation-fixtures.R, hs_*_validation_fixture() helpers in R/validation-fixtures.R.

Independent hand-built-MME anchors (non-circular)

The supplied-variance fixtures above have a known weakness: they pin numbers the package solver itself produced (hs_solve_henderson_mme_reference()), so a wrong solver would still pass. Two independent anchors close that circularity by building the MME by hand in lambda form (parameterised by the variance ratio λ = σ²ₑ/σ²ₐ rather than the variances themselves),

lambda = sigma_e2 / sigma_a2
C   = [[ X'X ,        X'Z          ],
       [ Z'X , Z'Z + Ainv * lambda ]]

solving it directly, and asserting the package solver matches that hand solve — never a value the solver produced:

diag(3) PEV/reliability anchor — a tiny intercept-only design with Ainv = I, anchoring the PEV (diag(solve(C)) * sigma_e2) and reliability (1 − PEV / (sigma_a2 · diag(A))) reference functions (test-pev-reliability-anchor.R);
twelve-animal pedigree anchor — the Mrode-style pedigree with genuine relatedness (off-diagonal Ainv entries, two-column X), anchoring fixed effects and EBVs (test-pedigree-mme-anchor.R).

These prove the solver agrees with an independent hand solve of the same MME. They do not yet pin the published Mrode textbook EBV digits — that is a flagged maintainer-confirm follow-up, not a silent gap.

Pedigree-inverse comparator (Mrode9 / nadiv)

Sparse Ainv construction is checked against an external reference: the engine’s pedigree_inverse() is compared to nadiv::makeAinv() on the Mrode example-9.1 pedigree shipped by nadiv, when both are available (test-mrode-validation.R). This is a pedigree-inverse comparator only, not a fitted-output claim.

The negative controls (tests of tests)

test-negative-control.R proves the comparisons above are discriminating. Each control reuses the same pure-R reference and the same tolerance as a real validation test, feeds it a deliberately wrong value, and asserts rejection:

a fixed-effect vector off by ±0.1, or an EBV off by 0.5, fails the 1e-12 / 1e-8 bands;
a REML log-likelihood off by 0.01 fails the 1e-10 band;
a gryphon VA off by 0.5 or h² off by 0.1 fails the 0.02 recovery band;
a variance estimate scaled by 1.5 fails the 5e-2 cross-check band;
a DGP statistic biased by 0.1 fails the 0.06 recovery band, and an EBV accuracy of 0.3 fails the 0.5 floor.

All negative controls are pure R, fast, and CI-runnable. If a real fixture ever became vacuous, its paired control would start failing.

What runs where: public CI vs locally

The split matters because most of the strong evidence needs a live engine that the public build host does not have.

Runs in public CI (no Julia):

the pure-R REML reference recovery checks (the small-N simulated DGP and, when enhancer is installed, the gryphon pure-R recovery);
the supplied-variance fixtures evaluated against the independent R references;
the two independent hand-built-MME anchors;
the full negative-control suite;
the pure-R leg of the independent-optimizer cross-check (convergence, finite positive estimates).

Runs only locally (needs a local Julia + HSquared.jl + JuliaCall): all live-engine tests are skip-guarded with hs_julia_bridge_available() and skip cleanly when the bridge is absent. These include:

the live-engine gryphon recovery via supplied A (fit_sparse_reml and fit_ai_reml);
the full engine DGP recovery study in data-raw/dgp-recovery-study.R;
the Julia pedigree_inverse() / henderson_mme() / dense-and-sparse REML agreement against the supplied-variance and Mrode9 fixtures;
the start-independence, sparse-vs-dense, and AI-REML-vs-sparse optimizer cross-checks.

This is why CI staying green is necessary but not sufficient: it confirms the R lane and the discriminating power of the comparisons, while the engine recovery evidence is established locally with the engine in place.

Honest boundaries

What the v0.1 promotion does not claim:

REML only. ML is not implemented. REML = FALSE is rejected; use REML = TRUE (the default). That is the same rule for the default fit and for engine = "validate". The supplied-variance paths (henderson_mme, snp_blup) solve at given variances rather than estimating them.
No nominally calibrated SEs/CIs. A validated, coverage-calibrated SE/interval — and confint()/vcov() — is deferred. summary() and the extractors may surface experimental, asymptotic REML SEs/intervals (heritability_interval(), variance_component_standard_errors(), heritability_standard_error(), repeatability_interval(), and multivariate covariance_standard_errors()), clearly labelled experimental and not coverage-calibrated. The named univariate-pedigree h2 delta/profile/bootstrap and sigma-a2 profile/bootstrap methods have directional-conservative evidence; sigma-a2 delta/Wald remains experimental-only. None is a nominal-coverage permission.
accuracy() is not a validated accuracy. It is the square root of reliability — a model-internal quantity on the dense validation path — not the realised cor(EBV, true) used as recovery evidence. PEV and reliability are reported on the dense validation path only and are labelled as such.
Covered set and remaining experimental set. Four opt-in models plus the default-routed t=2 unstructured multivariate path are covered at validation scale — each passed a pre-declared recovery gate plus a same-estimand REML comparator, mirroring the twin engine’s covered gate: (1) the common-environment two-effect leg (target = "two_effect" with common_env()); (2) the arbitrary-N independent (1|g) multi-effect model (target = "multi_effect"); (3) random regression k = 2 (target = "random_regression"); (4) the direct–maternal correlated 2×2 G model (target = "direct_maternal"); and (5) t=2 unstructured multivariate cbind() on the default path. Covered at validation scale means the experimental label is retained, dense, and not production — the warranty is validation-scale only, and intervals are asymptotic/uncalibrated. The remaining opt-in models — repeatability (target = "repeatability"), the maternal-genetic two-effect leg (target = "two_effect" with maternal_genetic()), single-step, and SNP-BLUP — stay experimental/partial. Narrow genomic GREML is covered at validation scale on the default route, but its production and interval debt remains. k≥3 and diagonal multivariate stay experimental.
Multivariate covered scope is t=2 unstructured. Known-truth recovery and full-unstructured sommer / blupf90+ comparators back that scope. They do not cover k≥3, diagonal genetic structure, or interval calibration. The optional diagonal-residual sommer check is not the unstructured-R0 claim.
No production / ASReml-parity / large-scale claim. The v0.1 default fit is validated only on small, clean pedigrees (the gryphon anchor and the simulated DGP study). Large, real, or heavily inbred pedigrees, and the engine half of boundary stability, are not yet validated — they remain post-v0.1 hardening. There is no ASReml / BLUPF90 / DMU / WOMBAT parity claim.

For the current, authoritative status of every capability — including which rows are planned, partial, and covered, and the exact claim boundary attached to each — call validation_status(). It is a status table only: it does not run checks, fit models, or promote anything.

For v0.1 the default-path covered rows are the univariate Gaussian animal-model fit and its two recovery anchors (the known-truth DGP study and the published gryphon estimate). Beyond v0.1, four opt-in models plus default-routed t=2 unstructured multivariate are also covered at validation scale (see the bullet above). Everything else is partial or planned.

status <- validation_status()
status[status$status == "covered", c("capability_label", "phase", "status")]
#> <hs_validation_status>
#>   validation: status table only; checks are run by tests and CI
#>   public claims: developer evidence table, not the user-facing list
#>   of what you may report; see Can I fit and report this?
#>   this table is not the full list of covered routes -- opt-in
#>   random_regression (k = 2) and direct_maternal are covered with no
#>   row yet; see the current-limits article
#>   shown below: `capability_label` (current wording); look rows up by
#>   the stable `capability` id that dated evidence records cite
#>                                                                                                                              capability_label
#>  two-effect / arbitrary-N independent-effect estimator (opt-in; covered: common-env + (1|g) iid / A2=I; experimental: maternal / A2=pedigree)
#>                                                                   experimental supplied-relationship estimator (opt-in: genomic, single-step)
#>                                                                                      experimental multivariate REML estimator (default route)
#>                                                                                  univariate Gaussian animal-model fit (default path, AI-REML)
#>                                                                                       external published-REML recovery (gryphon, R reference)
#>                                                                                     known-truth DGP variance-component recovery (R reference)
#>    phase  status
#>  Phase 2 covered
#>  Phase 5 covered
#>  Phase 3 covered
#>  Phase 1 covered
#>  Phase 1 covered
#>  Phase 1 covered

To actually fit the validated model end-to-end, including installing the HSquared.jl engine from GitHub, see the README Installation section and the Getting started article (Fit (needs the Julia engine)). This article is about the evidence, not the workflow.