Model status

Status ledger, not a green checklist. What you may report is on Can I fit and report this?. Covered-at-validation-scale rows remain experimental packaging (experimental 0.9.0 release; public_covered_count is 7; CRAN availability is tracked separately). Julia engine-covered ≠ R covered: a bounded factor-analytic G route and a bounded genetic GLLVM route are opt-in partial on R; single-step also stays opt-in partial.

This page separates what exists from what is planned.

Exists now

The default hsquared() call fits the v0.1 univariate Gaussian animal model (y ~ fixed + animal(1 | id, pedigree = ped), REML) through the HSquared.jl engine, validated by known-truth recovery, the published gryphon REML anchor, and sommer agreement (see the V0.1 contract). Fitting requires a local Julia + HSquared.jl. The remaining items below are the surrounding surface: parser, validation atoms, data container, extractors, and advanced opt-in engine controls.

• Fits the v0.1 univariate Gaussian animal model by default and returns variance components, heritability, breeding values, fixed effects, fitted values, residuals, and fit_diagnostics() (with an at_boundary flag).
• R package scaffold and CI.
• Team operating memory and claim registers.
• hs_control() with default engine = "fit", engine = "validate" for a no-fit preview, and engine = "julia" for advanced engine targets.
• hs_control() stores planned CPU-thread, CUDA, AMDGPU, Metal, and oneAPI backend names as control metadata.
• backend_info() reports the planned backend vocabulary and marks backend execution unavailable.
• formula_status() reports parsed, reserved, and planned grammar terms without implying fitted support.
• validation_status() reports current validation atoms and planned comparator lanes without running checks or implying fitted support.
• data_status() reports hs_data() component presence, ID overlap, and pedigree, expression-feature, genotype-column, marker-map/genotype-marker alignment, annotation-feature, and environment-key diagnostics without implying fitted support.
• model_spec() previews the parsed v0.1 animal-model contract, including fixed-effect design columns, sparse animal-effect design dimensions, normalized pedigree ordering, and Julia target metadata. It does not fit a model.
• Narrow univariate Gaussian REML genomic() (supplied Ginv or a markers matrix) is covered at validation scale through explicit engine = "julia", target = "genomic"; it is not a default route. single_step() (supplied Hinv, or constructed Hinv from pedigree + markers), permanent(), common_env(), maternal_genetic(), and the supplied-relationship primaries relmat() (dense K) and precision() (supplied inverse Q) remain explicit engine = "julia" targets; some covered at validation scale, some experimental (see below). relmat()/precision() are REML-only, experimental, and not covered (they do not change the covered-model count).
• markers(), marker_scan(), and qtl_scan() reserve planned genomic/QTL formula vocabulary and currently error as not implemented.
• maternal_env(), paternal_genetic(), paternal_env(), cytoplasmic(), imprinting(), dominance(), and epistasis() reserve standard quantitative-genetic and inheritance-kernel vocabulary and currently error as not implemented.
• animal() as an inert formula marker.
• A v0.1 parser for animal(1 | id, pedigree = ped), and for animal(1 | id) when data is an hs_data() bundle with a pedigree component.
• A tested internal R-to-Julia payload shape with y, X, sparse Z, method, family, encoded IDs, normalized pedigree metadata, and Julia target metadata.
• An experimental opt-in control = hs_control(engine = "julia") path that can send the tiny payload into a sibling HSquared.jl checkout, marshal sparse Z through CSC slots, and normalize the returned Julia result into the internal hsquared_fit contract. When the sibling Julia package exposes dense validation extractors, this path also enriches the result with PEV/reliability fields.
• An explicit supplied-variance bridge target, engine_control = list(target = "henderson_mme", variance_components = ...), that calls Julia henderson_mme() for tiny validation examples. It does not estimate variance components and does not provide a log-likelihood. When the sibling Julia checkout exposes applicable dense validation extractors, it can attach PEV/reliability fields.
• An internal tiny validation fixture that pins R payload ordering, sparse Z, and live Julia pedigree_inverse() agreement for a three-animal pedigree when the local Julia bridge is available.
• An optional Mrode9/nadiv pedigree-Ainv comparator test when nadiv and a sibling HSquared.jl checkout are available.
• An internal supplied-variance Henderson MME validation fixture that compares an independent R reference solve with Julia henderson_mme() for fixed effects, EBVs, fitted values, h2, and optional dense validation-path PEV/reliability when a sibling HSquared.jl checkout is available.
• A tiny supplied-variance likelihood fixture that compares Julia dense REML and sparse REML evaluators, plus an ML hand-check target, when a sibling HSquared.jl checkout is available.
• A Mrode-style supplied-variance output fixture that pins Ainv, fixed effects, EBVs, fitted values, PEV, reliability, h2, ML log-likelihood, and dense/sparse REML log-likelihood against independent R references and the sibling Julia engine when available.
• The average-information REML optimizer (HSquared.fit_ai_reml()) that the default engine = "fit" uses to estimate the variance components. It is also reachable explicitly via control = hs_control(engine = "julia", engine_control = list(target = "ai_reml")), records variance_components_source = "estimated_ai_reml", recovers known truth in the DGP recovery study (near-unbiased), and recovers the published gryphon REML estimate within the maintainer-signed-off comparator band.
• The sparse NelderMead REML optimizer (HSquared.fit_sparse_reml()), reachable via control = hs_control(engine = "julia", engine_control = list(target = "sparse_reml")). It records variance_components_source = "estimated_sparse_reml", agrees with the average-information and dense optimizers, and recovers the published gryphon anchor.
• The first fitted-object/extractor contract over internal hsquared_fit objects and mocked Julia result fields, including variance components, heritability, EBV/BLUP aliases, PEV, reliability, accuracy, fixed effects, random effects, log-likelihood, AIC, prediction, fitted values, residuals, summaries, coef(), nobs(), and fit_diagnostics() for convergence, optimizer, target, iteration, and dense-validation-path metadata.
• Marker/QTL/eQTL extractor names: marker_effects(), marker_variance_explained(), qtl_table(), gwas_table(), eqtl_table(), and lod_scores(). marker_effects() returns the per-marker effects of an opt-in SNP-BLUP fit, and marker_variance_explained() returns descriptive fitted-marker shares for that same path (see above). gwas_table(scan) and lod_scores(scan) work for an already-computed hs_gwas object, but they are uncalibrated and not map-annotated. Fit-level QTL/GWAS/eQTL tables still require future hsquared_fit result fields.
• hs_data() as a lightweight input container with ID maps for phenotype, pedigree, genotype, expression, marker, annotation, and environment inputs.
• summary(hs_data(...)) reports ID overlap and mismatch counts for phenotype, pedigree, genotype, and expression components.
• summary(hs_data(...)) reports pedigree coverage, founder, parent-link, and conservative warning counts when a pedigree component is supplied.
• summary(hs_data(...)) reports marker-map and genotype-marker alignment diagnostics when marker or genotype components are supplied.
• hs_data() validates supplied marker maps for marker ID, chromosome, and non-negative numeric position columns. This is metadata validation only.
• When both genotypes and markers are supplied, hs_data() checks that genotype marker columns match marker-map IDs exactly.
• hs_data(expression = ..., annotation = ..., annotation_id = ...) checks expression feature columns against annotation rows and reports duplicate or unmatched annotation keys. This is metadata validation only.
• hs_data(environment = ..., environment_id = ...) checks environment metadata coverage against phenotype records and reports duplicate or unmatched environment keys. This is metadata validation only.
• model_spec() and hsquared() can use an hs_data() object for the v0.1 animal-model parser, reading model variables from phenotypes and using the bundle pedigree by default when animal(1 | id) omits pedigree =.
• Local tests for accepted syntax, rejected future syntax, and pedigree/data ID checks.

Covered at validation scale (opt-in unless noted)

These models are REML only, Julia-owned (R surfaces them), dense/validation-scale (n ≤ ~1000), and not production. Each has cleared a pre-declared bias/MCSE recovery gate and a same-estimand REML comparator (mirroring a covered gate in the HSquared.jl twin). Several remain explicit opt-in control = hs_control(engine = "julia", engine_control = list(target = ...)) targets. Multivariate cbind() is the exception: it routes on the default engine = "fit" call. Narrow genomic GREML is covered at validation scale only through its explicit target. Under engine = "julia", the default target "fit_animal_model" and target = "repeatability" enforce this scale as an engine-side dense-cell cap (default 1e6 = nobs^2 + nanimals^2; raise it with engine_control = list(max_dense_cells = ...), or switch to a sparse target). The default engine = "fit" path routes to the sparse-capable AI-REML fitter and enforces no such cap, and neither do the other routes below, so for all of those the n ≤ ~1000 stays a recommendation only.

• Multivariate Gaussian animal model: cbind(trait1, trait2) ~ animal(1 | id, pedigree = ped), routed on the default path (no target needed; target = "multivariate" still accepted). Returns G/R covariance and correlation matrices, per-trait h2, and cross-trait EBVs; missing trait cells may be NA. Covered claim is t=2 unstructured G0/R0 (2026-09-02); k≥3 and genetic_structure = "diagonal" stay experimental. R-lane evidence includes cold-start recovery, one full-unstructured sommer REML comparator, a published Mrode supplied-covariance BLUP/MME anchor, and a MCMCglmm Bayesian agreement probe. The Mrode anchor is supplied-covariance ; it checks fixed effects and BLUPs given known G0/R0, not the estimated covariances. No published textbook anchor for estimated multivariate G0/R0 exists, so that part of the evidence is comparator- and recovery-based rather than textbook-anchored. R↔︎engine element-wise parity is verified locally at pre-declared tolerances, not in CI. MCMCglmm is not REML parity. Intervals remain uncalibrated.
• Genomic GREML: genomic(1 | id, markers = M) or genomic(1 | id, Ginv = Q) through control = hs_control(engine = "julia", engine_control = list(target = "genomic")). Covered at validation scale, experimental, and opt-in; it is not an ordinary default route. The experimental 0.9.0 release label is retained; public_covered_count stays 7. Marker input means sample-frequency VanRaden1 with K_lambda = G + 0.01I; supplied Ginv construction is unknown. heritability() returns the genomic variance-component ratio sigma_g2 / (sigma_g2 + sigma_e2) on the declared relationship scale, not generally pedigree or population h². Interval/SE accessors are unavailable. The nine-cell recovery campaign is confirmatory history and is not claimed PASS here. Single-step stays opt-in partial; SNP-BLUP stays opt-in experimental.
• Common environment: animal(1 | id) + common_env(1 | group), target = "two_effect" (additive animal genetic + IID common environment, A2 = I). The common-environment leg is covered (mirrors the twin V3-TWOEFFECT-REML: a 48-seed bias/MCSE gate PASS + a blupf90+ same-estimand comparator ~1e-5, sommer cross-check ~2e-5); its h²/c² interval is asymptotic and not coverage-calibrated. The maternal genetic leg of the same two_effect target stays experimental (listed below).
• Arbitrary-N independent (1 | g): animal(1 | id, pedigree = ped) + (1 | g1) + (1 | g2) + ..., target = "multi_effect": one pedigree animal component plus any number of INDEPENDENT i.i.d. effects. Covered (mirrors the twin V3-NEFFECT-REML: a 48-seed bias/MCSE gate PASS + a sommer same- estimand comparator, with exact live R–Julia parity). The animal-block ratio is narrow-sense h²; the other blocks are variance-explained proportions, not heritabilities. INDEPENDENT effects only; random slopes (x | g) and correlated (x || g) terms are rejected.
• Random regression (reaction norm), k = 2: animal(rr(t, order = 2) | id, pedigree = ped), target = "random_regression": a linear reaction norm returning the 2×2 genetic covariance K_g (intercept × slope), residual variance, EBVs, and an h²(t) trajectory. Covered at k = 2 (mirrors the twin random-regression REML validation row: a 48-seed bias/MCSE gate PASS
  • a sommer leg() same-estimand comparator AGREE, exact live parity). rr_heritability(), not heritability(), is the reaction-norm heritability accessor: it returns h²(t) as a curve, never a bare scalar. heritability() on a random-regression fit errors and names rr_heritability(), because a single scalar h² is not defined for this model. Point estimate only (no interval on the curve); the residual is homogeneous and there is no permanent-environment term, so h²(t) can overstate for repeated-records designs. k ≥ 3 is experimental.
• Direct–maternal correlated model (2×2 G): animal(1 | id, pedigree = ped) + maternal_genetic(1 | dam), target = "direct_maternal", the first correlated random-effect model (estimates the direct–maternal covariance σ_dm). Covered at validation scale (mirrors the twin V4-DIRECT-MATERNAL: a pre-declared 48-seed bias/MCSE gate PASSED + a sommer covm() same-estimand REML comparator AGREE), still opt-in and dense/validation-scale (n ≤ ~1000). heritability() returns the Willham labelled triple (direct h²_d, maternal m², total h²_T, r_am), never a bare scalar; a negative r_am is real and expected.

Opt-in and experimental (not the default)

These models are not yet comparator/known-truth-validated at the bar above; each mirrors a partial gate in the HSquared.jl twin. They are reached through the opt-in control = hs_control(engine = "julia", ...) path.

Narrow genomic GREML is covered at validation scale through its explicit target = "genomic" route. SNP-BLUP and single-step remain explicit opt-in experimental targets.

• Repeatability / permanent environment: animal(1 | id) + permanent(1 | id), target = "repeatability" (needs repeated records per individual). Dense, validation-scale (n ≤ ~1000); the engine enforces this as a dense-cell cap (nobs^2 + nanimals^2, default 1e6), raisable via engine_control = list(max_dense_cells = ...), and exceeding it raises an hsquared_error rather than running unboundedly (hsquared#217).
• Maternal genetic (INDEPENDENT two-effect leg): animal(1 | id) + maternal_genetic(1 | dam), target = "two_effect" (additive + pedigree maternal effect, treated as independent). Uses the same estimator as the covered common-environment leg with exact live parity, but its own recovery gate and comparator on the maternal-A2 design are still owed, so it stays experimental. (The correlated direct–maternal model above is a different, covered target.)
• Single-step: single_step(1 | id, Hinv = Hinv) on a supplied single-step inverse, target = "single_step", or single_step(1 | id, pedigree = ped, markers = M), target = "single_step_construct", where the engine builds Hinv from the pedigree and genotyped-subset marker matrix.
• SNP-BLUP / RR-BLUP marker effects: genomic(1 | id, markers = M), target = "snp_blup". Returns per-marker effects via marker_effects() plus per-individual genomic breeding values. variance_components is optional: with supplied c(sigma_g2 = ..., sigma_e2 = ...) it solves at those variances; when omitted it estimates sigma_g2/sigma_e2 by REML from the markers. Mirrors the twin V2-SNPBLUP gate.

Current limits

How many of each term can a formula carry today? The default grammar is narrow, while the opt-in multi-effect route accepts any number of independent random intercepts. These are the structural limits:

• Exactly one primary effect. A formula carries exactly one of animal(1 | id, pedigree = ped) (a pedigree additive effect) or a supplied-relationship primary (genomic(1 | id, ...) / single_step(1 | id, Hinv = ...)). Two primary effects is an error.
• At most one additional random effect on the default grammar, and only alongside an animal() primary: one of permanent(1 | id), common_env(1 | group), or maternal_genetic(1 | dam). The opt-in target = "multi_effect" path lifts this to an arbitrary number of INDEPENDENT (1 | g) intercepts (covered at validation scale; see the opt-in sections above). Correlated multi-term structures other than the covered direct–maternal (2×2 G) target remain on the roadmap.
• Default grammar accepts random intercepts only. Ordinary 1 | group terms are parsed where the target permits them; raw random slopes (x | id) are rejected. The opt-in animal(rr(t, order = 2) | id, ...) route is the separate bounded random-regression path. On the default path a bare lme4-style (... | group) term is rejected with a pointer to the named effects above (it is never silently absorbed into the fixed effects); the opt-in target = "multi_effect" path is where bare (1 | g) intercepts are accepted.
• One numeric response, or a cbind() multi-trait response. A multi-trait cbind(...) Gaussian response with an animal() term routes to the multivariate fitter on the default path; the explicit engine = "julia", target = "multivariate" spelling still works. The route is default; t=2 unstructured is covered at validation scale (experimental label retained).
• The default route is Gaussian with an identity link. Non-Gaussian Poisson/binomial paths exist only as explicit experimental Julia targets; they are not default-route or reporting permissions.
• REML only. REML = FALSE (ML) is rejected; use REML = TRUE (the default). That is the same rule for the default fit and for engine = "validate". The supplied-variance paths (henderson_mme, and snp_blup when variances are supplied) solve at given variances rather than estimating them (snp_blup with variances omitted estimates them by REML).
• Fixed effects are unrestricted otherwise: any number of ordinary fixed-effect terms (sex + age + ...) on the right-hand side is fine.

General combinations of random slopes and correlated random terms, and formula-level structured covariance models remain on the roadmap. The opt-in independent multi-effect, k=2 random-regression, and bounded four-trait FA expert-control routes have the separate limits described on this page.

Not implemented yet

• General default model fitting beyond the v0.1 univariate Gaussian animal model and the default-routed t=2 unstructured multivariate cbind() path. Genomic GREML, repeatability, SNP-BLUP, single-step, and non-Gaussian (poisson(log)/binomial(logit), Laplace marginal likelihood or a hybrid variational-plus-Laplace objective, experimental conditional three-field h2 outputs) models fit only through explicit engine = "julia" targets; some covered at validation scale, some experimental (see the two sections above).
• ML estimation (REML = FALSE is rejected on the default fit and on engine = "validate"; only REML is implemented).
• R-side Ainv construction (the engine builds Ainv in Julia).
• Estimated variance components, heritability, EBVs, or BLUPs as a default for any model other than the v0.1 univariate Gaussian animal model or the default-routed t=2 unstructured multivariate path (the explicit models above return them; covered ones are validated at validation scale only).
• Log-likelihood or information criteria for supplied-variance Henderson MME bridge results.
• Production sparse PEV or reliability.
• Full Mrode animal-model fit-output validation with estimated variance components.
• Broad ASReml, DMU, WOMBAT, or additional BLUPF90 validation beyond the scoped comparator runs cited above.
• File-backed genotype/omics loading or streaming computation.
• Automatic model construction from genotype, marker, expression, annotation, or environment components in hs_data().
• Environmental random effects or multi-environment model terms from environment metadata.
• Automatic feature annotation joins, eQTL scans, or omics models from expression and annotation metadata.
• Allele coding, marker imputation, PLINK/VCF parsing, or marker scanning from marker maps.
• QTL, GWAS, or eQTL result generation, and marker-scan / marker-imputation results (the opt-in SNP-BLUP path above does return marker effects via marker_effects()).
• Production genomic or single-step workflows; genomic GREML is covered at validation scale only through its explicit target, and single-step remains opt-in partial (explicit engine = "julia" target; not R covered).
• Factor-analytic and structured cov() G matrices stay planned in the R formula (cov = fa(K) is not activated). A bounded four-trait, one-factor Gaussian pedigree FA fit is available through explicit expert controls and remains partial. Julia V4-FA engine status is not a broader R claim. The unstructured t=2 multivariate G/R path is covered at validation scale on the default call; see above.
• Paternal, dominance, epistasis, custom relationship/precision, QTL-style, selfing, clonal, haplodiploid, polyploid, cytoplasmic, imprinting, and broader GLLVM-style models. A bounded Poisson-log, three-trait, two-factor pedigree genetic GLLVM fit is now available through explicit expert controls and remains partial; see Genetic GLLVM opt-in. (Permanent environment, common environment, and maternal genetic effects fit through the opt-in targets above; common environment is covered at validation scale, the others experimental.)
• GPU execution.

Comparator targets

sommer and pedigreemm are already in use as v0.1 comparators: sommer is the signed-off V1-COMPARATORS agreement check on the gryphon anchor, and pedigreemm provides a one-sided REML-log-likelihood floor. The broader long-term comparator set additionally includes ASReml, MCMCglmm, BLUPF90, DMU, WOMBAT, JWAS.jl, XSim.jl, AGHmatrix, nadiv, drmTMB, gllvmTMB, DRM.jl, and GLLVM.jl.

Performance and coverage claims are evidence-gated. Public pages may call a feature working only after code, tests, documentation, and validation evidence exist.