## Abstract

Genomic selection (GS) has become a tool for selecting candidates in plant and animal breeding programs. In the case of quantitative traits, it is common to assume that the distribution of the response variable can be approximated by a normal distribution. However, it is known that the selection process leads to skewed distributions. There is vast statistical literature on skewed distributions, but the skew normal distribution is of particular interest in this research. This distribution includes a third parameter that drives the skewness, so that it generalizes the normal distribution. We propose an extension of the Bayesian whole-genome regression to skew normal distribution data in the context of GS applications, where usually the number of predictors vastly exceeds the sample size. However, it can also be applied when the number of predictors is smaller than the sample size. We used a stochastic representation of a skew normal random variable, which allows the implementation of standard Markov Chain Monte Carlo (MCMC) techniques to efficiently fit the proposed model. The predictive ability and goodness of fit of the proposed model were evaluated using simulated and real data, and the results were compared to those obtained by the Bayesian Ridge Regression model. Results indicate that the proposed model has a better fit and is as good as the conventional Bayesian Ridge Regression model for prediction, based on the DIC criterion and cross-validation, respectively. A computing program coded in the R statistical package and C programming language to fit the proposed model is available as supplementary material.

- Genomic Selection
- data augmentation
- asymmetric distributions
- GBLUP
- Ridge regression
- GenPred
- Shared Data Resources

In genetic studies of plants or animals, it is common to find quantitative traits whose distribution is not normal; this is because the data are obtained from multiple sources or contain isolated observations (Li *et al.* 2015). Landfors *et al.* (2011) noted that it is often necessary to normalize the data to remove variation introduced during the experiment’s development. However, such standard normalization techniques are not always able to remove bias because a large number of observations are positively or negatively affected by some treatment. In addition, suitable transformation for the data may be difficult to find, and may bring further complications when estimating and interpreting the results obtained (Fernandes *et al.* 2007). To avoid transformations, different methods have been developed that are flexible enough to represent the data and reduce unrealistic assumptions (Arellano-Valle *et al.* 2005). In the genomic selection framework (GS; Meuwissen *et al.* 2001), dense molecular marker genotypes and phenotypes are used to predict the genetic values of candidates for selection. The availability of high density molecular markers of many agricultural species, together with promising results from simulations (*e.g.*, Meuwissen *et al.* 2001) and empirical studies in plants (de los Campos *et al.* 2009; de los Campos *et al.* 2010; Crossa *et al.* 2010, 2011) and animals (*e.g.*, VanRaden 2008; Weigel *et al.* 2009), are prompting the adoption of GS in several breeding programs. The parametric model for GS explains phenotypes (; *i = 1,…,n*) as functions of marker genotypes () using a linear model of the form: , where is the number of phenotypic records, is the number of markers, represents the number of copies of a bi-allelic marker (*e.g.*, an SNP), and is the additive effect of the reference allele at the *j ^{th}* marker, . In matrix notation, the model is expressed as , where , and are vectors of phenotypes, marker effects and Gaussian zero mean errors, respectively, and is a matrix of marker genotypes of dimensions

*n*×

*p*. However, when the data are not normal, the normal regression methods generate inconsistent estimates with the natural distribution of the data and, therefore, the estimation of the conditional mean given the covariates is also inconsistent (Bianco

*et al.*2005).

With dense molecular markers, the number of markers exceeds the number of records in the reference population (* p > >n*) and, therefore, penalized regression estimation methods and their Bayesian counterparts are commonly used. Penalized estimation methods produce regression parameter estimates that are often equivalent to posterior modes. The literature on GS offers a long list of Bayesian models whose main difference is the prior distributions assigned to marker effects, which leads to what is known as the Bayesian Alphabet (Gianola 2013; de los Campos

*et al.*2013). The above-mentioned models assume that the response follows a normal distribution. Several phenotypic traits have distributions that are skewed, for example, female flowering, male flowering, the interval between male and female flowering, categorical measurements of diseases such as ordinal scale, counting data, etc. In these cases, either a normalizing transformation for the response variable (

*e.g.*, using Box-Cox transformation) or a model that deals with skew responses may be used. Varona

*et al.*(2008) proposed using linear mixed models with asymmetric distributions in the residuals to tackle the problem in the context of animal breeding when pedigree information is available. Nascimento

*et al.*(2017) proposed the Regularized Quantile Regression as a way to overcome the issue of non-symmetric distributions when marker information is available.

If a population is selected based on one trait and another trait of interest results that exceeds (does not exceed) some threshold, then the conditional distribution of , for a fixed leads to a distribution that is skewed (Arnold and Beaver 2000), such as the skew-normal (SN) distribution, which is of particular interest in this research. This distribution is a generalization of the normal distribution (Azzalini 1985) with a shape parameter added to adopt skewed forms. It has the advantage of being mathematically tractable and shares properties with the normal distribution; for example, the density of the SN is unimodal (Genton 2004). Varona *et al.* (2008) argues that the asymmetric distributions observed for the phenotypes are the result of environmental factors and that data can be modeled using non-symmetric residual distributions.

Based on the previous considerations and motivated by the fact that a great deal of traits in plant breeding have skew normal distributed, such as flowering time in most crop species, as well categorical traits such as diseases (ordinal, binomial, or counting data), in this study we propose a general Bayesian genomic regression model for skew-normal phenotypic traits with skew-normal random errors. The model uses a stochastic representation of the response variable (Arnold and Beaver 2000) in order to ease computations and it also works in the case that when *n*>*p*. It should point out, however, that the aim of the paper is not to describe and study the causes of the skew distribution but rather we assume that the skew data are given and thus the objective is to propose a robust statistical model that deals with the skew-normal distribution of the phenotypic and residuals.

The structure of this paper is as follows. In section 2, we present the statistical models and describe the latent variable model used in the regression. In section 3, we describe a simulation experiment that is performed to evaluate the predictive power of the proposed model. In section 4, we present an application with real data; section 5 includes the discussion and concluding remarks.

## MATERIALS AND METHODS

### Statistical models

In this section we introduce the statistical models to be used in the manuscript. We begin by giving a brief review of the skew normal model. Then we introduce the concept of data augmentation and we use this concept in order to generate a skew normal random variable. After that we introduce the “centered parameterization” in the skew normal model, regression with skew random errors. Finally, we present the pior, posterior and full conditional distributions for the regression model with skew normal residuals.

#### Skew-normal model:

A continuous random variable is said to follow the skew-normal law with shape parameter , denoted by if its density function is:(1)where and denote the density and cumulative distribution functions of a standard normal random variable, respectively. The subscript indicates the use of “direct parametrization” (Azzalini 1985). Note that the skew normal distribution reduces to the normal case when .

The mean and variance of are given by:..The coefficient of skewness of is:.If is a random variable defined by , then is said to have a skew-normal distribution with location (), scale (), and shape () parameters, and is denoted as . The density function of is given by:.It can be shown that the coefficient of skewness of corresponds to the skewness coefficient of .

#### Hidden truncation:

Let *V* and *W* be two random variables whose joint distribution is given as follows:,where denotes a bivariate random variable with mean and variance-covariance matrix and the random variable *U* is defined as follows:then , with (Arnold and Beaver 2000; Hea-Jung 2005; Liseo and Parisi 2013). The above representation allows writing an augmented likelihood function (Hea-Jung 2005; Liseo and Parisi 2013), “as if” we had observed the latent variable . The conditional distribution of is and , which is a truncated normal random variable with location parameter 0, scale parameter 1, lower truncation bound 0 and upper truncation bound . Therefore, the joint distribution of and is , that is:.(2)Note that the density function of can be obtained by integrating with respect to *z*; that is, .

Estimating the parameters in the direct parametrization is troublesome, so “centered parametrization” is more appropriate for parameter estimation and interpretation (Azzalini 1985; Pewsey 2000; Azevedo *et al.* 2011, among others). If and , then and so that . Figure 1 shows the density function of *U* for several values of the shape parameter and the corresponding values of . The random variable:,(3)is said to be a skew normal random variable with parameters , and , where is Pearson’s skewness coefficient given by , and the range of is (-0.99527, 0.99527). In this case, , . The usual notation is .

If we consider the following transformations:,(4)it can be shown, using Jacobians (Casella and Berger, 2002, Chapter 2), that the joint density of and is given by:(5)where . Note that the density function of can be obtained by integrating with respect to ; that is, , with . This representation is very convenient, because it allows us to write an augmented likelihood function (Hea-Jung 2005; Liseo and Parisi 2013), “as if” we had observed the latent value . The density function of the skew normal random variable under “centered parametrization” is a complicated function that was given by Azevedo *et al.* (2011):,where , , with , .

#### Regression with skew normal random errors:

Azzalini and Capitanio (1999) and Rusell and González (2002) proposed a simple linear regression model where the error terms are independent and identically distributed as . The proposed model is:;from the properties of the skew normal distribution, it follows that . The model can be easily extended to include more covariates; that is:.Azzalini and Capitanio (1999) and Rusell and González (2002) used the maximum likelihood method to estimate the parameters in the model. These ideas can be extended to the case of errors that are independent and identically distributed as .

#### Bayesian regression with skew normal random errors (BSN):

Let Then, the likelihood function is given by:.Let and the prior distribution of and a set of hyper-parameters that index the prior distributions. Then, by Bayes’ theorem, the joint posterior distribution of is as follows:Neither the joint posterior distribution nor the full conditional distributions of the parameters of interest have a closed form; therefore, the implementation of this model within the Bayesian framework is computing intensive. We propose using hidden truncation together with two standard MCMC techniques in Bayesian analysis: (i) Gibbs Sampling (Geman and Geman 1984) and (ii) Random Walk Metropolis Sampling to alleviate some of the computing burden.

#### Prior, posterior and full conditional distributions:

Consider the joint distribution of and Z given in (5). In the regression context, we set ; then the augmented likelihood function is:(6)In order to fully specify the Bayesian model, prior distributions for the unknown parameters must be defined. Let ; for based on the following transformation , where the prior for is denoted as , and depending on the hyper-parameters and , it can lead to a rich variety of shapes, just as in the case of the Beta distribution. For the intercept, the prior distribution is ; for the scale parameter, a scaled inverted chi-squared prior distribution was used, that is, and finally (see Sorensen and Gianola, 2002, p. 85, for details about the parametrization used in this paper). The joint prior distribution is.(7)By combining equations 6 and 7 throught the Bayes’ theorem, the posterior distribution of is given by:(8)The full conditional distributions of the parameters are given in Appendix A, which are the inputs for the Gibbs and the Random Walk Metropolis Samplers. In Appendix B, some pragmatical rules to elicitate values for the hyper-parameters , and , are given. The rules for setting , are based on those given by de los Campos *et al.* (2013) and Pérez and de los Campos (2014). In this paper, we set the hyper-parameters as follows: , ,, , where is the sample phenotypic variance and . We set in order to reduceshrinkage and because in practice it mimics a non informative but proper distribution. To sample from the full conditionals of and , we implemented a Random Walk Metropolis Sampler whose parameters are tuned so that the acceptation rate is about 0.23 (see Appendix A for details).

The BSN can be re-parametrized by replacing with ; if the prior distribution of marker effects is normal with mean 0 and variance , then the prior of is , which leads to a G-BLUP model (see de los Campos *et al.* 2013, for details about G-BLUP) but with skew normal residuals, that is, or, in matrix notation, , which is a skew linear mixed model, a particular case of the model proposed by Arellano-Valle *et al.* (2005, 2007), who relaxed all normality assumptions in a standard mixed model.

#### Bayesian ridge regression With random normal errors (BRR):

Regression with random normal errors is a special case of the proposed model when . The model is widely used in the GS selection literature (*e.g.*, de los Campos *et al.* 2013). In the GS context, the model is given by:,where are independent and identically distributed as .

The prior distributions for the unknown are: , for the scale parameter, a scaled inverted chi-squared distribution, that is, and finally .

The joint prior distribution is.(9)Thus, the posterior distribution of isThe required full conditional distributions of the parameters for implementing a Gibbs sampler can be found elsewhere (*e.g.*, de los Campos *et al.* 2013). We set the hyper-parameters using the same rules as in the BSN model. The BRR model can be fitted easily using the BGLR statistical package (Pérez and de los Campos 2014).

### Monte Carlo Simulation

In this section, we use simulated data using marker genotypes from a wheat dataset made publicly available by Crossa *et al.* (2010). The dataset includes genotypic information for 599 wheat lines which were genotyped for 1279 DArT markers coded as 0 and 1. We simulated the phenotypes using the following additive genetic model:,(9)where , with , , which leads to different degrees of skewness. The intercept parameter, was set equal to 3; 10 marker effects were sampled from a normal distribution with mean 0 and variance 0.5/10 (Pérez and de los Campos 2014), and the rest were set equal to 0, that is:The idea here is to verify, through simulation, whether the proposed model works satisfactorily. We therefore obtained point estimates for , , and . We also fitted the Bayesian Ridge Regression model and compared the estimates of regression coefficients, predictions and estimates of genetic values of both models. Let be the vector of posterior means for regression coefficients. Pearson’s correlation between the observed () and predicted values () is a goodness-of-fit measure; Pearson’s correlation between the “true” genetic values () and the predicted values () is a measure of how well the genetic values are estimated; finally, Pearson’s correlation between the “true” marker effects () and the estimated effects () is a measure that indicates how good a model is at uncovering marker effects (de los Campos *et al.* 2009). We also computed the effective number of parameters () and deviance information criterion () for the two fitted models (see Spiegelhalter *et al.* 2002, for more details).

The algorithm used in this simulation experiment is described briefly below.

Set , , and .

Simulate the phenotypes using equation (9).

Fit the regression model with skew normal random errors and obtain point estimates for , , and , that is, , , and . The point estimates correspond to the posterior means of the posterior distribution of the parameters of interest.

Fit the Bayesian Ridge Regression model and obtain point estimates for , , , that is, , , .

Compute the correlation between observed and predicted phenotypes, “true” and predicted genetic values, and “true” and estimated regression coefficients with both regression models.

Compute the effective number of parameters () and deviance information criterion () for the two fitted models.

Repeat steps 1 to 5 one hundred times and obtain the averages of correlations, intercept (), and .

### Application to real data

This dataset is from the Drought Tolerance Maize (DTMA) project of CIMMYT’s Global Maize Program (http://www.cimmyt.org). The dataset comes from a large study aimed at detecting chromosomal regions affecting drought tolerance. The genotypic data consist of information from 300 tropical inbred lines that were genotyped using 1,152 SNPs (Single Nucleotide Polymorphisms). The analyzed trait is Gray Leaf Spot (GLS) caused by the fungus *Cercospora zeae-maydis*, which was evaluated at three different sites, Kakamega (Kenya), San Pedro Lagunillas (Mexico) and Santa Catalina (Colombia) (see Supporting information). Crossa *et al.* (2011) analyzed a subset of these data; the response variable was transformed using Box-Cox transformation (Box and Cox 1964). Figure 2 shows density plots for GLS rating at the three sites. Kernel density was estimated using a Gaussian kernel, and the bandwidth for the kernel was estimated according to Venables and Ripley (2002). Figure 2 also shows the sample skewness index, , where , is the sample mean and is the sample standard deviation (see Joanes and Gill 1998); in the three cases, the distribution is skewed to the right, so most of the distribution is concentrated around small values of the response variable.

We propose using the regression model with skew normal random errors to predict disease resistance. We fitted two models: (1) the standard Bayesian Ridge Regression, where the errors , where “NIID” stands for “normally, independent and identically distributed”; and (2) the proposed model with skew normal random errors. The Bayesian Ridge Regression was fitted using the BGLR package (Pérez and de los Campos 2014), whereas the proposed model was fitted using the algorithm described in Appendix A. The models were first fitted using the full data, and subsequently 100 random partitions with 80% of observations in the training set and 20% of observations in the testing set were generated. The two models were fitted for each of these random partitions; then the phenotypes of the individuals in the testing set were predicted and the ability of each model to make predictions was evaluated using Pearson’s correlation between observed and predicted values. Inferences for each fit were based on 100,000 samples obtained after discarding 50,000 samples that were taken as burn-in. Convergence was checked by inspecting trace plots of the parameters.

### Data and program availability

The data and programs are available as File S1 which corresponds to a compressed zip folder. The zip folder also contains a description of the data and commands to read it into the R statistical software.

## RESULTS

### Monte Carlo Simulation

Table 1 shows point estimates for , , and for the BSN and BRR models for different values of . It also shows an estimate of , a regularization parameter that is widely used in Bayesian Ridge Regression. Higher values of the parameter are associated with more shrinkage; note that the estimates of are very similar in both models, so small values of could be associated with more precise estimates of . It is also clear from this table that the point estimates for and are very close to the real values used in the simulation. The correlation between observed and predicted values and the mean squared error is quite similar for both models and there is no clear winner. Finally, the algorithm is not able to estimate precisely the parameter for distributions that are slightly asymmetric.

Table 2 shows the effective number of parameters , the deviance information criterion (*DIC*), the correlation between “true” and estimated marker effects and the correlation between “true” and estimated signals. The table also shows that in general the and the *DIC* (small is better) favored the BSN model. The correlation between “true” and estimated marker effects is slightly better for BSN and the difference between the two models becomes clearer as increases. The same pattern is observed for the correlations between true and estimated genetic signals.

### Application to real data

#### Full data:

Table 3 shows estimates of the posterior means of parameters , and , as well as the effective number of parameters and the deviance information criterion (*DIC*). From Table 3 it is clear that the estimation of marker effects is more precise for the BSN model than for the BRR model; the and the *DIC* also favored the BSN model. The estimated parameter also supports the assumption that the skew normal random error is correct, and that the point estimate is not around 0, except in the case of San Pedro Lagunillas.

Figure 3 shows scatterplots of the predicted GLS using the BSN and BRR models. As expected, Pearson’s correlation between both predictions is very high (higher than 0.95). That implies that even when the data are skewed, if a BRR model is fitted in order to obtain candidates for selection, we can expect to obtain about the same individuals. Two models were fitted for each site by BRR and BSN.

#### Cross-validation:

Figure 4 shows scatterplots for Pearson’s correlation between observed and predicted values for individuals in the testing set obtained after fitting the BSN and BRR models for the three locations. When the correlations are higher for BSN than for BRR, this is represented by a filled circle, and by an open circle otherwise. The figure also shows the number of times Pearson’s correlation is higher for the BSN than for the BRR model. From this figure, it is clear that the BSN model predicts slightly better than the BRR model. Figure 5 shows a scatterplot for the mean squared errors in the testing set for the three locations. When the MSE in BSN is smaller than the MSE in BRR, this is represented by an open circle and by a filled circle otherwise. The number of times that the MSE in BRR is greater than the MSE in BSN is also shown in the plots. From this figure, it is clear that in general, the MSE for BRR is greater than the MSE for BSN. Table 4 shows the average Pearson’s correlation and mean squared error (MSE) between observed and predicted values in the testing set. The averages and the standard deviations are very similar for both models and the differences between the models are non-significant, but the figures suggest that the BSN model predicts slightly better than the BRR model.

## DISCUSSION AND CONCLUSIONS

We have proposed a Bayesian regression model for skewed responses with applications when *p > >n* in the GS context, but it can also be employed in other cases and, of course, when *p < n*. In addition, to generalize linear whole genome regression models for various discrete distributions (ordinal, binomial, etc.), this study further completes the Bayesian toolbox for whole genome regression. The proposed model uses a stochastic representation of a skew normal random variable in order to facilitate the computations; it also allows using standard MCMC techniques to fit the proposed model. Results of the simulation and of applications with real data suggest that the proposed model fits the data better and also predicts slightly better than the standard Ridge Regression model. The Ridge Regression model is a particular case of our model when . On the other hand, our results also suggest that BRR is a very robust model, although in the simulations data we already knew that it was the wrong model to fit; still, the predictive power of the model was very good. Although the conventional Bayesian whole-genome regression is robust, it does not correctly deal with skew phenotypic data, and this can decrease its genomic-enabled prediction accuracy and its goodness of fit to the data. Thus, the advantages of the proposed Bayesian whole-genome regression compensate its complexity and possible increases in computational time as compared to the conventional Bayesian ridge regression. The model proposed in this study is conceptually and operationally different, and presumably simpler than the skew-normal linear mixed model of Arellano-Valle *et al.* (2005) that uses a multivariate skew-normal distribution in order to relax normality.

Despite the fact that skewness is a major concern for breeding data analyses and may often be a result of uneven sampling of “high” and “low” performing individuals, selection, environmental effects, etc., the theoretical developments presented in this study are also applicable to many other areas of research in agronomy and in agriculture in general. For example, most crop flowering time data are indeed skewed, as well as categorical data representing different types of diseases as those presented in this research. So, skewness in phenotypical response can be the result of an artificial phenomena, the aim of this study was to propose a statistical model that will be more appropriate to deal with that problem.

Results of this study can be compared to results of two other studies, Crossa *et al.* (2011) and González-Camacho *et al.* (2012). Crossa *et al.* (2011) included one site in Mexico (San Pedro Lagunillas) that was also analyzed by transforming the original GLS ordinal scale using Box-Cox transformation; the prediction accuracy of different models (*e.g.*, Bayesian Lasso and Reproducing Kernel Hilbert Spaces) ranged from 0.416 to 0.462. Although strict comparison with the results obtained in this study is not possible because other random cross-validations were generated, the prediction accuracies of BSN (0.5489) and BRR (0.5450) models were higher than those previously reported by Crossa *et al.* (2011) for the same site.

Stochastic representation can be used to extend Reproducing Kernel Hilbert Space (de los Campos *et al.* 2010) models that in many empirical studies have led to more accurate predictions than Bayesian Ridge Regression models and Bayesian LASSO, among others (*e.g.*, Pérez-Rodríguez *et al.* 2012), so this is a topic for future research. Further studies to extend the model proposed in this study to include genotype × environment interaction should not be complicated. The proposed model can also be extended by assigning different prior distributions to the marker effects, for example, to induce variable selection. This could potentially lead to a new Bayesian alphabet with skew normal random errors.

## Acknowledgments

We thank the researchers of the Global Maize Program of the International Maize and Wheat Improvement Center (CIMMYT) who carried out the maize trials and provided the phenotypic and genotypic data analyzed in this article.

Authors’ contributions: PPR, RAP, SPE, CVC, JSE and JC conceived the idea and designed the study; RAP, JC and PPR developed the statistical models; RAP wrote the R scripts; RAP and PPR carried out the statistical analysis; RAP, JC and PPR wrote the first version of the manuscript; JC, SPE, CVC, and JSE contributed critical comments and to the writing of the final manuscript.

## APPENDIX A: FULL CONDITIONAL DISTRIBUTIONS

The full conditional distributions of the unknown parameters in the Bayesian Genomic Regression Model with Skew Normal Random Errors are given below. The derivation uses results from Bayesian linear models (*e.g.*, Sorensen and Gianola, 2002), as well as results from Rusell and González (2002) and Liseo and Parisi (2013).

The joint posterior distribution is as follows:where and , **.**

### 1. Intercept

where . Note that the right-hand side of the last expression is the kernel of a normal distribution with mean and variance with .

### 2. Regression coefficients

Here we obtain the conditional distribution of each of the elements of the vector , *i.e.*, , .where , is the *i*-th row of without the *j*-th column and is the vector without the *j*-th element. After some algebraic manipulations, it can be shown that the right-hand side of the above expression corresponds to the kernel of a normal distribution with mean and variance , with , .

### 3. Latent variables

where . After simplifying some terms, the right-hand side of the above equation can be written as:.The above expression is the kernel of a truncated normal distribution with location parameter and scale parameter , lower truncation bound 0 and upper truncation bound . In this work, we used the R library truncnorm (Trautmann *et al.*, 2014) to sample from this distribution.

### 4. Residual variance ()

Note that is a very complex function of and its kernel does not correspond to any known univariate density function, so we have to sample this parameter by Metropolis algorithm or another MCMC technique. Following Hea-Jung (2005), we considered the Random Walk Metropolis Algorithm with a de-constraint transformation to sample . Given the fact that let so that has support in . Note that the density function of can be obtained by using the transformation method (Casella and Berger, 2002, Chapter 2), which is given by:.In the Random Walk Metropolis Algorithm, we generated by choosing a proposal transition kernel to add noise to the current state. Assuming that the actual value of is , and that we wanted to update its value so that in the next iteration we had , we followed steps a-c below.

a. Sample , , where .

b. Sample , .

c. If then ; otherwise .

Once has been obtained, compute . The parameter can be modified so that we have an optimal acceptation rate (Roberts *et al.*, 1997).

### 5. Variance of markers ()

### 6. Correlation coefficient ()

Note that is a very complex function of and its kernel does not correspond to any known univariate density function, so we must obtain samples using the Metropolis algorithm or another MCMC technique. Here we propose using Fisher’s (1915) transformation of defined as , so that has support in . Note that the density function of can be obtained using the transformation method (Casella and Berger, 2002, Chapter 2), which is given by.In the Random Walk Metropolis Algorithm, we generated by choosing a proposal transition kernel to add noise to the current state. Assuming that the actual value of is and that we wanted to update its value so that in the next iteration we had , we followed steps d-f below.

d. Sample , , where .

e. Sample , .

f. If , then ; otherwise .

Once has been obtained, compute . The parameter can be modified so that we have an optimal acceptation rate (Roberts *et al.*, 1997).

The samples from the posterior distribution can be obtained using the Gibbs Sampler (Geman and Geman, 1984) and the Random Walk Metropolis algorithm. In the algorithm, we sampled each of the fully conditional distributions until we obtained a sample of the desired size. We implemented the algorithm in the R Statistical package (R Core Team, 2016). In order to speed up computations, the routines that sample from , were implemented in C programming language (Kernighan and Ritchie, 1988), a shared library was generated and then the compiled routines were used in R. The R script and the C source code are available upon request from the first author.

## APPENDIX B: SETTING THE HYPER-PARAMETERS FOR THE PRIOR DISTRIBUTIONS

The hyper-parameters can be set using a set of default rules similar to those used in the BGLR software (de los Campos *et al.*, 2013; Pérez and de los Campos 2014). With these rules we assigned proper but weakly informative prior distributions so that we partitioned the total variance of the phenotypes into two components: (1) the error and (2) the linear predictor, that is:

where . A priori, the total genetic variance is and the a priori average genetic variance is:(B.2)where . So from (B.1) and (B.2), the partition of the phenotypic variance is given by:(B.3)where .

### Setting the hyper-parameters for

In the parametrization used in this work, and . We set =5 so that the prior distribution has a finite mean. From equation (B.2),,(B.4)thus, if we replace the left-hand side of (B.4) with the mode of then:.(B.5)From (B.5) and solving for , we obtain . From the definition of heritability, , so that ; then:(B.6)Once we set we can set , and we only need to compute the phenotypic variance , and set as the proportion of the variance that is explained a priori by the markers. By default we set .

### Setting the hyper-parameters for

In the parametrization used in this work, and ,

we set and .

### Setting the hyper-parameters for

We set so that the prior assigned to the intercept is effectively a flat one.

### Setting the hyper-parameters for

The density function of the prior assigned to is:Figure B.1 shows the density function of for different values of the parameters . Setting different values for those hyper-parameters could lead to a rich variety of shapes, just as in the Beta distribution. Note that if we set , then we obtain a prior that corresponds to the one used in this work.

## Footnotes

*Communicating editor: J.-L. Jannink*Supplemental material is available online at www.g3journal.org/lookup/suppl/doi:10.1534/g3.117.300406/-/DC1.

- Received November 2, 2017.
- Accepted March 23, 2018.

- Copyright © 2018 Pérez-Rodríguez
*et al.*

This is an open-access article distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.