Inheritance systems roadmap • hsquared

This article shows how hsquared is expected to grow beyond the ordinary diploid additive animal model.

The core design rule is simple:

inheritance system -> relationship matrix K or precision matrix Q
model term         -> design matrix Z
engine target      -> variance component and BLUP/EBV calculation

That rule keeps animal breeding, plant breeding, evolutionary ecology, and genomic work on one modelling ladder. It also keeps the public API from turning into a bag of special cases.

What works today

The default fit is the v0.1 additive animal model:

ped <- data.frame(
  id   = c("s1", "s2", "d1", "d2", "o1", "o2", "o3", "o4"),
  sire = c(NA, NA, NA, NA, "s1", "s1", "s2", "s2"),
  dam  = c(NA, NA, NA, NA, "d1", "d2", "d1", "d2")
)
pheno <- data.frame(
  id  = ped$id,
  sex = c("m", "m", "f", "f", "m", "f", "m", "f"),
  age = c(4, 5, 4, 5, 1, 1, 1, 1),
  y   = c(72, 75, 65, 68, 31, 33, 35, 30)
)

fit <- hsquared(
  y ~ sex + age + animal(1 | id, pedigree = ped),
  data = pheno,
  family = gaussian(),
  REML = TRUE
)

The following standard extensions fit only through opt-in experimental engine targets:

# permanent environment / repeatability
y ~ animal(1 | id, pedigree = ped) + permanent(1 | id)

# IID common environment
y ~ animal(1 | id, pedigree = ped) + common_env(1 | litter)

# independent direct additive + maternal genetic effect
y ~ animal(1 | id, pedigree = ped) + maternal_genetic(1 | dam)

Those are not unusual-inheritance models. They are useful stepping stones: extra random effects, pedigree-linked second effects, and clear extractor boundaries.

Planned named effects

Some future effects already have reserved formula markers:

y ~ animal(1 | id, pedigree = ped) +
  maternal_env(1 | dam)

y ~ animal(1 | id, pedigree = ped) +
  paternal_genetic(1 | sire, pedigree = ped) +
  paternal_env(1 | sire)

y ~ animal(1 | id, pedigree = ped) +
  cytoplasmic(1 | maternal_line)

y ~ animal(1 | id, pedigree = ped) +
  imprinting(1 | id, pedigree = ped, parent = "maternal")

y ~ animal(1 | id, pedigree = ped) +
  dominance(1 | id, pedigree = ped) +
  epistasis(1 | id, pedigree = ped)

These examples do not fit yet. Today the parser rejects them as planned, not implemented.

The supplied-relationship primaries relmat(1 | id, K = K_custom) and precision(1 | id, Q = Q_custom) are the exception: they do fit through the opt-in engine = "julia" path (the same supplied-relationship-inverse engine as genomic()). They are REML-only, experimental, and not covered — supply a dense positive-definite K (or its inverse Q), and the matrix is carried as provenance, not estimated.

Supplying your own relationship matrix

Three spellings take a user matrix, and they are not interchangeable:

Term	Argument	What you pass	Inverse
`relmat(1 \\| id, K = K)`	`K`	the relationship / covariance itself	R inverts it (`Kinv = solve(K)`)
`precision(1 \\| id, Q = Q)`	`Q` (or `Kinv`)	the precision	used as given
`genomic(1 \\| id, Ginv = Ginv)`	`Ginv`	the genomic relationship inverse	used as given

K, Q, and Ginv must be square, symmetric, finite, and positive definite, with unique non-missing row and column names equal to the ids. Rows may be in any order; they are matched by name, not by position. The fitted coefficient from relmat() is Va / (Va + Ve) on the supplied kernel, not a pedigree heritability unless K is the pedigree A.

A monozygotic (MZ) twin, clone, or repeated-genotype kernel is singular by construction: identical individuals have identical rows, so the smallest eigenvalue is zero and the parser stops with `K` must be positive definite. This package does not accept that singular K and does not invent a new decomposition. The AE workaround is a small ridge, K + 1e-6 * diag(n). Variance Va and the REML log-likelihood match the unridged kernel; the printed residual is E - 1e-6 * Va, so printed h2 = Va / (Va + Ve) is slightly inflated unless you add 1e-6 * Va back to the denominator. A larger eps moves residual and printed h2 further while leaving Va unchanged. An ACE twin model (relmat() plus a pair / common-environment term) is not available from R; collapsing each MZ pair to one genotype is the design alternative. This is not an ACE estimator.

# One MZ pair (K = 1) and one DZ pair (K = 0.5). Unrelated between pairs.
twin_dat <- data.frame(
  id = c("mz1", "mz2", "dz1", "dz2"),
  zyg = c("MZ", "MZ", "DZ", "DZ"),
  y = c(1.10, 1.02, 0.40, 0.55)
)
K <- diag(nrow(twin_dat))
rownames(K) <- colnames(K) <- twin_dat$id
K["mz1", "mz2"] <- K["mz2", "mz1"] <- 1
K["dz1", "dz2"] <- K["dz2", "dz1"] <- 0.5
# The raw MZ kernel is singular. Ridge it before fitting.
K_fit <- K + 1e-6 * diag(nrow(K))
fit_k <- hsquared(
  y ~ relmat(1 | id, K = K_fit),
  data = twin_dat,
  control = hs_control(
    engine = "julia",
    engine_control = list(target = "relmat")
  )
)
variance_components(fit_k)
heritability(fit_k) # Va / (Va + Ve) on the ridged kernel

The Julia engine can fit a twin ACE model outside R: HSquared.fit_multi_effect() takes the kernel and the pair effect as two independent random effects, each a (Z, Ainv) pair. This route is experimental. On the full twinbmi data (n = 8,542) it returned converged = false (HSquared.jl#426), so treat it as unvalidated for real twin data. In R, relmat() cannot yet be combined with (1 | pair) or common_env() (hsquared#296), so there is no supported R route for ACE from a custom kernel. The snippet below is a smoke test of the call on small simulated data (150 pairs). It is run from Julia, not from R, and is shown unevaluated here. The kernel is ridged for the same reason as above, so the printed E again carries the 1e-6 * Va shift.

using HSquared, LinearAlgebra, SparseArrays, Random

Random.seed!(1)
npair = 150
zyg = [i <= 75 ? :MZ : :DZ for i in 1:npair]
n = 2 * npair
pair = repeat(1:npair, inner = 2)
K = Matrix{Float64}(I, n, n)
for p in 1:npair
    i = 2p - 1
    K[i, i+1] = K[i+1, i] = zyg[p] == :MZ ? 1.0 : 0.5
end
Kridge = K + 1e-6 * I
Zpair = sparse(1:n, pair, 1.0, n, npair)
a = cholesky(Symmetric(Kridge)).L * randn(n) * sqrt(0.6)
c = Zpair * (randn(npair) * sqrt(0.2))
y = 1.0 .+ a .+ c .+ randn(n) * sqrt(0.4)
X = ones(n, 1)
Ainv_K = sparse(inv(Symmetric(Kridge)))
effects = [(sparse(1.0I, n, n), Ainv_K), (Zpair, sparse(1.0I, npair, npair))]
fit = HSquared.fit_multi_effect(y, X, effects; method = :auto, verbose = false)
println(fit.variance_components)
println(fit.converged, " ", fit.loglik)

fit.variance_components.sigmas holds A and C in the order the effects were listed, and fit.variance_components.sigma_e2 holds E. The data above are simulated, so the numbers are only a smoke test of the call, not a twin analysis. In a real MZ/DZ study K has 1 between MZ co-twins and 0.5 between DZ co-twins, and Zpair indexes the shared pair environment.

Planned inheritance systems

Other systems are roadmap examples rather than exported syntax:

# planned future shapes, not implemented
animal(1 | id, pedigree = ped, inheritance = selfing(rate = s))
animal(1 | id, pedigree = ped, inheritance = clonal())
animal(1 | id, pedigree = ped, inheritance = haplodiploid())
animal(1 | id, pedigree = ped, inheritance = polyploid(ploidy = 4))
animal(1 | id, pedigree = ped, inheritance = cytoplasmic())

These names are useful design placeholders, but the package does not currently export selfing(), clonal(), haplodiploid(), polyploid(), or an inheritance = grammar for animal().

Why kernels first

Most inheritance systems become mixed-model terms only after the relationship or precision structure is clear.

System	Kernel question	Public risk if rushed
selfing / partial selfing	How does inbreeding change additive covariance through the pedigree?	A standard diploid A matrix can give the wrong expectation.
clonal / asexual	Which records share genotypes, ramets, environments, or clone lines?	Clone identity can be confounded with common environment.
haplodiploid	How do sex-specific inheritance and relatedness enter A or Q?	Parent-offspring and sib covariances are not the ordinary diploid ones.
polyploid	What ploidy, double reduction, and dosage assumptions define relatedness?	Diploid shortcuts can misstate plant-breeding relationships.
cytoplasmic	What is the maternal-line transmission structure?	It can be mistaken for maternal environment.
imprinting	Which parent-of-origin effect is modelled?	Maternal and paternal allele-origin effects can be conflated.
dominance / epistasis	Which non-additive relationship matrix is used?	Additive h2 can be over-interpreted as total genetic variance.
custom K / Q	Who constructed the matrix, and on what scale?	Scale, ordering, and positive-definiteness mistakes become silent model errors.

The engine should prefer sparse precision matrices when possible. If a system cannot supply a well-defined K or Q, the formula should stay planned.

Validation gates

An inheritance system should not move beyond planned until all of these are true:

mathematical definition of the estimand and relationship/precision matrix;
tiny deterministic pedigree or genotype example;
ID ordering and missing-parent behaviour pinned by tests;
positive-definite or positive-semidefinite checks where relevant;
comparator check against an established package or a published worked example when one exists;
recovery simulation for the variance component and predicted effects;
extractor and summary wording that separates additive, non-additive, maternal, paternal, cytoplasmic, and environmental interpretations.

For dominance and epistasis, nadiv is a natural comparison target. For diploid/autopolyploid A/G/H relationship construction, AGHmatrix is a natural comparison target. For breeding-program simulations, XSim.jl is a natural future partner.

User-facing rule

The common path should stay easy:

y ~ fixed + animal(1 | id, pedigree = ped)

Specialist inheritance should be explicit:

# planned
y ~ fixed +
  animal(1 | id, pedigree = ped) +
  cytoplasmic(1 | maternal_line)

The user should never need to guess whether “maternal” means maternal genetic, maternal environmental, cytoplasmic, parent-of-origin, or dam identity. Names should stay longer until the biology is unambiguous.

Background anchors

Useful software anchors include nadiv for non-additive genetic relatedness matrices and AGHmatrix for diploid and autopolyploid A/G/H relationship construction. The package-level roadmap also tracks custom relationship and precision matrices because many plant, animal, and ecological designs will arrive as a known kernel before they arrive as a named inheritance helper.

For now, these sources motivate the roadmap. They are not evidence of current support for selfing, clonal, haplodiploid, polyploid, cytoplasmic, imprinting, dominance, epistasis, or custom-kernel models.