Covariate Adjustment in Experiments with Matched Pairs

Covariate Adjustment in Experiments with Matched Pairs

\hypersetup{pageanchor=false}

KEYWORDS: Experiment, matched pairs, covariate adjustment, randomized controlled trial, treatment assignment, LASSO

JEL classification codes: C12, C14

\thispagestyle{empty} \setcounter{page}{1} \hypersetup{pageanchor=true}

Introduction

This paper studies inference on the average treatment effect in experiments in which treatment status is determined according to “matched pairs.” By a “matched pairs” design, 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 surveys of various field experiments: donner2000design, glennerster2013running, and 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. bai2021inference develop methods for inference on the average treatment effect in such experiments based on the difference-in-means estimator. In this paper, we pursue the goal of improving upon the precision of this estimator by exploiting observed, baseline covariates that are not used in determining treatment status.

To this end, we study a broad class of estimators for the average treatment effect based on a “doubly robust” moment condition. The estimators in this framework are distinguished via different “working models” for the conditional expectations of potential outcomes under treatment and control given the observed, baseline covariates. Importantly, because of the double-robustness, these “working models” need not be correctly specified in order for the resulting estimator to be consistent. In this way, the framework permits us to study both finite-dimensional and high-dimensional forms of covariate adjustment without imposing unreasonable restrictions on the conditional expectations themselves. Under high-level conditions on the “working models” and their corresponding estimators and a requirement that pairs are formed so that units within pairs are suitably “close” in terms of the baseline covariates, we derive the limiting distribution of the covariate-adjusted estimator of the average treatment effect. We further construct an estimator for the variance of the limiting distribution and provide conditions under which it is consistent for this quantity.

Using our general framework, we first consider finite-dimensional, linear adjustments. For this class of estimators, our main findings are summarized as follows. First, we find that estimators with such adjustments are not guaranteed to be weakly more efficient than the unadjusted difference-in-means estimator. This finding echoes similar findings by yang2001efficiency and tsiatis2008covariate in settings in which treatment is determined by i.i.d.\ coin flips, and freedman2008regression in a finite population setting in which treatment is determined according to complete randomization. See negi2021revisiting for a succinct treatment of that literature. Moreover, we find that this phenomenon persists even if the adjustments are interacted with treatment. In fact, doing so leads to no changes in precision. In this sense, our results diverge from those in settings with complete randomization and treated fraction one half, where adjustments based on the uninteracted and interacted linear adjustments both guarantee gains in precision. Last, we show that estimators with both uninteracted and interacted linear adjustments with pair fixed effects are guaranteed to be weakly more efficient than the unadjusted difference-in-means estimator.

We then use our framework to consider high-dimensional adjustments based on $\ell_1$ penalization. Specifically, we first obtain an intermediate estimator by using the LASSO to estimate the “working model” for the relevant conditional expectations. When the treatment is determined according to “matched pairs,” however, this estimator need not be more precise than the unadjusted difference-in-means estimator. Therefore, following cohen2020no-harm, we consider, in an additional step, an estimator based on the finite-dimensional, linear adjustment described above that uses the predicted values for the “working model” as the covariates and includes fixed effects for each of the pairs. We show that the resulting estimator improves upon both the intermediate estimator and the unadjusted difference-in-means estimator in terms of precision. Moreover, we provide conditions under which the refitted adjustments attain the relevant efficiency bound derived by armstrong2022asymptotic.

Concurrent with our paper, C23 considers covariate adjustment in experiments in which units are grouped into tuples with possibly more than two units, rather than pairs. Both our paper and C23 find that finite-dimensional, linear regression adjustments with pair fixed effects are guaranteed to improve precision relative to the unadjusted difference-in-means estimator, and show that such adjustments are indeed optimal among all linear adjustments. However, C23 does not pursue more general forms of covariate adjustments, including the regularized adjustments described above. Such results permit us to study nonparametric adjustments as well as high-dimensional adjustments using covariates whose dimension diverges rapidly with the sample size.

The remainder of our 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. In Section (ref), we introduce our general class of estimators based on a “doubly robust” moment condition. Under certain high-level conditions on the “working models” and their corresponding estimators, we derive the limiting behavior of the covariate-adjusted estimator. In Section (ref), we use our general framework to study a variety of estimators with finite-dimensional, linear covariate adjustment. In Section (ref), we use our general framework to study covariate adjustment based on the regularized regression. In Section (ref), we examine the finite-sample behavior of tests based on these different estimators via a small simulation study. We find that covariate adjustment can lead to considerable gains in precision. Finally, in Section (ref), we apply our methods to reanalyze data from an experiment using a “matched pairs” design to study the effect of macroinsurance on microenterprise. Proofs of all results and some details for simulations are given in the Online Supplement.

Setup and Notation

Let $Y_i \in \mathbf R$ denote the (observed) outcome of interest for the $i$th unit, $D_i \in \{0,1\}$ be an indicator for whether the $i$th unit is treated, and $X_i \in \mathbf R^{k_x}$ and $W_i \in \mathbf R^{k_w}$ denote observed, baseline covariates for the $i$th unit; $X_i$ and $W_i$ will be distinguished below through the feature that only the former will be used in determining treatment assignment. Further denote by $Y_i(1)$ the potential outcome of the $i$th unit if treated and by $Y_i(0)$ the potential outcome of the $i$th unit if not treated. The (observed) outcome and potential outcomes are related to treatment status by the relationship

For a random variable indexed by $i$, $A_i$, it will be useful to denote by $A^{(n)}$ the random vector $(A_1, \ldots, A_{2n})$. Denote by $P_n$ the distribution of the observed data $Z^{(n)}$, where $Z_i = (Y_i, D_i,X_i,W_i)$, and by $Q_n$ the distribution of $U^{(n)}$, where $U_i = (Y_i(1),Y_i(0),X_i,W_i)$. Note that $P_n$ is determined by (ref), $Q_n$, and the mechanism for determining treatment assignment. We assume throughout that $U^{(n)}$ consists of $2n$ i.i.d.\ observations, i.e., $Q_n = Q^{2n}$, where $Q$ is the marginal distribution of $U_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 the treatment assignment. Our object of interest is the average effect of the treatment on the outcome of interest, which may be expressed in terms of this notation as

We now describe our assumptions on $Q$. We restrict $Q$ to satisfy the following mild requirement:

Next, we describe our assumptions on the mechanism determining treatment assignment. In order to describe these assumptions more formally, we require some further notation to define the relevant pairs of units. The $n$ pairs may be represented by the sets $$\{\pi(2j-1), \pi(2j)\} \text{ for } j = 1, \ldots, n~,$$ where $\pi = \pi_n(X^{(n)})$ is a permutation of $2n$ elements. Because of its possible dependence on $X^{(n)}$, $\pi$ encompasses a broad variety of different ways of pairing the $2n$ units according to the observed, baseline covariates $X^{(n)}$. Given such a $\pi$, we assume that treatment status is assigned as described in the following assumption:

Following bai2021inference, our analysis will additionally require some discipline on the way in which pairs are formed. Let $\|\cdot\|_2$ denote the Euclidean norm. We will require that units in each pair are “close” in the sense described by the following assumption:

It will at times be convenient to require further that units in consecutive pairs are also “close” in terms of their baseline covariates. One may view this requirement, which is formalized in the following assumption, as “pairing the pairs“ so that they are “close” in terms of their baseline covariates.

bai2021inference provide results to facilitate constructing pairs satisfying Assumptions (ref)--(ref) under weak assumptions on $Q$. In particular, given pairs satisfying Assumption (ref), it is frequently possible to “re-order” them so that Assumption (ref) is satisfied. See Theorem 4.3 in bai2021inference for further details. As in bai2021inference, we highlight the fact that Assumption (ref) will only be used to enable consistent estimation of relevant variances.

Main Results

To accommodate various forms of covariate-adjusted estimators of $\Delta(Q)$ in a single framework, it is useful to note that it follows from Assumption (ref) that for any $d \in \{0,1\}$ and any function $m_{d, n}: \mathbf R^{k_x} \times \mathbf R^{k_w} \rightarrow \mathbf R$ such that $E[|m_{d, n}(X_i,W_i)|] < \infty$,

We note that (ref) is just the augmented inverse propensity score weighted moment for $E[Y_i(d)]$ in which the propensity score is $1/2$ and the conditional mean model is $m_{d,n}(X_i,W_i)$. Such a moment is also “doubly robust.” As the propensity score for the “matched pairs” design is exactly one half, we do not require the conditional mean model to be correctly specified, i.e., $m_{d, n}(X_i,W_i) = E[Y_i(d)|X_i,W_i]$. See, for instance, robins1995analysis. Intuitively, $m_{d, n}$ is the “working model” which researchers use to estimate $E[Y_i(d) | X_i, W_i]$, and can be arbitrarily misspecified because of (ref). Although $m_{d, n}$ will be identical across $n \geq 1$ for the examples in Section (ref), the notation permits $m_{d, n}$ to depend on the sample size $n$ in anticipation of the high-dimensional results in Section (ref). Based on the moment condition in (ref), our proposed estimator of $\Delta(Q)$ is given by

where, for $d \in \{0,1\}$,

and $\hat m_{d, n}$ is a suitable estimator of the “working model” $m_{d, n}$ in (ref).

By some simple algebra, we have\footnote{We thank the referee for this excellent point.}

where

It means our regression adjusted estimator can be viewed as a difference-in-means estimator, but with the “adjusted” outcome $\tilde Y_i$.

We require some new discipline on the behavior of $m_{d, n}$ for $d \in \{0, 1\}$ and $n \geq 1$:

Assumption (ref)(a) is an assumption to rule out degenerate situations. Assumption (ref)(b) is a mild uniform integrability assumption on the “working models.” If $m_{d, n}(\cdot) \equiv m_d(\cdot)$ for $d \in \{0, 1\}$, then it is satisfied as long as $E[m_d^2(X_i, W_i)] < \infty$. Assumption (ref)(c) ensures that units that are “close” in terms of the observed covariates are also “close” in terms of potential outcomes, uniformly across $n \geq 1$.

Theorem (ref) below establishes the limit in distribution of $\hat \Delta_n$. We note that the theorem depends on high-level conditions on $m_{d, n}(\cdot)$ and $\hat m_{d, n}(\cdot)$. In the sequel, these conditions will be verified in several examples.

In order to facilitate the use of Theorem (ref) for inference about $\Delta(Q)$, we next provide a consistent estimator of $\sigma_n(Q)$. Define

where $\tilde Y_i$ is defined in (ref). The variance estimator is given by

The variance estimator in (ref), in particular its component $\hat \lambda_n$, is analogous to the “pairs of pairs” variance estimator in bai2021inference. Such a variance estimator has also been used in abadie2008estimation in a related setting. Note that it can be shown similarly as in Remark 3.9 of bai2021inference that $\hat \sigma_n^2$ in (ref) is nonnegative.

Theorem (ref) below establishes the consistency of this estimator and its implications for inference about $\Delta(Q)$. In the statement of the theorem, we make use of the following notation: for any scalars $a$ and $b$, $[a \pm b]$ is understood to be $[a - b, a + b]$.

Linear Adjustments

In this section, we consider linearly covariate-adjusted estimators of $\Delta(Q)$ based on a set of regressors generated by $X_i \in \mathbf R^{k_x}$ and $W_i \in \mathbf R^{k_w}$. To this end, define $\psi_i = \psi(X_i, W_i)$, where $\psi: \mathbf R^{k_x} \times \mathbf R^{k_w} \to \mathbf R^p$. We impose the following assumptions on the function $\psi$:

Assumption (ref) is analogous to Assumption (ref). Note, in particular, that Assumption (ref)(a) rules out situations where $\psi_i$ is a function of $X_i$ only. See Remark (ref) for a discussion of the behavior of the covariate-adjusted estimators in such situations.

Linear Adjustments without Pair Fixed Effects

Consider the following linear regression model:

Let $\hat \alpha_n^{\rm naive}$, $\hat \Delta_n^{\rm naive}$, and $\hat \beta_n^{\rm naive}$ denote the OLS estimators of $\alpha$, $\Delta$, and $\beta$ in (ref). We call these estimators na\"ive because the corresponding regression adjustment is subject to Freedman's critique and can lead to an adjusted estimator that is less efficient than the simple difference-in-means estimator $\hat \Delta_n^{\rm unadj}$.

It follows from direct calculation that \[ \hat \Delta_n^{\rm naive} = \frac{1}{n} \sum_{1 \leq i \leq 2n} (Y_i - \psi_i' \hat \beta_n^{\rm naive}) (2 D_i - 1)~. \] Therefore, $\hat \Delta_n^{\rm naive}$ satisfies (ref)--(ref) with \[ \hat m_{d, n}(X_i, W_i) = \psi_i' \hat \beta_n^{\rm naive}~. \]

Theorem (ref) establishes (ref) and (ref) for a suitable choice of $m_{d, n}(X_i, W_i)$ for $d \in \{0, 1\}$ and, as a result, the limiting distribution of $\hat \Delta_n^{\rm naive}$ and the validity of the variance estimator.

Linear Adjustments with Pair Fixed Effects

Remark (ref) implies that in “matched pairs” designs, including interaction terms in the linear regression does not lead to an estimator with lower limiting variance than the one based on the linear regression without interaction terms. It is therefore interesting to study whether there exists a linearly covariate-adjusted estimator with lower limiting variance than the ones based on (ref) and (ref) as well as the difference-in-means estimator. To that end, consider instead the following linear regression model:

Let $\hat \Delta_n^{\rm pfe}$, $\hat \beta_n^{\rm pfe}$, and $\hat \gamma_{j, n}$, $1 \leq j \leq n$ denote the OLS estimators of $\Delta$, $\beta$, $\theta_j$, $1 \leq j \leq n$ in (ref), where “pfe” stands for pair fixed effect. It follows from the Frisch-Waugh-Lovell theorem that \[ \hat \Delta_n^{\rm pfe} = \frac{1}{n} \sum_{1 \leq i \leq 2n} (Y_i - \psi_i' \hat \beta_n^{\rm pfe}) (2 D_i - 1)~. \] Therefore, $\hat \Delta_n^{\rm pfe}$ satisfies (ref)--(ref) with \[ \hat m_{d, n}(X_i, W_i) = \psi_i' \hat \beta_n^{\rm pfe}~. \]

Theorem (ref) establishes (ref) and (ref) for a suitable choice of $m_{d, n}(X_i, W_i), d \in \{0, 1\}$ and, as a result, the limiting distribution of $\hat \Delta_n^{\rm pfe}$ and the validity of the variance estimator.

Regularized Adjustments

In this section, we study covariate adjustments based on the $\ell_1$-regularized linear regression. Such settings can arise if the covariates $W_i$ are high-dimensional or if the dimension of $W_i$ is fixed but the regressors include many sieve basis functions of $X_i$ and $W_i$. To accommodate situations where the dimension of $W_i$ increases with $n$, we add a subscript and denote it by $W_{n, i}$ instead. Let $k_{w, n}$ denote the dimension of $W_{n, i}$. For $n \geq 1$, let $\psi_{n,i} = \psi_n(X_i,W_{n,i})$, where $\psi_n: \mathbf R^{k_x} \times \mathbf R^{k_{w, n}} \to \mathbf R^{p_n}$ and $p_n$ will be permitted below to be possibly much larger than $n$.

In what follows, we consider a two-step method in the spirit of cohen2020no-harm. In the first step, an intermediate estimator, $\hat \Delta_n^{\rm r}$, is obtained using (ref) with a “working model” obtained through a $\ell_1$-regularized linear regression adjustments $m_{d,n}(X_i,W_i)$ for $d \in \{0,1\}$. As explained further below in Theorem (ref), when $m_{d,n}(X_i,W_i)$ is approximately correctly specified, such an estimator is optimal in the sense that it minimizes the limiting variance in Theorem (ref). When this is not the case, however, for reasons analogous to those put forward in Remark (ref), it needs not to have a limiting variance weakly smaller than the unadjusted difference-in-means estimator. In a second step, we therefore consider an estimator by refitting a version of (ref) in which the covariates $\psi_i$ are replaced by the regularized estimates of $m_{d,n}(X_i,W_i)$ for $d \in \{0,1\}$. The resulting estimator, $\hat \Delta_n^{\rm refit}$, has the limiting variance weakly smaller than that of the intermediate estimator and thus remains optimal under approximately correct specification in the same sense. Moreover, it has limiting variance weakly smaller than the unadjusted difference-in-means estimator. WDTT16 also consider high-dimensional regression adjustments in randomized experiments using LASSO. We differ from their work by considering the “matched pairs” design, and more importantly, discussing when and how regularized adjustments can improve estimation efficiency upon the difference-in-means estimator.

Before proceeding, we introduce some additional notation that will be required in our formal description of the methods. We denote by $\psi_{n,i,l}$ the $l$th components of $\psi_{n,i}$. For a vector $a \in \mathbf R^k$ and $0 \leq p \leq \infty$, recall that $$\|a\|_p = \Big ( \sum_{1 \leq l \leq k} |a_l|^p \Big )^{1/p}~,$$ where it is understood that $\|a\|_0 = \sum_{1 \leq l \leq k} I \{a_k \neq 0\}$ and $\|a\|_\infty = \sup_{1 \leq l \leq k} |a_l|$. Using this notation, we further define $$\Xi_n = \sup_{(x,w) \times \mathrm{supp}(X_i) \times \mathrm{supp}(W_i) }\|\psi_{n,i}(x,w)\|_\infty ~.$$

For $d \in \{0, 1\}$, define

where $\lambda_{d, n}^{\rm r}$ is a penalty parameter that will be disciplined by the assumptions below, $\hat{\Omega}_n(d) = \operatorname*{diag}(\hat{\omega}_1(d),\cdots,\break\hat{\omega}_{p_n}(d))$ is a diagonal matrix, and $\hat{\omega}_{n,l}(d)$ is the penalty loading for the $l$th regressor. Let $\hat \Delta_n^{\rm r}$ denote the estimator in (ref) with $\hat m_{d, n}(X_i,W_{n,i}) = \psi_{n, i}' \hat \beta_{d, n}^{\rm r}$ for $d \in \{0, 1\}$.

We now proceed with the statement of our assumptions. The first assumption collects a variety of moment conditions that will be used in our formal analysis:

Our analysis will, as before, also require some discipline on the way in which pairs are formed. For this purpose, Assumption (ref) will suffice, but we will need an additional Lipshitz-like condition:

We next specify our restrictions on the penalty parameter $\lambda_{d, n}^{\rm r}$.

We note that Assumption (ref)(b) permits $p_n$ to be much greater than $n$. It also requires sparsity in the sense that $s_n = o(\sqrt{n})$.

Finally, as is common in the analysis of $\ell_1$-penalized regression, we require a “restricted eigenvalue” condition. This assumption permits us to apply bickel2009simultaneous and establish the error bounds for $|\hat \alpha_{d,n}^{\rm r} - \alpha_{d,n}^{\rm r}|+||\hat \beta_{d,n}^{\rm r} - \beta_{d,n}^{\rm r}||_1$ and $\frac{1}{n}\sum_{1\leq i\leq 2n}I\{D_i=d\}\left(\hat \alpha_{d,n}^{\rm r} - \alpha_{d,n}^{\rm r}+\psi_{n,i}'(\hat \beta_{d,n}^{\rm r} - \beta_{d,n}^{\rm r})\right)^2$.

Using these assumptions, the following theorem characterizes the behavior of $\hat \Delta_n^{\rm r}$:

Theorem (ref) below shows the “refit” estimator for the ATE is weakly more efficient than both $\hat \Delta_n^{\rm unadj}$ and $\hat \Delta_n^{\rm r}$. To state the results, define $\Gamma_{n,i} = (\psi_{n, i}' \beta_{1, n}^{\rm r}, \psi_{n, i}' \beta_{0, n}^{\rm r})'$, $\hat \Gamma_{n,i} = (\psi_{n, i}' \hat \beta_{1, n}^{\rm r}, \psi_{n, i}'\hat \beta_{0, n}^{\rm r})$, and $\hat \Delta_n^{\rm refit}$ as the estimator in (ref) with $\psi_i$ replaced by $\hat \Gamma_{n,i}$. Note that $\hat \Delta_n^{\rm refit}$ remains numerically the same if we include the intercept $\hat \alpha_{d,n}^{\rm r}$ in the definition of $\hat \Gamma_{n,i}$. Following Remark (ref), $\hat \Delta_n^{\rm refit}$ is the intercept in the linear regression of $(D_{\pi(2j-1)}-D_{\pi(2j-1)})(Y_{\pi(2j-1)} - Y_{\pi(2j)})$ on constant and $(D_{\pi(2j-1)}-D_{\pi(2j-1)})(\hat \Gamma_{n,\pi(2j-1)} - \hat \Gamma_{n,\pi(2j)})$. Replacing $\hat \Gamma_{n,i}$ by $\hat \Gamma_{n,i} + (\hat \alpha_{1,n}^{\rm r},\hat \alpha_{0,n}^{\rm r})'$ will not change the regression estimators.

The following assumption will be employed to control $\Gamma_{n,i}$ in our subsequent analysis:

The following theorem characterizes the behavior of $\hat \Delta_n^{\rm refit}$:

Simulations

In this section, we conduct Monte Carlo experiments to assess the finite-sample performance of the inference methods proposed in the paper. In all cases, we follow bai2021inference to consider tests of the hypothesis that

with $\Delta_{0}=0$ at nominal level $\alpha=0.05$.

Data Generating Processes

We generate potential outcomes for $d\in\{0,1\}$ and $1\leq i\leq2n$ by the equation

where $\mu_{d},m_{d}\left(X_{i},W_{i}\right),\sigma_{d}\left(X_{i},W_{i}\right)$, and $\epsilon_{d,i}$ are specified in each model as follows. In each of the specifications, ($X_{i},W_{i},\epsilon_{0,i},\epsilon_{1,i}$) are i.i.d. across $i$. The number of pairs $n$ is equal to 100 and 200. The number of replications is 10,000.

It is worth noting that Models 1, 2, 3, 4, 7, 10, 12, and 13 imply homogeneous treatment effects because $m_{1}\left(X_{i},W_{i}\right) = m_{0}\left(X_{i},W_{i}\right)$. Among them, $E[Y_i(a)|X_i,W_i] - E[Y_i(a)|X_i]$ is linear in $W_i$ in Models 1, 2, and 12. Models 5, 8, 11, and 14 have heterogeneous but homoscedastic treatment effects. In Models 6, 9, and 15, however, the implied treatment effects are both heterogeneous and heteroscedastic. Models 12-15 contain high-dimensional covariates.

We follow bai2021inference to match pairs. Specifically, if $\text{dim}\left(X_i\right)=1$, we match pairs by sorting $X_i, i = 1, \ldots, 2n$. If $\text{dim}\left(X_i\right)>1$, we match pairs by the permutation $\pi$ calculated using the {\tt R} package nbpMatching. For more details, see bai2021inference. After matching the pairs, we flip coins to randomly select one unit within each pair for treatment and another for control.

Estimation and Inference

We set $\mu_0=0$ and $\mu_{1}=\Delta$, where $\Delta=0$ and $1/4$ are used to illustrate the size and power, respectively. Rejection probabilities in percentage points are presented. To further illustrate the efficiency gains obtained by regression adjustments, in Figure (ref), we plot the average standard error reduction in percentage relative to the standard error of the estimator without adjustments for various estimation methods.

Specifically, we consider the following adjusted estimators.

See Section (ref) in the Online Supplement for the regressors used in the regularized adjustments.

For Models 1-11, we examine the performance of estimators (i)-(v). For Models 12-15, we assess the performance among estimators (i) and (v) in high-dimensional settings. Note that the adjustments are misspecified for almost all the models. The only exception is Model 1, for which the linear adjustment in $W_i$ is correctly specified because $m_d(X_i,W_i)$ is just a linear function of $W_i$.

Simulation Results

Tables (ref) and (ref) report rejection probabilities at the 0.05 level and power of the different methods for Models 1--11 when $n$ is 100 and 200, respectively. Several patterns emerge. First, for all the estimators, the rejection rates under $H_0$ are close to the nominal level even when $n=100$ and with misspecified adjustments. This result is expected because all the estimators take into account the dependence structure arising in the “matched pairs” design, consistent with the findings in bai2021inference.

Second, in terms of power, “pfe” is higher than “unadj”, “na\"ive”, and “na\"ive2” for all eleven models, as predicted by our theory. This finding confirms that “pfe” is the optimal linear adjustment and will not degrade the precision of the ATE estimator. In contrast, we observe that “na\"ive” and “na\"ive2” in Model 3 are even less powerful than the unadjusted estimator “unadj”. Figure (ref) further confirms that these two methods inflate the estimation standard error. This result echoes Freedman's critique freedman2008regression that careless regression adjustments may degrade the estimation precision. Our “pfe” addresses this issue because it has been proven to be weakly more efficient than the unadjusted estimator.

Third, the improvement of power for “pfe” is mainly due to the reduction of estimation standard errors, which can be more than 50% as shown in Figure (ref) for Models 4--9. This means that the length of the confidence interval of the “pfe” estimator is just half of that for the “unadj” estimator. Note the standard error of the “unadj” estimator is the one proposed by bai2021inference, which has already been adjusted to account for the cross-sectional dependence created in pair matching. The extra 50% reduction is therefore produced purely by the regression adjustment. For Models 10-11, the reduction of standard errors achieved by “pfe” is more than 40% as well. For Model 1, the linear regression is correctly specified so that all three methods achieve the global minimum asymptotic variance and maximum power. For Model 2, $m_d(X_i,W_i) - E[m_d(X_i,W_i)|X_i] = \gamma (W_i - E[W_i|X_i])$ so that the linear adjustment $\gamma W_i$ satisfies the conditions in Theorem (ref). Therefore, “pfe”, as the best linear adjustment, is also the best adjustment globally, achieving the global minimum asymptotic variance and maximum power. In contrast, “na\"ive” and “na\"ive2” are not the best linear adjustment and therefore less powerful than “pfe” because of the omitted pair dummies.

Finally, the “refit” method has the best power for most models as they automatically achieve the global minimum asymptotic variance when the dimension of $W_i$ is fixed.

Tables (ref) and (ref) report the size and power for the “refit” adjustments when both $W_i$ and $X_i$ are high-dimensional. We see that the size under the null is close to the nominal 5% while the power for the adjusted estimator is higher than the unadjusted one. Figure (ref) further illustrates the reduction of the standard error is more than 30% for all high-dimensional models.

\newcolumntype{L}{>{\arraybackslash}X} \newcolumntype{C}{>{\arraybackslash}X}

Empirical Illustration

In this section, we revisit the randomized experiment with a matched pairs design conducted in groh2016macroinsurance. In the paper, they examined the impact of macroinsurance on microenterprises. Here, we apply the covariate adjustment methods developed in this paper to their data and reinvestigate the average effect of macroinsurance on three outcome variables: the microenterprises' monthly profits, revenues, and investment.

The subjects in the experiment are microenterprise owners, who were the clients of the largest microfinance institution in Egypt. In the randomization, after an exact match of gender and the institution's branch code, those clients were grouped into pairs by applying an optimal greedy algorithm to additional 13 matching variables. Within each pair, a macroinsurance product was then offered to one randomly assigned client, and the other acted as a control. Based on the pair identities and all the matching variables, we re-order the pairs in our sample according to the procedure described in Section 5.1 of Jiang2022QTE. The resulting sample contains 2824 microenterprise owners, that is, 1412 pairs of them.\footnote{See groh2016macroinsurance and Jiang2022QTE for more details.}

Table (ref) reports the ATEs with the standard errors (in parentheses) estimated by different methods. Among them, “GM” corresponds to the method used in groh2016macroinsurance.\footnote{groh2016macroinsurance estimated the effect by regression with regressors including some baseline variables, a dummy for missing observations, and dummies for the pairs. Specifically, for profits and revenues, the regressors are the baseline value for the outcome of interest, a dummy for missing observations, and pair dummies; for investment, the regressors only include pair dummies. The standard errors for the “GM” ATE estimate are calculated by the usual heteroskedastity-consistent estimator. The “GM” results in Table (ref) were obtained by applying the Stata code provided by groh2016macroinsurance.} The description of other methods is similar to that in Section (ref).\footnote{Specifically:

} The results in this table prompt the following observations.

First, aligning with our theoretical and simulation findings, we observe that the standard errors associated with the covariate-adjusted ATEs, particularly those for the “naïve2” and “pfe” estimates, are generally lower compared to the ATE estimate without any adjustment. This pattern is consistent across nearly all the outcome variables. To illustrate, when examining the revenue outcome, the standard errors for the “pfe” estimates are 10.2% smaller than those for the unadjusted ATE estimate.

Second, the standard errors of the “refit” estimates are consistently smaller than those of the unadjusted ATE estimate across all the outcome variables. For example, when profits are the outcome variable, the “refit” estimates exhibit standard errors 7.5% smaller than those of the unadjusted ATE estimate. Moreover, compared with those of the “pfe” estimates, the standard errors of “refit” are slightly smaller.

Conclusion

This paper considers covariate adjustment for the estimation of average treatment effect in “matched pairs” designs when covariates other than the matching variables are available. When the dimension of these covariates is low, we suggest estimating the average treatment effect by a linear regression of the outcome on treatment status and covariates, controlling for pair fixed effects. We show that this estimator is no worse than the simple difference-in-means estimator in terms of efficiency. When the dimension of these covariates is high, we suggest a two-step estimation procedure: in the first step, we run $\ell_1$-regularized regressions of outcome on covariates for the treated and control groups separately and obtain the fitted values for both potential outcomes, and in the second step, we estimate the average treatment effect by refitting a linear regression of outcome on treatment status and regularized adjustments from the first step, controlling for the pair fixed effects. We show that the final estimator is no worse than the simple difference-in-means estimator in terms of efficiency. When the conditional mean models are approximately correctly specified, this estimator further achieves the minimum variance as if all relevant covariates are used to form pairs in the experiment design stage. We take the choice of variables to use in forming pairs as given and focus on how to obtain more efficient estimators of the average treatment effect in the analysis stage. Our paper is therefore silent on the important question of how to choose the relevant matching variables in the design stage. This topic is left for future research.