Exact Bias Correction for Linear Adjustment of Randomized Controlled Trials

Randomized Controlled Trials (RCTs) are popular in empirical economics [AP08, DGK07, Gle17, LR11]. When estimating average treatment effects, adjustment for pretreatment covariates with linear regression is a commonly recommended practice because it can reduce the variability of estimates. However, adjusting for covariates remains somewhat controversial, in large part because of an influential critique from David Freedman [Fre08a, Fre08b].

Scholars have worked to address Freedman’s critiques and to understand the extent that they can and do matter for empirical practice in empirical economics. Using Freedman’s own framework, [Lin13a] showed that arguments 1 and 2 were resolved by small modifications to practice. Freedman’s efficiency result may be addressed through simple modifications to the regression specification, namely including treatment by covariate interactions [Bli73, Oax73]. Then it can be shown that the adjusted estimator is never less asymptotically efficient than the unadjusted estimator. Regarding argument 2, Lin proves that robust standard errors ([Whi80, Hub67, Eic67], see also [SA12]) are asymptotically conservative in Freedman’s setting, guaranteeing the validity of large-sample inference. On argument 3, [Lin13a] (see also [Lin13b]) notes that the leading term of the bias is in fact estimable and can be shown to be small in a real-world empirical example. However, the small-sample bias of the regression estimator was not yet fully resolved.

Since [Lin13a], there have been notable papers that have proposed unbiased regression-type estimators for experimental data. [MSY13] demonstrate that if the regression model is fully saturated (see also [AI17] and [Imb10]), then the associated effect estimate is unbiased conditional on the event that treatment is not collinear with any covariate stratum. This approach cannot generally be used without coarsening continuous covariates. [AM13] proposes the use of auxiliary data, demonstrating that the suitable use of hold-out samples ensures the finite-sample unbiasedness of the associated regression estimator, but the paper does not consider efficiency properties. More recently, [WGB18] extended [AM13] to propose an innovative but computationally expensive split-sample approach for completely randomized experiments.

The primary contribution of this paper is to resolve Freedman’s third theoretical argument by proposing finite-sample-exact, closed-form bias corrections without adding any new assumptions. Our idea builds on [Lin13a]’s proposal to estimate the leading term of the bias, but further develops a novel finite-sample exact bias correction encompassing all higher-order terms [Fre08a]. We derive these bias corrections for both the noninteracted and interacted linear regression estimators. We prove that the estimators have the same limiting distributions as the non-bias-adjusted estimators, implying that they could replace existing estimators in instances where bias is a prevailing concern (e.g., trials that may be aggregated in meta-analysis). We further provide a numerical illustration demonstrating these properties.
Finally, we remind the readers that the practice of debiasing estimators is not uncontroversial. [TE93]²²2We thank Winston Lin for suggesting this reference. has warned that the bias correction could be dangerous in practice due to its high variability³³3As will be shown in the simulations, when the performance is measured by the Root Mean Squared Error (RMSE), there is no clear dominance among the estimators: in some cases the RMSEs of the debiased estimators are strictly smaller than those of other estimators, and in other cases larger.. In real-world decision making processes, people may express different preferences for different statistical properties (i.e. unbiasedness or low Mean Squared Error)⁴⁴4See [WGB18] for an anecdotal example of a policy-maker favoring unbiasedness.. Our results shall imply that in large samples the additional variation caused by the bias correction is negligible, but for small samples, in some cases, we find it important to account for the sampling variability of the additional terms. To address this problem, we propose a simple modification to the standard error estimation procedure. Such modification, based on recomputing OLS residuals using the debiased estimators, is shown to work well on our simulated datasets. We make recommendations for practice in the Simulation section.
The organization of the paper is as follows: Section 2 includes the model setup and assumptions; Section 3 considers a characterization of bias terms of the OLS estimators; Section 4 proposes the bias corrections; Section 5 includes simulation results with both simulated datasets and a real world dataset. In the appendix one can find the proofs for theorems in Section 3 and Section 4, and more simulation results.

We follow the setting of [Fre08a] and [Lin13a], which assume a Neyman [SNDS90] model with covariates. There are $n$ subjects indexed by $i=1,...,n$. For each subject we observe an outcome $Y_{i}$ and a column vector of covariates $\textbf{z}_{i}=(z_{i1},z_{i2},...,z_{iK})\in\mathbb{R}^{K}$. The dimension of the covariates, $K$, does not change with the sample size.
Each subject has two potential outcomes $a_{i}$ and $b_{i}$ (cf., the stable-unit-treatment-value-assumption [Rub90]). We observe $Y_{i}=a_{i}$ if $i$ is chosen for treatment arm $A$ (treated group) and $Y_{i}=b_{i}$ if i is chosen for arm $B$ (control group). Let $T_{i}$ be the dummy variable for treatment arm $A$. Thus the observed outcome for $i$ is $Y_{i}=a_{i}T_{i}+b_{i}(1-T_{i})$.
The experiment is assumed to be completely randomized: $n_{A}$ out of $n$ subjects are randomly assigned to arm $A$ and the remaining $n_{B}=n-n_{A}$ subjects to arm $B$. Random assignment is the only source of randomness in the model. We do not assume a superpopulation: the $n$ subjects are the population of interest.
We introduce some notation. let $n$ be the population size, $n_{A}$ and $n_{B}$ be the number of subjects in treatment arms $A$ and $B$, respectively. Let $[A]=\{i\mid T_{i}=1\}$ denote the set of individuals who are chosen for arm $A$ and similarly $[B]=\{i\mid T_{i}=0\}$. Let $\bar{x}=\frac{1}{N}\sum_{i=1}^{N}x_{i}$, $\bar{x}_{A}=\frac{1}{n_{A}}\sum_{i\in[A]}x_{i}$ and $\bar{x}_{B}=\frac{1}{n_{B}}\sum_{i\in[B]}x_{i}$ denote the population average, group $A$ average, and group $B$ average, respectively, of possibly a vector-valued variable $x$. The average treatment effect (ATE) can be written in this notation as:

and the difference-in-means estimator:

Simiarly we can write $\frac{1}{n}\sum_{i=1}^{n}\mathbf{z}_{i}\mathbf{z}_{i}^{\prime}=\overline{\mathbf{z}\mathbf{z}^{\prime}}$ for $\mathbf{z}_{i}\in\mathbb{R}^{K}$ and $\frac{1}{n}\sum_{i=1}^{n}a_{i}\mathbf{z}_{i}=\overline{a\mathbf{z}}$ for $a_{i}\in\mathbb{R}$ and $\mathbb{z}\in\mathbb{R}^{K}$.
We make the following assumptions throughout the paper, which are standard in the literature. (cf., [Fre08b, Fre08a, Lin13a]).

All assumptions are employed regularly in the literature. They are used to derive consistency and asymptotic normality for the estimator. Assumption 3 requires each arm receives a nontrivial fraction of subjects over the asymptotic sequence of the models. Assumption 4 is without loss of generality: in practice, researchers can just demean each covariate and apply our method.
We remind the readers of the definitions of our two OLS regression adjusted ATE estimators. The first estimator comes from a noninteracted OLS regression where one regresses observed outcome $Y$ on the treatment indicator $T$ and demeaned pretreatment covariates $Z$. The coefficient estimate for $T$ is the noninteracted OLS regression adjusted ATE estimator. The second estimator comes from a fully interacted OLS regression where one regresses observed outcome $Y$ on the treatment indicator $T$, demeaned pretreatment covariates $Z$, and interaction terms of the treatment indicators and demeaned pretreatment covariates. The coefficient estimate for $T$ is the interacted OLS regression adjusted ATE estimator.
Finally, we prepare some notation for the sections below. Let $a^{*}_{i}$ and $b^{*}_{i}$ be the centered potential outcomes, namely, $a^{*}_{i}=a_{i}-\bar{a}$ and $b^{*}_{i}=b_{i}-\bar{b}$. Let $\widehat{D}=\overline{\mathbf{z}\mathbf{z}^{\prime}}-{p}_{A}\overline{\mathbf{z}}_{A}\overline{\mathbf{z}}_{A}^{\prime}-{p}_{B}\overline{\mathbf{z}}_{B}\overline{\mathbf{z}}_{B}^{\prime}$, $\widehat{N}={p}_{A}(\overline{a\mathbf{z}}_{A}-\overline{a}_{A}\overline{\mathbf{z}}_{A})+{p}_{B}(\overline{b\mathbf{z}}_{B}-\overline{b}_{B}\overline{\mathbf{z}}_{B})$, $D=\overline{\mathbf{z}\mathbf{z}^{\prime}}$ and $N={p}_{A}\overline{a\mathbf{z}}+{p}_{B}\overline{b\mathbf{z}}$. With this notation the regression coefficients estimators of the pretreatment covariates in the noninteracted case can be written as $\widehat{Q}={\widehat{D}}^{-1}\widehat{N}$, and the population coefficients $Q=D^{-1}N$. Denote the (rescaled) leverage of $i$th data point as $h_{i}=z_{i}^{\prime}D^{-1}z_{i}$.
Further define $\widehat{D}_{A}=\overline{\mathbf{z}\mathbf{z}^{\prime}}_{A}-\overline{\mathbf{z}}_{A}\overline{\mathbf{z}}_{A}^{\prime}$, $\widehat{N}_{A}=\overline{a\mathbf{z}}_{A}-\overline{a}_{A}\overline{\mathbf{z}}_{A}$, $\widehat{D}_{B}=\overline{\mathbf{z}\mathbf{z}^{\prime}}_{B}-\overline{\mathbf{z}}_{B}\overline{\mathbf{z}}_{B}^{\prime}$, $\widehat{N}_{B}=\overline{b\mathbf{z}}_{B}-\overline{b}_{B}\overline{\mathbf{z}}_{B}$, $N_{A}=\overline{a\mathbf{z}}$ and $N_{B}=\overline{b\mathbf{z}}$. Thus the regression coefficients estimator of the pretreatment covariates in the interactive case can be written as $\widehat{Q}_{A}=\widehat{D}^{-1}_{A}\widehat{N}_{A}$ and $\hat{Q}_{B}=\widehat{D}^{-1}_{B}\widehat{N}_{B}$. Their population counterparts are $Q_{A}=D^{-1}N_{A}$ and $Q_{B}=D^{-1}N_{B}$.

As shown in [Lin13a], the OLS regression adjusted ATE estimator can be written as:

for the noninteracted case and

for the interacted case, where $\widehat{Q}$, $\widehat{Q}_{A}$ and $\widehat{Q}_{B}$ are the OLS coefficients in front of the covariates.

From the coefficient decomposition one can directly characterize the bias term of the ATE estimator. Note that the bias terms are of order $O_{p}(n^{-1})$.

Note that this result implies that the bias terms of the interactive ATE estimator are also of order $O_{p}(n^{-1})$.

Having established the decomposition, we now derive estimators of each bias term for use as bias corrections. We show that these bias estimates are (i) exactly unbiased and (ii) have estimation error of $O_{p}(n^{-1})$. It follows that use of this bias correction with an adjusted estimator yields a finite-sample unbiased estimator with the limit distribution of the adjusted estimator. We remind the readers that $h_{i}=z_{i}^{\prime}D^{-1}z_{i}$, the (rescaled) leverage of $i$th data point as defined in Section 2.

We again begin with the noninteracted case.

The results are derived analogously in the interacted case.

In this section we apply our estimators on several datasets.

We briefly comment on variance estimation and confidence interval construction. We showed in Section 4 that our debiased estimators have the same asymptotic distributions as those of the OLS estimators. This implies that for large samples we can recenter the OLS confidence intervals with our debiased estimators and expect the same coverage probabilities. For small samples, however, we find it important to account for the sampling variability of the additional terms. Indeed, in one of the simulations below, we found that a naive recentering procedure may lead to severe undercoverage. To address this problem, we propose a simple procedure shown to work well on our simulated datasets.⁵⁵5Another way is to directly estimate the variances of the additional terms, but this may be cumbersome to do. In this procedure, one first runs the OLS regression and computes the debiased estimator. Then one replaces the OLS treatment coefficient with the debiased estimated coefficients and recomputes the OLS residuals, keeping all other coefficients the same. Finally, one computes the variances and constructs confidence intervals for the debiased estimator using the same formula as for the OLS estimators. In the simulations below, such procedures will be denoted by BC, which stands for bias-corrected. In Appendix B, one can find a more detailed comparison of this new procedure with standard ones. In practice, we recommend researchers to use our debiased estimators with this procedure, the BC-HC2 heteroskedasticity-robust variance estimator with a Satterthwaite adjustment for inference⁶⁶6For a discussion of the Satterthwaite adjustment, see [Sat46],[BM02], [Lin13a] and [Imb10]..