This document follows a genetic evaluation in the order it actually happens, and at each step it says three things: what the matrices are and what they mean, what the algorithm does with them, and which function runs it. The formulas are the ones the package computes, not the textbook version they descend from; where the implementation makes a choice, the choice is stated.
The data for every example is simulated here, so each result has its truth beside it.
s <- simulate_breeding(n_founders = 60, n_generations = 3,
offspring_per_generation = 150, h2 = 0.4,
n_markers = 600, seed = 1)
str(s, max.level = 1)
#> List of 4
#> $ pedigree :'data.frame': 510 obs. of 3 variables:
#> $ data :'data.frame': 510 obs. of 3 variables:
#> $ tbv : Named num [1:510] -0.396 0.116 -0.528 1.009 0.208 ...
#> ..- attr(*, "names")= chr [1:510] "F001" "F002" "F003" "F004" ...
#> $ genotypes:List of 21. The model
Every fit in this package is a linear mixed model,
y = X b + Z u + e, u ~ N(0, D), e ~ N(0, R)where y holds the observations, b the fixed
effects, u all the random effects stacked, and
X and Z their incidence matrices. What
separates one model from another is not the algebra but the shape of
D and R: an animal model puts the pedigree
relationship matrix inside D, a repeatability model adds a
second block, a multi-trait model gives R a full covariance
between traits, a reaction norm gives each animal a vector of
coefficients instead of a scalar. Most of the literature calls that
covariance G; this document reserves G for the
genomic relationship matrix of section 6 and writes D for
the covariance of the random effects throughout.
The package writes all of that in one formula. Fixed terms are bare,
cov() marks a fixed covariate, and the marked terms are
random: animal(), maternal(),
sire() carry a relationship matrix, pe() and
random() are independent, rn() is a random
regression, indirect() is the associative effect of pen
mates.
pe() is worth its full definition, because it is the
term the repeatability model is named after. It is the permanent
environment of a subject with repeated records: what all the records of
that subject share and is NOT additive genetic, so the non-additive
genetic effects land there too (Mrode and Pocrnic, 2023, Eqn 5.1). The
repeatability itself is (var(animal) + var(pe)) / total,
the two share cells of the printed table added
together.
The unit of layout is the covariance group, not the
term. Two terms that share group = go into ONE
covariance matrix with the covariance between them estimated:
Var(u_g) = C_g (x) K_gwhere C_g is the small covariance among the terms of the
group and K_g is the relationship matrix among their
levels; D is the block diagonal of these, one block per
group. This is why direct-maternal, the reaction norm and the
associative model have no dedicated fitter: they are the same engine
with a different incidence and the same penalty.
2. The pedigree: A, F, and A inverse
The additive relationship matrix A holds twice the
kinship between every pair of animals. It is never built. What the
equations need is its inverse, and Henderson (1976) showed that the
inverse has a closed form: for each animal i with Mendelian
sampling variance d_i,
A^-1 = sum_i (1/d_i) v_i v_i', v_i = +1 at i, -1/2 at sire, -1/2 at damso A^-1 is built by adding one small outer product per
animal, sparse from the start, never through a dense A. The
Mendelian variance depends on how much of the parents’ ancestry is
already shared:
d_i = 1/2 - 1/4 (F_sire + F_dam) both parents known
d_i = 3/4 - 1/4 F_known one parent known
d_i = 1 neitherF is Wright’s (1922) inbreeding coefficient, computed by
the Meuwissen and Luo (1992) trace: for animal i, walk the
ancestors in decreasing index order accumulating l_j
contributions, and F_i = sum_j l_j^2 d_j - 1. A max-heap
over the ancestors is what makes a single pass correct: every descendant
of an ancestor has a larger index and has therefore already deposited
its share. This trace and the earlier strategy of Quaas (1976) share the
point: the diagonal of A comes out of the decomposition
while A^-1 is being built, instead of storing
A to read it back. Mrode and Pocrnic (2023, ch. 3) carry
both rules through a small pedigree by hand, which is where an
implementation gets checked against printed numbers.
p <- pedigree(s$pedigree)
head(p, 3)
#> Pedigree with 3 animals, 3 founder(s)
#> Mean F 0, maximum 0, 0 animal(s) with F > 0
#> id sire dam F
#> 1 F001 NA NA 0
#> 2 F002 NA NA 0
#> 3 F003 NA NA 0
ai <- a_inverse(s$pedigree) # i, j, x, n: the lower triangle in triplets
c(animals = ai$n, nonzeros = length(ai$x), mean_F = round(mean(p$F), 4))
#> animals nonzeros mean_F
#> 510.0000 1842.0000 0.0063Two details in the implementation are worth knowing because they are
silent when wrong. A parent that is cited but has no line of its own is
an error, not a new founder. And when sire and dam are the same animal,
the coefficients are merged BEFORE the outer product, which puts
1/d on that parent’s diagonal; emitting the terms one by
one into a half-stored triangle would give 3/4 d^-1 and no
error message.
Metafounders. When the base population is not one
homogeneous unrelated pool (different lines, different years, animals
whose parents are simply unknown), Legarra et al. (2015) replace the
unknown parents with metafounders: virtual base animals with a
self-relationship gamma that measures how related that base
actually was. The implementation reduces to one line: a metafounder row
is preloaded with F = gamma - 1 and d = gamma,
and because every formula above reads only the parents’ F,
the whole generalization propagates through the Mendelian variances, the
trace and A^-1 with no special case. gamma
must be strictly inside (0, 2), and the two ends are
refused for different reasons. The upper one is Mendelian: the offspring
of two copies of the same metafounder has
d = 1/2 - 1/4 (2 (gamma - 1)) = 1 - gamma/2, which vanishes
at gamma = 2, a base pool completely homozygous by descent,
where every base animal is the same animal. The lower one is that a base
pool needs a positive self-relationship to be a pool at all, and the
package’s error message names both. Zero is a limit with a meaning
rather than a degeneracy: gamma -> 0 is complete
heterozygosity by descent, an unrelated base, and the recursion
collapses onto the classical unknown parent of Meuwissen and Luo (1992)
(Legarra et al., 2015; Mrode and Pocrnic, 2023, sec. 3.8).
ped_mf <- s$pedigree
ped_mf$sire[ped_mf$sire == "0"] <- "L1"
ped_mf$dam[ped_mf$dam == "0"] <- "L1"
p_mf <- pedigree(ped_mf, metafounders = "L1", gamma = 0.7)
round(mean(p_mf$F[p_mf$id %in% s$pedigree$id]), 4) # base relatedness now shows up as F
#> [1] 0.3541Sires and maternal grandsires. A sire model often
has only the sire and the maternal grandsire of each animal. Then
u_i = u_s / 2 + u_k / 4 + m_i, with the maternal granddam
unknown: the grandsire path weighs 1/4, the Mendelian variance is
11/16 - F_s/4 - F_k/16, and the inverse gets
(1, -1/2, -1/4) in the rules above (Mrode and Pocrnic,
2023, secs. 3.6 and 3.7). Read as a dam, the grandsire would enter with
1/2, and the A^-1 would still be symmetric and positive
definite, so nothing downstream could catch it. The format is therefore
declared, with sire_mgs() or
pedigree(type = "sire_mgs"), and a third column named like
a grandsire column is refused until it is.
3. The mixed model equations
Henderson’s equations (Henderson, 1950) solve for b and
u at once. Writing W = [X Z], the package
assembles
C_s = W'W + blockdiag_g [ kron(s2e * C_g^-1, K_g^-1) ], rhs = W'yand solves C_s sol = rhs. Two conventions matter. The
system is scaled by the residual variance (the penalty carries
s2e * C_g^-1, not C_g^-1), so the solution is
the same BLUP but the matrix is in units of s2e. And inside
a term, the column of Z for coefficient a and
level l is a * n_levels + l: the level
varies fastest. Reversing that convention does not break
anything visibly: the model still assembles and converges, and estimates
something else.
W'W is accumulated one record at a time. Each record
touches a handful of columns and contributes a small clique among them,
so the cost is the sum of squared row-lengths, not anything quadratic in
the number of columns. The penalty is added afterwards into the same
triplet list, and the triplet-to-sparse conversion sums
duplicates, which is not a convenience but a requirement: in
A^-1 the mate-mate term arrives once per offspring, and a
build where the last write wins would silently discard part of the
matrix.
The restricted likelihood, from the MME:
-2logL = (n - p) log s2e + log|D| + log|C_s| + (e'e + u'D^-1 u) / s2ee = y - Xb - Zu are the residuals at the MME solution.
At that solution e'e + u'D^-1 u equals
y'y - sol'rhs, the form the MME suggest, but the package
forms it from the residuals. y'y - sol'rhs subtracts two
numbers of the order of n (mean y)^2 and inherits the error of the
solution at first order. So a constant added to y, which the intercept
absorbs and REML ignores, used to move the -2logL: 1.1e-8 relative at
1000 standard deviations and 1.1e-6 at 1e4. The residual form is
stationary in the solution and does not carry the mean.
D is the covariance of the random effects of section 1,
and it carries the same s2e scaling as C_s:
per group the code accumulates
nl log|C_g / s2e| + dim log|K_g|. With
weights =, a - sum_i log w_i term joins the
expression, the Jacobian of the weighted rescaling, constant in theta,
so it moves the value and not the optimum.
The package computes the same quantity a second way, through the
dense V = Z D Z' + R, and the tests demand the two agree.
They share nothing but the model, so agreement is evidence; a single
path checked against itself would be none.
d <- s$data
ev <- eval_internal(y ~ cg + animal(id), d, s$pedigree, theta = c(0.4, 0.6))
c(mme = ev$neg2logl, dense_V = ev$neg2logl_V, difference = abs(ev$neg2logl - ev$neg2logl_V))
#> mme dense_V difference
#> 4.681071e+02 4.681071e+02 1.597300e-114. Solving: sparse Cholesky and the selected inverse
C_s is factored as L L' by an up-looking
sparse Cholesky, preceded by a minimum degree ordering (George and Liu,
1989) that reduces fill-in and by a symbolic analysis that computes the
pattern of L once. On a clique, minimum degree would cost
the cube of its size rather than scale with the nonzeros, so a node
whose degree passes 80 percent of the nodes still alive leaves the graph
and goes to the end of the order; at the start, a node is set aside when
its degree passes both 10 sqrt(n), the dense-row rule of
AMD (Amestoy, Davis and Duff, 1996), and 0.8 of the 99th percentile of
the degrees, provided the set aside is at most a quarter of the graph.
Ordering never changes a result, only speed and fill. The symbolic
factorization is cached across REML iterations, which is sound only
because the pattern of C_s does not depend on the variance
components: the assembly pushes the same slots every time, explicit
zeros included.
What the dense nodes leave behind is a completely full block at the
end of L: k columns that cost about
k^3 / 3 flops per factorization, which dominates the fit
once k is large. The fit reports it as
fit$dense_block = c(dense = k, columns = n); in a single
step without APY, k is at least the number of genotyped
animals. That block is factored by scalar C++, not by BLAS, and nothing
in the package runs in parallel.
Prediction error variances need the diagonal of C_s^-1,
and the full inverse of a sparse matrix is dense, and for a pedigree of
any size unaffordable. The identity of Takahashi, Fagan and Chin (1973)
gives exactly the elements of the inverse that sit on the pattern of
L, which is all the diagonal needs:
ch <- sparse_chol(ai)
si <- selected_inverse(ai)
c(logdet = round(ch$logdet, 4), dense_block = ch$dense_block)
#> logdet dense_block
#> 312.6306 69.0000A position outside the pattern is not zero; it is unknown. The function returns only what it computed.
5. Estimating the variance components
REML maximizes the likelihood of error contrasts instead of the
likelihood of y. What that removes is not a bias in
b (the fixed effects come out as the same GLS solution
under either criterion) but the downward bias maximum likelihood carries
IN THE VARIANCE COMPONENTS, whose cause is that those fixed effects were
estimated from the same records. Mrode and Pocrnic (2023, ch. 17) write
the restricted likelihood as three terms: a weighted sum of squares of
the residuals, log|V|, and log|X'V^-1X|, the
last of them a penalty for having estimated b and the one
maximum likelihood leaves out (Patterson and Thompson, 1971). The
n - p in the exponent of section 3 is that same correction
seen from the side of the mixed model equations.
The package maximizes that restricted likelihood by average information REML (Gilmour, Thompson and Cullis, 1995; the animal-breeding derivation is in Mrode and Pocrnic, 2023, ch. 17), which needs three ingredients per iteration.
The score. For a covariance group with matrix
C_g, define the quadratic of the solution and the trace
from the selected inverse,
Q[a,b] = u_a' K^-1 u_b
T[a,b] = s2e * tr( K^-1 [C_s^-1]_{a,b} )
M_g = nl * C_g^-1 - C_g^-1 (Q + T) C_g^-1and the derivative of -2logL with respect to the element
(i,j) of C_g is M_g[i,i] on the
diagonal and M_g[i,j] + M_g[j,i] off it, because a
symmetric element appears twice. For the residual,
d(-2logL)/ds2e = (n - p - q + trpen)/s2e - e'e/s2e^2where q is the number of random columns and
trpen the trace of the penalty against C_s^-1.
The - q + trpen is not decoration: the random effects
absorb degrees of freedom too, and a version without it still converges,
to the wrong place, reporting convergence. That is the worst failure
mode available, which is why the score is checked against central finite
differences on every parameter.
The step. The average information matrix is the
average of the observed and expected information, cheap because it needs
no second derivatives, and the step is
theta <- theta - (AI + lambda diag(AI))^-1 score. This
lambda is the Levenberg-Marquardt damping (Levenberg, 1944;
Marquardt, 1963), and it is a third meaning of the letter: not the
s2e/s2a ratio that scales K^-1 inside the
mixed model equations of section 3, and not the ridge that
apy_inverse() calls lambda. The damping starts
at 1e-3, multiplies by ten whenever a candidate is rejected
(inadmissible covariance, singular damped matrix, or a likelihood that
got worse) and divides by ten when one is accepted. Before the AI loop,
a few EM rounds run: EM (Dempster, Laird and Rubin, 1977) is monotone
and stays inside the cone, so it costs little and lands the search in
the right basin.
Stopping. Convergence is judged on the relative
change of the components,
sqrt(sum(dtheta^2)/sum(theta^2)) < tol, never on an
absolute threshold on the score, whose norm grows with the number of
records. A small step is not an optimum, though, so
converged also requires the Newton decrement
g' AI^-1 g, about twice the remaining gap in -2logL, to
fall under 2e-4, components resting at a boundary excluded (in
model_mt() and model_ar1(), after checks of
the Karush-Kuhn-Tucker condition and of the boundary face on the
likelihood, a second walk from inside the cone and a polish along the
face); it is reported as newton_dec. At
maxiter with the decrement already under tolerance the
likelihood is flat along some direction, and the fit declares
convergence and says so. Standard errors come from the inverse of the AI
matrix, SE_k = sqrt(2 [AI^-1]_kk). When AI is singular at
the end (an eigenvalue below 1e-8 of the largest; in
model_mt() and model_ar1() the matrix is first
scaled by the standard deviation of each trait, so the verdict does not
depend on the unit of a trait), the fit names the components involved
and their standard errors are NaN: the data do not separate them, the
reported point depends on start=, and only combinations are
estimable (for the indirect-effect model, the identifiability question
of Cantet and Cappa, 2008).
fit <- model(y ~ cg + animal(id), d, s$pedigree)
fit
#> AI-REML fit of 'y'
#> converged in 7 iteration(s), relDelta 1.83e-09, Newton decrement 7.26e-19, 0.008 s
#> -2logL 468.1000386
#> 510 record(s), 513 column(s) in the equations
#> fixed column(s) removed for linear dependence: cg=g2
#>
#> component estimate std_error share correlation
#> var(animal) 0.400864 0.0956138 0.4024
#> var(residual) 0.595274 0.0688334 0.5976
#>
#> fixed effects (implicit intercept; dropped columns are zero; compare by contrast):
#> intercept cg=g1 cg=g3
#> 10.00302128 0.03647810 0.00221236Two rules protect the result. A covariance that stops being
positive-definite ends the fit rather than being nudged back, and the
starting point comes from var(y) unless start=
gives one, so no run inherits a neighbour’s answer.
6. Genomics: from markers to H
G. VanRaden (2008) builds the genomic relationship
matrix from centred marker dosages. With p_j the allele
frequency at marker j and Z = M - 2p,
G = Z Z' / (2 sum_j p_j (1 - p_j))Missing genotypes are imputed with the marker mean and never with
zero: zero is a valid genotype, and confusing the two changes both the
frequencies and every element of G with no error to notice.
Markers that are fixed or entirely missing are dropped and counted.
A22, and the trap. The single step needs the
pedigree relationships among the genotyped animals,
A22. That is the inverse of the genotyped block of the
inverse, not the genotyped block of A^-1, two matrices of
the same shape, both symmetric, only one correct. The package computes
it by a sparse Schur complement,
A22^-1 = B22 - B21 B11^-1 B12where B11 is the non-genotyped block.
B11^-1 is never formed: one sparse factorization and one
triangular solve per genotyped animal do the same work.
H. Genomic and pedigree information describe the
same animals on two different scales, so G is first brought
to the scale of A22 by an affine adjustment that matches
both the mean diagonal and the mean off-diagonal, then blended:
b = (mean_diag(A22) - mean_off(A22)) / (mean_diag(G) - mean_off(G))
a = mean_off(A22) - b * mean_off(G)
G* = (1 - w) (a + b G) + w A22The blend w (5 percent by default) keeps G*
invertible and admits that markers are not the whole truth. Then the
single-step inverse of Aguilar et al. (2010) and Christensen and Lund
(2010) says the whole relationship matrix H has an inverse
that differs from A^-1 only in the genotyped block:
H^-1 = A^-1 + [ 0 0 ]
[ 0 G*^-1 - A22^-1 ]This is the identity that makes the single step practical: nothing
dense is needed for the non-genotyped animals, and the correction is a
small block added to a sparse matrix. It is also self-checking: with
blend = 1, G* becomes A22
exactly, the correction block vanishes, and H^-1 must equal
A^-1 through the whole fit. That collapse is one of the
package’s tests. Mrode and Pocrnic (2023, ch. 12) work the same identity
through a small herd, number by number.
gen <- list(ids = s$genotypes$ids, m = s$genotypes$m)
f_ss <- model(y ~ cg + animal(id), d, s$pedigree, genotypes = gen)
f_id <- model(y ~ cg + animal(id), d, s$pedigree, genotypes = gen, blend = 1)
c(single_step = f_ss$neg2logl, blend_one = f_id$neg2logl, pedigree = fit$neg2logl)
#> single_step blend_one pedigree
#> 480.9321 468.1000 468.1000Inverting G* without paying the cube. Two approximations are available, and they are the same idea at different granularity. APY (Misztal, Legarra and Aguilar, 2014; Mrode and Pocrnic, 2023, sec. 12.4) splits the genotyped into a core and the rest; every non-core animal is regressed on the core, and what is left is a scalar Mendelian residual:
P = Gcc^-1 Gcn
m_j = g_jj - g_jc P_j
G_APY^-1 = [ Gcc^-1 + P Mnn^-1 P' -P Mnn^-1 ]
[ -Mnn^-1 P' Mnn^-1 ]The Vecchia recursion (Vecchia, 1988) generalizes it: each animal
conditions on its own k nearest relatives among the
previous ones rather than on one shared core. Henderson’s
A^-1 is this same recursion with the parents as the
conditioning set, exact there, because a pedigree is Markovian. Schafer,
Katzfuss and Owhadi (2021) prove the factor is optimal in
Kullback-Leibler divergence given the pattern, and that enlarging the
sets never hurts.
An honest limit: both approximations cut arithmetic, not memory.
G*^-1 is a dense matrix of the genotyped animals on all
three paths. What H^-1 keeps of it is only the nonzeros, so
with APY the dense block of the factor is at least the core size plus
one rather than the whole genotyped set.
How large a core? The dimensionality of the genomic information
bounds it: Pocrnic et al. (2016) take as many animals as eigenvalues of
G explain 98% of its trace, drawn at random.
apy_core_select() applies that rule, reports the cost of
one factorization with the core against the exact one, and takes
size = and include = (proven sires, say) to
override the draw; apy_core = "auto" calls it inside the
fit, which then keeps the core in fit$apy. Choose the core
once and pass it to every fit.
gid <- gen$ids
f_apy <- model(y ~ cg + animal(id), d, s$pedigree, genotypes = gen,
apy_core = gid[seq(1, length(gid), by = 3)])
f_vec <- model(y ~ cg + animal(id), d, s$pedigree, genotypes = gen, vecchia_k = 60)
round(c(apy = cor(ebv(f_ss, "animal"), ebv(f_apy, "animal")),
vecchia = cor(ebv(f_ss, "animal"), ebv(f_vec, "animal"))), 4)
#> apy vecchia
#> 0.9753 0.9675Markers as equations. The other route never forms
G at all. Liu et al. (2014) write the equivalent model with
the marker effects as unknowns, u_g = Z g + p_g, and the
joint precision turns the Schur term into blocks of A^-1,
so most of the system is the same sparse matrix the pedigree already
gives. snp_blup() solves it by preconditioned conjugate
gradients, which need only matrix-vector products, and every application
of A22^-1 uses the matrix-free identity of Masuda et
al. (2017):
A22^-1 v = A^22 v - A^21 (A^11)^-1 A^12 vone sparse factorization of the non-genotyped block, one triangular solve per iteration. Nothing dense is ever built. The variance components are given rather than estimated: this is the solver of routine evaluation, where components come from a periodic REML run.
7. What the formula can say
With the covariance group as the unit, several classical models are incidences rather than fitters.
Direct-maternal (Willham, 1972): the animal’s own genes and its dam’s maternal ability in one 2x2 matrix, with the direct-maternal covariance, usually negative, and the number the model exists for.
Reaction norm (Kirkpatrick, Lofsvold and Bulmer, 1990; in
animal models, Schaeffer and Dekkers, 1994): each animal gets a vector
of coefficients on an environmental gradient through
rn(id, base = ...), with legendre() building
the orthogonal basis. The genetic variance becomes a function of the
environment, which h2_curve() reads, evaluating every
random regression at the point, a permanent environment on the same
basis included; h2() refuses a reaction norm, because its
heritability is a curve and not a number.
Indirect genetic effects (Griffing, 1967; Muir and
Schinckel, 2002; Bijma, Muir and Van Arendonk, 2007): the incidence of a
record marks the animal’s PEN MATES, so a phenotype carries its own
genes plus the associative effects of those it lives with. Direct and
indirect share a group, so the fit returns var(animal),
var(indirect) and the covariance between them, and the SIGN
of that covariance is what separates heritable competition from
heritable co-operation.
What decides the response is not that covariance on its own but the total breeding value, into which the group size enters twice:
TBV = A_D + (n - 1) A_S
s2_TBV = s2_AD + 2 (n - 1) s_ADS + (n - 1)^2 s2_AS
s2_P = s2_AD + (n - 1) [1 + (n - 2) r] s2_AS + 2 (n - 1) r s_ADS + (the other components)
T2 = s2_TBV / s2_PBecause (n - 1) multiplies the covariance and
(n - 1)^2 the associative variance, the group size sets
both the magnitude of the response and, under strong heritable
competition, its direction: the response need not follow the selection
pressure (Bijma, Muir and Van Arendonk, 2007; Mrode and Pocrnic, 2023,
ch. 9). The phenotypic variance of a record carries the indirect effects
of its n - 1 mates, with r the relationship
between them, so heritability here is not var(animal) over
the plain sum of the components, and the share column of
the printed fit is left blank. h2(fit, n = , r = ) divides
by s2_P, and t2(fit, n, r) returns
s2_TBV, s2_P, T2 and the direct heritability
with delta-method standard errors.
The book fixes n = 3; real pens are unequal, and with a
coefficient of 1 per mate the social sum grows with
n_i - 1.
indirect(id, pen = "pen", dilution = d) scales every mate’s
entry to (n_i - 1)^(-d): d = 0 (the default)
is the book’s sum, d = 1 the mate mean, and the regime in
between is the dilution of Bijma (2010), chosen empirically by fitting a
small grid of d and comparing -2logL, which
stays comparable because only the incidence changes. Under dilution the
(n - 1) of the TBV and of the covariance term of
s2_P becomes (n - 1)^(1 - d), and the one that
multiplies s2_AS in s2_P becomes
(n - 1)^(1 - 2d) (Bijma, 2010); h2() and
t2() read d from the term. The residual
inherits the same heterogeneity,
var(e_i) = s2_ED + (n_i - 1)^(1 - 2d) s2_ES, the same
d diluting the social environmental deviation;
indirect_residual() estimates the ratio
s2_ES / s2_ED by profile REML over weighted fits, and its
documentation says which of the two the data actually pins.
associative_matrix(dilution = d) writes the same residual
as an exact covariance for a kernel() term.
Whether the three genetic components can be told apart is a property
of the design, not of the fitter (Cantet and Cappa, 2008). With pens of
one size built from two full-sib families, the data identify only two
combinations, the TBV variance and s2_AD - 2 s_ADS + s2_AS:
fits from different starts land on different components with the same
-2logL, and the fit says SINGULAR and names them. Pens of different
sizes, or relatives spread across pens, separate them.
Multi-trait: cbind() on the left, a full
residual covariance between traits, and missingness handled by the
pattern of each record rather than by dropping it.
AR(1)/CAR(1) residual (Wade and Quaas, 1993): for repeated
measures along time, the residual correlation decays as
rho^|dt|. The inverse of that correlation matrix is
tridiagonal (the classical Markov result), so the model costs no more
than the iid one, and rho comes out with its standard
error. On a grid, rho^k is the correlation of records k
steps apart and rho may be negative; it is negative only
between records an odd number of steps apart. The step is the largest
number that divides every gap (Pinheiro and Bates, 2000,
corAR1; the AR1 of ASReml for equally spaced
levels). Off a grid the correlation is rho^dt with
rho in [0, 1) (corCAR1; EXP of
ASReml), since a continuous-time Markov process has no negative
correlation. A few records whose time is off the grid of the others can
therefore turn a negative correlation into a wrong positive or null one;
model_ar1() checks for that after the fit and stops naming
those records.
The textbook treatment of these models is Mrode and Pocrnic (2023): the repeatability model in ch. 5, maternal traits in ch. 8, social interaction in ch. 9, random regression in ch. 10.
model(weight ~ cg + animal(id, group = "g") + maternal(dam, group = "g"), d, ped)
model(y ~ cg + rn(id, base = c("phi0", "phi1")) + pe(id), d2, ped)
model(y ~ cg + animal(id, group = "g") + indirect(id, pen = "pen", group = "g"), d, ped)
model_mt(cbind(t1, t2) ~ cg + animal(id), d, ped)
model_ar1(y ~ cg + animal(id) + pe(id), d, ped, subject = "id", time = "day")8. Reading the evaluation
The solution of the MME gives the breeding values directly. Their reliability comes from the prediction error variance, which is where the selected inverse pays off, and the accuracy of an animal’s EBV is
acc = sqrt( 1 - PEV / ((1 + F) sigma2_a) )with the inbreeding coefficient in the denominator because an inbred
animal has more additive variance to predict:
var(u_i) = a_ii s2a = (1 + F_i) s2a (Mrode and Pocrnic,
2023, ch. 3). The identity underneath,
PEV = var(a - a^) = C22 s2e, is Henderson (1975).
accuracy() divides by the prior the animal actually has.
For a pedigree animal that is 1 + F. For a genotyped animal
in a single step it is the diagonal of H, which in that
block is the diagonal of G*, and that is not
1 + F_ped: on the population simulated at the top of this
vignette, every animal genotyped, the two diagonals differ by up to 0.18
and an individual accuracy by up to 0.067. The herd mean does not move;
the individual does. A fit with genotypes= carries the
diagonal of G*, and accuracy() uses it.
Heritability is read off theta, and both halves of the
ratio have to be picked rather than assumed. theta is a
NAMED vector, and its first entry is the animal variance only when the
formula happens to put that term first:
random(lot) + animal(id) puts var(random)
there. It also carries the rho(residual) of
model_ar1(), a correlation that belongs in no variance
denominator. The covariances of a group = are another
matter. A direct-maternal covariance is part of the phenotypic variance
of a record, with coefficient 1 (Willham, 1972); in the associative
model the coefficients are those of s2_P in section 7,
which depend on the pen size. h2() takes the numerator by
name and divides by the phenotypic variance of the trait, covariances
included, and the share column of the printed fit uses the
same denominator, left blank where that denominator needs the pen
size.
The fixed effects come with the same solution vector:
fit$b, named term=level, also reachable as
coef(fit, effects = "fixed"). The parametrization is the
implicit intercept, with linearly dependent columns dropped
(dropped_x) and a dropped level holding the zero solution,
while the book’s examples zero one level per factor and fit no
intercept. The two conventions agree only on CONTRASTS, differences
between levels of the same factor, never on the raw values; contrasts
are invariant to the parametrization, individual entries are not.
e <- ebv(fit)
a <- accuracy(fit, s$pedigree)
round(c(cor_with_truth = cor(e[names(s$tbv)], s$tbv),
mean_accuracy = mean(a), h2 = h2(fit)), 3)
#> cor_with_truth mean_accuracy h2
#> 0.719 0.736 0.402From there the decisions: selection_index() combines
traits with economic weights, and rank_drift() answers what
an updated evaluation actually changed.
9. The Bayesian alternative
The same models can be sampled instead of maximized.
gibbs() draws all locations in one block through the same
sparse Cholesky, then draws each covariance from its conjugate
inverse-Wishart and the residual from a scaled inverse chi-square. The
prior is a choice, prior =: "jeffreys", the
default, is p(s2) ~ 1/s2; "flat" is uniform on
the variance; "uniform_sd" is uniform on the standard
deviation, the one Gelman (2006) recommends for a variance with few
levels; and c(df =, scale =) is the proper scaled inverse
chi-square of the BLUPF90 Gibbs programs. The improper ones can leave
the posterior of a weakly identified component with spurious mass near
zero (Hobert and Casella, 1996). For a covariance group the default is
worse: |C|^-(d+1)/2 grows without bound as C
nears a singular matrix, so the posterior has infinite mass at a
correlation of +-1 even when the data identify every component, and a
chain that gets there stays. With direct and indirect effects in one
group on pens of 2 to 8 (1600 records, 30 data sets) that happened in 4
of 30 under "jeffreys" and in none under
"flat", whose posterior means landed beside REML;
gibbs() warns when a chain sits there. On a Gaussian trait
the flat prior or a proper one avoids that boundary; with
family = "probit" only a proper one does, because the
binary likelihood stays positive as the variances grow and the flat
prior leaves that posterior improper too. The chain also has no SINGULAR
verdict of its own, so a design that only the likelihood fails to
identify runs without one (pens all of 4 holding two full sibs from each
of two families, or records on unrelated founders only: the message of
the chain was empty in 29 and in 30 of 30 data sets, where
model() on the same formula says SINGULAR in all 60). Where
REML gives a point and a standard error from the information matrix, the
sampler gives the whole posterior, worth the cost when a component sits
near a boundary, where the normal approximation behind the standard
error is least trustworthy. summary() of the chain gives
the posterior quantiles and the heritability computed draw by draw. The
sampler goes back to Geman and Geman (1984); Wang, Rutledge and Gianola
(1993, 1994) brought it to variance components in animal models, and
Mrode and Pocrnic (2023, ch. 18) set out the same sampler and the
conjugate draws it uses. With family = "probit" it samples
a binary trait on its liability, each record’s liability drawn from a
normal truncated at the threshold and the residual fixed at 1 (Albert
and Chib, 1993; Sorensen et al., 1995), which avoids the Laplace
approximation. It does not rescue every design: with
indirect() and the components estimated on binary records
(3200 records, 30 data sets) the chains reached a bulk ESS of 400 in
none, and the Laplace minimum sat at |r| > 0.99 in 29,
where REML on the liability itself recovers the components.
g <- gibbs(y ~ cg + animal(id), d, s$pedigree, n_iter = 2000, burnin = 400, thin = 4)
data.frame(component = colnames(g$samples),
posterior_mean = round(unname(g$mean), 3),
reml = round(unname(fit$theta), 3))
#> component posterior_mean reml
#> 1 var(animal) 0.411 0.401
#> 2 var(residual) 0.593 0.59510. A declared covariance: dominance, epistasis, breed origin
Sections 2 and 6 derived the covariance of a random term from the
pedigree or from the markers. kernel(id, K = ) removes the
derivation: the term’s covariance matrix is DECLARED, symmetric
positive-definite, its rownames naming the levels, and every animal in
K gets an equation whether or not it has a record. Nothing
else changes (the same kron(C^-1, K^-1) penalty, the same
score, the same selected inverse), so any model whose only novelty is
the covariance structure is a constructor plus a marker, not a new
fitter (Mrode and Pocrnic, 2023, chs 13-14).
The constructors cover the classical cases.
dominance_matrix() builds Cockerham’s (1954) D from the
pedigree, d_xy = 0.25 (a_sf a_dm + a_sd a_fm) with diagonal
1, the covariance of the dominance deviations, fitted BESIDE the
additive term so the two partition the genetic variance. The genomic
versions are g_dominance() (Vitezica, Varona and Legarra,
2013) and g_epistasis(), the Hadamard square of
G with its diagonal averaged to one;
genomic_inbreeding() reads f = 1 - het/N per
animal, and fitted as a cov() its coefficient in
fit$b is the inbreeding depression.
For a multibreed pedigree, partial_a() splits
A by breed origin (Garcia-Cortes and Toro, 2006): one
partial matrix per founder breed and one per segregating pair, built by
the tabular recursion with the breed contribution on the diagonal. A
purebred contributes nothing to the other breed’s matrix, and its row
there is entirely zero, which kernel() reads as “no
equation for this level, keep its records”; so each matrix enters as its
own term and the breed-specific solutions sum to the combined breeding
value.
D <- dominance_matrix(s$pedigree)
f10 <- model(y ~ cg + animal(id) + kernel(id, K = D), d, s$pedigree)
round(f10$theta, 3)
#> var(animal) var(kernel) var(residual)
#> 0.400 0.434 0.164Read that split with care. This data has ONE record per animal and
was simulated with no dominance, and D has a unit diagonal,
so var(kernel) and var(residual) are separated
only by the sib covariances, a weakly identified partition, and the
likelihood is nearly flat between them. The estimate above is one point
on that plateau, not a discovery of dominance. The hands-on vignette
plants a real dominance signal on repeated records, where the partition
is estimable and the fit recovers it.
11. The threshold model: a categorical trait on a liability
A calving score has ordered categories, not a variance. The threshold
model puts a Gaussian LIABILITY l = Xb + Zu + e underneath
and maps it to the observed category through thresholds:
y = j when t_{j-1} < l <= t_j. Only
signs and orderings of l are observable, so the residual is
FIXED at 1 (probit identifiability) and the thresholds replace the
intercept. model_threshold() fits it by Fisher scoring on
the Gianola and Foulley (1983) system, with the variance components on
the liability scale GIVEN via start=: the Gaussian REML of
section 5 has no likelihood for a categorical trait, and the package
says so rather than pretending. With estimate = TRUE the
components are estimated instead, by approximate marginal maximum
likelihood: the locations are integrated by the Laplace approximation
around the mode, and the estimate is the MINIMUM of that Laplace -2logL,
the number the fit reports, with standard errors from its numerical
Hessian. (The EM step of Foulley, Im, Gianola and Hoeschele, 1987,
s2 = (u' K^-1 u + tr(K^-1 C^uu)) / q, was the estimator
before; its fixed point is not that minimum, because the step ignores
how the weights move with s2.) Its known weakness is a binary trait with
few records per level of the random effect, where the Laplace
approximation underestimates the variance (Tempelman, 1998); on 80 sires
with 50 daughters each and a planted variance of 0.15 it averaged 0.161
(SE 0.007) over 10 replicates, and gibbs(family = "probit")
avoids that approximation. cbind(quant, binary) fits the
joint quantitative-binary analysis of Foulley, Gianola and Thompson
(1983), with the PEV of u1 and of the ranking value u2, and
predict() the probability of Eqn 15.25; Mrode and Pocrnic
(2023, ch. 15) carry both through worked examples, which are the
package’s gates.
The scale matters when comparing heritabilities: a 0/1 trait analysed
linearly reports an OBSERVED-scale h2, related to the liability one by
the Dempster and Lerner (1950) factor z^2 / (p(1-p)) at
incidence p.
d11 <- d
d11$cat <- as.integer(cut(d11$y, quantile(d11$y, c(0, 0.6, 0.85, 1)),
include.lowest = TRUE))
f11 <- model_threshold(cat ~ cg + animal(id), d11, s$pedigree, start = 0.35,
verbose = FALSE)
round(f11$thresholds, 3)
#> t1 t2
#> 0.249 1.165
round(h2_observed(0.35 / 1.35, 0.4), 3) # what that liability h2 shows on a 0/1 scale
#> [1] 0.16112. Time until failure
Longevity is a time, and part of it is unobserved: a cow still alive
at culling day has a LOWER BOUND on her length of productive life, not a
missing value. The Weibull frailty model (Kachman, 1999; Mrode and
Pocrnic, 2023, ch. 16) writes the hazard as
h(t) = h0(t) exp(x'b + z'u) with a Weibull baseline; a
complete record contributes its density, a censored one its survival
S(t), and that difference is where the censoring
information lives: model_survival() therefore refuses to
run without censor=. Solutions are LOG RELATIVE RISKS
(exp() gives the relative risk of failing), the Weibull
rho is estimated by Newton alongside the effects, and the
frailty variance comes from the minimum of the Laplace -2logL of the
marginal likelihood (one or several components, indirect()
included), or is given, as the book does.
A covariate that changes during a life is constant by pieces, and
each piece is an elementary record (entry, stop] of the
same subject, with the failure only on the last piece. The
cumulative hazard of a piece is
Lambda0(stop) - Lambda0(entry), so a life split into pieces
with the same covariates has exactly the likelihood of the unsplit
record. The pieces are checked: no overlap, the event only at the end,
the frailty constant within the subject, and gaps refused unless
gaps = "allow". A first entry above zero is left
truncation, accepted and flagged, because with a frailty its conditional
likelihood is naive.
cows <- data.frame(
cow = as.character(8:19),
herd = as.character(c(1,1,1,1,1,1,2,2,2,2,2,2)),
ysp = as.character(c(3,4,1,2,3,1,4,1,2,3,4,2)),
code = c(0,1,0,1,1,1,1,1,0,1,0,1),
lpl = c(40,47,22,28,50,33,49,29,23,37,35,30))
peds <- data.frame(
animal = as.character(1:19),
sire = c(rep("0", 7), "1","1","4","4","5","5","1","1","5","5","4","4"),
dam = c(rep("0", 7), "2","3","2","9","3","8","6","7","14","6","7","3"))
f12 <- model_survival(lpl ~ herd + ysp + animal(cow), cows, peds, censor = "code",
rho = 1, lambda = 1, sigma2 = 0.4, verbose = FALSE)
round(coef(f12, "fixed"), 3) # Example 16.1: the p.291 solutions of the book
#> herd=2 ysp=2 ysp=3 ysp=4
#> -1.631 -2.346 -3.149 -2.982References
The complete list is in the package README. The works this
walkthrough leans on directly: Henderson (1950) for the MME; Henderson
(1976) for A^-1; Wright (1922) and Meuwissen and Luo (1992)
for inbreeding; Patterson and Thompson (1971) and Gilmour, Thompson and
Cullis (1995) for REML and average information; Legarra et al. (2015)
for metafounders; VanRaden (2008), Aguilar et al. (2010) and Christensen
and Lund (2010) for the single step; Misztal, Legarra and Aguilar
(2014), Pocrnic et al. (2016) and Schafer, Katzfuss and Owhadi (2021)
for the approximate inverses; Liu et al. (2014), Masuda et al. (2017)
and Vandenplas et al. (2018) for markers as equations; George and Liu
(1989) and Amestoy, Davis and Duff (1996) for the ordering; Takahashi,
Fagan and Chin (1973) for the selected inverse; Willham (1972),
Kirkpatrick et al. (1990), Schaeffer and Dekkers (1994), Griffing
(1967), Muir and Schinckel (2002), Bijma, Muir and Van Arendonk (2007),
Bijma (2010), Cantet and Cappa (2008) and Wade and Quaas (1993) for the
models of section 7; Geman and Geman (1984) and Wang, Rutledge and
Gianola (1993, 1994) for the sampler, Gelman (2006) and Hobert and
Casella (1996) for its priors, and Albert and Chib (1993) and Sorensen
et al. (1995) for its probit mode; Cockerham (1954), Vitezica et
al. (2013) and Garcia-Cortes and Toro (2006) for the declared
covariances of section 10; Gianola and Foulley (1983), Foulley, Gianola
and Thompson (1983), Foulley, Im, Gianola and Hoeschele (1987),
Tempelman (1998) and Dempster and Lerner (1950) for the threshold model;
Kachman (1999) for survival. Mrode and Pocrnic (2023) is the textbook
behind the whole sequence, chapter by chapter.