Dual-Orthogonal Arrays for Order-of-Addition Two-Level Factorial Experiments

Abstract

A new class of orthogonal arrays called dual-orthogonal arrays is proposed in this paper to design order-of-addition two-level factorial experiments in which both component addition orders and component levels can be varied over treatments. Dual-orthogonal arrays can be viewed as an optimal combination of order-of-addition orthogonal arrays and two-level orthogonal arrays. Based on these two different concepts of orthogonality, when a compound model is used to fit the observed data, both pairwise order effects and component main effects can be estimated with optimal efficiency. Under the assumption of normality, these two kinds of parametric effects can also be inferred independently. A three-drug combination study is first used to show that dual-orthogonal arrays can be practical for real-world studies. Both combinatorial and computational methods are then introduced to construct dual-orthogonal arrays. Additionally, a design catalogue is generated for future work.

KEYWORDS:

• Interchange algorithm
• Main-effect plan
• Optimal design
• Screening design

1 Introduction

Screening experiments are frequently conducted in experimental studies. Under the sparsity-of-effects assumption, relatively few active effects with a significant impact on the responses can be initially identified from screening experiments. After leaving out those inactive effects, follow-up experiments can be designed to explore the treatment–response relationships more thoroughly such that experimental resources can be utilized more efficiently. This sequential experimentation strategy has been widely used in drug discovery, quality engineering and many other fields. The reader is referred to Dean and Lewis (2006) for an excellent introduction to screening experiments in diverse disciplines. Recently, a new class of experiments called order-of-addition experiments has received increasing attention from medical researchers, quality engineers and statisticians. The most remarkable feature of these kinds of experiments is that the responses may depend on material addition orders and/or process execution orders. Conventionally, these materials and processes are called the components. By deliberately varying component addition orders, several trials can be conducted to screen out active effects for the purpose of understanding how the responses depend on the addition orders of these components. Therefore, screening experiments may also play a key role in addressing order-of-addition problems.

Van Nostrand (1995) introduced the concept of pairwise order factors to connect order-of-addition designs and two-level factorial designs. Based on this key concept, Voelkel (2019) defined order-of-addition orthogonal arrays to design order-of-addition experiments. Peng, Mukerjee, and Lin (2019) further showed that the overall mean and main effects of all pairwise order factors can be estimated with optimal efficiency if the pairwise order matrix is an order-of-addition orthogonal array of strength two. On the basis of these seminal papers, several optimal and near-optimal order-of-addition designs have been constructed for real-world studies. Recent proposals include those by Chen, Mukerjee, and Lin (2020), Winker, Chen, and Lin (2020), Zhao, Lin, and Liu (2021), Tsai (2022), Wang and Mee (2022), and Zhao, Lin, and Liu (2022). These approaches only take component addition orders into account while component levels are assumed to be fixed throughout the experiments. Sometimes, changing component levels, in addition to changing component addition orders, may have a significant impact on the responses. If screening experiments are capable of evaluating the potential impact of varying component addition orders and component levels, then experimental costs would be reduced considerably. Voelkel (2019) noted that process variables can be incorporated into order-of-addition orthogonal arrays. Rios and Lin (2022) and Rios, Winker, and Lin (2022) addressed this research question under the framework of mixture experiments in which component levels of every treatment must satisfy the sum-to-one constraint. This physical constraint may not be appropriate for nonmixture experiments. Xiao et al. (2022) proposed an active learning procedure to identify optimal order-level combinations. However, their approach seems unable to screen out easy-to-interpret active effects. To solve these practical problems, the existing techniques for multifactorial experiments can be applied such that component levels can also be determined systematically. These kinds of experiments are called the order-of-addition factorial experiments because both component addition orders and component levels can be varied over treatments. By combining order-of-addition orthogonal arrays and two-level orthogonal arrays, a new class of orthogonal arrays is introduced in this paper to design such experiments. Because these orthogonal arrays have two different orthogonal properties, they are called dual-orthogonal arrays.

The remainder of this paper is organized as follows. Section 2 introduces a three-drug combination experiment as a motivating problem. Section 3 presents a compound model to fit the observed data. Dual-orthogonal arrays are defined in Section 4 and related optimality results are established. The three-drug combination experiment is analyzed in Section 5 to show that dual-orthogonal arrays can be practical. Both combinatorial and computational methods are proposed in Section 6 to construct dual-orthogonal arrays. Additionally, a design catalogue is generated for future work. Concluding remarks are given in the final section.

2 Three-Drug Combination Experiment

Combination therapies of several drugs have been widely used in medical treatments because they can enhance drug efficacy and reduce drug resistance and toxicity. Several trials are often conducted to explore optimal treatment orders and dose levels when developing new combination therapies. Efficient experimental designs can thus be used to extract as much information as possible from the observed data. To study whether or not different treatment orders of chemotherapeutics and their dose levels have a significant impact on lymphoma cell growth, Wang, Xu, and Ding (2020) conducted a series of in vitro experiments. In the first experiment, three FDA-approved chemotherapeutics for clinical trials, namely, paclitaxel, doxorubicin and mitoxantrone, denoted by 1, 2, and 3, respectively, were sequentially added into lymphoma cell cultures in different orders. In addition to varying addition orders, chemotherapeutics 1 and 2 were set at two different dose levels $\frac{1}{6}{\text{IC}}_{50}$ and $\frac{1}{8}{\text{IC}}_{50}$, denoted by + and –, to study their dosage effects, where ${\text{IC}}_{50}$ represents the half maximal inhibitory concentration. The dose level of chemotherapeutic 3 was held constant throughout the experiment. A 24-treatment design consisting of all order-level combinations of six addition orders and four level combinations was used to design this experiment. After conducting each treatment, the response was then measured by a tumor inhibition index. All observed data are collected in Table 1.

Table 1 Three-drug combination experiment.

Treatment	Time 1	Time 2	Time 3	y
1^‡	(1, +)	(2, +)	3	39.91
2^‡	(1, +)	3	(2, +)	44.38
3^†	(2, +)	(1, +)	3	17.08
4^‡	(2, +)	3	(1, +)	20.88
5^†	3	(1, +)	(2, +)	34.68
6^†	3	(2, +)	(1, +)	37.37
7^‡	(1, –)	(2,–)	3	30.00
8^†	(1, –)	3	(2,–)	47.18
9^‡	(2,–)	(1, –)	3	25.10
10^†	(2,–)	3	(1, –)	33.60
11^‡	3	(1, –)	(2,–)	35.04
12^†	3	(2,–)	(1, –)	35.04
13^†	(1, +)	(2,–)	3	44.33
14^‡	(1, +)	3	(2,–)	38.18
15^†	(2,–)	(1, +)	3	22.26
16^‡	(2,–)	3	(1, +)	31.40
17^†	3	(1, +)	(2,–)	38.91
18^‡	3	(2,–)	(1, +)	42.30
19^†	(1, –)	(2, +)	3	44.87
20^†	(1, –)	3	(2, +)	43.93
21^‡	(2, +)	(1, –)	3	26.02
22^†	(2, +)	3	(1, –)	22.56
23^‡	3	(1, –)	(2, +)	31.15
24^‡	3	(2, +)	(1, –)	37.19

^† : Fraction 1;

^‡ : Fraction 2.

From an experimental design perspective, the 24-treatment design in Table 1 can be viewed as the full design. Because the full design can provide the most comprehensive information regarding the treatment–response relationship of interest, it is a natural candidate for the experiment. However, when the number of components is moderately large, conducting all order-level combinations can be expensive and/or time-consuming. To reduce experimental costs, conducting a subset of the full design is more practical. After adopting this design strategy, an immediate question is how to select a subset of treatments from the full design such that the fitted model can be estimated with optimal efficiency? The three-drug combination experiment will be revisited in Section 5 to address this practical problem.

3 A Compound Model

Suppose that an n-treatment order-of-addition factorial experiment is conducted to study m components, where components 1 to u have two levels and the remaining components have constant levels. Let ${\mathit{t}}_h$ denote the hth treatment given by $${\mathit{t}}_h=\left[\left(c_{h1},l_{h1}\right),\left(c_{h2},l_{h2}\right),\dots ,\left(c_{hm},l_{hm}\right)\right].$$

The ith ordered pair $(c_{hi},l_{hi})$ in ${\mathit{t}}_h$ consists of component c_hi to be added in the ith step and its level l_hi, where $h=1,2,\dots ,n$ and $i=1,2,\dots ,m$. If component c_hi is set at the high level, then l_hi would be +. Otherwise, it would be –. Specifically, if component c_hi is held at a constant level throughout the experiment, then l_hi would be omitted from ${\mathit{t}}_h$ for notational simplicity. To conduct the treatment ${\mathit{t}}_h$, a practitioner must add the m components sequentially according to the addition order and level combination specified by ${\mathit{t}}_h$. After conducting the treatment ${\mathit{t}}_h$, the following linear model is proposed to fit the observed data: (1) $$y_h=\mu +\sum _{i=1}^{m-1}\sum _{j=i+1}^m{\gamma }_{ij}{z}_{h,ij}+\sum _{i=1}^u{\beta }_i{x}_{h,i}+{ϵ}_h,$$(1) where μ represents the overall mean, γ_ij denotes the pairwise order effect of components i and j, $z_{h,ij}$ represents the pairwise order factor given by $$z_{h,ij}=\{\begin{array}{cc}+1& \mathrm{ }\text{if}\mathrm{ }\text{component}\mathrm{ }i\mathrm{ }\text{is}\mathrm{ }\text{added}\mathrm{ }\text{before}\mathrm{ }\text{component}\mathrm{ }j\mathrm{ }\text{in}\mathrm{ }{\mathit{t}}_h;\\ -1& \mathrm{ }\text{if}\mathrm{ }\text{component}\mathrm{ }i\mathrm{ }\text{is}\mathrm{ }\text{added}\mathrm{ }\text{after}\mathrm{ }\text{component}\mathrm{ }j\mathrm{ }\text{in}\mathrm{ }{\mathit{t}}_h,\end{array}$$

β_i denotes the main effect of component i, $x_{h,i}$ represents the component level indicator variable given by $$x_{h,i}=\{\begin{array}{cc}+1& \mathrm{ }\text{if}\mathrm{ }\text{component}\mathrm{ }i\mathrm{ }\text{is}\mathrm{ }\text{set}\mathrm{ }\text{at}\mathrm{ }\text{the}\mathrm{ }\text{high}\mathrm{ }\text{level}\mathrm{ }\text{in}\mathrm{ }{\mathit{t}}_h;\\ -1& \mathrm{ }\text{if}\mathrm{ }\text{component}\mathrm{ }i\mathrm{ }\text{is}\mathrm{ }\text{set}\mathrm{ }\text{at}\mathrm{ }\text{the}\mathrm{ }\text{low}\mathrm{ }\text{level}\mathrm{ }\text{in}\mathrm{ }{\mathit{t}}_h,\end{array}$$and ϵ_h denotes the error term. All error terms are assumed to be independent and identically distributed normal random variables with zero mean and constant variance. Actually, the compound model in (1) is a simple combination of the pairwise order model for order-of-addition experiments and the main effects model for two-level factorial experiments. It has $p=1+q+u$ easy-to-interpret parameters in total, where $q=m(m-1)/2$.

Formally, an n-treatment order-of-addition two-level factorial design can be denoted by $$\mathit{T}=\left\{{\mathit{t}}_1,{\mathit{t}}_2,\dots ,{\mathit{t}}_n\right\}.$$

Based on the compound model in (1), after conducting all n treatments, the $n\times 1$ response vector y can be expressed in the following matrix form: $$\mathit{y}=\mu \mathbf{1}+\mathit{Z\gamma }+\mathit{X\beta }+\mathit{ϵ},$$where $\mathbf{1}$ denotes the $n\times 1$ vector of ones, Z represents the n × q pairwise order matrix, $\mathit{\gamma }$ denotes the $q\times 1$ vector of pairwise order effects, X represents the n × u component level matrix, $\mathit{\beta }$ denotes the $u\times 1$ vector of component main effects, and $\mathit{ϵ}$ represents the $n\times 1$ vector of error terms. When the least squares method is used to estimate the compound model in (1), the moment matrix M has the form $$\mathit{M}=\frac{1}{n}\left[\begin{array}{ccc}\hfill n& \hfill {\mathbf{1}}^{\top }\mathit{Z}& \hfill {\mathbf{1}}^{\top }\mathit{X}\\ \hfill {\mathit{Z}}^{\top }\mathbf{1}& \hfill {\mathit{Z}}^{\top }\mathit{Z}& \hfill {\mathit{Z}}^{\top }\mathit{X}\\ \hfill {\mathit{X}}^{\top }\mathbf{1}& \hfill {\mathit{X}}^{\top }\mathit{Z}& \hfill {\mathit{X}}^{\top }\mathit{X}\end{array}\right].$$

The moment matrix M plays a key role in characterizing optimal designs because its inverse is proportional to the variance-covariance matrix of the least squares estimators of all parameters in (1). Related design issues will be addressed in the next section.

4 Dual-Orthogonal Arrays

An n-treatment design is said to be pairwise order balanced if, for every pair of i and j, component i is added before component j in $n/2$ treatments and component i is added after component j in the remaining $n/2$ treatments. Given such a design, the two numbers + 1 and –1 occur equally often in every column of its pairwise order matrix Z, with the result that the least squares estimator of an arbitrary pairwise order effect ${\widehat{\gamma }}_{ij}$ is uncorrelated with that of the overall mean $\widehat{\mu }$. Similarly, an n-treatment design is said to be component level balanced if, for every i, component i is set at the high level in $n/2$ treatments and component i is set at the low level in the remaining $n/2$ treatments. Given such a design, the two numbers + 1 and –1 occur equally often in every column of its component level matrix X, with the result that the least squares estimator of an arbitrary component main effect ${\widehat{\beta }}_i$ is uncorrelated with that of the overall mean $\widehat{\mu }$. Additionally, the pairwise order factors are said to be orthogonal to the component level indicator variables if the inner product of an arbitrary column of Z and an arbitrary column of X is equal to zero. Given such a design, ${\widehat{\gamma }}_{ij}$ is uncorrelated with ${\widehat{\beta }}_i$ for every pair of ${\widehat{\gamma }}_{ij}$ and ${\widehat{\beta }}_i$. After imposing these symmetric properties, the moment matrix M can be expressed as a block diagonal matrix given by $$\mathit{M}=\frac{1}{n}\left[\begin{array}{ccc}n& {\mathbf{0}}^{\top }& {\mathbf{0}}^{\top }\\ \mathbf{0}& {\mathit{Z}}^{\top }\mathit{Z}& \mathbf{0}\\ \mathbf{0}& \mathbf{0}& {\mathit{X}}^{\top }\mathit{X}\end{array}\right].$$

By definition, an n × u matrix, whose entries are either + 1 or –1, is called a two-level orthogonal array of strength two, denoted by OA $(n,2^u,2)$, if, for any two columns, the ordered pairs $(-1,-1),(+1,-1)$, $(-1,+1)$, and $(+1,+1)$ occur equally often. Based on the theoretical results developed by Kiefer (1975), two-level orthogonal arrays of strength two are universally optimal for estimating the overall mean and main effects of all experimental factors. Although pairwise order factors can be used to convert an order-of-addition design to a two-level design, as noted in Lin and Peng (2019), the pairwise order matrix cannot be a two-level orthogonal array of strength two due to the transitive property of pairwise order factors. Therefore, Voelkel (2019) proposed order-of-addition orthogonal arrays as a generalization of two-level orthogonal arrays to design order-of-addition experiments. By definition, an n × q pairwise order matrix, where $q=m(m-1)/2$, is called an order-of-addition orthogonal array of strength two, denoted by OofA-OA $(n,m,2)$, if, for any two columns, the frequencies of the ordered pairs $(-1,-1),(+1,-1),(-1,+1)$, and $(+1,+1)$ are proportional to those of the full order-of-addition design. Peng, Mukerjee, and Lin (2019) further showed that an order-of-addition design is optimal for estimating the overall mean and main effects of all pairwise order factors if its pairwise order matrix is an order-of-addition orthogonal array of strength two. Next, these two different concepts of orthogonality will be combined to design order-of-addition two-level factorial experiments.

4.1 Definition

Suppose that an n-treatment order-of-addition factorial experiment is conducted to study m components, where components 1 to u have two levels and the remaining components have constant levels. When the compound model in (1) is used, a matrix pair $(\mathit{Z},\mathit{X})$ can be determined accordingly.

Definition 1.

A matrix pair $(\mathit{Z},\mathit{X})$ is called a dual-orthogonal array of strength two, denoted by DOA $(n,m,2^u,2)$, if it satisfies the following conditions.

1. The pairwise order matrix Z is an OofA-OA $(n,m,2)$.

2. The component level matrix X is an OA $(n,2^u,2)$.

3. The ordered pairs $(-1,-1),(+1,-1)$, $(-1,+1)$, and $(+1,+1)$ occur equally often in any $n\times 2$ matrix consisting of an arbitrary column vector of Z and an arbitrary column vector of X.

Specifically, when u = 1, that is, there is only one two-level component, condition (B) of Definition 1 requires the $n\times 1$ component level matrix X to be a two-level orthogonal array of strength one, denoted by OA $(n,2^1,1)$.

Example 1.

A fractional design called Fraction 1 is determined by the 12 treatments marked with $†$ in Table 1. Their pairwise order factors and component level indicator variables are given in the supplementary materials. It is straightforward to see that the pairwise order matrix Z is an OofA-OA $(12,3,2)$ and the component level matrix X is an OA $(12,2^2,2)$. Both conditions (A) and (B) of Definition 1 are satisfied. The ordered pairs $(-1,-1),(+1,-1),(-1,+1)$, and $(+1,+1)$ occur equally often in any 12 × 2 matrix consisting of an arbitrary column vector of Z and an arbitrary column vector of X. Therefore, condition (C) of Definition 1 is also satisfied. By definition, the matrix pair $(\mathit{Z},\mathit{X})$ determined by Fraction 1 is a DOA $(12,3,2^2,2)$. Similarly, another fractional design called Fraction 2 is determined by the 12 treatments marked with ‡ in Table 1. The matrix pair $(\mathit{Z},\mathit{X})$ determined by Fraction 2 is also a DOA $(12,3,2^2,2)$. Additionally, the full design in Table 1 determines a DOA $(24,3,2^2,2)$.

4.2 Optimality

As shown in Example 1, dual-orthogonal arrays of strength two do exist. Next, $\varphi$-optimal designs are introduced before presenting related optimality results. Peng, Mukerjee, and Lin (2019) noted that a design is said to be $\varphi$-optimal if its moment matrix M has maximum $\varphi (\mathit{M})$ for every concave and signed permutation invariant optimality criterion $\varphi (·)$. Many alphabetic-optimality criteria with statistically meaningful justifications, such as A-, D- and E-optimality criteria, are included in this important class of optimality criteria.

Theorem 1.

An order-of-addition two-level factorial design is $\varphi$-optimal for estimating the compound model in (1) if the corresponding matrix pair $(\mathit{Z},\mathit{X})$ is a DOA $(n,m,2^u,2)$.

The proof of Theorem 1 is presented in the supplementary materials. Because many statistical software and libraries provide built-in functions to compute the determinant of a matrix, the D-optimality criterion will be used later to search for dual-orthogonal arrays. Specifically, a design is said to be D-optimal if its moment matrix has maximum determinant among all competing designs.

5 Three-Drug Combination Experiment Revisited

The three-drug combination experiment is analyzed to illustrate dual-orthogonal arrays.

5.1 Main Effects Analysis Using Full Data

Because the full design in Table 1 determines a DOA $(24,3,2^2,2)$, when the compound model in (1) is used, the moment matrix ${\mathit{M}}_0$ is a block diagonal matrix given by $${\mathit{M}}_0=\frac{1}{24}\left[\begin{array}{cccccc}\hfill 24& \hfill 0& \hfill 0& \hfill 0& \hfill 0& \hfill 0\\ \hfill 0& \hfill 24& \hfill 8& \hfill -8& \hfill 0& \hfill 0\\ \hfill 0& \hfill 8& \hfill 24& \hfill 8& \hfill 0& \hfill 0\\ \hfill 0& \hfill -8& \hfill 8& \hfill 24& \hfill 0& \hfill 0\\ \hfill 0& \hfill 0& \hfill 0& \hfill 0& \hfill 24& \hfill 0\\ \hfill 0& \hfill 0& \hfill 0& \hfill 0& \hfill 0& \hfill 24\end{array}\right].$$

The block diagonal structure ensures that the model sum of squares can be further decomposed as (2) $$\begin{array}{c}\text{SSM}\left({\gamma }_{12},{\gamma }_{13},{\gamma }_{23},{\beta }_1,{\beta }_2\right)=\text{SSM}\left({\gamma }_{12},{\gamma }_{13},{\gamma }_{23}\right)\\ + \text{SSM}\left({\beta }_1,{\beta }_2\right),\end{array}$$(2) where $\text{SSM}(·)$ denotes the model sum of squares when all parameters in the parentheses are included in the fitted model. Based on this decomposition, the null hypothesis $H_0:{\gamma }_{12}={\gamma }_{13}={\gamma }_{23}=0$ can be tested independently of the null hypothesis $H_0:{\beta }_1={\beta }_2=0$. After conducting the standard F-tests, their p-values are equal to 0.004 and 0.780, respectively. The pairwise order effects are statistically significant at the 5% level and the component main effects are not significant. These findings are similar to those discovered by Wang, Xu, and Ding (2020). To further identify active pairwise order effects, seven candidate models consisting of all combinations of these three pairwise order effects are evaluated using Mallows’s C_p criterion. Their adjusted coefficients of determination are also reported for reference purposes. As shown in Table 2, where the two symbols · and ∘ are used to indicate that a pairwise order effect is included and not included in a candidate model, Mallows’s C_p ranks the pairwise order model M5 as the optimal model. This model also has the largest adjusted R-squared, where the least squares estimates of γ₁₂ and γ₂₃ are equal to 4.030 and –3.129, respectively. Based on these two estimates, one might conclude that chemotherapeutic 1 should be administered before chemotherapeutic 2 and chemotherapeutic 3 should also be administered before chemotherapeutic 2. Accordingly, the two treatment orders 132 and 312 are expected to yield better mean efficacy of tumor cell inhibition. These two orders will be further evaluated in Section 5.3.

Table 2 Goodness-of-fit measures for analyzing pairwise order effects.

		Parameter			Full design		Fraction 1		Fraction 2
Model	p	γ12	γ13	γ23	Cp	$R_{\text{adj}}^2$	Cp	$R_{\text{adj}}^2$	Cp	$R_{\text{adj}}^2$
M1	2	·	∘	∘	5.409	0.340	1.047	0.531	5.911	0.101
M2	2	∘	·	∘	19.742	– 0.032	12.735	– 0.074	8.950	– 0.095
M3	2	∘	∘	·	8.722	0.254	8.710	0.135	1.680	0.375
M4	3	·	·	∘	7.065	0.318	2.816	0.493	7.803	0.009
M5	3	·	∘	·	2.381	0.446	2.001	0.540	2.845	0.365
M6	3	∘	·	·	6.808	0.325	8.443	0.169	2.110	0.418
M7	4	·	·	·	4.000	0.429	4.000	0.482	4.000	0.354

5.2 Main Effects Analysis Using Fractional Data

As noted in Example 1, Fraction 1 determines a DOA $(12,3,2^2,2)$. When the compound model in (1) is used, the moment matrix ${\mathit{M}}_1$ is a block diagonal matrix given by $${\mathit{M}}_1=\frac{1}{12}\left[\begin{array}{cccccc}\hfill 12& \hfill 0& \hfill 0& \hfill 0& \hfill 0& \hfill 0\\ \hfill 0& \hfill 12& \hfill 4& \hfill -4& \hfill 0& \hfill 0\\ \hfill 0& \hfill 4& \hfill 12& \hfill 4& \hfill 0& \hfill 0\\ \hfill 0& \hfill -4& \hfill 4& \hfill 12& \hfill 0& \hfill 0\\ \hfill 0& \hfill 0& \hfill 0& \hfill 0& \hfill 12& \hfill 0\\ \hfill 0& \hfill 0& \hfill 0& \hfill 0& \hfill 0& \hfill 12\end{array}\right].$$

The decomposition in (2) can also be applied when analyzing the observed data of Fraction 1. The p-values of the standard F-tests for evaluating $H_0:{\gamma }_{12}={\gamma }_{13}={\gamma }_{23}=0$ and $H_0:{\beta }_1={\beta }_2=0$ are equal to 0.049 and 0.332, respectively. After ranking all candidate models in Table 2, the pairwise order model M5 is still identified as the optimal model. Because Fraction 2 also determines a DOA $(12,3,2^2,2)$, its moment matrix, denoted by ${\mathit{M}}_2$, is identical to ${\mathit{M}}_1$, with the result that the decomposition in (2) can also be applied. After conducting the standard F-tests, the p-values for evaluating $H_0:{\gamma }_{12}={\gamma }_{13}={\gamma }_{23}=0$ and $H_0:{\beta }_1={\beta }_2=0$ are equal to 0.097 and 0.308, respectively. As shown in Table 2, goodness-of-fit measures for the pairwise order models M3, M5, M6, and M7 are very close. Overall, M5 still appears to be a promising model. In this example, the fractional and full data provide either the same or very similar results. Most interestingly, the moment matrices of these two fractional designs are identical to the moment matrix of the full design, that is, ${\mathit{M}}_1={\mathit{M}}_2={\mathit{M}}_0$. These two fractional designs have the same symmetric properties as the full design under the compound model in (1). Furthermore, based on Theorem 1, all parameters in (1) can be estimated with optimal efficiency using these two fractional designs. Therefore, conducting such fractional designs can be a cost-effective alternative under resource constraints. The problem raised in Section 2 has been addressed.

5.3 Optimal Treatments

In some drug combination studies, changing dose levels of chemotherapeutics may have a significant impact on a tumor inhibition index. This phenomenon is not observed in the three-drug combination experiment probably due to the fact that the dose levels are chosen too conservatively such that their true dosage effects are undetectable. Therefore, dosage effects will not be taken into account in the subsequent analysis. To seek a better fit, four candidate models are generated by adding all combinations of estimable two-factor interaction effects of pairwise order factors to M5. The reader is referred to Mee (2020) for a detailed discussion regarding estimable interaction effects of pairwise order factors. These candidate models are evaluated according to their adjusted coefficients of determination. The two most promising models both indicate that the treatment order 132 may yield better mean efficacy than that of 312. Related analysis results are reported in the supplementary materials. Consequently, the treatments with the form $${\mathit{t}}^{\star }=\left\{\left(1,\pm \right),3,\left(2,\pm \right)\right\}$$are declared to be optimal. In particular, chemotherapeutics 1 and 2 can be set at their most economical dose levels. The three-drug combination experiment was also analyzed by other researchers. Wang, Xu, and Ding (2020) did not recommend a treatment after analyzing the three-drug combination experiment because their conclusions were drawn based on two other follow-up experiments. Xiao et al. (2022) identified ${\mathit{t}}_8=\left\{(1,-),3,(2,-)\right\}$ in Table 1 as optimal using a 15-treatment sequential design. Most interestingly, their proposal ${\mathit{t}}_8$ is included in ${\mathit{t}}^{\star }$. Experimental designs determined by dual-orthogonal arrays of strength two are model-based one-shot designs, so they are less flexible than active learning procedures. However, these designs are $\varphi$-optimal for estimating the compound model in (1) in which all parametric effects are easy to interpret. As shown in this example, dual-orthogonal arrays can be practical for real-world studies. Several construction methods will be introduced in the next section to generate more dual-orthogonal arrays.

6 Construction Methods

Let ⊗ denote the Kronecker product. A simple but important approach is introduced first to construct dual-orthogonal arrays of strength two.

Theorem 2.

The matrix pair $({\mathit{Z}}^{\ast }\otimes {\mathbf{1}}_{n_2},{\mathbf{1}}_{n_1}\otimes {\mathit{X}}^{\ast })$ is a DOA $(n_1\times n_2,m,2^u,2)$ if ${\mathit{Z}}^{\ast }$ is an OofA-OA $(n_1,m,2)$ and ${\mathit{X}}^{\ast }$ is an OA $(n_2,2^u,2)$.

The proof of Theorem 2 is omitted in this paper since it is straightforward. Note that the 24 treatments in Table 1 were generated using the same technique as that described in Theorem 2. These kinds of dual-orthogonal arrays can be quite large when n₁ and n₂ are moderately large. Two additional methods will be proposed in this section to generate dual-orthogonal arrays with more economical run sizes.

6.1 Combinatorial Method

Voelkel (2019) noted that the ordered pairs $(-1,-1),(+1,-1),(-1,+1)$, and $(+1,+1)$ would occur equally often in two columns of an order-of-addition orthogonal array of strength two if their pairwise order factors have no common component. Suppose that an OofA-OA $(n,m+2u,2)$, denoted by W, has been generated for $m+2u$ components. To satisfy condition (A) of Definition 1, an OofA-OA $(n,m,2)$, denoted by Z, can be generated for m components by deleting those columns corresponding to extra pairwise order factors from W. Next, an OA $(n,2^u,2)$, denoted by X, can be obtained by properly assigning u extra pairwise order factors to component level indicator variables such that both conditions (B) and (C) of Definition 1 are satisfied. Specifically, Algorithm 1 is designed to implement this combinatorial approach.

Algorithm 1

Generate DOAs using extra pairwise order factors

1: input m, n, u and W

2: for $h=1,2,\dots ,n$ do

3: for $i=1,2,\dots ,m-1$ do

4: for $j=i+1,i+2,\dots ,m$ do

5: $z_{h,ij}\leftarrow w_{h,ij}$

6: end for

7: end for

8: for $i=1,2,\dots ,u$ do

9: $x_{h,i}\leftarrow w_{h,m+2i-1m+2i}$

10: end for

11: end for

12: output $(\mathit{Z},\mathit{X})$

Example 2.

Wang, Xu, and Ding (2020) reported another order-of-addition two-level factorial experiment for studying five chemotherapeutics, where chemotherapeutics 1 to 4 had constant levels and chemotherapeutic 5 had two levels. In this five-drug combination experiment, 60 treatments were conducted. Suppose that a smaller design is required for the experiment due to resource constraints. First, an OofA-OA $(24,7,2)$, denoted by W, is generated by calling the R function optFederov in the package AlgDesign developed by Wheeler (2022). The design matrix of this OofA-OA $(24,7,2)$ is presented in Table 3. Next, pairwise order factors involving components 1 to 5 are used to generate an OofA-OA $(24,5,2)$, denoted by Z. Lastly, set $x_{h,5}=w_{h,67}$ for all h to get an OA $(24,2^1,1)$, denoted by X. The matrix pair $(\mathit{Z},\mathit{X})$ is a DOA $(24,5,2^1,2)$. Its design matrix is presented in Table 3. To highlight the relationship between these two design matrices, components 6 and 7 in Table 3 are displayed in boldface.

Table 3 Design matrices of an OofA-OA $(24,7,2)$ and a DOA $(24,5,2^1,2)$.

OofA-OA design matrix							DOA design matrix
2	1	3	5	7	6	4	2	1	3	$(5,-)$	4
2	1	6	7	5	3	4	2	1	$(5,+)$	3	4
2	3	4	7	6	1	5	2	3	4	1	$(5,-)$
2	6	4	5	3	1	7	2	4	$(5,+)$	3	1
3	1	2	5	4	6	7	3	1	2	$(5,+)$	4
3	1	6	4	5	2	7	3	1	4	$(5,+)$	2
3	6	7	5	4	1	2	3	$(5,+)$	4	1	2
4	1	5	3	6	7	2	4	1	$(5,+)$	3	2
4	1	7	6	3	5	2	4	1	3	$(5,-)$	2
4	3	2	7	5	1	6	4	3	2	$(5,-)$	1
4	6	2	5	7	1	3	4	2	$(5,+)$	1	3
5	1	2	7	4	3	6	$(5,-)$	1	2	4	3
5	1	3	4	7	2	6	$(5,-)$	1	3	4	2
5	2	6	4	7	1	3	$(5,+)$	2	4	1	3
5	4	6	2	3	1	7	$(5,+)$	4	2	3	1
5	7	6	3	4	1	2	$(5,-)$	3	4	1	2
6	1	3	2	7	4	5	1	3	2	4	$(5,+)$
6	1	4	7	2	3	5	1	4	2	3	$(5,+)$
6	3	5	7	2	1	4	3	$(5,+)$	2	1	4
7	1	4	5	2	6	3	1	4	$(5,-)$	2	3
7	1	6	2	5	4	3	1	2	$(5,-)$	4	3
7	2	3	4	5	1	6	2	3	4	$(5,-)$	1
7	4	3	2	6	1	5	4	3	2	1	$(5,-)$
7	5	3	6	2	1	4	$(5,-)$	3	2	1	4

An OofA-OA $(n,m+2u,2)$ must be used as a template when Algorithm 1 is used to generate a DOA $(n,m,2^u,2)$. As far as I know, small order-of-addition orthogonal arrays of strength two for more than eight components are rare in the existing literature. Currently, Algorithm 1 is unable to generate dual-orthogonal arrays of strength two with n < 200 and u > 2 due to the lack of suitable templates. To generate more dual-orthogonal arrays, a computational method will be used in the next section.

6.2 Computational Method

Given an OofA-OA $(n,m,2)$ as in condition (A) of Definition 1, Algorithm 2 is designed to search for an OA $(n,2^u,2)$ such that both conditions (B) and (C) are satisfied. First, an OA $(n,2^1,1)$, denoted by ${\mathit{x}}_i$, is randomly generated in line 7 to determine initial levels of component i over all treatments. Next, a new proposal ${\mathit{w}}_i$ is generated in line 11 by randomly switching a pair of + 1 and –1 in ${\mathit{x}}_i$. This simple operation was also used by Nguyen (1996) and Xu (2002) for different purposes. It ensures that the new proposal ${\mathit{w}}_i$ is still an OA $(n,2^1,1)$. The procedure would yield an OA $(n,2^i,2)$ if the D-efficiency achieves the theoretical optimum in line 14. This D-efficiency upper bound is justified in the supplementary materials. To evaluate the relative D-efficiency of the current and new proposals, a computationally efficient method, which is also justified in the supplementary materials, is applied in line 17. The current proposal ${\mathit{x}}_i$ would be replaced with the new proposal ${\mathit{w}}_i$ if a more D-efficient component level matrix is generated by ${\mathit{w}}_i$. To prevent the search procedure from being trapped in a local optimum, $n(n-1)/2$ new proposals will be generated to improve the current proposal. Specifically, the pre-specified upper limit $n(n-1)/2$ is determined based on my experience in implementing Algorithm 2. After considering $n(n-1)/2$ new proposals, as noted in line 24, if the theoretical optimum is not yet achieved, another random initial will be considered. Given an OA $(n,2^{i-1},2)$, Algorithm 2 may or may not succeed to get an OA $(n,2^i,2)$ by incorporating an OA $(n,2^1,1)$. The search procedure will be terminated, as noted in line 5, when an OA $(n,2^u,2)$ is obtained or v random initials have been proposed for adding an OA $(n,2^1,1)$ to the current component level matrix. As a greedy local search procedure, Algorithm 2 is not guaranteed to obtain an OA $(n,2^u,2)$ for every single implementation. It can only guarantee to return an OA $(n,2^i,2)$ with $i\le u$ within an acceptable computation time. Practitioners must implement Algorithm 2 several times with a sufficiently large number of random initials v to get an OA $(n,2^u,2)$ satisfying conditions (B) and (C) of Definition 1. If the D-efficiency upper bounds in lines 14 and 17 are replaced with $b\times n^i$, where $0<b<1$, Algorithm 2 would return a near OA $(n,2^u,2)$ much faster.

Algorithm 2

Determine orthogonal component level matrices of DOAs

1: input m, n, u, v and Z

2: Generate V according to (7) of Peng, Mukerjee, and Lin (2019)

3: $\mathit{H}\leftarrow \mathit{I}-\frac{3(m-1)}{n(m+1)}\mathit{Z}{\mathit{Z}}^{\top }+\frac{3}{n(m+1)}\mathit{ZV}{\mathit{Z}}^{\top }$

4: initialize i, j, k and ${\mathit{X}}_i$

5: while i < u and j < v do

6: $i\leftarrow i+1$

7: Generate an OA $(n,2^1,1)$ randomly, denoted by ${\mathit{x}}_i$

8: ${\mathit{X}}_i\leftarrow \left[{\mathit{X}}_{i-1},{\mathit{x}}_i\right]$

9: $j\leftarrow j+1$

10: if $k<n(n-1)/2$ then

11: Generate ${\mathit{w}}_i$ by randomly switching a pair of + 1 and –1 in ${\mathit{x}}_i$

12: ${\mathit{W}}_i\leftarrow \left[{\mathit{X}}_{i-1},{\mathit{w}}_i\right]$

13: $k\leftarrow k+1$

14: if $\det ({\mathit{W}}_i^{\top }\mathit{H}{\mathit{W}}_i)=n^i$ then

15: ${\mathit{X}}_i\leftarrow {\mathit{W}}_i$

16: initialize j and k

17: else if $\det ({\mathit{X}}_i^{\top }\mathit{H}{\mathit{X}}_i)<\det ({\mathit{W}}_i^{\top }\mathit{H}{\mathit{W}}_i)$ and $\det ({\mathit{W}}_i^{\top }\mathit{H}{\mathit{W}}_i)<n^i$ then

18: ${\mathit{x}}_i\leftarrow {\mathit{w}}_i$

19: initialize k

20: goto 10

21: else

22: goto 10

23: end if

24: else if $k\ge n(n-1)/2$ then

25: initialize k

26: go to 7

27: end if

28: end while

29: $\mathit{X}\leftarrow {\mathit{X}}_i$

30: output X

Example 3.

Suppose that all chemotherapeutics have two dose levels in the five-drug combination experiment in Example 2. Algorithm 1 needs an OofA-OA $(n,15,2)$ as a template to generate a DOA $(n,5,2^5,2)$. In this case, the run size n must be larger than or equal to 108. By combining optFederov and Algorithm 2, a hybrid procedure is used to generate a more economical design. First, an OofA-OA $(24,5,2)$, denoted by Z, is generated by executing the R function optFederov to satisfy condition (A) of Definition 1. Next, Algorithm 2 is implemented to search for an OA $(24,2^5,2)$, denoted by X, such that both conditions (B) and (C) of Definition 1 are satisfied. By definition, the matrix pair $(\mathit{Z},\mathit{X})$ is a DOA $(24,5,2^5,2)$. Its design matrix is presented in Table 4.

Table 4 Design matrix of a DOA $(24,5,2^5,2)$.

DOA design matrix
(1, +)	(2, +)	(4, –)	(3, +)	(5, –)
(1, +)	(2,–)	(5, –)	(3, –)	(4, +)
(1, –)	(3, –)	(4, +)	(5, –)	(2, +)
(1, –)	(5, –)	(4, –)	(3, +)	(2,–)
(2, +)	(1, +)	(3, +)	(5, +)	(4, –)
(2, +)	(1, –)	(4, –)	(5, –)	(3, –)
(2,–)	(3, –)	(4, –)	(1, –)	(5, +)
(2, +)	(4, +)	(5, +)	(3, +)	(1, –)
(2,–)	(5, +)	(1, +)	(4, +)	(3, –)
(3, +)	(1, –)	(4, –)	(2,–)	(5, +)
(3, +)	(1, +)	(5, +)	(2, +)	(4, +)
(3, +)	(2,–)	(5, –)	(4, –)	(1, –)
(3, +)	(4, +)	(2, +)	(5, –)	(1, +)
(3, –)	(5, +)	(2,–)	(1, –)	(4, +)
(4, +)	(1, –)	(3, –)	(2, +)	(5, +)
(4, +)	(1, +)	(5, +)	(3, +)	(2,–)
(4, +)	(2,–)	(3, –)	(1, +)	(5, –)
(4, –)	(3, –)	(5, –)	(1, +)	(2,–)
(4, –)	(5, +)	(2,–)	(1, +)	(3, +)
(5, –)	(1, –)	(2,–)	(3, +)	(4, +)
(5, +)	(3, –)	(1, +)	(4, –)	(2, +)
(5, –)	(3, –)	(2, +)	(4, –)	(1, +)
(5, +)	(4, –)	(1, –)	(2, +)	(3, –)
(5, –)	(4, +)	(2, +)	(3, +)	(1, –)

As shown in Example 3, dual-orthogonal arrays of strength two can be generated using a hybrid procedure. To satisfy condition (A) of Definition 1, the R function optFederov is first executed to get an OofA-OA $(n,m,2)$, denoted by Z. Because optFederov cannot guarantee to achieve the global optimum, the moment matrix must be calculated to confirm whether or not the pairwise order matrix Z is really an OofA-OA $(n,m,2)$. Additionally, optFederov seems unable to generate an OofA-OA $(n,m,2)$ for m > 7. The order-of-addition orthogonal arrays of strength two provided by Tsai (2022) are used for the cases of m = 8 and 9. Given an OofA-OA $(n,m,2)$, denoted by Z, as in condition (A) of Definition 1, Algorithm 2 is implemented to get an OA $(n,2^u,2)$, denoted by X, such that both conditions (B) and (C) of Definition 1 are satisfied. After pairing Z and X, the matrix pair $(\mathit{Z},\mathit{X})$ is a DOA $(n,m,2^u,2)$. A series of dual-orthogonal arrays of strength two is generated using this hybrid procedure. The design parameters are summarized in Table 5 and the corresponding design matrices are collected in the supplementary materials. In Table 5, for every m, the number of dual-orthogonal arrays generated by the hybrid procedure is denoted by $n_{\text{doa}}$. Except for the case of m = 9, several designs are provided for each m, where $n_{\min }$ and $n_{\max }$ denote the run sizes of the smallest and largest designs in the catalogue. Researchers can choose order-of-addition two-level factorial designs with the most suitable run sizes within $n_{\min }$ and $n_{\max }$ for their experiments. All designs in Table 5 allow us to consider u = m two-level components. To design experiments with u < m, practitioners just need to skip those unused component levels. Some dual-orthogonal arrays of strength two with $n<n_{\min }$ and u < m are also generated using the hybrid procedure. Their design matrices are collected in the catalogue. Additionally, the hybrid procedure can be used to generate dual-orthogonal arrays of strength two with $n>n_{\max }$. However, such searches will be computationally expensive. More discussions regarding $n_{\min }$ and $n_{\max }$ can be found in the supplementary materials.

Table 5 Parameters of order-of-addition two-level factorial designs in the catalogue.

m	u	$n_{\text{doa}}$	$n_{\min }$	$n_{\max }$
3	3	3	12	36
4	4	3	24	48
5	5	9	24	120
6	6	57	36	720
7	7	209	48	5040
8	8	3	96	192
9	9	1	96	96

7 Concluding Remarks

Dual-orthogonal arrays are proposed in this paper to design order-of-addition two-level factorial experiments. An order-of-addition two-level factorial design is $\varphi$-optimal for estimating the compound model in (1) if the corresponding matrix pair $(\mathit{Z},\mathit{X})$ is a dual-orthogonal array of strength two. Nonetheless, the compound model in (1) only takes main effects into account. As shown in Section 5.3, adding interaction effects to the fitted model can be practical for exploring optimal treatments. Under the framework of order-of-addition factorial experiments, three types of interaction effects can be studied: (I) interaction effects between pairwise order factors, (II) interaction effects between component level indicator variables and (III) interaction effects between pairwise order factors and component level indicator variables. The treatment-response relationships of interest would be better understood if these three types of interaction effects are identified properly. Mee (2020) established a theoretical framework to study all Type I interaction effects. When the expanded model with all Type II interaction effects is used, relative D-efficiencies of some small designs in the catalogue are reported in the supplementary materials. Although these designs appear to be highly D-efficient, they do not seem to be optimal. In some real-world studies, Type III interaction effects may also play a key role in identifying optimal treatments. Although pairwise order effects and component main effects can be estimated separately from an order-of-addition experiment and a two-level factorial experiment, Type III interaction effects are not estimable using all observed data, because component levels are fixed while conducting the order-of-addition experiment and the component addition orders are fixed while conducting the two-level factorial experiment. If an order-of-addition two-level factorial design determined by a dual-orthogonal array of strength two is used, then some Type III interaction effects would be estimable. Because both component addition orders and component levels are varied over treatments, experimental designs determined by dual-orthogonal arrays of strength two enable us to explore more interaction effects after identifying active main effects. However, the relative D-efficiency of a dual-orthogonal array of strength two for estimating an expanded model with all Type I or Type III interaction effects is unavailable because currently there is no theoretical D-efficiency upper bound.

Based on subject matter knowledge and/or prior information, when some potentially active interaction effects are specified by researchers, a two-stage procedure can be used to obtain a highly D-efficient design for estimating the requirement set consisting of all main effects and the user-specified interaction effects. First, the hybrid procedure in Section 6.2 is run multiple times to get several dual-orthogonal arrays of strength two. Next, the most D-efficient nonsingular design for estimating the requirement set is chosen for the experiment. The elitist design is not only $\varphi$-optimal for estimating the compound model in (1) but also highly D-efficient for estimating the requirement set. Given a requirement set, a highly D-efficient design can also be generated by replacing the determinants in lines 14 and 17 of Algorithm 2 with the determinants of the expanded moment matrices. Because Algorithm 2 updates the model matrix corresponding to the main effects in each iteration, the model matrix corresponding to the interaction effects can also be updated accordingly. The determinant of the expanded moment matrix can then be calculated to evaluate whether or not a new proposal should be accepted. By repeating this operation several times, the D-efficiency for estimating the requirement set can be improved. Nonetheless, a suitable stopping criterion cannot be determined because a theoretical D-efficiency upper bound does not currently exist. On the other hand, dual-orthogonal arrays of strength three or more also seem to be promising tools to construct optimal designs for estimating different types of interaction effects. Because order-of-addition orthogonal arrays and two-level orthogonal arrays can have strength three or more, both conditions (A) and (B) of Definition 1 can be extended easily. However, a statistically well-justified extension for condition (C) of Definition 1 is unavailable. Additionally, when designing order-of-addition factorial experiments, some components may have three or more levels. These multi-level components enable us to explore quadratic, cubic and higher-order component main effects. Extending dual-orthogonal arrays to three or more levels can also be practical. These interesting topics may be worth pursuing in future research.

Supplementary Materials

Supplementary materials of this article include the following sections:

(S1)Fractional three-drug combination experiments

(S2)Proof of Theorem 1

(S3)Optimal treatment determination

(S4)Fast evaluation of relative D-efficiency

(S5) $n_{\min }$ and $n_{\max }$

(S6)Relative D-efficiency for an expanded model

(S7)Design catalogue

(S8)Code

Acknowledgments

I would like to thank the editor, associate editor and two anonymous referees for their valuable comments and constructive suggestions. I am also grateful to Computer and Information Networking Center of National Taiwan University for the support of high-performance computing facilities.

Additional information

Funding

Shin-Fu Tsai’s research was supported in part by National Science and Technology Council of Taiwan (Grant Numbers NSTC 109-2118-M-002-004-MY2 and 111-2118-M-002-005-MY2).