Properties of designs

This section shows the structural and statistical properties of the orthogonal arrays, conference designs and D-efficient designs generated by the Orthogonal Array package. The properties of the arrays and designs are calculated using the array_link object or functions from the package.

Definitions of arrays and designs

Before introducing the structural and statistical properties, we define orthogonal arrays, conference designs and D-efficient designs:

An orthogonal array (OA) of strength \({t}\), \({N}\) runs and \({n}\) factors at \({s}\) levels is an \({N}\times {n}\) array of symbols \(0, \ldots,({s}-1)\), such that for every subset of \({t}\) columns, every \({t}\)-tuple occurs equally often [Rao47]. The set of all strength-\({t}\) OAs with \({N}\) runs and \({n}\) factors at \({s}\) levels is denoted by \({\operatorname{OA}({N}; {t}; {s}^{n})}\). If \({s=2}\), the OA is called a two-level OA and the set of all strength-\({t}\) two-level OAs with \({N}\) runs and \({n}\) factors is denoted as \({\operatorname{OA}({N}; {t}; {2}^{n})}\).

For \({N}\) even, a conference design [SEG19] \(C\) is an \({N}\times {n}\) array which satisfies \({C}^{T}C = (n-1) I_{n}\), with \({C}_{ii} = 0\) and \({C}_{ij} \in \{-1,1\}\), for \({i} \neq {j}\) and \({i}, {j} = 1, \ldots, n\). A \({N}\times {N}\) conference design \(E\) such that \(E{E}^{T} = (n-1) I_{n}\) is called a conference matrix; see [EN95], [CD06] and [XLB12].

A D-optimal design [DAT07] (\(X\)) is an \({N}\times {n}\) array which maximizes the \(D\)-efficiency, defined as \({(\operatorname{det}({X}^{T}_{M}{X}^{\phantom{T}}_{M})^{1/p})/N}\), for a given \({N}\times {p}\) model matrix \({X}_{M}\) (for details see Model matrices). The Orthogonal Array package uses a coordinate-exchange algorithm to generate designs that optimize the \(D\)-efficiency. Since there is no guarantee that the resulting designs have the largest possible \(D\)-efficiency, we refer to them as D-efficient designs in this documentation. An orthogonal array is called D-optimal orthogonal array if it provides the largest \(D\)-efficiency among all comparable orthogonal arrays.

Structural properties of an array

The OApackage can calculate the rank of an array, defined as the maximum number of linearly independent column or row vectors in the array. The rank of an array is useful for several other functions in the package. For two-level arrays, the OApackage can also check if the arrays are foldover arrays. A two-level array is called a foldover array if half of its runs are mirror images of the other half, in the sense that the factor levels are changed from 0 to 1 and from 1 to 0.

For example, to calculate the rank of a two-level orthogonal array and determine whether the array is a foldover array, one can use the methods array_link::rank() and array_link::foldover():

Calculate rank of array and test for foldover

>>> array = oapackage.exampleArray(1) # Select an example two-level orthogonal array
>>> array.showarray() # Show the two-level orthogonal array
array:
  0   0   0   0   0
  0   0   0   0   0
  0   0   0   1   1
  0   0   1   0   1
  0   1   0   1   0
  0   1   1   0   0
  0   1   1   1   1
  0   1   1   1   1
  1   0   0   1   1
  1   0   1   0   1
  1   0   1   1   0
  1   0   1   1   0
  1   1   0   0   1
  1   1   0   0   1
  1   1   0   1   0
  1   1   1   0   0
 >>> print(array.rank()) # Calculate the rank of the array
 5
 >>> print(array.foldover()) # Determine if the array is foldover
 False

Other structural properties such as whether an array involves two levels or is symetric can be found in the documentation of array_link, which shows the full set of methods available.

Model matrices

For orthogonal arrays and conference designs, the OApackage can calculate different model matrices. The model matrices available depend on the type of array and design.

An example on how to generate an interaction model matrix for a two-level orthogonal array is shown below.

Calculate interaction effects model matrix

>>> array=oapackage.exampleArray(0,1)
exampleArray 0: array in OA(8,2, 2^2)
>>> array.showarray()
array:
  0   0
  0   0
  0   1
  0   1
  1   0
  1   0
  1   1
  1   1
>>> M=oapackage.array2modelmatrix(array, 'i')
>>> print(M)
[[ 1. -1. -1.  1.]
 [ 1. -1. -1.  1.]
 [ 1. -1.  1. -1.]
 [ 1. -1.  1. -1.]
 [ 1.  1. -1. -1.]
 [ 1.  1. -1. -1.]
 [ 1.  1.  1.  1.]
 [ 1.  1.  1.  1.]]

Statistical properties of orthogonal arrays

Orthogonal arrays are commonly evaluated in terms of their generalized wordlength pattern [XW01] (GWLP). Two-level OAs are also commonly evaluated in terms of their \({J}_{k}\)-characteristics and \(F\)-vectors [DT99]. The OApackage can calculate all these statistical criteria: array_link::GWLP(), array_link::Fvalues(), array_link::Jcharacteristics().

The following example shows how to calculate the GWLP, \({F}_{k}\)-values and \({J}_{k}\)-characteristics from an array_link object:

Calculate GWLP and F-values

>>> al=oapackage.exampleArray(1,1) # Select an example array
exampleArray 1: array 3 in OA(16, 2, 2^5)
>>> gwlp = al.GWLP() # Calculate its generalized word length pattern
>>> print('GWLP: %s'% str(gwlp) )
GWLP: (1.0, 0.0, 0.0, 1.0, 1.0, 0.0)
>>> print('F3-value: %s' % str(al.Fvalues(3))) # Calculate the F_3 values
F3-value: (4, 6)
>>> print('F4-value: %s' % str(al.Fvalues(4))) # Calculate the F_3 values
F4-value: (1, 4)
>>> print('J3-characteristics: %s' % str(al.Jcharacteristics(3))) # Calculate the J_3-characteristics
J3-characteristics: (-8, -8, 0, 0, 0, -8, 0, -8, 0, 0)

We now briefly mention some technical details of the \({J}_{k}\)-characteristics, the \({F}_{k}\)-values and the GWLP.

Optimality criteria for D-efficient designs

In [ES17], D-efficient designs for the model including the intercept, all main effects and all two-factor interactions are generated. The OApackage provides functionality to compute the optimality criteria used to generate the D-efficient designs in [ES17]. Moreover, the package can calculate the well-known \(A\)- and \(E\)-optimality criteria from the literature on Optimal Experimental Design [DAT07]. The functions to perform the calulcations are array_link::Defficiency(), array_link::DsEfficiency(), array_link::Aefficiency(), array_link::Eefficiency().

The following example shows how to calculate the \(D\)-, \({D}_{s}\)-, \(A\)- and \(E\)-efficiency for a design that permits the estimation of the interaction model.

Calculate optimality criteria for D-efficient designs

# Select an array that can estimate the interaction model
>>> al = oapackage.exampleArray(11, 1)
exampleArray 11: D-optimal array in OA(44, 2^8)
>>> print('D-efficiency: %.4f' % al.Defficiency())
D-efficiency: 0.8879
>>> print('Ds-efficiency (Eendebak and Schoen, 2017): %.4f' % al.DsEfficiency())
Ds-efficiency (Eendebak and Schoen, 2017): 0.8059
>>> print('A-efficiency for the interaction model: %.4f' % al.Aefficiency())
A-efficiency for the interaction model: 0.7906
>>> print('E-efficiency for the interaction model: %.4f' % al.Eefficiency())
E-efficiency for the interaction model: 0.3602

Projection capacities

Other relevant statistical criteria to evaluate a two-level design with \(N\) runs and \(k\) factors include the so-called projection estimation capacity (PEC) and projection information capacity (PIC) [LST07]. These criteria focus on the projections of the two-level design onto a smaller number of factors. More specifically, the PEC and PIC summarize the performance of all the \(N\)-run subdesigns with \(l \leq k\) factors in terms of the capacity to estimate the interaction model and the \(D\)-efficiency for this model, respectively.

The PEC and PIC are based on the so-called PEC and PIC sequences, which are formally defined as follows. Let \(PEC_{l}\) denote the proportion of \(N\)-run \(l\)-factor subdesigns that permit the estimation of the interaction model in \(l\) factors, that is, the model including the intercept, all \(l\) main effects and all \(l(l-1)/2\) two-factor interactions. The PEC sequence is the vector \((PEC_{1}, PEC_{2}, \ldots, PEC_{k})\). Now, let \(PIC_{l}\) denote the average \(D\)-efficiency for the interaction model in \(l\) factors accross all \(N\)-run \(l\)-factor subdesigns. The PIC sequence is the vector \((PIC_{1}, PIC_{2}, \ldots, PIC_{k})\). The OApackage can calculate the PEC and PIC sequences of two-level designs with PECsequence() and PICsequence(), respectively.

The following example shows how to compute the PEC and PIC sequences of a two-level orthogonal array using the OApackage.

Calculate the PEC and PIC sequences

>>> al=oapackage.exampleArray(1,1)
exampleArray 1: array 3 in OA(16, 2, 2^5)
>>> PEC = al.PECsequence()
>>> print('PEC sequence: %s'% ','.join(['%.2f' % x for x in PEC]) )
PEC sequence: 1.00,1.00,1.00,0.80,0.00
>>> PIC = al.PICsequence()
>>> print('PIC sequence: %s'% ','.join(['%.2f' % x for x in PIC]) )
PIC sequence: 1.00,1.00,0.95,0.66,0.00

Properties of conference designs

In [SEG19], it is shown that the \(F_4\) vector is useful for classifying definitive screening designs [XLB12] that are generated by folding over a conference design. To calculate the \(F_4\) vector, we first need to compute the \(J_4\)-characteristics of the conference design. The calculations for the \(J_k\)-characteristics of conference designs are similar as for orthogonal arrays; see Statistical properties of orthogonal arrays. Consider a definitive screening design constructed from an \(N\)-run conference design with at least four factors. The \(F_4\) vector of this design collects the frequencies of the \(J_4\)-characteristics of \({2N − 8\lambda}\) for \({\lambda = 1,\ldots,N/4}\) when \(N\) is a multiple of 4, or \({\lambda = 1, \ldots, (N − 2)/4}\) when \(N\) is an odd multiple of 2.

Note: the \(J_4\)-values of a definitive screening design generated by folding over a conference design are twice the value of the \(J_4\)-values of the conference design. The \(F_4\) vector of a conference design and the corresponding definitive screening design are equal.

Calculate the F4 vector for a conference design

>>> import oapackage
>>> array=oapackage.exampleArray(47,1)
exampleArray 47: third conference design in C(20,8)
>>> F4=array.FvaluesConference(4)
>>> print(F4)
(0, 2, 4, 51, 13)
>>> definitive_screening_design = oapackage.conference2DSD(array)

The individual \(J_k\)-characteristics can be calculated with the method Jcharacteristics_conference(). For conference designs, we can calculate the projection statistics using conferenceProjectionStatistics().

Calculate projection statistics for conference designs

>>> array = oapackage.exampleArray(46, 1)
exampleArray 46: second conference design in C(20,8)
>>> pec, pic, ppc = oapackage.conference.conferenceProjectionStatistics(array)
>>> print('Projection estimation capacity for 4 columns: %.3f'  % pec)
Projection estimation capacity for 4 columns: 0.986
>>> J3 = oapackage.Jcharacteristics_conference(array, number_of_columns = 3)