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.
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.705303e-124. 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 in tiles of 64 x 64 on the threads that
br_threads() sets (OpenMP). Each output number is written
by one thread and every sum runs in a fixed order, so the factor is the
same bit for bit with 1 thread or with 16.
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(), 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() 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.
fit <- model(y ~ cg + animal(id), d, s$pedigree)
fit
#> AI-REML fit of 'y'
#> converged in 7 iteration(s), relDelta 1.55e-09, Newton decrement 5.17e-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. 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 ]An honest limit: the approximation cuts arithmetic, not memory.
G*^-1 is a dense matrix of the genotyped animals on both
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)])
round(c(apy = cor(ebv(f_ss, "animal"), ebv(f_apy, "animal"))), 4)
#> apy
#> 0.9753Markers 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.
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.
The textbook treatment of these models is Mrode and Pocrnic (2023): the repeatability model in ch. 5, maternal traits in ch. 8, random regression in ch. 10.
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. 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). 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.
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" 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; c(df =, scale =) is the proper scaled inverse
chi-square of the BLUPF90 Gibbs programs; and "half_t" is
the hierarchical prior of Huang and Wand (2013), each standard deviation
half-t with its density flat at zero, each correlation uniform, proper
and conjugate, its scale read from the design. 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
"jeffreys" 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;
gibbs() warns when a chain sits there. Left unset, each
kind of variance gets the prior its validation chose: the half-t on a
scalar term (no chain stuck at a small true variance, where
"jeffreys" held 35 of 360 chains at zero),
"flat" on a covariance group (the half-t removed the trap
too but pulled a direct-maternal correlation of -0.8 toward zero, by
0.43 posterior SD, where "flat" stayed at +0.04), and
"flat" on the residual. On a Gaussian trait the default
also runs the interweaving step of Yu and Meng (2011) on every estimated
term (asis): after the usual draw of a covariance, the
effects are rewritten in the non-centred parametrization and the
covariance is drawn again by Metropolis-Hastings, which leaves the
posterior unchanged. Against the same chains without it the bulk ESS per
iteration rose 5.2 to 10.6 times on direct-maternal and
rn() groups, and 2.1 times (a planted variance of 0.05) and
3.5 times (a planted 0) on a scalar term alone. With
family = "probit" the prior has to be proper.
"jeffreys" leaves the posterior improper at a variance of
zero, and where the random effects can reproduce the category of every
record, as with one binary record per animal, the binary likelihood
stays positive as the variances grow and "flat" and
"uniform_sd" leave it improper too (where a level holds
records of both categories the likelihood falls and they need not);
there the default is half_t(df = 4) with a normal prior on
the fixed effects, and "jeffreys", "flat" and
"uniform_sd" are refused unless
allow_improper = TRUE. The chain also has no SINGULAR
verdict of its own, so a design that only the likelihood fails to
identify runs without one (records on unrelated founders only: the
message of the chain was empty in 30 of 30 simulated data sets, where
model() on the same formula says SINGULAR in all 30). 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). Its default is
half_t(df = 4) with a normal prior on the fixed effects
(prior_fixed = 1). The default adds the PX-DA move of Wang
and Roy (2018) (px = TRUE) and the interweaving step on
each covariance group (asis = TRUE). gibbs()
refuses that step under "jeffreys", because its target, the
conditional of chol(C) given the non-centred effects, has
no finite integral.
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.398 0.401
#> 2 var(residual) 0.606 0.595References
The complete list is in REFERENCES.md. 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) for the approximate inverse; 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) and
Schaeffer and Dekkers (1994) for the models of section 7; Geman and
Geman (1984) and Wang, Rutledge and Gianola (1993, 1994) for the
sampler, Gelman (2006), Hobert and Casella (1996) and Huang and Wand
(2013) for its priors, Yu and Meng (2011) for the interweaving step, and
Albert and Chib (1993) and Sorensen et al. (1995) for its probit mode.
Mrode and Pocrnic (2023) is the textbook behind the whole sequence,
chapter by chapter.