A New Design-Based Variance Estimator for Finely Stratified Experiments

A New Design-Based Variance Estimator for Finely Stratified Experiments

KEYWORDS: Experiments, Finite Population, Average Treatment Effect, Matched Pairs, Stratification

JEL classification codes: C12, C31, C35, C36

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

Introduction

This paper considers the problem of design-based inference on the average treatment effect in finely stratified experiments. Here, by “design-based” we mean that there is no sampling uncertainty and the only randomness stems from the variation in treatment assignment itself; by “finely stratified” we mean that treatment is assigned in the following fashion: units are first stratified into groups of a fixed size $k$ according to baseline covariates and then, within each group, a fixed number $\ell < k$ are assigned uniformly at random to treatment and the remainder to control. A prominent special case of this framework is a matched pairs design, in which $k = 2$ and $\ell = 1$. While our primary focus is on experiments where treatment is assigned at the level of the individual, in Remark (ref) we explain how our analysis can accommodate experiments where the treatment is assigned at the cluster level.

In this setting, we first establish, by way of motivation, a result that shows under mild conditions that inference using the usual difference-in-means estimator requires an estimator of its variance that is at least asymptotically upward-biased. After having motivated the importance of such estimators, we review the standard argument for why constructing such estimators is particularly challenging for the case in which $\ell = 1$ or $k - \ell = 1$; in particular, in this case the sample variance of the outcomes under treatment (if $\ell = 1$) or control (if $k - \ell = 1$) computed in each stratum are identically zero. With this in mind, we propose a novel estimator that remains well-defined even in this challenging case, and demonstrate that our estimator is upward-biased under minimal assumptions. The construction of our variance estimator involves a pairing of the strata, and its bias depends on the corresponding differences of the average treatment effects in these paired strata. As a result, if the strata are paired in such a way that ensures “similar” strata are adjacent, then the bias of our estimator will naturally be small.

We then compare our variance estimator with a canonical variance estimator that has been proposed previously for the case in which $\ell = 1$ or $k - \ell = 1$. The estimator, proposed by imai2008variance, is based on the sample variance obtained by viewing the difference-in-means estimates for each stratum as independent observations. Existing results demonstrate that this estimator is also upward-biased for the variance of the difference-in-means estimator, but we demonstrate that the bias of our estimator is lower (and as a consequence leads to more precise inferences) whenever we can ensure that the “adjacent" strata used to construct our estimator are sufficiently similar. To further distinguish our estimator from the estimator proposed in imai2008variance, as well as a second estimator proposed by fogarty2018mitigating which modifies the Imai estimator by partialling out the average covariate values in each stratum, we introduce a framework motivated by a thought experiment in which the finite population can be modeled as being drawn once in an i.i.d.\ fashion from a well-behaved probability distribution. In this framework, we show that the limiting bias of our estimator is strictly smaller than that of imai2008variance unless the treatment effects are homogeneous, and strictly smaller than that of fogarty2018mitigating unless the conditional average treatment effect is linear in the covariates. As a result, confidence intervals based on our variance estimator will be strictly shorter than those based on the other variance estimators unless these exceptional restrictions hold.

The literature on stratified block randomization dates back to at least fisher1935design. For the case in which $\ell = 1$ or $k - \ell = 1$, the variance estimator in imai2008variance has been extended to different settings by imai2009essential, pashley2021insights, and zhu2024design-based. fogarty2018regression-assisted study regression adjustment for matched pair designs, and liu2020regression-adjusted studies regression adjustment for stratified designs primarily for the case where $\min\{\ell, k - \ell\} > 1$. ding2017paradox considers randomization inference for matched pairs, alongside other randomization schemes. We note that all of these papers employ a design-based framework as we do in this paper. As a consequence, their analyses differ from work which studies inference in finely stratified experiments from a superpopulation perspective bai2022optimality, bai2022inference, bai2024inference-b, bai2024covariate, bai2024inference, jiang2024bootstrap, cytrynbaum2024covariate, cytrynbaum2021optimal, bai2025efficiency; see also abadie2008estimation, who consider superpopulation estimation of the variance in matched pair designs, conditional on the covariates, in an alternative sampling framework where pairs are sampled instead of the units.

The remainder of the paper is organized as follows. In Section (ref), we describe our setup and notation. Section (ref) contains our main results. There we first highlight the importance of upward-biased variance estimators, and then review existing proposals before introducing our novel variance estimator. In Section (ref), we compare our variance estimators with existing ones in a particular limiting thought experiment. We examine the finite-sample behavior of confidence intervals based on all these estimators through a simulation study in Section (ref).

Setup and Notation

Consider an experiment consisting of $i \in \{1, \ldots, n\}$ units. For the $i$th unit, let $Y_{i} \in \mathbf{R}$ denote their observed outcome, $D_i \in \{0, 1 \}$ denote their received treatment, $X_i \in \mathbf{R}^{p}$ denote their observed, baseline covariates, and $Y_i(d)$ denote their potential outcome under treatment $d \in \{0,1\}$. As usual, the observed outcomes are related to the potential outcomes via the relationship

In what follows, it will be convenient to use the following shorthand notation: for a generic random vector $A_i$ indexed by $i$, let $A^{(n)} := (A_i: 1 \leq i \leq n)$. For later use, we also define $W_i := (Y_i(1),Y_i(0),X_i)'$.

In the design-based framework that we maintain throughout the paper, the potential outcomes and covariates of the units in the experiment are modeled as nonrandom quantities, with the only source of randomness arising from the treatment assignment mechanism. Our analysis concerns “finely stratified” designs in which the covariates are used to stratify units in the experiment into groups, i.e., strata, of fixed size $k$, and then, within each group, $\ell < k$ units are chosen uniformly at random to assign to treatment and the remainder to control; of particular interest to us will be the special case when $\min\{\ell, k - \ell\} = 1$. For ease of exposition, we assume throughout that $n = mk$, so that $m$ denotes the number of strata. Formally, we model the strata as a partition of $\{1, \dots, n\}$: $$\Lambda_n := \{ \lambda_j \subset \{1, \dots, n\}: 1 \leq j \leq m \}~,$$ with $|\lambda_j| = k$. It is worth emphasizing that we have suppressed in the notation the fact that each $\lambda_j$ (and therefore $\Lambda_n$ itself) can depend on $X^{(n)}$. Since we maintain that $k$ is fixed throughout the paper, when we write that $n \to \infty$, it should be understood that $m \to \infty$. Using this notation, the (joint) distribution of treatment assignment is characterized by the following assumption:

In other words, $\ell$ out of $k$ units in each stratum are treated uniformly at random, independently across strata.

Our parameter of interest is the finite population average treatment effect, given by $$\Delta_n := \bar Y_n(1) - \bar Y_n(0)~,$$ where, for $d \in \{0,1\}$, $$\bar Y_n(d) := \frac{1}{n} \sum_{1 \leq i \leq n} Y_i(d)~.$$ A natural estimator of $\Delta_n$ is given by the usual difference-in-means estimator, i.e., $$\hat \Delta_n := \frac{1}{n(1)} \sum_{1 \leq i \leq n} Y_i D_i - \frac{1}{n(0)} \sum_{1 \leq i \leq n} Y_i (1 - D_i)~,$$ where $n(1) = \sum_{1 \le i \le n} D_i = n \eta$, with $\eta := \ell/k$ the proportion of the $n$ units that are treated, and $n(0) = \sum_{1 \le i \le n}(1 - D_i) = n(1 - \eta)$. Note that because we have assumed that $n = mk$, both $n \eta = m \ell$ and $n ( 1 - \eta) = m (k - \ell)$ are integers.

Main Results

In this section, we begin, by way of motivation, with a result that highlights the importance of variance estimators that are at least asymptotically upward-biased, defined precisely in the statement of Theorem (ref) below. Under mild regularity conditions, we show, in particular, that valid inference for $\Delta_n$ based on a standard $t$-statistic requires a variance estimator that is asymptotically upward-biased for $\operatorname{Var}[\hat{\Delta}_n]$. Of course, a simple, sufficient condition for an estimator to be asymptotically upward-biased is that it is upward-biased. In Section (ref), we review some prior proposals for constructing variance estimators in finely stratified experiments that are upward-biased under only the assumption that $D^{(n)}$ satisfies Assumption (ref). We further discuss in Remark (ref) some improvements upon these estimators that typically result in estimators that are only asymptotically upward-biased under assumptions stronger than Assumption (ref) and therefore result in valid inference less generally than their upward-biased counterparts. Finally, in Section (ref), we present our main results. We introduce a novel estimator of $\operatorname{Var}[\hat{\Delta}_n]$ for finely stratified experiments and show that, like the estimators that we review in Section (ref), it is upward-biased for $\operatorname{Var}[\hat{\Delta}_n]$ under only the assumption that $D^{(n)}$ satisfies Assumption (ref). In contrast to these other estimators, however, we further argue that the magnitude of the bias of our estimator depends explicitly on the quality of the groupings of the experimental units into strata, in a sense made precise by Theorem (ref) below.

Motivating Result

The following theorem formalizes the sense in which an asymptotically upward-biased estimator for $\operatorname{Var}[\hat{\Delta}_n]$ is required for valid inference on $\Delta_n$.

Note that because $\tilde{V}_n$ is defined as an estimator of $\operatorname{Var}[\hat \Delta_n]$ and not the asymptotic variance of $\sqrt{n}(\hat \Delta_n - \Delta_n)$, the conditions in the theorem are all stated with a scaling by $n$. When (ref) holds, we say that $\tilde{V}_n$ is asymptotically upward-biased for $\operatorname{Var}[\hat{\Delta}_n]$. A simple sufficient condition for (ref) is, of course, that

When (ref) holds, we say that $\tilde{V}_n$ is upward-biased for $\operatorname{Var}[\hat{\Delta}_n]$. Theorem (ref) thus demonstrates that, under mild regularity conditions, confidence intervals for $\Delta_n$ of the form $[\hat \Delta_n \pm \sqrt{\tilde{V}_n}\cdot z_{1 - \alpha/2}]$ are valid in the sense of having the right limiting coverage probability, if and only if $\tilde{V}_n$ is asymptotically upward-biased for $\operatorname{Var}[\hat{\Delta}_n]$. These assumptions include weak restrictions on the finite population (i.e., (ref)), a non-degeneracy condition (i.e., (ref)), and a requirement that $n\cdot\tilde{V}_n$ is consistent for its expectation (i.e., (ref)). This last condition will typically hold under some further restrictions on finite population moments; see, e.g., Theorem (ref) below. Finally, we note that it is immediate from this result that the length of the interval $[\hat \Delta_n \pm \sqrt{\tilde{V}_n}\cdot z_{1 - \alpha/2}]$ will necessarily be (asymptotically, weakly) shorter whenever it is constructed using a variance estimator $\tilde{V}_n$ with correspondingly smaller asymptotic upwards bias. As a consequence, variance estimators with smaller bias lead to (asymptotically) more precise inferences.

The estimators discussed in Sections (ref) and (ref) below can be shown to satisfy (ref), and therefore (ref), whenever $D^{(n)}$ satisfies Assumption (ref). Such estimators therefore lead to confidence intervals for $\Delta_n$ that are valid whenever the hypotheses of Theorem (ref) are satisfied. It is worth emphasizing, however, that establishing (ref) for other estimators may in general require additional assumptions beyond the hypotheses of Theorem (ref). In Remarks (ref) and (ref) we discuss some examples of such estimators, and we illustrate using our simulations in Section (ref) that their resulting confidence intervals may lead to under-coverage when these additional conditions are not satisfied.

Review of Some Upwards Biased Estimators of $\operatorname{Var}[\hat{\Delta}_n]$

In this section, we review some upward-biased estimators of $\operatorname{Var}[\hat{\Delta}_n]$. To this end, recall that it can be shown using standard arguments imbens2015causal that

where \[S^2_j(d) := \frac{1}{k-1}\sum_{i \in \lambda_j}\left(Y_i(d)- \bar Y_{j, n}(d)\right)^2~, \hspace{5mm} S_{j, \Delta}^2 := \frac{1}{k-1}\sum_{i \in \lambda_j}\left(Y_i(1) - Y_i(0) - \Delta_{j,n}\right)^2~,\] with \[\bar Y_{j, n}(d) := \frac{1}{k}\sum_{i \in \lambda_j}Y_i(d)~, \hspace{5mm} \Delta_{j,n} := \frac{1}{k}\sum_{i \in \lambda_j}\left(Y_i(1) - Y_i(0)\right)~.\] In settings where $\min\{\ell, k - \ell\} > 1$, the construction of an upward-biased variance estimator is straightforward: an unbiased estimator of $S^2_j(d)$ for $d \in \{0, 1\}$ is given by the sample variance of the outcomes for units within stratum $j$ assigned to treatment $d$, which we denote by $\hat{S}^2_j(d)$ imbens2015causal, pashley2021insights. A simple estimator of $\operatorname{Var}[\hat{\Delta}_n]$ is thus given by

and its bias is $\frac{1}{nm}\sum_{1 \le j \le m}S_{j, \Delta}^2 \ge 0$. Note that this estimation strategy exploits the lower bound $S_{j, \Delta}^2 \ge 0$; we return to this observation in Remark (ref) below.

This estimation strategy fails, however, when there is exactly one treated or control unit in a stratum, i.e., when $\min\{\ell, k - \ell\} = 1$, since in this case the corresponding sample variance is identically zero. In these settings, a canonical estimator considered in the literature instead computes the variance of the stratum-level average treatment effect estimates imai2008variance: \[\hat{V}_n^{\rm IM} := \frac{1}{m (m - 1)}\sum_{1 \leq j \leq m}\left(\hat{\Delta}_{j,n} - \hat{\Delta}_n\right)^2~,\] where \[\hat{\Delta}_{j,n} := \frac{1}{\ell}\sum_{i \in \lambda_j}Y_iD_i - \frac{1}{k - \ell}\sum_{i \in \lambda_j}Y_i(1-D_i) \] is the difference in means in the $j$th stratum. This variance estimator serves as the basic scaffolding for many of the recent estimators proposed in the literature on inference in stratified randomized experiments: it is a special case of the small-block variance estimator proposed in pashley2021insights when all strata are the same size (see also zhu2024design-based), a special case of the pair-cluster variance estimator of de2024level in the setting of an individual-level randomized experiment, and a special case of the regression-based estimator due to fogarty2018mitigating, which we introduce in Section (ref). Note that the bias of $\hat{V}^{\rm IM}_n$ is given by imbens2015causal,fogarty2018mitigating

so that $\hat{V}^{\rm IM}_n$ is an upward-biased estimator of $\operatorname{Var}[\hat{\Delta}_n]$. Moreover, under mild assumptions, it is straightforward to establish using Chebyshev's inequality that $\hat V_n^{\rm IM}$ is consistent for its expectation in the sense of (ref), and thus can be used to construct valid confidence intervals in the sense of (ref).

A Novel Paired-Strata Estimator of $\operatorname{Var}[\hat{\Delta}_n]$

Building on our earlier work on super-population approaches to inference for finely stratified designs bai2022inference,bai2024inference-b, bai2024covariate,bai2024inference,bai2025efficiency, we now propose a novel estimator of $\operatorname{Var}[\hat{\Delta}_n]$, primarily for settings where $\min\{\ell, k - \ell\} = 1$. Our estimator is constructed by pairing together strata; for simplicity, we pair adjacent strata together, so the pairs are given by $\{(\lambda_{2j - 1}, \lambda_{2j}): 1 \le j \le \lfloor m/2 \rfloor\}$, but we emphasize that this pairing should be understood as being obtained after permuting the indices of the strata $\{\lambda_j: 1 \leq j \leq m\}$. In practice, it will be desirable to do this in a way that ensures “similar” strata are adjacent. With this is mind, our estimator is defined as follows:

where \[ \hat \tau_n^2 = \frac{1}{m} \sum_{1 \leq j \leq m} \hat \Delta_{j, n}^2\] and \[ \hat \kappa_n = \frac{2}{m} \sum_{1 \leq j \leq \lfloor \frac{m}{2} \rfloor} \hat \Delta_{2j - 1, n} \hat \Delta_{2j, n}~. \] Note that it follows from straightforward algebraic manipulation that $\hat{V}_n \ge 0$. The construction of $\hat V_n$ can be motivated using a specific limiting thought experiment that we introduce in Section (ref) below. The following theorem shows that $\hat V_n$ is an upward-biased estimator of $\operatorname{Var}[\hat{\Delta}_n]$ when $D^{(n)}$ satisfies Assumption (ref) and, under mild additional restrictions, is consistent for its expectation.

Theorem (ref) justifies the use of $\hat V_n$ for inference about $\Delta_n$ and further demonstrates that the magnitude of its bias depends on the differences of the stratum-level average treatment effects across pairs of strata, and hence the bias decreases whenever the pairs of strata are formed in a way that increases their homogeneity in this sense. The bias will therefore be small not only in the extreme instance in which the stratum-level average treatment effects are similar across all strata, but also when stratum-level average treatment effects vary across strata and are only similar for adjacent strata. Corollary (ref) documents the relative biases of $\hat{V}_n$ versus $\hat{V}^{\rm IM}_n$.

We note that, a sufficient condition for (ref) to hold is that the left-hand side is non-negative: we expect this to be the case whenever “similar" strata are paired together. We argue that such a condition holds quite generally for large $m$ in Section (ref) through additional restrictions on $W^{(n)}$ and $\Lambda_n$ that are motivated by a particular limiting thought experiment.

Comparison with Other Variance Estimators

In this section, we provide a framework that permits us to further discriminate among the variance estimators of $\operatorname{Var}[\hat{\Delta}_n]$ discussed in Section (ref). This framework is formalized by restrictions on the sequence of finite populations and experimental designs $(W^{(n)},\Lambda_n)$ specified in the assumptions below. These restrictions are motivated by a limiting thought experiment in which $W^{(n)}$ is itself a realization of an i.i.d.\ sample from a fixed probability distribution, and the strata $\Lambda_n$ are formed in such a way that units with similar observable characteristics are grouped together, but, as explained in Remark (ref) below, the restrictions can also be shown to hold in other contexts as well.

For the analysis of the estimator proposed by fogarty2018mitigating, introduced below, we need to further impose the following assumption:

Assumption (ref)(a) states that the finite population “moments” of the outcomes should converge to well-defined limit quantities; the condition holds almost surely by the strong law of large numbers if $W^{(n)}$ is modeled as an i.i.d.\ sample from the probability distribution $Q$. Assumption (ref)(b) states that the average products of the outcomes in stratum-pairs should converge to well-defined limit quantities in such a way that their covariate values are also being matched in the limit; this condition can also be justified if $W^{(n)}$ is modeled as an i.i.d.\ sample from $Q$ and the stratification $\Lambda_n$ satisfies the property that

following similar arguments to those used in Lemma C.2 in bai2024inference. In words, (ref) requires that the average of the squared distances of the covariate values in a paired stratum converges to zero asymptotically. This property is guaranteed by specific matching algorithms under appropriate assumptions; see, for instance, bai2022inference and cytrynbaum2021optimal for details. Assumption (ref) introduces similar conditions to Assumption (ref) but also requires that these conditions hold for the covariate vectors $X^{(n)}$. Note in particular that $E_Q[\tilde X \tilde Y(d)] = E_Q[\tilde X E_Q[\tilde Y(d) | \tilde X]]$ by the law of iterated expectations.

Using Assumption (ref), we have the following expression for the limit of $\operatorname{Var}[\hat{\Delta}_n]$ after normalizing appropriately:

We can now motivate the construction of $\hat V_n$ in Section (ref) through Theorem (ref). For the reasons outlined in Remark (ref), we focus on estimating

which is an upper bound for $V$. The construction of $\hat{V}_n$ can be motivated following arguments similar to those used in our earlier work on the analysis of finely stratified experiments from a super-population perspective bai2022inference,bai2024inference-b,bai2024inference,bai2025efficiency. First, by expanding the square, it can be shown under Assumption (ref) that

It thus remains to construct a consistent estimator of the final term and subtract it from $\hat{\tau}^2_n$. For this purpose, it can also be shown under Assumption (ref) that

To understand why it is natural to expect $\hat \kappa_n$ to satisfy (ref) under Assumption (ref)(b), note that $\hat \kappa_n$ is the average of products of the differences in means in a pair of strata. Each pair of differences in means are independent conditional on the stratification. Furthermore, for a sufficiently good stratification, each of them has conditional expectation close to $E_Q[\tilde Y(1) - \tilde Y(0) | \tilde X]$ for the same value of $\tilde X$. For this reason, we expect $\hat \kappa_n$ to converge in the desired way.

In the remainder of the section we compare the asymptotic variance $V$ to the (appropriately scaled) limit of $\hat{V}_n$ and to the limits of the variance estimators proposed in imai2008variance (i.e., $\hat{V}^{\rm IM}_n$), as well as the variance estimator proposed by fogarty2018mitigating, which we now describe. Let $R$ be an $m \times L$ matrix for $L < m$ with rank $L$ and $H_R = R (R' R)^{-1} R'$ denote the projection matrix onto its column space. For each $R$, consider the estimator given by

where $\hat{\delta}_n=(\hat \Delta_{1, n},\dots,\hat \Delta_{m, n})'$. Note that $\hat{V}_n^{\rm F}(\iota_m) = \hat{V}_n^{\rm IM}$, where $\iota_m$ is an $m \times 1$ vector of ones, so that $\hat{V}_n^{\rm F}(R)$ is indeed a generalization of $\hat{V}^{\rm IM}_n$. fogarty2018mitigating focuses on the estimator given by $\hat V_n^{\rm F} = \hat V_n^{\rm F}(Q_2)$, where \[ Q_2 =

\] for $\bar X_{j, n} = \frac{1}{k} \sum_{i \in \lambda_j} X_i$ and $\mu_{X, n} = \frac{1}{n} \sum_{1 \leq i \leq n} X_i$. In particular, fogarty2018mitigating argues that $\hat{V}_n^{\rm F}$ can be less conservative than $\hat{V}_n^{\rm IM}$ whenever the covariates are predictive of the treatment effects in an appropriate sense. The bias for $\hat V_n^{\rm F}$ is fogarty2018mitigating \[ E[\hat V_n^{\rm F}] - \operatorname{Var}[\hat \Delta_n] = \frac{1}{m^2} \delta_n' (\operatorname{diag}(I-H_{Q_2}))^{-1/2}(I-H_{Q_2})(\operatorname{diag}(I-H_{Q_2}))^{-1/2} \delta_n \geq 0~, \] where $\delta_n = (\Delta_{1, n}, \dots, \Delta_{m, n})'$, so $\hat V_n^{\rm F}$ is also an upward-biased estimator for $\operatorname{Var}[\hat \Delta_n]$.

As discussed in Section (ref), it is possible to establish under weak assumptions on the sequence of population “moments” that each of these variances is consistent for their expectations in the sense of (ref). Accordingly, in what follows it suffices to study the limit of the expectations of these estimators under Assumption (ref), and compare them to $V$.

Theorem (ref) demonstrates that $\hat{V}_n$ converges to our desired upper bound $V^{\rm obs}$, but that in general $\hat{V}^{\rm IM}_n$ and $\hat{V}^{\rm F}_n$ do not. Specifically, $\hat{V}^{\rm IM}_n$ only attains $V^{\rm obs}$ when treatment effects are sufficiently homogeneous, in the sense that $E_Q[\tilde{Y}(1) - \tilde{Y}(0)|\tilde{X}]$ is constant; $\hat{V}^{\rm F}_n$ only attains $V^{\rm obs}$ when the conditional average treatment effect $E_Q[\tilde{Y}(1) - \tilde{Y}(0)|\tilde X]$ is linear and therefore equal to the best linear predictor of $\tilde{Y}(1)-\tilde{Y}(0)$ given $1$ and $\tilde{X}$. As a consequence, Theorem (ref) demonstrates that, although confidence intervals constructed based on these variance estimators all cover $\Delta_n$ with probability at least $1 - \alpha$ in the limit, the confidence interval based on $\hat{V}_n$ is the least conservative whenever Assumptions (ref) (and (ref)) are satisfied. We further illustrate this phenomenon via simulation in Section (ref).

Simulations

In this section we examine the finite-$n$ behavior of several confidence intervals of the form described in (ref), constructed using different variance estimators. The variance estimators we consider are our proposed estimator $\hat{V}_n$, the imai2008variance estimator $\hat{V}^{\rm IM}_n$, the fogarty2018mitigating estimator $\hat{V}^{\rm F}_n$, and an alternative variance estimator $\hat{V}^{\rm alt}_n$ which can also be motivated using Theorem (ref) in our limiting thought experiment. Importantly, however, $\hat{V}^{\rm alt}_n$ is not necessarily upward-biased under Assumption (ref) alone. As a consequence, we will demonstrate that the resulting confidence interval will not appropriately cover $\Delta_n$ whenever the limiting thought experiment described in Assumption (ref) fails to hold. To construct $\hat{V}^{\rm alt}_n$, consider the following decomposition:

Under Assumption (ref)(a), straightforward consistent estimators of $\operatorname{Var}_Q[\tilde{Y}(d)]$ and $E_Q[\tilde{Y}(d)]^2$ are given by $\hat{\sigma}^2_n(d)$ and $\hat{\mu}_n(d)^2$, where

Under Assumption (ref)(b), an estimator of $E_Q[E_Q[\tilde{Y}(d)|\tilde{X}]^2]$ is \[ \hat \varsigma_n(d) = \frac{2}{m} \sum\limits_{1 \leq j \leq \lfloor\frac{m}{2}\rfloor} {\sum_{i \in \lambda_{2j}, i' \in \lambda_{2j - 1}}} Y_iY_{i'}I\{D_i = D_{i'} = d\}~.\] An alternative estimator of $V^{\rm obs}$ under Assumption (ref) is thus given by

\[\hat{V}^{\rm alt}_n := \frac{1}{n}\left(\frac{\hat{\sigma}^2_n(1) + \hat{\mu}_n(1)^2 - \hat{\varsigma}_n(1)}{\eta} + \frac{\hat{\sigma}^2_n(0) + \hat{\mu}_n(0)^2 - \hat{\varsigma}_n(0)}{1 - \eta} \right)~.\] We remark that, in the case of matched pairs, $\hat{V}_n^{\rm alt}$ also coincides with the estimator that would be obtained if we applied the estimation strategy in (ref) to the paired strata directly; see Remark (ref) for further discussion.

We generate our population $((Y_i(1), Y_i(0), X_i): 1 \le i \le n)$ using i.i.d.\ draws from a probability distribution. The potential outcomes are generated according to the equation:

where $((X_i, \epsilon_{0i}, \epsilon_{1i}): 1\leq i \leq n)$ are i.i.d., $(X_i, \epsilon_{0i}, \epsilon_{1i})$ are independent with $X_i \sim U[0,1]$, $\epsilon_{di} \sim N(0, 1)$, the parameters $\mu_d$ are given by $\mu_1 = 0.25$, $\mu_0 = 0$, and $\mu_d(\cdot)$ are specified as

For a population of size $n = 1000$, the potential outcomes and covariates are generated as i.i.d.\ draws once and then fixed in repeated samples. When considering populations of size $n \in \{100,250,500,750\}$, we draw a subsample from our population of size $1000$ once and then fix this in repeated samples.

In each Monte Carlo iteration, we assign units to treatment using a matched pairs design (i.e., a finely stratified design with $\ell = 1$ and $k = 2$), under two alternative matching methods, which we call “good match” and “bad match.” For “good match,” we sort units in increasing values of $X_i$ and then pair adjacent units. For “bad match,” we sort units in increasing values of $X_i$ and then pair them such that the unit with the smallest value of $X$ is matched with the unit with the largest value of $X$, the second smallest with the second largest, and so on. In both cases, when re-arranging the strata to construct our paired strata $\{(\lambda_{2j-1}, \lambda_{2j}): 1 \le j \le \lfloor m/2 \rfloor\}$, we sort strata in increasing values of their covariate averages and then pair adjacent pairs. As a result, the “good match” design minimizes the average squared distances of the covariate values in paired strata and can be shown to satisfy the conditions of our limiting thought experiment in Assumptions (ref) and (ref). In contrast, the “bad match” design is constructed specifically to fail the conditions of the limiting thought experiment.

Tables (ref) and (ref) report the coverage and average length of $95\%$ confidence intervals computed across 5000 Monte Carlo replications. When matches are good, we see that all estimators have appropriate (albeit conservative) coverage. In Table (ref) with good matches, the average length produced by $\hat{V}^{\rm IM}$ is largest while the average lengths produced by the other three estimators are comparable; this is in line with the theoretical results presented in Theorem (ref) given that the conditional average treatment effect is linear under Model 1. In Table (ref) with good matches, the average lengths produced by $\hat{V}_n$ and $\hat{V}_n^{\rm alt}$ are shortest, again in line with the theoretical results presented in Theorem (ref).

When matches are bad, the conditions in Theorem (ref) no longer apply, and instead we must rely on the asymptotic validity provided by Theorem (ref) alone. With this in mind we see that in both models, the estimators $\hat{V}_n^{\rm IM}$, $\hat{V}_n^{\rm F}$ and $\hat{V}_n$ all provide appropriate coverage; this is in line with our theoretical results given that all three estimators are upward-biased under Assumption (ref). The average lengths of the confidence intervals based on these estimators are comparable. On the other hand, $\hat{V}_n^{\rm alt}$, which is not guaranteed to be upward-biased in general, undercovers even in large populations.

Recommendations for Empirical Practice

Based on our theoretical results as well as the simulation study above, we conclude with some recommendations for empirical practice. Overall, we argue that our novel estimator $\hat{V}_n$ is well suited for performing design-based inference on the average treatment effect, particularly in settings with $\min\{\ell, k - \ell\} = 1$. Since the magnitude of the bias of this estimator depends explicitly on the quality of the pairing of the strata, we recommend that practitioners form the strata pairs by applying an optimal non-bipartite matching algorithm greevy2004optimal to the stratum-level covariate means $\bar{X}_j = \frac{1}{k}\sum_{i \in \lambda_j}X_i$; this algorithm is available in the {\tt R} package {\tt nbpMatching}.

An interesting feature of our results relative to our findings in prior work on the analysis of randomized experiments from a super-population perspective bai2022inference,bai2024inference-b,bai2024inference,bai2025efficiency is that, by appealing directly to Theorem (ref), asymptotically valid inference on $\Delta_n$ is possible under minimal assumptions on the nature of the matching procedure. In contrast, inference on the average treatment effect from a super-population perspective seems to require an assumption like (ref). However, we note that inferences based on the design-based estimator $\hat{V}_n$ are not guaranteed to be valid if we view the sample as being drawn from a larger (finite or super-) population. In such cases, it can be shown that the super-population variance estimators proposed in bai2022inference bai2024inference, bai2025efficiency remain valid in the limiting thought experiment of Section (ref).