A posterior sample answers the questions a breeder asks in probabilities: how likely is a genotype to be among the best, how stable is it, and what yield can be expected in a cell that was never grown. This vignette goes through bam_superiority(), bam_ranks(), bam_stability(), bam_predict() and bam_winners(), and through the model choice tools bam_compare() and bam_ppc().
2 Probabilities of superior performance and stability
Dias et al. (2022) select the best \(s\) fraction of \(G\) genotypes, that is \(\lceil sG \rceil\) of them, in every draw. Write \(\bar{\eta}_{i\cdot}\) for the mean of genotype \(i\) over the environments and \(\sigma^2_{GE,i} = \frac{1}{E}\sum_j(\eta_{ij} - \bar{\eta}_{i\cdot} - \bar{\eta}_{\cdot j} +
\bar{\eta}_{\cdot\cdot})^2\) for its interaction variance. Then
the probability of superior performance is the share of draws in which \(\bar{\eta}_{i\cdot}\) is among the best \(\lceil sG \rceil\);
the probability of superior stability is the share in which \(\sigma^2_{GE,i}\) is among the smallest \(\lceil sG \rceil\);
the joint probability is their product;
the conditional probability is the probability of superior performance within each environment;
the pairwise probability is \(P(\bar{\eta}_{i\cdot} > \bar{\eta}_{i'\cdot})\).
The tile of genotype \(i\) and rank \(r\) shows \(P(\text{rank}_i = r)\). A rank interval that spans several places is the honest answer to “which is second?”. Ranks within one environment come from .env:
bam_ranks(fit = fit , higher_better =TRUE , .env ="12" )$summary[1:3]
bam_stability() gives the posterior of two stability measures for every genotype: the interaction variance of Shukla (1972), and the AMMI stability value of Purchase et al. (2000), \[\mathrm{ASV}_i = \sqrt{\left(\frac{\lambda_1}{\lambda_2}\,\lambda_1^{1/2}\alpha_{i1}\right)^2 +
(\lambda_2^{1/2}\alpha_{i2})^2}.\] The weight \(\lambda_1/\lambda_2\) becomes unstable when \(\lambda_2\) is drawn near zero, so it is computed from the posterior means of \(\lambda_1\) and \(\lambda_2\), and the scores from each draw. The yield stability index (YSI) adds the mean rank of the yield to the mean rank of the ASV; small is good.
st <-bam_stability(fit = fit , higher_better =TRUE )st
bayes_ammi() stores the deviance information criterion (DIC; Spiegelhalter et al. 2002) and the widely applicable information criterion (WAIC; Watanabe 2010) of the observed plots. Both are computed while the chains run, so no draws need to be kept. bam_compare() puts several fits side by side:
Comparison of Bayesian models (smaller is better)
model terms variances DIC pD WAIC p_waic se
<char> <int> <char> <num> <num> <num> <num> <num>
1: K2 2 per environment 11736.67 86.5514 11743.04 83.1810 42.0135
2: K3 3 per environment 11736.29 97.0824 11744.16 93.1994 41.9472
3: K1 1 per environment 11751.68 68.9559 11756.55 66.9555 41.7350
lambda_last_lower lambda_last_mean dWAIC dWAIC_se
<num> <num> <num> <num>
1: 738.9639 2101.743 0.0000 0.0000
2: 3.9045 1135.815 1.1235 2.9675
3: 4233.4691 5259.183 13.5094 6.5130
Sources: Spiegelhalter et al. (2002); Watanabe (2010); Vehtari et al. (2017).
Full references: bam_cite(x).
plot(x = cp)
WAIC estimates the deviance of new plots of the same genotypes and environments (Gelman et al. 2014). se is the standard error of each WAIC; dWAIC_se is the standard error of the difference from the best model, computed from the plot-by-plot differences, and is the one to use when comparing (Vehtari et al. 2017). Here K2 has the smallest WAIC; K3 is 1.1 above it, 0.4 standard errors of the difference, and K1 is 13.5 above (2.1 standard errors). The column lambda_last_lower adds the posterior view: the HPD interval of the last singular value of K2 reaches down to 740, so a term whose interval nearly reaches zero improves prediction little and is weakly supported. Prefer the simplest model among those that the standard errors do not separate.
7.2 Posterior predictive checks
bam_ppc() simulates replicated trials from the posterior and compares statistics of the data with their replicated distribution (Gelman et al. 1996). A posterior predictive p-value near 0 or 1 marks a feature of the data the model does not reproduce.
bayesplot::mcmc_areas(x = fit$draws, pars =c("tau[4]", "tau[5]", "tau[6]"), prob =0.95)
9 Interpretation
bam_interpret(x = sp , higher_better =TRUE , top =5 )
Probabilities of superior performance and stability
Selecting 2 genotype(s): the highest joint probability of superior
performance and stability is that of 6 (0.13).
Findings:
- Probability of superior performance, largest first: 6 0.89, 5 0.63, 4 0.47;
under 0.005: 9, 3, 1, 2, 7, 8.
- Probability of superior stability, largest first: 7 0.72, 3 0.72, 2 0.24, 6
0.15, 5 0.10, 9 0.06, 1 0.02; under 0.005: 8, 4.
- Most likely superior in performance: 6; in stability: 7.
Recommendations:
- For a wide recommendation consider 6 and 5 (joint probability of at least
0.05); the probabilities by environment (component conditional) show where
each genotype is superior.
Cautions:
- The joint probability multiplies the two probabilities as in Dias et al.
(2022), which treats them as independent.
Sources:
- Dias, K. O. G., dos Santos, J. P. R., Krause, M. D., Piepho, H.-P.,
Guimaraes, L. J. M., Pastina, M. M. and Garcia, A. A. F. (2022). Leveraging
probability concepts for cultivar recommendation in multi-environment
trials. Theoretical and Applied Genetics. doi:10.1007/s00122-022-04041-y
bam_interpret(x = st , higher_better =TRUE , top =5 )
Bayesian stability
Best yield stability index (expected ranks of mean and ASV): 6, 5, 3, 7 and
9; most probably the most stable: 7 (0.56).
Findings:
- ASV with HPD intervals: 6 214 (99.8 to 329); 5 226 (94.4 to 339); 3 130
(8.34 to 231); 7 103 (3.48 to 205); 9 253 (133 to 367); 4 406 (297 to 512);
2 197 (77.4 to 314); 1 303 (179 to 421); 8 375 (245 to 490).
Recommendations:
- Consider 6, 5 and 3: they combine a good mean with a small AMMI stability
value.
Cautions:
- The ASV uses the first two terms; with one term it is the absolute first
score.
Sources:
- Purchase, J. L., Hatting, H. and van Deventer, C. S. (2000). Genotype x
environment interaction of winter wheat (Triticum aestivum L.) in South
Africa: II. Stability analysis of yield performance. South African Journal
of Plant and Soil, 17(3), 101-107.
- Shukla, G. K. (1972). Some statistical aspects of partitioning
genotype-environmental components of variability. Heredity, 29(2), 237-245.
doi:10.1038/hdy.1972.87
- Farshadfar, E. (2008). Incorporation of AMMI stability value and grain
yield in a single non-parametric index (GSI) in bread wheat. Pakistan
Journal of Biological Sciences, 11(14), 1791-1796.
10 References
Dias, K. O. G., dos Santos, J. P. R., Krause, M. D., Piepho, H.-P., Guimaraes, L. J. M., Pastina, M. M. and Garcia, A. A. F. (2022). Leveraging probability concepts for cultivar recommendation in multi-environment trials. Theoretical and Applied Genetics. doi:10.1007/s00122-022-04041-y
Shukla, G. K. (1972). Some statistical aspects of partitioning genotype-environmental components of variability. Heredity, 29(2), 237–245.
Purchase, J. L., Hatting, H. and van Deventer, C. S. (2000). Genotype x environment interaction of winter wheat in South Africa: II. Stability analysis of yield performance. South African Journal of Plant and Soil, 17(3), 101–107.
Spiegelhalter, D. J., Best, N. G., Carlin, B. P. and van der Linde, A. (2002). Bayesian measures of model complexity and fit. Journal of the Royal Statistical Society B, 64(4), 583–639.
Watanabe, S. (2010). Asymptotic equivalence of Bayes cross validation and widely applicable information criterion in singular learning theory. Journal of Machine Learning Research, 11, 3571–3594.
Gelman, A., Hwang, J. and Vehtari, A. (2014). Understanding predictive information criteria for Bayesian models. Statistics and Computing, 24(6), 997–1016.
Vehtari, A., Gelman, A. and Gabry, J. (2017). Practical Bayesian model evaluation using leave-one-out cross-validation and WAIC. Statistics and Computing, 27(5), 1413–1432.
Gelman, A., Meng, X.-L. and Stern, H. (1996). Posterior predictive assessment of model fitness via realized discrepancies. Statistica Sinica, 6(4), 733–807.
Gabry, J., Simpson, D., Vehtari, A., Betancourt, M. and Gelman, A. (2019). Visualization in Bayesian workflow. Journal of the Royal Statistical Society A, 182(2), 389–402.
Vehtari, A., Gelman, A., Simpson, D., Carpenter, B. and Burkner, P.-C. (2021). Rank-normalization, folding, and localization: an improved R-hat for assessing convergence of MCMC. Bayesian Analysis, 16(2), 667–718.
Crossa, J., Perez-Elizalde, S., Jarquin, D., Cotes, J. M., Viele, K., Liu, G. and Cornelius, P. L. (2011). Bayesian estimation of the additive main effects and multiplicative interaction model. Crop Science, 51(4), 1458–1469. doi:10.2135/cropsci2010.06.0343
