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.
63,573 characters · 13 sections · 57 citation commands
Inference in Experiments with Matched Pairs and Imperfect Compliance
KEYWORDS: Matched pairs, randomized controlled trial, experiments, noncompliance, imperfect compliance
JEL classification codes: C12, C14
\thispagestyle{empty} \setcounter{page}{1}
This paper studies inference for the local average treatment effect in randomized controlled trials with imperfect compliance, when treatment status is determined according to a “matched pairs" design. By “matched pairs,” we mean that units are sampled i.i.d.\ from the population of interest, paired according to observed, baseline covariates and finally, within each pair, one unit is selected at random for treatment. This method is used routinely in all parts of the sciences. Indeed, commands to facilitate its implementation are included in popular software packages, such as sampsi in Stata. References to a variety of specific examples can be found, for instance, in the following textbook treatments of randomized experiments: glennerster2014running, riach2002field, rosenberger2015randomization. See also bruhn2009pursuit, who, based on a survey of selected development economists, report that 56% of researchers have used such a design at some point. In many such experiments, compliance may be imperfect: some recent examples of experiments featuring both matched pairs and imperfect compliance are groh2016macroinsurance and resnjanskij2021can. Under weak assumptions that ensure pairs are formed so that units within pairs are suitably “close” in terms of observed, baseline covariates, we derive a variety of results pertaining to inference about the local average treatment effect in such experiments.
We first study the behavior of the usual Wald (i.e., two-stage least squares) estimator of the local average treatment effect. When all observed, baseline covariates are used in forming pairs, we find that the estimator is efficient among all estimators for the local average treatment effect in the sense that it achieves the lower bound on the limiting variance developed in bai2023efficient over a broad class of treatment assignment schemes that hold the marginal probability of treatment assignment equal to one half, thereby including matched pairs, in particular, as a special case. On the other hand, we find that the conventional heteroskedasticity-robust estimator obtained from a two-stage least squares regression (with and without pair fixed effects) is conservative in that its limit in probability is always weakly larger than the limiting variance, and strictly larger unless treatment effect heterogeneity is constrained in a particular fashion. As a result, we provide an alternative estimator of the limiting variance and show that it is consistent. In our simulation study, we find that tests using the Wald estimator together with conventional heteroskedasticity-robust estimators may as a consequence have worse power when compared to tests using the Wald estimator together with our estimator of its limiting variance.
Next, we analyze the behavior of covariate adjusted Wald estimators for settings in which there are additional observed baseline covariates that were not used when pairing units. We first derive the limiting behavior of a class of covariate-adjusted estimators indexed by different “working models” for the conditional expectations of the outcome and treatment take-up with respect to the covariates. Importantly, these working models need not be correctly specified for the true conditional expectations in order for this estimator to remain consistent for the local average treatment effect, but we show that the limiting variance of the estimator is minimized, in particular, when they are correctly specified. We then specialize these results to the case where the working models are linear as a function of the covariates, and derive the form of the optimal linear working model. We further show that the resulting optimal linear covariate adjusted estimator can be computed using a two-stage least squares estimator using treatment assignment as an instrument for treatment status in an instrumental variables regression of the outcome on the following quantities: treatment status; (functions of) the observed baseline covariates; and pair fixed effects. We emphasize, however, that the resulting heteroskedasticity robust variance estimator is not guaranteed to be consistent, and so we also provide a suitable consistent estimator of the asymptotic variance. In our simulation study, we find that tests using our optimal linear covariate adjusted estimator together with our consistent estimator of its limiting variance have better power than tests constructed using the unadjusted Wald estimator or tests constructed using sub-optimal linear covariate adjustments.
Our paper builds upon the analysis of bai2022mp, who analyzed the behavior of the difference-in-means estimator of the average treatment effect in context of experiments with matched pairs and perfect compliance. Our covariate-adjusted estimator is inspired by the analysis in bai2023covariate, who studied the use of additional, observed, baseline covariates in experiments with matched pairs and perfect compliance to improve the precision with which we can estimate the average treatment effect cytrynbaum2023covariate. We emphasize, however, that none of the aforementioned papers permit imperfect compliance, which, as argued by athey2017econometrics, is one of the most common complications in even the most well designed experiments. We note, however, that imperfect compliance has been studied in the context of randomized controlled trials with other treatment assignment schemes, such as stratified block randomization: see, for example, ansel2018ols, bugni2021inference, and jiang2022improving. We also emphasize that all of these papers, like ours, carry out their analysis in a superpopulation sampling framework. In this way, our analysis differs from the analysis of experiments with a finite population sampling framework athey2017econometrics, chaisemartin2022at, ding2017bridging, ren2021model.
The remainder of the paper is organized as follows. In Section (ref) we describe our setup and notation. In particular, there we describe the precise sense in which we require that units in each pair are “close” in terms of their baseline covariates. Our main results concerning the Wald estimator are contained in Section (ref). In Section (ref), we develop results pertaining to our covariate-adjusted estimator that exploits additional observed, baseline covariates not used in pairing units. In Section (ref), we examine the practical relevance of our theoretical results via a small simulation study. In Section (ref), we provide a brief empirical illustration of our proposed tests using data from an experiment in groh2016macroinsurance. Finally, we conclude in Section (ref) with some recommendations for empirical practice guided by both our theoretical results and our simulation study. As explained further in that section, we do not recommend the use of the Wald estimator with the conventional heteroskedasticity-robust estimator of its limiting variance because it is conservative in the sense described above; we instead encourage the use of the Wald estimator with our consistent estimator of its limiting variance because it is asymptotically exact, and, as a result, can be considerably more powerful. When there are additional, observed, baseline covariates that are not used when forming pairs, we recommend the use of our covariate-adjusted Wald estimator with our consistent estimator of its limiting variance. Proofs of all results are provided in the Appendix.
Let $Y_{i} \in \mathbf{R}$ denote the (observed) outcome of the $i$th unit, $A_i \in \{0, 1\}$ be an indicator for whether or not unit $i$ is assigned to treatment, $D_i \in \{0, 1 \}$ be an indicator for whether or not unit $i$ decides to take up treatment, $X_i \in \mathbf{R}^{k_x}$ denote observed, baseline covariates for the $i$th unit which are used for matching, and $W_i \in \mathbf{R}^{k_w}$ denote observed, baseline covariates for the $i$th unit which will be used when we consider covariate adjustment in Section (ref). In contrast to the setting considered in bai2022mp, we allow for imperfect compliance, i.e. for $D_i \ne A_i$. Further denote by $Y_i(d)$ the potential outcome of the $i$th unit if they make treatment decision $d \in \{0, 1\}$, and by $D_i(a)$ the potential treatment decision of the $i$th unit if assigned to treatment $a \in \{0, 1\}$. The observed treatment decision and potential treatment decision are related to treatment assignment via the usual relationship
and the observed outcome and potential outcome are related to treatment decision via the relationship
We will also often make use of the following alternative representation for the observed outcome, which is numerically equivalent to ((ref)):
where
for $a \in \{0, 1\}$. In words, $\tilde Y_i(a)$ represents the “intention-to-treat" potential outcome for unit $i$ when assigned to treatment $a \in \{0, 1\}$.
Following angristimbens1994late, each participant in the experiment can be categorized into one of four types: units for which $D_i(1) = 1$ and $D_i(0) = 0$ are referred to as compliers, units for which $D_i(1) = 1$ and $D_i(0) = 1$ are referred to as always takers, units for which $D_i(1) = 0$ and $D_i(0) = 0$ are referred to as never takers, and finally units for which $D_i(1) = 0$ and $D_i(0) = 1$ are referred to as defiers. We use the notation
below to indicate whether or not unit $i$ is a complier.
Throughout the paper we will study inference on samples with $2n$ observations, so that $n$ indexes the number of pairs of observations. For a random variable indexed by $i$, say for example $A_i$, it will be useful to denote by $A^{(n)}$ the random vector $(A_1, ..., A_{2n})$. Denote by $P_n$ the distribution of the observed data $(Y^{(n)}, D^{(n)}, A^{(n)}, X^{(n)}, W^{(n)})$, and by $Q_n$ the distribution of $(Y^{(n)}(1), Y^{(n)}(0), D^{(n)}(1), D^{(n)}(0), X^{(n)}, W^{(n)})$. Note that $P_n$ is jointly determined by ((ref)), ((ref)), $Q_n$, and the mechanism for determining treatment assignment. We assume that our sample consists of $2n$ i.i.d.\ observations i.e. that $Q_n = Q^{2n}$, where $Q$ is the marginal distribution of $(Y_i(1), Y_i(0), D_i(1), D_i(0), X_i, W_i)$. We therefore state our assumptions below in terms of assumptions on $Q$ and the mechanism for determining treatment assignment. Indeed, we will not make reference to $P_n$ in the sequel and all operations are understood to be under $Q$ and the mechanism for determining treatment assignment.
Our object of interest is the local average treatment effect, which may be expressed in our notation as
For a pre-specified choice of $\Delta_0$, the testing problem of interest is
at level $\alpha \in (0, 1)$.
We begin by describing our primary assumptions on the data generating process.
Assumption (ref)(a)--(b) are mild restrictions imposed to rule out degenerate situations and to permit the application of suitable laws of large numbers and central limit theorems. Assumption (ref)(c) is a smoothness requirement that ensures that units that are “close" in terms of their baseline covariates are suitably comparable. Similar smoothness requirements are also considered in bai2022mp, and generally play a key role in establishing the asymptotic exactness of our proposed tests. Assumptions (ref)(d)--(e) are the standard “monotonicity" and “relevance" conditions of angristimbens1994late which ensure that the probability limit of the Wald estimator which we define in Section (ref) can be interpreted as the local average treatment effect $\Delta(Q)$.
Next, we describe our assumptions on the mechanism determining treatment assignment. Following the notation in bai2022mp, the $n$ pairs can be represented by the sets
where $\pi = \pi_n\left(X^{(n)}\right)$ is a permutation of $2n$ elements. Given such a $\pi$, we assume that treatment status is assigned as described in the following assumption:
We further require that the units in each pair be “close" in terms of their baseline covariates in the following sense:
We will also sometimes require that the distances between units in adjacent pairs be “close" in terms of their baseline covariates:
bai2022mp, bai2023inference, and cytrynbaum2023designing provide several examples of pairing algorithms which satisfy Assumptions (ref)--(ref). The simplest such example is when $X_i \in \mathbf{R}$, in which case we can order units from smallest to largest according to $X_i$ and pair adjacent units. It then follows from Theorem 4.1 in bai2022mp that Assumptions (ref)--(ref) are satisfied as long as $E[X_i^2] < \infty$.
In this section, we study the asymptotic properties of the standard Wald estimator (i.e., the two-stage least squares estimator of $Y_i$ on $D_i$ using $A_i$ as an instrument) of $\Delta(Q)$ under a matched pairs design. In order to introduce this estimator, define
Using this notation, the Wald estimator is defined as
Note that this estimator may be obtained as the ratio of the estimator of the coefficient of $A_i$ in an ordinary least squares regression of $Y_i$ on a constant and $A_i$ (the “reduced form") to the estimator of the coefficient of $A_i$ in an ordinary least squares regression of $D_i$ on a constant and $A_i$ (the “first stage"). Theorem (ref) establishes the limiting distribution of $\hat{\Delta}_n$ under a matched pairs design.
We derive Theorem (ref) by reproducing arguments in Lemma S.1.4 in bai2022mp, but replacing the usual potential outcomes $Y_i(a)$ with the transformed outcomes $Y_i^\ast(a)$\footnote{We note that although a version of Theorem (ref) can be obtained by directly appealing to more recent results in bai2023efficient, the regularity conditions introduced in this paper are in fact weaker and, in our view, easier to interpret than the general conditions considered there.}. As a result, the the numerator of $\nu^2$ corresponds exactly to the limiting variance obtained in bai2022mp with the usual potential outcomes $Y_i(a)$ replaced with $Y_i^\ast(a)$. In particular, when there is perfect compliance, so that $D_i = A_i$, $P \{C_i = 1\} = 1$, and $D_i(a) = a$, the limiting variance we obtain in Theorem (ref) corresponds exactly to the limiting variance derived in bai2022mp. It can be shown that our expression for $\nu^2$ attains the efficiency bound derived in bai2023efficient over a broad class of treatment assignments which include matched pairs as a special case frolich2007nonparametric, hong2010semiparametric
In this section, we construct a consistent variance estimator for the limiting variance $\nu^2$, and then contrast this to the asymptotic behavior of standard regression-based variance estimators. As noted in the discussion following Theorem (ref), the expression for $\nu^2$ corresponds exactly to the limiting variance obtained in bai2022mp with the usual potential outcomes $Y_i(a)$ replaced with the transformed outcomes $Y_i^\ast(a)$. We thus follow the variance construction from bai2022mp, but with a feasible version of $Y_i^\ast(a)$ defined as
This strategy leads to the following variance estimator:
where
Note that the construction of the numerator of $\hat{\nu}^2_n$ can be motivated using a similar intuition to what has been previously discussed in bai2022mp: to consistently estimate \[ E[(E[Y^\ast_i(1)|X_i] - E[Y^\ast_i(0)|X_i])^2]~, \] ideally we would like access to four different units with similar values of $X_i$, of which two are treated. However, because each pair only contains two units, we need to average across “pairs of pairs" of units, where two pairs are grouped together so that they are “close” in terms of $X_i$. We establish the following consistency result for $\hat{\nu}^2_n$:
Assumption (ref)(a) implies that $\nu^2 > 0$. We therefore immediately obtain the following corollary which establishes the asymptotic exactness of a $t$-test for the null hypothesis (ref) constructed using the variance estimator $\hat{\nu}^2_n$:
Next, we consider the limiting behavior of the usual heteroskedasticity-robust variance estimator obtained from a two-stage least squares regression of $Y_i$ on a constant and $D_i$, using $A_i$ as an instrument, which we denote by $\hat \omega_n^2$. Theorem (ref) derives the limit in probability of $\hat \omega_n^2$.
Finally, we consider the heteroskedasticity-robust variance estimator obtained from a two-stage least squares regression with pair fixed effects, which is a specification commonly used in practice: see groh2016macroinsurance, and resnjanskij2021can for examples. Specifically, let $\hat \alpha_n$ be the two-stage least squares estimator for $\alpha$ in the linear regression
using $A_i$ as an instrument for $D_i$. Let $\hat \omega^2_{n,\rm pfe}$ denote the usual heteroskedasticity-robust variance estimator (HC0) for $\hat \alpha_n$ and let $\hat \omega^2_{n,\rm pfe, HC1}$ denote the robust variance estimator with a degree-of-freedom adjustment (HC1). In this particular case,
because the number of regressors in (ref) is $n + 1$.
From Theorems (ref) and (ref) we obtain that neither cluster-robust standard error is consistent for $\nu^2$ unless the baseline covariates are irrelevant for the transformed potential outcomes $Y_i^*(a)$ in an appropriate sense. In Section (ref), we illustrate that this conservativeness translates into a lack of power relative to tests constructed using our consistent variance estimator $\hat{\nu}_n^2$.
In this section, we consider a generalization of the estimator $\hat \Delta_n$ defined in Section (ref) that allows for covariate adjustment using the additional, observed, baseline covariates $W^{(n)}$ that were not used when forming pairs. In Section (ref), we derive general results on covariate adjustment using arbitrarily-specified working models of the conditional expectations of the outcome and treatment take-up with respect to the covariates. In Section (ref), we show how a careful application of these earlier results using an optimal linear covariate-adjusted estimator can ensure an improvement over the unadjusted Wald estimator $\hat{\Delta}_n$ in terms of precision.
Following bai2023covariate, note that it can be shown under Assumption (ref) that for any $a \in \{0, 1 \}$, $m_{a, \tilde Y}: \mathbb{R}^{k_x} \times \mathbb{R}^{k_w} \to \mathbb{R}$, and $m_{a, D}: \mathbb{R}^{k_x} \times \mathbb{R}^{k_w} \to \mathbb{R}$ such that $E[|m_{a, \tilde Y}(X_i, W_i)|] < \infty$, $E[|m_{a, D}(X_i, W_i)|] < \infty$,
Here, we view $m_{a, D}$ and $m_{a, \tilde Y}$ as “working models” of $E[\tilde{Y}_i(a)|X_i,W_i]$ and $E[D_i(a)|X_i,W_i]$, respectively, but we emphasize that (ref)--(ref) hold even if these models are incorrectly specified. These two moment conditions motivate a covariate-adjusted estimator defined as \[ \hat \Delta_n^{\rm adj} = \frac{\hat \psi_n^{\rm adj}(1) - \hat \psi_n^{\rm adj}(0)}{\hat \phi_n^{\rm adj}(1) - \hat \phi_n^{\rm adj}(0)}~, \] where
and $\hat m_{a, \tilde Y}$ and $\hat m_{a, D}$ are suitable estimators of $m_{a, \tilde Y}$ and $m_{a, D}$. Note that if we set $\hat m_{a,\tilde{Y}} = \hat m_{a, D} = 0$, then $\hat{\Delta}_n^{\rm adj}$ simplifies to $\hat{\Delta}_n$. Let \[m_{a, \tilde Y D}(X_i,W_i) = m_{a, \tilde Y}(X_i,W_i) - \Delta(Q) m_{a, D}(X_i,W_i)~.\] Theorem (ref) below establishes the limiting distribution of $\hat{\Delta}_n^{\rm adj}$ for a matched-pairs design under the following high-level assumption on the working models:
In Section (ref) we provide low-level sufficient conditions for Assumptions (ref) and (ref)-(ref) for the special case in which $m_{a,\tilde{Y}}$ and $m_{a,D}$ correspond to linear working models that are optimal in the sense of minimizing the limiting variance of $\hat \Delta_n^{\rm adj}$ among all linear working models. We derive Theorem (ref) by reproducing arguments in Theorem 3.1 in bai2023covariate, but replacing the working model $m_a$ with the transformed working model $m_{a, \tilde YD}$. As a result, the then numerator $\nu_{\rm adj}^2$ corresponds exactly to the limiting variance obtained in bai2023covariate with the working models in bai2023covariate replaced with $m_{a, \tilde YD}$. As expected, $\nu^2_{\rm adj} = \nu^2$ when the working models are set to zero. In general, $\nu^2_{\rm adj}$ is not guaranteed to be weakly smaller than $\nu^2$ for all choices of working models, but we note that $\nu_{\rm adj}^2$ is minimized (and thus smaller than $\nu^2$) when $\nu_{1, \rm adj}^2 = 0$, i.e., when the working models satisfy
with probability one. This property holds, in particular, when $m_{a,\tilde{Y}}$ and $m_{a,D}$ are correctly specified. Moreover, under correct specification, it can be shown that $\nu^2_{\rm adj}$ coincides with the efficiency bound derived in bai2023efficient.
Next, we construct a consistent variance estimator for the limiting variance $\nu_{\rm adj}^2$. As noted in the discussion following Theorem (ref), the numerator for $\nu_{\rm adj}^2$ corresponds exactly to the limiting variance obtained in bai2023covariate with the usual potential outcomes $Y_i(a)$ replaced with the transformed outcomes $Y_i^\ast(a)$ and working models in bai2023covariate replaced with $m_{a, \tilde YD}$. We thus follow the variance construction from bai2023covariate, but with a feasible version of $Y_i^\ast(a)$ defined as \[\hat{Y}_i = Y_i - \hat{\Delta}_nD_i~,\] and a feasible version of $m_{a, \tilde{Y}D}(X_i,W_i) = m_{a, \tilde Y}(X_i, W_i) - \Delta(Q) m_{a, D}(X_i, W_i)$ defined as \[ \hat m_{a, \tilde Y D}(X_i, W_i) = \hat m_{a, \tilde Y}(X_i, W_i) - \hat{\Delta}_n \hat m_{a, D}(X_i, W_i)~. \] This strategy leads to the following variance estimator:
where
The following theorem establishes the consistency of $\hat{\nu}_{n, \rm adj}^2$ for $ \nu_{\rm adj}^2$:
We now consider a simple setting of practical interest which specializes Theorem (ref) to the case in which $m_{a,\tilde{Y}}$ and $m_{a,D}$ are the optimal linear working models in a sense to be made formal below. To that end, let $\zeta_i = \zeta(X_i, W_i)$ be a user-specified transformation of the baseline characteristics $(X_i, W_i)$ given by some function $\zeta: \mathbb R^{k_x} \times \mathbb R^{k_w} \to \mathbb R^p$. Let $\hat m_{a, \tilde Y}(X_i, W_i) = \zeta_i' \hat \beta_n^Y$ and $\hat m_{a, D}(X_i, W_i) = \zeta_i' \hat \beta_n^D$ for $a \in \{0, 1\}$, where $\hat \beta_n^Y$ and $\hat \beta_n^D$ are the ordinary least squares estimators of $\beta^Y$ and $\beta^D$ in the following two linear regressions with pair fixed effects:
Lemma (ref) in the appendix shows that the adjusted estimator $\hat{\Delta}_n^{\rm adj}$ obtained from this choice of $\hat m_{a, \tilde Y}$ and $\hat m_{a, D}$ could alternatively be computed as the regression coefficient $\hat{\alpha}_n^{\rm IV}$ in the following two-stage least squares regression with $A_i$ as an instrument for $D_i$:
Versions of this estimator have been used in, for instance, groh2016macroinsurance and resnjanskij2021can. It is also a natural counterpart for the OLS estimator with pair fixed effects which is widely used in settings with perfect compliance; see, for instance, glennerster2014running. In order to analyze the large-sample behavior of this estimator, we introduce the following assumption:
Theorem (ref) shows that Assumption (ref) provides low-level sufficient conditions for Assumption (ref) and (ref)--(ref), and further verifies the optimality of these working models among all linear working models:
We conclude by noting that, since the unadjusted estimator $\hat{\Delta}_n$ is contained in the class of covariate adjusted estimators for which $m_{a,D}$ and $m_{a,\tilde Y}$ are linear in $\zeta_i$ (by setting both working models to zero), Theorem (ref) establishes that the optimal linear adjusted estimator is guaranteed to be (weakly) more efficient than the unadjusted Wald estimator, regardless of whether or not the linear working models are correctly specified.
In this section, we examine the finite-sample behavior of the estimation and inference procedures introduced in Section (ref). Following the simulation designs in bai2022mp, the potential outcomes and take-up decisions are generated according to the equations:
where $\mu_d$, $m_d(X_i)$, $\sigma_d(X_i)$, and $\epsilon_{d, i}$ are specified in each model below. We consider the following model specifications:
For each model, let $\Delta_0$ denote the value of the LATE when $\mu_1 = 0$.\footnote{For each model, $\Delta_0$ is computed numerically with a large sample. The values are -0.0000203726, 0.0859858425, 0.0903371248 for Models 1-3.} Because $\dim(X_i) = 1$, we construct pairs by sorting units according to $X_i$ and matching adjacent units. Table (ref) reports the rejection probabilities from 5,000 Monte Carlo replications for testing the null hypothesis (ref) against the alternative implied by setting $\mu_1 = 1/2$, using $t$-tests constructed using either the regression-based variance $\hat{\omega}_n^2$, the regression-based variance with pair fixed effects and HC1 correction $\hat{\omega}_{n, \rm pfe, HC1}^2$, or the consistent estimator $\hat{\nu}^2_n$. As expected given our theoretical results, tests based on $\hat \omega_n^2$ and $\hat{\omega}_{n, \rm pfe, HC1}^2$ are generally conservative, which leads to a loss of power under the alternative relative to the test based on our consistent variance estimator $\hat \nu^2_n$. Figure (ref) displays the power curves for testing the null hypothesis (ref) for Models 1--3 with $\mu_1 \in [-1, 1]$ and $2n = 200$. We emphasize that, although $\hat{\omega}_{n, \rm pfe, HC1}^2$ is less conservative than $\hat{\omega}_n^2$ for these simulation designs, our theoretical results in Section (ref) establish that this is not guaranteed to be the case in general.
In this section, we examine the finite-sample behavior of the estimation and inference procedures introduced in Section (ref). Following the simulation designs in bai2023covariate, the potential outcomes and take-up decisions are generated according to the equations:
where $\mu_d$, $m_d(X_i, W_i)$, $\sigma_d(X_i, W_i)$, and $\epsilon_{d, i}$ are specified in each model below. We consider the following model specifications:
As in Section (ref), for each model, let $\Delta_0$ denote the value of the LATE when $\mu_1 = 0$.\footnote{For each model, $\Delta_0$ is computed numerically with a large sample. The values are -0.0007846080, -0.0005474909, -0.0013187170, 0.0224019752 for Models 1-4.} We construct pairs by sorting units according to $X_i$ and matching adjacent units. Table (ref) reports the rejection probabilities from 5,000 Monte Carlo replications for testing the null hypothesis (ref) against the alternative implied by setting $\mu_1 = 1/2$, using $t$-tests constructed using three different linear covariate-adjusted estimators: the unadjusted estimator $\hat{\Delta}_n$ (denoted by unadj), the estimator obtained from a two-stage least squares regression of $Y_i$ on a constant, $D_i$, and covariates $W_i$, but without pair fixed effects (denoted by naive), and the optimal linear estimator described in (ref) with $\zeta_i = W_i$ (denoted by pfe). Note that each of these estimators corresponds to a special case of our regression adjusted estimator $\hat{\Delta}_n^{\rm adj}$ for three different working models.
As expected given our theoretical results in Section (ref), all three tests maintain exact size under the null hypothesis. Power under the alternative hypothesis increases for all model specifications and sample sizes as we move from the unadjusted estimator to the naive linear adjusted estimator and from the naive estimator to the optimal linear adjusted estimator. Table (ref) reports the bias and root mean-squared error (RMSE) for the same set of models and sample sizes. Here, we find that RMSE decreases as we move from the unadjusted estimator to the optimal linear estimator, with no qualitative differences in bias.
In this section, we illustrate our findings by revisiting the empirical application in groh2016macroinsurance. groh2016macroinsurance designed a matched-pair experiment in Egypt to study the effect on microenterprises of acquiring insurance against macroeconomic shocks.\footnote{To most closely align the dataset with our theoretical results, we made the following modifications to the dataset: (1) for each outcome variable, we drop pairs if at least one of the individuals in that pair has a missing outcome variable, (2) we drop pairs if at least one of the individuals in that pair is missing treatment assignment (the eligibility of purchasing the insurance), treatment status (whether a company actually purchased the insurance), or any baseline covariates, (3) we keep only pairs with exactly two individuals (there were 39 pairs with only one individual and one “block" with 16 individuals), (4) if necessary, we drop one pair from the end of the resulting dataset to ensure that the sample size is divisible by 4. (5) to construct the pairs of pairs when computing $\hat{\nu}_{n}$ and $\hat{\nu}_{n,adj}$, we use the R package nbpMatching to match pairs of pairs such that the conditions in Theorem 4.3 of bai2022mp are satisfied. Modifications (1)-(5) result in an average sample size reduction of 96 observations (3.29% of total sample size) across outcomes.} The eligibility to purchase macroinsurance was offered to companies in the treatment group. The take-up rate of purchasing in the treatment group was 37%: among 1481 companies in the treatment group, 548 of them purchased insurance. We also note among 1480 companies in the control group, 5 of them purchased insurance as well.
Table (ref) reports the estimated LATEs for a collection of outcomes, using both the unadjusted Wald estimator $\hat{\Delta}_n$ as well as the linearly adjusted estimator $\hat{\Delta}_n^{\rm adj}$ which uses the same covariates as those considered in the analysis in groh2016macroinsurance.\footnote{We note however that we exclude the female dummy and branchid dummies, which were used in the original regression specifications in groh2016macroinsurance, since these are perfectly collinear with our pair fixed effects.} For the unadjusted estimator $\hat{\Delta}_n$ we report the standard errors obtained from the regression-based variance estimators $\hat{\omega}_n^2$ and $\hat{\omega}_{n, \rm pfe, \rm HC1}^2$ as well as the standard errors obtained from our consistent variance estimators $\hat{\nu}^2_n$. For the adjusted estimator $\hat{\Delta}_n^{\rm adj}$ we report the standard errors obtained from our consistent variance estimator $\hat{\nu}^{2}_{n, \rm adj}$. Our findings are consistent with the theoretical results presented in Sections (ref) and (ref): for the unadjusted estimates, standard errors constructed from $\hat{\nu}^2_n$ are smaller than those constructed from $\hat{\omega}^2_n$ and comparable to those constructed from $\hat{\omega}^2_{n, \rm pfe, HC1}$. Given Theorem (ref), this suggests that there is limited heterogeneity in $E[Y_i^{\ast}(1) - Y_i^{\ast}(0)|X_i]$ in this application. For the adjusted estimates, we find that the standard errors constructed from $\hat{\nu}^2_{n, \rm adj}$ are smaller than those constructed from $\hat{\nu}^2_n$ across all outcomes. However, point estimates for profits and the aggregate index change in such a way that these are no longer significant at the $10\%$ level.
Based on our theoretical results as well as the simulation study above, we conclude with some recommendations for practitioners when conducting inference about the local average treatment effect in matched-pairs experiments. Our findings are that the standard Wald estimator is generally consistent and asymptotically normal under matched-pair designs, but its limiting variance is smaller than what would be obtained under i.i.d.\ assignment. It follows that inferences using standard heteroskedasticty-robust estimators of the variance will typically be conservative. We therefore recommend that practitioners use our consistent variance estimator $\hat{\nu}^2_n$ instead. When considering covariate adjustment, our findings are that the two-stage least squares estimator with pair fixed effects leads to an estimator that is optimal in the sense of having smallest limiting variance among all linearly-adjusted estimators. An important caveat, however, is that the usual heteroskedasticty-robust estimator of the variance is not consistent for its variance. As a result, we recommend that practitioners use our consistent variance estimator $\hat{\nu}^2_{n, \rm adj}$ instead.