Semi-parametric genomic-enabled prediction of genetic values using reproducing kernel Hilbert spaces methods

Summary

Prediction of genetic values is a central problem in quantitative genetics. Over many decades, such predictions have been successfully accomplished using information on phenotypic records and family structure usually represented with a pedigree. Dense molecular markers are now available in the genome of humans, plants and animals, and this information can be used to enhance the prediction of genetic values. However, the incorporation of dense molecular marker data into models poses many statistical and computational challenges, such as how models can cope with the genetic complexity of multi-factorial traits and with the curse of dimensionality that arises when the number of markers exceeds the number of data points. Reproducing kernel Hilbert spaces regressions can be used to address some of these challenges. The methodology allows regressions on almost any type of prediction sets (covariates, graphs, strings, images, etc.) and has important computational advantages relative to many parametric approaches. Moreover, some parametric models appear as special cases. This article provides an overview of the methodology, a discussion of the problem of kernel choice with a focus on genetic applications, algorithms for kernel selection and an assessment of the proposed methods using a collection of 599 wheat lines evaluated for grain yield in four mega environments.

1. Introduction

Prediction of genetic values is relevant in plant and animal breeding, as well as for assessing the probability of disease in medicine. Standard genetic models view phenotypic outcomes (y_i; i=1, …, n) as the sum of a genetic signal (g_i) and of a residual (∊_i), that is: y_i=g_i+∊_i. The statistical learning problem consists of uncovering genetic signal from noisy data, and predictions (ĝ_i) are constructed using phenotypic records and some type of knowledge about the genetic background of individuals.

Family structure, usually represented as a pedigree, and phenotypic records have been used for the prediction of genetic values in plants and animals over several decades (e.g. Fisher, 1918; Wright, 1921; Henderson, 1975). In pedigree-based models (P), a genealogy is used to derive the expected degree of resemblance between relatives, measured as Cov(g_i, g _i′), and this provides a means for smoothing phenotypic records.

Dense molecular marker panels are now available in humans and in many plant and animal species. Unlike pedigree data, genetic markers allow follow-up of Mendelian segregation; a term that in additive models and in the absence of inbreeding accounts for 50% of the genetic variability. However, incorporating molecular markers into models poses several statistical and computational challenges such as how models can cope with the genetic complexity of multi-factorial traits (e.g. Gianola & de los Campos, 2008), and with the curse of dimensionality that arises when a large number of markers is considered. Parametric and semi-parametric methods address these two issues in different ways.

In parametric regression models for dense molecular markers (e.g. Meuwissen et al., 2001), g_i is a parametric regression on marker covariates, x_ik with k=1, …, p indexing markers. The linear model takes the form: , where β_k is the regression of y_i on x_ik. Often, p>>n and some shrinkage estimation method such as ridge regression (Hoerl & Kennard, 1970a, 1970b) or LASSO (Least Absolute Shrinkage and Selection Operator, Tibshirani, 1996), or their Bayesian counterparts, are used to estimate marker effects. Among the latter, those using marker-specific shrinkage such as the Bayesian LASSO of Park & Casella (2008) or methods BayesA or BayesB of Meuwissen et al. (2001) are the most commonly used. In linear regressions, dominance and epistasis may be accommodated by adding appropriate interactions between marker covariates to the model; however, the number of predictor variables is extremely large and modelling interactions is only feasible to a limited degree.

Reproducing kernel Hilbert spaces (RKHS) regressions have been proposed for semi-parametric regression on marker genotypes, e.g. Gianola et al. (2006) and Gianola & van Kaam (2008) . In RKHS, markers are used to build a covariance structure among genetic values; for example, Cov(g_i, g _i′)∝K(x_i, x_i′), where x_i, x_i′ are vectors of marker genotypes and K(., .), the reproducing kernel (RK), is some positive definite (PD) function (de los Campos et al., 2009a). This semi-parametric approach has several attractive features: (a) the methodology can be used with almost any type of information set (e.g. covariates, strings, images and graphs). This is particularly important because techniques for characterizing genomes change rapidly; (b) some parametric methods for genomic selection (GS) appear as special cases and (c) computations are performed in an n-dimensional space. This provides RKHS methods with a great computational advantage relative to some parametric methods, especially when p>>n.

This article discusses and evaluates the use of RKHS regressions for genomic-enabled prediction of genetic values of complex traits. Section 2 gives a brief review of RKHS regressions. A special focus is placed on the problem of kernel choice. We discuss cases where a genetic model (e.g. additive infinitesimal) is used to choose the kernel and others where the RK is chosen based on its properties (e.g. predictive ability). Section 3 presents an application to an extensive plant breeding data set where some of the methods discussed in Section 2 are evaluated. Concluding remarks are provided in Section 4.

2. RKHS regression

RKHS methods have been used in many areas of application such as spatial statistics (e.g. ‘Kriging’; Cressie, 1993), scatter-plot smoothing (e.g. smoothing splines; Wahba, 1990) and classification problems (e.g. support vector machines; Vapnik, 1998), just to mention a few. Estimates in RKHS regressions can be motivated as solutions to a penalized optimization problem in an RKHS or as posterior modes in a certain class of Bayesian models. A brief description of RKHS estimates in the context of penalized estimation is given first in section 2(i), with its Bayesian interpretation introduced later in section 2(ii). A representation of RKHS regressions that uses orthogonal basis functions is given in section 2(iii). This section ends in 2(iv) with a discussion of the problem of kernel choice.

(i) Penalized estimation in RKHS

A standard problem in statistical learning consists of extracting signal from noisy data. The learning task can be described as follows (Vapnik, 1998): given data {(y_i, t_i)}_{i =1}ⁿ, originating from some functional dependency, infer this dependency. The pattern relating input, t_i∊T, and output, y_i∊Y, variables can be described with an unknown function, g, whose evaluations are g_i=g(t_i). For example, t_i may be a vector of marker genotypes, t_i=x_i and g may be a function assigning a genetic value to each genotype. Inferring g requires defining a collection (or space) of functions from which an element, ĝ, will be chosen via a criterion (e.g. a penalized residual sum of squares or a posterior density) for comparing candidate functions. Specifically, in RKHS, estimates are obtained by solving the following optimization problem:

(1)

where g∊H denotes that the optimization problem is performed within the space of functions H, a RKHS; l(g, y) is a loss function (e.g. some measure of goodness of fit); λ is a parameter controlling trade-offs between goodness of fit and model complexity; and ||g||_H ² is the square of the norm of g on H, a measure of model complexity. A technical discussion of RKHS of real-valued functions can be found in Wahba (1990) ; here, we introduce some elements that are needed to understand how ĝ is obtained.

Hilbert spaces are complete linear spaces endowed with a norm that is the square root of the inner product in the space. The Hilbert spaces that are relevant for our discussion are RKHS of real-valued functions, here denoted as H. An important result, known as the Moore–Aronszajn theorem (Aronszajn, 1950), states that each RKHS is uniquely associated with a PD function that is a function, K(t_i, t _i′), satisfying for all sequences, {α_i}, with α_i≠0 for some i. This function, K(t_i, t _i′), also known as the RK, provides basis functions and an inner product (therefore a norm) to H. Therefore, choosing K(t_i, t _i′) amounts to selecting H; the space of functions where (1) is solved.

Using that duality, Kimeldorf & Wahba (1971) showed that the finite-dimensional solution of (1) admits a linear representation , or in matrix notation, , where K={K(t_i, t _i′)} is an n×n matrix whose entries are the evaluations of the RK at pairs of pint in T. Further, in this finite-dimensional setting, . Using this in (1) and setting l(g, y) to be a residual sum of squares, one obtains: , where is the solution of

(2)

and y={y_i} is a data vector. The first-order conditions of (2) lead to . Further, since K=K′ and K⁻¹ exists, pre-multiplication by K⁻¹ yields, . Therefore, the estimated conditional expectation function is , where P(λ, K)=K(K+λI)⁻¹ is a smoother or influence matrix.

The input information, t_i∊T, enters into the objective function and on the solution only through K. This allows using RKHS for regression with any class of information sets (vectors, graphs, images, etc.) where a PD function can be evaluated; the choice of kernel becomes the key element of model specification.

(ii) Bayesian interpretation

From a Bayesian perspective, can be viewed as a posterior mode in the following model: ; P(ε, α|σ_∊², σ_g ²)=N(ε|0, Iσ_∊²)N(α|0, K⁻¹σ_g ²). The relationship between RKHS regressions and Gaussian processes was first noted by Kimeldorf & Wahba (1970) and has been revisited by many authors (e.g. Harville, 1983; Speed, 1991). Following de los Campos et al. (2009a), one can change variables in the above model, with g=K α, yielding

(3)

Thus, from a Bayesian perspective, the evaluations of functions can be viewed as Gaussian processes satisfying Cov(g_i, g_{i ′})∝K(t_i, t _i′). The fully Bayesian RKHS regression assumes unknown variance parameters, and the model becomes

(4)

where p(σ_∊², σ_g ²) is a (proper) prior density assigned to variance parameters.

(iii) Representation using orthogonal random variables

Representing model (4) with orthogonal random variables simplifies computations greatly and provides additional insights into the nature of the RKHS regressions. To this end, we make use of the eigenvalue (EV) decomposition (e.g. Golub & Van Loan, 1996) of the kernel matrix K=ΛΨΛ′ , where Λ is a matrix of eigenvectors satisfying Λ′Λ=I; Ψ=Diag{Ψ_j}, Ψ₁⩾Ψ₂⩾ …⩾Ψ_n>0, is a diagonal matrix whose non-zero entries are the EVs of K; and j=1, …, n, indexes eigenvectors (i.e. columns of Λ) and the associated EV. Using these, (4) becomes

(5)

To see the equivalence of (4) and (5), note that Λδ is multivariate normal because so is δ. Further, E(Λδ)=ΛE(δ)=0 and Cov(Λδ)=ΛΨΛ′σ_g ²=Kσ_g ². Therefore, equations (4) and (5) are two parameterizations of the same probability model. However, equation (5) is much more computationally convenient, as discussed next.

The joint posterior distribution of (5) does not have a closed form; however, draws can be obtained using a Gibbs sampler. Sampling regression coefficients from the corresponding fully conditional distribution, p(δ|y, σ_∊², σ_g ²), is usually the most computationally demanding step. From standard results of Bayesian linear models, one can show that , where ELSE denotes everything else other than δ, C=[Λ′Λ+σ_∊²σ_g ⁻² Ψ⁻¹]=Diag{1+σ_∊²σ_g ⁻²Ψ_j ⁻¹} and . This simplification occurs because Λ′Λ=I and Ψ=Diag{Ψ_j}. The fully conditional distribution of δ is multivariate normal, and the (co)variance matrix, σ_∊² C⁻¹, is diagonal; therefore . Moreover, p(δ_j|ELSE) is normal, centred at [1+σ_∊²σ_g ⁻²Ψ _j ⁻¹]⁻¹y _.j and with variance σ_∊²[1+σ_∊²σ_g ⁻²Ψ _j ⁻¹]⁻¹. Here, y _.j=λ′_j y, where λ_j is the jth eigenvector (i.e. the jth column of Λ). Note that model unknowns are not required for computing y _.j, implying that these quantities remain constant across iterations of a sampler. The only quantities that need to be updated are [1+σ_∊²σ_g ⁻²Ψ_j ⁻¹] and σ_∊²[1+σ_∊²σ_g ⁻²Ψ_j ⁻¹]. If model (5) is extended to include other effects (e.g. an intercept or some fixed effects), the right-hand side of the mixed model equations associated to p(δ|ELSE) will need to be updated at each iteration of the sampler; however, the matrix of coefficients remains diagonal and this simplifies computations greatly (see Appendix).

In equation (5), the conditional expectation function is a linear combination of eigenvectors: . The EV are usually sorted such that Ψ₁⩾Ψ₂⩾ …⩾Ψ_n>0. The prior precision variance of regression coefficients is proportional to the EV, that is, Var(δ_j)∝Ψ_j. Therefore, the extent of shrinkage increases as j does. For most RKs, the decay of the EV will be such that for the first EV [1+σ_∊²σ_g ⁻²Ψ_j ⁻¹] is close to one, yielding negligible shrinkage of the corresponding regression coefficients. Therefore, linear combinations of the first eigenvectors can then be seen as components of g that are (essentially) not penalized.

(iv) Choosing the RK

The RK is a central element of model specification in RKHS. Kernels can be chosen so as to represent a parametric model, or based on their ability of predicting future observations. Examples of these two approaches are discussed next.

The standard additive infinitesimal model of quantitative genetics (e.g. Fisher, 1918; Henderson, 1975), is an example of a model-driven kernel (e.g. de los Campos et al., 2009a). Here, the information set (a pedigree) consists of a directed acyclic graph and K(t_i, t _i′) gives the expected degree of resemblance between relatives under an additive infinitesimal model. Another example of an RKHS regression with a model-derived kernel is the case where K is chosen to be a marker-based estimate of a kinship matrix (usually denoted as G, cf, Ritland, 1996; Lynch & Ritland, 1999; Eding & Meuwissen, 2001; Van Raden, 2007; Hayes & Goddard, 2008). An example of a (co)variance structure derived from a quantitative trait locus (QTL)-model is given in Fernando & Grossman (1989) .

Ridge regression and its Bayesian counterpart (Bayesian ridge regression (BRR)) can also be represented using (4) or (5). A BRR is defined by y=Xβ+ε and p(ε, β, σ_∊², σ_β²)=N(ε|0, Iσ_∊²)N(β|0, Iσ_β²)p(σ_∊², σ_β²). To see how a BRR constitutes a special case of (5), one can make use of the singular value decomposition (e.g. Golub & Van Loan, 1996) of X=UDV′. Here, U (n×n) and V (p×n) are matrices whose columns are orthogonal, and D=Diag{ξ_j} is a diagonal matrix whose non-null entries are the singular values of X. Using this in the data equation, we obtain y=UDV′β+ε=Uδ+ε, where δ=DV′β. The distribution of δ is multivariate normal because so is that of β. Further, E(δ)=DV′E(β)=0 and Cov(δ)=DV′VD′σ_β²=DD′σ_β²; thus, δ~N[0, Diag{ξ_j ²}σ_β²]. Therefore, a BRR can be equivalently represented using (5) with Λ=U and Ψ=Diag{ξ_j ²}. Note that using Λ=U and Ψ=Diag{ξ_j ²} in (5) implies K=UDD′U′=UDV′VD′U′=XX′ in (4). Habier Fernando & Dekkers (2009) argue that as the number of markers increases, XX′ approaches the numerator relationship matrix, A. From this perspective, XX′ can also be viewed just as another choice for an estimate of a kinship matrix. However, the derivation of the argument follows the standard treatment of quantitative genetic models where genotypes are random and marker effects are fixed, whereas in BRR, the opposite is true (see Gianola et al., 2009 for further discussion).

In the examples given above, the RK was defined in such a manner that it represents a parametric model. An appeal of using parametric models is that estimates can be interpreted in terms of the theory used for deriving K. For example, if K=A then σ_g ² is interpretable as an additive genetic variance and σ_g ²(σ_g ²+σ_∊²)⁻¹ can be interpreted as the heritability of the trait. However, these models may not be optimal from a predictive perspective. Another approach (e.g. Shawe-Taylor & Cristianini, 2004) views RKs as smoothers, with the choice of kernel based on their predictive ability or some other criterion. Moreover, the choice of the kernel may become a task of the algorithm.

For example, one can index a Gaussian kernel with a bandwidth parameter, θ, so that K(t_i,t _i′|θ)=exp{−θd(t_i,t _i′)}. Here, d(t_i,t _i′) is some distance function and θ controls how fast the covariance function drops as points get further apart as measured by d(t_i, t _i′). The bandwidth parameter may be chosen by cross-validation (CV) or with Bayesian methods (e.g. Mallick et al., 2005). However, when θ is treated as uncertain in a Bayesian model with Markov chain Monte Carlo (MCMC) methods, the computational burden increases markedly because the RK must be computed every time that a new sample of θ becomes available. It is computationally easier to evaluate model performance over a grid of values of θ; this is illustrated in section 3.

The (co)variance structure implied by a Gaussian kernel is not derived from any mechanistic consideration; therefore, no specific interpretation can be attached to the bandwidth parameter. However, using results for infinitesimal models under epistasis one could argue that a high degree of epistatic interaction between additive infinitesimal effects may induce a highly local (co)variance pattern in the same way that large values of θ do. This argument is revisited later in this section.

The decay of the EV controls, to a certain extent, the shrinkage of estimates of δ and, with this, the trade-offs between goodness of fit and model complexity. Transformations of EV (indexed with unknown parameters) can also be used to generate a family of kernels. One such example is the diffusion kernel K_α=ΛDiag{exp(αΨ_j)}Λ′ (e.g. Kondor & Lafferty, 2002). Here, α>0 is used to control the decay of EV. In this case, the bandwidth parameter can be interpreted as a quantity characterizing the diffusion of signal (e.g. heat) along edges of a graph, with smaller values being associated with more diffusion.

A third way of generating families of kernels is to use closure properties of PD functions (Shawe-Taylor & Cristianini, 2004). For example, linear combinations of PD functions, , with σ_{g .}²⩾0, are PD as well. From a Bayesian perspective, and are interpretable as variance parameters. To see this, consider extending (4) to two random effects so that: g=g₁+g₂ and, p(g₁, g₂|, )=N(g₁|0, K₁ )N(g₂|0, K₂ ). It follows that g~N(0, K₁ +K₂ ), or equivalently , where and . Therefore, fitting an RKHS with two random effects is equivalent to using in (4). Extending this argument to r kernels one obtains: , where . For example, one can obtain a sequence of kernels, {K_r}, by evaluating a Gaussian kernel over a grid of values of a bandwidth parameter {θ_r}. The variance parameters, , associated with each kernel in the sequence can be viewed as weights. Inferring these variances amounts to inferring a kernel, , which can be seen as an approximation to an optimal kernel. We refer to this approach as kernel selection via kernel averaging (KA); an example of this is given in section 3.

The Haddamard (or Schur) product of PD functions is also PD, that is, if K ₁(t_i, t _i′) and K ₂(t_i, t _i′) are PD, so is K(t_i, t _i′)=K ₁(t_i, t _i′)K ₂(t_i, t _i′); in matrix notation, this is usually denoted as K=K₁#K₂. From a genetic perspective, this formulation can be used to accommodate non-additive infinitesimal effects (e.g. Cockerham, 1954; Kempthorne, 1954). For example, under random mating, linkage equilibrium and in the absence of selection, K=A#A={a(i, i′)²} gives the expected degree of resemblance between relatives under an infinitesimal model for additive×additive interactions. For epistatic interaction between infinitesimal additive effects of qth order, the expected (co)variance structure is, K={a(i, i′)^{q +1}}. Therefore, for q⩾1 and i≠i′, the prior correlation,

decreases, i.e. the kernel becomes increasingly local, as the degree of epistatic interaction increases, producing an effect similar to that of a bandwidth parameter of a Gaussian kernel.

3. Application to plant breeding data

Some of the methods discussed in the previous section were evaluated using a data set consisting of a collection of historical wheat lines from the Global Wheat Breeding Programme of CIMMYT (International Maize and Wheat Improvement Center). In plant breeding programmes, lines are selected based on their expected performance and collecting phenotypic records is expensive. An important question is whether phenotypes collected on ancestor lines, together with pedigrees and markers, can be used to predict performance of lines for which phenotypic records are not available yet. If so, breeding programmes could perform several rounds of selection based on marker data only; with phenotypes measured every few generations. The reduction in generation interval attainable by selection based on markers may increase the rate of genetic progress and, at the same time, the cost of phenotyping would be reduced (e.g. Bernardo & Yu, 2007; Heffner et al., 2009). Thus, assessing the ability of a model to predict future outcomes is central in breeding programmes.

The study presented in this section attempted to evaluate: (a) how much could be gained in predictive ability by incorporating marker information into a pedigree-based model, (b) how sensitive these results are with respect to the choice of kernel, (c) whether or not Bayesian KA is effective for selecting kernels and (d) how RKHS performs relative to a parametric regression model, the Bayesian LASSO (BL; Park & Casella, 2008).

(i) Materials and methods

The data comprise family, marker and phenotypic information of 599 wheat lines that were evaluated for grain yield (GY) in four environments. Single-trait models were fitted to data from each environment. Marker information consisted of genotypes for 1447 Diversity Array Technology (DArT) markers, generated by Triticarte Pty. Ltd (Canberra, Australia; http://www.triticarte.com.au). Pedigree information was used to compute additive relationships between lines (i.e. twice the kinship coefficient; Wright, 1921) using the Browse application of the International Crop Information System, as described in McLaren et al. (2005) .

A sequence of models was fitted to the entire data set and in a CV setting. Figure 1 gives a summary of the models considered. In all environments, phenotypes were represented using equation y_i=μ+g_i+∊_i, where y_i (i=1, …, 599) is the phenotype of the ith line; μ is an effect common to all lines; g_i is the genetic value of the ith line; and ∊_i is a line-specific residual. Phenotypes were standardized to a unit variance in each of the environments. Residuals were assumed to follow a normal distribution , where σ_∊² is the residual variance. The conditional distribution of the data was , where g=(g ₁, …, g_n)′. Models differed on how g_i was modelled.

Fig. 1. Alternative models for prediction of genetic values. Phenotypic records (y) were always the sum of a genetic signal (g) and a vector of Gaussian residuals (ε). Models differed on how g was represented, as described in the figure. BL, Bayesian LASSO; RKHS, reproducing kernel Hilbert spaces regression; λ, LASSO regularization parameter; θ, RKHS bandwidth parameter; σ_·², variance parameter; KA, kernel averaging; N(., .), normal density; DE(.), double-exponential density.

In a standard infinitesimal additive model (P, standing for pedigree-model), genetic values are g=a with p(a|σ_a ²)=N(0, Aσ_a ²), where σ_a ² is the additive genetic variance and A={a(i,i′)}, as before, is the numerator relationship matrix among lines computed from the pedigree. This is a RKHS with K=A.

For marker-based models (M), two alternatives were considered: BL and RKHS regression. In the BL, genetic values were a linear function of marker covariates, g=Xβ, where X is an incidence matrix with marker genotypes codes and β=(β₁, …, β_p)′, the vector of regression coefficients, was inferred using the BL of Park & Casella (2008) . Following de los Campos et al. (2009b), the prior density of the regularization parameter of the BL, here denoted as , was , which is flat over a fairly wide range. This model is denoted as M_BL.

In marker-based RKHS regressions (M_K) g=f_θ, where f_θ=(f _θ,1, …, f _θ,n)′ was assigned a Gaussian prior with null mean and (co)variance matrix Cov(f_θ)∝K_θ={exp(−θk ⁻¹d_ii′)}. Here, θ is a bandwidth parameter, d_{ii ′}=||x_i−x_i′||² is the square Euclidean distance between marker codes x_i=(x_{i 1}, …, x_ip)′ and x_{i ′}=(x_{i′ 1}, …, x_i′p)′, and . Models were fitted over a grid of values of θ and are denoted as M_{k ,θ}. The optimal value of the bandwidth parameter is expected to change with many factors such as: (a) distance function; (b) number of markers, allelic frequency and coding of markers, all factors affecting the distribution of observed distances and (c) genetic architecture of the trait, a factor affecting the expected prior correlation of genetic values (see section 2(iv)). We generated a grid of values, θ∊{0·1, 0·25, 0·5, 0·75, 1, 2, 3, 5, 7, 10}, that for this data set allowed exploring a wide variety of kernels. Figure 2 gives a histogram of the evaluations of the kernel for two extreme values of the bandwidth parameter; θ=0·25 gives very high prior correlations, while θ=7 gives a kernel matrix with very low correlations in the off-diagonal.

Fig. 2. Histogram of the evaluations of Gaussian kernel K(i,i′)=exp{−θk ⁻¹d_{ii ′}} by value of the bandwidth parameter (θ=0·25 left and θ=7, right). Here, d_{ii ′}=||x_i−x_i′||² is the squared Euclidean distance between marker codes x_i=(x_{i 1}, …, x_ip)′ and x_i′=(x_{i′ 1}, …, x_i′p)′ , and .

A model where g was the sum of two components: g=f_0·25+f₇, with p(f_0·25, f₇|, )=N(f_0·25|0,K_0·25 )N(f₇|0, K₇ ) was fitted as well. This model is referred to as M_KA, standing for marker-based model with ‘kernel-averaging’. Note that K_0·25 and K₇ provide very different kernels (see Fig. 2). With more extreme values of the bandwidth parameter, marker information is virtually lost. Indeed, choosing θ=0 gives a kernel matrix full of ones and θ→∞ gives K_θ→I, and averaging these two kernels gives a resulting (co)variance structure that does not use marker information at all.

Finally, a sequence of models including pedigree and marker data (PM) was obtained by setting g=a+Xβ, denoted as PM_BL; g=a+f_θ, θ={0·1, 0·25, 0·5, 0·75, 1, 2, 3, 5, 7, 10}, denoted as PM_{k ,θ}; and, g=a+f_0·25+f₇, denoted as PM_KA.

In all models, variance parameters were treated as unknown and assigned identical independent scaled inverse chi-square prior distributions with three degrees of freedom and scale parameters equal to 1, p(σ_·²)=χ⁻²(σ_·²|df=3, S=1). Samples from posterior distributions for each of the models were obtained with a Gibbs sampler (see de los Campos et al., 2009b, for the case of M_BL and PM_BL, and the Appendix for RKHS models). Inferences were based on all 35 000 samples obtained after discarding 2000 samples as burn-in. The distribution of prediction errors was estimated using a 10-fold CV (e.g. Hastie et al., 2009).

(ii) Results

Figure 3 shows the posterior means of the residual variance in M_{k ,θ} and PM_{k ,θ} versus values of the bandwidth parameter θ obtained when models were fitted to the entire data. Each panel in Fig. 3 corresponds to one environment and the horizontal lines give the posterior means of the residual variance from P and PM_KA. Table A1 of the Appendix gives estimates of the posterior means and of the posterior standard deviations of the residual variance from each of the 25 models, by environment. The posterior means of the residual variances indicate that models M and PM fitted the data better than P, and PM_KA gave almost always better fit than M_{k ,θ} and PM_{k ,θ}. In all environments, the posterior mean of the residual variance decreased monotonically with θ; this was expected because K_θ becomes increasingly local as the bandwidth parameter increases. In environments 2, 3 and 4, the slopes of the curves relating the posterior mean of residual variance to θ were gentler for PM_{k ,θ} than for M_{k ,θ}. This occurs, because in PM_{k ,θ}, the regression function has two components, one of which, the regression on the pedigree, is not a function of the bandwidth parameter. Models M_BL and PM_BL did not fit the training data as well as most of the RKHS counterparts, with a posterior mean of the residual variance that was close to that of M_{k ,0·1} and PM_{k ,0·5}, respectively (see Table A1 of the Appendix).

Fig. 3. Estimated posterior mean of the residual variance versus values of the bandwidth parameter, θ, by environment and model. K _θ is a marker-based RKHS regression with bandwidth parameter θ; Pedigree & Markers K _θ uses pedigree and markers, here, θ is the value of the bandwidth parameter for markers. Pedigree & Markers K _0·25+K ₇ uses pedigree and markers with kernel averaging (KA) E1–E4: environments where the lines were evaluated.

The contribution of a, that is, the regression on the pedigree, to the conditional expectation function, g, can be assessed via the posterior mean of σ_a ² (see Fig. A1 in the Appendix). The posterior mean of σ_a ² was larger in P models than in their PM counterparts; this was expected, because in P the regression on the pedigree is the only component of the conditional expectation function that contributes to phenotypic variance. Within PM_{k ,θ}, the posterior mean of σ_a ² was minimum at intermediate values of the bandwidth parameters. At extreme values of θ, the RK may not represent the types of patterns present in the data and, thus, the estimated conditional expectation function would depend more strongly on the regression on the pedigree (large values of σ_a ²).

Plots in Fig. 4 give the estimated mean-squared error (MSE) between CV predictions and observations versus values of the bandwidth parameter (x-axis), by environment and model. The predictive MSE of the P and PM_KA models are displayed as horizontal dashed lines, and values of those for the BL (both in M_BL and PM_BL) are shown at the bottom of the panels. Table A2 in the Appendix gives the estimated MSE by model and environment, respectively.

Fig. 4. Estimated MSE between CV predictions (yHat) and observations (y) versus values of the bandwidth parameter, θ, by environment and model. K _θ is a marker-based RKHS regression with bandwidth parameter θ; Pedigree & Markers K _θ uses pedigree and markers, here, θ is the value of the bandwidth parameter for markers. Pedigree & Markers K _0·25+K ₇ uses pedigree and markers with KA. E1–E4: environments where the lines were evaluated.

Overall, models including marker information had better predictive ability than pedigree-based models. For example, relative to P, using PM_KA yielded decreases in MSE between CV predictions of observations of 20·4, 8·8, 7·0 and 11·0% for E1 through E4, respectively (Table A2 in the Appendix). Thus, it appears that sizable gains in predictive ability can be attained by considering markers and pedigrees jointly, as in PM_KA. These results are in agreement with some empirical studies (e.g. Corrada Bravo et al., 2009; de los Campos et al., 2009b) that provided evidence of a gain in predictive ability by jointly considering markers and pedigree information. However, marker density in this study was relatively low; as marker density increases it is expected that the relative importance of considering pedigree information will decrease (e.g. Calus & Veerkamp, 2007).

As shown in Fig. 4, the value of the bandwidth parameter that gave the best predictive ability was in the range (2,4), except for environment E2 in which values of θ near one performed slightly better. The value of the bandwidth parameter that was optimal from the perspective of predictive ability was similar in M and PM models (Fig. 4 and Table A2 in the Appendix). However, the difference between the predictive ability of PM_{k ,θ} and M_{k ,θ} models was larger for extreme values of θ, indicating that PM models are more robust than M models with respect to the choice of θ. Again, this occurs because PM_{k ,θ} involves some form of KA (between the RK evaluated in the pedigree, A, and the one evaluated in marker genotypes, K_θ).

In all environments, KA had an estimated PMSE that was either close or lower than the one obtained with any specific value of the bandwidth parameter (Fig. 4 and Table A2 in the appendix). This was observed both in models with and without pedigree. These results suggest that KA can be an effective way of choosing the RK. Finally, PM_KA had higher predictive ability than PM_BL; this suggests a superiority of semi-parametric methods. However, PM_BL outperformed PM_{k ,θ} for extreme values of the bandwidth parameter, illustrating, again, the importance of kernel selection. Moreover, the superiority of RKHS methods may not generalize to other traits or populations.

Using data from US Jersey sires (n=1446) genotyped with the BovineSNP50 BeadChip (42 552 Single-nucleotide polymorphisms (SNPs)) de los Campos et al. (2010) compared the predictive ability of several RKHS models for predicted transmitting abilities of milk production, protein content and daughter pregnancy rate. Models evaluated in that study were: (a) BRR, i.e. K=XX′; (b) a Gaussian kernel evaluated over a grid of values of the bandwidth parameter, i.e. K_θ; (c) KA using the two most extreme kernels in the sequence {K_θ}; and (d) a model where K was a marker-based estimate of a kinship matrix, i.e. K=G. Results in that study are in agreement with findings reported here in that using KA gave predictive ability similar to that achieved with best performing kernel in the sequence {K_θ}. The comparison between KA, BRR and using K=G yielded mixed results: for milk yield all models performed similarly; however, for protein content BRR and G outperformed KA and the opposite was observed for daughter fertility, illustrating that the optimal choice of kernel may be trait dependent.

4. Concluding remarks

Incorporating molecular markers into models for prediction of genetic values poses important statistical and computational challenges. Ideally, models for dense molecular markers should be: (a) able to cope with the curse of dimensionality; (b) flexible enough to capture the complexity of quantitative traits and (c) amenable for computations. RKHS regressions can be used to address some of these challenges.

Coping with the curse of dimensionality and with complexity. In RKHS, the curse of dimensionality is controlled by defining a notion of smoothness of the unknown function with respect to pairs of points in input space, Cov[g(t_i), g(t _i′)]∝K(t_i, t _i′). The choice of RK becomes a central element of model specification in RKHS regressions.

As a framework, RKHS is flexible enough to accommodate many non-parametric and some parametric methods, including some classical choices such as the infinitesimal model. The frontier between parametric and non-parametric methods becomes fuzzy; models are better thought as decision rules (i.e. maps from data to estimates) and best evaluated based on performance. Predictive ability appears as a natural choice for evaluating model performance from a breeding perspective.

From a non-parametric perspective, kernels are chosen based on their properties (e.g. predictive ability). To a certain extent, this choice can be made a task of the algorithm. KA offers a computationally convenient method for kernel selection, and results on this study, as well as those of de los Campos et al. (2010), suggests that KA is an effective strategy for kernel selection.

Computational considerations. RK Hilbert spaces methods offer enormous computational advantages relative to most of the parametric methods for regression on molecular markers. This occurs for two reasons: (a) the model can be represented in terms of n unknowns and (b) factorizations such as EV or Singular value decompositions can be used to arrive at highly efficient algorithms. Unfortunately, these benefits cannot be exploited in linear models, y=Xβ+ε, with marker-specific prior precision variances of effects such as BayesA or Bayesian LASSO. This provides RKHS with a great computational advantage relative to those methods, especially when p>>n.

Contribution of marker genotypes to prediction of genetic values. Unlike pedigrees, molecular markers allow tracing Mendelian segregation; potentially, this should allow better predictions of genetic values. Results from this study confirm this expectation. Overall, PM models outperformed P models. Further, most RKHS regression yielded better predictions than those attained with the Bayesian LASSO. However, this did not occur for every RK, indicating that the choice of the kernel is one of the main challenges when applying kernel-based methods. As stated, our results as well as those of de los Campos et al. (2010) suggest that KA provides an effective way of choosing a kernel.

Future challenges. In the kernels used in this study all SNPs contributed equally to the RK. As the number of available markers increases, a high number is expected to be located in regions of the genome that are not associated with genetic variability of a quantitative trait. Ideally, the RK should weight each marker based on some measure of its contribution to genetic variance. In linear models such as the Bayesian LASSO, or methods Bayes A or Bayes B, the prior variances of marker effects, which are marker specific, act as weights assigned to each of the markers (e.g. de los Campos et al., 2009b).

In RKHS models, one could think of kernels where the contribution of each marker to the kernel is weighted according to some measure of its contribution to genetic variance. For example, one could derive weighted estimates of kinship in which each marker obtains a differential contribution. Alternatively, with a Gaussian kernel, one could think of attaching a bandwidth parameter to each marker. For example, one could use , where θ_k and d(x_ik, x_{i ′k}) are a bandwidth parameter and a distance function associated with the kth marker.

An approach similar to that above-described was evaluated by Long et al. (2010) who used radial-basis functions evaluated on principal components (as opposed to individual markers) derived from marker genotypes. Results of that study indicate that the use of input-specific bandwidth parameters may improve predictive ability relative to a model based on a single bandwidth parameter. However, inferring these weights (or bandwidth parameters) poses several statistical challenges when p>>n. This occurs because the kernel must be re-computed every time the bandwidth parameters are updated. A natural alternative is to use two-step procedures, with a first step in which an approximation to the weights (or bandwidth parameters) is employed (e.g. with some form of single-marker regression) and a second step where genetic values are inferred. Irrespective of whether single or two-step approaches are used, the development and evaluation of algorithms for computing weighted kernels seem to constitute a central area of research for the application of RKHS to genomic models.