Table of contents
\[
\newcommand{\vy}{\mathbf{y}}
\newcommand{\vx}{\mathbf{x}}
\newcommand{\vb}{\mathbf{b}}
\newcommand{\ve}{\mathbf{e}}
\newcommand{\vone}{\mathbf{1}}
\newcommand{\vbeta}{\boldsymbol{\beta}}
\newcommand{\vmu}{\boldsymbol{\mu}}
\newcommand{\mX}{\mathbf{X}}
\newcommand{\mZ}{\mathbf{Z}}
\newcommand{\mV}{\mathbf{V}}
\newcommand{\mI}{\mathbf{I}}
\newcommand{\mJ}{\mathbf{J}}
\newcommand{\mSigma}{\boldsymbol{\Sigma}}
\newcommand{\mY}{\mathbf{Y}}
\newcommand{\mK}{\mathbf{K}}
\newcommand{\mH}{\mathbf{H}}
\newcommand{\va}{\boldsymbol{\alpha}}
\newcommand{\vg}{\boldsymbol{\gamma}}
\newcommand{\vt}{\boldsymbol{\tau}}
\newcommand{\vd}{\boldsymbol{\delta}}
\newcommand{\vtheta}{\boldsymbol{\theta}}
\newcommand{\vv}{\mathbf{v}}
\newcommand{\E}{\mathrm{E}}
\newcommand{\Var}{\mathrm{Var}}
\newcommand{\Cov}{\mathrm{Cov}}
\newcommand{\tr}{\operatorname{tr}}
\]
Abstract
The sites regression (SREG) model puts the genotype main effects into the bilinear terms, so its biplot shows which genotypes do best where. This vignette fits it with the first year as prior, then shows how the weight of the earlier trial changes the posterior and how WAIC chooses the number of terms.
The SREG model with the first year as prior
<- bayes_ammi (.data = cultivo2009.y = y.gen = entry.env = site.rep = rep.nIter = 3000 .prior_data = cultivo2008.n0 = NULL .model = "SREG" .k = 2 .burnin = 1000 .thin = 1 .chains = 2 .priors = bys_priors (mu_mean = 0 mu_var = 1e+8 tau_mean = 0 tau_var = 1e+8 delta_mean = 0 delta_var = 1e+8 lambda_mean = 0 lambda_var = 1e+8 nu = 1 s2 = NULL .level = 0.95 .seed = 2012 .align = TRUE .waic = TRUE
Bayesian SREG with 2 multiplicative term(s), an earlier trial as prior (n0 = 3): 12 genotypes, 25 environments, 3 replicates
2 chain(s) x 3000 kept draws (iterations 3000 after a burn-in of 1000, thinning 1); 2.0 s
Posterior means, standard deviations and 95% HPD intervals
parameter mean sd lower upper ess rhat
<char> <num> <num> <num> <num> <num> <num>
1: mu 3848.46 30.35 3790.6 3907.58 6000 1.0002
2: sigma2 1686438.64 101459.21 1483150.1 1881754.16 5552 1.0007
3: delta[1] 128.59 146.43 -142.6 426.81 6000 1.0002
4: delta[2] -533.71 147.60 -829.2 -256.99 6069 1.0000
5: delta[3] 866.34 149.15 564.0 1154.05 6203 1.0003
6: delta[4] -380.47 147.99 -666.8 -90.91 5551 0.9999
7: delta[5] 1535.74 149.97 1239.4 1824.21 6219 1.0000
8: delta[6] -95.27 149.14 -376.5 196.74 5845 0.9999
9: delta[7] -1758.38 147.83 -2037.1 -1461.39 6313 0.9998
10: delta[8] 350.94 150.17 64.3 650.00 6000 1.0011
11: delta[9] -1471.15 147.61 -1756.5 -1180.45 6000 1.0011
12: delta[10] -745.08 147.81 -1029.4 -447.52 6167 0.9999
13: delta[11] -1203.17 149.42 -1494.8 -912.56 5763 1.0002
14: delta[12] -1897.55 148.98 -2192.3 -1610.74 6000 0.9999
15: delta[13] -969.30 149.85 -1261.3 -688.02 6000 0.9999
16: delta[14] -550.54 149.61 -847.5 -266.79 5646 0.9999
17: delta[15] -2214.80 148.81 -2497.1 -1914.30 6000 1.0005
18: delta[16] 91.07 151.50 -209.2 381.44 6000 0.9999
19: delta[17] 161.26 150.95 -131.9 459.36 6404 1.0002
20: delta[18] 1845.05 148.78 1561.0 2134.90 6000 1.0001
21: delta[19] -867.08 149.69 -1147.8 -567.86 6000 1.0003
22: delta[20] -1111.97 150.32 -1387.3 -799.39 6000 1.0000
23: delta[21] -1625.33 148.66 -1916.0 -1333.87 5454 1.0005
24: delta[22] 4918.80 148.58 4633.1 5218.63 5814 1.0001
25: delta[23] 2864.00 151.26 2561.8 3152.60 5649 1.0003
26: delta[24] 910.64 152.01 609.9 1205.86 6234 1.0002
27: delta[25] 1751.38 154.86 1432.6 2035.94 6000 1.0001
28: lambda[1] 6250.53 579.32 5139.4 7388.08 3888 1.0004
29: lambda[2] 3268.10 577.56 2222.5 4508.37 4130 0.9999
parameter mean sd lower upper ess rhat
<char> <num> <num> <num> <num> <num> <num>
Scores of the bilinear terms: x$summary (blocks alpha, gamma); biplot coordinates: x$scores
Sources: Perez-Elizalde et al. (2012); Crossa et al. (2011); Yan et al. (2000);
Hoff (2009); Spiegelhalter et al. (2002); Watanabe (2010). Full references:
bys_cite(x).
<- bayes_gge_biplot (fit = fsburnin = 0 thin = 1 pb = 0.05 plot_stable = FALSE plot_unstable = TRUE ncolors = 5 .which = "both" .grid = 100 $ biplot
On an SREG biplot the first axis reflects mostly the genotype means, so genotypes far out on it are above or below average rather than interacting.
The weight of the earlier trial
.n0 sets how many observations per cell the earlier trial counts for. By default it is its number of replicates (3).
<- c (0.5 , 3 , 12 )<- :: rbindlist (l = lapply (X = wFUN = function (n0) {<- bayes_ammi (.data = cultivo2009.y = y.gen = entry.env = site.rep = rep.nIter = 1500 .prior_data = cultivo2008.n0 = n0.model = "AMMI" .k = 2 .burnin = 500 .thin = 1 .chains = 1 .priors = bys_priors (mu_mean = 0 mu_var = 1e+8 tau_mean = 0 tau_var = 1e+8 delta_mean = 0 delta_var = 1e+8 lambda_mean = 0 lambda_var = 1e+8 nu = 1 s2 = NULL .level = 0.95 .seed = 4 .align = TRUE .waic = TRUE :: data.table (n0 = n0mu = f$ summary[parameter == "mu" , mean]alpha1 = f$ summary[block == "tau" , mean][1 ]lambda1 = f$ summary[parameter == "lambda[1]" , mean]
n0 mu alpha1 lambda1
<num> <num> <num> <num>
1: 0.5 3690.593 484.3266 3602.049
2: 3.0 3849.168 136.0094 3578.433
3: 12.0 3980.321 -154.0709 4484.151
As \(n_0\) grows the posterior moves from the second year (grand mean 3627) towards the first (4069). With \(n_0 = 3\) it sits halfway. A small \(n_0\) is the choice when the earlier trial was grown under different conditions.
How many terms?
WAIC compares fits with different numbers of terms (Watanabe 2010; Vehtari et al. 2017):
<- lapply (X = c (1 , 2 , 3 , 5 )FUN = function (K) {bayes_ammi (.data = cultivo2009.y = y.gen = entry.env = site.rep = rep.nIter = 1500 .prior_data = cultivo2008.n0 = NULL .model = "AMMI" .k = K.burnin = 500 .thin = 1 .chains = 1 .priors = bys_priors (mu_mean = 0 mu_var = 1e+8 tau_mean = 0 tau_var = 1e+8 delta_mean = 0 delta_var = 1e+8 lambda_mean = 0 lambda_var = 1e+8 nu = 1 s2 = NULL .level = 0.95 .seed = 5 .align = TRUE .waic = TRUE names (fits) <- paste0 ("K" , c (1 , 2 , 3 , 5 ))<- do.call (what = bys_compare, args = fits)
Comparison of Bayesian models (smaller is better)
model terms prior DIC pD WAIC p_waic
<char> <int> <char> <num> <num> <num> <num>
1: K1 1 earlier trial (n0 = 3) 15489.01 26.7424 15493.10 29.7240
2: K2 2 earlier trial (n0 = 3) 15493.61 37.4422 15499.51 41.2472
3: K3 3 earlier trial (n0 = 3) 15504.81 44.3373 15511.25 48.0848
4: K5 5 earlier trial (n0 = 3) 15518.21 56.6405 15524.42 58.8980
se lambda_last_lower lambda_last_mean dWAIC dWAIC_se
<num> <num> <num> <num> <num>
1: 50.1988 1863.926 3023.724 0.0000 0.0000
2: 51.5078 2249.549 3069.854 6.4048 7.3493
3: 51.6190 1456.046 2366.735 18.1503 8.8759
4: 52.0605 1093.993 1900.412 31.3123 10.7548
Sources: Spiegelhalter et al. (2002); Watanabe (2010); Vehtari et al. (2017).
Full references: bys_cite(x).
WAIC prefers K1; K2 lie within two standard errors of the difference (dWAIC_se), so the data do not separate them. With the first year as a strong prior, the singular values of every term are pulled towards that year’s, and few terms suffice for prediction. The paper chose three terms with the Bayes factor, which compares the marginal likelihoods of the models rather than their predictive accuracy.
References
Perez-Elizalde, S., Jarquin, D. and Crossa, J. (2012). A general Bayesian estimation method of linear-bilinear models applied to plant breeding trials with genotype x environment interaction. Journal of Agricultural, Biological, and Environmental Statistics , 17, 15–37. doi:10.1007/s13253-011-0063-9
Yan, W., Hunt, L. A., Sheng, Q. and Szlavnics, Z. (2000). Cultivar evaluation and mega-environment investigation based on the GGE biplot. Crop Science , 40(3), 597–605.
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.
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.
