Cross-Validation Without Doing Cross-Validation in Genome-Enabled Prediction

Cross-validation of methods is an essential component of genome-enabled prediction of complex traits. We develop formulae for computing the predictions that would be obtained when one or several cases are removed in the training process, to become members of testing sets, but by running the model using all observations only once. Prediction methods to which the developments apply include least squares, best linear unbiased prediction (BLUP) of markers, or genomic BLUP, reproducing kernels Hilbert spaces regression with single or multiple kernel matrices, and any member of a suite of linear regression methods known as “Bayesian alphabet.” The approach used for Bayesian models is based on importance sampling of posterior draws. Proof of concept is provided by applying the formulae to a wheat data set representing 599 inbred lines genotyped for 1279 markers, and the target trait was grain yield. The data set was used to evaluate predictive mean-squared error, impact of alternative layouts on maximum likelihood estimates of regularization parameters, model complexity, and residual degrees of freedom stemming from various strengths of regularization, as well as two forms of importance sampling. Our results will facilitate carrying out extensive cross-validation without model retraining for most machines employed in genome-assisted prediction of quantitative traits.

overall results are combined; this is known as a K-fold CV (Hastie et al. 2009). A loss function, such as mean-squared error (MSE) or predictive correlation is computed to gauge the various predictive machines compared. However, the process must be repeated a number of times, with folds reconstructed at random (whenever possible) to obtain measures of CV uncertainty (e.g., Okut et al. 2011;Tusell et al. 2014). CV is computationally taxing, especially when Bayesian prediction models with a massive number of genomic covariates and implemented via Markov chain Monte Carlo (MCMC) are involved in the comparison.
Stylized formulae (e.g., Daetwyler et al. 2008) suggest that the expected predictive correlation ("accuracy") in genome-enabled prediction is proportional to training sample size (n). On intuitive grounds, more genetic variability ought to be spanned as a training sample grows, unless additional cases bring redundant information. With larger n, it is more likely that genomic patterns appearing in a testing set are encountered in model training. Although the formulae do not always fit real data well (Chesnais et al. 2016), it has been observed that a larger n tends to be associated with larger predictive correlations (Utz et al. 2000;Erbe et al. 2010).
Arguably, there is no better expectation than what is provided by a CV conducted under environmental circumstances similar to those under which the prediction machine is going to be applied. When n is small, the largest possible training set sample size is attained in a leaveone-out (LOO) CV, e.g., Ober et al. (2015) with about 200 lines of Drosophila melanogaster. In LOO CV, n 2 1 cases are used for model training, to then predict the single out-of-sample case. Model training involves n implementations, each consisting of a training sample of size n 2 1 and a testing set of size 1.
It is not widely recognized that it is feasible to obtain CV results by running the model only once, which is well known for least-squares regression (e.g., Seber and Lee 2003). Here, we show that this idea extends to other prediction machines, such as ridge regression (Hoerl and Kennard 1970), genome-based best linear unbiased prediction (GBLUP; Van Raden 2008), and reproducing kernel Hilbert spaces regression (RKHS; Gianola et al. 2006;Gianola and van Kaam 2008). It is also shown that the concept can be applied in a MCMC context to any Bayesian hierarchical model, e.g., members of the "Bayesian alphabet" (Meuwissen et al. 2001;Gianola et al. 2009;Gianola 2013). This manuscript reviews available results for least-squares based CV, and shows how CV without actually doing CV can be conducted for ridge regression, BLUP of marker effects, GBLUP, and RKHS for any given kernel matrix. It is described how importance sampling can be used to produce Bayesian CV by running MCMC only once, which has great advantage in view of the intensiveness of MCMC computations. Illustrations are given by using a well characterized data set containing wheat grain yield as phenotype and 1279 binary markers as regressors, and the paper concludes with a Discussion. Most technical results are presented in a series of Appendices, to facilitate reading the main body of the manuscript.

CROSS-VALIDATION WITH ORDINARY LEAST-SQUARES Setting
A linear model used for regressing phenotypes on additive codes at biallelic marker loci (21, 0, 1 for aa, Aa and AA, respectively), such as in a genome-wide association study, is y ¼ Xbþe: (1) Here, y is the n · 1 vector of phenotypic measurements, X ¼ fx ij g is the n · p marker matrix, and x ij is the number of copies of a reference allele at locus j observed in individual i; b is a p · 1 vector of fixed regressions on marker codes, known as allelic substitution effects. Phenotypes and markers are often centered; if an intercept is fitted, the model is expanded by adding b 0 as an effect common to all phenotypes. The residual vector is assumed to follow the distribution e $ Nð0; Is 2 e Þ; where I is an n · n identity matrix, and s 2 e is a variance parameter common to all residuals.
The basic principles set here carry to the other prediction methods discussed in this paper. In this section, we assume rankðXÞ ¼ p , n so that ordinary least-squares (OLS) or maximum likelihood under the normality assumption above can be used. The OLS estimator of b iŝ b ¼ ðX9XÞ 21 X9y, and the fitted residual for datum i isê i ¼ y i 2 x i 9b, where x i 9 is the i th row of X: Assuming the model holds, EðbÞ ¼ b, so the estimator is unbiased. A review of the pertinent principles is given in Appendix A, from which we extract results.
It is shown in Appendix A that the uncertainty of a prediction, as measured by variance, increases with p (model complexity), and decreases with the size of the testing set, n test . Two crucial matters in genome-enabled prediction must be underlined. First, if the model is unnecessarily complex, prediction accuracy (in the MSE sense) degrades unless the increase in variance is compensated by a reduction in prediction bias. Second, if the training set is made large at the expense of the size of the testing set, prediction mean squared error will be larger than otherwise. The formulae of Daetwyler et al. (2008) suggest that expected prediction accuracy, as measured by predictive correlation (not necessarily a good metric; González-Recio et al. 2014), increases with n. However, the variability of the predictions would increase, as found by Erbe et al. (2010) in an empirical study of Holstein progeny tests with alternative CV layouts. Should one aim at a higher expected predictive correlation or at a more stable set of predictions at the expense of the former? This question does not have a simple answer.
Leave-one-out (LOO) cross-validation LOO is often used when n is small and there is concern about the limited size of the training folds. Let X ½2i be X with its i th row ðx i 9Þ removed, so that its order is ðn 2 1Þ · p: Since are p · p matrices. Likewise, if y ½2i is y with its i th element removed, the OLS right-hand sides in LOO are Making use of (81) in Appendix B, the least-squares estimator of b formed with the i th observation deleted from the model is expressible asb Employing Appendix C, the estimator can be written in the form where h ii ¼ x i 9ðX9XÞ 21 x i andê i ¼ y i 2 x i 9b is the fitted residual using all n observations in the analysis; the fitted LOO residual is Hence, the LOO estimator and prediction error can be computed directly from the analysis carried out with the entire data set: no need for n implementations. Making use of (6), the realized LOO CV mean squared error of prediction is and the expected mean squared error of prediction is given by where D ¼ fð12h ii Þ 22 g is an n · n diagonal matrix. As shown in Appendix A the LOO expected PMSE gives an upper bound for the expected squared error in least-squares based CV. The extent of overstatement of the error depends on the marker matrix X (via the h9s) and on the prediction biases d i : Hence, LOO CV represents a conservative approach, with the larger variance of the prediction resulting from the smallest possible testing set size ðn test ¼ 1Þ If the prediction is unbiased, the d 2 terms vanish, and it is clear that observations with h 2 values closer to 1 contribute more to squared prediction error than those with smaller values, as the model is close to overfitting the former type of observations.

Leave-d-out cross-validation
The preceding analysis suggests that reallocation of observations from training into testing sets is expected to reduce PMSE relative to the LOO scheme. Most prediction-oriented analyses use K 2 fold CV, where K is chosen arbitrarily (e.g., K ¼ 2; 5; 10) as mentioned earlier; the decision of the number of folds is usually guided by the number of samples available. Here, we address this type of scheme generically by removing d out of the n observations available for training, and declaring the former as members of the testing set. Let X ½2d be X with d of its rows removed, and y ½2d be the data vector without the d corresponding data points. As shown in Appendix D, where H d ¼ X ½d ðX9XÞ 21 X9 ½d ; note that ðI2H d Þ 21 does not always exist but can be replaced by a generalized inverse, andb ½2d will be invariant to the latter if p , n 2 d. Predictions are invariant with respect to the generalized inverse used. The similarity with (5) is clear: ðI2H d Þ 21 appears in lieu of 1 2 h ii ; X9 ½d of order p · d contains (in columns) the rows of X being removed, andê ½d ¼ y ½d 2 X ½db are the residuals corresponding to the left out d 2 plet obtained when fitting the model to the entire data set. The error of prediction of the d phenotypes entering into the testing set is The mean-squared error of prediction of the d observations left out (testing set) becomes The prediction bias obtained by averaging over all possible data sets is where m d is a d · 1 vector of true means of the distribution of observations in the testing sets. After algebra, and E y ½d ;yj;X;X ½d ½PMSEðdÞ ¼ Observe that the term in brackets is a matrix counterpart of (76) in Appendix A, with H d playing the role of h ii in the expression. The two terms in the equation above represent the contributions and bias and (co) variance to expected squared prediction error. The next section illustrates how the preceding logic carries to regression models with shrinkage of estimated allelic substitution effects ðbÞ:

CROSS-VALIDATION WITH SHRINKAGE OF REGRESSION COEFFICIENTS
BLUP of markers (ridge regression) Assume again that phenotypes and markers are centered. Marker effects b are now treated as random variables and assigned the normal Nð0; Is 2 b Þ distribution, where s 2 b is a variance component. The BLUP of b (Henderson 1975) is given by is a shrinkage factor taken as known. BLUP has the same mathematical form as the ridge regression estimator (Hoerl and Kennard 1970), developed mainly for tempering problems caused by colinearity among columns of X in regression models where p , n; and with all regression coefficients likelihood-identified. The solution vector b r can also be assigned a Bayesian interpretation as a posterior expectation in a linear model with Gaussian residuals, and Nð0; Is 2 b Þ used as prior distribution, with variance components known (Dempfle 1977;Gianola and Fernando 1986). A fourth view of b r is as a penalized maximum likelihood estimator under an L 2 penalty (Hastie et al. 2009). Irrespective of its interpretation, b r provides a "point statistic" of b for the n , p situation. In BLUP, or in Bayesian inference, it is not a requirement that the regression coefficients are likelihood identified. There is one formula with four interpretations (Robinson 1991).
Given a testing set with marker genotype matrix X test ; the point prediction of yet to be observed phenotypes is X test b r : We consider LOO CV because subsequent developments assume that removal of a single case has a mild effect on s 2 e and s 2 b . This assumption is reasonable for animal and plant breeding data sets where n is large, so removing a single observation should have a minor impact on, say, maximum likelihood estimates of variance components. If l is kept constant, it is shown in Appendix E that is the residual from ridge regression BLUP applied to the entire sample. A similar expression for leave-d-out cross-validation using the same set of variance components is also in Appendix E. If d=n is smaller and n is reasonably large, the error resulting from using variance components estimated from the entire data set should be small.
The error of predicting phenotype i is now y i 2 x9 i b r ½2i and is expressible as The first term is the average squared prediction bias, and the second is the prediction error variance. As s 2 b /0, (18) tends to the leastsquares PMSEð1Þ:

Genomic BLUP
Once marker effects are estimated as b r , a representation of genomic BLUP (GBLUP) for n individuals is the n · 1 vector g ¼ Xb r with its i th element beingĝ i ¼ x9 i b r . In GBLUP, "genomic relationship matrices" are taken as proportional to XX9 (where X often has centered columns); various genomic relationship matrices are in, e.g., Van Raden (2008), Astle and Balding (2009) and Rincent et al. (2014). Using (16) LOO GBLUP (i.e., excluding case i from the training sample) iŝ The formula above requires finding b r ; given l; the procedure entails solving p equations on p unknowns and finding the inverse of C is impossible or extremely taxing when p is large. A simpler alternative based on the well-known equivalence between BLUP of marker effects and of additive genotypic value is used here. If g ¼ Xb is a vector of marked additive genetic values and b $ Nð0; Is 2 b Þ; then g $ Nð0; XX9s 2 b Þ: Many genomic relationship matrices are expressible as G ¼ XX9 c for some constant c; so that g $ Nð0; Gs 2 g Þ and s 2 g ¼ cs 2 b is called "genomic variance" or "marked additive genetic variance" if X encodes additive effects; clearly, there is no loss of generality if c ¼ 1 is used, thus preserving the l employed for BLUP of marker effects. The model for the "signal" g becomes Letting C ¼ I þ G 21 l, then BLUPðgÞ ¼ĝ ¼ fĝ i g ¼ C 21 y is GBLUP using all data points. Appendix F shows how LOO GBLUP and d-out GBLUP can be calculated indirectly from elements or blocks of C; and elements of y:

RKHS regression
In RKHS regression (Gianola et al. 2006;Gianola and van Kaam 2008), input variables, e.g., marker codes, can be transformed nonlinearly, potentially capturing both additive and nonadditive genetic effects (Gianola et al. 2014a(Gianola et al. , 2014b, as further expounded by Jiang and Reif (2015) and Martini et al. (2016). When a pedigree or a genomic relationship matrix is used as kernel, RKHS yields pedigree-BLUP and GBLUP, respectively, as special cases (Gianola and de los Campos 2008;de los Campos et al. 2009de los Campos et al. , 2010. The standard RKHS model is with g ¼ Ka (and therefore g i ¼ k i 9a); K is an n · n positive (semi)definite symmetric matrix so that K ¼ K9; a ¼ K 21 g when the inverse is unique, and a $ Nð0; K 21 s 2 a Þ: BLUPðaÞ can be obtained by solving the system with solution (since K9K ¼ K 2 and K is invertible) The BLUP of g under a RKHS model is where K stands for "kernel," and l K ¼ s 2 e s 2 a : Putting C 21 K ¼ ðI þ K 21 l RKHS Þ 21 ; the RKHS solutionĝ K ¼ C 21 K y has the same form as BLUPðgÞ ¼ĝ ¼ C 21 y; as given in the preceding section. Using Appendix F, it follows that and PMSE K ðdÞ ¼ẽ 9 d;Kẽd;K d : The previous expressions reduce to the LOO CV situation by setting d ¼ 1:

BAYESIAN CROSS-VALIDATION Setting
Many Bayesian linear regression on markers models have been proposed for genome-assisted prediction of quantitative traits (e.g., Meuwissen et al. 2001;Heslot et al. 2012;de los Campos et al. 2013;Gianola 2013). All such models pose the same specification for the Bayesian sampling model (a linear regression), but differ in the prior distribution assigned to allelic substitution effects. Implementation is often via MCMC, where computations are intensive even in the absence of CV; shortcuts and approximations are not without pitfalls. Is it possible to do CV by running an MCMC implementation only once? What follows applies both to LOO and d 2 out CV situations as well as to any member of the Bayesian alphabet Gianola 2013) Suppose some Bayesian model has been run with MCMC, leading to S samples collected from a distribution with posterior density pðujy; HÞ; here, u are all unknowns to be inferred and H denotes hyper-parameters arrived at in a typically subjective manner, e.g., arbitrary variances in a four-component mixture distribution assigned to substitution effects (MacLeod et al. 2014). In CV, the data set is partitioned into y ¼ ðy test ; y train Þ, training and testing sets are chosen according to the problem in question, and Bayesian learning is based on the posterior distribution ½ujy train ; H. Predictions are derived from the predictive distribution of the testing set data the preceding assumes that y test is independent of y train given u , a standard assumption in genome-enabled prediction. The point predictor chosen most often is the expected value of the predictive distribution In the context of sampling, representation (29) (29). We describe Bayesian LOO CV, but extension to a testing set of size d is straightforward. In LOO, the data set is partitioned as y ¼ ðy i ; y 2i Þ; i ¼ 1; 2; . . . ; n; where y i is the predictand and y 2i is the vector containing all other phenotypes in the data set. A brute force process involves running the Bayesian model n times, producing the posterior distributions ½ujy 2i ; H; i ¼ 1; 2; . . . ; n: Since LOO CV is computationally formidable in an MCMC context, procedures based on drawing samples from ½ujy; H and converting these into realizations from ½ujy 2i ; H can be useful (e.g., Gelfand et al. 1992;Gelfand 1996;Vehtari and Lampinen 2002). Use of importance sampling, and of sampling importance resampling (SIR), algorithms for this purpose is discussed next. Cantet et al. (1992) and Matos et al. (1993) present early applications of importance sampling to animal breeding.

Importance sampling
We seek to estimate the mean of the predictive distribution of the left-out data point Eðy i y 2i ; HÞ: Since e i $ Nð0; s 2 e Þ is independent of y 2i , one has Here, w i ðbÞ ¼ pðbjy 2i ;HÞ pðbjy;HÞ is called an "importance sampling" weight Albert 2009). Expression (32) implies that the posterior mean of b in a training sample can be expressed as the ratio of the posterior means of wðbÞb, and of wðbÞ taken under a Bayesian run using the entire data set. It is shown in Appendix F that, given draws b ðsÞ ; s 2ðsÞ e ðs ¼ 1; 2; . . . ; SÞ from the full-posterior distribution, the posterior expectation can in equation (32) be estimated aŝ where By making reference to (31), it turns out that a Monte Carlo estimate of the mean of the predictive distribution of datum i in the Bayesian LOO CV is given bŷ This type of estimator holds for any Bayesian linear regression model irrespective of the prior adopted, i.e., it is valid for any member of the "Bayesian alphabet" Gianola 2013). In d 2 out CV, the prediction iŝ where where for j being a member of the d 2 plet of observations forming the testing set. The importance sampling weights are the reciprocal of conditional likelihoods; this specific mathematical representation can produce imprecise estimates of posterior expectations, especially if the posterior distribution with all data has much thinner tails than the posterior based on the training set. Vehtari and Lampinen (2002) calculate the "effective sample size" for a LOO CV as If all weights are equal over samples, the weight assigned to any draw is S 21 ; and the variance of the weights is 0; yielding S eff ;i ¼ S; on the other hand, S eff ;i can be much smaller than S if the variance among weights is large, e.g., when some weights are much larger than others. The SIR algorithm described by Rubin (1988), Smith and Gelfand (1992) and Albert (2009) can be used to supplement importance sampling; SIR can be viewed as a weighted bootstrap. Let the sampled values and the (normalized) importance sampling weights be b ðsÞ and w i;s ; respectively, for i ¼ 1; 2; . . . ; n and s ¼ 1; 2; . . . ; S: Then, obtain a resample of size S by sampling with replacement over b ð1Þ ; b ð2Þ ; . . . ; b ðSÞ with unequal probabilities proportional to w i;1 ; w i;2 ; . . . ; w i;S ; respectively, obtaining realizations b ð1Þ rep ; b ð2Þ rep ; . . . ; b ðSÞ rep : Finally, average realizations for estimating Eðbjy 2i ; HÞ in (31).
The special case of "Bayesian GBLUP" The term "Bayesian GBLUP" is unfortunate but has become entrenched in animal and plant breeding. It refers to a linear model that exploits genetic or genomic similarity matrix among individuals (as in GBLUP), but where the two variance components are unknown and learned in a Bayesian manner. Prior distributions typically assigned to variances are scale inverted chi-square processes with known scale and degrees of freedom parameters (e.g., Pérez and de los Campos 2014). The model is y ¼ g þ e; with e $ Nð0; Is 2 e Þ; the hierarchical prior is gjs 2 g $ Nð0; Gs 2 g Þ; s 2 g S 2 g ; n g and s 2 e S 2 e ; n e , where the hyperparameters are H ¼ ðS 2 g ; n g ; S 2 e ; n e Þ. A Gibbs sampler iteratively loops over the conditional distributions Above, ELSE denotes the data plus H and all other parameters other than the ones being sampled in a specific conditional posterior distribution; x 22 df indicates the reciprocal of a draw from a central chisquare distribution on df degrees of freedom. The samples of g are drawn from a multivariate normal distribution of order n: Its mean,ĝ; is GBLUP computed at the current state of the variance ratio, which varies at random from iteration to iteration; the covariance matrix is  for each of seven CV settings, each replicated 100 times at random. Training sets had size 599 -d; d = 10, 50, 100, 200, 300, 400, and 500.
The current GBLUP is calculated asĝ ¼ C 21 y ¼ ðI þ G 21 lÞ 21 y; in this representation, G 21 must exist. If it does not, a representation of GBLUP that holds iŝ where B is an n · n matrix of regression coefficients of g on y. Likewise Once the Gibbs sampler has been run and burn-in iterations discarded, S samples become available for posterior processing, with sample s consisting of fg with the product of the normal densities taken over members of the testing set. Sampling weights may be very unstable if d is large.

Data availability
The authors state that all data necessary for confirming the conclusions presented in the article are represented fully within the article.

EVALUATION OF METHODOLOGY
The methods described here were evaluated using the wheat data available in package BGLR (Pérez and de los Campos 2014). This data set is well curated and has also been used by, e.g., Crossa et al. (2010), Gianola et al. (2011) and Long et al. (2011). The data originated in trials conducted by the International Maize and Wheat Improvement Center (CIMMYT), Mexico. There are 599 wheat inbred lines, each genotyped with 1279 DArT (Diversity Array Technology) markers and planted in four environments; the target trait was grain yield in environment 1. Sample size was then n ¼ 599 with p ¼ 1279 being the number of markers. These DArT markers are binary ð0; 1Þ and denote presence or absence of an allele at a marker locus in a given line. There is no information on chromosomal location of markers. The objective of the analysis was to illustrate concepts, as opposed to investigate a specific genetic hypothesis. The data set of moderate size allowed extensive replication and reanalysis under various settings.
LOO vs. leave-d-out CV: ordinary least-squares The linear model had an intercept and regressions on markers 301 through 500 in the data file; markers were chosen arbitrarily. Here, p ¼ 201 and n ¼ 599; ensuring unique OLS estimates of substitution effects, i.e., there was no rank deficiency in X. Seven CV layouts were constructed in which testing sets of sizes 2, 3,. . ., 7, or 8 lines were randomly drawn (without replacement) from the 599 inbred lines. Training set sizes decreased accordingly, e.g., for d ¼ 7; training sample size was 599 -7= 592. Larger sizes of testing sets were not considered because ðI 2 H d Þ in (9) became singular as d increased beyond that point. The training-testing process was repeated 300 times at random, to obtain an empirical distribution of prediction mean squared errors.
For LOO CV, regression coefficients were calculated using (5), and the predictive mean squared error was computed as in (7). For the leave-d-out CV, regressions and PMSE were computed with (9) and (11), respectively. Figure 1 shows that the median PMSE for leave-d-out CV was always smaller than the LOO PMSE (horizontal line), although it tended toward the latter as d increased, possibly due to the increasingly smaller training sample size. PMSE in LOO was 1.12, while it ranged from 0.80 to 1.04 for testing sets containing two or more lines. An increase in testing set size at the expense of some decrease in training sample size produced slightly more accurate but less variable predictions (less spread in the distribution of PMSE); this trend can be seen in the box plots depicted in Figure 1. Differences were small but LOO was always less accurate.

BLUP of marker effects
The developments for ridge regression or BLUP of marker effects depend on assuming that allocation of observations into testing sets, with a concomitant decrease in training set size, does not affect the ratio of variance components appreciably.
First, we examined consequences of removing each of the 599 lines at a time on maximum likelihood estimates (MLE) of marker ðs 2 b Þ and residual ðs 2 e Þ variances. The model was as in (1), without an intercept (phenotypes were centered), assuming b $ Nð0; Is 2 b Þ and e $ Nð0; Is 2 e Þ were independently distributed, and with all 1279 markers used when forming XX9: An eigen-decomposition of XX9 coupled with the R function optim (G. de los Campos, personal communication) was used for computing MLE. It was assumed that convergence was always to a global maximum, as it was not practical to monitor the 599 implementations for convergence in each case. MLE of l ¼ With all lines used in the analysis, MLEðlÞ ¼ 189:9: Figure 2 displays the 599 estimates of l, and the resulting empirical cumulative distribution function when LOO was used. Removal of a single line produced MLE of l ranging from 174.5 to 195.6 (corresponding to estimates of s 2 b spanning the range 5:1 2 5:7 · 10 23 ); the 5 and 95 percentiles of the distribution of the LOO estimates of l were 185.9 and 192.7, respectively. Model complexity (Ruppert et al. 2003;Gianola 2013) was gauged by evaluating the "effective number of parameters" as p eff ¼ tr½XðX9X þ IlÞ 21 X9; the "effective degrees of freedom" are n e ¼ n 2 p eff : For the entire data set with n ¼ 599 and p ¼ 1279; variation of l from 174.5 to 195.6 was equivalent to reducing p eff from 164.2 to 155.3, with n e ranging from 435.8 to 443.7. These metrics confirm that the impact of removing a single line from the training process was fairly homogeneous across lines.
Next, we excluded d ¼ 10, 50, 100, 200, 300, 400, and 500 lines from the analysis, while keeping p ¼ 1279 constant. The preceding was done by sampling with replacement the appropriate number of rows from the entire data matrix ðy; XÞ; and removing these rows from the analysis; the procedure was repeated 100 times at random for each value of d; to obtain an empirical distribution of the MLE. must be exerted on marker effects to learn signal from 1279 markers as sample size decreases. (b) MLE of variance components have a finite sample size bias, which might be upwards for l here; bias cannot be measured, so the preceding is conjectural.
In short, it appears that keeping l constant in a LOO setting is reasonable; however, the estimated variance ratio was sensitive with respect to variation in training set size when 100 or more lines, i.e., 15% or more of the total number of cases, were removed for model training.
To assess the impact on PMSE of number of lines removed from training and allocated to testing, 300 testing sets for each of d ¼ 10; 20; . . . ; 80; 90 lines were formed by randomly sampling (with replacement) from the 599 lines. The regression model was trained on the remaining lines using the entire battery of 1279 markers with l ¼ 190. Figure 5 shows the distribution of PMSE for each of the layouts. A comparison with Figure 1 shows that the PMSEs for BLUP were smaller than for OLS; this was expected because, even though training set sizes were smaller than those used for OLS, BLUP predictions with l ¼ 190 are more stable and the model was more complex, since 1279 markers were fitted jointly. As testing set size increased, median PMSE was 0.68 ðd ¼ 10Þ 0.70 ðd ¼ 20Þ and 0.72-0.73 for the other testing set sizes. For LOO, PMSE was 0.72. As in the case of OLS, the distribution of PMSE over replicates became narrower as d grew. As anticipated, decreases in training set size produced a mild deterioration in accuracy of prediction (in an MSE sense) but generated a markedly less variable CV distribution. Testing sets of about 10% of all lines produced a distribution of PMSE with a similar spread to what was obtained with larger testing sets without sacrificing much in mean accuracy. We corroborated that attaining the largest possible training sample is not optimal if done at the expense of testing set size, because predictions are more variable.
Testing sets of size d ¼ 100 to d ¼ 500 (in increments of 50 lines) were evaluated as in the preceding case, again using 300 replicates for each setting and with l ¼ 190: Comparison of Figure 6 with Figure 5 indicates that a marked deterioration in PMSE ensued, which may be  due to insufficient regularization or overfitting from the decrease in training set size. For example, a testing set of size 500 implies that the model with 1279 markers was trained using only 99 inbred lines. In this case, stronger regularization (shrinkage of regression coefficients toward 0) may be needed than what is effected by l ¼ 190: To examine whether overfitting or insufficient regularization was the source of degradation in PMSE, the effective residual degrees of freedom were calculated for each of 70 combinations of d ¼ ð50; 100; . . . ; 450; 500Þ and l ¼ ð100; 150; . . . ; 350; 400Þ; each combination was replicated 50 times by sampling with replacement from ðy; XÞ and extracting the appropriate number of rows: The n e values were averaged over the 50 replications. Figure 7 displays n e plotted against training set size (n 2 d ¼ 99, 149,. . .,499, 549): the impact of variations in l on n e was amplified in absolute and relative terms as training set size increased. For instance, for n 2 d ¼ 99, each observation in the training set contributed 0.48 and 0.74 residual degrees of freedom when l varied from 100 to 400; when n 2 d ¼ 549; the corresponding contributions were 0.64 and 0.82. Figure 8 shows how n e varied with l for each of the training set sizes.
Overfitting did not seem to be the cause of degradation in PMSE because the models "preserved" a reasonable number of degrees of freedom in each case considered. These results reinforce the point that, in shrinkage-based methods, such as GBLUP or any member of the "Bayesian alphabet" Gianola 2013), there is an interplay between sample size, model complexity, and strength of regularization. The effective number of parameters in the training process is given by n 2 d 2 n e ; and here it varied from 25.64 ðn 2 d ¼ 99; l ¼ 400Þ to 198.73 ðn 2 d ¼ 549; l ¼ 100Þ: Even though p ¼1279 markers were fitted, the model was not able to estimate beyond about 200 linear combinations of marker effects. This illustrates that the "prior" (i.e., the distribution assigned to marker effects) matters greatly when n ( p . In other words, there were $ 1079 estimates of marker effects that are statistical artifacts from regularization, and which should not be construed as sensible estimates of marker locus effects, as pointed out by Gianola (2013). Bayesian learning would gradually improve over time if n would grow faster than p; which seems unlikely given a tendency toward overmodeling as sequence and postgenomic data accrue.
Genomic BLUP Standard GBLUP,ĝðlÞ; of genotypic values of the 599 lines ðĝ i ; i ¼ 1; 2; . . . ; 599Þ was computed with l ¼ 190. Subsequently, g ½2i ðlÞ was obtained for each of the 599 lines, i.e., the GBLUP of all lines after removing line i in the training process. Euclidean distances betweenĝðlÞ andĝ ½2i ðlÞ were calculated as this metric measures the extent to which removal of line i influences model training. The minimum and maximum absolute distances were 3 · 10 24 and 1.082, respectively, and the coefficient of variation of distances was about 80%. An observation was deemed influential when d i ðlÞ≧0.83, the 99-percentile of the empirical distribution. Figure 9 (top panel) shows a scatter plot of the d i ðlÞ; influential lines (28, 440, 461, 503, 559, and 580) correspond to points on top of the horizontal line. The relationship between the phenotype of the line excluded in the LOO CV is shown in the bottom panel: larger phenotypes tended to be associated with larger distances.  d = 99, 149, 199, 249, 299, 349, 399, 449, 499, and 549) at selected values of the regularization parameter (l = 100,150, 200, 250, 300, 350, and 400). Values are averages of 50 random replications.
Using data from all lines, GBLUP isĝ ¼ C 21 y; the influence of phenotype of line j ¼ 1; 2; . . . ; 599 on GBLUP of line i is element ij of the matrix @ĝ @y9 ¼ C 21 : Observe that is as a matrix of n · n regression coefficients; its i th row contains the regression of the genotype of line i on the n phenotypes. A measure of overall influence (leverage) of line j is the average of values (or absolute values) of elements in column j of C 21 : Clearly, leverages depend on relatedness structure and on l but not on phenotypes. Figure 10 depicts plots of LOESS regressions (Cleveland 1979) of Euclidean distance between GBLUP calculated with all lines, and LOO GBLUP on two measures of leverage: the average of absolute values of C 21 over rows for each line (leverage 1), and the average of elements of C 21 over rows, by line (leverage 2). LOESS fits (span parameter equal to 0.50) indicated that leverage 1 informs about the impact of removing a specific line in LOO: the larger the leverage 1 of a line, the larger the effect of its removal from the training process.

RKHS
We built a kernel matrix K with typical element (ranging between 0 and 1) for i; j ¼ 1; 2:::; 599: Here x i is the 1279 · 1 vector of marker genotypes in line i; h 1 ; h 2 , and h 3 are bandwidth parameters tuned to establish "global," "regional" and "local" similarities between individuals (as h increases, similarity decreases); w 1 ; w 2 , and w 3 are weights assigned to the three sources of similarity, such that 0 , w i , 1 and P 3 i¼1 w i ¼ 1: We arbitrarily chose h 1 ¼ 1 2 ; h 2 ¼ 2, and h 3 ¼ 4; and w 1 ¼ 0:5; w 2 ¼ 0:30, and w 3 ¼ 0:20. From K we created three additional kernels by placing w i ¼ 1 for i ¼ 1; 2; 3; leading to matrices K 1 ; K 2 , and K 3 : Mean off-diagonal elements of the four kernel matrices were 0.73 ðK 1 Þ, 0.29 ðK 2 Þ, 0.09 ðK 3 Þ, and 0.47 ðKÞ; these values can be interpreted as correlations between pairs of individuals. Hence, K 1 and K 3 produced the highest and lowest degrees of correlation, respectively; complexity of the models increases from kernel 1 to kernel 3 because the fit to the data increases with h (de los .
The four no-intercept RKHS models had the basic form where K was either as in (48) or K i ¼ 1; 2; 3: Variance components were estimated by maximum likelihood, producing as estimates of l ¼ s 2 e s 2 a : l 1 ¼ 0:16; l 2 ¼ 0:30, l 3 ¼ 0:32, and l K ¼ 0:21: The effective number of parameters was calculated (e.g., kernel 1) as p 1 ¼ trðI þ K 21 l 1 Þ 21 ; yielding 224.5, 319.2, and 376.3 for kernel matrices 1, 2, and 3, respectively; for K the effective number of parameters fitted was 330.3. As expected, model complexity increased as the model became more "local." We fitted the four RKHS models to all 599 lines, and conducted a LOO CV for each model. In the fitting process, the corresponding regularization parameter was employed; e.g., for K 2 ; l 2 ¼ 0:30. For each of the models, RKHS predictions of genotypic values were calculated aŝ g i ðl i ; h; wÞ ¼ À I þ K 21 i l i Á 21 y; i ¼ 1; 2; 3; K: The implicit dependence of predictions on bandwidths ðhÞ and weights ðwÞ is indicated in the notation above, but not used hereinafter. In LOO (line j left out in the training process), predictions and where c jj i is the j th diagonal element of ðI þ K 21 i l i Þ 21 ; and Predictive MSEs were 0.6795 ðK 1 Þ; 0.6446 ðK 2 Þ; 0.6555 ðK 3 Þ, and 0.6439 ðKÞ; predictive correlations were 0.566, 0.597, 0.591, and 0.598. Differences between kernels with respect to the criteria used were nil, but the model combining three kernels conveying differing degrees of locality had a marginally smaller MSE and a slightly larger correlation. LOO prediction errors plotted against line phenotypes are shown in Figure 11 for the four kernels used. Prediction errors were larger in absolute value for lines with lowest and highest grain yields, suggesting that the model may benefit by accounting for possibly heterogeneous residual variances. K 2 and K 3 captured some substructure in the distribution of fitted residuals. The more global kernel ðK 1 Þ; arguably capturing mostly additive effects, did not suggest any substructure, which reemerged when the three kernels were combined into K: The preceding exercise illustrates that predictive correlations and PMSE do not fully describe the performance of a prediction machine.
Bayesian GBLUP with known variance components Bayesian GBLUP with known variances has a closed form solution: using all data, the posterior distribution is gjy; s 2 e ; s 2 g $ Nðĝ; C 21 s 2 e Þ; whereĝ ¼ C 21 y and C 21 ¼ ðI þ G 21 lÞ 21 . Set G ¼ XX9; l ¼ 190, and s 2 e ¼ 0:54; X is centered. This problem was attacked (with Monte Carlo error) by drawing independent samples from the 599-variate normal posterior distribution; no MCMC is needed. Using (34), the importance sampling weights for LOO are w i;s ¼ where g ðsÞ i is sample s for line i; g ðsÞ is drawn from Nðĝ; C 21 0:54Þ. The importance sampling weight becomes Observe that the likelihood that y i confers to g ðsÞ i is inversely proportional to w i;s : Hence, samples in which the phenotype removed ðy i Þ confers little likelihood to the g i receive more weight.
We took S ¼15,000 independent samples from Nðĝ; C 21 0:54Þ: The effective number of weights per line, calculated with (38) ranged from 76.4 to 14,983.5; the median (mean) weight was 10,991.0 (9789.2), and the first and third quartiles of the distribution were 6826.1 and 14,983.5, respectively. On average, 1.54 independent samples were required for drawing an effectively independent LOO posterior sample. Figure 12  illustrates the variability among some arbitrarily selected lines of the mean number of importance weights. A phenotype having a small mean importance weight would be "surprising" with respect to the model. However, there are theoretical and numerical issues with the weights used here, a point retaken in the discussion. Figure 13 shows that GBLUP using the entire sample of size n fitted closely. Correlations between predictors were 0.91 in all cases. Meansquared errors were 0.40 (GBLUP, entire sample), 0.73 (Bayes LOO), and 0.73 (LOO GBLUP). The correlation between predictions and phenotypes was 0.81 for GBLUP (entire sample), 0.52 for LOO GBLUP, and 0.51 for Bayes LOO.
The importance sampling scheme worked well. Ionides (2008) suggested a modification of the importance ratio assigned to sample s and data point i; r i;s ¼ p 21 ðy i b ðsÞ ; s 2ðsÞ e Þ, as follows where r i is the average (over samples) importance sample ratio for line i in the context of the wheat data set. Ionides (2008) argued that truncation of large importance sampling weights (TIS) would be less sensitive with respect to the importance sampling density function used (the posterior distribution using all data points in our case) than IS. We converted the IS weights into normalized TIS weights using rule (55) applied to the normalized IS weights. As mentioned earlier, the effective number of normalized IS weights ranged between 76.4 and 14983.5; for normalized TIS weights, the effective number spanned from 691.9 to 13,970. TIS produced more stable weights than standard IS. However, truncation of weights introduces a bias, which may affect predictive performance adversely (Vehtari et al. 2016). However, it may be that TIS made the weights "too homogeneous," thus creating a bias toward the posterior distribution obtained with all data points. If all weights were equal to a constant, IS or TIS would retrieve the full posterior distribution.

DISCUSSION
Cross-validation (CV) has become an important tool for calibrating prediction machines in genome-enabled prediction (Meuwissen et al. 2001), and is often preferred over resampling methods such as bootstrapping. It gives a means for comparing and calibrating methods and training sets (e.g., Isidro et al. 2015). Typically, CV requires dividing data into training and testing folds, and the models must be run many times. An extreme form of CV is LOO; here if the sample has size n; each observation is removed in the training process, and labeled as a testing set of size 1. Hence, n different models must be run to complete a LOO CV. Our paper presented statistical methodology aimed at enabling extensive CV in genome-enabled prediction using a suite of methods. These were OLS, BLUP on markers, GBLUP, RKHS, and Bayesian procedures. Formulae were derived that enable arriving at the predictions that would be obtained if one or more cases were to be excluded from the training process and declared as members of the testing set. In the cases of OLS, BLUP, GBLUP, and RKHS, and assuming that the ratio of variance components do not change appreciably from those that apply to the entire sample, the formulae are exact. The deterministic formulae can also be applied in a multiple-kernel or multiple-random factors setting, given the variance components. For example, consider the bikernel RKHS regression (e.g., de los Campos et al. 2010;Tusell et al. 2014) Here, g P ¼ K P a P is genetic signal captured by pedigree, and g M ¼ K M a M is genetic signal captured by markers; a P and a M are unknown and independently distributed RKHS regression coefficients, and K P and K M are positive-definite similarity matrices. Suppose that K P ¼ A, i.e., the numerator relationship matrix based on the assumption of additive inheritance, and that K M is a Gaussian Figure 12 Effective number of importance samples for LOO Bayesian GBLUP (given the variances) for selected lines; 15,000 independent samples were drawn from the posterior distribution of genotypic values.
kernel such as those employed earlier in the paper. The standard assumption for the RKHS regression coefficients is where s 2 aA and s 2 aM are variance components associated with kernels A and K M ; respectively. Hence, g P $ Nð0; As 2 aA Þ and g M $ Nð0; K M s 2 aM Þ: The BLUP of g P and g M can be found by solving the linear system . Hence, the solutions satisfŷ Applying the logic leading to (102) for the LOO situation, the preceding equations can be written as where c ii A ðc ii M Þ are the i th diagonal elements of C 21 A ðC 21 M Þ; respectively; g i;P andĝ i;M are the corresponding solutions to the system of equations (58). The prediction of the left-out phenotype isỹ i ¼g i;P þg i;M : The situation is different for Bayesian models solved by sampling methods such as MCMC. Typically, there is no closed form solution, except in some stylized situations, so sampling must be used. The Bayesian model must be run with the entire data set and posterior samples weighted, e.g., via importance sampling, to convert realizations into draws pertaining to the posterior distribution that would result from using just the CV training set. Unfortunately, importance weights can be extremely variable. In (34), it can be seen that weights are proportional to the reciprocal of likelihoods evaluated at the posterior samples, so if a data point (or a vector of data points) "left out" confers a tiny likelihood to the realized value of the parameter sampled, then the sample is assigned a very large weight; on the other hand, if the likelihood is large, the weight is small. This phenomenon produces a large variance among weights, which we corroborated empirically. Another view at the issue at stake is as follows: the importance weight is w i ðbÞ ¼ pðbjy 2i ; HÞ pðbjy; HÞ ; hence, if the posterior obtained with all data points has much thinner tails than the posterior density constructed by excluding one or more cases, the weights can "blow up." The preceding problem would be exacerbated by including more than one observation in the testing set. Vehtari et al. (2016) examined seven data sets, and used "brute force" LOO MSE as a gold standard to examine the performance of various forms of importance sampling. TIS gave a better performance than standard importance sampling (IS) weights in two out of seven comparisons, with the standard method being better in three of the data sets; there were two ties. Vehtari et al. (2016) suggested another method called Pareto smoothed importance sampling (PSIS) that was better than IS in four of the data sets (two ties), and better than TIS in three of the analyses (two ties). Calculation of PSIS is involved, requiring several steps in our context: (a) compute IS ratios; (b) for each of the 599 lines, fit a generalized Pareto distribution to the 20% largest values found in (a); (c) replace the largest M IS ratios by expected values of the order statistics of the fitted generalized Pareto distribution, where M ¼ 0:20 · S; (d) for each line, truncate the new ratios. Clearly, this procedure does not lend itself to large scale genomic data, and the results of Vehtari et al. (2016), obtained with small data sets and simple models, are not conclusive enough. Additional research is needed to examine whether TIS, IS, or PSIS are better for genome-enabled prediction.
It is known (e.g., Henderson 1973Henderson , 1984Searle 1974) that the best predictor, that is, the function of the data with the smallest squared prediction error (MSE) under conceptual repeated sampling (i.e., infinite number of repetitions over the joint distribution of predictands and predictors) is Eðgjy; parametersÞ. This property requires knowledge of the form of the joint distribution, and of its parameters. In the setting of the case study, under multivariate normality, and with a zero-mean model and known variance components, GBLUP is the best predictor. However, in CV, the property outlined above does not hold. One reason is that the data set represents a single realization of the conceptual scheme. Another reason is that incidence and similarity matrices change at random in CV, plus parameters are estimated from the data at hand. For example, if datum i is removed from the analysis, the training model genomic relationship matrix becomes G ½2i;2i ; whereas, if observation j is removed, the matrix used is G ½2j;2j : Further, the entire data set is used in the CV, so yet-to-be observed data points appear in the training process at some point. The setting of best prediction requires that the structure of the data remains constant over repeated sampling, with the only items changing being the realized values of the data ðyÞ, and of the unobserved genotypic values ðgÞ: The CV setting differs from the idealized scheme, and expectations based on theory may not always provide the best effective guidance in predictive inference.
In conclusion, CV appears to be the best gauge for calibrating prediction machines. Results presented here provide the basis for conducting extensive cross-validation from results of a single run with all data. Future research should evaluate importance sampling schemes for more complex Bayesian models, e.g., those using thick-tailed processes or mixtures as prior distributions. An important challenge is to make the procedures developed here computationally cost-effective, so that software for routine use can be developed.

ACKNOWLEDGMENTS
The Editor and two anonymous reviewers suggested useful comments concerning the structure of the manuscript. We are grateful to the Institut Pasteur de Montevideo, Uruguay, for providing office space and facilities to D.G. from January-March 2016. Part of this work was done while D.G. was a Hans Fischer Fellow at the Institute of Advanced Study, Technical University of Munich-München, Germany; their support is gratefully acknowledged. Research was partially supported by a United States Department of Agriculture (USDA) Hatch Grant (142-PRJ63CV) to D.G., and by the Wisconsin Agriculture Experiment Station.

LITERATURE CITED APPENDIX A: PREDICTION WITH LEAST-SQUARES
Following Seber and Lee (2003), consider out-of-sample prediction, and suppose that the distribution of residuals in left-out data (testing set) of size n test is as in a training set of size n, with the training phenotypes represented as y; residuals in training and testing sets are assumed to be mutually independent. In genome-enabled prediction, (1) is merely an instrumental model that may bear little resemblance with the state of nature, so we let the true distribution be y test $ Nðm; Is 2 e Þ; allowing for the almost certain possibility that m 6 ¼ X test b, where X test is the marker matrix in the testing set. In quantitative genetics, m is an unknown function of quantitative trait locus (QTL) genotypes and of their effects; the latter may not be additive, and dominance or epistasis may enter into the picture without contributing detectable variance.
Model training yieldsb, and the point predictor of y test is X testb . The mean squared error of prediction ðPMSEÞ, conditionally on training data ðy; XÞ and on X test but averaged with respect of all possible testing sets of size, n test , is EðPMSEjy; X; X test Þ ¼ 1 n test E y test jy;X;Xtest h À y test 2X testb Á 9 À y test 2 X testb Á i : Define the conditional prediction bias asd ¼ E y test jy;X;Xtest À y test 2 X testb jy Á ¼ m 2 X testb : Standard results on expectation of quadratic forms applied to (63) Here, the expected prediction bias is after assuming for convenience of representation that Eðy test Þ ¼ EðyÞ ¼ m; i.e., that testing and training sets have the same size and unknown true mean m; note that @X testb @y9 ¼ X test ðX9XÞ 21 X9 is a matrix of "influences" describing how variation in training phenotypes affects predictions.
Using (68) in (67), and this in (66), the expected prediction mean-squared error averaged over an infinite number of training and testing sets, but with fixed marker genotype matrices is If X test ¼ X; that is, in the special case where the same genotypes appear in the training and testing sets, increases with p (model complexity) and decreases with n test (equal to n here). The impact of increasing complexity on bias is difficult to discuss in the absence of knowledge of QTL and their effects. When p/n, the training data set fits better and better and the model increasingly copies both signal and error: predictions become increasingly poorer. Note that Var½ðy 2 XbÞ ¼ MMs 2 e ¼ ðI 2 HÞs 2 e ; where M ¼ I 2 XðX9XÞ 21 X9 ¼ I 2 H. Applying this result in (8) where 1 n ðd9DdÞ is the average squared prediction bias, and s 2 e n tr½DM is the average prediction error variance. Note that and that Employing (74) and (75) in (73), Examine the difference ðD MSE Þ in expected prediction mean-squared error between the LOO procedure, as given in (76), and that from the "distinct testing set" layout given in (70), and set n ¼ n test . One obtains Since p ¼ trðXðX9XÞ 21 X9Þ ¼ trðHÞ ¼ P n i¼1 h ii ; The h ii are bounded between 0 and 1, so the two terms in (78) are positive.