Nothing here is checked against itself. Each piece answers to an independent path:
| what | against what |
|---|---|
| A^-1 and F | tabular A by the classic recursion; A^-1 A = I
|
| sparse Cholesky, selected inverse |
solve() and the package’s own dense path |
| -2logL of the sparse MME | dense V form, a path with nothing in common |
| analytic score | central finite differences, ALL parameters |
| full fit | recovery of the components used to simulate the data |
| A22^-1 | inverse of the tabular-A block, plus the trap gate (block 22 of A^-1 != A22^-1) |
| single step | identity: blend 1 forces H^-1 == A^-1 through the whole fit |
| multi-trait | V form, finite differences on all parameters, and the collapse: both missing == row removed |
| APY | identity: the output is the exact inverse of the G that APY implies; core = everyone == exact |
| fixed-effect solutions | the published fixed effects of Examples 4.1, 5.1, 5.2 and 8.1 (by contrast), plus GLS and dense-MME rebuilds in plain R at 1e-6 |
| h2() and the share column | the denominators rebuilt by hand: Willham’s for direct-maternal; 4
var(sire) in a sire model and f var(sire) / 4 of the known maternal
grandsires with sire(mgs =); the residual at the harmonic
mean weight; NA with a warning when the SINGULAR verdict holds the
numerator (tests/testthat/test-estrutura-e6.R, simulated
data at given components) |
| structure of the data before the fit | a random term nested in a fixed effect or a single level stops in
every fitter that estimates and warns at given components; a
relationship term whose records all carry the unknown code 0 stops in
every fitter, in the ids rule; repeated records without
pe() and few levels warn
(tests/testthat/test-estrutura-e6.R; simulated data
only) |
| multi-trait, a level missing for one trait | sparse and dense V routes agree to 7.5e-12 after the per-trait drop; the components land where the single-trait fits land |
| Gibbs priors: the half-t, and the defaults | exactness, independent of N
(validation/gibbs_half_t_exactness.R,
tests/testthat/test-gibbs-meia-t.R): with no data the two
conditionals of the half-t of Huang and Wand (2013) recover the half-t
marginal of each standard deviation and the
(1 - r^2)^(nu/2 - 1) density of each correlation for d = 1,
2, 3 and nu = 2, 4 (KS sqrt(n) D at most 1.72 in 400 tests
of 1e5 draws, threshold 3), and the shape (nu + nl)/2 fails
the same test (at least 14.9 in 100 seeds); y, 10 y and 0.1 y and a
covariate of rn() times 10 or 0.1 give the same draws
rescaled (relative error at most 2.3e-13). On a direct-maternal group
(validation/gibbs_group_prior.R, 2650 animals over 4
generations, 30 replicates per design) with a planted correlation of
-0.8, the half-t moved the posterior mean of the correlation +0.428
posterior SD from REML (SE 0.032; flat +0.052), toward 0, and covered
var(maternal) in 25 of 30; at -0.3 and on rn()
with d = 2 and 3 it stayed inside or near flat’s distance from REML.
Scalar terms (validation/gibbs_scalar_prior.R, 2000
records, 5, 20 or 200 levels, planted 0 or 0.05, 30 data sets per cell):
no half-t chain stuck below 1e-6 var(y) at 0.05 (0 of 360;
"jeffreys" 35 of 360, and every chain at 0), coverage 171
of 180 ("jeffreys" 160 of 180), and slow mixing at a true 0
with 200 levels (median bulk ESS 116 to 128); at the size of a herd
(mode padrao, 2100 animals, 30 data sets per design,
"uniform_sd" on the residual) the half-t landed beside
REML, mean minus REML +0.04 and +0.03 posterior SD with one record per
animal and +0.09 / -0.07 / +0.09 with two records and pe(),
coverage 0.90 to 1.00 and a bulk ESS of 400 in 30 of 30. On scalar terms
with 5 and 20 levels (mode escalar_flat) no chain was stuck under either
residual prior and the coverage of 0.05 stayed at 113 of 120. These
numbers set the defaults: the half-t on scalar terms,
"flat" on covariance groups ("jeffreys" holds
chains at the boundary, and the half-t pulled the -0.8 correlation; the
flat runs are in validation/gibbs_group_prior.R),
"flat" on the residual (by a rule written before those
runs, whose coverage clause was read after the run as that of the scalar
term; the residual’s own coverage was 109 against 110 of 120 at a
planted 0.05), and the ASIS step on every estimated term of a Gaussian
trait (of a binary one, on the covariance groups). With
"uniform_sd" on the residual and the rest of the default on
the direct-maternal and rn() designs (mode padrao, 30
replicates each, 2 chains of 16 000 and of 8000): no chain of 240 held,
coverage 0.900 to 1.000, the -0.8 correlation +0.038 posterior SD (SE
0.009) from REML, a bulk ESS of 400 and an R-hat under 1.01 on every
component of every data set; the default itself, "flat" on
the residual, was not run on these two designs. The ASIS step, against
the same chains without it, raised the smallest bulk ESS per iteration
5.18 to 10.63 times there, and on a scalar term alone (mode asis) 2.12
times at a planted 0.05 and 3.47 at 0, with the posterior unchanged; in
that one-term model it costs as much as the rest of an iteration, and
per second it lost with 5 levels (0.82 and 0.93 times at 0) and on
random() with 20 levels at 0.05 (0.94).
tests/testthat/test-gibbs-padrao.R checks each part against
the same configuration named, and
tests/testthat/test-gibbs.R that the default lands where
REML lands. Binary trait
(validation/gibbs_probit_odegard.R,
validation/gibbs_probit_default.R, 30 data sets per design
unless stated): with one binary record per animal and 80 fixed classes
(10 data sets), every prior with flat fixed effects drifted to
var(animal) of about 1500 in 8 of 10, and the default,
half_t(df = 4) with prior_fixed = 1 and
px, left 4 of 10 at medians of 15 to 62, covering the
planted 0.25 in 3 of 10 (without px, the same posterior, 3
of 10 at 11 to 49: chains of 4 x 25 000 that do not cover it); in a sire
model with many small fixed classes prior_fixed = 1 without
px put var(sire) at 0.169 against 0.1 planted
(coverage 0.70; flat fixed effects 0.225 and 0.40); with large classes
at 2% incidence it recovered it (0.108, coverage 0.97). Simulated only:
no part of this row is validated on real data |
| sire / maternal-grandsire pedigree | exact: the pedigree expanded with a dummy dam per animal; the A^-1
printed for Example 15.2; the mixed pedigree (dam =)
against the same expansion only where the dam is missing, A to 1e-12 and
the fit identical; at scale, 20 replicates of 36 000 records recover
var(sire) = va / 4 within 1.6%
(validation/sire_mgs_recovery.R) |
| APY core by eigenvalues | the count by two routes (eigenvalues of G, singular values of Z);
"auto" == the same core passed by hand; the Lanczos
estimate against the exact count on both sides of the Gram matrix,
within 2% and 4 standard errors + 2 (20 probe seeds: at most 1.04% and
3.8 standard errors, validation/apy_core_lanczos_se.R); the
quadrature error, paired with the same probes on the exact eigenvectors
and averaged over v from 80 to 99.5%, under 0.1%; the standard error
against the spread of 20 independent probe sets, ratio within 0.7 to 1.4
(0.82 to 1.31 over eight draws at the gated levels); a probe that stops
at an invariant subspace gives the count of its probes on the exact
eigenvectors, and the test fails with that branch removed; at 10 000 x
10 000 over ten seeds, mean 5957.3 against 5953 exact (z = 1.1) and a
paired quadrature error of +2.3 counts, 0.04%
(validation/apy_core_lanczos.R) |
| h_inverse() | the single-step formula rebuilt in R, exact and with APY |
| Gibbs, the whole chain | with flat priors the posterior of the components is the REML
likelihood (Harville, 1974): four chains match a 2-D quadrature of the
-2logL at z = 0.15 and 0.08 (validation/gibbs_harville.R);
mcmc_diagnostics() equals posterior 1.7.0 to 1e-10 (R-hat,
bulk/tail ESS, MCSE), and its verdict is calibrated on simulated chains:
never not mixed on stationary ones; not mixed
83% and 100% of the time with one of 4 chains shifted 0.5 and 1 sd;
burn-in too short every time for a 50 sd start with 4
chains (validation/mcmc_diagnostics_calibration.R). At 20
000 levels and 4 chains the diagnostics take 0.22 to 0.23 s (3 runs;
validation/mcmc_diagnostics_scale.R, iid levels); serial
and parallel chains identical |
| threads | the tiled tail, the dense inverses, the products and the selected
inverse against chol(), solve() and the
formulas; the same bits with 1 and 4 threads, and the whole suite passes
at both |
| metafounders with genotypes | H(Gamma) rebuilt from the dense A(Gamma) and G05 at 1e-8;
snp_blup() == model() at the same theta; 20
simulated replicates of two bases reduce the trend bias from -0.141 to
-0.081 (validation/hgamma_recovery.R) |
estimate_gamma() |
one pseudo-EM step is the metafounder block of H(Gamma) built
densely; GLS shows its published bias; "ml" against the
identity A(gamma) = (1 - gamma/2) A + gamma 11’ and the maximum of the
dense likelihood; the engine version against the previous R version:
pseudo-EM and GLS within 1e-10 with the same iteration count, with
genotypes in the last generation only, with genotyped parents and with
every animal genotyped; ML summaries within 1e-12 and gamma within 1e-6
with A22^-1 dense, sparse and with no non-genotyped animal; A22^-1 X
with 70 columns against the dense A22; a metafounder with no genotyped
descendant gives the Gamma of the problem without it |
snp_effects() |
G*^-1 u read back from the exported H^-1, ((H^-1 - A^-1) u)_g +
A22^-1 u_g, with ungenotyped ancestors, on each route (exact and APY)
with and without metafounders, and, on the exact route, on a sire /
maternal-grandsire pedigree; APY also against the APY inverse built
densely in R; with metafounders, the marker effects of
snp_blup() at the same theta; the previous R version within
1e-9; 1 vs 4 threads identical; a different blend and other genotypes
are refused |
| A22 by Colleau, the sparse Schur | H^-1 against the formula with A22 from the dense A, with inbreeding
and with a sire-MGS pedigree; a22_inverse() against
solve(); the APY route with core = everyone gives the exact
H^-1 to 1e-12; with two metafounders and a full Gamma, core = everyone
gives the exact H(Gamma) to 1e-10 and a partial core the APY inverse
written in R on the dense G* to 1e-8 |
pegs() |
at fixed variances the exact multivariate ridge, a dense
mk x mk solve; the structures are identities where they
must be; recovery with 2000 animals x 2000 markers x 3 traits; the Julia
reference gives the same numbers on the same data |
sire(sire, mgs =) |
-2logL and BLUP against the dense GLS with the incidence built by hand |
snp_blup(), any structure |
the dense single-step solve built in R with H^-1 from
G* without the affine step: breeding values and marker
effects to 1e-6 of their SD with one animal term, direct-maternal in one
group, direct and maternal in separate groups and a reaction norm;
rpg near 1 falls back to the pedigree BLUP; a planted QTL
comes out on top |
| genotype storage | double, integer and raw matrices give the same G, H^-1 (exact and
APY), APY core (both routes), Gamma (pseudo-EM, GLS, ML), fit,
backsolved SNP effects (exact and APY) and ssSNPBLUP, bit for bit, with
missing values in the matrix; read_blupf90_snp() returns
the matrix that was written, subsets by id |
| residuals and leverage | the dense R P y, GLS fixed part and
diag(R P R) in plain R for model() (weights,
pe, rn, mgs, single step) and model_mt(). Not validated on
real data |
| fixed-effect design (sparse) | the sparse rank keeps the columns of the dense Gram-Schmidt it
replaced on 150 random designs (crossed, nested and disconnected
classes, levels with no used record, collinear and nested covariates,
weights), on edge columns (tiny, NaN, Inf, nearly collinear), and on
unscaled polynomials and near copies. In near-singular designs every
decision is certified or the rank falls back to that Gram-Schmidt, and
it never keeps more columns than records. dropped_x of
model() (with and without weights) and
model_mt() (per trait) equals a dense reference rebuilt in
R. The bits are the same with 1 and 4 threads. With the rest of the
engine unchanged, a check found the results bit-identical to the dense
build over the 3 085 values returned by the fits of every test file but
this change’s own (validation/fixed_design_identity.R).
Scale (memory, not exactness): 200 000 records x 701 candidate columns,
fixed part 145-153 MB of peak working set (232 MB committed), three sets
of 3 replicates, with dropped_x checked against the table
(validation/fixed_design_memory.R) |
| missing labels and values | the rule has two implementations, the engine’s and R’s. One battery
of values (numeric NA, NaN, Inf,
-Inf; blank, tab and line-break text, ” NaN “) gives the
same types, the same count and the same first row through both, in every
fitter. A missing trait drops only its record, with the -2logL of the
data without that row |
unit, shift and boundaries in model_mt()
|
y, 100 y and 1e4 y give the
same SINGULAR verdict and the same standard errors up to the factor in
model() and model_mt(); one trait in other
units gives the same fit up to the scale factors and the Jacobian; a
constant added to y leaves the REML -2logL alone; a
covariance boundary that is not the optimum is not certified, and a
certified point on a flat face is polished to the face optimum
(tests/testthat/test-unidade-partida-fronteira.R). At 24
000 records the -2logL agrees to 6e-15 (relative) under a shift of up to
1000 standard deviations
(validation/unit_shift_start_scale.R). In simulated sweeps
of 368 fits from four starts per data set, run during development and
not kept in validation/, no certified point sits more than
9.4e-5 above the best of its starts, and one known pair of starts still
certifies points 1.7e-4 apart. The real fits that motivated the change
were not rerun |
| ids of the data against the pedigree (founders, unknown 0, pruning, spellings, 2^53, blocks) |
tests/testthat/test-estrutura-e1.R (20 gates);
validation/pedigree_ids_scale.R: pruning and founders exact
(difference 0 in -2logL, EBV and accuracy) at 20 200 records with a 125
200-line pedigree; the rule alone takes 0.03 s, and the fit with the
national pedigree takes the time of the fit pruned by hand (median of 3:
11.97 s against 10.90 s, replicates 11.1-14.5 against 10.8-12.4 s on a
shared machine); simulated only |
The tests in tests/testthat run these comparisons on
every build, so a change that breaks one of the identities cannot pass
quietly. simulate_breeding() is what they are built on: it
generates the pedigree, the phenotypes and the genotypes together, so
every check has the truth beside it.
The works cited here by author and year are listed in full in the references.