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.
96,355 characters · 15 sections · 78 citation commands
Inference in Cluster Randomized Trials with Matched Pairs
KEYWORDS: Experiment, matched pairs, cluster-level randomization, randomized controlled trial, treatment assignment
JEL classification codes: C12, C14
\thispagestyle{empty} \setcounter{page}{1}
This paper studies the problem of inference in cluster randomized experiments where treatment status is determined according to a “matched pairs” design. Here, by a cluster randomized experiment, we mean one in which treatment is assigned at the level of the cluster; by a “matched pairs” design, we mean that the sample of clusters is paired according to baseline, cluster-level covariates and, within each pair, one cluster is selected at random for treatment. Cluster matched pair designs feature prominently in all parts of the sciences: examples in economics include angrist2009effects, fryer2014injecting, banerjee2015miracle, crepon2015estimating, bruhn2016impact, glewwe2016better, fryer2018pupil and romero2020outsourcing.
Following recent work in bugni2024inference, we develop our results in a sampling framework where clusters are realized as a random sample from a population of clusters. Importantly, in this framework cluster sizes are modeled as random and “non-ignorable," meaning that “large" clusters and “small" clusters may be heterogeneous, and, in particular, the effects of the treatment may vary across clusters of differing sizes. The framework additionally allows for the possibility of two-stage sampling, in which a subset of units is sampled from the set of units within each sampled cluster.
We first study the large-sample behavior of a weighted difference-in-means estimator under two distinct sets of assumptions on the matching procedure. Specifically, we distinguish between settings where the matching procedure does or does not match on a function of cluster size. For both cases, we establish conditions under which our estimator is asymptotically normal and derive simple, closed-form expressions for the asymptotic variance. Using these results, we establish formally that employing cluster size as a matching variable in addition to baseline covariates delivers a weak (and often strict) improvement in asymptotic efficiency relative to matching on baseline covariates alone, and in fact achieves full efficiency in a broad class of experimental designs: see Remark (ref) for further discussion. We then propose a variance estimator which is consistent for either asymptotic variance depending on the nature of the matching procedure. Combining these results establishes the asymptotic exactness of tests based on our estimators.
We then consider the asymptotic properties of two commonly recommended inference procedures based on linear regressions of the individual-level outcomes on a constant and cluster-level treatment. The first inference procedure clusters at the level of treatment assignment. The second inference procedure clusters at the level of assignment pairs, as recently recommended in de2019level. We establish that both procedures are generally conservative in our framework.
Next, we study the behavior of a randomization test which permutes the treatment status for clusters within pairs. We establish the finite-sample validity of such a test for testing a certain null hypothesis related to the equality of potential outcome distributions under treatment and control, and then establish asymptotic validity for testing null hypotheses about the size-weighted average treatment effect. We emphasize, however, that the latter result relies heavily on our choice of test statistic, which is studentized using our novel variance estimator. In simulations, we find that this randomization test controls size more reliably than any of the other inference procedures we consider in the paper, while delivering comparable power.
Finally, we derive large-sample results for a covariate-adjusted version of our estimator, which is designed to improve precision by exploiting additional baseline covariates which were not used for treatment assignment. As discussed in bai2024covariate and cytrynbaum2023covariate, standard covariate adjustments based on a regression using treatment-covariate interactions negi2021revisiting are not guaranteed to improve efficiency when treatment assignment is not completely randomized. For this reason, we consider a modified version of the estimator developed in bai2024covariate for individual-level matched pair experiments. Our results show that our covariate-adjusted estimator is guaranteed to improve asymptotic efficiency relative to the unadjusted estimator.
The analysis of data from cluster randomized experiments and data from experiments with matched pairs has received considerable attention donner2000design, athey2017econometrics, hayes2017cluster, but most recent work has focused on only one of these two features at a time. Recent work on the analysis of cluster randomized experiments includes middleton2015unbiased, su2021model, schochet2021design, and wang2022model bugni2024inference. We note in particular that both middleton2015unbiased and su2021model discuss the benefits of using cluster size as a covariate in regression adjustment in the context of completely randomized experiments. Recent work on the analysis of matched pairs experiments includes jiang2020bootstrap, cytrynbaum2021designing, bai2024inference, and bai2022optimality bai2022inference. Two papers which focus specifically on the analysis of cluster randomized experiments with matched pairs are imai2009essential and de2019level. Both papers maintain a finite-population perspective, where the primary source of uncertainty is “design-based," stemming from the randomness in treatment assignment. In such a framework, both papers study the finite and large-sample behavior of difference-in-means type estimators and propose corresponding variance estimators which are shown to be conservative. In contrast, our paper maintains a “super-population" sampling framework and proposes a novel variance estimator which is shown to be asymptotically exact in our setting. In Appendix (ref), we repeat some of the simulation exercises we consider in the main text in a design-based framework. There we illustrate that our estimator may have benefits in the design-based framework as well.
The remainder of the paper is organized as follows. In Section (ref) we describe our setup and notation. Section (ref) presents our main results. Section (ref) studies the finite-sample behavior of our proposed tests via a simulation study. We conclude with recommendations for empirical practice in Section (ref).
In this section we introduce the notation and assumptions which are common to both matching procedures considered in Section (ref). We broadly follow the setup and notation developed in bugni2024inference. Let $Y_{i,g} \in \mathbf{R}$ denote the (observed) outcome of interest for the $i$th unit in the $g$th cluster, $D_g \in \{0, 1\}$ denote the treatment received by the $g$th cluster, $X_g \in \mathbf{R}^k$ the observed, baseline covariates for the $g$th cluster, and $N_g \in \mathbf{Z}_+$ the size of the $g$th cluster. In what follows we sometimes refer to the vector $(X_g, N_g)$ as $W_g$. Further denote by $Y_{i,g}(d)$ the potential outcome of the $i$th unit in cluster $g$, when all units in the $g$th cluster receive treatment $d \in \{0, 1\}$. As usual, the observed outcome and potential outcomes are related to treatment assignment by the relationship
In addition, define $\mathcal{M}_g$ to be the (possibly random) subset of $\{1, 2, \ldots, N_g\}$ corresponding to the observations within the $g$th cluster that are sampled by the researcher. We emphasize that a realization of $\mathcal{M}_g$ is a set whose cardinality we denote by $|\mathcal{M}_g|$, whereas a realization of $N_g$ is a positive integer. For example, in the event that all observations in a cluster are sampled, $\mathcal{M}_g = \{1, \ldots, N_g\}$ and $|\mathcal{M}_g| = N_g$. We assume throughout that our sample consists of $2G$ clusters and denote by $P_G$ the distribution of the observed data $$Z^{(G)} := (((Y_{i,g} : i \in \mathcal{M}_g), D_g, X_g, N_g) : 1 \leq g \leq 2G )~,$$ and by $ Q_G$ the distribution of $$(((Y_{i,g}(1),Y_{i,g}(0) : 1 \le i \le N_g), \mathcal{M}_g, X_g, N_g) : 1 \leq g \leq 2G )~.$$ Note that $P_G$ is determined jointly by (ref) together with the distribution of $D^{(G)} := (D_g : 1 \leq g \leq 2G)$ and $Q_G$, so we will state our assumptions below in terms of these two quantities.
We now describe some preliminary assumptions on $Q_G$ that we maintain throughout the paper. In order to do so, it is useful to introduce some further notation. To this end, for $d \in \{0,1\}$, define $$\bar Y_g(d) := \frac{1}{|\mathcal{M}_g|} \sum_{i \in \mathcal{M}_g} Y_{i,g}(d)~.$$ Further define $R_G(\mathcal{M}_g^{(G)}, X^{(G)}, N^{(G)})$ to be the distribution of $$((Y_{i,g}(1), Y_{i,g}(0) : 1 \le i \le N_g) : 1 \leq g \leq 2G) ~\big |~ \mathcal{M}_g^{(G)}, X^{(G)}, N^{(G)}~,$$ where $\mathcal{M}_g^{(G)} := (\mathcal{M}_g : 1 \leq g \leq 2G)$, $X^{(G)} := (X_g : 1 \leq g \leq 2G)$ and $N^{(G)} := (N_g : 1 \leq g \leq 2G)$. Note that $Q_G$ is completely determined by $R_G(\mathcal{M}_g^{(G)}, X^{(G)}, N^{(G)})$ and the distribution of $(\mathcal{M}_g^{(G)}, X^{(G)}, N^{(G)})$. The following assumption states our main requirements on $Q_G$ using this notation.
For completeness, we reproduce some of the observations from bugni2024inference regarding these assumptions. Assumptions (ref)(a)--(b) formalize the idea that our sample consists of an i.i.d\ sample of clusters whose cluster sizes are random and potentially related to the potential outcomes. As shown in bugni2024inference, an important implication of Assumptions (ref)(a)--(b) for our purposes is that
is an i.i.d.\ sequence of random vectors. Assumptions (ref)(c)--(d) impose some mild regularity on the (conditional) moments of the distribution of cluster sizes and potential outcomes, in order to permit the application of relevant laws of large numbers and central limit theorems. Note that Assumption (ref)(c) does not rule out the possibility of observing arbitrarily large clusters, but does place restrictions on the heterogeneity of cluster sizes. For instance, two consequences of Assumptions (ref)(a) and (c) are that \[\frac{\sum_{1\le g \le G}N_g^2}{\sum_{1\le g\le G}N_g} = O_P(1)~,\] and \[\frac{\max_{1\le g \le G}N_g^2}{\sum_{1 \le g \le G}N_g} \xrightarrow{P} 0~,\] which mirror heterogeneity restrictions imposed in the analysis of clustered data when cluster sizes are modeled as non-random hansen2019asymptotic. We use Assumption (ref)(c) extensively when establishing asymptotic normality in Theorems (ref) and (ref); recent work by sasaki2022non and chiang2023cluster, however, suggests that one may be able to sometimes obtain asymptotic normality even when $E[N_g^2] = \infty$, provided that certain delicate conditions about the tail behavior of $N_g$ are satisfied. When the tails of the distribution of $N_g$ are so heavy that asymptotic normality fails, it may be possible to extend the recent work on subsampling based inference in chiang2023cluster to our setting, but we leave this extension for future work.
Assumptions (ref)(e)--(f) impose high-level restrictions on the two-stage sampling procedure. Assumption (ref)(e) allows the subset of observations sampled by the experimenter to depend on $X_g$ and $N_g$, but rules out dependence on the potential outcomes within the cluster itself. Assumption (ref)(f) is a high-level assumption which guarantees that we can extrapolate from the observations that are sampled to the observations that are not sampled. It can be shown that Assumptions (ref)(e)--(f) are satisfied if $\mathcal{M}_g$ is drawn as a random sample without replacement from $\{1, 2, \ldots, N_g\}$ in an appropriate sense bugni2024inference.
Our object of interest is the size-weighted cluster-level average treatment effect, which may be expressed in our notation as \[\Delta(Q_G) =E\left[ \frac{N_g}{E[N_g]}\left(\frac{1}{N_g}\sum_{1 \le i \le N_g}(Y_{i,g}(1) - Y_{i,g}(0))\right)\right] = E\left[ \frac{1}{E[N_g]}\sum_{1 \le i \le N_g}(Y_{i,g}(1) - Y_{i,g}(0))\right]~.\] This parameter, which weights the cluster-level average treatment effects proportional to cluster size, can be thought of as the average treatment effect where individuals are the unit of interest. Note that Assumptions (ref)(a)--(b) imply that we may express $\Delta(Q_G)$ as a function of $R$ and the common distribution of $(\mathcal{M}_g, X_g, N_g)$. In particular, this implies that $\Delta(Q_G)$ does not depend on $G$. Accordingly, in what follows we simply denote $\Delta = \Delta(Q_G)$.
In Sections (ref)--(ref), we study the asymptotic behavior of the following size-weighted difference-in-means estimator:
where \[\hat{\mu}_G(d) := \frac{1}{N(d)}\sum_{1 \le g \le 2G}I\{D_g = d\}\frac{N_g}{|\mathcal{M}_g|}\sum_{i\in \mathcal{M}_g}Y_{i,g}~,\] with \[N(d) := \sum_{1 \le g \le 2G}N_gI\{D_g = d\}~.\] Note that this estimator may be obtained as the estimator of the coefficient of $D_g$ in a weighted least squares regression of $Y_{i,g}$ on a constant and $D_g$ with weights equal to $N_g/|\mathcal{M}_g|$. In the special case that all observations in each cluster are sampled, so that $\mathcal{M}_g = \{1, 2, \ldots, N_g\}$ for all $1 \le g \le G$ with probability one, this estimator collapses to the standard difference-in-means estimator. However, it is important to note that outside of this special case, the standard difference-in-means estimator is not consistent for the size-weighted average treatment effect $\Delta$, and is instead consistent for an “$|\mathcal{M}_g|$-weighted" treatment effect; see bugni2024inference for details. In Section (ref) we consider a covariate-adjusted modification of $\hat{\Delta}_G$ which is designed to incorporate additional baseline covariates which were not used for treatment assignment.
In this section, we consider the asymptotic behavior of $\hat{\Delta}_G$ for two distinct types of cluster-matched pair designs. Section (ref) studies a setting where cluster size is not used as a matching variable when forming pairs. Section (ref) considers the setting where we do allow for pairs to be matched based on cluster size in an appropriate sense made formal below.
In this section, we consider a setting where cluster size is not used as a matching variable. First, we describe our formal assumptions on the mechanism determining treatment assignment. The $G$ pairs of matched clusters may be represented by the sets \[\{\pi(2j - 1), \pi(2j)\} \text{ for } j = 1, \dots, G,\] where $\pi = \pi_G(X^{(G)})$ is a permutation of $2G$ elements, and the right-hand side of this equality emphasizes that, since the permutation represents the result of the matching procedure, it is in fact a function of the cluster-level covariates $X^{(G)}$. Given such a $\pi$, we assume that treatment status is assigned as follows:
Assumption (ref) states that, after pairs are formed according to the baseline covariates, which cluster is treated in a pair is determined by a coin flip independently of all other variables. We further require that the clusters in each pair be “close" in terms of their baseline covariates in the following sense:
bai2022inference provide results which facilitate the construction of pairs which satisfy Assumption (ref). For instance, if $\mathrm{dim}(X_g) = 1$, then by simply pairing clusters by ordering them from smallest to largest according to $X_g$ and then pairing adjacent clusters, it follows from Theorem 4.1 in bai2022inference that Assumption (ref) is satisfied if $E[X_g^2] < \infty$. When $\dim(X_g) > 1$ and a suitable matching procedure is used (for instance the {\tt nbpmatching} package in R), it follows from the discussion in Appendix (ref) that Assumption (ref) is satisfied when $E[\|X_g\|^d] < \infty$ for $d \geq \mathrm{dim}(X_g) + 1 $ .
Next, we state the additional assumptions on $Q_G$ we require beyond those stated in Assumption (ref):
Assumption (ref)(a) is a smoothness requirement analogous to Assumption 2.1(c) in bai2022inference that ensures that units within clusters which are “close" in terms of their baseline covariates are suitably comparable. If $X_g$ is discrete and clusters are matched perfectly in that the distance between pairs in Assumption (ref) is zero, Assumption (ref)(a) is not needed. Assumption (ref)(b) imposes an additional restriction on the distribution of cluster sizes beyond what is stated in Assumption (ref)(c). Under these assumptions, we obtain the following result:
The proof of Theorem (ref) proceeds by studying the joint distribution of the random numerators and denominators of $\hat{\mu}_G(d)$ for $d \in \{0, 1\}$ using techniques similar to those used in bai2022inference, carefully taking into consideration the potential dependence between cluster sizes and outcomes, and then applying the Delta method. Remarkably, the resulting asymptotic variance we obtain in Theorem (ref) corresponds exactly to the asymptotic variance of the difference-in-means estimator for matched pairs designs with individual-level assignment bai2022inference, but with transformed cluster-level potential outcomes given by $\tilde{Y}_g(d)$. Accordingly, our result collapses exactly to theirs when $P\{N_g = 1\} = 1$.
In this section, we repeat the exercise in Section (ref) in a setting where the assignment mechanism matches on baseline characteristics and (some function of) cluster size in an appropriate sense to be made formal below. Recall the definition $W_g = (X_g, N_g)$, and let $W^{(G)}:=(W_g: 1 \le g \le 2G)$. First, we describe how to modify our assumptions on the mechanism determining treatment assignment. The $G$ pairs of clusters are still represented by the sets \[\{\pi(2j - 1), \pi(2j)\} \text{ for } j = 1, \dots, G~,\] however, now we allow the permutation $\pi = \pi_G(W^{(G)})$ which determines the pairing to depend on cluster sizes as well as $X^{(G)}$. Given such a $\pi$, we now assume that treatment status is assigned as follows:
We also require some modifications on our regularity conditions for how pairs are formed and our smoothness requirements on the potential outcomes; we provide further discussion in Remark (ref) below:
Under our modified matching procedure and regularity conditions, we obtain the following analog to Theorem (ref):
Note that the asymptotic variance $\nu^2$ has exactly the same form as $\omega^2$ from Section (ref), with the only difference being that the final term of the expression conditions on both cluster characteristics $X_g$ and cluster size $N_g$.
In this section, we construct variance estimators for the asymptotic variances $\omega^2$ and $\nu^2$ obtained in Section (ref). In fact, we propose a single variance estimator that is consistent for both $\omega^2$ and $\nu^2$ depending on the nature of the matching procedure. As noted in the discussion following Theorem (ref), the expressions for $\omega^2$ and $\nu^2$ correspond exactly to the asymptotic variance obtained in bai2022inference with the individual-level outcome replaced by a cluster-level transformed outcome. We thus follow the variance construction from bai2022inference, but replace the individual outcomes with feasible versions of these transformed outcomes. To that end, consider the observed adjusted outcome defined as:
where
We then propose the following variance estimator:
where
Note that the construction of $\hat{v}^2_G$ can be motivated using the same intuition as the variance estimators studied in bai2022inference and bai2024inference: to consistently estimate quantities like (for instance) $E[E[\tilde{Y}_g(1)|X_g]^2]$ which appear in $\omega^2$, ideally we would like to average over the products of the average outcomes of two treated clusters with similar values of covariates. By construction, however, only one cluster in each pair is treated, and our solution is to instead average across “pairs of pairs" of clusters. As a consequence, we will additionally require that the matching algorithm satisfy the condition that “pairs of pairs" of clusters are sufficiently close in terms of their baseline covariates/cluster size, as formalized in the following two assumptions:
As noted in bai2022inference, given pairs which satisfy Assumptions (ref) or (ref), it is possible to reorder the pairs so that Assumptions (ref) or (ref) are satisfied. We then obtain the following two consistency results for the estimator $\hat{v}_G^2$:
By combining Theorems (ref)--(ref) with Theorem (ref), asymptotically exact tests and confidence intervals can be constructed using a $t$-statistic studentized by $\hat{v}_G$. Next, we derive the limits in probability of two commonly recommended variance estimators obtained from a (weighted) linear regression of the individual-level outcomes $Y_{i,g}$ on a constant and cluster-level treatment $D_g$. The first variance estimator we consider, which we denote by $\hat{\omega}^2_{\rm CR,G}$, is simply the cluster-robust variance estimator of the coefficient of $D_g$ as defined in equation (ref) in the appendix. Theorem (ref) derives the limit in probability of $\hat{\omega}^2_{\rm CR,G}$ under a matched pair design which matches on baseline covariates as defined in Section (ref), and shows that it is generally too large relative to $\omega^2$.
The next variance estimator we consider, which we denote by $\hat{\omega}^2_{\rm PCVE,G}$, is the variance estimator of the coefficient of $D_g$ obtained from clustering on the assignment pairs of clusters as defined in equation (ref) in the appendix. de2019level call this the pair-cluster variance estimator (PCVE)\footnote{We emphasize, however, that de2019level propose their variance estimator in a finite population “design-based" inferential framework, which is distinct from the superpopulation framework we consider here. In Appendix (ref) we repeat some of the simulation exercises we consider in Section (ref) in a design-based framework. There we illustrate that our estimator may have benefits in the design-based framework as well.}. Theorem (ref) derives the limit in probability of $\hat{\omega}^2_{\rm PCVE,G}$ in the special case where $N_g = n$ for $g = 1,\ldots, 2G$ for some fixed $n$ and $|\mathcal{M}_g| = N_g$, and shows that it is generally too large relative to $\omega^2$.
Although we do not derive the limit in probability of $\hat{\omega}^2_{\rm PCVE, G}$ in the general case, our simulation evidence in Section (ref) suggests that the limit of $\hat{\omega}^2_{\rm PCVE, G}$ remains conservative, and that the conditions under which it is consistent for $\omega^2$ are the same as those in equation (ref). From Theorems (ref) and (ref) we obtain that neither cluster-robust standard error is consistent for $\omega^2$ unless the baseline covariates are irrelevant for the potential outcomes in an appropriate sense. In particular, equation (ref) holds when the average treatment difference for the sampled units in a cluster are homogeneous, in the sense that $\bar{Y}_g(1) - \bar{Y}_g(0)$ is constant. We further note that the conditions under which $\hat{\omega}^2_{\rm CR,G}$ and $\hat{\omega}^2_{\rm PCVE,G}$ are consistent for $\omega^2$ are exactly analogous to the conditions under which bai2022inference derive (in the setting of an individual-level matched pairs experiment) that the two-sample $t$-test and matched pairs $t$-test are asymptotically exact, respectively.
In this section, we study the properties of a randomization test based on the idea of permuting the treatment assignments for clusters within pairs. In Section (ref) we present some finite-sample properties of our proposed test, and in Section (ref) we establish its large sample validity for testing the null hypothesis $H_0: \Delta = 0$.
First, we define the test. In words, the randomization test constructs its critical value from the empirical distribution of the test statistic obtained by permuting the treatment assignments within pairs. In practice, such a distribution can be approximated by randomly permuting the treatment status of clusters within the same pair: for each pair of clusters, the treatment status of the two clusters remains the same with probability one-half and is flipped otherwise. The test statistic is then calculated based on these permuted treatment assignments and the critical value is determined by the $1 - \alpha$ quantile of resulting distribution of all such permutation statistics. Formally, denote by $\mathbf{H}_G$ the group of all permutations on $2G$ elements and by $\mathbf{H}_G(\pi)$ the subgroup that only permutes elements within pairs defined by $\pi$: \[\mathbf{H}_G(\pi) = \{h \in \mathbf{H}_G: \{\pi(2j-1), \pi(2j)\} = \{h(\pi(2j-1)), h(\pi(2j))\} \text{ for } 1 \le j \le G\}~.\] Define the action of $h \in \mathbf{H}_G(\pi)$ on $Z^{(G)}$ as follows: \[hZ^{(G)} = \{((Y_{i,g}: i \in \mathcal{M}_g), D_{h(g)}, X_g, N_g): 1 \le g \le 2G\}~.\] The randomization test we consider is then given by \[\phi_G^{\rm rand}(Z^{(G)}) = I\{T_G(Z^{(G)}) > \hat{R}^{-1}_G(1 - \alpha)\}~,\] where \[\hat{R}_G(t) = \frac{1}{|\mathbf{H}_G(\pi)|}\sum_{h \in \mathbf{H}_G(\pi)}I\{T_G(hZ^{(G)}) \le t\}~,\] with \[T_G(Z^{(G)}) = \left|\frac{\sqrt{G}\hat{\Delta}_G}{\hat{v}_G}\right|~.\]
In this section we present some finite-sample properties of the proposed test. Consider testing the null hypothesis that the distribution of potential outcomes within a cluster are equal across treatment and control conditional on observable characteristics and cluster size:
Note (ref) is stronger than the statement that the average treatment effect $\Delta = 0$. As a consequence, we are able to establish the following result on the finite sample validity of our randomization test for testing (ref):
In this section, we establish the large-sample validity of the randomization test $\phi_G^{\rm rand}$ for testing the null hypothesis
Note (ref) is implied by (ref). In Remark (ref) we describe how to modify the test for testing non-zero null hypotheses.
Theorems (ref) and (ref) highlight that the randomization test $\phi^{\rm rand}_G(Z^{(G)})$ is asymptotically valid for testing (ref) while additionally retaining the finite-sample validity described in Section (ref) under the null hypothesis (ref). In Section (ref) we illustrate the benefit of this additional robustness on the small-sample behavior of $\phi^{\rm rand}_G(Z^{(G)})$ relative to tests constructed using Gaussian critical values. We note that, unlike for the null hypothesis considered in Section (ref), the choice of test statistic $T_G$ is crucial for establishing Theorem (ref). Similar observations have been made in related contexts in janssen1997studentized, chung2013exact, bugni2018inference and bai2022inference.
In this section, we consider a linearly covariate-adjusted modification of $\hat{\Delta}_G$ that is designed to improve precision by exploiting additional observed baseline covariates that were not used for treatment assignment. To that end, we consider a setting in which we observe two sets of baseline covariates, $X_g$ and $C_g$, where $X_g \in \mathbf R^k$ denotes the original set of baseline covariates used for treatment assignment, and $C_g \in \mathbf R^{\ell}$ denotes the covariates in addition to $X_g$ that were not used for treatment assignment. Note that $C_g$ could also include cluster-level aggregates of individual-level outcomes, including intracluster means and quantiles. Before proceeding, we note that for the remainder of Section (ref), Assumption (ref) should be understood to hold with $(X_g, C_g)$ in place of $X_g$.
Our primary focus will be on settings in which the cluster size $N_g$ is used in determining the pairs. We note that similar results continue to hold under suitable modifications of our assumptions when $N_g$ is not used in determining pairs by simply replacing $W_g$ with $X_g$ throughout. As in Section (ref), let $\pi = \pi_G(W^{(G)})$ denote the permutation that determines the pairs. We then assume that treatment status is assigned as follows:
We consider a linearly covariate-adjusted estimator of $\Delta$ based on a set of regressors generated by $W_g$ and $C_g$; define $\psi_g = \psi(W_g, C_g)$, where $\psi: (\mathbf R^k \times \mathbf Z_+) \times \mathbf R^\ell \to \mathbf R^p$. We impose the following assumptions on $\psi$:
Assumption (ref)(a) implies that none of the components of $\psi_g$ can be perfectly predicted only by $W_g$. Assumptions (ref)(b)--(c) form the counterpart to Assumption (ref)(d), and Assumption (ref)(d) is the counterpart to Assumption (ref).
As discussed in bai2024covariate and cytrynbaum2023covariate, standard covariate adjustments based on a regression using treatment-covariate interactions negi2021revisiting are not guaranteed to improve efficiency when treatment assignment is not completely randomized. For this reason, we consider a modified version of the adjusted estimator developed in bai2024covariate for individual-level matched pair experiments. Let $\hat \beta_G$ denote the OLS estimator of the slope coefficient in the linear regression of $\left ( \frac{1}{2G} \sum_{1 \leq g \leq 2G} N_g \right ) (\hat Y_{\pi(2g - 1)} - \hat Y_{\pi(2g)})(D_{\pi(2g - 1)} - D_{\pi(2g)})$ on a constant and $(\psi_{\pi(2g - 1)} - \psi_{\pi(2g)})(D_{\pi(2g - 1)} - D_{\pi(2g)})$. We then define our covariate-adjusted estimator as
where \[ \bar \psi_G = \frac{1}{2G} \sum_{1 \leq g \leq 2G} \psi_g~. \] Theorem (ref) derives the limiting distribution of $\hat \Delta_G^{\rm adj}$, and, importantly, it shows that the limiting variance of $\hat \Delta_G^{\rm adj}$ is no larger than that of $\hat \Delta_G$ in (ref) and is strictly smaller unless $\psi$ is “irrelevant” for $\tilde Y_g(1) + \tilde Y_g(0)$ after “controlling” for $W_g$, in the sense made precise below.
Note that the asymptotic variance $\varsigma^2$ has the same form as the variance $\nu^2$, but with new transformed outcomes $Y^*_g(d)$ which can be expressed as covariate-adjusted versions of the original transformed outcomes $\tilde{Y}_g(d)$. Exploiting this observation is what allows us to establish that $\varsigma^2 = \nu^2 - \kappa^2$. As a consequence, we find that the asymptotic variance of $\hat{\Delta}^{\rm adj}_G$ is lower than that of $\hat{\Delta}_G$ whenever the adjustment is appropriately “relevant," in the sense that $\kappa^2 \ne 0$.
For variance estimation, define
We then propose the following variance estimator:
where
The following theorem establishes the consistency of the variance estimator:
In this section, we examine the finite-sample behavior of the estimation and inference procedures considered in Sections (ref)-(ref). We further compare these procedures to tests and confidence intervals constructed using the standard cluster-robust variance estimator (CR) and the pair cluster variance estimator (PCVE) proposed in de2019level. For $d \in \{0, 1\}$, $1 \le g \le 2G$, the potential outcomes are generated according to the equation \[Y_{i,g}(d) = \mu_d(X_g, C_g) + 2\epsilon_{d,i,g}~.\] Where, in each specification, $(X_g, C_g)$, $g = 1, \ldots, 2G$ are i.i.d.\ with $X_g, C_g\sim\ Beta(2,4)$, and $(\epsilon_{0,i,g}, \epsilon_{1,i,g})$, $g = 1, \ldots, 2G$, $i = 1, \ldots, N_g$ are i.i.d.\ with $\epsilon_{0,i,g}, \epsilon_{1,i,g} \sim N(0,1)$ independently. Note that $C_g$ are additional cluster level covariates which are used to determine the cluster size $N_g$, but are not used directly for matching. Throughout Section (ref) we assume that we observe the entire cluster, that is, we assume $\mathcal{M}_g = \{1, 2, \ldots, N_g\}$; in Appendix (ref) we repeat the simulation exercise in Section (ref) for other choices of $\mathcal{M}_g$. We consider the following two specifications for $\mu_d$:
Note that Model 1 satisfies the homogeneity condition in (ref) whereas Model 2 does not. In both cases, $N_g$, $g=1,\dots,2G$ are i.i.d.\ with $N_g \sim\ Binomial(R, C_g)+(500-R)$, where $R$ determines the difference in maximum and minimum cluster sizes. In particular $R$ satisfies the property that $N_g \in [N_{min}, N_{max}]$ with $N_{max}-N_{min}=R$ and we consider $R \in \{49,149,249,349,449\}$ with $N_{max}=500$ fixed. For each model and distribution of cluster sizes, we consider two alternative pair-matching procedures. First, we consider a design which matches clusters using $X_g$ only. To construct these pairs, we sort the clusters according to $X_g$ and pair adjacent clusters. Next, we consider a design which matches clusters using both $X_g$ and $N_g$. To construct these pairs, we match the clusters according to their Mahalanobis distance using the non-bipartite matching algorithm from the R package {\tt nbpMatching}.
Tables (ref)--(ref) report the coverage and average length of $95\%$ confidence intervals constructed using our variance estimator as well as the CR and PCVE estimators. For Model 1 in Table (ref), we find that, in accordance with Theorems (ref)--(ref), the CR variance estimator is extremely conservative, whereas our proposed variance estimator (denoted $\hat{v}_G^2$) and the PCVE variance estimator have exact coverage asymptotically. This feature translates to significantly smaller confidence intervals: on average the confidence intervals constructed using $\hat{v}_G^2$ or PCVE are almost half the length of those constructed using CR when $G \ge 50$. However, the confidence intervals constructed using $\hat{v}_G^2$ or PCVE undercover when $G < 50$. We find similar results when matching on both $X_g$ and $N_g$ in Table (ref). Comparing across Tables (ref) and (ref) we find that, in line with the discussions following Theorems (ref) and (ref), matching on $N_g$ in addition to $X_g$ results in a large reduction in the average length of confidence intervals constructed using $\hat{v}_G^2$ (or PCVE), but no change in the average length of confidence intervals constructed using CR.
Moving to Model 2 in Tables (ref) and (ref), here we find that confidence intervals constructed using CR continue to be conservative, but now the confidence intervals constructed using PCVE are also conservative, and numerically very similar to those constructed using CR. In contrast, the confidence intervals constructed using $\hat{v}_G^2$ remain exact asymptotically. Once again this translates to smaller confidence intervals for $\hat{v}_G^2$: on average the confidence intervals constructed using $\hat{v}_G^2$ are approximately $25\%$ smaller than those constructed using CR or PCVE when $G \ge 50$. However, once again we find that the confidence intervals constructed using $\hat{v}^2_G$ can undercover when $G < 50$, with the size of the distortion growing as a function of the cluster size heterogeneity.
Next, to further address the small-sample coverage distortions observed in Tables (ref)-(ref), we study the size and power of $0.05$-level hypothesis tests conducted using our proposed randomization test, as well as standard $t$-tests constructed using the CR and PCVE estimators, in Tables (ref)--(ref) below.\footnote{Here we move to studying the properties of hypothesis tests instead of confidence intervals to avoid having to perform test-inversion for our randomization test, but we expect that similar results would continue to hold for confidence intervals as well.} In Table (ref) we find that tests based on the CR variance estimator are extremely conservative, and this translates to having essentially no power against our chosen alternative. Tests based on the PCVE estimator produce non-trivial power, but also size-distortions in small samples. In contrast, since Model 1 satisfies the null hypothesis considered in (ref), our randomization test is valid in finite samples by construction, and displays comparable power to the PCVE-based test even when the latter does not control size. When moving to Model 2 in Table (ref) we are only guaranteed that the randomization test is asymptotically valid, but we find that the test is still able to control size in small samples as long as cluster-size heterogeneity is not too large. Importantly, in such cases, both the CR and PCVE-based tests also fail to control size. Finally, the randomization test displays favorable power relative to both the CR and PCVE-based tests throughout Table (ref) except for some cases when $G = 12$.
In this section, we examine the finite-sample behavior of the covariate-adjusted estimator considered in Section (ref). We consider the following modification of Model 2: let $C_g = (C_{1g}, C_{2g})$,
with $X_g \sim U[0,1]$ generated independently of all other variables, and modify the distribution of $N_g$ so that $N_g \sim Binomial(R, 1 - C_{2g}) + (500 - R)$.
Tables (ref) and (ref) report the coverage and average length of $95\%$ confidence intervals constructed using our variance estimators when matching using $X_g$ and both $X_g$ and $N_g$, respectively, for $\hat{\Delta}_G$ versus $\hat{\Delta}^{\rm adj}_G$ with $\psi_g = C_g$. In accordance with Theorem (ref), we find that for moderate to large samples ($G \ge 50$), covariate adjustment leads to smaller average CI lengths.
Based on our theoretical results as well as the simulation study above, we conclude with some recommendations for practitioners when conducting inference for cluster matched pair designs. The methods in this paper are primarily tailored for inference in a super-population framework; as explained in bai2024primer, such a sampling framework may be viewed as an approximation to a regime where a small fraction of the total population of clusters is sampled. Simulation evidence in Appendix (ref), however, suggests that our methods compare favorably against existing methods even in finite-population settings. Formal results in a finite population framework can be established by following the general strategy presented in Appendix A.1 in bai2024primer.
Our recommendations depend on whether the number of clusters is moderately large (e.g., at least 50 pairs) or small (e.g., less than 50 pairs). If the number of clusters is moderately large, then our recommendation is that practitioners should employ either the covariate-adjusted tests based on the covariate-adjusted estimator $\hat{\Delta}^{\rm adj}_G$ defined in Section (ref) paired with its corresponding variance estimator $\mathring\varsigma_G^2$ and a normal critical value or the unadjusted tests based on the unadjusted estimator $\hat{\Delta}_G$ introduced in Section (ref) paired with its corresponding variance estimator $\hat{v}^2_G$ and a normal critical value.
If, on the other hand, the number of clusters is small, then we recommend instead that practitioners use the randomization test based on the un-adjusted estimator $\hat{\Delta}_G$ paired with its corresponding variance estimator $\hat{v}^2_G$ outlined in Section (ref). In our simulations, this test controlled size more reliably than any of the other inference procedures we considered in the paper, while delivering comparable power. Note that by modifying the test as in Remark (ref), the test could also be inverted to construct confidence intervals if desired.
In general, all of our results crucially hinge on the assumption that clusters in a pair are sufficiently “close” (Assumptions (ref) and (ref)), and such a condition becomes difficult to satisfy as the dimension of $X_g$ increases. For this reason, we recommend that practitioners construct their pairs using a small subset of the baseline covariates that they believe have the highest explanatory power (including possibly cluster size itself). The experimental data can then be analyzed by using either the un-adjusted or adjusted estimators we propose in this paper.