Revisiting Randomization with the Cube Method

Revisiting Randomization with the Cube Method

\smallbreak

{\let\number\value{footnote}\relax\footnote{{

We would like to thank Riccardo d'Adamo, Isaiah Andrews, Yuehao Bai, Stéphane Bonhomme, Guillaume Chauvet, Russell Davidson, Max Farrell, Xavier d'Haultfoeuille, Marc Henry, Keisuke Hirano, Xinran Li, John List, Kirill Ponomarev, Pauline Rossi, Jean Rubin, Felix Schleef, Azeem Shaikh, Sami Stouli, Max Tabord-Meehan, Yves Tillé, Panos Toulis and many other colleagues in workshops and conferences whose comments and suggestions helped improve this paper. In particular, we would like to thank everybody at the 2023 Bristol Econometric Study Group, the 2023 Advances with Field Experiments Conference, the 18th IZA & 5th IZA/CREST Conference, the 2023 European Winter Meeting of the Econometric Society, the AMSE Big Data and Econometrics Seminar, the 2024 RCEA International Conference in Economics, Econometrics, and Finance, the 2024 BSE Summer Forum Workshop on Microeconometrics and Policy Evaluation, the 2024 AFSE Annual Conference, the 2024 Annual Conference of the International Association for Applied Econometrics, the Econometrics Workshop at UChicago, the PhD Workshop at Penn State.}}}

\setcounter{footnote}{0} \justifying \@startsection{section}{2}{0mm}{-1.5\baselineskip}{1\baselineskip}{\normalfont}{Introduction}

Prior to running RCTs, researchers often have access to rich information on experimental units through baseline surveys or administrative data. Researchers would, ideally, like to create treatment and control groups that balance available covariates. How to randomize using available covariates is still a matter of debate, especially when the number of covariates is large.

We here introduce a randomization method based on the cube algorithm proposed by deville_efficient_2004 for survey sampling. Units are iteratively assigned to treatment or control groups. At each step, one unit is allocated, under the constraint that means of covariates in the two groups are equal. Assignment probabilities of unallocated units are updated at each step to ensure that empirical means of available covariates are the same in control and treatment.

The cube method aims at balancing chosen moments, such as mean, variance, or cross-moments. Still, it does not impose balancing constraints on moments considered of little relevance (for instance, high-order moments). In contrast, stratification and matched-pairs designs, two prominent methods proposed by bai_inference_2022,bai_optimality_2022 and cytrynbaum_optimal_2023, aim at reducing the distance between groups' joint probability distributions of covariates to achieve balance of all moments. Crucially, we show that reducing the distance between joint probability distributions becomes infeasible in finite samples for more than a few covariates. Stratification may even backfire, i.e., increasing imbalances, when stratifying naively on “too many” covariates. The cube method greatly mitigates this curse of dimensionality. Figure (ref) illustrates how imbalances when randomizing with the cube method are much reduced compared to other methods as the number of covariates used for balancing $p$ grows, for a fixed sample size $n$.

The set of balanced covariates has a direct impact on the precision of estimated treatment effects. Maximal precision gains are achieved by selecting the set of covariates that are the most correlated with potential outcomes. But relevant covariates are only known after the RCT was run and must be chosen under some uncertainty. Furthermore, RCTs often estimate several treatment effects and corresponding outcomes could be correlated with different covariates, notably pretreatment outcomes. Balancing on a larger set of covariates increases the likelihood of selecting the most relevant ones. The cube method allows increasing the number of covariates, while being robust to the addition of irrelevant ones. Figure (ref) displays the precision of the average treatment effect estimator in a situation in which only one covariate, out of 30, correlates with potential outcomes. The cube method allows balancing on all 30 covariates, increasing precision, while not being penalized by adding 29 irrelevant covariates. The cube method stands out as other methods suffer from precision loss when balancing on irrelevant covariates.

This paper contributes to three streams of research. First, we extend the cube method, developed first for survey sampling by deville_efficient_2004. The cube method is routinely used for sampling by national statistical institutes tille_ten_2011. Our technical contribution is to extend the scope of the cube method beyond sampling for estimating treatment effects in RCTs. Moving from sampling to RCTs implies redefining balancing constraints and deriving asymptotic properties for estimators of the population average treatment effect (PATE) and the sample average treatment effect (SATE). For the PATE, we show that the semiparametric efficiency bound in hahn_role_1998 is attained for large $n$ and fixed $p$ under a linearity assumption of the conditional expectation of potential outcomes. We also provide valid inference strategies for the estimators. As is the case for the other methods achieving the bound, precision improves when the share of the variance of potential outcomes explained by the covariates is larger. We thus formally motivate the interest of using a method that balances almost exactly a large number of covariates $p$ for a given sample size $n$.

Our second contribution is to compare the balancing performance of existing randomization methods when $p$, the number of covariates used for balancing, increases. Asymptotic properties of randomization methods as the number of units $n$ gets large are well studied. In contrast, their asymptotic behavior when $p$ grows has been seldom studied yet, to the best of our knowledge. We show that observed patterns in Figure (ref) are de facto generic for bounded covariates. A key finding is the identification of three distinct regimes where stratification exhibits different behaviors. For a small number of covariates (i.e., $p\ll \ln(n)$), stratification improves balancing compared to complete randomization. When $p\approx \ln(n)$, there is a critical regime where the balancing quality deteriorates quickly. Last, for $p \gg \ln(n)$, because of the small strata issue, stratification is similar to a coin toss and worse than complete randomization. In sum, stratification exhibits different balancing properties when $p$ varies. We also derive upper and lower bounds for imbalances in matched-pair designs, showing that contrary to stratification, it always performs better than complete randomization but that imbalances increase rather quickly for $p\gg \ln(n)$. In sharp contrast, the cube method grants the balance of (selected moments of) covariates with no critical change when $p\ll n$, allowing balancing on a larger set of covariates.

Last, we contribute to the the literature that establishes pros and cons of randomization methods to guide experimenters when choosing a randomization method bruhn_pursuit_2009, athey_chapter_2017,bai_primer_2024. We briefly review the recent use of covariate-balancing designs in RCTs, and discuss practical implications of using the cube method for the publication process. Indeed, as explained above, the cube method removes most of the bad luck that may arise from sampling errors, as the most unfavorable samples have a null probability of being selected. Avoiding large imbalances has implications for the publication process. snyder_examining_2024 show that balance checks generates p-hacking and/or publication bias.

The remainder of this paper is structured as follows. Section (ref) introduces the potential outcome framework and covariate balancing. Section (ref) presents the cube algorithm and its application to RCTs. Section (ref) gives the balancing properties of the cube method, compares imbalances to other methods, and provides novel asymptotic expressions for the variance of average treatment effect estimators. We then specify two ways of performing inference based on asymptotic normality and the randomization mechanism. Section (ref) uses simulated and experimental data to show our precision gains and how the cube method might be less constrained by the curse of dimensionality. Finally, Section (ref) reviews current practices in RCTs and discusses practical considerations of the cube method.

\@startsection{section}{2}{0mm}{-1.5\baselineskip}{1\baselineskip}{\normalfont}{Setup} This section presents the potential outcomes framework, provides assumptions on the data-generating process, and formally defines covariate balancing.

\@startsection{subsection}{2}{0mm}{-1.2\baselineskip}{1\baselineskip}{\normalfont}{Data Generating Process and Assignment Design}

We consider the standard Neyman-Rubin framework of potential outcomes where $Y_i(0)$ is the outcome of unit $i$ when untreated and $Y_i(1)$ is the outcome when treated. We consider $X_i$ a vector of $p$ covariates. According to the literature, we assume i.i.d.ness and the existence of second-order moments for all these variables.

The researcher observes $(X_1,\ldots,X_n)$ for a finite sample of size $n$. She wants to randomly allocate these $n$ units to treatment according to a design $\Pi$, i.e., a distribution on the set of the possible treatment allocations $\{0,1\}^n$. If the design $\Pi$ does not depend on the potential outcomes, it balances potential outcomes in the treatment and control groups in average, avoiding selection bias. The design $\Pi$ could depend on $(X_1,...,X_n)$. For instance, the treatment probability of a unit $i$ could depend on $X_i$ for various reasons, such as efficiency, cost of the treatment depending on $X_i$, or subpopulations of particular interest. In the following, $D_i$ is the dummy variable indicating if $i$ is treated or untreated. Researchers have to choose not only each individual selection probability $\mathbb{P}_{\Pi}(D_i=1|X_1,...,X_n)$ but the full design $\Pi$ that determines $\mathbb{P}_{\Pi}\left(\cap_{i=1,...,n}D_i=d_i|X_1,...,X_n\right)$ for any potential allocations $(d_i)_{i=1,...,n}\in \{0,1\}^n$. A major issue is exploiting the knowledge of $(X_1,...,X_n)$ to define a “good” design $\Pi$ to go beyond the balancing of potential outcomes in average. To study this question, let us formulate the assumption on the class of design we consider in the following.

Assumption (ref) is usual and necessary in the literature on treatment effects estimation. Equation (ref) means that assignment is independent of the unknown potential outcomes, conditional on the auxiliary information $X$. Equation (ref) specifies that the assignment probability of unit $i$ could depend on $X_i$ but not on $X_j$ for $j\neq i$. It also states that the propensity score $p(X_i)$ fulfills a common support condition. \\ In what follows, we denote the propensity score $p(X_i)$ as $\pi_i$. Randomization methods are often presented under the assumption that $\pi_i=1/2$ for all $i$. However, the cube method handles heterogeneous assignment probabilities, which are of particular interest in RCTs, for at least three reasons. First, in view to minimize the variance of the estimator of the ATE, the optimal assignment probabilities corresponding to the so-called Neyman allocation are $\pi_i= V(Y (1)|X)^{1/2} \left(V(Y (1)|X)^{1/2} + V(Y (0)|X)^{1/2}\right)^{-1}$. Second, even if $V(Y (1)|X) = V(Y (0)|X)$ and if the researcher wants to minimize the variance of estimates of $E(Y (1) - Y (0))$, she could also adapt some assignment probabilities to heterogeneous costs $c_1(X)$, $c_0(X)$ of treatment and control to fulfill a budget constraint. Third, the researcher could be interested in treatment effects on some subpopulations $X = x$ that would not be precisely estimated if using a constant assignment probability. More generally, in an armed bandit perspective, researchers may adapt assignment probabilities with respect to what they learned to maximize some objective, explore treatment effects on some subpopulations, and minimize regret. Our proposition of design accommodates any propensity score type, offering complete flexibility to researchers concerning its definition.\\ Assumption (ref) ensures that for any variable $W$ such that $(D_1,...,D_n)\perp \!\!\! \perp (W_1,...,W_n)|X_1,...,X_n$ we have:

Equation (ref) is true for $W_i=X_i$ and $W_i=(Y_i(0),Y_i(1))$. \@startsection{subsection}{2}{0mm}{-1.2\baselineskip}{1\baselineskip}{\normalfont}{Parameters of interest and estimators} After the experiment, the researcher observes $Y_i=Y_i(1)\times D_i+Y_i(0)\times(1-D_i)$. She will thus never observe both potential outcomes for the same unit. researchers are generally interested in estimating the sample and population average treatment effects given by

and

respectively.\footnote{In some cases, they are interested in similar parameters for some subpopulations: $\frac{1}{\sum_{i=1}^n\mathds{1}\{X_i\in \mathcal{X}\}}\sum_{i=1}^n(Y_i(1)-Y_i(0))\mathds{1}\{X_i\in\mathcal{X}\}$ or $\mathbb{E}\left[Y_i(1)-Y_i(0)|X_i\in \mathcal{X}\right]$. Estimators of these quantities are defined by restricting the sample to units such that $X_i\in\mathcal{X}$ and the asymptotic properties of these estimators follow from a straightforward adaptation of what is presented below.}

In this paper, we will focus on the Horvitz-Thompson estimator (HT) and the Hájek estimator (H), which are of central interest in RCTs. The Horvitz-Thompson estimator is

which is unbiased under Assumption (ref) and (ref) for both the SATE and the PATE and is the difference between the inverse probability weighting estimators on the treated and the control group.

The Hájek estimator is

and corresponds as well to the inverse probability weighting OLS estimator $$\widehat{\theta}_H=\arg\min_{\theta}\min_a\sum_{i=1}^nw_i\left(Y_i-a-\theta D_i\right)^2$$ for $w_i=\frac{1}{\pi_i}$ if $D_i=1$ and $w_i=\frac{1}{1-\pi_i}$ if $D_i=0$. Let $n_T$ denote the number of treated units and $n_C$ the number of control units. When $\pi_i$ is constant, $\hat{\theta}_{H}=\frac{1}{n_T}\sum_{i: D_i=1}Y_i-\frac{1}{n_C}\sum_{i:D_i=0} Y_i$ is the difference between the average on the treated group and the control group whereas $\hat{\theta}_{HT}=\frac{1}{\mathbb{E}(n_T)}\sum_{i: D_i=1} Y_i-\frac{1}{\mathbb{E}(n_C)}\sum_{i: D_i=0}Y_i$ is a slight modification of this difference of averages. If, additionally, $n_T$ and $n_C$ are fixed, both estimators are identical to the difference-in-means estimator.

\@startsection{subsection}{2}{0mm}{-1.2\baselineskip}{1\baselineskip}{\normalfont}{Balancing Constraints}

Randomization methods generate control and treament groups that are balanced on average. A more stringent requirement is to generate groups that are exactly balanced:

Equation (ref) describes equality between the estimated weighted averages in the treatment and control groups. A perfectly balanced assignment eliminates any allocation to the treatment that does not balance perfectly the covariates between treatment and control groups.\\ A common practice in experiments is to form treatment and control groups of fixed sizes, $n_T$, and $n_C=n-n_T$, respectively. This is equivalent to satisfying the following constraint for any possibe allocation $(d_1,...,d_n)$:

In that case we also have: $n_C=\sum_{i=1}^n(1-d_i)=\sum_{i=1}^n(1-\pi_i)=\mathbb{E}(n_C)$. As recommended by deville_efficient_2004, we also balance a constant ($X_{ji}=1$), for the treatment and control groups:

Under such assignment $\widehat{\theta}_{HT}$ defined in ((ref)) is equal to $\widehat{\theta}_H$ defined in ((ref)). Notice that we can rewrite (ref), (ref), (ref) as

with $Z_{1i}= (1, \frac{\pi_i}{1-\pi_i},\pi_i,\frac{X_i'}{1-\pi_i})'$ and $Z_{0i}= ( \frac{1-\pi_i}{\pi_i},1,1-\pi_i,\frac{X_i'}{\pi_i})'=\frac{1-\pi_i}{\pi_i}Z_{1i}$. If assignment probabilities are homogeneous (i.e., $\pi_i=\pi$), the balancing covariates are reduced to $Z_{1i}=Z_{0i}=(1,X_i')$ due to perfect multicollinearity, but this is no more the case if the $\pi_i$ are heterogeneous.\\ It is worth noticing that exact balance is not always attainable: for instance, if $n=101$ and $\pi_i=1/2$. Imposing (ref) implies $n_T=50.5$, which is simply impossible. But statistical analysis ensures that balancing up to a $o_p\left(\frac{1}{\sqrt{n}}\right)$ is sufficient to take full advantage of the auxiliary information $X_i, \pi_i$. We here propose an "almost" exactly balancing design $\Pi$ such that for $(D_i)_{i=1,...,n}\sim \Pi$,

As we will show below, the cube method always achieves almost exact balancing. As a consequence precision gains are asymptotically equivalent from those achieved under exact balancing.

\@startsection{section}{2}{0mm}{-1.5\baselineskip}{1\baselineskip}{\normalfont}{The Cube Method} deville_efficient_2004 first introduced the cube method to produce samples balanced to the population. The cube method consists of an algorithm in two steps: the flight and landing phases. The technique gets its name from the graphical representation of a sampling problem. Equation (ref) ensures that balancing treatment and control groups in an experimental setting for some covariates is equivalent to balancing the treatment group to the entire sample. Let us consider the $n$-cube $C=[0,1]^n$. Each vertex of $C$ (from $2^n$ possibilities) represents a possible allocation: for instance, $(1,1,...,1)$ corresponds to the situation where all units are allocated to treatment, $(1,0,1,0,...,1,0)$ corresponds to the case where the treatment group is $\{i: i \text{ odd}\}$. A sampling design $\Pi$ corresponds to how a vertex is selected. Recall that we consider a framework where researchers impose that Equation (ref) holds for $\Pi$ and a vector $(\pi_i)_{i=1,...,n}$.

We will first describe the cube algorithm without balancing constraints before moving to the more interesting case where the balancing constraints in Equation (ref) are considered. Whatever the set of balancing constraints, the cube method is a discrete martingale that moves in (at most) $n$ steps from the interior point $\bm{\pi}(0)=(\pi_i)_{i=0}^n$ to $\bm{\pi}(n)=(D_i)_{i=0}^n$ a vertex of $C$. Let us consider the case without constraints. At the first step, one chooses a random direction for $\bm{\pi}(1)-\bm{\pi}(0)$ and a step size such that $\bm{\pi}(1)$ belongs to a facet of $C$ and that $\mathbb{E}[\boldsymbol{\pi}(1)|\boldsymbol{\pi}(0)]=\boldsymbol{\pi}(0)$. After this step, because $\bm{\pi}(1)$ belongs to a facet of $C$, one component $i_0$ of $\bm{\pi}(1)$ is equal to 0 or 1, selecting $D_{i_0}=\pi_{i_0}(1)$ one has thus assigned a first unit to either treatment or control. Because a facet of a $n$-cube is a $(n-1)$-cube, one can then repeat the process in a $(n-1)$-cube, and so on, until landing in a vertex of $C$. At the final step $n$, one will have $(D_i)_{i=1,...,n}=\bm{\pi}(n)\in\{0,1\}^n$ and $\mathbb{E}[D_i]=\pi_i$ (i.e., every unit is allocated to the treatment group with the probability specified by the researcher). These successive steps are the flight phase and for the cube method without balancing constraint, allocation $(D_i)_{i=1,...,n}$ is always determined at the end of this phase. Figure (ref) illustrates graphically the method. \\

In Figure (ref), all vertices of the $n$-cube can be selected, meaning that all individuals could be allocated to the control group. We now consider that the researcher wants to allocate a fixed number $n_T$ of units to the treatment and $n_c$ units to the control. This can be achieved with the cube method as soon as $\sum_{i=1}^n\pi_i=n_T$. The condition that exactly $n_T$ units are assigned to the treatment can be expressed as a balancing constraint. Indeed, because $n_T=\sum_iD_i$ and $\sum_i \pi_i=n_T$, the fixed size condition is equivalent to $\sum_{i}\frac{Z_iD_i}{\pi_i}=\sum_{i}Z_i$ for $Z_i=\pi_i$. Let $K$ the set of vectors $s$ in the $n$-cube $C$ such that $\sum_i s_i=n_T$. $K$ is a closed convex set and its extreme points are vertices of $C$, that is the set of allocations respecting the fixed-size constraints. $K$ is contained in an affine subspace of dimension $n-1$ of direction $V:=\{v: \sum_{i=1}^n v_i=0\}$, we have $K=C\cap \left\{\bm{\pi}(0)+v: \sum_i v_i=0\right\}$. The cube method selects randomly an element of $V$ for the direction of $\bm{\pi}(1)-\bm{\pi}(0)$ and fixes the step size such that $\bm{\pi}(1)$ is a border point of $K$ and that $E(\bm{\pi}(1)|\bm{\pi}(0))=\bm{\pi}(0)$. After this first step, $\bm{\pi}(1)$ belongs to a facet of $C$ and a unit $i_1$ is assigned either to the treatment either to the control group. Units $i\neq i_0$ remain unassigned and we have $\sum_{i:i\neq i_0}\pi_i(1)=n_T-D_{i_0}$. We can then replicate the first step after replacing $n_T$ by $n_T-D_{i_1}$ the sample $\{1,...,n\}$ by $\{1,...,n\}\backslash\{i_0\}$ and to allocate a second unit and to update assignment probability as $\bm{\pi}(2)$. At step $n-1$, $\bm{\pi}(n-1)$ belongs to the extreme points of $K$, this ends the flight phase. If, for instance, $\pi_i=1/2$ and $n$ is even, the extreme points of $K$ are some vertices of $C$, so the assignment is achieved. Now imagine that one has 101 units to assign with equal probability to the treatment and control groups. exact balancing on the two group sizes is not possible: 101 is an odd integer and it is not feasible to assign 50.5 units to the treatment. A popular solution is to consider $\pi_i=50/101$ or $\pi_i=51/101$ and to sample randomly 50 (or 51) elements among the 101 units. However, this strategy does not accommodate easily with heterogeneous probabilities of assignment and does not generalize to take into account many balancing constraints. With the cube method described above, for each step $t$ of the flight phase we have $\sum_i \pi_i(t)=50.5$ and the extreme points of $K$ are not anymore vertices of $C$. In that case, at the end of the flight phase, $n-1=100$ units are assigned at the end of the flight phase with $(n-1)/2=50$ units to the treatment and $(n-1)/2=50$ units to the control. The cube method can be completed with a last phase that randomly assigns to the treatment of the control the remaining unit ensuring that $n_T=50$ or $51$ and $E(n_T)=50.5$. In that case, the sizes of treatment and control groups are not exactly fixed but almost fixed (in fact as fixed as possible as soon as we respect the initial assignment probabilities $\pi_i=1/2$). This second phase is called landing phase. These two phases, the flight phase and the landing phase, can be generalized to the case where the researcher wants to impose several balancing constraints and heterogeneous probabilities of treatment.

Let us describe the cube method with an arbitrary number of balancing constraints defined by $q$ variables $Z_i$. A point $\mathbf{s}\in C$ will satisfy an equation analog to (ref) if

Let $A_i=\frac{Z_i}{\pi_i}$ and $A=(A_1,...,A_n)$ the matrix of size $q\times n$. Then (ref) is equivalent to

$K=C\cap Q$ is, therefore, the $(n-q)$-polytope that contains all the points in $C$ such that (ref) holds. At the first step, one chooses a random direction in $v\in \ker(A)$ and we select the unique $\lambda>0$ such that $\bm{\pi}(1):=\bm{\pi}(0)+\lambda v$ is on a facet of $K$ and that $\mathbb{E}[\boldsymbol{\pi}(1)|\boldsymbol{\pi}(0)]=\boldsymbol{\pi}(0)$. Because any facet of $K$ is the intersection of a facet of $C$ with $Q$, a component $i_0$ of $\boldsymbol{\pi}(1)$ is 0 or 1 and defining $D_{i_0}=\pi_{i_0}(1)$ one has assigned a first unit. Next, one applies a similar step for the facet of $K$ instead of $K$ and $\bm{\pi}(1)$ as a starting point instead of $\bm{\pi}(0)$. After $n-q$ steps, one has reached a vertex of $K$. This process corresponds to the flight phase in deville_efficient_2004. If this vertex of $K$ is also a vertex of $C$, the flight phase allocates every unit, and the two groups are exactly balanced (see Figures (ref) and (ref)). But in many cases, the vertex of $K$ is not a vertex of $C$, and there remain at most $q$ units to assign during the landing phase deville_efficient_2004 (see Figure (ref)).

Say that at the end of the flight phase, one has not assigned $r\leq q$ units and let $\boldsymbol{\pi^\ast}=\boldsymbol{\pi}(n-q)$ be the updated treatment probabilities at this stage. The landing phase of the cube method assigns the $r$ missing units such that $\mathbb{E}[D_i | \boldsymbol{\pi^\ast}]=\boldsymbol{\pi^\ast}$. grafstrom_doubly_2013 describe two methods for the landing phase (these are also the options used in sampling packages): (i) Linear programming: one considers all the $2^r$ allocations for these units and assigns probabilities to each allocation to minimize a cost function and satisfy $\mathbb{E}[D_i | \boldsymbol{\pi^\ast}]=\boldsymbol{\pi^\ast}$. Sampling probabilities are chosen to minimize $$\mathbb{E}\left(\sum_{i\notin S}Z'_i(D_i-\pi_i^{\ast})M\sum_{i\notin S}Z_i(D_i-\pi_i^{\ast})\big|W\right),$$ where $S$ is the set of units allocated at the flight phase, $W=(S, (D_i)_{i\in S}, (\pi_i^{\ast})_{i\notin S}, (Z_{i})_{i=1,...,n})$, and $M$ is a symmetric positive-definite matrix $q\times q.$ Common choices for $M$ are the identity matrix or the inverse-covariance matrix. After solving this minimization problem, the researcher randomly draws an allocation using these probabilities. (ii) Suppression of variables: if $r>20$, solving a linear problem becomes computationally difficult. In that case, at the end of the flight phase, one can drop a covariate (i.e., a constraint) and continue with the flight phase. One can thus successively drop variables until attaining a vertex of $C$. This method, however, implies that the researcher has to define an order to drop the covariates, ideally from the least to the most important.

\@startsection{section}{2}{0mm}{-1.5\baselineskip}{1\baselineskip}{\normalfont}{Statistical Properties of the Cube Method} This Section shows how the cube method allows obtaining an “almost-exact balance” between the treatment and control groups and relates this balance to gains in precision for treatment effect estimators.

\@startsection{subsection}{2}{0mm}{-1.2\baselineskip}{1\baselineskip}{\normalfont}{Balancing Approximations}

\@startsection{subsubsection}{3}{0mm}{-0.8\baselineskip}{0.4\baselineskip}{\normalfont}{Balance Checks}

As explained above, designing an allocation mechanism that always produces exactly-balanced groups is generally impossible. However, we here prove that the cube method is successful, under certain conditions, in creating almost-exactly-balanced samples in the sense of Equation (ref).

To check balance properties after allocating individuals according to the design $\Pi$, researchers are interested in computing the difference

Because $\mathbb{P}_{\Pi}(D_i=1|(X_{i'})_{i'=1,...,n})=\pi_i$, we have $\mathbb{E}\left(\Delta^{\Pi}_{j,n}\right)=0$ and under weak conditions on $\Pi$, we have

where $\mathbb{V}(\Delta^{\Pi}_{j,n})$ is an asymptotic variance depending on $\Pi$ and the distribution of $X$.

For the so-called baseline balance tests, researchers often consider the $t$-statistic

where $\widehat{\mathbb{V}}(\Delta^{\Pi}_{j,n})$ is a consistent estimator of the asymptotic variance of $\Delta^{\Pi}_{j,n}$ to test the null hypothesis of exact balance. $t^{\Pi}_{j,n}$ is then associated to a $p$-value $p^{\Pi}_{j,n}$ which take values between 0 and 1. As explained by snyder_examining_2024, when creating balance tests for RCTs, small $p$-values (below $0.15$) are considered problematic and are usually underreported.

Let us first consider a naive mechanism that does not use baseline information to assign units. Such situations correspond to the case where the design $\Pi$ is a coin toss, i.e., a design where each unit $i$ is allocated to the treatment independently of the allocation of other units: $$\mathbb{P}_{\Pi}\left(\bigcap_{i=1}^n\{D_i=d_i\}\big|(X_{i})_{i=1}^n\right)=\prod_{i=1}^n\pi_i^{d_i}(1-\pi_i)^{1-d_i}.$$ A coin toss does not balance any variable nor group sizes.

When $\pi_i=\frac{n_T}{n}$ for any $i$, another popular design is sampling without replacement of $n_T$ treated units, also known as complete randomization: $$\mathbb{P}_{\Pi}\left(\bigcap_{i=1}^n\{D_i=d_i\}\big|(X_{i})_{i=1}^n\right)=\binom{n}{n_T}^{-1}\mathds{1}\left\{\sum_{i=1}^nd_i=n_T\right\}.$$ Complete randomization only balances constant variables. In that case, the sample of treated and control groups are fixed, and the design is also balanced on the constant $\big(\sum_{i=1}^n D_i=n_T$, $\sum_{i=1}^n (1-D_i)=n-n_T$ and $\sum_{i=1}^n\frac{D_i}{\pi_i}=n \big)$.

Under such assignments and Assumptions (ref) and (ref), and more generally for any design $\Pi$ such that (ref) holds with $\mathbb{V}(\Delta_{j,n}^{\Pi})>0$, we have $\Delta_{j,n}^{\Pi}=O_p\left(\frac{1}{\sqrt{n}}\right)$, $t_{j,n}^{\Pi}\stackrel{d}{\longrightarrow} \mathcal{N}\left(0,1\right)$ and $p_{j,n}^{\Pi}\stackrel{d}{\longrightarrow}\mathcal{U}(0,1)$. This result means that if one randomizes naively, control and treatment groups will present imbalances with a strictly-positive probability. Moreover, for a confidence level of $100(1-\alpha)\%$, there exists always $100\alpha\%$ chance of obtaining significant differences. If an researcher evaluates the balance of 10 independent covariates at the $85\%$ confidence level snyder_examining_2024, there is more than $80\%$ chance of having at least one significant difference. This magnitude questions the mere implementation of such widely used tests. Even if a multiple F-test with a confidence level of $85\%$ mitigates this rejection rate, the null hypothesis of simultaneously balanced covariates is rejected by construction with a 15% chance.

The cube method ensures that these tests are unnecessary since we can balance control and treatment groups in any covariate $(X_j)_{j=1,...,p}$. This is achieved because $\mathbb{V}(\Delta^{\Pi}_{j,n})=0$ for any $j=1,...,p$ in (ref). Performing these tests would not make sense since we never reject the null hypothesis by construction. However, one might report them if the editor worries about researchers randomizing badly. Usual balancing strategies are stratified or matched-pairs designs. These methods ensure $\mathbb{V}(\Delta^{\Pi}_{j,n})=0$ if the covariates $(X_j)_{j=1,...,p}$ are all discrete but will always generate imbalances for continuous ones since the researcher needs to discretize or aggregate them before randomizing.

The following proposition explains how the balancing approximations are satisfied with the cube method. Because the number $q$ of balancing constraints in Equation (ref) could be large with the cube method, we are also explicit on how $q$ affects balancing approximations to allow us to consider a framework where $q$ tends to $\infty$.

Proposition (ref) shows that, as $n$ grows, the cube method ensures the balancing Equation (ref) as soon as the second-order moments of $X$ exist. Furthermore, if moments of order $r>2$ exist for $X$, (ref) holds as soon as $q=O\left(n^{\frac{1}{2}-\frac{1}{r}}\right)$. $q$ can even be $o\left(\sqrt{\frac{n}{\ln(n)}}\right)$ if the covariates $X$ are all sub-Gaussian or $o(\sqrt{n})$ if they are bounded. This means that with probability tending to one, the $p$-values of balance tests tend to 1. Balance is thus never rejected for large $n$ contrary to randomization under a design $\Pi$ such that (ref) does not hold.

\@startsection{subsubsection}{3}{0mm}{-0.8\baselineskip}{0.4\baselineskip}{\normalfont}{Comparison with Other Methods}

We here compare the balancing properties of the cube method with other randomization methods. For the sake of simplicity, we fix $\pi_i=1/2$. We check imbalances by using the Horvitz-Thompson estimators for the average difference between the control and treatment groups $B_{n,p}(X)=\frac{2}{n}\sum_{i=1}^nX_iD_i-X_i(1-D_i)$ and looking at their squared Euclidean norm $||B_{n,p}(X)||^2=\frac{4}{n^2}\sum_{j=1}^p\left(\sum_{i=1}^nX_{ji}D_i-X_{ji}(1-D_i)\right)^2.$

Assumption (ref) imposes mild conditions over the baseline covariates. In particular, the components of the vector $X_i$ are not assumed to be independent. Figure (ref) illustrates our main results for a simple case where this assumption holds: $(X_{ji})_{j=1,\ldots,p,i=1\ldots,n}$ are independent and follow a uniform distribution on $[0,1].$ In this case, we have $V(X_{ji})=1/12$, $\mathbb{E}(X_{ji}^2)=1/3$, and $\underline{C}=\overline{C}=1$.

A coin toss ensures that treatment and control groups are balanced on average (i.e., $\mathbb{E}[B_{n,p}]$=0). However, we can still have imbalances between groups for a given allocation. Additionally, a coin toss will often generate different sizes between the treatment and control groups, meaning that it fails to balance on a constant. Proposition (ref) in Appendix (ref) shows that for a coin toss, $\frac{4\underline{C}}{3}\frac{p}{n}\leq\mathbb{E}[||B_{n,p}(X)||^2]\leq\frac{4\overline{C}}{3}\frac{p}{n}$. For the case illustrated in Figure (ref), we thus have $\mathbb{E}[||B_{n,p}(X)||^2]=\frac{4p}{3n}$.

Complete randomization improves upon the coin toss procedure. In this method, the researcher fixes the group sizes. Since we here assume $\pi_i=1/2$, the researcher randomly chooses an allocation among those having an equal number of treated and untreated units. The researcher still does not use any information on the baseline covariates to refine the randomization process but manages to reduce the imbalances due to different group sizes. Indeed, under Assumption (ref) and complete randomization $\frac{\underline{C}}{3}\frac{p}{n}\leq\mathbb{E}[||B_{n,p}(X)||^2]\leq\frac{\overline{C}}{3}\frac{p}{n}$, so bounds reduce by four, relative to the coin toss. For the example in Figure (ref), we have $\mathbb{E}[||B_{n,p}(X)||^2]=\frac{p}{3n}$. This result clearly shows the advantages of using designs with fixed sample sizes such as complete randomization, but also the cube method with the constraints in Section (ref) or matched-pairs design.\footnote{One can show that the gains from fixed group sizes are not present if one uses a difference-in-means estimator instead. Moreover, in the case of $n$ even and $\pi_i=1/2$, this estimator is equivalent to the Horvitz-Thompson for complete randomization, matched-pairs design, and the cube method. However, it can lead to more precise estimates for coin tosses or stratified designs.}

In sharp contrast with naïve methods, covariate-adaptive randomization uses baseline information to improve balance between treatment and control. Stratified designs are the most used and studied covariate-adaptive method. Stratification has a long tradition in RCTs fisher_design_1935,higgins_improving_2016. Stratified designs are the most popular assignment mechanisms used in RCTs as they are simple to grasp and can produce balanced samples. This method consists of using one or several baseline variables to create blocks or strata and then using complete randomization inside each stratum. A common practice in experiments is to block on gender, meaning that randomization is performed independently amongst male and female units, generating the same proportion of men and women in each treatment arm. When using dummy variables to define the strata, stratified or blocked randomization allows almost exact balancing of the variables used to create them. athey_chapter_2017 recommend balancing on small strata since this method generates substantial precision gains. However, stratified designs do not come without any limitations. Notably, the type and number of covariates that one wants to balance can impose some difficulties. Facing continuous covariates, such as income or grades, makes it impossible to stratify without the researcher deciding how to create the strata. Assumption (ref) imposes continuous covariates. We thus will focus on two ways of generating (possibly-)small strata, discretization and matched pairs.

First, the researcher can discretize continuous variables using $\ell$-quantiles for each covariate, generating thus $\ell^p$ strata. Stratifying will produce balance gains as long as the number of units remains large compared to the number of strata. In particular, we show in Proposition (ref) that whenever $n\ell^{-p}\to \infty$, stratified designs through discretization outperform complete randomization. Discretizing baseline covariates, however, does not ensure fixed sizes for each stratum. In particular, if the number of strata is big compared to the sample size, there is a big chance of having some strata with only one unit. In the limit case where $n\ell^{-p}\to 0$, every non-empty strata has one unit with probability one, and stratifying through discretization performs strictly worse than complete randomization and approximates a coin toss. We also show that, in both limit cases, imbalances grow at a rate of $p/n$. We observe this behavior in Figure (ref) for $\ell=2,4$. In this example, the stratified designs perform better than complete randomization whenever $l^p\leq n/2$ and similarly to a coin toss for $l^p\geq 32n$. Balancing deterioration can thus occur quite rapidly when stratifying is done by discretizing many continuous variables. It is worth noticing that this issue is not exclusive to continuous variables, as it arises when stratifying using many categorical variables. In that case, the number of strata equates to the product of covariate support cardinalities.

To eliminate the issue of single-unit strata, researchers may use a more sophisticated way of creating their strata: matched-pair designs. Following greevy_optimal_2004,bai_inference_2022,bai_optimality_2022, the researcher can create $n/2$ strata of two units to minimize the average intra-strata distance. The researcher thus creates pairs of two units that resemble each other. After constructing these strata, the researcher randomly allocates one to treatment. By doing so, she creates control and treatment groups that are very similar. We show in Proposition (ref) that this design always outperforms complete randomization. However, this type of strata construction works by trying to have a similar joint distribution of $X$ between treatment and control. This approach to balancing implies that for large $p$, it becomes more difficult to find pairs of units close to each other. In particular, we show that under this design and Assumption (ref), $\mathbb{E}[||B_{n,p}(X)||^2]\geq \frac{p}{n}\left(\frac{1}{3}-\sqrt{\frac{2\ln(n-1)+4\ln\overline{C}}{p}}\right)$. This result entails that the number of balancing covariates $p$ is large relative to $\ln n$, balance gains shrink, and imbalances increase at the rate of $p/n$. Figure (ref) illustrates this effect. Indeed, when $p$ becomes larger than $\ln 500\approx6$, the matched-pairs design performs better than complete randomization, but its relative gains quickly reduce.

For all these randomization methods but matched pairs, imbalances grow at the rate of $p/n$. For matched pairs, this rate holds whenever $\ln n=o(p).$ We now show that the cube method is less concerned by this curse of dimensionality since imbalances grow at a rate $p^2/n^2$. If $p$ remains smaller than $n$, then this result implies a much slower balance deterioration than for the methods described above. Proposition (ref) gives this upper bound for $\mathbb{E}[||B_{n,p}(X)||^2||]$ when using the cube method. This proposition is also stated and proved in Appendix (ref).

The upper bound depends on the matrix $M$, described in equation Equation (7) in deville_efficient_2004, used during the landing phase. Notably, one can take $M$ the identity matrix, and we have $\frac{\lambda_{max}(M)}{\lambda_{min}(M)}=1$. Then, we see that the cube method outperforms other methods that grow at a rate of $p/n$. This is clearly illustrated in Figure (ref), where we see that imbalances increase only very lightly on the number of covariates when randomizing with the cube method. The main difference between the cube and other designs is that it balances selected moments of the covariates instead of balancing the whole joint distribution of $X$. It thus reduces the burden of balancing a higher number of covariates. Balancing moments can also be achieved through other methods. In particular, we can perform re-randomization such that the stopping criterion requires balancing moments of $X$ or perform a Gram-Schmidt walk design harshaw_balancing_2024.

Re-randomization is another method that allows obtaining balance between covariates that has gained focus in the last decades morgan_rerandomization_2012,li_asymptotic_2018,imbens_experimental_2011. The main idea of re-randomization is to completely randomize repeatedly until the obtention of balanced groups. Some researchers perform re-randomization without prespecifying it. This repetition affects treatment probabilities in an unknown manner, which induces invalid inference bruhn_pursuit_2009,athey_chapter_2017. There are, however, several ways of performing re-randomization that allow valid inferences to some extent. Most of them rely on the simulation of the distribution under the re-randomization procedure used to assign units. This implies that researcher should draw a large number $N$ of balanced samples. One can keep randomizing until $||B_{n,p}(X)||^2\leq 4\frac{(p+1)^2}{n^2}$ to achieve the upper bound of Proposition (ref) (with $M=Id$). However, this upper bound is not sharp and to compare re-randomization with the cube method, we counted how many times an researcher should sample with naive randomization to get $N=1000$ samples that are balanced as well as $B^{\ast 2}=\mathbb{E}\left(||B_{n,p}(X)||^2\right)$, where the previous expectation is computed for the cube method through simulations. Under the design used in Figure (ref), and for complete randomization, $12/4\times n\times||B_{n,p}(X)||^2$ converges in distribution to $\chi^2(p)$. Next, the probability to achieve balancing as good as the cube is $F_{\chi^2(p)}(12/4*n*B^{\ast2})$. To have $N=1000$ samples balanced as well as the cube, researchers thus have to sample approximately $1000/F_{\chi^2(p)}(3*n*B^{\ast2})$. For $p=3$, the researcher have to sample more than $10^6$ samples and for $p=10$ this is more than $9,98\times 10^{11}$.

The probability of getting a sample that has the same properties as the cube method becomes small very quickly, so it becomes demanding computationally, in particular, if one wants several allocations to perform randomization-based inference.

harshaw_balancing_2024 recently developed the Gram-Schmidt walk design to obtain balanced groups in RCTs. As the authors consider a tradeoff between balance and robustness, they impose a choice of parameter $\phi\in(0,1]$. For $\phi=1$, the algorithm from harshaw_balancing_2024 reduces to a coin toss and next $\frac{4\underline{C}p}{3n}\leq \mathbb{E}\left(||B_{n,p}(X)||^2\right)\leq \frac{4\overline{C}p}{3n}$ under Assumption (ref). For $\phi\in(0;1)$, we conjecture that $\underline{K}_1 \phi \frac{p}{n}+\underline{K}_0 (1-\phi)\frac{p^2}{n^2}\leq\mathbb{E}\left(||B_{n,p}(X)||^2\right)\leq \overline{K}_1 \phi \frac{p}{n}+\overline{K}_0 (1-\phi)\frac{p^2}{n^2}$ for some constant $\underline{K}_0,\underline{K}_1,\overline{K}_0,\overline{K}_1$. The balance of the Gram-Schmidt method increases as $\phi$ tends to zero. However, theoretical results in harshaw_balancing_2024 and implementation of the Julia package only hold for $\phi$ positive. The choice of $\phi$ is thus critical but difficult to justify. On a side note, the cube method does not require any (subjective) parameter choice.

\@startsection{subsection}{2}{0mm}{-1.2\baselineskip}{1\baselineskip}{\normalfont}{Variance Reduction}

The balance between covariates in the control and treatment groups is also beneficial if these variables are related to the potential outcomes. In this case, using the cube method will also reduce the variance of the Horvitz-Thompson and Hájek estimators.

Assumption (ref) states that potential outcomes are linearly related to observable covariates. However, we allow heterogeneity in treatment effects by specifying different equations for control and treatment groups.

This conjecture establishes that as $n$ increases, the cube method tends to Poisson sampling. As $n$ goes to infinity, the dependence between the assignment of a finite number of individuals disappears. We draw this conjecture from results in deville_variance_2005 and simulations that confirm it.

To have a benchmark for the gains in variance decline, we compare the cube method with a coin toss.

Proposition (ref) shows the gain in asymptotic variance from balancing covariates using the cube method. The reduction is more substantial when $X$ explains more of the potential outcomes. Estimates of the ATE are thus more precise when using the cube method. This reduction can represent significantly lower costs when conducting an RCT. Notice that under the same set of assumptions, $V_0^\ast$ corresponds to the semiparametric efficiency bound in hahn_role_1998. Simulations in Sections (ref) illustrate these gains.

\@startsection{subsection}{2}{0mm}{-1.2\baselineskip}{1\baselineskip}{\normalfont}{Inference} This section provides properties of the cube algorithm and methods to perform inference. We elicit two main techniques of conducting inference, one based on the asymptotic properties of the HT estimator and the other based on the randomization mechanism.

\@startsection{subsubsection}{3}{0mm}{-0.8\baselineskip}{0.4\baselineskip}{\normalfont}{Asymptotics-based Inference}

Some methods, such as re-randomization, alter the inclusion probabilities in a manner that is unclear to the researcher imbens_experimental_2011. When the criterion for selection is known and behaves in a known way, such as the Mahalanobis distance, one can perform conservative inference. However, balance is imperfect for numerous covariates. Since the cube method assigns treatment only once, we can perform asymptotic-based inference. We here give the asymptotic properties and propose an easy way to construct exact confidence intervals.

To construct a confidence interval, one would like to estimate either $V_0$ or $V_0^\ast$. Estimating $V_0/n$ is impossible without making assumptions on the relation between $\varepsilon_i(1)$ and $\varepsilon_i(0)$. This issue is common in RCTs. We can, nonetheless, easily construct an unbiased estimator $\widehat{V}$ for $V_0^\ast/n$. Let $\widehat{\beta}_d$ and $\widehat{\varepsilon}_i(d)$ be the estimated coefficients and residuals, respectively, of a regression of $Y_i(d)$ on $Z_{di}$, for $d\in\{0,1\}$. We then have

\\ with $\widehat{\Omega}=\frac{1}{n-1}\sum_{i=1}^n\left(Z_{1i}\widehat{\beta_1}-Z_{0i}\widehat{\beta_0}-\frac{1}{n}\sum_{i'=1}^n\left(Z_{1i'}\widehat{\beta_1}-Z_{0i'}\widehat{\beta_0}\right)\right)^2.$

Then, we can test the weak hypothesis

and construct the confidence interval based on

\\ In Section (ref), we perform simulations that confirm the exact coverage rate of this confidence interval when $n$ is big enough.

\@startsection{subsubsection}{3}{0mm}{-0.8\baselineskip}{0.4\baselineskip}{\normalfont}{Randomized-based Inference}

We here study the properties of randomization-based inference when permuting treatment status while satisfying balancing constraints. For these tests, we consider the stronger null hypothesis:

Notice that testing this hypothesis, under Assumptions (ref) and (ref) is equivalent to testing $(Y_i)_{i=1}^n\perp \!\!\! \perp (D_i)_{i=1}^n | X_1,\ldots X_n$ (Proof in Appendix (ref)).\\ To explain the test, we introduce some new notation. Let $G_n$ be the set of all possible $2^n$ assignments. Then, we can define the set of assignments $G_n^{cube}\subseteq G_n$ satisfying the constraints imposed by the cube method. That is, with Assumptions (ref) and (ref),

$$G_n^{cube}=\left\{g\in G_n : \Delta_{j,n}=o_p\left(\frac{q}{\sqrt{n}}\right) \text{ for } 1\leq j\leq p\right\}.$$\\ We note $\mathbf{P_n}=(Y_i,D_i, X_i)_{i=1}^n$ the observed values, and $\mathbf{P_n^{(g)}}=(Y_i,D_i^{(g)}, X_i)_{i=1}^n$, the new data where we have reassigned treatment according to $g\in G_n^{cube}$. For computational facility, we can replace $G_n^{cube}$ by $G_n^B=\{g_1,\ldots,g_B\}$, such that $g_1$ is the assignment really obtained and $(g_i)_{i=1}^B$ are drawn independently from a uniform distribution on $G_ n^{cube}$.\\ Then, for a given test statistic $T_n(\mathbf{P_n})$,we consider the test

$$\phi^{rand}(\mathbf{P_n})=\mathbbm{1}\left\{T_n(\mathbf{P_n})>c_n(\mathbf{P_n},1-\alpha)\right\}$$ with $$c_n(\mathbf{P_n},1-\alpha)=\inf\left\{t\in \mathbb{R} : \frac{1}{B}\sum_{g\in G_n^{B}}\mathbbm{1}\{T_n(\mathbf{P_n^{(g)}}) \leq t\}\geq 1-\alpha\right\}.$$

Proposition (ref) indicates that if $T_n(\mathbf{P_n}) > c_n(\mathbf{P_n},1-\alpha)$, we reject the null hypothesis (ref) at the $\alpha$ level. The proof is similar to previous results on other covariate-adaptive assignment mechanisms heckman_analyzing_2010,heckman_inference_2011,lee_multiple_2014,bai_inference_2022, but it is presented for completeness. This proposition ensures that we can compute Fisher's $p$-values by comparing our test statistic with those produced by other assignments made by the Cube method.

\@startsection{section}{2}{0mm}{-1.5\baselineskip}{1\baselineskip}{\normalfont}{Simulations} This section compares the cube method to other randomization methods by performing Monte Carlo simulations. We are interested in examining the impact of introducing new covariates in the variance of treatment effect estimates. For this purpose, we evaluate different randomization methods using one example from data following a simple DGP in the spirit of Figure (ref) and another using data-driven methods from an empirical application. We use packages in R for randomization available here: https://rdrr.io/cran/BalancedSampling/. Interestingly enough, the cube method is not computationally demanding. While running simulations included in the present paper, we systematically found the cube method to have an order of magnitude faster than the algorithms we used to implement other presented methods.

\@startsection{subsection}{2}{0mm}{-1.2\baselineskip}{1\baselineskip}{\normalfont}{Simple DGP} For $k=1,\ldots,K$ the number of iterations, $j=1,\ldots,p$ the number of covariates, and $i=1,\ldots,n$, the number of observations, we independently draw $X_{jik}\sim\mathcal{U}(0,1)$ and $\varepsilon_{ik}(d)\sim\mathcal{N}(0,1) $, for $d=0,1$. We then generate the potential outcomes $Y_{ik}(0)=1+ (X_{ik}-1/2)'\beta_0+\varepsilon_{ik}(0)$ and $Y_{ik}(1)=1+X_{ik}'\beta_1+(X_{ik}-1/2)'A(X_{ik}-1/2)+\varepsilon_{ik}(1),$ with $A=(1/20)\times(\mathbbm{11'}-\operatorname{diag}(1))$. Notice that, in this example, $\theta_0^\ast=0.$

We consider $n=500$, $p=30$, $\beta_0=(1,\boldsymbol{0}')'$, $\beta_1=2\beta_0$, so only one covariate and noise explain variations in the individual treatment effect. We assume that the researcher knows she should always balance this covariate. Still, she does not have previous information about the (un)informativeness of the 29 other covariates. In these simulations, the researcher has to choose which simulation method she uses and how many covariates to include. That choice corresponds to an assignment design $\Pi$ and generates treatment statuses $D_{ik}^\Pi$. We estimate the PATE using the HT estimator $\widehat{\theta}^\Pi_{HT,k}$. To evaluate the precision entailed by the assignment design, we perform $K=5,000$ simulations and compute the standard deviation of the estimator over the simulations. Since the PATE is null, this is equivalent to estimating the root mean square error (RMSE).

Figure (ref) shows the RMSE of the HT estimator by number of covariates and randomization method. The simulations show that the cube method is always competitive. Since including more covariates deteriorates balancing only very lightly, precision gains are maintained even when $p=30$. This behavior is not present for other randomization methods. Indeed, stratification using the median (quartile) leads to worse precision than complete randomization as soon as $p>6$ ($p>3)$ and converges to a coin toss for $p>12$ ($p>6)$. Moreover, using a matched-pairs design improves from complete randomization but, for $p>3$, underperforms compared to the cube method: when $p$ increases, precision for matched-pairs design worsens, whereas it remains the same when using the cube. By allowing an abundant set of covariates, the cube method improves the exploitation of balancing gains, even when the researcher chooses to balance covariates that are not explicative of potential outcomes. This behavior could arise if the researcher is interested in several treatment outcomes and collects their pre-treatment values. Then, she would ideally want to balance them all, even if only one covariate is explicative of one outcome. As described through these simulations, the cube method ensures precision gains for a particular outcome variable, even when balancing another 29 baseline variables.

\@startsection{subsection}{2}{0mm}{-1.2\baselineskip}{1\baselineskip}{\normalfont}{Empirical Data} We further illustrate the properties of the cube method by using experimental data from gerber_one_2020. The authors investigate how informing potential voters about the closeness of an election affects their beliefs and voting behavior. Since the experimental data only represents one of many possible samplings, we proceed by generating a superpopulation. We create a large dictionary with baseline outcomes, covariates, and demographics. We consider all possible interactions and second-order polynomials. We thus generate a dataset of 6,424 observations and 7,381 covariates (hereon denoted by $X$), with 3,193 individuals in the treatment group. We consider beliefs about the closeness of the election as the main outcome $Y$. We ran two lasso regressions separately for treated and control units to train two models, $f_1$ and $f_0$. We then estimate $s_1^2=\widehat{\mathbb{V}}(Y-f_1(X)|D=1)$ and $s_0^2=\widehat{\mathbb{V}}(Y-f_0(X)|D=0)$. To generate the superpopulation we draw $N=50,000$ individuals, with replacement and we generate $Y_i(1)=f_1(X_i)+\varepsilon_i(1)$ and $Y_i(0)=f_0(X_i)+\varepsilon_i(0)$ for $i=1,\ldots,N$ with $(\varepsilon_i(1);\varepsilon_i(0))\sim\mathcal{N}\left((0~,~0) , (s_1^2 ~ ~0.5s_1s_0~,~0.5s_1s_0 ~~ s_0^2 )\right).$ We thus obtain a superpopulation $(X_i, Y_i(1), Y_i(0))_{i})_{i=1,\ldots,N}$.

We then run $K=10,000$ Monte Carlo simulations, where for every iteration, we draw $n\in(100,256,500,1000)$\footnote{We select 256 because exact inference methods for matched pair designs require a sample size divisible by four bai_inference_2022.} individuals, allocate them according to five treatment allocation methods: complete randomization, stratified randomization using median values for continuous variables, matched-pairs design using the Mahalanobis distance when balancing multiple covariates, and the cube method with the two first moments per variable. For stratified designs, matched-pairs design, and the cube method, we balance between 1 and 12 covariates. When balancing only one, we use the pre-treatment value of $Y$. For the 12 covariates, we consider five pre-treatment outcomes and seven baseline covariates. We always prioritize pre-treatment outcomes as they are likely the most explicative variable for their post-treatment counterpart. We set $\pi_i=\frac{1}{2}$. For complete randomization, matched pairs, and the cube method, we compute the HT estimator. For these methods, we compute confidence intervals using, respectively, White standard errors, Equation (14) in bai_optimality_2022, and Equation (ref) in Section (ref) above. For stratification, since $\pi_i=\frac{1}{2}$, we use an OLS regression with strata fixed-effects, which gives consistent estimators and exact inference as shown by bugni_inference_2018.

Table (ref) reports estimators of the effective sample size $\operatorname{ESS}=\frac{\mathbb{V}(\widehat{\theta}_{HT}^\Pi)}{\mathbb{V}(\widehat{\theta}_{HT}^{\operatorname{CR}})}\times n$ of each design, based on the variance of the estimators across the $K$ iterations. The ESS indicates, for every allocation design, the experimental sample size required to estimate the treatment effect with the same precision as with complete randomization. We see that almost every covariate-adaptive method does better than complete randomization, as they allow to reduce the sample size, often substantially. The only exception is stratification with many covariates. In general, if $n<2^p$, stratification becomes worse than complete randomization. This phenomenon is due to the small strata issue and worsens with the strata fixed-effects estimator. Across different $n$ sizes, we see that the matched-pairs design does better than the cube method when we balance a few covariates (up to three, in general). However, once balancing more covariates, the cube method becomes more efficient. As expected, these relative gains to the matched-pairs design are more apparent for smaller $n$ since the curse of dimensionality is more stringent. Notice there are clear gains from using the cube method even when allowing for non-linearities in the imputation of $Y(1)$ and $Y(0)$. These results and those in Section (ref), along with coverage rates displayed in Table (ref) in the Appendix, constitute suggestive evidence for a possible relaxation of Assumption (ref).

Tables (ref)-(ref) in the Appendix show additional results for each design: standard deviation and bias of the HT estimator, confidence coverage rate, and power for testing a null PATE. In particular, Table (ref) verifies that the coverage rate is exact for $n$ large enough.

\@startsection{section}{2}{0mm}{-1.5\baselineskip}{1\baselineskip}{\normalfont}{Discussion}

Although RCTs are now common in economics, there are still some aspects of randomization that can be further clarified. We here discuss two of these, in light of the contribution of the cube method. We first show that, despite being widely available, pre-treatment information is still underused. The cube method could improve practices. We also discuss how pre-analysis plans are often vague about which covariates should be included. In a second section, we explain how the cube method may help suppress a publication bias arising from imbalances.

\@startsection{subsection}{2}{0mm}{-1.2\baselineskip}{1\baselineskip}{\normalfont}{On the use of balancing methods in RCTs}

bai_optimality_2022 indicates that among 5,000 RCTs in the AEA RCT Registry, more than 800 are stratified (i.e., about 16%). We complement this insight by gathering information from 104 randomized controlled trials (RCTs) published in top-5 and AEA journals\footnote{American Economic Review, AEJ: Applied Economics, AEJ: Economic Policy, AEJ: Macroeconomics, AEJ: Microeconomics, Econometrica, Journal of Political Economy, Quarterly Journal of Economics, and Review of Economic Studies} between 2019 and 2023. Specifically, we examined their collection of baseline information, baseline outcomes, and randomization methods. Figure (ref) summarizes these details. Our findings indicate a lack of consensus among RCTs regarding the method used for allocating individuals to treatment. Most published papers (54%) employ a stratified design, followed by completely randomized designs (34%). A minority of researchers utilize alternative methods such as matching or re-randomization. However, there is substantial agreement regarding the collection of baseline data, with 90% of papers gathering information before treatment allocation. Nevertheless, this information is not always utilized during the randomization process, as only 46% of the papers leverage it to achieve covariate balance. The remaining studies collect baseline data for balance tests, covariate adjustment in regression, and/or heterogeneity analysis of treatment effects. If outcomes of interest are relatively stable over time, researchers should be interested in balancing their pre-treatment values, as they are highly likely to be correlated with potential outcomes. In our sample, 65% of researchers collect these variables, yet only 23% incorporate them into the allocation design, indicating an area for improvement in experimental design and inference. When various outcomes are considered in RCT, the curse of dimensionality arising in stratification, matched pair design, or re-randomization (where computational time could become prohibitive) may prevent researchers from balancing on a large set of pre-treatment outcomes and sociodemographic covariates. In view of results in Section (ref), the cube method could greatly improve experimental randomization by allowing balancing on more variables than other methods.

\@startsection{subsection}{2}{0mm}{-1.2\baselineskip}{1\baselineskip}{\normalfont}{“Balancing checks”: cube randomization and publication bias} Researchers routinely provide “balance tables”, i.e., a comparison of the moments of available covariates between control and treatment. Due to bad luck, the researcher can expect a certain proportion of significant imbalances. And the likelihood of imbalances increases with the number of covariates used. There seems to be a gray area around the reporting of balance checks. In particular, pre-analysis plans often report no clear justification regarding the choice of covariates to include in balance tables.

Researchers are not very comfortable with reporting substantial unbalances. snyder_examining_2024 analyzing a large set of balance checks find that the editorial process removes an ample part of studies reporting imbalances, perhaps as much as 30%. Imbalances increase the risk of rejection by journals and, even worse, may provide incentives to engage in p-hacking (e.g., removing some covariates from the balance table). Our Proposition (ref) ensures that cube randomization definitively solve these issues making balance checks superfluous.

\@startsection{subsection}{2}{0mm}{-1.2\baselineskip}{1\baselineskip}{\normalfont}{Future research and limitations}

To conclude this discussion, let us mention some limitations of our results and agenda for future research. We use an assumption of linearity in the conditional expectation of potential outcomes to establish some asymptotic results. Based on simulations, we conjecture that our results hold without linearity assumption, but we do not have a formal proof at the current stage.