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.
102,889 characters · 13 sections · 125 citation commands
Inference for Two-stage Experiments under Covariate-Adaptive Randomization
KEYWORDS: Randomized controlled trials, two-stage randomization, matched pairs, stratified block randomization, causal inference under interference
JEL classification codes: C13, C21 \hypersetup{pageanchor=false} \thispagestyle{empty} \hypersetup{pageanchor=true} \setcounter{page}{1}
This paper considers the problem of inference in two-stage randomized experiments under covariate-adaptive randomization. Here, a two-stage randomized experiment refers to a design where clusters (e.g., households, schools, or graph partitions) are initially randomly assigned to either a control or treatment group. Subsequently, random assignment of units within each treated cluster to either treatment or control is carried out based on a pre-determined treated fraction. Covariate-adaptive randomization refers to randomization schemes that first stratify according to baseline covariates and then assign treatment status so as to achieve “balance” within each stratum. Two-stage randomized experiments are widely used in social science (see for example Duflo2003,Haushofer2016,McKenzie2021), and discussed by statisticians (see for example hudgens2008), as a general approach to causal inference with interference; that is, when one individual's treatment status affects outcomes of other individuals. Moreover, practitioners often use covariate information to design more efficient two-stage experiments Duflo2003,Ichino2012,Beuermann2015,Muralidharan2015,Hidrobo2016,Rogers2018,Kinnan2020,Banerjee2021,Malani2021. However, to the best of my knowledge, there has not yet been any formal analysis on covariate-adaptive randomization in two-stage randomized experiments. Accordingly, this paper establishes general results about estimation and inference for two-stage designs under covariate-adaptive randomization. Subsequently, I propose and examine the optimality of two-stage designs with “matched tuples”, i.e. a generalized matched-pair design (see Bai-optimal and matched-tuple).
This paper examines covariate-adaptive randomization for two-stage experiments within a comprehensive framework that encompasses matched tuples designs, stratified block randomization, and complete randomization as special cases. The framework relies on finely stratified randomization (see cytrybaum2021 and bai2024efficiency), which involves grouping clusters into homogeneous strata of size $k$ and then assigning treatment entirely at random within each stratum.\footnote{The terms “cluster” and “stratum” are both used in the literature to describe groupings of units, which can lead to confusion. Here, a cluster is defined as a pre-determined group of units (e.g., households, schools, or graph partitions), and a stratum as a group of clusters that share similar baseline cluster-level covariates.} Within this framework, I propose a set of difference-in-“average of averages” estimators and analyze the statistical inference for four parameters of interest: equally-weighted and size-weighted primary effects, and equally-weighted and size-weighted spillover effects, under the assumption of homogeneous partial interference, where interference is confined within clusters. I establish conditions under which these four estimators are asymptotically normal and construct consistent estimators of their corresponding asymptotic variances. These results collectively validate the asymptotic validity of tests based on these estimators.
This paper then considers the asymptotic properties of a commonly recommended inference procedure based on a linear regression with cluster-robust standard errors. My findings suggest that the corresponding $t$-test is generally valid but conservative. I also demonstrate that in the first stage of cluster-level assignment, covariate information about clusters is important for both designing efficient experiments and consistently estimating variances under covariate-adaptive randomization. However, in the second stage of unit-level assignment, while individual-level covariate information is useful for improving efficiency, it is not required for the proposed inference method. Specifically, I show that consistent variance estimators can be constructed using only the cluster-level covariates from the first stage design, regardless of the use of individual-level covariates in the second stage.
Next, I apply the results to study optimal use of covariate information in two-stage designs. Here, by “optimal”, I mean designs that achieve the minimum asymptotic variances within the class of designs considered in the paper. For all estimands of interest, the designs in the first and second stage affect the efficiency independently. Thus, I am able to identify optimal designs in the first and second stage separately and use them together as the optimal two-stage design. My result shows that, at each stage, the asymptotically optimal design is a matched tuples design where clusters or units are matched based on an index function (similar to Bai-optimal) that is specific to the given estimator. In a simulation study, the results demonstrate that properly designed two-stage experiments utilizing the optimality results outperform other designs. However, the efficiency gain achieved through proper second-stage randomization is significantly lower compared to the first stage under my simulation specifications.
In the empirical literature, it is common to match or stratify on a small set of covariates expected to be most predictive of outcomes and to adjust for other pre-treatment covariates ex-post. Building on mp-cluster and bai2023, I propose a covariate-adjusted version of my estimator and discuss the conditions under which this estimator enhances asymptotic efficiency compared to the unadjusted version.
Finally, this paper evaluates the proposed inference method against various regression-based methods commonly used in empirical literature in a simulation study and empirical application. The simulation study confirms the asymptotic exactness of the inference results and highlights that statistical inference based on various ordinary least squares regressions could either be too conservative or invalid. Specifically, my result verifies that the commonly used regression with cluster-robust standard errors is conservative, while the other regression-based methods examined in the paper, such as regressions with strata fixed effects or heteroskedasticity-robust standard errors, are generally invalid. In the empirical application, I demonstrate the proposed inference method based on the experiment conducted in foos2017 and compare it with regression-based methods. The empirical findings are consistent with the results of the simulation study.
The analysis of data from two-stage randomized experiments and experiments under covariate-adaptive randomization has received considerable attention, but most work has focused on only one of these two features at a time. Previous work on the analysis of two-stage randomized experiments includes Hirano2010, liu2014, Rigdon2015, Baird2018, basse2018, basse2019, Imai2021, imai2022, VAZQUEZBARE2022 and Gonzalo-two-stage. Recent work on the analysis of covariate-adaptive experiments includes Bugni2018, yichong2021, bai-inference, Bai-optimal, mp-cluster, matched-tuple, yichong2022, bai2023, cytrybaum2021 and bai-attrition. In fact, both basse2018 and Imai2021 applied their inference methods, which do not account for covariate information, to two-stage experiments under covariate-adaptive randomization.\footnote{basse2018 analyzes the empirical application from Rogers2018, whose design involves stratification on school, grade, and prior-year absences. Imai2021 analyzes the empirical application from Kinnan2020, whose design involves matching villages (clusters) and households into small blocks.} My framework of analysis follows closely Bugni2022, in which they formalize cluster randomized experiments in a super population framework.
This paper contributes to the methodology for a growing number of empirical papers using two-stage experiments with covariate-adaptive randomization. For instance, Muralidharan2015, Hidrobo2016, foos2017, Rogers2018 and Banerjee2021 conducted two-stage randomized experiments that stratify clusters or units into a small number of large strata according to their baseline covariates, typically known as stratified block randomization. Duflo2003, Ichino2012, Beuermann2015, Kinnan2020 and Malani2021 conducted two-stage randomized experiments in which clusters or units are matched into small strata according to their baseline covariates, commonly known as matched pairs, matched triplets or matched tuples designs.
The rest of the paper is organized as follows. Section (ref) describes the setup and notation. Section (ref) presents the main results. Section (ref) discusses the optimality of matched tuples designs. Section (ref) introduces the covariate-adjusted estimator. Section (ref) examines the finite sample behavior of various experimental designs through simulations. Section (ref) illustrates the proposed inference methods in an empirical application based on the experiment conducted in foos2017. Finally, I conclude with recommendations for empirical practice in Section (ref).
Let $Y_{i,g}$ and $X_{i,g}$ denote the observed outcome and individual baseline covariates of the $i$th unit in the $g$th cluster, respectively. Denote by $Z_{i,g}$ the indicator for whether the $i$th unit in the $g$th cluster is treated or not. Let $C_g$ denote the observed baseline covariates for the $g$th cluster, $N_g$ denote the size of the $g$th cluster, $H_g$ denote the target fraction of units treated in the $g$th cluster, and $G$ the number of observed clusters. In addition, define $\mathcal{M}_g$ as the (possibly random) subset of $\{1,...,N_g\}$ corresponding to the observations within the $g$th cluster that are sampled by the researcher. Let $M_g=|\mathcal{M}_g|$ denote the number of units in set $\mathcal{M}_g$. In other words, the researcher randomly assigns treatments to all $N_g$ units in the $g$th cluster but only observes or conducts analysis on a subset of units sampled from the $g$th cluster Beuermann2015,Muralidharan2015,Haushofer2016,Hidrobo2016,Aramburu2019,Haushofer2019,Banerjee2021,Malani2021. Denote by $P_G$ the distribution of the observed data
This paper considers a setup where units are partitioned into a large number of clusters. In this context, the paper studies a two-stage randomized experiment with binary treatment in both stages. In the first stage, a fraction of $\pi_1$ clusters are randomly assigned to the treatment group, while the remaining clusters are assigned to the control group with no treated units. Then, conditional on the assignment in the first stage, a fraction of $\pi_2$ individuals from treated clusters are assigned to the treatment group, while the remaining units are assigned to the control group. Such a binary design is widely used in empirical literature Duflo2003,Ichino2012,Haushofer2016,foos2017,Haushofer2019. Moreover, while some experiments have multiple treated fractions, researchers often analyze them as binary designs Beuermann2015,basse2018,Imai2021.
The two-stage experiment closely resembles the split-plot design (see shi2022rerandomization,zhao2022splitplot), where \(H_g\) represents the whole-plot (cluster-level) randomization and \(Z_{i,g}\) represents the subplot (within-cluster) randomization. In split-plot designs, \(H_g\) and \(Z_{i,g}\) are usually assumed to be independent and represent two binary factors of treatment. However, in two-stage designs, \(H_g\) represents the intended treated fraction and thus does not correspond to a real treatment; it is correlated with \(Z_{i,g}\) through the relation \(H_g = \sum_{1 \leq i \leq N_g} Z_{i,g} / N_g\).\footnote{Strictly speaking, the equality holds up to a finite sample error, i.e. $\lfloor H_g N_g \rfloor = \sum_{1 \leq i \leq N_g} Z_{i,g}$.} This distinction indicates that it could be a promising direction for future research to develop a general framework that allows dependence between the first-stage and second-stage randomizations, encompassing both split-plot and two-stage designs.
In this section, I provide assumptions on the interference structure that assume no interference across clusters and exchangeable/homogeneous interference within clusters. Let $Y_{i,g}(\mathbf{z}, n)$ denote the potential outcome of the $i$th unit in the $g$th cluster, where $n$ denotes the cluster size and $\mathbf{z}$ denotes a realized vector of assignment for all units in all clusters, i.e., $\mathbf{z}=((z_{i,g}: 1\leq i\leq n): 1\leq g\leq G)$, where $z_{i,g} \in \{0,1\}$ denotes a realized assignment for the $i$th unit in the $g$th cluster. Following previous work hudgens2008, basse2018, basse2019, Forastiere2021, Imai2021, I assume the following about potential outcomes.
Under Assumption (ref), potential outcomes can be simplified as $Y_{i,g}(z,n,n_1)$ where $n_1$ denotes the number of treated units in the cluster. Following this notation, we define
to be the potential outcome under the individual treatment status $z \in \{0,1\}$ and the cluster target treated fractions $h \in \mathcal{H} \subseteq [0, 1]$, where $\mathcal{H}$ is a pre-determined set of treated fractions.\footnote{For example, when the cluster size is $3$ and the target treated fraction is $0.5$, there will be one treated unit in the cluster. Other rounding approaches, like the ceiling function, to handle fractional numbers of treated units can also be easily accommodated.} As mentioned before, this paper considers binary treatments, i.e. $\mathcal{H} = \{0, \pi_2\}$, throughout the paper.\footnote{Extending the designs to accommodate multiple treatment fractions is technically straightforward. Related work can be found in Bugni2019.} Furthermore, the (observed) outcome and potential outcomes are related to treatment assignment by the relationship $Y_{i,g} = Y_{i,g}(Z_{i,g}, H_{g})$. Denote by $Q_G$ the distribution of
The distribution $P_G$ of observed data and its sampling procedure can be described in three steps. First, $\{(\mathcal{M}_g, C_g, N_g): 1\leq g \leq G\}$ are i.i.d samples from a population distribution. Second, potential outcomes and baseline individual covariates are sampled from a conditional distribution $R_G(\mathcal{M}^{(G)}, C^{(G)}, N^{(G)})$, which is defined as follows:
Finally, $P_G$ is jointly determined by the relationship $Y_{i,g} = Y_{i,g}(Z_{i,g}, H_{g})$ together with the assignment mechanism, which will be described in Section (ref), and $Q_G$, which is described in the first two steps. Note that \(A^{(G)}\) denotes the vector \((A_1, \dots, A_G)\) for any random variable \(A\), and \(X_g\) represents the vector \((X_{i,g}: 1 \leq i \leq N_g)\). The following assumption states my requirements on $Q_G$ using this notation.
The sampling procedure of a cluster randomized experiment used in this paper closely follows that formalized by mp-cluster and Bugni2022. Assumption (ref) is essentially the same as Assumption 2.2 in Bugni2022, which formalizes the sampling procedure of i.i.d. clusters (Assumptions (ref) (a)-(b)) and imposes mild regularity conditions (Assumptions (ref) (c)-(d)). Furthermore, Assumption (ref) (e) accommodates a second-stage sampling process within a given cluster that may depend on cluster-level and individual-level covariates as well as cluster sizes. This flexibility permits $\mathcal{M}_g$ to be potentially determined through stratified sampling, as discussed in cytrybaum2021. Finally, Assumption (ref) (f) is a high-level assumption that ensures the extrapolation from the observations that are sampled to those that are not sampled.
In the context of the sampling framework described above, this paper considers four parameters of interest, including primary and spillover effects that are equally or (cluster) size-weighted. For different choices of (possibly random) weights $\omega_g$, $1\leq g\leq G$ satisfying $E[\omega_g]=1$, we define the average primary effects and spillover effects under general weights as follows.
Denote by $\theta^{P}_1(Q_G)$ and $\theta^{S}_1(Q_G)$ the equally-weighted cluster-level average primary and spillover effects with $\omega_g=1$, and $\theta^{P}_2(Q_G)$ and $\theta^{S}_2(Q_G)$ the size-weighted cluster-level average primary and spillover effects with $\omega_g=N_g/E[N_g]$. The consideration of weighted estimands is motivated by the non-ignorability of cluster sizes. According to Bugni2022, cluster sizes are considered ignorable if the individual-level average treatment effect is independent of the cluster size. Formally, this is expressed as:
for all \(1\leq g \leq G\) and \(z \in \{0, 1\}\). Cluster sizes are non-ignorable whenever ((ref)) is not satisfied. When cluster sizes are non-ignorable, different weights can lead to distinct parameters. The selection between these two types of estimands—equally weighted or size-weighted—depends on the analytical focus: whether the primary interest is on the clusters themselves or the individuals within these clusters. For instance, in assessing the impact of an educational program on students' academic performance, if policymakers are concerned with improvements at the school level, equally weighted estimands are appropriate. Conversely, if the focus is on student-level outcomes, then size-weighted estimands become relevant.
The primary effects $\theta^{P}_1(Q_G)$ and $\theta^{P}_2(Q_G)$ are the differences in the averaged potential outcomes of treated units from treated clusters and control units from control clusters. In contrast, the spillover effects $\theta^{S}_1(Q_G)$ and $\theta^{S}_2(Q_G)$ are the differences in the averaged potential outcomes of control units from treated clusters and control units from control clusters. In many empirical settings, the estimation and comparison of primary and spillover effects play a crucial role in addressing important research questions Duflo2003.
In summary, the formulas for the four parameters of interest are listed in Table (ref). These estimands have been proposed and studied in previous literature hudgens2008, toulis2013, basse2018, Imai2021, but mostly in a finite population framework. This paper adopts the terminology “primary” and “spillover” effects from basse2018, which are respectively referred to as “total” and “indirect” effects in hudgens2008. Previous works on interference have also studied other estimands, such as direct effects and overall effects hudgens2008,Wager2021,Imai2021, but I do not explore these estimands further in this paper.
For estimating the four parameters of interest, I propose the following estimators analogous to the difference-in-“average of averages” estimator in Bugni2022:
where $G_T = \sum_{1\leq g \leq G} I\{H_g = \pi_2\}$, $G_C = \sum_{1\leq g \leq G} I\{H_g = 0\}$, and $N_T = \sum_{1\leq g \leq G} I\{H_g = \pi_2 \} N_g$, $N_C = \sum_{1\leq g \leq G} I\{H_g = 0 \} N_g$ and
where $M_g^z = \sum_{i\in\mathcal{M}_g} I\{ Z_{i,g} = z \}$ with $z \in \{0,1\}$.
By definition, the “first/individual average” $\bar{Y}_{g}^{1}$ from the primary effect estimator is taken over all treated units within the $g$-th cluster if the cluster is treated, and all control units within the $g$-th cluster if the cluster is assigned to control. When it comes to estimating spillover effects, the “first/individual average” $\bar{Y}_{g}^{0}$ is taken over all control units within the $g$-th cluster if the cluster is treated, and all control units within the $g$-th cluster if the cluster is assigned to control. Then, the “second/cluster average” is a cluster-level average of $\bar{Y}_{g}^{1}$ or $\bar{Y}_{g}^{0}$ taken within groups of treated and untreated clusters as featured in a usual difference-in-means estimator.
The proposed estimators can be obtained from ordinary least squares regressions using different weighting schemes. Let \( L_{i,g} = I\{H_g = \pi_2\}(1 - Z_{i,g}) \) denote the indicator for untreated units within treated clusters. Consider the following linear model for an ordinary least squares regression:
Note that the estimators \(\hat{\theta}_1^P\) and \(\hat{\theta}_1^S\) may be obtained by estimating coefficients \(\beta_1\) and \(\beta_2\) from a weighted least squares regression of equation ((ref)) using weights \(1/M_g\). Similarly, \(\hat{\theta}_2^P\) and \(\hat{\theta}_2^S\) may be derived using weights \(N_g/M_g\) (see Appendix (ref) for formal derivations). Moreover, the unweighted least squares regression produces the “sample” size-weighted estimators. Taking \(\beta_1\) as an example:
where \( M_1 = \sum_{1 \leq g \leq G} I\{H_g = \pi_2\} M_g \), and \( M_0 = \sum_{1 \leq g \leq G} I\{H_g = 0\} M_g \). These “sample” size-weighted estimators are identical to \(\hat{\theta}_2^P\) and \(\hat{\theta}_2^S\) when outcomes of all units from each cluster are observed or the number of observed units is proportional to the cluster size, i.e., \( M_g / N_g = c \) for \( 0 < c \leq 1 \).
My estimators are closely related to those studied in previous methodological literature. For example, equally-weighted estimators $\hat{\theta}^P_1$ and $\hat{\theta}^S_1$ are identical to the household-weighted estimators from basse2018, which are closely related to the estimators in hudgens2008. $\hat{\theta}^P_1$ and $\hat{\theta}^S_1$ may also be obtained through the “household-level regression” proposed in basse2018, which is equivalent to running two separate ordinary least squares regressions of $\bar Y_{g}^1$ on a constant and $I\{H_g=\pi_2\}$, and $\bar Y_{g}^0$ on a constant and $I\{H_g=0\}$. Size-weighted estimators $\hat{\theta}^P_2$ and $\hat{\theta}^S_2$ are closely related to the individual-weighted estimator proposed by basse2018. In previous studies such as basse2018, VAZQUEZBARE2022, and Gonzalo-two-stage, researchers have investigated estimators obtained through a widely used saturated regression in multi-treatment experiments, similar to the least squares regression described by equation ((ref)).
In empirical literature, various regression estimators are used for estimating primary and spillover effects. One widely used estimator is described in equation ((ref)) Haushofer2016,Haushofer2019. Another estimator that produces the same set of estimators is through the alternative regression $Y_{i,g} = a + b_1 Z_{i,g} + b_2 I\{H_g=\pi_2\} + u_{i,g}$ Duflo2003,Ichino2012, where the estimators are related to those from ((ref)) as follows: $\hat\beta_1 = \hat b_1 + \hat b_2$ and $\hat \beta_2 = \hat b_2$. Some empirical works use either or both of the two separate regressions: $Y_{i,g} = \alpha + \beta_1 Z_{i,g} + \epsilon_{i,g}$ and $Y_{i,g} = \alpha + \beta_2 L_{i,g} + \epsilon_{i,g}$ Beuermann2015,Hidrobo2016,Aramburu2019. In many cases, estimators obtained from regressions with fixed effects are reported along with those without fixed effects Ichino2012. Section (ref) will examine the validity of statistical tests based on regressions with and without fixed effects.
In this section, I investigate the asymptotic properties of the estimators presented in Section (ref) within a finely stratified randomization framework. Specifically, in the first stage, clusters are partitioned into a large number of small strata of a fixed size, with the assignment mechanism being a completely randomized design (also known as a permuted block design) independently applied within each stratum. Formally, consider \( n \) strata of size \( k \) (each stratum consisting of \( k \) clusters), formed by matching clusters according to a function \( S:\text{supp}((C_g,N_g)) \rightarrow \mathbf{R}^\ell \). Denote by \( S^{(G)}=(S_1, \dots, S_G) \) the vector of variables used for matching, where \( S_g = S(C_g, N_g) \). Within each stratum, \( l \) clusters are randomly selected and assigned to the treatment group.\footnote{Extending the setup to a more general framework with varying stratum sizes and heterogeneous treatment fractions is indeed possible; see Section 3.2 of cytrybaum2021.} Specifically, \( G = nk \) and \( \pi_1 = l/k \), where \( 0 < l < k \), and \( l \) and \( k \) are mutually prime. Furthermore, I consider a second-stage stratification on units from a given cluster. Denote by \( B_g = (B_{i,g}: 1 \leq i \leq N_g) \) the vector of strata on units in the \( g \)th cluster, constructed from observed baseline covariates \( X_{i,g} \) for the \( i \)th unit using a function \( B:\text{supp}(X_{i,g}) \rightarrow \mathcal{B}_g \).\footnote{Asymptotics are not considered in the second-stage design; thus, the second stage could employ finely stratified designs like matched-pair, or those with coarse stratification such as stratified block randomization.}
To start with, I describe my assumptions on the treatment assignment mechanism in the first stage. Formally, let
denote $n$ sets each consisting of $k$ elements that form a partition of $\{1,\dots, G\}$.
I assume treatment status is assigned to clusters as follows:
Assumption (ref) formally describes the assignment mechanism of a two-stage experiment with finely stratified randomization in the first stage. Further, units in each pair are required to be “close” in terms of their stratification variable $S_g$ in the following sense:
The validity of the variance estimators relies on the following condition that the distances between units in adjacent blocks are considered “close” in relation to their baseline covariates:
The next step is to formalize the assumption of independence between the first and second stage designs. To begin with, I utilize the notation $\{Z_{i,g}(h): h \in \mathcal{H}\}$, representing the “potential treatment” for various treated fractions $h\in\mathcal{H}$, and relate the (observed) individual treatment indicator and potential individual-level treatment indicator as follows:
The underlying motivation for this “potential outcome style” notation becomes evident when considering that in two-stage experiments, the realized treatment assignment in the first stage is almost always correlated with that in the second stage (e.g., $H_g = \frac{1}{N_g}\sum_{1\leq i \leq N_g} Z_{i,g}$). Yet, the “potential” individual-level treatment assignment, for any specified target treated fraction, can be independent of the cluster-level assignment of that target treated fraction. This is similar to the classic potential outcome model, where treatment assignment is independent of potential outcomes but likely correlates with observed outcomes.
Then, my requirements on the treatment assignment mechanism for the second stage are summarized in the following assumption:
Assumption (ref) (a) rules out any confounders between the first-stage and second-stage treatment assignments, which is typically satisfied in most two-stage experiments. Assumption (ref) (b) is analogous to Assumption (ref) (a). Assumption (ref) (c) requires that the marginal assignment probability for each stratum and the realized treated fraction in the observed subset of units both equal the intended treated fraction $h$, up to a finite sample error that diminishes as cluster size increases.\footnote{In the proof of the main results, I only need \( E[Z_{i,g}(h) \mid B_g] = \frac{1}{M_g} \sum_{i \in \mathcal{M}_g} Z_{i,g}(h) \) to hold for unbiasedness. However, in practice, these two quantities, along with the treated fraction for the entire cluster, need to align with the intended treated fraction \( h \) so that they are consistent with the notations of the potential outcomes \( Y_{i,g}(z,h) \).} An example of this could be (individual-level) stratified block randomization, where the treated fraction remains constant across all strata, with observed units drawn from a random subset of these strata.
Finally, I impose the following assumption on $Q_G$ in addition to Assumption (ref):
Assumption (ref)(a) is a smoothness requirement analogous to Assumption 3(ii) in Bai-optimal ensuring that units within clusters which are “close” in terms of their baseline covariates are suitably comparable. Assumption (ref)(b) imposes an additional restriction on the distribution of cluster sizes beyond what is stated in Assumption (ref)(c).
The following theorem derives the asymptotic behavior of estimators for equally-weighted and size-weighted effects.\footnote{Throughout the paper, $V_1(1)$ and $V_2(1)$ denote the variances of primary effects, while $V_1(0)$ and $V_2(0)$ represent the variances of spillover effects. In other words, the notation $z\in \{0, 1\}$ (as in $V_1(z)$) represents the individual's own treatment status.}
Theorem (ref) implies that covariate information is important to establish asymptotically exact inference for the four estimands of interest under covariate-adaptive randomization. Many empirical studies rely on statistical inference based on the regression in equation ((ref)) with HC2 cluster-robust standard errors. While this procedure is also proposed in basse2018 and Gonzalo-two-stage, the regression coefficients it produces generally do not provide consistent estimates for the estimands in Table (ref). As discussed in Section (ref), if all units in each cluster are sampled ($N_g = M_g$) or the number of sampled units is proportional to cluster size ($M_g/N_g=c$ for $0< c < 1$), this procedure yields consistent point estimates for size-weighted effects but may still be conservative (see Appendix (ref)). Therefore, I aim to develop asymptotically exact inference methods based on my theoretical results.
To begin with, I introduce consistent variance estimators for the asymptotic variances from Theorem (ref). To estimate $V_1(z)$ and $V_2(z)$, I follow the construction of “pairs of pairs” in bai-inference and matched-tuple, and replace the individual outcomes with the averaged outcomes $\bar Y_g^{z}$ (as defined in Section (ref)) and adjusted averaged outcomes $\tilde{Y}_g^z$ , respectively. The definition of the adjusted average outcomes is given as follows:
where $G_g = \sum_{1\leq j \leq G} I\{H_g = H_j\}$. Here, I present the construction of variance estimator $\hat V_1(z)$ for $V_1(z)$. Similarly, $\hat V_2(z)$ can be constructed by simply replacing $\bar Y_{g}^z$ with $\tilde Y_{g}^z$, and thus details are omitted. Let $\hat \Gamma^z_n(h) = \frac{1}{n k(h)} \sum_{1\leq g \leq G} \bar Y_g^z I \{H_g = h\}$ where $k(h) = \sum_{i\in\lambda_j} I\{H_i=h\}$ denotes the number of units under assignment $H_i = h$ in the $j$-th strata. In the setup of binary treatment, it becomes that $k(\pi_2)=l$ and $k(0)=k-l$. Finally, my estimator for $V_1(z)$ is then given by
with
where
Based on the variance estimators, I propose the “adjusted” $t$-test with the aforementioned variance estimators as my method of inference throughout the rest of the paper. As an example, the “adjusted” $t$-test for equally-weighted primary effect, i.e. $H_0: \theta_1^P(Q_G) = \theta_{0}$, is given by
where $z_{1-\frac{\alpha}{2}}$ represents $1-\frac{\alpha}{2}$ quantile of a standard normal random variable.
The subsequent analysis yields consistency results for the estimators $\hat V_1(z)$ and $\hat V_2(z)$ and validity results for the adjusted $t$-test:
Note that the variance estimator $\hat V_1(z)$ (or $\hat V_2(z)$) depends on the assignment mechanism in the first stage through the strata indicator $S_g$, but not on the assignment mechanism in the second stage. This means that valid statistical inference based on $\phi_G(V^{(G)})$ does not require knowledge of the assignment mechanism in the second stage. We can see this by observing that the first term in equations ((ref)), which is the only term affected by the second-stage design, can be consistently estimated by the first term in equation ((ref)). My approach leverages the cluster-level averaged outcomes and benefits from large samples of clusters, without explicitly modeling intra-cluster correlations as done in the previous literature Gonzalo-two-stage.
In this section, I introduce two optimality results related to two-stage randomized experiments, as discussed in Sections (ref). The first result provides insights into the optimal design for the initial stage, while the second addresses the optimal design for the second stage, taking into account additional assumptions about the assignment mechanism and covariance among unit outcomes within clusters. These findings indicate that particular finely stratified designs maximize statistical precision when estimating parameters outlined in Table (ref).
First, I present a result that identifies the optimal functions for matching in the first-stage, targeting various parameters of interest.
A direct implication of Theorem (ref) is that it characterizes the optimal functions to match on within the class of finely stratified designs. These functions are referred to as “index function” in Bai-optimal. As noted in Remark (ref), when $S_g$ is categorical, finely stratified designs correspond to stratified block randomization, which implies that the optimal finely stratified designs is also asymptotically optimal among all large strata designs described in Appendix (ref). It is important to note that when discussing the optimal design for the first stage, we are comparing different first-stage designs for any fixed second-stage design (and vice versa for the second-stage design optimality).
The subsequent discussion examines the optimality of finely stratified designs in the second stage of the experiment. The second-stage randomization is formalized in the following assumptions.
Additionally, I assume that the covariance of outcomes between any pairs of units within a cluster is homogeneous. In other words, the covariance does not depend on the individual-level covariates of units in the same cluster. Formally, the assumption is stated as follows:
Assumption (ref) is a weaker assumption than assuming that outcomes of units are independent and identically distributed (i.i.d) within a cluster, as it only requires conditional independence between individual covariates and the covariance of outcomes. It is analogous to the standard homoscedasticity assumption, which assumes constant variance of errors in a regression model, except that it is a statement about covariance instead of variance. Under these two additional assumptions I obtain the following optimality result:
Though practitioners may not have knowledge of the index functions in Theorem (ref) and (ref), optimal stratification can be determined in some special cases. For instance, in experiments where the first-stage design uses only a univariate covariate $C_g$ Ichino2012, and practitioners expect a monotonic relationship between $S_g$ and $C_g$, the optimal stratification is to order the units by $C_g$ and group adjacent units. Similar results apply to the second-stage design. In more general cases where monotonicity does not hold or the baseline covariates are multivariate, a suitable matching algorithm bai-inference, cytrybaum2021 that directly matches on vectors of covariates can be asymptotically as efficient if the sample size is sufficiently large. In cases where the sample size is not sufficiently large, McKenzie2009 and Bai-optimal suggest matching on the baseline outcome, when available. If none of the aforementioned options is available, matching in a sub-optimal way can still be effective, as both Bai-optimal and simulation results from Section (ref) demonstrate that matching units sub-optimally can be more effective than completely randomized designs or some sub-optimal stratified block randomization designs. In this case, it could be useful to consider the recommendations in Remarks (ref) and (ref) for the choice of covariates.
In the empirical literature, it is common to match or stratify on a small set of covariates expected to be most predictive of outcomes, and to adjust for additional pre-treatment covariates ex-post. Consequently, this section introduces a linearly covariate-adjusted modification of $\hat\theta_2^P$, the size-weighted primary effect estimator. Adjusted estimators for other estimands follow a similar methodology and are thus omitted for brevity.
To begin, I introduce a new set of baseline covariates \(L_g\) that were not used for treatment assignment. These covariates \(L_g\) may include cluster-level aggregates of individual-level outcomes, such as intracluster means and quantiles. For the remainder of Section (ref), the assumptions specified in Section (ref) are modified such that \(C_g\) is replaced by \((C_g, L_g)\) throughout. In particular, references to Assumption (ref) should now be considered to include \((C_g, L_g)\) instead of \(C_g\). Following this, the treatment status is assigned as follows:
I consider a linearly covariate-adjusted estimator based on a set of regressors generated by $C_g, N_g, L_g$. To this end, define $\psi_g = \psi( C_g, N_g, L_g)$, where $\psi: \text{supp}((C_g,N_g,L_g)) \to \mathbf R^p$. We impose the following assumptions on $\psi$:
I extend the covariate-adjusted estimator from mp-cluster to accommodate the finely stratified design with a general treatment fraction $\pi_1$, as discussed in this paper. Let $\hat\mu_{1,j}$ represent the averaged value of $\tilde Y_g^1 \bar N_G$ among treated clusters within the $j$-th tuple, i.e., $g \in \lambda_j$. Similarly, $\hat\mu_{0,j}$ denotes the corresponding value for control clusters. Additionally, $\hat\psi_{1,j}$ and $\hat\psi_{0,j}$ refer to the averaged values of $\psi_g$ for treated and control clusters, respectively. Formaly, define
where $\bar N_G = \sum_{1\leq g \leq G} N_g/G$. Then, I define the linear adjustment coefficient $\hat\beta_2^P$ as the ordinary least squares (OLS) estimator of the slope coefficient in the linear regression of $\hat\mu_{1,j} - \hat\mu_{0,j}$ on a constant and $\hat\psi_{1,j} - \hat\psi_{0,j}$. Finally, I introduce my covariate-adjusted estimator for the size-weighted primary treatment effect as follows:
where
The following theorem derives the asymptotic behavior of my covariate-adjusted estimator for $\theta_2^P$, and, importantly, it shows that the limiting variance of $\hat \theta_2^{P, adj}$ is no larger than that of $\hat \theta_2^{P}$ in Theorem (ref) and can be strictly smaller.
For variance estimation, I employ the same methodology as \(\hat{V}_2(z)\) but with a modification: \(\tilde{Y}_g^z\) is replaced by \(\mathring{Y}_g^z = \tilde{Y}_g^z - \frac{(\psi_g - \bar{\psi}_G)' \hat\beta_2^P}{ \frac{1}{G}\sum_{1 \leq g \leq G} N_g }\). The consistency of this variance estimator follows from combining the arguments used to establish Theorem (ref) and those used to establish Theorem 3.2 in bai2023.
In this section, I illustrate the results presented in Section (ref) with a simulation study. To begin with, potential outcomes are generated according to the equation:
for $(z,h) \in \{(0,0),(0,\pi_2),(1,\pi_2)\}$, where
All simulations are performed with a sample of $200$ clusters, in which all units are sampled, i.e. $N_g=M_g$.
This section examines the performance of optimal matched tuples designs and several other designs via comparison of their MSEs (Mean Squared Errors). For simplicity, the parameters are given as follows: $\alpha_{z,h} = \beta_{z,h} = 1, \gamma = 1/100$ for all $(z,h)\in \{(0,0),(0,\pi_2),(1,\pi_2)\}$. This model configuration is referred to as “homogeneous model” since treatment effects are fully captured by $\mu_{z,h}$ and thus are homogeneously additive in this setting. A more complicated “heterogeneous model” will be introduced later. According to Theorem (ref), the optimal index functions for equally-weighted and size-weighted effects in the first stage are
In the second stage, the optimal finely stratified design matches on $X_{1,i,g}/(X_{2,i,g}+0.1)$ according to Theorem (ref). This section considers the following experimental designs for both stages:
Table (ref) shows the ratio of the MSE of each design relative to the MSE of the design with completely randomized assignments (C) in both stages, computed across 1000 Monte Carlo iterations. The rows indicate first-stage designs, and columns indicate second-stage designs. The lowest values in each row are marked in bold. In all designs, treatment effects are set to zero by assigning $\mu_{z,h} = 0$ for all $(z,h)\in {(0,0),(0,\pi_2),(1,\pi_2)}$, and the treated fraction is set to $1/2$ in both stages. As expected from Theorem (ref) and (ref), the matched-tuples design with complete matching (MT-C) outperforms the other designs in the first stage for all parameters of interest while remaining optimal in the second stage for many cases. However, it is noticeable that the assignment mechanism in the first stage has a greater effect on statistical precision than the second stage.
In this section, the focus shifts from optimality to studying the finite sample properties of different tests for the following null hypotheses of interest:
against the alternative hypotheses:
In Table (ref), the six assignment mechanisms with covariate-adaptive randomization (Design 2-7 in Section (ref)) for the first and second stages are considered, resulting in a total of 36 different designs. Hypothesis tests are performed at a significance level of 0.05, and rejection probabilities under the null and alternative hypotheses are computed from 1000 Monte Carlo iterations in each case. Tests are constructed as “adjusted $t$-tests” using the asymptotic results from Theorem (ref)-(ref). For stratified designs in the first stage (S-2, S-4 and S-4O), tests for equally- and size-weighted effects are performed using the variance estimators $\hat V_3(z)$ and $\hat V_4(z)$ (see ((ref)) and ((ref)) in Appendix (ref)). For matched tuples designs in the first stage (MT-A, MT-B and MT-C), tests for equally- and size-weighted effects are performed using the variance estimators $\hat V_1(z)$ and $\hat V_2(z)$. The results show that the rejection probabilities are universally around 0.05 under the null hypothesis, which verifies the validity of tests based on my asymptotic results across all the designs. Under the alternative hypotheses, the rejection probabilities vary substantially across the first-stage designs while remaining relatively stable across the second-stage designs. \textbf{MT-C} stands out as the most powerful design for the first-stage. These findings are consistent with previous section.
Next, the validity of commonly used regression-based inference methods in the empirical literature is tested. These methods are tested under both the “homogeneous model” from the previous simulation study in Section (ref) and a “heterogeneous model” in which two parameters are modified as follows: $\alpha_{1,\pi_2} = \beta_{1,\pi_2} = 2$, $\alpha_{0,\pi_2} = \beta_{0,\pi_2} = 0.5$, and $\alpha_{0,0} = \beta_{0,0} = 1$. The key difference between the two models is whether the conditional expectations of potential outcomes are identical or different across different exposures $(z,h)$. Four commonly used regression methods are considered in this study:
Note that due to full sampling, i.e. $N_g = M_g$, regressions without fixed effects (“OLS robust” and “OLS cluster”) output the same estimators as the size-weighted estimators $\hat \theta_2^P$ and $\hat \theta_2^S$. Most of the previous empirical analysis on covariate-adaptive two-stage experiments report cluster-robust standard errors in their main results, which could either be “OLS cluster” basse2018 or “OLS with group fixed effects (clustered)” Duflo2003,Ichino2012. For brevity, Table (ref) includes only six designs: those with either S-4O or MT-C in the first stage, and C, S-4O, or MT-C in the second stage. The table reveals that test results can be either conservative or invalid across different regression methods and designs. For stratified designs in the first stage, methods based on “robust” standard errors tend to over-reject, while methods based on “clustered” standard errors tend to under-reject. For matched tuples designs, “OLS cluster” is conservative, and the remaining methods could be invalid as they may over-reject the null hypothesis under some model specifications and parameters of interest. Similar results can also be found in the previous literature on covariate-adaptive randomization. For example, matched-tuple demonstrated that inferences based on OLS regressions with strata fixed effects could be invalid. On the other hand, dechaisemartin2022 documented that in cluster randomized experiments, $t$-test based on clustered standard errors tend to over-reject the null hypothesis when strata fixed effects are included, and under-reject otherwise. Therefore, it can be concluded that, with the exception of “OLS cluster” being conservative, the other three inference methods based on regression are generally invalid.
In this section, the inference methods introduced in Section (ref) are illustrated using data collected in foos2017. The experiment conducted by foos2017 is a randomly assigned spillover experiment in the United Kingdom designed to identify social influence within heterogeneous and homogeneous partisan households. The study first stratified $5190$ two-voter households into three blocks based on the latest recorded party preference of the experimental subject\footnote{Before assigning treatments, the researchers randomly selected one individual per household to potentially receive treatments, whom they mark as “experimental subjects”. In other words, the second-stage assignment is a complete randomization. Specifically, this two-stage design corresponds to “S3-C” (using the notation from the simulation section).}: “Labour” supporter,“rival party” supporter and those who were “unattached” to a party. Then experimental subjects or equivalently their households were randomly assigned to three groups: high partisan intensity treatment, low partisan intensity and control\footnote{The empirical treatment fractions for “Labour” supporters are 0.217, 0.217, and 0.566 for the high-intensity, low-intensity, and control groups, respectively. For “rival party” supporters, the corresponding fractions are 0.222, 0.215, and 0.563. For “unattached” individuals, they are 0.208, 0.226, and 0.566.}. Experimental subjects allocated to treatment groups were called by telephone and encouraged to vote in the PCC election on November 15, 2012. The “high partisan intensity” was formulated in a strongly partisan tone, explicitly mentioning the Labour Party and policies multiple time, while the “low partisan intensity” treatment message avoided all statements about party competition.
In the original analysis of foos2017, their main focus was on analyzing treatment effects conditional on a wide range of pre-treatment covariates. That said, in the final column of Table 1 in foos2017, they report estimators for (unconditional) primary and spillover effects, which are based on calculations of averages over separate experimental subjects and unassigned subjects. In contrast, my estimators do not distinguish experimental subjects from unassigned subjects and take averages solely based on treatment or spillover status. Another difference in my analysis is that estimators are calculated by pooling the two treatment arms, i.e. high and low partisan intensity, to maintain consistency with the setup of the paper\footnote{ Specifically, treated households effectively received a “random treatment”: high partisan intensity with some probability and low partisan intensity with the complementary probability. The pooled treatment still follows a complete randomization design within each stratum and therefore satisfies all assumptions related to treatment assignment. } . In contrast, foos2017 provide separate estimates for each treatment arm.
Table (ref) compares point estimates of treatment effect on turnout percentage and confidence intervals obtained from the four regression methods listed in Section (ref) with those based on my theoretical results, namely “adjusted $t$-test”. Since cluster (household) size is fixed, equally-weighted and size-weighted estimators and estimands collapse into one. Moreover, full sampling ($N_g = M_g = 2$) makes the point estimates of “adjusted $t$-test” and “OLS robust/cluster” equivalent. In the simulation study, it is found that `OLS robust” and “OLS fe robust” tend to over-reject the null hypothesis, which is consistent with the empirical results in Table (ref) that they both have narrower confidence interval than the “adjusted $t$-test”. Furthermore, “OLS cluster” and “OLS fe cluster” are shown to be conservative in the simulation study, and accordingly, they both have wider confidence intervals than the “adjusted $t$-test” in Table (ref). Therefore, the empirical findings are consistent with the simulation study in Table (ref).
Based on the theoretical results and the supporting simulation study, I conclude with the following recommendations for empirical practice, particularly in conducting inference about the parameters of interest, as listed in Table (ref). In scenarios where sizes of all strata are considerably large, such as more than 50 clusters as exemplified in simulation S-4, we advise practitioners to utilize $\hat V_3(1)$ and $\hat V_3(0)$, as defined in ((ref)), for estimating the equally-weighted primary effect $\theta_1^P$ and the spillover effect $\theta_1^S$. Similarly, $\hat V_4(1)$ and $\hat V_4(0)$, as detailed in ((ref)), should be employed for the size-weighted primary effect $\theta_2^P$ and the spillover effect $\theta_2^S$. However, when it is unclear whether the strata size is sufficiently large, or more commonly, when the experimental design involves a matched-tuples design with only one or two observations per treatment arm, we recommend the application of $\hat V_1(1), \hat V_1(0)$ and $\hat V_2(1), \hat V_2(0)$ as indicated in ((ref)) for the corresponding equally-weighted and size-weighted effects.
The results of this study have shown that tests based on the regression specified in equation ((ref)) with HC2 cluster-robust standard errors are valid but potentially conservative, which would result in a loss of power relative to our proposed test. Further, it's critical to note that regressions using strata fixed effects or heteroskedasticity-robust standard errors have generally been found invalid in the simulation study.
Based on the optimality results for the first-stage design, I recommend selecting cluster-level covariates for matching according to the parameters of interest, as elaborated in Remark (ref), while adhering to the established guidelines from previous studies McKenzie2009, bai-inference, Bai-optimal, cytrybaum2021. For the second stage, it is advisable to first evaluate the impact of the design on efficiency, as detailed in Remark (ref), and then assess whether the benefits of second-stage randomization outweigh its costs. Should this be the case, implementing a finely stratified second-stage randomization is recommended, taking into account intra-cluster correlation, as discussed in Remark (ref).