Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
92,655 characters · 18 sections · 68 citation commands
Covariate Adjustment in Stratified Experiments
{Keywords: Matched Pairs, Analysis of Covariance, Blocking, Robust Standard Error, Treatment Effects.} \\
{JEL Codes: C10, C14, C90}
\onehalfspacing
This paper studies covariate adjusted estimation of the average treatment effect (ATE) in stratified experiments. Researchers often make use of both stratified treatment assignment and ex-post covariate adjustment to improve the precision of experimental estimates. Indeed, out of a survey of over $50$ experimental papers published in the AER and AEJ between 2018-2023, we found that $57\%$ use stratified randomization, and $80\%$ used some form of ex-post covariate adjustment. An influential paper by lin2013 showed in a design-based setting that the regression estimator with full treatment-covariate interactions is always asymptotically weakly more efficient than difference of means estimation for completely randomized designs. negi2021 extended these results to estimation of the ATE using data sampled from a superpopulation. However, questions remain about the interaction between stratification and regression adjustment and the implications of combining these methods for both estimator efficiency and the power and validity of inference methods. To study these questions, we work in the stratified randomization framework of cytrynbaum2023, which includes matched tuples designs (e.g. matched pairs), coarse stratification, and complete randomization as special cases.
We show that the lin2013 interacted regression adjustment is generically inefficient in the family of linearly adjusted estimators, with asymptotic efficiency only in the limiting case of complete randomization. Motivated by this finding, we characterize the efficient linear covariate adjustment for a given stratified design, providing several new estimators that achieve the optimal variance.
Our first result derives the optimal linear adjustment coefficient for a given stratification. We show that asymptotically the interacted regression estimator uses the wrong objective function, minimizing a marginal variance objective that is totally insensitive to the stratification. By contrast, the optimal adjustment coefficient minimizes a mean-conditional variance objective, conditional on the covariates used to stratify. Intuitively, the efficient covariate adjustment is tailored to the stratification, ignoring fluctuations of the estimator that are predictable by the stratification covariates. Section (ref) draws an interesting connection with partially linear regression (robinson88), showing that efficient linear adjustment of a stratified design is asymptotically equivalent to doubly-robust semiparametric adjustment of an iid\ design. Intuitively, stratification contributes the nonparametric component of the semiparametric adjustment function.
Our second set of results develops feasible versions of the optimal linear adjustment derived in Section (ref). First, we show that if the conditional expectation of the adjustment covariates is linear in a known set of transformations of the stratification variables, then adding the latter to the interacted regression restores optimality. Next we relax this assumption, providing four different regression estimators that are asymptotically efficient under weak conditions. For matched pairs experiments or in settings with limited treatment effect heterogeneity, the non-interacted regression with a full set of pair fixed effects is asymptotically efficient. More generally, we show asymptotic optimality of within-stratum (inconsistently) partialled versions of the Lin and tyranny-of-the-minority estimators (lin2013). We also define a “group OLS” estimator, extending a proposal of imbens2015 for matched pairs experiments to a larger class of designs. We show that this group OLS estimator is also asymptotically optimal.
Our final contribution is to develop novel asymptotically exact inference methods for covariate adjusted estimation under stratified designs. Confidence intervals based on the usual heteroskedasticity robust variance estimator are known to be conservative in stratified experiments (bai2021inference). By contrast, the coverage probabilities of our proposed confidence intervals converge to the specified nominal level, with no overcoverage. Our approach applies to a generic family of linear covariate adjustments and randomization schemes, including as special cases non-interacted regression adjustment, the lin2013 interacted regression, and all of the other estimators considered in this paper. Simulations and an empirical application to the experiment in baysan2022 suggest that the usual robust confidence intervals can substantially overcover in stratified experiments, while our confidence intervals have close to nominal coverage.
We present several extensions of our main results in the appendix. In the first, we consider estimation and inference in stratified experiments with noncompliance. As a simple corollary of our results on ATE estimation, we characterize the optimal linearly adjusted Wald estimator for the LATE (imbens1994), construct a feasible implementation of the efficient adjustment, and provide asymptotically exact inference methods. We also study efficient linear adjustment for finely stratified designs with non-constant treatment proportions, as in cytrynbaum2023, and briefly consider the problem of efficient nonlinear adjustment.
There has been significant interest in treatment effect estimation under different experimental designs in the recent literature. Some papers studying covariate adjustment under stratified randomization include bugni2018inference, fogarty2018b, liu2020, lu2022, ma2020, reluga2022, wang2021, ye2022, zhu2022, and chang2023. These works differ from our paper in at least one of the following ways: (1) studying inference on the sample average treatment effect (SATE) rather than the ATE in a superpopulation, (2) restricting to coarse stratification (stratum size going to infinity), or (3) proving weak efficiency gains but not optimality. In a finite population setting, zhu2022 shows asymptotic efficiency of a projection-based estimator numerically equivalent to the “partialled Lin” approach considered in Section (ref). In the same setting, lu2022 prove efficiency of a tyranny-of-the-minority style regression similar but not equivalent to one the considered in Section (ref). Both papers give conservative inference on the SATE, while we provide asymptotically exact inference on the ATE using a generalized pairs-of-pairs (abadie2008) style approach. Remarks (ref) and (ref) in Section (ref) below provide a detailed comparison.
Relative to the above papers, the superpopulation framework considered here creates some new technical challenges. For example, as pointed out in bai2021inference, matching units into data-dependent strata post-sampling produces a complicated dependence structure between the treatment assignments and random covariates. We deal with this using a tight-matching condition (Equation (ref)) and martingale CLT analysis similar to cytrynbaum2022local. This setting also has analytical advantages, which allow us to establish new conceptual results. For example, the population level characterization of the optimal adjustment coefficent in Section (ref) allows us to give explicit necessary and sufficient conditions for the efficiency of several commonly used regression estimators. The efficiency of interacted regression under a “rich covariates” condition, as well as the equivalence between optimal linear adjustment of stratified designs and doubly-robust semiparametric adjustment appear to be new observations in this literature. To the best of our knowledge, we give the first asymptotically exact inference on the ATE for general covariate adjusted estimators under finely stratified randomization.
Independently, bai2023adjustment study covariate adjustment under matched pairs randomization in a superpopulation framework. They also find that regression adjustment without pair fixed effects may be inefficient, while adding pair fixed effects restores efficiency. Relative to our work, they additionally study regularized regression adjustment under high-dimensional asymptotics, which we do not consider. By contrast, we study more general forms of stratification, allowing coarse and fine stratification with arbitrary treatment proportions $\propfn \not = 1/2$. For such designs, the strata fixed effects estimator may still be inefficient. To fix this, we introduce novel forms of linear adjustment that are efficient under general stratified designs.
The rest of the paper is organized as follows. In Section (ref) we define notation and introduce the family of stratified designs that we will consider throughout the paper. Section (ref) gives our main results, characterizing optimal covariate adjustment and constructing efficient estimators. Section (ref) provides asymptotically exact inference on the ATE for generic linearly adjusted estimators. In Sections (ref) and (ref), we study the finite sample properties of our method, including both simulations and an empirical application to the experiment in baysan2022. Section (ref) concludes with some recommendations for practitioners.
For a binary treatment $d \in \{0, 1\}$, let $Y_i(1)$, $Y_i(0)$ denote the treated and control potential outcomes, respectively. For treatment assignment $\Di$, let $Y_i = Y_i(\Di) = \Di Y_i(1) + (1-\Di) Y_i(0)$ be the observed outcome. Let $X_i$ denote covariates. Consider data $(X_i, Y_i(1), Y_i(0))_{i=1}^n$ sampled i.i.d. from a superpopulation of interest. We are interested in estimating the average treatment effect in this population, $\ate = E[Y(1) - Y(0)]$. After sampling units $i = 1, \dots, n$, treatments $\Dn$ are assigned by stratified randomization. In particular, we use the “local randomization” framework introduced in cytrynbaum2022local.
Experiment Timing: Suppose that the experimenter does the following
We are agnostic about the exact time at which the covariates are observed, subject to the constraints above. For example, it could be that only $\psi(X)$ is observed at the design stage, while the full vector $X$ is collected later with the outcomes, and the experimenter chooses to adjust for $h(X) \subseteq X$. Alternatively, the full vector $X$ could be observed at the design stage, but the experimenter chooses to only stratify on $\psi(X)$, and adjusts for $h(X) \subseteq X$ at step (3). We may or may not have $\psi(X) \sub h(X)$.\footnote{Our asymptotic framework lets $h(X)$, $\psi(X)$ be fixed as $n \to \infty$.} \\
Consider the unadjusted estimator given by the coefficient $\est$ on $D$ in the regression $Y \sim 1 + D$. Before discussing covariate adjustment, we first state a helpful variance decomposition for $\est$ that will be used extensively below. Let $\catefn(X) = E[Y(1) - Y(0) | X]$ denote the conditional average treatment effect (CATE) and $\hk_d(X) = \var(Y(d) | X)$ the heteroskedasticity function. Define the balance function
We often denote $\balancefn = \balancefn(X; \propfn)$ in what follows. cytrynbaum2022local shows that if $\Dn \sim \localdesigncond(\psi, \propfn)$ then $\rootn(\est - \ate) \convwprocess \normal(0, V)$ with
The variance $V$ is in fact the hahn1998 semiparametric variance bound\footnote{armstrong2022 shows that this variance bound also holds for stratified designs.} for the $\ate$ (with covariates $\psi(X)$), providing a formal sense in which stratification does nonparametric regression adjustment “by design.” The middle term is the most important for our analysis below. For example, in this notation the difference in asymptotic efficiency between stratifications $\psi_1$ and $\psi_2$ (for fixed $\propfn$) is simply $E[\var(\balancefn | \psi_1)] - E[\var(\balancefn | \psi_2)]$. Note also that $E[\var(\balancefn | \psi)] \leq \var(\balancefn)$ for any $\psi$, showing how stratification removes the variance due to fluctuations that are predictable by $\psi(X)$. \\
Moving beyond the difference of means estimator $\est$, suppose that at the analysis stage, the experimenter has access to covariates $h(X)$, which may strictly contain $\psi(X)$. One may try to further improve the efficiency of $\ate$ estimation by regression adjustment using these covariates, either using standard the regression $Y \sim 1 + D + h$ or the regression $Y \sim 1 + D + h + Dh$ (with de-meaned covariates) studied in lin2013. We study the interaction between covariate adjustment and stratification in Section (ref) below, characterizing the optimal linear adjustment.
In this section, we begin by studying the efficiency of commonly used covariate-adjusted estimators of the $\ate$ under stratified randomization. lin2013 showed that in a completely randomized experiment, equivalent to $\Dn \sim \localdesigncond(\psi, \propfn)$ with $\psi=1$, regression adjustment with full treatment-covariate interactions is asymptotically weakly more efficient than difference of means estimation. negi2021 extended this result to ATE estimation in the superpopulation framework that we use in this paper. Interestingly, we show that this result is atypical. For a general stratified experiment with $\psi \not = 1$, lin2013 style regression adjustment may be strictly inefficient relative to difference of means. The problem is that the interacted regression solves the wrong optimization problem, minimizing a marginal variance objective when, due to the stratification, it should instead minimize a mean-conditional variance objective, conditional on the stratification variables $\psi$. In fact, the Lin estimator is totally insensitive to the stratification, estimating the same adjustment coefficient for any stratified design $\Dn \sim \localdesigncond(\psi, \propfn)$. Because of this, interacted regression is generically sub-optimal and in some cases can even be strictly less efficient than difference of means. Before proceeding, we state our main assumption.
Now we are ready to define the Lin estimator and state our first result. Denote $\hi = h(X_i)$ and de-meaned covariates $\hitilde = \hi - \en[\hi]$, with $\en[\hi] \equiv n \inv \sum_{i=1}^n \hi$. The Lin estimator $\estlin$ is the coefficient on $\Di$ in the interacted regression
Define the within treatment arm covariate means $\bar h_1 = \en[\hi \Di] / \en[\Di]$ and $\bar h_0 = \en[\hi (1-\Di)] / \en[1-\Di]$. The Lin estimator $\estlin$ can be related to the difference of means estimator $\est$ as
Here, the adjustment coefficient $\coefflin$ is $\coefflin = (1-\propfn)(\wh a_1 + \wh a_0) + \propfn \wh a_0$, where $\wh a_0$ and $\wh a_1$ are the coefficients on $\hitilde$ and $\Di \hitilde$ in Equation (ref). The following theorem characterizes the asymptotic properties of this estimator under stratified designs.
The variance $V$ differs from the variance of the unadjusted estimator only in the middle term, which changes from $E[\var(\balancefn | \psi)]$ in the unadjusted case to $E[\var(\balancefn - \coefflinpop'h | \psi)]$ for the interacted regression. Crucially, the second statement of Theorem (ref) shows that the adjustment coefficient $\coefflinpop$ attempts to minimize a marginal variance, instead of the mean-conditional variance that shows up in $V$ above. Because of this, the estimator may be inefficient for general stratifications $\psi \not = 1$, since in general \[ \coefflinpop = \argmin_{\gamma \in \mr^{\dimh}} \var(\balancefn - \gamma'h) \not = \argmin_{\gamma \in \mr^{\dimh}} E[\var(\balancefn - \gamma'h | \psi)] \equiv \gamma^*. \]
Observe that the Lin estimator is completely insensitive to the experimental design, estimating the same adjustment coefficient $\coefflinpop = \argmin_{\gamma} \var(\balancefn - \gamma'h)$ for any stratification variables $\psi(X)$. The following example shows that this can lead to strict inefficiency relative to difference of means estimation.
An important special case occurs when the design is completely randomized ($\psi = 1$) or if the covariates and stratification variables are independent $h(X) \indep \psi(X)$. In this case, the Lin estimator is weakly more efficient than difference of means since we have \[ E[\var(\balancefn - \coefflinpop'h | \psi)] = \var(\balancefn - \coefflinpop'h) = \min_{\gamma} \var(\balancefn - \gamma'h) \leq \var(\balancefn). \]
An analogue of Theorem (ref) also holds for the non-interacted regression estimator $Y_i \sim 1 + \Di + \hi$ under stratified designs $\Dn \sim \localdesigncond(\psi, \propfn)$. The non-interacted estimator is known to be inefficient relative to difference of means even for completely randomized experiments unless $\propfn = 1/2$ or treatment effects are homogeneous. For completeness, we give asymptotic theory and optimality conditions for this estimator under stratified randomization in Section (ref) in the appendix.
We noted above that the Lin estimator $\estlin$ can be written in the canonical form $\estlin = \est - \coefflin'(\hbarone - \hbarzero)$. In fact, most commonly used adjusted estimators can be written in the standard form $\est_{adj} = \est - \wh \gamma'(\hbarone - \hbarzero)$ for some $\wh \gamma$, up to order $\Op(n \inv)$ factors. The following theorem describes the asymptotic properties of general covariate-adjusted estimators $\est_{adj}$ of this form. To avoid carrying around factors of $\propfn$ in our variance expressions, in what follows we scale adjusted estimators by the normalization constant $\propconstant = \sqrt{\propfn(1-\propfn)}$.
We define a linearly-adjusted estimator to be asymptotically efficient if it globally minimizes the asymptotic variance $V(\gamma)$ in the previous theorem.
Note that efficiency is defined relative to a design $\Dn \sim \localdesigncond(\psi, \propfn)$ and covariates $h(X)$. Setting $\gamma = 0$ recovers unadjusted estimation, so any optimal estimator is in particular weakly more efficient than difference of means.
Optimal Adjustment Coefficient. If $E[\var(h | \psi)] \succ 0$, then the optimization in Definition (ref) is solved uniquely by a mean-conditional OLS coefficient
Intuitively, fine stratification makes treatment-control imbalances in the covariates $h(X)$ and the potential outcomes $Y(d)$ that are predictable by $\psi$ small enough that they do not contribute to first-order asymptotic variance. Because of this, the optimal covariate-adjusted estimator $\est - \gamma^*(\hbarone - \hbarzero)\propconstant$ ignores such fluctuations, minimizing the mean-conditional variance objective $E[\var(\balancefn - \gamma'h | \psi)]$, instead of the marginal variance $\var(\balancefn - \gamma'h)$ targeted by the Lin estimator.
Optimal Covariates. Intuitively, the form of the variance in Equation (ref) suggests adjusting for variables $h$ that contain predictive information not already contained in $\psi$. The (unknown) optimal covariates are $h^* = \balancefn$. In this case, $\gamma^* = 1$ makes the middle variance term identically zero, and $\estadj$ achieves the armstrong2022 semiparametric variance bound.
Sample Average Treatment Effect. Theorem (ref) may be extended to covariate-adjusted estimation of the sample average treatment effect $\sate = \en[Y_i(1) - Y_i(0)]$. Defining the conditional treatment effect variance $\hkte(X) = \var(Y(1) - Y(0) | X)$, one can show that $\rootn(\est - \sate) \convwprocess \normal(0, V_S(\gamma))$ with
In particular, the optimal adjustment for estimating the ATE and the SATE are the same, with $\gamma^*_{SATE} = \gamma^*$.
Before continuing, we briefly mention some extensions to the framework above that are studied in detail in Appendices (ref)-(ref).
Experiments with Noncompliance. In settings with noncompliance, we may instead consider estimation and inference on the local average treatment effect (LATE) of imbens1994. As a simple application of our main results, Section (ref) characterizes the optimal linear adjustment for estimating the LATE, constructs feasible efficient estimators, and provides asymptotically exact inference on the LATE under stratified randomization with ex-post covariate adjustment.
Varying Treatment Proportions. cytrynbaum2023 extends Definition (ref) to allow fine stratification with non-constant assignment propensity $\propfn(\psi)$. Section (ref) in the appendix characterizes the optimal adjustment coefficient for such designs and derives a feasible efficient estimator.
Nonlinear Adjustment. In some settings, it may be more natural to use nonlinear or nonparametric covariate adjustment to improve efficiency, for example in experiments with binary outcomes. Section (ref) in the appendix characterizes the optimal adjustment over a general function space $\mc H$ for finely stratified designs with varying propensity $\propfn(\psi)$. Feasible estimation of the optimal nonlinear adjustment is an interesting problem that we leave for future work.
This section shows that optimal linear adjustment of a stratified design is as efficient as semiparametric partially linear regression adjustment in an experiment with iid treatments, with adjustment function that is linear in $h(X)$ and nonparametric in $\psi(X)$. This suggests that experimenters stratify on a small set of covariates expected to be most predictive of outcomes at design-time, and (efficiently) adjust for the remaining covariates ex-post. See below for a more detailed discussion of stratification vs.\ adjustment.
The main result of this section shows first-order asymptotic equivalence of the following (design, estimator) pairs \[ (\Dn \sim \localdesigncond(\psi, \propfn), \text{optimal linear}) \iff (\Di \simiid \bern(\propfn), \text{optimal semiparametric}). \]
To define the latter, consider the within-arm partially linear regression models
for $d \in \{0, 1\}$. Define the partially linear adjustment function $F_d(x) = g_d^*(\psi(x)) + h(x)'\gamma_d^*$ and consider a robins95 style augmented inverse propensity weighting (AIPW) estimator
The next theorem shows that optimal linear adjustment of the design $\Dn \sim \localdesigncond(\psi, \propfn)$ is asymptotically equivalent to optimal semiparametric adjustment with nonparametric $\psi(X)$ and linear $h(X)$ components.
The limiting variance $V^*$ is the same as the optimal linearly adjusted variance $V(\gamma^*)$ from Definition (ref). Intuitively, stratification contributes the nonlinear component of the optimal model $F_d(x)$ above, while optimal adjustment contributes the linear component. The optimal adjustment coefficient $\gamma^* = \sqrt{\frac{1-\propfn}{\propfn}} \gamma_1 + \sqrt{\frac{\propfn}{1-\propfn}} \gamma_0$, for partially linear coefficients $\gamma_1, \gamma_0$ defined in Equation (ref) above.
Stratification vs.\ Regression Adjustment. Theorem (ref) shows that stratification provides nonparametric control over the fluctuations of the outcomes predictable by $\psi(X)$, while (linear) adjustment only provides linear control. In first-order asympotics, this suggests that we stratify on all available covariates, since the variance $V^*$ above is minimized by setting $\psi(X) = X$. However, this may perform poorly in finite samples due to a curse of dimensionality for stratification as $\dim(\psi)$ increases. For example, cytrynbaum2023 shows the variance convergence rate $n \var(\est) = V + \Op(n^{-2/(\dim(\psi)+1)})$ for the variance $V$ in Equation (ref), which may be slow even for moderate $\dim(\psi)$. Intuitively, this suggests stratifying on a small set\footnote{It is difficult to give concrete guidance for choosing $\dim(\psi)$, since the relevant quantities such as $E[\var(\balancefn | \psi)]$ are not estimable at design-time, before we have outcome data. The rate above suggests $\dim(\psi) = o(\log n)$ to achieve the variance $V$ in Equation (ref).} of covariates $\psi(X)$ expected to be most predictive of outcomes at design time, and planning to optimally adjust for less predictive covariates $h(X)$ ex-post.
The next two sections show how to construct linearly-adjusted estimators for the design $\Dn \sim \localdesigncond(\psi, \propfn)$ that achieve the optimal variance $V^*$.
This section provides a “rich covariates” style condition on the relationship between adjustment covariates and stratification variables under which a simple parametric correction of the Lin estimator is fully efficient. The basic idea is to include rich functions $z(\psi)$ of the stratification variables in the adjustment set alongside the additional covariates we would like to adjust for ex-post. The main result of this section shows that including $z(\psi)$ as covariates forces the Lin estimator to solve the mean-conditional variance minimization problem of Definition (ref), restoring asymptotic optimality. An analogous result holds for the non-interacted regression estimator $Y \sim 1 + D + h$ if $\propfn = 1/2$. As a simple application of this section's results, Example (ref) shows that for coarsely stratified designs the Lin estimator with leave-one-out strata indicators is efficient.
Consider adjusting for covariates $h(X) = (w(X), z(\psi))$. The main assumption of this section requires that the conditional mean $E[w | \psi]$ is well-approximated by known transformations $z(\psi)$ of the stratification variables.
Our next theorem shows that adding such transformations $z(\psi)$ to the adjustment set recovers full efficiency for the Lin estimator.
In practice, Theorem (ref) suggests including flexible functions $z(\psi)$ of the stratification variables in the adjustment set. The proof is given in Section (ref) of the supplement. The following corollary follows shows that if $\propfn = 1/2$ (matched pairs) or if treatment effect heterogeneity is limited then the non-interacted regression $Y \sim 1 + D + w + z$ with rich strata controls $z(\psi)$ is also asymptotically efficient.
The condition $E[\cov(Y(1)-Y(0), w | \psi)] = 0$ limits the explanatory power of covariates $w$ for treatment effect heterogeneity, conditional on the stratification variables.
The next example uses Theorem (ref) to show that including leave-one-out strata indicators as covariates in the Lin estimator restores asymptotic efficiency for coarsely stratified designs.
Leaving behind the rich covariates Assumption (ref), the next section provides new adjusted estimators that are fully efficient for any design in the class $\Dn \sim \localdesigncond(\psi, \propfn)$ under weak conditions.
In this section, we study several adjusted estimators that are asymptotically efficient under weak conditions for any stratified design $\Dn \sim \localdesigncond(\psi, \propfn)$. For matched pairs designs, or in settings with limited treatment effect heterogeneity, the non-interacted regression including treatment, covariates, and pair fixed effects is efficient. More generally, we show that the following estimators are efficient under weak assumptions.
The main new condition we impose in this section is that the adjustment covariates are not collinear, conditionally on the stratification variables. This guarantees that the optimal adjustment coefficient $\gamma^*$ is unique with $\gamma^* = E[\var(h | \psi)] \inv E[\cov(h, \balancefn | \psi)]$, as discussed in Section (ref).
Note that this assumption rules out adjustment for functions $h(\psi)$ of the stratification variables. To see why it is necessary, consider that, for example, in a regression with full strata fixed effects $Y \sim D + h + z^n$, covariates $\hi = h(\psii)$ would be asymptotically collinear with the strata fixed effects $z^n = (\one(i \in \group_j))_{j=1}^{n/k}$. More intuitively, the problem is that $h(\psi)$ has too little residual variation within local regions of $\psi(X)$ space defining the fine strata. We noted earlier that $\est - \alpha'(\bar h_1 - \bar h_0) = \est - \op(\negrootn)$ for any $\alpha \in \mr^{d_h}$, so such adjustment cannot improve first-order efficiency. Nevertheless, one may still wish to adjust for $h(\psi)$ to correct finite sample imbalances not controlled by the design. Adjustment for such variables needs to be handled slightly differently, and we construct modified efficient estimators for this purpose in Section (ref) below.
Recall that for $\propfn = a/k$, a finely stratified design $\Dn \sim \localdesigncond(\psi, \propfn)$ partitions the experimental units $\{1, \dots, n\}$ into $n/k$ disjoint groups $\group$. Define the fixed effects estimator $\estnaivefe$ by the least squares equation
The next theorem shows that $\estnaivefe$ is fully efficient in the case of matched pairs or if treatment effect heterogeneity is limited, but may be inefficient in general.
See Section (ref) for the proof. Asymptotically exact inference for the $\ate$ using $\estnaivefe$ is available using the tools in Section (ref).
In the rest of this section, we develop estimators that are fully efficient for any finely stratified design, without imposing any assumptions on treatment effect heterogeneity or treatment proportions.
First, we define a partialled version of the Lin estimator. Let $\group(i)$ denote the group that unit $i$ belongs to and define the within-group partialled covariates \[ \hicheck = \hi - \frac{1}{k} \sum_{j \in \group(i)} \hj. \] For example, if $k = 2$ this is just the within-pair covariate difference $\hipartial = (1/2)(\hi - h_{m(i)})$, where $i$ is matched to $m(i)$. We can think of $\hicheck$ as an inconsistent but approximately unbiased signal for the non-parametrically residualized covariate $\hi - E[\hi | \psii]$. Next, we use these partialled covariates in the Lin regression
Define the partialled Lin estimator $\estlinpartial$ to be the coefficient on $\Di$ in this regression. For reference, similarly to the Lin regression we may write this in the standard form $\estlinpartial = \est - \coefflinpartial'(\hbarone - \hbarzero)\propconstant$ with adjustment coefficient $\coefflinpartial = (\wh a_1 + \wh a_0) \sqrt{\frac{1-\propfn}{\propfn}} + \wh a_0 \sqrt{\frac{\propfn}{1-\propfn}}$, where $\wh a_0$ and $\wh a_1$ are coefficients on $\hicheck$ and $\Di \hicheck$. \\
Our main result in Theorem (ref) below shows that the partialled Lin estimator $\estlinpartial$ is asymptotically efficient in the sense of Definition (ref), with $\coefflinpartial \convp \gamma^*$ for the optimal adjustment coefficient $\gamma^*$.
Next, we generalize an estimator proposed by imbens2015 for covariate adjustment in matched pairs experiments to more general stratified designs. For each group of units $\group = 1, \dots, n/k$ in the design $\Dn \sim \localdesigncond(\psi, \propfn)$, define the within-group difference of means of outcomes and covariates \[ \yg = \frac{1}{k} \sum_{i \in \group} \frac{Y_i\Di}{\propfn} - \frac{1}{k} \sum_{i \in \group} \frac{Y_i(1-\Di)}{1-\propfn} \quad \text{and} \quad \hg = \frac{1}{k} \sum_{i \in \group} \frac{\hi \Di}{\propfn} - \frac{1}{k} \sum_{i \in \group} \frac{\hi (1-\Di)}{1-\propfn}. \] For any group-indexed array $(x_g)_g$, denote $\eg[x_g] = \frac{k}{n} \sum_{\group} x_g$. Define the Group OLS estimator $\estgroupols$ by the regression
with $\eg[(1, \hg) e_g] = 0$. For motivation, note that if $h = 0$ then this becomes $\yg = \estgroupols + e_g$ and $\estgroupols$ is just the unadjusted estimator $\estgroupols = \bar Y_1 - \bar Y_0$. More generally, the adjusted version can be written $\estgroupols = \eg[y_g] - \coeffgroupols'\eg[\hg] = \est - \coeffgroupols'(\bar h_1 - \bar h_0)$ with adjustment coefficient $\coeffgroupols = \var_g(h_g)\inv \cov_g(h_g, y_g)$. The estimators $\estgroupols$ and $\estlinpartial$ are numerically identical for the case of matched pairs, but not for $\propfn \not = 1/2$. The main result of this section shows that $\estgroupols$ is asymptotically equivalent to the partialled Lin estimator $\estlinpartial$, and both are asymptotically optimal.
Finally, we define tyranny-of-the-minority (ToM) adjustment, extending lin2013. To do so, define the adjustment coefficient
Define the ToM estimator $\esttom = \est - \coefftom'(\hbarone - \hbarzero)\propconstant$. The main difference between the ToM and Partialled Lin adjustment coefficients is that $\coefftom$ estimates the conditional variance $E[\var(h |\psi)]$ only once, using the sample variance $\var_n(\hicheck)$ for the full experimental sample. By contrast, partialled Lin estimates this term separately in each treatment arm, using $\var_n(\hicheck | \Di=1)$ and $\var_n(\hicheck | \Di=0)$. Because of this, we expect $\esttom$ to be more stable than $\estlinpartial$ in small experiments.
The main result of this section shows that all three estimators above are asymptotically equivalent and efficient in the sense of Definition (ref).
Methods for asymptotically exact inference on the $\ate$ using these estimators are discussed in Section (ref) below. Our simulations and empirical results show that the partialled Lin, Group OLS, and ToM estimators behave very similarly in finite samples.
In this section, we provide modified versions of the previous estimators that allow further adjustment for covariates $z(\psi)$ that are functions of the stratification variables. As discussed above, this cannot improve first-order efficiency but may improve finite sample performance by correcting for any remaining imbalances in $\psi$ not controlled by the stratification.
Denote $\zi = z(\psii)$. For each estimator $\est_k$ above with $k \in \{FE, PL, G, TM\}$, we define a modified estimator of the form $\wh \tau_k = \est_k - \wh \alpha_k'(\bar z_1 - \bar z_0) \propconstant$. For the fixed effects estimator, define $\estnaivefez$ to be the coefficient on $\Di$ in the regression $Y_i \sim (1, \Di, \hicheck, \zi)$. For the partialled Lin estimator, define $\estlinpartialz$ to be the coefficient on $\Di$ in the regression $Y_i \sim (1, \hicheck, \zi) + \Di(1, \hicheck, \zi)$. Define the modified ToM estimator to be as in Equation (ref), with $(\hicheck, \zi)$ in place of $\hicheck$. Finally, define the modified group OLS estimator $\estgroupolsz = \estgroupols - \coeffgroupolsadjust'(\bar z_1 - \bar z_0) \propconstant$, with $\coeffgroupolsadjust = \coefflinpartialadjust$. Our next theorem shows that these estimators are asymptotically equivalent to the original versions of each estimator that do not adjust for $z(\psi)$. However, the simulations in Sections (ref) and (ref) show that they may perform better in small experiments.
From the second statement of the theorem, we can interpret the modified estimators as taking a conservative approach that ignores stratification on $\psi$ and adjusts for imbalances in $z(\psi)$ as if the experiment were completely randomized.
In this section, we provide asymptotically exact confidence intervals for the ATE in stratified experiments using generic linearly adjusted estimators. Overcoverage is known to be a problem for inference based on the usual Eicker-Huber-White (EHW) variance estimator in stratified experiments. For example, bai2021inference shows that the EHW variance estimators for $Y \sim 1 + D + h$ and the fixed effects regression $Y \sim D + h + z^n$ are asymptotically conservative for matched pairs designs if $h = 0$. To the best of our knowledge, we give the first asymptotically exact inference methods for covariate-adjusted ($h \not = 0$) $\ate$ estimation under general stratified designs. Our main inference result applies to any estimator of the form $\est - \wh \gamma'(\hbarone - \hbarzero)\propconstant + \op(\negrootn)$. In particular, this enables asymptotically exact inference on the $\ate$ using any of the estimators in this paper. Our confidence intervals are shorter than those produced by EHW in the simulations and empirical application below, taking full advantage of the efficiency gains from both stratification and covariate adjustment.
To define our inference methods, consider such an estimator $\est(\wh \gamma) = \est - \wh \gamma'(\hbarone - \hbarzero)\propconstant$ with $\wh \gamma \convp \gamma$. Define the augmented potential outcomes $Y_i^a(d) = Y_i(d) - \propconstant \wh \gamma' \hi$ for $d \in \{0, 1\}$ and the augmented outcome $Y^a_i = Y_i - \propconstant \wh \gamma' \hi$. Then apparently
Our strategy is to apply the inference results of cytrynbaum2023 for difference of means estimation $\est = \bar Y_1 - \bar Y_0$ to the difference of augmented potential outcomes $\bar Y^a_1 - \bar Y^a_0$. To do so, let $\groupset_n$ denote the set of groups in Definition (ref). For each $\group \in \groupset_n$ define the group centroid $\bar \psi_{\group} = |\group|\inv \sum_{i \in \group} \psii$. Let $\groupmatching: \groupset_n \to \groupset_n$ be a bijective matching between groups satisfying $\groupmatching(\group) \not = \group$, $\groupmatching^2 = \identity$, and the homogeneity condition
In practice, $\groupmatching$ is obtained by simply matching the group centroids $\bar \psi_{\group}$ into pairs using the derigs1988 matching algorithm. Let $\groupsetnu_n = \{\group \cup \groupmatching(\group): \group \in \groupset_n\}$ be the unions of paired groups formed by this matching. Define $a(\group) = \sum_{i \in \group} \Di$ and $k(\group) = |\group|$. Finally, define the variance estimator components
Next, define the variance estimator
Our inference strategy begins with the sample variance of adjusted estimator, which is consistent for the asymptotic variance of $\estadj$ under an iid design, but too large under stratified designs. We correct this sample variance using the estimators above, which measure how well the stratification variables predict augmented outcomes in local regions of the covariate space. This section's main result shows that $\varest$ is consistent for the limiting variance of Theorem (ref), enabling asymptotically exact inference on the $\ate$ using adjusted estimators.
By Theorem (ref) and our previous asymptotic results in Theorem (ref), the confidence interval $\wh C = [\est(\wh \gamma) \pm \varest^{1/2} c_{1-\alpha/2} / \rootn]$ with $c_{\alpha} = \Phi \inv(\alpha)$ is asymptotically exact in the sense that $P (\ate \in \wh C ) = 1-\alpha + o(1)$.
In this section, we use simulations to test the finite sample performance of the estimators studied above. We consider quadratic outcome models of the form \[ Y_i(d) = \psii'Q_d\psii + \psii'L_d + c_d \cdot u(X_i) + \residual_i^d \quad \quad E[\residual_i^d | X_i] = 0 \] for $d \in \{0,1\}$. The component $u_i = u(X_i)$ represents covariate signal that is independent of the stratification variables $\psi(X_i)$. After implementing the design $\Dn \sim \localdesigncond(\psi, \propfn)$, we receive access to scalar covariates $\hi$ that are correlated with both $\psii$ and $Y_i(d)$. In particular, suppose that $\hi = \psii'Q_h \psii + \psii'L_h + u_i$ with $E[u_i | \psii] = 0$. In the following simulations, we let $\psii \sim N(0, I_{m})$, $u_i \sim N(0, 1)$, and $\epsilon_i^d \sim N(0, 1/10)$ with $(\psii, \ui, \epsilon_i^d)$ jointly independent. We use treatment proportions $\propfn = 2/3$ unless otherwise specified. With $m \equiv \dim(\psi)$, let $A \in \mr^{m \times m}$ have $A_{ij} = 1$ for $i\not=j$ and $A_{ii} = 0$. We simulate the following DGP's:
We begin by comparing the efficiency properties of different linearly adjusted estimators. Unadj refers to simple difference of means (unadjusted). The Lin estimator is studied in Theorem (ref). Naive refers to the non-interacted regression $Y \sim (1, D, h)$, (Theorem (ref)). FE refers to the fixed effects estimator (Theorem (ref)) and Plin the partialled Lin estimator (Theorem (ref)). GO refers to Group OLS and \textbf{ToM} refers to Tyranny-of-the-Minority estimation (Theorem (ref)). \textbf{Strata Controls} refer to modified versions of each of the previous estimators that further adjust for parametric strata controls $z(\psi)$, as discussed in Section (ref). In our simulations, we set $z(\psi) = \psi$. \textbf{Ad} refers to an adaptive\footnote{This estimator is pointwise asymptotically equivalent to $\estlinpartial$. Issues with post model-selection inference (e.g.\ Leeb2005) appear to be less worrying here, since even under the fixed alternative $\gamma^* \not = \coefflinpop$, the Lin estimator is still $\rootn$-consistent and asymptotically unbiased.} estimator that sets $\estadj = \estlin$ if $\wh V(\coefflin) \leq \wh V(\coefflinpartial)$ and $\estadj = \estlinpartial$ otherwise,\footnote{ We could also use a cross-fit version of $\wh V(\gamma)$ to reduce bias. However, the in-sample criterion performed quite well in our simulations.} including parametric controls $z(\psi)=\psi$ in both cases.\\
Table (ref) studies finite sample efficiency. We present the mean squared error (MSE) ratio, relative to unadjusted estimation, for each of the adjusted estimators above. The bottom line of the table reports the excess risk $R_k$ of each estimator $k$ relative to the optimal estimator. To define this, let $\mse_{k, s}$ be the relative MSE of estimator $k$ in simulation $s$. Then we set $R_k = (1/S)\sum_s (\mse_{k, s} - \min_j \mse_{j, s})$, averaging over all simulations in the table. All results are calculated using $2000$ Monte Carlo repetitions.
In models 1, 2, and 3, both Naive and Lin style linear adjustment are strictly inefficient relative to unadjusted estimation. These models have marginal covariance $\cov(Y(d), h) > 0$ but conditional covariance $E[\cov(Y(d), h | \psi)] < 0$, conditional on the stratification variables. Because of this, the optimal adjustment coefficient $\gamma^* < 0$, while the Naive and Lin regressions estimate positive adjustment coefficients $\gamma_N, \gamma_L > 0$, leading to even worse performance than unadjusted estimation in some cases. For Models 4 and 5, the Naive and Lin methods are competitive with the generic efficient methods from Section (ref). This is because in these cases we made it so that $\cov(Y(d), h) \approx E[\cov(Y(d), h | \psi)]$, so that “by chance” $\gamma^*$ is close to $\gamma_N$ and $\gamma_L$. However, the parametric coefficients $\gamma_N$ and $\gamma_L$ are estimated more precisely than the semiparametric object $\gamma^* = E[\var(h | \psi)] \inv E[\cov(h, \balancefn | \psi)]$. For Model 6, Lin with $z(\psi) = \psi$ controls is (approximately) optimal by Theorem (ref), since $E[w | \psi]$ is (approximately) linear in $\psi$.
Summarizing our findings, the Lin, Plin, and Naive estimators with parametric Strata Controls $z(\psi)=\psi$ had low excess risk across specifications, while the Ad estimator was the most efficient overall. The Naive and \textbf{Lin} estimators without strata controls or within-stratum partialling had large MSE. The \textbf{Plin}, \textbf{GO}, and \textbf{ToM} estimators had similar MSE across model specifications. These generic methods performed the best in regimes with large $n$, small $\dim(\psi)$, and nonlinear $E[h | \psi]$. In these cases, the gap $\gamma_L - \gamma^*$ between the sub-optimal Lin coefficient and optimal coefficient $\gamma^*$ dominates the additional variability $\var(\wh \gamma^*) > \var( \wh \gamma_L)$ required to estimate $\gamma^*$ (this variability increases with $\dim(\psi)$). For example, \textbf{Plin} with $z(\psi)$ controls performs the best when $(n, \dim(\psi)) = (1200, 2)$, but \textbf{Lin} with $z(\psi)$ controls is much better when $\dim(\psi) = 5$. The \textbf{Ad} estimator used a variance pre-test to choose between \textbf{Plin} and \textbf{Lin} (including $z(\psi)$ controls), allowing it to perform well in both regimes.
Table (ref) reports finite sample efficiency and coverage properties of the asymptotically exact inference methods developed in Section (ref). We let $n=1200$ and $\dim(\psi) = 5$. The first panel shows $\%$ change in confidence interval length relative to unadjusted estimation. All confidence intervals are computed using the method in Theorem (ref). We see that the relative efficiency of different estimators are reflected by our inference methods. In particular, asymptotically exact inference allows researchers to report shorter confidence intervals when a more efficient adjustment method is used. In the second panel, we show coverage probabilities for our asymptotically exact confidence interval across a range of linearly adjusted estimators. The final panel shows coverage probabilities for confidence intervals based on the usual HC2 variance estimator, where applicable. The HC2-based confidence intervals significantly overcover.
In this section we apply our methods to the experiment in baysan2022,\footnote{The data is available from baysan2022data.} who estimates the effect of a political information campaign on support for a 2017 Turkish referendum removing checks and balances on executive power. The campaign was administered by the opposition Republican People's Party (CHP), who opposed the referendum. Randomization was performed at the neighborhood level, stratified on quartiles of CHP vote share in the previous 2015 elections. The main outcome is the “No” vote share in the 2017 referendum.\footnote{baysan2022 estimates effects of the campaign on vote share at both the ballot box and neighborhood level. We focus on the neighborhood level effects.} Due to the cost of administering the campaign, $\propfn = 2/11$ out of $n=550$ total neighborhoods were treated. In the original analysis, baysan2022 performed non-interacted covariate adjustment (Theorem (ref)) for $h(X) = $ number of registered voters, number of valid votes, number of votes for the CHP in 2015, CHP vote share in 2015, voter turnout, and CHP vote share quartile fixed effects.
In the first block of Table (ref), we replicate the neighborhood-level analysis of baysan2022. $\estadj$ is the point estimate from each adjustment strategy, SE is the asymptotically exact standard error from Section (ref), and EHW is the usual robust standard error (HC2). Estimates in the “strata controls $z(\psi)$” section include quartile fixed effects, while the leftmost section does not. The results in Section (ref) show that Lin adjustment with quartile fixed effects is efficient in this case, and indeed this has the smallest estimated standard error. The generic efficient estimators have slightly larger SE. The asymptotically exact standard errors from Section (ref) are generally similar to or smaller than EHW, except for the Lin, FE, and Plin estimators with $z(\psi)$ controls. However, our simulation also showed that EHW standard errors may severely undercover in these cases.\footnote{We also note that bai2023tuples have found the EHW standard error from a linear regression with block fixed effects to be potentially invalid in a related problem.}
Overall, changing the adjustment method did not have an economically meaningful effect on the conclusions of the study, and we recover the null effect of baysan2022 in all cases. The covariate $h_k = $ “CHP vote share in 2015” is highly predictive of $Y = $ “CHP vote share in 2017,” so adjusting for this variable ex-post provides a modest variance reduction even after stratifying on 2015 vote share quartiles. However, the estimated optimal coefficient $\gamma^*_k \approx 0.27$ and Lin coefficient $\gamma_{L,k} \approx 0.31$ are quite similar, so (inefficient) Lin adjustment still performs quite well. The other covariates such as $h_j = $ “voter turnout” are very weak predictors of outcomes, so changing the adjustment coefficient on these variables doesn't matter much.
Next, we ask how each estimator would have performed in the experiment in baysan2022 under counterfactual randomization procedures, such as fine stratification.\footnote{Algorithms and inference methods for fine stratification with $\propfn \not =1/2$ have only been developed recently, e.g.\ bai2020pairs and cytrynbaum2023.} To do so, we follow the nonparametric imputation strategy in bai2020pairs, defining potential outcomes $\wh Y_i(d) = Y_i$ if $\Di=d$ and matching imputation $\wh Y_i(d) = Y_{j(i)}(d)$ with $j(i) = \argmin_{j: \Dj=d} |X_i - X_j|_2$ if $\Di \not = d$. We let the matching variables $X_i$ include all controls used in the analysis of baysan2022. Given the imputed data $(X_i, \wh Y_i(0), \wh Y_i(1))_{i=1}^n$, we do the following simulation exercise: (1) draw treatment assignments $\Dn \sim \localdesigncond(\psi, \propfn)$, (2) reveal outcomes $\wh Y_i = \wh Y_i(\Di)$ and (3) form each estimator $\estadj$. We report average point estimates and standard errors over $N=2000$ Monte Carlo repetitions of this procedure.
The first block of Table (ref) uses this imputation procedure to reproduce the empirical results in Table (ref), stratifying by quartiles of CHP vote share and adjusting for exactly the same covariates. The standard errors are very similar to those in the empirical analysis, which provides some validation for this imputation exercise. In the second block of Table (ref), we simulate a design with fine stratification on 2015 CHP vote share, rather than just stratifying by quartiles of the vote share as in baysan2022. We used a matched $11$-tuples design, letting $\Dn \sim \localdesigncond(\psi, \propfn)$ for $\propfn=2/11$ and $\psi = (\text{2015 CHP vote share})$. Covariates $h(X)$ are as above, with $z(\psi) = \psi$. In the third block, we simulate a matched pairs design $\Dn \sim \localdesigncond(\psi, 1/2)$. Note that $\propfn=1/2$ was infeasible in the original experiment due to the high cost of treatment. The last block uses the design $\Dn \sim \localdesigncond(\psi_{alt}, \propfn)$ for $\psi_{alt} = (\text{CHP vote share}, \text{Num.\ of registered voters}, \text{Num.\ of valid votes})$, $\propfn=2/11$, and covariates $h = \text{Turnout}$.
We make some brief observations about this simulation exercise. First, note that the Naive and Lin adjustment are strictly less efficient than unadjusted estimation under simulated fine stratification, consistent with Theorems (ref) and Section (ref). Lin and partialled Lin with $z(\psi)$ controls are the most efficient. Adjustment for extra covariates $h$ doesn't significantly improve efficiency relative to the baseline efficiency gain from finely stratifying on $\psi$ and adjusting for $z(\psi)$ ex-post. Using a matched pairs design $\propfn = 1/2$ improves efficiency, though the improvement is small considering that this design would require providing the information campaign to $175$ extra neighborhoods. Finally, fine stratification on $\psi_{alt}$ significantly reduces efficiency. This is because the extra covariates are not very predictive of outcomes, but stratifying on these covariates force us to use worse matches on the important covariate $\psi = (\text{2015 CHP vote share})$.
Stratified randomization and covariate adjustment are both commonly used in the design and analysis of experiments. In general, experimenters should stratify on a few variables $\psi(X)$ expected to be most predictive of outcomes at design-time, and plan to adjust for imbalances in the remaining covariates $h(X)$ ex-post, as discussed in Section (ref). Our analysis showed that under stratified randomization, the usual regression adjusted estimators can be inefficient. Motivated by this, we provide feasible alternatives that are asymptotically optimal in the class of linearly adjusted estimators. We conclude by giving some recommendations for empirical practice based on the theory, simulations, and empirical results above.
We recommend that applied researchers use either (1) the Lin estimator with parametric strata controls $z(\psi)$ (e.g.\ $z(\psi) = \psi)$ or (2) the partialled Lin estimator with parametric controls $z(\psi)$, since these estimators performed the best across our simulations and empirical application. Lin with parametric controls $z(\psi)$ is efficient under a rich covariates condition (Section (ref)), while partialled Lin is generically efficient (Section (ref)). Both estimators are robust to treatment effect heterogeneity, while the strata fixed effects estimator (Theorem (ref)) is not unless $\propfn = 1/2$.
In our simulations, partialled Lin had good finite sample performance in regimes where $n$ was large relative to $\dim(\psi)$, especially when $E[h|\psi]$ was very nonlinear. Lin with $z(\psi)=\psi$ controls performed better when $\dim(\psi)$ was large relative to $n$, or if $E[h | \psi]$ was approximately linear. To decide which regime we are in, we suggest model selection using a variance pre-test, choosing Lin if $\wh V(\wh \gamma_L) \leq \wh V(\coefflinpartial)$ and partialled Lin otherwise. This adaptive estimator (Ad in Section 6) was efficient in both regimes and had good coverage properties. We leave a more general study of such post model-selection estimators in this context to future work.
Regardless of the adjustment strategy, we recommend using the asymptotically exact confidence intervals provided in Section (ref). Our simulations showed close to nominal coverage for these confidence intervals across all considered estimators. By contrast, confidence intervals based on the HC2 robust variance estimator often had significant overcoverage.