Changelog

Unreleased

• Authorship: added co-authors Yefeng Yang (ORCID 0000-0002-8610-4016) and Szymon Drobniak (ORCID 0000-0001-8101-6247) in Project.toml, CITATION.cff, and Documenter authors; Shinichi Nakagawa remains primary contact. No version or capability change.

• Identifiability warning #416: the fit-free rank pre-check that fit_two_effect_reml and fit_multi_effect_reml already ran now also runs in fit_sparse_multi_effect_aireml, fit_multi_effect_mc_reml, fit_repeatability_reml (one record per animal makes sigma_pe2 and sigma_e2 inseparable) and fit_direct_maternal_reml. The numerical rank check warns of possible weak identification. Nonzero covariance columns are now normalized before SVD, so the unchanged 1e-8 relative tolerance applies to their directions rather than their raw scales; this changes numerical warning decisions in the existing dense routes as well. Rescaling invariance assumes the computed columns remain finite and nonzero; overflow or underflow during covariance construction can lose information, and nonfinite columns skip the check. No warning is not evidence of identification or practical precision. A warning does not prove structural non-identifiability or identify which components are confounded. The check only logs: estimates, converged and result fields are unchanged. The sparse and matrix-free routes skip it above 1e7 stacked cells.

• Genetic GLLVM (#50, owner decision 25 Oct 2026): fit_gllvm_laplace_reml, the GeneticGLLVMFit fit object and a new genetic_gllvm_result_payload(fit) are now exported so the R lane can bridge them. Both are labelled partial and experimental (the payload carries status = "partial", experimental = true, scale = "latent"). The payload reports only rotation-invariant functionals of G_lat (never raw loadings) and no heritability. GeneticGLLVMFit gains a trailing families field (per-trait response family tokens). No capability status, row count or version changes.

• Changed #413: the dense fit_two_effect_reml and fit_repeatability_reml fitters (and two_effect_ratio_interval, repeatability_interval, which call them) default to scale-aware starting variances set from the response variance: sigma_e2 = var(y)/2 and the two random-effect variances var(y)/4 each, the same rule as initial = :auto on the sparse multi-effect fitter. The default start was (1, 1, 1), which is blind to the unit of y: on a response with variance near 1e7 (milk yield in kg x 1000) NelderMead stopped at a wrong near-zero-heritability point with converged = true, while y / 1000 gave the correct REML answer. The optimum now no longer depends on the unit of the response. A constant response falls back to (1, 1, 1). An explicit initial is used as given, and passing (1, 1, 1) explicitly with var(y) far from 1 still emits the unit-scale start warning (#447). The default start changes the NelderMead path, so estimates on unit-variance data can differ from earlier releases in the last digits and the reported optimum can differ where the old start stopped early; the sim/ recovery harnesses pass the unit start explicitly and are unchanged. converged still means only that the simplex tolerance was met, not that the REML score is small (superseded by #424 below).

• Changed #424: the eleven NelderMead fitters (fit_variance_components, fit_sparse_reml, fit_two_effect_reml, fit_multi_effect_reml, fit_direct_maternal_reml, fit_repeatability_reml, fit_random_regression_reml, fit_multivariate_reml, fit_multivariate_repeatability_reml, fit_gllvm_laplace_reml, and the NelderMead branches of fit_laplace_reml) no longer take converged from Optim.converged alone. One shared helper in src/likelihood.jl also requires, at the returned optimum, a numerically zero gradient (Newton decrement at most 1e-3 log-likelihood units) and curvature in every direction (smallest eigenvalue of the finite-difference Hessian in correlation form above 1e-3). Coordinates already at a boundary (a variance share at or below 1e-6, a covariance that lost rank or a variance in it, an FA uniqueness at its floor, a #327 safety rail) are left out of the check, and the rotation directions of free loadings are set aside. A flat REML ridge, such as Va and Ve with one record per unrelated animal, or a rank-1 factor-analytic G_lat for two traits in fit_gllvm_laplace_reml, now reportsconverged = false, and so does a simplex that stopped short of the optimum. The check costs2d^2 + 4dextra objective evaluations ford checked parameters. Point estimates of a single fit are unchanged; the multivariate FA two-start search now prefers a start that passes the check, and bootstrap refits that fail it count as not converged.

None of the entries below change any capability status, row count, or package version. public_covered_count stays 7, and no entry below touches Project.toml. The package version is 0.9.0; its released section is recorded below.

• Part of #408: new variance_component_uncertainty(fit; level = 0.95) returns the variance-component covariance, its standard errors, the h² standard error and the logit-delta h² interval of a converged REML AnimalModelFit from one average-information matrix. Every field is bit-identical to the standalone call, which each build and factorise their own matrix, so the bundle saves two of three factorisations. The R bridge does not call it yet. No estimand, interval or coverage claim changes.

• #411 and #409 (no result change; outputs are bitwise identical to before on the test fixtures). The refitting forms of repeatability_interval, two_effect_ratio_interval, direct_maternal_interval and multi_effect_ratio_interval now reuse the dense A = inv(Ainv) (and dense X, Z, y) that their fit already built, instead of inverting each Ainv a second time for the finite-difference information. The profile intervals (heritability_interval(fit; method = :profile) and variance_component_interval) keep one mixed-model-equation workspace per interval: the cross products and the sparsity pattern are formed once, each evaluation refills the coefficient matrix in place and refactors it with the cached symbolic analysis (cholesky!), where before every evaluation re-formed the cross products, concatenated the blocks and reordered. The number of evaluations is unchanged.

• Part of #411: repeatability_interval, two_effect_ratio_interval and direct_maternal_interval each gain a method that takes an already-fitted result (repeatability_interval(fit, y, X, Z, Ainv; ...), and likewise for the other two) and computes only the finite-difference information at fit.variance_components, so a caller that has just fitted no longer repeats the whole dense Nelder-Mead fit inside the interval. The existing refitting methods now call the same code and return identical values. repeatability_interval also accepts max_dense_cells and forwards it to its fit; before, a cap raised for fit_repeatability_reml was ignored by the interval, which refused any model over the default 1,000,000 cells. Still open: the R bridge calls these functions with the refitting form, the sparse multi_effect_* route is not yet used for the repeatability and direct-maternal bridge targets, and each interval still forms one dense inv(Ainv).

• reliability's denominator diag(inv(Ainv)) comes from 1 .+ F when Julia builds Ainv from pedigree rows: pedigree_inverse already computes the inbreeding coefficients F (Meuwissen & Luo) for Henderson's rules, so the spec now carries them (animal_model_spec(...; relationship_diag)) and method = :auto reads them instead of a selected inverse of Ainv. The payload-v2 route carries this diagonal when it builds a pedigree block in Julia; callers that pass only Ainv, including the current R formula bridge, still obtain it through the selected inverse. Values are unchanged to 1e-10; explicit :selinv and :dense keep their own paths. Szymon Drobniak.

• Fixed #350: prediction_error_variance and reliability default to method = :auto, which takes the sparse Takahashi selected-inverse path for a sparse Ainv and the dense path for a dense genomic Ginv (:selinv and :dense stay explicit; the recursion over a dense factor is 67x slower than inv), and reliability reads diag(inv(Ainv)) through the selected inverse of Ainv instead of forming the dense inverse, so result_payload and breeding_values_plot_data no longer form any dense inverse for a sparse Ainv. The original q = 5,000 timing was removed because its pedigree generator, call script, and hardware were not recorded; it was not a reproducible package benchmark. Numerics are unchanged to machine precision. Szymon Drobniak (finding). The selected-inverse path refuses a numerically singular coefficient matrix (e.g. a duplicated column of X) with a named ArgumentError, using a relative-pivot test (L_ii² / C_ii) that is invariant to the units of a covariate: a well-posed fit with a covariate stored at magnitude 10⁶ or more is not refused, and neither is a degenerate-but-well-posed boundary fit (the guard's floor is eps-level, 1e-15; a duplicated column of X measures 7.4e-18).

• Fixed #370/#371: K-effect post-fit uncertainty can call multi_effect_uncertainty once and reuse one covariance calculation for component standard errors, ratio standard errors, and sum-ratio intervals. This applies only to the multi-effect K-effect route; the univariate animal-model route is unchanged. Performance figures from an undistributed great-tit dataset are intentionally omitted here because readers cannot reproduce them from this repository.

• Fixed #334: tools/write_validation_status_page.jl no longer stamps a <!-- regenerated: <timestamp> --> comment into docs/src/validation-status.md. The stamp came from now(UTC) on every call and was rewritten unconditionally, so a plain docs/make.jl build dirtied the tracked page even when the generated table body was unchanged. The BEGIN/END generated-block markers are unchanged. Docs-build hygiene only; no table content changed.

• Clarified #333: fit_snp_blup's docstring no longer states the GBLUP↔SNP-BLUP equivalence unconditionally. The equivalence is exact for the unregularized G; 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. Docstring-only; no arithmetic changed; the per-cell measurement is in issue #333; docs/src/genomic-models.md gives the reason, not the numbers.

• Fixed hsquared#210 (Julia half): the Willham total-heritability docstrings and the stored convention string in direct_maternal_interval wrote σ_P for a variance. Changed every such user-visible string in src/ to σ²_P, consistently, with no change to any arithmetic (src/likelihood.jl, direct_heritability/maternal_ratio/ total_heritability docstrings and the returned convention field).

• Fixed hsquared#212 (Julia half): fit_payload_v2 accepts and forwards initial/iterations on the :multi_effect (dense) and :direct_maternal arms, which previously hardcoded the engine call and discarded both; fit_gblup_reml, fit_single_step_reml and fit_metafounder_single_step_reml gain an iterations keyword. Defaults unchanged (nothing omits the keyword entirely). The R half of #212 is separate and does not land here.

• Fixed #331: covariance_structure_lrt's df counts identified parameters — _mv_nparams removes the r(r−1)/2 rotational indeterminacy of Λ for :lowrank and :factor_analytic, so a Ledermann-valid comparison is no longer refused with df = 0. The function now always uses the plain χ²df reference and never the 50:50 chi-bar mixture: a :factor_analytic null is a regular submanifold and reports boundary = false, while a :lowrank null stays boundary = true with the naive tail and its direction explicitly unknown — the word "conservative" is gone. A new reference field names which distribution produced pvalue. Behaviour change: a structured comparison with df = 1 (reachable at :lowrank, t = 2, rank = 1) now reports twice its former p-value.

• Fixed hsquared#214/#217 (Julia half): max_dense_cells is a keyword on the dense-validation fitters, and fit_repeatability_reml now runs the dense-cell guard.

• Fixed #327: NonGaussianFit gains a boundary field; nongaussian_three_field_payload refuses boundary-riding fits; an opt-in restart_check (with restart_estimate) is added; the :nbinom and :gaussian branches gain the same log-unit rail the gamma branch already has. No capability-status or count changes. Defaults are unchanged for the seven families that already had a bound; for :gaussian and :nbinom the new rail applies on the default path, so a fit whose optimum lies outside exp(log(initial) ± 8) now stops on the rail and sets boundary = true where it previously searched unbounded.

• Fixed #343 (follow-on to hsquared#212/#337): fit_payload_v2's :multi_effect dispatch now forwards initial/iterations on the opt-in scale_method = :auto route too (previously forwarded on :dense only), via fit_multi_effect's own kwargs... to whichever engine :auto selects. Default calls on either route are unaffected (nothing omits the keyword entirely).

• Clarified #347 (follow-on to #327): nongaussian_three_field_payload's boundary-refusal ArgumentError, and the corresponding fit_laplace_reml docstring paragraph, no longer advise restart_check = true as a way to clear a flagged boundary — that check is one-directional and can only make boundary stricter, never clear it. Both now name the actual lever: a different initial, which recentres the log-scale search bracket/rail. Inside the payload builder the family gate now runs before the boundary gate, so the message's ± 6 bracket is exact for every family that can reach it (the jointly-estimated families stop on a ±8-log-unit rail instead); a boundary-flagged fit of an unsupported family is now refused by the family message rather than the boundary message. Wording plus that one refusal-order change; no numerical behaviour change.

0.9.0 (experimental)

This release retained the experimental package status and did not authorize a new covered-capability count. It consolidated the Julia engine and R-bridge contracts present at the 0.9.0 tag; subsequent development remains under Unreleased above.

0.8.0 (experimental)

Experimental numbered bump after both 0.8 engine pillars are covered (FA V4-FA + single-step V2-SSHINV). Version tracks the engine pillar pair, not R-public count. public_covered_count stays 7. Experimental label retained. Not production, not Julia General. R FA stays planned. R single_step() stays opt-in partial.

0.7.0 (experimental)

Experimental numbered bump after R G10. Version tracks R-public covered capability (public_covered_count 6 → 7). Engine V4-MV-REML was already covered; this is not an engine covered flip. Experimental label retained. Not production, not Julia General. k≥3 / diagonal multivariate and interval calibration stay out of the covered claim.

0.6.0 (experimental; historical)

This historical experimental marker followed the R-public unstructured multivariate route, which changed public_covered_count from 5 to 6. V4-MV-REML was already engine-covered, so this was not an engine covered-row promotion. It did not establish production readiness, a Julia General registration, interval calibration, or a 0.9 release. The dated coordination record is 2026-09-02.

0.5.0 (experimental)

Experimental numbering change (0.0.1 → 0.5.0). The engine remains experimental. This is not a covered flip, not a production-ready release, and not a Julia General registration. public_covered_count stays 5. Install from the GitHub URL until a General PR merges.

• Corrected the Julia #44/R-bridge status wording after hsquared PR #96. The R twin already has an opt-in target = "nongaussian" bridge for Poisson and Binomial LA/VA fits plus the PR #95 normalizer fixture; the remaining bridge gap is per-record varying-trial formula/bridge activation and deeper validation/comparator/calibration evidence, not wholesale absence of R non-Gaussian bridge fitting.

• Recorded hsquared PR #95 consumption of the test/fixtures/non_gaussian_parity/ target in Julia status ledgers. This banks R Julia-free NonGaussianFit normalizer parity only; per-record varying-trial activation, interval calibration, external comparator evidence, and covered-status promotion remain open.

• Added test/fixtures/non_gaussian_parity/, a Julia-native #44 bridge payload target for nongaussian_result_payload(::NonGaussianFit). The fixture pins Poisson Laplace and per-record Binomial variational payloads for R normalizer tests; it is not per-record varying-trial R activation, interval calibration, or external comparator evidence.

• Added pedigree normalization and sparse Ainv construction utilities.

• Added low-level animal-model specification validation.

• Added dense Gaussian ML/REML log-likelihood evaluation at supplied variance components.

• Added sparse REML log-likelihood evaluation at supplied variance components using the Henderson MME determinant identity.

• Added experimental fit_sparse_reml() and fit_animal_model(...; target = :sparse_reml) dispatch for REML-only sparse validation optimization.

• Added an experimental sparse prediction-error-variance / reliability path: prediction_error_variance(...; method = :selinv) and reliability(...; method = :selinv) use a Takahashi selected inverse of the sparse Henderson MME coefficient matrix (kernel adapted from DRM.jl, MIT). The default stays :dense and matches it to machine precision on tiny fixtures.

• Added an experimental average-information REML estimator: fit_ai_reml() and fit_animal_model(...; target = :ai_reml) estimate the two variance components by AI/Newton steps on the sparse Henderson MME (score from the selected inverse, AI matrix from working-variate re-solves). Validated to recover the same optimum as the NelderMead optimizers; REML-only, Gaussian, experimental.

• Added genomic_relationship_matrix() — the VanRaden (2008) genomic relationship matrix G from a 0/1/2 (or dosage) marker matrix. Construction utility only (Phase 2 start); no genomic fitting yet.

• Added genomic_relationship_inverse() — the ridge-regularized dense inverse inv(G + ridge·I) of a genomic relationship matrix, intended for later GBLUP use. Construction utility only; not wired into model fitting, and no single-step (H-matrix) blending.

• Added fit_gblup() — genomic BLUP at supplied variance components: feeds a genomic Ginv into the existing Henderson MME (Ginv in the Ainv slot). Reproduces an independent dense MME and reproduces pedigree BLUP when G = A; experimental, supplied-variance only, no variance-component estimation, no single-step, no external comparator.

• Added fit_snp_blup() and centered_markers() — SNP-BLUP / RR-BLUP marker effects via the existing Henderson MME (Z = W centered markers, identity prior, σ²_marker = σ²_g/k). gebv = W·â equals GBLUP to machine precision (the GBLUP↔SNP-BLUP equivalence, validated via the marginal V); experimental, supplied-variance only, unweighted VanRaden method-1.

• Added test/fixtures/genomic_gblup_snpblup_target/ — a Julia-native #49 comparator target for supplied-variance GBLUP/SNP-BLUP with a positive-definite supplied-frequency VanRaden G, Ginv, beta, GEBVs, marker effects, metadata, and a generator. This pins a future external comparator target; it is not external comparator evidence.

• Added single_marker_scan() — a Phase-5 fixed-effect single-marker screening utility that residualizes y and centered marker dosages against X, then reports marker effects, supplied-variance standard errors, Wald z-scores, and chi-square statistics plus approximate two-sided Gaussian/Wald p-values, Bonferroni-adjusted p-values, Benjamini-Hochberg q-values, and fixed-effect known-variance LOD-equivalent scores. marker_manhattan_data() prepares plot-ready marker IDs, chromosome/position metadata, cumulative plot positions, raw p-values, and -log10(p) values, including overloads for already-validated HSMarkerMapSpec / HSData marker metadata. marker_region_data() prepares one-chromosome or coordinate-window data slices from the same row-aligned scan fields, preserving scan indices and optional marker-variance proportions for future regional plot/fine-mapping front ends without activating them. marker_significance_summary() reports nominal returned-marker-set raw, Bonferroni, and BH significance flags/counts plus top-marker provenance from the same scan fields; it is not a calibrated genome-wide threshold workflow. marker_effects() prepares deterministic sorted marker-effect summaries from the same scan fields, with optional metadata alignment. marker_scan_table() prepares row-aligned scan tables from the same direct scan fields in original scan order, with allele variances, marker-variance contributions, optional total-variance proportions, optional mixed/LOCO fields when present, and optional marker-map metadata alignment. gwas_table(), qtl_table(), and eqtl_table() are semantic wrappers over those already-computed direct scan tables; they add an analysis label plus optional trait or expression-feature metadata, but do not run GWAS, interval mapping, or expression-wide eQTL scans. marker_variance_explained() prepares deterministic marker-level variance-contribution summaries as 2p(1-p) * effect^2, with optional total-variance proportions and metadata alignment. marker_qq_data() prepares sorted observed/expected QQ plot data from the same direct scan output. marker_genomic_inflation() computes a genomic-control-style lambda_GC diagnostic from returned chi-square statistics. This is direct Julia engine tooling only: no mixed-model GWAS/QTL, no interval-mapping or mixed-model LOD workflow, no marker-file parser, no plotting backend, no regional_plot()or fine-mapping activation, no p-value calibration or calibrated/correlated-marker genome-wide threshold workflow, no calibrated PVE/model R² claim, no Rmarker_scan() syntax, no bridge payload change, and no comparator parity.

• Added opt-in sim/phase5_marker_scan_recovery.jl for seeded marker-signal recovery of the direct fixed, supplied-variance mixed, and supplied LOCO marker-scan helpers outside CI. Default seed 20260614 passes all three cases on a half-sib simulated design. This is internal recovery smoke evidence, not broad multiple-testing calibration, QTL/eQTL validation, public R syntax, bridge payload change, or comparator parity.

• Added mixed_model_marker_scan() — a dense, supplied-variance GLS marker-screening helper that forms V = sigma_a2 * Z * A * Z' + sigma_e2 * I from supplied variance components and a supplied relationship precision, then runs marker-by-marker Wald tests conditional on X. Tests pin reduction to single_marker_scan() when the random-effect design contributes zero covariance and agreement with an independent GLS calculation. Direct Julia engine tooling only: no marker-scan variance-component estimation, no LOCO, no sparse production scan, no calibrated p-values, no calibrated PVE/model R² claim, no plotting backend, no R marker_scan() syntax, no bridge payload change, and no comparator parity. The shared marker_effects(), marker_variance_explained(), Manhattan, QQ, and inflation helpers work over the returned scan fields.

• Added loco_relationship_precisions() and loco_mixed_model_marker_scan() — dense validation-scale leave-one-group-out marker-screening helpers. The construction helper drops each marker group, builds a VanRaden relationship matrix from the remaining markers, and applies the existing ridge-regularized inverse. The scan helper selects the matching precision before running the dense supplied-variance GLS scan. Tests pin explicit LOCO precision identities and marker-wise agreement with separate mixed_model_marker_scan() calls. Direct Julia engine tooling only: no public LOCO defaults, no marker-scan variance-component estimation, no sparse production scan, no calibrated p-values, no calibrated PVE/model R² claim, no plotting backend, no R marker_scan() syntax, no bridge payload change, and no comparator parity. The same effect-summary, marker-variance, and diagnostic helpers remain direct Julia-only.

• Internal: deduped the numerator-relationship recursion into _numerator_relationship(pedigree[, rows]) (one source shared by inbreeding_coefficients and single-step A₂₂ construction); removed the test-only duplicate. No public API or behavior change.

• Documented and tested genomic reliability/PEV/accuracy: for a GBLUP fit the existing extractors use the regularized genomic self-relationship diag(inv(Ginv)) = diag(G) + ridge as the denominator (the ridge perturbs the reported reliability/accuracy), and method = :selinv PEV matches the dense diagonal. No logic change.

• Added an internal single-step H-inverse construction utility _single_step_Hinv — H⁻¹ = A⁻¹ + scatter(τG⁻¹ − ωA₂₂⁻¹) on the genotyped rows (with the subtle A₂₂⁻¹ = inv(A[g,g]), not (A⁻¹)[g,g]). Property-checked (reduction, locality, symmetry, distinctness, scattered rows, singular-G guard); unexported, not wired into fitting, blending/τ/ω/ridge defaults not comparator-validated.

• Validated genomic REML variance-component estimation: the existing REML optimizers (fit_ai_reml, fit_sparse_reml) estimate σ²g/σ²e on a genomic Ginv spec — AI and NelderMead reach the same optimum and fit_gblup at the estimate reproduces the REML breeding values. Experimental; no external comparator; no new code (reuses the Phase-1 optimizers).

• Added a "Genomic models" documentation page with runnable examples for G/Ginv, GBLUP, genomic REML, SNP-BLUP, and the single-step H⁻¹ utility, and an explicit experimental / not-yet-R-wired / no-external-comparator boundary.

• Added experimental heritability uncertainty: variance_component_covariance, variance_component_standard_errors, heritability_standard_error, and heritability_interval — a logit-transform delta interval (always in (0,1)) built on the REML AI matrix, with a self-contained standard-normal quantile. Asymptotic / REML-only; wide and unreliable at small n.

• Added repeatability_mme() — the first Phase-3 slice: a supplied-variance Henderson solve of the two-random-effect repeatability / permanent-environment animal model (additive a + permanent-environment pe) via a block-diagonal relationship precision. Matches an independent marginal-GLS BLUP and reduces to the animal model as σ²pe → 0. Experimental, supplied-variance only — no REML estimation of the three components, no R-facing model-spec yet.

• Added fit_repeatability_reml() — REML estimation of (σ²a, σ²pe, σ²e) for the repeatability / permanent-environment model (dense two-random-effect REML loglik + NelderMead), returning the repeatability coefficient t and h². Deterministic checks (loglik reduces to the animal-model REML at σ²pe=0; BLUPs match repeatability_mme; optimum beats a grid) + a one-off seeded recovery. Experimental, dense/validation-scale, REML-only; no committed recovery test (suite kept RNG-free), no intervals, no R-facing model-spec.

• Added two_effect_mme() — the general supplied-variance kernel for two independent random effects (each with its own incidence, relationship inverse, and variance), covering common-environment and maternal-environment models; repeatability_mme is now its Z2=Z1, A2=I special case. Validated against an independent marginal-GLS BLUP. Experimental; no correlated direct–maternal genetic, no estimation, no R-facing model-spec.

• Added fit_two_effect_reml() — REML estimation of the two-effect-model variances (common-environment c², maternal variance, etc.) via the dense two-effect REML loglik; fit_repeatability_reml is the Z2=Z1, A2=I reduction. Experimental, dense/validation-scale, REML-only; no committed recovery test, no intervals, no R-facing model-spec.

• Added a "Standard QG models" documentation page with runnable examples for repeatability (MME + REML, repeatability coefficient), common-environment, and maternal-environment models, with the experimental / not-yet-R-wired boundary.

• Validation status (v0.1 gate): corrected the V1-AI-REML evidence to cite the committed finite-difference REML Hessian check (V1-HERIT-CI) instead of an uncommitted "250-animal observed-information" claim, and recorded the R-lane external validation — V1-MRODE-FIT and V1-COMPARATORS move to covered_external (the engine recovers the published gryphon REML estimate, Wilson 2010 h²=0.470, via supplied A_gryphon, and agrees with sommer, within the maintainer-signed-off band). Status rows + tests only; no engine change.

• Reconciled ROADMAP.md with merged reality: Phases 1-3 engine utilities (genomic prediction, single-step H-inverse, AI-REML, the Takahashi selected inverse, repeatability / two-effect models) are now listed as landed experimental, validation-scale utilities instead of "not implemented".

• Added multivariate_mme() — the first Phase-4 (multivariate Gaussian) slice: a supplied-covariance Henderson solve of the balanced multi-trait animal model (Kronecker MME with genetic precision Ainv⊗G0⁻¹ and residual I⊗R0⁻¹), returning per-trait fixed effects and EBVs. Validated against an independent loop-built MME, an independent marginal-GLS BLUP, the univariate animal model at t=1, and t independent single-trait fits at diagonal G0/R0. Experimental, dense/validation-scale; balanced data with a shared design only; no G0/R0 estimation, no missing records, no R-facing model-spec.

• Added genetic_correlation() — converts a covariance matrix (e.g. G0, R0) to the corresponding correlation matrix; also extracts the genetic correlation from a multivariate_mme result.

• Extended multivariate_mme() to unbalanced / missing-trait records: a missing or NaN entry in Y marks an unobserved trait, handled via the per-individual residual precision inv(R0[Sᵢ, Sᵢ]); balanced data reduces to the previous fast path. Validated against an independent loop-built MME and a marginal-GLS BLUP with per-individual residual blocks.

• Added fit_multivariate_reml() — estimates the multi-trait genetic/residual covariances G0, R0 by dense REML (log-Cholesky-parameterized Nelder–Mead; handles missing records), returning the estimates, correlations, per-trait heritabilities, and breeding values at the estimate. Validated by the t = 1 reduction to the univariate REML (the loglik equals sparse_reml_loglik exactly), grid-beating, and EBV consistency with multivariate_mme. Experimental, dense/validation-scale; the opt-in recovery harness accepts --seed or explicit --seeds lists and prints summaries outside CI, but is not broadly multi-seed calibrated and has no external-comparator parity yet.

• Added copy-returning Julia-side multivariate result accessors: variance_components, fixed_effects, breeding_values/EBV/BLUP, and REML-only heritability. These wrap existing multivariate_mme / fit_multivariate_reml fields and do not change result_payload() or the R bridge contract.

• Added test/fixtures/phase4_multitrait_parity/, a deterministic two-trait CSV fixture for R-lane sommer/ASReml/BLUPF90 parity work. It records a Julia REML target (G0, R0, beta, EBVs, h², and loglik) and CI checks fast self-consistency at the stored target covariances. Its README and docs/dev-log/decisions/2026-06-14-multitrait-comparator-protocol.md define the comparator protocol; it is not external comparator evidence.

• Added opt-in sim/phase4_multivariate_reml_recovery.jl for seeded two-trait known-truth recovery of the unstructured multivariate REML estimator outside CI. Default seed 20260616 passed with relative errors G = 0.174500 and R = 0.131056 against thresholds 0.25 and 0.20. This is internal recovery evidence, not multi-seed calibration or external comparator parity.

• Updated the validation-status documentation table to include the current Phase 2, Phase 3, Phase 4, and Phase 4B validation rows and their claim boundaries.

• Hardened the multivariate engine after a 7-lens adversarial review (each finding verified by running Julia): multivariate_mme / fit_multivariate_reml now reject non-finite observed phenotypes (Inf — only missing/NaN mark an unobserved trait), non-finite X/Z/Ainv, and empty-trait columns with clear errors (previously a silent all-NaN result or an opaque SingularException); and fit_multivariate_reml's loglik is now the full REML log-likelihood including the (N−p')·log(2π) constant, on the same scale as the other loglik functions (safe for LRT/AIC). Regression-tested.

• Added Phase-4B structured genetic covariance support for the dense multivariate REML engine: diagonal_covariance(), lowrank_covariance(), and factor_analytic_covariance() build diag(σ²), ΛΛ', and ΛΛ' + Ψ, and fit_multivariate_reml(...; genetic_structure = :diagonal | :lowrank | :factor_analytic, rank = K) estimates constrained genetic covariance structures while keeping residual R0 unstructured. Validated by deterministic constructor identities, structure metadata, copy-returning genetic_structure()/genetic_loadings()/genetic_uniqueness() accessors, loglik equality to the existing evaluator, PSD/PD covariance checks, constrained loglik ≤ the unstructured fit, and deterministic sign-canonicalization of returned loading columns. The 2026-06-14-loading-rotation-identifiability decision note records this as a sign-only metadata convention and defers full loading rotation/interpretation. Experimental, dense/validation-scale; no R-facing covariance-structure syntax, bridge change, full loading rotation/interpretation convention, or external comparator yet.

• Added opt-in Phase-4B recovery harness sim/phase4b_structured_covariance_recovery.jl for seeded low-rank and factor-analytic covariance recovery outside CI. It now accepts explicit --seeds lists and prints per-case summaries. This strengthens internal recovery tooling but does not add R-facing syntax, bridge payload fields, broad multi-seed calibration, or external comparator parity.

• Added docs/dev-log/decisions/2026-06-14-multivariate-recovery-calibration-protocol.md, a shared protocol for the unstructured and structured multivariate recovery harnesses. It defines the minimum seed counts, run-plan fields, summaries, and claim gate required before any broad multi-seed calibration claim.

• Executed that recovery calibration protocol outside CI on the predeclared seed lists and recorded negative calibration evidence in docs/dev-log/recovery-checkpoints/: unstructured passed 6/10, factor-analytic passed 8/10, and low-rank passed 9/10, with all fits converged. No broad multi-seed calibration claim is made.

• Added deterministic sim/summarize_recovery_calibration.jl tooling to parse committed recovery harness logs and regenerate the Markdown case summary, failure-mode table, and failed-seed list without rerunning stochastic simulations.

• Added docs/dev-log/recovery-checkpoints/2026-06-14-multivariate-recovery-calibration-failure-modes.md, which records that unstructured failures were 3 G-only plus 1 G+R, factor-analytic failures were 1 G-only plus 1 G+R, and low-rank had 1 R-only failure.

• Documented a throttled local command form for opt-in recovery/calibration harnesses (JULIA_NUM_THREADS=1, OPENBLAS_NUM_THREADS=1, OMP_NUM_THREADS=1, VECLIB_MAXIMUM_THREADS=1, nice -n 15) so long stochastic runs do not monopolize interactive workstations.

• Added docs/dev-log/decisions/2026-06-14-calibration-failure-response.md, which records that the failed predeclared calibration run cannot be rescued by silently relaxing thresholds, dropping failed seeds, or rerunning until the pass count improves.

• Added a "Multivariate models" documentation page with a runnable balanced two-trait example and the experimental / not-yet-R-wired boundary.

• Expanded planned backend marker/control vocabulary to include threaded CPU, AMDGPU, Metal, and oneAPI markers alongside CPU, CUDA, and auto metadata.

• Added backend_info() typed status diagnostics for planned backend rows with execution marked unavailable.

• Added planned genomic/QTL model-term vocabulary reservations: genomic(), single_step(), markers(), marker_scan(), and qtl_scan().

• Added planned standard quantitative-genetic model-term vocabulary reservations: permanent(), common_env(), maternal_genetic(), maternal_env(), paternal_genetic(), paternal_env(), cytoplasmic(), imprinting(), dominance(), epistasis(), relmat(), and HSquared.precision() in direct Julia code.

• Added a Documenter model-spec grammar page mirroring the R twin's status separation for parsed, reserved, and planned syntax.

• Added formula_status() grammar diagnostics and a Documenter status table mirroring the R twin's parsed/reserved/planned grammar rows.

• Added validation_status() diagnostics for covered, external, partial, and planned validation rows.

• Added max_dense_cells guards for the temporary dense validation path.

• Added experimental dense variance-component optimization for validated low-level animal-model specs.

• Added experimental low-level variance-component, fixed-effect, EBV/BLUP, fitted-value, and heritability extractors for the dense spec path.

• Switched breeding_values(fit) to the Henderson MME solve at the fit's variance components.

• Switched fitted_values(fit) to the same Henderson MME solve at the fit's variance components.

• Added experimental dense prediction-error-variance and reliability extractors for the dense spec path.

• Extended validation-scale prediction-error-variance and reliability extractors to supplied-variance HendersonMMEResult objects.

• Extended supplied-variance variance_components() and heritability() extractors to HendersonMMEResult objects.

• Added EBV() and BLUP() aliases for breeding_values(), plus accuracy() as a checked square-root transformation of reliability.

• Added fit_diagnostics() as metadata-only extraction for low-level AnimalModelFit and supplied-variance HendersonMMEResult objects.

• Added experimental direct payload fit_animal_model(y, X, Z, Ainv; ...) target for bridge-shaped inputs.

• Added explicit fit_animal_model(...; target = :henderson_mme, variance_components = ...) dispatch for supplied-variance Henderson MME solving.

• Added henderson_mme() for sparse Henderson mixed-model-equation solving at supplied variance components.

• Added a shared R/Julia Henderson mixed-model-equation validation fixture for the supplied-variance output path.

• Added result_payload() with field names aligned to the R hsquared_fit extractor contract.

• Added HSData, HSDataIDMap, and id_map() as an in-memory mirror of the R hs_data() input-container contract.

• Added HSData marker-map metadata validation and genotype-marker alignment checks.

• Added data_status(::HSData) diagnostics mirroring the R twin's data_status() surface for component presence, ID-overlap counts, pedigree status, genotype status, marker-alignment status, expression status, annotation-feature status, and environment-key status. Diagnostic only; no bridge payload, raw-pedigree Ainv construction, genotype parsing, relationship construction, annotation or environment covariate joins, automatic expression-feature joins, environmental model terms, marker scan, genomic fitting, eQTL/omics fitting, or QTL/eQTL claim.

• Recorded the R twin's hs_data() genotype-status diagnostics from hsquared head f067cd9 and mirrored them in Julia HSData as metadata diagnostics only.

• Recorded the R twin's hs_data() expression-status diagnostics from hsquared head 06cdf59 and mirrored them in Julia HSData as metadata diagnostics only.

• Recorded the R twin's hs_data() annotation-feature diagnostics from hsquared head 87888d9 and mirrored them in Julia HSData as metadata diagnostics only.

• Recorded the R twin's hs_data() environment-key diagnostics from hsquared head e7fbb31 and mirrored them in Julia HSData as metadata diagnostics only.

• Added sparse_csc_matrix() for R Matrix::dgCMatrix slot marshalling.

• Recorded the R twin's PEV/reliability bridge extractor contract while keeping Julia result_payload() fields unchanged.

• Recorded the R twin's tiny/local bridge enrichment of PEV/reliability from exported Julia extractors, still without widening base result_payload().

• Recorded the R twin's supplied-variance target = "henderson_mme" enrichment from prediction_error_variance(mme) and reliability(mme), still without widening base result_payload().

• Recorded the R twin's opt-in supplied-variance Henderson MME bridge target, with explicit no-log-likelihood/no-variance-estimation boundary.

• Recorded the R twin's sparse Z bridge marshalling handoff.

• Recorded the R twin's optional nadiv::Mrode9 pedigree-Ainv comparator evidence.

• Recorded the R twin's model_spec() preview surface for the v0.1 formula-to-bridge payload.

• Recorded the R twin's hs_data() parser integration for the v0.1 formula-to-bridge payload without changing the Julia payload shape.

• Recorded the R twin's animal(1 | id) shorthand for data = hs_data(..., pedigree = ped) as R-side formula ergonomics only; the explicit animal(1 | id, pedigree = ped) contract and Julia payload shape are unchanged.

• Added DocumenterVitepress documentation scaffold.

• Added audience and comparator programme notes.

• Added genomics/QTL/eQTL/GLLVM/GPU/HPC strategic roadmap.

• Added a backend and algorithm roadmap page for CPU, threads, CUDA, AMDGPU, Metal, oneAPI, AI-REML, Takahashi selected inversion, Woodbury paths, APY, and claim gates.

• Kept high-level fitting entry points as honest placeholders.