fit <-
bayes_ammi(
.data = Maiz
, .y = y
, .gen = entry
, .env = site
, .rep = rep
, .nIter = 3000
, .k = 2
, .burnin = 500
, .thin = 1
, .chains = 2
, .priors = bam_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 = 11
, .align = TRUE
, .model = "AMMI"
, .het = FALSE
, .waic = TRUE
)Posterior distributions of stability indices
2026-11-05
Source:vignettes/bayesammi_11_stability_indices.qmd
\[ \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}} \]
1 Abstract
Dozens of stability indices are in use, each computed once from the cell means of a trial. Computed on every posterior draw of the cell means instead, each index gets a posterior mean, an HPD interval and, for every genotype, a probability of being the best. bam_indices() does this for five classical indices and bam_ammi_indices() for six AMMI-based ones. The cell means come from any fit: AMMI or GGE (bayes_ammi()) or a stability model (bayes_stability()).
2 Classical indices
For one draw of the cell means \(\eta_{ij}\), with interaction \(\epsilon_{ij} = \eta_{ij} - \bar{\eta}_{i\cdot} - \bar{\eta}_{\cdot j} + \bar{\eta}_{\cdot\cdot}\):
| Index | Formula | Better |
|---|---|---|
| Wricke’s ecovalence | \(W_i = \sum_j \epsilon_{ij}^2\) | smaller |
| Shukla’s stability variance | \(\frac{G}{(G-2)(E-1)}W_i - \frac{\sum_k W_k}{(G-1)(G-2)(E-1)}\) | smaller |
| Lin and Binns’ superiority | \(\sum_j (\eta_{ij} - M_j)^2 / (2E)\), \(M_j\) the best in environment \(j\) | smaller |
| Annicchiarico’s reliability | \(\bar{Z}_i - z_{1-p} s_{Z_i}\), \(Z_{ij} = 100\,\eta_{ij}/\bar{\eta}_{\cdot j}\) | larger |
| Eskridge’s safety-first | \(\bar{\eta}_{i\cdot} - z_{1-\alpha} s_i\), \(s_i^2\) the variance across environments | larger |
ind <-
bam_indices(
fit = fit
, higher_better = TRUE
, .prob = 0.25
, .alpha = 0.05
)
indPosterior stability indices (larger values of the trait are better)
Posterior means; rank_overall is the mean of the expected ranks over the indices (1 = best)
Gen wricke shukla lin_binns annicchiarico safety_first rank_overall
<char> <num> <num> <num> <num> <num> <num>
1: 6 2464000 135300 118900 104.50 3049 2.570
2: 5 2855000 161800 152200 102.30 2942 3.192
3: 3 1129000 44930 392500 96.20 2930 3.213
4: 7 1112000 43790 508100 92.99 2736 4.241
5: 2 2225000 119100 612900 90.25 2732 5.370
6: 9 3395000 198200 308900 94.29 2552 5.383
7: 4 7851000 499800 245700 96.98 2579 6.125
8: 1 4611000 280600 672000 88.58 2775 6.532
9: 8 7841000 499100 1159000 80.39 2550 8.374
Most probably the best on each index:
index Gen prob_best mean lower upper
<fctr> <char> <num> <num> <num> <num>
1: wricke 7 0.453 1112000.0 66.68 2938000.0
2: shukla 7 0.453 43790.0 -36900.00 167000.0
3: lin_binns 6 0.705 118900.0 27980.00 226900.0
4: annicchiarico 6 0.709 104.5 99.48 109.5
5: safety_first 6 0.542 3049.0 2742.00 3354.0
All HPD intervals and probabilities: x$indices
Sources: Wricke (1962); Shukla (1972); Lin and Binns (1988); Annicchiarico
(1992); Eskridge (1990). Full references: bam_cite(x).
plot(x = ind)
The indices measure different things. Wricke’s and Shukla’s indices measure stability alone and rank genotypes in the same order. Lin and Binns’ index measures distance from the best, so it mixes mean and stability. The rank correlation between Wricke’s and Lin and Binns’ posterior means is 0.08. Shukla’s estimator can be negative, as in the classical analysis; 5 of its 9 HPD intervals reach below zero, a sign that the stability variance of those genotypes cannot be told from zero.
Because the indices are computed from the cell means, without plot error, they describe the stability of the expected performance. The HPD intervals show how well 20 environments determine them; they are wide.
3 AMMI-based indices
With coordinates \(\mathrm{PC}_{ik} = \lambda_k^{1/2}\alpha_{ik}\) and the share \(\theta_k = \lambda_k^2/\sum_m \lambda_m^2\) of each fitted term:
| Index | Formula | Source |
|---|---|---|
| ASV | \(\sqrt{(\lambda_1^2/\lambda_2^2\,\mathrm{PC}_{i1})^2 + \mathrm{PC}_{i2}^2}\) | Purchase et al. (2000) |
| SIPC | \(\sum_k \lvert\mathrm{PC}_{ik}\rvert\) | Sneller et al. (1997) |
| EV | \(\sum_k \alpha_{ik}^2 / K\) | Zobel (1994) |
| Za | \(\sum_k \theta_k \lvert\alpha_{ik}\rvert\) | Zali et al. (2012) |
| WAAS | \(\sum_k \theta_k \lvert\mathrm{PC}_{ik}\rvert\) | Olivoto et al. (2019) |
| WAASY | mean and WAAS rescaled to 0-100 and weighted | Olivoto et al. (2019) |
ai <-
bam_ammi_indices(
fit = fit
, higher_better = TRUE
, .weight_y = 50
)
aiPosterior AMMI-based stability indices (2 term(s); WAASY with weight 50 on the mean)
Posterior means; rank_overall is the mean of the expected ranks over the indices (1 = best)
Gen ASV SIPC EV Za WAAS WAASY rank_overall
<char> <num> <num> <num> <num> <num> <num> <num>
1: 3 130.10 19.92 0.05095 0.1646 11.650 67.600 2.504
2: 7 99.51 22.02 0.08482 0.1422 9.614 63.650 2.658
3: 6 210.90 27.92 0.07374 0.2580 18.510 79.290 3.737
4: 2 192.50 28.22 0.08125 0.2425 17.210 48.010 4.273
5: 5 219.60 32.11 0.10290 0.2768 19.640 73.780 4.364
6: 9 252.50 32.34 0.09467 0.3057 22.000 54.950 5.142
7: 1 296.40 37.70 0.12270 0.3590 25.850 33.960 6.476
8: 4 402.40 43.50 0.15210 0.4679 34.200 48.210 7.565
9: 8 376.60 53.51 0.23690 0.4737 33.680 5.992 8.281
Most probably the best on each index:
index Gen prob_best mean lower upper
<fctr> <char> <num> <num> <num> <num>
1: ASV 7 0.583 99.51000 5.678e+00 204.1000
2: SIPC 3 0.395 19.92000 5.479e+00 36.0500
3: EV 3 0.385 0.05095 3.757e-06 0.1475
4: Za 7 0.526 0.14220 1.018e-02 0.2696
5: WAAS 7 0.551 9.61400 4.464e-01 18.5300
6: WAASY 6 0.566 79.29000 6.224e+01 100.0000
All HPD intervals and probabilities: x$indices
Sources: Purchase et al. (2000); Sneller et al. (1997); Zobel (1994); Zali et
al. (2012); Olivoto et al. (2019). Full references: bam_cite(x).
plot(x = ai)
WAASY gives the mean and the stability equal weights (.weight_y = 50). A breeder who values the mean more raises the weight:
data.table::rbindlist(
l = lapply(
X = c(30, 50, 70)
, FUN = function(wy) {
a <-
bam_ammi_indices(
fit = fit
, higher_better = TRUE
, .weight_y = wy
)
a$indices[index == "WAASY"][order(-mean)][1:3, list(weight = wy, Gen, mean, prob_best)]
}
)
) weight Gen mean prob_best
<num> <char> <num> <num>
1: 30 7 74.80483 0.319666667
2: 30 3 74.70340 0.296166667
3: 30 6 72.28332 0.216500000
4: 50 6 79.29195 0.565500000
5: 50 5 73.78358 0.297333333
6: 50 3 67.60028 0.087833333
7: 70 6 86.30058 0.696333333
8: 70 5 80.16992 0.292833333
9: 70 4 63.66496 0.0066666674 Indices from a stability model
The same classical indices can be computed from the cell means of a bayes_stability() fit:
fs <-
bayes_stability(
.data = Maiz
, .y = y
, .gen = entry
, .env = site
, .rep = rep
, .model = "ER"
, .nIter = 2000
, .burnin = 500
, .thin = 1
, .chains = 1
, .priors = bam_stab_priors(
mu_var = NULL
, g_var = NULL
, h_var = NULL
, b_var = 100
, nu = 1
, s2 = NULL
, nu_d = 1
, s2_d = NULL
)
, .level = 0.95
, .alpha = 0.05
, .seed = 11
, .waic = TRUE
)
bam_indices(
fit = fs
, higher_better = TRUE
, .prob = 0.25
, .alpha = 0.05
)$means Gen wricke shukla lin_binns annicchiarico safety_first rank_overall
<char> <num> <num> <num> <num> <num> <num>
1: 6 1822522 86285.70 156276.3 105.29231 3155.480 1.8692
2: 3 1607810 71756.31 489402.1 95.20404 2910.646 3.4322
3: 5 3901674 226980.19 226178.4 100.68529 2741.897 4.0675
4: 9 2996713 165742.18 409454.0 94.96091 2612.972 4.6945
5: 7 2905543 159572.83 674857.3 91.22647 2895.690 4.9025
6: 2 2475026 130440.08 719596.2 89.47825 2649.955 5.3462
7: 4 6580280 408239.21 274473.6 97.27475 2493.226 5.6621
8: 1 6387406 395187.58 772477.1 85.52653 2488.491 7.3905
9: 8 10736649 689497.31 1399778.8 80.73232 3017.037 7.63535 Interpretation
bam_interpret(
x = ind
, higher_better = TRUE
, top = 5
)Posterior stability indices
Best on average over the indices (mean expected rank): 6, 5, 3, 7 and 2.
Findings:
- Most probably the best on wricke: 7 (probability 0.45).
- Most probably the best on shukla: 7 (probability 0.45).
- Most probably the best on lin_binns: 6 (probability 0.71).
- Most probably the best on annicchiarico: 6 (probability 0.71).
- Most probably the best on safety_first: 6 (probability 0.54).
- The indices disagree on the best genotype: they weigh mean and stability
differently.
Recommendations:
- Prefer genotypes that rank well on several indices and whose HPD intervals
separate from the others; overlapping intervals mean the trial cannot tell
them apart.
Cautions:
- Stability indices describe the environments of this trial; confirm over
years before recommending.
Sources:
- Wricke, G. (1962). Uber eine Methode zur Erfassung der okologischen
Streubreite in Feldversuchen. Zeitschrift fur Pflanzenzuchtung, 47, 92-96.
- 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
- Lin, C. S. and Binns, M. R. (1988). A superiority measure of cultivar
performance for cultivar x location data. Canadian Journal of Plant
Science, 68(1), 193-198. doi:10.4141/cjps88-018
- Annicchiarico, P. (1992). Cultivar adaptation and recommendation from
alfalfa trials in Northern Italy. Journal of Genetics and Breeding, 46,
269-278.
- Eskridge, K. M. (1990). Selection of stable cultivars using a safety-first
rule. Crop Science, 30, 369-374.6 References
Wricke, G. (1962). Uber eine Methode zur Erfassung der okologischen Streubreite in Feldversuchen. Zeitschrift fur Pflanzenzuchtung, 47, 92–96.
Shukla, G. K. (1972). Some statistical aspects of partitioning genotype-environmental components of variability. Heredity, 29(2), 237–245.
Lin, C. S. and Binns, M. R. (1988). A superiority measure of cultivar performance for cultivar x location data. Canadian Journal of Plant Science, 68(1), 193–198. doi:10.4141/cjps88-018
Annicchiarico, P. (1992). Cultivar adaptation and recommendation from alfalfa trials in Northern Italy. Journal of Genetics and Breeding, 46, 269–278.
Eskridge, K. M. (1990). Selection of stable cultivars using a safety-first rule. Crop Science, 30, 369–374.
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.
Sneller, C. H., Kilgore-Norquest, L. and Dombek, D. (1997). Repeatability of yield stability statistics in soybean. Crop Science, 37(2), 383–390.
Zobel, R. W. (1994). Stress resistance and root systems. In Proceedings of the Workshop on Adaptation of Plants to Soil Stress, INTSORMIL Publication 94-2, 80–99. University of Nebraska-Lincoln.
Zali, H., Farshadfar, E., Sabaghpour, S. H. and Karimizadeh, R. (2012). Evaluation of genotype x environment interaction in chickpea using measures of stability from AMMI model. Annals of Biological Research, 3(7), 3126–3136.
Olivoto, T., Lucio, A. D. C., da Silva, J. A. G., Marchioro, V. S., de Souza, V. Q. and Jost, E. (2019). Mean performance and stability in multi-environment trials I: combining features of AMMI and BLUP techniques. Agronomy Journal, 111(6), 2949–2960. doi:10.2134/agronj2019.03.0220
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.