Model 1: threshold Bayes A

where, λ is the underlying liability variable vector for y, μ is the population mean, 1 is a column vector (n×1) of ones; b = {b _j } corresponds to the vector for the regression coefficient estimates of the p markers or SNP assumed normally and independently distributed a priori as $N(0,{\sigma }_j^2)$, where ${\sigma }_j^2$ is an unknown variance associated with marker j. The prior distribution of ${\sigma }_j^2$ is assumed to be distributed as the scaled inverse chi-square ${\sigma }_j^2~{\upsilon }_j{s}_j^2{\chi }_{{\upsilon }_j}^{-1}$, with υ _j = 4 and $s_j^2=0.002$. Elements of the incidence matrix X, of order n × p, may be set up as for different additive, dominant or epistatic models. In the more practical scenario, it takes values -1, 0 or 1 for marker genotypes aa, Aa and AA, respectively. The residuals (e) are assumed to be distributed as $N(\mathbf{0},{\sigma }_e^2)$, with residual variance ${\sigma }_e^2=1$, as stated above. As in a regular threshold model, two parameters have to be set fixed (e.g. threshold and the residual variance are set to zero and one, respectively) since these parameters are not identifiable in a liability model.

Model 2: threshold Bayesian LASSO

where λ is the vector of liabilities for all individuals, μ is the population mean, 1 is a column vector (n × 1) of ones; $\stackrel{\mathbf{\wedge }}{\mathbf{\beta }}$ are the LASSO estimates with their respective incidence matrix X as described for model TBA. As a modeling choice, e was considered the vector of independently and identically distributed residuals, as $\mathbf{e}\sim N(0,{\sigma }_e^2)$. In accordance with tradition, we fixed the threshold to be 0 and the residual variance to be 1 as described for model TBA; alternate choices result in the same model.

where ${\sigma }_e^2$ is the residual variance, and γ is a parameter controlling the shrinkage of the distribution. Inferences about γ may be done in different ways [16]. To follow the Bayesian specifications, a gamma prior is proposed here for γ ² , with known rate (r) and shape (δ) hyper-parameters, as described by de los Campos et al. [17]. Samples from posterior distributions of those estimates may be drawn from the Gibbs sampling algorithm described in de los Campos et al. [17], with the corresponding data augmentation algorithm for liabilities, as described for TBA.

Model 3: gradient boosting

Each predictor ( h _m (y; X) for m ∈ (1, M)) is applied consecutively to the residual from the committee formed by the previous ones. This algorithm can be calculated using importance sampling learning ensembles as follows:

(Initialization): Given data (y, X), let the prediction of phenotypes be F₀ = μ, with μ being the population mean.

Then, for m in {1 to M}, with M being large, calculate the loss function (L) for $\left(y_i,F_{m-1}({\mathbf{x}}_i)+h(y_i;{\mathbf{x}}_i,j_m)\right)$

where j _m is the SNP (only one SNP is selected at each iteration) that minimizes ${\displaystyle \sum _{i=1}^{n}L\left(y_i,F_{m-1}({\mathbf{x}}_i)+h(y_i;{\mathbf{x}}_i,j_m)\right)}$ at iteration m, h(y _i ; x _i , j _m ) is the prediction of the observation using SNP j at the current iteration, F _{m -1}(x _i ) is the updated prediction at the previous iteration and L(·) is a given loss function. The updated prediction at each iteration m may be expressed as F _m (x _i ) = F _{m -1}(x _i )+v·h(y _i ; x _i , j _m ) with v being some shrinkage factor that, without loss of generality, can be assumed constant and small (0< v <1), but it may be optimized to balance predictive ability and computation time.

Therefore, after the initialization, the algorithm flows as follows:

Step 1: Compute residuals as ${\mathbf{r}}_m=\mathbf{y}-{\displaystyle \sum _{i=0}^{m-1}v\cdot F_{m-1}({\mathbf{x}}_i)}$, and fit the weak learner for each SNP j (j ∈{1,..., p}) to current residuals, where ν was set to 0.01.

Step 2: Select SNP j, where $j=\arg {\min }_j{\displaystyle \sum _{i=1}^{n}L\left(y_i,F_{m-1}({\mathbf{x}}_i)+h(y_i;{\mathbf{x}}_i,j_m)\right)}$, i.e. the SNP minimizing the loss function.

Step 3. Update predictions as F _m (x _i ) = F _{m -1}(x _i )+ν·h(y _i ; x _i , j _m ), (i∈{1,..., n}), where h(y _i ; x _i , j _m ) is the estimate for individual i obtained by regressing the current residual (r _i ) at iteration m on its genotype for the SNP selected in step 2.

Step 4: Increase the iteration index m by 1, and repeat steps 2-4 until a convergence criterion is reached.

Here, we used ordinary least square regression as predictor h(y; X) and two different loss functions: the L ₂ loss function (L2B), which is a quadratic error term in the form (y _i -F _m (y _i ; x _i , j _m ))², and a pseudo-Huber loss function (LhB) in the form $\log \left[\cosh \left(y_i-F_m(y_i;{\mathbf{x}}_i,j_m)\right)\right]$. The pseudo Huber loss function is a priori more appealing for discrete traits because it is continuous, differentiable, greater than or equal to the logit loss function and overcomes the disadvantage of the squared loss by becoming more linear when (y _i -F _m (y _i ; x _i , j _m )) tends to infinite. The choice of the number of iterations, M, is a model comparison problem which may be overcome in many different ways [12, 20]. Here, a cross-validation design was used as described in González-Recio et al. [8]. More details on the gradient boosting can be found in Freund and Schaphire [21], Friedman [12] and González-Recio et al. [8].

Model 4: Random Forest

Random Forest can be viewed as a machine learning ensemble algorithm and was first proposed by Breiman [11]. It is massively non-parametric, robust to over-fitting and able to capture complex interaction structures in the data, which may alleviate the problems of analyzing genome-wide data. This algorithm constructs many decision trees on bootstrapped samples of the data set, averaging each estimate to make final predictions. This strategy, called bagging [22], reduces error prediction by a factor of the number of trees.

A RF algorithm aimed at genome-wide prediction is described next, in a more extensive manner than the previous methods, as this is the first time that this algorithm is used in a genomic breeding value prediction context:

Each tree (h _t (y; X) for t∈(1, T)) is distinct from any other in the ensemble as it is constructed from n samples from the original data set selected at random with replacement, and at each node only a small group of SNP are randomly selected to create the splitting rule. Each tree is grown to the largest extent possible until all the terminal nodes are maximally homogeneous. Then, c_t is some shrinkage factor averaging the trees. The trees are independent identically distributed random vectors, each of them casting a unit vote for the most popular outcome of the disease at a given combination of SNP genotypes.

Each tree minimizes the average loss function of the bootstrapped data, and is constructed using a heuristic approach as follows:

1. First, bootstrapped samples from the whole data set are drawn with replacement so that realization (y _i , x _i ) may appear several times or not at all in the bootstrapped set Ψ^(t) t = (1,..., T).

with L(y, h _t (X)) being a certain loss function. i.e. SNP j is the one that minimizes a given loss function at the current node, and is selected in this step. The algorithm takes a fresh look at the data that have arrived at each node and evaluate all possible splits. Many loss functions can be chosen (e.g. logit function, squared loss function, misclassification rate, entropy, Gini index, ...). The behavior of a given loss function may depend on the nature of the problem. The squared loss function is popular for continuous response variables, and the logit function for categorical responses.

3. Split the node in two child nodes according to SNP j genotype that one individual may or may not have (e.g. individuals with the risk allele will pass to a child node, and the remaining animals will pass to the other child node).

4. Repeat steps 2-3 until a minimum node size is reached (usually <5). The predicted value of the genotype x _i is the majority vote for the outcome at the terminal nodes (for regression problems, it is the average phenotype of the individuals in the node).

Finally, a large amount of trees are constructed repeating steps 1-4 to grow a random forest. The forest may be stopped when the generalization error averaged across the out of bag samples (see section below) have converged. Convergence may be visually tested but it may also be determined using traditional methods for convergence testing of Monte Carlo Markov chains.

Final predictions can be made by averaging the values predicted at each tree to obtain a probability of being susceptible. In a naïve 0 = non-susceptible/1 = susceptible scenario, individuals with probability <0.5 may be considered as non-susceptible. To predict observations of new individuals, their marker genotypes are passed down each tree, and the estimate of the corresponding terminal nodes is assigned to the new individual in each tree. The predictions of each tree in the RF algorithm are averaged for each animal to compute the final prediction.

There are two main aspects that can be tuned in random forest: the first one is the number of SNP or covariates sampled at random for each node (mtry). Generalized cross-validation strategies can be used to optimize mtry. In high dimensional problems such as GWAS, Goldstein et al. [23] have suggested mtry to be fixed to >0.1 p. The algorithm may speed up for smaller mtry values. Nonetheless, cross-validation can be used to determine the best value of mtry for each trait, although at an expense of increasing computation time. Genetic background may influence the behavior of this tuning parameter. The second aspect is the criterion to select the best SNP to split the node. As commented above, different criteria may be used and the best choice may depend on the nature of the problem. Entropy theory seems the most appealing to evaluate genomic information on discrete traits (as concluded from pilot studies, results not shown). Other loss functions such as the L₁-loss function or the misclassification rate could be implemented in an easy manner. Without loss of generality we show how to implement the entropy theory in the node splitting decision. The information gain (IG) for each covariate s drawn at random in a given node was calculated as described in Long et al. [9]:

where $N_k=N_k^{+}+N_k^{-}$, and $H(\mathrm{Pr}(y))=-{\displaystyle \sum _{y\in A}\mathrm{Pr}(y){\log }_2\mathrm{Pr}(y)}$ is the entropy of the probability distribution of y, and A is the set of all states that y can take ({0,1}). The SNP covariate with the highest IG at each node is used to split the node into two new child nodes, each one containing the individuals from the parent node with the risk or the non-risk allele, respectively.

There are two features involved in the RF algorithm that deserve further attention: the out of bag samples, and the variable importance.

Data sets

Simulated and field data sets were used for the model comparisons. Description of these data is given next.

Predictive ability

Performance of the models was based on predictive ability to correctly predict genetic susceptibility in the testing sets. The true genetic susceptibilities of individuals in the simulated data set are known. However, true genetic merits are unknown in the field data case. Therefore, predictive ability was evaluated in a different manner in the field data, as described below.

Simulated data set

Bayesian regression showed larger Pearson correlations than ensemble algorithms in the scenario with a larger number of QTL. Method BTL achieved the largest Pearson correlation (0.38), followed by TBA and LhB (0.33). Method RF showed the smallest Pearson correlation (0.22) in this simulated scenario and the largest AUC (0.72). This suggests that RF ranked individuals less accurately than other methods when a large number of QTL affects additively the trait, but the method is more accurate than other methods at discerning between healthy and affected individuals.

It must be pointed out that the simulated scenarios are purely additive and other more realistic scenarios with a more complex interaction between genes and biological pathways might provide different results.

Field data set

The three data sets had a disease occurrence of 50%. The relative predictive importance obtained with RF for each SNP covariate $x_j^l$ in each line is plotted in Figure 1. Many more SNP were identified as predictors of SH in line A than in line B and C, suggesting that many more genomic regions may be associated to SH in line A than in line B or C. Lines B and C showed few genomic regions with a large relative importance variable associated to the genetic resistance to SH. Thirty seven, four and six SNP had a larger relative variable importance than 50% in lines A, B and C, respectively. The odds ratio of SNP with VI > 50% ranged from 1.41 to 2.17 in line A, from 2.56 to 3.03 in line B and from 1.86 to 2.50 in line C, suggesting a considerable risk of being susceptible to SH of those animals carrying the unfavorable alleles. The SNP with the largest importance estimate (VI = 100%) in line C had also the maximum VI in line B, but had a VI < 21% in line A. These results suggest that the genetic variants presented in line B and C in this genomic region provide a relatively larger predictive ability of SH than genetic variants in the same genomic region in line A. The relative VI of the most important SNP in line A was lower than 2% in lines B and C, although other SNP in LD with those may have been detected in these lines. Fifty, 44 and 48 markers with VI greater than 99.5 percentile were found in lines A, B and C, respectively. Most represented chromosomes were SSC4, SSC7, SSC14 and SSC17 in line A, SSC1, SSC2, SSC6 and X chromosome in line B, and SSC8 in line C. Validation of these results and conclusions about their role in genetic or biological pathways should be performed on different populations and studies.

None of the methods was clearly preferred in the crossbred (line C), where similar phi correlations were found for RF, TBA and boosting, with larger robustness for LhB at α < 0.50. No differences were found between RF, TBA and LhB to correctly detect most extreme animals in the crossbred line. The Huber loss function was more robust than the squared sum of residuals at analyzing binary traits, in accordance with its resemblance with the L₁ loss function.

The true genetic architecture of SH is obviously unknown and no conclusions on its relationship with the performance of the different methods can be extracted. There was no clear relationship between the preferred method and the number of relevant genomic regions identified in each line (Figure 1). The choice of the model to be used in genome-wide prediction of traits like SH may depend on the interest of the breeder. For instance, detection of susceptible animals was done more accurately in line A using RF, whereas the Bayesian regressions were preferred in line B. Thus, a different method may be desired depending on the objective of the breeding program. The model with higher sensitivity would be preferred in a breeding program aiming at eradicating a given disease or trait. In a specifity+sensitivity scenario, RF was the best method at α = 0.05, 0.25 and 0.50, and it also showed the larger AUC values, regardless of the line.

Results showed that RF had the lowest risk, among methods used here, of misclassifying animals for low-medium heritability discrete traits in all lines, although all methods had considerable misclassification risks at α = 0.50. However, in a disease resistance genome-assisted prediction context, for instance, we are mainly interested in correctly detecting the most susceptible or resistant animals (lower α values), and RF seemed to perform slightly better than the Bayesian regressions to detect susceptibility to SH in this population, mainly in line A. Note that the threshold versions presented here incorporate n liability variables to be estimated in the model, increasing the parameterization of the models, and therefore hampering their predictive ability.

Results from the analyses of the crossbred line were not conclusive, as different behaviors between methods were found for different α values. This may be explained by the larger genetic heterogeneity expected in line C which may not be captured with only 5000 markers.

A small number of animals was used in the testing set and only punctual estimates are given here. This may be important at low α levels with a smaller number of records. Uncertainty about these estimates may be reported from their posterior densities [34] in the case of Bayesian methods and using bootstrap or cross-validation strategies in the case of this version of RF [11]. Uncertainties are not reported in this study because this data set aims at serving just as an example of three different models applied to discrete traits in a genome-assisted prediction context without overloading the discussion. Furthermore, the preferred model may be case-specific.

The misclassification rate and the logit function were also used as splitting criteria in RF but with poorer predictive ability (results not shown). Here, hyperparameters were set as fixed, although it is possible to assign them a prior distribution for their estimation [35]. Nonetheless, a minor improvement on predictive ability is expected if the ad-hoc choice of the parameters is within a sensible range of values.