Skip to contents

A worked animal-model analysis (gryphon)

Source: vignettes/articles/gryphon-worked-example.Rmd
gryphon-worked-example.Rmd

Intervals: asymptotic / experimental · not coverage-calibrated. This is a teaching walk of the default univariate animal-model path on gryphon. Whiskers, when shown, are not coverage-calibrated 95% CIs for publication.

This article walks one univariate animal model end to end on a real teaching dataset — fit, heritability with its experimental uncertainty, breeding values and accuracy, and a diagnostic plot — so you can see the whole hsquared extractor surface in one place. It complements the broad Fitting models tour and the sommer / pedigreemm benchmark (which checks the same fit against external packages).

Two honesty notes first:

  • Fitting needs a local engine. The default hsquared() fit marshals to a local Julia + HSquared.jl + JuliaCall. The code below is shown but not executed at build time; the heritability value quoted is the published/recovered estimate (see the benchmark article for the executed agreement).
  • gryphon uses a supplied relationship matrix. The open Wilson et al. 2010 gryphon pedigree has no ancestral loops (parents precede offspring; no self-parents). This article still uses the supplied A_gryphon matrix as the signed-off recovery path. Most of your own datasets will use the ordinary animal(1 | id, pedigree = ped) path shown in Fitting models.

The data and question

The gryphon birth-weight dataset (enhancer::DT_gryphon / A_gryphon, from Wilson et al. 2010, An ecologist’s guide to the animal model) is the canonical teaching example. The question is the classic one: how heritable is birth weight? — i.e. what fraction of phenotypic variance is additive genetic, h² = σ²a / (σ²a + σ²e).

library(hsquared)
utils::data("DT_gryphon", package = "enhancer")
str(DT_gryphon[, c("ANIMAL", "BWT")])

Fit the animal model

The univariate Gaussian animal model is BWT ~ 1 + animal, by REML. With a raw pedigree you would write hsquared(BWT ~ animal(1 | id, pedigree = ped), data = dat); for gryphon the relationship information is supplied via the engine path.

# Requires a local Julia + HSquared.jl + JuliaCall.
fit <- hsquared(BWT ~ animal(1 | id, pedigree = ped), data = dat)
fit

The first fit in a session takes about 30 s while Julia starts and compiles. Later fits take under a second for a few hundred animals. The call prints no progress, so the first one can look frozen.

Read the result

The published REML estimates are σ²a = 3.3954, σ²e = 3.8286, giving h² = 0.470 — birth weight is moderately heritable. summary() reports the point estimates and, when a local engine provides them, the experimental uncertainty surfaces (asymptotic, REML-only — see Validation evidence):

summary(fit)
heritability(fit)          # point estimate (≈ 0.470)
heritability_interval(fit) # experimental logit delta-method CI, if available
heritability_standard_error(fit)
variance_component_standard_errors(fit)

The interval and standard errors are reported as point-estimate-plus-bounds only; they are not coverage-calibrated and should not be read as validated inference.

Breeding values and accuracy 

ebv <- breeding_values(fit)   # one predicted additive genetic value per animal
head(ebv)
accuracy(fit)                 # sqrt(reliability), when reliability is present

A diagnostic look 

plot(fit, type = "variance")   # components, with experimental +/- 1.96 SE whiskers if present
plot(fit, type = "residuals")  # residuals vs fitted

Interpreting h² = 0.47

A heritability near one-half means roughly half the variation in birth weight among these animals is additive genetic — there is appreciable raw material for a response to selection, and the predicted breeding values rank animals by genetic merit for the trait. As always, h² is a property of this population in this environment, not a fixed constant of the trait.

Honest boundaries

  • The fit requires a local engine; this article uses the supplied A_gryphon relationship matrix as the signed-off recovery path. The open gryphon pedigree has no ancestral loops. Your own analyses will typically use animal(1 | id, pedigree = ped).
  • The CI/SE surfaces are experimental (engine row V1-HERIT-CI, partial): asymptotic, REML-only, unreliable at small n. They are reported, not validated. The Wald SE and logit-delta CI are not valid when additive variance is near 0. The likelihood-ratio test of V_A = 0 uses a 50:50 mixture of chi^2_0 and chi^2_1 (Self and Liang 1987). R anova() does not run it; Julia nested_lrt(..., df = 1, boundary_df = 1) already does.
  • gryphon is teaching/simulated data; the published numbers are an external anchor, not a known-truth simulation. For bias/coverage and external-package agreement see Validation evidence and the benchmark.

See also