Predicting complex quantitative traits with Bayesian neural networks: a case study with Jersey cows and wheat
 Daniel Gianola^{1, 2, 3},
 Hayrettin Okut^{1, 4}Email author,
 Kent A Weigel^{2} and
 Guilherme JM Rosa^{1, 3}
DOI: 10.1186/147121561287
© Gianola et al; licensee BioMed Central Ltd. 2011
Received: 22 July 2011
Accepted: 7 October 2011
Published: 7 October 2011
Abstract
Background
In the study of associations between genomic data and complex phenotypes there may be relationships that are not amenable to parametric statistical modeling. Such associations have been investigated mainly using singlemarker and Bayesian linear regression models that differ in their distributions, but that assume additive inheritance while ignoring interactions and nonlinearity. When interactions have been included in the model, their effects have entered linearly. There is a growing interest in nonparametric methods for predicting quantitative traits based on reproducing kernel Hilbert spaces regressions on markers and radial basis functions. Artificial neural networks (ANN) provide an alternative, because these act as universal approximators of complex functions and can capture nonlinear relationships between predictors and responses, with the interplay among variables learned adaptively. ANNs are interesting candidates for analysis of traits affected by cryptic forms of gene action.
Results
We investigated various Bayesian ANN architectures using for predicting phenotypes in two data sets consisting of milk production in Jersey cows and yield of inbred lines of wheat. For the Jerseys, predictor variables were derived from pedigree and molecular marker (35,798 single nucleotide polymorphisms, SNPS) information on 297 individually cows. The wheat data represented 599 lines, each genotyped with 1,279 markers. The ability of predicting fat, milk and protein yield was low when using pedigrees, but it was better when SNPs were employed, irrespective of the ANN trained. Predictive ability was even better in wheat because the trait was a mean, as opposed to an individual phenotype in cows. Nonlinear neural networks outperformed a linear model in predictive ability in both data sets, but more clearly in wheat.
Conclusion
Results suggest that neural networks may be useful for predicting complex traits using highdimensional genomic information, a situation where the number of unknowns exceeds sample size. ANNs can capture nonlinearities, adaptively. This may be useful when prediction of phenotypes is crucial.
Background
Challenges in the study of associations between genomic variables (e.g., molecular markers) and complex phenotypes include the possible existence of cryptic relationships that may not be amenable to parametric statistical modeling, as well as the high dimensionality of the data, illustrated by the growing number of single nucleotide polymorphisms, now close to 10 million in humans http://www.genome.gov/11511175. These associations have been investigated primarily using naïve singlemarker regressions and, more recently, with Bayesian linear regression models of various types [1–3] but that assume additive inheritance almost invariably, while typically ignoring interactions and nonlinearity. Taking into account these phenomena may enhance the ability of predicting outcomes, and this is relevant in genomeassisted management of livestock and plants and in individualized medicine.
There has been a growing interest in the use of nonparametric methods for prediction of quantitative traits based on reproducing kernel Hilbert spaces regressions on markers [2, 4–7] and radial basis functions models [8] or related approaches [9]. Artificial neural networks (ANN) provide an interesting alternative because these learning machines can act as universal approximators of complex functions [10, 11]. ANNs can capture nonlinear relationships between predictors and responses and learn about functional forms in an adaptive manner, because a series of transformations called activation functions are driven by parameters. ANNs can be viewed as a computer based system composed of many processing elements (neurons) operating in parallel [12], and also as a schematic of Kolmogorov's theorem for representation of multivariate functions [13]. An ANN is determined by the network structure, represented by the number of layers and of neurons, by the strength of the connections (akin to nonparametric regression coefficients) between inputs, neurons and outputs, and by the type of processing performed at each neuron, represented by a linear or nonlinear transformation: the activation function. Neural networks have the potential of accommodating complex relationships between input and response variables, as well as of difficult to model interactions among inputs. For these reasons, ANNs are interesting candidates for the analysis of complex traits affected by cryptic forms of gene X gene interaction, and many algorithms for training (fitting) such networks are now available [14].
In this study we investigated the performance of several ANN architectures using Bayesian regularization (a method for coping with the "small n, large p" problem that arises in statistical models including a massive number of explanatory variables) when predicting milk production traits in a sample of Jersey cows or mean grain yield in hundreds of inbred wheat lines. The architectures considered differed in terms of number of neurons and activation functions used, and the input (predictor) variables were derived from pedigree and molecular marker information on the corresponding samples. The paper begins with a brief account of Bayesian regularized neural networks, of their connection with linear random regression models often used in quantitative genetics, and of how Bayesian regularization is made. Subsequently, it is shown how a neural network treatment of genomic data can enhance predictive ability over and above that using pedigree information (in Jerseys) or linear Bayesian regression on markers (in both cows and wheat), which is representative of a standard approach in quantitative genomics.
Methods
For clarity of presentation the methodology is presented first, as the main objective of the paper was to cast neural networks in a quantitative genetics predictive context. Subsequently, a description of the two sets of data used to illustrate how the Bayesian neural networks were run is provided. As stated, the first data set consisted of milk, protein and fat yield in dairy cows. The second set represented 599 lines of wheat, with mean grain yield as target trait.
Excursus: FeedForward Neural Networks
where e_{ i } ~ (0, σ^{2}) and σ^{2} is a variance parameter. If g(.) is a linear or identity activation function, the model is a linear regression on the adaptive covariates f_{ k }(w'_{ k } p_{ i }); if, further, f_{ k }(.), is also linear, the regression model is entirely linear. The term "adaptive" means that the covariates are functions of unknown parameters, the {w_{ kj }} connection strengths, so the networks can "learn" the relationship between explanatory variables and phenotypes, as opposed to posing it arbitrarily, as it is the case in standard regression models. In this manner, this type of neural network can also be viewed as a regression model, but with the extent of nonlinearity dictated by the type of activation functions used. Since the number of parameters increases linearly with the number of neurons, and the number of predictors given by the length of p (e.g., the number of markers) can amply exceed sample size, it is necessary to treat the connection strengths as random effects in which case the Bayesian connection is immediate [15, 16]. This approach is called "Bayesian regularization".
Fisher's infinitesimal model viewed as a neural network
 I)
t = u + e = Cz σ_{ u } + e = Cu* + e,
where z is a vector of independent standard normal deviates, u* = z σ_{ u } ~ (0, I σ^{2}_{ u }) and e ~ (0, I σ^{2}) is a residual vector with σ^{2} interpretable as environmental variance.
II) t = AA^{1}u + e = Au** + e,
where u** = A^{1}u ~ (0, A^{1} σ^{2}_{ u }), and
III) t = A^{1}Au + e = A^{1}u*** + e,
where u*** = Au ~ (0, A^{3}σ^{2}u).
Here, a bias parameter b is included for the sake of generality. Hence, the additive model can be viewed as a singleneuron network regression on either elements of the Cholesky decomposition of the numerator relationship matrix, on the relationships themselves or on the elements of the inverse of A, with the strengths of the connections represented by the corresponding entries of u*, u** and u***, respectively.
Here, the inputs are entries a_{ ij } of the relationship matrix, connecting individual i to all other individuals in the genealogy; the ${u}_{j}^{**\left[k\right]}$ coefficient is the connection strength for input j in neuron k; b_{ k } is the bias parameter associated with neuron k; g_{ k } is an activation function peculiar to neuron k; w_{ k } is the connection strength between the activated emission from neuron k and the output layer, b is the outer layer bias parameter and g(.) is the outer activation function, which may be linear or nonlinear, although it is typically taken as linear for quantitative responses. The nonlinear transformations modify the connection strengths between additive relationships and phenotypes in an adaptive manner, underlining the potential for an improvement in predictive ability.
Given the availability of dense markers in humans and animals, an alternative or complementary source of input that can be used in equation (2) consists of the elements of a markerbased relationship matrix, as in [17]; in this case the a_{ ij } coefficients are replaced by g_{ ij }, i.e., elements of some genome or markerderived relationship matrix G. As noted by [2], when G is proportional to XX', where X is the incidence matrix of a linear regression model on markers, this is equivalent to Bayesian ridge regression. Of course, nothing precludes using both pedigreederived and markerderived inputs in the construction of a neural network.
Bayesian regularization
The objective in ANNs is to arrive at some configuration that fits the training data well but that it also has a reasonable ability of predicting yet to be seen observations. This can be achieved by placing constraints on the size of the network connection strengths, e.g., via shrinkage, and the process is known as regularization. A natural way of attaining this compromise between goodness of fit and predictive ability is by means of Bayesian methods [2, 11, 15, 18]. In this section, an approach used often for Bayesian regularization in neural networks [18, 19] is presented along the lines of the hierarchical models employed by quantitative geneticists [15].
which does not have closed form, because of nonlinearity. Recall that b can be set to 0 provided the observations are suitably centered.
and E_{ w } = w'w. It follows that maximizing $L\left(w\phantom{\rule{2.77695pt}{0ex}}D,{\sigma}^{2},{\sigma}_{w}^{2},M\right)$ is equivalent to minimizing F(α, β). This function is often referred to as a "penalized" sum of squares, and it embeds a compromise between goodness of fit, given by the sum of squares of the residuals E_{ D }, and the degree of model complexity, given by the sum of squares of the network weights E_{ w }. The value w = w^{MAP} that maximizes $L\left(\mathsf{\text{w}}\phantom{\rule{2.77695pt}{0ex}}D,{\sigma}^{2},{\sigma}_{w}^{2},M\right)$ is the mode of the conditional (given the variances) posterior density of the connection strengths; MAP stands for "maximum a posteriori".
If the additive infinitesimal model is represented as a neural network, the coefficient of heritability is given by ${h}^{2}=\frac{1}{2\alpha}\u2215\left(\frac{1}{2\alpha}+\frac{1}{2\beta}\right)=\frac{\beta}{\alpha +\beta}$. As it can be seen in equation (6), if α<<β, the fitting or training algorithm places more weight on goodness of fit. If α>>β, the algorithm emphasizes reduction in the values of w (shrinkage), which produces a less wiggly output function [20].
Given α and β, the w = w^{MAP} estimates can be found via any nonlinear maximization algorithm as in, e.g., the threshold and survival analysis models of quantitative genetics [21].
Tuning parameters α and β
A standard procedure used in neural networks (and in the software employed here) infers α and β by maximizing the marginal likelihood of the data in equation (5); this corresponds to what is often known as empirical Bayes. Because (5) does not have a closed form (except in linear neural networks), the marginal likelihood is approximated using a Laplacian integration done in the vicinity of the current value w = w^{MAP}, which depends in turn on the values of the tuning parameters at which the expansion is made. This type of approach for nonlinear mixed models has been used in animal breeding for almost two decades [22].
These expressions, as well as (7), are similar to those that arise in maximum likelihood estimation of variance components [24–26], which vary depending on the algorithm used. Since β is a positive parameter, it must be that $n>m2{\alpha}_{MAP}tr{H}_{MAP}^{1}$. The quantity $\gamma =m2{\alpha}_{MAP}tr{H}_{MAP}^{1}$ is called "effective number of parameters" in the network [20] and its value ranges from 0 (or 1, if an overall intercept b is fitted) to m, the total number of connection coefficients and bias parameters in the network. If γ is close to n overfitting results, leading to poor generalization. It follows that the computations are similar to those used in the linear and nonlinear models employed by quantitative geneticists, the salient aspect being that a neural network can be strongly nonlinear.
More details on computing procedures for neural networks are in [12, 14, 18, 23, 27, 28]. Typically, an algorithm proceeds as follows: 1) initialize α, β and all elements in w. 2) Take one step of the LevenbergMarquardt algorithm to minimize the objective function F(α, β) and find the current value of w. 3) Compute γ, the effective number of parameters, using the GaussNewton approximation to the Hessian matrix in the LevenbergMarquardt training algorithm. 4) Compute updates α_{ new } and β_{ new }; and 5) iterate steps 24 until convergence.
Neural Network Architectures Evaluated and Implementation
where the $\u0109$ coefficients are the estimated linear regressions of the traits on the activated emissions fired by each of the two neurons.
MATLAB's neural networks toolbox [29] was used for fitting the architectures studied, using Bayesian regularization in all cases. As mentioned earlier, two combinations of activation functions were tried: 1) the hyperbolic tangent sigmoid function for activating emissions from each neuron in the hidden layer, plus a linear activation function from the hidden to the output layer, and 2) a linear activation throughout, this corresponding to a linear model with random regression coefficients. To avoid spurious effects caused by starting values in each iterative sequence, the networks were trained 20 times in the Jersey data and 50 times in the wheat data set, for each of the architectures. In Jerseys, each run randomly allocated 60% of the animals to a training set, 20% to a validation set and 20% to a testing set; results reported are the average of the 20 runs for each of the configurations. In wheat, the records were randomly partitioned into a training (480 lines) and a testing (119 lines) set. Each of the 50 random repeats matched exactly those of [28], to provide a comparison with the predictive ability of the Bayesian Lasso and of support vector regression models used with the wheat data set by [30].
The neural networks were fitted to data in the training set, with the α and β parameters, connection strengths and biases modified iteratively. In the Jersey data, as parameters changed in the course of training, the predictive ability of the network was gauged in parallel in the validation set, which was expected to be similar in structure to the testing set, because they were randomly constructed. The same was done with the wheat data, except that there was no "intermediate" validation set. Once the mean squared error of prediction reached an optimal level, training stopped, and this led to the best estimates of the network coefficients. This estimated network was then used for predicting the testing set; predictive correlations (Pearson) and meansquared errors were evaluated.
MATLAB uses the LevenbergMarquardt algorithm (based on linearization) for computing the posterior modes in Bayesian regularization, and backpropagation is employed to minimize the penalized residual sum of squares. The maximum number of iterations (called epochs) in backpropagation was set to 1000, and iteration stopped earlier if the gradient of the objective function was below a suitable level or when there were obvious problems with the algorithm [28, 29, 31]. Each of these settings corresponds to the default values provided by MATLAB.
Jersey cows data
where p_{ ij } = a_{ ij } (additive relationship between cows i and j) or g_{ ij } (genomederived relationships). Thus, for each cow the input vector p_{ i } had order 297 × 1.
The matrix Z = ME contains "centered" codes, such that the mean of the values in any of its columns is null; Z can be used as in incidence matrix in marker assisted regression models [17, 32, 33]. Then $Z{Z}^{\prime}=G\times 2\sum _{i=1}^{35,798}{\mu}_{i}\left(1{\mu}_{i}\right)$ is interpretable as a covariance matrix, analogous to $A{\sigma}_{u}^{2}$ in the infinitesimal model. The term $2\sum _{i=1}^{35,798}{\mu}_{i}\left(1{\mu}_{i}\right)$ holds under linkage equilibrium only, and cannot be construed as additive genetic variance of marker effects in the classical sense of [33]; its relationship to additive genetic variance in a finite locus or infinitesimal model is tenuous [16, 34].
Wheat lines data
There were 599 wheat lines, each genotyped with 1279 DArT markers (Diversity Array Technology) generated by Triticarte Pty. Ltd. (Canberra, Australia; http://www.triticarte.com.au). The DArT markers may take on two values, denoted by their presence or absence. In this data set, the overall mean frequency of the allele coded as "1" was 0.5607, with a minimum of 0.0083 and a maximum of 0.9866. Markers with a minor allele frequency lower than 0.05 were removed. Missing genotypes at locus j of line i were imputed using samples from the marginal distribution of marker genotypes, that is, ${x}_{ij}~Bernoulli\left({\widehat{p}}_{j}\right)$, where ${\widehat{p}}_{j}$ is the estimated allele frequency computed from the nonmissing genotypes [34]. The phenotype studied was average grain yield of each line. The data came from the International Maize and Wheat improvement Center, Mexico, and it can be downloaded from R package BLR http://cran.rproject.org/web/packages/BLR/index.html; more information can be found in [30, 35, 36]. The wheat data was partitioned randomly into a training set (480 lines) and a test set (119 lines), exactly as in [30].
Results
Degree of complexity
Effective number of parameters (± standard errors), by trait, in Jerseys.^{1}
Network  Fat yield (pedigree)  Fat yield (genomic)  Milk yield (pedigree)  Milk yield (genomic)  Protein yield (pedigree)  Protein yield (genomic) 

Linear  123 ± 5.6  166 ± 2.0  124 ± 7.6  162 ± 2.9  118 ± 8.5  151 ± 4.5 
1 neuron  91 ± 4.9  142 ± 2.0  93 ± 5.8  166 ± 2.0  91 ± 10.3  144 ± 2.5 
2 neurons  104 ± 5.8  128 ± 7.6  122 ± 6.5  145 ± 7.8  114 ± 8.0  136 ± 8.0 
3 neurons  107 ± 5.8  132 ± 5.7  123 ± 5.1  129 ± 6.0  126 ± 6.9  141 ± 4.9 
4 neurons  108 ± 5.8  129 ± 4.7  112 ± 4.7  131 ± 5.8  129 ± 5.4  138 ± 6.0 
5 neurons  106 ± 4.9  127 ± 4.9  118 ± 4.8  132 ± 5.4  131 ± 4.9  138 ± 5.6 
6 neurons  114 ± 3.3  128 ± 7.5  122 ± 5.1  132 ± 5.6  136 ± 4.6  137 ± 5.0 
Effective number of parameters, predictive correlations, and mean squared errors of prediction: wheat.^{1}
ANN architectures  Linear  1 neuron  2 neurons  3 neurons  4 neurons 

Criterion  
Effective number of parameters  299 ± 5.5  260 ± 6.1  253 ± 5.9  238 ± 5.5  220 ± 2.8 
Correlations in testing set  0.48 ± 0.03  0.54 ± 0.03  056 ± 0.02  0.57 ± 0.02  0.59 ± 0.02 
Mean squared error in testing set  0.99 ± 0.04  0.77 ± 0.03  0.74 ± 0.03  0.71 ± 0.02  0.72 ± 0.02 
The effective number of parameters behaved differentially with respect to model architecture and this depended on the input variables used. When using pedigrees in the Jersey data, the hyperbolic tangent activation function in the 1neuron model reduced γ drastically, relative to the linear model (1 neuron with linear activation throughout). Then, an increment in number of neurons from 2 to 6 increased model complexity relative to that of the 1 neuron model with nonlinear activation, but not beyond that attained with the linear model, save for protein yield. For this trait, γ was 118 for the linear model, and ranged from 126 to 136 in models with 3 through 6 neurons. When genomic relationships were used as inputs, γ was largest for the linear model for fat and protein yield, and for the 1neuron model with a nonlinear activation function in the case of milk yield. In wheat, the effective number of parameters decreased as architectures became more complex. There was large variation among runs in effective number of parameters for both data sets, but there was not a clear pattern in the variability.
Predictive performance
Prediction mean squared errors (± standard errors) by trait: Jerseys.^{1}
Network  Fat yield (pedigree)  Fat yield (genomic)  Milk yield (pedigree)  Milk yield (genomic)  Protein yield (pedigree)  Protein yield (genomic) 

Linear  1.19 ± 0.07  0.86 ± 0.05  1.09 ± 0.05  0.88 ± 0.04  1.00 ± 0.04  0.75 ± 0.07 
1 neuron  1.01 ± 0.04  0.74 ± 0.03  0.99 ± 0.04  0.81 ± 0.03  0.97 ± 0.04  0.71 ± 0.04 
2 neurons  0.93 ± 0.05  0.70 ± 0.03  0.96 ± 0.05  0.76 ± 0.04  1.02 ± 0.04  0.72 ± 0.04 
3 neurons  0.92 ± 0.04  0.71 ± 0.03  0.98 ± 0.02  0.78 ± 0.04  0.96 ± 0.06  0.80 ± 0.04 
4 neurons  0.99 ± 0.04  0.84 ± 0.04  0.98 ± 0.04  0.72 ± 0.04  0.90 ± 0.06  0.70 ± 0.03 
5 neurons  0.99 ± 0.04  0.86 ± 0.04  1.00 ± 0.05  0.80 ± 0.04  0.93 ± 0.04  0.77 ± 0.04 
6 neurons  0.95 ± 0.03  0.77 ± 0.04  1.02 ± 0.05  0.79 ± 0.03  0.95 ± 0.03  0.76 ± 0.05 
Correlation coefficients (± standard errors) in the Jersey testing data set, by trait.^{1}
Pedigree relationships  Genomic relationships  

Network  Fat yield  Milk yield  Protein yield  Fat yield  Milk yield  Protein yield 
Linear  0.11 ± 0.04  0.07 ± 0.03  0.02 ± 0.02  0.43 ± 0.02  0.42 ± 0.03  0.44 ± 0.02 
1 neuron  0.23 ± 0.03  0.10 ± 0.03  0.09 ± 0.02  0.51 ± 0.02  0.45 ± 0.02  0.44 ± 0.02 
2 neurons  0.22 ± 0.03  0.08 ± 0.01  0.08 ± 0.03  0.49 ± 0.02  0.46 ± 0.03  0.51 ± 0.02 
3 neurons  0.22 ± 0.02  0.13 ± 0.02  0.10 ± 0.03  0.53 ± 0.02  0.52 ± 0.02  0.47 ± 0.02 
4 neurons  0.20 ± 0.02  0.09 ± 0.02  0.14 ± 0.02  0.45 ± 0.03  0.52 ± 0.02  0.47 ± 0.03 
5 neurons  0.23 ± 0.02  0.13 ± 0.02  0.15 ± 0.02  0.42 ± 0.03  0.50 ± 0.02  0.47 ± 0.02 
6 neurons  0.27 ± 0.02  0.10 ± 0.03  0.11 ± 0.02  0.48 ± 0.04  0.54 ± 0.02  0.50 ± 0.03 
The predictive correlations in wheat (Table 2) ranged from 0.48 with the linear ANN (equivalent to Bayesian ridge regression) to 0.59 for the nonlinear ANN with 4 neurons. Clear and significant differences between linear and nonlinear architectures were observed. The improvements over the linear ANN were 11.2, 14.3, 15.8 and 18.6% in predictive correlation for 1, 2, 3 and 4 neurons in the hidden layer, respectively. Mean squared errors were also 2329% smaller than in the linear ANN.
In the Jerseys, the large variability in mean squared error among runs (Table 3) precludes making strong statements about differences among architectures. However, predictive correlations (Table 4) were clearly larger for the nonlinear ANN. For fat yield, the results with pedigrees employed as input suggest that a nonlinear, adaptive use of additive relationships (as done in all networks with the hyperbolic tangent activation function) can improve predictive performance beyond that of the infinitesimal model. Further, use of genomic relationships led to more reliable prediction of phenotypes than use of pedigree information as measured by the predictive correlations in Table 4. The relative increase in strength of association, as measured by the correlation, is much larger than the ones that have been reported, e.g., in dairy cattle [32, 37], when prediction of breeding values of bulls was made from genomic information, as opposed to from pedigrees. Our result is encouraging and suggests that genomic data may also play an important role in prediction of individual outcomes (as opposed to breeding value), an issue of relevance in medicine [4].
Shrinkage
Discussion
Models for prediction of fat, milk and protein yield in cows using pedigree and genomic relationship information as inputs, and wheat yield using molecular markers as predictor variables were studied. This was done using Bayesian regularized neural networks, and predictions were benchmarked against those from a linear neural network, which is a Bayesian ridge regression model. In the wheat data, the comparison was supplemented with results obtained by our group using RKHS or support vector methods. Different network architectures were explored by varying the number of neurons, and using a hyperbolic tangent sigmoid activation function in the hidden layer plus a linear activation function in the output layer. This combination has been shown to work well when extrapolating beyond the range of the training data [36]. The choice of number of neurons can be based on crossvalidation, as in the present data, or on standard Bayesian metrics for model comparison [11, 15].
The LevenbergMarquardt algorithm, as implemented in MATLAB, was adopted to optimize weights and biases, as this procedure has been effective elsewhere [38]. Bayesian regularization was adopted to avoid overfitting and to improve generalization, and crossvalidation was used to assess predictive ability, as in [28, 39].
Because Bayesian neural networks reduce the effective number of weights relative to what would be obtained without regularization, this helps to prevent overfitting [40]. For the networks we examined, even though the total number of parameters, e.g., in Jerseys, ranged from 300 to 1795, the effective number of parameters varied from only 91 to 136 when pedigree relationships were used, and from 127 to 166 when genomic relationships were used as inputs, illustrating the extent of regularization. There were differences in predictive abilities of different architectures but the small sample used dictates a cautious interpretation. Nevertheless, the results seem to support networks with at least 2 neurons, which has been observed in several studies [20, 28, 41–43]. This suggests that linear models based on pedigree or on genomic relationships may not provide an adequate approximation to genetic signals resulting from complex genetic systems. Because highly parameterized models are penalized in the Bayesian approach, we were able to explore complex architectures. However, there was evidence of overfitting in the Jersey training set, where correlations between observed and predicted data in the training set were always larger than 0.90, sometimes exceeding 0.95. This explains why correlations were much lower in the testing set, which is consistent with what was observed in other studies with neural networks [42]. Although more parameters in a model can lead to smaller error in the training data, it cannot be overemphasized that this is not representative of prediction error in an independent data set, as shown by [43] working with human stature.
Our results with ANN for wheat are at least as good as those obtained with the same data in two other studies. Crossa et al. [35] found crossvalidation correlations with the following values when various methods were used: pedigree information (BLUP), 0.45; pedigreebased reproducing kernel Hilbert spaces regression (RKHS), 0.60; RKHS with both pedigree and markers, 0.61; Bayesian Lasso with markers, 0.46; Bayesian Lasso with markers and pedigree, 0.54, and Bayesian ridge regression on markers, 0.49. Long et al., [30] compared the Bayesian Lasso with four supportvector regression models consisting of the combination of two kernels and two loss functions. The predictive correlation for wheat yield (average of 50 random repeats of the crossvalidation exercise) was 0.52 for the Bayesian Lasso, and ranged between 0.50 and 0.58 for the support vector implementations. Hence, it seems clear, at least for wheat yield in this data set, that the nonparametric methods can outperform a strong learner, the Bayesian Lasso, and that the neural networks are competitive with the highly regarded support vector methods [11].
Thus, the so defined breeding value of individual i depends on the values of the input covariates observed on this individual, on all connection strengths and bias parameters from inputs to neurons in the middle layer (the u's and the b's), and on all connection strengths from the middle layer to the output layer (the w's). In order to obtain an estimate of breeding value the unknown quantities would replaced by the corresponding maximum a posteriori (MAP) estimates or, say, by the estimate of their posterior expectation if a Markov chain Monte Carlo scheme is applied [44].
Another issue is that of assessing the importance of an input relative to that of others. For example, in a linear regression model on markers, one could use a point estimate of the substitution effect or its "studentized" value (i.e., the point estimate divided by the corresponding posterior standard deviation), or some measure that involves estimates of substitution effects and of allelic frequencies. A discussion of some measures of relative importance of inputs in an ANN is in [28, 43], for example, the ratio between the absolute value of the estimate of a given connection strength, and the sum of the absolute values of all coefficients in the network.
Conclusion
Nonlinear neural networks tended to outperform benchmark linear models in predictive ability, and clearly so in the wheat data. Bayesian regularization allowed estimation of all connection strengths even when n<<p, and the effective number of parameters was much smaller than the corresponding nominal number. Although the study was based on small samples, and the differences found may be reflective of random variation, especially in the Jersey data, our results suggest that the neural networks may be useful for predicting complex traits using highdimensional genomic information, a situation where the number of coefficients that need to be estimated exceeds sample size. Neural networks have the ability of capturing nonlinearities, and do so adaptively, which may be useful in the study of quantitative traits under complex gene action, and particularly when prediction of outcomes is crucial, such as in personalized medicine.
In summary, predictive ability seemed to be enhanced by use of Bayesian neural networks. Due to small sample sizes no claim is made about superiority of any specific nonlinear architecture. As it has been observed in many studies, the superiority of one predictive model over another depends on the species, trait and environment, and the same will surely hold for ANNs.
Abbreviations
 ANN:

artificial neural network
 BR:

Bayesian regularization
 BRANN:

Bayesian regularization artificial neural network
 LASSO:

Least absolute shrinkage and selection operator
 MAP:

Maximum a posterior
 NN:

Neural network
 RKHS:

Reproducing kernel Hilbert spaces regression
 SNP:

Single nucleotide polymorphism
 SSE:

Sum of squared error.
Declarations
Acknowledgements
Research was supported by the Wisconsin Agriculture Experiment Station and by grants from Aviagen, Ltd., Newbridge, Scotland, and Igenity/Merial, Duluth, Georgia, USA.
Authors’ Affiliations
References
 Meuwissen THE, Hayes BJ, Goddard ME: Prediction of total genetic value using genomewide dense marker maps. Genetics. 2001, 157: 18191829.PubMed CentralPubMed
 de los Campos G, Gianola D, Rosa GJM: Reproducing kernel Hilbert spaces regression: A general framework for genetic evaluation. Journal of Animal Science. 2009, 87: 18831887. 10.2527/jas.20081259.View ArticlePubMed
 de los Campos G, Gianola D, Allison DB: Predicting genetic predisposition in humans: the promise of whole genome markers. Nature Reviews Genetics. 2010, 11: 880886. 10.1038/nrg2898.View ArticlePubMed
 de los Campos G, Gianola D, Rosa GJM, Weigel KA, Crossa J: Semiparametric genomicenabled prediction of genetic values using reproducing kernel Hilbert spaces methods. Genetics Research. 2010, 92: 295308. 10.1017/S0016672310000285.View ArticlePubMed
 Gianola D, Fernando RL, Stella A: Genomic assisted prediction of genetic value with semiparametric procedures. Genetics. 2006, 173: 17611776. 10.1534/genetics.105.049510.PubMed CentralView ArticlePubMed
 Gianola D, van Kaam JBCHM: Reproducing kernel hilbert spaces regression methods for genomic assisted prediction of quantitative traits. Genetics. 2008, 178: 22892303. 10.1534/genetics.107.084285.PubMed CentralView ArticlePubMed
 Gianola D, de los Campos G: Inferring genetic values for quantitative traits nonparametrically. Genetics Research. 2008, 90: 525540. 10.1017/S0016672308009890.View ArticlePubMed
 Long N, Gianola D, Rosa GMJ, Weigel KA, Kranis A, GonzálezRecio O: Radial basis function regression methods for predicting quantitative traits using SNP markers. Genetics Research. 2010, 92 (3): 20925. 10.1017/S0016672310000157.View ArticlePubMed
 Ober U, Erbe M, Long N, Porcu E, Schlather M, Simianer H: Predicting genetic values: a kernelbased best linear unbiased prediction with genomic data. Genetics. 2011, 188 (3): 695708. 10.1534/genetics.111.128694.PubMed CentralView ArticlePubMed
 Alados I, Mellado JA, Ramos F, AladosArboledas L: Estimating UV Erythema1 irradiance by means of neural networks. Photochemistry and Photobiology. 2004, 80: 351358. 10.1562/20040312RA111.1.View ArticlePubMed
 Bishop CM: Pattern Recognition and Machine Learning. 2006, Singapore: Springer
 Lamontagne L, Marchand M: Advances in Artificial Intelligence. 2006, Berlin: SpringerView Article
 Pereira BDB, Rao CR: Data Mining using Neural Networks: A Guide for Statisticians. 2009, [http://www.po.ufrj.br/basilio/publicacoes/livros/2009_datamining_using_neural_networks.pdf]
 Lampinen J, Vehtari A: Bayesian approach for neural networks review and case studies. Neural Networks. 2001, 14: 257274. 10.1016/S08936080(00)000988.View ArticlePubMed
 Sorensen D, Gianola D: Likelihood, Bayesian and MCMC methods in quantitative genetics. 2002, New York: SpringerView Article
 Gianola D, de los Campos G, Hill WG, Manfredi E, Fernando R: Additive genetic variability and the Bayesian alphabet. Genetics. 2009, 183: 347363. 10.1534/genetics.109.103952.PubMed CentralView ArticlePubMed
 Van Raden PM: Efficient methods to compute genomic predictions. J Dairy Sci. 2008, 91: 44144423. 10.3168/jds.20070980.View Article
 MacKay DJC: Baysian interpolation. Neural Computation. 1992, 4: 415447. 10.1162/neco.1992.4.3.415.View Article
 Titterington DM: Bayesian methods for neural networks and related models. Statistical Science. 2004, 19: 128139. 10.1214/088342304000000099.View Article
 Foresee FD, Hagan MT: GaussNewton approximation to Bayesian learning. Proc IEEE Int Conf Neural Networks. 1997, 19301935.
 Gianola D: Inferences from mixed models in quantitative genetics. Handbook of Statistical Genetics. Edited by: Balding DJ, Bishop M, Cannings C. 2007, West Sussex UK: John Wiley & Sons, Third
 Tempelman RJ, Gianola D: Marginal maximum likelihood estimation of variance components in Poisson mixed models using Laplace integration. Genetics, Selection, Evolution. 1993, 25: 305319. 10.1186/12979686254305.PubMed CentralView Article
 Xu M, Zengi G, Xu X, Huang G, Jiang R, Sun W: Application of Bayesian regularized BP neural network model for trend analysis, acidity and chemical composition of precipitation in North. Water, Air, and Soil Pollution. 2006, 172: 167184. 10.1007/s1127000590688.View Article
 Smith SP, Graser HU: Estimating variance components in a class of mixed models by restricted maximum likelihood. J Dairy Sci. 1986, 69: 1165
 Graser HU, Smith SP, Tier B: A derivativefree approach for estimating variance components in animal models by restricted maximum likelihood. J Anim Sci. 1987, 64: 1362
 Hassami M, Anctil F, Viau AA: Selection of an artificial neural network model for the postcalibration of weather radar rainfall estimation. Journal of Data Science. 2004, 220: 107124.
 MacKay JCD: Information Theory, Inference and Learning Algorithms. 2008, Cambridge; Cambridge University Press
 Okut H, Gianola D, Rosa GJM, Weigel KA: Prediction of body mass index in mice using dense molecular markers and a regularized neural network. Genetics Research. 2011, 93: 189201. 10.1017/S0016672310000662.View ArticlePubMed
 Beal MH, Hagan MT, Demuth HB: Neural Network Toolbox' 6 User's Guide. 2010, The MathWorks, Inc
 Long N, Gianola D, Rosa GJM, Weigel KA: Application of support vector regressions to genomeassisted prediction of quantitative traits. Theoretical and Applied Genetics. 2011, (under review)
 Haykin S: Neural Networks: Comprehensive Foundation. 2008, New York USA: PrenticeHall
 Habier D, Fernando RL, Dekkers JCM: The impact of genetic relationship information on genomeassisted breeding values. Genetics. 2007, 177 (4): 23892397.PubMed CentralPubMed
 Van Raden PM, Tooker ME, Cole JB: Can you believe those genomic evaluations for young bulls?. Journal of Dairy Science. 2009, 92 (ESuppl 1): 314
 Falconer DS, McKay TFC: Introduction to Quantitative Genetics. 1996, Malaysia: Longmans Green
 Crossa J, de los Campos G, Perez P, Gianola D, Dreisigacker S, Burgueño J, Araus JL, Makumb D, Yan J, Singh R, Arief V, Banzinger M, Braun HJ: Prediction of genetic values for quantitative traits in plant breeding using pedigree and molecular markers. Genetics. 2010, 186: 713724. 10.1534/genetics.110.118521.PubMed CentralView ArticlePubMed
 Perez P, de los Campos G, Crossa J, Gianola D: Genomicenabled prediction based on molecular markers and pedigree using the Bayesian Linear Regression package in R. The Plant Genome. 2010, 3: 106116. 10.3835/plantgenome2010.04.0005.PubMed CentralView ArticlePubMed
 Hayes BJ, Bowman BJ, Chamberlain AJ, Goddard ME: Invited review: Genomic selection in dairy cattle: Progress and challenges. J Dairy Sci. 2009, 92: 433443. 10.3168/jds.20081646.View ArticlePubMed
 Maier HR, Dandy CG: Neural networks for the prediction and forecasting of water resources variables: a review of modelling issues and applications. Environmental Modelling & Software. 2000, 15: 101124. 10.1016/S13648152(99)000079.View Article
 Demuth H, Beale M, Hagan M: Neural Network Toolbox™ 6 User's Guide. 2009, The MathWorks, Inc
 Fernandez M, Caballero J: Ensembles of Bayesianregularized genetic neural networks for modeling of acetylcholinesterase inhibition by huprines. Chem Biol Drug Des. 2006, 68: 201212. 10.1111/j.17470285.2006.00435.x.View ArticlePubMed
 Winkler DA, Burden FR: Modelling bloodbrain barrier partitioning using Bayesian neural nets. Journal of Molecular Graphics and Modelling. 2004, 22: 499505. 10.1016/j.jmgm.2004.03.010.View ArticlePubMed
 Joseph H, Huang WL, Dickman M: Neural network modelling of coastal algal blooms. Ecol Model. 2003, 159: 179201. 10.1016/S03043800(02)002818.View Article
 Sorich MJ, Miners JO, Ross AM, Winker DA, Burden FR, Smith PA: Comparison of linear and nonlinear classification algorithms for the prediction of drug and chemical metabolism by human UDPGlucuronosyltransferase isoforms. J Chem Inf Comput Sci. 2003, 43: 20192024. 10.1021/ci034108k.View ArticlePubMed
 Makowski R, Pajewski NM, Klimentidis YC, Vazquez AI, Duarte CW, Allison DA, de los Campos G: Beyond missing heritability: prediction of complex traits. PLOS Genetics. 2011, 7: 19.
Copyright
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.