Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
88,985 characters · 16 sections · 70 citation commands
Optimality of Matched-Pair Designs in Randomized Controlled Trials
}
Keywords: Matched-pair design, baseline outcome, stratified randomization, experiment, randomized controlled trial
JEL classification codes: C12, C13, C14, C90
\thispagestyle{empty} \setcounter{page}{1}
This paper studies the optimality of matched-pair designs in randomized controlled trials (RCTs). Matched-pair designs are examples of stratified randomization, in which the researcher partitions a set of units into strata (groups) based on their observed covariates and assigns a fraction of units in each stratum to treatment. A matched-pair design is a stratified randomization scheme with two units in each stratum.
Stratified randomization is prevalent in economics. Among the 5,000 RCTs in the AEA RCT Registry, more than 800 are stratified. The schemes in these papers, however, differ vastly in terms of the covariates used to stratify and how fine the strata are. Among these 800 RCTs, around 50 use matched-pair designs. Moreover, 56% of the researchers interviewed in bruhn2009pursuit have used matched-pair designs at some point in their research. Yet, despite the frequency with which applied researchers make decisions about how to stratify, there are few general econometric results on whether matched-pair designs lead to better precision of estimators of treatment effects than other stratified randomization schemes and the best way to pair units.
I derive the exact form of the stratified randomization scheme that has the maximum statistical precision for estimating the average treatment effect (ATE). The optimal scheme is a matched-pair design. In an important special case, the optimal design is to order the units according to the baseline values of the primary outcome variable of interest and then pair the adjacent units. When I simulate this simple design using data from 10 recent papers in the American Economic Journal: Applied Economics, I find it lowers the standard error of the difference-in-means estimator by 10% on average, and by up to 34%, relative to the designs actually used in those studies. I also find some more complicated stratifications with strata of four units according to multiple covariates could further lower both the MSE and the standard error. Based on these findings, I make practical recommendations across a wide range of empirical settings.
In Section (ref), I study settings where the treated fractions are identically $\frac{1}{2}$ across strata. In such settings, a common estimator for the ATE is the difference in the means of the treated and control groups. The properties of the difference-in-means estimator, however, vary substantially with how the researcher stratifies. To begin, consider the thought experiment where we know the distributions of the potential outcomes. Let $Y(1)$ denote the potential outcome if a unit is treated and let $Y(0)$ denote the potential outcome if it is not treated. Let $X$ denote the observed, baseline covariates. I define an index function $E[Y(1) + Y(0) | X]$, the expected sum of the potential outcomes given the covariates. My first result shows the mean-squared error (MSE) of the difference-in-means estimator is minimized by a matched-pair design, where units are ordered according to this index function and paired adjacently. My optimality result holds at any sample size and without any distributional assumption beyond the existence of moments. In particular, my result does not rely on restrictions on treatment effects heterogeneity.
I describe a special case where the optimal stratification is feasible even without knowing the index function. Suppose $X$ contains a single covariate. Further suppose both $Y(1)$ and $Y(0)$ are higher in expectation when $X$ is higher, so that $E[Y(1) + Y(0) | X]$ is increasing in $X$. In this case, pairing units according to $X$ is optimal. An important example in empirical practice is when $X$ is the baseline value of the primary outcome variable of interest. For instance, in angrist2009effects, the primary outcome variable of interest is a test score and the treatment is an educational program, so we expect a higher baseline test score ($X$) implies a higher endline test score ($Y(1)$ and $Y(0)$) in expectation.
If researchers are unsure about the monotonicity condition, or if multiple covariates are available, then the optimal stratification is generally unknown because the index function is generally unknown. As such, Section (ref) studies several feasible procedures. With multiple covariates, I study pairing units to minimize the (Mahalanobis) distances of the covariates. In settings with auxiliary data, such as data from pilot studies, I propose several matched-pair designs in which the index function is approximated by a proxy based on the auxiliary data.
In Section (ref), to compare the performance of these practical procedures, I study the asymptotic properties of the difference-in-means estimator. I show that relative to not stratifying, pairing according to any function of the covariates can only reduce the limiting variance of the difference-in-means estimator. Moreover, the limiting variance is lower if the stratifying variables explain a larger proportion of the variation in $Y(1) + Y(0)$.
In Section (ref), I conduct a simulation study using data from a systematically selected set of 10 RCTs from recent issues of the American Economic Journal: Applied Economics. Relative to the original stratifications used in those 10 papers, if the researchers had just paired the units according to their baseline outcomes, then the MSE of the difference-in-means estimator would be 24% smaller on average and 56% smaller in some cases. The standard error of the difference-in-means estimator would be 10% smaller on average and 34% smaller in some cases.
Among all methods in the simulation, pairing units to minimize the sum of the squared Mahalanobis distances of the covariates usually leads to the smallest MSEs. When the number of covariates is large, however, the standard error could be even larger than the standard error when pairing according to the baseline outcome alone. Intuitively, this is because the quality of the variance estimator is lower when the curse of dimensionality is more severe. An alternative that balances the MSE and the standard error is to match units into sets of four, instead of pairs, to minimize the sum of the squared Mahalanobis distance of the covariates. Such a method has both smaller MSEs and standard errors than pairing according to the baseline outcome alone while being computationally more intensive.
I conclude with recommendations for empirical practice in Section (ref). I recommend different stratifications based on the availability of auxiliary datasets and whether one main outcome of interest clearly dominates the others. All of my recommended procedures are defined by pairing units or matching units into sets of four according to all or a subset of the available covariates.
\paragraph{Related Literature} This paper is most closely related to barrios2013optimal and tabord-meehan2020stratification. barrios2013optimal studies minimizing the variance of the difference-in-means estimator. He is the first to show pairing units according to my index function is optimal among all matched-pair designs, albeit under the assumption of homogeneous treatment effects. My optimality result holds among all stratified randomization schemes and with heterogeneous treatment effects. tabord-meehan2020stratification studies optimality within a class of stratification trees. Because the number of strata is fixed in his asymptotic framework, he can optimize over the treated fraction in each stratum. In a matched-pair design, the number of strata is half of the sample size and hence not fixed as the sample size increases, so matched-pair designs are precluded in his framework. In Section (ref) of the supplement, I elaborate on the comparison between the two papers and further note that combining our procedures is straightforward.
The following papers also study matched-pair designs: greevy2004optimal study pairing units to minimize the sum of the squared Mahalanobis distances of the covariates. imai2008variance studies matched-pair designs, focusing on the sample ATE. The inference methods in this paper build on and extend those in bai2021inference. In addition, inference under matched-pair designs has also been studied in abadie2008estimation, who assume a different sampling framework, fogarty2018mitigating,fogarty2018regression-assisted, who provides conservative estimators for the limiting variance, and de_chaisemartin2021at in a finite-population setting.
Let $Y_i$ denote the observed outcome of interest for the $i$th unit, let $D_i$ denote the treatment status for the $i$th unit, and let $X_i$ denote the observed, baseline covariates for the $i$th unit. Further denote by $Y_i(1)$ the potential outcome of the $i$th unit if treated and by $Y_i(0)$ if not treated. As usual, the observed outcome is related to the potential outcomes and treatment status by the relationship \[ Y_i = Y_i(1) D_i + Y_i(0) (1 - D_i)~. \] For ease of exposition, I assume the sample size is even and denote it by $2n$. I assume $((Y_i(1), Y_i(0), X_i): 1 \leq i \leq 2n)$ is an i.i.d. sequence of random vectors. Note the potential outcomes and the covariates are drawn from a population and hence are random instead of fixed. For any random vector indexed by $i$, $A_i$, define $A^{(n)} = (A_1, \dots, A_{2n})'$. The main parameter of interest is the average treatment effect (ATE): \[ \theta = E[Y_i(1) - Y_i(0)]~. \]
In stratified randomization, I first partition the set of units into strata. Formally, I define a stratification $\lambda = \{\lambda_s: 1 \leq s \leq S\}$ as a partition of $\{1, \dots, 2n\}$:
Let $\Lambda_n$ denote the set of all stratifications of $2n$ units. Define $n_s = |\lambda_s|$ and $\tau_s$ as the treated fraction in stratum $\lambda_s$. A matched-pair design is simply a stratified randomization scheme with $S = n$ and $n_s = 2$ for $1 \leq s \leq S$. I define $\Lambda_n^{\rm pair} \subseteq \Lambda_n$ as the set of all matched-pair designs for $2n$ units.
I make the following assumption on the treatment assignment scheme:
Assumption (ref) implies
In other words, treatment status and potential outcomes are conditionally independent given the covariates. Assumption (ref) also implies $n_s$ has to be even because a unit cannot be cut in half. Note the distribution of the vector of treatment status $D^{(n)}$ depends on $\lambda$. Most results below can be extended to settings where $\tau_s, 1 \leq s \leq S$ are identical but not $\frac{1}{2}$, or where they are additionally allowed to vary across subpopulations. See Remark (ref) for details.
For all treatment assignment schemes in the main text, I estimate the ATE by the difference in the means of the treated and control groups. Formally, for $d \in \{0, 1\}$, define \[ \hat \mu_n(d) = \frac{1}{n} \sum_{1 \leq i \leq 2n: D_i = d} Y_i~. \] The difference-in-means estimator is defined as \[ \hat \theta_n = \hat \mu_n(1) - \hat \mu_n(0)~. \] The difference-in-means estimator is widely used because it is simple and transparent. Under Assumption (ref), it coincides with the OLS estimator for the coefficient in the linear regression of the outcome on treatment status and strata fixed effects and the OLS estimator from the fully saturated version of that regression, both of which are also widely used in analyses of RCTs. See, for example, duflo2007using, glennerster2013running, and crepon2015estimating.
This section studies the optimal stratification. To preview the results, define the index function
I show the optimal stratification is given by ordering the units according to $g_i = g(X_i)$ and then pairing the adjacent units. In the special case where $X_i$ is a scalar and $E[Y_i(1) | X_i = x]$ and $E[Y_i(0) | X_i = x]$ are both weakly increasing (or both weakly decreasing) in $x$, the optimal stratification is given by ordering the units according to $X_i$ and then pairing the adjacent units.
The analysis in this section is conditional on $X^{(n)}$. In this section only, instead of the population ATE, I focus on the ATE conditional on $X^{(n)}$: \[ \theta_n = \frac{1}{2n} \sum_{1 \leq i \leq 2n} E[Y_i(1) - Y_i(0) | X_i]~. \] Focusing on $\theta_n$ simplifies the discussion. Moreover, conditional on a fixed sample with covariates $X^{(n)}$, I can only hope to be unbiased for $\theta_n$ instead of $\theta$. The conclusions of the theorems in this section are the same regardless of whether the parameter of interest is $\theta_n$ or $\theta$.
My objective function is the MSE of $\hat \theta_n$ for $\theta_n$ conditional on $X^{(n)}$ under a stratification $\lambda \in \Lambda_n$: \[ \operatorname{MSE} (\lambda | X^{(n)}) = E_\lambda[(\hat \theta_n - \theta_n)^2 | X^{(n)}]~. \] Here, the notation $E_\lambda$ indicates the distribution of the vector of treatment status $D^{(n)}$ depends on the stratification. I consider minimizing the conditional MSE over the set of all stratifications:
In what follows, I derive the optimal stratification as the solution to (ref). I emphasize that by a simple bias-variance decomposition, one can show (ref) is equivalent to the problem where $\theta_n$ is replaced by $\theta$, so focusing on $\theta_n$ in this section is genuinely without loss of generality.
Solving (ref) involves two intermediate results, each carrying additional insights into the problem. To describe the first intermediate result, I define the ex-ante bias of $\hat \theta_n$ for $\theta_n$ conditional on $X^{(n)}$ as \[ \operatorname{Bias}_{n, \lambda}^{\rm ante}(\hat \theta_n | X^{(n)}) = E_\lambda[\hat \theta_n | X^{(n)}] - \theta_n~, \] and the ex-post bias of $\hat \theta_n$ for $\theta_n$ conditional on $X^{(n)}$ and $D^{(n)}$ as \[ \operatorname{Bias}_n^{\rm post}(\hat \theta_n | X^{(n)}, D^{(n)}) = E[\hat \theta_n | X^{(n)}, D^{(n)}] - \theta_n~. \] Here, ex-ante bias refers to the bias conditional only on the covariates, before treatment status is realized; ex-post bias refers to the bias conditional on both the covariates and treatment status, after treatment status is realized. Note in the definition of the ex-post bias, the $\lambda$ subscript does not appear because $D^{(n)}$ is already given. Note from the definition of the difference-in-means estimator that
By Assumption (ref), the marginal treatment probability of each unit satisfies $E_\lambda[D_i | X^{(n)}] = \frac{1}{2}$, and together with the conditional independence assumption in (ref), they imply \[ E_\lambda[\hat \theta_n | X^{(n)}] = \theta_n~. \] Therefore, the ex-ante bias is identically zero across $\lambda \in \Lambda_n$, which is not surprising because the ex-ante bias should be zero if we run an experiment. By the law of iterated expectations, \[ E_\lambda[\operatorname{Bias}_n^{\rm post}(\hat \theta_n | X^{(n)}, D^{(n)}) | X^{(n)}] = \operatorname{Bias}_{n, \lambda}^{\rm ante}(\hat \theta_n | X^{(n)}) = 0~, \] so the mean of the ex-post bias over the distribution of treatment status equals the ex-ante bias, which is zero.
The first intermediate result is a decomposition of the conditional MSE in (ref). Because $E_\lambda[\hat \theta_n - \theta_n | X^{(n)}] = 0$, by the law of total variance,
For any $\lambda \in \Lambda_n$, the first term on the right-hand side of (ref) equals
which is identical across all $\lambda \in \Lambda_n$. Note I used the conditional independence assumption in (ref), the facts that $\theta_n$ is a constant given $X^{(n)}$, that $D_i(1 - D_i) = 0$ for $1 \leq i \leq 2n$, and that $E_\lambda[D_i | X^{(n)}] = \frac{1}{2}$. Hence, (ref) is further equivalent to minimizing the second term on the right-hand side of (ref), which is the variance of the ex-post bias: \[ \operatorname{Var}_\lambda[\operatorname{Bias}_n^{\rm post}(\hat \theta_n | X^{(n)}, D^{(n)}) | X^{(n)}]~. \] The discussion so far leads to my first intermediate result:
Next, I describe the second intermediate result in solving (ref). The result states any stratification is a convex combination of matched-pair designs. Formally, for $\lambda, \lambda' \in \Lambda_n^{\rm pair}$ and $\delta \in [0, 1]$, define $\delta \lambda \oplus (1 - \delta) \lambda'$ as the randomization between $\lambda$ and $\lambda'$ such that $\lambda$ is implemented with probability $\delta$. Define the convex hull formed by all convex combinations of any finite number of matched-pair designs as
In other words, a member of the convex hull is the “mixing” of $J$ matched-pair designs, where $J$ is finite.
For example, suppose $2n = 4$. Then, four stratifications are possible:
$\lambda^0$ puts four units in the same stratum. $\lambda^1$ pairs 1 and 2 together and 3 and 4 together. $\lambda^2$ and $\lambda^3$ are defined similarly. Note that implementing each of the three matched-pair designs with probability $1 / 3$ is equivalent to implementing $\lambda^0$, in the sense that the distributions of $(D_1, D_2, D_3, D_4)$ are the same under the two implementations. Indeed, under $\lambda^1$, $(D_1, D_2, D_3, D_4)$ takes the following four values each with probability $1/4$: $(1, 0, 1, 0)$, $(1, 0, 0, 1)$, $(0, 1, 1, 0)$, $(0, 1, 0, 1)$. Similarly, under $\lambda^2$, it takes the following four values each with probability $1/4$: $(1, 0, 0, 1)$, $(1, 1, 0, 0)$, $(0, 0, 1, 1)$, $(0, 1, 1, 0)$. Under $\lambda^3$, it takes the following four values each with probability $1/4$: $(1, 0, 1, 0)$, $(1, 1, 0, 0)$, $(0, 1, 0, 1)$, $(0, 0, 1, 1)$. Accordingly, under $\frac{1}{3} \lambda^1 \oplus \frac{1}{3} \lambda^2 \oplus \frac{1}{3} \lambda^3$, it takes the following six values each with probability $1/6$: $(1, 1, 0, 0)$, $(1, 0, 1, 0)$, $(1, 0, 0, 1)$, $(0, 1, 1, 0)$, $(0, 1, 0, 1)$, $(0, 0, 1, 1)$. This distribution is the same as that of $(D_1, D_2, D_3, D_4)$ under $\lambda^0$, where two out of four units are treated uniformly at random. As a result, $\lambda^0 \in \mathrm{co}(\{\lambda^1, \lambda^2, \lambda^3\})$, meaning $\lambda^0$ can be written as a convex combination of the three matched-pair designs.
I show in Section (ref) of the supplement that the result above holds in general and summarize it into the following lemma:
Combining Lemmas (ref)--(ref) to minimize the MSE as in (ref) is now straightforward. To state the result, I need an equivalent notation for matched-pair designs. Recall that a permutation of $\{1, \dots, 2n\}$ is a function that maps $\{1, \dots, 2n\}$ onto itself. Let $\Pi_n$ denote the group of all permutations of $\{1, \dots, 2n\}$. A matched-pair design is a stratified randomization scheme with \[ \lambda = \{\{\pi(2s - 1), \pi(2s)\}: 1 \leq s \leq n\}~, \] where $\pi \in \Pi_n$. Recall the definition of the index function $g$ in (ref) and order the units by defining $\pi^g \in \Pi_n$ that satisfies $g_{\pi^g(1)} \leq \dots \leq g_{\pi^g(2n)}$. Define the stratification
The stratification in (ref) is given by ordering the units according to $g_i$ and then pairing the adjacent units. I now show it minimizes the MSE as in (ref).
For each $\lambda \in \Lambda_n$, define $V(\lambda)$ as the objective in Lemma (ref). Recall $g^{(n)} = (g_1, \ldots, g_n)'$. Then,
where the first equality follows from the definition of the ex-post bias, the second equality follows from (ref) and the fact that $\theta_n$ is a constant given $X^{(n)}$, and the last two equalities follow by inspection. Recall the variance of $D_i$ is $\frac{1}{4}$. Also recall that for a matched-pair design, the covariance between treatment status of the two units in a pair is $- \frac{1}{4}$, and that of units across pairs is 0. Therefore, for any $\lambda = \{\{\pi(1), \pi(2)\}, \allowbreak \dots, \{\pi(2n - 1), \pi(2n)\}\} \in \Lambda_n^{\rm pair}$, \[ V(\lambda) = \frac{1}{4n^2}\sum_{1 \leq s \leq n} (g_{\pi(2s-1)} - g_{\pi(2s)})^2~. \] Therefore, $V(\lambda)$ is proportional to the sum of squared distances of $g$ within each pair. By Lemma (ref) in the supplement, which is a simple consequence of the Hardy-Littlewood-P\'olya rearrangement inequality, $V(\lambda^g(X^{(n)})) \leq V(\lambda)$ for any $\lambda \in \Lambda_n^{\rm pair}$. Therefore, $\lambda^g(X^{(n)})$ minimizes the MSE among $\Lambda_n^{\rm pair}$, the set of all matched-pair designs.
To conclude $\lambda^g(X^{(n)})$ is optimal among the set of all stratifications $\Lambda_n$, note each stratification is a mixing of matched-pair designs, and no “mixed strategy” has a better payoff than the optimal “pure strategy.” Formally, by Lemma (ref), any $\lambda \in \Lambda_n$ can be written as \[ \lambda = \bigoplus_{1 \leq j \leq J} \delta_j \lambda^j~, \] where $\lambda^j \in \Lambda_n^{\rm pair}$, $\delta_j \geq 0$ for $1 \leq j \leq J$, and $\sum_{1 \leq j \leq J} \delta_j = 1$. As a result,
where the equality follows from the definition of the MSE, the first inequality follows because any weighted average of a set of numbers is weakly larger than the minimum across them, and the last inequality follows because $\lambda^g(X^{(n)})$ minimizes $\operatorname{MSE}(\lambda | X^{(n)})$ across $\Lambda_n^{\rm pair}$. Therefore, I have established my main theorem on the optimal stratification:
Note from (ref) that the index function $g_i$ is a scalar regardless of the dimension of $X_i$. Moreover, the optimal stratification depends not on the values but merely on the ordering of $g_i$. For instance, if $X_i$ is univariate and $g(x)$ is monotonic in $x$, then the optimal stratification in (ref) is given by ordering the units by $X_i$ and then pairing the adjacent units. This scenario arises in many settings, especially if $X_i$ is the baseline value of the primary outcome variable of interest, which is collected in the baseline survey before treatment is assigned. For instance, angrist2009effects study the effect of an educational program on test scores. In their paper, $X_i$ is the baseline test score, so we expect both $E[Y_i(1) | X_i]$ and $E[Y_i(0) | X_i]$ are weakly increasing in $X_i$. I record this result as a theorem. Let $\pi^X \in \Pi_n$ be such that $X_{\pi^X(1)} \leq \dots \leq X_{\pi^X(2n)}$.
The optimal stratification in Theorem (ref) depends on the index function $g$, which is generally unknown, so the optimal stratification is also generally unknown. Therefore, researchers often need to approximate the index function with some proxies, possibly with the help of auxiliary data. This section studies a wide range of feasible stratification methods. Some procedures are based on data from pilot experiments, which are smaller-scale copies of the main experiment run on the same population. Depending on the availability of a pilot experiment and its sample size, different procedures are available. I switch the parameter of interest back to the population ATE $\theta$, recalling that all results in the previous section hold for both $\theta_n$ and $\theta$.
According to Theorem (ref), if $X_i$ is univariate and the index function $g(x)$ is monotonic in $x$, then the optimal stratification is given by pairing units according to $X_i$. A prominent example is where $X_i$ is the baseline value of the primary outcome variable of interest, and $E[Y_i(1) | X_i = x]$ and $E[Y_i(0) | X_i = x]$ are both weakly increasing or both weakly decreasing in $x$.
Even if the monotonicity condition fails, units can still be paired according to their baseline outcomes. Theorem (ref) and Remark (ref) below study the limiting variance of the difference-in-means estimator. They reveal that if we need to choose a single covariate to pair on, the smallest limiting variance is attained by pairing units according to a covariate that explains the largest proportion of the variation in the potential outcomes. bruhn2009pursuit note the baseline outcome is often such a covariate. Simulation evidence in Section (ref) further shows pairing units according to the baseline outcome performs better than the status-quo methods in terms of both the MSE and the standard error of the difference-in-means estimator.
Regardless of whether the baseline outcome is available, if $X_i$ is multivariate, then researchers can also pair units to minimize the sum of the squared Mahalanobis distances of the covariates:
Here, $\hat \Sigma_n$ is the sample variance matrix of $X$. (ref) is simply the squared Euclidean distance if $\hat \Sigma_n$ is the identity matrix, and $\hat \Sigma_n^{-1}$ serves as a scale normalization because different covariates may be measured in different units or have different standard deviations. Note \[ d(x_1, x_2) = \|\hat \Sigma_n^{-1/2} (x_1 - x_2)\|^2~, \] where $\hat \Sigma_n^{-1/2}$ is the square root of $\hat \Sigma_n$. So the Mahalanobis distance between $x_1$ and $x_2$ equals the Euclidean distance between $\hat \Sigma_n^{-1/2} x_1$ and $\hat \Sigma_n^{-1/2} x_2$.
When the baseline outcome is unavailable but a large amount of auxiliary data is available, I can calculate a sample counterpart of (ref). In general, the auxiliary data needs to come from pilot experiments, but in one special case, even observational data suffices. If the conditional ATEs are homogeneous, meaning
then the ordering of $g_i$ is the same as that of $E[Y_i(0) | X_i]$. Suppose we have an observational dataset where the distribution of $(Y_i(0), X_i)$ is the same as that in the main experiment. As an example, suppose in an RCT to study the effect of educational program on test scores, the researcher has administrative data on the test scores of the previous cohort, and the distributions of $(Y_i(0), X_i)$ are the same across the two cohorts. Then, they can estimate $E[Y_i(0) | X_i = x]$ by a nonparametric regression using the data for the previous cohort and pair the units in the current cohort according to the predicted values in the regression. A key requirement is that the estimator for $E[Y_i(0) | X_i = x]$ is consistent in the sense of Assumption (ref) and Theorem (ref) below. Then, as the sample sizes of the auxiliary data and the main experiment both increase, the limiting variance of $\hat \theta_n$ when units are paired according to the predicted values of $E[Y_i(0) | X_i]$ is the same as that under the optimal stratification in (ref).
Next, I consider settings with data from a pilot experiment. Let $m$ denote the sample size of the pilot experiment. I assume the pilot units are drawn from the same population as the main experiment.
I start by investigating settings where the sample size of the pilot experiment is large. Formally, in the asymptotic framework, I allow both $m$ and $n$ to go to infinity. I pair units according to a suitable estimator $\tilde g_m$ of the index function $g$, where $\tilde g_m$ comes from a nonparametric regression using the pilot data. Again, a key requirement is that $\tilde g_m$ is consistent for $g$ in the sense of Assumption (ref) and Theorem (ref) below. Then, as the sample sizes of the auxiliary data and the main experiment both increase, the limiting variance of $\hat \theta_n$ when units are paired according to $\tilde g_m$ is the same as that under the optimal stratification in (ref).
If the pilot data is imperfect in the sense that it does not come from the same population as the main experiment, or if the estimation method for constructing $\tilde g_m$ is not flexible enough, then $\tilde g_m$ may not converge to $g$ but instead to another function $h$. In that case, the limiting variance of $\hat \theta_n$ is different from that under the optimal stratification in (ref) and depends on $h$, but it is still smaller than that under no stratification. See Theorem (ref) and Remark (ref) for details.
In practice, even if pilot data is available, its sample size is often small. In those settings, pairing units according to $\tilde g_m$ generally does not ensure efficiency, unlike in settings with large pilots. We may be concerned that $\tilde g_m$ is a poor approximation of the index function $g$, and as a result, if units are paired according to $\tilde g_m$, then both the conditional MSE and the limiting variance of $\hat \theta_n$ are large.
Researchers could of course ignore the information in the small pilot and implement the procedures in Section (ref). If they would like to incorporate information from the pilot experiment, they can consider the following procedure. For $d \in \{0, 1\}$, let $\tilde \beta_m(d)$ denote the OLS estimators of the linear regression coefficients among the treated or untreated units in the pilot experiment and let $\tilde \Omega_m(d)$ denote the variance estimators in OLS assuming homoskedasticity (see Section (ref) of the supplement for details). Further define
I pair the units to minimize the sum of the following distances of the covariates:
To shed some light on the behavior of such a minimization problem, I consider two extreme cases. If $\tilde \Omega_m = 0$, which means $\tilde \beta_m$ is very precise, then the solution is given by pairing units according to $\tilde g_m = X_i' \tilde \beta_m$. If $\tilde \Omega_m$ is large, which means $\tilde \beta_m$ is very imprecise, then the second term on the right-hand side of (ref) dominates the first term, so the solution is close to a paired matching weighted by $\tilde \Omega_m$. Therefore, the solution can be viewed as penalizing the pairing according to $\tilde g_m(x) = x' \tilde \beta_m$, with the penalization determined by the variance estimator $\tilde \Omega_m$. I refer to the solution as the penalized matched-pair design. In Section (ref) of the supplement, I show it is optimal in a Bayesian framework.
Each matched-pair design discussed in this section has a counterpart where units are matched into sets of four instead of pairs. Specifically, I first pair the units and then pair the pairs using the midpoints of all pairs, as in Section 4 of bai2021inference. Such a design is also discussed by athey2017econometrics. It often increases the MSE relative to its paired version but often improves inference for the ATE, especially with multiple covariates. In particular, simulation evidence in Section (ref) shows that, with multiple covariates, the test with matched sets of four usually has the correct size but the test with matched pairs often severely underrejects. I refer interested readers to Section (ref) of the supplement for a detailed discussion.
A frequent concern in experiments is attrition, meaning units in the baseline survey may drop out in the follow-up survey, so their covariates are available but outcomes are not. I emphasize that even when a unit attrites, the entire pair may not need to be dropped. If attrition happens, then I redefine the difference-in-means estimator using only non-attritors. If attrition is independent of treatment status conditional on the covariates, then this estimator is consistent for the ATE for non-attritors. I refer interested readers to Section (ref) of the supplement for details. The case with differential attrition, as in the setting of lee2009training, is an interesting topic for future work.
Another related question that frequently arises in the design of experiments is that some studies are implemented in multiple waves. Although a full-length discussion of such settings is beyond the scope of the paper, a possible solution is to implement the procedures discussed in this section repeatedly. For instance, in the first wave, researchers could pair the units according to their baseline outcomes. In the second wave, they could use the data from the first wave as pilot data, and implement the pilot-based procedures discussed earlier in this section. They can repeatedly implement the pilot-based procedures in the following waves. In Section (ref) of the supplement, I discuss how to pool the data from multiple waves for estimation and inference.
The optimality result in Section (ref) pinpoints the optimal stratification but is silent on how the feasible procedures in Section (ref) compare with each other. To make such a comparison, this section studies the asymptotic properties of the difference-in-means estimator. I also provide inference methods for the ATE under different stratifications. The main difficulty in deriving the theoretical results is that under matched-pair designs, treatment status across units is heavily dependent; in fact, treatment status of the two units in a pair is perfectly correlated. I extend the results in bai2021inference by allowing units to be paired according to functions of the covariates instead of the covariates themselves, and furthermore allowing the function to be random and dependent on auxiliary data. To begin, I make the following mild moment restriction on the distributions of potential outcomes:
I provide general results when units are paired according to a measurable function $h$ that maps from the support of $X_i$ into $\mathbf R$. The results can be easily specialized to the procedures in Section (ref). Let $\pi^h \in \Pi_n$ be such that $h_{\pi^h(1)} \leq \dots \leq h_{\pi^h(2n)}$ and define the stratification that pairs units according to $h$ as \[ \lambda^h(X^{(n)}) = \{\{\pi^h(2s - 1), \pi^h(2s)\}: 1 \leq s \leq n\}~. \]
To describe the requirements on $h$, define $\mathbf H$ to be the set of all measurable functions mapping from the support of $X_i$ into $\mathbf R$ such that the following three conditions hold:
(a) is a mild restriction to rule out degenerate situations and to permit the application of suitable laws of large numbers and central limit theorems, and (c) is another mild moment restriction to ensure the pairs are “close” in the limit. Some restrictive primitive conditions for (b) are provided in Section (ref) of the supplement. I assume $h$ lies in the set $\mathbf H$:
The next theorem establishes the limiting distribution of $\hat \theta_n$ when units are paired according to $h$, where $h$ satisfies Assumption (ref).
Next, I consider inference for the ATE when units are paired according to $h \in \mathbf H$. For any prespecified $\theta_0 \in \mathbf R$, I am interested in testing
at level $\alpha \in (0, 1)$. To do so, it suffices to provide a consistent estimator for the limiting variance $\varsigma_h^2$ in (ref). To describe such an estimator, for $d \in \{0, 1\}$, define the variance estimator among units with $D = d$ as \[ \hat \sigma_n^2(d) = \frac{1}{n} \sum_{1 \leq i \leq 2n: D_i = d} (Y_i - \hat \mu_n(d))^2~. \] In addition, define
and
The calculation in (ref) of the supplement shows $\hat \varsigma_{h, n}^2$ is nonnegative. The correction term $\hat \rho_n$ is constructed by averaging the product of the sum of the outcomes of adjacent pairs of pairs, as in bai2021inference.
The following theorem shows the variance estimator in (ref) is consistent for the limiting variance in (ref).
The discussion in the last subsection applies to settings where units are paired according to a fixed function $h \in \mathbf H$ or a random function $\tilde h_m$ such that $\tilde h_m \in \mathbf H$ with probability one. Such settings are most relevant when the pilot sample size $m$ is small. Next, I consider settings where $\tilde h_m$ converges to a fixed function $h \in \mathbf H$ in a suitable sense as $m \to \infty$. Let $Q_X$ denote the marginal distribution of $X_i$.
Assumption (ref) is commonly referred to as the $L^2$-consistency of the $\tilde h_m$ for $h$. When the dimension of $X_i$ is fixed and suitable smoothness conditions hold, $L^2$-consistency is satisfied by series and sieves estimators newey1997convergence,chen2007large and kernel estimators li2007nonparametric. In some high-dimensional settings, when the dimension of $X_i$ increases with $n$ at suitable rates, it is satisfied by the LASSO estimator buhlmann2011statistics,belloni2014inference, regression trees and random forests gyorfi2002distribution-free,wager2015adaptive, neural nets white1990connectionist,farrell2018deep, and support vector machines steinwart2008support. The results therein are either exactly as stated in Assumption (ref) or one of the following:
It is straightforward to see (a) implies Assumption (ref). Furthermore, (b) implies Assumption (ref) by Markov's inequality.
The next theorem shows that if $\tilde h_m$ is $L^2$-consistent for $h$, then as the sample sizes of both the pilot and main experiments increase, the limiting variance of $\hat \theta_n$ when units are paired according to $\tilde h_m$ is the same as that when units are paired according to $h$.
We briefly comment on inference when units are paired according to multiple covariates. For the settings in (ref) and (ref), the variance estimators are slightly more complicated than that in (ref) because the distances in (ref) and (ref) cannot be written as distances between two scalars, but the correction term is similar in spirit to (ref). I defer the discussion to Section (ref) of the supplement. In addition, note combining data from both the pilot and main experiments for estimation and inference is possible. I defer the discussion to Section (ref) of the supplement.
When units are paired using multiple covariates, the simulation evidence in Section (ref) shows that when the sample size is not large enough relative to the number of covariates, the size of the test is often strictly smaller than the nominal level. The reason is that the asymptotic results rely on the assumption that units are “close,” in the sense that a suitable normalization of the sum of distances between the covariates within each pair is close to zero. When the sample size is not large enough relative to the number of covariates, the procedures here suffer from the curse of dimensionality, so the units paired together are not close enough in terms of their covariates, and hence, the asymptotic results do not approximate the finite-sample distribution of $\hat \theta_n$ very well. The problem is mitigated by matching units into sets of four instead of pairs. Specifically, I first pair the units and then pair the pairs using the midpoints of all pairs, as in Section 4 of bai2021inference. In Section (ref) of the supplement, I propose a valid test for (ref) when units are matched into sets of four. Simulation evidence in Section (ref) shows the size of my proposed test is close to the nominal level in finite sample.
In this section, I examine the performance of the practical procedures in Section (ref) and the inference methods in Section (ref) via a simulation study calibrated to a systematically selected set of 10 RCTs from recent issues of the American Economic Journal: Applied Economics. I focus on settings with small or no pilots, because they are the most common settings in practice. I searched the 11 issues from October 2018 to April 2021 and collected 28 papers running RCTs. I exclude 11 papers for which the treatment is assigned at the cluster level instead of the unit level. I further exclude four papers with a network/spillover structure. I also exclude one paper for which the sample size is too small (less than 20). Finally, I exclude two papers for which the data is confidential. I end up with 10 papers, which are listed in Table (ref). For each paper, I list whether the baseline outcome is available, the original randomization method, and the number of additional covariates besides the baseline outcome in the main regression specification of the paper. The full details of the data are available in Section (ref) of the supplement.
For each paper, I denote the sample size by $2n$. I use the original sample except for barrera-osorio2019medium-, where the original data contains 15,759 observations and 16 covariates, so one replication in the simulation takes almost six hours. For barrera-osorio2019medium- only, I take half of the observations as the population to reduce the computational time for one replication to about an hour, which is about the same as that using the next largest dataset. I begin by imputing the unobserved potential outcomes. For the $i$th unit, I denote the original data by $(Y_i^\ast, D_i^\ast, X_{1i}^\ast, X_{2i}^\ast)$, where $Y_i^\ast$ denotes its observed outcome, $D_i^\ast$ denotes its treatment status, $X_{1i}^\ast$ denotes its baseline outcome if available, and $X_{2i}^\ast$ denotes the other covariates in the main regression specification of the paper. Let $Y_i^\ast(1), Y_i^\ast(0)$ denote the potential outcomes for the $i$th unit. For the $i$th unit, $Y_i^\ast(D_i)$ is observed, and I construct $Y_i^\ast(1 - D_i)$ according to the following models:
In Model 1, treatment effects are homogeneous, so (ref) holds. In Models 2 and 3, treatment effects are heterogeneous. Note the baseline outcome predicts the potential outcomes better in Model 3 than in Model 2.
For each replication, I simulate new data $((Y_i(1), Y_i(0), X_{i1}, X_{2i}): 1 \leq i \leq 2n)$ by drawing $2n$ units from the empirical distribution of $((Y_i^\ast(1), Y_i^\ast(0), X_{1i}^\ast, X_{2i}^\ast): 1 \leq i \leq 2n)$, so each unit in the original data is drawn with equal probability, with replacement.
For each paper, I implement several stratifications. Note the covariates may have been selected ex-post by authors on the basis of predictive power, while ideally I would like to include only the covariates specified in the pre-analysis plans. Unfortunately, only one paper includes such information in the AEA RCT Registry, and the pre-registered covariates are the same as those in the regression analysis. I stratify on the baseline outcome whenever it is available. To stratify based on data from pilot experiments, I reached out to the authors of all 10 papers to request pilot data. 9 out of 10 replied, 8 of whom said they did not run a pilot, and the last said they ran a pilot, but the data was lost. Therefore, I simulate pilot data by drawing with replacement from the empirical distribution at a sample size of $\lfloor 0.2 \cdot (2n) \rfloor$. To study the setting with a small and fixed pilot, I fix the pilot data throughout all replications. I also consider several stratifications with matched sets of four, as in athey2017econometrics. Specifically, strata are constructed by first pairing the units and then pairing the pairs according to their midpoints, as in Section 4 of bai2021inference. The complete list of stratifications are:
The original stratifications are listed in Table (ref). I do not consider re-randomization in lee2021poverty, because the exact implementation is unclear from the original paper, and inference under re-randomization is complicated. See also Remark (ref) for a comparison between re-randomization and matched-pair designs. Section (ref) of the supplement contains the results for several additional stratifications. Although it is interesting to investigate the performance of regression adjustment with the stratifications in Origin, inference with regression adjustment under stratified randomization is still an open question.
I consider the following inference methods:
For matched sets of four, athey2017econometrics propose a test in a sampling framework different from ours. I show in Section (ref) of the supplement that because of the differences in sampling frameworks, the test in athey2017econometrics does not control size in my setting unless the conditional ATEs are homogeneous. Therefore, I defer these simulation results to Section (ref) of the supplement.
For each paper, each model, and each stratification, across $1,000$ replications, I calculate three metrics of performance: (1) the MSE of estimating $\theta$ using $\hat \theta_n$, reflecting the precision of the estimator; (2) the average rejection probability of testing (ref) for $\theta_0 = \theta$, reflecting the size of the test; and (3) the average standard error, which directly determines the length of the confidence interval for the ATE.
The main results of the simulation study are summarized in Table (ref). I only report the summary statistics across all papers and models and defer the raw numbers to Section (ref) of the supplement. In particular, for each stratification, I report the average and $[\min, \max]$ across all papers and models of
Rows are labeled according to the stratifications.
In this subsection, I discuss major takeaways about statistical precision (specifically, MSE) from Table (ref). I focus on five questions that are particularly relevant to empirical practice.
First, how much statistical precision are researchers leaving on the table with their current stratification methods? To answer this question, I compare the MSEs under the stratifications used in the original paper (Origin) and the MSEs when pairing according to the baseline outcome (MP base). Relative to the original stratifications used in those 10 papers, if the researchers had just paired the units according to their baseline outcomes, the MSE would be 24% smaller on average and 56% smaller in some cases. In fact, in many models, the MSE under original stratification is almost the same as the MSE when not stratifying. As a result, if researchers had paired the units according to their baseline outcomes, they could have reached the same statistical precision with a much smaller sample size.
Second, how much statistical precision would researchers leave on the table by pairing units according to the baseline outcome rather than using more complicated feasible methods? Note MP X usually has the smallest MSEs across all methods. On average, the MSE under MP base is about 39% larger than that under MP X, and 10% larger than that under MS X. As a result, researchers indeed sacrifice some statistical precision by pairing units according to the baseline outcome along instead of pairing or matching into sets of four according to all covariates. Note, however, that MP base picks up about half of the difference between the MSEs under the best feasible method (MP X) and the status-quo methods (Origin).
Third, what is the value of collecting just the baseline outcome rather than other covariates? To answer this question, I compare the MSEs under MP base and MP X2. Note the MSE under MP base is on average only 11% larger than and sometimes almost the same as that under MP X2. As a result, the statistical precision of pairing according to the baseline outcome alone is comparable to the statistical precision of pairing according to all other covariates. Note the number of other covariates is close to or larger than 10 in most cases, so the per-covariate return for collecting all of them is limited relative to collecting the baseline outcome alone.
Fourth, are matched sets of four better than matched pairs in terms of the MSE? To answer this question, I compare the MSEs of the MS methods and the MP methods. The MSE under MS base is 4% larger than that under MP base, and the MSE under MS X is 27% larger than that under MP X. Therefore, the statistical precision is higher with matched pairs. The difference is pronounced when I match according to multiple covariates but tiny when I match only according to the baseline outcome.
Fifth, does the best pairing method based on a small pilot dominate pairing according to $X$? In other words, what is the value of having a small pilot? First, note the MSE under MP pen is 19% smaller than that under MP pilot, so the penalized matched-pair design indeed has better precision than the naïve plug-in procedure. Meanwhile, the MSE under MP pen is almost the same as that under MP X, and even in the most favorable case, it is only about 8% smaller. Therefore, the return for a small pilot in terms of the MSE is negligible.
I also study the performance of regression-adjusted estimators in stratifications with one stratum (None-reg). With one stratum, the regression-adjusted estimator usually has slightly smaller MSEs than the difference-in-means estimator. In almost all cases, however, the MSE is larger than those under all methods with matched pairs or matched sets of four, regardless of whether all or only a subset of the covariates in the regression adjustment are used in the matching. Section (ref) of the supplement contains the results for the regression-adjusted estimator in lin2013agnostic, which additionally includes the interactions of treatment status and covariates. The results are qualitatively similar to those for None-reg. Therefore, most of the gains in precision from matching according to the covariates cannot be retrieved by controlling for the same covariates via ex-post regression adjustment.
Next, I discuss the properties of the inference methods. I start with the size of the tests. For matched pairs, note both the adjusted $t$-test and the “matched pairs” $t$-test control size well across all papers and models. When I pair units according to the baseline outcome (MP base), the size of the adjusted $t$-test is almost always close to 5%. The size of the “matched pairs” $t$-test is also close to 5%. The underrejection phenomenon for the “matched pairs” $t$-test in Remark (ref) is very mild, reflecting the treatment effects heterogeneity is not very large. Several relatively noticeable cases include: for Model 2 of paper 5, the size of the “matched pairs” $t$-test is 3.9% but the size of the adjusted $t$-test is 5.4%; for Model 3 of paper 5, the size of the “matched pairs” $t$-test is 2.8% but the size of the adjusted $t$-test is 3.2%.
When I pair units according to multiple covariates (MP X, MP X2, and MP pen), the “matched pairs” $t$-test is still conservative except in Model 1, for the same reason as mentioned in Remark (ref). At the same time, the adjusted $t$-test also becomes conservative---its size is often smaller than 5%. The reason is that the asymptotic results in this paper rely on the assumption that units are “close,” in the sense that a suitable normalization of the sum of the distances between the covariates within each pair is close to zero. When pairing according to multiple covariates, however, the procedures suffer from the curse of dimensionality, so the units paired together are not close enough in terms of their covariates. Therefore, the asymptotic results do not approximate the finite-sample distribution of $\hat \theta_n$ very well, and my variance estimator does not approximate the actual variance of $\hat \theta_n$ very well.
When matching according to multiple covariates, the conservativeness of the tests is somewhat alleviated by matching units into sets of four instead of pairs. The size of the test under MS X is close to 5% even when the size of the test under MP X is much smaller than 5%. On the other hand, such a difference is virtually nonexistent when I match only according to the baseline outcome---the size of the test under MP base and that under MS base are both close to 5%. Our current asymptotic framework cannot explain the difference in size because the variance estimators for both MP and MS methods are consistent for the limiting variance. The exact reason for the difference is an interesting topic for future work.
I now turn to the standard errors. The findings are mostly similar to those for the MSEs, though with some important exceptions. The standard error under MP base is 10% smaller on average and 34% smaller in some cases than that under Origin. At the same time, although the MSE under MP X is smaller than that under MS X, the standard error of MP X is larger than that under MS X. In fact, the standard error under MP X is about the same as that under MP base, and the standard error under MP X2 is often larger than that under MP base. Therefore, although pairing according to multiple covariates is desirable for the MSE, it is often not the best choice for inference, because the standard error is too large and the size of the test could be strictly smaller than the nominal level. By matching units into sets of four instead of pairs according to the same set of covariates, researchers could lower the standard error and bring the size close to the nominal level. Note, however, that the MSE will increase, as discussed in Section (ref).
I emphasize the validity of my tests relies on the assumptions on the sampling framework. In this paper, I assume units are drawn from a superpopulation, and the potential outcomes and the covariates are random. de_chaisemartin2021at, on the other hand, study a finite-population setting in which the potential outcomes and the covariates are fixed. Such a setting is particularly relevant if we have a convenience sample instead of a random sample drawn from a large population. de_chaisemartin2021at show in these settings that if the number of pairs is small, then the tests in my paper and bai2021inference may not control size, and the “matched pairs” $t$-test in imbens2015causal could become preferable.
Finally, I consider settings with multiple outcomes. I take the example of abel2019bridging, where a primary outcome, a secondary outcome, and the baseline outcomes of both are available. The paper studies job-searching behaviors. The primary outcome is the search hours and the secondary outcome is the number of applications sent. I study the estimation of the ATE of the secondary outcome. The missing potential outcomes are imputed as in Model 1, assuming the treatment effect is zero for everyone, and Model 3, using the nearest neighbor in terms of the baseline value of the secondary outcome. I consider the following stratifications:
In light of the results in Section (ref), I only consider the adjusted $t$-tests. For each stratification, I calculate the MSE, the size of testing (ref) with $\theta_0 = \theta$, and the average standard error. The results are displayed in Table (ref). Rows are labeled according to the stratifications. As expected, because the secondary outcome is of interest, for both models, stratifying on the baseline value of the secondary outcome produces smaller MSEs than stratifying on the baseline value of the primary outcome. For both models, the MSEs when stratifying on the baseline value of the primary outcome are close to the MSEs with one stratum, reflecting the baseline value of the primary outcome is a poor predictor of the secondary outcome. The smallest MSE is attained by pairing according to the baseline values of both outcomes, and the second smallest MSE is attained by matching units into sets of four according to both baseline outcomes. In all models, the size of the test is close to the nominal level. The test under MP 1+2 slightly underrejects for the same reason as in Section (ref), although the problem is mild because I only match on two variables. The standard errors are ranked in the same way as the MSEs except for that of MP 1+2. In both models, MS 1+2 produces the smallest standard errors across all methods.
Based on the theoretical results, in settings with large pilots, I recommend researchers to pair units according to the estimated index function from nonparametric regressions. If the conditional ATEs are homogeneous and researchers have access to a large observational dataset from the same population as that of the main experiment, then I recommend pairing according to predicted values from nonparametric regressions of the outcome on the covariates in the observational dataset. In what follows, I focus on settings with small or no pilots because these are the most common settings in practice.
A simple approach that researchers can take, assuming there is only one primary outcome of interest and its baseline value is available, is to pair units according to the baseline outcome. Indeed, if the baseline outcome is the only available covariate and the index function is monotonic in it, then my theoretical results show pairing units according to the baseline outcome is optimal at any sample size. The simulation results in Section (ref) also show pairing according to the baseline outcome improves upon the status-quo methods in terms of both the MSE and the standard error of the difference-in-means estimator. Further, this approach has the advantage of simplicity.
When multiple covariates are available, an attractive alternative is to match units into pairs or sets of four according to the baseline outcome and other covariates. Unless the number of covariates is very small, I recommend matched sets of four over pairs because the standard error is usually smaller with matched sets of four. In my simulation study, when I use the baseline outcome together with all the covariates that the authors control for in their regressions, matching units into sets of four leads to smaller MSEs and standard errors than pairing on the baseline outcome alone. Note that the good performance of this design could be due to the fact that the authors selected the covariates with the best predictive power ex-post, something that is not feasible at the time of randomization. Nevertheless, forming sets of four according to the baseline outcome and other covariates is an attractive alternative to pairing according to the baseline outcome alone.
In my simulation study, when multiple outcomes are of interest, pairing on one of them may not improve the MSE of the other outcomes. In those settings, researchers could consider matching units into sets of four to minimize the sum of the squared Mahalanobis distances of the baseline values of all outcomes of interest and perhaps some additional covariates.
A further question is whether pilot experiments are worth running for the sole purpose of improving the precision for estimating the ATE. My simulation results only show minor gains in statistical precision when using pilot-based stratifications instead of matching directly on the covariates. Therefore, although pilot experiments are essential for other aspects of the design of the main experiment, they are not as helpful in improving the precision of the estimator for the ATE.
Another natural question is whether one can retrieve the gains in precision from matched pairs or sets of four units by controlling for the same set of covariates through ex-post regression adjustment. Our simulation results show the answer is negative, although with one stratum, regression adjustment slightly lowers the MSE and the standard error relative to the unadjusted difference-in-means estimator. I further note the difference-in-means is unbiased for the ATE in finite sample under all stratifications considered in this paper, while the regression-adjusted estimators are only consistent for the ATE asymptotically. A very interesting direction for future work is to combine regression adjustment with stratifications defined by matched pairs or matched sets of four units.
For inference, researchers can use the test with the variance estimator in (ref). They can also use the test in Theorem 10.1 of imbens2015causal, which is valid albeit sometimes conservative. Finally, I emphasize that my framework assumes units are drawn from a superpopulation and the potential outcomes and the covariates are random. If we have a convenience sample instead of a random sample drawn from a large population, and the sample size is small, then de_chaisemartin2021at show the tests in this paper and bai2021inference may not control size, and the test in imbens2015causal could become preferable.