Efficient strategies for leave-one-out cross validation for genomic best linear unbiased prediction

1A random multiple-regression model that simultaneously fit all allele substitution effects for additive markers or haplotypes as uncorrelated random effects was proposed for Best Linear Unbiased Prediction (BLUP), using whole-genome data. Breeding values are defined as the sum of the effects of all the markers or haplotypes, and their estimates are widely used for prediction of the merit of selection candidates. Estimates of marker or haplotype effects are used to predict breeding values of individuals that were not present in a previous analysis commonly referred to as training. An alternative earlier published approach to use marker or haplotype information fits breeding values as random effects based on covariances defined by a “genomic relationship matrix” computed from genotypes. These two models have been shown to be equivalent in terms of predicting breeding values and we refer to them here as marker effect models (MEM) or breeding value models (BVM), the latter often known as Genomic Best Linear Unbiased Prediction (GBLUP).

2Cross validation is often used to quantify the predictive ability of a statistical model. In k-fold cross validation, the whole dataset is partitioned into k parts with k analyses, where one part is omitted for training with validation on the omitted part. Leave-one-out cross validation (LOOCV) is a special case of k-fold cross validation with k=n, the number of observations. When the dataset is small, leave-one-out cross validation is appealing as the size of the training set is maximized. However, naive application of LOOCV is computationally intensive, requiring n analyses.

3We show below how LOOCV can be performed using either the MEM or BVM with little more effort than is required for a single analysis with n observations.

4Use of the MEM is more efficient when the number n of individuals is larger than the number p of markers, because for this model the mixed model equations are of order p plus the number of other effects. When n<p, estimated breeding values can be obtained more efficiently by solving the mixed model equations for the BVM of order n plus the number of other effects. We deal with the special case where the only other effect is a general mean and phenotypes have been pre-corrected for other nuisance variables. Efficient strategies for LOOCV using this special case for MEM when n≥p and BVM when p≥n are shown below.

5The MEM for GBLUP can be written as

6where y, a n×1 vector for phenotypes, has been pre-corrected for all fixed effects other than μ, the overall mean, X is the n×p matrix of marker covariates, is a p×1 random vector of the allele substitution effects and e is a n×1 random vector of residuals. Often it is assumed that marker effects are identically and independently distributed (iid) random variables with null means and variances. Thus, under the usual assumption that the residuals are iid with null means and variances, E(y)=1 μ. When MEM is used, LOOCV can be performed by using a well-known strategy used in least-squares regression to compute the predicted residual sum of square (PRESS) statistic.

7BLUP of can be obtained by solving the mixed model equations

8where X ∗= [ 1 X],, D is a diagonal matrix whose elements are 0 followed by a p vector of 1s and.

9Now, BLUP for, where observation j is left out, can be obtained as

10where is X ∗ with the jth row removed and y −j is y with the jth element removed.

11Suppose is the jth row of X ∗, then from the matrix inverse lemma,

12where the quadratic is the jth diagonal element of H=X ∗(X ∗′ X ∗+D λ)−1 X ∗′.

13Using (3) in (4), the prediction residual for the jth observation can be written as

14These prediction errors can be squared and accumulated over n realizations to compute PRESS defined as. The accuracy of genomic prediction is often quantified as the correlation between the predicted and observed values of y j, and that correlation can be estimated from the values of, which can be computed efficiently as, using the observed values of y j. When a specific group of individuals is of interest, prediction accuracies and PRESS can also be calculate using for individuals in that group.

15When n<p, the genomic prediction of the breeding value can be obtained more efficiently by solving the mixed model equations for the BVM:

16where u=X β,, Z is the identity matrix of order n and other variables are as in the MEM. Further, in both models E(y)=1 μ, and. These two models are said to be equivalent, and linear functions predicted from one model are identical to corresponding predictions from the other model. Two efficient strategies for LOOCV using the BVM are shown below.

17The mixed model equations for this model are:

18where. Due to the relative order of the coefficient matrices for the MEM and the BVM, when n<p, is more efficiently obtained as. Similarly, can be obtained more efficiently as. Using these two equalities, the formula for becomes:

19where the quadratic is the jth diagonal element of C=Z ∗(Z ∗′ Z ∗+G λ)−1 Z ∗′.

20Another efficient strategy for BVM is shown here. First we consider the situation where y has been pre-corrected for μ in addition to nuisance effects so that E(y)=0 and we define. Now matrix Q is constructed by augmenting the covariance matrix of y with one leading row and column as

21To obtain the prediction error for observation j, the second row and column of Q are permuted with row and column j+1. In this manner Q has its rows and columns symmetrically permuted as P j′Q P j =W, where the permutation matrix P j is obtained by permuting the second row of the n order identity matrix with row j+1. So the permuted matrix is:

22where we will define the leading 2×2 matrix as, and the other partitions as, and C=V −j,−j, where −j denotes that the jth element, row or column has been removed. Defining W 11 as the top left or leading 2×2 sub-matrix in W −1 corresponding to the position of A in W, and using partitioned inverse-matrix identities, the inverse of W 11 can be written as,

23Now V j,−j in element (2,1) of the above inverse matrix is the vector of covariances between y j and y −j and is the inverse of the covariance matrix of y −j. Thus, is the Best Linear Predictor (BLP) of y j given y −j, and element (2,1) of (10) is the prediction error of y j. The element (2,2) in (10) is the prediction error variance (PEV) for y j, where. PEV can also be used to calculate theoretical reliability for individual i as, and characterizing the distributions of reliability for all the individuals in a dataset has a number of practical applications. Note this allows us to obtain the PEV of every individual and the distribution of these values provide information as to the robustness of genomic predictions across the population of individuals represented in the dataset. This PEV is determined by the genomic variance-covariance matrix and does not depend on y. Two different datasets could generate the same PRESS statistic but with different distributions of PEV.

24Now, because the permutation matrix P j is orthogonal,, and the elements of W 11 that are of interest in terms of predicting individual j can be obtained directly from Q −1 as

25It follows that, which is the off-diagonal element of the inverse of the 2×2 matrix W 11, can be written in terms of Q −1 as

26where q i,j is the element from row i and column j of Q −1. Thus, once Q −1 is computed, for all j can be computed using (12), and these values can be used to compute PRESS as. To estimate the correlation between the predicted and observed values of y j, the value of is efficiently obtained as the difference.

27Now we consider the situation without pre-correcting y for μ, where E(y)=1 μ. Now the mixed model (7) contains both fixed and random effects. Note that the mixed model equations that correspond to this mixed effects model can be derived by treating μ as “random” with null mean and large variance. So, let

28for sufficiently large value of. Then under this assumption, E(y)=0 and, and thus is the BLP from the random effects rather than mixed effects model of y j given y −j. This BLP obtained from the model with random μ will be numerically very close to the BLUP obtained from the mixed model with fixed μ. The Q matrix corresponding to the BLP with random μ is constructed as and prediction residuals are obtained as (12).

29Phenotypes y and genotypes X at 5 markers for 3 individuals are in Table 1. Assume and the overall mean μ is the only fixed effect. In LOOCV strategy for MEM and strategy I for BVM, the diagonal elements of H for MEM and C for BVM, which are in the denominators of (6) and (9), are in Table 2. The numerators of (6) and (9) are obtained by solving the MME (2) and (8). Then prediction errors are calculated as in (6) and (9) and shown in Table 4. In LOOCV strategy II for BVM, the Q matrix (Table 3) is constructed using, which is sufficiently large relative to for μ to be indistinguishable from a fixed effect with a flat prior. The prediction errors are calculated as (12) and shown in Table 4. The MEM strategy and BVM strategy I gave identical prediction errors and identical PRESS for this numerical example were numerically very close to those from the BVM strategy II.

30Two datasets were simulated using XSim, where 1,000 offspring were sampled from random mating of 100 parents for 10 non-overlapping generations, to compare the computational efficiencies for naive and efficient strategies using BVM or MEM for LOOCV in GBLUP. Dataset I was simulated with 1,000 observations and 10,000 SNP markers for a p≫n scenario. Dataset II was simulated with 1,000 observations and 100 markers for a n≫p scenario. The processor used in the analyses was a 1.4 GHz Intel Core i5 with 4 GB of memory.

31For dataset II, efficient MEM is 99 times faster than the naive application (2.979 s versus 0.030 s) (Table 5). All strategies implemented in Julia, a scientific programming language, gave virtually identical prediction accuracies defined as the correlation between y and for each dataset. For dataset I, efficient BVM is 786 times faster than the naive application (3.107 s versus 2,442.59 s) (Table 5).

32In genomic prediction, the candidates to be predicted are often offspring that are genotyped but not yet phenotyped. In this situation, LOOCV using all individuals in the training dataset will provide an upper bound for the accuracy of prediction, because ancestors in the training dataset with large numbers of descendants have more accurate predictions than descendants. A better estimate of the accuracy of prediction can be obtained by applying LOOCV to only terminal offspring in the training dataset.

33Efficient strategies for LOOCV in GBLUP are presented in this paper. LOOCV strategy I and II for BVM are more efficient when p≫n. LOOCV strategy for MEM is more efficient when n≫p. The accuracy of genomic prediction is often quantified as the correlation between the predicted and observed values of y j, and this correlation can be estimated efficiently using LOOCV strategies. Compared to naive application of LOOCV, which is computationally intensive, LOOCV can be implemented efficiently.