Sparse inverse of the dominance relationship matrix, by sire x dam subclasses
Source:R/nonadditive.R
dominance_inverse.RdThe method of Hoeschele & VanRaden (1991), generalized to inbreeding, overlapping
generations, unknown parents and repeated matings. The dominance deviation of an animal
with both parents known is split into the effect of its sire x dam subclass and a
deviation within it, d_i = h_SD + delta_i, with Var(h) = (F / 4) s2d,
F_{SD,XY} = a_SX a_DY + a_SY a_DX, and Var(delta_i) = Delta_i s2d,
Delta_i = 1 - [(1 + F_S)(1 + F_D) + 4 F_i^2] / 4; an animal without both parents has
no subclass and Delta_i = 1. Then D = W (F / 4) W' + diag(Delta) is the Cockerham D
of dominance_matrix() as an identity, inbreeding included, and when every
Delta_i > 0 (see below) the precision of the augmented vector [d; h],
Q = [[Delta^-1, -Delta^-1 W], [-W' Delta^-1, 4 F^-1 + W' Delta^-1 W]], exists and is
sparse. With kernel(id, Kinv = dominance_inverse(ped)) the component is the same s2d
as with kernel(id, K = dominance_matrix(ped)) and the subclass levels are latent
effects, like the ancestors without records in A^-1, so the likelihood is the same.
The tests check, at fixed components, -2logL, score and average information in
model() and model_mt(), -2logL and score in model_ar1(), and in gibbs() the
starting point and, by Monte Carlo, the posterior mean of the effects against the BLUP.
The REML optimum is the same; the path to it need not be (the EM warm-up of model()
moves differently with the latent levels).
Usage
dominance_inverse(
ped,
animals = NULL,
route = c("auto", "dense", "sparse"),
max_pairs = 5e+06
)Arguments
- ped
data.frame with animal, sire and dam, as in
model(); unknown parent 0 or NA- animals
the animals that get a dominance level, by default all of the pedigree. Pass the animals with records: the marginal of D on them is exact and the precision carries only their subclasses. A fit refuses a record with an observed response of an animal left out here, which would otherwise be dropped from the analysis without notice
- route
"auto","dense"or"sparse", see above; all give the same marginal on the animals- max_pairs
ceiling on the pairs the sparse closure may create, a guard on memory: a pair costs a few hundred bytes while the closure is built (an estimate from the layout, about 320, not measured), so the default allows on the order of 1.5 GB. Values above
2^32 - 1, the most pairs the closure can index, act as2^32 - 1
Value
the lower triangle of Q in triplets, list(i, j, x, n, id) as a_inverse()
returns, plus: type, per level, "animal", "subclass" (the sire x dam pair of an
animal in Q) or "ancestral" (a pair of the sparse closure kept as a latent level);
prior, the prior variance of each level at unit s2d (1 for an animal, F_cc / 4
for a pair), which accuracy() divides by; logdet, log|K| of the augmented
covariance by the closed form, which a fit uses instead of factoring Q again while
the triplets and logdet still match signature (any change to either and the fit
factors Q); start_scale, the geometric mean of the
eigenvalues of D[animals, animals] (of D itself without animals =), which the
multi-trait, AR(1) and Gibbs fitters use for their starting values so that they start
where kernel(K = D[animals, animals]) starts; route, "dense", "sparse", or
"none" when no animal of Q has both parents (then Q = I); decision, why that
route: "given" (route =), "compared" (both costed), "closure", "memory",
"ordering_work" and "ordering_cost" (the early exits above, all to the dense
route), "dense_too_big" (the sparse route is the only one) or "no_subclass";
n_subclasses; n_pairs;
pairs, the two animals of each pair level in pedigree order, one row per level named
by its label (a level h of pair x, y has Var(h) = (a_xx a_yy + a_xy^2) / 4 and
Cov(h, h') = (a_xx' a_yy' + a_xy' a_yx') / 4 at unit s2d);
closure (pairs created and processed, and whether it was abandoned); flops, the
sum of squared column counts of the Cholesky factor of the reference model under each
route, and cost, the same with the sparse part of the factor weighted 8 times the
dense final block, which is what "auto" compares (NA when not computed, and the
sparse ones Inf when its ordering was abandoned); work, the list entries read by
the ordering of the sparse candidate (ordering, up to where it stopped) and by that
of the dense one (dense), and the ceiling of the former, the latter plus 0.05
times the dense cost; memory, the estimated peak bytes of building Q by the dense
route; entries, the triplets each route assembles for the pair block of Q, the full
triangle for the dense one and, counted before assembly, the sparse candidate's (NA
when not counted), which the memory exit compares; seconds, the wall time of the
pair closure, of the choice of route and of building Q; signature, a hash of the
triplets and of logdet; and left_out, the pedigree animals without a level. A
subclass is labelled "sire x dam" (the parents of its first animal in pedigree
order) and an ancestral pair "a x b" in pedigree order. In ebv() the "animal"
levels are the dominance deviations d-hat and the pair levels the h-hat;
k_level_type in the fit and the type column of solutions() tell them apart. An
animal of the pedigree left out of Q has no d-hat; its expected dominance is the
h-hat of its sire x dam pair when that pair is a level (a full sib is in Q, or the
sparse route kept the pair), and is not reported otherwise.
Details
Two exact routes build F^-1. "dense" forms F of the subclasses from the additive
relationships between the parents (Colleau) and inverts it, which costs n^3 / 3 in
the number of subclasses: it wins with few subclasses of many offspring (litters).
"sparse" uses the recursion of the paper on pairs of animals, which is the row of
(I - P) (x) (I - P) with residual variance d_x d_y (2 d_x^2 for the own pair of an
animal), d_x the Mendelian coefficient of A, over the ancestral closure of the
subclasses, and integrates out exactly every pair without animals that has at most one
child: it wins when each subclass has few offspring (dairy). The paper derives the
recursion for non-inbred animals and prunes pairs by a rule it does not prove; the
residual variances under inbreeding and the pruning rule here are this package's
generalization, checked against dominance_matrix() on inbred pedigrees (sire x
daughter, self-fertilization, reciprocal crosses, animals used as sire and as dam,
overlapping generations) and not taken from the paper.
"auto" compares the symbolic Cholesky cost of the joint mixed-model equations of a
reference model, the additive animal effect with A^-1 plus this precision with one
record per animal, and takes the cheaper. The cost is the sum of squared column counts
of the factor, with the sparse columns weighted 8 times the dense final block, which
the package factors in register-blocked tiles; the weight is a fixed constant, so the
route, and with it the levels of the term, depend neither on the machine nor on the
number of threads. The reference model is a proxy: fixed effects, permanent
environment and repeated records are not in it. The dense candidate is costed first
and bounds the work spent on the sparse one, which is dropped for the dense route
(decision says which exit) when its pair closure holds more pairs than the dense
block has entries or more memory than the dense route needs at its peak; when its pair
block of Q would be assembled from more triplets than the full triangle of the dense
block (entries: both routes assemble Q in triplets, and A^-1 and the animal part,
common to both, stay out of the count); or when the minimum-degree ordering of its
equations reads more list entries than the ordering of the dense candidate plus 0.05
times the dense cost (work): an ordering of a candidate with heavy fill is far slower
per operation than the tiled factorization, and the fit would pay that same ordering
before its first iteration. Every exit is a deterministic count, never a timing. The
choice is skipped when the dense route does not fit, which leaves the sparse one, and
when route = is given.
validation/dominance_hv91_routes.R times the choice and both routes (build plus one
evaluation of the mixed-model equations at fixed components) on eight pedigrees of
3,230 to 60,100 animals, from litters to one offspring per pair, with three replicates
except two single runs that its output marks. The decision took 0.02 to 1.9 s. In five
of the seven pedigrees where both routes were timed, the route chosen was the cheaper one
for a fit (build plus ten evaluations): with litters the sparse route cost 2 times more
(litters of 8, four generations) and 140 times more (litters of 6, six generations, 9,280
animals: 4.0 s against 572 s, the sparse route a single run), with two offspring per pair
1.3 times, and with one offspring per pair the dense route cost 11 and 5.6 times more
(3,230 and 15,040 animals; on the latter 233 s against 42 s, the dense route a single
run). On 25,100 animals with 300 records in distinct subclasses the two cost the same
(2.1 and 2.2 s) and "auto" took the sparse one. On 60,100 animals with 800 records
after five generations the ordering of the sparse candidate passed the work ceiling and
"auto" took the dense route, which cost 1.9 times the sparse one (18.1 s against
9.5 s). In the litter pedigree of 12,280 animals the sparse route is not timed: run once
on its own, it was still building Q after 608 s, past 2 GB of memory, when it was
stopped (forcing route = "sparse" on large litters can cost that much); "auto" takes
the dense one there. It is a heuristic over a proxy model, not a guarantee: where the
two routes cost about the same it may take either, and the work ceiling can send a
cheaper sparse route to the dense one, as in the 60,100-animal pedigree.
The dense route needs at most about 18 n_sub^2 + 8 n_parents^2 bytes at its peak
(memory), and is not offered above 4 GB. validation/dominance_hv91_memory.R
measures the rise of a fresh R process above where it stood, three runs each: 72 MB for
the dense route with 2,400 subclasses, against 119 MB estimated; and in two litter
pedigrees where "auto" picks the dense route, 79 against 50 MB on 14,730 animals and
112 against 88 MB on 19,530, so the decision adds 24 to 29 MB and stays under the
estimate of the dense peak (85 and 152 MB).
A fit holds the precision in R (the triplets, 16 bytes per entry), as the penalty of the
term in the engine (12 bytes per entry; the engine's input copy is released once the
equations are set up), and builds the factor and the selective inverse of the
mixed-model equations on top. One evaluation at fixed components with 2,000 subclasses
and 1.8 million entries in Q rose 209 MB, about 116 bytes per entry of Q (the same
script). The 4 GB ceiling above guards only the construction; check the memory before
fitting a precision with tens of millions of entries. The long steps (the closure, the
orderings, the factorizations) check for an interrupt between columns or panels: a time
limit set with setTimeLimit() stops them, which the tests check; Ctrl-C goes through
the same check and was not tested interactively.
WHEN IT DOES NOT EXIST. Delta_i <= 0 leaves nothing for the within-subclass deviation,
and the function stops, naming the animal. With self-fertilization it happens in the
second generation (Delta = -0.125), with full-sib mating in the sixth (parents with
F = 0.59). Two cases. When the subclass has two or more animals in Q, D itself is not
positive-definite on them (two full sibs have D_ij = 1 - Delta, an eigenvalue of
Delta), so no route fits it. When the animal is alone in its subclass the Cockerham D
may still be positive-definite (the selfing example is, smallest eigenvalue 0.5) and
only the subclass representation is missing: if that animal has no record, leave it out
with animals = (its Delta enters nothing else); otherwise use
kernel(id, K = dominance_matrix(ped)).
WHAT IT DOES NOT FIX. The Cockerham D with unit diagonal is the classic approximation
under inbreeding (the full dominance covariance needs further identity components); this
is the inverse of THAT matrix, with the same approximation as the dense route. A sire /
maternal-grandsire pedigree has no sire x dam subclass and is refused; a parent code
with no line of its own (a metafounder or unknown-parent group label) is refused as in
pedigree(), since the Cockerham D is defined on a base of unrelated, non-inbred
founders.
References
Hoeschele, I. & VanRaden, P.M. (1991). Rapid inversion of dominance relationship matrices for noninbred populations by including sire by dam subclass effects. Journal of Dairy Science 74:557-569.
Cockerham, C.C. (1954). An extension of the concept of partitioning hereditary variance for analysis of covariances among relatives when epistasis is present. Genetics 39:859-882.
Examples
ped <- data.frame(animal = as.character(1:8),
sire = c("0", "0", "0", "1", "1", "3", "3", "4"),
dam = c("0", "0", "0", "2", "2", "2", "2", "6"))
Dinv <- dominance_inverse(ped)
table(Dinv$type)
#>
#> animal subclass
#> 8 3
Q <- matrix(0, Dinv$n, Dinv$n)
Q[cbind(Dinv$i, Dinv$j)] <- Dinv$x
Q[cbind(Dinv$j, Dinv$i)] <- Dinv$x
# the animal block of Q^-1 is the dense D
max(abs(solve(Q)[1:8, 1:8] - dominance_matrix(ped)))
#> [1] 5.551115e-17