Package {hsquared}

Title:	'R' Interface for 'Julia'-Backed Quantitative-Genetic Models
Version:	0.9.0
Description:	Provides an interface for quantitative-genetic models in 'R', backed by the 'HSquared.jl' engine for 'Julia'. Fits Gaussian pedigree animal models by restricted maximum likelihood, using average-information restricted maximum likelihood (AI-REML; Gilmour et al. (1995) <doi:10.2307/2533274>) for the univariate model. Returns heritability, variance components, breeding values, and covariance and correlation matrices for multivariate models. Formula and pedigree validation can run without 'Julia'. Model fitting requires a local 'Julia' and 'HSquared.jl'; use hs_control(engine = "validate") to preview the model contract without fitting. What you may report is listed on the package article 'Can I fit and report this?', not in validation_status(), which is a developer evidence table. Model-specific restrictions are documented in that article. Uncertainty intervals are not coverage-calibrated.
License:	MIT + file LICENSE
URL:	https://itchyshin.github.io/hsquared/, https://github.com/itchyshin/hsquared, https://github.com/itchyshin/HSquared.jl
BugReports:	https://github.com/itchyshin/hsquared/issues
Encoding:	UTF-8
Language:	en-US
Imports:	ggplot2, lifecycle, Matrix, stats
Suggests:	enhancer, JuliaCall, knitr, lme4, nadiv, pedigreemm, rmarkdown, sommer, testthat (≥ 3.0.0), withr
VignetteBuilder:	knitr
Config/testthat/edition:	3
Config/roxygen2/version:	8.1.0
NeedsCompilation:	no
Packaged:	2026-10-09 12:31:53 UTC; z3437171
Author:	Shinichi Nakagawa [aut, cre, cph], Yefeng Yang [aut, cph], Szymon Drobniak [aut, cph]
Maintainer:	Shinichi Nakagawa <itchyshin@gmail.com>
Repository:	CRAN
Date/Publication:	2026-10-21 13:10:26 UTC

hsquared: R Interface for Julia-Backed Quantitative-Genetic Models

Description

The hsquared package is the R-facing interface for heritability, breeding-value, G-matrix, and inheritance-structured quantitative-genetic models backed by the HSquared.jl Julia engine.

Details

v0.1 fits the univariate Gaussian animal model y ~ fixed + animal(1 | id, pedigree = ped) by REML (average-information) through the HSquared.jl engine: the default hsquared() call fits when a local Julia and HSquared.jl are available, and otherwise errors with install guidance. A cbind() multivariate Gaussian response also routes on that default path; multivariate is covered at validation scale and stays experimental. Narrow genomic GREML is an explicit experimental target (target = "genomic") and is covered at validation scale. The covered common-environment two-effect, independent multi-effect, direct-maternal, and k = 2 random-regression models remain opt-in. Repeatability, maternal two-effect, and non-Gaussian (poisson(log)/binomial(logit), Laplace marginal likelihood or variational ELBO) models are opt-in and experimental; single-step stays opt-in partial. Factor-analytic G is engine-covered on the Julia engine (V4-FA) and planned, not fitted, on the R formula (Julia engine-covered is not R covered).

Current limitations

This package is an experimental 0.9.0 release; it is not production or 1.0.0, and CRAN availability is tracked separately. Default hsquared() fitting requires a local Julia installation and HSquared.jl; use hs_control() with engine = "validate" to check the model contract without fitting. Report point estimates only for covered routes; uncertainty intervals are experimental and not coverage-calibrated.

What you may report is listed on Can I fit and report this? – not in validation_status(). That article lists every covered route and its reporting scope. validation_status() is a developer evidence table of validation atoms; it is not the user-facing list. The model status and validation evidence articles are published on the package website rather than as installed vignettes.

Author(s)

Maintainer: Shinichi Nakagawa itchyshin@gmail.com (ORCID) [copyright holder]

Authors:

• Shinichi Nakagawa itchyshin@gmail.com (ORCID) [copyright holder]

• Yefeng Yang (ORCID) [copyright holder]

• Szymon Drobniak (ORCID) [copyright holder]

See Also

Useful links:

• https://itchyshin.github.io/hsquared/

• https://github.com/itchyshin/hsquared

• https://github.com/itchyshin/HSquared.jl

• Report bugs at https://github.com/itchyshin/hsquared/issues

Animal-model formula marker

Description

Usage

animal(formula, pedigree = NULL, ...)

Arguments

Details

animal() is the additive-genetic term in an hsquared() formula. Write animal(1 | id, pedigree = ped), or animal(1 | id) when data is an hs_data() object with a pedigree component. Fitting happens in hsquared(), not here: calling animal() on its own is a syntax marker and returns NULL.

Value

NULL, invisibly. The call is interpreted by hsquared() when it appears inside a model formula.

Examples

# n = 4 is a syntax demo, not a number for a paper.
ped <- data.frame(
  id = c("sire", "dam", "off1", "off2"),
  sire = c(NA, NA, "sire", "sire"),
  dam = c(NA, NA, "dam", "dam")
)
dat <- data.frame(
  id = c("sire", "dam", "off1", "off2"),
  sex = c("m", "f", "m", "f"),
  weight = c(42, 38, 40, 37)
)
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
)

Inspect planned compute backends

Description

Usage

backend_info(control = hs_control())

Arguments

Details

backend_info() reports which backend names are accepted by hs_control() and whether any of them are execution-ready from the R package. In the current package state, backend names are control metadata only: they are selectable but not dispatched.

Value

A data frame of backend status records with class "hs_backend_info".

Extract breeding values

Description

Usage

breeding_values(object, ...)

EBV(object, ...)

BLUP(object, ...)

Arguments

Details

breeding_values() is part of the planned v0.1 fitted-object contract. It works for hsquared_fit objects that contain a Julia result. EBV() and BLUP() are aliases for applied quantitative-genetic workflows.

Value

Breeding value results for hsquared_fit objects.

Extract common-environment effects

Description

Usage

common_env_effects(object, ...)

Arguments

Details

common_env_effects() returns the predicted common-environment effects of the opt-in, experimental two-effect model (target = "two_effect").

Value

Common-environment effect results for two-effect hsquared_fit objects.

Extract the common-environment variance ratio (c2)

Description

Usage

common_env_proportion(object, ...)

Arguments

Details

common_env_proportion() returns the estimated common-environment variance ratio c2 = sigma_c2 / (sigma_a2 + sigma_c2 + sigma_e2) from the opt-in, experimental two-effect model (target = "two_effect" with a common_env() term).

Falconer fence: c2 is a variance ratio (the proportion of phenotypic variance from the shared common-environment effect), not a heritability. In the same fit, heritability() reports narrow-sense h2 = sigma_a2 / (sigma_a2 + sigma_c2 + sigma_e2) within this two-effect model. The returned data frame carries this note as an "interpretation" attribute. For an interval on c2, see common_env_proportion_interval().

Value

A one-row data frame with term and estimate, plus an "interpretation" attribute, for two-effect hsquared_fit objects.

Extract an experimental common-environment / maternal variance-ratio interval

Description

Usage

common_env_proportion_interval(object, ...)

maternal_proportion_interval(object, ...)

Arguments

Details

common_env_proportion_interval() and maternal_proportion_interval() return an experimental large-sample (logit delta-method) confidence interval for the second variance ratio (c2 / m2, ratio2) of the opt-in two-effect model, available only when the fit contains it. On the same two-effect fit, heritability_interval() returns the matching interval for the direct heritability (h2, ratio1).

This mirrors the engine row V3-TWOEFFECT-REML: the interval is the asymptotic delta-method CI built from the two-effect REML observed information (the finite-difference Hessian of the two-effect REML log-likelihood at the optimum). It is asymptotic, delta-method, REML only, and NOT coverage-calibrated – on small samples the REML surface is flat and the interval is unreliable (the parametric bootstrap is the only finite-sample- aware path). No calibrated coverage is claimed.

Boundary honesty: when the ratio's variance component sits on the boundary (sigma -> 0) it is flagged (boundary = TRUE) and lower/upper are NA, never a spuriously tight CI. lower_clamped/upper_clamped flag when a bound reaches the numerical ⁠(0, 1)⁠ rails.

Falconer fence: c2 / m2 is a variance ratio, not a heritability (see common_env_proportion() / maternal_proportion()).

Value

A one-row data frame with estimate (the ratio), lower, upper, level, se, lower_clamped, upper_clamped, and boundary, plus an "interpretation" attribute, for two-effect hsquared_fit objects that contain it.

Extract experimental multivariate covariance standard errors

Description

Usage

covariance_standard_errors(object, ...)

Arguments

Details

covariance_standard_errors() returns experimental large-sample (delta-method) standard errors for the multivariate genetic/residual covariance and correlation matrices and per-trait h^2, for an opt-in unstructured multivariate fit, when the engine returned them.

Heavy caveats (engine multivariate REML validation row, partial): the strict per-seed recovery gate is still a non-pass (7/12 unstructured seeds in the updated study), but a 12-seed bias/MCSE study (twin HSquared.jl#78) shows no detectable bias – all six covariance parameters have ⁠|bias| <= 2*MCSE⁠ (largest 0.84*MCSE), a low-power non-rejection consistent with an unbiased estimator (not a proof), with EBV accuracy ~ 0.90 in both traits; a cold-start replication (#79) reaches the same optimum on all 12 seeds, so it is not a warm-start artifact. The per-seed gate failures reflect sampling variance of the estimated G at this design, not a detected bias. These SEs remain asymptotic, REML-only, unstructured-only (the engine refuses structured / factor-analytic fits, whose loadings are rotation-nonidentified), omitted at a flat/boundary optimum, not coverage-calibrated, with no external comparator, and not a validated capability.

The returned genetic_correlation SE is on the correlation (r) scale. A naive symmetric interval built as ⁠estimate +/- 1.96 * se⁠ on that scale can leave ⁠(-1, 1)⁠ for strong correlations. The engine instead builds its own genetic-correlation interval on the Fisher-z scale (z = atanh(r), se_z = se_r / (1 - r^2), endpoints tanh(z -+ q * se_z)), which always lands inside ⁠(-1, 1)⁠. That range-respecting endpoint is a separate property from calibration: the Fisher-z transform does not make these SEs coverage-calibrated (see above).

Value

A named list of standard-error matrices genetic_covariance, residual_covariance, genetic_correlation, residual_correlation, and a per-trait heritability SE vector, for hsquared_fit objects that contain them.

Likelihood-ratio test for genetic covariance structure

Description

Usage

covariance_structure_lrt(constrained, full, ...)

Arguments

Details

covariance_structure_lrt(constrained, full) is an experimental nested likelihood-ratio test comparing two multivariate fits on the same data: a constrained genetic structure (currently genetic_structure = "diagonal") against the full "unstructured" fit. The statistic is 2 * (logLik(full) - logLik(constrained)) on df = n_genetic_params(full) - n_genetic_params(constrained).

For diagonal-vs-unstructured the null (off-diagonal genetic covariances = 0) is interior, so the chi^2 reference is exact (boundary = FALSE). Structures whose null lies on a rank/PSD boundary (low-rank / factor-analytic) would need a chi^2-mixture correction and are gated out of the R bridge for now.

It mirrors the engine multivariate REML validation row (partial): asymptotic, REML-only, dense validation-scale, with the multivariate recovery calibration not yet passed – a reported test, not a validated one. Both fits must be on the same response, fixed effects, and pedigree.

Value

A one-row data frame with statistic, df, pvalue, boundary, and the constrained / full genetic-structure labels.

Inspect hsquared data-container status

Description

Usage

data_status(data)

Arguments

Details

data_status() gives a direct user-facing view of the checks stored in an hs_data() object. It reports component presence, ID overlap diagnostics, pedigree diagnostics, expression-feature diagnostics, genotype-column diagnostics, and marker-map/genotype-marker alignment diagnostics when those inputs are supplied. When annotation_id is supplied, it reports expression-feature annotation coverage diagnostics. When environment_id is supplied, it also reports environment metadata coverage diagnostics. It does not fit models, build genomic relationship matrices, add eQTL terms, or add environment-effect terms.

Value

An "hs_data_status" object.

Examples

ped <- data.frame(
  id = c("a", "b", "c", "d"),
  sire = c(NA, NA, "a", "a"),
  dam = c(NA, NA, "b", "c")
)
dat <- data.frame(
  y = c(1, 2, 3),
  sex = c("f", "m", "f"),
  id = c("a", "c", "d")
)
bundle <- hs_data(phenotypes = dat, pedigree = ped)
data_status(bundle)

Direct-maternal correlated model extractors

Description

Usage

direct_heritability(object, ...)

direct_variance(object, ...)

partner_variance(object, ...)

direct_maternal_covariance(object, ...)

total_heritability(object, ...)

Arguments

Details

These extractors summarize an opt-in, covered at validation scale direct-maternal correlated model (target = "direct_maternal"), fitted with animal(1 | id, pedigree = ped) + maternal_genetic(1 | dam). Validation-scale coverage concerns point estimates within the documented scope; this is not a production-scale or calibrated-interval claim. The model estimates a 2x2 genetic covariance matrix G_dm between the direct additive and maternal additive effects, plus a residual variance.

Willham fence: the direct-maternal model (Willham 1963, 1972) distinguishes the direct narrow-sense heritability ⁠h2_d = sigma_ad / sigma^2_P⁠ from the maternal variance ratio ⁠m2 = sigma_am / sigma^2_P⁠, the Willham total (selection-response) heritability ⁠h2_T = (sigma_ad + 1.5*sigma_dm + 0.5*sigma_am) / sigma^2_P⁠, and the genetic covariance sigma_dm. ⁠sigma^2_P = sigma_ad + sigma_am + sigma_dm + sigma_e2 = Var(y_i)⁠ for a non-inbred base (coefficient 1 on sigma_dm because 2 * A[i,dam] = 2 * (1/2) = 1). h2 is denominator-dependent under maternal effects; compare (co)variance components, not h2 values, across software (ASReml/BLUPF90/WOMBAT/sommer/MCMCglmm all leave sigma^2_P to the user). A negative genetic correlation r_am is real and biologically expected in many livestock traits; it reflects an antagonistic direct-maternal relationship and does not by itself indicate a model failure.

Note: the 2x2 G_dm formulation is due to Willham (1963, 1972), not Falconer (1965). Falconer's single-m model fixes ⁠r_am = +/-1⁠ and is a special case (Bijma 2011).

direct_heritability() returns the direct narrow-sense heritability with an interpretation attribute. total_heritability() returns Willham's selection-response heritability h2_T. direct_variance(), partner_variance(), and direct_maternal_covariance() return the raw variance-component scalars. genetic_correlation() on a direct-maternal fit returns the between-effect correlation r_am.

Value

• direct_heritability(): a one-row data frame with term and estimate, plus an "interpretation" attribute.

• total_heritability(): a one-row data frame with term = "total_willham" and estimate, plus an "interpretation" attribute.

• direct_variance(), partner_variance(), direct_maternal_covariance(): single numerics (the raw variance component).

See Also

variance_components(), genetic_correlation(), maternal_effects(), total_heritability()

Examples

# Synthetic result object: demonstrate extraction, not model fitting.
# Phenotypic variance is 0.30 + 0.15 - 0.10 + 0.55 = 0.90.
fit_dm <- structure(list(
  spec = list(target = "direct_maternal"),
  result = list(
    direct_variance = 0.30, partner_variance = 0.15, covariance = -0.10,
    variance_components = data.frame(
      component = c("direct", "maternal", "covariance", "residual"),
      estimate = c(0.30, 0.15, -0.10, 0.55)),
    heritability = data.frame(term = "direct", estimate = 0.30 / 0.90),
    genetic_correlation = data.frame(
      term_1 = "direct", term_2 = "maternal",
      estimate = -0.10 / sqrt(0.30 * 0.15)),
    maternal_effects = data.frame(id = "dam1", value = 0.1),
    converged = TRUE)
), class = "hsquared_fit")

direct_heritability(fit_dm)    # h2_d (direct narrow-sense)
total_heritability(fit_dm)     # h2_T (Willham selection-response)
genetic_correlation(fit_dm)    # r_am (may be negative)
direct_variance(fit_dm)        # sigma_ad
partner_variance(fit_dm)       # sigma_am
direct_maternal_covariance(fit_dm)  # sigma_dm
variance_components(fit_dm)    # all four components
maternal_effects(fit_dm)       # maternal EBVs

Reserved factor-analytic and G-matrix extractors

Description

Usage

genetic_loadings(object, ...)

specific_variance(object, ...)

latent_breeding_values(object, ...)

Arguments

Details

These extractor names are reserved for future factor-analytic G-matrix results. The current package can report invariant covariance and correlation matrices from multivariate fits, but it does not yet expose interpreted loadings, uniqueness/specific variance, or latent breeding values. Loading columns (genetic_loadings()) and latent breeding values (latent_breeding_values()) are rotation-nonunique until a rotation or constraint policy is validated. Future hsquared_fit methods reserve effect and rotation controls, but these controls currently error rather than implying that loading axes are interpretable. The rotation-invariant genetic eigenstructure and evolvability geometry are available now via eigen_G() and the g_matrix_geometry family.

specific_variance() (Psi, the factor-analytic specific/unique variances) is a different case: for ⁠G = Lambda Lambda' + Psi⁠, rotating ⁠Lambda -> Lambda Q⁠ with ⁠QQ' = I⁠ leaves ⁠Lambda Q Q' Lambda' = Lambda Lambda'⁠ unchanged, so Psi is rotation-invariant and identified, not rotation-nonunique. The engine payload for structured multivariate fits carries Psi as genetic_uniqueness, explicitly marked identified, alongside the excluded, rotation-nonidentified loadings. This extractor still errors on the R surface, but for a different reason: genetic_structure = "factor_analytic" (and "lowrank") are planned on the R surface and are not yet activated on the R-to-Julia bridge, not because Psi is unidentified.

Value

These reserved extractors currently error for hsquared_fit objects.

Inspect fitted-model diagnostics

Description

Usage

fit_diagnostics(object, ...)

Arguments

Details

fit_diagnostics() returns a compact diagnostics table for an hsquared_fit object. It is an inspection helper over the current result payload: it does not refit the model, rerun validation checks, or promote an experimental bridge target to production support.

Two rows share the word "boundary" but report unrelated things (hsquared#230):

• at_boundary / at_boundary_condition derive from the fitted estimates: a variance component at or near zero, or an inadmissible negative estimate (hs_fit_boundary_class()). Emitted for the Gaussian animal-model / genomic / single-step routes that report one primary variance-share layout; the ⁠_condition⁠ row is dropped (not NA) when the fit is interior, since there is nothing extra to say.

• search_boundary / search_boundary_condition mirror the non-Gaussian bridge's optimizer state: result$boundary (documented as fit$result$boundary in ?hs_control) reports whether HSquared.jl's single-variance Brent search stopped on the rail of its log-scale bracket (⁠log(initial$sigma_a2) +/- 6⁠, HSquared.jl#327) or, under restart_check = TRUE, whether the two-start gap fence fired; either way the estimate is a function of the search, not the data. Only fits whose result carries a boundary field (the non-Gaussian route) gain this pair; a Gaussian engine = "fit" fit gains neither row (not even NA). Unlike at_boundary_condition, search_boundary_condition is always emitted alongside a present search_boundary flag, reading "interior" rather than disappearing when the flag is FALSE – the search-bracket rail is a binary optimizer state with no benign-vs- inadmissible split to only sometimes report. In practice a search_boundary = TRUE fit is refused before it reaches an hsquared_fit object at all, so a returned fit's search_boundary is always FALSE; TRUE is reachable only on a hand-built result object.

Value

A data frame with metric and value columns and class "hs_fit_diagnostics".

Extract fixed effects

Description

Usage

fixef(object, ...)

Arguments

Details

fixef() is part of the planned v0.1 fitted-object contract for hsquared_fit objects.

Value

Fixed-effect estimates for hsquared_fit objects.

Inspect formula grammar status

Description

Usage

formula_status()

Details

formula_status() reports which pieces of the planned hsquared() formula language are parsed today, reserved as syntax markers, or still roadmap-only. It is a status table, not a model-fitting helper.

The printed header is derived from the rows being printed, so it cannot lag the table. Covered-versus-experimental fences live in ⁠$current_behavior⁠ and are not dumped into the default print.

Value

A data frame of formula grammar records with class "hs_formula_status".

Examples

formula_status()

G-matrix geometry and evolvability

Description

Usage

evolvability(object, direction, ...)

respondability(object, direction, ...)

conditional_evolvability(object, direction, ...)

autonomy(object, direction, ...)

mean_evolvability(object, ...)

g_max(object, ...)

variance_along_gradient(object, direction, normalize = TRUE, ...)

eigen_G(object, ...)

Arguments

Details

These extractors summarise the genetic variance-covariance matrix G of an multivariate hsquared_fit (a cbind() response, which fits on the default path) through its rotation-invariant geometry, following Hansen & Houle (2008). They are defined on G itself (not on factor loadings), so they are well defined for any multivariate fit – unstructured, diagonal, or (when bridged) low-rank / factor-analytic – and do not depend on a loading rotation convention.

• eigen_G() returns the genetic eigenstructure: values (the variance along each genetic principal axis, descending) and vectors (the genetic principal components, sign-canonicalised so the largest-magnitude element of each is positive).

• g_max() returns the leading genetic axis: its eigenvalue and eigenvector (the direction of maximum evolvability).

• mean_evolvability() is the average evolvability over random selection directions, tr(G) / t.

• evolvability(fit, direction) is ⁠e(beta) = beta'G beta⁠ (unit beta): the additive genetic variance available to directional selection along direction.

• variance_along_gradient(fit, direction, normalize) is ⁠beta'G beta⁠ for a unit direction (normalize = TRUE, equal to evolvability()) or for the raw direction (normalize = FALSE).

• respondability(fit, direction) is ⁠||G beta||⁠: the length of the response.

• conditional_evolvability(fit, direction) is ⁠1 / (beta' G^-1 beta)⁠: evolvability when all other directions are held under stabilising selection (requires a positive-definite G).

• autonomy(fit, direction) is conditional_evolvability / evolvability in ⁠[0, 1]⁠ (requires a positive-definite G).

These are computed in R from the fitted genetic_covariance(fit) and match the engine's evolvability.jl definitions (verified by a skip-guarded live parity test). They are descriptive geometry of the estimated G; they carry the same experimental, REML-only, not-coverage-calibrated status as the multivariate fit itself and report no standard errors.

Value

A scalar for the directional metrics; a list for eigen_G() (values, vectors) and g_max() (eigenvalue, eigenvector).

References

Hansen, T. F., & Houle, D. (2008). Measuring and comparing evolvability and constraint in multivariate characters. Journal of Evolutionary Biology, 21(5), 1201-1219.

Extract a supplied metafounder Gamma matrix

Description

Usage

gamma_matrix(object, ...)

Arguments

Details

gamma_matrix() returns the supplied metafounder Gamma matrix carried by an experimental metafounder or H^Gamma single-step hsquared_fit object. It is provenance for the fitted relationship, not an estimated parameter.

Value

A numeric matrix with metafounder labels when the fitted object contains a supplied Gamma payload.

Genomic and QTL formula markers

Description

Usage

genomic(formula, G = NULL, Ginv = NULL, markers = NULL, ...)

single_step(
  formula,
  H = NULL,
  Hinv = NULL,
  pedigree = NULL,
  markers = NULL,
  group = NULL,
  Gamma = NULL,
  tau = 1,
  omega = 1,
  blend_weight = 0,
  ridge = 0,
  ...
)

markers(M, model = c("random", "fixed", "scan"), ...)

marker_scan(M, map = NULL, ...)

qtl_scan(position, genotype_probs = NULL, ...)

Arguments

Details

These functions provide readable formula vocabulary for genomic, single-step, marker-effect, GWAS, and QTL/eQTL models. Called directly they are inert (they return NULL); they take meaning only inside an hsquared() formula. genomic() fits a narrow Gaussian REML random-intercept model only through the explicit ⁠engine = "julia", target = "genomic"⁠ path. single_step() also remains opt-in and experimental. The reported genomic coefficient is sigma_g2 / (sigma_g2 + sigma_e2) on the declared relationship scale, not generally an average marginal phenotypic-variance fraction or pedigree/population heritability. Accepted genomic forms are genomic(1 | id, Ginv = Ginv) or genomic(1 | id, markers = M) (GREML, or SNP-BLUP via target = "snp_blup"), and single_step(1 | id, Hinv = Hinv) (a precomputed inverse) or single_step(1 | id, pedigree = ped, markers = M) (the engine constructs H^-1 from the pedigree + genotyped-subset markers via target = "single_step_construct"). Adding group plus supplied Gamma fits a supplied-Gamma, validation-scale H^Gamma path through target = "metafounder_single_step"; Gamma is supplied, not estimated. When data is an hs_data() container that bundles a pedigree and genotypes, single_step(1 | id) resolves both from the bundle (the animal(1 | id) precedent), so neither ⁠pedigree =⁠ nor ⁠markers =⁠ is required; explicit arguments override the bundle. The remaining markers (markers(), marker_scan(), qtl_scan()) are still inert syntax reservations that the parser rejects with a planned-not-implemented message.

Value

NULL, invisibly. Calls are interpreted by hsquared() when they appear inside model formulas.

Post-fit relatedness-corrected marker scan (GWAS)

Description

Usage

gwas(
  object,
  markers,
  marker_ids = NULL,
  method = c("mixed", "single", "loco"),
  marker_groups = NULL,
  genome_wide = FALSE,
  n_permutations = 1000L,
  seed = 1L,
  ...
)

Arguments

Details

gwas() runs a dense, supplied-variance, relatedness-corrected mixed-model (GLS) Wald marker scan on a fitted Gaussian animal model, reusing the fit's estimated variance components ⁠(sigma^2_a, sigma^2_e)⁠ and pedigree relationship. It is an experimental, validation-scale screen that surfaces the Julia-owned HSquared.mixed_model_marker_scan().

The default p_value/bonferroni_p/bh_qvalue columns are NOT genome-wide calibrated. They are marker-by-marker Wald (nominal) p-values plus deterministic Bonferroni and Benjamini-Hochberg adjustments over the supplied marker set only. Do not report genome-wide significance from those columns.

Genome-wide calibration is available via genome_wide = TRUE (which requires method = "single", the validated fixed-effect scope). It adds a genome_wide_p column: the exact per-dataset add-one permutation p-value (the marker is declared genome-wide significant when genome_wide_p <= alpha, alpha = 0.05). For each analysis the permutation null is rebuilt from the analysed phenotype (y permuted conditional on X, re-scanned n_permutations times), and the genome-wide p is ⁠(1 + #{null max >= observed})/(n_permutations + 1)⁠. This is the HSquared.jl genome_wide_marker_scan engine call, whose family-wise type-I control is validated at validation and production scale (the per-dataset add-one REBUILD gate; the anti-conservative (1-alpha) quantile rule is NOT used). The result carries a calibration attribute (method permutation_addone, empirical_type1 = NA because the per-dataset rule's validity is by construction and externally validated, named in validation_reference). SCOPE: fixed-effect / intercept-only; the relatedness-corrected mixed-model genome-wide null is NOT yet calibrated, so genome_wide = TRUE is rejected for method = "mixed"/"loco".

A leave-one-group-out (LOCO) scan is available with method = "loco" and a marker_groups argument: when a marker is tested, the genomic relationship correction is built from the markers not in that marker's group (e.g. its chromosome), so the marker's own signal does not leak into the background relationship. The LOCO relationship is genomic (VanRaden) while the reused variance components are pedigree-estimated – a scale mismatch that keeps this validation-scale (see method).

Value

An hs_gwas data frame with one row per marker: marker, effect, se, z, chisq, p_value, bonferroni_p, bh_qvalue, lod (and, when genome_wide = TRUE, genome_wide_p), carrying a scan_method attribute (and, when genome_wide = TRUE, a calibration attribute). Its print() method restates the uncalibrated-significance caveat for the nominal columns (and, for method = "single", the absence of any relatedness correction; for method = "loco", the genomic-vs-pedigree scale mismatch).

Extract heritability estimates

Description

Usage

heritability(object, ...)

Arguments

Details

heritability() is part of the planned v0.1 fitted-object contract. It works for hsquared_fit objects that contain a Julia result.

Falconer fence for the opt-in two-effect model (target = "two_effect"): the reported number is the narrow-sense direct heritability h2 = sigma_a2 / (sigma_a2 + sigma_2 + sigma_e2) within that model (the additive-genetic variance is divided by the total phenotypic variance, which now includes the second component sigma_2). The second component's variance ratio (common-environment c2 or maternal m2) is not a heritability and is returned separately by common_env_proportion() / maternal_proportion().

Willham fence for the opt-in direct-maternal correlated model (target = "direct_maternal"): heritability() returns the labelled triple – direct h2_d, maternal m2, Willham total h2_T, and r_am – as a data frame (Willham 1963, 1972). ⁠sigma^2_P = sigma_ad + sigma_am + sigma_dm + sigma_e2 = Var(y_i)⁠ (coefficient 1 on sigma_dm). A warning is issued because h2 is denominator-dependent under maternal effects and h2_T < h2_d is expected when r_am < 0. Use direct_heritability() or total_heritability() for targeted accessors without the warning.

Reaction-norm fence for the opt-in random-regression model (target = "random_regression"): heritability is a trajectory h2(t) over the covariate, not a scalar, so heritability() errors on such a fit and names the implemented accessor. Use rr_heritability() for the h2(t) curve.

A non-converged fit still returns the engine number so you can inspect it, but warns: that number is not an estimate. A near-zero value from a failed fit is not evidence that heritability is zero. Use fit_diagnostics() before reading any number. logLik() already refuses a non-converged fit; heritability() keeps the value and shouts instead.

Value

Heritability results for hsquared_fit objects.

Extract an experimental heritability confidence interval

Description

Usage

heritability_interval(object, ...)

Arguments

Details

heritability_interval() returns an experimental large-sample confidence interval for h^2. It is available only when an hsquared_fit object contains the interval field, which the default Gaussian animal-model fit (engine = "fit") populates from the engine's HSquared.heritability_interval() when a local Julia engine is present and the estimate is interior to ⁠(0, 1)⁠. On the opt-in two-effect fit it returns the direct-heritability ratio interval (ratio1), and on the opt-in multi-effect fit (target = "multi_effect", K >= 3 blocks) it returns the ANIMAL block's ratio interval (the animal additive variance over the total phenotypic variance); the other blocks' variance-ratio intervals are surfaced separately in fit$result$variance_ratio_intervals.

The interval leg is a REML-only, asymptotic (logit delta-method or profile) approximation returned by the engine. The 2000-rep C1 coverage confirm (HSquared.jl DRAC job 47925485) places the univariate pedigree h² interval at the directional-conservative claim level under doc-34 §4: delta over-covers (worst Ĉ 0.969), profile is in-band (worst Ĉ 0.950), and bootstrap mildly under-covers versus nominal 0.95 at the governing interior cell (worst Ĉ 0.924 at h²=0.5, still at least 0.90). Those figures were measured on one interpretable small design (q = 120), at interior ⁠h² in {0.3, 0.5, 0.7}⁠, at the 0.95 level, in DRAC job 47925485; unbalanced designs were not in the confirm grid, so these numbers say nothing about coverage under imbalance. That is not coverage-calibrated at nominal and not "never under-covers." Profile and bootstrap h² interval legs were measured in the same confirm but are not separately surfaced by this extractor; only the engine-returned interval is shown. It is reported as a point estimate plus bounds, not a validated (coverage-calibrated) capability, and remains unreliable near the ⁠h² → 0⁠ boundary. Genomic fits still error (scale-labelled interval not validated). The underlying estimators V3-TWOEFFECT-REML / V3-NEFFECT-REML are covered, but this interval is not.

Value

A one-row data frame with estimate, lower, upper, level, se (NA for the profile method), and method, for hsquared_fit objects that contain it.

Create hsquared control options

Description

hs_control() records execution and storage controls for hsquared model calls. The default engine = "fit" fits the validated v0.1 Gaussian animal model through the HSquared.jl engine (average-information REML); engine = "validate" parses and validates the contract without fitting, returning the validated spec invisibly; and engine = "julia" exposes advanced engine controls.

Usage

hs_control(
  engine = c("fit", "validate", "julia"),
  backend = c("auto", "cpu", "threads", "cuda", "amdgpu", "metal", "oneapi"),
  accelerator = c("auto", "none", "gpu", "cuda", "amdgpu", "metal", "oneapi"),
  precision = c("float64", "float32"),
  save = c("minimal", "full", "tiny"),
  engine_control = list()
)

Arguments

Value

An object of class "hs_control".

Create an hsquared data container

Description

Usage

hs_data(
  phenotypes,
  pedigree = NULL,
  genotypes = NULL,
  markers = NULL,
  expression = NULL,
  annotation = NULL,
  annotation_id = NULL,
  environment = NULL,
  environment_id = NULL,
  id = "id"
)

Arguments

Details

hs_data() collects phenotype, pedigree, genotype, marker, expression, annotation, and environment inputs into one checked container. It is a lightweight data-contract object for future genomic, QTL/eQTL, and multi-omics workflows. It does not fit models. The v0.1 parser can use an hs_data object directly as data, reading model variables from phenotypes and making named components such as pedigree available to formula terms.

summary(hs_data(...)) reports ID overlap diagnostics, pedigree diagnostics, and, when expression, genotype, or marker components are supplied, expression-feature, genotype-column, marker-map, and genotype-marker alignment diagnostics. When annotation_id is supplied, it reports expression-feature annotation coverage diagnostics. When environment_id is supplied, it also reports environment metadata coverage diagnostics.

Value

An hs_data object.

Examples

ped <- data.frame(
  id = c("a", "b", "c", "d"),
  sire = c(NA, NA, "a", "a"),
  dam = c(NA, NA, "b", "c")
)
dat <- data.frame(
  y = c(1, 2, 3),
  sex = c("f", "m", "f"),
  id = c("a", "c", "d")
)
bundle <- hs_data(phenotypes = dat, pedigree = ped)
bundle
summary(bundle)

Forest plot of a known-truth recovery study

Description

Usage

hs_recovery_forest(data)

Arguments

Details

Visualizes a bias +/- 2MCSE recovery study (e.g. data-raw/multivariate-recovery-study.R): each target's bias with its +/- 2MCSE interval and a zero-bias reference line. Targets whose interval covers zero show "no detectable bias".

Value

A ggplot object.

Fit a quantitative-genetic model

Description

The default path fits a univariate Gaussian animal model by REML: y ~ sex + age + animal(1 | id, pedigree = ped). Fitting needs a local Julia and HSquared.jl. Preview the same call with control = hs_control(engine = "validate").

Usage

hsquared(
  formula,
  data,
  family = stats::gaussian(),
  REML = TRUE,
  control = hs_control(),
  ...
)

Arguments

Details

What you may report is listed on Can I fit and report this? – not in validation_status(). Other routes live in formula_status() and on that limits page.

Value

A "hsquared_fit" object from a successful default-path fit. Missing Julia produces an install-guidance error. When engine = "validate", the validated model spec is returned invisibly as a named list (not the classed object from model_spec()).

ggplot2 visualizations for hsquared results

Description

Usage

## S3 method for class 'hsquared_fit'
autoplot(
  object,
  type = c("variance", "breeding_values", "g_matrix", "g_geometry", "reaction_norm",
    "rr_eigenfunctions", "rr_surface"),
  ...
)

## S3 method for class 'hs_gwas'
autoplot(object, type = c("manhattan", "qq"), ...)

Arguments

Details

autoplot() methods render the quantitative-genetic results an hsquared_fit (or a gwas() scan) carries as ggplot2 objects, in the style of the brms/bayesplot ecosystem and consistent with the sister packages drmTMB/gllvmTMB. They are uncertainty-first: where the fit carries the experimental standard errors / reliabilities, the figures show them (clearly labelled experimental and asymptotic).

Available types for autoplot.hsquared_fit():

• "variance" (default) – a horizontal forest of the variance components and per-trait h^2, each with approximate 95% intervals (⁠+/- 1.96 * SE⁠) when the fit carries the experimental standard errors.

• "breeding_values" – a sorted caterpillar of the estimated breeding values, with ⁠+/- 1.96 * sqrt(PEV)⁠ bands when prediction error variances are available, faceted by trait for multivariate fits.

• "g_matrix" – a rotation-invariant genetic-correlation heatmap of the estimated G for multivariate fits (correlations are invariant to the factor rotation; raw loadings are never plotted – the ratified cross-lane convention). Off-diagonal cells involving a low-h^2 trait are flagged as imprecise (threshold low_h2, default 0.1).

• "g_geometry" – a scree of the rotation-invariant genetic eigenstructure (eigenvalues = variance per genetic axis, with percent variance explained) for multivariate fits. Axis directions / loadings are never drawn (rotation-arbitrary; span-ambiguous under repeated eigenvalues).

• "reaction_norm" – for random-regression fits, the genetic-variance and heritability trajectories across the covariate (faceted). The heritability trajectory carries the same caveat as rr_heritability(): with a homogeneous residual and no permanent-environment term it can overstate h^2(t) for repeated-records designs.

• "rr_eigenfunctions" – for random-regression fits, the rotation-invariant eigenfunctions psi_j(t) of K_g as covariate functions (faceted by axis, labelled by percent genetic variance). Signs are arbitrary and the curves are span-ambiguous under repeated eigenvalues (do not over-read).

• "rr_surface" – for random-regression fits, the genetic covariance surface ⁠S(s, t) = phi(s)' K_g phi(t)⁠ over the covariate grid as a heatmap (pass correlation = TRUE for the genetic-correlation surface). Supplied-K_g descriptive.

autoplot() on a gwas() scan (hs_gwas) draws type = "manhattan" (default) or type = "qq" (observed vs expected -log10(p) with a y = x null and the genomic-inflation lambda_GC as a diagnostic). Both carry the EXPERIMENTAL, NOT-genome-wide-calibrated caveat. The Julia lane has a fixed-panel calibration smoke harness, but R significance thresholds are not activated.

The figure helpers are deliberately modular (each takes a tidy data frame and returns a ggplot) so they can be factored into a shared visualization package later.

Value

A ggplot object.

Block unsupported likelihood-inference helpers

Description

These methods intentionally fail with explicit scope messages. Validated confint() / vcov() / profile() / anova() surfaces are not implemented. Experimental interval extractors (heritability_interval(), variance_component_standard_errors()) exist when the engine returned those fields; they are not coverage-calibrated. Point estimates remain variance_components(), heritability(), and fit_diagnostics().

Usage

## S3 method for class 'hsquared_fit'
confint(object, parm, level = 0.95, ...)

## S3 method for class 'hsquared_fit'
vcov(object, ...)

## S3 method for class 'hsquared_fit'
profile(fitted, ...)

## S3 method for class 'hsquared_fit'
anova(object, ...)

Arguments

Value

These functions always error.

Extract planned marker, QTL, GWAS, and eQTL results

Description

Usage

marker_effects(object, ...)

marker_variance_explained(object, ...)

qtl_table(object, ...)

gwas_table(object, ...)

eqtl_table(object, ...)

lod_scores(object, ...)

Arguments

Details

These extractor names cover genomic, QTL, GWAS, and eQTL fitted results. They return values only when an hsquared_fit object contains the corresponding result field. marker_effects() and marker_variance_explained() are populated by the opt-in SNP-BLUP path (target = "snp_blup"). The variance-explained table is a descriptive fitted-marker share, computed as effect squared times centered marker variance and normalized across markers; it is not a marker-scan p-value, QTL statistic, or causal decomposition under linkage disequilibrium. The gwas_table() and lod_scores() methods for an already-computed hs_gwas object expose the current uncalibrated marker-scan table and marker-level LOD scores. Fit-level QTL/GWAS/eQTL result tables and map-annotated / calibrated scan outputs remain reserved for future results.

Value

The requested marker or scan result for hsquared_fit objects that contain the corresponding field.

Extract maternal genetic effects

Description

Usage

maternal_effects(object, ...)

Arguments

Details

maternal_effects() returns the predicted maternal genetic effects of the opt-in, experimental maternal two-effect model (target = "two_effect" with a maternal_genetic() term).

Value

Maternal genetic effect results for maternal hsquared_fit objects.

Extract the maternal variance ratio (m2)

Description

Usage

maternal_proportion(object, ...)

Arguments

Details

maternal_proportion() returns the estimated maternal variance ratio m2 = sigma_m2 / (sigma_a2 + sigma_m2 + sigma_e2) from the opt-in, experimental maternal two-effect model (target = "two_effect" with a maternal_genetic() term).

Falconer fence: m2 is a variance ratio (the proportion of phenotypic variance from the maternal (dam) genetic effect), not a heritability. In the same fit, heritability() reports the narrow-sense direct h2 = sigma_a2 / (sigma_a2 + sigma_m2 + sigma_e2) within this two-effect model; the direct and maternal genetic effects are modelled as uncorrelated here (no direct-maternal genetic covariance). The returned data frame carries this note as an "interpretation" attribute. For an interval on m2, see maternal_proportion_interval().

Value

A one-row data frame with term and estimate, plus an "interpretation" attribute, for maternal two-effect hsquared_fit objects.

Reserved metafounder effect extractor

Description

Usage

metafounder_effects(object, ...)

Arguments

Details

metafounder_effects() reserves the extractor name for future fitted metafounder solutions. Current experimental metafounder and H^Gamma fits expose supplied Gamma and group-assignment provenance via gamma_matrix() and metafounder_groups(), but the engine does not yet return explicit combined-system metafounder effects for extraction.

Value

This extractor currently errors for all objects.

Extract supplied metafounder group assignments

Description

Usage

metafounder_groups(object, ...)

Arguments

Details

metafounder_groups() returns the ID-keyed metafounder group assignments carried by an experimental metafounder or H^Gamma single-step hsquared_fit object. These assignments are provenance for the supplied relationship, not estimated grouping parameters.

Value

A data frame with columns id, metafounder_group, and is_metafounder. Rows are animals, keyed by id; is_metafounder flags animals whose unknown-parent slot carries a supplied group label, not metafounders themselves. The metafounders are the pseudo-populations that are the rows and columns of gamma_matrix().

Inspect a parsed hsquared model specification

Description

Usage

model_spec(formula, data, family = stats::gaussian(), REML = TRUE, ...)

Arguments

Details

model_spec() validates the narrow v0.1 animal-model grammar and returns the R-side model specification that would be used to build the Julia bridge payload. It is a preview and diagnostics helper: it does not fit a model.

Value

An "hs_model_spec" object.

Extract multivariate covariance and correlation matrices

Description

Usage

genetic_covariance(object, ...)

G_matrix(object, ...)

residual_covariance(object, ...)

R_matrix(object, ...)

genetic_correlation(object, ...)

residual_correlation(object, ...)

Arguments

Details

These extractors return the genetic (G) and residual (R) covariance or correlation matrices from multivariate hsquared_fit objects (a cbind() response, which fits on the default path). G_matrix() is an applied-workflow alias for genetic_covariance(), and R_matrix() is an alias for residual_covariance(). Use them after checking fit_diagnostics() because likelihood-based summaries are intentionally blocked when a multivariate fit has not converged.

Value

A numeric matrix for hsquared_fit objects that contain the requested multivariate result field.

Examples

# Synthetic result object: demonstrate extraction, not model fitting.
# No Julia installation is needed for this example.
G <- matrix(c(0.4, 0.1, 0.1, 0.3), 2,
            dimnames = list(c("weight", "length"), c("weight", "length")))
R <- diag(c(0.6, 0.7))
dimnames(R) <- dimnames(G)
fit_mv <- structure(list(
  spec = list(target = "multivariate"),
  result = list(genetic_covariance = G, residual_covariance = R,
                genetic_correlation = cov2cor(G),
                residual_correlation = cov2cor(R),
                heritability = data.frame(trait = colnames(G),
                                          estimate = diag(G) / diag(G + R)),
                converged = TRUE)
), class = "hsquared_fit")

genetic_covariance(fit_mv)
G_matrix(fit_mv)
residual_covariance(fit_mv)
R_matrix(fit_mv)
genetic_correlation(fit_mv)
residual_correlation(fit_mv)
heritability(fit_mv)

Extract permanent-environment effects

Description

Usage

permanent_effects(object, ...)

Arguments

Details

permanent_effects() returns the predicted permanent-environment effects of the opt-in, experimental repeatability model.

Separating sigma_a2 from sigma_pe2 needs repeated records per individual; see ?hs_control (target = "repeatability") for that identifiability requirement.

Value

Permanent-environment effect results for repeatability hsquared_fit objects.

Diagnostic plots for a fitted animal model

Description

A base-graphics diagnostic plot for hsquared_fit objects. Two panels are available:

Usage

## S3 method for class 'hsquared_fit'
plot(x, type = c("variance", "residuals"), ...)

Arguments

Details

• type = "variance" (default) plots the estimated variance components as points. When the fit carries the experimental variance-component standard errors (see variance_component_standard_errors()), it adds approximate ⁠+/- 1.96 * SE⁠ whiskers and labels the panel experimental; those intervals are asymptotic, REML-only, and not coverage-calibrated. The governing coverage study found under-coverage for the sigma^2_a delta/Wald interval (worst interior coverage 0.897 at h^2 = 0.5), so no conservative direction is claimed for the plotted variance-component SE.

• type = "residuals" plots residuals against fitted values (with a zero reference line), when the fit carries fitted values and a response.

No plotting dependency is added; this uses base graphics only.

Value

x, invisibly.

Extract prediction error variances

Description

Usage

prediction_error_variance(object, ...)

Arguments

Details

prediction_error_variance() is part of the planned v0.1 fitted-object contract. It returns values only when an hsquared_fit object contains a Julia result field for prediction error variances.

The engine defines PEV as the animal-block diagonal of the inverse of the precision-scaled Henderson mixed-model-equation coefficient matrix, so diag(C^-1) is directly the prediction error variance in variance units, with no separate sigma_e2 factor to apply.

Value

Prediction error variances for hsquared_fit objects.

Quantitative-genetic formula markers

Description

Usage

permanent(formula, ...)

common_env(formula, ...)

maternal_genetic(formula, pedigree = NULL, ...)

maternal_env(formula, ...)

paternal_genetic(formula, pedigree = NULL, ...)

paternal_env(formula, ...)

cytoplasmic(formula, ...)

imprinting(formula, pedigree = NULL, parent = c("maternal", "paternal"), ...)

dominance(formula, pedigree = NULL, D = NULL, Dinv = NULL, ...)

epistasis(formula, pedigree = NULL, E = NULL, Einv = NULL, ...)

relmat(formula, K = NULL, Kinv = NULL, Q = NULL, ...)

precision(formula, Q = NULL, ...)

group(formula, ...)

unknown_parent_group(formula, ...)

metafounder(formula, pedigree = NULL, group = NULL, Gamma = NULL, ...)

inbreeding(formula, ...)

Arguments

Details

These functions provide readable formula vocabulary for standard quantitative-genetic, parental, inheritance, and custom-kernel models. Called directly they are inert (they return NULL); they take meaning only inside an hsquared() formula. permanent(), common_env(), and maternal_genetic() now fit through an opt-in, experimental engine path (engine = "julia", REML-only, not the default, mirroring a partial validation gate) as the second random effect alongside animal(). metafounder() now fits through an opt-in, supplied-variance, experimental engine path (engine = "julia", target = "metafounder"); Gamma and the variance components are supplied, not estimated. relmat() and precision() now fit through an opt-in, EXPERIMENTAL supplied-relationship primary engine path (engine = "julia", REML-only, not the default, mirroring a partial validation gate): relmat(1 | id, K = K) supplies a dense symmetric positive-definite relationship/covariance matrix K (the parser marshals the inverse Kinv = solve(K); relmat() also accepts a directly supplied inverse via ⁠Kinv =⁠/⁠Q =⁠, no solve), and precision(1 | id, Q = Q) supplies the precision (inverse) directly (⁠Kinv =⁠ also accepted); both fit through the SAME supplied-relationship-inverse path as genomic()'s Ginv, and the supplied matrix is provenance, not an estimate. The remaining markers (paternal/maternal-environment, dominance, epistasis, cytoplasmic, imprinting, genetic groups / unknown-parent-groups, and inbreeding) are still inert syntax reservations that the parser rejects with a planned-not-implemented message.

Some markers use generic names (e.g. group(), inbreeding()). They are formula-only tokens detected by call head and are not meant to be called directly, so attaching hsquared alongside a package that exports a same-named function (e.g. pedigreemm::inbreeding()) may print a masking message; this is expected and harmless.

Value

NULL, invisibly. Calls are interpreted by hsquared() when they appear inside model formulas.

Random-regression (reaction-norm) extractors

Description

Usage

rr_covariance(object, ...)

random_coefficients(object, ...)

rr_genetic_variance(object, at = NULL, n = 25L, ...)

rr_heritability(object, at = NULL, n = 25L, ...)

rr_correlation(object, at = NULL, n = 25L, ...)

rr_eigenfunctions(object, at = NULL, n = 25L, ...)

Arguments

Details

These extractors summarize an opt-in, experimental random-regression (reaction-norm) fit (target = "random_regression"), fitted with an animal(rr(covariate, order = k) | id, pedigree = ped) term. The model estimates a ⁠k x k⁠ genetic covariance matrix K_g among an animal's normalized-Legendre random-regression coefficients plus a single homogeneous residual variance.

• rr_covariance() returns the estimated ⁠k x k⁠ coefficient genetic covariance matrix K_g.

• random_coefficients() returns the per-animal predicted Legendre coefficients (long format: id, coefficient, value).

• rr_genetic_variance() returns the additive genetic variance trajectory ⁠v_g(t) = phi(t)' K_g phi(t)⁠ across covariate points.

• rr_heritability() returns the heritability trajectory h^2(t) = v_g(t) / (v_g(t) + sigma_e^2). Because the residual is homogeneous and there is no permanent-environment term yet, this can OVERSTATE h^2(t) for repeated-records designs (test-day, growth curves).

• rr_correlation() returns the genetic correlation surface among the covariate points.

• rr_eigenfunctions() returns the eigen-decomposition of the coefficient covariance K_g evaluated as functions of the covariate: the eigenvalues (additive genetic variance carried by each axis), the proportion of genetic variance explained, the sign-canonicalized eigen-coefficients, and the eigenfunctions ⁠psi_j(t) = phi(t)' v_j⁠ over the covariate grid. These are rotation-invariant functionals of K_g (mirroring the engine's rr_eigenfunctions); no raw, rotation-arbitrary loadings are returned.

The trajectories are computed in R from the estimated K_g and the recorded covariate standardization range; at is supplied on the ORIGINAL covariate scale (defaulting to a grid over the fitted range) and re-standardized to ⁠[-1, 1]⁠ internally, matching the Julia engine's basis convention.

Value

rr_covariance() returns a numeric matrix; random_coefficients() returns a data frame; rr_genetic_variance() and rr_heritability() return data frames with covariate and value columns. rr_correlation() returns a numeric matrix. rr_eigenfunctions() returns a list containing eigenvalues, variance proportions, eigen-coefficients and evaluated eigenfunctions.

Out-of-range at

at must lie inside the fitted covariate range (object$result$random_regression$lower to ⁠$upper⁠); a value outside that range errors with a message naming the fitted range, instead of silently returning the nearest endpoint's value. This matches the Julia engine (HSquared.legendre_basis), which refuses a standardized covariate outside ⁠[-1, 1]⁠ rather than clamping it (#213).

Examples

# Synthetic result object: demonstrate extraction, not model fitting.
# The supplied covariance describes two Legendre coefficients.
fit_rr <- structure(list(
  spec = list(target = "random_regression"),
  result = list(coefficient_covariance = diag(c(1, 0.2)),
                residual_variance = 0.5,
                random_coefficients = data.frame(
                  id = rep("animal1", 2),
                  coefficient = c("legendre0", "legendre1"),
                  value = c(0.1, 0.2)),
                random_regression = list(
                  covariate = "age", lower = 1, upper = 5, order = 2),
                converged = TRUE)
), class = "hsquared_fit")

rr_covariance(fit_rr)
random_coefficients(fit_rr)
rr_genetic_variance(fit_rr)
rr_heritability(fit_rr)
rr_correlation(fit_rr, at = c(1, 3, 5))
rr_eigenfunctions(fit_rr)

Extract random effects

Description

Usage

ranef(object, ...)

Arguments

Details

ranef() is part of the planned v0.1 fitted-object contract for hsquared_fit objects.

Value

Random-effect estimates for hsquared_fit objects.

Objects exported from other packages

Description

These objects are imported from other packages. Follow the links below to see their documentation.

ggplot2

autoplot()

Extract reliability and accuracy estimates

Description

Usage

reliability(object, ...)

accuracy(object, ...)

Arguments

Details

reliability() is part of the planned v0.1 fitted-object contract. It returns values only when an hsquared_fit object contains a Julia result field for reliability estimates. accuracy() returns the square root of reliability for hsquared_fit objects.

The engine defines reliability as 1 - PEV_i / (sigma_a2 * A_ii), using the dense relationship matrix A = inv(Ainv) implied by the supplied precision. For a genomic spec (Ainv = Ginv) this A_ii is 'diag(inv(Ginv)) = diag(G)

• ridge⁠(the regularized genomic self-relationship, often not equal to 1), so the genomic ridge perturbs the reported reliability and accuracy. Reliability values are **not clipped**;⁠accuracy()' rejects a reliability value below 0 or above 1 rather than silently clipping it.

Value

Reliability estimates for hsquared_fit objects.

Extract repeatability estimates

Description

Usage

repeatability(object, ...)

Arguments

Details

repeatability() reports the repeatability R = (Va + Vpe) / Vp of the opt-in, experimental repeatability (permanent-environment) model. It works for hsquared_fit objects fitted with engine_control = list(target = "repeatability").

Separating sigma_a2 from sigma_pe2 needs repeated records per individual; see ?hs_control (target = "repeatability") for that identifiability requirement.

Value

Repeatability results for repeatability hsquared_fit objects.

Extract an experimental repeatability confidence interval

Description

Usage

repeatability_interval(object, ...)

Arguments

Details

repeatability_interval() returns an experimental large-sample (logit delta-method) confidence interval for the repeatability coefficient t = (Va + Vpe) / Vp of the opt-in repeatability (permanent-environment) model, available only when the fit contains it.

It mirrors the engine row V3-REPEAT-REML (partial): the engine's repeatability REML estimator and this interval are engine-internal self-consistency tested (recovery of t and interval bracketing / range / level-nesting / point-estimate match on seeded fixtures), but there is no external comparator, no h^2 interval, and no deep-pedigree validation. It is asymptotic, REML-only, unreliable at small n or near the (0, 1) boundary (where the engine throws and the field is omitted), and not a validated capability.

Value

A one-row data frame with estimate (the repeatability t), lower, upper, level, and se, for hsquared_fit objects that contain it.

Response-scale prediction helpers

Description

predict(), fitted(), and residuals() are part of the planned v0.1 fitted-object contract for univariate hsquared_fit objects. They are univariate-only: the multivariate model (a cbind() response) fits multiple traits jointly and is intentionally out of v0.1 response-scale scope, so these methods stop with a scope message pointing to breeding_values(), genetic_covariance(), and residual_covariance().

Usage

## S3 method for class 'hsquared_fit'
predict(object, ...)

## S3 method for class 'hsquared_fit'
fitted(object, ...)

## S3 method for class 'hsquared_fit'
residuals(object, ...)

Arguments

Value

Response-scale predictions, fitted values, or residuals for univariate hsquared_fit objects.

hsquared ggplot2 theme

Description

Usage

theme_hsquared(base_size = 12)

Arguments

Details

A light, publication-oriented ggplot2 theme shared by the hsquared autoplot() figures. Exported so users can restyle or extend the figures.

Value

A ggplot2 theme object.

Inspect validation evidence status

Description

Usage

validation_status()

Details

validation_status() reports the current validation atoms and planned comparator lanes for hsquared. It is a developer evidence table, not the user-facing list of what you may report. For that list, read Can I fit and report this?. This function does not run validation checks, fit models, or promote any capability to working status.

Value

A data frame of validation status records with class "hs_validation_status".

This table is not the full list of covered routes

The rows here are validation atoms, grouped by evidence chain rather than by public modelling route, so the covered row count is not the number of covered routes. Two capabilities that count as publicly covered have no row of their own yet:

• the opt-in reaction-norm model, k = 2: animal(rr(x, order = 2) | id, pedigree = ped) with engine_control = list(target = "random_regression");

• the opt-in correlated direct-maternal 2x2 G model: animal(1 | id, pedigree = ped) + maternal_genetic(1 | dam) with engine_control = list(target = "direct_maternal").

Both are covered at validation scale, opt-in, REML-only, and dense (n <= 1000 or so). Their scope and reporting limits are recorded in the repository capability ledger (docs/design/capability-status.md) and, for readers, in Can I fit and report this?. Read this table together with that enumeration: an absent row is not an absent capability, and this table alone does not tell you everything you may report.

Bridge production fences (0.9)

The R-to-Julia bridge is live at validation scale for covered default and opt-in targets, but not production-complete. Infrastructure capabilities (⁠R-to-Julia bridge payload⁠, ⁠opt-in experimental Julia engine⁠) stay partial even when individual targets are covered. Targets split into: default live (engine = "fit"), opt-in live (engine_control$target), PATH_ONLY smoke (C1-ext intervals — not a fit bridge), and blocked planned grammar (for example genetic_structure = "factor_analytic"). payload_v2 routing is limited to direct_maternal and multi_effect. See docs/design/45-bridge-production-fences-DRAFT.md.

Capability identifiers versus labels

capability is a stable identifier. Dated evidence records – comparator runs under ⁠docs/dev-log/comparator-runs/⁠ and check-log entries – cite these strings verbatim to name the row they report against, so an identifier is not rewritten once evidence points at it. Look rows up by capability.

capability_label is the current wording, and is what reader-facing surfaces print. The two differ only where an identifier's wording has been overtaken by a later change; today that is the multivariate row, whose cbind() route became default rather than opt-in.

Examples

validation_status()

Extract experimental variance-component and heritability standard errors

Description

Usage

variance_component_standard_errors(object, ...)

heritability_standard_error(object, ...)

Arguments

Details

variance_component_standard_errors() and heritability_standard_error() return experimental large-sample (delta-method) standard errors derived from the REML average-information matrix. They are available only when an hsquared_fit object contains them; the default Gaussian animal-model fit populates them from the engine when a local Julia engine is present and the AI matrix is invertible. When the engine could not compute them, hsquared() warns once at fit time and keeps the engine's reason in attr(fit, "bridge_errors"), a character vector named by the engine function that failed (HSquared.jl#351); the field is then absent, not NA.

These mirror the engine row V1-HERIT-CI (partial): asymptotic, REML-only, and unreliable at small n or near a variance-component boundary (where the AI matrix is ill-conditioned and the fields are omitted). The 2000-rep C1 coverage confirm (job 47925485) measured the additive-variance delta/Wald interval implied by these SEs to under-cover (0.897 at nominal 0.95, h²=0.5), placing it at experimental-only: the SE is a point-estimate reference only, not a calibrated and not a conservative interval. (A profile additive-variance interval is directional-conservative by the same run — including the h²=0.3 cell 0.963 — but is not surfaced by this delta-method extractor. Additive-variance bootstrap worst-cell Ĉ 0.918 is the same DC mild-under branch as h² bootstrap.) Not coverage-calibrated, not a validated capability.

Value

variance_component_standard_errors() returns a data frame with component and se. heritability_standard_error() returns a one-row data frame with term and se, matching heritability() so the two merge on term (hsquared#236). Only for hsquared_fit objects that contain them. The fit object may still store heritability_se as a bare numeric internally; this extractor is the public shape.

Extract variance components

Description

Usage

variance_components(object, ...)

Arguments

Details

variance_components() is part of the planned v0.1 fitted-object contract. It works for hsquared_fit objects that contain a Julia result.

Value

Variance component results for hsquared_fit objects.