Genome-enabled prediction of genetic values using radial basis function neural networks

The availability of high density panels of molecular markers has prompted the adoption of genomic selection (GS) methods in animal and plant breeding. In GS, parametric, semi-parametric and non-parametric regressions models are used for predicting quantitative traits. This article shows how to use neural networks with radial basis functions (RBFs) for prediction with dense molecular markers. We illustrate the use of the linear Bayesian LASSO regression model and of two non-linear regression models, reproducing kernel Hilbert spaces (RKHS) regression and radial basis function neural networks (RBFNN) on simulated data and real maize lines genotyped with 55,000 markers and evaluated for several trait–environment combinations. The empirical results of this study indicated that the three models showed similar overall prediction accuracy, with a slight and consistent superiority of RKHS and RBFNN over the additive Bayesian LASSO model. Results from the simulated data indicate that RKHS and RBFNN models captured epistatic effects; however, adding non-signal (redundant) predictors (interaction between markers) can adversely affect the predictive accuracy of the non-linear regression models.


Introduction
The availability of high density panels of molecular markers has catalyzed the adoption of genomic selection (GS) methods in animal and plant breeding (Meuwissen et al. 2001); empirical evidence has demonstrated a superiority of marker-based models over pedigree-based models for predicting complex traits (e.g., VanRaden 2008;Hayes et al. 2009;de los Campos et al. 2009a;Crossa et al. 2010Crossa et al. , 2011. However, most applications of GS use additive linear regression models, and there may be still opportunities for increasing prediction accuracy even further by capturing non-additive sources of genetic variability such as dominance or epistasis. Evidence of epistatic effects in plant traits is vast (Holland 2001(Holland , 2006. For instance, Dudley (2008) found the presence of epistasis in oil, protein, and starch in different crosses of maize lines, and Dudley and Johnson (2010) reported that adding two locus interactions to the model increases prediction power. Despite this, experiments performed in maize have not provided sizable estimates of epistatic variance components (Hallauer and Miranda 1981), perhaps reflecting the fact that even highly epistatic systems generate a great deal of additive variance (e.g., Hill et al. 2008). Also, there is a lack of well-established methods for incorporating epistasis in the prediction of complex traits in plant breeding programs (Hallauer and Miranda 1981;Bernardo 2002).
Interactions between marker alleles at two or more loci can be accommodated in a linear model using appropriate contrasts. However, this is feasible only when the number of markers (p) is moderate. In GS, however, p is usually large, making parametric modeling of complex epistatic interactions unfeasible. An alternative is to use semiparametric regressions (e.g., Gianola et al. 2006), such as kernel-based methods (e.g., Wahba 1978;Gianola et al. 2006; Gianola and van Kaam 2008) or neural networks (NNs) , with the expectation that such procedures can capture complex higher order interaction patterns. The use of reproducing kernel Hilbert spaces (RKHS) for prediction of complex traits was first proposed by Gianola et al. (2006), and empirical studies have demonstrated good prediction accuracy in plant (e.g., Crossa et al. 2010Crossa et al. , 2011 and chicken data (Gonzalez-Recio et al. 2008;Long et al. 2010). However, a potential limitation of RKHS regressions is that the basis functions used for regression must be defined a priori.
In NN, the basis functions are inferred from the data, giving NN great potential for capturing complex interactions between predictor variables (Hastie et al. 2009). Such flexibility comes at the price of a substantial increase in computational demands and the risk of overfitting the training data. Radial basis function neural networks (RBFNNs) are a particular class of NN that have features that make them attractive for applications in GS. First, it has been shown that RBFNNs have universal approximation properties (e.g., Park and Sandberg 1991). Second, RBFNN combines, in a single framework, features of NNs and of RKHS, and both approaches have been widely shown to be promising for predicting phenotypes of complex traits. Further, algorithms exist [e.g., the orthogonal least-squares method proposed by Chen et al. (1991)] that make the computational burden of fitting a RBFNN much less than that of a comparable standard NN.
The RBFNNs have been applied as a prediction and classification tool in many different domains (Jayawardena and Fernando 1998;Takasaki and Kawamura 2007;Zheng et al. 2011). However, they have not been evaluated in the context of genomic selection. In this article, we (1) review the concepts of RBFNN, (2) discuss the connection between these methods and RKHS regressions, and (3) compare the predictive performance of RBFNN with that of RKHS and of an additive linear regression model (Bayesian LASSO). We also illustrate the use of these models on simulated and real maize lines genotyped with high density markers and evaluated for several trait-environment combinations.

Simulated data sets
This data set was simulated by Zhang and Xu (2005) and has a sample size of 600 individuals. The genome has a single chromosome (1,800 cM long) and 121 evenly spaced markers with a 15 cM per marker interval. The authors simulated 9 main QTL effects and 13 interactions between different QTL effects; all QTL effects overlapped with markers. Each QTL had a contribution to phenotypic variance that varied from 0.5-20 %. Models were fitted to two simulated data sets, including the 121 evenly spaced marker covariates indicating the genotype of the jth marker, and the 121(121 ? 1)/2 = 7,381 marker 9 marker first order interactions.

Maize data sets
The maize data represent 21 trait-environment combinations measured in 300 tropical inbred lines genotyped with 55,000 SNPs each. First, we considered eight trait-environment combinations including four traits [grain yield (GY), female flowering (FFL) or days to silking, male flowering time (MFL) or days to anthesis, and anthesissilking interval (ASI)], each evaluated under severe drought stress (SS) and in well-watered (WW) environments. This data set was previously used by Crossa et al. (2010) for the assessment of prediction performance of the BL and RKHS methods, but using only 1,148 SNPs.
Second, the 300 maize lines were evaluated in 9 international environments for gray leaf spot (GLS), a disease caused by the fungus Cercospora zeae-maydis, which is pandemic in Africa. Now recognized as one of the most significant yield-limiting diseases of maize worldwide, GLS is associated with the rapid adoption of conservation agriculture techniques. The 9 environments for GLS had appreciable levels of disease infection. Third, grain yields (GY) of these 300 maize lines were also measured in a large number of relatively high yielding environments (GY-HI) and low yielding environments (GY-LO). Finally, phenotypes for northern corn leaf blight (NCLB), a disease caused by the fungus Exserohilum turcicum, were taken from disease trials evaluated in two environments.

Linear and non-linear regressions on marker genotypes
In GS, phenotypes ðy i ; i ¼ 1; . . .; nÞ are regressed on p marker covariates using a regression function that maps from marker genotypes x ij 2 f0; 1; 2g onto the real line, that is f ðx i1 ; . . .; x ip Þ. Methods differ on (a) how f ðx i1 ; . . .; x ip Þ is structured (e.g., linear vs. non-linear functions of marker genotypes) and (b) how the parameters are estimated. In all models, the response variable was described as the sum of an effect common to all lines (l), a genetic value f(x i ), and a model residual " e i , that is, Residuals were assumed to be independent draws from a normal distribution with null mean and variance equal to r 2 e ni , where n i is defined below. Models differed in how marker information was used to describe f(x i ). Phenotypes were standardized within trait-by-environment combination; therefore, in all cases the response was " y i ¼ 1 SDÂn i P n i k¼1 y ik , where n i is the number of replicates available for the ith trait-by-environment combination, and SD is the sample standard deviation of the within trait-by-environment line means.

Linear model
In linear additive models for GS (e.g., Meuwissen et al. 2001 where b 0 is an intercept and fb j g p j¼1 are marker effects. In practice, the number of markers can vastly exceed the number of records; therefore, shrinkage estimation procedures are commonly used to estimate marker effects. This approach has been used successfully for predicting genetic values in plants and animals. However, the additive specification may not be optimal if dominance or epistasis effects make a sizeable contribution to total genetic variance. As stated, the linear additive model can be extended to accommodate dominance or epistasis by adding the appropriate effects. However, when p is large, modeling complex epistatic patterns using interactions is not feasible. Here, genetic values are represented using a linear regression on marker genotypes and marker effects were estimated using the Bayesian LASSO of Park and Casella (2008), as implemented in the BLR package of R (de los Campos and . Further details about this model and about the algorithms implemented in BLR can be found in de los Campos et al. (2009a) and Pérez et al. (2010). These articles also provide general guidelines for choosing hyper-parameters which were followed here to determine (1) the prior scale (S), (2) the degrees of freedom (df) of the scaled-inverse Chi-square distribution assigned to the residual variance, and (3) the shape (s) and rate (r) parameter of the gamma distribution assigned to the regularization parameter. In our implementation, df was set equal to 4 and the scale was set to 1, this gives a prior density with a prior expectation equal to 0.5 (i.e., one half of the sample variance of the standardized phenotypes) and it is relatively flat around its mode. Pérez et al. (2010) also provide guidelines for choosing the rate and shape parameters of the BL and the proposed approach is to choose these hyper-parameters so that the prior has a mode that is relatively flat in the neighborhood of where MSx is the average (across subjects) sum of squares of marker genotypes. This quantity varies across data sets. Here, we set the rate and shape parameters to 1 9 10 -4 and 0.6, respectively; these values give a prior that has a mode close to 30 and is flat in a relatively wide range of values fork:

Reproducing kernel Hilbert spaces (RKHS) regression
The RKHS model has been suggested as an alternative to the linear model. Its proponents (e.g., Gianola et al. 2006) have argued that such a procedure may capture complex interaction patterns that may not be accounted for in the linear model, and simulations as well as empirical evidence have hinted a superiority of this approach over linear models for predicting phenotypes for some traits (e.g., de los Campos et al. 2009bCrossa et al. 2010). In a RKHS model, the regression function takes the following form: where x i ¼ ðx i1 ; . . .; x ip Þ 0 and x i 0 ¼ ðx i 0 1 ; . . .; x i 0 p Þ 0 are vectors of marker genotypes, a i 0 are regression coefficients, and Kðx i ; x i 0 Þ is a positive definite function (the reproducing kernel, RK) evaluated in a pair of lines which are denoted by i and i 0 . This can be, for example, a Gaussian kernel, where h is a bandwidth parameter and kx i À x i 0 k is the Euclidean distance between the vectors of marker genotypes in lines i and i'. The RK provides a set of n basis functions, fKðx i ; x i 0 Þg n i¼1 , which are non-linear on marker genotypes; however, the regression function is simply a linear combination of the basis functions provided by the RK. To prevent over-fitting, the vector of regression coefficients, ða 1 ; . . .; a n Þ, is estimated using shrinkage estimation procedures such as penalized or Bayesian regressions. Clearly, the set of basis functions is defined a priori via the choice of kernel, and an inappropriate selection may limit the ability of RKHS to capture complex patterns.
As stated above, in this model the regression function is linear on the RK. We used a Gaussian kernel, together with a strategy of kernel averaging (KA, de los Campos et al. 2010), for implicit selection of optimal values of the bandwidth parameter. In particular, we defined three extreme kernels: is a standardized squared Euclidean distance, V j is the sample variance of the jth marker, q 05 is the 5th percentile of d 2 ii 0 , and h 1 ¼ 5; h 2 ¼ 1; h 3 ¼ 1=5 are values of the bandwidth parameter, such that K 1 ðx i ; x i 0 ; h 1 Þ gives extremely local basis functions and K 3 ðx i ; x i 0 ; h 3 Þ gives basis functions with a much wider span. Figure 4 (Appendix 1) gives a histogram (for the ASI-SS maize data set) of the off-diagonal entries of the three kernels. K 1 has very small off-diagonal values, K 2 gives off-diagonal values concentrated in the [0.2, 0.6] interval and K 3 gives off-diagonal values concentrated in the [0.7, 0.9] interval.
Kernel averaging was implemented using Bayesian methods, as described by . The joint prior distribution of this Bayesian RKHS regression has eight hyper-parameters; the prior scale (S) and degrees of freedom (df) of the scaled-inverse Chi-square distribution assigned to the residual variance, and those of the distributions assigned to the variances associated with each of the three RK (the scale and the degrees of freedom hyperparameters). In our implementation, we set the df = 4, because this gives relatively un-informative priors, and chose the scale parameters so that (1) the prior expectation of the residual variance was one half of the sample variance of the standardized phenotypes (in our case S ¼ ðdf À 2Þ=2 ¼ 1) and (2) the prior expectation of the variance of each of the kernels was 1/6 of the sample variance of standardized phenotypes (in our case S ¼ ðdf À 2Þ=6 ¼ 1=3).

Single hidden layer neural network
In a NN, the basis functions are inferred from the data, which give NN great flexibility in terms of capturing complex patterns. The rest of this section gives an overview of these procedures. We begin by reviewing a standard single hidden layer NN for a continuous response with the RBFNN introduced subsequently.
A graphical representation of a single hidden layer neural network is given in Fig. 1. This NN can be thought of as a twostage regression (e.g., Hastie et al. 2009). In the first stage (hidden layer), M data-derived basis functions, fz mi g i¼n;m¼M i¼1;m¼1 , are inferred; in the second stage (the output layer), the response is regressed on the basis functions (inferred in the hidden layer) using a non-linear procedure (Fig. 1).
In the hidden layer, one data-derived predictor (or basis function) is inferred at each of M neurons. These dataderived predictors are formed by first inferring a score (u mi ), which is a linear combination of the input variables (marker genotypes, in our case), and then transforming this score using a non-linear activation function, uðÁÞ, that is Fig. 1 Graphical representation of a single hidden layer feedforward neural network (NN). In the hidden layer, input variables x i ¼ ðx i1 ; . . .; x ip Þ (j ¼ 1; . . .; p markers) are combined using a linear function, u mi ¼ w m0 þ P p j¼1 x ij w mj (m = 1,…,M), and subsequently transformed using a non-linear activation function, u m ðÁÞ, yielding a set of M (M = number of neurons) inferred scores, z mi ¼ u m ðu mi Þ. These scores are used in the output layer as basis functions to regress the response using the linear activation function on the data-derived predictors y i ¼ uðw 0 þ P M m¼1 z mi w m Þ þ e i ; uðÁÞ could be either an identity or any other function where w m0 is an intercept (also referred to as 'bias' term), and w m ¼ fw mj g m¼M;j¼p m¼1;j¼1 is a vector of regression coefficients (the socalled 'weights').
Subsequently, in the output layer, phenotypes are regressed on the data-derived features, fz mi g i¼n;m¼M i¼1;m¼1 , according to y i ¼ uðw 0 þ P M m¼1 z mi w m Þ þ e i , where uðÁÞ is usually a linear activation function and e i is a model residual. For a continuous outcome, uðÁÞ, may simply be an identity link, so that y i ¼ w 0 þ P M m¼1 z mi w m þ e i . Model specification in NN refers to the choice of architecture (i.e., the number of hidden layers and of neurons per hidden layer) and the type of activation function.
The activation function is a monotonic map from a score defined in the real line onto the interval [0, 1] (for a sigmoid function) or onto the interval [-1, 1] (for a hyperbolic tangent function). For example, the sigmoid function is z mi ¼ u m ðu mi ; hÞ ¼ 1 1þexpðÀh u mi Þ , where h is a parameter controlling the shape of the activation function. The use of data-derived predictors and activation functions, together with the possibility of using multiple neurons and layers, gives NN great flexibility for capturing complex interaction patterns between predictors; however, the computational burden can be extremely high and over-fitting may occur.

Radial basis function neural network
The RBFNN was first proposed by Broomhead and Lowe (1988) and Poggio and Girosi (1990), who applied regularization theory to solve ill-conditioned problems in the approximation/interpolation of a function. Figure 2 gives a graphical display of a single hidden layer RBFNN with M neurons (M B n). The output layer is exactly as that shown for NN in Fig. 1; the main difference between the standard NN and an RBFNN is how the hidden layer is structured, that is, how the basis functions are inferred. In an RBFNN, the basis functions consist of a pre-determined number of radial basis functions (RBFs), each of which is indexed by parameters (e.g., centroid; see below for further explanation) to be estimated from the data.
A radial basis function, wðÁÞ, is a map of pairs of vectors, fx i ; cg, onto the real line, with the peculiarity that the map depends only on the Euclidean distance between the two vectors (input vector, x i , and centroid vector, c), that is, wðx i ; cÞ ¼ wðkx i À ckÞ. The Gaussian kernel is a particular case of this. The illustration in Fig. 2 uses a Gaussian RBF; however, the methodology can be applied using other RBFs, such as splines, multi-quadrics, etc. For a given set of centroids fc 1 ; . . .; c M g (M vectors each of order p), the set of parameters involved in an RBFNN (the weights of the output layer, x ¼ fw 0 ; :w 1 :; w M g) can include a large  Fig. 2 Graphical representation of a single hidden layer (Gaussian) radial basis function neural network (RBFNN). In the hidden layer, information from input variables ðx i1 ; . . .; x ip Þ (j ¼ 1; . . .; p markers) is first summarized by means of the Euclidean distance between each of the input vectors {x i } with respect to M (data-inferred) (M = number of neurons) centers {c m }, that is u mi ¼ h m jjx i À c m jj 2 . These distances are then transformed using the Gaussian function, z mi ¼ expðÀu mi Þ, yielding M data-derived scores. These scores are used in the output layer as basis functions for the linear regression, y i ¼ wðw 0 þ P M m¼1 z mi w m Þ þ e i ; wðÁÞ is usually an identity function Theor Appl Genet (2012) 125:759-771 763 number of unknowns; thus, shrinkage estimation methods may be needed. The regularization approach for solving a learning (approximation/interpolation) problem is to search for a function f ðx i ; xÞ that approximates the training set of response data; this function has input vectors, x i 2 R p (the domain of the function), responses y i 2 R; ði ¼ 1; . . .; nÞ, and the weight vector x. In other words, we need to find the functional U½f ðx i ; xÞ that minimizes the cost function H½f ðx i ; xÞ (Poggio and Girosi 1990;Kecman 2001), where P n i¼1 ðy i À f ðx i ; xÞÞ 2 is a residual sum of squares between the response y i and the approximating function f ðx i ; xÞ (i.e., a measure of goodness-of-fit); k is a small, positive number (the Lagrange multiplier), also called the regularization parameter, that controls the trade-off between fitness and model complexity; U½f ðx i ; xÞ is a measure of complexity of f ðx i ; xÞ and a penalty function also called a stabilizer that enforces the smoothness of f ðx i ; xÞ. The regularization parameter k, which is commonly proportional to the extent of noise in data, determines the influence of this stabilizer and controls the trade-off between the two terms of H½f ðx i ; xÞ. The stabilizer (or penalty) function U½f ðx i ; xÞ can take several forms (i.e., spline, multi-quadric, radial basis, Gaussian, etc.).
When U½f ðx i ; xÞ takes a symmetrical radial form, a particular regularized solution that minimizes H½f ðx i ; xÞ is given by the linear combination of the Gaussian RBFs (Poggio and Girosi 1990;Kecman 2001): where w 0 is the intercept, w m are the weights of the linear layer, c m are the centers of the RBFs and w m ðkx i À ckÞ ¼ exp½Àhkx i À c m k 2 are Gaussian RBFs that depend only on the Euclidean norm of the difference vector x i À c m . The weights (w m ), the centroids (c m ), and h are estimated in such a way that the fit between f ðx i ; xÞ and the desired response is optimum.

Estimating the parameters of the RBFNN
To estimate the parameters of a RBFNN, the weights w m of the linear output layer are determined using the ordinary least-squares method, once the RBF (Gaussian in this case) w m ðÁÞ (0 \ m B M), their corresponding centers, and the bandwidth h of the RBF have been assigned values. Several methods are available for selecting the centers (Haykin 1994); in this study, the centroids were selected using the orthogonalization least-squares procedure proposed by Chen et al. (1991). This method sequentially selects the centers of the RBF such that each new selected center is orthogonal to the previous ones. The selected centers maximize the decrease in the mean squared error of the RBFNN, and the algorithm stops when the number of centers attains a desired precision, or when the number of centers is equal to the number of input vectors, that is, when M = n.

Relationship between RBFNN and RKHS
The RBFNN is closely related to RKHS regression. In particular, if in Fig. 2 we let the activation function of the output layer be the identity function wðw 0 þ P M m¼1 z mi w m Þ ¼ w 0 þ P M m¼1 z mi w m and the number of neurons be equal to n, with c m ¼ x i 0 , then the structure of the conditional expectation function of the RKHS regression and the structure of the RBFNNs are exactly the same. In the RBFNN, the strategy is to select a set of basis functions by estimating centers (c m ), and each center defines a basis function. Typically, the number of centers (or neurons, in this case) is much smaller than the number of data points. The strategy in RKHS regression is different: a large set of basis functions is offered to the algorithm (at least n, one per data point, and more, when kernel averaging is used; see Eq. [1]), but the contribution of each of these basis functions to the conditional expectation (i.e., the a's) is estimated using shrinkage estimation procedures. In the statistical learning literature, this is known as 'automatic knot selection' (Ruppert et al. 2003) and is the strategy used by the smoothing spline (Wahba 1990). Arguably, the performance of an RBFNN could be improved if a shrinkage estimation procedure was used, instead of the least-squares method of Chen et al. (1991), but the latter is computationally simpler.

Model comparison
The predictive ability of the additive Bayesian LASSO linear model, RKHS, and the RBFNN was evaluated. A total of 50 independent random partitions of each of the 23 data sets into training (90 % of the data points) and testing (10 % of the data points) were generated. For each of these partitions, models were fitted to the training set data, and prediction accuracy was evaluated in the testing data set. Accuracy was assessed by means of Pearson's correlation between predictions and observations and by the predictive mean squared error (PMSE ¼ n À1 tst P n test i¼1 ð" y i À" y i Þ 2 , where" y i is the predicted value), both evaluated in testing data sets of size n tst .
The number of times a given model had a higher correlation (or smaller PMSE) than another was counted and represented in a graph, to produce a visual assessment of the ''winner'' models in terms of correlation and PMSE.

Results
The average (across 50 training-testing partitions) correlations between predictions and observations obtained with the simulated and real data sets are given in Table 1. Results for PMSE are given in Table 2 (Appendix 2). Given the similarity of results for correlations and PMSE, here we will concentrate on correlations only.

Simulated data sets
The analysis involving 121 marker covariates showed a marked superiority of RKHS (correlation 0.757) and of RBFNN (correlation 0.770) over the Bayesian LASSO (correlation 0.643). Here, RBFNN outperformed RKHS slightly. These results confirm that RKHS and RBFNN are able to capture patterns (perhaps generated by epistatic effects) that cannot be detected by a linear model for additive effects.
However, when marker main effects and two-marker interactions were fitted, the performance of the linear model increased markedly (correlation 0.797) and that of the semi-parametric procedures decreased (average correlation 0.547, for both RKHS and RBFNN). These results indicate that, for non-linear models, the information on interaction between predictors incorporated into the input space becomes redundant in the feature space. On the other hand, for the linear model, the information on the marker 9 marker interaction incorporated in the input space is useful to predict the feature space. The linear model was able to detect this via estimates of regression coefficients which weight the contribution of each marker to the estimated conditional expectation. On the other hand, in both RBFNN and RKHS, each marker gets a similar weight in the basis function or kernel, and the effect of adding non-signal covariates reduces method performance.

Maize data sets
Overall, that is, averaged across training-testing partitions, the three methods performed similarly, with only a slight superiority of RKHS (average correlation 0.553) over RBFNN (average correlation 0.547) and the linear model (average correlation 0.542) (Table 1). Similarly, RKHS had a slightly smaller average PMSE (0.645) than RBFNN (0.656) and BL (0.658) ( Table 2). Figure 3a-b (and Fig. 5a-b in Appendix 2) gives the correlations (and PMSE) obtained with RKHS versus BL, and RBFNN versus BL. In these figures, each dot represents the estimated correlations (and PMSE) for each of the two methods included in the plot and corresponds to one of the 1,050 analyses (21 trait-environment combinations 9 50 training-testing partitions) conducted. A point above the 45°line represents an analysis where the method whose predictive correlation (and PMSE) is given on the vertical axis outperformed the one whose correlation (and PMSE) is given on the horizontal axis. Although there is a slight overall superiority of RKHS and RBFNN over the linear model (they outperformed the linear model 66 and 60 % of the times, respectively; see Table 1), the average performance across traits and environments was similar.

Flowering (FFM, MFL, ASI)
For traits FFL and MFL (Table 1), the three models achieved high prediction accuracy (correlations over 0.75), whereas for ASI they achieved moderate correlations. These results are in agreement with those reported by Crossa et al. (2010) for these traits. For FFL and MFL, the predictive accuracy obtained under well-watered conditions was higher and more stable (across partitions) than that obtained under drought stress. For these traits, we observed, in general, a slight superiority (1-3 % in the correlation) of RKHS or RBFNN over the additive Bayesian LASSO.

Grain yield
For yield traits (Table 1) we obtained moderately high correlations in well-watered (GY-WW) and high-yield environments (GY-HI), and a lower predictive performance under drought stress (GY-SS) and low-yield environments (GY-LO). These results highlight the difficulties of predicting performance under stress conditions and reinforce the importance of having a precise phenotypic system for controlling local plot-to-plot variability in field trials under restrictive conditions. The analysis of GY traits showed slightly better prediction of BL and RKHS over RBFNN.

Gray leaf spot
Estimated predictive correlations ranged from 0.220 to 0.596, depending mostly on environment. Although there were some differences across models, their ranking was not clear; the BL, RKHS and RBFNN methods were best in 4, 3 and 2 of the 9 environments, respectively.

Northern corn leaf blight
The estimated predictive correlations for these trait-environment combinations were moderate to high, and in the two environments we observed better performance of the semi-parametric procedures: RKHS was best in environment 1 and RBFNN was best in environment 2.

Discussion and conclusions
Our empirical results, in which 21 maize data sets represented different traits and environments, indicated that the three models considered had a very similar overall prediction accuracy, with a slight superiority of RKHS and RBFNN over the additive Bayesian LASSO model. In general, these results are similar and sometimes slightly better than other findings using similar data sets. The sample size (300 maize lines) may be a limiting factor for obtaining better discrimination between the predictive accuracy of these models. Results from the simulated data suggest that, for non-linear models, introducing interactions between predictors (markers) in the input space may not be necessary for predicting the feature space; however, this interaction information in the input space is necessary (but feasible to be incorporated in real situations) when the feature space is predicted by means of a linear model. These results were confirmed when using real data.
The simulated data example not only shows that RKHS or RBFNN can capture epistatic patterns, but also indicates that adding non-signal predictors (as might happen using 55 K, 100 K or denser platforms) can adversely affect the predictive accuracy of these models, because in the current formulations of RKHS and RBFNN all markers are equally weighted. Possible ways to overcome this problem would be to (1) introduce unknown marker weights in the kernel, which could be computationally challenging; (2) use arbitrary weights or pre-selecting markers based on an adhoc procedure (e.g., single marker regression or information gain); or (3) obtain haplotypes and examine their prediction accuracy. This is an issue that requires further study.
Given the hundreds of thousands of markers, including all pair-wise (or higher order) interactions among markers in linear models becomes a difficult and almost impossible problem to solve. As pointed out initially by Gianola et al. (2006), and subsequently corroborated by Long et al. (2010), non-parametric models do not impose strong assumptions on the phenotype-genotype relationship and allow capturing interactions among loci. The results of these real data sets, comprising maize trials conducted to measure several traits under a wide range of environmental conditions, agreed with previous findings in animal breeding and with simulated experiments in the sense that sometimes a non-parametric treatment of markers may account for epistatic effects that are not captured by linear additive regression models.
The two kernel models considered, RBFNN and RKHS, had some similarities and displayed good predictive abilities in several trait-environment combinations. While RKHS with kernel averaging is robust for any combination of traits and environments, the two non-parametric models, RBFNN and RKHS, seem to be useful for predicting quantitative traits with complex underlying gene action under varying types of interaction with different environmental conditions. While the additive linear model seems to be robust when hundreds of non-signal predictors are included in the model, the degraded performance of RKHS and RBFNN when a large number of non-signal markers are added to the model requires further investigation, along the previously described lines.
Although parametrically estimating all possible regression coefficients in a linear model is not feasible for large p, it is possible to make further improvements on the accuracy of the RKHS and RBFNN models by introducing differential weights in SNPs, as shown by Long et al. (2011) for RBFs. Further, the output layer of the RBFNN used in this study does not use a regularized regression but, rather, ordinary least squares. Using a shrinkage regression model for the output layer of the RBFNN may offer an extra increase in accuracy. This needs further investigation in the context of genomic prediction.
Open Access This article is distributed under the terms of the Creative Commons Attribution License which permits any use, distribution, and reproduction in any medium, provided the original author(s) and the source are credited.  Fig. 4 Histograms of the offdiagonal entries of each of the three kernels used (K 1 , K 2 , K 3 ) in the RKHS model for the ASI-SS maize data set In a when the best model is RKHS, this is represented by a white circle; when the best model is BL, this is represented by a black circle. In b when best model is RBFNN, this is represented by a white circle; when the best model is BL, this is represented by a black circle