Genomic Models

Experimental · opt-in · not the default fit

Engine APIs on this page are validation-scale and dense. The R genomic GREML surface stays partial/held. This page is not public activation and does not change public_covered_count.

HSquared.jl builds the genomic relationship engine on top of the same Henderson mixed-model-equation machinery as the pedigree animal model. The functions below are engine APIs (the Julia package's own functions). A narrow R-twin public-activation candidate now routes Gaussian REML models with one genomic random intercept from raw markers or a supplied Ginv into the existing supplied-precision solver. For marker input, the internal bridge construction freezes sample-frequency, unweighted VanRaden method 1 with K_lambda = G + 0.01I, and records engine-owned ID, marker, kernel, and precision fingerprints. The deterministic cross-twin fixture in test/fixtures/genomic_public_activation_target/ is construction and route- identity evidence. A later fail-closed boundary candidate matched an independent base-R oracle on all 240 sealed holdouts and corrected 30 classifications with 0 losses, but the conjunctive gate was still BOUNDARY_HOLDOUT_FAIL: one cell was 5.99x slower at p95 than the old path, above the frozen 3x cap. The 240 seeds are spent, the nine-cell recovery campaign did not launch, and this is not broad recovery, production-scale evidence, or permission to activate the R default. The R surface therefore remains partial/held with no Rose/G10 promotion basis. single_step() remains a separate experimental surface. The dense paths here are validation-scale only — they do not gain the sparse selected-inversion advantage.

Twin R article: Genomic prediction (genomic(), SNP-BLUP, single-step syntax and reporting fences).

Genomic relationship matrix G and its inverse

genomic_relationship_matrix builds the VanRaden (2008) G from a 0/1/2 (or dosage) marker matrix. genomic_relationship_inverse adds a ridge before inverting: a VanRaden G is rank-deficient — column-centering puts the all-ones vector in its null space (rank(G) ≤ n − 1), so a ridge is required.

using HSquared, LinearAlgebra

M = [0.0 1 2; 2 1 0; 1 1 1; 0 2 1]   # 4 individuals x 3 markers
G = genomic_relationship_matrix(M)
Ginv = genomic_relationship_inverse(G; ridge = 0.01)
round.(Ginv; digits = 3)

4×4 Matrix{Float64}:
 42.166  41.534  -8.15   24.45
 41.534  41.852  -8.307  24.92
 -8.15   -8.307  91.772  24.685
 24.45   24.92   24.685  25.945

GBLUP

fit_gblup solves the genomic animal model by placing Ginv in the same relationship slot the pedigree animal model uses, then reusing henderson_mme. At supplied variance components it returns the usual result object and works with every extractor (fixed_effects, breeding_values, heritability, …).

y = [10.0, 12.0, 11.0, 13.0]
X = ones(4, 1)
Z = Matrix(1.0I, 4, 4)
fit = fit_gblup(y, X, Z, Ginv, 1.0, 2.0)
(beta = fixed_effects(fit), gebv = round.(breeding_values(fit).values; digits = 4))

(beta = [11.5], gebv = [-0.5163, 0.3055, -0.1054, 0.3162])

The genomic variance components can be estimated by REML, by running the existing optimizers on a genomic spec:

M6 = [0.0 1 2; 2 1 0; 1 1 1; 0 2 2; 1 0 2; 2 1 1]
y6 = [10.0, 12.0, 11.0, 9.0, 13.0, 10.5]
Ginv6 = genomic_relationship_inverse(genomic_relationship_matrix(M6); ridge = 0.05)
spec = animal_model_spec(y6, ones(6, 1), Matrix(1.0I, 6, 6), Ginv6)
est = fit_animal_model(spec; target = :ai_reml)
est.variance_components

(sigma_a2 = 1.5583914976021154, sigma_e2 = 0.500615205536898)

SNP-BLUP and the GBLUP↔SNP-BLUP equivalence

fit_snp_blup fits marker effects: the centered markers W are the random-effect design, with an identity prior and per-marker variance σ²g / k (k = 2 Σ p(1 − p)). The implied breeding values gebv = W·â equal the GBLUP breeding values for the same data — the classic equivalence (verified to machine precision via the marginal covariance, the singular-G-safe route).

That equivalence applies when the two routes imply the same covariance. The activation candidate instead fits K_lambda = G + 0.01I; adding the ridge changes the covariance kernel, so this evidence does not establish an unqualified exact ridge-regularized GBLUP–SNP-BLUP equivalence.

snp = fit_snp_blup(y, X, M, 1.0, 2.0)
(marker_effects = round.(snp.marker_effects; digits = 4), gebv = round.(snp.gebv; digits = 4))

(marker_effects = [-0.0022, 0.4133, -0.4112], gebv = [-0.5128, 0.3051, -0.1039, 0.3116])

The random block is deliberately labelled marker_effects, not breeding values: on a SNP-BLUP spec the random effects are marker effects, and reusing the breeding_values / EBV vocabulary there would mislabel them.

Single-step H⁻¹

An internal helper HSquared._single_step_Hinv assembles the single-step relationship inverse

over the genotyped animals, where A₂₂⁻¹ = inv(A[g, g]) is the inverse of the submatrix of A (not the submatrix of A⁻¹ — the two differ). It is a validation-scale construction helper. The exported single_step_inverse, fit_single_step, and fit_single_step_reml wrappers expose the same dense relationship-precision path for tests and bridge targets. Its blending / τ / ω / ridge knobs are not comparator-validated. Inputs must be finite and symmetric; A and A⁻¹ and the returned precision must be positive definite; genotyped row indices must be unique; and require τ > 0, ω ≥ 0, blend weight in [0, 1], and non-negative ridge.

The supplied-Γ metafounder variant uses the same update with A replaced by the animal block of A^Γ:

metafounder_single_step_inverse, fit_metafounder_single_step, and fit_metafounder_single_step_reml are dense, validation-scale bridge primitives. Γ is supplied, not estimated. At Γ = 0, the helpers reduce to the ordinary pedigree single-step path. They do not add R-facing formula syntax or external BLUPF90 evidence by themselves.

Validation boundary

Engine evidence now available:

• VanRaden G on a hand-computed fixture (symmetric, PSD, pinned entries);

• regularized Ginv (defining identity (G + ridge·I)·Ginv ≈ I, ridge/PD guards);

• GBLUP against an independent dense MME (~1e-15) and reproducing pedigree BLUP when G = A (~1e-30);

• SNP-BLUP gebv = W·â equal to GBLUP (~1e-16) via the marginal V, for both n < m and n > m;

• genomic reliability / PEV / accuracy from the diag(inv(Ginv)) = diag(G) + ridge denominator, with selinv PEV matching the dense diagonal;

• genomic REML: AI-REML and NelderMead reach the same optimum, and a seeded simulation recovers the variance components; the supplied-Ginv estimator is covered through the preregistered 48-seed recovery gate and the historical same-precision blupf90+ comparator;

• the v0.7 activation fixture freezes the sample-frequency VanRaden-1 construction, K_lambda = G + 0.01I, Q_lambda = inv(K_lambda), provenance fingerprints, and marker-fed versus supplied-precision fit identity across repeated records, a nonconstant fixed effect, and unphenotyped genotyped IDs;

• the #49 genomic GBLUP/SNP-BLUP target is mirrored and consumed in the R twin (hsquared PR #84) by recomputing supplied-frequency G, Ginv, GBLUP MME, and SNP-BLUP route agreement without live Julia;

• single-step H⁻¹ reduction (H⁻¹ = A⁻¹ when G = A₂₂), locality, symmetry, and the A₂₂⁻¹ ≠ (A⁻¹)[g,g] distinctness guard;

• supplied-Γ H^Γ construction: reduction to ordinary single-step at Γ = 0 and equality to the manually built A^Γ + ordinary single-step path.

Still planned / coordinated:

• component/allocation profiling on already-open discovery data, output-equivalent performance work, a new preregistration, and a fresh untouched holdout before end-to-end recovery can resume; the narrow R genomic() candidate stays held and single_step() remains separate;

• a broader construction oracle beyond the completed independent base-R reconstruction of the exact candidate kernel. A fresh hash-pinned blupf90+ run now links the exact candidate precision to the supplied-precision estimator, but it is one point-estimate comparison rather than recovery or production evidence;

• sparse / APY G and GPU acceleration of the dense products;

• comparator-validated single-step blending defaults.