\section{Introduction}\label{sec:intro}


The setting we consider is the two-arm non-sequential randomized study with a continuous endpoint (e.g. a pill-placebo double-blind clinical trial assessing blood glucose improvement) that seeks inference for an additive treatment effect. There are $2n$ individuals which are sample-size balanced between the two groups. Each subject $i$ is placed into the treatment or control group and this information is encoded in $w_i \in \braces{-1, +1}$ where -1 denotes assignment to control and +1 denotes assignment to treatment. The vector $\bv{w} := \bracks{w_1~\ldots~w_{2n}}$ is called an \textit{assignment}, \textit{allocation}, or a \textit{randomization}. The collection of the legal $\bv{w}$ vectors and the probability distribution on the legal vectors in the experimental setting is known as an experimental \textit{design}, \textit{strategy}, \textit{algorithm} or \textit{procedure}.

There are $p$ covariates with the $j$th's value for subject $i$ denoted $x_{ij}$ that are assumed known in advance of the randomization (the setting of known measurements is sometimes called \textit{offline} to distinguish it from the \textit{sequential} setting, the latter not being the setting considered in this paper). After any randomization, the values of the subjects' covariates in each group are approximately the same; and it is this fact that is the main reason randomization is employed when seeking causal inference \citep{Cornfield1959}. Thus, \textit{imbalance} (of which there are many principled metrics) in the covariate values between the treatment subjects and the control subjects should be small. Upon completion, the sample responses $y_i$, the result of an unknown process we call the \emph{response model}, are assessed. These sample responses are used to compute a causal estimate of the population average treatment effect $\beta_T$.

Our contribution is a new experimental design that provides (1) lower error in the measurement of the average treatment effect and (2) is robust when the response model is linear and/or non-linear. If the response model is linear, the best experimental design is one that optimizes the treatment assignment for minimal covariate imbalance which we call \emph{imbalance-optimizing designs}. If the response model is nonlinear, a good strategy is to optimize the treatment assignment via creating binary matches (i.e. nonbipartite pairings of subjects) with small intramatch covariate distance. If the response model is a combination of both linear and nonlinear, then both criteria (minimize covariate imbalance and matching) become important. Our hybrid design does both: first by creating optimal nonbipartite matches on covariate distance and then optimizing intramatch assignments to provide small covariate imbalance across the whole sample. The latter can be accomplished by many methods; we investigate rerandomization and a pairwise greedy-switching heuristic (the subject of our previous work) where this heuristic now operates on the differences of the covariate values within each match in contrast to the raw subject covariate values. 

To illustrate the advantages of our hybrid experimental design that both matches and imbalance-optimizes, consider linear and nonlinear response models, respectively $y_i = x_i + \frac{1}{2}\beta_T w_i + \mathcal{E}_i$ and $y = x^2 + \frac{1}{2} \beta_T w_i + \mathcal{E}$ where $\beta_T$ we set to 1, $x_i$'s are drawn iid from $U(0, 3)$ and the $\mathcal{E}_i$'s are drawn iid from $\normnot{0}{0.1^2}$. We use a small sample size of $2n = 40$ and consider three designs: (G) pure imbalance optimization via the greedy pair switching design which will be detailed in Section~\ref{sec:methodology}, (M) optimal non-bipartite matching and (MG) optimal non-bipartite matching followed by imbalance optimization via the greedy pair switching design. We simulate the above system 1,000 times and measure covariate imbalance using the metric log$_{10}(|\bar{\bv{x}}_T - \bar{\bv{x}}_C|)$ and squared error between the estimate $\ybar_T - \ybar_C$ and $\beta_T$. We average over the simulations to estimate both means and display the results in Table~\ref{tab:basic_demo}.

\begin{table}[h]
\centering
\begin{tabular}{cccc|ccc}
&&Average & Average & \multicolumn{3}{c}{p value for the average} \\
Response && Log$_{10}$  & Squared & \multicolumn{3}{c}{squared error comparison} \\
 Model & Design & Imbalance  &  Error & G & M & MG \\ \hline
Linear & G  & -3.90 & 0.00099 & - & $\approx 0$ & 0.74  \\
Linear & M  & -1.87 & 0.00135 & - & - & $\approx 0$  \\
Linear & MG  &  -4.55 & 0.00105 & - & - & - \\ \hline
Nonlinear & G &  - &0.04350 & - & $\approx 0$  & $\approx 0$ \\
Nonlinear & M &  - &0.00614 & - & - & 0.08\\
Nonlinear & MG &  - & 0.00273 & - & - & - \\
\end{tabular}
\caption{Results of the illustrative simulation explained in the text.}
\label{tab:basic_demo}
\end{table}

In the linear model, the pure imbalance-optimizing design (G) provides better MSE performance than matching (M) and employing the hybrid design combining matching with imbalance-optimization (MG) does not penalize performance (i.e. G and MG have MSE's that are statistically equal). When the model is non-linear, the covariate imbalance is less relevant, and  thus (M) provides much better MSE performance than (G) and the hybrid design (MG) outperforms (M) although not statistically significantly here (we will see in the results section that this advantage is likely only for very small sample sizes). Thus the hybrid (MG) design is robust and adaptive as it has the best performance in both linear and nonlinear settings.

When examining the imbalance, M achieves small covariate imbalance but not when compared to G / MG which is expected as that is not the former's primary objective in contrast to the latters' primary objective. MG has smaller covariate imbalance than G as it can combine the reduction in imbalance from M with the largely independent reduction in imbalance from G, the pure imbalance-optimizing procedure. The nearly order of magnitude difference in imbalance between MG and G does not translate to a significantly lower MSE as there is diminishing returns for imbalance optimization; the theoretical reasons for this will become clear in the following section.


\section{Methodology}\label{sec:methodology}

\subsection{Setup and Assumptions}\label{subsec:setup}

We denote the responses $\bv{y} = \bracks{y_1, \ldots, y_{2n}}^\top$ where the number of subjects $n$ is assumed even. The \emph{assignment} or \emph{allocation} vector is $\bv{w} = \bracks{w_1, \ldots, w_n}^\top$ whose entries are either +1 (the subject received $T$) or -1 (the subject received a $C$) and $\bv{w} \in \braces{-1,+1}^{2n}$. We define a \emph{design} $D$ as $W_D$, a discrete uniform random variable with support $\mathbb{W}_D \subseteq \braces{-1,+1}^{2n}$ and number of vectors $\abss{\mathbb{W}_D} \leq 2^{2n}$. We restrict the designs we consider to those that have the mirror property \citet[Assumption 2.2]{Kapelner2020}, \citet[Assumption 1b]{Kallus2018}, \citet{Nordin2020} i.e. for any $\bv{w}$, if treatment assignments and control assignments were switched then that resulting assignment would also be in the design: $\bv{w} \in \mathbb{W} \Rightarrow \bv{1} - \bv{w} \in \mathbb{W}$. We further restrict the designs we consider to be \emph{forced balance procedures} where all allocations have the same number of treated and control subjects \citep[Chapter 3.3]{Rosenberger2016}. This is a minor restriction is a standard design known as the balanced complete randomization design (BCRD) or completely random design, denoted $\mathbb{W}_{BCRD} = \braces{\bv{w}\,:\, \sum w_i = 0}$ which has $\binom{2n}{n}$ assignments. If a design further restricted this set, i.e. $\mathbb{W}_D \subset \mathbb{W}_{BCRD}$, we call $D$ a \emph{restricted design}.

We assume the only source of the randomness in the response is thus the treatment assignments $\bv{w}$. This assumption on the source of randomness is termed the \emph{randomization model} \citep[Chapter 6.3]{Rosenberger2016}, the \textit{Fisher model} or the \textit{Neyman model} whereby \qu{the $2n$ subjects are the population of interest; they are not assumed to be randomly drawn from a superpopulation} \citep[page 297]{Lin2013}.

As in \citet{Kallus2018}, we employ the simple \emph{differences-in-means} estimator, 

\vspace{-0.25cm}\begin{eqnarray}\label{eq:estimator}
\hat{\beta} := \frac{\bv{w}^\top \bv{y}}{2n} = \frac{1}{2}(\Ybar_T - \Ybar_C)
\end{eqnarray} 

which is population average treatment estimator regardless of the response model and the effect of the treatment \citep[Example 4.7]{Lehmann1998}.

\subsection{Previous Literature}

Our design idea stems from a combination of two long-standing streams of research: first, the pure-imbalance optimizing design literature and second, the binary-matching bias-reducing literature.

For the pure-imbalance optimizing design, we first reference the rerandomization design R. Here, one begins with vectors from BCRD and discards those whose covariate imbalance is beyond a desired threshold. This idea that goes back to Fisher, the father of randomized experiments who was aware of imbalanced allocations and ironically warned against pure unrestricted randomization, the procedure that he is famous for. 

Here, $\mathbb{W}_{R(a)} := \braces{\bv{w}\,:\, d_M(\bv{w}) \leq a, \bv{w} \in \mathbb{W}_{BCRD}}$ where 

\vspace{-0.25cm}\begin{eqnarray}\label{eq:mahal}
d_M(\bv{w}) := n (\bar{\bv{x}}_T - \bar{\bv{x}}_C)^\top \hat{\Sigma}_X^{-1} (\bar{\bv{x}}_T - \bar{\bv{x}}_C),
\end{eqnarray}

a popular metric of covariate imbalance among the treatment and control groups known as the Mahalanobis distance, $a$ is an upper-bound acceptability threshold of the covariate imbalance, $\bar{\bv{x}}_T$ and $\bar{\bv{x}}_C$ are the averages across the $p$ covariates in the treatment and control groups respectively and $\hat{\Sigma}_X$ is the sample variance-covariance matrix of all $2n$ subjects. Since $a < \infty$ implies that $\mathbb{W}_{R(a)} \subset \mathbb{W}_{BCRD}$, rerandomization is a restricted design by construction. Assuming multivariate-normally distributed covariates and an additive treatment effect but no assumption on the nature of the response model, the multiplicative reduction in MSE of the difference-in-averages treatment estimator when using the R design is $\eta(p, a) R^2$ where 

\vspace{-0.25cm}\begin{eqnarray}\label{eq:rerand_reduction}
\eta(p, a) := 1 - \frac{2}{p}  \frac{\gamma(p/2 + 1, a/2)}{\gamma(p/2, a/2)},
\end{eqnarray} 

$\gamma(\cdot, \cdot)$ denotes the lower incomplete gamma function and $R^2$ is variance explained in an OLS linear regression of the response on the covariates \citep{Morgan2012}.

This reduction in treatment estimator error comes from two sources (1) a reduction of the imbalance in the $p$ covariates which is at their control and (2) the degree of linearity of the covariates in the response as measured by $R^2$. 

As for their first source of estimator error reduction, lowering the imbalance in rerandomization is only slightly more impressive than BCRD. There are many ideas for achieving better reduction in covariate imbalance. First, \citet[Section 3.3]{Kallus2018} conjectured using a heuristic argument that the optimal reduction (i.e. for the best possible vector) is exponential $2^{-\Omega(n)}$. Numerical optimization such as in \citet{Bertsimas2015} employ heuristics to approximate the optimal vector. Some heuristics come with theoretical guarantees, e.g. the greedy pair-switching design of \citet{Krieger2019} (henceforth denoted KAK19) uses a greedy heuristic that begins with a draw from BCRD, considers switching every T/C pair and retains the one with the greatest reduction in the imbalance metric and stops switching when it is no longer fruitful. The resulting vectors have a provably very low imbalance $O_p(n^{-(1 + 2/p)})$ which can then be further enhanced by generating many $\bv{w}$'s and retaining only the best of those akin to rerandomization. Moreover, since there is a provably low number of switches, the degree of randomness is nearly that of BCRD. We denote this restricted design as G.

Better rates than offered by G are unlikely to matter as the imbalance reduction multiple (Equation~\ref{eq:rerand_reduction}) has diminishing returns in $a$. After a procedure such as KAK19, the performance of the estimator is much more sensitive to the second source of error, i.e. the strength of the covariates' affect on the response (as gauged by $R^2$). If the $R^2$ is low (in the case of a non-linear response model and/or a high noise setting), covariate imbalance reduction will not be too fruitful. Moreover, \citet[Section 2.3.3]{Kallus2018} proves that minimizing $a$ in rerandomization is fruitful for minimax estimator error reduction only in the case of a linear response model. \citet{Kallus2018} explains the reason why this advantage is only limited to linear response models: even if $d_M(\bv{w}) = 0$ (i.e. the first moments are matched perfectly), the covariate distributions in the two arms could be very different. In a nonlinear response model, the estimator performance can depend on more subtle differences in the covariates' distributions (e.g. their tail behavior).

The second stream of research, binary matching, implicitly attempts to equalize the covariate distributions in both arms. Here, the $2n$ subjects' indices are first organized into a \textit{matched pair structure} $\mathcal{M}_d$, a set whose elements are \textit{pairs} (a set of two subject indices). The structure $\mathcal{M}_d$ is created by first assuming a distance function $d(\bv{x}_r, \bv{x}_s) : \mathbb{R}^p \times \mathbb{R}^p \rightarrow \mathbb{R}_{\geq 0}$ whose inputs are subject covariates. A $2n \times 2n$ nonnegative symmetric matrix is then computed consisting of distances for all pairs of subjects. The optimal match structure $\mathcal{M}_d^\star$ would be one of the $\binom{2n}{n}$ sets that minimizes the sum of the resulting intramatch distances i.e. $\sum_{\braces{r,s} \in \mathcal{M}} d(\bv{x}_r, \bv{x}_s)$. The algorithm that produces the optimal match structure is known as \emph{optimal nonbipartite matching}. This algorithm has been shown to be reducible to a polynomial-time algorithm (see \citealp{Lu2011} for details and history). To create an assignment vector $\bv{w}$, the individual assignments within each of the $n$ pairs are randomly allocated to T/C or C/T via $n$ iid Bernoulli draws. We denote the optimal binary matching design with the Mahalanobis distance metric (Equation~\ref{eq:mahal}) as M and its allocation space by $\mathbb{W}_{M} = \{\bv{w}\,:\,w_r = -w_s, \braces{r,s} \in \mathcal{M}_{d_M}^\star\} \subset \mathbb{W}_{BCRD}$. 

Can binary matching play a role in experimental design? This M design was investigated by \citet{Greevy2004} who reported better imbalance and higher power using randomization tests. However, matching in order to have robust estimation under a nonlinear response model has a long literature in observational studies and we recommend \citet{Stuart2010} for a broad overview. (Without the ability to manipulate $\bv{w}$, the matching procedure is quite different as it is limited to fixed $w_i$'s and is called \emph{bipartite matching}). \citet{Rubin1979} was the first to show this robustness and recommended that $d$ be specified as the Mahalanobis distance of Equation~\ref{eq:mahal}. 

Further, \citet[Section 2.3.2]{Kallus2018} proves that binary matching via $d$ is the minimax variance minimizing experimental design for response functions that are Lipschitz continuous with respect to metric $d$. It is reasonable to conclude in a real-world experimental setting (for instance in a clinical trial in medicine) that the response would be continuous with derivatives not changing too quickly (i.e. satisfying the Lipschitz conditions). Thus, to reduce the variance of the treatment estimator in nonlinear response models, we elect to use a binary matching design.

Our contribution is relatively simple: we first compute $\mathcal{M}_{d_M}^\star$ as part of an M design and this binary pairing structure will provide robustness to nonlinear response models. Then instead of assigning the specific pairs' subjects via Bernoulli draws as explained above, we further restrict the pair assignments to insist on small imbalance by using R (resulting in the matching-then-rerandomization design we term MR) or by using G (resulting the matching-then-greedy-pair-switching design we term MG).  We turn to the details of MR and MG now.






\subsection{Our Algorithms}

First, $\mathcal{M}_{d_M}^\star$ is computed by running the optimal non-bipartite matching algorithm implemented by \citet{Beck2016} in the \texttt{R} package \texttt{nbpMatching}.

For R after M, $\mathbb{W}_{MR(a)} := \braces{\bv{w}\,:\, d_M(\bv{w}) \leq a, \bv{w} \in \mathbb{W}_M}$ i.e. we draw many different assignments from $\mathbb{W}_M$ and retain those whose Mahalanobis distances meet the threshold $a$ thereby finding small imbalances subject to the match structure.

For G after M, we cannot simply input $\bv{w} \in \mathbb{W}_M$ directly into KAK19 as this would violate the $\mathcal{M}_{d_M}^\star$ structure, breaking the binary pairings. Instead, we temporarily view the matched pairs as individual subjects (of which there are $n$) as the input into KAK19. This means we temporarily view the intramatch covariate vector differences as the subjects' covariates which are inputted into KAK19. Put another way, the T arm now refers to a matched pair assigned T/C where the covariate difference vector is computed by taking the first subject's covariates minus the second subject's covariates. The C arm would be the opposite. In one iteration of KAK19, all \textit{pair of pairs switches} are considered: from (T/C)/(C/T) to (C/T)/(T/C) or vice versa. One of these will result in the largest reduction in imbalance as measured among the pair of pairs which is algebraically the same imbalance as measured among the entire $\bv{w}$ vector for the $2n$ subjects. This best pair of pairs switch is retained and the algorithm stops when an iteration cannot reduce the imbalance any longer. This scheme thereby preserves $\mathcal{M}_{d_M}^\star$.


Under our designs, to evaluate hypotheses concerning the value of $\beta_T$ we can use the randomization test and to construct approximate confidence intervals for $\beta_T$ we can invert the randomization test. For details, see \citet[Section 2.2]{Morgan2012} and KAK19, Section 6; both provide a careful step-by-step explanation and a literature review.

One can consider an alternative implementation of the G in MG as follows. After M, each iteration of G can perform one switch of one pair (instead of a switch of a pair of pairs). Here, we now have fewer switch possibilities, i.e. $n$ and not $\binom{n}{2}$, and we conjecture this procedure would not perform as well as the one proposed above. 

One can also consider a GM (or RM) procedure i.e. to first imbalance-optimize via G (or R) and then employ binary matching. In that scenario, there will be a different $\mathcal{M}$ for each initial imbalance optimized vector based on a bipartite matching. The gain here would be that the imbalance-optimization would be slightly better but it comes at a cost: bipartite matches are much more constrained than nonbipartite matches and thus we conjecture that these designs will have much lower robustness to nonlinear models. 

\subsection{Theoretical Investigation of MG and MR}\label{sec:theory}

\subsubsection{Imbalance}

The key idea of our designs is to retain the match structure while optimizing imbalance. Instead of employing G or R using $\mathbb{W}_{BCRD}$ as a base, we now employ G or R using $\mathbb{W}_{M}$ as a base. We are not aware of theoretical results for the imbalance after optimal nonbipartite matching when $p > 1$ and thus we focus on the case when $p = 1$. 

In this univariate case, consider the order statistics $X_{(1)}, \ldots, X_{(2n)}$. Optimal matching creates $\mathcal{M} = \braces{\braces{2n, 2n-1}, \ldots, \braces{2,1}}$, i.e. pairing the first largest covariate with the second largest covariate value, then pairing the third largest with the fourth largest, etc. We will consider the case where $X_1,\ldots,X_{2n}$ are iid $\uniform{0}{1}$. Define $\Delta_i = X_{(2i)} - X_{(2i-1)}$ for $i = 1, 2, \ldots, n$. By \citet[Theorem 6.6c]{Dasgupta2011}, $\bracks{\Delta_1,\ldots,\Delta_n}$ has the same distribution as $\oneover{S}\bracks{E_1,\ldots,E_n}$, where $E_1,\ldots,E_n,\ldots,E_{2n+1}$ are iid $\exponential{1}$ and $S := \sum_{i=1}^{2n+1} E_i = O_p(n)$. Since the denominator $S$ is common to all entries of the vector, applying the greedy pair-switching of KAK19 to $\bracks{\Delta_1,\ldots,\Delta_n}$ is the same as (in terms of distribution) applying the algorithm to $\bracks{E_1,\ldots,E_n}$ and then dividing by $S$. 

The imbalance of assignments in G is $O_p(n^{-3})$ by KAK19 (Theorem 1), a distribution-free result. When dividing by $S$, which is $O_p(n)$, the resultant rate of our MG design is $O_p(n^{-4})$.


The imbalance of assignments in R is less than $a$ by construction. However, if $a$ is too small, then practical computational constraints dictate that all BCRD assignments one generates will have imbalance greater than $a$, resulting in the inability to locate any assignments in $\mathbb{W}_{R(a)}$. KAK19 reframes the discussion of R's imbalance by making this computational constraint explicit. Let $n_R$ denote the number of BCRD assignments one searches through for each resultant assignment in R. For every collection of $n_R$, the assignment that provides the minimum imbalance in that collection is retained. These retained vectors then have imbalance which is order $n_R^{-1} O_p(n^{-1/2})$ where $O_p(n^{-1/2})$ is the imbalance of BCRD. It then follows by the same argument above that when dividing by $S$, the resultant imbalance of MR is then $n_R^{-1} O_p(n^{-3/2})$. 

Extending the theoretical results for imbalance in MG beyond the assumption of uniformly distributed covariates is difficult. The main result needed to prove that the order of the imbalance of the assignments in G is $O_p(n^{-3})$ is to show that there exists a pair of covariate values to switch $X_i, X_j$ where $w_i = 1$ and $w_j = -1$ that result in a necessary reduction in imbalance. The analogous result here would require the locating of $\Delta_i, \Delta_j$ where $w_i = 1$, $w_{i+1} = -1$, $w_j = -1$ and $w_{j+1} = 1$ that results in the same reduction in imbalance after switching both pairs. Showing the necessary reduction in imbalance assumed that the $X_i$'s are continous and iid. We were able to use this result in this paper when the $X_i$'s were uniformly distributed, but if the $X_i$'s are normally distributed, it is much more complicated as the $\Delta_i$'s corresponding to covariate values in a tail are stochastically larger than the $\Delta_i$'s corresponding to covariate values in the center. The imbalance of MG in the normal case is shown by simulation to be lower than G (see Figure~\ref{fig:mse_density_plots_normal} in the Supporting Information).





\subsubsection{Degree of Randomness}

We now turn to degree of randomness in MG and MR. We reiterate that in experimental design, it is important for the assignments to be highly random as this randomness provides insurance against a large variance in the treatment effect estimator due to unobserved covariates \citep{Kapelner2020}. Both G and R are less random than BCRD since they are forms of restricted randomization. In the case of G, KAK19 (Theorem 2 and Proposition 1) shows that G's assignments are asymptotically as random as BCRD according to three degree-of-randomness metrics: (a) the pairwise entropy metric of assignments, (b) the standard error metric of probability of pairwise assignments and (c) the accidental bias metric of \citet[Section 5]{Efron1971} defined as the maximum eigenvalue of the variance covariance matrix of $W$, the random variable producing assignments in any strategy (see Section~\ref{subsec:setup}). The reason for G's high degree of randomness is that the number of pairwise switches made in G is small, being only $O_p(\sqrt{n})$. Thus, since G starts with BCRD, and its assignments are not changed significantly, the resulting $\mathbb{W}_G$ is very large. In the case of R, we directly limit the number of designs based on the cutoff $a$. As, $a$ decreases, $\mathbb{W}_R$ shrinks relative to $\mathbb{W}_{BCRD}$, making R less random. Since $a$ can never be too small due to computational limitations, it follows that randomness is not substantially deteriorated in R. Thus, MG is as random as M asymptotically and MR-in-practice is likely as random as M asymptotically as well. And M is highly random according to all three degree-of-randomness metrics, only slightly less random than BCRD. 

\section{Simulation Results}\label{sec:simulations}

To illustrate our main results we turn to simulations. All simulations herein were performed with \texttt{GreedyExperimentalDesign}, an \texttt{R} package available on \texttt{CRAN} whose core is implemented in \texttt{Java} for speed. We measure mean squared error (MSE) of the difference-in-averages estimator for an additive average treatment effect when the response model depends on underlying covariate(s). 

We conjecture that (a) the MSE performance of MR/MG will be similar to imbalance-optimizing designs G/R when the response model is purely linear, (b) the MSE performance of MR/MG will be similar to the pure nonbipartite matching design M when the response model is purely nonlinear and (c) the MSE performance of MR/MG will be better than either exclusively imbalance-optimizing G/R or exclusively nonbipartite matching M when the response model is a hybrid of a strong linear component and strong nonlinear component.

To demonstrate the conjecture we consider $n=100$, two covariates and five response models that take the form $\bv{Y} = \beta_T\bv{w} + f(x_1, x_2) + \bv{\errorrv}$ where $x_1, x_2$ are independent draws from $U(-\sqrt{3},+\sqrt{3})$ so that their variance is one and the $n$ components of $\bv{\errorrv}$ are iid draws from $\normnot{0}{0.5^2}$. The five models are (Z) a null zero model $f = 0$ for calibration and ensuring the simulation is working (L) a purely linear model $f = 3x_1 + 3x_2$, (LsNL) a mostly linear but slightly nonlinear model $f = 3x_1 + 3x_2 + x_1^2$, (LNL) a linear and nonlinear model $f = 3x_1 + 3x_2 + x_1^2 + x_2^2 + x_1 x_2$ and (NL) a purely nonlinear model $f = x_1^2 + x_2^2 + x_1 x_2$. 

We then consider the six designs of the previous section all of which feature an equiprobable draw from a set of mirrored allocations: (BCRD) balanced and completely randomized, (R) rerandomization of \citet{Morgan2012} where only the top 1\% of BCRD vectors are retained, (G) the greedy pair switching of KAK19, (M) optimal nonbipartite matching, and our two new hybrid approaches of (MG) optimal nonbipartite matching followed by greedy pair switching on the within-pair covariate differences and (MR) optimal nonbipartite matching followed by rerandomization on the within-pair covariate differences where the top 1\% of pairings are retained. Thus, we have $5 \times 6$ = 30 sample size-model-design simulation settings.

For each sample size-model-design simulation setting, we simulate 50 different covariate and noise value settings, $x_1, x_2, \mathcal{E}_1, \ldots, \mathcal{E}_n$. Within each specific covariate-noise realization setting, we draw 500 unique allocation vectors from each of the six design strategies. We define MSE as $\expesubnostr{X_1, X_2, \bv{\errorrv}}{\cexpesubnostr{W_D}{(\hat{\beta}_T - \beta_T)^2}{X_1, X_2}}$ for each sample size-model-design setting. This MSE is estimated by first taking an arithmetic average of the 500 squared errors within the covariate-noise setting and then averaging these over the 50 covariate-noise settings. Estimates for matching designs (M and MR) usually are computed by the average paired differences. This estimate is algebraically equivalent to the naive difference-in-averages estimate and thus estimates are computed uniformly for all six designs.


\begin{figure}[h]
\centering
\includegraphics[width=6.5in]{fig1.pdf}
\caption{Density plots for MSE (on a log scale) over the covariate distributions, noise and assignments when $n=100$ for the five models (the vertically stacked panes) and six designs (depicted as different colors) explained in the text. The vertical lines indicate the MSE estimate (explained in the text) by design. Lines directly adjacent to one another are the result of a slight jitter. Statistically significant differences are tabulated in Table~\ref{tab:all_pvals}.}
\label{fig:mse_density_plots}
\end{figure}

Figure~\ref{fig:mse_density_plots} provides density plots for MSE when $n=100$ for all five models and six designs and Table~\ref{tab:all_pvals} tabulates the statistically significance of each design comparison. The null response model (Z) serves as a check on the integrity of our simulation: since the covariates are disjointed from the response, all designs perform equally as seen in the top pane of Figure~\ref{fig:mse_density_plots} and for this reason it was omitted from Table~\ref{tab:all_pvals}. For all other response models where the covariates inform the response, the figure and tables illustrates the main thrust of this paper's contribution: (1) Our designs, MR and MG, are always the best performers (or tied for the best performance) indicating robustness to linear, nonlinear and mixtures of linear and nonlinear response models. Other observations support this main thrust. (2) BCRD is the worst performer (or tied for worst performance) which is expected given that it does not allocate subjects with the benefit of knowing the covariate values. (3) In the exclusively linear model (L), the pure imbalance optimizing designs perform the best i.e. G, R, MR, MG all perform equally with M lagging. (4) In the exclusively nonlinear model (NL), the matching designs perform the best i.e. M, MR, MG all perform equally. (5) In the linear-nonlinear models (LsNL, LNL), our matching-plus-imbalance-optimizing designs MR and MG are equally the best, the matching M is second best and the pure imbalance optimizing designs G and R lag behind M. (6) R lags behind G because its imbalance is a much smaller rate (see KAK19 for details why). MR also lags behind MG for the same reason but the difference is not appreciable in these simulations.

\begin{table}[h]
\centering

\begin{subtable}{.5\linewidth}{
\begin{tabular}{rlllll}
 & R & G & M & MR & MG \\ 
  \hline
BCRD & *** & *** & *** & *** & *** \\ 
  R &  &  &  &  &  \\ 
  G &  &  & ** &  &  \\ 
  M &  &  &  & ** & ** \\ 
  MR &  &  &  &  &  \\ 
\end{tabular}
}
\caption{Model: L}
\label{tab:pvals1}
\end{subtable}~
\begin{subtable}{.5\linewidth}{
\begin{tabular}{rlllll}
 & R & G & M & MR & MG \\ 
  \hline
BCRD & *** & *** & *** & *** & *** \\ 
  R &  &  & *** & *** & *** \\ 
  G &  &  & ** & *** & *** \\ 
  M &  &  &  & * & * \\ 
  MR &  &  &  &  &  \\ 
\end{tabular}
}
\caption{Model: LsNL}
\label{tab:pvals2}
\end{subtable}\\
~\\~\\

\begin{subtable}{.5\linewidth}{
\begin{tabular}{rlllll}
 & R & G & M & MR & MG \\ 
  \hline
BCRD & *** & *** & *** & *** & *** \\ 
  R &  &  & *** & *** & *** \\ 
  G &  &  & *** & *** & *** \\ 
  M &  &  &  & * & * \\ 
  MR &  &  &  &  &  \\ 
\end{tabular}
}
\caption{Model: LNL}
\label{tab:pvals3}
\end{subtable}~
\begin{subtable}{.5\linewidth}{
\begin{tabular}{rlllll}
 & R & G & M & MR & MG \\ 
  \hline
BCRD &  & *** & *** & *** & *** \\ 
  R &  & *** & *** & *** & *** \\ 
  G &  &  & *** & *** & *** \\ 
  M &  &  &  &  &  \\ 
  MR &  &  &  &  &  \\ 
\end{tabular}
}
\caption{Model: NL}
\label{tab:pvals4}
\end{subtable}

\caption{Tukey-Kramer comparisons of the MSE for each design pair by model. * indicates $p < 0.05$, ** indicates $p < 0.01$ and *** indicates $p < 0.001$. The Z model is not displayed as none of its design pairs' MSEs tested statistically significant.}
\label{tab:all_pvals}
\end{table}


We also repeated the simulation for three other sample sizes ($n=32$, $n=132$ and $n=200$) and they all have qualitatively similar results. The tabulation of MSE estimates for all sample sizes, designs and response models can be found in Tables~\ref{fig:ols_n_32}-\ref{fig:ols_n_200} in the Supporting Information. The entire simulation was also repeated for the case of normally distributed covariates and the results were not substantially different (see Figure~\ref{fig:mse_density_plots_normal} in the Supporting Information).

\section{Concluding Remarks}\label{sec:conclusion}



We proposed a randomized experimental design that is a hybrid between two well studied strategies: those that optimize covariate imbalance and those that use binary matching. The former is known to be minimax optimal when the response model is linear in the covariates and the latter is known to be minimax optimal when the response model is continuous with limits on its derivative with respect to all $p$ covariates. Additionally, we have theory that shows that the fusing of both strategies fortuitously enhances the covariate imbalance all while retaining a degree of randomization close to that of the classic completely random designs. 

Because we fuse both types of design together without sacrificing the performance of either, our designs provides very low MSE performance when estimating a population average treatment effect in the purely linear case (because we optimize covariate imbalance), the purely nonlinear case (as we our assignments are optimal nonbipartite matched pairs) and response models that feature both linear and nonlinear components. Thus we expect our design to be very powerful when employed in real-world settings such as clinical trials and Internet-based experimentation.

There are many extensions to this work. First, there are different types of matching such as ratio matching \citep[Section 3.1.2]{Stuart2010}, matching in more than two arms (ibid, Section 6.1.4), matching with discrete covariates (the Mahalanobis distance function is not the most appropriate here as discussed in \citealt{Gu1993}) and matching with some covariates being more important than others. As some matches may not be acceptable, one can caliper match; this would result in a subset of the subjects being matched and a subset being unmatched. In this case, the G algorithm would have to be redesigned but the R procedure would be fairly straightforward. Further, one may wish to use the OLS estimator instead of the classical difference-in-averages estimator (Equation~\ref{eq:estimator}). One can argue that restricted randomization can leave you vulnerable to high MSE estimation due to imbalance in unobserved covariates \citep[Sections 2.2.1, 2.2.3, 2.3.1 and 2.3.3]{Kapelner2020}. Therein, they proved that restricted designs with a high degree of randomness, such as our hybrid designs in this work (see Section~\ref{sec:theory}), are highly random and thus are not susceptible. Further, \citet[Equation 14]{Kapelner2020} showed that the OLS estimator improves asymptotic MSE by an order of magnitude in sample size. However, imbalance-optimizing designs (such as R and G) still provide improvement in the linear case. Covariate adjustment cannot aid with the nonlinear component of the response model and thus our hybrid designs will perform very well in conjunction with the OLS estimator on nonlinear models. And the OLS estimator is actually preferred in matching designs \citep{Rubin1979}. Future work can also adapt this design to the sequential setting, where individuals in the experiment arrive one by one and must be assigned quickly to an arm. We can apply the on-the-fly matching procedure in \citet{Kapelner2014} with the rerandomization approach in \citet{Zhou2018} in the same vein as our hybrid procedure. 

\subsection*{Replication}

All figures and tables can be reproduced by running the \texttt{R} code found at \url{https://github.com/kapelner/GreedyExperimentalDesign/blob/master/hybrid_paper}.

\subsection*{Acknowledgements}

This research was supported by Grant No 2018112 from the United States-Israel Binational Science Foundation (BSF).

\bibliographystyle{apa}
