- Email: [email protected]
Optimal designs for second-degree Kronecker model mixture experiments
Journal of Statistical Planning and Inference 123 (2004) 117 – 131
www.elsevier.com/locate/jspi
Optimal designs for second-degree Kronecker model mixture experiments Thomas Klein∗ Institut fur Mathematik, Universitat Augsburg, D-86135 Augsburg, Germany Received 11 August 2001; accepted 26 January 2003
Abstract The paper investigates optimal designs in the second-degree Kronecker model for mixture experiments. Three groups of novel results are presented: (i) characterization of feasible weighted centroid designs for a maximum parameter system, (ii) derivations of D-, A-, and E-optimal weighted centroid designs, and (iii) numerically p -optimal weighted centroid designs. Results on a quadratic subspace of invariant symmetric matrices containing the information matrices involved in the design problem serve as a main tool throughout the analysis. c 2003 Elsevier B.V. All rights reserved. MSC: 62K05 Keywords: Designs for mixtures; Design optimality; D-, A-, E-optimality; Kiefer’s p -criteria; Weighted centroid designs; Second-order Kronecker model; Quadratic subspace; Invariant symmetric matrices
1. Introduction A mixture experiment is an experiment in which the m factors t1 ; : : : ; tm are nonnegative and subject to the simplex restriction t1 + · · · + tm = 1, that is, the factors represent relative proportions of m ingredients blended in a mixture. Thus, the experimental conditions of a mixture m experiment are points in the probability simplex Tm = {t = (t1 ; : : : ; tm ) ∈ [0; 1]m : i=1 ti = 1}. Cornell (1990) lists numerous examples and applications of mixture experiments. A particular polynomial regression model for mixture experiments suggested by Draper and Pukelsheim (1998) is the second-degree Kronecker model 2
f : Tm → Rm ; ∗
t = (t1 ; : : : ; tm ) → t ⊗ t = (ti tj ) i; j=1; :::; m
ordered lexicographically
Tel.: +49-821-5982202; fax: +49-821-5982280. E-mail address: [email protected] (T. Klein).
c 2003 Elsevier B.V. All rights reserved. 0378-3758/$ - see front matter doi:10.1016/S0378-3758(03)00145-9
;
118
T. Klein / Journal of Statistical Planning and Inference 123 (2004) 117 – 131
with model equation E[Yt ] = f(t) =
m
ii ti2 +
i=1
m
( ij + ji )ti tj :
(1)
i; j=1 i¡j
Here Yt , the response under experimental condition t ∈ Tm , is taken to be a real-valued 2 random variable, while = ( 11 ; 12 ; : : : ; mm ) ∈ Rm is the vector of unknown parameters. All observations taken in an experiment are assumed to be uncorrelated and to have common unknown variance 2 ∈ (0; ∞). The paper determines optimal designs for a maximum subsystem of parameters in the Kronecker model, where Kiefer’s p -functions serve as optimality criteria. These results complement Galil’s and Kiefer’s (1977) paper on p -optimal designs in an alternative second-degree mixture model, the Sche78e model E[Yt ] = fS (t) =
n i=1
i ti +
m
ij ti tj
(2)
i; j=1 i¡j m+1
with parameter vector = (1 ; : : : ; m ; 12 ; 13 ; : : : ; m−1; m ) ∈ R( 2 ) , introduced by ScheGHe (1958, 1963). Both the Kronecker and the ScheGHe model are based on the same space of regression polynomials, but diGer in their choice of representing this space. Draper and Pukelsheim (1998) and Prescott et al. (2002) put forward several advantages of the Kronecker model, e.g. the homogeneity of regression terms and reduced ill-conditioning of information matrices. Interestingly, both models share the same invariance properties, see the appendix. An outline of the paper is as follows: Section 2 formulates the parameter system K of interest and the associated design problem. Based on the completeness Theorem 2.4 in Draper et al. (2000), the design problem is reduced to the class of weighted centroid designs. In Section 3 we present a characterization of weighted centroid designs under which K is estimable. In Section 4 we derive D-, A-, and E-optimal weighted centroid designs, while Section 5 presents numerical results on p -optimal weighted centroid designs. 2. The design problem 2
Since the Kronecker model’s full parameter vector ∈ Rm is not estimable, we consider a maximum parameter subsystem instead, that is, a parameter system K where the range R(K) coincides with the span of the regression range X := {f(t) : t ∈ Tm }. This formalizes the idea of estimating as many parameters as possible, see Klein (2001). A canonical maximum parameter subsystem is as follows. We denote the canonical unit vectors in Rm by e1 ; : : : ; em and set eij := ei ⊗ej for i; j= m 1; : : : ; m. The canonical unit vectors in R( 2 ) are denoted by Eij with 1 6 i ¡ j 6 m in m2 ×s lexicographic order. We write s := ( m+1 2 ), and deIne the matrix K := (K1 ; K2 ) ∈ R
T. Klein / Journal of Statistical Planning and Inference 123 (2004) 117 – 131
119
by K1 :=
m i=1
eii ei ;
K2 :=
m 1 (eij + eji )Eij : 2 i; j=1
(3)
i¡j
Then the matrix K satisIes R(K) = span X, and the maximum parameter subsystem considered in the following is: ( ii )16i6m 2 ∈ Rs for all ∈ Rm : K = 1 ( 2 ( ij + ji ))16i¡j6m An experimental design for a mixture experiment is a probability measure with Inite support on Tm . The statistical properties of a design are reJected by its moment matrix M () := Tm f(t)f(t) d. Here, the matrix M () contains the fourth moments of . Let Is denote the identity matrix in Rs×s . The information matrix 2
CK (M ()) := min{LM ()L |L ∈ Rs×m ;
LK = Is };
captures the amount of information the design contains on K , see Pukelsheim (1993, Section 3.10). As a consequence of K being a maximum parameter system, information matrices for K take a particularly simple form. The following lemma quotes this result from Klein (2001). Lemma 2.1. De9ne the matrix L˜ := (K K)−1 K = (K1 ; 2K2 ) ∈ Rs×m with K1 ; K2 from ˜ L˜ for all moment matrices M , that is, information (3). Then we have CK (M ) = LM matrices for K are linear transformations of moment matrices. 2
A scalar measure for the amount of information inherent to CK (M ()) is provided by Kiefer’s p -functions or matrix means p (C) = ((1=s)traceC p )1=p for all p ∈ [ − ∞; 1] and C from the set PD(s) of positive deInite s × s matrices. The family of matrix means includes the popular T-, D-, A-, and E-criteria, corresponding to p = 1, 0, −1, and −∞, respectively. The problem of Inding a design with maximum information on the parameter subsystem K can now be formulated as Maximize
p (CK (M ())) with ∈ T
subject to
CK (M ()) ∈ PD(s);
(4)
where T denotes the set of all designs on Tm . The side condition CK (M ()) ∈ PD(s) is equivalent to the existence of an unbiased linear estimator for K under , see Theorem 3.15 in Pukelsheim (1993). In this case, the design is called feasible for K . Any design solving problem (4) for a Ixed p ∈ [ − ∞; 1] is called p -optimal for K in T. Due to a result of Draper et al. (2000) on weighted centroid designs, the problem can be substantially reduced. Denition 2.2. For m ¿ 2 and j ∈ {1; : : : ; m}, the jth elementary centroid design j is the uniform distribution on the centroids of depth j, that is, on all points taking the
120
T. Klein / Journal of Statistical Planning and Inference 123 (2004) 117 – 131
j form 1=j i=1 meki ∈ Tm with 1 6 k1 ¡ k2 ¡ · · · ¡ kj−1 ¡ kj 6 m. A convex combination ( ) := j=1 j j with = ( 1 ; : : : ; m ) ∈ Tm is called a weighted centroid design with weight vector . We denote the set of all weighted centroid designs by (Tm ). Weighted centroid designs are exchangeable, that is, they are invariant under permutation of the mixture ingredients. This property leads to H-invariance of information matrices, HCK (M ())H = CK (M ())
for all H ∈ H;
∈ (Tm );
(5)
where H is the matrix group deIned in (A.1). The appendix summarizes some results on H-invariant information matrices which we utilize in the following sections. Theorems 6.4 and 7.4 in Draper and Pukelsheim (1999) and Theorem 2.4 in Draper et al. (2000) show that the set (Tm ) is an essentially complete class of designs relative to the Kiefer ordering of moment matrices in the Kronecker model. See Pukelsheim (1993, Chapter 14) for an introduction to the Kiefer ordering. Since the group H in (5) is a subgroup of the orthogonal group, all matrix means p are H-invariant and the target function of design problem (4) is isotonic relative to the Kiefer ordering. Therefore, the set (Tm ) is also essentially complete relative to the target function, that is, for any given design ∈ T there is a weighted centroid design () satisfying (p ◦ CK ◦ M )() ¿ (p ◦ CK ◦ M )() for all p ∈ [ − ∞; 1]. As a consequence, we may restrict the set of competing designs in problem (4) to (Tm ), thus obtaining a mere allocation problem for the weight vector ∈ Tm .
3. Feasible weighted centroid designs First we Ind an explicit description of feasible weighted centroid designs, that is, weighted centroid designs ( ) whose information matrix for K is regular. Let Cj := CK (M (j )) denote the information matrix of the elementary centroid design j , for 1 6 j 6 m. From the linearity of the information matrix mapping CK in Lemma 2.1 we get, for every ∈ Tm , CK (M (( ))) = j Cj with J( ) := {j = 1; : : : ; m | j ¿ 0}; (6) j∈J( )
whence R(CK (M (( )))) =
R(Cj ):
(7)
j∈J( )
Eq. (7) suggests studying the ranges of the information matrices C1 ; : : : ; Cm of the elementary centroid designs. This is done in Lemma A.4. By utilizing this result, the following theorem characterizes feasible weighted centroid designs for K . The essential ingredient is the indicator set J( ) from (6), collecting those subscripts j whose weight j is nonzero.
T. Klein / Journal of Statistical Planning and Inference 123 (2004) 117 – 131
121
Theorem 3.1. Let = ( 1 ; : : : ; m ) ∈ Tm be a weight vector of a weighted centroid design ( ) and set J( ) := {j | j ¿ 0} as in (6). Then the following assertions hold: (i) For m ∈ {2; 3}, the weighted centroid design ( ) is feasible for K if and only if J( ) ⊇ {1; 2}. (ii) For m ¿ 4, the weighted centroid design ( ) is feasible for K if and only if |J( ) ∩ {1; : : : ; m − 2}| ¿ 2
or
|J( ) ∩ {2; : : : ; m − 1}| ¿ 2:
Proof. We only prove (ii). Using Lemma A.4 and applying the dimension formula to the right-hand side of (7), we see that CK (M (( ))) is regular if and only if the set J( ) satisIes the condition given in (ii). For four or more ingredients, m ¿ 4, we combine our Theorem 3.1 with Theorem 3.2 from Draper et al. (2000). Corollary 3.2. Suppose m ¿ 4 and p ∈ [ − ∞; 1]. Then there is a weighted centroid design ( ) which is p -optimal for K in T and whose weight vector ∈ Tm satis9es, for some j ∈ {2; : : : ; m − 2}, one of the six conditions (i) J( ) = {1; j};
(iii) J( ) = {1; j; j + 1};
(v) J( ) = {j; j + 1; m};
(ii) J( ) = {j; j+1};
(iv) J( ) = {1; j; m};
(vi) J( ) = {1; j; j+1; m}:
4. Optimal weighted centroid designs We begin the investigation of p -optimal weighted centroid designs for K by adapting a general Equivalence Theorem to our speciIc problem. Theorem 4.1. Let ∈ Tm be the weight vector of a weighted centroid design ( ) which is feasible for K , and let J( ) be the set of active indices from (6). Furthermore, let C := CK (M (( ))) and p ∈ (−∞; 1]. Then ( ) is p -optimal for K in T if and only if = trace C p for all j ∈ J( ); p−1 trace Cj C 6 trace C p otherwise: Proof. We start from Theorem 7.19 in Pukelsheim (1993), that ( ) is p -optimal for K if and only if there is a generalized inverse G of M := M (( )) with trace M (())GKC p+1 K G 6 trace C p
for all ∈ Tm :
(8)
122
T. Klein / Journal of Statistical Planning and Inference 123 (2004) 117 – 131
With C = (K K)−1 K MK(K K)−1 , M = K(K K)−1 K M , and M (()) = KCK (M (( )))K from Lemma 2.1, the left-hand side is trace M (())GKC p+1 K G = trace (K GMK(K K)−1 ) CK (M (()))(K GMK(K K)−1 )C p−1 :
(9)
Due to the feasibility of ( ) we have R(K) = R(M ). Hence K = MZ for some Z, and so K GMK(K K)−1 = Z MGMK(K K)−1 = Z MK(K K)−1 = I( m ) . Now the right-hand 2 side of (9) simpliIes to trace CK (M (()))C p−1 , and (8) turns into trace CK (M (())) C p−1 6 trace C p for all ∈ Tm . According to equation (6) we can write the left-hand m p−1 , giving trace Cj C p−1 6 trace C p for all 1 6 j 6 m. Fiside as j=1 j trace Cj C nally, equality must hold for any j ∈ J( ). The proof is complete. The conditions for p -optimality provided above can be explicitly solved for the weight vector ∈ Tm in the special cases of D- and A-optimality (that is, p = 0 and p = −1) in Theorems 4.3 and 4.4, involving only weighted centroid designs with the Irst and second weight positive. The following theorem shows that this property implies uniqueness of the optimizers. Theorem 4.2. Let p ∈ (−∞; 1) and ( ) be a weighted centroid design that is p optimal for K in T , with ∈ Tm . (i) If J( ) = {1; 2}, then ( ) is the unique solution of problem (4). (ii) If J( ) = {1; 2; 3}, then there is no further exchangeable design Q ∈ T that is p -optimal for K in T . If there is a non-exchangeable p -optimal design for K , then all its support points are centroids of depths 1; 2, or 3. Proof. We only prove the Irst assertion. The proof of (ii) parallels that of (i) but excludes step (IV) below. In steps (II) and (III) of the proof we utilize some facts on fourth moments of exchangeable designs taken from Draper et al. (2000). That 4 is, an exchangeable design Q ∈ T possesses Ive distinct fourth moments & ( ) Q = t1 d , Q 4 Q &22 () Q = t12 t22 d , Q &211 () Q = t12 t2 t3 d , Q and &1111 () Q = t1 t2 t3 t4 d , Q &31 () Q = t13 t2 d , where the moments &211 () Q and &1111 () Q only occur for m ¿ 3 and m ¿ 4. The identity &31 () Q = &22 () Q is equivalent to Q being a weighted centroid design. The fourth moments of a weighted centroid design () are convex combinations of the fourth m moments of the elementary centroid designs, e.g. &4 (()) = j=1 j &4 (j ) for all =(1 ; : : : ; m ) ∈ Tm . Draper et al. (2000, p. 581) list the fourth moments of 1 ; : : : ; m . (I) First we show that M := M (( )) is the unique p -optimal moment matrix for K in M (T). Theorem 8.14 in Pukelsheim (1993) states that the moment matrix M () of any solution ∈ T to design problem (4) with p ∈ (−∞; 1) satisIes MGK = K;
(10)
where G is a generalized inverse of M (( )) satisfying optimality condition (8) for ( ). From our Theorem 4.1 it follows that any generalized inverse enjoys this property. We choose the Moore–Penrose inverse M + = K(K K)−1 C −1 (K K)−1 K with C :=
T. Klein / Journal of Statistical Planning and Inference 123 (2004) 117 – 131
123
CK (M (( ))). Inserting M + and using M () = KCK (M ())K , equation (10) becomes K = KCK (M ())C −1 . Post-multiplication with CK yields M = M (). (II) Suppose () with ∈ Tm is a further weighted centroid design which is p -optimal for K in (Tm ). According to part (I) of the proof, M = M (()) holds, that is, ( ) and () have the same fourth moments. We thus obtain m 1 1 2 = &4 (()) = j ; + 3 8m j m m j=1
m
j−1 2 = &31 (()) = j ; 8m(m − 1) j 3 m(m − 1) j=2
0 = &211 (()) =
m j=3
0 = &1111 (()) =
(j − 1)(j − 2) j ; j 3 m(m − 1)(m − 2)
m j=4
(j − 1)(j − 2)(j − 3) j : j 3 m(m − 1)(m − 2)(m − 3)
The equations for &211 (()) and &1111 (()) entail 3 = · · · = m = 0. Hence, the equations for &4 (()) and &31 (()) imply = . (III) Let Q ∈ T be an exchangeable design that is p -optimal for K in T. Then the Q Q Thus, Q is a weighted fourth moments of Q and ( ) coincide, implying &31 ()=& 22 (). centroid design, and so Q = ( ) as shown in (II). (IV) Now we drop exchangeability of the design ∈ T, but assume supp ⊆ supp ( ). From (I) we have CK (M ()) = C. Due to rank M (()) = |supp ( )| = s, the set {f(t) | t ∈ supp ( )} is linearly independent. Theorem 8.7 in Pukelsheim (1993) yields f ({f(t)}) = ( )f ({f(t)}) for all t ∈ supp ( ), and so = ( ). (V) Finally, let ∈ T be an arbitrary design which is p -optimal for K in T. Then, the exchangeable design Q := 1=m! R∈Perm(m) R obtained by averaging over the group Perm(m) is also optimal, and supp ⊆ supp . Q According to part (III), =() Q holds, which implies supp ⊆ supp ( ). Hence, part (IV) establishes = ( ). We may now derive optimal weighted centroid designs for the T-, D-, A- and E-criteria (that is, p = 1; 0; −1; −∞). Since the T-criterion is linear, the elementary centroid design 1 is the unique maximizer of 1 ◦ CK ◦ M in (Tm ), with optimum value (T = 1=s. However, 1 is not feasible for K due to Theorem 3.1. The cases of D- and A-optimality are covered by the following theorems. Since the D-criterion is invariant under reparameterizations, Theorem 4.3 applies to the Kronecker model as well as to the ScheGHe model. Theorem 4.3. In the second-degree mixture model with m ¿ 2 ingredients, the unique D-optimal design for K is 2 m−1 1 + 2 : ( (D) ) = m+1 m+1 The maximum of the D-criterion is (D = (1=s) 2−2(m−1)=(m+1) .
124
T. Klein / Journal of Statistical Planning and Inference 123 (2004) 117 – 131
Proof. Within the ScheGHe model (2), the design ( (D) ) is D-optimal for the full parameter vector, see Kiefer (1961). Since the D-criterion is invariant under linear re-parameterization of the space of regression polynomials (see e.g. GaGke, 1981), this design is also D-optimal for K within the Kronecker model. Uniqueness follows from Theorem 4.2. Theorem 4.4. In the second-degree Kronecker model with m ¿ 2 ingredients, the unique A-optimal design for K is √ √ 2(m − 1)[2(m − 1) − m + 3] 2(m − 1) m + 3 − m − 3 (A) ( ) = 1 + 2 : 4m2 − 9m + 1 4m2 − 9m + 1 The maximum of the A-criterion is (A =
√ 2(m + 1)[4m2 − 7m + 7 − 4(m − 1) m + 3] : m(4m2 − 9m + 1)2
Proof. Uniqueness of the proposed optimal design again follows from Theorem 4.2. The actual optimality proof consists of two parts: First, deriving an optimality candidate, and second, verifying its optimality. We concentrate on the Irst part. Let = ( 1 ; 2 ; 0; : : : ; 0) ∈ Tm be a weight vector with J( ) = {1; 2} and suppose ( ) is A-optimal for K in T. With the abbreviation C( ) := CK (M (( ))), Theorem 4.1 thus implies = trace C( )−1 for j ∈ {1; 2}; −2 (11) trace Cj C( ) 6 trace C( )−1 otherwise: The proof of Lemma A.4 already demonstrated how to compute an inverse matrix in the space Sym (s; H) by solving a system of linear equations. By the same approach we obtain the blocks of C( )−1 in the partitioning suggested by Lemma A.1, namely m m Im ; (C( )−1 )21 = − V1 ; (C( )−1 )11 = 1 2 1 (C( )−1 )22 =
m[4(m − 1) 1 + 2 ] m I( m ) + W2 : 2 2 1 2 4 1
Lemma A.2 yields the blocks of Cj C( )−2 for j = 1; : : : ; m. Using trace U2 = trace W2 = trace W3 = 0, we can now evaluate both sides of the optimality condition (11), trace C( )−1 =
m2 [4(m − 1)2 1 + (m + 3) 2 ] ; 4 1 2
trace Cj C( )−2 =
(12)
m2 32(m − 1)2 (j − 1) 12 4j 3 12 22
+ 16(m − 1)(j − 1)(j − 2) 1 2 + (j − 2)2 (m + j + 2) 22 :
(13)
T. Klein / Journal of Statistical Planning and Inference 123 (2004) 117 – 131
125
Equating (12) and (13) for j = 1 and j = 2 and inserting 2 = 1 − 1 , we obtain a quadratic equation for 1 . Its unique solution in (0; 1) is the weight 1(A) given in the assertion. By construction, the weight vector (A) = ( 1(A) ; 1 − 1(A) ; 0; : : : ; 0) satisIes the two equations in condition (11). The remaining m − 2 inequalities still need to be established. To this end, we compute trace Cj C( (A) )−2 − trace Cj+1 C( (A) )−2 m2 (j − 1) m(j 3 + 11j 2 + 12j + 4) = 3 4j (j + 1)3 12(A) √ + 8 m + 3(j 3 − 3j 2 − 5j − 2) − 2(j 3 − 17j 2 − 22j − 8)
for j = 2; : : : ; m − 1. Both terms on the right-hand side are found to be positive, whence trace C2 C( (A) )−2 ¿ · · · ¿ trace Cm C( (A) )−2 . Therefore, ( (A) ) is indeed A-optimal for K in T. The approach of the above proof, that is, solving the optimality condition from Theorem 4.1 for the weights 1 ; 2 , does not carry over to E-optimality. In fact, there are only partial results on E-optimal weighted centroid designs, quoted from Klein (2001) in the following two lemmas. Lemma 4.5. In the second-degree Kronecker model with m = 2 or m = 4 ingredients, the weighted centroid designs 5 6 for m = 2; 11 1 + 11 2 (E) ( ) = 13 12 9 for m = 4 46 1 + 23 2 + 46 3 are E-optimal for K in T, and there is no further E-optimal weighted centroid 1 1 design for K . The maxima of the E-criterion are (E = 11 for m = 2, and (E = 46 for m = 4. Lemma 4.6. In the second-degree Kronecker model for mixture experiments with m ¿ 5 ingredients, there is no weighted centroid design ( ) with J( ) ⊆ {1; 2; 3} which is E-optimal for K in T. As is seen in the following section, all numerically p -optimal weight vectors (p) , with p ∈ [ − 50; 1], satisfy J( (p) ) ⊆ {1; 2; 3}. However Lemma 4.6 says that this condition does not hold for E-optimal weight vectors when m ¿ 5. This is in line with Galil’s and Kiefer’s (1977) numerical results that p -optimal weight vectors for the full parameter vector of the ScheGHe model use the Irst three plus a further weight. 5. Numerical results The numerical computation of p -optimal weighted centroid designs for K divides into two steps, that is, Inding a local optimality candidate and examining a global
126
T. Klein / Journal of Statistical Planning and Inference 123 (2004) 117 – 131
Fig. 1. p -Optimal weights for K , with m = 4 and m = 10.
optimality condition for this candidate. As the Irst step requires greater eGort when the target function is not diGerentiable we omit the case of E-optimality. In order to Ind an optimality candidate we use Lemma A.3 to compute the eigenvalues of CK (M (( ))) as functions of ∈ Tm . As gradient methods for maximizing the target function of design problem (4) only work on open sets, each of the 6(m − 3) cases from Corollary 3.2 is treated separately. For the global optimality part we examine the condition given in Theorem 4.1, where the spectral decomposition of information matrices from Lemma A.3 is utilized once more. All p -optimal weighted centroid designs with p within the investigated range [ − 50; 1] are covered by the two cases J( ) = {1; 2} and J( ) = {1; 2; 3}. Since we are able to explicitly prove or disprove optimality of a candidate, the numerical algorithm can treat these two cases Irst and does not have to run through the whole list of cases mentioned above. Furthermore, Theorem 4.2 provides uniqueness of optimal designs ( ) with J( ) = {1; 2}; {1; 2; 3} within the set (Tm ). This justiIes the notation (p) for the p -optimal weight vector.
T. Klein / Journal of Statistical Planning and Inference 123 (2004) 117 – 131
(p)
(p)
127
(p)
Fig. 2. p -Optimal weight vectors in barycentric coordinates. The curve p → ( 1 ; 2 ; 3 ) in barycentric coordinates. For p ∈ (p2 ; 1] the curve runs along the 1 - 2 -edge of the simplex, for p ¡ p2 it extends into the interior. Open dots mark the optimal weight vectors for the T-, D-, A-, E-, and −50 -criteria.
Fig. 1 exempliIes our numerical results by listing the optimal weights 1(p) , 2(p) , 3(p) and the corresponding optimal values (p := (p ◦ CK ◦ M ◦ )( (p) ) for selected values of p ∈ [−∞; 1] and for m ∈ {4; 10}. The qualitative behavior of the functions p → i(p) with i = 1; 2; 3 and m ¿ 3 is captured by three critical values 0 ¿ p1 ¿ p2 ¿ p3 . For p ∈ (p2 ; 1) only two positive weights 1(p) and 2(p) are present, while for p ¡ p2 we also have 3(p) ¿ 0. For p descending from 1 to p2 , the weight 1(p) decreases when m 6 4. In the case m ¿ 5 a local minimum of 1(p) is found at p=p1 . Letting p descend from p2 , the weight 2(p) decreases while 3(p) increases. Interestingly, the Irst weight decreases but remains nearly constant. For m ¿ 7 a further change in monotonicity is observed at p = p3 where 1(p) starts a slight increase. Fig. 2 displays the curve p → ( 1(p) ; 2(p) ; 3(p) ) in barycentric coordinates for m = 3; : : : ; 10. Appendix A. H-Invariant information matrices The probability simplex Tm is invariant under permutation, that is, renumbering the m ingredients. Let Sm be the symmetric group of order m, and let R, = i=1 e,(i) ei be the matrix associated with , ∈ Sm . Then the group Perm(m) of m × m permutation matrices acts on Tm according to (R, ; (t1 ; : : : ; tm ) ) → (t,−1 (1) ; : : : ; t,−1 (m) ) and on T according to (R; ) → R = ◦ R−1 . If a design satisIes = R for all R ∈ Perm(m), we call it exchangeable. We deIne the group R, 0 : , ∈ Sm with H := H, = 0 S,
128
T. Klein / Journal of Statistical Planning and Inference 123 (2004) 117 – 131
S, :=
m i; j=1 i¡j
E(,(i); ,( j))↑ Eij ∈ Perm
m
for all , ∈ Sm
2
(A.1)
where (,(i); ,(j))↑ denotes the pair of indices ,(i), ,(j) in ascending order. The group H is a subgroup of Perm(s), and is isomorphic to Sm . As for the Kronecker model, the group H acts on an information matrix C for K through C → HCH . The equivariance property CK (M (R, )) = H, CK (M ())H, for , ∈ Sm , ∈ T links the action of H on CK (M (T)) to the action of Perm(m) on T. Within the ScheGHe model (2), we have MS (R, ) = H, MS ()H, for all moment matrices MS (). Therefore, we end up with the same invariance properties for both the Kronecker and the ScheGHe model. Since H is a subgroup of the orthogonal group, the space Sym (s; H) := {C ∈ Sym (s) : HCH = C for all H ∈ H} of H-invariant symmetric matrices is a quadratic subspace of the space Sym (s) of symmetric s × s matrices, that is, a subspace closed under formation of powers C n with n ∈ N. Basic facts about quadratic subspaces in experimental design can be found in Pukelsheim (1993, Chapters 13 and 14). The particular quadratic subspace Sym (s; H) is analyzed in Klein (2001, 2003). We brieJy summarize the results of this analysis. Lemma A.1. Let 1m denote the vector (1; : : : ; 1) ∈ Rm , and de9ne U2 := 1m 1m − Im ∈ Sym (m) and m m m m Eij (ei + ej ) ; V2 := Eij ek ∈ R( 2 )×m ; V1 := i; j=1 i¡j
W2 :=
W3 :=
i; j=1 k=1 i¡j k∈{i; j}
m
m
i; j=1 i¡j
k;‘=1 k¡‘ |{i; j}∩{k;‘}|=1
m
m
i; j=1 i¡j
k;‘=1 k¡‘ {i; j}∩{k;‘}=?
; Eij Ek‘
Eij Ek‘ ∈ Sym
m
2
:
Then any matrix C ∈ Sym (s; H) can be uniquely represented in the form aIm + bU2 cV1 + dV2 C = cV1 + dV2 eI m + fW2 + gW3 2
with coe=cients a; : : : ; g ∈ R. The terms containing V2 , W2 and W3 only occur for m ¿ 3 or m ¿ 4, respectively. In particular, dim Sym (s; H) = 4; 6; 7
for m = 2; m = 3 and m ¿ 4:
T. Klein / Journal of Statistical Planning and Inference 123 (2004) 117 – 131
129
When computing powers or inverses of matrices in Sym (s; H), products of the matrix blocks introduced in Lemma A.1 have to be evaluated. The following Lemma A.2 lists a multiplication table for Im , U2 , V1 , V2 , I( m ) , W2 , and W3 . 2
m
( 2 )×m , and W2 ; Lemma A.2.
any m ¿ 2, the matrices U2 ∈ Sym (m), V1 ; V2 ∈ R mFor W3 ∈ Sym ( 2 ) from Lemma A.1 satisfy the following equations:
(i) U22 V1 V1 (ii) V1 U2 W2 V1 W3 V1 (iii) V1 V1
= = = = = =
(m − 1)Im + (m − 2)U2 ; (m − 1)Im + U2 ; V1 + 2V2 ; (m − 2)V1 + 2V2 ; (m − 3)V2 ; 2I m + W2 ; 2
V1 V2 V2 V2 V 2 U2 W2 V 2 W 3 V2 V1 V2
= = = = = =
V2 V1 = (m − 2)U2 ;
m−1 Im + m−2 U2 : 2 2 (m − 2)V1 + (m − 3)V2 ; (m − 2)V1 + 2(m − 3)V2 ;
m−2 V1 + m−3 V2 : 2 2 V2 V1 = W2 + 2W3 ;
V2 V2 = (m − 2)I m + (m − 3)W2 + (m − 4)W3 ; W22 W32
2
= 2(m − 2)I m + (m − 2)W2 + 4W3 ; 2
m + m−3 W2 + m−4 W3 ; I = m−2 2 2 2 2
W2 W3 = W3 W2 = (m − 3)W2 + 2(m − 4)W3 : A spectral analysis of H-invariant symmetric matrices can be derived which comprises descriptions of eigenvalues as well as eigenspaces. This expands earlier Indings of Galil and Kiefer (1977), who already provided the eigenvalue part of the analysis. We summarize the results on eigenvalues in the following lemma. Our results on eigenspaces are not used explicitly in this paper, whence the reader is referred to Klein (2001, 2003) for the description of eigenspaces. Lemma A.3. Let a; : : : ; g ∈ R be the coe=cients of a matrix C ∈ Sym (s; H) w.r.t. the basis given in Lemma A.1 with d; f and g occurring only when m ¿ 3 or m ¿ 4, respectively. Furthermore, de9ne
2 D1 := a + (m − 1)b − e − 2(m − 2)f − m−2 g 2 +2(m − 1)[2c + (m − 2)d]2 ; D2 := [a − b − e − (m − 4)f + (m − 3)g]2 + 4(m − 2)(c − d)2 : Then, in the case m ¿ 4, the matrix C has the eigenvalues 51 := e − 2f + g;
√ 1 a + (m − 1)b + e + 2(m − 2)f + m−2 g ± D1 ; 2 2 √ 1 := [a − b + e + (m − 4)f − (m − 3)g ± D2 ] 2
52; 3 := 54; 5
130
T. Klein / Journal of Statistical Planning and Inference 123 (2004) 117 – 131
with multiplicities m(m − 3)=2; 1, and m − 1, respectively. In the case m = 2 only the eigenvalues 52 ; 53 ; 54 occur, whereas for m = 3 there are four eigenvalues 52 ; : : : ; 55 . Information matrices of weighted centroid designs are H-invariant, see equation (5). According to (6), we focus on the information matrices C1 ; : : : ; Cm of elementary centroid designs. For j = 1; : : : ; m we block partition C11; j C21; j Cj = C21; j C22; j and apply Lemma 2.1 to the representation of the moment matrices M (j ) given by Draper et al. (2000). We obtain 1 1 j−1 C11; j = 3 Im + 3 U2 ; j m m−1 j m C21; j =
2 j−1 j−2 2 j−1 V1 + 3 V2 ; j3 m m − 1 j m m−1 m−2
C22; j =
4 j−1 j−2 4 j−1 j−2 j−3 4 j−1 W2 + 3 W3 ; I m + 3 j m m−1 m−2 j m m−1 m−2 m−3 j3 m m − 1 2
where the terms containing V2 ; W2 and W3 only occur for m ¿ 3 and m ¿ 4, respectively. Motivated by equation (7), the following lemma describes the ranges of C1 ; : : : ; Cm , and intersections of these ranges. Lemma A.4. Let 1m denote the vector (1; : : : ; 1) ∈ Tm . Then the following results hold for the information matrices C1 ; : : : ; Cm of elementary centroid designs: (i) R(C1 ) = Rm × {0( m ) }, and rank C1 = m, for all m ¿ 2. 2 D (ii) R(Cj ) = R 2D1; j , and rank Cj = ( m2 ), for m ¿ 4 and j ∈ {2; : : : ; m − 2}, where 2; j
D1; j := (m − 3)[(m − 2)V1 + (j − 2)V2 ]; D2; j := (m − 2)(m − 3)I m + (j − 2)(m − 3)W2 + (j − 2)(j − 3)W3 : (iii) R(Cm−1 ) = R
U2 2V2
(iv) R(Cm ) = span
2
, and rank Cm−1 = m, for all m ¿ 3.
1m 2 1 m 2
, and rank Cm = 1, for all m ¿ 2.
(v) For all m ¿ 2 and i; j ∈ {1; : : : ; m}, i = j, we have {0m } × R(V1 )⊥ for i; j ∈ {2; : : : ; m − 2}; R(Ci ) ∩ R(Cj ) = {0s } otherwise: In particular, dim(R(Ci ) ∩ R(Cj )) =
m (m − 3) for all i; j ∈ {2; : : : ; m − 2}, i = j. 2
T. Klein / Journal of Statistical Planning and Inference 123 (2004) 117 – 131
131
Proof. As a prototype we establish (ii), the other proofs running quite similar. Due to the quadratic subspace property of Sym (s; H), the inverse of D2; j is a linear combination of the matrices I( m ) , W2 , W3 . Hence we may set I( m ) = (aI( m ) + bW2 + cW3 )D2; j , 2 2 2 for some a; b; c ∈ R. By means of the multiplication table in Lemma A.2, the right-hand side expands into a linear combination of I( m ) ; W2 ; W3 . Comparing coeScients results 2 in a linear equation system for a; b; c, with unique solution 1 a∗ = [(m − 3)(m − 4)j 2 − (m − 4)(m − 5)j + 4]; 6 1 b∗ = − (j − 2)[(m − 4)j + 2]; 6 2 c∗ = (j − 1)(j − 2); 6 where 6 = j(j − 1)(m − j)(m − j − 1)(m − 2)(m − 3). Hence, D2; j is nonsingular with inverse a∗ I( m ) + b∗ W2 + c∗ W3 . 2 C Now we adopt the block notation for Cj . The matrix C21; j is a scalar multiple of 22; j D1; j , and Lemma A.2 implies 2D2; j D1; j C11; j = c1 N31 for some constant c1 ∈ R: C21; j 2D2; j D This establishes R(Cj ) = R 2D1; j , and rankCj = ( m2 ). 2; j
References Cornell, J.A., 1990. Experiments with Mixtures. Designs, Models, and the Analysis of Mixture Data, Second Edition. Wiley, New York. Draper, N.R., Pukelsheim, F., 1998. Mixture models based on homogeneous polynomials. J. Statist. Plann. Inference 71, 303–311. Draper, N.R., Pukelsheim, F., 1999. Kiefer ordering of simplex designs for Irst- and second-degree mixture models. J. Statist. Plann. Inference 79, 325–348. Draper, N.R., Heiligers, B., Pukelsheim, F., 2000. Kiefer-ordering of simplex designs for second-degree mixture models with four or more ingredients. Ann. Statist. 28, 578–590. GaGke, N., 1981. Some classes of optimality criteria and optimal designs for complete two-way layouts. Ann. Statist. 9, 893–898. Galil, Z., Kiefer, J., 1977. Comparison of simplex designs for quadratic mixture models. Technometrics 19, 445–453. Kiefer, J., 1961. Optimum designs in regression problems, II. Ann. Math. Statist. 32, 298–325. Klein, T., 2001. Optimale VersuchsplVane im Kronecker-Modell zweiten Grades fVur Mischungsexperimente. WiWner, Augsburg. Klein, T., 2003. Invariant symmetric block matrices for the design of mixture experiments. Linear Algebra Appl., in press. Prescott, P., Dean, A.M., Draper, N.R., Lewis, S.M., 2002. Mixture experiments: Ill-conditioning and quadratic model speciIcation. Technometrics 44, 260–268. Pukelsheim, F., 1993. Optimal Design of Experiments. Wiley, New York. ScheGHe, H., 1958. Experiments with mixtures. J. Roy. Statist. Soc. Ser. B 20, 344–360. ScheGHe, H., 1963. The simplex-centroid design for experiments with mixtures. J. Roy. Statist. Soc. Ser. B 25, 235–251.
www.elsevier.com/locate/jspi
Optimal designs for second-degree Kronecker model mixture experiments Thomas Klein∗ Institut fur Mathematik, Universitat Augsburg, D-86135 Augsburg, Germany Received 11 August 2001; accepted 26 January 2003
Abstract The paper investigates optimal designs in the second-degree Kronecker model for mixture experiments. Three groups of novel results are presented: (i) characterization of feasible weighted centroid designs for a maximum parameter system, (ii) derivations of D-, A-, and E-optimal weighted centroid designs, and (iii) numerically p -optimal weighted centroid designs. Results on a quadratic subspace of invariant symmetric matrices containing the information matrices involved in the design problem serve as a main tool throughout the analysis. c 2003 Elsevier B.V. All rights reserved. MSC: 62K05 Keywords: Designs for mixtures; Design optimality; D-, A-, E-optimality; Kiefer’s p -criteria; Weighted centroid designs; Second-order Kronecker model; Quadratic subspace; Invariant symmetric matrices
1. Introduction A mixture experiment is an experiment in which the m factors t1 ; : : : ; tm are nonnegative and subject to the simplex restriction t1 + · · · + tm = 1, that is, the factors represent relative proportions of m ingredients blended in a mixture. Thus, the experimental conditions of a mixture m experiment are points in the probability simplex Tm = {t = (t1 ; : : : ; tm ) ∈ [0; 1]m : i=1 ti = 1}. Cornell (1990) lists numerous examples and applications of mixture experiments. A particular polynomial regression model for mixture experiments suggested by Draper and Pukelsheim (1998) is the second-degree Kronecker model 2
f : Tm → Rm ; ∗
t = (t1 ; : : : ; tm ) → t ⊗ t = (ti tj ) i; j=1; :::; m
ordered lexicographically
Tel.: +49-821-5982202; fax: +49-821-5982280. E-mail address: [email protected] (T. Klein).
c 2003 Elsevier B.V. All rights reserved. 0378-3758/$ - see front matter doi:10.1016/S0378-3758(03)00145-9
;
118
T. Klein / Journal of Statistical Planning and Inference 123 (2004) 117 – 131
with model equation E[Yt ] = f(t) =
m
ii ti2 +
i=1
m
( ij + ji )ti tj :
(1)
i; j=1 i¡j
Here Yt , the response under experimental condition t ∈ Tm , is taken to be a real-valued 2 random variable, while = ( 11 ; 12 ; : : : ; mm ) ∈ Rm is the vector of unknown parameters. All observations taken in an experiment are assumed to be uncorrelated and to have common unknown variance 2 ∈ (0; ∞). The paper determines optimal designs for a maximum subsystem of parameters in the Kronecker model, where Kiefer’s p -functions serve as optimality criteria. These results complement Galil’s and Kiefer’s (1977) paper on p -optimal designs in an alternative second-degree mixture model, the Sche78e model E[Yt ] = fS (t) =
n i=1
i ti +
m
ij ti tj
(2)
i; j=1 i¡j m+1
with parameter vector = (1 ; : : : ; m ; 12 ; 13 ; : : : ; m−1; m ) ∈ R( 2 ) , introduced by ScheGHe (1958, 1963). Both the Kronecker and the ScheGHe model are based on the same space of regression polynomials, but diGer in their choice of representing this space. Draper and Pukelsheim (1998) and Prescott et al. (2002) put forward several advantages of the Kronecker model, e.g. the homogeneity of regression terms and reduced ill-conditioning of information matrices. Interestingly, both models share the same invariance properties, see the appendix. An outline of the paper is as follows: Section 2 formulates the parameter system K of interest and the associated design problem. Based on the completeness Theorem 2.4 in Draper et al. (2000), the design problem is reduced to the class of weighted centroid designs. In Section 3 we present a characterization of weighted centroid designs under which K is estimable. In Section 4 we derive D-, A-, and E-optimal weighted centroid designs, while Section 5 presents numerical results on p -optimal weighted centroid designs. 2. The design problem 2
Since the Kronecker model’s full parameter vector ∈ Rm is not estimable, we consider a maximum parameter subsystem instead, that is, a parameter system K where the range R(K) coincides with the span of the regression range X := {f(t) : t ∈ Tm }. This formalizes the idea of estimating as many parameters as possible, see Klein (2001). A canonical maximum parameter subsystem is as follows. We denote the canonical unit vectors in Rm by e1 ; : : : ; em and set eij := ei ⊗ej for i; j= m 1; : : : ; m. The canonical unit vectors in R( 2 ) are denoted by Eij with 1 6 i ¡ j 6 m in m2 ×s lexicographic order. We write s := ( m+1 2 ), and deIne the matrix K := (K1 ; K2 ) ∈ R
T. Klein / Journal of Statistical Planning and Inference 123 (2004) 117 – 131
119
by K1 :=
m i=1
eii ei ;
K2 :=
m 1 (eij + eji )Eij : 2 i; j=1
(3)
i¡j
Then the matrix K satisIes R(K) = span X, and the maximum parameter subsystem considered in the following is: ( ii )16i6m 2 ∈ Rs for all ∈ Rm : K = 1 ( 2 ( ij + ji ))16i¡j6m An experimental design for a mixture experiment is a probability measure with Inite support on Tm . The statistical properties of a design are reJected by its moment matrix M () := Tm f(t)f(t) d. Here, the matrix M () contains the fourth moments of . Let Is denote the identity matrix in Rs×s . The information matrix 2
CK (M ()) := min{LM ()L |L ∈ Rs×m ;
LK = Is };
captures the amount of information the design contains on K , see Pukelsheim (1993, Section 3.10). As a consequence of K being a maximum parameter system, information matrices for K take a particularly simple form. The following lemma quotes this result from Klein (2001). Lemma 2.1. De9ne the matrix L˜ := (K K)−1 K = (K1 ; 2K2 ) ∈ Rs×m with K1 ; K2 from ˜ L˜ for all moment matrices M , that is, information (3). Then we have CK (M ) = LM matrices for K are linear transformations of moment matrices. 2
A scalar measure for the amount of information inherent to CK (M ()) is provided by Kiefer’s p -functions or matrix means p (C) = ((1=s)traceC p )1=p for all p ∈ [ − ∞; 1] and C from the set PD(s) of positive deInite s × s matrices. The family of matrix means includes the popular T-, D-, A-, and E-criteria, corresponding to p = 1, 0, −1, and −∞, respectively. The problem of Inding a design with maximum information on the parameter subsystem K can now be formulated as Maximize
p (CK (M ())) with ∈ T
subject to
CK (M ()) ∈ PD(s);
(4)
where T denotes the set of all designs on Tm . The side condition CK (M ()) ∈ PD(s) is equivalent to the existence of an unbiased linear estimator for K under , see Theorem 3.15 in Pukelsheim (1993). In this case, the design is called feasible for K . Any design solving problem (4) for a Ixed p ∈ [ − ∞; 1] is called p -optimal for K in T. Due to a result of Draper et al. (2000) on weighted centroid designs, the problem can be substantially reduced. Denition 2.2. For m ¿ 2 and j ∈ {1; : : : ; m}, the jth elementary centroid design j is the uniform distribution on the centroids of depth j, that is, on all points taking the
120
T. Klein / Journal of Statistical Planning and Inference 123 (2004) 117 – 131
j form 1=j i=1 meki ∈ Tm with 1 6 k1 ¡ k2 ¡ · · · ¡ kj−1 ¡ kj 6 m. A convex combination ( ) := j=1 j j with = ( 1 ; : : : ; m ) ∈ Tm is called a weighted centroid design with weight vector . We denote the set of all weighted centroid designs by (Tm ). Weighted centroid designs are exchangeable, that is, they are invariant under permutation of the mixture ingredients. This property leads to H-invariance of information matrices, HCK (M ())H = CK (M ())
for all H ∈ H;
∈ (Tm );
(5)
where H is the matrix group deIned in (A.1). The appendix summarizes some results on H-invariant information matrices which we utilize in the following sections. Theorems 6.4 and 7.4 in Draper and Pukelsheim (1999) and Theorem 2.4 in Draper et al. (2000) show that the set (Tm ) is an essentially complete class of designs relative to the Kiefer ordering of moment matrices in the Kronecker model. See Pukelsheim (1993, Chapter 14) for an introduction to the Kiefer ordering. Since the group H in (5) is a subgroup of the orthogonal group, all matrix means p are H-invariant and the target function of design problem (4) is isotonic relative to the Kiefer ordering. Therefore, the set (Tm ) is also essentially complete relative to the target function, that is, for any given design ∈ T there is a weighted centroid design () satisfying (p ◦ CK ◦ M )() ¿ (p ◦ CK ◦ M )() for all p ∈ [ − ∞; 1]. As a consequence, we may restrict the set of competing designs in problem (4) to (Tm ), thus obtaining a mere allocation problem for the weight vector ∈ Tm .
3. Feasible weighted centroid designs First we Ind an explicit description of feasible weighted centroid designs, that is, weighted centroid designs ( ) whose information matrix for K is regular. Let Cj := CK (M (j )) denote the information matrix of the elementary centroid design j , for 1 6 j 6 m. From the linearity of the information matrix mapping CK in Lemma 2.1 we get, for every ∈ Tm , CK (M (( ))) = j Cj with J( ) := {j = 1; : : : ; m | j ¿ 0}; (6) j∈J( )
whence R(CK (M (( )))) =
R(Cj ):
(7)
j∈J( )
Eq. (7) suggests studying the ranges of the information matrices C1 ; : : : ; Cm of the elementary centroid designs. This is done in Lemma A.4. By utilizing this result, the following theorem characterizes feasible weighted centroid designs for K . The essential ingredient is the indicator set J( ) from (6), collecting those subscripts j whose weight j is nonzero.
T. Klein / Journal of Statistical Planning and Inference 123 (2004) 117 – 131
121
Theorem 3.1. Let = ( 1 ; : : : ; m ) ∈ Tm be a weight vector of a weighted centroid design ( ) and set J( ) := {j | j ¿ 0} as in (6). Then the following assertions hold: (i) For m ∈ {2; 3}, the weighted centroid design ( ) is feasible for K if and only if J( ) ⊇ {1; 2}. (ii) For m ¿ 4, the weighted centroid design ( ) is feasible for K if and only if |J( ) ∩ {1; : : : ; m − 2}| ¿ 2
or
|J( ) ∩ {2; : : : ; m − 1}| ¿ 2:
Proof. We only prove (ii). Using Lemma A.4 and applying the dimension formula to the right-hand side of (7), we see that CK (M (( ))) is regular if and only if the set J( ) satisIes the condition given in (ii). For four or more ingredients, m ¿ 4, we combine our Theorem 3.1 with Theorem 3.2 from Draper et al. (2000). Corollary 3.2. Suppose m ¿ 4 and p ∈ [ − ∞; 1]. Then there is a weighted centroid design ( ) which is p -optimal for K in T and whose weight vector ∈ Tm satis9es, for some j ∈ {2; : : : ; m − 2}, one of the six conditions (i) J( ) = {1; j};
(iii) J( ) = {1; j; j + 1};
(v) J( ) = {j; j + 1; m};
(ii) J( ) = {j; j+1};
(iv) J( ) = {1; j; m};
(vi) J( ) = {1; j; j+1; m}:
4. Optimal weighted centroid designs We begin the investigation of p -optimal weighted centroid designs for K by adapting a general Equivalence Theorem to our speciIc problem. Theorem 4.1. Let ∈ Tm be the weight vector of a weighted centroid design ( ) which is feasible for K , and let J( ) be the set of active indices from (6). Furthermore, let C := CK (M (( ))) and p ∈ (−∞; 1]. Then ( ) is p -optimal for K in T if and only if = trace C p for all j ∈ J( ); p−1 trace Cj C 6 trace C p otherwise: Proof. We start from Theorem 7.19 in Pukelsheim (1993), that ( ) is p -optimal for K if and only if there is a generalized inverse G of M := M (( )) with trace M (())GKC p+1 K G 6 trace C p
for all ∈ Tm :
(8)
122
T. Klein / Journal of Statistical Planning and Inference 123 (2004) 117 – 131
With C = (K K)−1 K MK(K K)−1 , M = K(K K)−1 K M , and M (()) = KCK (M (( )))K from Lemma 2.1, the left-hand side is trace M (())GKC p+1 K G = trace (K GMK(K K)−1 ) CK (M (()))(K GMK(K K)−1 )C p−1 :
(9)
Due to the feasibility of ( ) we have R(K) = R(M ). Hence K = MZ for some Z, and so K GMK(K K)−1 = Z MGMK(K K)−1 = Z MK(K K)−1 = I( m ) . Now the right-hand 2 side of (9) simpliIes to trace CK (M (()))C p−1 , and (8) turns into trace CK (M (())) C p−1 6 trace C p for all ∈ Tm . According to equation (6) we can write the left-hand m p−1 , giving trace Cj C p−1 6 trace C p for all 1 6 j 6 m. Fiside as j=1 j trace Cj C nally, equality must hold for any j ∈ J( ). The proof is complete. The conditions for p -optimality provided above can be explicitly solved for the weight vector ∈ Tm in the special cases of D- and A-optimality (that is, p = 0 and p = −1) in Theorems 4.3 and 4.4, involving only weighted centroid designs with the Irst and second weight positive. The following theorem shows that this property implies uniqueness of the optimizers. Theorem 4.2. Let p ∈ (−∞; 1) and ( ) be a weighted centroid design that is p optimal for K in T , with ∈ Tm . (i) If J( ) = {1; 2}, then ( ) is the unique solution of problem (4). (ii) If J( ) = {1; 2; 3}, then there is no further exchangeable design Q ∈ T that is p -optimal for K in T . If there is a non-exchangeable p -optimal design for K , then all its support points are centroids of depths 1; 2, or 3. Proof. We only prove the Irst assertion. The proof of (ii) parallels that of (i) but excludes step (IV) below. In steps (II) and (III) of the proof we utilize some facts on fourth moments of exchangeable designs taken from Draper et al. (2000). That 4 is, an exchangeable design Q ∈ T possesses Ive distinct fourth moments & ( ) Q = t1 d , Q 4 Q &22 () Q = t12 t22 d , Q &211 () Q = t12 t2 t3 d , Q and &1111 () Q = t1 t2 t3 t4 d , Q &31 () Q = t13 t2 d , where the moments &211 () Q and &1111 () Q only occur for m ¿ 3 and m ¿ 4. The identity &31 () Q = &22 () Q is equivalent to Q being a weighted centroid design. The fourth moments of a weighted centroid design () are convex combinations of the fourth m moments of the elementary centroid designs, e.g. &4 (()) = j=1 j &4 (j ) for all =(1 ; : : : ; m ) ∈ Tm . Draper et al. (2000, p. 581) list the fourth moments of 1 ; : : : ; m . (I) First we show that M := M (( )) is the unique p -optimal moment matrix for K in M (T). Theorem 8.14 in Pukelsheim (1993) states that the moment matrix M () of any solution ∈ T to design problem (4) with p ∈ (−∞; 1) satisIes MGK = K;
(10)
where G is a generalized inverse of M (( )) satisfying optimality condition (8) for ( ). From our Theorem 4.1 it follows that any generalized inverse enjoys this property. We choose the Moore–Penrose inverse M + = K(K K)−1 C −1 (K K)−1 K with C :=
T. Klein / Journal of Statistical Planning and Inference 123 (2004) 117 – 131
123
CK (M (( ))). Inserting M + and using M () = KCK (M ())K , equation (10) becomes K = KCK (M ())C −1 . Post-multiplication with CK yields M = M (). (II) Suppose () with ∈ Tm is a further weighted centroid design which is p -optimal for K in (Tm ). According to part (I) of the proof, M = M (()) holds, that is, ( ) and () have the same fourth moments. We thus obtain m 1 1 2 = &4 (()) = j ; + 3 8m j m m j=1
m
j−1 2 = &31 (()) = j ; 8m(m − 1) j 3 m(m − 1) j=2
0 = &211 (()) =
m j=3
0 = &1111 (()) =
(j − 1)(j − 2) j ; j 3 m(m − 1)(m − 2)
m j=4
(j − 1)(j − 2)(j − 3) j : j 3 m(m − 1)(m − 2)(m − 3)
The equations for &211 (()) and &1111 (()) entail 3 = · · · = m = 0. Hence, the equations for &4 (()) and &31 (()) imply = . (III) Let Q ∈ T be an exchangeable design that is p -optimal for K in T. Then the Q Q Thus, Q is a weighted fourth moments of Q and ( ) coincide, implying &31 ()=& 22 (). centroid design, and so Q = ( ) as shown in (II). (IV) Now we drop exchangeability of the design ∈ T, but assume supp ⊆ supp ( ). From (I) we have CK (M ()) = C. Due to rank M (()) = |supp ( )| = s, the set {f(t) | t ∈ supp ( )} is linearly independent. Theorem 8.7 in Pukelsheim (1993) yields f ({f(t)}) = ( )f ({f(t)}) for all t ∈ supp ( ), and so = ( ). (V) Finally, let ∈ T be an arbitrary design which is p -optimal for K in T. Then, the exchangeable design Q := 1=m! R∈Perm(m) R obtained by averaging over the group Perm(m) is also optimal, and supp ⊆ supp . Q According to part (III), =() Q holds, which implies supp ⊆ supp ( ). Hence, part (IV) establishes = ( ). We may now derive optimal weighted centroid designs for the T-, D-, A- and E-criteria (that is, p = 1; 0; −1; −∞). Since the T-criterion is linear, the elementary centroid design 1 is the unique maximizer of 1 ◦ CK ◦ M in (Tm ), with optimum value (T = 1=s. However, 1 is not feasible for K due to Theorem 3.1. The cases of D- and A-optimality are covered by the following theorems. Since the D-criterion is invariant under reparameterizations, Theorem 4.3 applies to the Kronecker model as well as to the ScheGHe model. Theorem 4.3. In the second-degree mixture model with m ¿ 2 ingredients, the unique D-optimal design for K is 2 m−1 1 + 2 : ( (D) ) = m+1 m+1 The maximum of the D-criterion is (D = (1=s) 2−2(m−1)=(m+1) .
124
T. Klein / Journal of Statistical Planning and Inference 123 (2004) 117 – 131
Proof. Within the ScheGHe model (2), the design ( (D) ) is D-optimal for the full parameter vector, see Kiefer (1961). Since the D-criterion is invariant under linear re-parameterization of the space of regression polynomials (see e.g. GaGke, 1981), this design is also D-optimal for K within the Kronecker model. Uniqueness follows from Theorem 4.2. Theorem 4.4. In the second-degree Kronecker model with m ¿ 2 ingredients, the unique A-optimal design for K is √ √ 2(m − 1)[2(m − 1) − m + 3] 2(m − 1) m + 3 − m − 3 (A) ( ) = 1 + 2 : 4m2 − 9m + 1 4m2 − 9m + 1 The maximum of the A-criterion is (A =
√ 2(m + 1)[4m2 − 7m + 7 − 4(m − 1) m + 3] : m(4m2 − 9m + 1)2
Proof. Uniqueness of the proposed optimal design again follows from Theorem 4.2. The actual optimality proof consists of two parts: First, deriving an optimality candidate, and second, verifying its optimality. We concentrate on the Irst part. Let = ( 1 ; 2 ; 0; : : : ; 0) ∈ Tm be a weight vector with J( ) = {1; 2} and suppose ( ) is A-optimal for K in T. With the abbreviation C( ) := CK (M (( ))), Theorem 4.1 thus implies = trace C( )−1 for j ∈ {1; 2}; −2 (11) trace Cj C( ) 6 trace C( )−1 otherwise: The proof of Lemma A.4 already demonstrated how to compute an inverse matrix in the space Sym (s; H) by solving a system of linear equations. By the same approach we obtain the blocks of C( )−1 in the partitioning suggested by Lemma A.1, namely m m Im ; (C( )−1 )21 = − V1 ; (C( )−1 )11 = 1 2 1 (C( )−1 )22 =
m[4(m − 1) 1 + 2 ] m I( m ) + W2 : 2 2 1 2 4 1
Lemma A.2 yields the blocks of Cj C( )−2 for j = 1; : : : ; m. Using trace U2 = trace W2 = trace W3 = 0, we can now evaluate both sides of the optimality condition (11), trace C( )−1 =
m2 [4(m − 1)2 1 + (m + 3) 2 ] ; 4 1 2
trace Cj C( )−2 =
(12)
m2 32(m − 1)2 (j − 1) 12 4j 3 12 22
+ 16(m − 1)(j − 1)(j − 2) 1 2 + (j − 2)2 (m + j + 2) 22 :
(13)
T. Klein / Journal of Statistical Planning and Inference 123 (2004) 117 – 131
125
Equating (12) and (13) for j = 1 and j = 2 and inserting 2 = 1 − 1 , we obtain a quadratic equation for 1 . Its unique solution in (0; 1) is the weight 1(A) given in the assertion. By construction, the weight vector (A) = ( 1(A) ; 1 − 1(A) ; 0; : : : ; 0) satisIes the two equations in condition (11). The remaining m − 2 inequalities still need to be established. To this end, we compute trace Cj C( (A) )−2 − trace Cj+1 C( (A) )−2 m2 (j − 1) m(j 3 + 11j 2 + 12j + 4) = 3 4j (j + 1)3 12(A) √ + 8 m + 3(j 3 − 3j 2 − 5j − 2) − 2(j 3 − 17j 2 − 22j − 8)
for j = 2; : : : ; m − 1. Both terms on the right-hand side are found to be positive, whence trace C2 C( (A) )−2 ¿ · · · ¿ trace Cm C( (A) )−2 . Therefore, ( (A) ) is indeed A-optimal for K in T. The approach of the above proof, that is, solving the optimality condition from Theorem 4.1 for the weights 1 ; 2 , does not carry over to E-optimality. In fact, there are only partial results on E-optimal weighted centroid designs, quoted from Klein (2001) in the following two lemmas. Lemma 4.5. In the second-degree Kronecker model with m = 2 or m = 4 ingredients, the weighted centroid designs 5 6 for m = 2; 11 1 + 11 2 (E) ( ) = 13 12 9 for m = 4 46 1 + 23 2 + 46 3 are E-optimal for K in T, and there is no further E-optimal weighted centroid 1 1 design for K . The maxima of the E-criterion are (E = 11 for m = 2, and (E = 46 for m = 4. Lemma 4.6. In the second-degree Kronecker model for mixture experiments with m ¿ 5 ingredients, there is no weighted centroid design ( ) with J( ) ⊆ {1; 2; 3} which is E-optimal for K in T. As is seen in the following section, all numerically p -optimal weight vectors (p) , with p ∈ [ − 50; 1], satisfy J( (p) ) ⊆ {1; 2; 3}. However Lemma 4.6 says that this condition does not hold for E-optimal weight vectors when m ¿ 5. This is in line with Galil’s and Kiefer’s (1977) numerical results that p -optimal weight vectors for the full parameter vector of the ScheGHe model use the Irst three plus a further weight. 5. Numerical results The numerical computation of p -optimal weighted centroid designs for K divides into two steps, that is, Inding a local optimality candidate and examining a global
126
T. Klein / Journal of Statistical Planning and Inference 123 (2004) 117 – 131
Fig. 1. p -Optimal weights for K , with m = 4 and m = 10.
optimality condition for this candidate. As the Irst step requires greater eGort when the target function is not diGerentiable we omit the case of E-optimality. In order to Ind an optimality candidate we use Lemma A.3 to compute the eigenvalues of CK (M (( ))) as functions of ∈ Tm . As gradient methods for maximizing the target function of design problem (4) only work on open sets, each of the 6(m − 3) cases from Corollary 3.2 is treated separately. For the global optimality part we examine the condition given in Theorem 4.1, where the spectral decomposition of information matrices from Lemma A.3 is utilized once more. All p -optimal weighted centroid designs with p within the investigated range [ − 50; 1] are covered by the two cases J( ) = {1; 2} and J( ) = {1; 2; 3}. Since we are able to explicitly prove or disprove optimality of a candidate, the numerical algorithm can treat these two cases Irst and does not have to run through the whole list of cases mentioned above. Furthermore, Theorem 4.2 provides uniqueness of optimal designs ( ) with J( ) = {1; 2}; {1; 2; 3} within the set (Tm ). This justiIes the notation (p) for the p -optimal weight vector.
T. Klein / Journal of Statistical Planning and Inference 123 (2004) 117 – 131
(p)
(p)
127
(p)
Fig. 2. p -Optimal weight vectors in barycentric coordinates. The curve p → ( 1 ; 2 ; 3 ) in barycentric coordinates. For p ∈ (p2 ; 1] the curve runs along the 1 - 2 -edge of the simplex, for p ¡ p2 it extends into the interior. Open dots mark the optimal weight vectors for the T-, D-, A-, E-, and −50 -criteria.
Fig. 1 exempliIes our numerical results by listing the optimal weights 1(p) , 2(p) , 3(p) and the corresponding optimal values (p := (p ◦ CK ◦ M ◦ )( (p) ) for selected values of p ∈ [−∞; 1] and for m ∈ {4; 10}. The qualitative behavior of the functions p → i(p) with i = 1; 2; 3 and m ¿ 3 is captured by three critical values 0 ¿ p1 ¿ p2 ¿ p3 . For p ∈ (p2 ; 1) only two positive weights 1(p) and 2(p) are present, while for p ¡ p2 we also have 3(p) ¿ 0. For p descending from 1 to p2 , the weight 1(p) decreases when m 6 4. In the case m ¿ 5 a local minimum of 1(p) is found at p=p1 . Letting p descend from p2 , the weight 2(p) decreases while 3(p) increases. Interestingly, the Irst weight decreases but remains nearly constant. For m ¿ 7 a further change in monotonicity is observed at p = p3 where 1(p) starts a slight increase. Fig. 2 displays the curve p → ( 1(p) ; 2(p) ; 3(p) ) in barycentric coordinates for m = 3; : : : ; 10. Appendix A. H-Invariant information matrices The probability simplex Tm is invariant under permutation, that is, renumbering the m ingredients. Let Sm be the symmetric group of order m, and let R, = i=1 e,(i) ei be the matrix associated with , ∈ Sm . Then the group Perm(m) of m × m permutation matrices acts on Tm according to (R, ; (t1 ; : : : ; tm ) ) → (t,−1 (1) ; : : : ; t,−1 (m) ) and on T according to (R; ) → R = ◦ R−1 . If a design satisIes = R for all R ∈ Perm(m), we call it exchangeable. We deIne the group R, 0 : , ∈ Sm with H := H, = 0 S,
128
T. Klein / Journal of Statistical Planning and Inference 123 (2004) 117 – 131
S, :=
m i; j=1 i¡j
E(,(i); ,( j))↑ Eij ∈ Perm
m
for all , ∈ Sm
2
(A.1)
where (,(i); ,(j))↑ denotes the pair of indices ,(i), ,(j) in ascending order. The group H is a subgroup of Perm(s), and is isomorphic to Sm . As for the Kronecker model, the group H acts on an information matrix C for K through C → HCH . The equivariance property CK (M (R, )) = H, CK (M ())H, for , ∈ Sm , ∈ T links the action of H on CK (M (T)) to the action of Perm(m) on T. Within the ScheGHe model (2), we have MS (R, ) = H, MS ()H, for all moment matrices MS (). Therefore, we end up with the same invariance properties for both the Kronecker and the ScheGHe model. Since H is a subgroup of the orthogonal group, the space Sym (s; H) := {C ∈ Sym (s) : HCH = C for all H ∈ H} of H-invariant symmetric matrices is a quadratic subspace of the space Sym (s) of symmetric s × s matrices, that is, a subspace closed under formation of powers C n with n ∈ N. Basic facts about quadratic subspaces in experimental design can be found in Pukelsheim (1993, Chapters 13 and 14). The particular quadratic subspace Sym (s; H) is analyzed in Klein (2001, 2003). We brieJy summarize the results of this analysis. Lemma A.1. Let 1m denote the vector (1; : : : ; 1) ∈ Rm , and de9ne U2 := 1m 1m − Im ∈ Sym (m) and m m m m Eij (ei + ej ) ; V2 := Eij ek ∈ R( 2 )×m ; V1 := i; j=1 i¡j
W2 :=
W3 :=
i; j=1 k=1 i¡j k∈{i; j}
m
m
i; j=1 i¡j
k;‘=1 k¡‘ |{i; j}∩{k;‘}|=1
m
m
i; j=1 i¡j
k;‘=1 k¡‘ {i; j}∩{k;‘}=?
; Eij Ek‘
Eij Ek‘ ∈ Sym
m
2
:
Then any matrix C ∈ Sym (s; H) can be uniquely represented in the form aIm + bU2 cV1 + dV2 C = cV1 + dV2 eI m + fW2 + gW3 2
with coe=cients a; : : : ; g ∈ R. The terms containing V2 , W2 and W3 only occur for m ¿ 3 or m ¿ 4, respectively. In particular, dim Sym (s; H) = 4; 6; 7
for m = 2; m = 3 and m ¿ 4:
T. Klein / Journal of Statistical Planning and Inference 123 (2004) 117 – 131
129
When computing powers or inverses of matrices in Sym (s; H), products of the matrix blocks introduced in Lemma A.1 have to be evaluated. The following Lemma A.2 lists a multiplication table for Im , U2 , V1 , V2 , I( m ) , W2 , and W3 . 2
m
( 2 )×m , and W2 ; Lemma A.2.
any m ¿ 2, the matrices U2 ∈ Sym (m), V1 ; V2 ∈ R mFor W3 ∈ Sym ( 2 ) from Lemma A.1 satisfy the following equations:
(i) U22 V1 V1 (ii) V1 U2 W2 V1 W3 V1 (iii) V1 V1
= = = = = =
(m − 1)Im + (m − 2)U2 ; (m − 1)Im + U2 ; V1 + 2V2 ; (m − 2)V1 + 2V2 ; (m − 3)V2 ; 2I m + W2 ; 2
V1 V2 V2 V2 V 2 U2 W2 V 2 W 3 V2 V1 V2
= = = = = =
V2 V1 = (m − 2)U2 ;
m−1 Im + m−2 U2 : 2 2 (m − 2)V1 + (m − 3)V2 ; (m − 2)V1 + 2(m − 3)V2 ;
m−2 V1 + m−3 V2 : 2 2 V2 V1 = W2 + 2W3 ;
V2 V2 = (m − 2)I m + (m − 3)W2 + (m − 4)W3 ; W22 W32
2
= 2(m − 2)I m + (m − 2)W2 + 4W3 ; 2
m + m−3 W2 + m−4 W3 ; I = m−2 2 2 2 2
W2 W3 = W3 W2 = (m − 3)W2 + 2(m − 4)W3 : A spectral analysis of H-invariant symmetric matrices can be derived which comprises descriptions of eigenvalues as well as eigenspaces. This expands earlier Indings of Galil and Kiefer (1977), who already provided the eigenvalue part of the analysis. We summarize the results on eigenvalues in the following lemma. Our results on eigenspaces are not used explicitly in this paper, whence the reader is referred to Klein (2001, 2003) for the description of eigenspaces. Lemma A.3. Let a; : : : ; g ∈ R be the coe=cients of a matrix C ∈ Sym (s; H) w.r.t. the basis given in Lemma A.1 with d; f and g occurring only when m ¿ 3 or m ¿ 4, respectively. Furthermore, de9ne
2 D1 := a + (m − 1)b − e − 2(m − 2)f − m−2 g 2 +2(m − 1)[2c + (m − 2)d]2 ; D2 := [a − b − e − (m − 4)f + (m − 3)g]2 + 4(m − 2)(c − d)2 : Then, in the case m ¿ 4, the matrix C has the eigenvalues 51 := e − 2f + g;
√ 1 a + (m − 1)b + e + 2(m − 2)f + m−2 g ± D1 ; 2 2 √ 1 := [a − b + e + (m − 4)f − (m − 3)g ± D2 ] 2
52; 3 := 54; 5
130
T. Klein / Journal of Statistical Planning and Inference 123 (2004) 117 – 131
with multiplicities m(m − 3)=2; 1, and m − 1, respectively. In the case m = 2 only the eigenvalues 52 ; 53 ; 54 occur, whereas for m = 3 there are four eigenvalues 52 ; : : : ; 55 . Information matrices of weighted centroid designs are H-invariant, see equation (5). According to (6), we focus on the information matrices C1 ; : : : ; Cm of elementary centroid designs. For j = 1; : : : ; m we block partition C11; j C21; j Cj = C21; j C22; j and apply Lemma 2.1 to the representation of the moment matrices M (j ) given by Draper et al. (2000). We obtain 1 1 j−1 C11; j = 3 Im + 3 U2 ; j m m−1 j m C21; j =
2 j−1 j−2 2 j−1 V1 + 3 V2 ; j3 m m − 1 j m m−1 m−2
C22; j =
4 j−1 j−2 4 j−1 j−2 j−3 4 j−1 W2 + 3 W3 ; I m + 3 j m m−1 m−2 j m m−1 m−2 m−3 j3 m m − 1 2
where the terms containing V2 ; W2 and W3 only occur for m ¿ 3 and m ¿ 4, respectively. Motivated by equation (7), the following lemma describes the ranges of C1 ; : : : ; Cm , and intersections of these ranges. Lemma A.4. Let 1m denote the vector (1; : : : ; 1) ∈ Tm . Then the following results hold for the information matrices C1 ; : : : ; Cm of elementary centroid designs: (i) R(C1 ) = Rm × {0( m ) }, and rank C1 = m, for all m ¿ 2. 2 D (ii) R(Cj ) = R 2D1; j , and rank Cj = ( m2 ), for m ¿ 4 and j ∈ {2; : : : ; m − 2}, where 2; j
D1; j := (m − 3)[(m − 2)V1 + (j − 2)V2 ]; D2; j := (m − 2)(m − 3)I m + (j − 2)(m − 3)W2 + (j − 2)(j − 3)W3 : (iii) R(Cm−1 ) = R
U2 2V2
(iv) R(Cm ) = span
2
, and rank Cm−1 = m, for all m ¿ 3.
1m 2 1 m 2
, and rank Cm = 1, for all m ¿ 2.
(v) For all m ¿ 2 and i; j ∈ {1; : : : ; m}, i = j, we have {0m } × R(V1 )⊥ for i; j ∈ {2; : : : ; m − 2}; R(Ci ) ∩ R(Cj ) = {0s } otherwise: In particular, dim(R(Ci ) ∩ R(Cj )) =
m (m − 3) for all i; j ∈ {2; : : : ; m − 2}, i = j. 2
T. Klein / Journal of Statistical Planning and Inference 123 (2004) 117 – 131
131
Proof. As a prototype we establish (ii), the other proofs running quite similar. Due to the quadratic subspace property of Sym (s; H), the inverse of D2; j is a linear combination of the matrices I( m ) , W2 , W3 . Hence we may set I( m ) = (aI( m ) + bW2 + cW3 )D2; j , 2 2 2 for some a; b; c ∈ R. By means of the multiplication table in Lemma A.2, the right-hand side expands into a linear combination of I( m ) ; W2 ; W3 . Comparing coeScients results 2 in a linear equation system for a; b; c, with unique solution 1 a∗ = [(m − 3)(m − 4)j 2 − (m − 4)(m − 5)j + 4]; 6 1 b∗ = − (j − 2)[(m − 4)j + 2]; 6 2 c∗ = (j − 1)(j − 2); 6 where 6 = j(j − 1)(m − j)(m − j − 1)(m − 2)(m − 3). Hence, D2; j is nonsingular with inverse a∗ I( m ) + b∗ W2 + c∗ W3 . 2 C Now we adopt the block notation for Cj . The matrix C21; j is a scalar multiple of 22; j D1; j , and Lemma A.2 implies 2D2; j D1; j C11; j = c1 N31 for some constant c1 ∈ R: C21; j 2D2; j D This establishes R(Cj ) = R 2D1; j , and rankCj = ( m2 ). 2; j
References Cornell, J.A., 1990. Experiments with Mixtures. Designs, Models, and the Analysis of Mixture Data, Second Edition. Wiley, New York. Draper, N.R., Pukelsheim, F., 1998. Mixture models based on homogeneous polynomials. J. Statist. Plann. Inference 71, 303–311. Draper, N.R., Pukelsheim, F., 1999. Kiefer ordering of simplex designs for Irst- and second-degree mixture models. J. Statist. Plann. Inference 79, 325–348. Draper, N.R., Heiligers, B., Pukelsheim, F., 2000. Kiefer-ordering of simplex designs for second-degree mixture models with four or more ingredients. Ann. Statist. 28, 578–590. GaGke, N., 1981. Some classes of optimality criteria and optimal designs for complete two-way layouts. Ann. Statist. 9, 893–898. Galil, Z., Kiefer, J., 1977. Comparison of simplex designs for quadratic mixture models. Technometrics 19, 445–453. Kiefer, J., 1961. Optimum designs in regression problems, II. Ann. Math. Statist. 32, 298–325. Klein, T., 2001. Optimale VersuchsplVane im Kronecker-Modell zweiten Grades fVur Mischungsexperimente. WiWner, Augsburg. Klein, T., 2003. Invariant symmetric block matrices for the design of mixture experiments. Linear Algebra Appl., in press. Prescott, P., Dean, A.M., Draper, N.R., Lewis, S.M., 2002. Mixture experiments: Ill-conditioning and quadratic model speciIcation. Technometrics 44, 260–268. Pukelsheim, F., 1993. Optimal Design of Experiments. Wiley, New York. ScheGHe, H., 1958. Experiments with mixtures. J. Roy. Statist. Soc. Ser. B 20, 344–360. ScheGHe, H., 1963. The simplex-centroid design for experiments with mixtures. J. Roy. Statist. Soc. Ser. B 25, 235–251.