- Email: [email protected]
New classes of D-optimal edge designs
Journal of Statistical Planning and Inference 136 (2006) 248 – 259 www.elsevier.com/locate/jspi
New classes of D-optimal edge designs C. Koukouvinos∗ , S. Stylianou Department of Mathematics, National Technical University of Athens, Zografou 15773, Athens, Greece Received 1 July 2003; accepted 2 June 2004 Available online 25 July 2004
Abstract Edge designs are screening experimental designs that allow a model-independent estimate of the set of relevant variables, thus providing more robustness than traditional designs. In this paper we construct new classes of D-optimal edge designs. This construction uses weighing matrices of order n and weight k together with permutation matrices of order n to obtain D-optimal edge designs. One linear and one quadratic simulated screening scenarios are studied and compared using linear regression and edge designs analysis. © 2004 Elsevier B.V. All rights reserved. MSC: primary 62K15; secondary 05B20 Keywords: Weighing matrix; Edge design; Construction; Screening experiment
1. Introduction Screening experiments are used for the correct identification of relevant variables out of a large number of variables that possibly influence some characteristic y. The aim of screening is to detect cheaply this subset of influential variables; then more accurate designs on lower-dimensional subspace are used to account for nonlinearities. Since the non-influential variables are not explored in subsequent measurements, an important requirement for screening designs is that they place measurements such that the region of interest is well explored. The number of trials N needed for the identification of the active factors should be kept as small as possible. The experimental expense ∗ Corresponding author.
E-mail address: [email protected] (C. Koukouvinos). 0378-3758/$ - see front matter © 2004 Elsevier B.V. All rights reserved. doi:10.1016/j.jspi.2004.06.003
C. Koukouvinos, S. Stylianou / Journal of Statistical Planning and Inference 136 (2006) 248 – 259
249
for empirical model building depends drastically on the number of variables in the experiment. For example, in the case of quadratic modelling, the number of experiments needed is O(n2 ). An acceptable reduction of runs to O(n), such that correct identification of relevance is retained, is of great importance since the complexity is significantly reduced from quadratic to linear order. This will have practical purpose in cases where n is rather large. When n is small, the experiment is usually performed with two-level fractional factorial designs suitable for optimal estimation of the coefficients {} in a model such as y = 0 +
n
j x j +
j =1
n n
ij xi xj .
(1)
i=1 j =i+1
To estimate all main effects i and all interactions ij , again a large number of runs (O(n2 )) is needed. To be accurate, note that N = (n2 + n + 2)/2 runs are needed to estimate all main and all interactions in an experiment involving n variables. For example, when n = 35, the number of runs needed is 631. Usually, screening for a large number of variables is frequently done by the use of linear models and orthogonal designs such as fractional factorial designs and Plackett–Burman designs Plackett and Burman, 1946. However, the final section explores by an example the risk of screening with linear models and models including second order interactions and quadratic terms. The aim of edge designs is that they can be evaluated in two independent ways, see Elster and Neumaier (1995). In the following, In denotes the identity matrix of order n, j T = (1, 1, . . . , 1) is a 1 × n vector with all entries equal to 1, and J = jj T . Further, x denotes the Euclidean norm, that is x2 = x12 + · · · + xn2 , and (x, y) = x1 y1 + · · · + xn yn is the standard inner product in Rn . Let Q = [−1, 1]n , and denote X the N × (n + 1) design matrix with Xij = fj (x i ). When x i ∈ Q is the ith setting of the variable vector x i = (x1i , . . . , xni )T and f is given by the assumed linear relation y i = 0 +
n j =1
j xji + εi ,
(2)
where y i denotes the ith measurement and εi the measurement error; the εi (i = 1, . . . , N ) are assumed to be independent and identically distributed. Details on screening designs and how to analyze unreplicated factorials can be found in Aboukalam and Al-Shiha, 2001; Box and Meyer, 1986. Edge designs can be used without the assumption of a specific model, and in any model the set of relevant factors will be found. It is known that the design matrix X of any design for linear model (2) satisfies det XT X N n+1 , with equality if and only if X T X = N I .
(3)
The Deff of a design is defined to be Deff = (D/D ∗ )1/(n+1) , where D is the determinant of XT X, and D ∗ = N n+1 is the upper bound given in (3). A design X is called D-optimal if it achieves this bound and thus it has Deff = 1.
250
C. Koukouvinos, S. Stylianou / Journal of Statistical Planning and Inference 136 (2006) 248 – 259
Elster and Neumaier, 1995 introduced a new class of screening designs which allow a model-independent test for active variables. This is achieved by arranging the measurements into a set E of pairs that within these pairs the coordinates differ in one component only. They called such pairs, edges since in the optimal case they are located at the edges of the cube and refer to designs consisting of a collection of edges as edge designs. They have shown that the simplest edge design, “one factor at a time” design, has very low Deff . An edge design is said to be minimal if it has just n edges, one for each variable. Elster and Neumaier (1995) conjecture that a minimal edge design cannot be D-optimal and they prove this result in the most natural class of minimal edge designs, namely those with N = 2n runs and n disjoint edges. Independent of any particular model, data collected with edge designs may be evaluated using the assumption that only a few, say p, of the n factors are active; i.e. contribute to the variability in the observations. This so-called factor sparsity assumption, mentioned e.g. by Lenth (1989), is very natural in screening experiments and implies that almost all differences zij :=yi − yj , (i, j ) ∈ E, consist of noise only. If we assume that the noise in the data is additive, normally distributed with zero mean and variance 2 , the n − p of the zi are normally distributed with zero mean and variance 22 . Because of the unknown number of outliers, the variance must be estimated in a robust way. For example, we can use the median estimate =
median{|zij | : (i, j ) ∈ E} , 21/2 × 0.675
(4)
which is consistent when p =0 (see Lenth, 1989), and hence can be expected to give reliable results when p >n. Outliers then determine active factors. In general, the use of edge designs guarantees that irrelevant variables are never treated as relevant, in contrast to classical designs, and there is only a very small chance that a relevant variable is not correctly recognized, which occurs when the two function values nearly agree on the corresponding edges. Thus, screening with edge designs is robust. A permutation matrix is the identity with its rows re-ordered. Permutation matrices are orthogonal and so if P is a permutation matrix then P T P = P P T = I. The product of permutation matrices is also a permutation matrix. A (0, 1, −1) matrix W = W (n, k) of order n satisfying W W T = W T W = kI n is called a weighing matrix of order n and weight k. A weighing matrix W = W (n, k) for which W T = −W is called a skew-weighing matrix. For more results on the existence and constructions of weighing matrices the reader should consider Craigen (1996) or Geramita and Seberry (1979). Elster and Neumaier (1995) constructed edge designs using conference matrices of order n. A conference matrix is actually a weighing matrix W =W (n, n−1) of order n and weight n−1. Conference matrices exist only for some odd values n while weighing matrices W (n, k) can be constructed for several values of n and k. Thus, the construction of edge designs from weighing matrices is a natural generalization of the construction from conference matrices. The next section presents the construction of new D-optimal edge designs using weighing matrices.
C. Koukouvinos, S. Stylianou / Journal of Statistical Planning and Inference 136 (2006) 248 – 259
251
2. Construction of new D-optimal edge designs 2.1. D-optimal edge designs from W (n, n − 2) The next theorem describes how a weighing matrix of order n and weight n − 2 can be used to construct D-optimal edge designs with 4n runs and 2n edges. Theorem 1. Let W = W (n, n − 2) be a weighing matrix of order n and weight n − 2. If there exist P, Q permutation matrices of order n such that W + P + Q is {1, −1} matrix and W T P = −P T W , then there exists a D-optimal edge design with 4n runs and 2n edges. Proof. Set
W +P +Q −W − P − Q , W +P −Q −W − P + Q
j j X= j j
X ∈ {1, −1}4n×(n+1) .
It is easy to see that the ith row together with the (2n+i)th row and the (n+i)th row together with the (3n + i)th row, i = 1, 2, . . . , n, define edges. Thus, there are two edges for each variable. Consequently, X is an edge design with 4n-runs and 2n-edges (two edges for each variable). To show that X is a D-optimal edge design we need to prove that XX T = 4nI n+1 . With simple calculations we have
4n 0 . X X= 0 4(W T W + P T P + QT Q + P T W + W T P )
T
Since, W T W =(n−2)In , P T P =QT Q=In , P T W =−W T P we have that X T X =4nI n+1 and the proof is completed. The following example illustrates how the above theorem works. Example 1. Let
1 0 −1 W =W (6, 4)= −1 0 −1 0 0 1 0 0 0 0 0 0 Q= 0 1 0 1 0 0 0 0 0
0 0 −1 −1 1 1 1 1 0 0 0 −1 −1 1 −1 −1 −1 0 0 0 0 0 0 1 0 1 0 . 0 0 0 0 0 0 1 0 0
1 1 0 1 −1 0
1 0 0 1 1 0 , P= −1 0 0 0 1 0
1 0 0 0 0 0
0 0 0 1 0 0
0 0 1 0 0 0
0 0 0 0 0 1
0 0 0 , 0 1 0
252
C. Koukouvinos, S. Stylianou / Journal of Statistical Planning and Inference 136 (2006) 248 – 259
The conditions stated in Theorem 1 are satisfied since W P T = −P W T , W + P + Q is {1, −1} matrix and P P T = P T P = QT Q = QQT = I . Following the proof of Theorem 1 we obtain the design matrix X, XT
1 1 1 = 1 ¯ 1 1 1
1 1 1¯ 1 1 1 1
1 1¯ 1 1 1 1 1
1 1¯ 1 1 1¯ 1 1¯
1 1 1¯ 1 1¯ 1¯ 1
1 1¯ 1¯ 1¯ 1 1 1
1 1¯ 1¯ 1¯ 1 1¯ 1¯
1 1¯ 1 1¯ 1¯ 1¯ 1¯
1 1 1¯ 1¯ 1¯ 1¯ 1¯
1 1 1¯ 1¯ 1 1¯ 1
1 1¯ 1 1¯ 1 1 1¯
1 1 1 1 1¯ 1¯ 1¯
1 1 1 1¯ 1¯ 1 1
1 1 1¯ 1 1 1 1¯
1 1¯ 1 1 1 1¯ 1
1 1¯ 1¯ 1 1¯ 1 1¯
1 1¯ 1¯ 1 1¯ 1¯ 1
1 1¯ 1¯ 1¯ 1¯ 1 1
1 1¯ 1¯ 1 1 1¯ 1¯
1 1¯ 1 1¯ 1¯ 1¯ 1
1 1 1¯ 1¯ 1¯ 1 1¯
1 1 1 1¯ 1 1¯ 1
1 1 1 1¯ 1 1 1¯
1 1 1 1 , 1 1¯ 1¯
where 1¯ stands for −1. This matrix satisfies equation (3) with equality since X T X = 24I7 , and thus this is a D-optimal (Det XT X = 247 = N n+1 , Deff = 1) edge design with 24 runs and 12 edges (two edges for each variable). The condition stated in Theorem 1 is trivially satisfied if the weighing matrix used is skew. This result is very useful for the choice of a proper weighing matrix and it is proved in Corollary 1. Corollary 1. Let W = W (n, n − 2) be a skew-weighing matrix of order n and weight n − 2. Then there exists a D-optimal edge design with 4n-runs and 2n-edges. Proof. Set P = In and Q = (qij ), where qij = 1 − (abs(wij ) + pij ). Using these matrices we have that W T P = W T In = W T = −W = −P T W and from Theorem 1 we have the result. Example 2. Let
0 1 −1 0 W = W (4, 2) = 0 1 −1 0 0 0 1 0 0 0 0 1 Q= . 1 0 0 0 0 1 0 0
0 1 −1 0 , 0 −1 1 0
1 0 P = 0 0
0 1 0 0
0 0 1 0
0 0 , 0 1
Since W is a skew weighing matrix the condition required by Corollary 1 is satisfied and thus W T P = −P T W , W + P + Q is {1, −1} matrix and P P T = P T P = QT Q = QQT = I . From Corollary 1 we obtain the design matrix X, 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1¯ 1 1¯ 1¯ 1 1¯ 1 1 1¯ 1¯ 1¯ 1¯ 1 1 1 T X = 1 1 1 1 1¯ 1¯ 1¯ 1¯ 1 1 1 1¯ 1¯ 1¯ 1¯ 1 . 1 1¯ 1 1 1¯ 1 1¯ 1¯ 1¯ 1¯ 1 1 1 1 1¯ 1¯ 1 1 1¯ 1 1¯ 1¯ 1 1¯ 1 1¯ 1¯ 1 1¯ 1 1 1¯
C. Koukouvinos, S. Stylianou / Journal of Statistical Planning and Inference 136 (2006) 248 – 259
253
Table 1 The existence of W (n, n − 2) for small values of n n 4 6
8
10
12
16
18
20
In Craigen (1996) see Table
52.30
52.29
52.29
52.30
52.28
52.29
52.29
52.28
This matrix satisfies Eq. (3) with equality since X T X = 16I5 and thus, this is a D-optimal (Det XT X = 165 = N n+1 , Deff = 1) edge design with 16 runs and 8 edges (two edges for each variable). Weighing matrices W (n, n − 2) exist in many cases, especially when n is even. An explicit table of the parameters of smaller instances for which the construction works is given in Table 1. If n is odd and a W (n, n − 2) exists, then n − 2 is a perfect square. In this case, the only possible value for n 20 is n = 11 which does not exist, see Craigen, 1996, Table 52.28. More weighing matrices with these parameters can be found in Craigen (1996). 2.2. D-optimal edge designs from W (n, k) The Theorems presented in this section use any weighing matrix and suitable permutation matrices to obtain D-optimal edge designs. Theorem 2. Let W = W (n, k) be a weighing matrix of order n and weight k. If there exist P1 , P2 , . . . , Pn−k permutation matrices such that W + P1 + P2 + . . . + Pn−k is a {1, −1} matrix, then there exists a D-optimal edge design with n2n−k edges, and the matrix n−k+1 )×(n+1) X ∈ {1, −1}(n2 is the design matrix. Proof. Set Z the 2n−k+1 × (n − k + 1) matrix which corresponds to a full 2n−k+1 factorial design and W P1 P A= .2 . .. Pn−k We define Y = (Z ⊗ I )A, where ⊗ is the Kronecker product and j X = ... Y . j
We shall show that X is the desirable (n2n−k+1 ) × (n + 1) matrix. This matrix consists of n2n−k edges. With simple calculations we have n−k+1
n2 0 T X X= . 0 Y TY
254
C. Koukouvinos, S. Stylianou / Journal of Statistical Planning and Inference 136 (2006) 248 – 259
Since, Y T Y = n2n−k+1 In we have that X T X = n2n−k+1 In+1 and the proof is completed. The next example shows how Theorem 2 can be used to give a D-optimal edge design with 192 runs and 96 edges. Example 3. Let
1 1 1 0 ¯ 1 1 W = W (12, 9) = ¯ 1 0 1 1 ¯ 1 0
0 0 0 1 0 0 P1 = 0 0 0 0 0 0 0 0 0 0 0 0 P2 = 0 0 0 0 0 1
0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0
0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0
1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0
1¯ 0 1 1 1 1 0 1¯ 1 1¯ 0 1
1 0 1¯ 1 1¯ 0 1 1 1 1¯ 0 1
0 1¯ 1 1¯ 0 1 1 1 1 0 1 1¯
1¯ 1 0 0 1 1¯ 1 1 1 1 1¯ 0
0 0 0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 1 0 0 0 1 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0
0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0
0 0 0 0 0 0 , 0 0 1 0 0 0 0 0 0 0 1 0 , 0 0 0 0 0 0
1 1 1 1¯ 1 0 0 1 1¯ 1¯ 0 1
1 1 1 1 0 1¯ 1 1¯ 0 0 1 1¯
0 1 1¯ 1 1 1 1¯ 1 0 0 1 1¯
1 1¯ 0 1 1 1 1 0 1¯ 1 1¯ 0
1¯ 1 0 0 1¯ 1 1 0 1¯ 1 1 1
1 0 1¯ 1¯ 1 0 0 1¯ 1 1 1 1
0 1¯ 1 1 0 1¯ , 1¯ 1 0 1 1 1
C. Koukouvinos, S. Stylianou / Journal of Statistical Planning and Inference 136 (2006) 248 – 259
0 0 0 0 0 0 P3 = 0 1 0 0 0 0
0 0 0 0 0 0 0 0 0 0 1 0
0 0 0 0 0 0 0 0 0 1 0 0
0 0 0 0 0 0 0 0 1 0 0 0
0 0 0 0 0 0 0 0 0 0 0 1
0 0 0 0 0 0 1 0 0 0 0 0
0 0 0 0 0 1 0 0 0 0 0 0
0 0 0 0 1 0 0 0 0 0 0 0
0 0 0 1 0 0 0 0 0 0 0 0
0 0 1 0 0 0 0 0 0 0 0 0
0 1 0 0 0 0 0 0 0 0 0 0
255
1 0 0 0 0 0 . 0 0 0 0 0 0
The conditions stated in Theorem 2 are satisfied since W + P1 + P2 + P3 is {1, −1} matrix and Pi PiT = PiT Pi = I , for i = 1, 2, 3. Following the procedure described in the proof of Theorem 2, the design matrix X is obtained. This matrix satisfies equation (3) with equality since XT X = 192I13 and thus, this is a D-optimal (Det X T X = 19213 = N n+1 , Deff = 1) edge design with 192 runs and 96 edges (eight edges for each variable). When the number of trials is large and experiments are expensive it is important to keep the number of runs as low as possible. With the following Theorem the number of trials is diminished to half, compared to Theorem 2. To achieve that, one more restriction on weighing matrix and on the first permutation matrix, is required. This restriction and the obtained result are described in Theorem 3. Theorem 3. Let W = W (n, k) be a weighing matrix of order n and weight k. If there exist P1 , P2 , . . . , Pn−k permutation matrices such that W + P1 + P2 + . . . + Pn−k is a {1, −1} matrix and W T P1 = −P1T W then there exists a D-optimal edge design with n2n−k−1 edges, n−k and the matrix X ∈ {1, −1}(n2 )×(n+1) is the design matrix. Proof. Set Z the 2n−k+1 × (n − k + 1) matrix which corresponds to a full 2n−k+1 factorial design. We select from Z the rows(runs) zi that satisfy zi1 zi2 = 1. These are 2n−k rows(runs) and suppose that these rows form the (2n−k ) × (n − k + 1) matrix H. Let
W P1 P2 .. .
A=
Pn−k and define Y = (H ⊗ I )A and
j X = ... j
Y
,
X ∈ {1, −1}(n2
n−k )×(n+1)
.
256
C. Koukouvinos, S. Stylianou / Journal of Statistical Planning and Inference 136 (2006) 248 – 259
This matrix consists of n2n−k−1 edges. The design matrix X satisfies X T X = n2n−k In+1 since Y T Y = n2n−k In . The importance of Theorem 3 is shown using Example 4. In this example we obtain a D-optimal edge design with 96 runs and 48 edges. Example 4. Let
1 1 0 ¯ 1 1 1 W = W (12, 9) = ¯ 1 0 1 ¯ 1 ¯ 1 0
0 0 1 0 0 0 P1 = 0 0 0 0 0 0 0 0 0 0 0 0 P2 = 0 0 0 0 0 1
0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0
1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0
1¯ 1¯ 1 0 1 1 1 1 0 1¯ 0 1
1¯ 1 0 1 0 1 1¯ 0 1¯ 1 1¯ 1
1 1¯ 0 1 0 1 1 1¯ 1 1 1¯ 0
1 0 1 1 1 1¯ 0 1 1¯ 0 1¯ 1¯
1 0 1¯ 1 1 0 0 1¯ 1¯ 1¯ 1 1
1¯ 1 1 1 1 0 0 1¯ 1 0 1 1¯
0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 1 0 0 0 0 0 0 0 0 0
0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0
0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 1 , 0 0 0 0 0 0 0 1 0 0 0 0 , 0 0 0 0 0 0
0 1¯ 1¯ 0 1 1 1¯ 1 0 1 1 1¯
0 1 1¯ 1 1¯ 1 1 1 0 1¯ 0 1¯
1¯ 1¯ 1¯ 1 0 1¯ 1¯ 0 1 1¯ 1¯ 0
0 1¯ 1 0 1¯ 1 1¯ 1¯ 1¯ 1¯ 0 1¯
1¯ 0 1¯ 1¯ 1 0 , 1 1¯ 1¯ 0 1¯ 1¯
C. Koukouvinos, S. Stylianou / Journal of Statistical Planning and Inference 136 (2006) 248 – 259
0 0 0 0 0 0 P3 = 0 1 0 0 0 0
0 0 0 0 0 0 0 0 1 0 0 0
0 0 1 0 0 0 0 0 0 0 0 0
0 0 0 1 0 0 0 0 0 0 0 0
1 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 1 0 0 0 0 0 0 0
0 1 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 1 0 0 0 0 0 0
0 0 0 0 0 0 1 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 1
0 0 0 0 0 0 0 0 0 0 1 0
257
0 0 0 0 0 0 . 0 0 0 1 0 0
It is easy to see that the conditions stated in Theorem 3 are satisfied since W +P1 +P2 +P3 is {1, −1} matrix, W T P1 =−P1T W and Pi PiT =PiT Pi =I for i =1, 2, 3. Following the proof of Theorem 3, the design matrix X is obtained. This matrix satisfies equation (3) with equality since XT X = 96I13 and thus this is a D-optimal (Det X T X = 9613 = N n+1 , Deff = 1) edge design with 96 runs and 48 edges (four edges for each variable). Remark 1. Although the number of variables in Examples 3 and 4 is the same, Example 4 gives an important reduction of runs. Moreover, the restriction required by Theorem 3 is not hard to be fulfilled (see for example, Corollary 2). Corollary 2. If W = W (n, k) is a skew-weighing matrix of order n and weight k then there exists a D-optimal edge design with n2n−k−1 edges and n2n−k runs. The design matrix is n−k X ∈ {1, −1}(n2 )×(n+1) . Proof. Set P1 = In and select any P2 , P3 , . . . , Pn−k permutation matrices such that W + P1 + P2 + · · · + Pn−k is a {1, −1} matrix. Then W T P1 = W T In = −W = −P1T W and thus the result is obtained using Theorem 3. Some weighing matrices W (n, k) are given in Craigen (1996). In the next section, we present an example to show that edge designs are very useful when nonlinearities appear. 3. Simulated screening scenarios The data used in these two simulated screening scenarios are obtained using first a noisy linear function h and then a noisy quadratic function f, where y = h(x1 , x2 , x3 , x4 , x5 , x6 ) = 10 + 6x3 + 16x4 + ε,
(5)
u = f (x1 , x2 , x3 , x4 , x5 , x6 ) = −3 + x3 + x4 + 3(x3 + 2x4 )2 + ε,
(6)
and in both cases xi ∈ [−1, 1] (i = 1, . . . , 6), ε ∼ N(0, 1).
258
C. Koukouvinos, S. Stylianou / Journal of Statistical Planning and Inference 136 (2006) 248 – 259
Table 2 Edge design analysis in the noisy linear function simulated case x3
x6
x5
x2
x1
x4
y1 − y13 11.81
y2 − y14 −3.01
y3 − y15 1.33
y4 − y16 0.02
y5 − y17 0.95
y6 − y18 31.28
Table 3 Edge design analysis in the noisy quadratic function simulated case x3
x6
x5
x2
x1
x4
u1 − u13 −20.47
u2 − u14 3.53
u3 − u15 0.34
u4 − u16 1.48
u5 − u17 0.78
u6 − u18 22.91
In both cases we used the design X as it is given in Example 1. The response vector y in the case of the noisy linear function is y T = (0.22, 30.86, 34.06, 0.31, 0.54, 19.76, 19.98, −11.84, −11.87, 20.89, 17.86, −1.15, −11.60, 33.87, 32.73, 0.29, −0.41, −11.52, 31.20, −11.96, −12.09, 21.46, 20.79, 31.62) while the response vector u in the case of the noisy quadratic function is uT = (0.39, 24.53, 25.97, 0.94, 0.54, −0.93, −1.12, 20.56, 22.87, 1.35, −1.11, 0.22, 20.86, 28.06, 26.31, 0.54, −0.24, 21.98, 26.16, 22.13, 22.89, −2.14, −1.15, 26.40). Case 1: Noisy linear function h and the corresponding response vector y:An analysis of the data with the software package SPSS, using linear regression revealed, as expected, two main effects (x3 and x4 ) and gave an estimated linear model h(x)=6.01x3 +16.09x4 +10.16+ε, with ε of mean zero and standard deviation = 0.89. An analysis of the edges (see Table 2) shows the same relevant variables (x3 and x4 ) with an absolute value of more than four times the robust estimation ˜ = 2.2 of as this is defined in (4). Note that it is enough to use only one edge for each variable since the results obtained from edges on the same variable are similar. Since relevant factors found by both methods are the same, we conclude that these are the correct active factors. In this case, screening with six columns of a Hadamard matrix of order 24 reveals the same relevant factors. Case 2: Noisy quadratic function f and the corresponding response vector u: An analysis of the data with the software package SPSS, using linear regression revealed, as expected, no main effects and gave an estimated linear model f (x) = 11.84 + ε, with ε of mean zero and standard deviation =12.5. An analysis of the edges (see Table 3) revealed two relevant variables (x3 and x4 ) with an absolute value of more than four times the robust estimation ˜ = 2.6 (using the six first edges) of as this is defined in (4). Note that it is enough to use only one edge for each variable since the results obtained from edges on the same variable are similar. Since relevant factors found by the two methods are different, we conclude that the union of the two subsets should contain the correct active factors but the expected model could not be linear. In this case, screening with six columns of a Hadamard matrix of order 24 reveals no relevant factors, and thus no further investigation is possible.
C. Koukouvinos, S. Stylianou / Journal of Statistical Planning and Inference 136 (2006) 248 – 259
259
Although the second case is extreme, it is easy to run into trouble of the same kind with a more realistic, nonlinear functional relationship. For these reasons, a robust screening design, such as the edge design, can be proved to be very useful. References Aboukalam, M.A.F., Al-Shiha, A.A., 2001. A robust analysis for unreplicated factorial experiments. Comput. Statist. Data Anal. 36, 31–46. Box, G.E.P., Meyer, R.D., 1986. An analysis for unreplicated fractional factorials. Technometrics 28, 11–18. Craigen, R., 1996. Weighing matrices and conference matrices. in: Colbourn, C.J., Dinitz, J.H. (Eds.), The CRC Handbook of Combinatorial Designs. CRC Press, Boca Raton, FL, pp. 496–504. Elster, C., Neumaier, A., 1995. Screening by conference designs. Biometrika 82, 589–602. Geramita, A.V., Seberry, J., 1979. Orthogonal Designs: Quadratic Forms and Hadamard Matrices. Marcel Dekker, New York–Basel. Lenth, R.S., 1989. Quick and easy analysis of unreplicated factorials. Technometrics 31, 469–473. Plackett, R.L., Burman, J.P., 1946. The design of optimum multifactorial experiments. Biometrika 33, 305–325.
New classes of D-optimal edge designs C. Koukouvinos∗ , S. Stylianou Department of Mathematics, National Technical University of Athens, Zografou 15773, Athens, Greece Received 1 July 2003; accepted 2 June 2004 Available online 25 July 2004
Abstract Edge designs are screening experimental designs that allow a model-independent estimate of the set of relevant variables, thus providing more robustness than traditional designs. In this paper we construct new classes of D-optimal edge designs. This construction uses weighing matrices of order n and weight k together with permutation matrices of order n to obtain D-optimal edge designs. One linear and one quadratic simulated screening scenarios are studied and compared using linear regression and edge designs analysis. © 2004 Elsevier B.V. All rights reserved. MSC: primary 62K15; secondary 05B20 Keywords: Weighing matrix; Edge design; Construction; Screening experiment
1. Introduction Screening experiments are used for the correct identification of relevant variables out of a large number of variables that possibly influence some characteristic y. The aim of screening is to detect cheaply this subset of influential variables; then more accurate designs on lower-dimensional subspace are used to account for nonlinearities. Since the non-influential variables are not explored in subsequent measurements, an important requirement for screening designs is that they place measurements such that the region of interest is well explored. The number of trials N needed for the identification of the active factors should be kept as small as possible. The experimental expense ∗ Corresponding author.
E-mail address: [email protected] (C. Koukouvinos). 0378-3758/$ - see front matter © 2004 Elsevier B.V. All rights reserved. doi:10.1016/j.jspi.2004.06.003
C. Koukouvinos, S. Stylianou / Journal of Statistical Planning and Inference 136 (2006) 248 – 259
249
for empirical model building depends drastically on the number of variables in the experiment. For example, in the case of quadratic modelling, the number of experiments needed is O(n2 ). An acceptable reduction of runs to O(n), such that correct identification of relevance is retained, is of great importance since the complexity is significantly reduced from quadratic to linear order. This will have practical purpose in cases where n is rather large. When n is small, the experiment is usually performed with two-level fractional factorial designs suitable for optimal estimation of the coefficients {} in a model such as y = 0 +
n
j x j +
j =1
n n
ij xi xj .
(1)
i=1 j =i+1
To estimate all main effects i and all interactions ij , again a large number of runs (O(n2 )) is needed. To be accurate, note that N = (n2 + n + 2)/2 runs are needed to estimate all main and all interactions in an experiment involving n variables. For example, when n = 35, the number of runs needed is 631. Usually, screening for a large number of variables is frequently done by the use of linear models and orthogonal designs such as fractional factorial designs and Plackett–Burman designs Plackett and Burman, 1946. However, the final section explores by an example the risk of screening with linear models and models including second order interactions and quadratic terms. The aim of edge designs is that they can be evaluated in two independent ways, see Elster and Neumaier (1995). In the following, In denotes the identity matrix of order n, j T = (1, 1, . . . , 1) is a 1 × n vector with all entries equal to 1, and J = jj T . Further, x denotes the Euclidean norm, that is x2 = x12 + · · · + xn2 , and (x, y) = x1 y1 + · · · + xn yn is the standard inner product in Rn . Let Q = [−1, 1]n , and denote X the N × (n + 1) design matrix with Xij = fj (x i ). When x i ∈ Q is the ith setting of the variable vector x i = (x1i , . . . , xni )T and f is given by the assumed linear relation y i = 0 +
n j =1
j xji + εi ,
(2)
where y i denotes the ith measurement and εi the measurement error; the εi (i = 1, . . . , N ) are assumed to be independent and identically distributed. Details on screening designs and how to analyze unreplicated factorials can be found in Aboukalam and Al-Shiha, 2001; Box and Meyer, 1986. Edge designs can be used without the assumption of a specific model, and in any model the set of relevant factors will be found. It is known that the design matrix X of any design for linear model (2) satisfies det XT X N n+1 , with equality if and only if X T X = N I .
(3)
The Deff of a design is defined to be Deff = (D/D ∗ )1/(n+1) , where D is the determinant of XT X, and D ∗ = N n+1 is the upper bound given in (3). A design X is called D-optimal if it achieves this bound and thus it has Deff = 1.
250
C. Koukouvinos, S. Stylianou / Journal of Statistical Planning and Inference 136 (2006) 248 – 259
Elster and Neumaier, 1995 introduced a new class of screening designs which allow a model-independent test for active variables. This is achieved by arranging the measurements into a set E of pairs that within these pairs the coordinates differ in one component only. They called such pairs, edges since in the optimal case they are located at the edges of the cube and refer to designs consisting of a collection of edges as edge designs. They have shown that the simplest edge design, “one factor at a time” design, has very low Deff . An edge design is said to be minimal if it has just n edges, one for each variable. Elster and Neumaier (1995) conjecture that a minimal edge design cannot be D-optimal and they prove this result in the most natural class of minimal edge designs, namely those with N = 2n runs and n disjoint edges. Independent of any particular model, data collected with edge designs may be evaluated using the assumption that only a few, say p, of the n factors are active; i.e. contribute to the variability in the observations. This so-called factor sparsity assumption, mentioned e.g. by Lenth (1989), is very natural in screening experiments and implies that almost all differences zij :=yi − yj , (i, j ) ∈ E, consist of noise only. If we assume that the noise in the data is additive, normally distributed with zero mean and variance 2 , the n − p of the zi are normally distributed with zero mean and variance 22 . Because of the unknown number of outliers, the variance must be estimated in a robust way. For example, we can use the median estimate =
median{|zij | : (i, j ) ∈ E} , 21/2 × 0.675
(4)
which is consistent when p =0 (see Lenth, 1989), and hence can be expected to give reliable results when p >n. Outliers then determine active factors. In general, the use of edge designs guarantees that irrelevant variables are never treated as relevant, in contrast to classical designs, and there is only a very small chance that a relevant variable is not correctly recognized, which occurs when the two function values nearly agree on the corresponding edges. Thus, screening with edge designs is robust. A permutation matrix is the identity with its rows re-ordered. Permutation matrices are orthogonal and so if P is a permutation matrix then P T P = P P T = I. The product of permutation matrices is also a permutation matrix. A (0, 1, −1) matrix W = W (n, k) of order n satisfying W W T = W T W = kI n is called a weighing matrix of order n and weight k. A weighing matrix W = W (n, k) for which W T = −W is called a skew-weighing matrix. For more results on the existence and constructions of weighing matrices the reader should consider Craigen (1996) or Geramita and Seberry (1979). Elster and Neumaier (1995) constructed edge designs using conference matrices of order n. A conference matrix is actually a weighing matrix W =W (n, n−1) of order n and weight n−1. Conference matrices exist only for some odd values n while weighing matrices W (n, k) can be constructed for several values of n and k. Thus, the construction of edge designs from weighing matrices is a natural generalization of the construction from conference matrices. The next section presents the construction of new D-optimal edge designs using weighing matrices.
C. Koukouvinos, S. Stylianou / Journal of Statistical Planning and Inference 136 (2006) 248 – 259
251
2. Construction of new D-optimal edge designs 2.1. D-optimal edge designs from W (n, n − 2) The next theorem describes how a weighing matrix of order n and weight n − 2 can be used to construct D-optimal edge designs with 4n runs and 2n edges. Theorem 1. Let W = W (n, n − 2) be a weighing matrix of order n and weight n − 2. If there exist P, Q permutation matrices of order n such that W + P + Q is {1, −1} matrix and W T P = −P T W , then there exists a D-optimal edge design with 4n runs and 2n edges. Proof. Set
W +P +Q −W − P − Q , W +P −Q −W − P + Q
j j X= j j
X ∈ {1, −1}4n×(n+1) .
It is easy to see that the ith row together with the (2n+i)th row and the (n+i)th row together with the (3n + i)th row, i = 1, 2, . . . , n, define edges. Thus, there are two edges for each variable. Consequently, X is an edge design with 4n-runs and 2n-edges (two edges for each variable). To show that X is a D-optimal edge design we need to prove that XX T = 4nI n+1 . With simple calculations we have
4n 0 . X X= 0 4(W T W + P T P + QT Q + P T W + W T P )
T
Since, W T W =(n−2)In , P T P =QT Q=In , P T W =−W T P we have that X T X =4nI n+1 and the proof is completed. The following example illustrates how the above theorem works. Example 1. Let
1 0 −1 W =W (6, 4)= −1 0 −1 0 0 1 0 0 0 0 0 0 Q= 0 1 0 1 0 0 0 0 0
0 0 −1 −1 1 1 1 1 0 0 0 −1 −1 1 −1 −1 −1 0 0 0 0 0 0 1 0 1 0 . 0 0 0 0 0 0 1 0 0
1 1 0 1 −1 0
1 0 0 1 1 0 , P= −1 0 0 0 1 0
1 0 0 0 0 0
0 0 0 1 0 0
0 0 1 0 0 0
0 0 0 0 0 1
0 0 0 , 0 1 0
252
C. Koukouvinos, S. Stylianou / Journal of Statistical Planning and Inference 136 (2006) 248 – 259
The conditions stated in Theorem 1 are satisfied since W P T = −P W T , W + P + Q is {1, −1} matrix and P P T = P T P = QT Q = QQT = I . Following the proof of Theorem 1 we obtain the design matrix X, XT
1 1 1 = 1 ¯ 1 1 1
1 1 1¯ 1 1 1 1
1 1¯ 1 1 1 1 1
1 1¯ 1 1 1¯ 1 1¯
1 1 1¯ 1 1¯ 1¯ 1
1 1¯ 1¯ 1¯ 1 1 1
1 1¯ 1¯ 1¯ 1 1¯ 1¯
1 1¯ 1 1¯ 1¯ 1¯ 1¯
1 1 1¯ 1¯ 1¯ 1¯ 1¯
1 1 1¯ 1¯ 1 1¯ 1
1 1¯ 1 1¯ 1 1 1¯
1 1 1 1 1¯ 1¯ 1¯
1 1 1 1¯ 1¯ 1 1
1 1 1¯ 1 1 1 1¯
1 1¯ 1 1 1 1¯ 1
1 1¯ 1¯ 1 1¯ 1 1¯
1 1¯ 1¯ 1 1¯ 1¯ 1
1 1¯ 1¯ 1¯ 1¯ 1 1
1 1¯ 1¯ 1 1 1¯ 1¯
1 1¯ 1 1¯ 1¯ 1¯ 1
1 1 1¯ 1¯ 1¯ 1 1¯
1 1 1 1¯ 1 1¯ 1
1 1 1 1¯ 1 1 1¯
1 1 1 1 , 1 1¯ 1¯
where 1¯ stands for −1. This matrix satisfies equation (3) with equality since X T X = 24I7 , and thus this is a D-optimal (Det XT X = 247 = N n+1 , Deff = 1) edge design with 24 runs and 12 edges (two edges for each variable). The condition stated in Theorem 1 is trivially satisfied if the weighing matrix used is skew. This result is very useful for the choice of a proper weighing matrix and it is proved in Corollary 1. Corollary 1. Let W = W (n, n − 2) be a skew-weighing matrix of order n and weight n − 2. Then there exists a D-optimal edge design with 4n-runs and 2n-edges. Proof. Set P = In and Q = (qij ), where qij = 1 − (abs(wij ) + pij ). Using these matrices we have that W T P = W T In = W T = −W = −P T W and from Theorem 1 we have the result. Example 2. Let
0 1 −1 0 W = W (4, 2) = 0 1 −1 0 0 0 1 0 0 0 0 1 Q= . 1 0 0 0 0 1 0 0
0 1 −1 0 , 0 −1 1 0
1 0 P = 0 0
0 1 0 0
0 0 1 0
0 0 , 0 1
Since W is a skew weighing matrix the condition required by Corollary 1 is satisfied and thus W T P = −P T W , W + P + Q is {1, −1} matrix and P P T = P T P = QT Q = QQT = I . From Corollary 1 we obtain the design matrix X, 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1¯ 1 1¯ 1¯ 1 1¯ 1 1 1¯ 1¯ 1¯ 1¯ 1 1 1 T X = 1 1 1 1 1¯ 1¯ 1¯ 1¯ 1 1 1 1¯ 1¯ 1¯ 1¯ 1 . 1 1¯ 1 1 1¯ 1 1¯ 1¯ 1¯ 1¯ 1 1 1 1 1¯ 1¯ 1 1 1¯ 1 1¯ 1¯ 1 1¯ 1 1¯ 1¯ 1 1¯ 1 1 1¯
C. Koukouvinos, S. Stylianou / Journal of Statistical Planning and Inference 136 (2006) 248 – 259
253
Table 1 The existence of W (n, n − 2) for small values of n n 4 6
8
10
12
16
18
20
In Craigen (1996) see Table
52.30
52.29
52.29
52.30
52.28
52.29
52.29
52.28
This matrix satisfies Eq. (3) with equality since X T X = 16I5 and thus, this is a D-optimal (Det XT X = 165 = N n+1 , Deff = 1) edge design with 16 runs and 8 edges (two edges for each variable). Weighing matrices W (n, n − 2) exist in many cases, especially when n is even. An explicit table of the parameters of smaller instances for which the construction works is given in Table 1. If n is odd and a W (n, n − 2) exists, then n − 2 is a perfect square. In this case, the only possible value for n 20 is n = 11 which does not exist, see Craigen, 1996, Table 52.28. More weighing matrices with these parameters can be found in Craigen (1996). 2.2. D-optimal edge designs from W (n, k) The Theorems presented in this section use any weighing matrix and suitable permutation matrices to obtain D-optimal edge designs. Theorem 2. Let W = W (n, k) be a weighing matrix of order n and weight k. If there exist P1 , P2 , . . . , Pn−k permutation matrices such that W + P1 + P2 + . . . + Pn−k is a {1, −1} matrix, then there exists a D-optimal edge design with n2n−k edges, and the matrix n−k+1 )×(n+1) X ∈ {1, −1}(n2 is the design matrix. Proof. Set Z the 2n−k+1 × (n − k + 1) matrix which corresponds to a full 2n−k+1 factorial design and W P1 P A= .2 . .. Pn−k We define Y = (Z ⊗ I )A, where ⊗ is the Kronecker product and j X = ... Y . j
We shall show that X is the desirable (n2n−k+1 ) × (n + 1) matrix. This matrix consists of n2n−k edges. With simple calculations we have n−k+1
n2 0 T X X= . 0 Y TY
254
C. Koukouvinos, S. Stylianou / Journal of Statistical Planning and Inference 136 (2006) 248 – 259
Since, Y T Y = n2n−k+1 In we have that X T X = n2n−k+1 In+1 and the proof is completed. The next example shows how Theorem 2 can be used to give a D-optimal edge design with 192 runs and 96 edges. Example 3. Let
1 1 1 0 ¯ 1 1 W = W (12, 9) = ¯ 1 0 1 1 ¯ 1 0
0 0 0 1 0 0 P1 = 0 0 0 0 0 0 0 0 0 0 0 0 P2 = 0 0 0 0 0 1
0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0
0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0
1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0
1¯ 0 1 1 1 1 0 1¯ 1 1¯ 0 1
1 0 1¯ 1 1¯ 0 1 1 1 1¯ 0 1
0 1¯ 1 1¯ 0 1 1 1 1 0 1 1¯
1¯ 1 0 0 1 1¯ 1 1 1 1 1¯ 0
0 0 0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 1 0 0 0 1 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0
0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0
0 0 0 0 0 0 , 0 0 1 0 0 0 0 0 0 0 1 0 , 0 0 0 0 0 0
1 1 1 1¯ 1 0 0 1 1¯ 1¯ 0 1
1 1 1 1 0 1¯ 1 1¯ 0 0 1 1¯
0 1 1¯ 1 1 1 1¯ 1 0 0 1 1¯
1 1¯ 0 1 1 1 1 0 1¯ 1 1¯ 0
1¯ 1 0 0 1¯ 1 1 0 1¯ 1 1 1
1 0 1¯ 1¯ 1 0 0 1¯ 1 1 1 1
0 1¯ 1 1 0 1¯ , 1¯ 1 0 1 1 1
C. Koukouvinos, S. Stylianou / Journal of Statistical Planning and Inference 136 (2006) 248 – 259
0 0 0 0 0 0 P3 = 0 1 0 0 0 0
0 0 0 0 0 0 0 0 0 0 1 0
0 0 0 0 0 0 0 0 0 1 0 0
0 0 0 0 0 0 0 0 1 0 0 0
0 0 0 0 0 0 0 0 0 0 0 1
0 0 0 0 0 0 1 0 0 0 0 0
0 0 0 0 0 1 0 0 0 0 0 0
0 0 0 0 1 0 0 0 0 0 0 0
0 0 0 1 0 0 0 0 0 0 0 0
0 0 1 0 0 0 0 0 0 0 0 0
0 1 0 0 0 0 0 0 0 0 0 0
255
1 0 0 0 0 0 . 0 0 0 0 0 0
The conditions stated in Theorem 2 are satisfied since W + P1 + P2 + P3 is {1, −1} matrix and Pi PiT = PiT Pi = I , for i = 1, 2, 3. Following the procedure described in the proof of Theorem 2, the design matrix X is obtained. This matrix satisfies equation (3) with equality since XT X = 192I13 and thus, this is a D-optimal (Det X T X = 19213 = N n+1 , Deff = 1) edge design with 192 runs and 96 edges (eight edges for each variable). When the number of trials is large and experiments are expensive it is important to keep the number of runs as low as possible. With the following Theorem the number of trials is diminished to half, compared to Theorem 2. To achieve that, one more restriction on weighing matrix and on the first permutation matrix, is required. This restriction and the obtained result are described in Theorem 3. Theorem 3. Let W = W (n, k) be a weighing matrix of order n and weight k. If there exist P1 , P2 , . . . , Pn−k permutation matrices such that W + P1 + P2 + . . . + Pn−k is a {1, −1} matrix and W T P1 = −P1T W then there exists a D-optimal edge design with n2n−k−1 edges, n−k and the matrix X ∈ {1, −1}(n2 )×(n+1) is the design matrix. Proof. Set Z the 2n−k+1 × (n − k + 1) matrix which corresponds to a full 2n−k+1 factorial design. We select from Z the rows(runs) zi that satisfy zi1 zi2 = 1. These are 2n−k rows(runs) and suppose that these rows form the (2n−k ) × (n − k + 1) matrix H. Let
W P1 P2 .. .
A=
Pn−k and define Y = (H ⊗ I )A and
j X = ... j
Y
,
X ∈ {1, −1}(n2
n−k )×(n+1)
.
256
C. Koukouvinos, S. Stylianou / Journal of Statistical Planning and Inference 136 (2006) 248 – 259
This matrix consists of n2n−k−1 edges. The design matrix X satisfies X T X = n2n−k In+1 since Y T Y = n2n−k In . The importance of Theorem 3 is shown using Example 4. In this example we obtain a D-optimal edge design with 96 runs and 48 edges. Example 4. Let
1 1 0 ¯ 1 1 1 W = W (12, 9) = ¯ 1 0 1 ¯ 1 ¯ 1 0
0 0 1 0 0 0 P1 = 0 0 0 0 0 0 0 0 0 0 0 0 P2 = 0 0 0 0 0 1
0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0
1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0
1¯ 1¯ 1 0 1 1 1 1 0 1¯ 0 1
1¯ 1 0 1 0 1 1¯ 0 1¯ 1 1¯ 1
1 1¯ 0 1 0 1 1 1¯ 1 1 1¯ 0
1 0 1 1 1 1¯ 0 1 1¯ 0 1¯ 1¯
1 0 1¯ 1 1 0 0 1¯ 1¯ 1¯ 1 1
1¯ 1 1 1 1 0 0 1¯ 1 0 1 1¯
0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 1 0 0 0 0 0 0 0 0 0
0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0
0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 1 , 0 0 0 0 0 0 0 1 0 0 0 0 , 0 0 0 0 0 0
0 1¯ 1¯ 0 1 1 1¯ 1 0 1 1 1¯
0 1 1¯ 1 1¯ 1 1 1 0 1¯ 0 1¯
1¯ 1¯ 1¯ 1 0 1¯ 1¯ 0 1 1¯ 1¯ 0
0 1¯ 1 0 1¯ 1 1¯ 1¯ 1¯ 1¯ 0 1¯
1¯ 0 1¯ 1¯ 1 0 , 1 1¯ 1¯ 0 1¯ 1¯
C. Koukouvinos, S. Stylianou / Journal of Statistical Planning and Inference 136 (2006) 248 – 259
0 0 0 0 0 0 P3 = 0 1 0 0 0 0
0 0 0 0 0 0 0 0 1 0 0 0
0 0 1 0 0 0 0 0 0 0 0 0
0 0 0 1 0 0 0 0 0 0 0 0
1 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 1 0 0 0 0 0 0 0
0 1 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 1 0 0 0 0 0 0
0 0 0 0 0 0 1 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 1
0 0 0 0 0 0 0 0 0 0 1 0
257
0 0 0 0 0 0 . 0 0 0 1 0 0
It is easy to see that the conditions stated in Theorem 3 are satisfied since W +P1 +P2 +P3 is {1, −1} matrix, W T P1 =−P1T W and Pi PiT =PiT Pi =I for i =1, 2, 3. Following the proof of Theorem 3, the design matrix X is obtained. This matrix satisfies equation (3) with equality since XT X = 96I13 and thus this is a D-optimal (Det X T X = 9613 = N n+1 , Deff = 1) edge design with 96 runs and 48 edges (four edges for each variable). Remark 1. Although the number of variables in Examples 3 and 4 is the same, Example 4 gives an important reduction of runs. Moreover, the restriction required by Theorem 3 is not hard to be fulfilled (see for example, Corollary 2). Corollary 2. If W = W (n, k) is a skew-weighing matrix of order n and weight k then there exists a D-optimal edge design with n2n−k−1 edges and n2n−k runs. The design matrix is n−k X ∈ {1, −1}(n2 )×(n+1) . Proof. Set P1 = In and select any P2 , P3 , . . . , Pn−k permutation matrices such that W + P1 + P2 + · · · + Pn−k is a {1, −1} matrix. Then W T P1 = W T In = −W = −P1T W and thus the result is obtained using Theorem 3. Some weighing matrices W (n, k) are given in Craigen (1996). In the next section, we present an example to show that edge designs are very useful when nonlinearities appear. 3. Simulated screening scenarios The data used in these two simulated screening scenarios are obtained using first a noisy linear function h and then a noisy quadratic function f, where y = h(x1 , x2 , x3 , x4 , x5 , x6 ) = 10 + 6x3 + 16x4 + ε,
(5)
u = f (x1 , x2 , x3 , x4 , x5 , x6 ) = −3 + x3 + x4 + 3(x3 + 2x4 )2 + ε,
(6)
and in both cases xi ∈ [−1, 1] (i = 1, . . . , 6), ε ∼ N(0, 1).
258
C. Koukouvinos, S. Stylianou / Journal of Statistical Planning and Inference 136 (2006) 248 – 259
Table 2 Edge design analysis in the noisy linear function simulated case x3
x6
x5
x2
x1
x4
y1 − y13 11.81
y2 − y14 −3.01
y3 − y15 1.33
y4 − y16 0.02
y5 − y17 0.95
y6 − y18 31.28
Table 3 Edge design analysis in the noisy quadratic function simulated case x3
x6
x5
x2
x1
x4
u1 − u13 −20.47
u2 − u14 3.53
u3 − u15 0.34
u4 − u16 1.48
u5 − u17 0.78
u6 − u18 22.91
In both cases we used the design X as it is given in Example 1. The response vector y in the case of the noisy linear function is y T = (0.22, 30.86, 34.06, 0.31, 0.54, 19.76, 19.98, −11.84, −11.87, 20.89, 17.86, −1.15, −11.60, 33.87, 32.73, 0.29, −0.41, −11.52, 31.20, −11.96, −12.09, 21.46, 20.79, 31.62) while the response vector u in the case of the noisy quadratic function is uT = (0.39, 24.53, 25.97, 0.94, 0.54, −0.93, −1.12, 20.56, 22.87, 1.35, −1.11, 0.22, 20.86, 28.06, 26.31, 0.54, −0.24, 21.98, 26.16, 22.13, 22.89, −2.14, −1.15, 26.40). Case 1: Noisy linear function h and the corresponding response vector y:An analysis of the data with the software package SPSS, using linear regression revealed, as expected, two main effects (x3 and x4 ) and gave an estimated linear model h(x)=6.01x3 +16.09x4 +10.16+ε, with ε of mean zero and standard deviation = 0.89. An analysis of the edges (see Table 2) shows the same relevant variables (x3 and x4 ) with an absolute value of more than four times the robust estimation ˜ = 2.2 of as this is defined in (4). Note that it is enough to use only one edge for each variable since the results obtained from edges on the same variable are similar. Since relevant factors found by both methods are the same, we conclude that these are the correct active factors. In this case, screening with six columns of a Hadamard matrix of order 24 reveals the same relevant factors. Case 2: Noisy quadratic function f and the corresponding response vector u: An analysis of the data with the software package SPSS, using linear regression revealed, as expected, no main effects and gave an estimated linear model f (x) = 11.84 + ε, with ε of mean zero and standard deviation =12.5. An analysis of the edges (see Table 3) revealed two relevant variables (x3 and x4 ) with an absolute value of more than four times the robust estimation ˜ = 2.6 (using the six first edges) of as this is defined in (4). Note that it is enough to use only one edge for each variable since the results obtained from edges on the same variable are similar. Since relevant factors found by the two methods are different, we conclude that the union of the two subsets should contain the correct active factors but the expected model could not be linear. In this case, screening with six columns of a Hadamard matrix of order 24 reveals no relevant factors, and thus no further investigation is possible.
C. Koukouvinos, S. Stylianou / Journal of Statistical Planning and Inference 136 (2006) 248 – 259
259
Although the second case is extreme, it is easy to run into trouble of the same kind with a more realistic, nonlinear functional relationship. For these reasons, a robust screening design, such as the edge design, can be proved to be very useful. References Aboukalam, M.A.F., Al-Shiha, A.A., 2001. A robust analysis for unreplicated factorial experiments. Comput. Statist. Data Anal. 36, 31–46. Box, G.E.P., Meyer, R.D., 1986. An analysis for unreplicated fractional factorials. Technometrics 28, 11–18. Craigen, R., 1996. Weighing matrices and conference matrices. in: Colbourn, C.J., Dinitz, J.H. (Eds.), The CRC Handbook of Combinatorial Designs. CRC Press, Boca Raton, FL, pp. 496–504. Elster, C., Neumaier, A., 1995. Screening by conference designs. Biometrika 82, 589–602. Geramita, A.V., Seberry, J., 1979. Orthogonal Designs: Quadratic Forms and Hadamard Matrices. Marcel Dekker, New York–Basel. Lenth, R.S., 1989. Quick and easy analysis of unreplicated factorials. Technometrics 31, 469–473. Plackett, R.L., Burman, J.P., 1946. The design of optimum multifactorial experiments. Biometrika 33, 305–325.