Fit the Bayesian Multi-Environment model and apply a random partition cross-validation

This example illustrates how to fit a model where there is only one environment and several traits, even though it may contradict the model. To do this, the WheatMadaToy data set is used. To load the data set, the data() function as shown below:

Then we proceed to define the model to be adjusted, as the data set only includes an environment where several traits were evaluated, the BME model is used, and for its implementation, first we need order the dataset as follows:

This is the most important step in the analysis, because if the dataset is not ordered by the identifiers, may cause conflicts between the analysis producing incorrect estimations. Then, the matrix design for the genetic effects should be generated, as shown below,

Then, we extract the phenotypic responses and it’s converted to matrix object as shown in the following command,

Finally, the model is adjusted, and for demonstration purposes only 1250 iterations are used to adjust the model.

GID	PH	FL	FE	NS	SY	NP
9	29.7776	-8.8882	-4.93900	1.04100	169.06	28.8025
11	3.2210	-7.1111	-0.36940	-3.88940	-107.19	58.2516
12	6.1670	-9.5337	-12.43680	2.58250	-160.54	17.1278
15	6.8117	4.6377	11.78860	-0.03378	235.70	-19.6571
20	-14.4480	3.2525	6.40780	-14.23460	131.87	42.2962
21	-13.2185	3.8902	0.09722	5.35680	164.06	36.8239

To know all the details about the output of the fitted model, we use the str() function, which returns the following information:

Where we can observe that the returned data corresponds to the observed values ($Y), the parameters provided for the model fit ($nIter, $burnIn, $thin, etc.), the predicted values ($yHat) and the estimations of the beta coefficients, random effects of lines and the genetic and residual covariances ($beta, $SD.beta, $b1, $SD.b1, $varTrait, $vare, etc.). Since breeders are interested in understanding the matrix of covariance and genetic correlation, these can be extracted with the following command,

To convert this covariance matrix into a correlation matrix it is suggested to use the following command,

Where we can observe that there is not a high correlation (greater than 0.5). Below is a sample of how we can obtain the residual covariance (correlation) matrix,

To convert the residual matrix into the correlation matrix it is suggested to use the following command,

Where we can observe that the residual of traits is not highly correlated (greater than 0.5). On the other hand, to extract the predicted values from the model it is necessary call the $yhat values of the object within the adjusted model, for demonstration purposes we will only extract the first 6 predictions for the 6 traits evaluated.

The software also allows plotting the observed values against the predicted values by trait, as shown below:

Cross validation with random partitions

On the other hand, the package also allows making cross-validation for the predictions for this we require a data.frame object with the phenotypes as illustrated below:

pheno <- data.frame(GID = phenoMada[, 1], Env = '', Response = phenoMada[, 3])

Once the object is generated, we will use the CV.RandomPart() function, to form the training and testing sets of each random partition of the cross validation, it is suggested to provide the number of random partitions, the percentage of the data to be used for the testing and a seed to guarantee a reproducible analysis.

CrossV <- CV.RandomPart(pheno, NPartitions = 4, PTesting = 0.2, set_seed = 123)

Since the partitions have been generated for cross validation, we will use the BME() function to fit the model, however, we will include the above object in the testingSet parameter to implement cross validation. In addition, we will use the summary() function, to show the resulting predictive capacity by trait evaluated.

pm <- BME(Y = Y, Z1 = Z.G, nIter = 1000, burnIn = 500, thin = 2, bs = 50,
          testingSet = CrossV, progressBar = FALSE)
summary(pm)
#>   Environment Trait Pearson SE_Pearson  MAAPE SE_MAAPE
#> 1                FE  0.4715     0.2096 0.7472   0.0714
#> 2                FL  0.1293     0.1624 0.7996   0.0739
#> 3                NP  0.4975     0.2827 0.7033   0.0334
#> 4                NS  0.6522     0.0950 0.6664   0.0382
#> 5                PH  0.5811     0.1315 0.7754   0.0757
#> 6                SY  0.1956     0.1446 0.7029   0.0613

We can observe that the results are given around the traits that have been evaluated, showing for each of them the predictive capacity obtained using the Pearson correlation average of all the partitions used in the cross validation; as well as the mean arctangent percentage error (MAAPE). From the results obtained, we can emphasize that the traits that have obtained a better predictive capacity were the NS and PH traits. While the SY and FL traits have low predictive capabilities. In addition, the package offers the ability to generate graphs of the results obtained, for it we will use the function that appears next to show a boxplot with the results obtained.

boxplot(pm)

In this figure we can observe that the PH trait has obtained the best predictive capacity on average, while the SY trait has obtained the lowest predictive capacity. It is also possible to graph the predictive capacity of the traits under study using the MAAPE, just specify it through the select parameter, as shown below,

boxplot(pm, select = 'MAAPE')

In this figure we can observe that the FL trait has been the one that has obtained the highest MAAPE index and therefore presents the worst prediction, but it is like the rest of them, so we cannot trust in the predictions of the analysis.

Fit the Bayesian Multi-Environment model and apply a random partition cross-validation

2019-02-02

Cross validation with random partitions