EconBase
← Back to paper

On the Efficiency of Finely Stratified Experiments

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.

92,858 characters · 11 sections · 102 citation commands

Rendered from LaTeX for readability, not typeset faithfully. Citation keys are highlighted; maths is left as source; figures, tables and equation environments are summarised rather than reproduced; unrecognised commands are greyed out so nothing is silently dropped. Email addresses are removed.

On the Efficiency of Highly Stratified Experiments

spacing{1.1} \begin{abstract} This paper studies the use of highly stratified designs for the efficient estimation of a large class of treatment effect parameters that arise in the analysis of experiments. By a “highly stratified” design, we mean experiments in which units are divided into blocks of a fixed size and a proportion within each block is assigned to a binary treatment uniformly at random. The class of parameters considered are those that can be expressed as the solution to a set of moment conditions constructed using a known function of the observed data. They include, among other things, average treatment effects, quantile treatment effects, and local average treatment effects as well as the counterparts to these quantities in experiments in which the unit is itself a cluster. In this setting, we establish three results. First, we show that under a highly stratified design, the na\"ive method of moments estimator achieves the same asymptotic variance as what could typically be attained under alternative treatment assignment mechanisms only through {\it ex post} covariate adjustment. Second, we argue that the na\"ive method of moments estimator under a highly stratified design is asymptotically efficient by deriving a lower bound on the asymptotic variance of regular estimators of the parameter of interest in the form of a convolution theorem. In this sense, highly stratified experiments are attractive because they lead to efficient estimators of treatment effect parameters “by design.” Finally, we strengthen this conclusion by establishing conditions under which a “fast-balancing” property of highly stratified designs is in fact necessary for the na\"ive method of moments estimator to attain the efficiency bound. \end{abstract}

KEYWORDS: Convolution theorem, Efficiency, Experiment, Experimental design, Highly stratified experiment, Matched pairs, Randomized controlled trial

JEL classification codes: C12, C14

\thispagestyle{empty} \setcounter{page}{1}

Introduction

This paper studies the use of highly stratified designs for the efficient estimation of a large class of treatment effect parameters that arise in the analysis of experiments. By a “highly stratified” design, we mean experiments in which units are divided into blocks of a fixed size based on their covariate values and a proportion within each group is assigned to a binary treatment uniformly at random. The canonical example of such a design is a matched pairs design: here the fixed size of the blocks equals two, and exactly one of the two units in each block is assigned to treatment at random, so that the marginal probability of treatment assignment is one half. More broadly, a “highly" stratified experiment is also intended to generalize the notion of a “finely" stratified experiment fogarty2018mitigating, for which in every block there is exactly one treated or control individual. On the other hand, our use of the term “highly stratified" in this context contrasts with what we could call a “coarsely stratified" design, in which units are divided into a small set of large blocks: see Example (ref) below for a formal definition. The class of parameters considered are those that can be expressed as the solution to a set of moment conditions constructed using a known function of the observed data. This class of parameters includes many causal parameters of interest: average treatment effects (ATEs), quantile treatment effects, and local average treatment effects as well as the counterparts to these quantities in experiments in which the unit is itself a cluster.

In the setting described above, we establish three results. First, we study the asymptotic properties of a na\"ive method of moments estimator under a highly stratified design. Here, by a na\"ive method of moments estimator, we mean an estimator constructed using a direct sample analog of the moment conditions. For example, in the case of the ATE, such an estimator is given by the Horvitz-Thompson estimator for the difference in means. We show that under a highly stratified design, the na\"ive method of moments estimator achieves the same asymptotic variance as what could typically be attained under alternative treatment assignment mechanisms only through {\it ex post} covariate adjustment using the same set of covariates. Such adjustment strategies frequently involve the nonparametric estimation of conditional expectations or similar quantities; see, for example, zhang2008improving, tsiatis2008covariate, jiang2022improving, jiang2022regression-adjusted and rafi2023efficient. We further illustrate that this feature of highly stratified experiments stems from the way in which highly stratified designs balance treatment status across covariate values, a property we define formally below and refer to as “fast-balancing.” Second, we derive a lower bound on the asymptotic variance of regular estimators of the parameter of interest in the form of a convolution theorem. This convolution theorem accommodates a large class of possible treatment assignment mechanisms, including covariate adaptive randomization efron1971forcing, wei1978adaptive,zelen1974randomization, pocock1975sequential,hu2012asymptotic,bugni2018inference,ye2022inference,ma2020statistical,ma2024new, re-randomization li2017general,li2018asymptotic,li2020rerandomization,li2020rerandomization-1, cytrynbaum2024finely, and highly stratified designs jiang2021bootstrap, bai2022inference, cytrynbaum2023designing. We show that the lower bound is attained by the na\"ive method of moments estimator under a highly stratified design. In this sense, the na\"ive method of moments estimator under a highly stratified design is asymptotically efficient. More succinctly, we say that highly stratified experiments lead to efficient estimators “by design.” Finally, we strengthen this conclusion by characterizing all regular asymptotically linear estimators for a large class of treatment assignment mechanisms and use this characterization to establish conditions under which the fast-balancing property of highly stratified experiments is in fact a necessary condition for the na\"ive method of moments estimator to attain the efficiency bound.

Together, these results demonstrate that highly stratified experiments lead to efficient estimators that prioritize transparency in that they preclude the researcher from “data snooping” associated with {\it ex post} nonparametric covariate adjustment. Importantly, concerns with this type of data snooping are not completely eliminated by typical pre-analysis plans because such adjustments involve choices, such as the choice of nonparametric estimator or tuning parameters, that are often not pre-registered prior to the experiment. The estimators are therefore attractive because they avoid performing nonparametric covariate adjustment in order to achieve efficiency and thereby remain “hands above the table” freedman2008regression, lin2013agnostic.

Our paper builds upon two strands of literature. The first strand of literature concerns the analysis of highly stratified experiments. Within this literature, our analysis is most closely related to bai2022inference, who derive the asymptotic behavior of the difference-in-means estimator of the ATE when treatment is assigned according to a matched pairs design, and cytrynbaum2023designing, who develops related results for an experimental design referred to as “local randomization” that permits the proportion of units assigned to treatment to vary with the baseline covariates. Beyond settings that study estimation of the ATE, bai2024inference-1 develops results for the analysis of different cluster-level average treatment effects and jiang2021bootstrap develop results analogous to those in bai2022inference for suitable estimators of the quantile treatment effect. This paper, like those just mentioned, operates in a “super-population” framework, in which the outcomes and covariates are assumed to be drawn as an i.i.d.\ sample from a population distribution. This is in contrast to an alternative strand of the literature that studies highly stratified experiments from the design-based perspective imai2008variance, imai2009essential, fogarty2018mitigating, fogarty2018regression-assisted, liu2020regression-adjusted, pashley2021insights, bai2025new. To our knowledge, our paper is the first to analyze the properties of highly stratified experiments in a general framework that accommodates any parameter that can be characterized as the solution to a set of moment conditions involving a known function of the observed data. In more recent work, cytrynbaum2024finely considers a similar framework to study highly stratified re-randomized experiments, which nest highly stratified experiments as a special case. However, his results assume that the moment functions which define the parameter of interest are continuous in a way that precludes parameters like the Quantile Treatment Effect. Moreover, we emphasize that none of the above papers formally establish the asymptotic efficiency of highly stratified experiments. The second strand of literature concerns bounds on the efficiency with which treatment effect parameters can be estimated in experiments. We note that, due to the potential for dependence in treatment assignments across individuals, we cannot immediately appeal to standard semi-parametric efficiency results van_der_vaart1998asymptotic, chen2008semiparametric. Two important recent papers in this literature studying efficiency bounds in the special case of estimating the ATE are armstrong2022asymptotic and rafi2023efficient. Even in this special case, their results differ from ours in important and empirically relevant ways; Remark (ref) provides an in-depth discussion of the connection between these results and ours. See also bai2022optimality for some finite-sample optimality properties of matched pairs designs for estimation of the ATE.

The remainder of this paper is organized as follows. In Section (ref), we describe our setup and notation. We emphasize, in particular, the way in which our framework can accommodate various treatment effect parameters of interest. Section (ref) derives the asymptotic behavior of the na\"ive method of moments estimator of our parameter of interest when treatment is assigned using a highly stratified design and studies estimation of the asymptotic variance. In Section (ref), we develop our lower bound on the asymptotic variance of regular estimators of these parameters and show that it is achieved by the the na\"ive method of moments estimator in a highly stratified design. In Section (ref), we characterize all regular asymptotically linear estimators and argue that the fast-balancing property is necessary for the na\"ive method of moments estimator to attain the efficiency bound. In Section (ref), we illustrate the practical relevance of our theoretical results through a simulation study. Finally, we conclude in Section (ref) with some recommendations for empirical practice guided by both these simulations and our theoretical results. Proofs of all results can be found in the supplementary material.

Setup and Motivation

Let $A_i \in \{0, 1\}$ denote the treatment status of the $i$th unit, and let $X_i \in \mathbf R^{d_x}$ denote their observed, baseline covariates. For $a \in \{0, 1\}$, let $R_i(a) \in \mathbf R^{d_r}$ denote a vector of potential responses. As we illustrate below, considering a vector of responses allows us to accommodate many parameters of interest. Let $R_i \in \mathbf R^{d_r}$ denote the vector of observed responses obtained from $R_i(a)$ once treatment is assigned. As usual, the observed responses and potential responses are related to treatment status by the relationship

equation[equation omitted — 66 chars of source]

We assume throughout that our sample consists of $n$ units. For any random vector indexed by $i$, for example $A_i$, we define $A^{(n)} = (A_1, \ldots, A_n)$. Let $P_n$ denote the distribution of the observed data $(R^{(n)}, A^{(n)}, X^{(n)})$, and $Q_n$ the distribution of $(R^{(n)}(1), R^{(n)}(0), X^{(n)})$. We assume $Q_n = Q^n$, where $Q$ is the marginal distribution of $(R_i(1), R_i(0), X_i)$. Given $Q_n$, $P_n$ is then determined by (ref) and the mechanism for determining treatment assignment. We assume that treatment assignment is performed such that a standard unconfoundedness assumption holds and such that the probability of assignment given $X_i$ is some known constant for every $1 \le i \le n$, as is often the case in most experiments:

assumptionTreatment status is assigned so that \begin{equation} (R^{(n)}(1), R^{(n)}(0)) \perp \!\!\! \perp A^{(n)} | X^{(n)} , \end{equation} and such that $P\{A_i = 1|X_i=x\} = \eta$, for some $\eta \in (0, 1)$ for all $1 \le i \le n$.

Assumption (ref) restricts the probability of assignment to be the fixed fraction $\eta$ across the entire experimental sample, but this restriction can be weakened so that $\eta$ is replaced by $\eta(X_i)$ for many of our subsequent results: see Remark (ref) for a discussion. Given Assumption (ref), it can be shown that $(X_i, A_i, R_i)$ are identically distributed for $1 \le i \le n$, and their marginal distribution does not change with $n$ (see Lemma A.7 in the Supplement). As a consequence, we denote the marginal distribution of $(X_i, A_i, R_i)$ by $P$. We consider parameters $\theta_0 \in \Theta \subset \mathbf R^{d_\theta}$ that can be defined as the solution to a set of moment equalities. As we show below, these parameters include a large class of causal parameters defined in terms of potential outcomes and potential treatments. Formally, let $m: \mathbf R^{d_x} \times \{0, 1\} \times \mathbf R^{d_r} \to \mathbf R^{d_\theta}$ be a known measurable function, then we consider parameters $\theta_0$ that uniquely solve the moment equality

equation[equation omitted — 72 chars of source]

We emphasize that $m(\cdot)$ is not a function of any unknown nuisance parameters, but may depend on the known value of $\eta$ in Assumption (ref). We present five examples of well-known parameters that can be described as (functions of) solutions to a set of moment conditions as in (ref).

example[Average Treatment Effect] Let $Y_i(a) = R_i(a)$ denote a scalar potential outcome for the $i$th unit under treatment $a \in \{0, 1\}$, and let $Y_i = R_i$ denote the observed outcome. Let $\theta_0 = E_Q[Y_i(1) - Y_i(0)]$ denote the average treatment effect (ATE). Under Assumption (ref), $\theta_0$ solves the moment condition in (ref) with \begin{equation} m(X_i, A_i, R_i, \theta) = \frac{Y_i A_i}{\eta} - \frac{Y_i (1 - A_i)}{1 - \eta} - \theta . \end{equation} For papers that consider estimators based on (ref), see hirano2001estimation and hirano2003efficient.
example[Quantile Treatment Effect] Let $Y_i(a) = R_i(a)$ denote a scalar potential outcome for the $i$th unit under treatment $a \in \{0, 1\}$, and let $Y_i = R_i$ denote the observed outcome. Let $\tau \in (0, 1)$ and $\theta_0 = (\theta_0(1), \theta_0(0))' = (q_{Y(1)}(\tau), q_{Y(0)}(\tau))'$, where \[ q_{Y(a)}(\tau) = \inf \{\lambda \in \mathbf R: Q \{Y_i(a) \leq \lambda \} \geq \tau\}~. \] In other words, $\theta_0$ is defined to be the vector of $\tau$th quantiles of the marginal distributions of $Y_i(1)$ and $Y_i(0)$. If we assume $q_{Y(a)}(\tau)$ is unique for $a \in \{0, 1\}$ in the sense that $Q\{Y(a) \leq q_{Y(a)}(\tau) + \epsilon\} > Q\{Y(a) \leq q_{Y(a)}(\tau)\}$ for all $\epsilon > 0$, then it follows from Assumption (ref) and Lemma 1 in firpo2007efficient that $\theta_0$ solves the moment condition in (ref) with \[ m(X_i, A_i, R_i, \theta) = \begin{pmatrix} \displaystyle \frac{A_i (\tau - I \{Y_i \leq \theta(1)\})}{\eta} \\ \displaystyle \frac{(1 - A_i) (\tau - I \{Y_i \leq \theta(0)\})}{1 - \eta} \end{pmatrix}~, \] for $\theta = (\theta^{(1)}, \theta^{(0)})'$. Note that the quantile treatment effect $q_{Y(1)}(\tau) - q_{Y(0)}(\tau)$ can then be defined as $h(\theta_0)$ where $h:\mathbf R^2 \to \mathbf R$ is given by $h(s,t) = s - t$.
example[Local Average Treatment Effect] Let $(\tilde{Y}_i(a), D_i(a)) = R_i(a)$ denote the vector of potential outcomes ($\tilde{Y}_i(a) \in \mathbf R$) and treatment take-up ($D_i(a) \in \{0, 1\}$) under treatment $a \in \{0, 1\}$, and let $(Y_i, D_i) = R_i$ denote the vector of observed outcomes and treatment take-up. Note here that $\tilde{Y}_i(a)$ corresponds to the potential outcome under assignment $a \in \{0,1\}$ and not to the potential outcome for a given take-up $D_i = d$. Suppose $E_Q[D_i(1) - D_i(0)] \ne 0$ and let \[ \theta_0 = \frac{E_Q[\tilde{Y}_i(1) - \tilde{Y}_i(0)]}{E_Q[D_i(1) - D_i(0)]}~. \] It then follows from Assumption (ref) that $\theta_0$ solves the moment condition in (ref) with \begin{equation} m(X_i, A_i, R_i, \theta) = \frac{Y_i A_i}{\eta} - \frac{Y_i (1 - A_i)}{1 - \eta} - \theta \left ( \frac{D_i A_i}{\eta} - \frac{D_i (1 - A_i)}{1 - \eta} \right ) . \end{equation} If we further assume instrument monotonicity (i.e., $P\{D_i(1) \ge D_i(0)\} = 1$) and instrument exclusion, then $\theta_0$ could be re-interpreted as the local average treatment effect (LATE) in the sense of imbens1994identification.
example[Weighted Average Treatment Effect] Let $Y_i(a) = R_i(a)$ denote a scalar potential outcome for the $i$th unit under treatment $a \in \{0, 1\}$, and let $Y_i = R_i$ denote the observed outcome. Let \[\theta_0 = E_Q\left[\frac{\omega(X_i)}{E_Q[\omega(X_i)]}\left(Y_i(1) - Y_i(0)\right)\right]~,\] for some known function $\omega: \mathbf R^{d_x} \to \mathbf R$. It then follows from Assumption (ref) that $\theta_0$ solves the moment condition in (ref) with \[ m(X_i, A_i, R_i, \theta) = \omega(X_i)\left(\frac{Y_i A_i}{\eta} - \frac{Y_i (1 - A_i)}{1 - \eta}\right) - \omega(X_i)\theta~.\] By defining $Y_i$ to be the average outcome in the $i$th cluster, $\theta_0$ defined in this way can accommodate the (cluster) size-weighted and equally-weighted average treatment effects considered in bugni2022inference and bai2024inference-1 in the context of cluster-level randomized controlled trials.
example[Log-Odds Ratio] Let $Y_i(a) = R_i(a) \in \{0, 1\}$ denote a binary potential outcome for the $i$th unit under treatment $a \in \{0, 1\}$, and let $Y_i = R_i$ denote the observed outcome. Suppose $0 < P\{Y_i(a) = 0\} < 1$ for $a \in \{0, 1\}$, and let $\theta_0 = (\theta_0(1), \theta_0(2))'$, where \[\theta_0(1) = \text{logit}(E_Q[Y_i(0)])~,\] \[\theta_0(2) = \text{logit}(E_Q[Y_i(1)]) - \text{logit}(E_Q[Y_i(0)])~,\] with $\text{logit}(z) = \log(\frac{z}{1-z})$, so that $\theta_0(2)$ denotes the log-odds ratio of treatment $1$ relative to treatment $0$. It follows from Assumption (ref) that $\theta_0$ solves the moment condition in (ref) with \[m(X_i, A_i, R_i, \theta) = \begin{pmatrix} \displaystyle 1 - A_i \\ \displaystyle A_i \end{pmatrix}\left(Y_i - \text{expit}(\theta(1) + \theta(2)A_i)\right)~,\] where $\text{expit}(z) = \frac{\exp(z)}{1 + \exp(z)}$. The log-odds ratio can then be defined as $h(\theta_0)$ where $h: \mathbf R^2 \to \mathbf R$ is given by $h(s,t) = t$. This parameter appears in, for example, zhang2008improving.

Additional examples could be obtained by considering combinations of Examples (ref)--(ref). For instance, combining the moment functions from Examples (ref) and (ref) would result in a weighted LATE parameter. Beyond these examples, certain treatment effect contrasts could also be related to the structural parameters in, for instance, an economic model of supply and demand: see, for example, the model estimated in casaburi2022using.

Throughout the rest of the paper we consider the asymptotic properties of the method of moments estimator $\hat{\theta}_n$ for $\theta_0$ which is constructed as a solution to the sample analogue of (ref):

equation[equation omitted — 103 chars of source]

Because $\hat{\theta}_n$ is constructed directly using the moment function $m(\cdot)$, we call $\hat{\theta}_n$ the na\"ive method of moments estimator. Note that $\hat \theta_n$ as defined in (ref) is closely related to standard estimators of the parameter $\theta_0$ in specific examples. For instance, in Example (ref), \[\hat{\theta}_n = \frac{1}{\eta}\sum_{1 \le i \le n}Y_iA_i - \frac{1}{1 - \eta}\sum_{1 \le i \le n}Y_i(1 - A_i)~,\] so that $\hat{\theta}_n$ is a Horvitz-Thompson analogue of the standard difference-in-means estimator for the ATE. In Example (ref), \[\hat{\theta}_n = \frac{\frac{1}{\eta}\sum_{1 \le i \le n}Y_iA_i - \frac{1}{1 - \eta}\sum_{1 \le i \le n}Y_i(1 - A_i)}{\frac{1}{\eta}\sum_{1 \le i \le n}D_iA_i - \frac{1}{1 - \eta}\sum_{1 \le i \le n}D_i(1 - A_i)}~,\] so that $\hat{\theta}_n$ is a Horvitz-Thompson analogue of the standard Wald estimator for the local average treatment effect.

Before proceeding, in the remainder of this section, we provide a more detailed summary of the main contributions of our paper. To this end, first note that if $A^{(n)}$ were assigned i.i.d., independently of $X^{(n)}$, then it can be shown under mild conditions on $m(\cdot)$ van_der_vaart1998asymptotic that the na\"ive method of moments estimator satisfies \[\sqrt{n}(\hat{\theta}_n - \theta_0) \overset{d}{\to} N(0, \mathbb{V})~,\] where

equation[equation omitted — 136 chars of source]

with $M = \frac{\partial}{\partial \theta'} E_P[m(X, A, R, \theta)] \Big |_{\theta = \theta_0}$. In Section (ref), we show that if we assign $A^{(n)}$ using a highly stratified design (see Assumption (ref) below for a formal definition) then, under appropriate assumptions so that we achieve “fast balance” of the treatment across covariate values (see Assumptions (ref), (ref) below), \[\sqrt{n}(\hat{\theta}_n - \theta_0) \overset{d}{\to} N(0, \mathbb{V}_*)~,\] where $\mathbb{V} \ge \mathbb{V}_*$ (see Theorem (ref)). Under i.i.d.\ assignment, the na\"ive method of moment estimator $\hat \theta_n$ cannot generally attain $\mathbb V_\ast$, but an estimator that attains $\mathbb{V}_*$ could instead be constructed by appropriately “augmenting” the moment function, and then considering an estimator which solves the augmented moment equation. For instance, if we consider the ATE in Example (ref), then it is straightforward to show that the following augmented moment function identifies $\theta_0$:

equation[equation omitted — 183 chars of source]

where $\mu_a(X_i) = E_Q[Y_i(a)|X_i]$. This choice of $m^*(\cdot)$ produces the well known doubly-robust moment condition for estimating the ATE robins1995analysis,hahn1998role. It can then be shown that an appropriately constructed two-step estimator, in which $\mu_1(\cdot)$ and $\mu_0(\cdot)$ are non-parametrically estimated in a first step, attains $\mathbb{V}_*$ tsiatis2008covariate,farrell2015robust,chernozhukov2017doubledebiasedneyman, rafi2023efficient. Intuitively, the estimator obtained from the augmented moment function $m^*(\cdot)$ performs nonparametric covariate adjustment by exploiting the information contained in $X^{(n)}$ that may not have been captured in the original moment function $m(\cdot)$. Similar nonparametric covariate adjustments based on augmented moment equations have been developed for other parameters of interest zhang2008improving, belloni2017program,jiang2022improving,jiang2022regression-adjusted. In this sense, we show that highly stratified designs can perform nonparametric covariate adjustment “by design” for the large class of parameters that can be expressed in terms of moment conditions of the form given in (ref), thus generalizing similar observations made in bai2022inference, bai2022optimality, and cytrynbaum2023designing in the special case of estimating the ATE. As we explain in the discussion following Theorem (ref), highly stratified experiments have this feature because they lead to fast-balancing of the treatment across covariate values, as defined formally in (ref) below.

Earlier work on efficient treatment effect estimation has noted that the variance $\mathbb{V}_*$ is in fact the efficiency bound for estimating $\theta_0$ under i.i.d.\ assignment cattaneo2010efficient. A natural follow-up question is whether or not $\mathbb{V}_*$ continues to be the efficiency bound for estimating $\theta_0$ under a highly stratified design, or more generally for complex experimental designs which induce dependence in the treatment assignments across individuals in the experiment. In Section (ref), we show that $\mathbb V_\ast$ continues to be the efficiency bound for estimating $\theta_0$ for a large class of treatment assignment mechanisms with a fixed marginal probability of treatment assignment, which includes highly stratified designs as a special case. We can thus conclude that, from the perspective of asymptotic efficiency, highly stratified designs are optimal experimental designs for a broad range of treatment effect estimation problems. In Section (ref) we build on this result and establish conditions under which efficient estimation of $\theta_0$ using the na\"ive method of moments estimator can be achieved only if the experimental design leads to fast-balancing of the treatment across covariate values. In this sense, we show that the fast-balancing property of highly stratified experiments is in fact a necessary condition for achieving efficient estimation of $\theta_0$ “by design.”

The Asymptotic Variance of Highly Stratified Experiments

In this section, we derive the asymptotic distribution of the method of moments estimator $\hat \theta_n$ when treatment is assigned by a highly stratified design over the baseline covariates $X^{(n)}$. This assignment mechanism uses the covariates $X^{(n)}$ to group units with similar covariate values into blocks of fixed size, and then assigns treatment completely at random within each block. In order to describe the assignment mechanism formally, we require some further notation to define the blocks of units. Let $\ell$ and $k$ be arbitrary positive integers with $\ell < k$ and set $\eta = \ell/k$. Here, $k$ is the total number of units in each block and $\ell$ is the number of treated units in each block. For simplicity, assume that $n$ is divisible by $k$. We then represent blocks of units using a partition of $\{1, \ldots, n\}$ given by \[\left\{\lambda_j = \lambda_j(X^{(n)}) \subseteq \{1, \ldots, n\}, 1 \le j \le n/k\right\}~,\] with $|\lambda_j| = k$. Because of its possible dependence on $X^{(n)}$, $\{\lambda_j: 1 \le j \le n/k\}$ encompasses a variety of different ways of blocking the $n$ units according to the observed, baseline covariates. We note, however, that our framework is not intended to reflect settings in which the blocks themselves are randomly sampled from the population of interest pashley2021insights. Given such a partition, we assume that treatment status is assigned as described in the following assumption:

assumptionTreatment status is assigned so that $(R^{(n)}(1), R^{(n)}(0)) \perp \!\!\! \perp A^{(n)} \big | X^{(n)}$ and, conditional on $X^{(n)}$, \[\{(A_i: i \in \lambda_j): 1 \le j \le n/k\}\] are i.i.d.\ and each uniformly distributed over all permutations of $(\underbrace{0, 0, \ldots, 0}_{k - \ell}, \underbrace{1, 1, \ldots, 1}_{\ell})$.

The assignment mechanism described in Assumptions (ref) generalizes the definition of a matched pairs design. In particular, we recover a matched pairs design if we set $(\ell, k) = (1, 2)$, with $\eta = 1/2$. Indeed, suppose $n$ is even and consider pairing the experimental units into $n / 2$ pairs, represented by the sets \[ \lambda_j = \{\pi(2j - 1), \pi(2j)\} \text{ for } j = 1, \ldots, n / 2~, \] where $\pi = \pi_n(X^{(n)})$ is a permutation of $n$ elements. Because of its possible dependence on $X^{(n)}$, $\pi$ encompasses a broad variety of ways of pairing the $n$ units according to the observed, baseline covariates $X^{(n)}$. Given such a $\pi$, we assume that treatment status is assigned so that Assumption (ref) holds and, conditional on $X^{(n)}$, $(A_{\pi(2j-1)}, A_{\pi(2j)}), j = 1, \ldots, n / 2$ are i.i.d.\ and each uniformly distributed over the values in $\{(0,1), (1,0)\}$. For some examples of such an assignment mechanism being used in practice, see, for instance, angrist2009effects, banerjee2015miracle, and bruhn2016impact.

remarkNote that Assumption (ref) generalizes matched pairs designs along two dimensions: first, it allows for treatment fractions other than $\eta = 1/2$. Second, it allows for choices of $\ell$ and $k$ which are not relatively prime. For instance, if we set $(\ell, k) = (2, 4)$, then $\eta = 1/2$ as in matched pairs, but now the assignment mechanism blocks units into groups of size $4$ and assigns two units to treatment, two units to control. Although Theorem (ref) below establishes that allowing for this level of flexibility has no effect on the asymptotic properties of our estimator, in our experience we have found that designs which employ these treatment “replicates” in each block can simplify the construction of variance estimators in practice; see Section (ref) for details, and imbens2011experimental for an early discussion.

Our analysis will require some discipline on the way in which the blocks are formed. In particular, we will require that the units in each block be close in terms of their baseline covariates in the sense described by the following assumption:

assumptionThe blocks used in determining treatment status satisfy \[ \frac{1}{n} \sum_{1 \leq j \leq n/k} \max_{i, i' \in \lambda_j} \|X_{i} - X_{i'}\|^2 \stackrel{P}{\rightarrow} 0~. \]

bai2022inference and cytrynbaum2023designing discuss blocking algorithms that satisfy Assumption (ref). When $X_i \in \mathbf R$ and $E_Q[X_i^2] < \infty$, a simple algorithm that satisfies Assumption (ref) is simply to order units from smallest to largest and then block adjacent units into blocks of size $k$. In the case of matched pairs, if $\mathrm{dim}(X_i) > 1$ and $E_Q[\|X_i\|^d] < \infty$ for $d \geq \mathrm{dim}(X_i) + 1$, then Assumption (ref) is satisfied by the nbpmatching algorithm in R that minimizes the sum of squared distances of $X$ within pairs: see Appendix A of bai2024inference-1 for details. Beyond the case of pairs, cytrynbaum2023designing demonstrates that the optimal blocking satisfies the following bound \[ \frac{1}{n} \sum_{1 \leq j \leq n/k} \max_{i, i' \in \lambda_j} \|X_{i} - X_{i'}\|^2 = O_P(n^{{2/d} - 2/(\rm{dim}(X_i)+1)})~, \] from which we can deduce that the rate of convergence depends on both the dimension of $X_i$ as well as the number of moments it possesses. Note further that Assumption (ref) can be satisfied even if $X_i$ contains discrete components, as may arise when summarizing a categorical variable numerically using, e.g., one-hot encoding.

Finally, we impose the following assumptions to derive the large-sample properties of $\hat{\theta}_n$. In what follows, when writing expectations and variances, we suppress the subscripts $P$ and $Q$ whenever doing so does not lead to confusion.

assumptionLet $m(\cdot) = (m_s(\cdot): 1 \le s \le d_{\theta})'$. The moment functions are such that \begin{enumerate}[\rm (a)] • For every $\epsilon > 0$, $\inf\limits_{\theta \in \Theta: \|\theta - \theta_0\| > \epsilon} \| E[m(X_i, A_i, R_i, \theta)] \| > 0$. • $E[m(X_i, A_i, R_i, \theta)]$ is differentiable at $\theta_0$ with a nonsingular derivative $M = \frac{\partial}{\partial \theta'} E[m(X, \allowbreak A, R, \theta)] \Big |_{\theta = \theta_0}$. • For $1 \leq s \leq d_\theta$, $E[((m_s(X, a, R(a), \theta) - m_s(X, a, R(a), \theta_0))^2] \to 0$ as $\theta \to \theta_0$ for $a \in \{0, 1\}$. • For $1 \leq s \leq d_\theta$, $\{m_s(x, a, r, \theta): \theta \in \Theta\}$ is pointwise measurable in the sense that there exists a countable set $\Theta^\ast$ such that for each $\theta \in \Theta$, there exists a sequence $\{\theta_m\} \subset \Theta^\ast$ such that $m_s(x, a, r, \theta_m) \to m_s(x, a, r, \theta)$ as $m \to \infty$ for all $x, a, r$. • (i) $\sup_{\theta \in \Theta} E[\|m(X, a, R(a), \theta)\|] < \infty$ for $a \in \{0, 1\}$. (ii) $\{m_s(x, 1, r, \theta): \theta \in \Theta^\ast\}$ and $\{m_s(x, 0, r, \theta): \theta \in \Theta^\ast\}$ are $Q$-Donsker for $1 \leq s \leq d_\theta$. • For $a \in \{0, 1\}$, $E[ m_s(X, a, R(a), \theta_0) | X = x]$ is $C$-Lipschitz for $1 \leq s \leq d_\theta$, for some constant $C < \infty$. \end{enumerate}

Assumption (ref)(a) is a standard assumption to ensure the solution to (ref) is “well separated.” It appears as a condition, for instance, in Theorem 5.9 in van_der_vaart1998asymptotic. Assumption (ref)(b) is a standard assumption used when deriving the properties of $Z$-estimators. See, for instance, Theorem 3.1 in newey1994large and Theorem 5.21 in van_der_vaart1998asymptotic. Because differentiability is imposed on their expectations instead of the moment functions themselves, the moment functions are allowed to be nonsmooth, as in Example (ref). Assumption (ref)(c) requires the moment function to be mean-square continuous in $\theta$. Assumption (ref)(d) is a standard condition to guarantee the measurability of the supremum of a suitable class of functions. In particular, it allows us to define expectations of suprema without invoking outer expectations. See Example 2.3.4 in van_der_vaart1996weak for details. Assumption (ref)(e) is a standard assumption which guarantees the existence of an integrable envelope and allows us to invoke a uniform law of large numbers and a uniform central limit theorem (see, for instance, page 81 of van_der_vaart1996weak for a definition of a Donsker class). In particular, this assumption can be verified for Examples (ref)--(ref). Assumption (ref)(f) is a common assumption which simplifies some arguments when studying highly stratified designs, and ensures units that are close in terms of the baseline covariates are also close in terms of their moments. Note that Assumption (ref)(f) could be dropped following the approximation arguments in Lemma C.5 of cytrynbaum2023designing; see also Examples (ref) and (ref) for further discussion.

The following theorem establishes the asymptotic variance of the na\"ive method of moments estimator when the treatment assignment mechanism is highly stratified in the sense of satisfying Assumptions (ref)--(ref). Its proof relies on a crucial technical result in han2021complex, which allows us to compare the empirical process that depends on the treatment assignments with the empirical process that only depends on i.i.d.\ quantities. As a consequence of us leveraging this result, Assumption (ref) is comparable to the typical assumptions imposed to study the properties of method of moments estimators using i.i.d.\ data.

theoremSuppose the treatment assignment mechanism satisfies Assumptions (ref)--(ref) and the moment functions satisfy Assumption (ref). Let $\hat{\theta}_n$ be defined as in (ref). Then, \begin{equation} \sqrt n(\hat \theta_n - \theta_0) = \frac{1}{\sqrt{n}} \sum_{1 \leq i \leq n} \psi^\ast(X_i, A_i, R_i, \theta_0) + o_P(1) . \end{equation} where \begin{align*} & \psi^\ast(X_i, A_i, R_i, \theta_0) \\ & = - M^{-1} \Big ( I\{A_i = 1\} (m(X_i, 1, R_i, \theta_0) - E[m(X_i, 1, R_i(1), \theta_0) | X_i]) \\ & + I\{A_i = 0\} (m(X_i, 0, R_i, \theta_0) - E[m(X_i, 0, R_i(0), \theta_0) | X_i]) \\ & + \eta E[m(X_i, 1, R_i(1), \theta_0) | X_i] + (1 - \eta) E[m(X_i, 0, R_i(0), \theta_0) | X_i] \Big ) . \end{align*} Further, we have that \begin{equation} \sqrt n(\hat \theta_n - \theta_0) \stackrel{d}{\to} N(0, \mathbb V_\ast) , \end{equation} where \begin{equation} \mathbb V_\ast = \operatorname{Var}[\psi^\ast(X_i, A_i, R_i, \theta_0)] . \end{equation}

In order to make Theorem (ref) useful for inference about $\theta_0$, we describe in Section (ref) an estimator $\hat{\mathbb{V}}_n$ of $\mathbb{V}_*$. We now sketch an argument of the proof of Theorem (ref) to highlight the fundamental role played by a “fast-balancing” property of highly stratified designs (see (ref) below). In the proof of Theorem (ref), we first establish that

equation[equation omitted — 147 chars of source]

To further establish (ref), it thus suffices to show that, under a highly stratified design,

equation[equation omitted — 194 chars of source]

To obtain this equivalence, consider the following decomposition of $m(\cdot)$:

align*[align* omitted — 450 chars of source]

Then the equivalence follows if we can show that

equation[equation omitted — 166 chars of source]

We call (ref) the fast-balancing condition for the function $m(\cdot)$. Intuitively, the fast-balancing condition imposes that the experimental design should balance the treatment across covariate values at a rate which is faster than sampling variation. To see why (ref) holds for a highly stratified design, let $\Omega(X_i) = E[m(X_i, 1,R_i(1),\theta_0) - m(X_i, 0,R_i(0),\theta_0)|X_i]$ and first note that, by Assumption (ref), \[E\bigg[\frac{1}{\sqrt{n}}\sum_{1 \le i \le n}(A_i - \eta)\Omega(X_i) \bigg \vert X^{(n)}\bigg] = 0~.\] Next, for $1 \leq s \leq d_\theta$, let $\Omega^{(s)}(X_i)$ denote the $s$th component of $\Omega(X_i)$. Then it can be shown using Assumption (ref) and (ref)(f) that for $1 \leq s \leq d_\theta$, \[\operatorname{Var}\bigg[\frac{1}{\sqrt{n}}\sum_{1 \le i \le n}(A_i - \eta)\Omega^{(s)}(X_i) \bigg \vert X^{(n)}\bigg] \le C^2\frac{\ell(k - \ell)}{k - 1}\left(\frac{1}{n}\sum_{1 \le j \le n/k} \max_{i, i' \in \lambda_j} \|X_i - X_{i'}\|^2\right)~,\] where $C$ denotes the Lipschitz constant in Assumption (ref)(f), and so the conditional variance converges in probability to zero under Assumption (ref). The fast-balancing condition (ref) then follows by an application of Chebyshev's inequality conditional on $X^{(n)}$ and the dominated convergence theorem. In Section (ref), we further argue that the fast-balancing condition (ref) is in fact a necessary condition which a given experimental design must satisfy to ensure (ref).

remarkNote it follows from (ref) that \begin{equation} \eta E_Q[m(X_i, 1, R_i(1), \theta_0)] + (1 - \eta) E_Q[m(X_i, 0, R_i(0), \theta_0)] = E_P[m(X_i, A_i, R_i, \theta_0)] = 0 , \end{equation} so that $E[\psi^\ast(X_i, A_i, R_i, \theta_0)] = 0$. It is further straightforward to show using Assumption (ref) that \begin{align} \mathbb V_\ast & = \operatorname{Var}[\psi^\ast(X_i, A_i,R_i, \theta_0)] \\ \nonumber & = M^{-1} \big ( E \big [ \eta \operatorname{Var}[m(X_i, 1, R_i(1), \theta_0) | X_i] + (1 - \eta) \operatorname{Var}[m(X_i, 0, R_i(0), \theta_0) | X_i] \big ] \\ \nonumber & + \operatorname{Var} \big [ \eta E[m(X_i, 1, R_i(1), \theta_0) | X_i] + (1 - \eta) E[m(X_i, 0, R_i(0), \theta_0) | X_i] \big ] \big ) (M^{-1})' \end{align} For instance, in the special case of the ATE (Example (ref)) we obtain that \begin{align} \nonumber \operatorname{Var}[\psi^\ast(X_i, A_i,R_i, \theta_0)] & = E\left[\frac{\operatorname{Var}[Y_i(1)|X_i]}{\eta} + \frac{\operatorname{Var}[Y_i(0)|X_i]}{1 - \eta} \right. \\ & + \left. \left(E[Y_i(1)- Y_i(0)|X_i] - E[Y_i(1) - Y_i(0)]\right)^2\right] , \end{align} which matches the asymptotic variance derived in bai2022inference for matched pairs. Theorem (ref) however accommodates a much larger class of parameters, including those introduced in Examples (ref)--(ref).
remarkBy comparing the variance expression in (ref) to the variance expression for $\mathbb{V}_*$, we obtain \begin{equation} \mathbb{V} - \mathbb{V}_* = \eta (1 - \eta) M^{-1} \operatorname{Var}[E[m(X_i, 1, R_i(1), \theta_0) - m(X_i, 0, R_i(0), \theta_0) | X_i]](M^{-1})^{\prime} , \end{equation} which is positive semidefinite. From this, we conclude that the asymptotic variance of the naive method of moments estimator $\hat{\theta}_n$ is lower in a highly stratified design compared to i.i.d.\ assignment. In Section (ref), we will further show that $\mathbb V_\ast$ is the lowest possible asymptotic variance among regular estimators for $\theta_0$ in a large class of treatment assignment mechanisms, including both i.i.d.\ assignment and highly stratified designs. When $d_\theta = 1$, we may express $\mathbb V - \mathbb V_\ast$ in terms of the “nonparametric $R^2$.” In particular, $\mathbb V - \mathbb V_\ast$ is proportional to $E[R_{g, X}^2(g_i, X_i) \operatorname{Var}[g_i]]$, where \begin{equation} R_{g, X}^2(g_i, X_i) = \frac{\operatorname{Var}[E[g_i|X_i]]}{\operatorname{Var}[g_i]} , \end{equation} and $g_i = m(X_i, 1, R_i(1), \theta_0) - m(X_i, 0, R_i(0), \theta_0)$. The quantity in (ref) measures how much of the variation in $g_i$ can be explained nonparametrically by $X_i$. See, for instance, chernozhukov2024long.
remarkNote that highly stratified designs include as a special case completely randomized designs, i.e., experiments in which a fixed proportion of the entire sample is assigned to treatment uniformly at random. To see this, consider, for instance, simply matching on an exogenously generated covariate. In this special case, we find from (ref) that $\mathbb{V}_{\ast} = \mathbb{V}$ as defined in Section (ref). In this way, we see that completely randomized experiments are asymptotically no more efficient than i.i.d.\ assignment.

Variance Estimation

In this subsection, we provide a consistent variance estimator for the asymptotic variance $\mathbb V_\ast$ in (ref). We suppose that a consistent estimator $\widehat M_n$ for $M$ is available, i.e., $\widehat M_n \xrightarrow{P} M$. In examples where $m$ is differentiable in $\theta$, including Examples (ref) and (ref)--(ref), the analog principle suggests that a natural estimator for $M$ is given by \[ \widehat M_n = \frac{1}{n} \sum_{1 \leq i \leq n} \frac{\partial}{\partial \theta'} m(X_i, A_i, R_i, \theta)\bigg|_{\theta = \hat \theta_n}~. \] In examples including Example (ref) where $m$ is nonsmooth in $\theta$, $M$ may consist of components that require nonparametric estimators. See, for instance, jiang2021bootstrap.

It then suffices to construct a consistent estimator for the “meat" in (ref). To motivate such an estimator, consider the expression in (ref) when $d_\theta = 1$. By the law of total variance, this middle component equals $\Sigma_1 + \Sigma_2$, where

align*[align* omitted — 681 chars of source]

For $a \in \{0, 1\}$, define \[ \hat \mu_n(a) = \frac{1}{\eta_a n} \sum_{1 \leq i \leq n} I \{A_i = a\} m(X_i, A_i, R_i, \hat \theta_n)~, \] where $\eta_1 = \eta$ and $\eta_0 = 1 - \eta$. The analog principle suggests that a natural estimator for $\Sigma_1$ is

align*[align* omitted — 354 chars of source]

To estimate $\Sigma_2$, we first define

align*[align* omitted — 226 chars of source]

and $\hat \varsigma_n(0, 1)$ similarly. Next, define \[ \hat \varsigma_n(1, 1) =

cases\frac{k}{n} \sum\limits_{1 \leq j \leq n / k} \frac{1}{\binom{\ell}{2}} \sum\limits_{i < i' \in \lambda_j: A_i = A_{i'} = 1} m(X_i, A_i, R_i, \hat \theta_n) \\ \times m(X_{i'}, A_{i'}, R_{i'}, \hat \theta_n)' & if \ell > 1 \\ \frac{2k}{n} \sum\limits_{1 \leq j \leq \frac{n}{2k}} {\sum\limits_{i \in \lambda_{2j}, i' \in \lambda_{2j - 1}: A_i = A_{i'} = 1}} m(X_i, A_i, R_i, \hat \theta_n) \\ \times m(X_{i'}, A_{i'}, R_{i'}, \hat \theta_n)' & if \ell = 1 .

\] Similarly, define \[ \hat \varsigma_n(0, 0) =

cases\frac{k}{n} \sum\limits_{1 \leq j \leq n / k} \frac{1}{\binom{k - \ell}{2}} \sum\limits_{i < i' \in \lambda_j: A_i = A_{i'} = 0} m(X_i, A_i, R_i, \hat \theta_n) \\ \times m(X_{i'}, A_{i'}, R_{i'}, \hat \theta_n)' & if k - \ell > 1 \\ \frac{2k}{n} \sum\limits_{1 \leq j \leq \frac{n}{2k}} {\sum\limits_{i \in \lambda_{2j}, i' \in \lambda_{2j - 1}: A_i = A_{i'} = 0}} m(X_i, A_i, R_i, \hat \theta_n) \\ \times m(X_{i'}, A_{i'}, R_{i'}, \hat \theta_n)' & if k - \ell = 1 .

\] Finally, define \[ \hat \Sigma_{2, n} = - \eta(1 - \eta) \big ( \hat \varsigma_n(1, 1) + \hat \varsigma_n(0, 0) - \hat \varsigma_n(1, 0) - \hat \varsigma_n(0, 1) - (\hat \mu_n(1) - \hat \mu_n(0)) (\hat \mu_n(1) - \hat \mu_n(0))' \big )~. \] The estimator $\hat \varsigma_n(1, 1)$ is constructed in one of two ways depending on the number of treated units in each block. If more than one unit in each block is treated, then we take the averages of all pairwise products of the treated units in each block, and average them across all blocks. We call this a “within block" estimator. If instead only one unit in each block is treated, then we take the product of two treated units in adjacent blocks. We call this a “between block" estimator, and note that similar constructions have been used previously in abadie2008estimation, bai2022inference, bai2024inference-1, and cytrynbaum2023designing. The estimator $\hat \varsigma_n(0, 0)$ is constructed similarly. A natural estimator for $\mathbb{V}_\ast$ is then given by \[\hat{\mathbb{V}}_n = \widehat{M}_n^{-1} \left(\hat{\Sigma}_{1,n} + \hat{\Sigma}_{2,n}\right) \left (\widehat{M}_n^{-1} \right )'~. \] Note that for specific choices of $m(\cdot)$, $\hat{\mathbb{V}}_n$ recovers estimators which have been studied in prior work on inference in highly stratified experiments. For instance, in the case of matched pairs with $(\ell, k) = (1,2)$ and $m(\cdot)$ as in Example (ref), so that $\theta_0$ is the ATE, $\hat{\mathbb{V}}_n$ exactly coincides with the estimator defined in equation (28) of bai2022inference.

In addition to Assumption (ref), we will now also require that the distances between units in adjacent blocks be “close" in terms of their baseline covariates:

assumptionThe blocks used in determining treatment status satisfy \[ \frac{1}{n} \sum_{1 \leq j \leq \lfloor n / 2 \rfloor} \max_{i \in \lambda_{2j - 1}, i' \in \lambda_{2j}} \|X_i - X_{i'}\|^2 \stackrel{P}{\to} 0~. \]

Note that given blocks which satisfy Assumption (ref), it is always possible to re-order the blocks such that the pairs $\{(\lambda_{2j-1}, \lambda_{2j})\}_{1 \le j \le \lfloor n/2 \rfloor}$ satisfy Assumption (ref), as long as we maintain the sufficient condition that $E[\|X_i\|^d] < \infty$ for $d \geq \mathrm{dim}(X_i) + 1$. This property could be achieved, for instance, by applying the {\tt nbpmatching} algorithm to the block-means $\{\bar{X}_j\}_{1 \le j \le n/k}$, where $\bar{X}_j := \frac{1}{k}\sum_{i \in \lambda_j}X_i$: see Lemma A.6 in cytrynbaum2023covariate for details.

To formally establish the consistency of $\hat{\mathbb V}_n$, we impose the following mild uniform integrability condition, as well as a Glivenko-Cantelli and uniform Lipschitz condition. All three conditions need only hold in an arbitrarily small neighborhood of $\theta_0$.

assumptionThere exists $\delta > 0$ such that \begin{enumerate}[(a)] • For $a \in \{0, 1\}$, \begin{equation*} \lim_{\lambda \to \infty} E \left [ \sup_{\theta \in \Theta: \|\theta - \theta_0\| < \delta} \|m(X_i, a, R_i(a), \theta)\|^2 I \left \{ \sup_{\theta \in \Theta} \|m(X_i, a, R_i(a), \theta)\| > \lambda \right \} \right ] = 0 . \end{equation*} • $\{E[m_s(X_i, a, R_i(a), \theta) | X_i = x]: \|\theta - \theta_0\| < \delta\}$ and $\{E[m_s(X_i, a, R_i(a), \theta) m(X_i, a, \allowbreak R_i(a), \theta)'| X_i = x]: \|\theta - \theta_0\| < \delta\}$ are $Q$-Glivenko Cantelli for $1 \leq s \leq d_\theta$. • For $a \in \{0, 1\}$, each component of $E[m(X, a, R(a), \theta) | X = x]$ and $E[m(X, a, R(a), \theta) \allowbreak m(X, a, R(a), \theta)' | X = x]$ is Lipschitz with a common Lipschitz constant across $\{\theta \in \Theta: \|\theta - \theta_0\| < \delta\}$. \end{enumerate}

Assumption (ref)(a) is a mild uniform integrability condition for the envelope function of the moment function in an arbitrarily small neighborhood of $\theta_0$. Assumption (ref)(b) is a mild condition that requires the conditional expectation of the moment functions to satisfy a uniform law of large numbers. See page 81 of van_der_vaart1996weak for a definition of a Glivenko-Cantelli class of functions. Assumption (ref)(c) strengthens Assumption (ref)(f) to hold uniformly in an arbitrarily small neighborhood of $\theta_0$.

The following theorem establishes the consistency of $\hat{\mathbb V}_n$ for $\mathbb V_\ast$:

theoremSuppose the treatment assignment mechanism satisfies Assumptions (ref), (ref), and (ref) and the moment functions satisfy Assumptions (ref) and (ref). Further suppose $\widehat M_n \xrightarrow{P} M$. Then, $\hat{\mathbb V}_n \xrightarrow{P} \mathbb V_\ast$.

An Efficiency Bound and the Necessity of “Fast-Balancing”

In Section (ref) we establish that $\mathbb{V}_{\ast}$ is the efficiency bound for a large class of experimental designs. As a consequence, we can conclude that highly stratified designs are asymptotically efficient “by design.” Building on this result, in Section (ref) we establish that a necessary condition for achieving the bound $\mathbb{V}_{\ast}$ when estimating $\theta_0$ using the na\"ive method of moments estimator is that the experimental design be fast-balancing, in the sense of (ref).

Efficiency Bound

An inspection of the asymptotic variance in (ref) reveals that $\mathbb V_\ast$ in fact coincides with the classical efficiency bound for estimating $\theta_0$ with i.i.d.\ assignment. For example, the variance derived in (ref) coincides with the efficiency bound derived in hahn1998role for estimating the ATE with a known marginal treatment probability $\eta$. Therefore, another way to interpret our result in Theorem (ref) is that the standard i.i.d.\ efficiency bound can be attained by a na\"ive method of moments estimator under a highly stratified design. On the other hand, because treatment status is not independent in a highly stratified design, a natural follow-up question is whether or not the efficiency bound for estimating $\theta_0$ changes relative to what can be obtained under i.i.d.\ assignment once we allow for more general assignment mechanisms. In this section, we show that $\mathbb{V}_*$ continues to be the efficiency bound for the class of parameters introduced in Section (ref), while allowing for a more general class of treatment assignment mechanisms. The main restriction on treatment assignment is given by Assumption (ref), which requires the marginal treatment probability to be known and equal to $\eta$. As mentioned earlier and explained in Remark (ref) below, it is possible to relax this requirement so that $\eta$ can be replaced by a known function $\eta(X_i)$. For a discussion of how our efficiency bound compares with other results in the literature, see Remark (ref).

We impose the following high-level assumption on the assignment mechanism:

assumptionThe treatment assignment mechanism is such that for any integrable function $\gamma: \mathbf R^{d_x} \to \mathbf R$, \[ \frac{1}{n} \sum_{1 \leq i \leq n} A_i \gamma(X_i) \stackrel{P}{\to} \eta E[\gamma(X_i)]~. \]

In words, Assumption (ref) requires that the assignment mechanism admits a law of large numbers for integrable functions of the covariate values. Examples (ref)--(ref) illustrate that the assumption holds for common treatment assignment mechanisms used in practice.

example[i.i.d.\ assignment] Let $A^{(n)}$ be assigned i.i.d., independently of $X^{(n)}$, such that $P\{A_i = 1\} = \eta$. Then it follows immediately by the law of large numbers that Assumption (ref) is satisfied.
example[Covariate-adaptive randomization (CAR)] Let $S: \mathbf R^{d_x} \to \mathcal S = \{1, \ldots, \allowbreak |\mathcal S|\}$ be a function that maps the covariates into a fixed, finite set of discrete strata. We call such a stratification “coarse", to distinguish it from highly stratified designs as defined in Section (ref). Define $S_i = S(X_i)$ and assume that treatment status is assigned so that \[ (R^{(n)}(1), R^{(n)}(0), X^{(n)}) \perp \!\!\! \perp A^{(n)} \big | S^{(n)}~, \] and that for $s \in \mathcal S$, \[ \frac{\sum_{1 \leq i \leq n} I \{S_i = s, A_i = 1\}}{\sum_{1 \leq i \leq n} I \{S_i = s\}} \stackrel{P}{\to} \eta~. \] This high-level assumption accommodates a large class of stratified assignment mechanisms, including stratified biased coin designs efron1971forcing, wei1978adaptive, minimization methods pocock1975sequential,hu2012asymptotic and stratified block randomization zelen1974randomization. It follows from Lemma C.4 in bugni2019inference that for any integrable function $\gamma(\cdot)$, \[ \frac{1}{n} \sum_{1 \leq i \leq n} A_i \gamma(X_i) \stackrel{P}{\to} \eta \sum_{s \in \mathcal S} P \{S_i = s\} E[\gamma(X_i)| S_i = s] = \eta E[\gamma(X_i)]~. \] Therefore, Assumption (ref) is satisfied.
example[CAR with general covariate features] ma2024new propose a family of covariate adaptive randomization procedures which assign treatment sequentially based on an imbalance metric defined by “feature maps” of (potentially continuous) covariates. It follows by Theorem 3.5 of their paper that Assumption (ref) is satisfied under appropriate conditions.
example[Matched pairs] Suppose $n$ is even and we assign treatment using a highly stratified design with $(\ell, k) = (1,2)$. As discussed at the beginning of Section (ref), such a design is also known as a matched pairs design. Assume that the pairing algorithm $\pi_n(X^{(n)})$ results in pairs that are close in the sense of Assumption (ref). It then follows from the same argument used to establish (ref) in the discussion following Theorem (ref) that for any Lipschitz integrable function $\gamma(\cdot)$, \[ \frac{1}{n} \sum_{1 \leq i \leq n} A_i \gamma(X_i) \stackrel{P}{\to} \frac{1}{2} E[\gamma(X_i)]~. \] By approximating integrable functions by Lipschitz integrable functions as in Lemma A.1 in hanneke2021universal, it can be shown that the convergence holds for any integrable function $\gamma(\cdot)$. Therefore, Assumption (ref) is satisfied.
example[Re-randomization] Re-randomization is an assignment mechanism in which researchers specify a balance criterion for the covariates, and then repeatedly generate assignments using a completely randomized design until an assignment is found which achieves an acceptable covariate distribution according to the balance criterion. The properties of re-randomization procedures have been studied in li2017general,li2020rerandomization, li2018asymptotic,li2020rerandomization-1, and cytrynbaum2024finely. It follows from Corollary 3.7 in cytrynbaum2024finely that Assumption (ref) holds for re-randomization designs, under appropriate assumptions.

We now present an efficiency bound for the parameter $\theta_0$ introduced in Section (ref). Formally, we characterize the bound via a convolution theorem that applies to all regular estimators of the parameter $\theta_0$. Following the definition on page 365 of van_der_vaart1998asymptotic, by a regular estimator, we mean an estimator whose asymptotic distribution is invariant to “local" perturbations of the data generating process: \[ \sqrt n(\tilde \theta_n - \theta(P_{t/\sqrt n, g})) \xrightarrow{P_{t/\sqrt n, g}} L ~, \] where $P_{t/\sqrt n,g}$ represents a “local" perturbation of the distribution $P$ along a “path” with “score" $g$. We leave the precise definition of regularity and related assumptions to Supplement A.2. In the paragraph following the statement of the theorem we provide some more details on the nature of our result.

theoremSuppose Assumptions (ref), (ref)(b), and (ref) hold, as well as Condition A.1 described in Supplement A.2. Further suppose $\mathbb V_\ast < \infty$. Let $\tilde{\theta}_n$ be any regular estimator of the parameter $\theta_0$ in the sense of (S.16) in Supplement A.2. Then, \[\sqrt{n}(\tilde{\theta}_n - \theta_0) \xrightarrow{d} L~,\] where \[L = N(0, \mathbb V_\ast) \ast B~, \] for $\mathbb V_\ast$ in (ref) and some fixed probability measure $B$ which is specific to the estimator $\tilde{\theta}_n$.

Given Theorem (ref) we call $\mathbb V_\ast = \operatorname{Var}[\psi^\ast(X_i, A_i,R_i, \theta_0)]$ the efficiency bound for $\theta_0$, since our result shows that this is the lowest asymptotic variance attainable by any regular estimator under our assumptions. Indeed, it follows from Anderson's lemma van_der_vaart1998asymptotic that the asymptotic loss of any regular estimator is bounded below by the loss under $N(0, \mathbb V_\ast)$ for any “bowl-shaped” loss function (including, in particular, square loss). We note that our assumptions on the assignment mechanism preclude us from immediately appealing to standard semi-parametric convolution theorems van_der_vaart1989asymptotic. Instead, we proceed by justifying an application of Theorem 3.1 in armstrong2022asymptotic combined with the convolution Theorem 3.11.2 in van_der_vaart1996weak to each $d_\theta$-dimensional parametric submodel separately, and then arguing that the supremum over all such submodels is attained by $\operatorname{Var}[\psi^{\ast}]$. A key observation is that in order to apply Theorem 3.1 in armstrong2022asymptotic to argue that the likelihood ratio process is locally asymptotically normal, the conditional information needs to settle down in the limit, which is guaranteed as long as Assumption (ref) is satisfied.

remarkFollowing similar arguments as those in Remark (ref), we can deduce that our efficiency bound agrees with well-known bounds for common parameters (like those presented in Examples (ref)--(ref)) in the setting of i.i.d.\ assignment. For example, we have noted in the case of the ATE (Example (ref)) that (ref) matches the efficiency bound under i.i.d.\ assignment derived in hahn1998role. See rafi2023efficient and armstrong2022asymptotic for related results in the context of stratified and response-adaptive experiments. Straightforward calculation also implies that, for the quantile treatment effect (Example (ref)), the efficiency bound is given by \begin{multline*} E \bigg [ \frac{1}{\eta} \frac{F_1 \big (\theta_0(1) | X_i \big ) \big (1 - F_1 \big (\theta_0(1) | X_i \big ) \big )}{f_1 \big (\theta_0(1) \big )^2} + \frac{1}{1- \eta} \frac{F_0 \big (\theta_0(0) | X_i \big ) \big (1 - F_0 \big (\theta_0(0) | X_i \big ) \big )}{f_0 \big (\theta_0(0) \big )^2} \\ + \bigg ( \frac{F_1 \big (\theta_0(1) | X_i \big ) - \tau}{f_1 \big (\theta_0(1) \big )} - \frac{F_0 \big (\theta_0(0) | X_i \big ) - \tau}{f_0 \big (\theta_0(0) \big )} \bigg)^2 \bigg ] , \end{multline*} which matches the efficiency bound under i.i.d.\ assignment derived in firpo2007efficient when the propensity score is set to $\eta$.
remarkThe efficiency bound in Theorem (ref) is attained by highly stratified experiments as in Theorem (ref) if no additional covariates are available for estimation beyond the set of covariates $X_i$ used in the design. In practice, researchers may consider adjusting for additional baseline covariates in order to improve efficiency. Suppose additional covariates $W^{(n)}$ are available and Assumption (ref) is modified such that \[ (R^{(n)}(1), R^{(n)}(0), W^{(n)}) \perp \!\!\! \perp A^{(n)} \big | X^{(n)}~. \] It can be shown that the efficiency bound, allowing for additional covariate adjustment based on $X_i$ and $W_i$, is \begin{equation} \mathbb V_\ast - \eta(1 - \eta) M^{-1} E[\operatorname{Var}[E[g_i|X_i,W_i]|X_i]] (M^{-1})' , \end{equation} where $g_i = m(X_i, 1, R_i(1), \theta_0) - m(X_i, 0, R_i(0), \theta_0)$. Then, as in Remark (ref), the potential gain in efficiency from exploiting $W_i$ in addition to $X_i$ is proportional to \[ E[R_{g, X, W}^2(g_i, X_i, W_i) \operatorname{Var}[g_i | X_i]]~, \] where \begin{equation*} R_{g, X, W}^2(g_i, X_i, W_i) = \frac{\operatorname{Var}[E[g_i|X_i,W_i]|X_i]}{\operatorname{Var}[g_i | X_i]} \end{equation*} is the nonparametric $R^2$ from regressing $g_i$ on $X_i$ and $W_i$ conditional on $X_i$. As a result, the scope for improving efficiency by adjusting for additional covariates is limited if $R_{g, X, W}^2$ is small. In the case of estimating the ATE, \[ g_i = \frac{Y_i(1)}{\eta} + \frac{Y_i(0)}{1 - \eta}~, \] so the scope for improvement depends on how much additional variation in the weighted potential outcomes can be explained by $W_i$ beyond $X_i$.
remarkHere, we comment on how Theorem (ref) relates to prior efficiency bounds in experiments with general assignment mechanisms. For the case of estimating the ATE, armstrong2022asymptotic derives an efficiency bound over a very large class of assignment mechanisms, including even response-adaptive designs, and shows that the bound is attained when units are assigned to treatment (control) with conditional probability proportional to the conditional standard deviation of the potential outcome under treatment (control). This type of assignment is sometimes referred to as the Neyman allocation. On the other hand, our results show that this bound may be quite loose whenever the assignment proportions are restricted to be anything not equal to the Neyman allocation, which is, of course, unknown. For example, the bound is not informative about what can be achieved if the assignment proportions were set to one half regardless of whether or not the conditional outcome variances across treatment and control are equal. Such settings frequently arise in practice due to logistical constraints or the absence of pilot data with which to estimate conditional variances of potential outcomes under treatment and control. Furthermore, as argued in cai2022performance, even if pilot data is available, these quantities may be estimated so poorly that exogenously constraining the assignment proportions to one half leads to more efficient estimates of the ATE in practice. Motivated by such concerns, rafi2023efficient derives an efficiency bound for the ATE over the class of “coarsely-stratified” assignment mechanisms studied in bugni2019inference, where the stratum-level assignment proportions are restricted {\it a priori} by the experimenter. This framework, however, rules out highly stratified designs. Finally, we once again emphasize that our analysis, unlike these other papers, applies to a general class of treatment effect parameters, including the ATE as a special case.
remarkAlthough we focus on the case where $\eta_i(X_i) = P\{A_i = 1|X_i\} = \eta$ is a constant, the proof of Theorem (ref) holds when $\eta_i(x) = \eta(x)$ for $1 \leq i \leq n$, where $\eta(x)$ is an arbitrary known and fixed function. In these settings, Lemma A.5 shows that the efficiency bound equals \begin{equation} \begin{split} \mathbb V_\ast & = \operatorname{Var}[\psi^\ast(X_i, A_i,R_i, \theta_0)] \\ & = M^{-1} \big ( E \big [ \eta(X_i) \operatorname{Var}[m(X_i, 1, R_i(1), \theta_0) | X_i] \\ & + (1 - \eta(X_i)) \operatorname{Var}[m(X_i, 0, R_i(0), \theta_0) | X_i] \big ] \\ & + \operatorname{Var} \big [ \eta(X_i) E[m(X_i, 1, R_i(1), \theta_0) | X_i] \\ & + (1 - \eta(X_i)) E[m(X_i, 0, R_i(0), \theta_0) | X_i] \big ] \big ) (M^{-1})' , \end{split} \end{equation} so that the only difference from (ref) is that $\eta$ is replaced by $\eta(X_i)$. If we additionally impose that $\eta(X_i) $ takes on a finite set of values $\{\eta_1, \dots, \eta_S\}$, then this bound could be achieved by separately implementing a highly stratified experiment over each set $\{i: \eta(X_i) = \eta_s\}$ for $1 \leq s \leq S$. In other words, separately within each stratum defined by the units for which $\eta(X_i) = \eta_s$, employ the assignment mechanism described in Assumptions (ref)--(ref) with $\ell/k = \eta_s$. For more general functions $\eta(\cdot)$, we conjecture one could employ the local randomization procedure in cytrynbaum2023designing.

The Necessity of “Fast-Balancing”

In this subsection, we provide conditions under which the fast-balancing condition described in (ref) is a necessary condition for efficient estimation of $\theta_0$ “by design.” As a supplement, we also provide necessary and sufficient conditions for an asymptotically linear estimator to be regular (in the sense of (S.16) in Supplement A.2) for a large class of treatment assignment mechanisms. Concretely, given an assignment mechanism, suppose $\tilde \theta_n$ is an asymptotically linear estimator for $\theta_0$ in the sense that

equation[equation omitted — 154 chars of source]

where $E[\psi(X_i, A_i, R_i, \theta_0)] = 0$ and $\operatorname{Var}[\psi(X_i, A_i, R_i, \theta_0)] < \infty$. The results in this section derive necessary and sufficient conditions for $\tilde \theta_n$ to be regular, and further demonstrate that in order for $\tilde{\theta}_n$ to be regular and efficient, either $\psi = \psi^\ast$ or a fast-balancing condition involving $\psi(\cdot)$ needs to be satisfied. Therefore, highly stratified designs are not only sufficient to guarantee efficiency when estimating $\theta_0$ using the na\"ive method of moments estimator, but their fast-balancing property is also necessary.

In order to study the behavior of the estimator under local alternatives, we will impose the following high-level assumption on the “imbalance” of the treatment assignments. To describe the assumption, let $\rho$ denote any metric that metrizes weak convergence.

assumptionThe treatment assignment mechanism is such that for any square-integrable function $\gamma: \mathbf R^{d_x} \to \mathbf R^{d_\theta}$ with $E[\gamma(X_i)] = 0$ , \[ \rho \bigg ( \mathcal L \Big ( \frac{1}{\sqrt n} \sum_{1 \leq i \leq n} (A_i - \eta) \gamma(X_i) \Big \vert X^{(n)} \Big ), ~N(0, V_\gamma^{\rm imb}) \bigg ) \xrightarrow{P} 0 \] for some deterministic variance $V_\gamma^{\rm imb}$, where $\mathcal L(\cdot | X^{(n)})$ denotes the conditional distribution given $X^{(n)}$.

In the following examples, we discuss Assumption (ref) in the context of some common treatment assignment mechanisms.

exampleRevisiting Example (ref), let $A^{(n)}$ be assigned i.i.d., independently of $X^{(n)}$, such that $P \{A_i = 1\} = \eta$. Then, by verifying the conditions of the Lindeberg-Feller CLT conditional on $X^{(n)}$, it can be shown that Assumption (ref) is satisfied with $V^{\rm imb}_\gamma = \eta (1 - \eta) \operatorname{Var}[\gamma(X_i)]$.
exampleRevisiting Example (ref), suppose treatment status is assigned using stratified block randomization, which is a special case of covariate-adaptive randomization where $A^{(n)}$ is such that \[\sum_{1 \leq i \leq n} A_i I \{S_i = s\} = \bigg\lfloor \eta \sum_{1 \leq i \leq n} I \{S_i = s\} \bigg\rfloor~,\] with all such assignments being drawn uniformly at random and independently across strata. It follows from Theorem 12.2.1 in lehmann2022testing combined with a subsequencing argument that Assumption (ref) is satisfied with $V_\gamma^{\rm imb} = \eta (1 - \eta) E[\operatorname{Var}[\gamma(X_i) | S_i]]$.
exampleRevisiting Example (ref), suppose treatment is assigned using the covariate adaptive randomization procedure described in ma2024new. Then it follows from Theorem 3.6 in their paper that, under appropriate assumptions, \[\frac{1}{\sqrt{n}}\sum_{1 \le i \le n}(A_i - \eta)\gamma(X_i) \xrightarrow{d} N(0, \tilde{V})~,\] for some variance $\tilde{V}$. Note, however, that this result is not conditional on $X^{(n)}$ and thus does not immediately imply Assumption (ref). We conjecture that a similar result could be established conditional on $X^{(n)}$ and thus Assumption (ref) would be satisfied.
exampleRevisiting Example (ref), suppose $n$ is even and we assign treatment using a matched pairs design. It then follows by arguing as in the discussion following Theorem (ref) that for any square-integrable Lipschitz function $\gamma(\cdot)$, \[ \operatorname{Var}\bigg[\frac{1}{\sqrt{n}} \sum_{1 \leq i \leq n}(A_i - \eta) \gamma(X_i)\bigg\vert X^{(n)}\bigg] \stackrel{P}{\to} 0~. \] Therefore, by Markov's inequality, Assumption (ref) is satisfied with $V_\gamma^{\rm imb} = 0$. By approximating square-integrable functions by square-integrable Lipschitz functions as in Lemma C.5 in cytrynbaum2023designing, it can be shown that the convergence holds for any square-integrable function $\gamma(\cdot)$.
exampleRevisiting Example (ref), we note that, following Corollary 3.7 in cytrynbaum2024finely, we do not expect Assumption (ref) to hold for re-randomization designs in general.

We are now ready to state a theorem that characterizes all regular asymptotically linear estimators and establishes the necessity of the fast-balancing condition for efficient estimation using $\hat{\theta}_n$.

theoremSuppose the treatment assignment mechanism satisfies Assumptions (ref) and (ref)--(ref). Suppose $\tilde \theta_n$ is an asymptotically linear estimator for $\theta_0$ in the sense of (ref). Then, $\tilde \theta_n$ is regular if and only if \[ \psi(x, a, r, \theta_0) = \psi^\ast(x, a, r, \theta_0) + \psi^\perp(x, a, \theta_0)~, \] for some function $\psi^\perp$ such that $E[\psi^\perp(X_i, A_i, \theta_0) | X_i] = \eta \psi^\perp(X_i, 1, \theta_0) + (1 - \eta) \psi^\perp(X_i, 0, \allowbreak \theta_0) = 0$. Furthermore, if $\tilde \theta_n$ is regular, it attains the efficiency bound if and only if \begin{equation} \frac{1}{\sqrt n} \sum_{1 \leq i \leq n} (A_i - \eta) E[\psi(X_i, 1, R_i(1), \theta_0) - \psi(X_i, 0, R_i(0), \theta_0) | X_i] = o_P(1) . \end{equation}

To further understand condition (ref), note that $E[\psi^\ast(X_i, 1, R_i(1)) | X_i] = E[\psi^\ast(X_i, 0, \allowbreak R_i(0)) | X_i]$, so

align*[align* omitted — 346 chars of source]

where the last equality follows from the fact that $E[\psi^\perp(X_i, A_i, \theta_0) | X_i] = 0$. As a result, (ref) holds either when the treatment assignment mechanism groups units with similar values of $\psi^\perp(X_i, 1, \theta_0) - \psi^\perp(X_i, 0, \theta_0)$, or when the estimator is based on the efficient influence function, so that $\psi^\perp = 0$. Recall that as an intermediate step in the proof of Theorem (ref), we showed in (ref) that for the na\"ive method of moments estimator $\hat \theta_n$, $\psi(\cdot) = - M^{-1} m(\cdot)$, so (ref) coincides with the fast-balancing condition in (ref). We can thus conclude from Theorem (ref) that the fast-balancing condition is necessary to achieve efficient estimation based on the na\"ive method of moments estimator, when the class of assignment mechanisms satisfy Assumption (ref).

exampleRevisiting Example (ref), recall $\hat \theta_n$ estimates the ATE based on the moment conditions in (ref). Direct calculation shows that for $\hat \theta_n$, \[ \psi^\perp(x, a, \theta_0) = (a - \eta) \Big ( \frac{\mu_1(x)}{\eta} + \frac{\mu_0(x)}{1 - \eta} \Big )~. \] As a result, Theorem (ref) demonstrates that $\hat{\theta}_n$ does not achieve the efficiency bound, unless \[ \frac{1}{\sqrt n} \sum_{1 \leq i \leq n} (A_i - \eta) \left ( \frac{\mu_1(X_i)}{\eta} + \frac{\mu_0(X_i)}{1 - \eta} \right ) = o_P(1)~, \] which is indeed the case in highly stratified experiments when the treatment assignment mechanism satisfies Assumption (ref).

Simulations

In this section, we illustrate the theoretical results in Sections (ref) and (ref) through a simulation study. Throughout this section, we set $\eta = 1/2$, and compare the mean-squared error (MSE), bias, the length of the confidence interval and its coverage rate for the following combinations of treatment assignment mechanisms and estimators:

enumerate[--] • i.i.d.\ treatment assignment and the naïve method of moments estimator • i.i.d.\ treatment assignment and covariate adjusted estimators • Matched pairs, i.e., a highly stratified design with $(\ell, k) = (1, 2)$, and the naïve method of moments estimator • Matched quadruplets, i.e., a highly stratified design with $(\ell, k) = (2, 4)$, and the naïve method of moments estimator

In Section (ref), we present the model specifications and estimators for estimating the ATE as in Example (ref). Supplement B contains the model specifications and estimators for estimating the LATE as in Example (ref). Section (ref) reports the simulation results for the MSE. The tables and discussion for bias and coverage are contained in Supplement B.

Average Treatment Effect

In this section, we present model specifications and estimators for estimating the ATE as in Example (ref). The Recall that in this case the moment function we consider is given by

equation*[equation* omitted — 112 chars of source]

with $R_i = Y_i$. For $a \in \{0, 1\}$ and $1\leq i \leq n$, the potential outcomes are generated according to the equation:

equation[equation omitted — 76 chars of source]

where $\mu_0(X_i) = \sum_{1 \leq l \leq 8} w_l \big ( X_{i, l} + \frac{1}{3} (X_{i, l}^2 - 1) \big )$ for $w = (2, 1, 1, 0.5, 0.01, 0.001, 0.0001, \allowbreak 0.00001)$, $\mu_1(X_i) = 0.2 + \mu_0(X_i)$, $\epsilon_i \sim N(0, 4)$, $(X_i, \epsilon_i)$, $1 \leq i \leq n$ are i.i.d., and for each $1 \leq i \leq n$, $(X_i, \epsilon_i)$ are independent. In each of the simulations to follow, when forming pairs/quadruplets or performing regression adjustment, we use only a subvector of the covariates $X_i$ consisting of the first $T$ covariates, for $T \in \{2, 4, 8\}$.

We consider the following three estimators for the ATE:

description• \begin{equation*} \hat\theta_n^{\rm unadj} = \frac{1}{n / 2} \sum_{1\leq i \leq n} ( Y_i A_i - Y_i(1-A_i) ) . \end{equation*} • \begin{equation*} \hat\theta_n^{\rm adj, 1} = \frac{1}{n} \sum_{1 \leq i \leq n} \big ( 2 A_i(Y_i - \hat\mu_1^Y(X_i)) - 2(1- A_i)(Y_i - \hat\mu_0^Y(X_i)) + \hat\mu_1^Y(X_i) - \hat\mu_0^Y(X_i) \big ) , \end{equation*} where $\hat\mu_a^Y(X_i)$ is the linear projection of $Y_i$ on $(1, (X_{i, l}, X_{i, l}^2: 1 \leq l \leq T))$ in the subsample with $A_i=a$. • \begin{equation*} \hat\theta_n^{\rm adj, 2} = \frac{1}{n} \sum_{1 \leq i \leq n} \big ( 2 A_i(Y_i - \hat\mu_1^Y(X_i)) - 2(1- A_i)(Y_i - \hat\mu_0^Y(X_i)) + \hat\mu_1^Y(X_i) - \hat\mu_0^Y(X_i) \big ) , \end{equation*} where $\hat\mu_a^Y(X_i)$ is the linear projection of $Y_i$ on on $(1, (X_{i, l}, X_{i, l}^2, X_{i, l} \allowbreak I\{X_{i, l} > \hat t_l\}: 1 \leq l \leq T))$ in the subsample with $A_i=a$, where $\hat t_l$ is the sample median of $X_{i, l}$, $1 \leq i \leq n$.

The first estimator $\hat\theta_n^{\rm unadj}$ is the naïve method of moments estimator given by the solution to (ref). The second and third estimators $\hat\theta_n^{\rm adj,1}$ and $\hat\theta_n^{\rm adj,2}$ are covariate-adjusted estimators which can be obtained as two-step method of moments estimators from solving the “augmented” moment equation (ref) described in the discussion at the end of Section (ref). $\hat\theta_n^{\rm adj,1}$ and $\hat\theta_n^{\rm adj,2}$ differ in the choice of basis functions used in the construction of the estimators $\hat{\mu}_a(x)$. Note that by the double-robustness property of the augmented estimating equation (ref), it can be shown that the adjusted estimators $\hat{\theta}_n^{\rm adj,1}$, $\hat{\theta}_n^{\rm adj,2}$ are consistent and asymptotically normal regardless of the choice of estimators $\hat{\mu}_a(x)$, but consistency of $\hat{\mu}_a(x)$ to $\mu_a(x)$ would ensure that $\hat{\theta}_n^{\rm adj,1}$, $\hat{\theta}_n^{\rm adj,2}$ are efficient under i.i.d.\ assignment robins1995analysis, tsiatis2008covariate, chernozhukov2017doubledebiasedneyman.

Simulation Results

We focus on the MSE and leave the discussion for bias and coverage to Supplement B. Table (ref) displays the ratio of the empirical MSE for each design/estimator pair relative to the MSE of the unadjusted estimator under i.i.d.\ assignment, computed across $4000$ Monte Carlo replications. As expected given our theoretical results, we find that the empirical MSEs of the na\"ive unadjusted estimator under a matched pairs/quads design closely match the empirical MSEs of the covariate adjusted estimators under i.i.d.\ assignment; this feature is particularly noteworthy given that the adjusted estimators are in fact correctly specified, and the correct specification would be unknowable in practice. We note that the MSE improvement of the unadjusted estimator with a matched pairs/quads design relative to i.i.d.\ assignment is typically worse with 2 or 8 covariates than with 4 covariates: this phenomenon stems from the fact that the first 4 covariates are much stronger predictors of the control outcome than the last 4 covariates, which are almost uninformative. Although we have found in prior work bai2024inference that the matched pairs design delivers a lower MSE than the matched quads design when the potential outcomes depend on the covariates linearly, we do not find that this is the case here with a nonlinear model. However, we do consistently find that the MSE with 8 covariates is smaller for the matched pairs design than the matched quads design. This comparison illustrates that, although as discussed in Remark (ref), our theoretical results imply matched pairs and matched quads designs are not distinguishable asymptotically, their finite-sample properties may differ, especially when the number of covariates is large. In particular, note Assumption (ref) is more stringent for matched quads, for which $k = 4$, than matched pairs, for which $k = 2$.

table[table omitted — 3,190 chars of source]

Recommendations for Empirical Practice

We conclude with some recommendations for empirical practice based on our theoretical results. Overall, our findings highlight the general benefit of highly stratified designs for designing efficient experiments: highly stratified experiments “automatically” perform fully-efficient covariate adjustment for a large class of interesting parameters. This finding generalizes similar observations made by bai2022inference, bai2022optimality and cytrynbaum2023designing for the special case of estimating the ATE.

Our simulation evidence suggests, however, that highly stratified experiments may produce less precise estimates than (correctly specified) covariate adjustment when the the dimension of $X_i$ is large relative to the sample size. For this reason, we recommend that practitioners construct their blocks using a subset of the baseline covariates that they believe have the highest explanatory power in terms of the nonparametric $R^2$ in (ref); the pre-treatment measure of the outcomes of interest, for example, is typically believed to be one such covariate bruhn2009pursuit. The experimental data can then be analyzed efficiently using an unadjusted method-of-moments estimator.

If one wishes to perform covariate adjustment with additional covariates beyond those used for blocking, then this can be done {\it ex-post}. As discussed in Remark (ref), the scope for improvement from covariate adjustment is limited by the nonparametric $R^2$ from the regression of the moment functions on the additional covariates, conditional on the ones used for matching; if one has already matched on the covariates with the highest explanatory power, then the potential gain in efficiency from adjusting for these additional covariates may be limited. We further caution that care must be taken to ensure that the adjustment is performed in such a way that it guarantees a gain in efficiency: see bai2024covariate and cytrynbaum2023covariate for related discussion. Recent work has developed such methods of covariate adjustment for specific parameters of interest bai2024covariate, bai2024inference-1, bai2025inference,cytrynbaum2023covariate, but we leave the development of a method of covariate adjustment which applies at the level of generality considered in this paper to future work.