Okutetal.GeneticsSelectionEvolution2013,45:34 Genetics http://www.gsejournal.org/content/45/1/34 Selection Evolution RESEARCH Open Access Predicting expected progeny difference for marbling score in Angus cattle using artificial neural networks and Bayesian regression models Hayrettin Okut1,2*, Xiao-Liao Wu1,3, Guilherme JM Rosa1,4, Stewart Bauck5, Brent W Woodward6, Robert D Schnabel7, Jeremy F Taylor7 and Daniel Gianola1,3,4 Abstract Background: Artificial neural networks (ANN) mimic the function of the human brain and are capable of performing massively parallel computations for data processing and knowledge representation.ANN can capture nonlinear relationships between predictors and responses and can adaptively learn complex functional forms, in particular, for situations where conventional regression models are ineffective. In a previous study, ANN with Bayesian regularization outperformed a benchmark linear model when predicting milk yield in dairy cattle or grain yield ofwheat. Althoughbreeding values rely on theassumption of additive inheritance, the predictive capabilities ofANN are ofinterest from the perspective of their potential to increase theaccuracyof prediction ofmolecular breeding values used for genomic selection.This motivated thepresentstudy, in which the aim was to investigate theaccuracy ofANN when predicting theexpectedprogeny difference (EPD) of marbling score inAngus cattle. Various ANN architectures were explored, which involved twotraining algorithms,two types of activation functions, and from 1to 4 neurons inhidden layers. For comparison, BayesCπ models were used to select a subset of optimal markers (referred to as feature selection), under theassumption of additive inheritance, and then the marker effects were estimatedusing BayesCπwith π set equalto zero. This procedure is referred to as BayesCpC and was implemented on a high-throughput computing cluster. Results: The ANN with Bayesian regularization method performed equally well for prediction ofEPD as BayesCpC, based onprediction accuracyand sum of squared errors.With the 3K-SNP panel, for example,prediction accuracy was 0.776 using BayesCpC,and ranged from 0.776to 0.807using BRANN. With the selected 700-SNP panel, prediction accuracy was 0.863for BayesCpC and ranged from 0.842to 0.858 for BRANN. However,prediction accuracy for the ANN with scaled conjugate gradient back-propagation was lower, ranging from 0.653 to 0.689with the3K-SNP panel, and from 0.743to 0.793 with theselected 700-SNPpanel. Conclusions: ANN with Bayesian regularization performed as well as linear Bayesian regression modelsin predicting additive genetic values, supporting theidea thatANN are useful as universal approximatorsof functions ofinterest inbreeding contexts. *Correspondence:[email protected] 1DepartmentofAnimalSciences,UniversityofWisconsin,Madison,WI53706, USA 2DepartmentofAnimalScience,BiometryandGeneticsBranch,Universityof YuzuncuYil,Van65080,Turkey Fulllistofauthorinformationisavailableattheendofthearticle ©2013Okutetal.;licenseeBioMedCentralLtd.ThisisanOpenAccessarticledistributedunderthetermsoftheCreative CommonsAttributionLicense(http://creativecommons.org/licenses/by/2.0),whichpermitsunrestricteduse,distribution,and reproductioninanymedium,providedtheoriginalworkisproperlycited. Okutetal.GeneticsSelectionEvolution2013,45:34 Page2of13 http://www.gsejournal.org/content/45/1/34 Background complexfunctionalforms,throughaseriesoftransforma- The availability of genome-wide dense marker panels for tions (i.e., activation functions) driven by parameters. many species of plants and animals has opened doors Multilayer feed-forward is the most common architecture for incorporating genomic information into practical used in ANN, which consists of neurons assembled into breeding programs, an approach known as genomic se- layers. The first layer is called the input layer (the left- lection[1].Itisnoweasytogenerategenome-widescans mostlayeroftheANN)thatacceptsdata(e.g.,SNPgeno- with more than one million SNPs (single nucleotide types) from sources external to the network, and the last polymorphisms) but these huge databases pose chal- layer (the right-most layer) is called the output layer that lengesincomputationalcapacity, dataanalysisandinter- contains output units of the network. In between these pretation of results for genomic selection [2]. For twolayersareso-calledhiddenlayers,becausetheir values example, even with an initial screening to reduce the arenotobservedinthetrainingset.Hence,anANNarchi- number of markers to less than ten thousand SNPs, it is tectureisspecifiedbythenumberoflayers,thenumberof still not feasible for most computational platforms to neuronsineachlayer,andthetypeofactivationfunctions evaluate all combinations of SNPassociations,evenwhen that are used. Neurons in each layer are connected to the low-dimension interactions are explored [3]. Hence, re- neurons from the previous and the subsequent layer duction of dimensionality and feature extraction arguably throughadaptablesynapticweights. play pivotal roles in current genomic studies [4]. The in- The network weights are evaluated in two steps: the tensive computation inherent in these problems has al- feed-forward and the back-propagation steps. In the teredthecourseofmethodologicaldevelopments, and the feed-forward step, information comes from the input same is true for genomic selection [5]. layer and each unit evaluates its activation function, In the genome-enabled prediction of genetic merit for which is transmitted to the units connected to the out- breeding purposes, parametric statistical methods tend put layer. The back-propagation step consists of running to make strong assumptions about functional forms and the whole network backwards to minimize the error the statistical distribution of marker effects. On the one function in the space of weights using the method of hand, ridge regression best linear unbiased prediction gradient descent. A set of weights that minimizes the assumes that all markers have an effect on the trait of error function is considered to be a solution of the interest, and that these effects share a common variance learning problemfor theANN. intheirdistribution.Thissimpleassumptionisobviously Determination of the number of neurons to be placed not true in reality. On the other hand, hierarchical linear in the hidden layer is a critical task in the design of Bayesian regression models, such as BayesA and BayesB ANN. A network with an insufficient number of neurons [1], allow marker effects to be estimated with differential maynotbecapableofcapturingcomplexpatterns.Incon- shrinkage. However, posterior inference, in particular for trast, a network with too many neurons will suffer from the variances of marker effects, depends heavily on the over-parameterization, leading to over-fitting and poor prior assumptions used in these models [6]. Hence, predictive ability of yet to be observed data. Two popular BayesCπ [7] was proposed to overcome some of the techniques to overcome the over-fitting problem in ANN above mentioned drawbacks. A BayesCπ model postu- models are Bayesian regularization and cross-validated lates an unknown probability π that a SNP has no effect early stopping (CVES). Bayesian regularization (BR) con- at all on the trait, while all non-zero SNP effects are as- strains (shrinks) the magnitude of the networks weights, sumedto berandomsamplesfrom anormaldistribution improvesgeneralizationandallowsbiasinparameteresti- with null mean and a common variance. Recently, there mates towards values that are considered to be plausible, has been interest in the use of non-parametric methods while reducing their variance; thus, there is a bias- for the prediction of quantitative traits, such as reprodu- variance trade-off [13]. Unlike other back-propagation cing kernel Hilbert space regressions [8,9], radial basis trainingalgorithmmethodsthatuseasinglesetofparam- functionmodels[10]andnon-parametricBayesianmodels eters,BRconsidersallpossiblevaluesofnetworkparame- with Dirichlet process priors [11]. These nonparametric tersweightedbytheprobabilityofeachsetofparameters. methodsmakeweakerassumptionsandcanbemoreflex- In practice, Bayesian regularized ANN can be computa- iblefordescribingcomplexrelationships[12]. tionally more robust and cost-effective than standard Artificial neural networks (ANN), also known as back-propagationnetsbecausetheycanreduceorelimin- neuro-computationalmodels,provide an appealing alter- atetheneedforlengthycross-validation. native for genome-enabled prediction of quantitative With cross-validated early stopping (CVES) methods, traits [13,14]. These machine learning methods can act the training data set is split into a new training and a asuniversalapproximatorsofcomplexfunctions[15]be- validationsetandagradientdescentalgorithmisapplied causetheyarecapableofcapturingnonlinearrelationships to the new training data set. The ANN performs an it- betweenpredictorsandresponsesandcanadaptivelylearn erative process, first using the training data set to Okutetal.GeneticsSelectionEvolution2013,45:34 Page3of13 http://www.gsejournal.org/content/45/1/34 calculate the weight and bias estimates, and then applies BeadChip [16] was extracted. After data quality control these parameter estimates in the validation data set to and screening, a total of 2421 polymorphic SNPs were calculate the prediction errors. The process iterates re- retained for analysis. The target variable to be predicted peatedly, substituting parameter estimates from the was EPD for marbling score, which had been computed training data set into the validation data set to find the by the American Angus Association using BLUP (best smallest possible average prediction errors for the valid- linear unbiased prediction) for each of these animals, ation data set. Training ceases when the error in the val- based upon their pedigree data and progeny carcass and idation data set increases in certain consecutive epochs ultrasound data [17]. In animal breeding, an EPD is de- (iterations) in order to avoid the problem of over-fitting fined as the predicted performance of a future offspring (the number of consecutive epochs is 6 by default in of an animal for a particular trait (marbling score), cal- MATLAB). The ANN parameter estimates with the best culated from measurement(s) of the animal’s own per- performanceinthe validationsetisthenusedonthetest- formance and/or the performance of its relatives under ingdatatoevaluatethepredictiveabilityofthenetwork. the assumption of additive inheritance. Hence, EPD rep- In a previous study, ANN with BR outperformed a resent a typical linear system since the EPD of an indi- benchmark linear model when predicting milk yield in vidual can be represented as a linear function of the dairy cattle or grain yield of wheat [14]. However, be- EPD of relatives and a residual term that reflects the fact cause breeding values are defined in terms of linear that an individual inherits a random sample of the genes functions based upon additive inheritance, the predictive present within its parents. The distribution of marbling performance of ANN relative to linear additive systems score EPD in the Angus sample was symmetric and is of some interest. This motivated the present study, in suggested a normal distribution (Figure 1), with mean which the aim was to investigate the accuracy of ANN and standard deviation estimated at 0.0265 and 0.254, for predicting expected progeny differences (EPD) for respectively [18]. marbling score in Angus cattle relative to hierarchical linear Bayesian regression models. Various ANN archi- Statisticalmethods tectures were explored, which involved two training Artificialneuralnetworks algorithms, two types of activation functions, and from 1 A schematic presentation of a multilayer perceptron to4neuronsinhiddenlayers. (MLP) feed-forward of an ANN is presented in Figure 2. This is a multi-layer feed-forward neural network with a Methods single middle (hidden) layer and four neurons. Training Datasets an ANN involves the minimization of an error function The data contained 3079 registered Angus bulls, geno- that depends on the network’s synaptic weights (the w’s typed with the Illumina BovineSNP50 BeadChip, from in Figure 2), and these weights are iteratively updated by which the SNP content for the Illumina Bovine3K a learning algorithm, to approximate the target variable. Mean 0.0265 Skewness 0.468 StdDev 0.254 Kurtosis 0.356 Range 1.920 Interquartile 0.330 Figure1Histogramanddensityplotofderegressedexpectedprogenydifferencesformarblingscorefor3079Anguscattle. Okutetal.GeneticsSelectionEvolution2013,45:34 Page4of13 http://www.gsejournal.org/content/45/1/34 b(1) 1 SNP1 X1 w1,1 ∑ f1 w2,1 b2(1) w1 SNP2 X2 . ∑ f2 SNP3 X3 . b3(1) w2 b(2) . . . ∑ f3 w3 ∑ g y . b(1) 4 SNP2421 X2421 w4,2421 w4 ∑ f4 Input layer Hidden layer, 4 neurons hyperbolic tangent Output layer, 1 neuron linear activation Figure2Schematicrepresentationofathree-layerfeed-forwardneuralnetwork.Genotypesfor2421(or700)SNPswereusedasinputs x={x|i=1,2,…,n},wherenisthenumberofindividualswithgenotypes;eachSNPwasconnectedtoupto4neuronsviacoefficientsw , j ij kj wherekdenotesneuronandjdenotesSNP;here,w isaweightfromahiddenlayerunitstotheoutputunit,f isanactivationfunctionapplied k k tohiddenlayerunits(e.g.,thehyperbolictangent),gisanactivationfunctionappliedtotheoutputlayerunit(e.g.,linear),b(1)andb2arebiases ofhiddenandoutputlayerunits,and^y isapredictedvalue. Theupdatingisusuallyaccomplishedbyback-propagating Xp bð1Þþ w x, where p=700 or 2421, and bð1Þ is the the error, which is essentially a non-linear least-squares k kj j k problem [19]. Back-propagation is a supervised learning j¼1 bias parameter defined in the hidden layer and this algorithmbasedonasuitableerrorfunction,the valuesof quantity is then transformed using some linear or which are determined by the target (i.e., marbling EPD) nonlinear activation function (f ) as: and the mapped (predicted) outputs of the network (i.e., k ! fitted values ofmarbling EPD).Weightsin anMLP archi- Xp tecture are determined by a back-propagation algorithm f bð1Þþ w x : ð1Þ k k kj j to minimize the sum of squares of errors using gradient- j¼1 descent methods [13]. During the training, weights and The above is the output of the hidden layer, which in biases in the ANN are successively adjusted based on the turn is delivered to the output layer (e.g., each neuron k input and the desired output. Each iteration of feed- in the hidden layer sums ANN input x after multiplying forward in a MLP constitutes two sweeps: forward acti- j them by the strengths of the respective connection vation to produce a desired output, and backward weights, w , and adds biases b and then computes its propagation of the computed error to update the values kj k outputasafunction ofthe sum),andiscollected as: for the weights and biases. The forward and backward ! sweeps are repeatedly performed until the ANN solution XK Xp agrees with the desired value to within a pre-specified bð2Þþ wkfk bð1Þþ wijxj ; ð2Þ tolerance [20,21]. k¼1 j¼1 (cid:1) (cid:3) Let x0j ¼ x1j;x2j::::::;xnj be a row vector that con- where wk is the weight specific to the kth neuron (k=1, tains SNP genotypes for all i = 1,…,n individuals, for 2,…,K), and b(2) is the bias parameter defined in the the jth SNP, where j = 1,…..,p, and p = 700 (referred output layer. Next, quantity (2) is activated with the fol- to as the 700-SNP panel) or p = 2421 (referred to as lowing function: " # the 3K-SNP panel). In an ANN, SNP genotypes are XK connected to each neuron in a single hidden layer via gð:Þ¼g w f ð:Þþbð2Þ : ð3Þ k k weights wkj, for k=1,…,K neurons, with each weight k¼1 defining a specific SNP-neuron connection with ap- The above yields the fitted marbling EPD, ^y , for each propriate biases (intercepts), bð1Þ;bð1Þ; ::::::;bð1Þ, each individual in the training set. Finally, the prediicted value 1 2 k pertaining to a specific neuron. The input into ofmarblingscoreEPD canbecomputedinthetestingset neuron k, prior to activation is expressed linearly as using a formula similar to (3). We used the hyperbolic Okutetal.GeneticsSelectionEvolution2013,45:34 Page5of13 http://www.gsejournal.org/content/45/1/34 tangent sigmoid and linear (identity) activation functions minimize the following function, equal to a penalized inthehiddenandoutputlayers,respectively. log-likelihood, We used BR and scaled conjugate gradient (SCG) F ¼βE ðDjw;MÞþαE ðwjMÞ; ð5Þ back-propagation as training algorithms. The basic idea D W of back-propagation algorithms is to adjust weights in the steepest descent direction (negative of the gradient), where E (w|M) is the sum of squares of network W such that the objective function decreases most rapidly weights, M is the ANN architecture, and α and β are [22]. However, in practice, this does not necessarily pro- positive regularization parameters that must be esti- duce the fastest convergence, although the function mated. The second term on the right hand side of (5), may decrease most rapidly along the negative of the known as weight decay, favors small values of w and de- gradient. Several conjugate gradient algorithms have creases the tendency of a model to over-fit the data. been proposed, in which a search is performed along Hence, training involves a tradeoff between model com- conjugate directions, generally leading to faster conver- plexityandgoodnessoffit.Largevaluesofαleadtopos- gence to a local function minimum than steepest des- terior densities of weights that are highly concentrated cent directions [22,23]. On the one hand, with the around zero, so that the weights effectively disappear, exception of SCG, the conjugate gradient algorithms the model discounts connections in the network [24,25], use linear searches at each iteration and thus, are com- and complex models are automatically self-penalized. putationally expensive as they require that the network From equation (5), if α << β, the fitting or training algo- response to all training inputs be computed several rithm places the most weight on goodness of fit. If α >> times for each search. On the other hand, the SCG al- β, emphasis is placed on reducing the magnitude of the gorithm combines a model-trust region approach with weights at the expense of goodness of fit, while produ- the conjugate gradient approach and uses a step size cinga smoothernetwork response[26]. scaling mechanism to avoid time-consuming linear In an empirical Bayesian framework, the “optimal” searches [23]. Hence, SCG can significantly reduce the weights are those that maximize the conditional poster- number of computations performed in each iteration, ior density P(w|D,α,β,M), which is equivalent to min- but may require more iterations to converge than do the imizing the regularized objective function F in equation other conjugate gradient algorithms [22]. SCG uses the (5). Minimization of F is identical to finding the (locally) CVES technique to prevent over-fitting, which is com- maximum a posteriori estimates of w, denoted wMP, monly used in neural networks and can improve their which minimize E using the back-propagation training D predictive ability (generalization). algorithms [24]. However, this is possible only if n > m, Training proceeds as follows. Let the data set be D= where misthenumber ofparameterstobeestimated. {y,{xi}i=1,…,n}, where xi is a vector of inputs (SNP geno- Bayestheoremyieldstheposteriordensityofαandβas: types) forindividualiandyisavectoroftarget variables (EPD). Once a set of weights, w, is assigned to the con- Pðα;βjD;MÞ¼PðDjα;β;MÞPðα;βjMÞ: nections in the networks, this defines a mapping from PðDjMÞ the input x to the output ^y . Let M denote a specific i i If the prior density P(α,β|M) is uniform, maximization network architecture (in terms of numbers of neurons of P(α,β|D,M) with respect to α and β is equivalent to and choice of activation functions), then the typical ob- maximization of P(D|α,β,M). Bayesian optimization of jective function used for training a neural network using the regularization parameters requires computation of CVESisthesum ofsquaredprediction errors(E ): D the Hessian matrix of the objective function F evaluated Xn at the optimum point wMP [27], but directly computing E ðDjw;MÞ¼ ð^y−yÞ2 ð4Þ D i i the Hessian matrix is not always necessary. As proposed i¼1 by MacKay [28], the Gauss-Newton approximation to the Hessian matrix can be used if the Levenberg-Marquardt forninput-target pairsdefining D. optimization algorithm is employed to locate the mini- Regularization produces shrinkage of parameter esti- mum of F [13,29,30]. The Levenberg–Marquardt training mates towards some plausible values, while at the same algorithm [31] achieves second-order trainingspeed with- time reducing their variance. In Bayesian models, re- out computing the Hessian matrix. Briefly, the Hessian gularization techniques involve imposing certain prior matrixisapproximatedas: distributions on the model parameters. In Bayesian ANN (e.g., BRANN), the objective function that was H¼J0J; ð7Þ specified in (4) has an additional term that penalizes large weights, in the hope of achieving a smoother map- where J is the Jacobian matrix that contains first deriva- ping. Gradient-based optimization is then used to tives of the network errors with respect to network Okutetal.GeneticsSelectionEvolution2013,45:34 Page6of13 http://www.gsejournal.org/content/45/1/34 parameters(theweightsand biases).Thegradientiscom- unknown and needs to be inferred, with the prior distri- putedas: bution ofπtaken tobeuniform between0and1, g¼J0e; ð8Þ πeUniformð0;1Þ: ð14Þ andnetworkparametersareupdatedas: A Bernoulli indicator variable, δ, is introduced to fa- j cilitate sampling from the mixtures of the SNP effects, wkþ1 ¼wk−ðJ0JþμIÞ(cid:2)J0e: ð9Þ asfollows: Here, μ is Levenberg's damping factor, which is ad- pðδijπÞ¼π1−δið1−πÞδi: justed at each iteration and guides the optimization process. If the reduction of error sum of squares is Hence, unconditionally,(cid:1) the va(cid:3)riable αj follows a rapid, a smaller value of μ is used to bring the algo- univariate-t distribution, t 0;S2α;υα , if δj=1, or is equal rithm closer to the Gauss–Newton algorithm. Alterna- to zero otherwise [6]. Posterior inference of unknown parameters in the Bayesian model via Markov chain tively, the damping factor is increased to give a step Monte Carlo (MCMC) implementation is described in to gradient descent direction if an iteration provides [7]. With a subset of, say k≤p, selected markers, the insufficient reduction in the error sum of squares [32]. Optimal values of regularization parameters, α and β, statistical model takes the same form as (12), replacing p with kforthenumberofmarkers. in BRANN can be calculated as: By assuming that all k of the selected SNPs (based on γ αMP ¼ ð10Þ the posterior model probability and including the fre- 2EwðwMPÞ quency of marker k during MCMC) have non-null effects on the quantitative trait, we define a BayesCπ and model with π=0, which was used for the statistical in- n−γ ference and model cross-validation subsequent to selec- βMP ¼ ; ð11Þ 2E ðwMPÞ tion of markers (referred to as post-selection hereafter). D So, posterior inference in BayesCπ with π = 0 was as for where 0 ≤ γ=m−2αMPtr(HMP)−1≤m the number of BayesC π, except thatπ was fixedatzeroandhence sam- effective parameters in the neural network, and m is plingoftheindicatorvectorδwasnolongerrelevant. the total number of parameters. Here HMP is the max- imum a posteriori estimate of H described in (7). Computationalimplementation MATLAB [31] was used tofit the ANN. Each neural net- BayesCπandBayesCwithπ=0 workconsistedofthreelayers(i.e.,input,hiddenandout- Considerthefollowinglinearmodel: put layers). The number of neurons in a single hidden layer variedfrom1to4.EachANNhad2421(or700)in- Xp y ¼μþ x α þe; ð12Þ puts (SNPs). Before processing, MATLAB automatically i ij j i rescaled all input and output variables using the j¼1 “mapminmax”functionsuchthattheyresidedintherange where y is the marbling score EPD for the ith animal; μ [−1, +1], to enhance numerical stability. Two combina- i is the overall mean; α is the substitution effect associ- tions of activation functions were used: (i) a set of hyper- j ated with the jth SNP (j=1,…,p); x is an indicator vari- bolic tangent sigmoidal activation functions from the ij able corresponding to the genotype(cid:1)at the(cid:3) jth SNP (0, 1, input layer to the hidden layer, plus a linear activation 2) for the ith individual, and eeN 0;σ2 is a residual functionfromthehiddenlayertotheoutputlayer,and(ii) i e term,whereσ2 istheresidualvariance. a set of linearactivation functions from the input layerto e The BayesCπ model [7] assumes that each SNP effect the hidden layer and from the hidden layer to the output is null with p(cid:1)roba(cid:3)bility π, or that it follows a normal dis- layer.Trainingwas stopped ifanyofthefollowingcriteria tribution,N 0;σ2 ,with probability 1-π,i.e.: were met: (i) a maximum number (1000) of epochs was α reached, (ii) performance had met a pre-specified (the (cid:5) (cid:4) (cid:1) (cid:3) performance function for feed-forward networks is the αi(cid:5)(cid:5)(cid:5)π;σ2α ¼eN00;σ2α wwiitthhpprroobbaabbiilliittyyð1−ππÞ: ð13Þ msueitaanblesqtuaragreet,eorrro(riv))ltehveelL,e(viiei)nbthereg-gMraadriqenutarwdtapsabraemlowetear exceeded1010. Here, σ2α is a variance common to all non-zero SNP For each ANN architecture, eight replicates were run. effects, which (cid:1)is assig(cid:3)ned a scaled inverse chi-square dis- Each replicate was independently initialized, in order to tribution, χ−2 vα;s2α . Furthermore, the value of π is eliminate spurious effects caused by the starting values, Okutetal.GeneticsSelectionEvolution2013,45:34 Page7of13 http://www.gsejournal.org/content/45/1/34 and to improve predictive ability. The results were pre- results than prediction using a smaller panel, yet the op- sented as averages across the eight replicates per ANN timal panel size may depend on many factors. In this architecture. study, we empirically chose the 700-SNP subset as the BayesCπ with π set equal to zero (referred to as the optimal panel. The fitting accuracy in the training set BayesCpC procedure) was implemented via a high- and predictive accuracy in the testing set using the opti- throughput computing pipeline to select SNPs for post- mal700-SNPsubset areillustrated inFigure3b. selection statistical inference and cross-validation [5]. This pipeline ran multiple chains for both feature selec- DeterminationofanoptimalANNarchitecture tion and cross-validation. A three-fold cross-validation The performance of the ANN architectures was exa- approachwasemployed,inwhichthewholedatasetwas mined based on the sum of squared errors (SSE) of divided into three approximately equal portions, with prediction with a 3K panel, averaged over eight inde- two-thirdsusedfortrainingandone-thirdusedfortest- pendent replicates of each ANN, for both BRANN and ing, and the portions used for training and testing were SCGANN. Each ANN had a distinct combination of rotated three times. Each cross-validation experiment training algorithm, transformation method, and number was randomly replicated eight times. Three parallel of neurons, but both, BRANN and SCGANN, had an MCMC chains were run for each feature-selection ana- input of 3K SNPs. The average SSE ranged from 13.936 lysis, and each consisted of 50 000 iterations after a to 16.882 for BRANN, and from 36.531 to 39.571 for burn-in of 5000 iterations, thinned every tenth iterate. SCGANN. Smaller SSE were produced by BRANN with MCMC sampling for each cross-validation consisted of nonlinear activation functions and from 2 to 4 neurons, 100 000 iterations, with a burn-in of 10 000 iterations, and by SCGANN with nonlinear activation functions thinnedeverytenthiterate. and from 1 to 4 neurons. There was no evidence that more complex networks (e.g., with more neurons or a Results non-linear transformation function) produced better DeterminationofanoptimalSNPpanelsize predictions than the linear model, as the ANN were The predictive performance of each ANN was examined similar in terms of their SSE. Possibly, this was because usingeitherthe3K-SNPpaneloranoptimalsubsetof700 marbling score EPD is estimated under an additive selectedSNPs.Thelatterwerederivedfromthe3K-panel, model of inheritance in which additive genetic merit has selected using the BayesCpC procedure with three-fold a linear relationship with SNP effects. Nevertheless, cross-validation.This was accomplishedbyexaminingthe BRANN performed as well as the linear models when predictionperformanceof varyingpanelsizesfrom100to predicting this linear system. Also, BRANN consistently 2400 SNPs in 100-SNP increments, and the optimal set produced a more accurate prediction of marbling score thatgavethebestpredictionincross-validationswasiden- EPD than did SCGANN. On average, SSE obtained from tified. The reason for not choosing the optimal subset BRANN were about 40% to 50% of those obtained from based on ANN models was because the selection tasks SCGANN (Figure 4). This is attributed to the use of with a grid of 24 candidate SNP-panels of varying sizes Bayesian regularization, which handles over-fitting better were too computationally intensive for BRANN. Never- than doestheSCG back-propagation. theless, the parallel-BayesCpC pipeline handled this task easily, because all jobs were submitted to run in parallel Predictiveperformanceusingthe3K-SNPpanel onaclusterofdedicatedcomputers. BRANN and BayesCpC performed very similarly with As shown in Figure 3a, the correlation between marb- the 3K-SNP panel and both methods yielded higher pre- ling score EPD and their fitted values in the training set diction accuracies than did SCGANN (Figure 5). On (referred to hereafter as fitting accuracy) increased al- average, the correlation in the training set was 0.941 most monotonically with panel size, until reaching a with BayesCpC. This correlation ranged from 0.942 to plateau at a panel size of 1400 SNPs. However, the 0.967 with BRANN, and from 0.796 to 0.897 with correlation between marbling score EPD and their pre- SCGANN. The average correlation (in the testing set) dicted values in the testing set (referred to hereafter as was 0.776 with BayesCpC, and ranged from 0.776 to predictive accuracy) reached its peak (0.863) with a 0.807 with BRANN, and from 0.653 to 0.689 with panel size of 700 SNPs, and decreased thereafter. The SCGANN. In general, these correlations increased decrease in prediction accuracy with > 700 selected slightly with the number of neurons, but no consistent SNPs possibly reflects over-fitting in the training set, pattern wasobserved(Figure 5). which, inthis case, happened much before the panel size Withthe3K-SNPpanel,thenumberofSNPs(i.e.,2421) exceeded the training set size (i.e., approximately 2000 exceeded the number of animals (~2000 animals) in the animals). Hence, with Bayesian regression models, pre- training set. This means that there were more parameters dictionusingmoreSNPsmaynotnecessarilyyieldbetter to be estimated than data points, even when all SNPs Okutetal.GeneticsSelectionEvolution2013,45:34 Page8of13 http://www.gsejournal.org/content/45/1/34 a Training Testing 0.96 0.94 0.92 0.9 0.88 s n 0.86 o ati 0.84 rel 0.82 r o 0.8 C 0.78 0.76 0.74 0.72 0.7 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 Number of Selected Markers X 100 b Figure3SelectionofanoptimalSNPpanelsizeforpredictingmarblingscoreexpectedprogenydifference(EPD)usingtheBayesCpC approach.(a)CorrelationsbetweenmarblingscoreEPDandtheirfitted(predicted)valuesinthetrainingandtestingsets,wheretheSNPpanel consistedof100,200,…,2400markers,respectively;(b)correlationbetweenEPD(inputvalues)andfittedvaluesinthetrainingset(upper)and correlationbetweenEPD(inputvalues)andpredictedvaluesinthetestingsetfortheoptimalsetof700SNPs(lower). entered the model linearly. This explains the better pre- the posterior mean (standard deviation) of π was 0.621 dictive performance of BRANN over SCGANN: Bayesian (0.0830), based on 2421 polymorphic SNPs; this means regularization imposes a penalty for model complexity that, on average, 918 SNPs were fitted in the model at and prevents over-fitting, while SCGANN captures both each MCMC iteration. Hence, over-fitting was not a signalandnoiseinthetrainingset,leadingtopoorpredic- concern. Posterior densities for π are shown in Figure 6 tion. The results illustrate that ANN with Bayesian and were computed using posterior samples obtained regularization can behave better than ANN with SCG from each of the three parallel chains and from all back-propagation. BRANN penalize the estimation of pooled samples of the three chains. The density plots for largeweightstoachieveasmoothermapping,whichleads the three parallel chains highly resembled each other toeffectivelymoreparsimoniousmodels. and had similar means, indicating that convergence and In the BayesCpC procedure, the BayesCπ model pos- length of the MCMC chains were appropriate for π. We tulates that a portion, π, of all SNPs have no effect on noticed that the estimate of π obtained from the 3K marblingEPD.Inahigh-densitySNPpanel,πistypically panel wassmallerthanthatbased on BovineSNP50geno- expected to be large, meaning that the portion of “sig- types, because the latter has typically been greater than nal” SNPs, 1 - π, is small, and the chance of over-fitting 0.95[5].Thus,itwouldseemthat,withalow-densitySNP is reduced. Using the Illumina Bovine3K SNP genotypes, panel, the interpretation of SNPs having putatively non- Okutetal.GeneticsSelectionEvolution2013,45:34 Page9of13 http://www.gsejournal.org/content/45/1/34 Figure4Sumofsquarederrors(y-axis)inthetestingsets,computedasaveragesfromeightreplicatesforBayesianneuralnetworks (BRANN)andANNcoupledwiththeSCGalgorithm(SCGANN),forANNwithdifferentnumbersofneurons. zero effects on the target traits should be taken with cau- thebestpredictionswhengeneralizedbeyondthetraining tion, because many could be distant from the functional set.Thisphenomenonisreferredtoaspoorgeneralization genes or quantitative trait loci. We also observed that the inmachinelearningmethods[12]. selected number of SNPs having non-zero effects based on BayesCπ in the training set did not correspond to the Predictiveperformanceusingthe700-SNPpanel number of SNPs in the optimal SNP panel for prediction Unlike the case with 3K-SNPs, prediction using the (918 vs 700 SNPs). We suspect that parameter π in a selected 700-SNPs was not challenged by over-fitting. BayesCπmodeldoesnotfullyinformonthesizeofanop- Hence, we observed a smaller difference in prediction timal SNP panel for prediction; a model that describes performance between BRANN and SCGANN because variation in a training set well does not necessarily yield regularization was not a decisive issue in the prediction. 1.2 Training Testing 1 0.8 0.6 0.4 0.2 0 ayes CpC Linear 1_Neuron _Neurons _Neurons _Neurons Linear 1_Neuron _Neurons _Neurons _Neurons B 2 3 4 2 3 4 BRANN SCGAN Figure5Correlationsbetweenmarblingscoreexpectedprogenydifferencesinthetraining(testing)setsandtheirfitted(predicted) valuesusingBayesCpCandBRANNandSCGANNwithdifferentnumbersofneuronsinthehiddenlayerandusingthe3K-SNPpanel. 1training=correlationsinthetrainingsets;testing=correlationsinthetestingsets;2BayesCpC=Bayesianregressionmodel,wheretheBayesCπ modelisusedforfeatureselectionandBayesCπwithπ=0isusedforpost-selectionstatisticalinferenceandcross-validation;BRANN=artificial neuralnetworkwithBayesianregularization;SCGANN=artificialneuralnetworkwithscaledconjugategradientback-propagation. Okutetal.GeneticsSelectionEvolution2013,45:34 Page10of13 http://www.gsejournal.org/content/45/1/34 5 5 4 4 sity 3 sity3 n n De 2 De2 1 1 0 0 0.3 0.4 0.5 0.6 0.7 0.8 0.9 0.3 0.4 0.5 0.6 0.7 0.8 0.9 Chain #1 Chain #2 5 5 4 4 sity 3 sity3 n n De 2 De2 1 1 0 0 0.3 0.4 0.5 0.6 0.7 0.8 0.9 0.3 0.4 0.5 0.6 0.7 0.8 0.9 Chain #3 pooled Figure6PosteriordensityforparameterπintheBayesCπmodelobtainedfromthreeparallelchainsandfrompooledsamplesofthe threechains. With the 700-SNP panel, the predictive performance of inheritance, and this smoothing may mask non-linear BayesCpC was slightly better than that of BRANN, variation in the response variable. A better way to possibly because the subset of SNPswasoptimal, at least analyze this trait would be to remove variation due to as selected by BayesCpC. As expected, BayesCπ and contemporary groups from the field data and then BRANN outperformed SCGANN in the prediction of analyze individual marblingphenotypes, but this was not marbling score EPD but differencesin performance were possible here because we did not have access to the raw smaller with the 700-SNP panel than with the 3K-SNP data, which were, in general, collected on progeny of the panel. Average prediction accuracy was 0.863 with genotyped bulls. Theaccuracy ofEPD varies between in- BayesCpC and ranged from 0.843 to 0.858with BRANN, dividuals, which suggests that the residuals may be whereas the average prediction accuracy with SCGANN heteroscedastic due to unequal prediction error vari- ranged from 0.743 to 0.793. For both BRANN and ances. An alternative is to use deregressed EPD as the SCGANN, the difference in predictive accuracy between target variable, for which the parent averages are re- the linear and non-linear activation functions was negli- moved and appropriate weights can be applied to ac- gible (Figure 7). count for heteroscedastic residuals [33]. Some reports suggest that training on deregressed EBVcould generate Discussion higher accuracy of genome-enhanced breeding values We investigated the predictive performance of ANN in than training on EBV [34]. However, the main purpose comparison to a hierarchical Bayesian regression model, of the present research was to investigate the predictive using 3K-SNP and 700-SNP panels. The 700 SNPs were performance of ANN in comparison with Bayesian re- preselected from the former based upon their power to gressionmodelsinalinearsystem.WeusedEPDinstead predict marbling EPD. Various ANN architectures to of deregressed proofs because correct deregression in- predict marbling score EPD were examined, in conjunc- volves matrix operations (the data available were incom- tionwithtwotrainingalgorithms,twotypesofactivation plete for correct deregression), or approximations. Both functions, and from 1 to 4 neurons in the hidden layer. EPD and deregressed EPD are heteroscedastic because We did not observe significant differences in predictive of an unequal amount of information, but the neural accuracy between linear and non-linear models, prob- network software useddidnot allowfor incorporation of ably because the relationship between marbling score heteroscedasticityina straightforwardmanner. EPD and SNP effects is theoretically linear. An EPD In genomic selection, the joint estimation of a genome- produces a smoothed data point based on additive wide SNP set is challenging due to the so-called “large p,