Choose the APY core by the eigenvalues of G

The core size is the number of eigenvalues of the genomic relationship matrix that explain variance of its trace, the rule of Pocrnic et al. (2016a): the dimensionality of the genomic information, which is what bounds how many animals the APY needs to condition on. The animals themselves are a random draw of that size, the standard practice (Pocrnic et al. 2016a,b; Bradford et al. 2017 found the definitions equivalent once the core is large enough). This is the count that OPTION snp_svd prints in preGSf90, followed by a random core. model(apy_core = "auto") and the other fitters call it with the defaults.

Usage

apy_core_select(
  genotypes,
  variance = 0.98,
  size = NULL,
  include = NULL,
  seed = 1L,
  method = c("auto", "exact", "lanczos"),
  probes = 30L,
  steps = 100L
)

Arguments

genotypes

list with ids and m (animals x markers, 0/1/2 coded); NA is imputed with the marker mean and monomorphic markers are left out, as in g_matrix()

variance

fraction of the trace of G the core must explain, strictly between 0 and 1; 0.98 is the value at which Pocrnic et al. (2016a) found realized accuracy to peak, 0.99 the conservative choice

size

a core size to use instead of the count (the eigenvalue table is still reported); NULL uses the count

include

ids that must be in the core, for example the proven sires (Masuda et al. 2016 found a small gain from sires and cows over a random core). The rest of the core is drawn from the other genotyped animals. If include is longer than the size, the core is include and a warning says so, never a silent truncation

seed

seed of the random draw (and of the Lanczos probes); the global random stream is restored afterwards, so a later gibbs() chain or simulation is not affected

method

"exact" (eigenvalues of the smaller Gram matrix), "lanczos" (stochastic Lanczos quadrature, no Gram matrix) or "auto" (exact up to min(n_animals, n_markers) = 4000)

probes, steps

random sign vectors and Lanczos steps of the "lanczos" route; the cost is 2 x steps passes over the genotype matrix with probes columns, plus the averaged rule, which solves ceiling(steps / 4) dense eigenproblems of order up to steps per probe in R, a cost that grows with probes x steps^4: with 30 probes 0.5 to 0.8 s at 100 steps, 7 to 12 s at 200 and 27 to 49 s at 300 (three runs in each of two sessions, validation/apy_core_lanczos_se.R)

Value

character vector with the ids of the core animals, in genotype order, of class breeding_apy_core, with attributes size, variance (the target), variance_explained (the fraction the chosen size reaches), eig (the counts at 90, 95, 98 and 99 percent, the preGSf90 table), eigenvalues (of G, decreasing), n_genotyped, n_markers, seed, method, and for "lanczos" count_se (the standard error of the count at variance over the random probes), probes and steps (eigenvalues is then NULL), and cost_ratio, the flops of one sparse factorization of the APY block relative to the exact dense one, (nj (nc + 1)^2 + (nc + 1)^3 / 3) / (n^3 / 3) with nc core and nj non-core animals

Details

The count is made on the raw VanRaden G, as in the source, not on the G* that the single-step then inverts (G scaled to A22 and blended with it): the A22 part has full rank and would inflate the count. The consequence is that the blended part is only partly covered by any core smaller than the genotyped set, so the APY breeding values stay a little away from the exact ones even when the core explains all of G. Pocrnic et al. (2016a) report correlations with the exact GEBV above 0.98 with the 98% core and above 0.99 with 99%; Cesarani et al. (2023), with a few thousand genotyped, found accuracy still growing up to 99%, which makes variance = 0.99 the conservative choice in small populations.

With method = "exact" the eigenvalues come from the smaller of the two Gram matrices of the centered genotypes (animals x animals or markers x markers), so the cost is cubic in min(n_animals, n_markers): about half a minute at 3000 with R's reference LAPACK, 12 minutes at 10 000, and at the 30 000 x 40 000 of a routine evaluation a 7 GB matrix before the decomposition starts. method = "lanczos" never forms it: stochastic Lanczos quadrature (Ubaru, Chen & Saad, 2017) runs steps Lanczos steps on G from each of probes random sign vectors, using only products with the genotype matrix (two passes per step, all probes at once, on br_threads() threads), and the Gauss nodes and weights of the tridiagonals estimate the spectral measure of G, where the count is read at the point its mass reaches variance. A single Gauss rule per probe is biased there by an amount that more probes do not reduce, because the Ritz values near the threshold fall in almost the same place for every probe; the size and the sign of that bias change with steps and variance. Each probe therefore contributes the average of the Gauss rules of its last ceiling(steps / 4) leading tridiagonals, whose nodes interlace, at no extra pass over the genotypes; a probe that stops early at an invariant subspace keeps its whole tridiagonal, whose rule is already exact. The count is an estimate; count_se is its standard error over the random probes, from the linearized count. On 1000 animals with 40 steps, eight independent draws of 600 probes, each split into 20 sets of 30, gave a spread of the 20 counts of 0.80 to 1.31 times the mean count_se at 90 to 99% (mean 1.01; validation/apy_core_lanczos_se.R). On 10 000 animals and 10 000 markers (validation/apy_core_lanczos.R, ten seeds) the exact count at 98% is 5953 and the estimates ran from 5934 to 5975, mean 5957.3 with a standard error of about 12; with the same probes on the exact eigenvectors the quadrature error was +2.3 counts (0.04%), against +21.4 (0.36%) for a single Gauss rule. "auto" takes the exact route up to min(n_animals, n_markers) = 4000 and Lanczos above. Choose the core once, keep the result and pass it to every fit, rather than "auto" in each of them.

References

Pocrnic, I., Lourenco, D.A.L., Masuda, Y., Legarra, A. & Misztal, I. (2016a). The dimensionality of genomic information and its effect on genomic prediction. Genetics 203:573-581.

Pocrnic, I., Lourenco, D.A.L., Masuda, Y. & Misztal, I. (2016b). Dimensionality of genomic information and performance of the Algorithm for Proven and Young for different livestock species. Genetics Selection Evolution 48:82.

Bradford, H.L. et al. (2017). Journal of Animal Breeding and Genetics 134:545-552.

Masuda, Y. et al. (2016). Journal of Dairy Science 99:1968-1974.

Cesarani, A. et al. (2023). Animal 17:100766.

Ubaru, S., Chen, J. & Saad, Y. (2017). Fast estimation of tr(f(A)) via stochastic Lanczos quadrature. SIAM Journal on Matrix Analysis and Applications 38:1075-1099.

Examples

set.seed(2)
m <- matrix(rbinom(60 * 400, 2, 0.4), 60, 400)
nucleo <- apy_core_select(list(ids = paste0("a", 1:60), m = m))
#> Warning: a core of 57 of 60 genotyped: one factorization with APY costs 104% of the exact one, a small gain for an approximation
nucleo
#> APY core: 57 of 60 genotyped, explaining 98.6% of G (400 markers, seed 1)
#> eigenvalues of G for 90/95/98/99% of its trace: 48 / 53 / 57 / 58 
#> one factorization with this core costs 104.3% of the exact one