Skip to contents

Animal-model formula marker

Source: R/animal.R
animal.Rd

[Experimental]

Usage

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

Arguments

formula

A random-effect expression. The v0.1 parser accepts 1 | id. The opt-in random-regression form is rr(x, order = k) | id; see rr().

pedigree

A pedigree data frame with individual, sire, and dam columns. Optional only when the enclosing hsquared() or model_spec() call uses data = hs_data(..., pedigree = ...). Unknown parents may be NA, 0, "0", or "". Other codes (".", "-9", "UNK", "*", the string "NA", -1) are treated as real IDs and error if they are not also rows. Every known parent must appear as an id row; phantom parents are 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. In v0.1 this route rejects any row with the same known sire and dam (selfing); no argument here reaches the engine's existing allow_selfing flag. The current workaround for a selfing or hermaphroditic pedigree is relmat(1 | id, K = A) with a hand-built or AGHmatrix-built relationship matrix A (an experimental route; see ?hs_control).

...

Reserved for future syntax such as cov =.

Value

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

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. The opt-in random-regression form is animal(rr(x, order = k) | id, pedigree = ped); see rr() for order (number of Legendre coefficients; order = 2 is intercept + slope) and the [-1, 1] mapping. Fitting happens in hsquared(), not here: calling animal() on its own is a syntax marker and returns NULL.

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")
)
#> Validated the v0.1 animal-model contract; no model was fitted (`engine = "validate"`). Julia fit target: `HSquared.fit_animal_model(y, X, Z, Ainv; method = :REML)`. The validated spec is returned invisibly as a named list; assign it to inspect the parsed contract (e.g. `spec$bridge$target`).
model_spec(
  weight ~ sex + animal(1 | id, pedigree = ped),
  data = dat
)
#> <hs_model_spec>
#>   response: weight
#>   family: gaussian(identity)
#>   method: REML
#>   observations: 4; fixed columns: 2; animal IDs: 4
#>   animal term: animal(1 | id, pedigree = ped); group: id
#>   random design: sparse Z with 4 nonzero entries
#>   Ainv: build_in_julia
#>   fitting: not run by `model_spec()`