Fast multi-trait genomic prediction by randomized Gauss-Seidel (PEGS)

The multivariate SNP-BLUP of Xavier and Habier (2022): for trait t, y_t = 1 mu_t + X b_t + e_t, the k effects of each marker jointly N(0, Sb) and the residuals N(0, I s2e_t) UNCORRELATED across traits. Instead of factoring the system of order p k, the k effects of one marker are solved together with the residual updated in place, sweeping the markers in a new random order at each pass; the variances are re-estimated inside the loop by pseudo-expectation, and Sb is bent back to positive-definite (Hayes and Hill, 1981) when it leaves the cone. It fits several traits at once in a fraction of the time of REML, and borrows strength across correlated traits.

Usage

pegs(
  data,
  traits,
  id,
  genotypes,
  environment = NULL,
  cov_structure = c("unstructured", "hcs", "xfa"),
  n_factors = 3L,
  maxiter = 1000L,
  tol = 1e-10,
  deflate_min = 0.75,
  estimate = TRUE,
  start = NULL
)

Arguments

data

data.frame with one row per animal and the traits in columns (traits =), or, with environment =, one row per record of one trait in long form (each animal at most once per environment); NA is a missing record

traits

names of the trait columns (wide form), or the single trait column (long). A factor or text trait is refused rather than read as its codes. There is no missing_code here: a missing record is NA, and a value that looks like a code (-999 held by two records or more, outside the range of the others) or lies more than 10 robust SDs from the median is warned about with the first row, as in model()

id

name of the animal column. A missing id (NA, NaN, Inf, -Inf, blank text or the text "NaN") is refused on every row, as a level column of model(); an id that is present but not among the genotyped animals leaves its record, counted in dropped

genotypes

list(ids, m) with dosages 0/1/2 (NA imputed by the marker mean); a repeated id, rownames that contradict ids or an animal below a 0.90 call rate are refused, as in model()

environment

name of the environment column for the long form: each level becomes a trait. A missing environment is refused on every row, by the same rule as id: it would become one more trait, or its records would vanish uncounted

cov_structure

"unstructured", "hcs" or "xfa"

n_factors

latent factors of "xfa"

maxiter

passes over the markers

tol

convergence on log10 of the sum of squared changes of the effects in a pass

deflate_min

smallest factor the bending may apply to the covariances

estimate

FALSE keeps start fixed: a pure multivariate ridge solve

start

list(Vb =, Ve =): the marker-effect covariance (k x k) and the residual variances; NULL starts from half the phenotypic variance

Value

a breeding_pegs object: mu, marker_effects (markers x traits), gebv (genotyped animals x traits), h2, Vb, Ve, Gcor (genetic correlations), bend, iters, converged, n_cov_params, n_records, dropped (records without genotype) and n_imputed.

Details

What it assumes, and where that fails. The residuals of different traits are independent: it is a method for the same trait in different ENVIRONMENTS (farms, locations, years), where each animal is recorded in one or a few of them. For traits recorded on the SAME animal and occasion the residual correlation has nowhere to go and leaks into the genetic covariance: on real data of number born and litter weight (the same litter) PEGS gave r_g = 0.83 against 0.33 +- 0.08 from the bivariate REML with the residual matrix estimated. Use model_mt() there. The pseudo-expectation variances are noisier than REML (standard errors 60 to 100 percent larger in the original paper), and so are the genetic correlations when the markers far outnumber the animals: in four simulated replicates (3 traits, r_g = 0.6, h2 = 0.5) the largest error of r_g ranged from 0.05 to 0.26 with 1000 animals and 5000 markers, and from 0.03 to 0.08 with 3000 animals, while the GEBV accuracy stayed at 0.73 to 0.80. The GEBV are where the method is strong.

cov_structure restricts Sb (Xavier et al., 2025): "hcs" is one common genetic correlation with a variance per trait (k + 1 parameters), "xfa" is n_factors latent factors with the diagonal kept, for many environments sparsely tested.

References

Xavier, A. & Habier, D. (2022). A new approach fits multivariate genomic prediction models efficiently. Genetics Selection Evolution 54:45.

Hayes, J.F. & Hill, W.G. (1981). Modification of estimates of parameters in the construction of genetic selection indices ('bending'). Biometrics 37:483-493.