Formula grammar roadmap • hsquared

This page is a map of the modelling language. It separates syntax that is parsed today from syntax that is reserved for later phases.

Parsed today

The current parser contract is intentionally narrow:

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

When users pass an hs_data() bundle with a pedigree component, the same contract can be written without repeating the pedigree argument:

bundle <- hs_data(phenotypes = dat, pedigree = ped)

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

The R side validates the formula, data, response, pedigree columns, and animal IDs, builds the internal bridge payload, and (by default) fits the v0.1 univariate Gaussian animal model by REML through HSquared.jl. Terms below that say parsed are recognized; only some of those also fit, and then only on a named target =. Reserved names error as planned, not implemented.

Unknown parents in pedigree = may be NA, 0, "0", or "". Other codes (".", "-9", "UNK") are not unknown-parent markers: they are treated as real IDs and error if they are not also rows. A parent that is not also an id row is a phantom parent and is rejected. Add founder rows (data.frame(id = setdiff(c(ped$sire, ped$dam), ped$id), sire = NA, dam = NA)) or recode those slots as NA. The R pedigree = path does not expose the Julia engine’s missing_values argument.

Standard quantitative-genetic extensions

These names are reserved so future models read like the biological question:

y ~ sex + age +
  animal(1 | id, pedigree = ped) +
  permanent(1 | id) +
  common_env(1 | litter)

y ~ sex + age +
  animal(1 | id, pedigree = ped) +
  maternal_genetic(1 | dam, pedigree = ped) +
  maternal_env(1 | dam) +
  paternal_genetic(1 | sire, pedigree = ped) +
  paternal_env(1 | sire)

Current status: permanent(), common_env(), and maternal_genetic() are parsed. They fit only with a named target = (repeatability, two-effect, or direct–maternal — see Fitting models). The default hsquared() call names that target and stops. The default grammar accepts at most one of those named extra effects beside animal(). More independent intercepts use bare (1 | group) terms with target = "multi_effect". A maternal-genetic plus permanent-environment combination is not available. maternal_env(), paternal_genetic(), and paternal_env() are reserved and error as planned, not implemented.

Inheritance and relationship kernels

Later phases should treat inheritance systems as relationship or precision kernels, not as one-off special cases:

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

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

y ~ relmat(1 | id, K = K_custom)
y ~ precision(1 | id, Q = Q_custom)

Current status: relmat() and precision() fit through the opt-in experimental engine = "julia" path (target = "relmat" / "precision"). They are not covered and not the default. K and Q must be square, symmetric, and positive definite, with unique row and column names equal to the ids (any row order; names are matched). dominance(), epistasis(), cytoplasmic(), and imprinting() remain syntax reservations. See Inheritance systems for a worked custom-K example, including singular monozygotic-twin kernels.

Genomic and marker terms

The genomic vocabulary follows the same rule: readable terms first, fitting only after the R-Julia contract and validation exist.

y ~ sex + batch + genomic(1 | id, Ginv = Ginv)

y ~ sex + batch +
  genomic(1 | id, Ginv = Ginv) +
  marker_scan(M, map = marker_map)

y ~ sex + batch + single_step(1 | id, Hinv = Hinv)

Current status: markers(), marker_scan(), and qtl_scan() are planned markers and do not fit yet. genomic() and single_step() now fit through an opt-in, experimental engine path (engine = "julia", target = "genomic" / "single_step" / "snp_blup"; REML-only or supplied-variance, not the default, mirroring a partial validation gate).

Multivariate and factor-analytic grammar

The first multivariate surface is the cbind() animal-model path, which routes on the default call:

cbind(weight, height) ~ sex + animal(1 | id, pedigree = ped)

The longer multivariate ladder is planned to connect ordinary G matrices and GLLVM-style latent genetic structure:

y ~ trait + trait:sex +
  animal(trait | id, pedigree = ped, cov = us()) +
  residual(trait | unit, cov = us())

y ~ trait + trait:sex +
  animal(trait | id, pedigree = ped, cov = lowrank(K = 2)) +
  residual(trait | unit, cov = diag())

y ~ trait + trait:sex +
  animal(trait | id, pedigree = ped, cov = fa(K = 2)) +
  residual(trait | unit, cov = fa(K = 1))

Planned covariance meanings:

cov = us(): full unstructured covariance.
cov = diag(): trait-specific variances only.
cov = lowrank(K): G = Lambda Lambda'.
cov = fa(K): G = Lambda Lambda' + Psi.

Current status: the cbind() multivariate animal-model path is parsed and routed by default, and fitted experimentally — the routing is default, the claim is still partial. The long-format animal(trait | id, cov = ...) grammar and residual structured covariance grammar are planned, not parsed or fitted.

Error rule

Unsupported terms fail early. There are two live shapes.

Reserved names such as markers(), dominance(), or maternal_env() stop as planned, not implemented, then name the nearest live path:

`markers()` is planned, not implemented. The live marker path is
`genomic(1 | id, markers = )` or `gwas(fit, markers)`. Run
`formula_status()` for the live list of which terms parse and which fit.

Terms that already parse but are not on the default path – common_env(), permanent(), maternal_genetic(), rr(...) – do not use that planned stop. They paste the next hsquared() call:

`common_env()` is not on the default path.

Closest working call:

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

That route is covered for point estimates only (validation-scale,
common-environment leg). Intervals are experimental.
See: current-limits article, formula_status().

That error is part of the user interface. It names the closest working call without pretending the default path already fitted the model.

For a compact table of the same status information from R, use:

formula_status()