Reference
HSquared.HSControl Type
HSControl(; backend = AutoBackend(), accelerator = :auto, precision = Float64,
save = :minimal, save_fitted = false, save_residuals = false,
save_design = false, save_factorization = false,
disk_cache = false)Store planned engine controls for future HSquared.jl model calls.
Phase 1 validates and records these controls so the R and Julia twins can agree on names before production backend dispatch exists. CPU is the trusted always-available path; accelerator names are future optional-extension metadata.
sourceHSquared.CPUBackend Type
CPUBackend()Select the trusted CPU execution path.
The CPU backend is the only always-available backend. Phase 1 records the backend choice only; it does not route production fitting through a backend dispatcher yet.
sourceHSquared.ThreadsBackend Type
ThreadsBackend()Select the planned multi-threaded CPU execution path.
This is metadata only. No threaded execution path is implemented yet.
sourceHSquared.CUDABackend Type
CUDABackend()Select the planned CUDA execution path.
CUDA is not a dependency and no GPU work is executed yet.
sourceHSquared.AMDGPUBackend Type
AMDGPUBackend()Select the planned AMDGPU/ROCm execution path.
AMDGPU.jl is not a dependency and no GPU work is executed yet.
sourceHSquared.MetalBackend Type
MetalBackend()Select the planned Apple Metal execution path.
Metal.jl is not a dependency and no GPU work is executed yet.
sourceHSquared.OneAPIBackend Type
OneAPIBackend()Select the planned Intel oneAPI execution path.
oneAPI.jl is not a dependency and no accelerator work is executed yet.
sourceHSquared.AutoBackend Type
AutoBackend()Let HSquared.jl choose an execution backend in a future implementation.
Until backend dispatch exists, AutoBackend() is control metadata only. CPU is the trusted fallback.
HSquared.BackendInfoRow Type
BackendInfoRowTyped status row returned by backend_info.
Rows describe accepted backend names and whether a backend is execution-ready. In Phase 1 all rows are selectable control metadata with execution_available == false and status == :planned.
HSquared.backend_info Function
backend_info(control = HSControl())Return typed status rows for planned compute backends.
The rows mirror the R twin's backend_info() shape: backend, accelerator, requested, selectable, execution_available, status, and note. In Phase 1 all rows are accepted control metadata, all have execution_available == false, and all have status == :planned.
HSquared.FormulaStatusRow Type
FormulaStatusRowTyped engine-term status row returned by formula_status.
Rows use the same columns as the R twin's formula_status(): term, category, phase, syntax_status, fitting_status, and current_behavior.
HSquared.formula_status Function
formula_status()Return a Julia-local diagnostic table for the engine-facing model-term vocabulary. The R package owns public formula parsing and its formula_status() table has a broader inventory; this table is not a mirror of R grammar status.
Julia does not parse formulas. syntax_status and fitting_status describe only the selected Julia marker or engine route represented by each row, not the R formula parser or its full status inventory. engine bridge marks the experimental low-level animal payload route; reserved and planned mark Julia placeholder or roadmap rows. Bounded FA and genetic GLLVM expert-control fits are separate routes and are not activated by cov = fa() or by a formula term. Consult validation_status() and the capability records for those fit boundaries.
HSquared.ValidationStatusRow Type
ValidationStatusRowTyped validation-status row returned by validation_status.
Each row records the validation item, current evidence, missing evidence, and the allowed public-claim boundary.
sourceHSquared.validation_status Function
validation_status()Return the current validation-evidence ladder for HSquared.jl.
This is a diagnostic table only. It does not run comparator packages, fit models, or promote any planned capability.
sourceHSquared.planned_model_terms Function
planned_model_terms()Return the planned model-term names reserved by HSquared.jl.
These names mirror the R twin's inert formula markers. They reserve formula vocabulary only; calling a term does not fit a model. Some corresponding direct Julia utilities are available through separate entry points.
sourceHSquared.planned_genomic_qtl_terms Function
planned_genomic_qtl_terms()Return the planned genomic, single-step, marker, and QTL term names.
sourceHSquared.planned_quantgen_terms Function
planned_quantgen_terms()Return the planned standard quantitative-genetic, parental, inheritance, and custom-kernel term names.
sourceHSquared.genomic Function
genomic(args...; kwargs...)Planned genomic relationship model term.
Mirrors the R genomic() formula marker. This function is intentionally not implemented yet.
HSquared.single_step Function
single_step(args...; kwargs...)Planned single-step pedigree/genomic relationship model term.
Mirrors the R single_step() formula marker. This function is intentionally not implemented yet.
HSquared.markers Function
markers(args...; kwargs...)Planned marker-effect model term.
Mirrors the R markers() formula marker. This function is intentionally not implemented yet.
HSquared.marker_scan Function
marker_scan(args...; kwargs...)Planned marker-scan / GWAS model term.
Mirrors the R marker_scan() formula marker. This function is intentionally not implemented yet.
HSquared.qtl_scan Function
qtl_scan(args...; kwargs...)Planned QTL interval-scan model term.
Mirrors the R qtl_scan() formula marker. This function is intentionally not implemented yet.
HSquared.permanent Function
permanent(args...; kwargs...)Planned permanent-environment model term.
Mirrors the R permanent() formula marker. This function is intentionally not implemented yet.
HSquared.common_env Function
common_env(args...; kwargs...)Planned common-environment model term.
Mirrors the R common_env() formula marker. This function is intentionally not implemented yet.
HSquared.maternal_genetic Function
maternal_genetic(args...; kwargs...)Planned maternal-genetic model term (engine-internal formula marker).
Mirrors the R maternal_genetic() formula marker. This stub is intentionally not implemented as a standalone engine path. Note: the correlated 2×2 direct– maternal capability IS covered at validation scale (opt-in) via fit_direct_maternal_reml / R target="direct_maternal" — this stub is only the engine-internal formula-reservation marker, not the fitting entry point.
HSquared.maternal_env Function
maternal_env(args...; kwargs...)Planned maternal-environment model term.
Mirrors the R maternal_env() formula marker. This function is intentionally not implemented yet.
HSquared.paternal_genetic Function
paternal_genetic(args...; kwargs...)Planned paternal-genetic model term.
Mirrors the R paternal_genetic() formula marker. This function is intentionally not implemented yet.
HSquared.paternal_env Function
paternal_env(args...; kwargs...)Planned paternal-environment model term.
Mirrors the R paternal_env() formula marker. This function is intentionally not implemented yet.
HSquared.cytoplasmic Function
cytoplasmic(args...; kwargs...)Planned cytoplasmic-inheritance model term.
Mirrors the R cytoplasmic() formula marker. This function is intentionally not implemented yet.
HSquared.imprinting Function
imprinting(args...; kwargs...)Planned parent-of-origin / imprinting model term.
Mirrors the R imprinting() formula marker. This function is intentionally not implemented yet.
HSquared.dominance Function
dominance(args...; kwargs...)Planned dominance relationship model term.
Mirrors the R dominance() formula marker. This function is intentionally not implemented yet.
HSquared.epistasis Function
epistasis(args...; kwargs...)Planned epistatic relationship model term.
Mirrors the R epistasis() formula marker. This function is intentionally not implemented yet.
HSquared.relmat Function
relmat(args...; kwargs...)Planned custom relationship-matrix model term.
Mirrors the R relmat() formula marker. This function is intentionally not implemented yet.
HSquared.precision Function
precision(args...; kwargs...)Planned custom precision-matrix model term.
Mirrors the R precision() formula marker. This function is intentionally not implemented yet.
HSquared.GaussianFamily Type
GaussianFamily()Marker for the Phase 1 Gaussian animal-model family.
sourceHSquared.HSData Type
HSDataIn-memory container for matched phenotypes, pedigree, genotypes, expression, marker annotation, and environmental data.
This is a conservative Phase 1 mirror of the R hs_data() contract. It stores the supplied objects and records exact-ID overlap diagnostics, but it does not construct relationship matrices, read file-backed data, or fit a model.
HSquared.HSDataIDMap Type
HSDataIDMapConservative ID-overlap summary for HSData.
IDs are matched exactly with Julia equality. The first slice deliberately does not coerce 1 to "1" or otherwise normalize values across sources.
HSquared.HSMarkerMapSpec Type
HSMarkerMapSpecValidated marker-map metadata stored by HSData.
This is metadata hygiene only. It does not parse genotype files, impute marker dosages, construct genomic relationships, or run marker scans.
sourceHSquared.HSGenotypeMarkerSpec Type
HSGenotypeMarkerSpecValidated alignment between genotype marker columns and a marker map.
sourceHSquared.HSAnnotationSpec Type
HSAnnotationSpecValidated expression-feature annotation metadata stored by HSData.
This is metadata hygiene only. It does not join annotation covariates into model matrices, fit eQTL or omics models, or run GLLVM workflows.
sourceHSquared.HSEnvironmentSpec Type
HSEnvironmentSpecValidated environment-key metadata stored by HSData.
This is metadata hygiene only. It does not add environmental model terms, join environment covariates into design matrices, or fit multi-environment models.
sourceHSquared.HSDataIDOverlapRow Type
HSDataIDOverlapRowOne ID-overlap count returned by data_status.
HSquared.HSDataPedigreeStatusRow Type
HSDataPedigreeStatusRowOne pedigree diagnostic returned by data_status.
HSquared.HSDataMarkerStatusRow Type
HSDataMarkerStatusRowOne marker-map or genotype-marker diagnostic returned by data_status.
HSquared.HSDataGenotypeStatusRow Type
HSDataGenotypeStatusRowOne genotype-component diagnostic returned by data_status.
HSquared.HSDataExpressionStatusRow Type
HSDataExpressionStatusRowOne expression-component diagnostic returned by data_status.
HSquared.HSDataAnnotationStatusRow Type
HSDataAnnotationStatusRowOne annotation-metadata diagnostic returned by data_status.
HSquared.HSDataEnvironmentStatusRow Type
HSDataEnvironmentStatusRowOne environment-metadata diagnostic returned by data_status.
HSquared.HSDataStatus Type
HSDataStatusDiagnostic container returned by data_status.
This mirrors the R twin's data_status() surface for component presence, ID-overlap counts, pedigree status, marker-map/genotype-marker alignment status, genotype status, expression status, expression-feature annotation status, and environment-key metadata status. It is diagnostic only and does not build model specifications, join covariates, or construct relationship matrices.
HSquared.id_map Function
id_map(data)Return the conservative ID-overlap map stored in an HSData.
HSquared.data_status Function
data_status(data::HSData)Return component, ID-overlap, marker-alignment, genotype, expression, annotation-feature, and environment-key diagnostics for an HSData object.
This mirrors the R twin's data_status() helper. It does not parse genotype files, construct genomic relationships, join expression, annotation, or environment covariates, build bridge payloads, or fit models.
HSquared.Pedigree Type
Pedigree(ids, sire, dam, original_order; missing_values = DEFAULT_UNKNOWN_PARENT_VALUES,
allow_selfing = false)Normalized pedigree with parent references encoded as row indices.
sire and dam use 0 for unknown parents. Rows are topologically sorted so known parents always occur before their offspring. Direct construction checks unique IDs, matching vector lengths, earlier-row parent indices, and that original_order is a permutation of the input rows. IDs cannot use values in missing_values, which defaults to missing, nothing, "", "0", and 0. Pass the same missing_values used to normalize custom parent markers. Same-parent sire/dam indices are rejected unless allow_selfing = true, matching normalize_pedigree. Prefer that function for raw parent labels.
HSquared.normalize_pedigree Function
normalize_pedigree(ids, sire, dam; missing_values = ..., allow_selfing = false)Validate, recode, and topologically sort an animal pedigree.
ids, sire, and dam must have the same length. Parent values must be either unknown-parent markers or values present in ids. By default, missing, nothing, "", "0", and 0 are treated as unknown-parent markers.
The returned Pedigree stores parents as integer row indices, with 0 for unknown parents.
allow_selfing = true permits self-fertilization (sire == dam), a Phase-3 non-standard-inheritance feature: the additive-relationship recursion and the Henderson Ainv rules already handle a self (a selfed offspring of a non-inbred parent has inbreeding F = 1/2, A_ii = 3/2). It is rejected by default to preserve the sexual-pedigree contract; an individual still cannot be its own parent.
HSquared.inbreeding_coefficients Function
inbreeding_coefficients(pedigree; max_relationship_cache = 10_000)
inbreeding_coefficients(ids, sire, dam; max_relationship_cache = 10_000)Return the inbreeding coefficient for each row of a normalized pedigree.
Inbreeding is computed by the Meuwissen & Luo (1992) method (_meuwissen_luo_inbreeding) without forming the dense relationship matrix. Runtime depends on the ancestry visited and heap operations. max_relationship_cache is accepted for signature compatibility but does not bound this path; it still governs the dense _numerator_relationship used by additive_relationship and single-step A₂₂, and as the inbreeding validation oracle.
HSquared.pedigree_inverse Function
pedigree_inverse(pedigree; max_relationship_cache = 10_000)
pedigree_inverse(ids, sire, dam; max_relationship_cache = 10_000)Construct the sparse inverse additive relationship matrix Ainv.
The implementation follows Henderson's direct inverse contribution pattern: each animal adds a scaled outer product of [1, -1/2, -1/2] over itself and its known parents. The matrix is returned as a SparseMatrixCSC{Float64,Int}.
HSquared.mendelian_sampling_variances Function
mendelian_sampling_variances(pedigree; max_relationship_cache = 10_000)
mendelian_sampling_variances(ids, sire, dam; max_relationship_cache = 10_000)Per-individual Mendelian sampling variances d_i — the diagonal of D in the decomposition A = T·D·Tᵀ of the additive relationship matrix (T unit lower-triangular). d_i is the variance of an individual's breeding value given its parents: 1 for a founder, 0.75 − 0.25·F_parent with one known parent, and 0.5 − 0.25·(F_sire + F_dam) with both known. Returned in pedigree.ids (topologically sorted) order; det(A) = ∏_i d_i.
The reciprocals (1/d_i) scale the Henderson pedigree_inverse contributions. The variances are also used in gene-dropping and within-family accuracy.
HSquared.additive_relationship Function
additive_relationship(pedigree; max_relationship_cache = 10_000)
additive_relationship(ids, sire, dam; max_relationship_cache = 10_000)Dense additive (numerator) relationship matrix A for a pedigree — the companion of the sparse pedigree_inverse (A = inv(Ainv)). A[i, i] = 1 + F_i and A[i, j] is twice the coancestry of i and j. Returned in pedigree.ids (topologically sorted) order.
Validation-scale and dense (bounded by max_relationship_cache); for production use the sparse pedigree_inverse. This is the additive companion of the exported dominance_relationship, cytoplasmic_relationship, and clonal_relationship matrices.
HSquared.dominance_relationship Function
dominance_relationship(pedigree)
dominance_relationship(ids, sire, dam)Dense dominance relationship matrix D (Cockerham). For animals x, y with both parents known — sires sx, sy, dams dx, dy — D[x, y] = ¼·(A[sx, sy]·A[dx, dy] + A[sx, dy]·A[dx, sy]), where A is the additive numerator relationship. The diagonal is set to 1, and any pair in which either animal has an unknown parent has D = 0. Returned in pedigree.ids order, aligned with pedigree_inverse.
Experimental Phase 3 primitive, validation-scale and dense — D is the relationship for a dominance random effect. The off-diagonal formula is general; the unit diagonal and the absence of dominance-inbreeding corrections hold for non-inbred parents (the standard textbook case). With unrelated, non-inbred parents, full sibs have D = 1/4; half sibs and parent–offspring have D = 0 under the corresponding unrelated-founder pedigree conditions.
HSquared.epistatic_relationship Function
epistatic_relationship(pedigree; kind = :additive_additive)
epistatic_relationship(ids, sire, dam; kind = :additive_additive)Dense epistatic relationship matrix as a Hadamard (element-wise) product of the additive A and dominance D relationship matrices (Henderson 1985):
kind = :additive_additive→A ∘ A,kind = :additive_dominance→A ∘ D,kind = :dominance_dominance→D ∘ D.
These give the orthogonal additive×additive, additive×dominance, and dominance×dominance epistatic relationship matrices used for epistatic variance components. Full sibs from unrelated, non-inbred parents have additive×additive = 1/4. Returned in pedigree.ids order. Experimental Phase 3 primitive, validation-scale and dense; inherits the dominance_relationship non-inbred-parent assumption. No public model-spec.
HSquared.metafounder_relationship Function
metafounder_relationship(pedigree, group_of, Gamma; max_relationship_cache = 10_000)
metafounder_relationship(ids, sire, dam, group_of, Gamma; ...)Dense metafounder-augmented additive relationship matrix A^Γ (animal × animal, in pedigree.ids order) for a SUPPLIED m×m symmetric positive-semidefinite metafounder covariance Γ and a metafounder assignment group_of. group_of is aligned to pedigree.ids like clonal_relationship's clone_of: entry i is the metafounder-group label of animal i's unknown parent slot(s) (both unknown parents of an animal share one group this slice); entries for animals with both parents known are ignored. Distinct labels resolve to Γ's rows 1..m in stable first-appearance order in the normalized pedigree.ids.
In the raw-array overload, group_of is supplied in the same original row order as ids; the method carries it through the pedigree's topological sort. Supply Γ in the group order after this sort, which can differ from the first group appearance in the caller's original rows; Γ is not itself reordered.
A^Γ is the existing tabular recursion with the leading Γ block seeded: a metafounder's self-relationship is Γ[i,i] (so its inbreeding F = Γ[i,i] − 1, which is NEGATIVE for Γ[i,i] < 1), and an animal's two unknown parents drawn from the same group m make it inbred with A_ii = 1 + Γ[m,m]/2 and two such founders related by Γ[m,m]. At Γ = 0 the result equals additive_relationship exactly.
Γ is SUPPLIED, NOT estimated. EXPERIMENTAL, validation-scale and dense (bounded by max_relationship_cache); no external (BLUPF90) comparator, no R-facing metafounder model-spec or bridge payload, and no covered claim. The separate single-step H^Γ utilities and supplied-variance wrappers are also experimental. Scale convention: Γ[m,m] is on the relationship scale (Γ[m,m] = 0 ⇒ classical founder), not an allele-frequency parameterization.
HSquared.metafounder_relationship_inverse Function
metafounder_relationship_inverse(pedigree, group_of, Gamma; max_relationship_cache = 10_000)
metafounder_relationship_inverse(ids, sire, dam, group_of, Gamma; ...)DESCRIPTIVE dense inverse of the ANIMAL block of A^Γ, i.e. inv(A^Γ_animals). This is distinct from metafounder_inverse: it is NOT the animal sub-block of the combined [metafounders; animals] precision (conflating the two corrupts either the MME or the descriptive inverse). At Γ = 0 it equals Matrix(pedigree_inverse(pedigree)). Validation-scale and dense; Γ supplied, not estimated. For the raw-array overload, group_of follows the original ids order and is reordered with the pedigree.
HSquared.metafounder_inverse Function
metafounder_inverse(pedigree, group_of, Gamma; max_relationship_cache = 10_000)
metafounder_inverse(ids, sire, dam, group_of, Gamma; ...)MME-ready sparse inverse of the COMBINED metafounder-augmented relationship over [metafounders 1..m; animals m+1..m+n] (dimension (m+n)×(m+n)): the metafounder block is seeded with inv(Γ) and each animal adds the existing Henderson [1, −½, −½]/d_k outer product over its (animal, sire, dam) combined indices, where a parent index may be a metafounder column. The Mendelian sampling variance d_k = ½ − ¼(F_s + F_d) uses the combined inbreeding (metafounder F = Γ − 1), so d_k may EXCEED ½ (heterozygote excess) — this is correct and not clamped.
Requires Γ symmetric POSITIVE DEFINITE (it inverts Γ). The numerical check uses s = max(1, maximum(abs, Γ)): symmetry tolerance 1e-10*s, and smallest eigenvalue strictly greater than 1e-12*s for this combined inverse. Thus even a well-conditioned Γ can be refused at a very small absolute scale. The descriptive relationship path instead permits eigenvalues down to -1e-10*s as roundoff. The leading m rows/columns are the metafounder solutions; the animal block of this combined inverse is not inv(A^Γ_animals) (see metafounder_relationship_inverse). A_combined · this ≈ I is the round-trip correctness anchor. Returned as a SparseMatrixCSC{Float64,Int}, mirroring pedigree_inverse. EXPERIMENTAL, validation-scale (still forms the dense A^Γ to get d_k, bounded by max_relationship_cache); not wired into henderson_mme this slice (the extra metafounder levels need the R bridge contract). The separate animal-only metafounder_animal_model uses the descriptive animal-block inverse; it does not solve the combined metafounder-effect system. For the raw-array overload, group_of follows the original ids order and is reordered with the pedigree. Γ supplied.
HSquared.metafounder_inbreeding Function
metafounder_inbreeding(pedigree, group_of, Gamma; max_relationship_cache = 10_000)
metafounder_inbreeding(ids, sire, dam, group_of, Gamma; ...)Metafounder-aware inbreeding coefficient F_i = A^Γ[i,i] − 1 for each animal (in pedigree.ids order). Under metafounders an animal can have F < 0 (heterozygote excess) when its parents' additive relationship is negative: F_i = A^Γ[sire,dam]/2 in the augmented pedigree. The negative value Γ[g,g] - 1 for Γ[g,g] < 1 belongs to the leading metafounder, rather than automatically to its animals. For example, founders assigned to one group with Γ[g,g] = 0.4 have animal F = 0.2, while the metafounder has F = -0.6. Γ supplied, not estimated. The raw-array overload takes group_of in the original ids order and returns coefficients in normalized pedigree.ids order.
HSquared.metafounder_animal_model Function
metafounder_animal_model(y, X, Z, pedigree, group_of, Gamma, sigma_a2, sigma_e2;
ids = pedigree.ids)Supplied-variance Gaussian animal-model BLUP under a metafounder-augmented relationship A^Γ (#53, Legarra et al. 2015). Builds the descriptive animal-only metafounder precision inv(A^Γ) via metafounder_relationship_inverse and solves the standard Henderson MME (henderson_mme) at supplied variance components, returning the HendersonMMEResult. At Γ = 0 this reduces EXACTLY to the classical animal model (metafounder_relationship_inverse → pedigree_inverse), so the fixed effects and EBVs match henderson_mme with pedigree_inverse.
SUPPLIED-variance and SUPPLIED-Γ only — neither Γ nor the variance components is estimated. This is the animal-only BLUP under A^Γ; the combined system with explicit metafounder effects (metafounder_inverse) is a separate path, and there is no R-facing model-spec or bridge payload.
HSquared.metafounder_single_step_inverse Function
metafounder_single_step_inverse(pedigree, group_of, Gamma, G, genotyped_rows; ...)
metafounder_single_step_inverse(ids, sire, dam, group_of, Gamma, G, genotyped_rows; ...)Dense validation-scale single-step relationship precision H^Γ⁻¹: build the animal-only metafounder relationship A^Γ from supplied Gamma, then apply the ordinary single-step update
over genotyped_rows. group_of is aligned to pedigree.ids; Gamma is supplied and not estimated. At Gamma = 0 this reduces to single_step_inverse with the classical pedigree relationship.
For the raw-array overload, group_of and genotyped_rows use the original ids row order. The wrapper maps them to the normalized pedigree order; G stays in the same order as the supplied genotyped_rows vector.
Experimental and dense/validation-scale. The blend_weight/tau/omega/ridge knobs inherit the ordinary single-step caveat: defaults are not comparator-validated.
HSquared.fit_metafounder_single_step Function
fit_metafounder_single_step(y, X, Z, pedigree, group_of, Gamma, G, genotyped_rows,
sigma_a2, sigma_e2; ...)Supplied-variance single-step GBLUP under the supplied-Γ metafounder relationship. This builds H^Γ⁻¹ with metafounder_single_step_inverse and delegates to fit_gblup. Dense/validation-scale; variance components and Gamma are inputs, not estimates.
HSquared.fit_metafounder_single_step_reml Function
fit_metafounder_single_step_reml(y, X, Z, pedigree, group_of, Gamma, G, genotyped_rows;
iterations = nothing, ...)REML-estimated variance components for the dense validation-scale H^Γ⁻¹ relationship precision. The metafounder covariance Gamma is supplied, not estimated. External comparator parity and R-facing formula syntax remain separate gates.
iterations (hsquared#212) forwards to fit_gblup_reml. Default nothing is byte-identical to the pre-#212-fix call (this keyword did not previously exist).
HSquared.cytoplasmic_relationship Function
cytoplasmic_relationship(pedigree)
cytoplasmic_relationship(ids, sire, dam)Dense cytoplasmic (maternal-lineage) relationship matrix C: C[i, j] = 1 if i and j share a maternal lineage (see maternal_lineage), else 0 (diagonal 1). This is the relationship structure of strictly maternally inherited factors (mitochondrial DNA, cytoplasm). Returned in pedigree.ids (topologically sorted) order, so it aligns with pedigree_inverse.
Experimental Phase 7 primitive — a relationship-construction helper for a non-standard inheritance system. C is a 0/1 same-lineage indicator (rank = number of maternal lineages), so it is singular whenever a lineage has more than one member: use it as the relationship for an i.i.d. cytoplasmic lineage random effect (a grouping), not as a matrix to invert.
HSquared.clonal_relationship Function
clonal_relationship(pedigree, clone_of)Dense additive relationship matrix for a pedigree containing clonal (asexual) ramets. clone_of is aligned to pedigree.ids: entry i is the id of the genet that individual i is a clonal copy of, or an unknown-parent marker (missing, "", "0", 0, …) if i is not a clone. Clonal ramets carry no new Mendelian variation, so each ramet is genetically identical to its genet and inherits the genet's whole row/column of the numerator relationship: C[i, j] = A[rep(i), rep(j)], where rep maps each individual to its ultimate genet (clone links are followed transitively). Returned in pedigree.ids order, aligned with pedigree_inverse.
Experimental Phase 3 non-standard-inheritance primitive. Because clonemates are identical rows, C is rank-deficient — it is a relationship matrix to use directly, not to invert. Record the genets sexually in pedigree; record ramets with unknown parents and mark them through clone_of. Ramets must be terminal: a row marked in clone_of is rejected if it occurs as a sire or dam in the sexual pedigree. Breeding through ramets requires a separate inheritance recursion.
HSquared.maternal_lineage Function
maternal_lineage(pedigree)
maternal_lineage(ids, sire, dam)Maternal-lineage label of each individual: the id of its earliest known maternal ancestor (the maternal founder reached by following dam links). Two individuals share a maternal lineage iff they carry the same label. Returned in pedigree.ids (topologically sorted) order, aligned with cytoplasmic_relationship and pedigree_inverse.
This is the inheritance pattern of strictly maternally transmitted factors (mitochondrial DNA / cytoplasm). An individual with no recorded dam is its own maternal founder.
sourceHSquared.single_step_inverse Function
single_step_inverse(Ainv, A, G, genotyped_rows; tau = 1.0, omega = 1.0,
blend_weight = 0.0, ridge = 0.0)Single-step genomic relationship inverse H⁻¹ (Aguilar et al. 2010; Christensen & Lund 2009): H⁻¹ = A⁻¹ + scatter(τ·Gʷ⁻¹ − ω·A₂₂⁻¹) over the genotyped rows, where A₂₂ = A[g, g] is the pedigree block among the genotyped animals (NOT (A⁻¹)[g, g]) and Gʷ = (1 − blend_weight)·G + blend_weight·A₂₂ is the optionally blended/ridged genomic block. With default tau = omega = 1 and ridge = 0, G = A₂₂ makes H⁻¹ reduce exactly to A⁻¹.
Ainv is the pedigree inverse A⁻¹, A the dense pedigree relationship matrix, G the genomic relationship in the supplied genotyped_rows sequence (the indices refer to rows of A and Ainv). The result is the relationship PRECISION for the single-step animal model — pass it where fit_gblup takes Ginv, or use fit_single_step.
Inputs must be finite and symmetric; A and Ainv must be positive definite, and genotyped_rows must contain unique valid indices. Require tau > 0, omega ≥ 0, blend_weight ∈ [0, 1], and ridge ≥ 0; the returned precision must also be positive definite.
Experimental, dense/validation-scale. The blend_weight/tau/omega/ridge knobs are not comparator-validated (defaults blend_weight = ridge = 0, tau = omega = 1); a singular raw G requires blending/ridging.
HSquared.genomic_relationship_matrix Function
genomic_relationship_matrix(markers; allele_frequencies = nothing, method = :vanraden1,
weights = nothing, backend = :cpu)VanRaden (2008) genomic relationship matrix G from a biallelic marker genotype matrix markers (rows = individuals, columns = markers; entries are the count of one allele, 0/1/2, or an imputed dosage in [0, 2]).
Allele frequencies are estimated from the columns (p_j = mean(markers[:, j]) / 2) unless supplied via allele_frequencies. Two standard constructions:
method = :vanraden1(default) — centered markers, single overall scale:G = Z Zᵀ / (2 Σ_j p_j(1 − p_j)),Z = markers − 2p.method = :vanraden2— per-marker standardized, equal marker weight:G = Zₛ Zₛᵀ / m,Zₛ[:, j] = (markers[:, j] − 2p_j) / √(2 p_j(1 − p_j)),mmarkers. Requires every marker polymorphic (0 < p_j < 1).
Optional per-marker weights (for weighted GBLUP / GWAS-informed prediction, method = :vanraden1 only) give G = Z·diag(w)·Zᵀ / Σ_j w_j·2 p_j(1 − p_j), which reduces exactly to the unweighted method-1 G when the weights are equal.
backend = :cpu (default) computes on the CPU as described. backend = :cuda (v0.7-G-D, opt-in) routes to the GPU twin gpu_genomic_relationship_matrix (the HSquaredCUDAExt extension) — the SAME estimand (it reuses the validated centered_markers verbatim) accelerated on-device; without a CUDA backend in scope it throws a MethodError asking you to load CUDA. The :cpu result is unchanged by the addition of this keyword.
CPU cross-products are symmetrized at floating-point precision before return, so the constructed relationship is exactly symmetric for downstream validation. This leaves the covariance estimand unchanged.
This is the Phase 2 genomic-relationship construction utility — it builds G only. Its regularized inverse is genomic_relationship_inverse, and the experimental supplied-variance GBLUP / SNP-BLUP fitting that consume G / Ginv are fit_gblup / fit_snp_blup. G is typically rank-deficient when there are fewer markers than individuals and needs regularization before inversion.
HSquared.genomic_relationship_inverse Function
genomic_relationship_inverse(G; ridge = 0.01, backend = :cpu)Regularized inverse of a genomic relationship matrix G (e.g. from genomic_relationship_matrix) — the genomic relationship inverse Ginv consumed by the experimental supplied-variance fit_gblup. A genomic G built from markers is usually rank-deficient (markers < individuals), so a ridge is added to the diagonal before inversion: inv(G + ridge·I).
This is a simple ridge regularization. Blending G with a pedigree A (single-step / H-matrix) exists as an internal validation utility and is not yet part of the public surface.
backend = :cpu (default) is unchanged; backend = :cuda (v0.7-G-D, opt-in) routes to the GPU twin gpu_genomic_relationship_inverse (MethodError without CUDA in scope).
HSquared.apy_genomic_relationship_inverse Function
apy_genomic_relationship_inverse(G, core; ridge = 0.0)APY (Algorithm for Proven & Young, Misztal et al. 2014) sparse-structured inverse of a genomic relationship matrix G — the large-genotyped-population scale path for Ginv (v0.8-S4/V8.5). It inverts only the core × core block and adds a DIAGONAL conditional-variance correction for the non-core animals, instead of the full O(n³) dense inverse: partition animals into core c and non-core n,
Mnn = Diagonal(gᵢᵢ − G[i,c]·Gcc⁻¹·G[c,i])(each non-core animal's conditional "Mendelian" variance);Ginv = [[Gcc⁻¹ + Gcc⁻¹Gcn·Mnn⁻¹·GncGcc⁻¹, −Gcc⁻¹Gcn·Mnn⁻¹]; [−Mnn⁻¹·GncGcc⁻¹, Mnn⁻¹]].
core is the vector of core-animal row indices into G. ridge is added to the core block before its Cholesky (a VanRaden G is rank-deficient). As the core set approaches the effective rank of G, Ginv_APY → inv(G); with core = 1:n (all animals) it equals the full regularized inv(G + ridge·I) EXACTLY. Returns the dense n × n Ginv (the APY STRUCTURE — the compute win is inverting Gcc (ncore³) + the diagonal Mnn, O(ncore³ + nnoncore·ncore²), not O(n³); the result is returned dense at validation scale — a sparse/on-device representation is future work).
EXPERIMENTAL, engine-internal, validation-scale, supplied-G; no core-selection algorithm (the caller supplies core), no R surface. Throws when a conditional variance is non-finite or is not positive relative to its subtraction scale (increase the core set or ridge).
HSquared.gpu_genomic_relationship_matrix Function
gpu_genomic_relationship_matrix(markers; allele_frequencies = nothing,
method = :vanraden1, weights = nothing,
precision = Float64)GPU-accelerated VanRaden genomic relationship matrix G — the device twin of genomic_relationship_matrix. STUB: the method is provided by the HSquaredCUDAExt package extension, which activates only when CUDA is loaded (using CUDA). Without a CUDA backend in scope, calling this throws a generic MethodError; load CUDA to activate the extension method.
It is a NUMERICAL ACCELERATION, not a new estimand: it reuses the validated CPU centering and validation (centered_markers) verbatim — same allele frequencies, same VanRaden scale k, same input guards — and runs only the dense W·Wᵀ / k GEMM (and the :vanraden2 / weighted variants) on the device, returning a CPU Matrix{Float64} identical (to floating-point tolerance) to the CPU result. The accepted keywords match genomic_relationship_matrix exactly, plus precision.
precision = Float64 (default) runs the GEMM in Float64 and matches the CPU result to BLAS round-off. precision = Float32 (v0.7-G-B) is an OPT-IN mixed-precision path: the centered marker matrix is downcast to Float32 for the O(n²·m) GEMM (where H100 Float32/TF32 throughput dominates) and the result is upcast to Float64 before the scaling — so the return is ALWAYS Matrix{Float64} (the drop-in contract holds), only the GEMM accumulation precision changes. Float32 is a speed-vs-accuracy trade: the Float32 G differs from the Float64 G by the GEMM round-off (absolute error ~Float32-eps scale, characterized on-device), NOT more accurate and NOT the default. The CPU↔GPU agreement test and the Float32 accuracy + speedup benchmark are opt-in cluster scripts (sim/drac/g1_gpu_genomic.jl, sim/drac/g_b_float32.jl), run on a CUDA device — never in CI. No performance or agreement claim holds until a committed run lands its artifact (Wave F execution model, doc 17).
HSquared.gpu_genomic_relationship_inverse Function
gpu_genomic_relationship_inverse(G; ridge = 0.01)GPU-accelerated regularized inverse of a genomic relationship matrix G — the device twin of genomic_relationship_inverse. STUB: the method is provided by the HSquaredCUDAExt package extension (loads on using CUDA); without a CUDA backend in scope this throws a MethodError.
Mirrors the CPU contract exactly: forms regularized = G + ridge·I, runs the dense Cholesky factorization + inverse on the device (CUSOLVER), and returns a CPU Matrix{Float64} equal (to tolerance) to inv(Symmetric(G + ridge·I)). A non positive-definite regularized G throws the same ArgumentError as the CPU twin ("increase ridge"). Float64 throughout; agreement/benchmark are opt-in cluster scripts (see gpu_genomic_relationship_matrix), never CI.
HSquared.centered_markers Function
centered_markers(markers; allele_frequencies = nothing)Center a biallelic marker matrix for SNP-BLUP / RR-BLUP, using the same allele-frequency centering and VanRaden scaling as genomic_relationship_matrix.
Returns a NamedTuple (W, p, k): W = markers − 2p is the centered marker matrix (each column sums to zero when frequencies are estimated from these markers), p are the per-marker allele frequencies (estimated from the columns unless supplied), and k = 2 Σ_j p_j(1 − p_j) is the VanRaden scale, so that genomic_relationship_matrix(markers) == W * Wᵀ / k. Computing W, p, and k together guarantees they share one p, which the GBLUP↔SNP-BLUP equivalence requires.
HSquared.single_marker_scan Function
single_marker_scan(y, X, markers; allele_frequencies = nothing,
sigma_e2 = nothing, marker_ids = nothing)Fixed-effect single-marker scan for biallelic marker dosages.
Each marker column is centered with centered_markers, residualized against the fixed-effect design X, and tested one marker at a time in the Gaussian linear model y = Xβ + marker * α + e. If sigma_e2 is omitted, the residual variance is the OLS residual mean square of y on X. A supplied value must be finite and positive. The returned NamedTuple contains marker_ids, effects, standard_errors, z_scores, chisq, p_values, bonferroni_p_values, bh_q_values, lod_scores, denominators, p, k, sigma_e2_used, and sigma_e2_source = :supplied | :estimated. p_values are approximate two-sided Gaussian/Wald p-values implied by sigma_e2_used. bonferroni_p_values and bh_q_values are deterministic Bonferroni and Benjamini-Hochberg adjustments over the returned marker set. lod_scores are known-variance fixed-effect LOD-equivalent scores, computed as chisq / (2log(10)).
This is a deterministic Phase 5 validation-scale utility. It is not a mixed model GWAS/QTL scan, does not account for relatedness or population structure, does not compute interval-mapping or mixed-model LOD scores or calibrated / correlated-marker multiple-testing workflows, and does not activate the R-facing marker_scan() formula term.
single_marker_scan(fit::AnimalModelFit, markers; allele_frequencies = nothing,
marker_ids = nothing)Run the fixed-effect single-marker scan of single_marker_scan on a fitted animal model, using the fit's y/X and its fitted residual variance σ²e. This is the relatedness-UNcorrected screen (no Z/Ainv); use mixed_model_marker_scan(::AnimalModelFit, ::AbstractMatrix) for the relatedness-corrected scan.
EXPERIMENTAL, dense/validation-scale; p-values NOT genome-wide calibrated (#48).
sourceHSquared.mixed_model_marker_scan Function
mixed_model_marker_scan(y, X, Z, Ainv, markers, sigma_a2, sigma_e2;
allele_frequencies = nothing, marker_ids = nothing)Supplied-variance mixed-model single-marker scan for biallelic marker dosages.
The helper forms the dense validation-scale marginal covariance V = sigma_a2 * Z * A * Z' + sigma_e2 * I, where A = inv(Ainv), then tests each centered marker as a fixed effect by generalized least squares conditional on the fixed-effect design X. The returned fields mirror single_marker_scan: marker effects, standard errors, Wald z-scores, chi-square statistics, approximate two-sided Gaussian/Wald p-values, Bonferroni-adjusted p-values, Benjamini-Hochberg q-values, LOD-equivalent scores, GLS denominators, marker IDs, allele frequencies, and the VanRaden scale, plus supplied variance components and target = :mixed_model_marker_scan.
This is an engine-internal, dense, supplied-variance Phase 5 utility. It is not variance-component estimation, LOCO, interval mapping, calibrated genome-wide testing, a plotting backend, a bridge payload, or activation of the R-facing marker_scan() formula term.
mixed_model_marker_scan(fit::AnimalModelFit, markers; allele_frequencies = nothing,
marker_ids = nothing)Run the dense supplied-variance mixed-model (relatedness-corrected GLS) marker scan of mixed_model_marker_scan on a fitted animal model, using the fit's y/X/Z, its relationship precision Ainv, and its fitted variance components (σ²a, σ²e). Equivalent to calling the explicit-argument method with those values pulled off fit.spec and fit.variance_components. The fit must have converged. Marker rows must be in the same observation order as fit.spec.y; this method has no observation-ID alignment information.
EXPERIMENTAL, dense/validation-scale. The Wald p-values are NOT genome-wide calibrated (see #48); this does not activate the R marker_scan() formula path or change any bridge payload.
HSquared.loco_relationship_precisions Function
loco_relationship_precisions(markers, marker_groups; allele_frequencies = nothing,
ridge = 0.01)Construct dense leave-one-group-out genomic relationship precisions from marker dosages.
marker_groups must provide one group label per marker column, for example a chromosome label. Labels must be present and non-empty before string conversion; the literal string "missing" is an ordinary supplied label. For each group, the helper drops that group's markers, builds a VanRaden genomic relationship matrix from the remaining markers using genomic_relationship_matrix, and returns its regularized dense inverse from genomic_relationship_inverse. The return value is a Dict{String, Matrix{Float64}} keyed by group label, ready to pass as the relationship_precisions argument to loco_mixed_model_marker_scan.
This is a validation-scale construction helper. It does not choose public LOCO defaults, estimate variance components, calibrate marker-scan p-values, parse marker files, or activate R-facing marker_scan() syntax.
HSquared.loco_mixed_model_marker_scan Function
loco_mixed_model_marker_scan(y, X, Z, relationship_precisions, marker_groups,
markers, sigma_a2, sigma_e2; ...)Leave-one-group-out supplied-variance mixed-model marker scan.
relationship_precisions is a dictionary or named tuple mapping each marker group label (for example a chromosome) to the relationship precision that should be used when testing markers in that group. The helper selects the matching precision for each marker, forms the dense validation-scale GLS covariance, and runs the same marker-by-marker Wald scan as mixed_model_marker_scan. Callers can build the dictionary with loco_relationship_precisions, or provide their own externally constructed matrices. Group labels and relationship keys must be present and non-empty before string conversion.
This helper only selects among supplied matrices; it does not choose public LOCO defaults, estimate variance components, calibrate p-values, run sparse production scans, or activate the R-facing marker_scan() formula term.
HSquared.marker_scan_result_payload Function
marker_scan_result_payload(scan)Return the bridge-facing marker-scan result payload.
This is the stable, boring NamedTuple shape for post-fit GWAS/QTL/eQTL bridge work. It preserves the row-aligned scan fields returned by single_marker_scan, mixed_model_marker_scan, or loco_mixed_model_marker_scan: marker IDs, effect estimates, standard errors, Wald statistics, chi-square values, nominal p-values, Bonferroni/BH-adjusted p-values, LOD-equivalent scores, denominators, allele frequencies, VanRaden scale, optional supplied variance components, and optional LOCO group metadata.
The payload does not run a scan, calibrate genome-wide thresholds, estimate marker-scan variance components, join marker maps, draw figures, or activate the R-facing marker_scan() formula term.
HSquared.genome_wide_threshold_from_null Function
genome_wide_threshold_from_null(null_max_statistics; alpha = 0.05, statistic = :chisq)Turn a SUPPLIED empirical null distribution of per-scan maximum statistics into a genome-wide significance threshold at level alpha: the (1 - alpha) empirical quantile of null_max_statistics. Returns (threshold, alpha, statistic, n_null). alpha must remain in (0, 1) after conversion to Float64; null statistics are converted to finite Float64 values.
The null distribution is the set of per-permutation (or per-null-simulation) maxima (the internal _scan_max_statistic extracts the max chi-square / max -log10 p of a scan); generate it with the opt-in sim/phase5_threshold_calibration.jl harness (phenotype permutation re-runs the marker scan, so the maximum is correlation/LD-aware). Because it uses the distribution of the MAXIMUM over the jointly-scanned markers, the threshold is correlation-aware; it may be less conservative than Bonferroni when markers are strongly redundant, but a finite empirical threshold is not guaranteed to sit below Bonferroni.
The (1 - alpha) type-7 quantile threshold and the add-one genome_wide_pvalue are two DIFFERENT estimators of the same level: they agree asymptotically, but at small n_null the quantile threshold is mildly anti-conservative (e.g. for n_null = 100, alpha = 0.05, an observed value exactly at the threshold has an add-one p of 6/101 ≈ 0.059, not 0.05). genome_wide_pvalue is the exact / conservative decision rule; use the quantile threshold for display and the add-one p-value for the formal accept/reject.
EXPERIMENTAL, validation-scale: this calibrates a threshold FROM a supplied null; it is not itself a production genome-wide-significance claim (that needs a realistic-design calibration — the #48 gate) and has no external comparator yet.
sourceHSquared.genome_wide_pvalue Function
genome_wide_pvalue(observed, null_max_statistics)Empirical genome-wide p-value for an observed per-scan maximum statistic against a SUPPLIED null distribution of maxima: (1 + #{null ≥ observed}) / (n_null + 1) (the standard add-one permutation p-value, so it is never 0). Null values are converted to finite Float64 values; comparisons retain the original Real observed, including values outside the Float64 range. The same null distribution feeds genome_wide_threshold_from_null. EXPERIMENTAL, validation-scale; not a production genome-wide-significance claim.
HSquared.genome_wide_marker_scan Function
genome_wide_marker_scan(y, X, markers; n_permutations = 1000, alpha = 0.05,
sigma_e2 = nothing, allele_frequencies = nothing,
marker_ids = nothing, rng = Random.default_rng())Fixed-effect single-marker scan with a per-dataset residual-permutation null. Runs single_marker_scan, then builds a PER-DATASET residual-permutation null of the per-scan MAXIMUM chi-square (permute the residuals of y on X, re-scan n_permutations times) and returns the add-one empirical genome-wide p-value (genome_wide_pvalue) for the observed maximum AND for each marker, plus the (1 - alpha) permutation threshold. sigma_e2 follows single_marker_scan: omitted means the OLS residual mean square of y on X. The returned fields include sigma_e2_used and sigma_e2_source. Permutation scans reuse that same residual variance.
The null is rebuilt for each dataset. Exact/conservative permutation interpretation requires exchangeability under the permitted permutations. The evidenced calibration scope is intercept-only X with exchangeable Gaussian residuals: validation-scale sim/phase5_qtl_addone_gate.jl and …_design_sweep.jl, plus the larger intercept-only sim/phase5_qtl_rebuild_production_gate.jl. The fixed-null reuse shortcut has a documented anti-conservative calibration result.
Supplied-covariate X is accepted as an experimental residual-permutation utility. General regression residuals need not be exchangeable; add-one counting alone does not establish exactness or family-wise type-I control for that design. Relatedness- corrected mixed-model genome-wide calibration remains separately unvalidated.
Returns the single_marker_scan fields plus genome_wide_p_values (per marker), genome_wide_threshold, genome_wide_p_min (the top marker's genome-wide p), n_permutations, alpha, null_max, a calibration NamedTuple, and target = :genome_wide_marker_scan.
HSquared.gwas_table Function
gwas_table(scan; trait = nothing, total_variance = nothing)
gwas_table(scan, marker_spec::HSMarkerMapSpec; ...)
gwas_table(scan, data::HSData; ...)Prepare a GWAS-labelled table from an already-computed direct marker scan.
This is a semantic wrapper around marker_scan_table: it preserves the row-aligned marker scan fields and adds analysis = :gwas. If trait is supplied, it is recorded as non-empty scalar table metadata. Existing scan statistics, optional marker-map metadata, and optional variance proportions are not recomputed.
This helper does not run a marker scan, estimate variance components, calibrate p-values, choose genome-wide thresholds, activate R-facing marker_scan() syntax, draw plots, or change the bridge payload.
HSquared.qtl_table Function
qtl_table(scan; trait = nothing, total_variance = nothing)
qtl_table(scan, marker_spec::HSMarkerMapSpec; ...)
qtl_table(scan, data::HSData; ...)Prepare a QTL-labelled table from an already-computed direct marker scan.
This is a semantic wrapper around marker_scan_table: it preserves the row-aligned marker scan fields and adds analysis = :qtl. If trait is supplied, it is recorded as non-empty scalar table metadata. The LOD-equivalent scores are the scan's existing fixed-effect known-variance scores; no interval mapping or mixed-model LOD workflow is performed here.
This helper does not run a marker scan, perform interval mapping, estimate variance components, calibrate p-values, choose QTL thresholds, activate R-facing marker_scan() / qtl_scan() syntax, draw plots, or change the bridge payload.
HSquared.eqtl_table Function
eqtl_table(scan; feature = nothing, total_variance = nothing)
eqtl_table(scan, marker_spec::HSMarkerMapSpec; ...)
eqtl_table(scan, data::HSData; ...)Prepare an eQTL-labelled table from an already-computed direct marker scan.
This is a semantic wrapper around marker_scan_table: it preserves the row-aligned marker scan fields and adds analysis = :eqtl. If feature is supplied, it is recorded as non-empty scalar table metadata for the expression feature, gene, or transcript represented by the scan.
This helper does not run expression-wide scans, classify cis/trans windows, join expression or annotation tables, estimate variance components, calibrate p-values, activate R-facing marker_scan() / eQTL syntax, draw plots, or change the bridge payload.
HSquared.marker_effects Function
marker_effects(scan; sort_by = :p_value, top_n = nothing,
decreasing = nothing)
marker_effects(scan, marker_spec::HSMarkerMapSpec; ...)
marker_effects(scan, data::HSData; ...)Prepare a deterministic marker-effect summary from a direct marker-scan result.
The helper expects the scan fields returned by single_marker_scan, mixed_model_marker_scan, or loco_mixed_model_marker_scan: marker effects, standard errors, Wald statistics, chi-square values, p-values, Bonferroni p-values, Benjamini-Hochberg q-values, LOD-equivalent scores, and denominators. It returns those fields sorted into a compact NamedTuple with scan_indices, sort_by, decreasing, top_n, and the scan target when available. Supported sort_by values are :p_value, :bonferroni_p_value, :bh_q_value, :chisq, :lod_score, :effect, and :abs_effect; p-value metrics sort ascending by default, while effect-size/statistic metrics sort descending by default.
When an already-validated HSMarkerMapSpec or HSData with marker metadata is supplied, marker IDs must match exactly and chromosome / position vectors are aligned to the returned summary order.
This is a direct Julia summary helper only. It does not calibrate p-values, choose thresholds, draw plots, activate R-facing marker_scan() syntax, or change the bridge payload.
HSquared.marker_variance_explained Function
marker_variance_explained(scan; total_variance = nothing,
sort_by = :marker_variance, top_n = nothing,
decreasing = nothing)
marker_variance_explained(scan, marker_spec::HSMarkerMapSpec; ...)
marker_variance_explained(scan, data::HSData; ...)Prepare deterministic marker-level variance-contribution summaries from a direct marker-scan result.
The helper expects scan marker_ids, marker effects, and allele frequencies p such as those returned by single_marker_scan, mixed_model_marker_scan, or loco_mixed_model_marker_scan. For each biallelic marker it computes 2p(1-p) * effect^2. If total_variance is supplied, it also returns proportion_variance_explained = marker_variance / total_variance. Supported sort_by values are :marker_variance, :proportion_variance_explained, :allele_variance, :effect, :abs_effect, and :p_value; p-value sorting requires a scan p_values field. Variance and effect-size metrics sort descending by default, while p-values sort ascending.
When an already-validated HSMarkerMapSpec or HSData with marker metadata is supplied, marker IDs must match exactly and chromosome / position vectors are aligned to the returned summary order.
This is a direct Julia summary helper only. It does not estimate marker-scan variance components, calibrate p-values, claim calibrated PVE/model R², choose thresholds, draw plots, activate R-facing marker_scan() syntax, or change the bridge payload.
HSquared.marker_scan_table Function
marker_scan_table(scan; total_variance = nothing)
marker_scan_table(scan, marker_spec::HSMarkerMapSpec; ...)
marker_scan_table(scan, data::HSData; ...)Prepare a deterministic row-aligned marker-scan table from a direct scan result.
The helper expects the scan fields returned by single_marker_scan, mixed_model_marker_scan, or loco_mixed_model_marker_scan: marker effects, standard errors, Wald statistics, chi-square values, p-values, Bonferroni p-values, Benjamini-Hochberg q-values, LOD-equivalent scores, denominators, and allele frequencies. It preserves the original scan order and returns those vectors with scan_indices, the scan target, allele variances, and marker-level variance contributions 2p(1-p) * effect^2. If total_variance is supplied, it also returns proportion_variance_explained = marker_variance / total_variance.
When an already-validated HSMarkerMapSpec or HSData with marker metadata is supplied, marker IDs must match exactly and chromosome / position vectors are aligned to the scan order.
This is a direct Julia table-preparation helper only. It does not sort markers, draw plots, calibrate p-values, choose thresholds, estimate marker-scan variance components, activate R-facing marker_scan() syntax, or change the bridge payload.
HSquared.marker_genomic_inflation Function
marker_genomic_inflation(scan; expected_median = 0.454936423119572)Compute a genomic-control-style inflation diagnostic from a direct marker-scan result.
The helper expects a chisq field such as the one returned by single_marker_scan, mixed_model_marker_scan, or loco_mixed_model_marker_scan. It returns a compact NamedTuple with lambda_gc = median(chisq) / expected_median, the observed median chi-square, the expected median, the number of markers, and the scan target when available. The default expected median is the 0.5 quantile of a one-degree-of-freedom chi-square distribution.
This is a diagnostic summary only. It does not calibrate p-values, correct the scan statistics, choose genome-wide thresholds, or activate R-facing marker_scan() syntax.
HSquared.marker_manhattan_data Function
marker_manhattan_data(scan; chromosomes = nothing, positions = nothing,
p_floor = floatmin(Float64), chromosome_gap = 1.0)
marker_manhattan_data(scan, marker_spec::HSMarkerMapSpec; ...)
marker_manhattan_data(scan, data::HSData; ...)Prepare plot-ready Manhattan data from a direct single_marker_scan result.
The helper returns marker IDs, raw p-values, -log10(p) values, chromosome labels, marker positions, cumulative plotting positions, and a deterministic plot order. If chromosome or position metadata is omitted, all markers are placed on chromosome "1" with sequential positions. Zero p-values are floored only for -log10 display values, and the floor is returned in the output. When a validated HSMarkerMapSpec or HSData with marker metadata is supplied, marker IDs must match the scan exactly; chromosome and position metadata are then aligned to the scan marker order, while chromosome display order follows the marker-map order.
This is a data-preparation helper only. It does not draw a plot, does not parse marker maps, does not run a mixed-model marker scan, and does not activate the R-facing marker_scan() formula term.
HSquared.marker_region_data Function
marker_region_data(scan; chromosomes, positions, chromosome, start = nothing,
stop = nothing, flank = 0, total_variance = nothing,
p_floor = floatmin(Float64))
marker_region_data(scan, marker_spec::HSMarkerMapSpec; ...)
marker_region_data(scan, data::HSData; ...)Prepare scan data for one chromosome or chromosome-window region.
The helper reuses marker_scan_table validation, then subsets the row-aligned marker-scan table by chromosome and optional coordinate bounds. When start and/or stop are supplied, a non-negative flank expands that window before filtering. Returned rows are ordered by position within the region, with original scan indices preserved.
This is a direct Julia data-preparation helper only. It does not draw a plot, choose thresholds, calibrate p-values, run interval mapping, activate gwas_table() / qtl_table() / eqtl_table(), activate R-facing marker_scan() syntax, or change the bridge payload.
HSquared.marker_significance_summary Function
marker_significance_summary(scan; alpha = 0.05)Summarize nominal marker hits from a direct marker-scan result.
The helper expects the scan fields returned by single_marker_scan, mixed_model_marker_scan, or loco_mixed_model_marker_scan: marker IDs, raw p-values, Bonferroni-adjusted p-values, Benjamini-Hochberg q-values, chi-square values, and LOD-equivalent scores. It returns per-marker significance flags, counts, marker IDs and scan indices for raw, Bonferroni, and BH summaries, plus the top marker by raw p-value.
The thresholds are nominal summaries over the markers already present in the scan: raw p <= alpha, Bonferroni raw-p threshold alpha / m, adjusted Bonferroni p_adj <= alpha, and BH q <= alpha. This helper does not calibrate correlated-marker genome-wide thresholds, estimate effective marker counts, correct p-values, choose public GWAS/QTL thresholds, draw plots, or activate R-facing marker_scan() syntax.
HSquared.marker_qq_data Function
marker_qq_data(scan; p_floor = floatmin(Float64))Prepare plot-ready QQ data from a direct single_marker_scan result.
The helper sorts observed p-values from smallest to largest and returns the paired expected uniform order-statistic p-values, raw and sorted p-values, marker IDs, sorted marker IDs, -log10 observed and expected values, the sort order, and the display p-value floor. Zero p-values are floored only for display values; raw p-values are preserved.
This is a data-preparation helper only. It does not draw a plot, estimate genomic inflation, calibrate p-values, run a mixed-model marker scan, or activate the R-facing marker_scan() formula term.
HSquared.sparse_csc_matrix Function
sparse_csc_matrix(nrow, ncol, colptr, rowval, nzval; index_base = :zero)Construct a Julia SparseMatrixCSC{Float64,Int} from compressed sparse column slots.
This is a bridge utility for R Matrix::dgCMatrix payloads. R stores column pointers and row indices with zero-based indexing, so index_base = :zero is the default. Use index_base = :one for Julia-style CSC slots.
HSquared.AnimalModelSpec Type
AnimalModelSpecBridge-ready low-level specification for the first Gaussian animal-model engine target.
The spec stores the response y, fixed-effect design X, animal-effect design Z, sparse relationship precision Ainv, encoded animal ids, family marker, fitting method, and optionally relationship_diag, a precomputed diag(inv(Ainv)) (1 .+ F for a pedigree Ainv; nothing when absent). It validates dimensions only; it does not fit a model.
HSquared.animal_model_spec Function
animal_model_spec(y, X, Z, Ainv; ids = nothing, family = GaussianFamily(),
method = :REML, relationship_diag = nothing)Validate and store the low-level inputs for the first Gaussian animal-model engine target.
This function checks dimensions, ID length, family, and method. It is the Julia-side mirror of the R parser contract, not a fitting routine.
relationship_diag, when supplied, must be diag(inv(Ainv)) in the row order of Ainv: for a pedigree precision that is 1 .+ inbreeding_coefficients(pedigree) of the same normalized pedigree that built Ainv. reliability then reads its denominator from it (under the default method = :auto) instead of a selected inverse of Ainv. It is checked for length, finiteness, and positivity only — not for consistency with Ainv — so pass it only when it comes from the same source. The R bridge attaches it automatically when it builds Ainv from pedigree rows.
HSquared.GaussianLikelihoodResult Type
GaussianLikelihoodResultResult from evaluating the Gaussian animal-model log-likelihood at supplied variance components.
sourceHSquared.gaussian_loglik Function
gaussian_loglik(spec, sigma_a2, sigma_e2; method = spec.method,
max_dense_cells = 1_000_000)Evaluate the Gaussian ML or REML log-likelihood at supplied variance components.
This Phase 1 evaluator is deliberately conservative: it forms dense matrices from the validated AnimalModelSpec so the likelihood can be tested before the production sparse solver lands. It does not optimize variance components and does not return a fitted model. max_dense_cells is a safety guard for this temporary dense path: if nobs^2 + nanimals^2 exceeds it, this throws an ArgumentError naming the observed cell count and the effective cap (hsquared#214) before any dense object is formed.
HSquared.sparse_reml_loglik Function
sparse_reml_loglik(spec, sigma_a2, sigma_e2)Evaluate the Gaussian REML log-likelihood at supplied positive variance components using the sparse Henderson mixed-model-equation identity.
This is a Phase 1 validation bridge toward the production sparse optimizer. It does not estimate variance components and it only evaluates REML.
sourceHSquared.sparse_multi_reml_loglik Function
sparse_multi_reml_loglik(y, X, effects, sigmas, sigma_e2) -> (loglik, beta, us)Evaluate the Gaussian REML log-likelihood, fixed effects, and per-block BLUPs of the general K-independent-random-effect model y = X·β + Σᵢ Zᵢ·uᵢ + e, uᵢ ~ N(0, σᵢ²·Aᵢ), e ~ N(0, σ_e²·I), at supplied positive variance components, using the sparse Henderson MME determinant identity — the K-block generalization of sparse_reml_loglik. effects is a vector of (Zᵢ, Ainvᵢ) pairs (same contract as multi_effect_mme).
The log-likelihood uses the package-wide LOGLIK_CONVENTION_FULL −0.5·[(n−p)·log(2π) + log|R| + log|G| + log|C| + y'Py] (identical to sparse_reml_loglik / fit_ai_reml). It therefore equals the dense fit_multi_effect_reml REML objective (_multi_effect_dense, LOGLIK_CONVENTION_OMIT_2PI) PLUS −0.5·(n−p)·log(2π); this offset is the only difference and is exact (#365). Fit results that call this path expose loglik_convention = :reml_full_constant and loglik_full_constant_offset = 0.0. Engine-internal, supplied-variance; it does not estimate.
HSquared.HendersonMMEResult Type
HendersonMMEResultResult from solving Henderson's mixed-model equations at supplied variance components.
This is a Phase 1 engine utility. It uses sparse design and relationship precision matrices, but it does not estimate variance components and is not a production sparse fitting claim by itself.
sourceHSquared.henderson_mme Function
henderson_mme(spec, sigma_a2, sigma_e2)Solve Henderson's mixed-model equations for fixed effects and animal-effect BLUPs/EBVs at supplied positive variance components.
This forms the sparse equation system [X'R^-1X X'R^-1Z; Z'R^-1X Z'R^-1Z + Ainv / sigma_a2] with R = sigma_e2 I. The coefficient matrix is symmetric positive definite by construction, so the solve uses cholesky(Symmetric(lhs)), the same factor sparse_reml_loglik and fit_ai_reml already use. A non-positive-definite coefficient matrix throws PosDefException. It is a supplied-variance solver and does not optimize variance components.
HSquared.solve_animal_model_pcg Function
solve_animal_model_pcg(spec, sigma_a2, sigma_e2; tol = 1e-10, maxiter = 1000,
preconditioner = :jacobi, matrix_free = false)Solve the supplied-variance Gaussian animal-model mixed-model equations by preconditioned conjugate gradient — the ITERATIVE companion of the direct henderson_mme factorization. It solves the SAME sparse symmetric positive-definite system C·[β; u] = rhs as henderson_mme iteratively, never forming a Cholesky factor. preconditioner = :jacobi (default) uses the diagonal preconditioner M⁻¹ = 1/diag(C); :ichol uses an incomplete-Cholesky IC(0) factor of the assembled C (a stronger, still-sparse preconditioner, with a diagonal-shift fallback on breakdown — requires matrix_free = false); :none is plain CG.
matrix_free = false (default) assembles C once (_sparse_mme_system) and applies it. matrix_free = true applies C·v directly from the sparse X, Z, Ainv matvecs (common = X·v_β + Z·v_u; top = X'·common/σe²; bottom = Z'·common/σe² + Ainv·v_u/σa²) and uses a matrix-free Jacobi diagonal — C is NEVER assembled. Both paths return the SAME solution (validated bit-for-bit close); the matrix-free path removes the C-assembly memory, the foundation for a future large-pedigree solver. Still no performance claim (no benchmark recorded).
Returns a NamedTuple:
beta— fixed effects (1:nfixedof the solution);breeding_values = (ids, values)— EBVs (pedigree.idsorder);iterations— CG iterations taken;relative_residual—‖rhs − C·x‖ / ‖rhs‖at the returned solution;converged— whetherrelative_residual ≤ tol;preconditioner.
EXPERIMENTAL. This is a CORRECTNESS primitive: it is validated to recover the direct henderson_mme solution (β and EBVs) to a tight tolerance on the tiny, Mrode9-shaped, and larger validation fixtures, and the Jacobi preconditioner is validated to reach the same solution (in no more iterations than plain CG). It makes NO performance / large-pedigree scaling claim — _sparse_mme_system still assembles C explicitly — and is not the default fit path. It is the iterative-solver foundation the production sparse path will build on. sigma_a2/sigma_e2 are SUPPLIED, not estimated.
HSquared.AnimalModelFit Type
AnimalModelFitExperimental low-level Gaussian animal-model fit object.
This is returned only for validated AnimalModelSpec inputs. It uses the current dense likelihood evaluator or sparse REML validation objective and a conservative optimizer path. The caller's original specification is retained when its method and converted relationship precision are unchanged. A method override or precision canonicalization is recorded in a validated specification.
HSquared.BreedingValues Type
BreedingValuesExperimental low-level container for animal-effect BLUPs/EBVs.
sourceHSquared.fit_sparse_reml Function
fit_sparse_reml(spec; initial = (sigma_a2 = 1.0, sigma_e2 = 1.0),
iterations = 1_000)Optimize the sparse Gaussian REML validation objective over positive variance components.
The optimizer works on log-variance parameters and uses sparse_reml_loglik as the objective. This is a Phase 1 validation path toward sparse fitting. It is REML-only, not AI-REML, not the default fitting path, and not a production sparse solver.
converged is the optimizer stopping flag for a finite returned result. It does not certify score stationarity or a closed-boundary optimum. Exact zero components are unavailable in this positive log-variance parameterization; nonfinite proposals are rejected and nonfinite final results raise an error.
HSquared.fit_ai_reml Function
fit_ai_reml(spec; initial = (sigma_a2 = 1.0, sigma_e2 = 1.0),
iterations = 100, tol = 1e-8)Estimate the Phase 1 Gaussian animal-model variance components by average-information (AI) REML.
Each iteration solves the sparse Henderson mixed-model equations, reads the variance-component score from the BLUP solution and the Takahashi selected inverse (the tr(Ainv * C^uu) term), forms the average-information matrix from two working-variate re-solves that reuse the same Cholesky factor, and takes an AI/Newton step with step-halving to keep the variance components positive.
REML-only and experimental: it is validated to recover the same optimum as the dense and sparse NelderMead optimizers, but is not yet checked against external comparators or hardened for boundary/large-pedigree cases. For a Gaussian linear mixed model whose covariance is linear in the variance components, average information is the arithmetic average of observed and expected REML information. Its working-variate quadratic need not equal the observed information. This Gaussian identity does not directly extend to Laplace-approximated non-Gaussian models.
An optional EM-REML warm-start (em_warmup, default 0) runs that many EM-REML iterations before the AI/Newton loop. The EM update is the closed form that zeroes the REML score, so it is monotone and stays strictly inside the parameter space — a few iterations hand the AI step a good, in-bounds start. On an identified problem this is optimum-invariant (the converged estimates are unchanged) and improves robustness to poor starting values; em_warmup = 0 is byte-identical to the pre-warm-start path. It does NOT change the σ²→0 / non-identified boundary behaviour (there the fit still returns converged = false), and on a non-identified surface the warm-start can shift the non-converged estimate. If a warm-start update makes the next mixed-model system non-positive-definite, the fitter returns the last factorable variances with converged = false and optimizer_status = "non_positive_definite" — read the flag, not the value.
When the selected-inverse score is cancellation-dominated and the pedigree has more than 512 animals (_AI_REML_BOUNDARY_TRACE_MAX_COLUMNS), AI-REML stops with optimizer_status = "boundary_score_unresolved" and returns the iterate. That iterate is not a certified REML optimum. Use fit_sparse_reml or raise em_warmup on the same spec when the optimum is at or near σ²a = 0 (HSquared.jl #407).
When the REML optimum is on the residual boundary (σ²e → 0), AI-REML often does not converge. The iterate can sit at a start-dependent point and the REML log-likelihood can get worse with more iterations. If the returned residual share is at or near zero, the fitter warns and still returns the iterate with converged = false. Use fit_sparse_reml on the same spec (HSquared.jl #412). This is a different failure from the q > 512 boundary_score_unresolved abort (#407).
A founders-only pedigree (A = I) with one record per animal does not identify σ²a separately from σ²e. The REML surface is then a ridge in σ²a + σ²e, and the reported h² follows the start ratio. The fitter warns and still returns the iterate; converged = true is not an identified heritability (HSquared.jl #422).
HSquared.fit_gblup Function
fit_gblup(y, X, Z, Ginv, sigma_a2, sigma_e2; ids = nothing, method = :REML)Genomic BLUP (GBLUP) at supplied variance components: solve the Gaussian animal model with a genomic relationship inverse Ginv in place of the pedigree Ainv, reusing the existing Henderson mixed-model-equation solve.
Ginv is the (regularized) inverse of a genomic relationship matrix, e.g. from genomic_relationship_inverse. With sample-estimated allele frequencies, a VanRaden G is rank-deficient (column-centering puts the all-ones vector in its null space, so rank(G) ≤ n − 1). Supplied frequencies do not imply that null vector. A singular G must be regularized before inversion; that is genomic_relationship_inverse's job, not this function's. sigma_a2 and sigma_e2 are the supplied genomic and residual variances; GBLUP here does not estimate them.
This is a thin convenience over animal_model_spec + henderson_mme: the genomic precision enters the same Ainv slot the pedigree animal model uses, so it returns the same HendersonMMEResult and works with every existing extractor (fixed_effects, breeding_values, heritability, …). It is experimental and engine-internal; the user-facing R genomic() model-spec mapping is coordinated separately and not part of this function. The dense Ginv path is validation-scale only — it does not gain the sparse selected-inversion advantage.
HSquared.fit_snp_blup Function
fit_snp_blup(y, X, markers, sigma_g2, sigma_e2; allele_frequencies = nothing,
ids = nothing)SNP-BLUP / RR-BLUP at supplied variance components: estimate marker effects with the existing Henderson MME, treating the centered markers as the random-effect design with an identity prior.
The centered markers W (from centered_markers) are the random-effect design Z, the relationship inverse is the identity I_m (markers a priori independent), and the per-marker variance is sigma_g2 / k with k = 2 Σ_j p_j(1 − p_j). Returns a NamedTuple (marker_effects, gebv, beta, k, p) where marker_effects = â, gebv = W·â are the implied genomic breeding values, and beta are the fixed effects.
gebv equals the GBLUP genomic breeding values for the same data and variances, exactly for G itself (the GBLUP↔SNP-BLUP equivalence) — but through fit_gblup's required Ginv = inv(G + ridge·I), the package's only GBLUP route, the two differ by O(ridge): about 1–2 % of sd(gebv) at the package default ridge = 0.01, because the ridge changes the covariance kernel. See issue #333 for the measured per-cell gap; Genomic models explains why the ridge breaks the equivalence but does not carry the measurement. The random block is deliberately labelled marker_effects (not breeding_values/EBV), because on this spec the random effects are marker effects, not animal breeding values. Experimental, supplied-variance only (no variance-component estimation); unweighted VanRaden method-1 / single identity prior only.
HSquared.fit_gblup_reml Function
fit_gblup_reml(y, X, Z, Ginv; initial = (sigma_a2 = 1.0, sigma_e2 = 1.0),
target = :ai_reml, iterations = nothing, ids = nothing)GBLUP with REML-estimated variance components: build the genomic animal-model spec (genomic precision Ginv in the Ainv slot) and estimate (sigma_a2, sigma_e2) by REML, instead of supplying them as in fit_gblup. A one-call convenience over animal_model_spec(...; method = :REML) + fit_animal_model; target = :ai_reml (default) or :sparse_reml. Returns the AnimalModelFit (estimated variance_components, genomic breeding_values, likelihood, converged, …).
iterations (hsquared#212), when supplied, forwards to whichever target optimizer is selected (fit_ai_reml's own default is 100; fit_sparse_reml's and fit_variance_components's is 1000) via fit_animal_model's kwargs... passthrough. Default nothing omits the keyword entirely, so the call is byte-identical to the pre-#212-fix behaviour (each target's own hardcoded default).
Experimental, dense/validation-scale (the dense Ginv path gains no sparse selected-inversion advantage); no external comparator, no R model-spec.
HSquared.fit_snp_blup_reml Function
fit_snp_blup_reml(y, X, markers; initial = (sigma_a2 = 1.0, sigma_e2 = 1.0),
allele_frequencies = nothing, target = :ai_reml, ids = nothing)SNP-BLUP / RR-BLUP with REML-estimated variance components: estimate the per-marker random-effect variance and residual variance by REML on the centered-marker spec (markers as the random-effect design with an identity prior), instead of supplying them as in fit_snp_blup. Returns a NamedTuple (marker_effects, gebv, beta, sigma_g2, sigma_e2, k, p, loglik, converged), where sigma_g2 = σ̂²_marker · k is the total genomic variance (k = 2 Σ_j p_j(1 − p_j)).
At the estimated variances this reproduces fit_snp_blup exactly (the REML wrapper just supplies the estimated variances). Experimental, dense/validation-scale, unweighted VanRaden method-1 / single identity prior; no external comparator.
HSquared.fit_single_step Function
fit_single_step(y, X, Z, Ainv, A, G, genotyped_rows, sigma_a2, sigma_e2;
tau = 1.0, omega = 1.0, blend_weight = 0.0, ridge = 0.0, ids = nothing)Single-step GBLUP at supplied variance components: build the single-step relationship inverse H⁻¹ (single_step_inverse) and solve the Gaussian animal model with H⁻¹ as the relationship precision (via fit_gblup). With default tau = omega = 1 and ridge = 0, G = A₂₂ reproduces the pedigree animal model exactly. Experimental, dense/validation-scale; supplied-variance.
HSquared.fit_single_step_reml Function
fit_single_step_reml(y, X, Z, Ainv, A, G, genotyped_rows;
tau = 1.0, omega = 1.0, blend_weight = 0.0, ridge = 0.0,
initial = (sigma_a2 = 1.0, sigma_e2 = 1.0), target = :ai_reml,
iterations = nothing, ids = nothing)Single-step GBLUP with REML-estimated variance components: build H⁻¹ (single_step_inverse) and estimate (sigma_a2, sigma_e2) by REML on the single-step spec (via fit_gblup_reml). With default tau = omega = 1 and ridge = 0, G = A₂₂ reproduces the pedigree-REML optimum. Experimental, dense/validation-scale.
iterations (hsquared#212) forwards to fit_gblup_reml (and from there to the target optimizer). Default nothing is byte-identical to the pre-#212-fix call (this keyword did not previously exist).
HSquared.repeatability_mme Function
repeatability_mme(y, X, Z, Ainv, sigma_a2, sigma_pe2, sigma_e2; ids = nothing)Supplied-variance Henderson solve of the repeatability / permanent-environment animal model with repeated records:
y = X·β + Z·a + Z·pe + e,
a ~ N(0, sigma_a2·A), pe ~ N(0, sigma_pe2·I), e ~ N(0, sigma_e2·I),where Z is the record→animal incidence (shared by the additive genetic effect a and the permanent-environment effect pe), and Ainv is the relationship inverse. The mixed-model equations carry a block-diagonal relationship precision blockdiag(Ainv/sigma_a2, I/sigma_pe2) for the stacked random effect [a; pe].
This is the first Phase-3 (standard quantitative-genetic) engine slice: a supplied-variance MME solve (it does not estimate the variance components), the analogue of henderson_mme for two random effects. Experimental and engine-internal; the R permanent() / repeatability model-spec mapping and REML estimation of the three variance components are coordinated separately and not part of this function. Returns a NamedTuple (beta, animal_effects, permanent_effects, variance_components). Identifiability of a vs pe requires repeated records (animals with more than one record).
HSquared.fit_repeatability_reml Function
fit_repeatability_reml(y, X, Z, Ainv; initial, iterations = 200, ids = nothing,
max_dense_cells = 1000000)Estimate the three variance components (sigma_a2, sigma_pe2, sigma_e2) of the repeatability / permanent-environment animal model by REML, by maximizing the dense two-random-effect REML log-likelihood over the log-variances (NelderMead).
Returns a NamedTuple with variance_components, the repeatability t = (sigma_a2 + sigma_pe2) / total, the heritability h² = sigma_a2 / total, beta, the a / pe BLUPs at the estimate, loglik, and converged, plus #365 convention fields (loglik_convention = :reml_omit_2pi, loglik_full_constant_offset = −½(n−p)log(2π), loglik_comparable_across_routes = false). Raw loglik is not comparable to the sparse repeatability route (fit_multi_effect / fit_sparse_multi_effect_aireml under scale_method = "auto") without comparable_loglik — AIC / LRT across those routes is silently wrong by exactly that offset.
Experimental and validation-scale: it forms the dense n×n marginal covariance, so it is for small problems, not production. REML-only. Uncertainty intervals for t / h² and the R model-spec mapping are not part of this function. Separating sigma_a2 from sigma_pe2 needs relationship contrast and replication; on small data the optimum can sit on a boundary (one variance → 0).
Unlike fit_variance_components/gaussian_loglik, this function takes raw y/X/Z/Ainv (no AnimalModelSpec) and previously formed the dense n×n marginal via inv(Symmetric(Matrix{Float64}(Ainv))) with no size guard at all (hsquared#217). max_dense_cells (default 1000000) now guards this path the same way the spec-based fitters are guarded, checked BEFORE the dense inverse is formed.
HSquared.comparable_loglik Function
comparable_loglik(fit) -> Float64Return fit.loglik on the package-wide LOGLIK_CONVENTION_FULL scale, using fit.loglik_full_constant_offset when present. Use this (not raw loglik) for AIC / BIC / LRT that span dense and sparse multi-effect / repeatability routes (#365). Throws if the fit does not carry convention metadata.
HSquared.reml_full_constant_offset Function
reml_full_constant_offset(n, p) -> Float64Additive term that converts an :reml_omit_2pi log-likelihood into the :reml_full_constant scale: −½(n−p)log(2π). Zero when n == p (degenerate REML). Exact identity behind HSquared.jl #365.
HSquared.LOGLIK_CONVENTION_FULL Constant
Package-wide full-constant REML loglik (−½[(n−p)log(2π) + …]). See #365.
HSquared.LOGLIK_CONVENTION_OMIT_2PI Constant
Dense multi-effect / two-effect / repeatability REML loglik (omits (n−p)log(2π)). See #365.
HSquared.two_effect_mme Function
two_effect_mme(y, X, Z1, Ainv1, Z2, Ainv2, sigma1, sigma2, sigma_e2;
ids1 = nothing, ids2 = nothing)Supplied-variance Henderson solve of a Gaussian model with two independent random effects:
y = X·β + Z1·u1 + Z2·u2 + e,
u1 ~ N(0, sigma1·A1), u2 ~ N(0, sigma2·A2), e ~ N(0, sigma_e2·I),with relationship inverses Ainv1 = A1⁻¹, Ainv2 = A2⁻¹. The mixed-model equations carry the block-diagonal precision blockdiag(Ainv1/sigma1, Ainv2/sigma2) for the stacked random effect [u1; u2].
This is the general engine kernel for the standard two-random-effect quantitative-genetic models: repeatability / permanent environment (Z2 = Z1, A2 = I; see repeatability_mme), common environment (Z2 = group incidence, A2 = I), and maternal-environment (Z2 = dam incidence, A2 = I). Experimental, supplied-variance, engine-internal — it does not estimate variances and does not cover correlated direct–maternal genetic effects (which need a 2×2 genetic covariance). Returns (beta, effect1, effect2, variance_components).
HSquared.fit_two_effect_reml Function
fit_two_effect_reml(y, X, Z1, Ainv1, Z2, Ainv2; initial, iterations = 200,
ids1 = nothing, ids2 = nothing)REML estimation of the variance components (sigma1, sigma2, sigma_e2) of the general two-independent-random-effect model (see two_effect_mme), by maximizing the dense two-effect REML log-likelihood (NelderMead). Covers common-environment (c² = ratio2) and maternal-environment variance estimation.
Returns a NamedTuple with variance_components, ratio1 = sigma1/total, ratio2 = sigma2/total, beta, the two BLUPs, loglik, and converged, plus #365 convention fields (loglik_convention = :reml_omit_2pi, loglik_full_constant_offset, loglik_comparable_across_routes = false). Raw loglik omits −½(n−p)log(2π); use comparable_loglik before AIC / LRT against sparse / gaussian_loglik routes.
Experimental, dense/validation-scale, REML-only; uncertainty intervals and the R model-spec mapping are not part of this function. On small data the optimum can sit on a boundary (a variance → 0).
converged is Optim.converged for NelderMead: the simplex objective spread fell below g_tol. It is not a score, information-rank, or multi-start check. When the two random effects are linearly dependent on the residual (the usual animal + dam case where each dam has one sire, so A = 0.5 I + 0.5 D among phenotyped animals), only two combinations of (Va, Vdam, Ve) are identified. Different starts can then return different h2 values at the same log-likelihood. Read ratio1 as one point on that ridge, not a unique heritability (HSquared.jl #416).
Default starts are (sigma1, sigma2, sigma_e2) = (1, 1, 1). REML is scale-equivariant, but NelderMead from that unit-scale start can stop at a wrong near-zero-h2 point when var(y) is far from 1, still with converged = true (HSquared.jl #413). The fitter warns in that case; pass an initial on the scale of the response (or rescale y) rather than reading the green flag as an identified optimum. Recorded recovery evidence used the unit-scale start, so the default is unchanged.
HSquared.two_effect_ratio_interval Function
two_effect_ratio_interval(y, X, Z1, Ainv1, Z2, Ainv2; level = 0.95,
which = :both, initial = ..., iterations = 200,
ids1 = nothing, ids2 = nothing, fd_step = 1e-4,
boundary_tol = 1e-6)Asymptotic delta-method confidence interval(s) for the variance ratios of the general two-effect REML model (fit_two_effect_reml):
ratio1 = σ1² / (σ1² + σ2² + σe²), ratio2 = σ2² / (σ1² + σ2² + σe²)(e.g. h² and c² / m² for a common-environment / maternal model). Fits by REML, forms the observed information as the central finite-difference Hessian of the two-effect REML log-likelihood (_two_effect_dense) at the optimum, and applies the delta method on the logit scale (so each interval lies in (0, 1)) — the same machinery as repeatability_interval and heritability_interval method = :delta.
which selects which ratio(s) to build an interval for (:both, :ratio1, or :ratio2); the unrequested ratio still reports its point estimate with a NaN interval. Returns a NamedTuple (ratio1, ratio2, level, converged) where each of ratio1/ratio2 is (estimate, lower, upper, se, lower_clamped, upper_clamped, boundary).
Boundary honesty: when a component sits on the variance boundary (σ → 0, i.e. σ_i / total ≤ boundary_tol) the ratio built on it is flagged boundary = true with a NaN interval — never a spuriously tight CI. A ratio built on a non-boundary numerator drops the degenerate component and uses the corresponding information sub-block conditional on that component held at its supplied value. The reported ratio and its gradient retain the full denominator, so e.g. ratio1 remains well-defined when only σ2² → 0; in that reduction it recovers heritability_interval on the underlying animal model. converged is carried from the REML fit.
Experimental, asymptotic, delta-method, REML only; the interval is a large-sample approximation and is NOT coverage-calibrated — on small samples the REML surface is flat and the interval is unreliable (the parametric bootstrap, bootstrap_variance_component_interval, is the only finite-sample-aware path). No calibrated coverage is claimed.
HSquared.multi_effect_mme Function
multi_effect_mme(y, X, effects, sigmas, sigma_e2; ids = nothing)Supplied-variance Henderson solve of the general model with an ARBITRARY number K of independent random effects — the K-block generalization of two_effect_mme:
y = X·β + Σ_{i=1}^{K} Z_i·u_i + e,
u_i ~ N(0, sigma_i·A_i), e ~ N(0, sigma_e2·I),where effects is a vector of (Z_i, Ainv_i) pairs (Z_i the n×q_i record→level incidence, Ainv_i the q_i×q_i relationship precision — pass I for a plain i.i.d. (1|group) effect), sigmas[i] the variance of effect i, and sigma_e2 > 0. The stacked random precision is blockdiag(Ainv_1/σ_1, …, Ainv_K/σ_K) over [u_1; …; u_K] and the random design is hcat(Z_1, …, Z_K), exactly as the two-effect kernel. ids may be nothing or a length-K vector of per-effect id vectors.
Returns (beta, effects = [(ids, values), …], variance_components = (sigmas, sigma_e2)). The K=2 case is byte-identical to two_effect_mme. Experimental, supplied-variance, engine-internal — it does not estimate variances and does not cover correlated effects (a 2×2 direct–maternal G needs kron(inv(G), Ainv), a different structure).
HSquared.fit_multi_effect_reml Function
fit_multi_effect_reml(y, X, effects; initial = nothing, iterations = 200,
ids = nothing, max_dense_cells = 1_000_000)REML estimation of the K+1 variance components (sigmas..., sigma_e2) of the general K-independent-random-effect model (see multi_effect_mme), by maximizing the dense REML log-likelihood (NelderMead over K+1 log-variances). effects is a vector of (Z_i, Ainv_i) pairs; initial, if supplied, is a length-K+1 vector of positive starting variances ([σ_1, …, σ_K, σ_e2]), otherwise all start at 1.
Returns a NamedTuple with variance_components = (sigmas, sigma_e2), per-effect ratios, beta, the K BLUPs (effects = [(ids, values), …]), loglik, converged, and per-component boundary flags (σ_i / total < 1e-6), plus #365 convention fields (loglik_convention = :reml_omit_2pi, loglik_full_constant_offset, loglik_comparable_across_routes = false). Raw loglik omits −½(n−p)log(2π) and is not comparable to fit_sparse_multi_effect_aireml / fit_multi_effect without comparable_loglik. Variance components are unaffected.
Reductions: the K=1 fit recovers the univariate animal-model REML optimum, and the K=2 fit is byte-identical to fit_two_effect_reml on identified data. EXPERIMENTAL, dense/validation-scale, REML-only, Gaussian, INDEPENDENT effects only (no correlated / direct–maternal covariance). The dense path forms an n×n V (guarded by max_dense_cells); it is an oracle for small K (≲4) and n (≲2000), NOT a production sparse estimator — that is the owed sparse AI-REML K-component path. On small/uninformative data the optimum can sit on a boundary (a variance → 0), reported via converged / boundary, never hidden. Non-identifiability (e.g. two A=I effects on the same grouping) can produce a flat ridge even when the optimizer reports converged = true. Optimizer convergence and boundary flags do not certify identifiability.
HSquared.fit_sparse_multi_effect_aireml Function
fit_sparse_multi_effect_aireml(y, X, effects; initial = nothing,
iterations = 100, tol = 1e-8, em_warmup = 0,
ids = nothing)Sparse average-information (AI) REML for the general K-independent-random-effect Gaussian animal model y = X·β + Σᵢ Zᵢ·uᵢ + e, uᵢ ~ N(0, σᵢ²·Aᵢ), e ~ N(0, σ_e²·I), estimating the K+1 variance components (σ₁²,…,σ_K², σ_e²). This is the sparse, scale-path generalization of the single-component fit_ai_reml; it is the production-shaped estimator behind the dense oracle fit_multi_effect_reml (which forms an n×n V and is guarded by max_dense_cells).
effects is a vector of (Zᵢ, Ainvᵢ) pairs (Zᵢ the n×qᵢ record→level sparse incidence, Ainvᵢ the qᵢ×qᵢ supplied relationship PRECISION — pass a sparse identity for a plain i.i.d. (1|group) effect), the same contract as multi_effect_mme. initial, if supplied, is a length-K+1 vector of positive starting variances [σ₁²,…,σ_K²,σ_e²] (default all 1), or :auto for a data-scaled start — half the response variance to the residual, the other half split evenly across the K effects. :auto uses the same objective from a different starting point; convergence and the endpoint must be checked for each fit. It pays off when the response is far from unit scale (measured 10 → 8 iterations, 0.67 s → 0.52 s, on an 11,856-record animal + permanent-environment fit with var(y) = 3.0) and costs iterations when it is not (12 against 10 on a unit-variance simulated fixture, test/). It is opt-in for a second reason: moving the optimizer path perturbs a converged estimate in its last bits and changes the reported iterations, and this package's recorded evidence — recovery checkpoints, pre-declared bias/MCSE gates — was run from (1,…,1).
Each iteration assembles the sparse Henderson MME coefficient matrix C = [X'R⁻¹X X'R⁻¹Z; Z'R⁻¹X Z'R⁻¹Z + blockdiag(Aᵢ⁻¹/σᵢ²)] (R = σ_e²·I), sparse-Cholesky factorizes it ONCE, reads each block's REML score ∂ℓ/∂σᵢ² = −0.5·(qᵢ − tr(Aᵢ⁻¹C^{uᵢuᵢ})/σᵢ² − uᵢ'Aᵢ⁻¹uᵢ/σᵢ²)/σᵢ² from the BLUP solution and the Takahashi selected inverse (selinv_block_traces, the tr(Aᵢ⁻¹C^{uᵢuᵢ}) terms — no dense inverse is formed), assembles the (K+1)×(K+1) average-information matrix from working-variate re-solves that reuse the same Cholesky factor, and takes an AI/Newton step with step-halving to keep every component positive. An interior fit reports convergence only when both the previous relative variance update and the evaluated-point score, scaled by the largest variance and residual degrees of freedom, are below tolerance. This is not a boundary KKT check. An optional EM-REML warm-start (em_warmup, default 0 = byte-identical to the pure AI path) hands the AI step a good in-bounds start.
REDUCTIONS: selected converged small-data K = 2 and K = 3 fixtures agree with dense fit_multi_effect_reml after aligning the REML likelihood constant; the K = 1 test reduces to fit_ai_reml (test/runtests.jl). These checks do not establish accuracy near a variance boundary. Returns a NamedTuple with variance_components = (sigmas, sigma_e2), per-effect ratios, beta, the K BLUPs (effects = [(ids, values), …]), loglik (LOGLIK_CONVENTION_FULL, identical to fit_ai_reml), loglik_convention = :reml_full_constant, loglik_full_constant_offset = 0.0, loglik_comparable_across_routes = true, converged, iterations, per-component boundary flags (σᵢ²/total < 1e-6), and estimator = :sparse_multi_effect_aireml. Raw dense fit_multi_effect_reml / fit_repeatability_reml loglik uses LOGLIK_CONVENTION_OMIT_2PI — convert with comparable_loglik before AIC / LRT across routes (#365).
EXPERIMENTAL, REML-only, Gaussian, INDEPENDENT effects only (no correlated / direct–maternal 2×2 G). The sparse machinery has selected dense reduction checks, but its scale/performance is NOT yet benchmarked (measure-first; sim/phase5_sparse_aireml_benchmark.jl is the opt-in scaffold) and it is NOT the public default fit path. On uninformative/non-identified data a component can approach the σ²→0 boundary; numerical score accuracy there is unresolved. When boundary-triggered streamed scores would exceed the column budget, the iterate takes an EM step instead of aborting (#441). If that EM step cannot stay finite and positive, status is boundary_score_unresolved and boundary is missing rather than a share-rule verdict from that iterate. Inspect converged before using a returned estimate.
HSquared.solve_multi_effect_pcg Function
solve_multi_effect_pcg(y, X, effects, sigmas, sigma_e2; tol = 1e-10, maxiter = 1000,
preconditioner = :jacobi, matrix_free = true, ids = nothing)Solve the supplied-variance K-INDEPENDENT-random-effect mixed-model equations by preconditioned conjugate gradient, matrix-free — the multi-effect generalization of solve_animal_model_pcg and the ITERATIVE companion of the direct multi-effect Cholesky in fit_sparse_multi_effect_aireml / sparse_multi_reml_loglik. It solves the SAME sparse SPD system C·[β; u₁; …; u_K] = rhs as the direct factorization (the coefficient matrix of _sparse_multi_lhs_rhs), never forming C.
effects is a vector of (Zᵢ, Aᵢ⁻¹) pairs (same contract as fit_sparse_multi_effect_aireml); sigmas are the K SUPPLIED random-effect variances and sigma_e2 the SUPPLIED residual variance (this is a supplied-variance SOLVE, not a REML fit — it does not estimate variance components).
matrix_free = true (default) applies C·v directly from the per-block sparse X, Zᵢ, Aᵢ⁻¹ matvecs (common = X·v_β + Σᵢ Zᵢ·v_{uᵢ}; top = X'·common/σe²; bottomᵢ = Zᵢ'·common/σe² + Aᵢ⁻¹·v_{uᵢ}/σᵢ²) with a matrix-free Jacobi diagonal — C is NEVER assembled, so the quadratic Cholesky fill-in of the environmental-group columns (the Phase 5 K≥2 scale bottleneck) is bypassed. matrix_free = false assembles C once (_sparse_multi_lhs_rhs) and applies it, enabling preconditioner = :ichol; both paths return the SAME solution.
preconditioner = :jacobi (default) uses M⁻¹ = 1/diag(C); :ichol (assembled only) an IC(0) factor with a diagonal-shift fallback; :none plain CG.
Returns a NamedTuple:
beta— fixed effects (1:pof the solution);effects— length-Kvector of(ids, values)per random block, in[u₁; …; u_K]order (blocki'sidsdefault to1:qᵢ, or the suppliedids[i]);iterations— CG iterations taken;relative_residual—‖rhs − C·x‖ / ‖rhs‖at the returned solution;converged— whetherrelative_residual ≤ tol;preconditioner,matrix_free.
EXPERIMENTAL — a CORRECTNESS primitive: validated to recover the direct multi-effect solve (β and per-block BLUPs) to a tight tolerance. It makes NO performance / large-q scaling claim (no benchmark is asserted here); it is the matrix-free iterative foundation the production sparse multi-effect path builds on. NOT the default fit path.
HSquared.mc_reml_block_traces Function
mc_reml_block_traces(X, effects, sigmas, sigma_e2; nprobe = 64, tol = 1e-9,
maxiter = 2000, seed = 0)MATRIX-FREE stochastic (Hutchinson) estimator of the K REML score-trace terms tr(Aᵢ⁻¹·C⁻¹[uᵢ,uᵢ]) of the K-INDEPENDENT-effect mixed-model equations — the building block for a matrix-free Monte-Carlo REML FIT (v0.8-S2 follow-on) at q → 10⁶, where the EXACT selinv_block_traces (a Takahashi selected inverse of the sparse Cholesky factor, used by fit_sparse_multi_effect_aireml) is factorization-limited by the K≥2 environmental-group Cholesky fill-in.
effects is the vector of (Zᵢ, Aᵢ⁻¹) pairs; sigmas/sigma_e2 are the SUPPLIED variance components at which the trace is evaluated (the same C as solve_multi_effect_pcg). For each block b: draw Rademacher probes z ∈ {±1}^{qᵦ}, embed into the full [β; u] space, solve C·x = ẑ MATRIX-FREE (the solve_multi_effect_pcg operator — C is never assembled or factorized), extract the b-block xᵦ = C⁻¹[uᵦ,uᵦ]·z, and accumulate the probe sample zᵀ·Aᵦ⁻¹·xᵦ. Since E[z zᵀ] = I, the probe MEAN is an UNBIASED estimate of the trace (E[zᵀ M z] = tr(M), M = Aᵦ⁻¹C⁻¹[uᵦ,uᵦ]).
Returns (traces, mcse) — each a length-K vector: the trace estimate and its Monte-Carlo standard error (probe SD / √nprobe; the honest noise band, ∝ 1/√nprobe). With one probe, mcse is undefined (NaN). Probe summaries are scaled to avoid intermediate overflow; unsupported nonfinite arithmetic throws an ArgumentError.
shared_probes (default false): with false, each block gets its own nprobe probes (nprobe·K solves). With true, one full-random probe over the whole random block yields all K block traces in nprobe solves. Both estimators are unbiased under exact solves. Shared probes can be cheaper, but cross-block terms contribute to their sampling variance, so they are not uniformly more precise at an equal solve budget. The returned mcse measures probe sampling variation and excludes PCG solve error.
EXPERIMENTAL — a stochastic-trace PRIMITIVE, NOT a fit. At validation scale the EXACT selinv_block_traces / fit_sparse_multi_effect_aireml is preferred (exact, no MC noise, few iterations); this exists ONLY for the large-q regime the direct factorization cannot reach. Deterministic given seed. No performance/scaling claim here (that is the pre-declared benchmark's job).
HSquared.fit_multi_effect_mc_reml Function
fit_multi_effect_mc_reml(y, X, effects; nprobe = 64, tol = 1e-4, iterations = 200,
seed = 0, pcg_tol = 1e-9, pcg_maxiter = 2000,
initial = nothing, ids = nothing)MATRIX-FREE Monte-Carlo REML fixed-point fit of the K-INDEPENDENT-effect Gaussian mixed model y = Xβ + Σᵢ Zᵢuᵢ + e, uᵢ ~ N(0, σᵢ²Aᵢ), e ~ N(0, σ²e I) — the FIT companion of the matrix-free solve solve_multi_effect_pcg and the estimator that lifts v0.8-S2 from a SOLVE to a FIT. It NEVER forms or factorizes the MME coefficient matrix C: each iteration is a matrix-free PCG solve plus the matrix-free Hutchinson trace mc_reml_block_traces, so it is not limited by the K≥2 environmental-group Cholesky fill-in that caps the direct/sparse AI-REML.
Uses REML score fixed-point updates with the form used by fit_sparse_multi_effect_aireml's warmup — σᵢ²_new = (ûᵢ'Aᵢ⁻¹ûᵢ + tr(Aᵢ⁻¹C⁻¹[uᵢ,uᵢ]))/qᵢ, σ²e_new = ê'ê/(n−p−Σqᵢ + Σ trᵢ/σᵢ²) — but with the trace terms MC-estimated. A FIXED probe seed is reused across iterations (correlated sampling), so the iteration is reproducible for a given seed. This is a stochastic REML fixed-point iteration, not an EM ascent guarantee; converged reports small relative changes in variance components when the Hutchinson trace_mcse is finite (nprobe > 1). It does not establish a unique or global optimum. optimizer_status is "fixed_point_relchange" for that stop, never "converged".
Returns a NamedTuple containing the fitted variance components, ratios, beta, random effects, and convergence diagnostics. trace_mcse is the Monte-Carlo standard error of the final trace estimate at the returned variance components. It is not a gradient, variance-component, or log-likelihood standard error, and it excludes PCG solve error. The accompanying trace_evaluation_variance_components records the variance components used for that diagnostic.
EXPERIMENTAL. The returned NamedTuple includes loglik; it is NaN unless compute_loglik = true, in which case matrix_free_reml_loglik supplies a stochastic REML log-likelihood estimate using SLQ for the log determinant (V8.1). At validation scale the EXACT fit_sparse_multi_effect_aireml is preferred (exact, no MC noise, fewer iterations); this exists for the large-q regime the direct factorization cannot reach. Current tests include a bounded exact AI-REML comparison (test/runtests.jl); they do not establish recovery rates. A pre-declared bias/MCSE recovery gate at scale + an external comparator are OWED before any covered claim.
HSquared.matrix_free_reml_loglik Function
matrix_free_reml_loglik(y, X, effects, sigmas, sigma_e2; slq_probes = 20, slq_steps = 40,
seed = 0, pcg_tol = 1e-9, pcg_maxiter = 2000)MATRIX-FREE Gaussian REML log-likelihood of the K-independent-effect model at the SUPPLIED variance components — the loglik companion of fit_multi_effect_mc_reml (v0.8 V8.1). It never forms or factorizes the MME coefficient matrix C: the only determinant term that needs C, log|C|, is estimated by stochastic Lanczos quadrature (SLQ) — slq_probes Rademacher probes, slq_steps matrix-free Lanczos steps each — while log|R| and log|G| are exact/cheap:
ℓ = −½[(n−p)·log(2π) + n·log(σ²e) + Σᵢ(qᵢ·log(σᵢ²) − log|Aᵢ⁻¹|) + log|C| + yᵀPy],
with log|Aᵢ⁻¹| from a one-time sparse Cholesky of each Aᵢ⁻¹ (the pedigree/identity precision — far cheaper than C; an O(q) Mendelian route is a future scale optimization). The log|C| term uses SLQ and yᵀPy uses a matrix-free PCG solve. It uses the SAME full-constant convention as sparse_multi_reml_loglik, but has both SLQ approximation error and PCG solve error.
Returns (loglik, loglik_mcse) — the estimate and the between-probe Monte-Carlo standard error for log|C|, scaled by ½ (∝ 1/√slq_probes). This MCSE does not measure finite-Lanczos quadrature bias or PCG solve error, so it does not alone quantify total likelihood error or support likelihood-ratio inference.
EXPERIMENTAL. The loglik is approximate and stochastic. At validation scale the exact sparse_multi_reml_loglik is preferred. Requires a positive-definite C; failed SLQ or nonfinite arithmetic throws an ArgumentError instead of reporting a numeric likelihood. With one probe, the probe MCSE remains undefined (NaN).
HSquared.matrix_free_reml_information Function
matrix_free_reml_information(y, X, effects, sigmas, sigma_e2; pcg_tol = 1e-11, pcg_maxiter = 5000)MATRIX-FREE average-information (AI) matrix of the K-independent-effect REML fit at the SUPPLIED variance components — the (K+1)×(K+1) information over [σ²₁, …, σ²_K, σ²e] used for asymptotic variance-component standard errors + ratio/h² intervals (v0.8 V8.3). It is EXACT (not stochastic): unlike the score's trace term, the AI matrix 0.5·WᵀPW is built purely from working-variate P-projections (W = [Z₁û₁/σ₁², …, Z_Kû_K/σ_K², ê/σ²e], with W ∈ ℝ^{n × (K+1)}, and P·w an MME re-solve), so it needs only matrix-free solves — NO Hutchinson trace. It reproduces the exact Cholesky-factor AI matrix (fit_sparse_multi_effect_aireml's information) to the PCG tolerance.
Returns the Symmetric (K+1)×(K+1) AI matrix. Nonfinite final information caused by unsupported arithmetic range throws an ArgumentError; finite information is not a conditioning or calibration guarantee.
HSquared.matrix_free_ratio_intervals Function
matrix_free_ratio_intervals(y, X, effects, sigmas, sigma_e2; level = 0.95,
boundary_tol = 1e-6, pcg_tol = 1e-11, pcg_maxiter = 5000)Asymptotic delta-method variance-ratio intervals for the matrix-free K-effect fit (v0.8 V8.3) — one per random effect (the animal-block ratio is narrow-sense h²; the others are variance-explained proportions). Builds the matrix-free AI matrix (matrix_free_reml_information) at the supplied estimate and applies the same logit delta-method (_ratio_delta_ci) as the exact multi_effect_ratio_interval. AI-based, ASYMPTOTIC, NOT coverage-calibrated; boundary components (σᵢ/total ≤ boundary_tol) are flagged. Returns a length-K vector of (estimate, lower, upper, se, lower_clamped, upper_clamped, boundary).
HSquared.fit_multi_effect Function
fit_multi_effect(y, X, effects; method = :auto, direct_max_n = 200_000, nprobe = 64,
verbose = true, kwargs...)Fit the K-INDEPENDENT-effect Gaussian mixed model, automatically choosing the solver by feasibility. This is the single entry point over the two multi-effect REML engines:
exact
fit_sparse_multi_effect_aireml— sparse AI-REML, exact gradient (a Cholesky selected inverse each iteration). The default where feasible; limited by the K≥2 environmental-group Cholesky fill-in at largeN.matrix-free
fit_multi_effect_mc_reml— Monte-Carlo REML fixed-point iteration, never forms or factorizesC(matrix-free solves + Hutchinson trace). Feasible where the exact path is fill-limited, with Monte-Carlo error in the trace estimates (trace_mcse). This is not a variance-component standard error.
method:
:auto(default) — route on feasibility::exactfor a single random effect (K = 1, the animal-model MME stays sparse-feasible to very largeq) OR whenN = p + Σqᵢ ≤ direct_max_n; otherwise:matrix_free, with a@infomessage (unlessverbose = false) noting the switch and that the trace estimates carry a Monte-Carlo standard error.:exact/:matrix_free— force the engine (may OOM / accept MC noise respectively).
The :auto direct_max_n threshold is a coarse, machine-agnostic heuristic calibrated to the measured K≥2 factorization cost (Phase 5 / v0.8-S2 benchmarks: the direct multi-effect Cholesky is already ~quadratic by q ≈ 50k and infeasible past ~10⁵–10⁶). It is deliberately conservative and ALWAYS overridable; a precise symbolic-fill predictor is future work. Both engines return loglik; it is exact for :exact, and NaN unless requested for :matrix_free, where it is a stochastic estimate. The matrix-free result also returns trace_mcse for trace-estimation error. The chosen engine's NamedTuple gains a dispatch field (:exact | :matrix_free) recording which ran.
EXPERIMENTAL. The multi-effect matrix-free result includes loglik (NaN unless requested, stochastic when enabled) and trace_mcse for trace-estimation error. These describe different quantities. nprobe, shared_probes, compute_loglik, slq_probes, slq_steps, and verbose are used only on the matrix-free route; supplying one of them on the :exact route is an error that names the unused keyword.
HSquared.multi_effect_ratio_interval Function
multi_effect_ratio_interval(y, X, effects; level = 0.95, which = :all,
initial = nothing, iterations = 200,
ids = nothing, fd_step = 1e-4,
boundary_tol = 1e-6)Asymptotic delta-method confidence interval(s) for the per-component variance ratios of the general K-independent-random-effect REML model (fit_multi_effect_reml) — the K-component generalization of two_effect_ratio_interval:
ratio_i = σ_i² / (Σ_{j=1}^{K} σ_j² + σ_e²), i = 1, …, K.Fits by REML, forms the observed information as the central finite-difference Hessian of the K-effect REML log-likelihood (_multi_effect_dense) over the K+1 variances at the optimum, and applies the delta method on the logit scale (so each interval lies in (0, 1)) — the SAME machinery as two_effect_ratio_interval / heritability_interval method = :delta, sharing the finite-difference information and per-ratio logit-delta helper.
which selects which component(s) to build an interval for: :all (default, all K) or an integer 1 ≤ i ≤ K (only component i; the others report their point estimate with a NaN interval). Returns a NamedTuple (ratios, level, converged) where ratios is a length-K vector of (estimate, lower, upper, se, lower_clamped, upper_clamped, boundary).
Boundary honesty (identical to the two-effect path): a component on the variance boundary (σ_i / total ≤ boundary_tol) is flaggedboundary = truewith aNaN interval — never a spuriously tight CI. A ratio built on a non-boundary numerator drops the degenerate component(s) and uses the corresponding information sub-block conditional on those components held at their supplied values. The reported ratio and its gradient retain the full denominator. At exact zero this gives the corresponding zero-component reduction; the approximation is not a near-boundary coverage guarantee.
Reductions: at K = 2 this matches two_effect_ratio_interval (same finite-difference information, same estimand); at K = 1 it matches heritability_interval method = :delta up to the finite-difference vs analytic-AI information difference (same estimand, different information estimator — a few-percent endpoint difference, not machine precision).
Experimental, asymptotic, delta-method, dense/validation-scale, REML only, Gaussian, INDEPENDENT effects only (no correlated / direct–maternal covariance); the dense path forms an n×n V. The interval is a large-sample approximation and is NOT coverage-calibrated — on small samples the REML surface is flat and the interval is unreliable (the parametric bootstrap, bootstrap_variance_component_interval, is the only finite-sample-aware path). No calibrated coverage is claimed.
HSquared.multi_effect_variance_component_covariance Function
multi_effect_variance_component_covariance(y, X, effects, sigmas, sigma_e2;
fd_step = 1e-4)Asymptotic covariance of the estimated [σ_1², …, σ_K², σ_e²] for a K-effect REML fit: the inverse of the observed information, formed as the central finite-difference Hessian of the REML log-likelihood at the estimate — the same machinery as multi_effect_ratio_interval and two_effect_ratio_interval, and the K-block analogue of variance_component_covariance.
Unlike multi_effect_ratio_interval, this differentiates the SPARSE sparse_multi_reml_loglik rather than the dense _multi_effect_dense, so it is usable on the same large problems fit_multi_effect(:auto) fits; it never densifies Aᵢ⁻¹.
Large-sample approximation, unreliable where the REML surface is flat. Throws rather than returning NaN when the information is not finite positive-definite (a flat or boundary optimum), or when a component sits so close to zero that the finite-difference step would take it non-positive. Experimental; REML only.
HSquared.multi_effect_variance_component_standard_errors Function
multi_effect_variance_component_standard_errors(y, X, effects, sigmas, sigma_e2;
fd_step = 1e-4)Asymptotic standard errors of [σ_1², …, σ_K²] and σ_e² for a K-effect REML fit, as a NamedTuple (sigmas, sigma_e2). See multi_effect_variance_component_covariance for the caveats.
HSquared.multi_effect_ratio_standard_errors Function
multi_effect_ratio_standard_errors(y, X, effects, sigmas, sigma_e2;
fd_step = 1e-4)Delta-method asymptotic standard errors of each variance RATIO ratioᵢ = σᵢ² / (Σⱼ σⱼ² + σ_e²) for a K-effect REML fit. For an animal block this is the standard error of narrow-sense h²; for any other block it is the standard error of that block's variance-explained proportion, which is NOT a heritability.
Uses multi_effect_variance_component_covariance and inherits its caveats and its refusals.
HSquared.multi_effect_sum_ratio_interval Function
multi_effect_sum_ratio_interval(y, X, effects, sigmas, sigma_e2;
which = 1:length(effects), level = 0.95,
fd_step = 1e-4, boundary_tol = 1e-6)Delta-method confidence interval for a SUMMED variance ratio
r = (Σ_{i ∈ which} σᵢ²) / (Σⱼ σⱼ² + σ_e²)of a K-effect REML fit, on the logit scale so the interval lies in (0, 1) — the same construction as _ratio_delta_ci's single-component case and as repeatability_interval, generalised to a sum of components.
With which = 1:2 on an animal + permanent-environment fit this is the REPEATABILITY coefficient t = (σ²_a + σ²_pe) / σ²_P. Unlike repeatability_interval, it takes the components from an already-computed fit and differentiates the SPARSE sparse_multi_reml_loglik, so it does not refit densely and carries no dense ceiling.
Returns a NamedTuple matching the single-ratio shape: estimate, lower, upper, se, lower_clamped, upper_clamped, boundary. Returns NaN endpoints with boundary = true rather than throwing when the ratio sits on a rail or the information is not positive definite. which must be a nonempty collection of unique integer component indices (Boolean values are invalid). Variances and fd_step must remain finite and positive after Float64 conversion; invalid inputs throw even when the ratio would lie on a rail. Asymptotic; REML only.
HSquared.multi_effect_uncertainty Function
multi_effect_uncertainty(y, X, effects, sigmas, sigma_e2;
which = 1:length(effects), level = 0.95,
fd_step = 1e-4, boundary_tol = 1e-6)Every asymptotic uncertainty product of a K-effect REML fit, from ONE information matrix. Returns a NamedTuple with
covariance— the(K+1)×(K+1)asymptotic covariance of[σ₁²,…,σ_K², σ_e²](multi_effect_variance_component_covariance);variance_component_se—(sigmas, sigma_e2), asmulti_effect_variance_component_standard_errors;ratio_se— the length-Kvector ofmulti_effect_ratio_standard_errors;sum_ratio_interval— the summed-ratio CI ofmulti_effect_sum_ratio_intervaloverwhich;level.
WHY IT EXISTS. Those three functions each build the covariance, and the covariance is a finite-difference Hessian of the REML log-likelihood — the single most expensive post-fit quantity in this package. A caller that wants all three (the R twin's target = "repeatability" sparse route wants exactly these three) was paying for three IDENTICAL information matrices. Measured on the great tit animal + permanent-environment fit: 4.40 s for the three separate calls, against 1.47 s for one covariance — and that 1.47 s is itself down from the pre-workspace cost (see _MultiREMLWorkspace).
The returned values are identical to calling the three functions separately; this only stops paying three times for one matrix. which has the same nonempty, unique integer selection contract as multi_effect_sum_ratio_interval. Invalid variances or finite-difference controls throw before information is computed.
Failure behaviour follows the individual functions. The covariance THROWS at a flat or boundary optimum, so this throws too, once, rather than reporting three separate failures for one cause. A summed ratio sitting on a 0/1 rail is different — the covariance is fine and only the logit interval is undefined, so sum_ratio_interval carries boundary = true with NaN endpoints while the standard errors are returned normally.
Experimental, REML only, asymptotic and NOT coverage-calibrated (the same fence as every component of it).
sourceHSquared.fit_direct_maternal_reml Function
fit_direct_maternal_reml(y, X, Zd, Zm, Ainv; initial = nothing,
iterations = 200, ids = nothing,
max_dense_cells = 1_000_000)REML estimation of the direct–maternal genetic model — one trait, one relationship A = Ainv⁻¹, a DIRECT additive effect (incidence Zd, record→animal) and a MATERNAL additive effect (incidence Zm, record→dam) with a 2×2 genetic covariance G_dm over [a_d; a_m] (Var = kron(G_dm, A)), plus residual σ²e. The 2×2 G_dm is estimated by dense REML over a log-Cholesky parameterization (so G_dm stays positive definite), reusing the multivariate _chol_params_to_cov machinery; the direct/maternal BLUPs come from the marginal GLS form.
Returns a NamedTuple with variance_components = (G_dm, sigma_ad, sigma_am, sigma_dm, sigma_e2), the direct–maternal genetic_correlation r_am, beta, direct_effects/maternal_effects (ids, values), loglik, and converged.
Reduction: with a diagonal G_dm (σ_dm = 0) _direct_maternal_dense equals the two-independent-effect dense model [(Zd, A), (Zm, A)]. COVERED at validation scale (opt-in, NOT the public default), dense (n ≤ ~1000), REML-only, Gaussian. This is the FIRST correlated random-effect structure (σ_dm ≠ 0), distinct from independent multi-effect and from multivariate G0-over-traits. Evidence: an R surface (target="direct_maternal" / maternal_genetic()), a sommer 4.4.5 covm() same-estimand REML comparator (AGREE, all entries ≤ 1.1e-2 rel.diff), and a PRE-DECLARED 48-seed bias/MCSE recovery gate (48/48 converged, all four |bias| ≤ 2·MCSE; see docs/dev-log/recovery-checkpoints/2026-07-01-direct-maternal-covered-evidence.md). INTERPRETATION FENCE (Willham): a negative r_am is real and expected; the direct heritability σ²_ad/σ²_P is NOT "the heritability" (the selection-relevant total additive variance involves σ_dm); callers must label direct-vs-total, never emit a bare h². On small/uninformative data or |r_am| → 1 the optimum can sit on a boundary (converged = false); identifiability generally needs designed data (the direct–maternal/PE/cytoplasmic confound).
HSquared.direct_maternal_interval Function
direct_maternal_interval(y, X, Zd, Zm, Ainv; level = 0.95, initial = nothing,
iterations = 200, ids = nothing, fd_step = 1e-4,
max_dense_cells = 1_000_000)Asymptotic delta-method standard errors and confidence intervals for the direct–maternal REML model (fit_direct_maternal_reml). Fits the model, forms the observed information as the central finite-difference Hessian of the dense REML log-likelihood over the four natural variance components θ = (σ²_ad, σ²_am, σ_dm, σ²e) at the optimum, inverts it for cov(θ̂), and delta-transforms to each reported quantity — the SAME finite-difference-Hessian + delta machinery as repeatability_interval.
Returns a NamedTuple with per-component (estimate, se, lower, upper) records for the variance components (sigma_ad, sigma_am, sigma_dm, sigma_e2), the direct–maternal genetic correlation r_am (Fisher-z interval, so it stays in (-1, 1)), and the Willham labelled triple direct_heritability (σ²_ad/σ²_P), maternal_ratio (σ²_am/σ²_P), and total_heritability (h²_T = (σ²_ad + 1.5·σ_dm + 0.5·σ²_am)/σ²_P, Willham (1972), with σ²_P = σ²_ad + σ²_am + σ_dm + σ²e — the SAME convention as the R total_heritability() surface).
INTERVALS ARE ASYMPTOTIC / UNCALIBRATED (normal-z Wald / delta on the observed REML information; Fisher-z for r_am), NOT coverage-calibrated — the same house convention as every other interval helper. Variance-component Wald bounds are NOT clamped (a lower bound may fall below 0 on a flat surface; that is honest and flagged by information_posdef). If the observed information is not positive definite (flat surface / boundary optimum, e.g. |r_am| → 1) the interval is undefined and an error is thrown, mirroring repeatability_interval. Internal / opt-in; this does NOT change the fitted result or any R-facing surface.
HSquared.NonGaussianFit Type
NonGaussianFitExperimental fitted-object container for the non-Gaussian animal model returned by fit_laplace_reml. Fields: variance_components, marginal_loglik, beta, breeding_values (the posterior-mode random effect), ids, converged, family, marginal, n_trials (the Binomial trials denominator for family = :binomial — a scalar common denominator or a per-record Vector{Int}; nothing for every other family), dispersion (the FIXED supplied overdispersion ρ for family = :beta_binomial; nothing for every other family — the negative-binomial theta is ESTIMATED and lives in variance_components, not here), boundary (true when the returned variance-component estimate sits at the search boundary — the bracket endpoint for the single-variance families, or the ±8-log-unit safety rail for the jointly-estimated families — meaning the point estimate is a function of the START VALUE, not the data; see #327), and restart_estimate (the second-start point estimate of sigma_a2 from an opt-in restart_check = true fit, nothing otherwise). Use the extractor functions breeding_values (→ BreedingValues(ids, values)), variance_components, and fixed_effects for the same access contract as AnimalModelFit; this is a distinct type so its extractors do not collide with the multivariate NamedTuple extractors.
HSquared.fit_laplace_reml Function
fit_laplace_reml(y, X, Z, Ainv; family = :gaussian, marginal = :laplace,
initial = nothing, ids = nothing, iterations = 200,
restart_check = false)Estimate the variance component(s) of the non-Gaussian animal model by maximising the integrated Laplace objective (marginal = :laplace) or the Gaussian variational objective (marginal = :variational) over the variance components. With fixed effects, non-Gaussian VA adds a Laplace beta correction and is not guaranteed to be a lower bound; see variational_marginal_loglik. family = :gaussian estimates (sigma_a2, sigma_e2) (NelderMead); family = :poisson, family = :bernoulli, and family = :binomial (which requires the n_trials keyword — a common scalar denominator OR a per-record integer vector of length length(y), the general cbind(successes, failures) GLMM) estimate the single sigma_a2 (Brent). family = :beta_binomial is the OVERdispersed logit-binomial (BetaBinomialResponse); it requires BOTH n_trials (scalar) and rho (the fixed overdispersion ρ ∈ (0,1)), estimates sigma_a2 (Brent) at that supplied fixed ρ, and is Laplace-only (marginal = :variational is rejected). A control this family does not use (rho outside :beta_binomial, n_trials outside :binomial/:beta_binomial, theta_init outside :nbinom/:gamma, initial.sigma_e2 outside :gaussian) is an error, not a silent drop. family = :bernoulli_probit is the binary threshold / liability-scale model (BernoulliProbitResponse, probit link Φ(η)); it estimates the single sigma_a2 (Brent) and is also Laplace-only (its variational expected information is response-dependent; marginal = :variational is rejected). Returns a NonGaussianFit with fields variance_components, marginal_loglik, beta, breeding_values, ids, converged, family, marginal, boundary, restart_estimate, and the extractor methods breeding_values(fit) / variance_components(fit) / fixed_effects(fit).
Every family's variance-component search is bounded, and converged = true alone does NOT mean the estimate is informative: a search that stops on its own boundary reports a point estimate that is a function of the START VALUE (initial), not the data (#327). fit.boundary reports this honestly for every family. Two mechanisms produce it: a Brent bracket endpoint at exp(log(sa0) ± 6) — the single-variance families (:poisson, :bernoulli, :binomial, :beta_binomial, :bernoulli_probit) and :ordered_probit with K = 2, where only σ²a is free — and a ±8-log-unit joint safety rail around the supplied start, for the jointly-estimated searches (:gamma, :ordered_probit with K ≥ 3, and now :gaussian and :nbinom, which had no rail at all before #327 and gained the one :gamma already used). The five-seed Bernoulli and Binomial recovery runs differ in their observed sigma_a2 estimates; that contrast does not establish the cause or a family-wide search-bound tendency (see sim/phase6_binomial_recovery.jl). nongaussian_three_field_payload (private) refuses a boundary = true fit.
When supplied, ids must contain exactly one unique identifier for each breeding value (size(Z, 2)).
restart_check = true (opt-in, doubles the cost of the fit) refits ONCE from a second start sa0 * exp(3.0) (hard-coded restart_check = false on that inner call, so there is no recursion) and compares the two sigma_a2 point estimates on the log scale: a gap > 0.01 sets boundary = true even when neither fit landed exactly on its own rail (the two-start fence itself was sized in the #327 wave-4 campaign comment — 27/27 truncated replicates detected, 0 false positives over 80 replicate pairs — but that measurement used :poisson only, a second start of `sa0
- 10
rather thansa0 * exp(3)`, and a relative rather than a log-scale gap; the
threshold shipped here is the same order of magnitude, not the measured configuration, and is unmeasured for the other eight families). The second estimate is exposed in restart_estimate (nothing when restart_check = false). This check is one-directional (#347): it can only turn boundary from false to true, never the reverse, so restart_check = true cannot clear an already-flagged boundary — it is a stricter detector, not a fix. The only lever that changes the point estimate itself is a different initial, which recentres the search bracket/rail.
EXPERIMENTAL, dense/validation-scale — the first fitted non-Gaussian step. For the Gaussian family the objective is the exact REML log-likelihood inside the ±8-log-unit box around the supplied initial that #327 added, so it recovers the same estimate as fit_sparse_reml whenever the optimum lies inside that box — which the default initial = (sigma_a2 = 1.0, sigma_e2 = 1.0) makes σ² ∈ [e⁻⁸, e⁸]. Outside it the search stops on the rail, returns a start-dependent estimate, and sets boundary = true; pass an initial on the scale of the data (or read fit.boundary) rather than assuming agreement. Exported as an experimental fitter; not the public default, not wired into the R formula path, no R model-spec, no external comparator.
HSquared.laplace_reml_interval Function
laplace_reml_interval(y, X, Z, Ainv; family = :poisson, marginal = :laplace,
level = 0.95, initial = nothing, n_trials = nothing,
iterations = 1000)Profile likelihood-ratio confidence interval for the single-variance-component non-Gaussian animal-model sigma_a2, by inverting the marginal LRT 2·(ℓ̂ − ℓ(sigma_a2)) ≤ χ²₁,level. Returns (sigma_a2, lower, upper, level, lower_clamped, upper_clamped, converged) — the *_clamped flags report whether an endpoint reached the search bound (the profile did not cross the χ² threshold within range) so a non-crossing endpoint is NOT a confidence limit. A point fit must converge and must not be flagged at its bounded search rail; otherwise the function throws instead of returning profile endpoints. Each profile evaluation must also converge and return a finite likelihood. For a returned interval, converged is true for the point fit.
Supports the single-variance-component families family = :poisson, family = :bernoulli, family = :bernoulli_probit, and family = :binomial (which requires n_trials — a scalar common denominator or a per-record integer vector, the same contract as fit_laplace_reml). The Gaussian two-component case needs nuisance profiling and is future work.
This is a profile-LIKELIHOOD-ratio interval, so marginal = :laplace is required: the variational :VA objective (including its Laplace beta correction when fixed effects are present) is not the marginal log-likelihood, so its differences are not χ²₁-calibrated — :variational throws rather than return an uncalibrated quantity dressed as a CI.
EXPERIMENTAL, asymptotic, single-component only — preliminary coverage CHARACTERIZATION only (sim/phase6_nongaussian_interval_coverage.jl, opt-in, validation-scale; observed coverage is cell-dependent with both under- and over-coverage, NOT a calibrated guarantee. Whether the interval is two-sided depends on where σ̂²a sits relative to the flat near-zero region of the profile, NOT on the family alone: a Binomial fit whose σ̂²a is clear of zero gives two interior LRT roots, but a small σ̂²a or the tested binary :bernoulli fixture can leave a flat profile so an endpoint reaches the search bound (flagged via *_clamped) — honest but not a confidence limit (the tested binary profile shape does not establish a general information limit or explain variance-estimation error). Reuses _profile_root.
HSquared.nongaussian_result_payload Function
nongaussian_result_payload(fit::NonGaussianFit)Bridge-ready, "boring" result payload (a NamedTuple of scalars / arrays / nested NamedTuples — Julia structs stay Julia-side) for a non-Gaussian animal-model fit from fit_laplace_reml. It mirrors the top-level shape of multivariate_result_payload (the univariate result_payload predates this convention and nests method/diagnostics differently) so the R twin can marshal one shape and the R non-Gaussian family-acceptance can fire.
Fields: engine, target = "nongaussian_reml", family ("gaussian"/"poisson"/"bernoulli"/"binomial"/"nbinom"/"beta_binomial"/"bernoulli_probit"), n_trials (the Binomial/beta-binomial trials denominator — a scalar common denominator or a per-record integer vector; nothing for other families — so a counts payload is self-describing on the data scale), dispersion (the FIXED overdispersion ρ for family = "beta_binomial"; nothing for every other family — the negative-binomial theta is ESTIMATED and rides in variance_components), method ("laplace"/"variational", resolved through the internal MarginalMethod dispatch from the stored marginal symbol), variance_components, fixed_effects, breeding_values = (ids, values), loglik, and converged.
The payload shape is deliberately family-uniform and therefore carries NO heritability field: a NonGaussianFit computes none, and h² is left to the consumer (it is derivable from variance_components for the Gaussian family, but on the liability scale the logit/log-link families have no residual-variance scale on which a single h² is defined here — surfacing one would be an unbacked claim). EXPERIMENTAL: the fitter is not the public default, not wired into the R formula path, and has no external comparator; the Bernoulli single-trial variance estimate has shown error in a five-seed recovery study; those seeds do not establish bias direction or cause. The R-facing method token ("laplace"/"variational") and family-acceptance shape are pending R-lane agreement before this is treated as a frozen contract.
HSquared.multivariate_mme Function
multivariate_mme(Y, X, Z, Ainv, G0, R0; ids = nothing, traits = nothing)Supplied-(co)variance Henderson solve of the multi-trait Gaussian animal model. For t traits, n records, and q related animals:
Y[i, :] = (X·B)[i, :] + (Z·U)[i, :] + E[i, :],
vec(Uᵀ) ~ N(0, A ⊗ G0), E[i, S_i] ~ N(0, R0[S_i, S_i]),where Y is n×t, the shared fixed-effect design X is n×p, the shared record→animal incidence Z is n×q, Ainv = A⁻¹ is the q×q relationship inverse, G0 is the t×t additive genetic covariance, and R0 is the t×t residual covariance. Records are ordered individual-major with trait fastest, so the mixed-model equations carry the genetic precision Ainv ⊗ G0⁻¹ on the random block and the residual precision R⁻¹ (block-diagonal over individuals) on the data:
[ X'R⁻¹X X'R⁻¹Z ] [vec(Bᵀ)] [ X'R⁻¹·y ]
[ Z'R⁻¹X Z'R⁻¹Z + Ainv⊗G0⁻¹ ] [vec(Uᵀ)] = [ Z'R⁻¹·y ]Unbalanced / missing-trait records are supported: an entry of Y that is missing or NaN marks an unobserved trait for that record. Such observations are dropped from the data and individual i's residual precision uses only the observed-trait submatrix inv(R0[S_i, S_i]) (where S_i is its observed-trait set). With every trait observed this reduces to the balanced model with residual precision I_n ⊗ R0⁻¹. A non-finite observed phenotype (Inf — only missing / NaN mark an unobserved trait), a non-finite X/Z/Ainv, or a trait with no observed records is rejected with a clear ArgumentError.
This is the Phase-4 (multivariate Gaussian) supplied-variance engine slice: an MME solve that does not estimate G0 / R0 (covariance-matrix estimation is a separate REML slice), the multi-trait analogue of henderson_mme. It assumes a fixed-effect / incidence design shared across traits; per-trait designs are not covered. Experimental, dense/validation-scale, supplied-covariance Julia helper. R's multivariate fitting routes are documented separately.
Returns a NamedTuple:
beta—p×tfixed effects (columnk= traitk);breeding_values—(ids, traits, values)withvaluestheq×tEBV matrix (columnk= traitk);genetic_covariance,residual_covariance— the suppliedG0,R0;genetic_correlation,residual_correlation— derivedt×tcorrelations;traits— the trait labels.
HSquared.fit_multivariate_reml Function
fit_multivariate_reml(Y, X, Z, Ainv; initial = nothing, iterations = 2000,
ids = nothing, traits = nothing,
genetic_structure = :unstructured, rank = nothing)Estimate the genetic and residual covariance matrices G0, R0 of the multi-trait animal model by dense REML. Inputs match multivariate_mme (balanced or with missing/NaN trait records); the marginal model is
y ~ N(X·β, V), V = Z_full·(A ⊗ G0)·Z_full' + R,with R block-diagonal over individuals (individual i's block R0[Sᵢ, Sᵢ]). The REML log-likelihood is maximized by Nelder–Mead. The default :unstructured fit uses a log-Cholesky parameterization for both G0 and R0; structured fits constrain only the additive genetic covariance and keep R0 unstructured positive definite. At the optimum the breeding values and fixed effects are computed from the marginal GLS BLUP form, which remains defined for low-rank genetic covariances.
The additive genetic covariance can be constrained with genetic_structure:
:unstructured(default) — full positive-definiteG0;:diagonal— independent trait-specific genetic variances;:lowrank—G0 = ΛΛ', requiringrank;:factor_analytic—G0 = ΛΛ' + Ψ, requiringrank. Fitted uniqueness usesψ_i = FA_UNIQUENESS_FLOOR + exp(θ_i)somin(ψ̂) ≥ 1e-4(Heywood interior bound in absolute trait-specific genetic-variance units). This frozen floor is not scale-equivariant near the bound: changing a trait's units bycchanges its variance byc²but does not change1e-4. Ledermann-saturated cells (ledermann_slack(t, K) ≤ 0) may still be fitted for diagnosis;require_fa_covered_flip_cellrefuses them as a covered-flip cell.
For structured fits, initial may supply G0, R0, loadings (traits × rank), and uniqueness for :factor_analytic; omitted fields use deterministic phenotypic-scale defaults. When genetic_structure = :factor_analytic and initial is omitted, the fitter runs both its standard and a deterministic balanced phenotypic-scale initialization. Away from the uniqueness floor, the second start preserves the initial genetic diagonal, halves its off-diagonal covariances, and retains the initial residual covariance; floor clamping can also change the genetic diagonal. This is a decomposition sensitivity check, not broad multistart protection. It returns the highest-likelihood converged fit when one exists; if neither converges, it returns the best finite fit with converged = false. fa_start_diagnostics records both outcomes, whether their fitted covariances disagree, and whether uniqueness is within 1% of the absolute floor. Supplying initial requests a single user-controlled start.
Returned genetic_loadings use a deterministic sign convention: for each factor column, the largest-absolute loading is non-negative. This removes arbitrary sign flips from returned metadata but does not impose a rotation or lower-triangular identification constraint; for rank > 1, loadings remain rotation-nonunique and should not be interpreted as uniquely identified factors.
The optimizer still searches the full raw loading vector (t * rank entries for :lowrank, plus t uniqueness parameters for :factor_analytic). It does not remove the K(K-1)/2 rotational degrees of freedom of Λ. Those flat directions do not change G0 = ΛΛ' (or ΛΛ' + Ψ). _mv_nparams subtracts them for LRT reporting only; converged and the search itself still use the unreduced parameter count (HSquared.jl #340).
Experimental, dense/validation-scale, REML-only, Gaussian. The REML estimator is validated by deterministic self-consistency checks (the t = 1 reduction recovers the univariate REML estimate; the multivariate REML log-likelihood is on the same full-constant scale as the univariate one; the optimum beats a coarse grid). Seeded fixtures and separately documented recovery campaigns concern specific trait/rank/design cells. They establish no broad recovery or inference calibration. For the unstructured multi-trait comparison, same-estimand external parity is discharged on ONE deterministic fixture (point-estimate) by sommer 4.4.5 and blupf90+ 2.60, with no ASReml or JWAS leg — treat multi-trait variance estimates as experimental. There is no published textbook anchor for estimated multivariate G0/R0: the R lane's Mrode Example 5.1 leg is a supplied-covariance BLUP/MME anchor, so it does not check the estimated covariances. The R twin reaches this fitter on its default cbind() route (hsquared MV-4); that claim stays partial, and R↔engine element-wise parity at the promotion fixture is discharged locally (hsquared A26) but is not CI-backed — the R lane's CI provisions no Julia, so the parity legs skip there.
Returns a NamedTuple:
genetic_covariance,residual_covariance— estimatedG0,R0(t×t);genetic_correlation,residual_correlation— derivedt×tcorrelations;heritability— per-traith² = diag(G0)/(diag(G0)+diag(R0));beta,breeding_values— fixed effects and EBVs at the estimate;genetic_structure,genetic_rank,genetic_loadings,genetic_uniqueness— structure metadata (nothingwhere not applicable);loglik— the full REML log-likelihood at the estimate, including the(N − p')·log(2π)constant, so it is on the same scale asgaussian_loglik/sparse_reml_loglik; comparisons still require matched data and compatible fixed-effect spaces;converged,iterations,traits;fa_start_diagnostics— for FA fits, per-start objective, termination, iteration count, objective range, covariance disagreement, whether a better nonconverged start exists, and distance to the uniqueness floor;nothingfor other genetic structures.
Non-finite or empty-trait inputs are rejected up front (see multivariate_mme), so the optimizer never returns plausible-looking covariances from Inf/NaN data; check converged for genuine boundary cases.
HSquared.multivariate_result_payload Function
multivariate_result_payload(result)Bridge-ready, "boring" result payload (a NamedTuple of scalars / arrays / nested NamedTuples — Julia structs stay Julia-side) for a multivariate REML fit from fit_multivariate_reml. It mirrors result_payload for the univariate animal model so the R twin can marshal one shape across traits.
Exposed only for the rotation-free genetic structures :unstructured and :diagonal. :lowrank / :factor_analytic are rejected on purpose: their loadings are rotation-nonidentified, so surfacing them across the bridge needs a rotation/interpretation convention first. The payload carries the estimated G0/R0, the per-trait genetic variances diag(G0), the genetic/residual correlations, the REML loglik, the genetic parameter count n_genetic_params (so a :diagonal-vs-:unstructured LRT df is just the difference of the two fits' counts), the fixed effects, breeding values, per-trait heritabilities, and converged. It deliberately omits genetic_loadings / genetic_uniqueness.
HSquared.multivariate_covariance_standard_errors Function
multivariate_covariance_standard_errors(fit, Y, X, Z, Ainv; fd_step = 1e-4)Asymptotic standard errors for an unstructured multivariate REML fit (fit_multivariate_reml). Pass the same Y, X, Z, Ainv used for the fit. Standard errors come from the observed-information matrix (a central finite-difference Hessian of the REML log-likelihood on the log-Cholesky parameterization, the scale the optimizer uses) propagated through the delta method to the genetic/residual covariances, the derived correlations, and the per-trait heritabilities.
Returns a NamedTuple of t×t SE matrices genetic_covariance, residual_covariance, genetic_correlation, residual_correlation (correlation-SE diagonals are 0), the length-t heritability SE vector, and the raw information matrix.
Requires converged = true and a finite positive representable difference step. Experimental, asymptotic, dense/validation-scale, REML-only. SEs are finite-difference approximations and are wide/unreliable at small n; not coverage-calibrated. Throws if the observed information is not finite positive-definite (a flat or boundary optimum). A fitted off-diagonal |r| ≥ 1 − 1e-6 is treated as a boundary even when the finite-difference Hessian appears PD (platform/roundoff). Structured/factor-analytic fits are not supported by this unstructured-coordinate SE implementation; rotation-invariant covariance quantities would need a separate identified parameterization and validation.
HSquared.covariance_structure_lrt Function
covariance_structure_lrt(constrained, full)Likelihood-ratio test comparing a nested constrained genetic-covariance fit against the full (less-constrained) fit, both from fit_multivariate_reml on the same data. Returns a NamedTuple with the LRT statistic = 2(ℓ_full − ℓ_constrained), the parameter-count difference df, the pvalue, a reference symbol naming which reference distribution produced it, a boundary flag, and a note. Both fits must report convergence, finite log-likelihoods, and equal trait dimensions. Only :diagonal, :lowrank, or :factor_analytic nested in :unstructured are accepted. The caller remains responsible for confirming that both fits use the same data, fixed-effect design, and relationship matrix. df is the nominal generic dimension difference: for :lowrank/:factor_analytic, _mv_nparams removes the r(r-1)/2 rotational indeterminacy of Λ (Λ and ΛQ give the same G). This count alone does not establish local identifiability at a particular loading matrix or regular LRT asymptotics.
Two kinds of structured null are distinguished. For neither of them is nested_lrt's chi-bar weighting invoked: this function always requests the plain, unmixed tail (boundary_df = 0), because nested_lrt's convex-cone (boundary_df ≥ 1) branches are for genuine variance-at-zero boundaries and neither structured null here is one:
No PSD-rank boundary flagged (
:diagonalor:factor_analyticnested in:unstructured,reference = :chisq,boundary = false): the plain χ²dftail is reported. For FA it is a regular asymptotic reference only at a locally identifiable, interior point where the covariance map has full local rank and the usual likelihood conditions hold. Rotation-quotient counting andΨ > 0do not prove those conditions. For example, att = 4, K = 1, four nonzero loadings yield Jacobian rank 8, whereas only two nonzero loadings leaveΨnonidentified even with positive slack.fit_multivariate_remluses the absolute genetic-variance floorψ_i = 1e-4 + exp(θ_i); an estimate close to the floor is another nonregular case. This helper checks neither local rank nor floor distance, soboundary = falseis not a regularity certificate. The 50:50 chi-bar correction for a single variance-at-zero boundary is not applied.PSD-boundary nulls (
:lowranknested in:unstructured,reference = :chisq_naive_boundary,boundary = true): a low-rank nullG = ΛΛ'(rankr < t) lies on the boundary of the PSD cone. The regular χ² reference is not justified by this helper; the relevant nonstandard boundary limit and its weights are not derived or computed here. The reported p-value is the naive χ²dftail, and its direction relative to the relevant limit law is not knowable here — never call it conservative.
Changed in #331: because boundary_df = 0 is now always requested, a structured comparison whose df is 1 no longer receives the 50:50 chi-bar mixture it received before. The :lowrank-in-:unstructured t = 2, rank = 1 case is the reachable instance: its reported p-value is now the full χ²₁ tail, exactly twice the previously reported value. The statistic and df are unchanged; only the reference distribution is.
Experimental, asymptotic, dense/validation-scale. Repeatability fits include a permanent-environment covariance. This animal-only parameter count does not cover them, so either input is rejected.
sourceHSquared.genetic_correlation Function
genetic_correlation(C::AbstractMatrix)
genetic_correlation(result)Convert a symmetric covariance matrix C (e.g. the additive genetic covariance G0 or residual covariance R0) into the corresponding correlation matrix D⁻¹ C D⁻¹, with D = diag(sqrt.(diag(C))). Diagonal entries are exactly 1.
The second method extracts genetic_correlation from a multivariate_mme result NamedTuple. Requires positive diagonal variances.
HSquared.genetic_structure Function
genetic_structure(result::NamedTuple)Return structured genetic covariance metadata from a multivariate REML result as (structure, rank).
HSquared.genetic_loadings Function
genetic_loadings(result::NamedTuple)Return a copy of the structured genetic loading matrix from a low-rank or factor-analytic multivariate REML result. Returns nothing when the fitted structure has no loading matrix.
HSquared.genetic_uniqueness Function
genetic_uniqueness(result::NamedTuple)Return a copy of the structured genetic uniqueness vector from a multivariate REML result. Returns nothing when the fitted structure has no uniqueness metadata.
HSquared.evolvability Function
evolvability(G, beta)Hansen & Houle (2008) evolvability e(β) = β̂ᵀ G β̂ — the additive genetic variance available in the (unit-normalized) selection-gradient direction β. G is a genetic covariance matrix (or a multivariate result, from which genetic_covariance is read). PSD-safe (works on a reduced-rank G). Along an eigenvector of G it equals that eigenvalue. Descriptive geometry, not a predicted response; see the module note for caveats.
HSquared.conditional_evolvability Function
conditional_evolvability(G, beta)Hansen & Houle (2008) conditional evolvability c(β) = 1 / (β̂ᵀ G⁻¹ β̂) — the genetic variance available in direction β when all other directions are held at their conditional optima (i.e. under a constraint). Requires G positive DEFINITE (it inverts G); singular or excessively ill-conditioned matrices throw. The condition estimate must be at most 1/sqrt(eps(Float64)). Along an eigenvector of G it equals that eigenvalue, and c(β) ≤ e(β) always. The unit direction is Euclidean in the supplied trait coordinates; values depend on trait scaling.
HSquared.respondability Function
respondability(G, beta)Hansen & Houle (2008) respondability r(β) = ‖G β̂‖ — the magnitude of the evolutionary response to a unit selection gradient in direction β. PSD-safe. Along an eigenvector of G it equals that (non-negative) eigenvalue.
HSquared.autonomy Function
autonomy(G, beta)Hansen & Houle (2008) autonomy a(β) = c(β) / e(β) ∈ (0, 1] — the fraction of the evolvability in direction β that is autonomous of genetic constraint from other traits. Requires G positive definite (via conditional_evolvability). Along an eigenvector of G it equals 1.
HSquared.variance_along_gradient Function
variance_along_gradient(G, beta; normalize = true)Additive genetic variance of the linear index βᵀ a, βᵀ G β. With normalize = true (default) β is scaled to unit length first, so this returns evolvability; with normalize = false it uses β as given (the raw genetic variance of the index for an arbitrary contrast). PSD-safe. Derived variances outside finite Float64 range throw ArgumentError.
HSquared.genetic_pca Function
genetic_pca(G)Genetic principal-component decomposition of G. Returns (values, vectors) with eigenvalues in DESCENDING order (the genetic variances along the genetic principal axes / "genetic lines of least resistance") and the corresponding eigenvectors as the columns of vectors, each deterministically sign-canonicalized (largest-magnitude entry positive). PSD-safe. Under repeated eigenvalues, individual eigenvectors within the repeated eigenspace are not unique. Interpret that eigenspace rather than an individual axis. This also applies to the leading axis when the largest eigenvalue is repeated. Individual axes can be unstable when eigenvalues are close. Eigenvectors and eigenvalues are expressed in the supplied trait coordinates; changing trait units can change both.
HSquared.g_max Function
g_max(G)Leading genetic principal axis of G: returns (eigenvalue, eigenvector) for the largest eigenvalue — g_max, the direction of maximum additive genetic variance ("genetic line of least resistance"). The eigenvector is sign-canonicalized. When the largest eigenvalue is repeated, g_max is not a unique direction; interpret the full leading eigenspace instead.
HSquared.mean_evolvability Function
mean_evolvability(G)Mean (unconditional) evolvability — the average of e(β) over uniformly random unit selection gradients, which has the exact closed form tr(G) / t (the mean eigenvalue of G). PSD-safe, deterministic. The population-averaged conditional evolvability and autonomy (Hansen & Houle "random skewers") have no comparable simple closed form and are left to future work.
HSquared.genetic_gllvm_descriptors Function
genetic_gllvm_descriptors(loadings; uniqueness = nothing)Rotation-invariant descriptors of a genetic-GLLVM latent layer with SUPPLIED traits × K loadings Λ (and optional positive traits-vector uniqueness Ψ). The implied among-trait genetic covariance is Σ_g = ΛΛ' (low-rank, uniqueness = nothing) or Σ_g = ΛΛ' + diag(Ψ) (factor-analytic). Returns a NamedTuple:
genetic_covariance—Σ_g;genetic_variances—diag(Σ_g);genetic_correlation— the correlation matrix ofΣ_g;communality—c²_t = (ΛΛ')_tt / Σ_g[t,t]∈[0,1], the per-trait fraction of genetic variance explained by the common latent factors (= 1whenΨis absent; the one genuinely new GLLVM descriptor);genetic_pca—(values, vectors)ofΣ_g(descending eigenvalues, sign-canonicalized eigenvectors);g_max— leading genetic principal axis ofΣ_g;rank/n_latent_factors— the latent-factor countK = size(Λ, 2).
DESCRIPTIVE, supplied-covariance only: Λ/Ψ are NOT estimated by this descriptor function, which computes no marginal / likelihood / fit or bridge payload. The genetic covariance, variances, correlations, and PCA eigenvalues are functionals of Σ_g; communality additionally depends on the supplied decomposition through Ψ. At fixed Ψ, an orthogonal rotation Λ → ΛQ leaves communality unchanged, but decompositions with the same Σ_g and different Ψ need not have the same communality. Raw loadings Λ are never returned (they are rotation-nonidentified). The genetic_pca eigenvectors are defined only up to sign, and a repeated eigenvalue identifies its eigenspace rather than a unique axis. Guards (dimension / positivity / rank) are delegated to lowrank_covariance and factor_analytic_covariance. The first foundation step of the genetic GLLVM (#50); the supplied-covariance latent objective and estimation were added in later slices below.
genetic_gllvm_descriptors(result::NamedTuple)Rotation-invariant genetic-GLLVM latent-structure descriptors for an ESTIMATED factor-analytic or low-rank multivariate REML fit (fit_multivariate_reml(...; genetic_structure = :factor_analytic | :lowrank, rank = K)). Reads the fit's estimated, rotation-invariant genetic covariance G = result.genetic_covariance and uniqueness Ψ (genetic_uniqueness; nothing ⇒ low-rank, Ψ = 0) — NEVER the rotation-nonidentified loadings — and returns the same NamedTuple as the supplied-loadings method, with communality = 1 − Ψ / diag(G) (the per-trait fraction of genetic variance from the common latent factors; = 1 for low-rank). For FA, this communality is conditional on the fitted G/Ψ decomposition; equal G matrices can yield different communalities. Synthetic results must carry consistent factor-count/loading dimensions and finite nonnegative uniqueness with Ψ <= diag(G). Loading values are checked when supplied, but descriptors are computed from G and Ψ, preserving the rotation-free output contract. Rejects the rotation-free :diagonal / :unstructured structures, which have no latent-factor interpretation.
HSquared.genetic_gllvm_gaussian_mme Function
genetic_gllvm_gaussian_mme(Y, X, Z, Ainv, loadings, R0; uniqueness = nothing, ids = nothing)Supplied-covariance Gaussian genetic-GLLVM latent solve (#50 slice 2). With a Gaussian response, the genetic-GLLVM latent layer η[i,t] = Σ_k Λ[t,k] g[i,k], g[·,k] ~ N(0, A) and, when supplied, trait-specific effects d[·,t] ~ N(0, Ψ[t,t] A) make the trait-level effect U = FΛ′ + D and the among-trait genetic covariance G_lat = ΛΛ' (+ diag Ψ) satisfy Cov(vec(U)) = G_lat ⊗ A — i.e. the Gaussian genetic GLLVM is EXACTLY the multivariate animal model at G0 = G_lat. This convenience builds G_lat from the SUPPLIED traits × K loadings Λ (+ optional positive uniqueness Ψ) and solves it through multivariate_mme, returning that solve (beta, breeding_values, genetic_covariance = G_lat, residual_covariance, genetic_correlation, residual_correlation, traits) augmented with the rotation-invariant latent_structure (genetic_gllvm_descriptors) and n_latent_factors = K.
G_lat must be positive definite for the multivariate genetic precision G0⁻¹ to exist: supply a positive uniqueness Ψ, or full-rank loadings (K ≥ traits). A pure low-rank G_lat (K < traits, no Ψ) is singular and is rejected by the solve. SUPPLIED-covariance only — Λ/Ψ/R0 are NOT estimated (that is slice 3), and only rotation-INVARIANT functionals of the latent structure are reported (never raw loadings). No R model-spec or bridge payload.
HSquared.genetic_pca_plot_data Function
genetic_pca_plot_data(G; n_axes = nothing)Plot-ready data for the rotation-invariant G-geometry figure (scree + biplot): delegates to genetic_pca and returns (eigenvalues, variance_explained, eigenvectors, loadings_scaled, axis_labels, rotation_invariant = true, is_eigenstructure_not_loadings = true). HARD CONTRACT: only the rotation-INVARIANT eigenstructure of G — eigenvalues, sign-canonicalized principal axes, and loadings_scaled = eigenvectors · √eigenvalue (biplot vector length) — is returned, NEVER raw factor-analytic loadings Λ. n_axes (≤ size(G,1)) limits the axes.
HSquared.genetic_correlation_plot_data Function
genetic_correlation_plot_data(G; traits = nothing, heritabilities = nothing)Plot-ready data for the genetic-correlation heatmap: delegates to genetic_correlation and returns (traits, genetic_correlations, heritabilities, rotation_invariant = true). Rotation-invariant (D⁻¹ G D⁻¹, unit diagonal, off-diagonals in [-1,1]). Optional heritabilities lets the drawing layer flag low-h² / imprecise cells.
HSquared.variance_components_plot_data Function
variance_components_plot_data(fit::AnimalModelFit; level = 0.95)Plot-ready data for the variance-component + heritability forest figure (plotting set B): tidy parallel vectors (term, estimate, lo, hi, panel, level, interval_method, interval_status, interval_reason, supplied = false) shaped to drop directly into the R hs_gg_forest contract. The variance-component rows (sigma_a2, sigma_e2) carry asymptotic estimate ± z·SE — NOT clamped, since an asymptotic CI can cross zero (surfaced, never hidden); the h2 row carries the logit-delta heritability_interval (always in (0,1)). lo/hi are NaN where the interval is unavailable (no fabricated whiskers); interval_reason records "ok" or the expected numerical/argument failure for each row. Unexpected exceptions, including interrupts and programming errors, are rethrown. interval_status is "experimental_asymptotic" (NOT coverage-calibrated) when any interval is present, else "none". interval_method is a coarse roll-up tag ("asymptotic_reml"): the VC-row whiskers are normal-Wald on the raw variance scale, the h2-row whisker is the logit-delta back-transform — both asymptotic, both from the REML AI matrix. supplied = false is the honest-status hinge — these are ESTIMATED, unlike the descriptive supplied-K_g/G plot-data sets. Intervals are REML-only; a non-REML fit degrades gracefully to points-only (lo/hi all NaN, interval_status = "none").
HSquared.breeding_values_plot_data Function
breeding_values_plot_data(fit::AnimalModelFit; trait = 1)Plot-ready data for the EBV "caterpillar" figure (plotting set B): tidy parallel vectors (id, trait, value, pev, pev_scale) shaped to drop directly into the R autoplot.R breeding-value plot (per the #93 R-twin alignment — this closes the last live-parity gap R flagged). value is the EBV (breeding_values), pev the prediction error variance (prediction_error_variance, sparse :selinv path since #350 for a sparse Ainv; a dense genomic Ginv keeps the dense path), and pev_scale = "validation" is the honest-status flag: the PEV is validated against the dense oracle at validation scale (V1-SELINV-PEV), NOT yet a production large-pedigree reliability claim. The R column convention is followed exactly (EBV as value). Univariate AnimalModelFit; trait is the (single) trait label. Plot-DATA only — no drawing backend, no estimation.
HSquared.hsquared_figure Function
hsquared_figure(data; kind = <inferred>, kwargs...)Draw an HSquared plotting-layer figure from a *_plot_data NamedTuple. STUB: the drawing method is provided by the HSquaredMakieExt package extension, which activates only when a Makie backend and AlgebraOfGraphics are loaded (using CairoMakie, AlgebraOfGraphics / using GLMakie, AlgebraOfGraphics). Without them, calling this throws a MethodError asking you to load the plotting weak dependencies.
Supported kinds (each consumes the matching preparer and renders the R-twin's honest-status drawing behaviors ON the figure — #93):
:variance_components←variance_components_plot_data— the VC + h² forest (AoG tabular/faceted scaffold, Makie annotations). Whiskers are RAW (never clamped); a whisker crossing the[0,1]boundary is annotated on the h² panel only (a variance-component whisker crossing 0 is expected/honest, not flagged);NaNwhiskers draw no whisker.:breeding_values←breeding_values_plot_data— the EBV caterpillar (AoG tabular scaffold, Makie uncertainty whiskers) (sorted EBV ±√PEV); thepev_scale = "validation"caveat is rendered in the subtitle.:g_geometry←genetic_pca_plot_data— the eigenvalue scree only, NEVER a loadings biplot (gated onis_eigenstructure_not_loadings); on a non-PDGa negative eigenvalue bar is drawn but the %-variance labels are suppressed and a "non-positive-definite G" note is added (a ">100% variance share" is the trap).:genetic_correlation←genetic_correlation_plot_data— the rotation-invariantD⁻¹GD⁻¹correlation heatmap (unit diagonal, off-diagonals in[-1,1]), NEVER raw loadings (gated onrotation_invariant); whenheritabilitiesare supplied, low-h² (imprecise) traits are flagged in the subtitle.:manhattan←marker_manhattan_data— chromosome-coloured scatter of cumulativeplot_positionsvs-log10(p), with a VISUAL-ONLY Bonferroni guide line; the subtitle states the p-values are nominal Wald and NOT genome-wide calibrated (#48), so the line is guidance only.:qq←marker_qq_data— observed vs expected-log10(p)with they = xuniform-null line; the subtitle carries the same NOT-genome-wide-calibrated caveat. λGC is intentionally not recomputed in the drawing layer (the preparer carries no χ²; recomputing would duplicate/srcnumerics and risk an uncalibrated read).:rr_variance←rr_genetic_variance_plot_data— the genetic-variance trajectoryv_g(t), plus anh²(t)panel ONLY when a residual was supplied; the subtitle flags it as supplied-K_gdescriptive and thath²(t)can overstate without a permanent-environment term.:rr_surface←rr_covariance_surface_plot_data— the genetic covariance/correlation surface heatmap (diverging RdBu centred at 0); the correlation surface uses a fixed(-1, 1)colorrange, the covariance surface a data-driven symmetric-about-0 range (a fixed(-1, 1)would clip).:rr_eigenfunctions←rr_eigenfunctions_plot_data— one covariance-function eigenfunction per line with a per-axis variance-explained legend; the subtitle states it is rotation-invariant, signs are arbitrary, and the span is ambiguous under repeated eigenvalues.
The honest-status caveat is always rendered on the figure (subtitle), sourced from the SAME flags the preparer carries — this is the drawing-layer half of the plotting layer's "subtitle drop is the only guardrail" defense (docs/design/13-plotting-layer.md). Drawing only — no estimation, no engine computation.
HSquared.legendre_basis Function
legendre_basis(t, order)Normalized Legendre basis vector φ(t) = [φ_0(t), …, φ_{order-1}(t)] at a standardized covariate t ∈ [-1, 1], where φ_n(t) = sqrt((2n+1)/2)·P_n(t) and P_n are the ordinary Legendre polynomials (Bonnet recurrence). The basis is orthonormal on [-1, 1] (∫_{-1}^{1} φ_m φ_n = δ_mn). order = k is the number of random-regression coefficients. Use standardize_covariate to map a raw covariate onto [-1, 1] first (this function throws if |t| > 1).
HSquared.legendre_design Function
legendre_design(ts, order)Build the n × order random-regression design matrix Φ whose row i is legendre_basis(ts[i], order) for standardized covariate points ts ∈ [-1, 1] (use standardize_covariate first). The columns are the normalized Legendre basis functions.
HSquared.random_regression_mme Function
random_regression_mme(y, X, Phi, Z, Ainv, K_g, sigma_e2; ids = nothing)Solve the SUPPLIED-covariance Henderson mixed-model equations for the polynomial random-regression animal model with a homogeneous scalar residual variance:
y_r = x_rᵀ β + φ(s_r)ᵀ a_{a(r)} + e_r, e_r ~ N(0, sigma_e2),
vec(a) ~ N(0, A ⊗ K_g), A = Ainv⁻¹,where Phi (n × k) holds the per-record basis rows φ(s_r)ᵀ (see legendre_design), Z (n × q) is the record→animal incidence, Ainv (q × q) the relationship precision, K_g (k × k, positive definite) the supplied genetic covariance among the k random-regression coefficients, and finite sigma_e2 > 0. The genetic precision block is Ainv ⊗ inv(K_g) (animal-outer, coefficient-fastest ordering, matching the random design W = face-splitting(Z, Phi)).
Returns (beta, random_coefficients = (ids, values), variance_components = (K_g, sigma_e2), basis = (ncoef,)), where values is the q × k matrix of per-animal coefficient vectors (row = animal, column = Legendre coefficient 0..k-1).
EXPERIMENTAL, dense/validation-scale, SUPPLIED-covariance: K_g/sigma_e2 are NOT estimated here — use fit_random_regression_reml to estimate them by REML. There is no R-facing model-spec or bridge payload, and homogeneous residual variance only.
HSquared.fit_random_regression_reml Function
fit_random_regression_reml(y, X, Phi, Z, Ainv; initial = nothing,
iterations = 2000, ids = nothing)Estimate the random-regression coefficient genetic covariance K_g (k×k) and the homogeneous residual variance σ²e of the polynomial random-regression animal model by dense REML. Inputs match random_regression_mme (Phi the n×k basis design from legendre_design, Z the n×q record→animal incidence, Ainv the q×q relationship precision), but here K_g/σ²e are ESTIMATED rather than supplied. The marginal model is
y ~ N(X·β, V), V = W·(A ⊗ K_g)·Wᵀ + σ²e·I,with A = Ainv⁻¹ and W = face-splitting(Z, Phi) (animal-outer, coefficient-fastest — see random_regression_mme). The REML log-likelihood is maximized by Nelder–Mead over a log-Cholesky parameterization of K_g and log σ²e (so K_g stays positive definite and σ²e positive). At the optimum the coefficient BLUPs and β come from the marginal GLS BLUP form.
initial may supply K_g and/or finite positive sigma_e2; omitted fields use phenotypic-scale defaults. Returns a NamedTuple:
variance_components = (K_g, sigma_e2)— estimatedk×kcovariance + residual;beta— fixed effects at the estimate;random_coefficients = (ids, values)—q×kper-animal coefficient BLUPs (row = animal, column = Legendre coefficient0..k-1);loglik— full REML log-likelihood at the estimate (same scale assparse_reml_loglik; thek = 1reduction withK_g = [2σ²a]equals the univariate REML log-likelihood);converged,iterations,basis = (ncoef = k,).
EXPERIMENTAL, dense/validation-scale, REML-only, Gaussian, homogeneous residual variance, no permanent-environment term. Validated by deterministic self-consistency: the k = 1 reduction recovers the univariate fit_sparse_reml optimum (K_g[1,1] = 2σ²a, equal σ²e, equal log-likelihood); the reported log-likelihood matches an independent marginal oracle and beats off-optimum points; the BLUPs reproduce random_regression_mme at the estimate. Known-truth K_g recovery and any WOMBAT/ASReml/JWAS comparator are not yet exercised. The R twin has a covered opt-in k=2 rr(covariate, order) surface, separate from the Julia engine row and not the default route; broader R orders and model structures remain open. The covered aim + reporting convention lock (k = 2, basis, K_g reporting, h²-curve, frozen (x|g), k ≥ 3) is docs/design/22-rr-convention-lock.md. As a dense GLS path, V's conditioning degrades as O(1/σ²e) toward the residual boundary, so a near-noiseless optimum (σ²e → 0) is conditioning-limited at this validation scale.
HSquared.standardize_covariate Function
standardize_covariate(a; lower = minimum(a), upper = maximum(a))Affinely map a raw covariate vector a (e.g. age/time) onto t ∈ [-1, 1] by t = 2(a - lower)/(upper - lower) - 1 (so lower → -1, upper → +1, midpoint → 0). Returns a Vector{Float64}. Bounds and covariates must be finite; bounds must be distinct after conversion to Float64. Throws if the standardized values cannot be represented as finite Float64 values.
HSquared.rr_genetic_variance Function
rr_genetic_variance(K_g, ts)Per-covariate additive genetic variance trajectory v_g(t_i) = φ(t_i)ᵀ K_g φ(t_i) for a supplied k×k random-regression coefficient genetic covariance K_g and standardized covariate points ts ∈ [-1, 1] (k = size(K_g, 1)). Returns (covariate = ts, values). Descriptive only — K_g is supplied, not estimated.
HSquared.rr_genetic_covariance_surface Function
rr_genetic_covariance_surface(K_g, ts)Genetic covariance surface G(t_i, t_j) = φ(t_i)ᵀ K_g φ(t_j) across the standardized covariate points ts, i.e. Φ K_g Φᵀ (m×m, symmetric, PSD whenever K_g is PSD). Returns (covariate = ts, values). Its diagonal is rr_genetic_variance — exactly for a positive-definite K_g; for a reduced-rank (PSD) K_g the surface diagonal may carry a tiny negative roundoff at a near-zero-variance point, whereas rr_genetic_variance clamps such values to 0.
HSquared.rr_genetic_correlation_surface Function
rr_genetic_correlation_surface(K_g, ts)Genetic correlation surface across the standardized covariate points: the correlation matrix of rr_genetic_covariance_surface (D⁻¹ G D⁻¹, unit diagonal). Reuses genetic_correlation, so it throws if any covariate point has non-positive genetic variance. Returns (covariate = ts, values).
HSquared.rr_heritability Function
rr_heritability(K_g, residual, ts)Heritability trajectory h²(t_i) = v_g(t_i) / (v_g(t_i) + σ²_e(t_i)) for a supplied coefficient genetic covariance K_g and supplied residual variance residual — either a finite positive scalar (homoscedastic) or a length-m finite positive vector (heteroscedastic across the covariate). K_g and σ²_e are SUPPLIED, not estimated. Returns (covariate = ts, values). The genetic variance trajectory must be representable as finite Float64 values; an overflowing trajectory is rejected.
The supplied residual is treated as the TOTAL non-additive-genetic phenotypic variance at each covariate point. In the canonical repeated-records random-regression model (test-day, growth curves) a permanent-environment effect is essentially always present, so to avoid OVERSTATING h²(t) you must supply v_pe(t) + σ²_e(t) here, not the residual alone (the permanent-environment term is a later slice).
HSquared.rr_eigenfunctions Function
rr_eigenfunctions(K_g, ts)Eigen-function (covariance-function) decomposition of a supplied k×k random-regression coefficient genetic covariance K_g (Kirkpatrick, Lofsvold & Bulmer 1990). Eigen-decomposes K_g = Σ_j λ_j v_j v_jᵀ (descending λ_j ≥ 0, via genetic_pca) and evaluates the corresponding eigenFUNCTIONS ψ_j(t) = φ(t)ᵀ v_j — the orthonormal genetic principal curves of the reaction norm — at the standardized covariate points ts ∈ [-1, 1] (k = size(K_g, 1)).
Returns (covariate, eigenvalues, eigen_coefficients, eigenfunctions, variance_explained):
eigenvalues—λ_jdescending (the genetic variance carried by each eigenfunction);eigen_coefficients—k×k, columnjis the Legendre-coefficient eigenvectorv_j, sign-canonicalized as ingenetic_pca;eigenfunctions—length(ts)×k, columnjisψ_jevaluated atts(= Φ v_j,Φ = legendre_design(ts, k));variance_explained—λ_j / Σλ(zeros ifK_gis the zero matrix).
The eigenfunctions are orthonormal on [-1, 1] (∫ ψ_i ψ_j = δ_ij) and the covariance surface reconstructs spectrally as Φ K_g Φᵀ = Σ_j λ_j ψ_j ψ_jᵀ. Rotation-invariant and DESCRIPTIVE — K_g is SUPPLIED, not estimated (like genetic_pca / evolvability on a G). Under repeated eigenvalues the individual eigenfunctions are span-ambiguous (as in genetic_pca); the eigenvalues, variance explained, and the spectral reconstruction remain well-defined.
HSquared.rr_eigenfunctions_plot_data Function
rr_eigenfunctions_plot_data(K_g, ts)Plot-ready data for the Kirkpatrick eigenfunction figure: delegates to rr_eigenfunctions and returns (covariate, eigenvalues, eigenfunctions, variance_explained, basis_order, rotation_invariant = true, supplied = true). DESCRIPTIVE, supplied-K_g, rotation-invariant — each curve carries its genetic variance share λ_j/Σλ; eigenfunction signs are arbitrary (do not over-read).
HSquared.rr_genetic_variance_plot_data Function
rr_genetic_variance_plot_data(K_g, ts; residual = nothing)Plot-ready data for the genetic-variance trajectory v_g(t) (and the heritability trajectory h²(t) when residual is supplied): delegates to rr_genetic_variance / rr_heritability. Returns (covariate, genetic_variance, heritability, basis_order, supplied = true). DESCRIPTIVE from a SUPPLIED K_g/residual — not a REML estimate, not a phenotypic variance.
HSquared.rr_covariance_surface_plot_data Function
rr_covariance_surface_plot_data(K_g, ts; correlation = false)Plot-ready data for the genetic covariance (or correlation) surface heatmap: delegates to rr_genetic_covariance_surface (or rr_genetic_correlation_surface when correlation = true). Returns (covariate, surface, is_correlation, supplied = true). DESCRIPTIVE, supplied-K_g; genetic, not phenotypic.
HSquared.diagonal_covariance Function
diagonal_covariance(variances)Build a diagonal trait covariance matrix from positive variance components. This is the direct Julia engine utility for the Phase-4B diag genetic covariance structure. It is a matrix builder only; it does not change the R formula contract.
HSquared.lowrank_covariance Function
lowrank_covariance(loadings)Build a low-rank covariance matrix ΛΛ' from a traits × rank loading matrix Λ. The result is positive semidefinite and may be singular when rank < traits.
HSquared.factor_analytic_covariance Function
factor_analytic_covariance(loadings, uniqueness)Build a factor-analytic covariance matrix ΛΛ' + Ψ, where Λ is a traits × rank loading matrix and Ψ is a positive diagonal uniqueness vector. This is the direct Julia engine utility for the Phase-4B fa(K) covariance structure.
HSquared.fit_variance_components Function
fit_variance_components(spec; initial = (sigma_a2 = 1.0, sigma_e2 = 1.0),
method = spec.method, iterations = 1_000,
max_dense_cells = 1_000_000)Optimize the dense Gaussian ML/REML objective over positive variance components.
The optimizer works on log-variance parameters and uses Optim.NelderMead(). This is an experimental Phase 1 path for tiny validation examples. It is not AI-REML and is not the production sparse solver. max_dense_cells guards this dense path (hsquared#214): if nobs^2 + nanimals^2 exceeds it, this throws an ArgumentError naming the observed cell count and the effective cap. This is also the path the R default animal() route reaches, via fit_animal_model(...; kwargs...) forwarding max_dense_cells here.
HSquared.variance_components Function
variance_components(fit)Return the additive and residual variance components from an experimental low-level AnimalModelFit.
variance_components(result::NamedTuple)Return the genetic and residual covariance matrices from a multivariate HSquared result as (genetic_covariance, residual_covariance). When the result is a multi-trait repeatability fit, also return permanent_covariance.
HSquared.variance_component_covariance Function
variance_component_covariance(fit)Asymptotic covariance of the estimated (sigma_a2, sigma_e2) for a REML AnimalModelFit: the inverse of the average-information matrix. This is a large-sample approximation and is unreliable on small samples, where the REML surface is flat and the matrix is ill-conditioned. Experimental; REML only.
HSquared.variance_component_standard_errors Function
variance_component_standard_errors(fit)Asymptotic standard errors of (sigma_a2, sigma_e2) for a REML fit, as a NamedTuple. See variance_component_covariance for the caveats.
HSquared.fixed_effects Function
fixed_effects(fit)Return the fixed-effect estimates from an experimental low-level AnimalModelFit.
fixed_effects(result::NamedTuple)Return the fixed-effect matrix from a multivariate HSquared result.
HSquared.fit_diagnostics Function
fit_diagnostics(fit)Return compact status metadata for an experimental low-level fit result.
This is an extractor over fields already stored on the result object. It does not refit a model, run an optimizer, compute PEV/reliability, or change the bridge-facing result_payload contract.
HSquared.breeding_values Function
breeding_values(fit)Return animal-effect BLUPs/EBVs for an experimental low-level AnimalModelFit.
The current implementation solves Henderson's mixed-model equations at the fit's variance components and returns the animal-effect block. Variance component estimation is still the experimental dense path; this only changes the EBV/BLUP extraction equation solve.
sourcebreeding_values(result::NamedTuple)Return the multivariate breeding-value metadata from a multivariate HSquared result as (ids, traits, values).
breeding_values(fit::GeneticGLLVMFit)Return the q × T trait genetic conditional modes on the link scale from a GeneticGLLVMFit (internal struct). These combine common factors and, for factor-analytic fits, trait-specific genetic modes. They are invariant to an orthogonal rotation of the common factors; they are not posterior means.
HSquared.EBV Function
EBV(fit)Alias for breeding_values, matching the R twin's applied quantitative-genetic extractor vocabulary.
HSquared.BLUP Function
BLUP(fit)Alias for breeding_values. For the Phase 1 animal-effect block, the returned values are the same animal BLUPs/EBVs as breeding_values.
HSquared.fitted_values Function
fitted_values(fit; include_random = true)Return fitted values for an experimental low-level AnimalModelFit.
The current implementation solves Henderson's mixed-model equations at the fit's variance components, then computes X * beta + Z * u from that supplied variance solution. Variance-component estimation is still the experimental dense path.
HSquared.heritability Function
heritability(fit)Return narrow-sense heritability for the univariate Gaussian animal model: sigma_a2 / (sigma_a2 + sigma_e2). Fixed-effect variance is not included in the denominator. For fit types that estimate permanent-environment, common-environment, or maternal variance, those fitted components are included in the corresponding heritability denominator. Components must be finite and nonnegative, with at least one positive. The ratio is scaled before summation to avoid overflow in the total variance.
heritability(result::NamedTuple)Return the per-trait heritability vector from a multivariate REML result.
sourceHSquared.heritability_standard_error Function
heritability_standard_error(fit)Delta-method asymptotic standard error of h² = sigma_a2 / (sigma_a2 + sigma_e2) for a REML fit, from variance_component_covariance. Asymptotic; see the caveats there.
HSquared.heritability_interval Function
heritability_interval(fit; level = 0.95, method = :delta)Experimental two-sided confidence interval for h² of a REML AnimalModelFit.
method = :delta (default) builds the interval on the logit scale (delta method) and back-transforms, so it always lies in (0, 1); it returns (heritability, lower, upper, level, se, lower_clamped, upper_clamped, method).
method = :profile inverts the REML likelihood-ratio statistic: it profiles the REML log-likelihood over the total variance at each fixed h² and reports the h² range where 2·(ℓmax − ℓprofile(h²)) ≤ χ²₁,level. Endpoints that reach the (0, 1) search bounds are clamped. It returns (heritability, lower, upper, level, lower_clamped, upper_clamped, method) (no se). The *_clamped flags report endpoints that reached the numerical search rails rather than a likelihood-ratio crossing. Nonfinite profile targets or an invalid point-estimate anchor throw an ArgumentError; they are not reported as numeric endpoints.
Both are large-sample approximations: on small samples the REML surface is flat, so the intervals are wide.
sourceHSquared.repeatability_interval Function
repeatability_interval(y, X, Z, Ainv; level = 0.95, initial = ..., iterations = 200,
ids = nothing, fd_step = 1e-4, boundary_tol = 1e-6)Asymptotic delta-method confidence interval for the repeatability t = (σ²a + σ²pe) / (σ²a + σ²pe + σ²e) of the repeatability / permanent-environment animal model. Fits by REML (fit_repeatability_reml), forms the observed information as the central finite-difference Hessian of the REML log-likelihood at the optimum, and applies the delta method to t on the logit scale (so the interval lies in (0, 1)). Returns (repeatability, lower, upper, level, se, converged, boundary, lower_clamped, upper_clamped).
Experimental, asymptotic. t is the well-identified summary of this model (the σ²a/σ²pe split is weakly identified, so a per-component SE is unreliable, but t is stable). A non-converged fit, a component too close to zero for a finite-difference step that stays positive, or a repeatability on the (0, 1) rail returns NaN limits with boundary = true rather than a number from outside the parameter space.
HSquared.prediction_error_variance Function
prediction_error_variance(fit; method = :auto)Return prediction error variances for animal-effect BLUPs/EBVs from an experimental low-level AnimalModelFit.
method = :auto (the default since #350) picks :selinv when Ainv is sparse (a pedigree precision) and :dense when it is dense (a genomic Ginv), because the Takahashi recursion over a dense factor is slower than a dense inverse. :selinv reads the random-effect diagonal of the mixed-model-equation coefficient-matrix inverse through the Takahashi selected inverse of the sparse coefficient matrix, whose cost is Θ(Σⱼ|L[:,j]|²): it tracks the fill-in of the factor, not nnz(L), so :selinv is always memory-safe at pedigree scale but not uniformly faster than :dense at moderate size on a high-fill pedigree. :dense forms and inverts the dense coefficient matrix; it is the validation oracle and agrees with :selinv to machine precision.
prediction_error_variance(result::HendersonMMEResult; method = :auto)Return prediction error variances for a supplied-variance Henderson MME result.
Same method choices as prediction_error_variance(::AnimalModelFit): :auto (default: :selinv for a sparse Ainv, :dense for a dense one), :selinv (sparse Takahashi selected inverse) or :dense (the dense coefficient-matrix inverse, validation oracle).
HSquared.reliability Function
reliability(fit; method = :auto, pev = nothing)Return animal-level reliability values for the Phase 1 univariate animal model.
Reliability is computed as 1 - PEV_i / (sigma_a2 * A_ii) with A_ii the diagonal of the relationship matrix A = inv(Ainv) implied by the supplied precision. For a pedigree Ainv this is 1 + F_i (F_i the inbreeding coefficient); for a genomic spec (Ainv = Ginv) it is diag(inv(Ginv)) = diag(G) + ridge (the regularized genomic self-relationship, often ≠ 1), so the ridge perturbs the reported reliability/accuracy and the same extractor yields genomic reliabilities. method selects both the PEV path and the A_ii path: :auto (default since #350) resolves to :selinv for a sparse Ainv and :dense for a dense one; :selinv reads diag(inv(Ainv)) through the sparse Takahashi selected inverse of Ainv (no dense A is formed); :dense forms inv(Ainv) densely (validation oracle). Values are not clipped; small examples can expose weakly informed animals directly.
A supplied pev must contain finite nonnegative values and each model animal ID exactly once. Values and relationship diagonals are aligned by ID; the returned order follows pev.ids. Ratios are evaluated without requiring the intermediate product sigma_a2 * A_ii to fit in Float64; a nonfinite ratio is rejected.
HSquared.accuracy Function
accuracy(fit; method = :dense)Return animal-level accuracy values as sqrt(reliability(fit; method)).
This is a validation-scale extractor over the existing reliability method. It does not add independent accuracy validation and it rejects non-finite or out-of-range reliability values instead of silently clipping them. method is forwarded to reliability unchanged (:dense is the tiny-validation reference; :selinv avoids forming the dense n×n relationship matrix — see reliability). The default stays :dense for byte-identical behavior; callers at large pedigree scale should pass method = :selinv explicitly.
HSquared.result_payload Function
result_payload(fit)Return a bridge-facing result payload with field names aligned to the R hsquared_fit contract.
This is an experimental low-level payload. It is intended to make the R-Julia result shape explicit before live bridge execution is widened beyond tiny validation paths.
The payload includes prediction_error_variance and reliability as standard fields (each a (ids, values) named tuple). The PEV is computed through the sparse Takahashi selected inverse (method = :selinv, cost Θ(Σⱼ|L[:,j]|²), tracking fill-in; no dense matrix is formed), which matches the dense MME inverse diagonal to machine precision for well-conditioned validation-scale fits (V1-SELINV-PEV). The R twin unpacks these top-level fields directly via hs_julia_id_values() (hsquared#21), so the opportunistic per-extractor enrichment is no longer required. The Henderson MME is solved once here and shared by breeding_values and fitted_values (no second factorization). The PEV is computed once here and reused by reliability (no second factorization). Since #350 the reliability denominator also reads the animal self-relationships diag(inv(Ainv)) through the sparse selected inverse of Ainv, so for a sparse (pedigree) Ainv no dense matrix is formed anywhere in this payload; a dense genomic Ginv keeps its dense path (:auto). This remains a validation-scale claim, not a production large-pedigree reliability claim.
HSquared.hsquared Function
hsquared(args...; kwargs...)Planned high-level Julia entry point for inheritance-aware quantitative-genetic models.
This function is intentionally not implemented in Phase 0.
sourceHSquared.fit_animal_model Function
fit_animal_model(spec; target = :variance_components, ...)Fit or solve the Phase 1 Gaussian animal-model engine target for a validated AnimalModelSpec.
The default target = :variance_components dispatches to fit_variance_components, the experimental dense validation optimizer. target = :sparse_reml dispatches to fit_sparse_reml, the experimental sparse REML validation optimizer. target = :ai_reml dispatches to fit_ai_reml, the experimental sparse average-information REML validation optimizer. target = :henderson_mme requires supplied variance_components and returns a HendersonMMEResult. The Henderson target solves mixed-model equations at supplied variance components; it does not estimate them and does not return log-likelihood, AIC, df, or optimizer diagnostics.
fit_animal_model(...)Fit a Gaussian animal model through one of the supported Phase 1 engine methods. Implemented overloads accept an AnimalModelSpec or y, X, Z, Ainv. Targets include :variance_components, :sparse_reml, :ai_reml, and :henderson_mme; see the method-specific documentation for their arguments and limits.
Unsupported call shapes reach the fail-closed fallback and throw Phase0NotImplementedError.
HSquared.Phase0NotImplementedError Type
Phase0NotImplementedError(operation)Legacy error type for public placeholders and unwired payload routes.
sourceHSquared.FA_UNIQUENESS_FLOOR Constant
FA_UNIQUENESS_FLOORHeywood interior bound for fitted factor-analytic uniqueness (1e-4). Each fitted uniqueness is parameterized as ψ_i = FA_UNIQUENESS_FLOOR + exp(θ_i), so min(ψ̂) ≥ 1e-4 in the absolute trait-specific genetic-variance units of the fit. The bound is not scale-equivariant near the floor.
HSquared.bootstrap_variance_component_interval Function
bootstrap_variance_component_interval(fit::AnimalModelFit; level = 0.95,
n_boot = 1000, estimator = :sparse_reml,
rng = Random.MersenneTwister(0x48324352), max_dense_cells = 1000000,
min_converged = 2, min_converged_rate = 0.0)Parametric (Gaussian) bootstrap percentile confidence intervals for sigma_a2, sigma_e2, and h² = σ²a/(σ²a+σ²e) of a fitted univariate Gaussian REML animal model — a cross-check on the asymptotic delta / profile-LRT intervals (heritability_interval, variance_component_interval).
Mechanism: at the fitted (β, σ²a, σ²e), simulate Gaussian responses over the SUPPLIED relationship — a* = chol(inv(Ainv)).L · randn · √σ²a, e* = randn · √σ²e, y* = Xβ + Za* + e* — refit each replicate with the SAME REML estimator (:sparse_reml → fit_sparse_reml; :ai_reml → fit_ai_reml), and take percentile endpoints from the interior replicate vectors via the in-package type-7 _empirical_upper_quantile (no Statistics dependency). A replicate is dropped only for an expected numerical failure (the same predicate as _likelihood_proposal_failure), a non-converged refit, or a non-finite/boundary variance. Any other exception is rethrown. The returned tuple reports n_converged plus n_dropped_error, n_dropped_boundary, and n_dropped_nonconverged. The percentile interval is therefore conditional on interior refits; the call throws when n_converged < min_converged (default 2) or when n_converged / n_boot < min_converged_rate (default 0).
Returns a NamedTuple: sigma_a2, sigma_e2, heritability (the point estimates from fit); sigma_a2_ci, sigma_e2_ci, heritability_ci (each (lower, upper)); level, n_boot, n_converged, the three drop counts, method = :parametric_bootstrap_percentile; and replicates (the per-component interior-replicate vectors).
The interval FUNCTION is deterministic: rng defaults to a fixed-seed MersenneTwister, so the result is reproducible at the call site (only opt-in sim harnesses vary the seed).
EXPERIMENTAL, REML-only, univariate-Gaussian, dense/validation-scale (it forms inv(Ainv) + chol(A), guarded by max_dense_cells: if nobs^2 + nanimals^2 exceeds it, this throws an ArgumentError naming the observed cell count and the effective cap, hsquared#214, before the dense inverse is formed). It is a PERCENTILE bootstrap (BCa is out of scope); its OWN coverage is NOT calibrated in CI (an opt-in coverage sim is deferred follow-up) — it is the cross-check the delta/profile interval debt names, not evidence that those intervals are correct. Rejects non-REML fits; multivariate / non-Gaussian bootstrap CIs are separate slices.
HSquared.fa_covered_flip_cell Function
fa_covered_flip_cell(t, K) -> BoolReturn true only when ledermann_slack(t, K) > 0. Saturated cells may still be fitted for diagnosis; they are refused as a covered-flip cell.
HSquared.fit_matrix_free_reml Function
fit_matrix_free_reml(spec; nprobe = 64, tol = 1e-4, iterations = 200, seed = 0,
pcg_tol = 1e-9, pcg_maxiter = 2000, initial = nothing,
shared_probes = false, compute_loglik = true)MATRIX-FREE Monte-Carlo REML fixed-point fit of the SINGLE-effect Gaussian animal model (F6) — the AnimalModelFit-returning single-effect face of fit_multi_effect_mc_reml, for the one regime where the exact sparse fit_ai_reml is measured-infeasible and the eigen-once rescue is unavailable: high fill-in AND n past the dense eigen cap.
The mechanism it avoids is specific. fit_ai_reml's per-iteration cost is dominated NOT by the sparse Cholesky (cheap even at high fill — 0.35 s at q=20 000) but by the Takahashi selected inverse that supplies its exact score trace tr(A⁻¹C^uu), which scales with nnz(L) and is recomputed every iteration (381 s per call at q=20 000, fill 471 — so the whole fit walls at ~25 min; docs/dev-log/recovery-checkpoints/2026-07-24-f0-adversarial-highfill-decision.md). This fitter replaces that exact trace with the matrix-free Hutchinson estimator mc_reml_block_traces and the MME solve with matrix-free PCG, so C is NEVER assembled or factorized during the fit and the fill-in never materializes.
The estimator is therefore STOCHASTIC. A fixed probe seed is reused across iterations, so the iteration is reproducible for a given seed. This wrapper returns an AnimalModelFit and does not expose the low-level trace_mcse field. Use fit_multi_effect_mc_reml to inspect trace Monte-Carlo error; that error is not a gradient or variance-component standard error and excludes PCG solve error. A relative-change stop is labelled optimizer_status = "fixed_point_relchange", not "converged".
compute_loglik = true (default) evaluates the EXACT sparse_reml_loglik ONCE at the returned estimate. That costs two sparse Choleskys and NO selected inverse — cheap relative to the fit it just replaced — so the returned log-likelihood is exact, not stochastic. Set compute_loglik = false (loglik NaN) to skip the final exact-likelihood factorizations. The fit still validates each supplied relationship precision as positive definite, which requires sparse Cholesky factorization.
Returns the standard AnimalModelFit (target = :matrix_free_reml), so the ordinary extractors and the result payload work unchanged.
EXPERIMENTAL, REML-only, two-component, Gaussian-only. Ledger row V1-MATFREE-REML (validation_status(), status partial) — read it before quoting any number from this fitter.
The pre-declared tail-scale known-truth recovery gate PASSED on 2026-09-01 at q = 25,000 (48/48 converged, mean relative error against generating truth 0.0016/0.0007 against a bound of 0.05; docs/dev-log/recovery-checkpoints/2026-09-01-f6-matfree-tail-recovery-result.md). That is recovery-to-truth for the frozen gate fixture at that one scale — not at any other. The external same-estimand REML comparator (ASReml-R 4.2.0.482, docs/dev-log/recovery-checkpoints/2026-07-28-asreml-matfree-comparator.md) ran at q = 2000, below the measured crossover, so an at-scale comparator in the high-fill tail is still OWED before any covered claim.
OPT-IN ONLY. On this branch the fitter is reachable only by calling it directly: fit_animal_model accepts :variance_components, :sparse_reml, :ai_reml, and :henderson_mme, and there is no :auto REML router, so nothing can select this estimator on a user's behalf. That is deliberate — the regime a route would serve is largely the one in which this fitter has not been measured — and it is PINNED IN CI by the testset V1-MATFREE-REML opt-in fence, so re-wiring a route in fails loudly rather than silently. That fence covers this animal-model target router only; fit_multi_effect's :auto is a different estimator and does route to a matrix-free engine.
HSquared.fit_multivariate_repeatability_reml Function
fit_multivariate_repeatability_reml(Y, X, Z, Ainv; initial = nothing,
iterations = 2000, ids = nothing,
traits = nothing)Dense REML for the multi-trait repeatability (animal + permanent environment) model — the engine half of hsquared #237.
Y = X B + Z u_a + Z u_pe + E
vec(U_aᵀ) ~ N(0, A ⊗ G0), vec(U_peᵀ) ~ N(0, I ⊗ P0)G0, P0, and R0 are unstructured t×t covariances (honest for t = 2). Per-trait heritability uses the PE denominator h²_k = G0[k,k] / (G0[k,k] + P0[k,k] + R0[k,k]); repeatability is t_k = (G0[k,k] + P0[k,k]) / (G0[k,k] + P0[k,k] + R0[k,k]).
Reductions: t = 1 recovers fit_repeatability_reml variance components; P0 = 0 recovers the fit_multivariate_reml log-likelihood. Experimental, dense/validation-scale, REML-only. No covered flip.
The R twin has a dedicated experimental matrix route for cbind() plus permanent(). The generic scalar reader does not transport matrix effects. Direct Julia calls and payload-v2 dispatch :multivariate_repeatability use Y plus pedigree and iid permanent-environment blocks.
HSquared.fit_payload_v2 Function
fit_payload_v2(payload; scale_method = :dense, initial = nothing,
iterations = nothing) → fit NamedTupleParse a payload-v2 request and run the dispatched estimator, returning the estimator's raw NamedTuple result. For a single-pedigree-block payload this is exactly the result of fit_animal_model(…).
initial and iterations (hsquared#212) are forwarded ONLY to the :multi_effect and :direct_maternal dispatch arms — the two engine calls this dispatcher previously hardcoded with neither. Default nothing for both reproduces the exact pre-#212-fix call (each underlying fitter's own default: fit_multi_effect_reml's initial = nothing / iterations = 200, fit_direct_maternal_reml's initial = nothing / iterations = 200). For :multi_effect, forwarding applies on BOTH scale_method routes (#343): the default :dense fitter fit_multi_effect_reml, and the opt-in :auto route's fit_multi_effect, which itself forwards kwargs... to whichever engine it selects (fit_sparse_multi_effect_aireml's initial = nothing / iterations = 100, or fit_multi_effect_mc_reml's initial = nothing / iterations = 200). A supplied initial, iterations, or non-default scale_method on any other wired arm (:animal, :two_effect, :multivariate_repeatability; also non-default scale_method on :direct_maternal) raises ArgumentError and names the unsupported keyword (HSquared.jl#436). :multivariate and :coefcov still raise Phase0NotImplementedError.
The :coefcov dispatch is a frozen slot: fit_payload_v2 raises Phase0NotImplementedError for it (§6: "no coefcov payload fitting route is currently wired").
HSquared.genetic_correlation_interval Function
genetic_correlation_interval(fit, Y, X, Z, Ainv; level = 0.95, method = :delta,
pairs = nothing, fd_step = 1e-4)Experimental two-sided confidence interval for each off-diagonal genetic correlation r_g[i,j] of an :unstructured multivariate REML fit. Returns a NamedTuple of parallel per-pair vectors over the t(t-1)/2 strict-lower trait pairs — (trait_i, trait_j, estimate, lower, upper, lower_clamped, upper_clamped, level, method, converged) — plus the symmetric t×t matrices lower_matrix/upper_matrix (NaN on the diagonal and off the returned pairs; NO fabricated diagonal whiskers, mirroring the seRG diagonal-zero convention). pairs optionally restricts to a subset of (i, j) index tuples (i > j).
method = :delta (the only method this slice): reuses the validated observed- information + delta-method multivariate_covariance_standard_errors seRG and builds the interval on the Fisher-z scale (z = atanh(r_g), se_z = se_rg/(1-r_g²), endpoints tanh(z ∓ q·se_z)) so it always lands in (-1, 1); lower_clamped/upper_clamped are false (delta never clamps). If the SE path throws (non-PD observed information / boundary optimum, e.g. r_g → ±1), a clear ArgumentError is propagated rather than a fabricated whisker. (method = :profile, a profile-LRT inversion with the (i,j) correlation pinned, is explicit follow-up.)
Structured fits (:diagonal/:lowrank/:factor_analytic) are rejected because this interval uses only the unstructured-coordinate SE path. Off-diagonal genetic correlations are 0 by construction under :diagonal; FA covariance correlations can be identified at regular points but their interval path has not been implemented or calibrated here.
EXPERIMENTAL, asymptotic, REML-only, and NOT coverage-calibrated — a Wald approximation on the Fisher-z scale. It does NOT extend the V4-MV-REML covered claim (which is the unstructured REML POINT ESTIMATE only); this is a new uncertainty surface landing as partial-quality. An opt-in coverage harness is future work.
HSquared.ledermann_slack Function
ledermann_slack(t, K)Ledermann slack for a t-trait rank-K factor-analytic covariance: (t − K)² − (t + K). The nominal parameter count after removing factor rotations is m = tK + t − K(K−1)/2; slack is twice the difference between the ambient covariance dimension and m. Half the slack is an actual covariance codimension only where the covariance map has this nominal rank. Strictly positive slack is required before a cell can be a covered-flip candidate; this dimension-count screen does not prove local or global identification. Slack ≤ 0 is Ledermann-saturated (S1 t=3 K=1 is the disclosure cell).
HSquared.multivariate_repeatability_result_payload Function
multivariate_repeatability_result_payload(result)Bridge-ready payload for fit_multivariate_repeatability_reml. Exposes the three covariance blocks R needs (G0, P0, R0), component names/dimensions, PE-aware per-trait heritability and repeatability, breeding values, and permanent-environment BLUPs.
target = "multivariate_repeatability_reml". Experimental; no covered flip.
HSquared.nested_lrt Function
nested_lrt(loglik_constrained, loglik_full; df, boundary_df = 0, label = "LRT")Fit-agnostic likelihood-ratio test helper. The statistic is 2·(loglik_full − loglik_constrained); the reference distribution depends on boundary_df, the number of constrained parameters lying ON a boundary of the full parameter space under the null:
boundary_df = 0(interior null): plainχ²_dftail, valid under regular local identification and ordinary likelihood asymptotics.boundary_df = 1(one parameter on its boundary, e.g. a variance fixed at 0): under the regular single variance-component boundary conditions, the 50:50 chi-bar-squared mixture½·χ²_df + ½·χ²_{df−1}(Self & Liang 1987; Stram & Lee 1994). The naiveχ²_dftail is larger than this mixture tail. A boundary count alone does not establish these regularity conditions.boundary_df ≥ 2: weights are not computed. The naiveχ²_dftail is returned with the legacy tokenmixture = :chisq_conservative. This token does not certify the reference law or conservativity for arbitrary models; those require an appropriate regular cone limit and its component degrees.
Log-likelihoods must be finite before and after Float64 conversion. A positive overflow of the likelihood-difference statistic has the chi-square tail limit zero.
Returns (statistic, df, boundary_df, pvalue, boundary, mixture, note). Asymptotic theory only — the single-boundary mixture is textbook, not recovery-calibrated.
HSquared.nongaussian_heritability Function
nongaussian_heritability(fit::NonGaussianFit; mu = nothing, n_trials = nothing, predictor_variance = 0.0)
nongaussian_heritability(sigma_a2, mu, family::ResponseFamily; predictor_variance = 0.0)Latent- and observation-scale heritability for a non-Gaussian animal model — the Nakagawa–Schielzeth (2017) / de Villemereuil (QGglmm) transform that fills the gap the family-uniform nongaussian_result_payload deliberately leaves (it carries NO heritability, since "reuse the Gaussian ratio" is wrong off the identity link).
Returns a self-describing NamedTuple: (family, sigma_a2, mu, latent_total_variance, h2_latent, h2_observation, var_distribution, var_link, converged, information_limited, caveat, method).
Latent/link scale h2_latent = V_A / (V_A + V_link + V_fixed): V_link = π²/3 (logit), σ²e (Gaussian), and 0 for the Poisson log link — which makes the Poisson latent h² DEGENERATE, returned as NaN (the precise reason the payload refuses a single h²). Observation/data scale uses the QGglmm decomposition V_A,obs = Ψ²·V_A (Ψ = E[g⁻¹′(η)], the average inverse-link derivative; by Stein's lemma the exact variance of the regression of the mean on the breeding value, so the theoretical population projection satisfies h2_observation ∈ [0,1] when observation variance is finite and positive, including zero at V_A = 0), integrating over the LINEAR-PREDICTOR distribution η ~ N(μ, V_A + V_fixed) (the π²/3 logit residual is NOT added to the integration variance) via the module's existing 20-node Gauss–Hermite for logit, which approximates the population moments, and the log-normal closed form for Poisson. Estimand: PROPORTION for Bernoulli/Binomial, COUNT for Poisson; Gaussian reduces to V_A/(V_A+σ²e) on both scales. Gaussian uses this conditional ratio on both scales and requires predictor_variance = 0; nonzero fixed-effect spread is rejected. Threshold families (:bernoulli_probit, :ordered_probit) report the liability scale: the latent scale IS the liability with V_link = 1 (Dempster–Lerner 1950), the selection-relevant primary heritability V_A/(V_A+1+V_fixed) — returned in h2_latent (independent of μ and the cutpoints). The BINARY :bernoulli_probit observed-0/1 scale IS computed (the QGglmm probit integration Ψ²V_A/[p̄(1−p̄)] = the Dempster–Lerner transform, verified equal); the ORDINAL (K>2) per-category observed scale needs the cutpoints and stays a follow-up (h2_observation = NaN).
mu (link-scale population mean) defaults to the fit's single intercept; with >1 fixed effect it is REQUIRED (and predictor_variance, the fixed-effect linear- predictor variance, is the NS "variance explained by fixed effects" term). This is additional fixed-predictor spread, not total predictor variance. V_η = V_A + V_fixed assumes zero covariance between genetic and fixed contributions and a normal predictor approximation. The function REFUSES a non-converged fit; sets information_limited = true for the tested single-trial Bernoulli condition; available five-seed evidence does not establish a population bias direction or cause; and returns h2_observation = NaN with a caveat for genuinely varying per-record n_trials (a single data-scale h² is ill-defined under varying denominators). Constant fit trial vectors use their common scalar trial count.
EXPERIMENTAL, dense/validation-scale; exact in its closed-form limbs and anchored to an independent quadrature oracle in test/runtests.jl, but its estimates inherit uncertainty in latent σ²a (especially for single-trial Bernoulli); it has a QGglmm external comparator for the logit + binary-probit observation scales (comparator/qgglmm_probit_observed/) but not yet for the other observation scales / MCMCglmm, nor a Fisher/Falconer sign-off — not the public default, not covered. Deliberately NOT added to nongaussian_result_payload (that shape stays family-uniform).
HSquared.parse_payload_v2 Function
parse_payload_v2(payload) → ParsedPayloadV2Parse a payload-v2 request object (a Dict, NamedTuple, or any property- accessible container from JuliaCall) and return a ParsedPayloadV2 struct carrying resolved engine inputs and the chosen dispatch symbol.
Implements the §4 back-compat alias: a v0.1 payload (top-level Z/Ainv, no payload_version or payload_version = 1L) is lifted to a single pedigree block with relmat_status = "build_in_julia".
Implements the §4 transition alias: a payload with legacy Z/Z2/effect2 but no random_effects is lifted to a two-block list.
The §2 grammar table governs which block types are accepted; unknown types raise an ArgumentError.
Independent-block count (HSquared.jl #352): 1 pedigree block → :animal; 2 independent blocks → :two_effect (animal + permanent environment when the second block is iid); 3 or more independent blocks → :multi_effect. Two effects are not rejected.
CONTRACT-ONLY (docs/design/21-payload-v2-multiblock-schema.md §6): no new estimator is added here.
sourceHSquared.require_fa_covered_flip_cell Function
require_fa_covered_flip_cell(t, K)Throw if (t, K) is Ledermann-saturated (slack ≤ 0). Does not block fitting; it is the honest refuse path for a covered-flip cell. Returns the positive slack.
HSquared.result_payload_v2 Function
result_payload_v2(fit, parsed::ParsedPayloadV2) → NamedTupleBuild the block-structured result payload (§5) from an estimator fit and the ParsedPayloadV2 produced by parse_payload_v2.
Single-pedigree-block fast path (§5): when dispatch == :animal, the fitter must return AnimalModelFit; the wrapper delegates to result_payload to preserve the complete flat v0.1 fields and current R extractors. Partial raw tuples are rejected rather than exposed as incomplete legacy results.
For multi-block fits the result carries:
variance_components.blocks— ordered list of per-block variance records.variance_components.residual— scalar σ²e for univariate fits; the experimental multivariate-repeatability route carries a trait covariance matrix.random_effects— ordered list of(name, ids, values)records.loglik,df,nobs,diagnostics,converged— top-level fields.
The experimental multivariate-repeatability extension labels its matrix columns with top-level traits. It is Julia-only until the R normalizer and parity tests support that extension (docs/design/21-payload-v2-multiblock-schema.md §5).
CONTRACT-ONLY (docs/design/21-payload-v2-multiblock-schema.md §5, §6).
sourceHSquared.structured_genetic_payload Function
structured_genetic_payload(result)Rotation-INVARIANT Julia result payload for a :lowrank or :factor_analytic multivariate REML fit from fit_multivariate_reml — the rotation-gated companion of multivariate_result_payload (which serves the rotation-free :unstructured / :diagonal structures).
The factor loadings of a low-rank / factor-analytic G are rotation-NONidentified (Λ and ΛQ give the same G), so they are NEVER surfaced. Instead this exposes rotation-INVARIANT summaries of the estimated genetic covariance G (the FA rotation convention, docs/dev-log/decisions/2026-06-19-fa-rotation-convention.md): the reconstructed G, per-trait genetic variances / correlations, the genetic eigenstructure (genetic_pca: descending eigenvalues + sign-canonicalized principal axes; when eigenvalues are repeated, interpret the eigenspace rather than individual axes), mean_evolvability, and — for :factor_analytic — the separate specific variances Ψ (genetic_uniqueness, a rotation-invariant fitted component; rotation invariance does not establish that the decomposition identifies Ψ). The payload reports genetic_uniqueness_identification = :not_assessed_by_fit; callers must not infer identification from convergence or positive Ledermann slack. The genetic_rank field is the requested factor count K, not an estimated rank of G or automatic rank selection. Trait labels and all trait-indexed fields follow the input Y column order. heritability is Gᵢᵢ/(Gᵢᵢ + Rᵢᵢ) on the relationship matrix reference scale, and mean_evolvability is coordinate- and unit-dependent. Fixed effects, breeding values, per-trait heritabilities, the REML loglik, and converged are rotation-invariant and carried through. It deliberately OMITS genetic_loadings; rotation_invariant = true and loadings_excluded = true make that self-describing.
HSquared.variance_component_interval Function
variance_component_interval(fit; level = 0.95, method = :profile)Profile likelihood-ratio confidence interval for the additive variance component sigma_a2 of a REML AnimalModelFit. It inverts 2·(ℓmax − ℓprofile(sigma_a2)) ≤ χ²₁,level while profiling the residual variance sigma_e2 as a NUISANCE at each candidate sigma_a2 — the variance-component companion of heritability_interval method = :profile (which profiles the total variance at fixed h²).
Returns (sigma_a2, lower, upper, level, lower_clamped, upper_clamped, method). The *_clamped flags report an endpoint that reached the (sigma_a2·1e-4, sigma_a2·1e4) search bound (the profile did not cross the χ² threshold within range), so a non-crossing endpoint is self-describing — on small samples the REML surface is flat and the interval clamps.
Applies unchanged to a genomic GBLUP / supplied-Ginv REML fit (fit_gblup_reml, or any method = :REML spec with Ginv in the Ainv slot): the profiler reads only the spec's precision through sparse_reml_loglik. On a genomic spec sigma_a2 is the GENOMIC additive variance and the interval is CONDITIONAL on the supplied Ginv (ridge + centering), so the implied A_ii = diag(inv(Ginv)) ≠ 1 — not a pedigree-scale sigma_a2. If a failed profile evaluation, including a non-PD precision, yields a nonfinite likelihood-ratio target, the interval throws an ArgumentError; this is not a clamped confidence endpoint.
Experimental, asymptotic, REML only; no coverage calibration.
sourceHSquared.gpu_fit_gblup Function
gpu_fit_gblup(y, X, Z, markers, sigma_a2, sigma_e2; ridge = 0.01,
method = :vanraden1, allele_frequencies = nothing, weights = nothing)DEVICE-RESIDENT genomic BLUP (GBLUP) at supplied variance components — the GPU twin of the CPU pipeline genomic_relationship_matrix → genomic_relationship_inverse → fit_gblup, run end-to-end on the device (v0.7 slice G-A). STUB: the method is provided by the HSquaredCUDAExt package extension (loads on using CUDA); without a CUDA backend in scope this throws a MethodError.
The genomic relationship G = W·Wᵀ/k (GEMM), its ridge-regularized inverse Ginv = inv(G + ridge·I) (Cholesky + solve), and the dense Henderson mixed-model equations C·[β; u] = rhs (C = [[X'X/σe² X'Z/σe²]; [Z'X/σe² Z'Z/σe² + Ginv/σa²]]) are built and solved on the device, keeping the q×q matrices G/Ginv/C ON-DEVICE across all three stages — only the marker matrix goes up and the fixed effects β + genomic breeding values u come back. This is a NUMERICAL ACCELERATION of the CPU GBLUP, NOT a new estimand: the marker centering/validation is the validated CPU centered_markers verbatim (same p, k, guards), and the returned (beta, breeding_values) equal (to floating-point tolerance) the CPU fit_gblup(y, X, Z, genomic_relationship_inverse(genomic_relationship_matrix(markers; method), ridge), sigma_a2, sigma_e2). sigma_a2/sigma_e2 are SUPPLIED, not estimated. Float64 throughout.
Returns a NamedTuple (beta, breeding_values = (ids, values), converged). The CPU↔GPU agreement test and the end-to-end benchmark are opt-in cluster scripts (sim/drac/g_a_device_gblup.jl), run on a CUDA device — never in CI. No performance or agreement claim holds until a committed run lands its artifact.