Algorithmic search for baseline minimum aberration designs

Journal of Statistical Planning and Inference 149 (2014) 172–182

Contents lists available at ScienceDirect

Journal of Statistical Planning and Inference journal homepage: www.elsevier.com/locate/jspi

Algorithmic search for baseline minimum aberration designs Pei Li a,1, Arden Miller b,n, Boxin Tang a a Department of Statistics and Actuarial Science Simon Fraser University, 8888 University Drive, Burnaby, British Columbia, Canada V5A 1S6 b Department of Statistics, University of Auckland, Private Bag 92019, Auckland 1142, New Zealand

a r t i c l e in f o

abstract

Article history: Received 6 November 2013 Received in revised form 17 February 2014 Accepted 17 February 2014 Available online 26 February 2014

This paper reviews the formulation of the K-aberration criterion for baseline two-level designs and the efficient complete search algorithm developed by Mukerjee and Tang (2012). An efficient incomplete search algorithm is proposed that can be used to find near optimal baseline designs in situations where the complete search algorithm is not feasible. Lower bounds for values of K2 and K3 are established. A catalogue of optimal (or near optimal) 20-run baseline designs is provided. & 2014 Elsevier B.V. All rights reserved.

Keywords: Baseline model Minimum aberration Orthogonal array Search algorithm Two-level design

1. Introduction Two level factorial and fractional factorial designs have been extensively studied since the pioneering work of Fisher and Yates. The study of optimal fractional factorial designs under the minimum aberration and related model robustness criteria has received significant attention over the last two decades, motivated by the wide applicability of these designs. Mukerjee and Wu (2006) provided a survey which concentrates on the so-called regular designs that arise through a defining equation. Recent work in this area has focused on nonregular designs and an excellent review is available in Xu et al. (2009). Both regular and nonregular two-level factorial designs have been of particular interest because of their popularity among practitioners. The most common approach to analysing data from a regular two-level fractional factorial design is to estimate a set of contrasts (main effects and interaction effects) as described by Box and Hunter (1961). For a 2m design these contrasts form an orthogonal set and thus the estimates are independent – in this paper we will refer to this as the orthogonal parameterization. For a regular fraction of a 2m design (i.e. a 2m
p design) any two effects are either orthogonal to each other or completely aliased with each other. The aberration criteria are often used to evaluate the extent of aliasing in a given 2m
p design. In general it is desirable to have as little aliasing as possible between low-order effects and the minimum aberration design achieves this goal. A less common but in some situations a more natural approach is the baseline parameterization. A baseline parameterization is natural if there is a null state or baseline level for each factor. For instance, in a toxicological study involving binary

n

Corresponding author. Tel.: þ64 9 373 7599x85053; fax: þ64 9 373 7018. E-mail addresses: [email protected] (P. Li), [email protected], [email protected] (A. Miller), [email protected] (B. Tang). 1 Present address: Syreon, 260 – 1401 West 8th Ave, Vancouver, BC, Canada V6H 1C9.

http://dx.doi.org/10.1016/j.jspi.2014.02.009 0378-3758 & 2014 Elsevier B.V. All rights reserved.

P. Li et al. / Journal of Statistical Planning and Inference 149 (2014) 172–182

173

factors that indicate the presence or the absence of a particular toxin, the state of absence can be regarded as a natural baseline level for each factor (Kerr, 2006). A baseline level need not strictly mean the zero level on some scale, but can also refer to a standard or control level. Thus, as noted by Banerjee and Mukerjee (2008), in an industrial experiment on possible quality improvement via a change in the settings of several machines, the current settings may constitute the baseline levels. Similarly, in an agricultural experiment to investigate possible improvement in productivity by changing the doses of several fertilizers, the currently used doses may be considered as baseline levels. Under the baseline parameterization, Mukerjee and Tang (2012) developed a design theory of optimality and robustness for estimating main effects. They first showed that orthogonal arrays of strength two are universally optimal for estimating main effects. If interactions cannot be safely ignored, consequently causing bias in the estimation of main effects, Mukerjee and Tang (2012) proposed an aberration criterion that aims at minimizing this bias in a sequential fashion. Through a complete search, they found all minimum aberration baseline designs for up to 16 runs. The present paper has two goals. The first is to provide an expository account of the theory of baseline designs in a manner that is friendly to practitioners. Such effort should help promote the applications of the ideas and methods from the baseline parameterization, which we think is more appropriate than the orthogonal parameterization in many practical situations. Sections 2 and 3 and the description of the complete search algorithm in Section 4 constitute a user-friendly review of the foundational work of Mukerjee and Tang (2012). The second goal is to find minimum aberration baseline designs for larger run sizes. For designs of 20 runs, we obtain minimum aberration baseline designs for up to 13 factors by a complete search of all possible designs. In Section 4, we develop a novel incomplete search algorithm that enabled us to find good designs under the K-aberration criterion for 14 or more factors. By comparing the results from the complete search to those from our incomplete search algorithm, we conclude that the algorithm performs quite well in identifying designs that are either optimal or close to optimal with respect to the K-aberration criterion. Section 5 contains a catalogue of the minimum K-aberration 20-run designs that were found using these search algorithms. In Section 6, we establish lower bounds for values of K2 and K3 for 20-run designs – the methodology we developed to do this can be easily extended to other run sizes. Comparing the values of K3 for the designs found by our incomplete search algorithm to these lower bounds provides further evidence that the algorithm performs well. Section 7 contains our concluding remarks. 2. The baseline parameterization and the minimum K-aberration criterion Consider an N-run two-level design in m factors where for each factor, one level has been designated as the baseline level (b) and the other as the test level (t). For the baseline parameterization of the linear model, the main effect of each factor is defined by a vector of length N which has a zero in the ith position if that factor was set to the baseline level for run i and a one if it was set to the test level. The interaction between two factors is represented by a vector of length N which has a zero in the ith position if either (or both) of the factors is at its baseline level and a one only if both factors are at their test levels. In general the interaction between g of the factors is a vector that has a zero in the ith position if any of the factors involved is at its baseline level and a one only if all of the factors are at their test levels. Note that any interaction vector can be easily obtained by multiplying the main effect vectors for the factors involved in an element-wise fashion – i.e. the ith entry for an interaction vector is simply the product of the ith entries for the relevant main effect vectors. The set of baseline effects for a 23 design is given in Table 1. Now let θ0 represent the intercept coefficient and θA ; θB ; …; θABC represent the coefficients for the baseline effects specified in Table 1. Then the expected values of the response for a baseline model applied to the 23 design are given in Table 2. The interpretation of the baseline coefficients can be surmised by examining this table. Main effects represent the difference in the expected response when one factor is changed from the baseline level to the test level while all other factors are held at the baseline level – e.g. θA ¼ μtbb
μbbb . A two-factor interaction is the difference between the expected response when two factors are set to their test levels and the expected response that would be obtained using an additive main effects model for these two factors while all other factors are held at the baseline level – e.g. θAB ¼ μttb
μtbb
μbtb þ μbbb . Notice that these interpretations are analogous to those for the coefficients of the orthogonal model with one vital difference: for the baseline model factors not involved in an effect are fixed at the baseline level whereas for the orthogonal model effects Table 1 Baseline effects for a 23 design. Factor

Baseline effects

A

B

C

A

B

C

AB

AC

BC

ABC

b t b t b t b t

b b t t b b t t

b b b b t t t t

0 1 0 1 0 1 0 1

0 0 1 1 0 0 1 1

0 0 0 0 1 1 1 1

0 0 0 1 0 0 0 1

0 0 0 0 0 1 0 1

0 0 0 0 0 0 1 1

0 0 0 0 0 0 0 1

174

P. Li et al. / Journal of Statistical Planning and Inference 149 (2014) 172–182

Table 2 Expected response for a 23 design using the baseline model. Factor

Expected response for the baseline model

A

B

C

b t b t b t b t

b b t t b b t t

b b b b t t t t

μbbb ¼ θ0 μtbb ¼ θ0 þ θA μbtb ¼ θ0 þ θB μttb ¼ θ0 þ θA þ θB þ θAB μbbt ¼ θ0 þ θC μtbt ¼ θ0 þ θA þ θC þ θAC μbtt ¼ θ0 þ θB þ θC þ θBC μttt ¼ θ0 þ θA þ θB þ θC þ θAB þ θAC þ θBC þ θABC

are defined as being averaged over the remaining factors. For example, the orthogonal model main effect of A is the difference in the expected response when A is changed from the low level to the high level averaged over the levels of all other factors. The question of under what circumstance might the baseline model be preferred to the (more commonly used) orthogonal model arises. Consider a situation where the baseline levels represent the current settings for the factors of a particular process and the test levels are being investigated as avenues of possible improvement. Now suppose that for reasons of continuity, the practitioners are reluctant to make wholesale changes to the process and thus wish to identify a small number of changes that will result in the largest improvement. In such a case, the baseline effects are more relevant than the orthogonal effects as they are defined in terms of holding the remaining factors at their baseline levels. We now turn our attention to the main effects only baseline model which will be the focus of the rest of this paper. The main effects only baseline model for m factors can be represented in matrix form as where ϵ  s2 IN :

Y ¼ Xθ þ ϵ

ð1Þ

The model matrix X ¼ ½1N Z is an N  ðm þ 1Þ matrix where the first column is all ones and each of the remaining columns is a main effect vector for one of the factors. The first element of the parameter vector θ represents the intercept for the baseline model and each of the remaining elements represents the main effect coefficient for one of the factors. The estimated parameter vector θ^ ¼ ½ðXt XÞ
1 Xt Y has covariance matrix s2 ðXt XÞ
1 . For example, if the design represented in Table 1 is used for a main effects only model, then we have 2 3 1 0 0 0 61 1 0 07 6 7 2 3 6 7 2 3 θ0 61 0 1 07 0:50
0:25
0:25
0:25 6 7 6 7 6 61 1 1 07 6 θA 7
0:25 0:50 0:00 0:00 7 7 6 7 7; ðXt XÞ
1 ¼ 6 X¼6 6 7 7; θ ¼ 6 6 7 4
0:25 61 0 0 17 0:00 0:50 0:00 5 4 θB 5 6 7 61 1 0 17 θC
0:25 0:00 0:00 0:50 6 7 6 7 41 0 1 15 1

1

1

1

If the main effects only baseline model is correct (i.e. all interactions are negligible), then θ^ ¼ ½ðXt XÞ
1 Xt Y is an unbiased estimator of θ. The question arises of whether this is the best possible design. Mukerjee and Tang (2012) show that for the model described in (1), if Z forms a two-symbol orthogonal array of strength two – i.e. for any set of two columns the four level combinations (0, 0), (0, 1), (1, 0) and (1, 1) occur with the same frequency – then the associated design is universally optimal among all N-run designs, for inference on θ. Therefore provided that all interactions are negligible, the 23 factorial design is optimal. However, if this is not the case then it is necessary to evaluate the impact of the active interactions on the estimates of the main effect coefficients. The following approach was developed by Mukerjee and Tang (2012). In the presence of active interactions, true model is given by Y ¼ Xθ þ Z2 θ2 þ⋯ þZm θm þϵ;

ð2Þ

where Zj is the model matrix for the j-factor interactions. By applying least square estimation for the main effects model to this true model, the expected value of θ^ is ^ ¼ ðXt XÞ
1 Xt ðXθ þZ2 θ2 þ ⋯ þZm θm Þ EðθÞ m

¼ θþ ∑ ðXt XÞ
1 Xt Zj θj : 2

The contribution to the bias in θ^ attributable to active j-factor interactions is ðXt XÞ
1 Xt Zj θj . The first row of the matrix ðXt XÞ
1 Xt Zj determines the bias in the intercept estimate and each of the remaining rows determines the bias in one of the main effect estimates. Since the intercept term is usually of little concern, we drop the first row and denote the matrix

P. Li et al. / Journal of Statistical Planning and Inference 149 (2014) 172–182

175

Table 3 Baseline effects for an eight-run one-factor-at-a-time design. Factor

Baseline effects

A

B

C

A

B

C

AB

AC

BC

ABC

b b t t b b b b

b b b b t t b b

b b b b b b t t

0 0 1 1 0 0 0 0

0 0 0 0 1 1 0 0

0 0 0 0 0 0 1 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

formed by the remaining m rows as Bj . Thus the bias in the estimated main effect coefficients induced by active j-factor interactions is determined by Bj θj . Going back to the 23 design represented in Table 1: 2

3

0 60 6 6 60 6 61 6 Z2 ¼ 6 60 6 60 6 6 40

0

0

0

1 0

07 7 7 07 7 07 7 7; 07 7 07 7 7 15

1

1

1

0 0 0

2 3 0 607 6 7 6 7 607 6 7 607 6 7 Z 3 ¼ 6 7; 607 6 7 607 6 7 6 7 405

2

0:5 6 B2 θ2 ¼ 4 0:5 0

0:5 0 0:5

3 θAB 76 θ 7 0:5 54 AC 5; θBC 0:5 0

32

2

3 0:25 6 7 B3 θ3 ¼ 4 0:25 5θABC : 0:25

1

Thus the bias in θ^ A is 0:5θAB þ 0:5θAC þ 0:25θABC the bias in θ^ B is 0:5θAB þ0:5θBC þ 0:25θABC and the bias in θ^ C is 0:5θAC þ0:5θBC þ 0:25θABC . It is possible to construct an eight-run design for three factors such that the main effects only baseline model can be estimated and that the main effect estimates are unbiased even if some or all of the interactions are active. This would be done by using a one-factor-at-a-time design – the runs and baseline effects for such a design are given in Table 3. It is easy to see that for this design both B2 and B3 are zero matrices. For this design, the covariance matrix for the estimated coefficients is given by 2

0:5

6
0:5 6 ðXt XÞ
1 s2 ¼ 6 4
0:5
0:5


0:5


0:5


0:5

3

1:0

0:5

0:5

1:0

0:5 7 7 2 7s : 0:5 5

0:5

0:5

1:0

Therefore, although this design results in unbiased main effect estimates in the presence of active interactions, the variances of those estimates are twice as large as for the 23 design. Clearly a trade-off has to be made between the bias caused by active interactions and the variance of the estimates. As the values of the interaction coefficients are not known, it is not possible to evaluate which design would produce main effect estimates with less mean square error. This suggests two possible approaches to finding optimal baseline designs. One approach would be to restrict the search to “no bias” maineffects-only designs and from these to select the design that optimized some measure of the overall precision of the estimated main effects. A second approach would be to restrict the search to main-effects-only designs that are optimal with respect to some measure of the overall precision of the estimated main effects and from these to select a design that minimizes the impact of active interactions. Mukerjee and Tang (2012) argued that given the well known empirical principles of effect hierarchy (interaction effects are less likely to be non-negligible than main effects) and effect sparsity (only a fraction of the interaction effects will be active) the second approach will often be preferable. In this paper the second approach is adopted following the minimum K-aberration procedure developed by Mukerjee and Tang (2012). First the set of designs is restricted to two-symbol orthogonal arrays of strength two – as mentioned earlier such designs are universally optimal for inference on θ for the main-effects-only baseline model. We use OA(N, m, 2, 2) to denote a two-symbol orthogonal array of strength two with N runs and m factors. Now the goal is to minimize the bias caused by B2 θ2 þB3 θ3 þ⋯ þBm θm . As the θj 's are not known, it is sensible to use a design that minimizes the overall size of the Bj 's. Mukerjee and Tang (2012) propose using K j ¼ trðBj Btj Þ as such a measure – note that this trace is equal to the sum of the squared elements of Bj . Applying the usual effect hierarchy principle (Wu and Hamada, 2009, pp. 172–3), which states that interactions of the same order are equally likely to be active and that lower order interactions are more likely to be

176

P. Li et al. / Journal of Statistical Planning and Inference 149 (2014) 172–182

active than higher order ones, we look for a design that sequentially minimizes K2, K3, …, Km. Formally we adopt the following definition: Definition 1. Consider two N-run designs d1 and d2, both given by orthogonal arrays of strength 2. Let s be the smallest integer s for which the values of Ks for d1 and d2 differ. If the value of K s for d1 is less than that for d2, then d1 is said to have less K-aberration than d2. A minimum K-aberration design is one such that no other design has less K-aberration than it. Note that a consequence of restricting designs to orthogonal arrays of strength 2 is that the run size of such designs is limited to multiples of four: N¼ 4t where t is a positive integer. 3. Isomorphism for baseline designs An interesting issue arises in considering isomorphism among designs when using baseline parameterization. The standard definition of isomorphism allows for the relabelling of levels but under baseline parameterization the control level and the test level are not interchangeable. As a result two designs which are isomorphic using the standard definition (and therefore have equivalent properties relative to the orthogonal parameterization) may have different properties relative to the baseline parameterization. Thus when considering baseline designs it is necessary to alter the definition of isomorphism as follows: Definition 2. Two designs are isomorphic if the array Z of one can be obtained from that of the other by row and column permutations. This definition is suitable for baseline designs since K2, …, Km remain unaltered under row and column permutations of Z but not necessarily under symbol interchange in one or more columns. In the following sections it will be necessary to also consider orthogonal arrays which are nonisomorphic using the standard definition. In order to distinguish between the two types of isomorphism we adopt the following definition: Definition 3. Two orthogonal arrays are combinatorially isomorphic if one can be obtained from the other by row and columns permutations and the renaming of symbols in one or more columns. Obviously, the two definitions are closely related with the definition of combinatorially isomorphic allowing label interchanges in addition to all of the operations allowed by the definition of isomorphic. As a result, any two designs that are isomorphic must also be combinatorially isomorphic but that the converse is not true. It follows that the complete set of nonisomorphic designs S 1 is at least as large as the complete set of combinatorially nonisomorphic designs S 2 and that each element of S 1 is combinatorially isomorphic to exactly one element of S 2 . Thus if S 2 is available for a particular scenario then S 1 can be generated as follows: 1. Select one design from S 2 and switch the symbols for a subset of the columns. Do this for every possible subset of columns to create a total of 2m designs. Some of these designs may be isomorphic to others, so eliminate designs as necessary to create the largest possible set of nonisomorphic designs. 2. Repeat step 1 for each of the other designs in S 2 . 3. The union of the sets of nonisomorphic designs created in steps 1 and 2 will be S 1 .

4. Identifying minimum K-aberration designs In this paper, we are interested in finding the minimum K-aberration design for N runs and m factors. The definition of K-aberration explicitly requires that the design be an orthogonal array of strength 2 denoted by OA(N, m, 2, 2) which restricts N to multiples of four. Mukerjee and Tang (2012) have provided a method of constructing the minimum Kaberration designs for m ¼ N
1 and m ¼N
2. First consider m¼N
1. Start with a Hadamard matrix that has a column of þ1's and delete this column. Then switch the signs of the elements of any column that has a first entry of þ 1 so that the entries of first row of the resulting matrix are all
1's. Then replace all
1's in the entire array with 0's. The resulting design will be the minimum K-aberration design for m ¼N
1. For m¼N
2 apply the same procedure except delete the column of þ1's plus any one additional column. However, extending this approach to m o N
2 does not appear to be possible – Mukerjee and Tang (2012) stated that it is difficult to extend the arguments underlying their procedures to general m. As no construction methods are currently known for such designs it is necessary to use search procedures. In order to have a complete search for the minimum K-aberration design for N runs and m factors, it is necessary and sufficient to evaluate a set of designs that contains all of the nonisomorphic baseline designs of that size which are OA(N, m, 2, 2)'s. Mukerjee and Tang (2012) suggested an efficient method of conducting a complete search in situations where the complete set of combinatorially nonisomorphic OA(N, m, 2, 2)'s is available. Suppose that for one such OA, the symbols are switched in a subset of the columns to create a new OA. Although the new array is combinatorially isomorphic to the original it may constitute a nonisomorphic baseline design. Now consider the set of all possible baseline designs that can be created in this manner (including the designs where no

P. Li et al. / Journal of Statistical Planning and Inference 149 (2014) 172–182

177

switching occurs). If there are p combinatorially nonisomorphic OA's in the initial set, then the generated set contains p  2m baseline designs. This new set must contain all of the nonisomorphic baseline designs and thus it contains the baseline design that sequentially minimizes the K2, K3, … sequence. However, it is not necessary to evaluate the full sequence for all p  2m designs in this set since the K2 value is not affected by symbol interchanges – see Mukerjee and Tang (2012). Therefore we can first obtain the value of K2 for the initial set of p OA's and only retain those that give the minimum value of K2. Suppose there are q such arrays, then for each of these generate new designs by symbol interchange and evaluate the values of K3, K4, … for this set of q  2m designs. It follows that the design(s) that sequentially minimize K3, K4, … must be the minimum K-aberration design(s). This procedure is summarized by the following algorithm: 1. Generate all of the combinatorially nonisomorphic OA(N, m, 2, 2)'s for given N and m using the symbols 1 and 0, and denote the total number of such arrays as p. 2. Compute K2 for each of the p arrays. Find the arrays with minimum K2, and let q be the number of such arrays. 3. For each of these minimum K2 arrays consider generating a new array by interchanging the 0's and 1's in a subset of the m columns. As there are 2m possible subsets (including the empty set and the complete set), 2m arrays are generated for each minimum K2 array creating q  2m arrays in total. 4. For each of these q  2m arrays, compute the sequence K3, K4, … and identify those arrays that sequentially minimize this sequence. These are the MA designs under baseline parameterization. This algorithm requires that K2 is evaluated p times and K3, K4, … are evaluated q  2m times. In situations where q  2m is too large to make this feasible, we propose a modified algorithm that is not complete but has a high probability of finding the minimum K-aberration design. Again, assume that the complete set of combinatorially nonisomorphic OA(N, m, 2, 2)'s is known. These arrays are generated and their K2 values are found. The subset of these arrays that have values of K2 equal to the minimum value is retained and used as the initial arrays for a series of systematic searches. In each case, the values of K3, K4, … are found for the initial array and for the set of all arrays that can be obtained by interchanging symbols in just one column of the initial array. If the array that sequentially minimizes the K3, K4, … sequence is not the initial array, it replaces the initial array as the current best array (CBA) and the process is repeated. When the point is reached where no further improvement can be made by interchanging symbols in just one column, the values of K3, K4, … are found for all of the arrays that can be created by interchanging the symbols in exactly two of the columns. If the array that sequentially minimizes the K3, K4, … sequence for this set has less K-aberration than the CBA, then it replaces the CBA and the search process starts again (i.e. go back to considering one column symbol interchanges). This procedure is used for each of the q “K2 minimum arrays” and the K3, K4, … sequences of the resulting arrays are compared to identify the one with the least K-aberration. A step by step description of this algorithm follows: 1. Generate all of the combinatorially nonisomorphic OA(N, m, 2, 2)'s for given N and m using the symbols 1 and 0, and denote the total number of such arrays as p. 2. Compute K2 for each of the p arrays. Find the arrays with minimum K2, and let q be the number of such arrays. 3. For each of these minimum K2 arrays, designate it as the CBA (current best array) and do the following search procedure: (a) Compute the sequence K3, K4, … for the CBA. (b) Generate m new arrays by interchanging the 0's and 1's in each column of the CBA, one column at a time and compute the K3, K4, … sequences for these arrays. If the CBA has less K-aberration than any of these designs go to step 3(c). Otherwise replace the CBA with the new array which has the least K-aberration and repeat step 3(a). If there is more than one such design, randomly choose one of these to replace the CBA. (c) For the CBA, create (m choose 2) new arrays by interchanging the 0's and 1's in all possible subsets of columns of size 2. Compute the K3, K4, … sequences for these arrays. If the CBA has less K-aberration than any of these designs go to step 4. Otherwise, replace the CBA with the new array which has the least K-aberration and go to back to step 3(a). If there is more than one such design, randomly choose one of these to replace the CBA. 4. For each of q arrays described in step 2, the search procedure in step 3 finds one array. Compare the K3, K4, … sequences for these arrays and the array that has the least aberration is the best K-aberration design by the search. For a more detailed description of the implementation of our search algorithm see Li (2011). Finally, note that this algorithm is an example of the local search algorithms described by Aarts and Lenstra (2003) which are heuristic but are known to perform surprisingly well. 5. A catalogue of minimum K-aberration 20-run designs Computer searches were conducted to identify the optimal 20-run baseline designs (with respect to K-aberration) for m ¼5 to m ¼19 factors. For m ¼5 to m ¼13, the complete search method described in the previous section was applied and thus the identified designs are minimum K-aberration designs. For m ¼14 to m ¼17, the second algorithm (incomplete search) was used and as a result these designs cannot be guaranteed to be minimum K-aberration designs. For m¼18 and

178

P. Li et al. / Journal of Statistical Planning and Inference 149 (2014) 172–182

m ¼19, the construction method of Mukerjee and Tang (2012) was applied and so these designs must be minimum K-aberration designs. Table 4 contains the K2, K3, … Km sequences for the best baseline designs (in terms of K-aberration) found by the computer searches. Tables 5, 6 and 7 contain the 11-factor, 13-factor and 19-factor designs respectively. The designs for the remaining numbers of factors can be obtained by taking subsets of columns from these designs. The following table summarizes the subsets to be used in each case. Number of factors (m)

Table

Columns

5 6 7 8 9 10 11 12 13 14 15 16 17 18 19

5 5 5 5 5 5 5 6 6 7 7 7 7 7 7

1–5 1–6 1–7 1–8 1–9 1–10 All 1–12 All 1–3, 5, 6, 8–14, 17, 18 1, 2, 6, 8–19 1–12, 15–18 1–17 1–18 All

In order to assess the effectiveness of the incomplete search algorithm, we apply it to the m ¼13 case and compare the results with information gathered from the complete search. For m¼13 there are 730 combinatorially nonisomorphic arrays of strength 2. Of these five produce baseline designs that have the minimum value of K2 which is 50.94. Thus both the complete search and the incomplete search are restricted to these five arrays – call these the base arrays. For the complete search, each base array generates 213 ¼8192 baseline designs by symbol permutations and thus a total of 40,960 designs are generated. Fig. 1 contains a histogram of the values of K3 for these designs. Note that the minimum K3 value is 90.47. For the incomplete search, each of the base arrays is used in turn as the starting array for the algorithm. If we do this the K3 values for the final designs found by the five searches are 91.60, 91.70, 91.51, 91.69 and 91.25. Therefore the overall best design found has K3 ¼91.25. This is the sixth smallest unique value of K3 from the set in Fig. 1 and corresponds to approximately the 0.06 percentile of that set. To get a better indication of the efficiency of the incomplete search algorithm, this process was repeated using randomly selected starting designs. Starting designs were restricted to designs which had a row of zeros, since we noted that such designs tend to have smaller K3 values. Note that all the arrays in our catalogue of nonisomorphic arrays of strength 2 had this property so this restriction is consistent with incomplete searches that were actually conducted. For each of the five base arrays, a starting design was generated by randomly selecting one of the rows and using symbol-switching to set that row to all zeros. The incomplete search algorithm was applied to each starting design and the minimum K3 design from the five final designs was selected. This procedure was repeated 200 times. The best possible value of K3 ¼90.47 was obtained in 14% of these searches. The highest K3 value obtained was 91.6 which corresponds to the 0.1 percentile of the K3 values in Fig. 1 – note the vertical dashed lined in this plot. Thus we conclude that although the incomplete search algorithm may not find the best design, it always produces very good designs.

Table 4 K2, K3, … Km sequences for the best baseline designs found by the computer searches – the values of Ki for 12r ir m are equal to 0 for all i and m Z 12. m

K2

K3

K4

K5

K6

K7

K8

K9

K10

K11

5 6 7 8 9 10 11 12 13 14 15 16 17 18 19

5.30 8.10 11.55 15.68 20.52 26.10 33.65 41.52 50.94 61.22 72.39 84.48 98.00 112.50 128.25

1.54 3.58 6.75 11.52 18.36 29.88 46.00 64.80 90.47 122.04 160.07 205.28 260.00 324.00 399.00

0.25 0.66 1.75 4.00 8.10 18.72 35.66 58.92 95.93 147.22 212.62 297.92 413.00 558.00 741.00

0.05 0.06 0.21 0.64 1.62 6.30 16.40 33.12 65.76 119.86 187.19 282.72 432.90 633.60 906.30

0.00 0.00 0.00 0.00 0.90 4.18 10.56 28.60 68.60 112.50 179.20 309.40 491.40 758.10

0.00 0.00 0.00 0.00 0.44 1.44 7.15 28.00 45.75 73.60 149.60 259.20 433.20

0.00 0.00 0.00 0.00 0.00 0.78 8.12 12.15 17.76 46.75 89.10 162.45

0.00 0.00 0.00 0.00 0.00 1.54 1.95 1.92 8.50 18.00 36.10

0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 1.62 3.61

0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00

P. Li et al. / Journal of Statistical Planning and Inference 149 (2014) 172–182

179

Table 5 Matrix for 11-factor 20-run design. 1

2

3

4

5

6

7

8

9

10

11

0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1

1 1 1 1 1 0 0 0 0 0 1 1 1 1 1 0 0 0 0 0

0 0 1 1 1 0 1 0 1 0 1 1 0 0 0 1 0 0 1 1

0 1 1 1 0 0 0 1 1 0 1 0 1 0 0 1 1 1 0 0

0 1 0 0 1 0 1 1 1 0 1 0 0 1 1 0 1 0 0 1

1 0 0 1 1 0 0 1 1 0 0 0 1 1 0 0 0 1 1 1

1 1 1 0 0 0 1 0 1 0 0 1 0 1 0 0 1 1 1 0

1 1 0 1 0 0 1 1 0 0 0 1 0 0 1 1 0 1 0 1

1 0 1 0 1 0 1 1 0 0 0 0 1 0 1 1 1 0 1 0

0 0 1 1 0 1 1 1 0 0 1 0 0 1 1 0 0 1 1 0

1 0 0 0 0 1 1 1 1 0 1 1 1 1 0 1 0 0 0 0

Table 6 Matrix for 13-factor 20-run design. 1

2

3

4

5

6

7

8

9

10

11

12

13

1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 1 1 1 1 1 0 0 0 0 0 1 1 1 1 1

1 0 0 1 1 0 0 1 1 0 0 0 1 1 0 0 0 1 1 1

1 1 1 0 0 0 0 1 0 1 0 0 1 1 0 1 1 1 0 0

1 0 0 1 0 1 0 0 1 1 1 0 1 0 1 0 1 1 0 0

0 0 0 1 1 0 1 1 0 1 1 0 0 1 1 1 0 1 0 0

1 0 1 1 0 1 0 0 0 1 1 0 0 1 0 1 0 0 1 1

0 0 1 0 1 0 1 0 1 1 1 0 1 1 0 0 1 0 0 1

0 1 1 1 0 0 1 0 1 0 1 0 1 0 0 1 0 1 1 0

1 0 1 0 1 0 0 1 1 0 1 0 0 0 1 1 1 0 1 0

0 0 1 1 0 1 1 1 0 0 0 0 1 1 1 0 1 0 1 0

1 1 0 0 1 1 1 0 0 0 1 0 0 1 0 0 1 1 1 0

0 1 1 1 1 0 0 0 0 1 0 0 0 0 1 0 1 1 1 1

6. Lower bounds on K2 and K3 In this section, we present results concerning the minimum values of K2 and K3 that can be obtained for designs of size N ¼20. First we need to introduce some notation. For an m-column design, let Ωs represent the collection of all s-tuples g1, g2, …, gs where 1 rg 1 o g 2 o⋯ o g s r m. Further let αðg 1 g 2 …g s Þ denote the number of rows that consist entirely of ones in the s-column subarray defined by taking columns g1, g2, …, gs from the original design. Mukerjee and Tang (2012) proved the following result:   1 ð3Þ K 2 ¼ mðm
1Þ þ 48=N 2 λ where λ ¼ ∑fαðg 1 g 2 g 3 Þ
ðN=8Þg2 : 4 Ω3 For N ¼ 20; fαðg 1 g 2 g 3 Þ
ðN=8Þg2 Z 0:25 and the equality clearly only holds when αðg 1 g 2 g 3 Þ ¼ 2 or 3. Thus for given m:     m 1 K 2 Z mðm
1Þ þ 48=400  0:25: 4 3 Clearly, the lower bounds obtained by this method are not necessarily sharp (i.e. attainable) as there is no guarantee that an m-column design exists with αðg 1 g 2 g 3 Þ ¼ 2 or 3 for every 3-tuple of columns in Ω3. In Table 8, the lower limits for K2 given

180

P. Li et al. / Journal of Statistical Planning and Inference 149 (2014) 172–182

Table 7 Matrix for 19-factor 20-run design. 1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0

1 1 1 1 1 0 0 0 0 0 1 1 1 1 1 0 0 0 0 0

1 0 0 0 0 1 1 1 1 0 1 1 1 1 0 1 0 0 0 0

0 1 1 0 0 1 1 0 0 1 1 1 0 0 1 1 1 0 0 0

1 1 1 0 0 0 0 1 0 1 0 0 1 1 0 1 1 1 0 0

1 0 0 1 0 1 0 0 1 1 1 0 1 0 1 0 1 1 0 0

0 1 0 1 1 1 0 0 1 0 0 1 0 1 0 1 1 1 0 0

0 0 0 1 1 0 1 1 0 1 0 1 1 0 1 1 0 1 0 0

0 0 1 0 1 0 1 1 1 0 1 0 0 1 1 0 1 1 0 0

0 0 1 1 0 0 1 0 1 1 0 1 1 1 0 0 1 0 1 0

0 0 1 1 1 1 0 1 0 0 1 0 1 0 0 1 1 0 1 0

0 1 0 0 1 1 1 0 0 1 1 0 1 1 0 0 0 1 1 0

0 1 0 1 0 0 0 1 1 1 1 0 0 1 1 1 0 0 1 0

0 1 1 0 0 1 0 1 1 0 0 1 1 0 1 0 0 1 1 0

1 0 1 0 1 0 0 0 1 1 1 1 0 0 0 1 0 1 1 0

1 0 1 1 0 1 1 0 0 0 0 0 0 1 1 1 0 1 1 0

1 0 0 0 1 1 0 1 0 1 0 1 0 1 1 0 1 0 1 0

1 1 0 1 0 0 1 1 0 0 1 1 0 0 0 0 1 1 1 0

1 1 0 0 1 0 1 0 1 0 0 0 1 0 1 1 1 0 1 0

10000

Frequency

8000

6000

4000

91.6

2000

0 90

95

100

105

Values of K3

Fig. 1. Histogram of K3 values for the 40,960 designs obtained by symbol permutations of the five base designs for m¼ 13.

by this method (in the row labelled “Limits A”) are compared to the values found by our search. These limits are good for low values of m but get progressively worse as m increases. Alternatively if the minimum value of K2 that is attainable for an (m
1)-column design is known, then a lower limit for an m-column design can be found which is generally tighter than the previously described limits. Let λm represent the value of λ for an m-column design and consider the values of λ for the m subarrays of dimensions N  ðm
1Þ for this design. ð2Þ ðmÞ Denote these values of λ by λð1Þ m
1 , λm
1 …λm
1 . Then λm ¼

ð1Þ ð2Þ ðmÞ λm
1 þ λm
1 þ⋯ þλm
1 m
3

since every element of Ω3 for the N  m array occurs in Ω3 for m
3 of the subarrays. Clearly, the best possible situation is to have each subarray have a value of λ that corresponds to the minimum value of K2 for m
1 columns and let this value be represented by λnm
1 . Thus we have λm Z

m n λ : m
3 m
1

P. Li et al. / Journal of Statistical Planning and Inference 149 (2014) 172–182

181

Table 8 Limits for K2. m

Search K2 Limits A Limits B

5

6

7

8

9

10

11

12

5.30 5.30

8.10 8.10 8.10

11.55 11.55 11.55

15.68 15.68 15.68

20.52 20.52 20.52

26.10 26.10 26.10

33.65 32.45 32.45

41.52 39.60 41.20

m

Search K2 Limits A Limits B

13

14

15

16

17

18

19

50.94 47.58 50.08

61.22 56.42 60.70

72.39 66.15 72.15

84.48 76.80 84.48

98.00 88.40 97.73

112.50 100.98 112.50

128.25 114.57 128.25

9

10

11

12

Table 9 Limits for K3. m 6

7 3.58 3.12

Search K3 Limits

8 6.75 6.69

11.52 11.52

18.36 18.36

29.88 27.77

46.00 43.24

64.80 64.40

m

Search K3 Limits

13

14

15

16

17

18

19

90.47 88.16

122.04 120.29

160.07 159.10

205.28 205.19

260.00 259.24

324.00 324.00

399.00 399.00

A lower bound of K2 for m columns can be obtained by substituting this expression for λ into (3). Again these bounds are not necessarily sharp since an N  m design where every N  ðm
1Þ subarray has the minimum possible value of K2 may not exist. The lower limits given by this “add-one-column” method are given in Table 8 in the row labelled “Limits B.” It is evident that this procedure generally produces bounds that are tighter than the previous method. For the present situation (N ¼20), lower bounds on K2 are not of practical importance since the search procedure has ensured that designs with the smallest possible values of K2 have been identified. However, it would be straight forward to adapt this methodology to find lower bounds for higher values of N where this may not be the case (i.e. situations where a complete set of combinatorially nonisomorphic OA is not available). In such situations, the bounds could be used to give an indication as to when an m-column design has been found with an optimal or near optimal value of K2. The add-one-column methodology can be extended to establish lower bounds for K3. In this case the following result from Mukerjee and Tang (2012) is needed: K3 ¼

4 N2

ð3T 13 þ T 23 Þ

T 13 ¼ ∑fαðg 1 g 2 g 3 Þg2 Ω3

where and

T 23 ¼ ∑½f2αðg 1 g 2 g 3 g 4 Þ
αðg 1 g 2 g 3 Þg2 þf2αðg 1 g 2 g 3 g 4 Þ
αðg 1 g 2 g 4 Þg2 Ω4

þ f2αðg 1 g 2 g 3 g 4 Þ
αðg 1 g 3 g 4 Þg2 þ f2αðg 1 g 2 g 3 g 4 Þ
αðg 2 g 3 g 4 Þg2 :

ð4Þ

Again we note that the best case scenario for an m-column array is to have every N  ðm
1Þ subarray give the minimum possible value of K3. Let T n13 and T n23 be the values of T13 and T23 for an N  ðm
1Þ that has this minimum value of K3. Then it is straightforward to show that for the m-column array T 13 Z

m Tn m
3 13

and

T 23 Z

m Tn : m
4 23

Substituting these expressions into (4) produces a lower bound for K3. Table 9 compares the limits obtained using this methodology to the values of K3 for the designs found by our search algorithms. The table indicates that generally the designs identified by the search algorithms have K3 values that are quite close to the calculated lower bounds. For m r13, these values are not of practical importance as a complete search was done which ensures that designs with the lowest possible values of K3 were found. However, for m Z14 the values of K3 for the designs found by the incomplete search

182

P. Li et al. / Journal of Statistical Planning and Inference 149 (2014) 172–182

algorithm are all reasonably close to the calculated bounds. Keeping in mind that the bounds are not necessarily sharp, this provides some additional evidence that the incomplete search algorithm has performed well. 7. Concluding remarks In this paper, we have provided a catalogue of baseline optimal designs for experiments of 20 runs for m ¼5–19 factors – baseline optimal designs for smaller run sizes are provided by Mukerjee and Tang (2012). A complete search algorithm was employed to find the minimum K-aberration designs for m r 13. For m Z14, we developed a search algorithm and demonstrate that it is very effective at finding designs that are near optimal in terms of K-aberration. Clearly, this algorithm can be easily extended to search for designs of larger run size. We have also developed methodology for establishing lower bounds on K2 and K3 which can also be easily extended for designs of more than 20 runs. Acknowledgements Boxin Tang's research is supported by the Natural Sciences and Engineering Research Council of Canada. The authors would like to thank two referees for their useful comments. References Aarts, E., Lenstra, J.K., 2003. Local Search in Combinatorial Optimization. Princeton University Press, New Jersey. Banerjee, T., Mukerjee, R., 2008. Optimal factorial designs for cDNA microarray experiments. Ann. Appl. Statist. 2, 366–385. Box, G.E.P., Hunter, K.P., 1961. The 2k
p fractional factorial designs. Technometrics 3, 311–351. Kerr, K.F., 2006. Efficient 2k factorial designs for blocks of size 2 with microarray applications. J. Qual. Technol. 38, 309–318. Li, P., 2011. Algorithmic Search of Baseline Minimum Aberration Designs (M.Sc. thesis). Department of Statistics and Actuarial Science, Simon Fraser University, Burnaby, British Columbia, Canada. Mukerjee, R., Tang, B., 2012. Optimal fractions of two-level factorials under a baseline parameterization. Biometrika 99, 71–84. Mukerjee, R., Wu, C.F.J., 2006. A Modern Theory of Factorial Designs. Springer, New York. Wu, C.F.J., Hamada, M.S., 2009. Experiments: Planning, Analysis and Optimization, 2nd ed. Wiley, Hoboken, New Jersey. Xu, H., Phoa, F.K.H., Wong, W.K., 2009. Recent developments in nonregular fractional factorial designs. Statist. Surveys 3, 18–46.