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

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.