- Email: [email protected]
Block designs for random treatment effects
Journal of Statistical Planning and Inference 137 (2007) 1446 – 1451 www.elsevier.com/locate/jspi
Block designs for random treatment effects Júlio S. de S. Bueno Filhoa , Steven G. Gilmourb,∗ a Departamento de Ciências Exatas, CP37 - UFLA, 37200-000 Lavras, MG - Brazil b School of Mathematical Sciences, Queen Mary, University of London, Mile End Road, London E1 4NS, UK
Received 13 September 2004; received in revised form 3 January 2006; accepted 2 February 2006 Available online 30 March 2006
Abstract A Bayesian design criterion for selection experiments in plant breeding is derived using a utility function that minimizes the risk of an incorrect selection. A prior distribution on the heritability parameter is used to complete the definition of the design optimality criterion. An example is given with evaluations of the criterion for different prior distributions on the heritability. Though coming from a genetic motivation this criterion should prove useful for any other types of experiments with random treatment effects. © 2006 Elsevier B.V. All rights reserved. Keywords: Genetic relatedness; Mixed model; Optimal design; REML
1. The problem In selection trials in plant and animal breeding the relative genetic merits of individuals must be estimated. The analysis of trials whose main objective is selection consists of fitting linear mixed models with genetic breeding values taken as random treatment effects, subject to a genetic covariance proportional to the pairwise relationship matrix (Henderson, 1975), and random block factor effects. REML estimation (Patterson and Thompson, 1971) is seen as the state of the art for the analysis of such a model, either in animal breeding (Henderson, 1985) or in plant breeding (Bernardo, 1996). Selection trials are usually carried out as incomplete block designs. Bueno Filho and Gilmour (2003) showed that information on genetical relatedness can improve the design of experiments, especially incomplete block experiments, in which the recovery of intertreatment information may play an important role. In the present paper, we provide a Bayesian justification for the optimal design criterion used there. Although it has a genetic motivation, this criterion should prove useful for other types of experiments with random treatment effects, since we show that it gives optimal designs for a broad class of utility functions. Section 2 details the model assumptions and gives our main result. A small example is used in Section 3 to introduce the benefits of using a prior distribution for the heritability parameter and results for this example are discussed in Section 4. Some general conclusions are drawn in Section 5.
∗ Corresponding author. Tel.: +44 20 7882 7833; fax: +44 20 8981 9587.
E-mail addresses: [email protected]fla.br (J.S. de S. Bueno Filho), [email protected] (S.G. Gilmour). 0378-3758/$ - see front matter © 2006 Elsevier B.V. All rights reserved. doi:10.1016/j.jspi.2006.02.002
J.S. de S. Bueno Filho, S.G. Gilmour / Journal of Statistical Planning and Inference 137 (2007) 1446 – 1451
1447
2. The model and design criterion Consider a model of the form y = X + Z + , where y is an n × 1 observable random vector, X is an n × b design matrix for nuisance parameters (blocks), is a b × 1 vector of unknown block means, Z is an n × t design matrix for treatments, is a t × 1 vector of unknown random treatment effects and is an n × 1 random vector of unknown errors. For the purposes of design construction it is convenient to treat block effects as fixed, so that E(y) = X. We assume that V () = 2 I and V () = 2 A, where 2 is a constant associated with the variance of treatment effects and A = {ij } is the matrix of the pairwise correlations between treatment effects. The elements of A are obtained directly from the knowledge of the relationship between individuals, known as identity-by-descent—see, for example, Lange (1997). Bueno Filho and Gilmour (2003) showed that (1/2 )V (ˆ) = , where −1 = r + A−1 − N (k )−1 N , k and r are diagonal matrices of the numbers of plots per block and the numbers of replications per treatment, respectively, = 2 /2 and N is the block-treatment concurrence matrix. Thus, the properties of must be studied to develop design criteria. Without loss of generality we will take 2 = 1 so that, in what follows, variances refer to scaled variances. Any useful optimality criterion for selection experiments should be based on discriminating between good and bad genetic material. Assume that we intend to select s out of t treatments for further experimentation. Let S be the set of the best s treatments. Then we are particularly interested in the comparisons ij = i − j , for i ∈ S and j ∈ / S. To ˆ minimize the probability of a particular mistake in selection, we must minimize P (ij < 0). Given unbiased estimators ˆ ij ). Since any such mistake leads to an incorrect selection, the design of the ’s, this is equivalent to minimizing V ( ˆ criterion should minimize i∈S j ∈S / V (ij ). The same criterion arises from quadratic loss utility on the comparisons of interest, i.e. U (Z, y, ) = −L(Z, y, ), where L(Z, y, ) =
(1)
ˆ ij − ij )2 . (
(2)
i∈S j ∈S /
We should choose Z to minimize L(Z, y, ), but this also depends on the data, y, and on the unknown parameters, . Integrating first over y|, gives ˆ ij ), Ey| {L(Z, y, )} = V ( i∈S j ∈S /
ˆ ij ) does not depend on and ˆ ij is an unbiased estimator of ij . since V ( Next we integrate over the prior distribution of to get ⎧ ⎫ ⎨ ⎬ ˆ ij ) E [Ey| {L(Z, y, )}] = E V ( ⎩ ⎭ i∈S j ∈S / ˆ ij )P (R = S), = V ( R∈ i∈R j ∈R /
where is the set of all possible sets of s treatments. If the prior is such that P (R = S) is constant, i.e. any treatment is equally likely to be among the best s, then E [Ey| {L(Z, y, )}] ∝
t−1 t i=1 j =i+1
This establishes the following result.
ˆ ij ). V (
1448
J.S. de S. Bueno Filho, S.G. Gilmour / Journal of Statistical Planning and Inference 137 (2007) 1446 – 1451
Theorem 1. Given any prior distribution for such that the prior probability that treatment i is among the best s treatments is constant ∀i ∈ {1, . . . , t} and given the utility function defined by Eqs. (1) and (2), maximizing the expected utility is equivalent to minimizing the criterion C(Z|, A) =
t−1 t
V (ˆi − ˆ j ),
(3)
i=1 j =i+1
assuming and A are known. This criterion is the one used by Bueno Filho and Gilmour (2003). In fact this result is even more general, since s need not be fixed and need not represent the number of treatments we will select. Relabel the above utility function Us (Z, y, ). Then our utility function can be any weighted average of Us (Z, y, ), for s=1, . . . , t. This is quite important in practice, since it is not usually the case that the number of treatments that will be selected from the experiment is fixed in advance of the experiment. Rather, an appropriate number to ensure enough variation in the selected population will be selected. This criterion is also not specific to the problem of selection, although that is our primary objective. It would also be a sensible utility function for ranking treatments and for estimating treatment effects. In fact, it corresponds to the commonly used quadratic loss utility function for estimating treatment differences. 3. Prior distribution for heritability Consider the following designs for comparing 6 treatments in 4 blocks of size 3: D1 D2 D3 D4 D5
: {1, 2, 3}, {1, 5, 6}, {2, 4, 6}, {3, 4, 5}, : {1, 3, 6}, {1, 4, 5}, {2, 3, 5}, {2, 4, 6}, : {1, 2, 5}, {1, 3, 4}, {2, 5, 6}, {3, 4, 6}, : {1, 3, 6}, {1, 4, 5}, {2, 3, 6}, {2, 4, 5}, : {1, 2, 3}, {1, 2, 3}, {4, 5, 6}, {4, 5, 6},
where Di denotes the design and braces denote blocks. If we consider the treatment labels to be exchangeable, the regular graph designs D1 and D2 are isomorphic, as are the partially balanced incomplete block designs D3 and D4 . The treatment concurrence matrices and the graphs for these designs were given by Bueno Filho and Gilmour (2003) who plotted the values of C(Z|, A) for each design for a range of different values and for various family structures, including one defined by ⎡1 ⎤ 0 0 0 18 81 2 ⎢0 1 0 1 0 1 ⎥ ⎢ 2 8 8⎥ ⎢ ⎥ ⎢0 0 1 1 1 0⎥ ⎢ ⎥ 2 8 8 A=⎢ ⎥. ⎢0 1 1 1 0 0⎥ 8 8 2 ⎢ ⎥ ⎢1 ⎥ ⎣ 8 0 18 0 21 0 ⎦ 1 8
1 8
0
0
0
1 2
With such structures, the treatment labels can no longer be considered to be exchangeable and so the five designs become genuinely different. The model with fixed treatment effects can be obtained by making → 0 (i.e. 2 → ∞), so the best designs in this situation are also the best designs for genetically unrelated treatments. For this situation D5 is the worst possible, being disconnected for fixed treatments. On the other hand, D5 is a good design for greater than approximately 0.4. This surprising result arises because the analysis with random treatment effects allows information on comparisons involving a particular treatment to be “borrowed” from comparisons involving related treatments. The use of point prior values for does not account for uncertainty in its value, and we may be tempted to identify a design as the best when it is probably not. Instead we can define a prior distribution on and find the optimal design
J.S. de S. Bueno Filho, S.G. Gilmour / Journal of Statistical Planning and Inference 137 (2007) 1446 – 1451
1449
averaged over this prior. The most usual way to describe the ratio of variances between treatments and experimental units in plant breeding is through the heritability, h2 = 2 /(2 + 2 ) = 1/(1 + ). This is a measure of the covariance between relatives as a proportion of the total phenotypic variance in experimental units and 0 < h2 < 1. Heritability values can be given prior probability distributions and prior elicitation is relatively straightforward because plant breeders have estimates of heritability from large amounts of previous data. Assuming some prior probability distribution on h2 , which implies a prior for it is possible to use our Bayesian criterion averaged over this prior. Thus, our utility function becomes E,h2 {U (Z, y, , h2 )} = Eh2 [E|h2 {U (Z, y, |h2 )}] = − Eh2 {C(Z|h2 , A)} 1 C(Z|h2 , A)p(h2 ) dh2 . = − 0
In general the evaluation of this utility function for any design might require Monte Carlo sampling from the prior distribution for h2 . However, in this small example it is possible to get analytical solutions through some lengthy, but straightforward, algebraic manipulations. A Beta prior distribution for h2 will be flexible enough in many practical situations and is easily interpreted. The Beta distribution with parameters and has probability density function 1 p(h2 ) = (h2 )( −1) (1 − h2 )( −1) /B( , ), where B( , ) = 0 u −1 (1 − u) −1 du. It has expectation /( + ) and, when used as a prior distribution, + represents the equivalent sample size of the prior information. The Beta distribution with parameters and is implied if we have independent priors for 2 and 2 which are, respectively, Inverse- 22 and Inverse- 22 (which also implies that / ∼ F2 ,2 ), although it is not necessary to assume this for the Beta prior for h2 to arise. The steps in a numerical evaluation of the criterion values are as follows:
1. for a known A matrix, write out the design matrix and the information matrix as a function of and hence as a function of h2 ; 2. write the criterion conditional on h2 , according to Eq. (3); 3. assign a Beta prior for h2 with and chosen to represent the prior mean heritability by /( + ) and the equivalent sample size by + ; 4. evaluate the integrated criterion value, taking the prior distribution of h2 as weights, 1 2 E,h2 {U (Z, y, , h )} = − C(Z|h2 , A)p(h2 ) dh2 . 0
In the next section, we present the results of following these steps to evaluate the Bayesian design criterion for each of the 5 designs considered, averaged over Beta priors for the heritabilities. 4. Results Tables 1 and 2 show some values of the Bayesian criterion calculated for 36 different Beta priors for the heritabilities. Moving down Table 1, the priors become more informative, while Table 2 represents priors with no uncertainty, and moving from left to right in each table the expected heritability becomes greater. The priors studied have many different shapes, including those that are highly skewed to the left and right, those that are bimodal, one which is uniform on [0, 1], those which tend towards normality, point priors and various others in between. It can be seen that for a wide range of flat priors the best design is in the class of regular graph designs (in fact D2 is clearly best). However, when the prior distribution gives a high probability to low heritabilities, we see that D5 can be the best. The tables show that in the case of little prior knowledge, it is safest to use a standard regular graph design. However, the prior knowledge that the heritability is low does not have to be very great before D5 becomes better, if only slightly better. With the prior having information equivalent to a sample size of 2, for example, D5 is the best for an average heritability of 0.1. For strong prior knowledge on heritability (equivalent information to an infinite sample size) regular graph designs are better just in situations with very high heritabilities (h2 > 0.7), situations that are unusual for many important economic traits.
1450
J.S. de S. Bueno Filho, S.G. Gilmour / Journal of Statistical Planning and Inference 137 (2007) 1446 – 1451
Table 1 Values of the design criterion and best designs (shown in bold) for Beta priors with different and Design
Beta parameters ( , ) (0.1,0.9)
(0.3,0.7)
(0.5,0.5)
(0.7,0.3)
(0.9,0.1)
D1 D2 D3 D4 D5
1.55 1.54 1.56 1.56 1.58
4.79 4.77 4.86 4.85 5.05
8.24 8.22 8.44 8.42 9.09
11.93 11.91 12.35 12.33 13.99
15.90 15.89 16.67 16.66 20.20
D1 D2 D3 D4 D5
(0.2,1.8) 1.48 1.47 1.48 1.48 1.46
(0.6,1.4) 4.61 4.59 4.63 4.62 4.61
(1,1) 7.99 7.96 8.1 8.08 8.3
(1.4,0.6) 11.69 11.66 12.01 11.98 12.98
(1.8,0.2) 15.78 15.77 16.49 16.47 19.53
D1 D2 D3 D4 D5
(0.5,4.5) 1.42 1.42 1.41 1.41 1.39
(1.5,3.5) 4.45 4.42 4.42 4.42 4.32
(2.5,2.5) 7.76 7.71 7.77 7.76 7.66
(3.5,1.5) 11.44 11.40 11.64 11.61 11.94
(4.5,0.5) 15.64 15.63 16.28 16.25 18.61
D1 D2 D3 D4 D5
(1,9) 1.40 1.39 1.39 1.39 1.37
(3,7) 4.38 4.35 4.34 4.34 4.24
(5,5) 7.65 7.6 7.63 7.62 7.44
(7,3) 11.32 11.28 11.46 11.43 11.49
(9,1) 15.57 15.55 16.16 16.13 18.03
D1 D2 D3 D4 D5
(5,45) 1.38 1.38 1.37 1.37 1.36
(15,35) 4.32 4.29 4.28 4.27 4.18
(25,25) 7.56 7.51 7.51 7.49 7.28
(35,15) 11.21 11.16 11.29 11.26 11.09
(45,5) 15.49 15.47 16.03 15.99 17.32
D1 D2 D3 D4 D5
(10,90) 1.38 1.37 1.37 1.37 1.36
(30,70) 4.31 4.28 4.27 4.27 4.17
(50,50) 7.55 7.49 7.49 7.48 7.26
(70,30) 11.19 11.14 11.27 11.23 11.05
(90,10) 15.48 15.46 16.01 15.97 17.21
Table 2 Values of the design criterion and best designs (shown in bold) for point priors with different h2 Design
D1 D2 D3 D4 D5
h2 0.1
0.3
0.5
0.7
0.9
1
1.38 1.37 1.37 1.37 1.36
4.31 4.28 4.26 4.26 4.16
7.53 7.48 7.47 7.46 7.24
11.18 11.12 11.24 11.21 11.00
15.46 15.45 16.00 15.95 17.09
18.00 18.00 19.00 19.00 24.00
J.S. de S. Bueno Filho, S.G. Gilmour / Journal of Statistical Planning and Inference 137 (2007) 1446 – 1451
1451
5. Final remarks For common situations of selection in plant breeding we are very likely to have either a sparse genetical relationship matrix, A, or discrete family structures. For some of these situations it will be possible to find optimal designs robust to assumptions on (assumptions on heritabilities). However, for other situations (complex pedigrees) the design choice may change dramatically at particular values using point prior estimates. Our Bayesian optimality design criterion results from a general framework that encompasses all cases, and involves just sampling from some reasonable distribution on heritabilities. For flat priors on heritabilities we have found that the good designs will usually be in the class of regular graph designs. Our results use utility functions and prior distributions to find optimal designs for likelihood-based analysis. It could be argued that we should use our prior information in the analysis of data. This is quite possible, of course, but is a different problem. If the results are to be made public, then the experiment has to be designed with a view to analyzing the data using a population of prior distributions, not just the experimenter’s. Our method makes use of the experimenter’s prior knowledge for the design, but recognizes that in practice the data will be analyzed by likelihood methods. Bayesian design for likelihood analysis is reviewed by Chaloner and Verdinelli (1996), who suggest that it should be used much more often in practice. Acknowledgements This work was conceived when the first author visited Queen Mary, University of London, with financial support from CAPES grant BEX-1248/00-6. It was continued when the second author visited Universidade Federal de Lavras, with financial support from the Engineering and Physical Sciences Research Council grant GR/S14009/01. We thank the referees for comments which helped us to improve the presentation of this work. References Bernardo, R., 1996. Best linear unbiased prediction of maize single-cross performance. Crop Sci. 36, 50–56. Bueno Filho, J.S. de S., Gilmour, S.G., 2003. Planning incomplete block experiments when treatments are genetically related. Biometrics 59, 375–381. Chaloner, K., Verdinelli, I., 1996. Bayesian experimental design: a review. Statist. Sci. 10, 270–304. Henderson, C.R., 1975. Best linear unbiased estimation and prediction under a selection model. Biometrics 31, 423–449. Henderson, C.R., 1985. Best linear unbiased prediction of nonadditive genetic merits in noninbred populations. J. Animal Sci. 60, 111–117. Lange, K., 1997. Mathematical and Statistical Methods for Genetic Analysis. Springer, New York. Patterson, H.D., Thompson, R., 1971. Recovery of inter-block information when block sizes are unequal. Biometrika 58, 545–554.
Block designs for random treatment effects Júlio S. de S. Bueno Filhoa , Steven G. Gilmourb,∗ a Departamento de Ciências Exatas, CP37 - UFLA, 37200-000 Lavras, MG - Brazil b School of Mathematical Sciences, Queen Mary, University of London, Mile End Road, London E1 4NS, UK
Received 13 September 2004; received in revised form 3 January 2006; accepted 2 February 2006 Available online 30 March 2006
Abstract A Bayesian design criterion for selection experiments in plant breeding is derived using a utility function that minimizes the risk of an incorrect selection. A prior distribution on the heritability parameter is used to complete the definition of the design optimality criterion. An example is given with evaluations of the criterion for different prior distributions on the heritability. Though coming from a genetic motivation this criterion should prove useful for any other types of experiments with random treatment effects. © 2006 Elsevier B.V. All rights reserved. Keywords: Genetic relatedness; Mixed model; Optimal design; REML
1. The problem In selection trials in plant and animal breeding the relative genetic merits of individuals must be estimated. The analysis of trials whose main objective is selection consists of fitting linear mixed models with genetic breeding values taken as random treatment effects, subject to a genetic covariance proportional to the pairwise relationship matrix (Henderson, 1975), and random block factor effects. REML estimation (Patterson and Thompson, 1971) is seen as the state of the art for the analysis of such a model, either in animal breeding (Henderson, 1985) or in plant breeding (Bernardo, 1996). Selection trials are usually carried out as incomplete block designs. Bueno Filho and Gilmour (2003) showed that information on genetical relatedness can improve the design of experiments, especially incomplete block experiments, in which the recovery of intertreatment information may play an important role. In the present paper, we provide a Bayesian justification for the optimal design criterion used there. Although it has a genetic motivation, this criterion should prove useful for other types of experiments with random treatment effects, since we show that it gives optimal designs for a broad class of utility functions. Section 2 details the model assumptions and gives our main result. A small example is used in Section 3 to introduce the benefits of using a prior distribution for the heritability parameter and results for this example are discussed in Section 4. Some general conclusions are drawn in Section 5.
∗ Corresponding author. Tel.: +44 20 7882 7833; fax: +44 20 8981 9587.
E-mail addresses: [email protected]fla.br (J.S. de S. Bueno Filho), [email protected] (S.G. Gilmour). 0378-3758/$ - see front matter © 2006 Elsevier B.V. All rights reserved. doi:10.1016/j.jspi.2006.02.002
J.S. de S. Bueno Filho, S.G. Gilmour / Journal of Statistical Planning and Inference 137 (2007) 1446 – 1451
1447
2. The model and design criterion Consider a model of the form y = X + Z + , where y is an n × 1 observable random vector, X is an n × b design matrix for nuisance parameters (blocks), is a b × 1 vector of unknown block means, Z is an n × t design matrix for treatments, is a t × 1 vector of unknown random treatment effects and is an n × 1 random vector of unknown errors. For the purposes of design construction it is convenient to treat block effects as fixed, so that E(y) = X. We assume that V () = 2 I and V () = 2 A, where 2 is a constant associated with the variance of treatment effects and A = {ij } is the matrix of the pairwise correlations between treatment effects. The elements of A are obtained directly from the knowledge of the relationship between individuals, known as identity-by-descent—see, for example, Lange (1997). Bueno Filho and Gilmour (2003) showed that (1/2 )V (ˆ) = , where −1 = r + A−1 − N (k )−1 N , k and r are diagonal matrices of the numbers of plots per block and the numbers of replications per treatment, respectively, = 2 /2 and N is the block-treatment concurrence matrix. Thus, the properties of must be studied to develop design criteria. Without loss of generality we will take 2 = 1 so that, in what follows, variances refer to scaled variances. Any useful optimality criterion for selection experiments should be based on discriminating between good and bad genetic material. Assume that we intend to select s out of t treatments for further experimentation. Let S be the set of the best s treatments. Then we are particularly interested in the comparisons ij = i − j , for i ∈ S and j ∈ / S. To ˆ minimize the probability of a particular mistake in selection, we must minimize P (ij < 0). Given unbiased estimators ˆ ij ). Since any such mistake leads to an incorrect selection, the design of the ’s, this is equivalent to minimizing V ( ˆ criterion should minimize i∈S j ∈S / V (ij ). The same criterion arises from quadratic loss utility on the comparisons of interest, i.e. U (Z, y, ) = −L(Z, y, ), where L(Z, y, ) =
(1)
ˆ ij − ij )2 . (
(2)
i∈S j ∈S /
We should choose Z to minimize L(Z, y, ), but this also depends on the data, y, and on the unknown parameters, . Integrating first over y|, gives ˆ ij ), Ey| {L(Z, y, )} = V ( i∈S j ∈S /
ˆ ij ) does not depend on and ˆ ij is an unbiased estimator of ij . since V ( Next we integrate over the prior distribution of to get ⎧ ⎫ ⎨ ⎬ ˆ ij ) E [Ey| {L(Z, y, )}] = E V ( ⎩ ⎭ i∈S j ∈S / ˆ ij )P (R = S), = V ( R∈ i∈R j ∈R /
where is the set of all possible sets of s treatments. If the prior is such that P (R = S) is constant, i.e. any treatment is equally likely to be among the best s, then E [Ey| {L(Z, y, )}] ∝
t−1 t i=1 j =i+1
This establishes the following result.
ˆ ij ). V (
1448
J.S. de S. Bueno Filho, S.G. Gilmour / Journal of Statistical Planning and Inference 137 (2007) 1446 – 1451
Theorem 1. Given any prior distribution for such that the prior probability that treatment i is among the best s treatments is constant ∀i ∈ {1, . . . , t} and given the utility function defined by Eqs. (1) and (2), maximizing the expected utility is equivalent to minimizing the criterion C(Z|, A) =
t−1 t
V (ˆi − ˆ j ),
(3)
i=1 j =i+1
assuming and A are known. This criterion is the one used by Bueno Filho and Gilmour (2003). In fact this result is even more general, since s need not be fixed and need not represent the number of treatments we will select. Relabel the above utility function Us (Z, y, ). Then our utility function can be any weighted average of Us (Z, y, ), for s=1, . . . , t. This is quite important in practice, since it is not usually the case that the number of treatments that will be selected from the experiment is fixed in advance of the experiment. Rather, an appropriate number to ensure enough variation in the selected population will be selected. This criterion is also not specific to the problem of selection, although that is our primary objective. It would also be a sensible utility function for ranking treatments and for estimating treatment effects. In fact, it corresponds to the commonly used quadratic loss utility function for estimating treatment differences. 3. Prior distribution for heritability Consider the following designs for comparing 6 treatments in 4 blocks of size 3: D1 D2 D3 D4 D5
: {1, 2, 3}, {1, 5, 6}, {2, 4, 6}, {3, 4, 5}, : {1, 3, 6}, {1, 4, 5}, {2, 3, 5}, {2, 4, 6}, : {1, 2, 5}, {1, 3, 4}, {2, 5, 6}, {3, 4, 6}, : {1, 3, 6}, {1, 4, 5}, {2, 3, 6}, {2, 4, 5}, : {1, 2, 3}, {1, 2, 3}, {4, 5, 6}, {4, 5, 6},
where Di denotes the design and braces denote blocks. If we consider the treatment labels to be exchangeable, the regular graph designs D1 and D2 are isomorphic, as are the partially balanced incomplete block designs D3 and D4 . The treatment concurrence matrices and the graphs for these designs were given by Bueno Filho and Gilmour (2003) who plotted the values of C(Z|, A) for each design for a range of different values and for various family structures, including one defined by ⎡1 ⎤ 0 0 0 18 81 2 ⎢0 1 0 1 0 1 ⎥ ⎢ 2 8 8⎥ ⎢ ⎥ ⎢0 0 1 1 1 0⎥ ⎢ ⎥ 2 8 8 A=⎢ ⎥. ⎢0 1 1 1 0 0⎥ 8 8 2 ⎢ ⎥ ⎢1 ⎥ ⎣ 8 0 18 0 21 0 ⎦ 1 8
1 8
0
0
0
1 2
With such structures, the treatment labels can no longer be considered to be exchangeable and so the five designs become genuinely different. The model with fixed treatment effects can be obtained by making → 0 (i.e. 2 → ∞), so the best designs in this situation are also the best designs for genetically unrelated treatments. For this situation D5 is the worst possible, being disconnected for fixed treatments. On the other hand, D5 is a good design for greater than approximately 0.4. This surprising result arises because the analysis with random treatment effects allows information on comparisons involving a particular treatment to be “borrowed” from comparisons involving related treatments. The use of point prior values for does not account for uncertainty in its value, and we may be tempted to identify a design as the best when it is probably not. Instead we can define a prior distribution on and find the optimal design
J.S. de S. Bueno Filho, S.G. Gilmour / Journal of Statistical Planning and Inference 137 (2007) 1446 – 1451
1449
averaged over this prior. The most usual way to describe the ratio of variances between treatments and experimental units in plant breeding is through the heritability, h2 = 2 /(2 + 2 ) = 1/(1 + ). This is a measure of the covariance between relatives as a proportion of the total phenotypic variance in experimental units and 0 < h2 < 1. Heritability values can be given prior probability distributions and prior elicitation is relatively straightforward because plant breeders have estimates of heritability from large amounts of previous data. Assuming some prior probability distribution on h2 , which implies a prior for it is possible to use our Bayesian criterion averaged over this prior. Thus, our utility function becomes E,h2 {U (Z, y, , h2 )} = Eh2 [E|h2 {U (Z, y, |h2 )}] = − Eh2 {C(Z|h2 , A)} 1 C(Z|h2 , A)p(h2 ) dh2 . = − 0
In general the evaluation of this utility function for any design might require Monte Carlo sampling from the prior distribution for h2 . However, in this small example it is possible to get analytical solutions through some lengthy, but straightforward, algebraic manipulations. A Beta prior distribution for h2 will be flexible enough in many practical situations and is easily interpreted. The Beta distribution with parameters and has probability density function 1 p(h2 ) = (h2 )( −1) (1 − h2 )( −1) /B( , ), where B( , ) = 0 u −1 (1 − u) −1 du. It has expectation /( + ) and, when used as a prior distribution, + represents the equivalent sample size of the prior information. The Beta distribution with parameters and is implied if we have independent priors for 2 and 2 which are, respectively, Inverse- 22 and Inverse- 22 (which also implies that / ∼ F2 ,2 ), although it is not necessary to assume this for the Beta prior for h2 to arise. The steps in a numerical evaluation of the criterion values are as follows:
1. for a known A matrix, write out the design matrix and the information matrix as a function of and hence as a function of h2 ; 2. write the criterion conditional on h2 , according to Eq. (3); 3. assign a Beta prior for h2 with and chosen to represent the prior mean heritability by /( + ) and the equivalent sample size by + ; 4. evaluate the integrated criterion value, taking the prior distribution of h2 as weights, 1 2 E,h2 {U (Z, y, , h )} = − C(Z|h2 , A)p(h2 ) dh2 . 0
In the next section, we present the results of following these steps to evaluate the Bayesian design criterion for each of the 5 designs considered, averaged over Beta priors for the heritabilities. 4. Results Tables 1 and 2 show some values of the Bayesian criterion calculated for 36 different Beta priors for the heritabilities. Moving down Table 1, the priors become more informative, while Table 2 represents priors with no uncertainty, and moving from left to right in each table the expected heritability becomes greater. The priors studied have many different shapes, including those that are highly skewed to the left and right, those that are bimodal, one which is uniform on [0, 1], those which tend towards normality, point priors and various others in between. It can be seen that for a wide range of flat priors the best design is in the class of regular graph designs (in fact D2 is clearly best). However, when the prior distribution gives a high probability to low heritabilities, we see that D5 can be the best. The tables show that in the case of little prior knowledge, it is safest to use a standard regular graph design. However, the prior knowledge that the heritability is low does not have to be very great before D5 becomes better, if only slightly better. With the prior having information equivalent to a sample size of 2, for example, D5 is the best for an average heritability of 0.1. For strong prior knowledge on heritability (equivalent information to an infinite sample size) regular graph designs are better just in situations with very high heritabilities (h2 > 0.7), situations that are unusual for many important economic traits.
1450
J.S. de S. Bueno Filho, S.G. Gilmour / Journal of Statistical Planning and Inference 137 (2007) 1446 – 1451
Table 1 Values of the design criterion and best designs (shown in bold) for Beta priors with different and Design
Beta parameters ( , ) (0.1,0.9)
(0.3,0.7)
(0.5,0.5)
(0.7,0.3)
(0.9,0.1)
D1 D2 D3 D4 D5
1.55 1.54 1.56 1.56 1.58
4.79 4.77 4.86 4.85 5.05
8.24 8.22 8.44 8.42 9.09
11.93 11.91 12.35 12.33 13.99
15.90 15.89 16.67 16.66 20.20
D1 D2 D3 D4 D5
(0.2,1.8) 1.48 1.47 1.48 1.48 1.46
(0.6,1.4) 4.61 4.59 4.63 4.62 4.61
(1,1) 7.99 7.96 8.1 8.08 8.3
(1.4,0.6) 11.69 11.66 12.01 11.98 12.98
(1.8,0.2) 15.78 15.77 16.49 16.47 19.53
D1 D2 D3 D4 D5
(0.5,4.5) 1.42 1.42 1.41 1.41 1.39
(1.5,3.5) 4.45 4.42 4.42 4.42 4.32
(2.5,2.5) 7.76 7.71 7.77 7.76 7.66
(3.5,1.5) 11.44 11.40 11.64 11.61 11.94
(4.5,0.5) 15.64 15.63 16.28 16.25 18.61
D1 D2 D3 D4 D5
(1,9) 1.40 1.39 1.39 1.39 1.37
(3,7) 4.38 4.35 4.34 4.34 4.24
(5,5) 7.65 7.6 7.63 7.62 7.44
(7,3) 11.32 11.28 11.46 11.43 11.49
(9,1) 15.57 15.55 16.16 16.13 18.03
D1 D2 D3 D4 D5
(5,45) 1.38 1.38 1.37 1.37 1.36
(15,35) 4.32 4.29 4.28 4.27 4.18
(25,25) 7.56 7.51 7.51 7.49 7.28
(35,15) 11.21 11.16 11.29 11.26 11.09
(45,5) 15.49 15.47 16.03 15.99 17.32
D1 D2 D3 D4 D5
(10,90) 1.38 1.37 1.37 1.37 1.36
(30,70) 4.31 4.28 4.27 4.27 4.17
(50,50) 7.55 7.49 7.49 7.48 7.26
(70,30) 11.19 11.14 11.27 11.23 11.05
(90,10) 15.48 15.46 16.01 15.97 17.21
Table 2 Values of the design criterion and best designs (shown in bold) for point priors with different h2 Design
D1 D2 D3 D4 D5
h2 0.1
0.3
0.5
0.7
0.9
1
1.38 1.37 1.37 1.37 1.36
4.31 4.28 4.26 4.26 4.16
7.53 7.48 7.47 7.46 7.24
11.18 11.12 11.24 11.21 11.00
15.46 15.45 16.00 15.95 17.09
18.00 18.00 19.00 19.00 24.00
J.S. de S. Bueno Filho, S.G. Gilmour / Journal of Statistical Planning and Inference 137 (2007) 1446 – 1451
1451
5. Final remarks For common situations of selection in plant breeding we are very likely to have either a sparse genetical relationship matrix, A, or discrete family structures. For some of these situations it will be possible to find optimal designs robust to assumptions on (assumptions on heritabilities). However, for other situations (complex pedigrees) the design choice may change dramatically at particular values using point prior estimates. Our Bayesian optimality design criterion results from a general framework that encompasses all cases, and involves just sampling from some reasonable distribution on heritabilities. For flat priors on heritabilities we have found that the good designs will usually be in the class of regular graph designs. Our results use utility functions and prior distributions to find optimal designs for likelihood-based analysis. It could be argued that we should use our prior information in the analysis of data. This is quite possible, of course, but is a different problem. If the results are to be made public, then the experiment has to be designed with a view to analyzing the data using a population of prior distributions, not just the experimenter’s. Our method makes use of the experimenter’s prior knowledge for the design, but recognizes that in practice the data will be analyzed by likelihood methods. Bayesian design for likelihood analysis is reviewed by Chaloner and Verdinelli (1996), who suggest that it should be used much more often in practice. Acknowledgements This work was conceived when the first author visited Queen Mary, University of London, with financial support from CAPES grant BEX-1248/00-6. It was continued when the second author visited Universidade Federal de Lavras, with financial support from the Engineering and Physical Sciences Research Council grant GR/S14009/01. We thank the referees for comments which helped us to improve the presentation of this work. References Bernardo, R., 1996. Best linear unbiased prediction of maize single-cross performance. Crop Sci. 36, 50–56. Bueno Filho, J.S. de S., Gilmour, S.G., 2003. Planning incomplete block experiments when treatments are genetically related. Biometrics 59, 375–381. Chaloner, K., Verdinelli, I., 1996. Bayesian experimental design: a review. Statist. Sci. 10, 270–304. Henderson, C.R., 1975. Best linear unbiased estimation and prediction under a selection model. Biometrics 31, 423–449. Henderson, C.R., 1985. Best linear unbiased prediction of nonadditive genetic merits in noninbred populations. J. Animal Sci. 60, 111–117. Lange, K., 1997. Mathematical and Statistical Methods for Genetic Analysis. Springer, New York. Patterson, H.D., Thompson, R., 1971. Recovery of inter-block information when block sizes are unequal. Biometrika 58, 545–554.