An R Package for Bayesian Analysis of Multi-environment and Multi-trait Multi-environment Data for Genome-Based Prediction

Evidence that genomic selection (GS) is a technology that is revolutionizing plant breeding continues to grow. However, it is very well documented that its success strongly depends on statistical models, which are used by GS to perform predictions of candidate genotypes that were not phenotyped. Because there is no universally better model for prediction and models for each type of response variable are needed (continuous, binary, ordinal, count, etc.), an active area of research aims to develop statistical models for the prediction of univariate and multivariate traits in GS. However, most of the models developed so far are for univariate and continuous (Gaussian) traits. Therefore, to overcome the lack of multivariate statistical models for genome-based prediction by improving the original version of the BMTME, we propose an improved Bayesian multi-trait and multi-environment (BMTME) R package for analyzing breeding data with multiple traits and multiple environments. We also introduce Bayesian multi-output regressor stacking (BMORS) functions that are considerably efficient in terms of computational resources. The package allows parameter estimation and evaluates the prediction performance of multi-trait and multi-environment data in a reliable, efficient and user-friendly way. We illustrate the use of the BMTME with real toy datasets to show all the facilities that the software offers the user. However, for large datasets, the BME() and BMTME() functions of the BMTME R package are very intense in terms of computing time; on the other hand, less intensive computing is required with BMORS functions BMORS() and BMORS_Env() that are also included in the BMTME package.

multienvironment multi-trait genome-based prediction and selection R-software multivariate analysis GenPred Shared data resources Genomic Prediction Genomic selection (GS) is a methodology used in plant breeding that was proposed by Meuwissen et al. (2001). It is a type of marker-assisted selection that consists of genotyping and phenotyping a training sample (reference population); with the help of statistical models, predictions of genomic estimated breeding values (GEBV) or phenotypic values of the testing sample (validation population) are obtained for which only genome-wide dense genetic marker data were available. GS does not depend on prior knowledge about a few, large-effect genes or QTL, since all markers are used simultaneously in the training of the statistical models. GS was first used in animal breeding (Hayes and Goddard 2010), but nowadays is being implemented in many crops, for example, maize (Crossa et al., 2014), cassava (de Oliveira et al., 2012), wheat (Rutkoski et al., 2011), sugar beet (Würschum et al., 2013), tomato (Yamamoto et al., 2016), rice (Spindel et al., 2015), apple (Kumar et al., 2012), pea (Burstin et al., 2015), cranberry (Covarrubias-Pazaran et al., 2018) and many others.
In recent years, an active area of research has begun to develop and improve existing statistical models for genomic selection (GS) due to the fact that successful GS implementation is strongly related to the accuracy of the predictions performed by statistical models. However, because there are no universally superior machines for prediction, many models have been proposed. For example, most of the proposed models are univariate and few are multivariate. Most of the univariate models are appropriate for continuous and Gaussian phenotypes, but there are several appropriate models for binary, ordinal and count traits. Some examples of implementations of models for non-Gaussian, non-continuous traits are unordered categorical (Heuer et al., 2016), binomial (Technow and Melchinger 2013) and ordinal categorical (Montesinos-López et al., 2015a,b). While multivariate models are used almost routinely nowadays, for the joint analysis of multiple-traits (e.g., Jia and Jannink 2012) as well as multiple-environments (e.g., Burgueño et al., 2012) and even multiple populations (e.g., Olson et al., 2012), there are few multivariate practical software programs for continuous and Gaussian phenotypes and there are scarcely any models and software for other types of response variables. To the best of our knowledge, almost none of the currently reported models consider mixed phenotypes including continuous, binary, ordinal, count, etc. traits. For this reason, it is clear that to increase the power of GS technology, it is of paramount importance to develop more models and improve the existing ones.
Multi-trait models in GS have been applied by many scientists. For example, Calus and Veerkamp (2011), Jia and Jannink (2012), Jiang et al. (2015), He et al. (2016), Schulthess et al. (2017), and Covarrubias-Pazaran et al. (2018) reported that multi-trait analysis outperforms unitrait analysis in terms of prediction accuracy and that the larger the correlation between traits, the larger the benefit of multi-trait analysis. The Multi-Trait Model (MTM) of de los Campos and Grüneberg (2016) is a mixed multi-trait Gaussian model under the Bayesian framework that uses a Gibbs sampler for inferences. Furthermore, Bayesian multi-output regressor stacking (BMORS) is a Bayesian version of the multi-trait regressor stacking method proposed by Spyromitros-Xioufis et al. (2012;2016). The training of BMORS has two stages: (1) a single univariate model is implemented using the GBLUP model, and (2) the resulting predictions are directly included by BMORS in an additional training stage. Thus, the concept of BMORS is that a second-stage model will correct the predictions of the first-stage model [using the predictions of the first-stage univariate GBLUP model (Spyromitros-Xioufis et al., 2012;2016)]. Montesinos-Lopez et al. (2016) were the first to develop a comprehensive theory for a Bayesian multi-trait multi-environment (BMTME) model for genome-based prediction. An improved version of BMTME allows general covariance matrices by using the matrix normal distribution that facilitates easy derivation of all full conditional distributions and permits a more efficient model in terms of time of implementation Montesinos-López et al. (2018a,b,c). In general, the matrix normal distribution model considerably improved in terms of implementation time over the time required by the original BMTME. Also, the Gibbs sampler for implementing the new BMTME model can be found in Montesinos-López et al. (2018a), and the priors of the model are given in detail in Montesinos-López et al. (2018b). Montesinos-López et al.
(2018a) provide the appropriate notations used for the matrix-variate normal distribution that is a generalization of the multivariate normal distributions to matrices. This plays a key role in building the BMTME model. The original software used by Montesinos-Lopez et al. (2016) to fit the BMTME was the first attempt to implement the multi-trait multi-environment theory when analyzing real data; however, the lack of the necessary optimization algorithms for efficiently applying the software made the original BMTME difficult to apply to real data.
It is also important to point out that even though the existing R statistical software for Bayesian analysis like 'stan' (https://mc-stan.org/) and 'JAGS' (https://en.wikipedia.org/wiki/Just_another_Gibbs_sampler) are very flexible for implementing Bayesian analysis, they are not user-friendly because the user needs a certain level of programming skills to correctly implement them (Stan Development Team 2018;Plummer 2018). These two software programs (stan and JAGS) also require more computational resources for their implementation since they are built not with conjugate priors. It is documented that multivariate analysis improves parameter estimation (Schulthess et al., 2017). For this reason, we agree with Castro et al. (2013) and Huang et al. (2015), who stated that multi-trait analysis is a powerful tool for clarifying the relationship and the effect of each studied variable and for building more efficient prediction models.
Due to the background of plant breeders, not only are new models needed, but the existing ones need to be improved. We also need reliable, efficient, user-friendly software in which breeders can implement the existing GS models. One popular R package in the context of genomic selection for continuous and ordinal data are the BGLR package of Pérez and de los Campos (2014) that was built under the Bayesian framework and is very flexible because it allows the use of a genomic relationship matrix (derived from marker or pedigree), and also allows implementing various methods like BayesA, BayesB, Bayes Lasso, Bayes Ridge and GBLUP and can deal with moderate datasets; however, it only allows the implementation of univaritate models. Therefore, to contribute to this requirement, we developed a Bayesian multi-trait and multienvironment (BMTME) R software that allows the implementation of multi-trait and multi-environment data for performing parameter estimates and evaluating the prediction performance of multiple traits that are studied in many environments. This BMTME package is different from existing ones [sommer (Covarrubias-Pazaran 2016), BGGE (Granato et al., 2018), ASREML (Gilmour et al., 1995) and MCMCglmm (Hadfield et al., 2010)] because it takes into account the genetic correlation between traits and between environments. The main difference of BMTME with sommer and ASREML is that our package was built under a Bayesian framework, while sommer and ASREML were based on a classical approach using restricted maximum likelihood. The difference between BGGE and our model is that our model is not only for multienvironment data but rather for multi-environment and multi-trait data simultaneously. On the other hand, the MCMCglmm package only allows a general covariance matrix for traits but not for environments, like the proposed BMTME package; however, it is important to point out that the MCMCglmm package allows modeling not only continuous responses but also binary, ordinal and counts.
The main objective of this research was to illustrate the application of the new BMTME with two real toy datasets; with these we show how to use the functions available in the BMTME package for implementing multienvironment (BME function), multi-trait and multi-environment data (BMTME function), as well as the Bayesian multi-output regressor stacking functions BMORS () and BMORS_ENV (). These two functions are very different to what the existing software [sommer (Covarrubias-Pazaran 2016), BGGE (Granato et al., 2018), ASREML (Gilmour et al., 1995) and MCMCglmm (Hadfield et al., 2010)] implements, since the theory behind this function is that of stacking methods. Stacking methods consist of training multiple learning algorithms for the same dataset and then combining the predictions to obtain the final predictions. In this study we used the initial BMTME of Montesinos-Lopez et al. (2016) but improved it by using the matrix variate normal distribution (Montesinos-López et al., 2018c) and the appropriate priors given by Montesinos-López et al. (2018a) and Montesinos-López et al. (2018b).

Statistical models
Multiple-environment Genomic Best Linear Unbiased Predictor (GBLUP) model: Since genotype · environment interaction is of paramount importance in plant breeding, the following univariate linear mixed model is often used for each trait: where y ij represents the response of the jth line in the ith environment (i ¼ 1; 2; . . . ; I, j ¼ 1; 2; . . . ; JÞ. E i represents the effect of the ith environment, g j represents the random genomic effect of the jth line, with g ¼ ðg 1 ; . . . ; g J Þ T $ Nð0; s 2 1 G g Þ; s 2 1 is a genomic variance, G g is of order J · J and represents the genomic relationship matrix (GRM) and is calculated (VanRaden 2008) as G g ¼ WW T p , where p denotes the number of markers and W is the matrix of markers of order J · p. The G g matrix is constructed using the observed similarity at the genomic level between lines, rather than the expected similarity based on pedigree. Further, gE ij is the random interaction term between the genomic effect of the jth line and the ith environment with gE ¼ ðgE 11 ; . . . ; gE IJ Þ T $ Nð0; s 2 2 I I 5GÞ, where s 2 2 is an interaction variance, and e ij is a random residual associated with the jth line in the ith environment distributed as Nð0; s 2 Þ where s 2 is the residual variance.
Bayesian multiple-trait multiple-environment (BMTME) model: The current BMTME model was implemented by Montesinos-López et al. (2018a,b,c). For a complete understanding of its description, first we provide the notations used for the matrix-variate normal distribution that plays a key role in building the BMTME model. Matrix-variate normal distribution is a generalization of the multivariate normal distribution to matrices. The (n·p) random matrix, M, has a matrixvariate normal distribution denoted as M $ NM n·p ðH; V; ΣÞ, if and only if, the (np·1) random vector vecðMÞ is distributed as multivariate normal as N np ðvecðHÞ; Σ5VÞ; therefore, NM n·p denotes the (n · p) dimensional matrix-variate normal distribution, H is a (n · p) location matrix, Σ is a (p · p) first covariance matrix, and V is a (n · n) second covariance matrix (Srivastava and Khatri 1979). vecð:Þ and 5 are the standard vector operator and Kronecker product, respectively. Unlike in a multivariate normal model where the data are concatenated into a single vector of length np, in a matrix-variate normal model, the data ðMÞ are in an n·p matrix where each column is a trait (Montesinos-López et al., 2018a). Therefore, the proposed BMTME model is defined as follows: where Y is of order n · L, with L the number of traits and n ¼ J · I, where J denotes the number of lines and I the number of environments, X is of order n · I, b is of order I · L, since b ¼ fb il g for i ¼ 1; ::; I and l ¼ 1; ::; L; Z 1 is of order n · J, b 1 is of order J · L and contains the genotype · trait interaction term since b 1 ¼ fgt jl g where gt jl is the effect of the genotype · trait interaction term for l ¼ 1; ::; J and for j ¼ 1; ::; L. Z 2 is of order n · IJ, b 2 is of order IJ · L and contains the genotype · environment · trait interaction, since b 2 ¼ fgEt jil g, where gEt jil is the effect of genotype · environment · trait interaction for j ¼ 1; ::; J, for i ¼ 1; ::; I and for l ¼ 1; ::; L. Vector b 1 is distributed under a matrix-variate normal distribution as NM J·L ð0; G g ; Σ t Þ; where G g is of order J · J and represents the Genomic Relationship Matrix (GRM) and is calculated using the VanRaden (2008) method as G g ¼ WW T p , where p denotes the number of markers and W the matrix of markers of order J · p; and Σ t is the unstructured genetic (co)variance matrix of traits of order where Σ E is an unstructured (co)variance matrix of order I · I and E is the matrix of residuals of order n · L with E $ NM n·L ð0; I n ; R e Þ, where R e is the unstructured residual (co)variance matrix of traits of order L · L, and G g is the genomic relationship matrix described above (Montesinos-López et al., 2018a). The BMTME model resulting from equation (2)  Next, we used the modified version of the Gibbs sampler of the original BMTME model proposed by Montesinos-Lopez et al. (2016) that was implemented in Montesinos-López et al. (2018a). It is important to point out that model (2) takes into account the genotype · environment terms in the (Z 2 b 2 Þ term and, for comparison purposes, we also ran the model in equation (2) but without the (Z 2 b 2 Þ term to study the effect on prediction performance with and without the genotype · environment term. The Gibbs sampler for implementing the BMTME model is found in Montesinos-López et al. Bayesian multi-output regressor stacking (BMORS): The proposed BMORS is a Bayesian version of the multi-trait regressor stacking method proposed by Spyromitros-Xioufis et al. (2012;2016). The training of BMORS consists of two stages. In the first stage, L single univariate models are implemented using the GBLUP model given in equation (1), but instead of using the resulting predictions directly as the final output, the BMORS includes an additional training stage where a second set of L meta-models are implemented for each of the L traits under study. Each meta-model is implemented with the following model: where the covariatesẐ 1ij ;Ẑ 2ij ; . . . ;Ẑ Lij represent the scaled predictions of each trait obtained with the GBLUP model in the first-stage analysis, and b 1 ; . . . ; b L are the beta coefficients for each covariate. The scaling of each prediction was performed by subtracting its mean (m lij ) and dividing by its corresponding standard deviation (s lij ), that is,Ẑ lij =ðŷ lij 2 m lij Þs 21 lij , for each l ¼ 1; . . . ; L. Therefore, the BMORS model contains as predictor information the scaled predictions of its response variables yielded by the first-stage models. In other words, the BMORS model is based on the idea that a secondstage model is able to correct the predictions of a first-stage model using information about the predictions of other first-stage models (Spyromitros-Xioufis et al., 2012;2016).

Real toy datasets
Mada dataset: This dataset was obtained from the study by Ben Hassen et al. (2018). The dataset is composed of a sample of 188 wheat lines evaluated for six traits. Each of the lines was evaluated in one environment. The lines were genotyped and 32,066 single nucleotide polymorphisms (SNPs) were obtained with a heterozygosity rate , 5% and a minor allele frequency (MAF) . 5%. A subset of the data were included in the package that includes 30 lines, and we named this dataset Mada. For more details, see the study by Ben Hassen et al. (2018). Raw markers are not included, and we provide the genomic relationship matrix (GRM) calculated according to the method of VanRaden (2008).
Maize dataset: This dataset was obtained from the study by Montesinos-Lopez et al. (2016). It consists of a sample of 309 maize lines evaluated for three traits: anthesis-silking interval (ASI), plant height (PH), and grain yield (GY). Each trait was evaluated in three optimal environments (Env1, Env2 and Env3). The lines were genotyped, 681,257 single nucleotide polymorphisms (SNPs) were obtained, and markers with more than 20% missing values were removed. After that, markers were imputed using observed allelic frequencies, and markers with MAF , 0.05 were removed, so that at the end of the quality control and imputation, 158,281 SNPs were still available for further analyses. To load this dataset in the package, we used only 30 lines, and we named this dataset Maize. For more details, see the study by Montesinos-Lopez et al. (2016).

Evaluation of prediction performance
We implemented cross-validation (CV) to evaluate the prediction performance. Two types of CV were implemented: K-fold cross-validation and random cross-validation.
K-fold cross-validation: Under this CV, the dataset was partitioned into K subsamples of equal size; each time K-1 of them were used for training (TRN) and the remaining one for testing (TST). In this CV, one observation cannot appear in more than one fold. In the design, some lines can be evaluated in some, but not all, target environments, which mimics a prediction problem faced by breeders in incomplete field trials. This CV strategy is exactly the same as the strategy denoted as CV2 that was proposed and implemented by Jarquín et al. (2017), where a certain portion of test lines (genotypes) in a certain portion of test environments is predicted, since some test lines that were evaluated in some test environments are assumed to be missing in others.
Random cross-validation: This CV strategy randomly splits the dataset into training (TRN) and testing data (TST). For each such split, the model is fitted to the TRN data, and predictive accuracy is assessed using the TST data. Since we used sampling with replacement, one observation may appear in more than one partition. The implemented CV mimics a prediction problem faced by breeders in incomplete field trials, where some lines may be evaluated in some, but not all, target environments. Since N ¼ J · I denotes the total number of records per each available trait, then to select lines in the TST dataset, we fixed the percentage of data to be used for TST [PTesting]. Then PTesting·N (lines) were chosen at random, and subsequently for each of these lines, one environment was randomly picked from I environments. The cells selected through this algorithm were allocated to the TST dataset, while the cells (ijÞ that were not selected were assigned to the TRN dataset. Lines were sampled with replacement if J , PTesting · N, and without replacement otherwise (López-Cruz et al., 2015). The metrics used to measure the prediction accuracy under both CV strategies were Pearson's correlation and the mean arctan absolute percentage error (MAAPE), which has the advantage that no zero estimates are produced when the response variable contains many zeros. They were calculated from each trait-environment combination for each of the testing sets and the average of all random partitions (folds) is reported as a measure of prediction performance.

Data availability
The data used in this study are included in the BMTME package, so once that package is installed, the datasets can be loaded into the R environment.
Installation of the BMTME package The aim of this section is to illustrate the use of the R BMTME package for analyzing multi-environment and multi-trait and multi-environment data from plant breeding programs. The BMTME package was built following the paper by Montesinos-Lopez et al. (2016) and implemented in the R statistical software (R Core Team 2018).
The development version of the BMTME package can be installed directly from the GitHub repository (https://github.com/frahik/ BMTME). In order to install the package, it is necessary to install the appropriate compilers; the installation process and the required tools depend heavily on the operating system. For example, in Windows it is necessary to install Rtools (https://cran.r-project.org/bin/windows/ Rtools/), and in modern versions of macOS, it is necessary to install XCode from App Store or the development tools for R from CRAN (https://cran.r-project.org/bin/macosx/tools/). In the case of Linux, it is necessary to install the C++ compilers included in your distribution, for example, g++ from GNU (https://www.gnu.org). Once the tools have been installed, use the following command to install the package within your R session: install.packages(9devtools9) devtools::install_github(9frahik/BMTME9) You can also find the package in the CRAN repository, and you can use the following command (see below) to install a version of the package from CRAN. This will avoid the need to install some dependencies manually and install the Rtools software using the following command: install.packages(9BMTME9) The R package BMTME is available at the following link: https:// cran.r-project.org/web/packages/BMTME/index.html.

RESULTS
The results are given in three main sections. The first section illustrates the use of the BME function for implementing multi-environment analysis, while the second and the third sections illustrate the use of the BMTME and BMORs functions for implementing multi-trait and multi-environment analyses.

The BME Function
This example illustrates how to fit a model when there is only one environment and several dependent variables. First, we load the library: library(BMTME) Then we load the Mada dataset: data("WheatMadaToy") Then we define the model to be adjusted; since the dataset only includes an environment where several dependent variables were evaluated, the BME model is used. To implement it, first we need to order the dataset as follows: phenoMada ,-(phenoMada[order(phenoMada$GID),]) rownames(phenoMada)=1:nrow(phenoMada) head(phenoMada) GID PH FL FE NS SY NP 1 9 29.7776 -8.8882 -4.93900 1.04100 169.06 28.8025 2 11 3. 2210 -7.1111 -0.36940 -3.88940 -107.19 58.2516 3 12 6.1670 -9.5337 -12.43680 2.58250 -160.54  This is a very important step in the analysis, because if the dataset is not ordered correctly, this may cause conflicts and produce incorrect estimations. Also, with the head() function we printed the phenotypic dataset, where the required format of the dataset requires a first column with the identifiers of the lines and then the names of all the traits. It is important to respect this format to be able to successfully implement the multi-environment (trait) datasets.
Then, the design matrix for the genetic effects should be generated, as shown below.
LG ,cholesky(genoMada) ZG ,model.matrix($0 + as.factor(phenoMada$GID)) Z.G ,-ZG % Ã % LG Then, we can extract the phenotypic responses that were converted to matrix object as shown in the following command: Y ,as.matrix(phenoMada[, -c(1)]) Finally, the model was adjusted, and 30,000 iterations were used to adjust the model. It is important to point out that bs is the block size for sampling from posterior distributions; we suggest using a value of at least 50 but less than 1000.
Next we used the names() function to identify all the available outputs of the fitted model.

names(fm)
[1] "Y" "nIter" "burnIn" "thin" "dfe" [6] "Se" "yHat" "SD.yHat" "beta" "SD.beta" [11] "b1" "SD.b1" "vare" "SD.vare" "varTrait" [16] "SD.varTrait" "NAvalues" Here we extracted the observed values ($Y), the predicted values ($yHat), the parameters provided for the model fit ($nIter, $burnIn, $thin, etc.), the estimates of the beta coefficients, random effects of lines and the genetic and residual covariances ($beta, $SD.beta, $b1, $SD.b1, $varTrait, $vare, etc. For demonstration purposes, we will only extract the first 6 predictions for the 6 evaluated traits. We also plotted the observed values against the predicted values for each trait, as follows (see Figure 1): plot(fm, trait = 9FL9) Since the code provided above is only appropriate for parameter estimation, now we provide the code required to evaluate the prediction accuracy using the BME() function. For this reason, first we built the random CV strategy with 10 random partitions, each with TRN = 80% and TST = 20%, using the following code: Here we see that the best prediction in terms of APC was found in trait PH (0.5612), while the worst was in trait SY (0.0242). However, in terms of MAAPE, the best prediction was observed in trait NP (0.7070), while the worst was found in trait FL (0.7751). With the boxplot(pm) function, we created a plot summarizing the predictions in terms of Pearson's correlation, but if users want this plot in MAAPE terms, they need to use the following code: boxplot(pm, select="MAAPE") ( Figure 2).

boxplot(pm)
It is important to point out that the BME function can be used with only 1 testing set that can be defined by the user, as shown in the following example: CrossV1 ,-sample(nrow(Y),15) pm ,-BME(Y = Y, Z1 = Z.G, nIter = 1250, burnIn = 500, thin = 2, bs = 50,testingSet = CrossV1) Next we summarize the prediction accuracy: Since only one training set and one testing set were used, the standard errors for both metrics appear with NaN, given that it is not possible to calculate the standard error because only one testing set is available.
The Bayesian Multi-Trait and Multi-Environment (BMTME) function This example illustrates how to fit a model with multiple traits and multiple environments. To do this, use the Maize dataset; first, load the data using the following function: data("MaizeToy") Next, order the dataset, rename the rows of the phenotypic dataset and print the first six observations of the data in order to see the structure required of the data, which consists of a first column that includes the lines, a second column that includes the environments and third, fourth and fifth columns that correspond to the traits under study.
We used the names() function to see all the things that can be extracted after fitting a model with the BMTME function.
names(fm) [1] "Y" "nIter" "burnIn" "thin" "dfe" "Se" "yHat" [8] "SD.yHat" "beta" "SD.beta" "b1" "b2" "vare" "SD.vare" [15] "varEnv" "SD.varEnv" "varTrait" "SD.varTrait" "NAvalues" We can extract the predicted and observed values, the random effects of lines, of lines·trait, lines·environment·trait, as well as the genetic covariances between traits and environments, and the residual covariance between traits. To extract the matrix of covariances between traits, we used the following: To convert this covariance matrix between traits into a correlation matrix, we used the following command: The observed and predicted values of each trait can be plotted using the following plot(), function ( Figure 3): We provide the corresponding code for the fivefold CV strategy and its implementation with the BMTME function: The summary information is given with the following code using the boxplot(pm, select="MAAPE", las = 2), where we added the parameter las = 2, to show the labels vertically and be able to distinguish the complete names (Figure 4). boxplot(pm, select="MAAPE", las = 2) The BMORS Function Since the BMORS function is only appropriate for evaluating prediction performance (but not for parameter estimation), we provide the required R script for evaluating the prediction performance of multiple trait and multiple environment data with a random CV strategy. Since we use the same dataset that was used to illustrate the BMTME function, we only provide the R code for building the predictor (ETA), the CV strategy and for implementing the BMORS function. To build the linear predictor, we used the following R code: ETA ,list(Env = list(X = Z.E, model = "BRR"), Gen = list(X = Z.G, model = "BRR"), EnvGen = list(X = Z.EG, model = "BRR")) Next, we provide the R code for implementing the random CV strategy: CrossValidation ,-CV.RandomPart(pheno, NPartitions = 10, PTesting = 0.2, set_seed = 123) Finally, the model will be implemented with 15,000 iterations, of which 10,000 will be used as burn-in to fit the model. Below we show the resulting predictive performance using the summary() function. The information = 9complete9 command inside the summary() function shows the prediction performance (in terms of Pearson's correlation and MAAPE) for all partitions implemented, but we only show the first 20 observations.  To create a graph with a summary of the predictions in terms of Pearson's correlation and in terms of MAAPE, we used the plot() function ( Figure 5). Because the names are composed of the evaluated traits and environments, we added the parameter las = 2 to show the labels in a vertical way and to distinguish the complete names of the trait-environment combinations. In addition, we used the par() function and the mar parameter to modify the margins of the graph. par(mar = c(6,4,2,1)) plot(pm, las = 2) Figure 5 shows that the lowest average Pearson's correlation obtained was observed in the ASI_KAK and ASI_KTI trait-environment combinations, while the highest average Pearson's correlation was obtained in the PH_KAK trait-environment combination. It is possible to create a boxplot with the results of the MAAPE, using the following command ( Figure 6): boxplot(pm, select ="MAAPE", las = 2) Figure 6 shows that the lowest MAAPE was for PH_KAK (best prediction), while the highest MAAPE was for the ASI_KAK traitenvironment combination (worst prediction).

The BMORS_ENV Function
This function is useful for predicting whole environments using the remaining environments as training. Next we provide the R code for evaluating the prediction performance of the same maize dataset, but using the KAK environment as training and the KTI and EBU environments as testing. Two important things to point out for using this function are: (a) that we provided not only the matrix of response variables, but also a data.frame that contains, in the first column, the names of the environments followed by information on all response variables, and (b) we did not create a separate file for specifying the training and testing individuals; we only specified in testingEnv which environments are used as testing; the remaining environments are used by default as training, as shown below.
dataset ,-phenoMaizeToy[, 2:5] pm ,-BMORS_Env(dataset, testingEnv = c(9KTI9, 9EBU9), ETA = ETA, covModel = 9BayesB9, nIter = 15000, burnIn = 10000, thin = 2, progressBar = TRUE, digits = 3) summary ( For this example, we specified that covModel = 9BayesB9, which means that the Bayesian BayesB model will be implemented for the second stage of the model where it is implemented (equation 3). In covModel, in addition to Bayesian Ridge regression (BRR) and BayesB, we can also implement BayesA, BayesC and Bayesian Lasso (BL); however, the BRR model is implemented by default.To create a graph with Pearson's correlation or the MAAPE index, we used the barplot() function, as shown below (Figure 7): barplot(pm) Figure 7 shows that the lowest Pearson's correlation obtained was in the ASI_EBU trait-environment combination, while the highest Pearson's correlation was obtained in the Yield_KTI trait-environment combination.

DISCUSSION
As mentioned in the introduction, we propose a Bayesian R package for implementing multi-environment and multi-trait and multi-environment analysis for parameter estimation and for evaluating prediction accuracy. We illustrate the four main functions [BME(), BMTME(), BMORS() and BMORS_Env()] of the BMTME package with real toy datasets starting from the type and preprocessing required to make correct use of each of these datasets for parameter estimation and for evaluating prediction performance. It is important to point out that one advantage of the BME and BMTME functions is that, in addition to being used to evaluate the prediction accuracy, they can also be used for parameter estimation, which allows estimating the random effects (lines, lines·environments for each trait) and variance-covariance matrices of genetic (for traits and environments) and residual (for traits) effects. The BMORS() and BMORS_Env() functions are not useful for obtaining parameter estimates of covariances between traits and environments because they implement univariate analysis at both stages. However, they have two important advantages: (a) they allow implementing even more complex predictors than the one specified in equation (1), which modifies the ETA list used to create the predictor, and (b) the computational resources required to implement it are much less than those needed by the BMTME() function for implementing multitrait and multi-environment data. This last point is observed in Figure 8 where the implementation time for the Mada and Maize datasets is reported. The figure shows that in the Mada dataset, the BMORS model was more than 15 times faster than the BMTME model (25.246/ 1.621= 15.572), while in the Maize dataset, the BMORS model was more than 37 times faster than the BMTME model (25.668/ 0.692= 37.093); these results were obtained with 10000 iterations of the Gibbs sampler. The BMTME R package provides very synthetic summaries (tables and plots) of the prediction accuracies, which are ready to be interpreted  and used to write the results in a manuscript. Additionally, we provide three types of cross-validations useful for breeders that are implemented in this R package, which is very simple to use and implement.
The main disadvantage of the BME() and BMTME() functions of the BMTME R package is that the computational resources required for their implementation are very demanding; fortunately, the parameter estimates involved are stabilized very quickly even with few iterations.
The toy examples used in this article are for illustration purposes and to help users follow, as easily as possible, the necessary steps for running the different processes. Comparing them with other software of similar type is not possible, as no similar software for simultaneously fitting multitrait multi-environment is currently available. For example, the MTM (de los Campos and Grüneberg 2016) is an efficient Bayesian multi-trait software but is not multi-trait and multi-environment. Future research on benchmarking the BME() and BMTME() functions of the BMTME package with other potential software to be developed in terms of computing time for processing large datasets should be performed. However, the BMORS() and BMORS_Env() functions that also belong to the BMTME R package are very efficient in terms of computational resources, which gives the user an alternative option for performing this type of analyses.
It is important to point out that the proposed BMTME package is different from existing multi-trait analysis software such as ASREML (Gilmour et al., 1995), sommer (Covarrubias-Pazaran 2016), BGGE (Granato et al., 2018) and MCMCglmm (Hadfield et al., 2010). In addition to taking into account variance-covariance matrices of traits (genetic and residual), it also takes into account the genetic covariance (correlation) between environments, which is estimated from the data. This can help improve parameter estimates and prediction accuracy when the degree of correlation between traits is moderate or high.
Multi-trait models are preferred over univariate-trait models because they have the following advantages: (a) they produce higher prediction accuracy because they have more information (direct or indirect) and better data connectedness (Colleau et al., 2009); (b) they improve index selection because optimal weight factors can be obtained for the total merit index (Colleau et al., 2009); and (c) they allow obtaining more precise genetic and residual covariances and incorporating them into expected breeding value (EBV) estimates for across-location, acrosscountry or across-region evaluations (Thompson and Meyer 1986;Schaeffer 2001).
Note that the two datasets used for illustrating the main functions of the BMTME R package are datasets with few lines (toy datasets) with the main intention that users interested in using the package can obtain results very quickly and practice using the software. However, the structure of the data are exactly the same as the structure of the data produced in plant breeding programs. The two toy datasets are included in the BMTME package to facilitate its implementation and allow users to practice using the R software.
To conclude, this paper presents the R package BMTME which allows the implementation of multi-trait, multi-trait and multienvironment analysis for estimating parameters (genetic correlation between traits and environments, residual correlation between traits, random effects of lines and lines·environments) and evaluating the prediction accuracies of many traits simultaneously. We illustrate the implementation of the main functions (BME, BMTE and BMORS) of the R package with two toy real datasets that are very common in plant breeding programs. We provide details of the characteristics that each of the datasets must have, and show how to build the CV strategies available in the package, how to prepare the data to implement the main functions of the BMTME package, how to extract the parameter estimates and how to obtain the summary and plots of prediction accuracies resulting from the implemented CV strategy. The computing time of the BME() and BMTE() functions of the BMTME R package for large datasets is significantly more demanding (in terms of time) than for the toy examples used in this study.