Exact Bias Correction for Linear Adjustment of Randomized Controlled Trials

Exact Bias Correction for Linear Adjustment of Randomized Controlled Trials

Introduction

Randomized Controlled Trials (RCTs) are popular in empirical economics angrist2008mostly,duflo2007using,glennerster2017practicalities,list2011field. 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 freedman2008A,freedman2008B.

Freedman argued that randomization does not justify the use of linear regression for completely randomized experiments. Freedman's theoretical arguments relied on three results proven under the randomization-based splawa1990application,imbens2015causal inferential paradigm:

Freedman's arguments were influential among scholars across multiple disciplines (e.g., dunning2012natural,BenRecht). Freedman's third argument garnered particular attention among social scientists. Notably, deaton2018understanding's critique of randomization in empirical economics argued that the bias introduced by regression undermines the gold standard argument for RCTs.

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, lin2013agnostic 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 blinder1973wage,oaxaca1973male. 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 (white1980heteroskedasticity,Hube67,eicker1967limit, see also samii2012equivalencies) are asymptotically conservative in Freedman's setting, guaranteeing the validity of large-sample inference. On argument 3, lin2013agnostic (see also lin2013essays) 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 lin2013agnostic, there have been notable papers that have proposed unbiased regression-type estimators for experimental data. miratrix2013adjusting demonstrate that if the regression model is fully saturated (see also athey2017econometrics and imbens2010better), 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. aronow2013class 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, wu2018loop extended aronow2013class 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 lin2013agnostic’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 freedman2008A. 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. tibshirani1993introduction\footnote{We thank Winston Lin for suggesting this reference.} has warned that the bias correction could be dangerous in practice due to its high variability\footnote{As 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)\footnote{See wu2018loop 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.

Setting, Assumptions and Notations

We follow the setting of freedman2008A and lin2013agnostic, which assume a Neyman splawa1990application 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 RUBIN1990279). 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_iT_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'=\overline{\mathbf{z}\mathbf{z}'}$ for $\mathbf{z}_i\in\mathbb{R}^K$ and $\frac{1}{n}\sum_{i=1}^na_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., freedman2008B,freedman2008A,lin2013agnostic).

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}'} - {p}_A \overline{\mathbf{z}}_A \overline{\mathbf{z}}_A'- {p}_B \overline{\mathbf{z}}_B \overline{\mathbf{z}}_B'$, $\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}'}$ 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'D^{-1}z_i$.\\ Further define $\widehat{D}_A= \overline{\mathbf{z}\mathbf{z}'}_A - \overline{\mathbf{z}}_A\overline{\mathbf{z}}_A'$, $\widehat{N}_A = \overline{a\mathbf{z}}_A-\overline{a}_A\overline{\mathbf{z}}_A$, $\widehat{D}_B= \overline{\mathbf{z}\mathbf{z}'}_B - \overline{\mathbf{z}}_B\overline{\mathbf{z}}_B'$, $\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$.

Bias Characterization

As shown in lin2013agnostic, 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.\\

A characterization of the bias terms is provided in this section. Note that both the noninteracted and interacted estimators can be written as sums of the difference-in-means estimator and a regression adjustment using group means and OLS coefficients. The bias comes from the regression adjustment terms, in particular from estimating the regression coefficients on the covariates. We shall characterize the bias terms of the coefficient estimator first. We start from the noninteracted case. Note from here on we shall assume for simplicity that all design matrices (i.e. $\hat{D}$ and $D$) are invertible. In case of noninvertible design matrices, our debiased procedure will still work after choosing a particular generalized inverse, and compute the ATE estimators according to the formulae above.

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})$.

Following the same steps as we did for the noninteracted estimator, we are able to derive analogous results for the interacted estimator.

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

Bias Corrections for Regression Components

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'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.

Simulations

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.\footnote{Another 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\footnote{For a discussion of the Satterthwaite adjustment, see satterthwaite1946approximate,bell2002bias, lin2013agnostic and imbens2010better.}.

Simulated Datasets

In this section we compare the performance of our debiased estimators with that of standard estimators using simulated datasets. We show the results of two simulation schemes here\footnote{In Appendix B, readers can find results for two more simulation schemes as well as some graphical information of the data generating processes.}. In each scheme, we first generate two dimensional covariates that are the quantiles of some prespecified distributions\footnote{For example, in scheme 1, with a sample of $N$ units, the covariates of $i$th unit are the $\frac{i}{N+1}$ quantiles of a Beta(0.5,0.5) distribution and a Triangle(0,1) distribution.}. We then compute the studentized leverage ratios and use those to impute potential outcomes. We consider three ways to impute potential outcomes. For all cases the average treatment effect is equal to 0. The experiment is a completely randomized experiment with 24 units and an inclusion probability of $\frac{1}{3}$ for the treatment arm. Table (ref) includes simulation details. Note that these schemes are designed specifically such that the finite sample bias is relatively large. \\ Table 2 and 3 report the simulation results for the two schemes. Our debiased estimators are exactly unbiased as expected. In terms of the root mean squared error (RMSE) the picture is less clear. There are cases where the debiased estimators dominated others (DGP 1.1 and DGP 1.3), and cases where the unadjusted estimator is the best (DGP 1.2 and DGP 2.2)\footnote{This is the artifact of the DGPs. Recall the variance formula for the difference-in-means estimator, for example from imbens2015causal.}. Note that DGP 1.3 and DGP 2.3 are constant effect models, in which the noninteracted OLS estimators are first order unbiased. However in DGP 1.3 we still observe a small, higher-order bias. \\ In terms of confidence interval coverages, observe in DGP 2.1, 2.2 and 2.3 that the original recentering intervals exhibit significant undercoverages. However, the procedure based on recomputing the OLS residuals with Satterthwaite adjustments works reasonably well. There is only one case DGP2.3, Non-Int, where the coverage is not very satisfactory. As shown in Table 6 and 8 in the appendix, the BC procedures do not add to the median confidence interval length significantly, although it tends to add to the average confidence interval length. The Satterthwaite could also add to the median (mean) confidence interval length: it typically increases the confidence interval length by at most 10 to 20 percent, compared with the Student-t adjustment.

Real Dataset

In this section we compare the performance of debiased estimators with that of standard OLS estimators on a real world dataset. We precisely follow lin2013agnostic's simulation setting. We generate our simulation from the experimental example of angrist2009incentives by simulating random assignments under the maintained hypothesis of no treatment effect. Because this setting assumes no effects, bias is expected to be negligible: freedman2008A notes that the leading term of the bias is greatest when treatment effects are heterogeneous. The simulation is thus not primarily meant to investigate bias, but rather the precision and coverage consequences of the use of our bias corrections in a real-world setting.\footnote{We thank Winston Lin for sharing the replication files.} \\ angrist2009incentives sought to measure the effects of support services and financial incentives on college students' academic achievement. The experiment randomly assigns eligible first-year undergraduate students into four groups. One treatment group was offered both support services and financial incentives. A second group was offered only support services and a third group only financial incentives. The control group was eligible only for standard university support services. As in lin2013agnostic, we only use the data for men in the services-and-incentives ($N=58$) and service-only ($N=99$) groups. The simulation datasets are generated assuming the treatment has no effect on any students. We replicate the experiments $10^7$ times, and each time randomly assign $58$ students to the services-and-incentives group and $99$ to the service-only group. The regression estimators estimate the treatment effects adjusting for high-school GPAs. The standard errors of the OLS estimators are estimated using the standard sandwich formulas.\\ Table 1 reports the simulation results from $10^7$ simulations. The first and second rows of the table show the means and standard deviations of the five estimators. All estimators are approximately unbiased after rounding and the variances of the debiased estimators are no larger than the standard estimators. The third row shows that all confidence intervals cover the true value of the ATE with approximately 95 percent probability. The fourth row reports the average length of the confidence intervals. On average, the intervals of regression-adjusted estimators are slightly narrower than that of the unadjusted estimator. (The width of the confidence intervals for the debiased estimators are mechanically identical to those for the standard estimators due to being constructed using the same SE estimators.)

Together, these results demonstrate that, in a real-world setting, our bias corrections can be effectively introduced without appreciably compromising the precision or coverage properties of regression adjusted estimators.

[4] [4]