Minimum-Size Mixed-Level Orthogonal Fractional Factorial Designs Generation: A SAS-Based Algorithm

Orthogonal fractional factorial designs (OFFDs) are frequently used in many fields of application, including medicine, engineering and agriculture. In this paper we present a software tool for the generation of minimum size OFFDs. The software, that has been implemented in SAS/IML, puts neither a restriction on the number of levels of each factor nor on the orthogonality constraints. The algorithm is based on the joint use of polynomial counting function and complex coding of levels and follows the approach presented in Fontana (2013).


Introduction
In this paper we present a software tool for generating minimum size orthogonal fractional factorial designs (OFFDs).The algorithm is based on the approach described in Fontana (2013).A preliminary version of it has been developed in the unpublished master's degree thesis Sampò (2011).
In this Section we give an overview of the use of OFFDs in practical applications and of the problems related to their generation.In Section 2 we briefly review the algebraic theory of OFFDs based on polynomial counting functions and strata.In Section 3 and Section 4 we describe the algorithm.Some applications of the algorithm are presented in Section 5. Finally, concluding remarks are given in Section 6. Section 1 and Section 2 are closely based on Section 1 and Section 2 of Fontana (2013).We include them here for completeness.
OFFDs are frequently used in many fields of application, including medicine, engineering and agriculture.They offer a valuable tool for dealing with problems where there are many factors involved and each run is expensive.They also keep the statistical analysis of the data quite simple.The literature on the subject is extremely rich.A non-exhaustive list of references, mainly related to the theory of the design of experiments, includes the fundamental paper of Bose (1947) and the following books: Raktoe, Hedayat, and Federer (1981), Collombier (1996), Dey and Mukerjee (1999), Wu and Hamada (2000), Mukerjee and Wu (2006) and Bailey (2008).
Orthogonal arrays (OAs) represent an important class of OFFDs, see, for example, Hedayat, Sloane, and Stufken (1999) and Schoen, Eendebak, and Nguyen (2010).Indeed an OA of appropriate strength can be used to solve the wide range of problems related to the study of the size of the main effects and interactions up to a given order of interest.
It is evident that in many real-life experiments, finding an OFFD with the smallest possible number of runs is of great importance.This is particularly true in the case where the cost of each experiment is high in terms of resources and/or time, such as in the study of the relationship between fuel consumption and the design parameters of a new car engine.
A large number of techniques are known to generate OFFDs and, in particular, OAs.For example: The case where all factors have the same number p of levels and p is a prime number or a power of a prime number, is commonly studied using Galois Fields and finite geometries.
Hadamard matrices are used for OAs where all factors have 2 levels and where strength is less than or equal to 3.
Difference schemes are a tool for constructing mixed OAs of strength 2.
From the above it is clear that there are several different methods covering different situations.When different methods are applied to certain problems the solutions that are found can be significantly different.For example minimum size OAs with eleven 2-level factors and strength 2 obtained using the Galois field GF (2) ≡ Z 2 have 16 runs while those obtained using Hadamard matrices have 12 runs.Thus the problem of finding minimum size OFFDs can be difficult for the non-expert user due to the difficulty of selecting the most appropriate method.
The joint use of polynomial indicator functions and complex coding of levels provides a general theory for mixed level OFFDs (see Pistone and Rogantin 2008).This theory puts neither a restriction on the number of levels of each factor nor on the orthogonality constraints.It also makes use of commutative algebra (see Pistone and Wynn 1996), and generalizes the approach to two-level designs as discussed in Fontana, Pistone, and Rogantin (2000).The definition of strata provided in Section 2 makes it possible to transform each OFFD into a solution of a homogeneous system of linear equations where the unknowns are positive integers.
The aim of this work is to develop a general software tool for minimum size OFFDs generation, where general means that it puts neither a restriction on the number of levels of each factor nor on the orthogonality constraints.

Algebraic generation of OAs
In this Section, for ease in reference, we report some relevant results of the algebraic theory of OFFDs from Fontana (2013).The interested reader can find further information, including the proofs of the propositions in Fontana (2013) itself and in Fontana et al. (2000), Pistone and Rogantin (2008), Fontana and Pistone (2010a) and Fontana and Pistone (2010b).

Fractions of a full factorial design
Let us consider an experiment which includes m factors D j , j = 1, . . ., m.Let us code the n j levels of the factor D j by the n j roots of the unity where ω gives the number of times it belongs to the multiset F.
We recall that the underlying set of elements F * is the subset of D that contains all the elements of D that appear in F at least once.We denote the number of elements of a fraction F by #F, with #F = ζ∈F * f * (ζ).
Example 2.1.Let us consider m = 1, n 1 = 3.We get In order to use polynomials to represent all the functions defined over D, including multiplicity functions, we define The function X j is called simple term or, by abuse of terminology, factor.
, the monomial function The function X α is called interaction term.
We observe that {X is a basis of all the complex functions defined over D. We use this basis to represent the counting function of a fraction according to Definition 2.2.
Definition 2.2.The counting function R of a fraction F is a complex polynomial defined over D so that for each ζ ∈ D, R(ζ) equals the number of appearances of ζ in the fraction.A 0/1 valued counting function is called an indicator function of a single replicate fraction F. We denote by c α the coefficients of the representation of R on D using the monomial basis With Proposition 2.1 from Pistone and Rogantin (2008), we link the orthogonality of two interaction terms with the coefficients of the polynomial representation of the counting function.
Proposition 2.1.If F is a fraction of a full factorial design D, R = α∈L c α X α is its counting function and [α − β] is the m-tuple made by the componentwise difference in the rings 3. the terms X α and X β are orthogonal on F, if and only if, c [α−β] = 0.
We now define projectivity and, in particular, its relation with OAs.Given I = {i 1 , . . ., i k } ⊂ {1, . . ., m} , i 1 < . . .< i k we define the projection π I as Definition 2.3.A fraction F factorially projects onto the I-factors, I = {i 1 , . . ., i k } ⊂ {1, . . ., m}, i 1 < . . .< i k , if the projection π I (F) is a multiple full factorial design, i.e., the multiset (D i 1 ×. ..×D i k , f * ) where the multiplicity function f * is constantly equal to a positive integer over It follows that F projects on the first factor and does not project on the second factor.
Definition 2.4.A fraction F is a mixed orthogonal array of strength t if it factorially projects onto any I-factors with #I = t.Proposition 2.2.A fraction factorially projects onto the I-factors, I = {i 1 , . . ., i k } ⊂ {1, . . ., m}, i 1 < . . .< i k , if and only if, all the coefficients of the counting function involving the I-factors only are 0. Proposition 2.2 can be immediately stated for mixed OAs.
Proposition 2.3.A fraction is an orthogonal array of strength t, if and only if, all the coefficients c α of the counting function up to the order t are 0.

Counting functions and strata
It follows from Propositions 2.1, 2.2 and 2.3 that the problem of finding OFFDs can be written as a polynomial system in which the indeterminates are the complex coefficients c α of the counting polynomial fraction.
Let us now introduce a different way to describe the full factorial design D and all its subsets.We consider the indicator functions 1 ζ of all the single points of D. The counting function R of a fraction F can be written as  Pistone and Rogantin (2008), we recall two basic properties which hold true for the full design D. Note that the notation gcd stands for greatest common divisor and lcm stands for least common multiple.
Proposition 2.4.Let X j be the simple term with level set 1.Over D, the term X r j takes all the values of Ω s j equally often, where s j = 1 if r = 0 and 2. Over D, the term X α takes all the values of Ω sα equally often, where s α = lcm(s 1 , . . ., s m ) and s i , that is determined according to the previous Item 1, corresponds to We refer to Ω sα as the level set of X α .Sometimes we also write s in place of s α to simplify the notation.Let us now define the strata that are associated with simple and interaction terms.
, where ω (sα) h ∈ Ω sα and s α is determined according to the previous Proposition 2.4.We use n α,h to denote the number of points of the fraction F that are in the stratum h , where s is determined by X α according to Proposition 2.4.
We now use a part of Proposition 3 from Pistone and Rogantin (2008) to obtain conditions on n α,h that make X α centered on the fraction F.
Proposition 2.6.Let X α be a term with level set Ω s on full design D and let P s (ζ) be the complex polynomial associated with the sequence (n α,h ) h=0,...,s−1 so that P s (ζ) = s−1 h=0 n α,h ζ h and Φ s the cyclotomic polynomial of the s-roots of the unity.1.Let s be prime.The term X α is centered on the fraction F, if and only if, its s levels appear equally often n α,0 = n α,1 = . . .= n α,s−1 .

Let s
The term X α is centered on the fraction F, if and only if, the remainder H s (ζ) = P s (ζ) mod Φ s (ζ), whose coefficients are integer linear combinations of n α,h , h = 0, . . ., s − 1, is identically zero.
If we recall that n α,h are related to the values of the counting function R of a fraction F by n α,h = ζ∈D α h y ζ , Proposition 2.6 allows us to express the condition X α is centered on F as integer linear combinations of the values R(ζ) of the counting function over the full design D. In Section 2.3, we will show the use of this property to generate fractional factorial designs.

Generation of fractions
We use strata to generate fractions that satisfy a given set of constraints on the coefficients of their counting functions.Formally, we give the following definition.
The set of all the fractions of D whose counting functions are C-compatible is denoted by OF (n 1 . . .n m , C).
For example let us consider OA(n, s m , t), i.e., OAs with n rows and m columns, where each column has s symbols, s prime and with strength t.Using Proposition 2.3 the coefficients of the corresponding counting functions must satisfy the conditions c α = 0 for all α ∈ C where C = {α ∈ L ≡ (Z s ) m : 0 < α ≤ t} and α is the number of non-null elements of α.It follows that OF (s m , C) = n OA(n, s m , t).Now using Proposition 2.6, we can express these conditions using strata.If we consider α ∈ C, we can write the condition c α = 0 as: .
We observe that the ((s − 1) × s) matrix B sα is the same for all α because s α is constant and equal to s.
To obtain all the conditions it is enough to vary α ∈ C. We therefore obtain the homogeneous system of linear equations AY = 0 where A is the (#C •(s−1)×s m ) matrix whose rows contain the values, over D, of the difference between two indicator functions of strata, 1 the s m column vector whose entries are the values of the counting function over D; 0 is the (#C • (s − 1)) column vector whose entries are all equal to 0.
Let us now consider the general case in which there are no restrictions on the number of levels and show our method for OA(n, 4 2 , 1).In this case the number of levels is a power of a prime, 4 = 2 2 .Using Proposition 2.3 the coefficients of the corresponding counting functions must satisfy the conditions c α = 0 for all α ∈ C where C = {α ∈ L ≡ Z 4 × Z 4 : α = 1}.Let us consider c 1,0 .From Item (1) of Proposition 2.4, X 1 takes the values in Ω s (1,0) where s (1,0) = 4. From Proposition 2.6, X 1 will be centered on F, if and only if, the remainder Lang 1965) and so we can compute the remainder Let us now consider c 2,0 .From Item (1) of Proposition 2.4, X 2 1 takes the values in Ω s (2,0) where s (2,0) = 2. From Proposition 2.6, X 2 1 will be centered on F, if and only if If we repeat the same procedure for all α, where α = 1, and if we recall that n α,h = ζ∈D α h y ζ , and so, for example, n (1,0),0 = y 1,1 + y 1,i + y 1,−1 + y 1,−i with i = √ −1, the OAs OA(n, 4 2 , 1) become the positive integer solutions of the following integer linear homogeneous system: AY = 0 where and 0 is the 10-rows column vector whose entries are all equal to 0.
It should be noted that the matrix A of the coefficients does not have full rank, e.g., the first and the fourth rows are equal.This aspect is discussed in Fontana and Pistone (2010b).In any case the solution method used here does not require a reduction to a full rank matrix.

The optimization problem
When the tests are very expensive or require a substantial amount of time, experimenters are interested in finding minimum size orthogonal fractional designs, i.e., fractional factorial designs that satisfy some orthogonality requirements and consist of the minimum number of points.
According to our formalization, the problem is equivalent to extracting one fraction F * from OF (n 1 . . .n m , C), such that the size of the fraction #F * is minimum.
The problem can be written as where A is the matrix built as explained in Section 2.3, 0 is the column vector that has the same numer of rows of A and whose entries are all equal to 0, 1 is the #D column vector whose entries are all equal to 1, 1 is its transpose and Y is a vector of positive integer numbers that contains the unknown counting function values, Y = [0, . . ., 0].

The algorithm and its SAS implementation
We now focus on the details of the algorithm that makes use of the theoretical results described in Section 2. The implementation used SAS software: in particular SAS/IML (SAS Institute, Inc. 2008a), for the construction of the constraints AY = 0 and SAS/OR (SAS Institute, Inc. 2008b), to solve the optimization problem as stated in Equation 1 of Section 2.4.
As described in the previous section, we state the problem of finding OFFDs that satisfy a set of conditions in terms of orthogonality as a homogeneous system of linear equations, AY = 0, in which the indeterminates are positive integers that contain the unknown counting values.For the sake of simplicity, we work with mixed level OAs of strength t, that is OA(n, n 1 • . . .• n m , t); the extension to mixed level OFFDs, regardless of the orthogonality constraints they must satisfy, is straightforward, as we will explain in Section 6.1.
With reference to Figure 1, the algorithm takes the levels n j for each factor, j = 1, . . ., m, and the required strength t as input; generates the constraint matrix A; performs the optimization phase; provides a minimum size OA as output.
Constraints generation is done by the following major steps: 1. We generate the full factorial Z n 1 × . . .× Z nm .It represents both the set L of multiindexes and the full factorial design with integer coding, i.e., the full factorial design where the levels of the factor D j are coded with the integer 0, . . ., n j − 1, j = 1, . . ., m.
We denote by α a generic element of L.
2. Let us denote the subset of L that contains all the α such that the norm of α is less than or equal to t, α = [0, . . ., 0] by C; for all α ∈ C ⊆ L we build the column vector X α (ζ) : ζ ∈ D and, coherently with Proposition 2.4, we also compute s α .
3. According to Proposition 2.6, for all s ∈ {2, . . ., s max }, where s max is the maximum of all the s α , α ∈ C, we compute the matrix B s .
Figure 1: The main modules of the algorithm.
(a) If s is prime, then Once the constraints AY = 0 are obtained, the optimization phase can be carried out in SAS or in another software tool like lp solve (Berkelaar, Eikland, and Peter 2004).
In the next section we explain in detail all the steps of the algorithm.

Full factorial design generation
We must build up the full factorial design To write the rows of L, we generalize the common conversion base algorithm (Cox, Little, and O'Shea 1997).We make a loop over the #L rows: the row index (from 0 to #L − 1) is the number k to be converted and r 1 . . .r m is the corresponding row of L. This is the algorithm: The notation [k] n indicates the residue of k modulo n, and int(k/n) is the integer part of the division k n .

Simple and interaction terms
From now on we use ω k,n j in place of ω to simplify the notation.We observe that the mapping τ n j is a group isomorphism between Ω n j with the complex product and Z n j with the sum modulo n j .It follows that τ n 1 × . . .× τ nm is a group isomorphism between D and L. Let us consider the term X r j and its values X r j (ω k,n j ) = ω r k,n j over D j , j = 1, . . ., m, r ∈ Z + .Using the isomorphism τ n j we can represent the vector [X r j (ω k,n j ) : ω r k,n j ∈ D j ] with the vector Now we consider the interaction term X α and its values over D, α ∈ L. From Proposition 2.4, X α takes all the values of Ω sα , where s α = lcm(s 1 , . . ., s m ) and s j , j = 1, . . ., m are determined as in Item 1 of the same Proposition.It follows that we can represent the vector [X α (ω k 1 ,n 1 , . . ., ω km,nm ) : ω k,n j ∈ D j , j = 1, . . ., m] with the vector Finally, by observing that the complex conjugate ω h,s of ω h,s ∈ Ω s is ω [s−h]s,s and that, using τ s , can be represented as [s − h] s , we build a matrix whose columns contain the integer representation of X α (ζ) for all α ∈ C and for all ζ ∈ D.
Furthermore we store all the values s α , α ∈ C in a vector and we calculate the maximum value s max of it, which will be used in the next steps of the algorithm.

Remainders
For each s, where s = 2, . . ., s max , there is one or more linear constraint over the coefficients n α,h .If s is prime we know from Proposition 2.6 that the conditions can be easily determined.For example for s = 5 the conditions correspond to the matrix B 5 : If s is not prime, then we have to generate the cyclotomic polynomial Φ s , the polynomial P s (ζ) and then to calculate the remainder H s of the division P s /Φ s .
From the fact that the roots of a cyclotomic polynomial are the primitive roots of unity (Lang 2002), we can calculate Φ s by dividing X s − 1 by the product of the cyclotomic polynomials Φ d i corresponding to the proper divisors d 1 , . . ., d s of s : .
We create a function that implements the polynomials division (Cox et al. 1997).Let f be the dividend polynomial and g the divisor, the pseudocode of the algorithm for finding q and r such that f = qg + r is: Input: g, f Output: q, r q = 0; r = f while r = 0 and LT (g) divides LT (r) do q = q + LT (r)/LT (g) r = r − (LT (r)/LT (g))g end while LT means leading term, i.e., the term of highest degree.To construct the cyclotomic polynomial Φ s we need three steps: 1. Find all proper divisors of s: d 1 , . . ., d k .From Proposition 2.6 we know that we are interested in calculating the remainder of the division of P s by the cyclotomic polynomial Φ s .In this case the dividend polynomial has symbolic coefficients n α,h , h = 0, . . ., s − 1:

Analyze each d
The division algorithm remains the same but at each step we have to store also these coefficients.We create an s×s matrix, where the column index is related to the degree and the row index to the polynomial coefficient.We explain this point considering s = 4.By definition P 4 (ζ) is and the corresponding matrix becomes )ζ and its related matrix Finally, by simply deleting the zero columns, we obtain the matrix B 4 :

Strata and constraints
The last step is to transform each matrix B s , s = 2, . . ., s max , that refers to n α,h , h = 0, . . ., s− 1 into a matrix A s that refers to y ζ , ζ ∈ D and then to stack all the A s into the matrix A.
In this way, it is possible to set up the optimization problem as stated in Equation 1 of Section 2.4.
We obtain The remainder matrix for s = 4 is By repeating the same procedure for all α ∈ C and stacking all the A s matrices, we obtain the A matrix, as reported at the end of Section 2.3.

The optimization problem
At this point our problem is to find the fraction that satisfies all the constraints while minimizing the number of points.We can use any optimization software that solves mixed-integer linear problems (MILP).In particular, we have constructed the set of constraints in such a way that can be directly used in SAS/OR and can also be exported for use with lp solve (Berkelaar et al. 2004).
SAS/OR procedures are used to optimize a linear function subject to linear and integer constraints.Specifically, we use the Optmodel procedure, which solves the general MILP of the form: where: A is an r × v matrix of technological coefficients; b is an r × 1 matrix of right-hand side constants; c is a v × 1 matrix of objective function coefficients; x is a v × 1 matrix of structural variables; l is a v × 1 matrix of lower bounds on x; u is a v × 1 matrix of upper bounds on x; S is a subset of the set of indices 1, . . ., v and identifies those variables that must take only integer values.
In our case: b is the r × 1 matrix whose entries are all equal to 0; c = 1 is the n 1 • . . .• n m × 1 matrix whose entries are all equal to 1; x is an n 1 • . . .• n m × 1 matrix that contains the unknown number of appearances of each point of D in the fraction F; l is the n 1 • . . .• n m × 1 matrix whose entries are all equal to 0; we do not need to set any value for the upper bounds u; S contains all n 1 • . . .• n m variables, because all of them must take integer values.
Finally, in order to exclude the empty fraction from the solutions, we add the following constraint: For further details about the Optmodel procedure see SAS Institute, Inc. (2008b).
As mentioned before, the constraints matrix A is stored in a dataset and can be exported to other optimization programs.Our algorithm exports the matrix A into a file that can be used by the well-known and widely-used open-source code lp solve (Berkelaar et al. 2004).

Results
We tested the algorithm by comparing our results with those reported in Schoen et al. (2010).We also ran other examples.The results are summarized in Table 1.The first column reports the factor set, i.e., the number of levels of each factor, #D is the cardinality of the full factorial design D, t and n are respectively the strength and the cardinality of the OA as generated by the algorithm.To deal with generic OFFD, it is enough, using Proposition 2.1, to modify C, the subset of the multi-indexes α ∈ L for which the coefficients of the counting function must be equal to zero and to leave the remaining part of the code unchanged.We illustrate this point with an example.Let us suppose that we consider 4 factors with 3 levels each and that we search for a minimum size OFFD such that all the main effects of the factors, D 1 , . . ., D 4 , and the interaction between factor D 2 and D 3 are estimable and orthogonal.
From Item 2 of Proposition 2.1, the main effect of the factor D i will be orthogonal to the main effect of the factor D j , i < j, if and only if, c α = 0 for α ∈ e(δ i 1 , δ i 2 , δ i 3 , δ i 4 ) + f (δ j 1 , δ j 2 , δ j 3 , δ j 4 ) with e, f = 1, 2 and δ k h = 0 if h = k and δ k h = 1 if h = k; the main effect of the factor D 1 (resp.D 4 ) will be orthogonal to the interaction between factor D 2 and D 3 , if and only if, c α = 0 for α ∈ {(e, f, g, 0)} (resp.α ∈ {(0, e, f, g)} ) with e, f, g = 1, 2; the main effect of the factor D 2 or of the factor D 3 will be orthogonal to the interaction between factor D 2 and D 3 , if and only if, c α = 0 for α ∈ {(0, e, f, 0)} with e, f = 1, 2.
It follows that C can be written as: ζ∈D y ζ 1 ζ with y ζ ≡ R(ζ) ∈ {0, 1, . . ., n, . ..}.The particular case in which R is an indicator function corresponds to y ζ ∈ {0, 1}.From Proposition 2.1 we obtain that the values of the counting function over D, y ζ , are related to the coefficients c α by determined taking the remainder H s of the division between the polynomial P s and the cyclotomic polynomial Φ s .4. Using all the X α (ζ) : ζ ∈ D , α ∈ C and the corresponding B sα and recalling that n α,h = ζ∈D α h y ζ , we build the constraints matrix A.
which represents the set of multi-indexes α.The total number of points in L is #L = n 1 •. ..•n m .The set of multi-indexes L can be represented by a matrix L with m columns and #L rows.
is not prime then we have to return to step 1 setting s = d i .3.WhenΦ d 1 , . . ., Φ d k are built, we obtain Φ s by dividing X s −1 by the product Φ d 1 •. ..•Φ d k .

Table 1 :
The results of some examples.Given the levels of each factor, n 1 , . . ., n m , and the strength t as input, the algorithm provides one of the minimum size mixed level OAs OA(n, n 1 • . . .• n m , t) as output.