Genomic prediction • hsquared

Experimental 0.9.0 release · explicit GREML target · count stays 7. Narrow genomic() GREML is covered at validation scale only through its explicit Julia target. SNP-BLUP and single-step stay opt-in and experimental. This page does not flip public_covered_count and does not claim the nine-cell recovery campaign as PASS.

This article shows the current genomic prediction surface in hsquared. It is for users who already have a genomic relationship inverse, a marker matrix, or a single-step relationship inverse and want to understand which path is live today.

markers = M must be coded as the biallelic allele count 0/1/2 (or an imputed dosage in [0, 2]); a matrix in the centered -1/0/1 convention (as produced by some external genomic-prediction tools) must be recoded (add 1) before it is passed in. A 0/1 panel of homozygous inbred lines (as in BGLR::wheat) must be doubled to 0/2; hsquared does not recode it, and a named hsquared_zero_one_markers warning is raised if every entry is in {0, 1}.

The genomic GREML paths are experimental, REML-only, and Julia-backed. Narrow genomic(1 | id, markers = M) / Ginv = Q requires engine = "julia", target = "genomic"; the ordinary engine = "fit" call remains the pedigree animal-model route. SNP-BLUP and single-step still need their own explicit targets. A repaired boundary candidate resolved all 240 sealed holdouts and corrected 30 classifications with no losses, but its preregistered runtime gate failed because one cell was 5.99 times slower at p95 than the old path, above the frozen 3 times cap. That history is confirmatory debt — not a claim that the nine-cell recovery campaign passed. All fitting paths require local Julia, JuliaCall, and a sibling HSquared.jl checkout.

Three current paths

The current R surface has three distinct genomic paths.

Path	Formula term	Target	Variances
GBLUP / genomic GREML from a supplied inverse	`genomic(1 \| id, Ginv = Ginv)`	`target = "genomic"`	estimated by the Julia AI-REML path
GBLUP / genomic GREML from raw markers	`genomic(1 \| id, markers = M)`	`target = "genomic"`	estimated after the frozen engine construction
SNP-BLUP / RR-BLUP marker effects	`genomic(1 \| id, markers = M)`	`target = "snp_blup"`	supplied by the user, or REML-estimated when omitted

There is also a single-step surface:

Path	Formula term	Target	Variances
supplied-H inverse single-step	`single_step(1 \| id, Hinv = Hinv)`	`target = "single_step"`	estimated by the Julia AI-REML path
constructed-H inverse single-step	`single_step(1 \| id, pedigree = ped, markers = M)`	`target = "single_step_construct"`	estimated after the engine builds Hinv
animal-only supplied-Gamma metafounder	`metafounder(1 \| id, pedigree = ped, group = mf_group, Gamma = Gamma)`	`target = "metafounder"`	supplied by the user (`sigma_a2`, `sigma_e2`)
supplied-Gamma H^Gamma single-step	`single_step(1 \| id, pedigree = ped, markers = M, group = mf_group, Gamma = Gamma)`	`target = "metafounder_single_step"`	estimated after the engine builds H^Gamma-inverse

The constructed single-step path is experimental and validation-scale. Metafounder support is also experimental and validation-scale: R validates the supplied metafounder labels and Gamma, then calls either the Julia-owned animal-only supplied-variance metafounder_animal_model() path or fit_metafounder_single_step_reml() for H^Gamma single-step. Gamma is supplied, not estimated, and there is no metafounder-specific extractor or external comparator evidence. APY, low-rank marker solvers, weighted genomic relationships, production-scale single-step construction, and Bayesian marker models are planned. A post-fit, relatedness-corrected marker scan is now available experimentally via gwas(fit, markers) (see Post-fit marker scan below) — its p-values are not genome-wide calibrated by default; an opt-in per-dataset permutation cutoff exists for the single-marker scope (see below). QTL/eQTL scans remain planned.

Supplied Ginv

Use this path when you already have a genomic relationship inverse whose row and column names match the IDs in the phenotype data.

ids <- paste0("id", seq_len(8))
pheno <- data.frame(
  id    = ids,
  batch = rep(c("A", "B"), each = 4),
  y     = c(72, 75, 65, 68, 31, 33, 35, 30)
)

set.seed(1)
M <- matrix(sample(0:2, 8 * 6, replace = TRUE), nrow = 8, ncol = 6)
rownames(M) <- ids

p <- colMeans(M) / 2
Z <- M - matrix(2 * p, nrow = 8, ncol = 6, byrow = TRUE)
k <- 2 * sum(p * (1 - p))
G <- tcrossprod(Z) / k + diag(0.01, 8)
dimnames(G) <- list(ids, ids)
Ginv <- solve(G)

fit_g <- hsquared(
  y ~ batch + genomic(1 | id, Ginv = Ginv),
  data = pheno,
  family = gaussian(),
  REML = TRUE,
  control = hs_control(
    engine = "julia",
    engine_control = list(target = "genomic")
  )
)

variance_components(fit_g)
heritability(fit_g)
breeding_values(fit_g)
fit_diagnostics(fit_g)

Read the random effect as a genomic breeding value: $\mathbf g \sim N(0, \sigma_g^2\mathbf Q^{-1})$ , with Ginv = Q supplied on the precision scale. hsquared does not infer how Q was constructed: method, allele-frequency source, ridge, denominator, and covariance-kernel fingerprint remain unknown; the ordered-ID and supplied-precision fingerprints are known.

The parser checks that every observed ID is present in the Ginv dimnames. The current path estimates the genomic and residual variance components by the Julia AI-REML engine. heritability(fit_g) retains term = "genomic" for compatibility but labels the coefficient genomic_variance_ratio on inverse_of_supplied_precision: sigma_g2 / (sigma_g2 + sigma_e2), conditional on the precision supplied by the user. It is not generally an average marginal phenotypic-variance fraction and is not pedigree-, founder-base-, population-, or universal narrow-sense heritability. Genomic heritability_interval() and heritability_standard_error() are unavailable until a scale-labelled interval contract is separately validated. The covered narrow GREML route remains validation-scale only: the sealed boundary candidate failed its runtime gate, the nine-cell campaign did not run, and production-scale genomic validation is not complete. Those are remaining debts, not a demotion of GREML; SNP-BLUP and both single-step routes remain opt-in partial.

Marker-built G

Use this path when you have an individual-by-marker dosage matrix M, with rownames matching the phenotype IDs.

fit_marker_g <- hsquared(
  y ~ batch + genomic(1 | id, markers = M),
  data = pheno,
  family = gaussian(),
  REML = TRUE,
  control = hs_control(
    engine = "julia",
    engine_control = list(target = "genomic")
  )
)

variance_components(fit_marker_g)
heritability(fit_marker_g)
breeding_values(fit_marker_g)

The frozen construction uses sample allele frequencies $p_j=\operatorname{mean}(M_{\cdot j})/2$ , unweighted VanRaden method 1, $G=WW'/k$ , and $K_\lambda=G+0.01I$ ; the REML solver receives $Q_\lambda=K_\lambda^{-1}$ . The fit carries engine-generated SHA-256 fingerprints for ID order, marker content, K_lambda, and Q_lambda. heritability(fit_marker_g) labels its coefficient genomic_variance_ratio = sigma_g2 / (sigma_g2 + sigma_e2) on K_lambda: the genomic variance-component ratio on the declared relationship scale. It is not generally an average marginal phenotypic-variance fraction and is not pedigree-, founder-base-, population-, or universal narrow-sense heritability. The fitted residual is the residual variance of the model with kernel G + 0.01 I; it equals the residual variance on G minus 0.01 * sigma_g2. Genomic variance and the REML log-likelihood are unchanged relative to an unridged G fit. To compare residual with GCTA / sommer / rrBLUP on G, add 0.01 * sigma_g2 back. The marker-built ridge is fixed; genomic() has no ridge argument.

This is convenient for small and validation-scale examples. It is not yet a production genotype pipeline: PLINK/VCF readers, imputation hooks, scaling/blending choices, APY, and large on-disk marker workflows remain future work.

SNP-BLUP marker effects

Use this path for per-marker effects. Supply variance_components when you want the effects at variance components you already trust (for example from a prior genomic GREML fit) — a supplied-variance solve. Omit variance_components to have hsquared() estimate sigma_g2/sigma_e2 by REML from the markers (fit_snp_blup_reml); the example below shows the supplied-variance form.

fit_snp <- hsquared(
  y ~ batch + genomic(1 | id, markers = M),
  data = pheno,
  family = gaussian(),
  REML = TRUE,
  control = hs_control(
    engine = "julia",
    engine_control = list(
      target = "snp_blup",
      variance_components = c(sigma_g2 = 1.0, sigma_e2 = 2.0)
    )
  )
)

marker_effects(fit_snp)
marker_variance_explained(fit_snp)
breeding_values(fit_snp)
heritability(fit_snp)
fit_diagnostics(fit_snp)

marker_effects(fit_snp) returns one effect per marker column. The genomic breeding values are the marker-derived predictions for each genotyped individual. marker_variance_explained(fit_snp) returns a descriptive table with each marker’s fitted contribution, computed as effect squared times the centered marker variance and normalized across markers. Treat it as a summary of the fitted SNP-BLUP, not as a marker-scan p-value, QTL signal, or causal decomposition under linkage disequilibrium. When you supply sigma_g2/sigma_e2 the fitted object records variance_components_source = "supplied"; if you omit them, hsquared() estimates the variances by REML and records variance_components_source = "estimated_snp_blup_reml".

Supplied Hinv single-step

Use this path when you already have a single-step inverse relationship matrix Hinv with dimnames matching the phenotype IDs.

fit_h <- hsquared(
  y ~ batch + single_step(1 | id, Hinv = Hinv),
  data = pheno,
  family = gaussian(),
  REML = TRUE,
  control = hs_control(
    engine = "julia",
    engine_control = list(target = "single_step")
  )
)

variance_components(fit_h)
heritability(fit_h)
breeding_values(fit_h)

Constructed Hinv single-step

Use this path when you have a pedigree for all animals and a marker matrix for the genotyped subset. Marker rows must be named by animal ID. Phenotyped animals do not all need marker rows; ungenotyped animals are the point of single-step.

fit_h_construct <- hsquared(
  y ~ batch + single_step(1 | id, pedigree = ped, markers = M),
  data = pheno,
  family = gaussian(),
  REML = TRUE,
  control = hs_control(
    engine = "julia",
    engine_control = list(target = "single_step_construct")
  )
)

variance_components(fit_h_construct)
heritability(fit_h_construct)
breeding_values(fit_h_construct)

The engine builds the pedigree inverse, dense pedigree relationship, genomic relationship for the genotyped subset, and the single-step precision H^-1. Construction knobs (tau, omega, blend_weight, ridge) exist for validation work but are not comparator-validated defaults.

Output checklist

For the current genomic paths, start with:

fit_diagnostics(fit_g)
variance_components(fit_g)
heritability(fit_g)
breeding_values(fit_g)

Do not call heritability_interval() or heritability_standard_error() on a genomic fit; those public surfaces are deliberately fenced. Inspect the relationship-scale and provenance columns returned by heritability().

For SNP-BLUP, also use:

marker_effects(fit_snp)
marker_variance_explained(fit_snp)

Post-fit marker scan (experimental)

gwas(fit, markers) 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 and pedigree relationship. Markers are one row per animal (in the fit’s pedigree order), one column per marker.

fit <- hsquared(y ~ x + animal(1 | id, pedigree = ped), data = dat)
scan <- gwas(fit, markers) # markers: animals x markers
scan # effect, se, z, chisq, p_value, bonferroni_p, bh_qvalue, lod
gwas_table(scan)
lod_scores(scan)

The default p-values are not genome-wide calibrated: they are nominal Wald p-values plus Bonferroni/Benjamini-Hochberg over the supplied markers only, with no realistic-LD production calibration and no external scan comparator. A scoped, opt-in genome-wide cutoff is available for the single-marker scan: gwas(fit, markers, method = "single", genome_wide = TRUE) adds a genome_wide_p column from an exact per-dataset add-one permutation rule (fixed-effect / intercept-only scope; rejected for method = "mixed"/"loco"; experimental, validation-scale, not the default). The relatedness-corrected mixed and leave-one-group-out loco (marker_groups =) scans stay uncalibrated. The Julia lane also has a fixed-marker-panel type-I smoke harness (HSquared.jl PR #134). Do not report genome-wide significance from the nominal p-value columns.

gwas_table(scan) and lod_scores(scan) are thin views of the already-computed hs_gwas result. Fit-level tabular gwas_table() / qtl_table() / eqtl_table() output names stay reserved for a planned map-annotated scan API:

qtl_table(fit_g)    # reserved (planned map-annotated API)
gwas_table(fit_g)   # reserved
eqtl_table(fit_g)   # reserved

Current evidence boundary

The current R genomic and single-step rows have different status: narrow genomic GREML is covered at validation scale, while supplied and constructed single-step routes remain opt-in partial. This page does not change the current public_covered_count of 7.

Covered in the R lane:

parser and payload checks for supplied Ginv, marker matrices, supplied Hinv, constructed single-step inputs, animal-only supplied-Gamma metafounder inputs, supplied-Gamma H^Gamma inputs, and SNP-BLUP variance inputs;
skip-guarded live bridge tests for genomic GREML, marker-built G, supplied Hinv, constructed single-step, animal-only supplied-Gamma metafounder, supplied-Gamma H^Gamma, and SNP-BLUP when Julia and the sibling engine are available;
an exact live identity gate showing the explicit marker route and the exact supplied-Q_lambda route return the same fitted quantities, plus an independent base-R reconstruction of the frozen symbolic inputs;
extractor checks for variance components, the scale-labelled genomic variance ratio, genomic breeding values, marker effects, and provenance diagnostics;
a Julia-native supplied-variance GBLUP/SNP-BLUP target fixture (HSquared.jl PR #140, commit 008ea4d) for future external comparator runs. This is a target fixture, not comparator evidence.

Still planned:

production relationship-matrix construction and scaling/blending controls;
component and allocation profiling on the already-open discovery data, output-equivalent performance work, a new preregistration, and a fresh untouched holdout block; the 240 opened holdout seeds cannot adjudicate a revised candidate, and the nine-cell recovery campaign remains held;
robustness studies for LD, population structure, imputation, base-frequency misspecification, real panels, and production genotype data; the sealed holdout used only independent HWE/no-LD markers;
production-scale and externally comparator-validated single-step construction;
external validation for animal-only metafounder and single-step H^Gamma;
APY and low-rank large-marker workflows;
PLINK/VCF/BCF readers and on-disk marker storage;
broader external construction comparators beyond the independent base-R calculation and the fresh hash-pinned exact-Q_lambda blupf90+ run;
marker scans, QTL, GWAS, and eQTL;
uncertainty intervals and reliability for production genomic predictions.

Background anchors

The current syntax follows standard genomic prediction ideas: marker-based prediction of total genetic value, genomic relationship matrices for GBLUP, and single-step blending of pedigree and genomic information. Useful anchors are Meuwissen, Hayes and Goddard (2001), VanRaden (2008), and Aguilar, Misztal, Johnson, Legarra, Tsuruta and Lawlor (2010).

For hsquared, those papers motivate the roadmap. The package exposes one narrow explicit genomic GREML target (covered at validation scale; experimental 0.9.0 release; count stays 7) while keeping SNP-BLUP, single-step, recovery-campaign PASS, interval, and broad production genomic claims behind explicit validation gates.