Can I fit and report this? • hsquared

This page answers one question: can I use hsquared for my analysis, and what am I allowed to write in the paper?

hsquared is an experimental package. Most of what follows is not about whether code runs; nearly all of it runs. It is about which numbers have enough evidence behind them to report.

The three-line version

Two default-routed covered paths exist: the univariate Gaussian animal model and t = 2 unstructured multivariate cbind() (validation-scale). Genomic GREML is also covered at validation scale, but only through its explicit engine = "julia", target = "genomic" route; ordinary-call / default-route genomic activation remains held (design-44 G5; BOUNDARY_HOLDOUT_FAIL). Report point estimates only within each stated scope.
Almost everything else fits through an opt-in engine target (some covered at validation scale, some experimental). k≥3 and diagonal multivariate stay experimental. Treat experimental fits as exploratory, or report them beside a comparator you ran yourself.
No standard error or interval in this package is nominally coverage-calibrated. Named univariate pedigree intervals have a narrower directional-conservative claim level, not a nominal-calibration claim. Do not report any interval as if nominal coverage had been established. The information matrix also differs by route: see Which information matrix?.

Three-rung validation ladder from planned to partial to covered. A side note marks internal self-agreement as insufficient for covered status.

Covered at validation scale · experimental · dense · REML-only still means experimental packaging. It is not 1.0 and not interval-calibrated. The ladder caption is dated to this source tree; the live count is validation_status() and this limits page after any covered-row update. This figure does not flip a row. public_covered_count stays 7.

What the status words mean

hsquared uses four words, in descending order of what you may claim. It deliberately never uses a bare “supported”, because that word hides the difference between runs and validated.

Word	The code	The evidence	What you may report
`covered`	Runs	Pre-declared recovery gate passed and an external same-estimand comparator agrees	Point estimates, within the stated scope
`partial`	Runs	Bridge-verified and partly evidenced; the chain has a named hole	Exploratory use; report only beside your own comparator
`experimental`	Runs	Labelling, not a status tier; every opt-in route carries it	Nothing on this package’s word alone
`planned`	Aborts as planned, not implemented	None	Nothing

This page is the user-facing list of what you may report. validation_status() is a developer evidence table of validation atoms; it is not a complete list of covered routes. Developers can inspect it in their own session:

validation_status()

Do I need Julia?

Yes, for fitting. This matters more here than in most R packages, because the default call is the one that needs it.

What you want	What you need
Parse and validate a model, inspect the design matrices and pedigree ordering	R only: `hs_control(engine = "validate")` and `model_spec()`
Actually fit and get heritability, variance components, EBVs	A local Julia, the `JuliaCall` package, and an `HSquared.jl` checkout

Without the engine, hsquared() still tells you whether your model is expressible, and errors with install guidance if you ask it to fit. See the Installation section of the README.

# no Julia required; same ped/dat as the README
hsquared(
  weight ~ sex + animal(1 | id, pedigree = ped),
  data = dat,
  control = hs_control(engine = "validate")
)

model_spec(
  weight ~ sex + animal(1 | id, pedigree = ped),
  data = dat
)

Route-by-route reporting permissions

The cards below are generated from validation_status(), a developer evidence table of validation atoms. They are not a complete list of covered routes — two covered routes have no row there yet and are listed on this page below the cards.

Status	Rows	What it means for you
`covered`	6	Pre-declared recovery gate passed and an external same-estimand comparator agrees. Point estimates are reportable within the stated scope.
`partial`	8	The code runs and is bridge-verified, but the evidence chain is incomplete. Exploratory use only.
`planned`	7	Not implemented. No reporting permission.

Those counts are rows in the evidence ledger, not user-facing models. Some rows are validation atoms (a pedigree-inverse check, a textbook fixture) or evidence for another row rather than a model you would fit. The cards below are the model routes.

Reader routes

Every route below runs. The question these cards answer is narrower and more useful: what may I put in a paper?

Univariate Gaussian animal model (the default path)

hsquared(y ~ sex + animal(1 | id, pedigree = ped), data = dat)

Can I fit it? Yes - implemented and covered on the default call.

Can I report the point estimate? Yes, within the stated scope - a pre-declared recovery gate passed and an external same-estimand comparator agrees.

Can I report an interval? Measured, but not nominally calibrated. Named univariate pedigree h2 delta/profile/bootstrap and sigma-a2 profile/bootstrap intervals have a directional-conservative claim level. Sigma-a2 delta/Wald is experimental-only. Do not report any of them as nominal 95% coverage.

Exact scope and caveat. One additive genetic effect, Gaussian response, REML only, fitted by average-information REML on the sparse mixed-model equations. Known-truth recovery is near-unbiased across an h2 grid of 0.2/0.4/0.6; the near-boundary cell h2 = 0.1 shows mild upward bias and 5% boundary pinning. ML is not implemented. REML = FALSE is rejected; use REML = TRUE (the default). That is the same rule for the default fit and for engine = "validate". The covered claim is this REML estimator, not ML. Named h2 delta/profile/bootstrap and sigma-a2 profile/bootstrap intervals are directional-conservative here; sigma-a2 delta/Wald remains experimental-only. None is nominally calibrated. This route has no dense-cell cap, but that is not a size-free claim: when additive variance is near 0 and the pedigree has more than 512 animals, the default fit can stop with converged = FALSE, optimizer_status boundary_score_unresolved, and an inflated residual. Check fit_diagnostics() after every large-pedigree fit.

Concrete fallback. This is the recommended route. If V_A is near 0 and the pedigree has more than 512 animals, retry with engine = “julia” and target = “sparse_reml” (see the sparse REML card). That fallback is experimental and opt-in; it does not change this route’s covered point-estimate claim. For a model with a second random effect, see the two-effect card below.

Animal model plus independent extra effects (common environment, arbitrary N)

hsquared(y ~ animal(1 | id, pedigree = ped) + common_env(1 | group),
         data = dat, control = hs_control(engine = "julia",
           engine_control = list(target = "two_effect")))
# more independent intercepts (hsquared#285):
hsquared(y ~ animal(1 | id, pedigree = ped) + (1 | nest) + (1 | year),
         data = dat, control = hs_control(engine = "julia",
           engine_control = list(target = "multi_effect")))

Can I fit it? Yes - implemented and covered at validation scale, behind an opt-in engine target.

Can I report the point estimate? Yes, within the stated scope - a pre-declared recovery gate passed and an external same-estimand comparator agrees.

Can I report an interval? No. Standard errors and intervals are asymptotic/delta-method, labelled experimental, and NOT coverage-calibrated. No route in this package currently carries an interval-reporting permission.

Exact scope and caveat. The common-environment leg (additive animal A plus i.i.d. environment, A2 = I) and the arbitrary-N generalization to independent i.i.d. effects are covered at validation scale on a pre-declared 48-seed recovery gate with an external same-estimand REML comparator. The default grammar accepts one named extra effect beside animal(); more independent intercepts use bare (1 | group) with target = "multi_effect". A maternal-genetic plus permanent-environment combination is not available. The animal-block ratio is narrow-sense h2; other blocks are variance-explained proportions, NOT heritabilities. Independent effects only - correlated, random-regression, and non-Gaussian structures are not this route.

Concrete fallback. The maternal-genetic leg (maternal_genetic(), A2 = pedigree A) runs on the same estimator but is still experimental: its own recovery gate and comparator are owed. Use the covered common-environment route, or treat a maternal fit as exploratory.

Sparse REML variance components (opt-in engine target)

hsquared(y ~ animal(1 | id, pedigree = ped), data = dat,
         control = hs_control(engine = "julia",
           engine_control = list(target = "sparse_reml")))

Can I fit it? Yes, but opt-in and experimental - the code runs; the evidence is incomplete.

Can I report the point estimate? Not on this package’s evidence alone. Report only beside a comparator you ran yourself, and say the route is experimental.

Exact scope and caveat. Reaches the same REML optimum as the default path and is cross-checked against an independent pure-R REML optimizer and the external pedigreemm package. It exists as an engine-target check, not as a separate user-facing model. It is also the existing fallback when the default AI-REML path stops unconverged on more than 512 animals with V_A near 0 (converged = FALSE, optimizer_status boundary_score_unresolved). That use is still experimental and opt-in; it does not newly cover this route.

Concrete fallback. Use the default hsquared() call, which is the covered route to the same estimates, except when that path stops on q > 512 animals with V_A near 0. Then this target is the existing fallback.

Repeatability / permanent environment (repeated records)

hsquared(y ~ animal(1 | id, pedigree = ped) + permanent(1 | id), data = dat,
         control = hs_control(engine = "julia",
           engine_control = list(target = "repeatability")))

Can I fit it? Yes, but opt-in and experimental - the code runs; the evidence is incomplete.

Can I report the point estimate? Not on this package’s evidence alone. Report only beside a comparator you ran yourself, and say the route is experimental.

Exact scope and caveat. REML-only repeatability optimizer returning additive, permanent-environment, and residual variances. Va and Vpe are only identifiable with repeated records per individual. No comparator and no known-truth recovery gate yet.

Concrete fallback. Fit the covered univariate animal model on individual means, or run your own comparator (for example sommer) alongside before reporting.

Genomic (GREML) and single-step relationship models

hsquared(y ~ genomic(1 | id, markers = M), data = dat,
         control = hs_control(engine = "julia",
           engine_control = list(target = "genomic")))

Can I fit it? Yes - implemented and covered at validation scale, behind an opt-in engine target.

Can I report the point estimate? Yes, within the stated scope - a pre-declared recovery gate passed and an external same-estimand comparator agrees.

Exact scope and caveat. 0.7 covered claim is opt-in genomic GREML only (target = “genomic”; genomic_variance_ratio on K_lambda; design-51/53). Accepts supplied Ginv or builds VanRaden1 G from markers. Ordinary-call / default-route activation remains held (design-44 G5; BOUNDARY_HOLDOUT_FAIL). Single-step / metafounder paths in this combined validation_status row remain exploratory (not the 0.7 covered claim). APY / production / intervals out of scope.

Concrete fallback. Report genomic_variance_ratio on the declared kernel with the no-anchor and scale disclosures. Do not treat single-step or default genomic routing as covered.

SNP-BLUP marker effects

hsquared(y ~ genomic(1 | id, markers = M), data = dat,
         control = hs_control(engine = "julia",
           engine_control = list(target = "snp_blup")))

Can I fit it? Yes, but opt-in and experimental - the code runs; the evidence is incomplete.

Can I report the point estimate? Not on this package’s evidence alone. Report only beside a comparator you ran yourself, and say the route is experimental.

Exact scope and caveat. VanRaden method-1 marker model returning per-marker effects, genomic breeding values, and fixed effects, with variances either supplied or REML-estimated. Weighted/Bayesian marker priors and JWAS/sommer/BLUPF90 parity are planned.

Concrete fallback. Use the genomic GREML route above for a breeding-value question, and treat marker effects as exploratory.

Multivariate Gaussian animal model (cbind response)

hsquared(cbind(y1, y2) ~ animal(1 | id, pedigree = ped), data = dat)
# routes on the default path; engine = "julia", target = "multivariate" still works

Can I fit it? Yes - implemented and covered on the default call.

Can I report the point estimate? Yes, within the stated scope - a pre-declared recovery gate passed and an external same-estimand comparator agrees.

Exact scope and caveat. REML-only, animal-model-only, dense/validation-scale. Returns G0/R0 covariance and correlation matrices, per-trait h2, and cross-trait EBVs. Covered claim is t = 2 unstructured G0/R0 only (2026-09-02 maintainer sign-off); k>=3 and genetic_structure = “diagonal” stay experimental. The R-lane evidence includes a 100-replicate t = 2 cold-start recovery study, one reproduced full-unstructured sommer comparator leg, a published Mrode-style supplied-variance anchor, and a Bayesian MCMCglmm agreement probe - which is NOT same-estimand REML parity. The Mrode anchor is SUPPLIED-covariance (fixed effects and BLUPs given known G0/R0), so it does NOT anchor the ESTIMATED G0/R0, and no published textbook anchor for estimated multivariate G0/R0 exists - that evidence is comparator- and recovery-based, not textbook-anchored. R<->engine element-wise parity is verified LOCALLY at pre-declared tolerances, NOT in CI (the CI runner has no Julia, so those checks skip). The engine row is covered; this R-public surface is now covered at validation scale; public_covered_count is 7. Since default routing the cbind route is selected on the default call.

Concrete fallback. For k>=3 or a diagonal genetic structure, treat the fit as experimental. For a single-trait question, use the covered univariate model.

Planned validation lanes and grammar

These rows are planned validation tasks in validation_status(), not a blanket statement that every related fit aborts. Check the partial routes on this page for bounded FA and genetic GLLVM fits. A planned row by itself gives no reporting permission.

Mrode fitted animal-model outputs (Phase 1)
ASReml comparison policy (Phase 1)
BLUPF90/DMU/WOMBAT comparison policy (Phase 1)
XSim simulation truth (Phase 5+)
genomic and QTL/eQTL validation (Phase 5+)
GLLVM-style multivariate validation (Phase 6)
CPU/GPU backend comparison (Phase 7+)

Before you report

Confirm this page still lists your route as reportable. Do not treat validation_status() as the complete user list.
Report point estimates for covered routes only.
Do not report a standard error or interval as nominally calibrated. The named univariate-pedigree directional-conservative methods are a narrower claim, not a nominal-coverage permission.
For a partial route, run an external comparator yourself and report both numbers.
Say in the methods that the fit came from an experimental package, and name the version.

Two covered routes without their own status row

validation_status() does not yet carry a separate row for two routes that docs/design/capability-status.md records as covered at validation scale:

Random regression / reaction norm, k = 2: animal(rr(age, order = 2) | id, pedigree = ped) via engine_control = list(target = "random_regression"). A pre-declared 48-seed recovery gate passed and a sommer leg() same-estimand comparator agrees on all K_g entries and the residual. k = 2 only. The residual is homogeneous and there is no permanent-environment term, so rr_heritability() can overstate h²(t) for repeated-records designs. K_g values are not comparable across basis-normalization conventions. Rescale the response to about unit variance before fitting; otherwise a large-scale trait can report converged = TRUE below the REML optimum and a wrong covariance function.
Direct–maternal correlated 2×2 G: animal(1 | id, pedigree = ped) + maternal_genetic(1 | dam) via engine_control = list(target = "direct_maternal"). A pre-declared 48-seed gate passed and a sommer covm() same-estimand comparator agrees. heritability() returns the labelled Willham triple: direct h²_d, maternal m², total h²_T, and r_am; never a bare scalar. Direct h² is not the total heritability, and a negative r_am is real and expected. Single relationship A only.

Both are opt-in, experimental, and REML-only. The engine-enforced dense rule on the default animal() route and target = "repeatability" is nobs^2 + nanimals^2 <= 1e6 (default max_dense_cells), counting every animal in the pedigree, not just records. Unlike those guarded routes, neither random regression nor direct-maternal currently accepts engine_control$max_dense_cells, so they have no engine-enforced cell cap. Read this page before reporting a number; validation_status() is a developer evidence table. The cards above are generated from the evidence table and will not mention them until those rows exist.

New partial FA and genetic GLLVM routes

The four-trait, one-factor Gaussian pedigree FA fit is available through an explicit Julia expert-control target. It returns a fitted genetic covariance, uniqueness estimates, genetic correlations, and diagnostics. The rank does not by itself identify every loading and uniqueness configuration. The tested fixture required a higher iteration limit. In a 200-seed Gaussian pedigree study with four traits and three records per animal, 141 fits converged and 110 met the declared covariance-recovery diagnostics, using two ordinary starts with 5000 iterations each. All 200 attempts were counted. Inspect convergence, start disagreement and floor diagnostics before interpreting a fit. See Multivariate Gaussian animal models for the exact call and the study’s limits.

The three-trait, two-factor Poisson-log pedigree genetic GLLVM fit is also an explicit expert-control route. It returns genetic covariance and correlations plus trait-level conditional genetic modes on the log-rate scale. See Genetic GLLVM opt-in for the exact call and limits.

Both routes are partial. Use their point estimates for exploration; for a paper, run a suitable comparator yourself. No interval, automatic factor rank, missing-response fit, or broader family claim follows from these cells.

Which information matrix?

SE and interval extractors do not all use the same information matrix. Swapping one matrix for the other changes the SE (the two are not the same estimand). The method itself is unchanged; this table only names what each route already computes.

Route	Information matrix	How it is built	Ratio SE / interval
Univariate animal model (default)	Average information	Analytic AI matrix from the REML fit	Delta method; logit scale for `heritability_interval()`
Repeatability / permanent environment	Observed information	Central finite-difference Hessian, `fd_step = 1e-4`	Logit-scale delta method for `repeatability_interval()`
Two-effect (`common_env` / maternal)	Observed information	Central finite-difference Hessian, `fd_step = 1e-4`	Logit-scale delta method for `heritability_interval()` and `common_env_proportion_interval()`
Multivariate `cbind()`	Observed information	Central finite-difference Hessian, `fd_step = 1e-4`	Delta method for `covariance_standard_errors()`; Fisher-z for genetic-correlation endpoints

See ?variance_component_standard_errors, ?heritability_interval, ?repeatability_interval, ?common_env_proportion_interval, and ?covariance_standard_errors. None of these intervals is nominally coverage-calibrated.

Testing V_A = 0 and boundary fits

Additive variance sits on the edge of the parameter space (V_A >= 0). The Wald SE and the logit-delta CI are not valid when V_A is near 0. The likelihood-ratio test of V_A = 0 uses a 50:50 mixture of chi^2_0 and chi^2_1 (Self and Liang 1987; Stram and Lee 1994), so a plain chi^2_1 p-value is twice too large. R anova() does not run that test. logLik() and AIC() are exposed for converged fits, so a reader can build a naive chi^2_1 LRT by hand – do not. The Julia engine already implements the mixture as HSquared.nested_lrt(loglik_constrained, loglik_full; df = 1, boundary_df = 1). That names an existing helper; it is not a new R test statistic.

Where to read next

Fitting models: the smallest honest workflow after you choose a route.
Visualizing models: diagnose a fitted result before extracting it.
Model status: the same ledger organised by model family.
Validation evidence: the evidence chain behind each row, atom by atom.
Function map and cheat sheet: which function to call at which step.
Formula grammar roadmap: what the parser reserves for later, and what it refuses today.
Progression & evidence: dated milestones that do not themselves change a reporting permission.