Data

The data set consists of potato field trials performed in years 2016–2018. The trials used in the study were conducted in experimental stations (sites) belonging to the Research Centre for Cultivar Testing (COBORU) (Table 1).

Table 1. Sites used in the 3-year variety trials conducted from 2016 to 2018

Most of the trials were organized in a randomized complete block design with three replicates. Only in 2018, at six sites, the trials were conducted according to a 1-resolvable design with three replicates and two blocks. On each plot, there were 60 plants planted. Plots that experienced unusual damage, such as boar damage, were excluded from further analyses. During the 3 years of study, in total, 22 varieties were tested. Since we were interested in assessing yielding stability, we choose to analyse varieties which were tested for 3 years. The list of varieties used in the present study with their country of origin is given in Table 2. During the 3-year period, most of the varieties were tested in all sites, only variety Everest was tested at five sites in years 2016 and 2017.

Table 2. Varieties used in the 3-year variety trials conducted from 2016 to 2018 and their country of origin

In Polish potato VCU and post-registration trials, one of the analysed characteristics is tuber yield, which is observed on plots. The observed tuber yields are expressed in decitonnes per hectare (dt/ha).

Statistical analysis

The observations were modelled according to a linear mixed model. For the sake of clarity, throughout the paper by environment we will mean a combination of year and location.

In the analysed data set, 1-resolvable design was used only in six out of 67 environments. Speed et al. (Reference Speed, Williams and Patterson1985) have shown that dropping the block effect for this design provides still a valid analysis. For this reason, we omitted in our model the effect of blocks nested within replicates and environments.

Let y _ikr be the potato tuber yield for the i-th variety (i  =  1, …, I) at the k-th environment (k  =  1, …, K) and in the r-th replicate (r  =  1, …, R). Then the model can be written as:

(1)$$y_{ikr} = \mu _i + u_k + v_{ik\;} + w_{kr} + e_{ikr}$$

where μ _i denotes the mean of the i-th variety. By u _k, v _ik , w _kr and e _ikr we denote in model (1) the random effects of environments, variety × environment interaction, of replicates nested within environment and of errors, respectively. Using vector notation, model (1) can be written as:

(2)$${\bi y} = {\bf X}{\bi \beta } + {\bf Z}_1{\bi u} + {\bf Z}_2{\bi v} + {\bf Z}_3{\bi w\;} + {\bi e}$$

where y is a vector of potato tuber yields and β is the vector of fixed variety means. By u, v, w, we denote in model (2) the random vectors of environment effects, of environment × variety interaction effects, and of replicates nested within environment effects, respectively. Further in model (2), by ${\bi e} = \lsqb {\bi {e}^{\prime}}_1\comma \;\ldots \comma \;\;{\bi {e}^{\prime}}_K\rsqb ^{\prime}$, where e_k  =  [e _1,k,1, …, e _I,k,R]^′ (k  =  1, …, K), we denote the random vector of errors. The matrices ${\bf X}$, ${\bf Z}_1$, ${\bf Z}_2$ and ${\bf Z}_3$ in model (2) are design matrices for β, u, v and w, respectively.

To fit the model, a hierarchical Bayesian linear mixed model was used. The model has three hierarchies or stages. In the first, it was assumed that observations are exchangeable samples from normal distribution:

$${\bi y}\vert {\bi \beta }\comma \;\;{\bi u}\comma \;\;{\bi v}\comma \;\;{\bi w}\comma \;{\bf R}_0\;\sim N\lpar {\bf X}{\bi \beta } + {\bf Z}_1{\bi u} + {\bf Z}_2{\bi v} + {\bf Z}_3{\bi w\;}\comma \;\oplus \sigma _k^2 {\bf I}_{n_k}\rpar $$

where ${\bf R}_0 = {\rm diag\lpar }\sigma _1^2 \comma \;\sigma _2^2 \comma \;\ldots \comma \;\sigma _K^2 \rpar$ is a matrix of residual variances, n _k is the number of plots in the k-th environment, and ${\oplus}$ denotes the direct sum of matrices. In case when the number of plots is equal in all trials, the above covariance matrix can be written as ${\bf R}_0\otimes {\bf I}$, where ⊗ denotes the Kronecker product.

In the second, the prior distributions are assigned to β, u, v and w; these correspond to the assumptions made on fixed and random effects in a Bayesian linear model. For the vector of variety means β a vague, large-variance Gaussian prior was assigned, i.e.

$${\bi \beta }\sim N\lpar 0\comma \;10^8{\bf I}\rpar $$

The vector of environment effects u was assigned a normal distribution with mean vector of zeros as:

$${\bi u}\vert \sigma _u^2 \sim N\lpar 0\comma \;\sigma _u^2 {\bf I}\rpar $$

whereas the vector of replicates nested within environments w was assigned a normal distribution with mean vector of zeros as:

$${\bi w}\vert \sigma _w^2 \sim N\lpar 0\comma \;\sigma _w^2 {\bf I}\rpar $$

For each trial (k  =  1, …, K), we assumed that ${\bi e}_k\vert \sigma _k^2 \sim N\lpar 0\comma \;\sigma _k^2 {\bf I}_{n_k}\rpar$. This implies that conditionally on $\sigma _1^2 \comma \;\ldots \comma \;\sigma _K^2$, the vector e follows

$${\bi e}\vert {\bf R}_0\;\sim N\;\lpar 0\comma \;\oplus \sigma _k^2 {\bf I}_{n_k}\rpar $$

The vector of environment × variety interaction v was assigned a normal distribution with mean vector of zeros and covariance matrix ${\bf G}_0\otimes {\bf I}$, i.e.

$${\bi v}\vert {\bf G}_0\;\sim N\lpar 0\comma \;{\bf G}_0\otimes {\bf I}\rpar .$$

Since we are interested in identifying the most stable varieties, by assuming different covariance matrices and modifying the list of random vectors in (2), we obtained different models. Our base model is Shukla model (SH) for which ${\bf G}_0 = {\rm diag\lpar }\sigma _{11}^2 \comma \;\ldots \comma \;\sigma _{II}^2 \rpar .$ Next, by removing from the model (2) the random vector of environment effects and assuming that ${\bf G}_0$ is completely unrestricted covariance matrix (i.e. the elements in the symmetric variance-covariance matrix ${\bf G}_0$ may take any value, as long ${\bf G}_0$), we obtained the environmental variance model (ENVIR; Piepho, Reference Piepho1999). Since in ENVIR model, the random vector of environment effects u is excluded, the random vector v is no longer vector of interaction effects rather a random vector of environment-variety effects. For both models, the diagonal elements of covariance matrix ${\bf G}_0$ are variances corresponding to the I varieties, which are interpreted as stability measures.

In the final or last stage, prior distributions were assigned to hyperparameters $\sigma _u^2$, $\;\sigma _w^2$, ${\bf R}_0{\bi \;}$and $\;{\bf G}_0$. We assumed that both $\sigma _u^2$ and $\;\sigma _w^2$ follow an inverse-Wishart distribution with 2.002 degrees of freedom and one-dimensional scale matrix equal to $\lpar 0.002/2.002\rpar {\bf I}$. For covariance matrix ${\bf G}_0$, as in Gelman et al. (Reference Gelman, Carlin, Stern, Dunson, Vehtari and Rubin2014), we assigned inverse-Wishart distribution prior distribution with ${\rm dim}{\bf G}_0 + 1$ degrees of belief. In our stability analysis, for covariance matrix ${\bf G}_0$ as the scale matrix, we used diagonal matrix of order ${\rm dim}{\bf G}_0$ with elements equal to empirical variances of each variety divided by ${\rm dim}{\bf G}_0$. Finally, for matrix ${\bf R}_0$, we assigned inverse-Wishart prior distribution with 2.002 degrees of belief. For matrix ${\bf R}_0$ as the scale matrix, we used the diagonal matrix of order K, with elements equal to (0.002/2.002). In the supplement of Longdon et al. (Reference Longdon, Hadfield, Webster, Obbard and Jiggins2011), it was pointed out that inverse-Wishart prior distribution defined in this way induces a marginal prior distribution on the variances that is equivalent to an inverse-γ distribution with shape and scale equal to 0.001.

Both SH and ENVIR models were fitted using R-package MCMCglmm (Hadfield, Reference Hadfield2012). Each model was run five times with different starting values. The length of each chain was set to 2 000 000 iterations with a burn-in period of 1 000 000 iterations and a thinning interval equal to 50. The convergence of chains was examined by visual inspection of trace plots and using the Gelman and Rubin (Reference Gelman and Rubin1992) convergence diagnostic (see also Cowles and Carlin, Reference Cowles and Carlin1996; Gelman et al., Reference Gelman, Carlin, Stern, Dunson, Vehtari and Rubin2014) implemented in the coda R-package (Plummer et al., Reference Plummer, Best, Cowles and Vines2006). For each model and chain, the deviance information criterion (Spiegelhalter et al., Reference Spiegelhalter, Best, Carlin and van der Linde2002; DIC; in smaller-is-better form; see also Sorensen and Gianola, Reference Sorensen and Gianola2002) was calculated. For clear presentation, only the results for the best model will be reported. For the best model, the estimates from all five chains were summarized using the summary function from coda R-package (Plummer et al., Reference Plummer, Best, Cowles and Vines2006).

Further, based on the obtained posterior mean yields and stability measure associated with the best model, all varieties were ranked. To combine the rankings, we applied Kang's (Reference Kang1988) non-parametric rank-sum stability procedure, where both posterior mean yield and posterior stability measure associated with the best model were used as a selection criterion. As in Kang (Reference Kang1988), the varieties with the lowest rank sum are the most desirable.

Further, for the best model, based on the variance components and the posterior variety means, coefficients of variation (CV) were calculated using the formula:

(3)$${\rm C}{\rm V}_i = \displaystyle{{{\hat{\tau }}_i} \over {{\hat{\mu }}_i}} \times 100\percnt \comma \;\quad i = 1\comma \;\ldots \comma \;I$$

where $\hat{\tau }_i$ is equal to $\sqrt {\hat{\sigma }_u^2 + \hat{\sigma }_{ii}^2 }$ for the SH model, $\sqrt {\hat{\sigma }_{ii}^2 }$ for the ENVIR model, and $\hat{\sigma }_{ii}^2$ is the estimate of the i-th diagonal element of ${\bf G}_0$ either in the SH model or the ENVIR model.

Since choosing which of the number of varieties for future sowing is a decision theoretic problem, optimum Bayesian choice of the best variety requires a utility function to be specified. In the present work, the critical-level approach to stability was used and for each variety, the utility function was defined in the following way: when the yield y is below the critical level γ, the utility is zero; above γ, the utility is equal to a positive constant c times y. Given the model parameters and assuming normality, the expected utility for variety i is equal to

(4)$$E\lpar {y_i\vert D\comma \;\gamma } \rpar = c\left[{\mu_i + \tau_i \times \varphi \left({\displaystyle{{\gamma -\mu_i} \over {\tau_i}}} \right)} \right]$$

where τ _i is equal to $\sqrt {\sigma _u^2 + \sigma _{ii}^2 }$ for the SH model, and $\sqrt {\sigma _{ii}^2 }$ for the ENVIR model,$\;\sigma _{ii}^2$ is the i-th diagonal element of ${\bf G}_0$ either in the SH model or the ENVIR model, φ( ⋅ ) is the density function of the standard normal distribution, and D, depending on the model, is a set of parameters of either the SH model or the ENVIR model. Given γ, the best variety is the one maximizing the posterior expected value of this expression. Since the parameters in the preceding equation are usually unknown, one has to use their estimates instead. From our observations, it follows that varieties with mean tuber yields above 42 t/ha are being recommended. For this reason, in the present work, we set the critical level γ equal to 42 t/ha. Moreover, from the discussions with farmers, it follows that the farmers' tuber yield is equal to 80% of the tuber yield obtained in the trials. For this reason, we set in (4) c equal to 0.8.

Next, for a risk analysis, the probabilities that the yields of varieties fall below a certain critical level in the environments were calculated. Since the probability of falling below a certain level depends on the variety mean and variance, it can be treated as a different stability measure which combines mean and variance in a meaningful way (Eskrigde, Reference Eskrigde1990; Eskridge and Mumm, Reference Eskridge and Mumm1992). We assumed that the tuber yield follows a normal distribution. The probability that the tuber yield of the i-th variety falls below a certain critical level in a randomly chosen environment is equal to (Eskrigde, Reference Eskrigde1990)

(5)$$p\lpar i\rpar = P\lpar y_{ik} \lt \delta \rpar = {\rm \Phi }\left({\displaystyle{{\delta -\mu_i} \over {\sqrt {\sigma_{ii}^2 } }}} \right)\comma \;\quad i = 1\comma \;\ldots \comma \;I$$

where μ _i is the mean for the i-th variety, $\sigma _{ii}^2$ is the variance of i-th variety across environments, Φ( ⋅ ) is the standard normal c.d.f.

Finally, by expressing p(i ₂) as a function of p(i ₁) (Piepho, Reference Piepho1996)

(6)$$p\;\lpar {i_2} \rpar = {\rm \Phi }\left({\displaystyle{{{\rm \Phi }^{{-}1}\lsqb {\,p\lpar i_1\rpar } \rsqb \sqrt {\sigma_{i_1i_1}^2 } + \mu_{i_1}-\mu_{i_2}} \over {\sqrt {\sigma_{i_2i_2}^2 } }}} \right)\comma \;\quad 0 \lt p\lpar i_1\rpar \lt 1$$

where Φ⁻¹( ⋅ ) is the inverse of the standard normal c.d.f., the relative risks for the best varieties from each country can be plotted. Since the parameters $\mu _{i_1}$, $\;\mu _{i_2}$, $\sigma _{i_1i_1}^2$, $\sigma _{i_2i_2}^2$ in the preceding equations are usually unknown, one has to use their estimates instead.