Optimal allocation of sample size for randomization-based inference from 2^K factorial designs

1 Introduction

Randomized 2 K factorial experiments are conducted to assess the marginal causal effects of K factors, each with two levels, along with their interactions on a response of interest. The two levels are often denoted as the “high level” and “low level” of the factor [1,2]. With K factors, there are 2 K unique treatment combinations to which units can be assigned. In the twentieth century, factorial designs have mostly been discussed in an industrial setting, whereas in recent times, there has been a lot of interest in their application in the social, behavioral and biomedical sciences and randomization-based inference from such designs (e.g., [3,4]).

Randomization-based inference is a useful methodology for drawing inference on causal effects of treatments in a finite population setting (e.g., [5,6]). A major advantage of randomization-based inference is that it applies even if the experimental units are not randomly sampled from a larger population, which is the case in most social science experiments [7,8]. The theory, methods, and applications of randomization-based inference for two-level factorial experiments with a completely randomized treatment assignment mechanism have been developed and discussed (e.g., [9,10]). Further, randomization-based inference from experiments with more general factorial structures and complex assignment mechanisms have been discussed in the study by Mukerjee et al. [11]. Connections between regression- and randomization-based causal inference from factorial experiments have been studied by Zhao and Ding [12].

Despite the growing literature in this area, most of the recent research on randomization-based inference of factorial experiments has been confined to the analysis side. On the design side, the main focus has been on rerandomization [3,13,14], which generalizes the idea of blocking by predefining an acceptable criterion for randomization based on covariate balance between treatment groups. There have also been extensions to fractional and incomplete factorial designs [15] and to the use of screening steps [16]. However, the distribution of the total number of experimental units into the 2 K treatment groups has not received much attention. Balanced designs that assign an equal number of units to the treatment groups are often the default choice, but it is unclear whether they are the “best” design under different conditions. One work that does discuss how to allocate units to optimize precision of factorial estimators from the randomization-based perspective is the one by Blackwell et al. [17]. That work explores the advantages of Neyman allocation [18,19] by extending the two-stage adaptive design proposed in the study by Hahn et al. [20] to multiple treatment designs. Dai et al. [21] similarly explores Neyman allocation but within sequential designs.

To motivate this problem, we consider an education experiment from [22], conducted to assess the causal effects of two different interventions, a student support program (SSP) and a student fellowship program (SFP), on the academic performance of freshmen. This is a 2 2 factorial experiment in which each unit (freshman) can receive only one of four treatment combinations: control (neither of the two), SSP only, SFP only, and student fellowship and support program (SFSP) that indicated both interventions applied to a unit. The units were divided into two blocks based on their sex. Table 1 shows the allocation of units within each block to the four treatment combinations.

Clearly, the design is unbalanced, with the highest number of units assigned to control and the fewest to the treatment SFSP. Such an allocation, among other reasons, could be motivated by the budget for experimentation. The question we investigate in this article is the following: Under what assumptions, conditions, and requirements will such an allocation be the best possible one (in terms of being able to precisely answer scientific questions of interest)?

The problem of finding optimal designs in the context of model-based inference has been extensively studied in the twentieth century (see, e.g., [23]). In such settings, optimal designs depend on a postulated outcome model that may be linear or nonlinear. For example, for binary responses, optimal designs based on logistic models are likely to differ from those based on probit models and depend on unknown model parameters [24–26]. We aim to develop optimal designs that are tied to model-free, randomization-based analysis for finite and super populations. In addition to being robust to model assumptions, our approach works for continuous as well as binary outcomes as long as the finite-population estimand is well-defined.

This article is organized as follows: The next section introduces basic notation and estimands for factorial experiments using the potential outcomes framework. In Section 3, we derive optimal allocations of the N units in a population to different treatment combinations under three commonly used optimality criteria for a completely randomized design (CRD). In addition to theoretical results for exact optimal allocations, we also provide numerical algorithms for obtaining integer solutions. In Section 4, we extend our results for CRDs to the setting of randomized block designs (RBDs). In Section 5, we derive optimality results under cost constraints. Two different applications of the proposed methodology, the motivating education experiment and an audit experiment conducted to assess discrimination, are described in Section 6. We conclude with a discussion, including opportunities for future work, in Section 7.

2 2 K factorial experiments under the potential outcomes framework

Here, we introduce some key definitions and notations from the study by Dasgupta et al. [9]. Consider a 2 K experiment with N units, in which the levels of each of the K factors are denoted by 0 and 1. Each treatment combination is of the form z j = ( z 1 … , z K ) , where z k ∈ { 0 , 1 } for k ∈ 1 , … , K . There are J = 2 K treatment combinations arranged in lexicographic order 1 , … , J , where treatment combination z j is such that j = 2 K − 1 z 1 + 2 K − 2 z 2 + … + z K + 1 . In other words, ( z 1 … z K ) is a binary representation of integer j − 1 . Thus, for example, in a 2 2 experiment, the four treatment combinations 00, 01, 10, and 11 are numbered as j = 1 , 2 , 3 , and 4, respectively, and the 8 th treatment combination in a 2 4 experiment is 0111. We will just refer to the treatment combination by its number ( j ) in notation below.

For i = 1 , … , N , under the stable unit treatment value assumption (SUTVA) [27], the i th unit has J = 2 K potential outcomes Y i ( 1 ) , … , Y i ( J ) corresponding to the J treatment combinations z 1 , … , z J . Let Y i denote the J × 1 vector of potential outcomes for unit i . For unit i , the unit-level main effect of factor k = 1 , … , K is defined as the difference between the averages of potential outcomes for unit i for which the levels of factor k are at levels 0 vs 1. Mathematically, it is a contrast of the form 2 − ( K − 1 ) λ k T Y i = 2 − ( K − 1 ) ∑ j = 1 J λ j k Y i ( j ) , where x T denotes the transpose of vector x , λ k is a J × 1 column-vector with coefficient λ j k such that λ j k = − 1 if the level of factor k in j th treatment combination is 0, and λ j k = 1 otherwise. For all k = 1 , … , K , λ k is a contrast vector, i.e., ∑ j = 1 J λ j k = 0 .

Proceeding along the lines of the study by Dasgupta et al. [9], for unit i , we can define K 2 two-factor interactions, K 3 three-factor interactions, and finally one K -factor interaction as contrasts of the form 2 − ( K − 1 ) λ T Y i , where the contrast vector λ for any interaction can be derived by element-wise multiplication of the contrast vectors of the corresponding main effects λ k , for factors involved in the interaction. Denoting the J − 1 = 2 K − 1 contrast vectors for the J − 1 unit-level factorial effects τ 1 i , … , τ J − 1 , i by λ 1 , … , λ J − 1 , respectively, we define the J × J matrix as follows:

where λ 0 is the J × 1 vector with all elements equal to one. We note that L is an orthogonal matrix with L L T = L T L = 2 K − 1 I J , where I J denotes the identity matrix of order J . For i = 1 , … , N , let τ i = ( 2 τ 0 i , τ 1 i , … , τ J − 1 , i ) T , where τ 0 i denotes the average of all potential outcomes for unit i . The linear transform between the vector of unit-level potential outcomes Y i and the vector of unit-level factorial effects τ i can be expressed as follows:

Having defined unit-level factorial effects, we now move to their population-level counterparts. Let Y ¯ = N − 1 ∑ i = 1 N Y i and τ = N − 1 ∑ i = 1 N τ i , respectively, denote the J × 1 vectors of average potential outcomes and the average factorial effects. Then, averaging equation (2) over i = 1 , … , N , the vector of population-level factorial effects is given as follows:

Note that the first element of τ is twice the average of all potential outcomes.

In a CRD, a pre-assigned number of units, N j , are randomly assigned to treatment j . The experiment generates an N × 1 vector of observed outcomes data from which the vector of factorial effects τ can be unbiasedly estimated. We examine the properties of an unbiased estimator of τ defined in the next section, with respect to its randomization distribution, and formulate the problem of optimally allocating the N units to the J treatment combinations in Section 3.

3 Optimal designs for completely randomized experiments

In a randomized experiment with N j units assigned to treatment combination j ∈ { 1 , … , J } , only one of the J potential outcomes is observed for unit i . This observed outcome is y i = Y i ( T i ) for i = 1 , … , N , where T i is the random treatment assignment variable for unit i taking value j if unit i receives treatment j . In a CRD, the joint probability distribution of ( T 1 , … , T N ) is

where 1 { A } denotes the indicator random variable for set A . Let y ¯ ( j ) = N j − 1 ∑ i = 1 N 1 { t i = j } Y i ( j ) denote the average response for treatment j . Let y ¯ denote the vector ( y ¯ ( 1 ) , y ¯ ( 2 ) , … , y ¯ ( J ) ) T of observed averages. Substituting y ¯ in place of Y ¯ in equation (3), we can unbiasedly estimate the vector of factorial effects as follows:

Lu [10] derived sampling properties of the estimator τ ^ with respect to its randomization distribution for the general case of unequal N 1 , … , N J . Lu [10] showed that τ ^ is an unbiased estimator of τ , and has the following finite-population covariance matrix:

where λ ˜ j represents the transpose of row j of the model matrix L defined in equation (1), τ i denotes the vector of unit-level factorial effects given by equation (2), τ denotes the vector of population-level factorial effects given by equation (3), and

denotes the variance of all N potential outcomes for treatment j with divisor N − 1 , where Y ¯ ( j ) = N − 1 ∑ i = 1 N Y i ( j ) .

In the spirit of classical optimal designs [23], we can define a design optimality criterion as a functional of the matrix V τ defined in equation (5). For example, the D-optimality criterion, which aims to minimize the determinant of the covariance matrix, or the A-optimality criterion, which aims to minimize the trace of the covariance matrix, or the E-optimality criterion, which aims to minimize the maximum eigenvalue of the covariance matrix, can be considered. However, the second term 1 ∕ ( N ( N − 1 ) ) ∑ i = 1 N ( τ i − τ ) ( τ i − τ ) T , which is a measure of heterogeneity of treatment effects, cannot be estimated from observed data, because none of the unit-level treatment effects τ i are estimable without additional assumptions due to the missing potential outcomes. Because 1 ∕ ( N ( N − 1 ) ) ∑ i = 1 N ( τ i − τ ) ( τ i − τ ) T is positive semi-definite, the first term of equation (5) can be considered as an upper bound of V τ , which is attained under specific restrictions on the potential outcomes (e.g., treatment effect homogeneity). Thus, we propose optimizing a functional of the first term of equation (5), which in turn is equivalent to optimizing a functional of the positive definite matrix

instead, where A = diag ( S 1 2 ∕ N 1 , … , S J 2 ∕ N J ) .

Another justification for using a functional of the matrix V ˜ as an optimality criterion comes from the super-population perspective mentioned in Section 2. Ding et al. [29] showed that the estimator τ ^ defined earlier is also an unbiased estimator of the super-population estimand τ SP . Furthermore, extending their argument for a single factor with two levels to the case of 2 K factorial designs, if V j 2 = Var [ Y i ( j ) ] , j = 1 , … , J , then a variance decomposition yields the following sampling variance of τ ^ ,

where Var SP denotes variance over the random sampling from the super population and random assignment, and V j 2 = E ( S j 2 ) represents expectation of S j 2 with respect to the distribution of the potential outcomes in the super population. This connection provides further motivation for the form of our optimization, but we focus on the finite-population setting going forward.

3.2 Computation of integer optimal designs using an integer programming approach

Many sources in the integer programming literature address constrained optimization under integer space constraints (see, e.g., [31–34]). We adopt the methods proposed in the study by Friedrich et al. [35] that are designed for settings very similar to the ones we consider. Friedrich et al. [35] considered the following integer programming problem:

(7) min N 1 , … , N J f ( N 1 , … , N J ) s.t. ∑ j = 1 J N j = N , l j ≤ N j ≤ u j ∀ j = 1 , 2 , … , J , N j ∈ Z + ∀ j = 1 , 2 , … , J ,

where Z + is the set of positive integers, ( l j , u j ) are the lower and upper bound constraints on N j , and f : R + J → R is a convex function. If f ( N 1 , … , N J ) is separable (i.e., can be expressed as ∑ j = 1 J f j ( N j ) ) then the greedy algorithm in Figure 1 finds the globally optimal integer solution of the minimization problem given in equation (7) under some regularity conditions that we show in Supplementary Material S2.

Figure 1

Greedy algorithm for separable functions.

From the proof of Theorem 1 in Supplementary Material S1, it follows that the A- and D-optimality criteria can, respectively, be expressed as ∑ j = 1 J ( S j 2 ∕ N j ) and ∑ j = 1 J ( log S j 2 ∕ N j ) . In Supplementary Material S2, we show that the conditions for global convergence of the greedy algorithm are met, and thus, convergence to the true integer optimal solution in the cases of A- and D-optimality are guaranteed if this algorithm is used. Hence, substitution of S j 2 ∕ N j and log S j 2 ∕ N j for f j ( N j ) in the algorithm described in Figure 1 leads to optimal integer solutions for the A-optimality criterion and the D-optimality criterion, respectively. In our implementation of this algorithm, we take l j = 2 and u j = N − 2 J for j = 1 , … , J to guarantee at least two units are assigned to each treatment combination, allowing variance estimation within each treatment group.

We provide a modified greedy algorithm described in Figure 2 for solving the E-optimality problem since the optimization problem cannot be written in a separable form as in the aforementioned A- or D-optimality criterion. In this algorithm, we take f j ( N j ) = S j 2 ∕ N j while keeping l j = 2 and u j = N − 2 J for j = 1 , … , J . In Section 4.2, the ability of this greedy algorithm to find E-optimal solutions is demonstrated empirically for blocked designs, discussed next.

Figure 2

Greedy algorithm for E-optimality.

4 Optimal allocation for factorial experiments with blocks

Consider a block-randomized 2 K factorial design with H blocks. That is, units are pre-assigned membership to one of h blocks based on some similarity metric (we do not consider how to form blocks here). Let M h denote the size of block h , h = 1 , … , H . In addition, let M h , j be the number of units in block h assigned to the treatment j ( j = 1 , … , J ). Finally, let b i ( h ) , i = 1 , … , N , be an indicator variable, taking the value 1 if unit i belongs to block h and 0 otherwise. Treatment assignment under a factorial RBD is equivalent to performing an independent CRD, as described in Section 3, within each block.

The population average treatment effect τ can be expressed as ∑ h M h τ h ∕ N , where τ h is the block-level vector of factorial effects and its estimator τ ^ is a weighted average of τ ^ h , an unbiased estimator of τ h defined in the same way as in equation (4) for block h . Extending equation (5) by noting the independence of the assignment to treatment across blocks, the covariance matrix of τ ^ h in a factorial RBD can thus be obtained as follows:

where λ ˜ j represents the transpose of row j of the model matrix L defined in equation (1) and S h , j 2 denotes the variance of all M h potential outcomes for units in block h under treatment j with divisor M h − 1 . The covariance matrix of τ ^ can then be expressed as Cov ∑ h M h τ ^ h ∕ N . Because the block-level treatment estimators τ ^ h are independent across blocks, we have

Writing S blk , j 2 = ∑ h = 1 H ( M h 2 ∕ N 2 ) ( S h , j 2 ∕ M h , j ) and proceeding along similar lines as in Section 3, we formulate a surrogate optimization problem to only optimize the first term, since the second term in the equation above is not identifiable. Then, choosing M h , j s to optimize some functional of Cov ( τ ^ ) is equivalent to optimizing a functional of the matrix

where A blk = diag ( S blk , 1 2 , … , S blk , J 2 ) .

5 Optimal allocation driven by cost constraints

So far, we have considered optimality criteria that are based on the covariance matrix of the estimated factorial effects, implicitly assuming that all treatment combinations are equally expensive (with respect to cost and/or time). However, such assumptions may not be true in many practical situations, and cost constraints can play an important role in determining optimal allocation. Thus, it is worthwhile to explore solutions to optimal allocation under cost constraints.

We consider the optimal allocation for 2 K factorial CRDs. Let the cost of assigning treatment combination j to one unit be C j > 0 , and the total available budget be C . In the new optimization problem, we replace the constraint ∑ j N j = N in the original problem described in Section 3.1 by the cost constraint ∑ j C j N j ≤ C . The new optimization problem is, therefore:

where V ˜ = ∑ j = 1 J S j 2 N j λ ˜ j λ j ˜ T and ψ ( V ˜ ) is a functional of V ˜ . A straightforward approach to incorporate this new constraint into our previous setting is to re-write the constraint as ∑ j N ˜ j ≤ C , where N ˜ j = N j C j is the total cost for the suggested allocation to treatment arm j . Note that the total sample size N is no longer fixed before the optimization. Under this one-to-one transformation N ˜ j = C j N j , the optimization problem in equation (10) is equivalent to minimizing the objective function over N ˜ j [36] and can be written as follows:

Because the optimal solution of the above optimization problem is attained at ∑ j N ˜ j = C , the inequality constraint can be replaced by the equality constraint. Then, proceeding along the lines of Theorem 1, one can obtain the cost for the optimal allocation to treatment arm j as N ˜ j ∝ S j ( C j ) , N ˜ j = C ∕ J and N ˜ j ∝ S j 2 C j as the A-, D-, and E-optimal solutions. These results are formalized in terms of the optimal proportion of the budget allocated to each treatment arm, which can be used to determine the number of units to assign to each treatment arm, in Theorem 3.

We use an example to demonstrate the applicability of Theorem 3. Consider a 2 2 factorial design, with C = 100 , per unit cost vector ( C 1 , … , C 4 ) = ( 0.1 , 4 , 4 , 9 ) . This setup represents many common scenarios where treatment arm 1 represents the control group 00 and involves a per-unit cost that is negligible compared to the ones with at least one active treatment. On the other hand, treatment arm 4 has both treatments at active level and involves the highest cost. Table 3 shows the A-, D-, and E- optimal proportions of total cost π j s for two different vectors of variances ( S 1 2 , … , S 4 2 ) . In one setting, we take the vector as ( 1 , 1 , 1 , 1 ) and in another, set it to ( 1 , 2 , 3 , 4 ) . For the sake of completeness, we also add a column of equal cost ( 1 , 1 , 1 , 1 ) , under which the p j ’s of Theorem 1 and π j ’s of Theorem 3 become identical, as explained in Remark 8. Thus, the optimal allocations in the first column of Table 3 can also be derived from Theorem 1 with N = C = 100 .

One can obtain the A-, D-, and E-optimal allocations of N j s by substituting the optimal π j s from Theorem 3 into N j = ( C π j ) ∕ C j . However, rounding these optimal N j s into nearest integers may lead to violation of the constraint ∑ j C j N j ≤ C . To avoid such possibilities, one can consider the optimal values of ⌊ ( C π j ) ∕ C j ⌋ as approximate integer solutions, where ⌊ x ⌋ denotes the largest integer contained in x .

Researchers may decide to impose an additional constraint on the optimization problem (5) that forces the sum of N j s to be exactly equal to a predetermined N . Such a problem would give optimal allocation under fixed N , unlike Theorem 3. However, imposing this additional constraint may force the set of feasible solutions to the optimization problem to be empty. For example, suppose for all j , C j > C ∕ N . Then, ∑ j C j N j > ∑ j ( C ∕ N ) N j , which exceeds the allowable cost C if the restriction ∑ j N j = N is imposed. Thus, additional conditions are necessary to guarantee that the feasible set is nonempty. Obtaining closed-form solutions under such conditions may not be straightforward, and one may need to rely on numerical methods to obtain such solutions.

6 Applications in real experiments

In this section, we demonstrate applications of the results and algorithms developed to two real-life experiments. First, we revisit the education example from [22] described in Section 1. Second, we discuss a pilot audit experiment reported in Libgober’s work [37] conducted to identify how perceptions of race, gender, and affluence affect access to lawyers and demonstrate how the proposed methodology can be used to design follow-up experiments in similar populations.