\[
\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
Crossa et al. (2011) fitted the Bayesian AMMI model with two bilinear terms to the maize trial Maiz. They ran two chains of 35000 iterations, discarded a burn-in of 5000 and thinned by 2. This vignette repeats that run and compares the posterior summaries with their Tables 1 and B1. The overall mean, the error variance, the main effects, the singular values and the scores of the first term agree within Monte Carlo error. The second scores differ, for a reason the paper does not discuss: the sign of a weak term.
The run of the paper
fit <-
bayes_ammi(
.data = Maiz
, .y = y
, .gen = entry
, .env = site
, .rep = rep
, .nIter = 30000
, .k = 2
, .burnin = 5000
, .thin = 2
, .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 = 2011
, .align = TRUE
, .model = "AMMI"
, .het = FALSE
, .waic = TRUE
)
fit$settings$seconds
Table 1
pub <- utils::read.csv(file = system.file("extdata", "crossa2011.csv", package = "bayesammi"))
s <- fit$summary
t1 <- merge(x = pub[pub$table == "1", ], y = s[, list(parameter, mean, sd, lower, upper)],
by = "parameter", suffixes = c(".paper", ".here"), sort = FALSE)
t1$dev_sd <- round(x = (t1$mean.here - t1$mean.paper) / t1$sd.paper, digits = 3)
num <- c("mean.paper", "mean.here", "sd.paper", "sd.here", "lower.paper", "lower.here", "upper.paper", "upper.here")
t1[num] <- lapply(X = t1[num], FUN = signif, digits = 5)
t1[, c("parameter", "mean.paper", "mean.here", "sd.paper", "sd.here", "dev_sd")]
parameter mean.paper mean.here sd.paper sd.here dev_sd
1 mu 4858.00 4857.600 32.71 32.539 -0.012
2 sigma2 763600.00 762820.000 43580.00 43631.000 -0.018
3 tau[1] -241.70 -241.050 91.78 92.348 0.007
4 tau[2] -255.00 -254.360 91.97 91.456 0.007
5 tau[3] -34.66 -35.607 91.75 91.862 -0.010
6 tau[4] 359.20 359.270 91.68 92.234 0.001
7 tau[5] 389.70 388.110 92.32 91.983 -0.017
8 tau[6] 472.50 471.850 91.47 92.025 -0.007
9 tau[7] -186.00 -185.970 92.19 92.082 0.000
10 tau[8] -574.90 -574.570 91.69 91.876 0.004
11 tau[9] 70.95 72.330 91.90 91.695 0.015
12 delta[1] -1245.00 -1244.600 141.80 142.300 0.003
13 delta[2] -647.30 -649.780 141.50 142.460 -0.018
14 delta[3] 246.60 246.630 141.60 141.440 0.000
15 delta[4] 371.60 371.080 141.60 142.390 -0.004
16 delta[5] 96.01 97.451 141.60 141.430 0.010
17 delta[6] 1448.00 1447.300 141.60 140.490 -0.005
18 delta[7] -1629.00 -1627.500 142.40 141.550 0.011
19 delta[8] -830.60 -831.110 142.50 142.130 -0.004
20 delta[9] 111.70 113.250 141.80 142.240 0.011
21 delta[10] -1921.00 -1920.500 141.80 141.840 0.003
22 delta[11] 449.20 449.080 143.00 143.990 -0.001
23 delta[12] 2659.00 2656.400 141.90 141.660 -0.018
24 delta[13] 1474.00 1473.400 142.60 141.550 -0.004
25 delta[14] 1194.00 1192.900 141.80 141.030 -0.008
26 delta[15] 189.30 189.090 141.60 141.600 -0.002
27 delta[16] 549.30 548.820 142.40 143.030 -0.003
28 delta[17] 17.62 18.970 141.70 142.970 0.010
29 delta[18] -305.50 -306.320 142.00 142.840 -0.006
30 delta[19] -2207.00 -2205.100 141.40 140.800 0.014
31 delta[20] -19.72 -19.467 140.80 141.200 0.002
32 lambda[1] 5483.00 5470.700 458.40 455.830 -0.027
33 lambda[2] 1627.00 1658.100 786.80 775.400 0.040
t1[, c("parameter", "lower.paper", "lower.here", "upper.paper", "upper.here")]
parameter lower.paper lower.here upper.paper upper.here
1 mu 4792.9000 4793.000 4921.0000 4920.100
2 sigma2 680570.0000 677110.000 849750.0000 847100.000
3 tau[1] -421.8000 -419.590 -61.8490 -56.881
4 tau[2] -433.0500 -435.610 -72.3490 -78.355
5 tau[3] -213.5800 -216.420 147.9200 144.300
6 tau[4] 180.6000 175.520 541.1600 535.570
7 tau[5] 206.3700 206.740 567.9200 566.120
8 tau[6] 296.9900 295.540 655.4900 654.760
9 tau[7] -363.9800 -369.400 -3.5244 -9.515
10 tau[8] -758.7400 -749.270 -401.0400 -392.460
11 tau[9] -107.7900 -107.690 247.4000 252.980
12 delta[1] -1518.8000 -1530.600 -963.8600 -974.240
13 delta[2] -925.3100 -925.220 -373.4200 -363.950
14 delta[3] -28.9430 -22.012 525.8100 529.360
15 delta[4] 93.7570 89.982 647.7600 649.060
16 delta[5] -186.8900 -180.950 365.2200 371.100
17 delta[6] 1170.0000 1178.100 1722.1000 1726.200
18 delta[7] -1906.1000 -1893.700 -1351.2000 -1337.500
19 delta[8] -1111.0000 -1111.500 -552.4400 -556.380
20 delta[9] -164.9300 -166.010 391.3200 392.030
21 delta[10] -2191.5000 -2197.800 -1637.0000 -1645.900
22 delta[11] 176.6500 168.460 738.1900 729.650
23 delta[12] 2387.1000 2385.300 2942.4000 2937.800
24 delta[13] 1188.2000 1194.900 1748.4000 1745.400
25 delta[14] 914.8800 912.030 1470.3000 1463.600
26 delta[15] -95.9140 -89.349 457.7000 462.400
27 delta[16] 265.9100 265.810 821.7100 826.410
28 delta[17] -265.4700 -254.970 290.3300 301.350
29 delta[18] -588.1600 -585.660 -34.3380 -28.091
30 delta[19] -2481.3000 -2481.600 -1928.5000 -1930.300
31 delta[20] -289.2900 -300.400 261.1700 250.560
32 lambda[1] 4573.8000 4581.700 6366.0000 6364.600
33 lambda[2] 2.0887 13.071 2822.7000 2833.800
The column dev_sd is the difference of the posterior means in units of the posterior standard deviation; its largest absolute value is 0.04. Both columns are Monte Carlo estimates, so they agree up to Monte Carlo error.
The error variance
The paper reports a prior scale of \(S_0 = 10^8\) with \(\nu = 1\), but its posterior of \(\sigma^2\) (mean 763600) is what a vague prior gives. Taken literally, with \(S^2_0 = 10^8\), the prior would add \(10^8/721\) to the posterior mean:
fit_s2 <-
bayes_ammi(
.data = Maiz
, .y = y
, .gen = entry
, .env = site
, .rep = rep
, .nIter = 3000
, .k = 2
, .burnin = 500
, .thin = 1
, .chains = 1
, .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 = 1e8
)
, .level = 0.95
, .seed = 1
, .align = TRUE
, .model = "AMMI"
, .het = FALSE
, .waic = TRUE
)
c(default = s[parameter == "sigma2", mean], s2_1e8 = fit_s2$summary[parameter == "sigma2", mean])
default s2_1e8
762815.3 922193.4
The package therefore defaults to the least-squares residual mean square as the scale, which reproduces Table 1.
Table B1: the first term
b1 <-
merge(
x = pub[pub$table == "B1" & grepl(pattern = ",1\\]$", x = pub$parameter), ]
, y = s[, list(parameter, mean, sd, lower, upper)]
, by = "parameter"
, suffixes = c(".paper", ".here")
, sort = FALSE
)
b1$dev_sd <- round(x = (b1$mean.here - b1$mean.paper) / b1$sd.paper, digits = 3)
b1[c("mean.here", "sd.here")] <- lapply(X = b1[c("mean.here", "sd.here")], FUN = round, digits = 4)
b1[, c("parameter", "mean.paper", "mean.here", "sd.paper", "sd.here", "dev_sd")]
parameter mean.paper mean.here sd.paper sd.here dev_sd
1 gamma[1,1] -0.2110 -0.2121 0.0751 0.0751 -0.015
2 gamma[2,1] 0.0169 0.0163 0.0757 0.0755 -0.008
3 gamma[3,1] -0.2427 -0.2426 0.0752 0.0748 0.001
4 gamma[4,1] 0.1389 0.1389 0.0734 0.0741 0.000
5 gamma[5,1] 0.0811 0.0808 0.0754 0.0751 -0.005
6 gamma[6,1] -0.0068 -0.0055 0.0766 0.0765 0.017
7 gamma[7,1] -0.0763 -0.0765 0.0753 0.0761 -0.003
8 gamma[8,1] -0.6353 -0.6340 0.0613 0.0614 0.021
9 gamma[9,1] 0.0359 0.0358 0.0752 0.0751 -0.002
10 gamma[10,1] -0.2052 -0.2056 0.0729 0.0730 -0.005
11 gamma[11,1] 0.3673 0.3676 0.0690 0.0690 0.005
12 gamma[12,1] 0.1510 0.1516 0.0754 0.0762 0.008
13 gamma[13,1] 0.2358 0.2356 0.0761 0.0764 -0.003
14 gamma[14,1] -0.0058 -0.0066 0.0743 0.0740 -0.011
15 gamma[15,1] 0.1950 0.1960 0.0755 0.0752 0.013
16 gamma[16,1] 0.0746 0.0747 0.0739 0.0740 0.002
17 gamma[17,1] -0.0008 -0.0004 0.0783 0.0781 0.006
18 gamma[18,1] 0.2042 0.2040 0.0769 0.0769 -0.003
19 gamma[19,1] -0.0826 -0.0836 0.0761 0.0770 -0.013
20 gamma[20,1] -0.0343 -0.0343 0.0760 0.0761 0.001
21 alpha[1,1] -0.3724 -0.3703 0.0742 0.0748 0.028
22 alpha[2,1] -0.2416 -0.2407 0.0756 0.0759 0.012
23 alpha[3,1] -0.1603 -0.1610 0.0753 0.0754 -0.009
24 alpha[4,1] 0.5016 0.5020 0.0647 0.0645 0.006
25 alpha[5,1] 0.2756 0.2759 0.0775 0.0773 0.003
26 alpha[6,1] 0.2649 0.2639 0.0727 0.0735 -0.013
27 alpha[7,1] -0.1152 -0.1158 0.0827 0.0820 -0.007
28 alpha[8,1] -0.4661 -0.4676 0.0776 0.0782 -0.019
29 alpha[9,1] 0.3135 0.3136 0.0741 0.0735 0.002
The largest absolute deviation is 0.028 posterior standard deviations.
Table B1: the second term
b2 <-
merge(
x = pub[pub$table == "B1" & grepl(pattern = ",2\\]$", x = pub$parameter), ]
, y = s[, list(parameter, mean, sd)]
, by = "parameter"
, suffixes = c(".paper", ".here")
, sort = FALSE
)
b2[c("mean.here", "sd.here")] <- lapply(X = b2[c("mean.here", "sd.here")], FUN = round, digits = 4)
b2[, c("parameter", "mean.paper", "mean.here", "sd.paper", "sd.here")]
parameter mean.paper mean.here sd.paper sd.here
1 gamma[1,2] 0.0840 -0.1314 0.2205 0.1979
2 gamma[2,2] -0.0562 0.0866 0.2060 0.1963
3 gamma[3,2] 0.0539 -0.0867 0.2187 0.2085
4 gamma[4,2] 0.0271 -0.0429 0.1869 0.1821
5 gamma[5,2] 0.0300 -0.0400 0.1952 0.1906
6 gamma[6,2] -0.0621 0.0973 0.2177 0.2032
7 gamma[7,2] 0.0331 -0.0524 0.1974 0.1910
8 gamma[8,2] -0.1302 0.2057 0.2395 0.1855
9 gamma[9,2] -0.0173 0.0230 0.1919 0.1899
10 gamma[10,2] 0.0184 -0.0340 0.1827 0.1792
11 gamma[11,2] -0.0189 0.0335 0.1773 0.1742
12 gamma[12,2] -0.0967 0.1561 0.2279 0.1934
13 gamma[13,2] -0.0762 0.1207 0.2350 0.2184
14 gamma[14,2] -0.0074 0.0114 0.1795 0.1768
15 gamma[15,2] 0.0884 -0.1396 0.2229 0.1959
16 gamma[16,2] -0.0014 0.0010 0.1822 0.1815
17 gamma[17,2] 0.1313 -0.2093 0.2619 0.2055
18 gamma[18,2] -0.1026 0.1661 0.2524 0.2159
19 gamma[19,2] 0.0246 -0.0398 0.2145 0.2127
20 gamma[20,2] 0.0784 -0.1252 0.2167 0.1916
21 alpha[1,2] -0.0841 0.1441 0.3087 0.2883
22 alpha[2,2] -0.1145 0.2020 0.2964 0.2433
23 alpha[3,2] -0.0568 0.0935 0.2620 0.2484
24 alpha[4,2] 0.0059 -0.0016 0.2222 0.2189
25 alpha[5,2] 0.1241 -0.2186 0.3260 0.2752
26 alpha[6,2] 0.0820 -0.1387 0.2619 0.2321
27 alpha[7,2] -0.1333 0.2362 0.3656 0.3086
28 alpha[8,2] 0.2483 -0.4364 0.4262 0.2385
29 alpha[9,2] -0.0716 0.1195 0.2805 0.2617
Here the means differ, mostly in sign, and the published means are smaller in absolute value. The second term is weak: the HPD interval of \(\lambda_2\) reaches down to near zero, both in the paper (2.09) and here (13.07). A chain without sign alignment then moves between the two mirror images \((\va_2, \vg_2)\) and \((-\va_2, -\vg_2)\). Its posterior means are averages over both images and shrink towards zero:
fit_free <-
bayes_ammi(
.data = Maiz
, .y = y
, .gen = entry
, .env = site
, .rep = rep
, .nIter = 30000
, .k = 2
, .burnin = 5000
, .thin = 2
, .chains = 1
, .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 = 2011
, .align = FALSE
, .model = "AMMI"
, .het = FALSE
, .waic = TRUE
)
d <- fit_free$draws[[1]][, "alpha[8,2]"]
c(share_positive = mean(d > 0), mean = mean(d), paper = pub$mean[pub$parameter == "alpha[8,2]"])
share_positive mean paper
0.51626667 0.01630027 0.24830000
The published value, 0.2483, is smaller in absolute value than the aligned posterior mean, -0.436, and has the opposite sign. A chain that spends a share \(q\) of its run in the image opposite to the least-squares orientation has a mean of about \((1 - 2q)\) times the aligned one. The published value corresponds to \(q \approx\) 0.78: the published chain spent most of its run in the other image. The aligned summaries in this package describe one orientation of the term, and the scores of genotypes then keep their meaning relative to each other.
plot(x = fit, which = "trace", pars = c("lambda[1]", "lambda[2]"))
plot(x = fit, which = "density", pars = c("lambda[1]", "lambda[2]"))
plot(x = fit, which = "biplot")
The coordinates are \(\lambda_k^{1/2}\alpha_{ik}\) and \(\lambda_k^{1/2}\gamma_{jk}\), as the paper defines them. The axes of its figures run from -1 to 1, so they show a rescaled version; the relative positions are the same.
bp <-
biplots(
model = fit
, burnin = 0
, thin = 0.25
, pb = c(0.05, 0.01)
, plot_stable = TRUE
, plot_unstable = TRUE
, ncolors = 5
, .which = "both"
, .grid = 100
)
bp$biplot
bp$biplot_data[axis1_significant == TRUE, list(byplot, name, x, axis1_lower, axis1_upper, group)]
byplot name x axis1_lower axis1_upper group
<char> <char> <num> <num> <num> <char>
1: Genotype 1 -27.40119 -39.009787 -17.0056562 negative
2: Genotype 2 -17.78439 -28.769313 -6.9059935 negative
3: Genotype 3 -11.92951 -22.387334 -0.5460779 negative
4: Genotype 4 37.11681 26.547850 46.4285760 positive
5: Genotype 5 20.33477 8.388511 30.8259011 positive
6: Genotype 6 19.51037 9.088838 30.3996106 positive
7: Genotype 8 -34.57251 -46.037088 -22.9895206 negative
8: Genotype 9 23.26336 11.850401 33.6434950 positive
9: Environment 1 -15.64906 -26.568920 -4.8800356 negative
10: Environment 3 -17.92078 -28.906052 -7.3859744 negative
11: Environment 8 -46.91103 -56.836455 -37.1519907 negative
12: Environment 10 -15.26108 -25.257391 -4.2795139 negative
13: Environment 11 27.22102 16.825997 37.3800310 positive
14: Environment 13 17.43228 6.882893 28.7666430 positive
15: Environment 15 14.44714 3.701040 25.1896927 positive
16: Environment 18 15.08558 3.669220 25.8574122 positive
On the first term, genotypes 4, 5, 6 and 9 form a positive group and genotypes 1 and 8 a negative group, as in the paper’s Figure 3. Sites S8 (negative) and S11, S13 and S18 (positive) stand out, as in its Figure 4. With the intervals of Table B1, genotypes 2 and 3 and some further sites also exclude zero at 95%; the paper’s figures show a subset of them.
References
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
