Optimizing the Allocation of Trials to Sub-regions in Multi-environment Crop Variety Testing

Abstract

New crop varieties are extensively tested in multi-environment trials in order to obtain a solid empirical basis for recommendations to farmers. When the target population of environments is large and heterogeneous, a division into sub-regions is often advantageous. When designing such trials, the question arises how to allocate trials to the different sub-regions. We consider a solution to this problem assuming a linear mixed model. We propose an analytical approach for computation of optimal designs for best linear unbiased prediction of genotype effects and their pairwise linear contrasts and illustrate the obtained results by a real data example from Indian nation-wide maize variety trials. It is shown that, except in simple cases such as a compound symmetry model, the optimal allocation depends on the variance–covariance structure for genotypic effects nested within sub-regions.

Appendices

Proof of Lemmas 1 and 2

To make use of the theoretical results that are available in the literature (see, e.g., Henderson 1975) for the prediction of random parameters, we will represent the model (1) as a particular case of the general LMM

$$\begin{aligned} \mathbf {Y}=\mathbf {X} {\varvec{\beta }} + \mathbf {Z} {\varvec{\zeta }} + {\varvec{\epsilon }} \end{aligned}$$

(29)

with design matrices $\mathbf {X}$ and $\mathbf {Z}$ for the fixed effects and the random effects, respectively. In (29), ${\varvec{\beta }}$ denotes the fixed effects and ${\varvec{\zeta }}$ are the random effects. The random effects and the observational errors ${\varvec{\epsilon }}$ are assumed to have zero mean and to be all uncorrelated with positive definite covariance matrices $\text{ Cov }\,({\varvec{\zeta }})=\mathbf {G}$ and $\text{ Cov }\,({\varvec{\epsilon }})=\mathbf {R}$, respectively. Random effects and observational errors are assumed to be uncorrelated.

To present model (1) in form (29), we follow the next steps:

$$\begin{aligned}&\mathbf {Y}_{ijk}=\mathbb {1}_L\,\mu _i+\mathbb {1}_L\,\alpha _{ik}+\mathbb {1}_L\,\lambda _{ij}+\mathbb {1}_L\,\gamma _{ijk}+\mathbf {b}_{ij}+{\varvec{\varepsilon }}_{ijk},\\&\quad i=1 \dots P,\quad k=1, \dots , K,\quad j=1, \dots , J_i, \end{aligned}$$

where $\mathbf {b}_{ij}=(b_{ij1}, \dots , b_{ijL})^\top$.

$$\begin{aligned} \mathbf {Y}_{ik}= & {} \mathbb {1}_{LJ_i}\,\mu _i+\mathbb {1}_{LJ_i}\,\alpha _{ik}+(\mathbb {I}_{J_i}\otimes \mathbb {1}_L)\,{\varvec{\lambda }}_{i}\\&+\,(\mathbb {I}_{J_i}\otimes \mathbb {1}_L)\,{\varvec{\gamma }}_{ik}+\mathbf {b}_{i}+{\varvec{\varepsilon }}_{ik},\\&\qquad i=1 \dots P,\quad k=1, \dots , K, \end{aligned}$$

where ${\varvec{\lambda }}_i=(\lambda _{i1}, \dots , \lambda _{iJ_i})^\top$ and ${\varvec{\gamma }}_{ik}=(\gamma _{i1k}, \dots , \gamma _{iJ_ik})^\top$.

$$\begin{aligned} \mathbf {Y}_{k}=\mathbf {F}{\varvec{\mu }}+\mathbf {F}{\varvec{\alpha }}_{k}+\mathbf {H}{\varvec{\lambda }}+\mathbf {H}{\varvec{\gamma }}_k+\mathbf {b}+{\varvec{\varepsilon }}_{k},\quad k=1, \dots , K, \end{aligned}$$

where $\mathbf {H}=(\mathbb {I}_J\otimes \mathbb {1}_L)$, ${\varvec{\mu }}=(\mu _1, \dots , \mu _P)^\top$, ${\varvec{\lambda }}=({\varvec{\lambda }}_1^\top , \dots , {\varvec{\lambda }}_P^\top )^\top$ and ${\varvec{\gamma }}_k=({\varvec{\gamma }}_{1k}^\top , \dots , {\varvec{\gamma }}_{Pk}^\top )^\top$.

$$\begin{aligned} \mathbf {Y}= & {} (\mathbb {1}_K\otimes \mathbf {F}){\varvec{\mu }}+(\mathbb {I}_K\otimes \mathbf {F}){\varvec{\alpha }}+(\mathbb {1}_K\otimes \mathbf {H}){\varvec{\lambda }}\\&+\,(\mathbb {I}_K\otimes \mathbf {H}){\varvec{\gamma }}+(\mathbb {1}_K\otimes \mathbb {I}_{LJ})\mathbf {b}+{\varvec{\varepsilon }}, \end{aligned}$$

where ${\varvec{\gamma }}=({\varvec{\gamma }}_1^\top , \dots , {\varvec{\gamma }}_K^\top )^\top$.

The latter equation may alternatively be written as

$$\begin{aligned} \mathbf {Y}=(\mathbb {1}_K\otimes \mathbf {F}){\varvec{\mu }}+(\mathbb {I}_K\otimes \mathbf {F}){\varvec{\alpha }}+\tilde{{\varvec{\varepsilon }}}, \end{aligned}$$

(30)

where $\tilde{{\varvec{\varepsilon }}}:=(\mathbb {1}_K\otimes \mathbf {H}){\varvec{\lambda }}+(\mathbb {I}_K\otimes \mathbf {H}){\varvec{\gamma }}+(\mathbb {1}_K\otimes \mathbb {I}_{LJ})\mathbf {b}+{\varvec{\varepsilon }}$. Model (30) is of form (29) with $\mathbf {X}=(\mathbb {1}_K\otimes \mathbf {F})$, $\mathbf {Z}=(\mathbb {I}_K\otimes \mathbf {F})$, $\mathbf {G}=\mathrm {Cov}({\varvec{\alpha }})=\sigma ^2\mathbb {I}_K\otimes \mathbf {D}$ and

$$\begin{aligned} \mathbf {R}= & {} \mathrm {Cov}(\tilde{{\varvec{\varepsilon }}})=\sigma ^2((v_1\mathbb {1}_K\mathbb {1}_K^\top +v_2\mathbb {I}_K)\otimes \mathbb {I}_J\otimes (\mathbb {1}_L\mathbb {1}_L^\top )\\&+\,v_3(\mathbb {1}_K\mathbb {1}_K^\top )\otimes \mathbb {I}_{LJ}+\mathbb {I}_{LJK}). \end{aligned}$$

According to Henderson (1975), the BLUP of the random effects ${\varvec{\zeta }}$ (which corresponds to ${\varvec{\alpha }}$ in our model (30)) is given by

$$\begin{aligned} \hat{{\varvec{\zeta }}}= & {} \left( \mathbf {Z}^\top \mathbf {R}^{-1}\mathbf {Z}+\mathbf {G}^{-1} -\mathbf {Z}^\top \mathbf {R}^{-1}\mathbf {X} (\mathbf {X}^\top \mathbf {R}^{-1}\mathbf {X})^{-}\mathbf {X}^\top \mathbf {R}^{-1}\mathbf {Z}\right) ^{-1}\nonumber \\&\cdot \left( \mathbf {Z}^\top \mathbf {R}^{-1} -\mathbf {Z}^\top \mathbf {R}^{-1}\mathbf {X} (\mathbf {X}^\top \mathbf {R}^{-1}\mathbf {X})^{-}\mathbf {X}^\top \mathbf {R}^{-1}\right) \mathbf {Y}. \end{aligned}$$

(31)

Using this formula, we obtain the BLUP for the genotype effects ${\varvec{\alpha }}$, which results in formula (2). The MSE matrix of the BLUP of the random effects ${\varvec{\zeta }}$ is given by

$$\begin{aligned} \mathrm {Cov}(\hat{{\varvec{\zeta }}}-{\varvec{\zeta }})=\left( \mathbf {Z}^\top \mathbf {R}^{-1}\mathbf {Z} +\mathbf {G}^{-1}-\mathbf {Z}^\top \mathbf {R}^{-1}\mathbf {X} (\mathbf {X}^\top \mathbf {R}^{-1}\mathbf {X})^{-}\mathbf {X}^\top \mathbf {R}^{-1}\mathbf {Z}\right) ^{-1}, \end{aligned}$$

(32)

where $\mathbf {A}^{-}$ denotes a generalized inverse of $\mathbf {A}$. By this formula, we obtain MSE matrix (4). Then, using the relation ${\varvec{\theta }}^{k,k'}={\varvec{\alpha }}_k-{\varvec{\alpha }}_{k'}=((\mathbf {e}_k-\mathbf {e}_{k'})^\top \otimes \mathbb {I}_P)\,{\varvec{\alpha }}$ between the genotype effects and their pairwise contrasts we obtain formulae (3) and (5).

Table 7 Optimal numbers of locations per sub-region and efficiency of balanced design compared to optimal designs with respect to the standard A-criterion in the FA model (late maturity) for different values of the total number of locations J and the error variance $\sigma ^2$.

Table 8 Optimal numbers of locations per sub-region and efficiency of balanced design compared to optimal designs with respect to the standard A-criterion in the FA model (medium maturity) for different values of the total number of locations J and the error variance $\sigma ^2$.

Sensitivity Analysis

1.1 Standard A-Criterion

We take values of the covariance matrix $\mathbf {V}$ from Tables 3, 4 and 5 in Kleinknecht et al. (2013) for late, medium and early maturity. Tables 7, 8 and 9 summarize the results for optimal designs for the standard A-criterion in FA model for late, medium, and early maturity, respectively.

Table 9 Optimal numbers of locations per sub-region and efficiency of balanced design compared to optimal designs with respect to the standard A-criterion in the FA model (early maturity) for different values of the total number of locations J and the error variance $\sigma ^2$.

Table 10 Optimal numbers of locations per sub-region and efficiencies of balanced and weighted designs compared to optimal designs with respect to the weighted A-criterion for the FA model (late maturity) for different values of the total number of locations J and the error variance $\sigma ^2$.

Table 11 Optimal numbers of locations per sub-region and efficiencies of balanced and weighted designs compared to optimal designs with respect to the weighted A-criterion for the FA model (medium maturity) for different values of the total number of locations J and the error variance $\sigma ^2$.

Table 12 Optimal numbers of locations per sub-region and efficiencies of balanced and weighted designs compared to optimal designs with respect to the weighted A-criterion for the FA model (early maturity) for different values of the total number of locations J and the error variance $\sigma ^2$.

Fig. 1

Efficiencies $\mathrm {Eff}_{a,P}$ (dark dots) $\mathrm {Eff}_{a,\ell }$ (big light dots) and $\mathrm {Eff}_{e,P}$ (small light dots) in dependence on the total number of allocations J for weighted A-criterion in the CS model for late (left panel), medium (middle panel) and early (right panel) maturity for $\sigma ^2=50$ (first row), $\sigma ^2=200$ (second row) and $\sigma ^2=400$ (third row)

1.2 Weighted A-Criterion

Tables 10, 11, and 12 summarize the results optimal designs with respect to the weighted A-criterion in the FA model for late, medium, and early maturity, respectively.

Figure 1 illustrates the behavior of efficiencies of balanced and weighted designs with respect to optimal approximate and exact designs ($\mathrm {Eff}_{a,P}$, $\mathrm {Eff}_{a,\ell }$ and $\mathrm {Eff}_{e,P}$ as in Sect. 4) in dependence on the total number of allocations J for weighted A-criterion in the CS model. For J, we considered all multiples of 5 between 15 and 200. The error variance is fixed at $\sigma ^2=50$, $\sigma ^2=200$ and $\sigma ^2=400$.