Adjustment with Many Regressors Under Covariate-Adaptive Randomizations

Adjustments with Many Regressors under Covariate-Adaptive Randomizations

Introduction

This paper studies linear regression adjustments (RAs) for the estimation and inference of the average treatment effect (ATE) under covariate-adaptive randomizations (CARs) when there are many regressors. CARs have recently seen growing use in a wide variety of randomized experiments in economic research.\footnote{See, for example, CCFNT16,greaney2016,jakiela2016,burchardi2019,anderson2021, and etc.} Under CARs, units are first stratified using some baseline covariates, and within each stratum, the treatment status is assigned independently to achieve the balance between the numbers of treated and control units. Given data from CARs, researchers often estimate and infer various treatment effect parameters via regressions with strata dummies and baseline covariates as controls, a practice known as RAs. However, F08b,F081 showed that the usual OLS regression with covariates can actually decrease the precision of the ATE estimator. L13 found that, under complete randomization, for a linear regression with covariates to guarantee efficiency improvement upon the simple difference-in-means estimator (i.e., “no-harm”), it must include a full set of interactions between treatment status and covariates. Under more complicated CARs, the “no-harm” regressions are performed stratum by stratum in a fully saturated fashion, as pointed out by BCS18 and YYS22.

In many economics applications, the “no-harm” linear regressions are usually coupled with high dimensional regressors, as noted by duflo2018machinistas.\footnote{Examples of such experiments include bursztyn2019,banerjee2020lack,dhar2022reshaping, among others.} Moreover, when the regressors are sieve bases with a growing dimension, causal estimators with “no-harm” RAs can potentially achieve the semiparametric efficiency bound, as shown by JLTZ22 and BJRSZ23. However, CJN18_ET,CJN18 demonstrated that the usual heteroskedasticity consistent inference method for linear regressions leads to misleading inference when there are many regressors. CJM18 found similar problems for a class of non-linear settings. Under CARs, this issue of dimensionality is further exacerbated with the “no-harm” RAs stratum by stratum.

In this paper, we consider the “many regressors” regime--an asymptotic regime where the number of regressors can diverge at most as fast as the sample size and discover a new efficiency trade-off of using “no-harm” RAs due to the dimensionality of regressors. On one hand, RAs can improve the estimation efficiency by incorporating information from covariates that are not used in the randomization; on the other hand, they can degrade estimation efficiency due to their estimation errors, which are not asymptotically negligible when the number of regressors grows at the same rate of the sample size. This trade-off poses potentially serious challenges for inference. First, ignoring the cost of RAs in the asymptotic variance estimation may lead to substantial over-rejection under the null. Second, a consistent variance estimator that accounts for the cost of RAs can potentially be larger than that of the simple estimator without RAs, even asymptotically, if the cost outweighs the benefit. Therefore, using a “no-harm” RA that incurs such costs can actually be harmful.

To resolve these issues, we derive the joint asymptotic distribution for estimators both with and without RAs under the “many regressors” asymptotic regime and propose a consistent estimator of the asymptotic covariance matrix, which accounts for both the benefit and cost of RAs. We then use this matrix to construct a new ATE estimator by optimally linearly combining the two estimators. We show that (1) the Wald test based on this new optimal linear combination estimator achieves the exact asymptotic size under the null and (2) the optimal linear combination estimator is weakly more efficient than estimators both with and without RAs. Furthermore, we conduct a local asymptotic power analysis and establish that the Wald test based on the optimal linear combination estimator is asymptotically the uniformly most powerful (UMP) test against two-sided local alternatives over a class of unbiased tests that only depend on estimators with and without RAs. This implies our new estimator is optimal locally among not only linear but also nonlinear combination estimators. These results collectively address the challenge of the efficiency trade-off unearthed in the paper.

We further consider an alternative asymptotic regime where the dimension of regressors is fixed or moderate. Under mild regularity conditions, we show that the same optimal linear combination estimator with the same covariance matrix estimator achieves the same properties (1) and (2) as above. Importantly, these two results hold even when the RAs are not approximately correctly specified.

Relation to the literature. Our results build on the work of CJN18, who studied the inference of linear coefficients in OLS regression with many regressors and independent or cluster-independent observations. We extend their analysis in three new directions. First, we consider observations generated by CARs, which introduce cross-sectional dependence among treatment assignments and outcomes. Second, we estimate the ATE in a two-step procedure by linearly combining intercepts from multiple OLS regressions with many regressors. Because the first step intercept estimators are asymptotically dependent and the second step of the procedure introduces additional estimation errors, the asymptotic distribution of our ATE estimator is not a simple application of the result by CJN18. Instead, we develop a new distributional theory for our ATE estimator using techniques from BCS17, the Yurinskii's coupling,\footnote{We avoid the mistake mentioned in Remark 2.1 of CMU22 by formulating our result such that the Gaussian approximation error $\delta_n$ is explicitly defined as a vanishing sequence with respect to the sample size.} and an anti-concentration inequality from CCK14. Third, we propose a new optimal linear combination estimator and discuss its efficiency improvement under the many regressors framework, which is novel in the literature. CJM18 also considered a two-step procedure, and their knife-edge order of the number of regressors is the square root of the sample size. Our paper differs from theirs because our first-step intercept estimators enter the second step linearly, resulting in a knife-edge rate of the same order as the sample size.

Our covariance matrix estimator builds on cross-fit estimators for both the variance of regression coefficients and the “variance component” with many regressors, as proposed by J22 and KSS2020. To address the dependence among the first-step intercept estimators, we also extend KSS2020 by proposing a new estimator for the “covariance component” with many regressors, involving coefficients from two separate linear regressions. Other recent studies on inference in linear regression with many regressors include AS23 and MS23.

This paper also relates to other works on RAs in randomized experiments. JLTZ22, JPTZ22, LTM20, MTL20, WDTT16, and YYS22 considered RAs for various causal parameters when the covariates are either fixed dimensional or high-dimensional but sparse. In these settings, the estimation errors of RAs are asymptotically negligible. However, as argued by LM21, the sparsity condition (1) may not hold in social science applications as it is unclear why the large majority of control coefficients should be very nearly zero, (2) is not invariant to linear transformations of the controls, and (3) depends on some tuning parameter (e.g., the penalty term in Lasso), which complicates small sample interpretation of the resulting inference. WZ23 further demonstrated that both Lasso and debiased Lasso can exhibit significant omitted variable biases, even when the coefficient vector is sparse and the sample size is sufficiently large compared to the number of controls. In such cases, the “long regression” approach often outperforms both Lasso and debiased Lasso, utilizing the modern high-dimensional OLS-based inference methods from CJN18_ET, CJN18, CJM18. Our approach does not need sparsity conditions, thus avoiding these issues. Under a finite-population framework with complete random sampling, LD21 and CMO23 proposed bias correction for the ATE while allowing the number of regressors to grow at a rate slower than the sample size. More recently, LYW23 introduced a debiased regression-adjusted estimator under complete random sampling, permitting the number of regressors to grow at the same rate as the sample size. In contrast, we adopt a super-population framework with general CARs. A more detailed comparison with LD21, CMO23, and LYW23 is provided in Remark (ref) below. T18 further considered “optimal" stratification under CARs with a pilot experiment.

Recent studies (Bai22, Bai22, C23optimal, C23optimal, and BLST23, BLST23) have highlighted the optimality of “finely stratified” experiments. However, this paper focuses on “coarse” CARs for two main reasons. First, CARs, including simple random sampling and stratified block randomization, are widely used in empirical research as demonstrated by works such as CCFNT16, jakiela2016, burchardi2019, anderson2021. Based on a survey of selected development economists, B09 reported that about 40% of researchers have used such a design at some point. Second, “fine stratifying” continuous covariates at the experimental design stage may not be applicable in some scenarios because (1) sometimes, all relevant covariates are discrete, as seen in dupas2018, where randomization is stratified by gender, occupation, and bank branch-- all discrete variables; (2) researchers may need to analyze either past experiments or those conducted by others using CARs. Our paper thus aligns with the literature on RAs under CARs by taking the randomization scheme as given and aims to develop more efficient estimators using covariate information at the data analysis stage. Therefore, it complements the strategy of “fine stratification” at the experimental design stage because our method applies at the data analysis stage.

The rest of this paper is organized as follows. Section (ref) introduces the setup. Section (ref) identifies the efficiency trade-off in RAs with many regressors. In Section (ref), we present the optimal linear combination estimator for ATE and its asymptotic properties. Section (ref) discusses the properties of the optimal linear combination estimator with fixed or moderate number of regressors. Sections (ref) and (ref) provide simulations and an empirical application, respectively. Proofs and additional figures are given in the Online Appendix.

Notation. For any positive integer $m$, let $0_m$, $1_m$ and $I_m$ be the $m \times 1$ vector of zeros, $m \times 1$ vector of ones, and the $m \times m$ identity matrix, respectively. Let $||\cdot||_2$ and $||\cdot||_{op}$ denote the $\ell_2$-norm for a vector and operator norm for a matrix, respectively. For a symmetric and positive semi-definite matrix $\Upsilon$, $\lambda_{\max}(\Upsilon)$ denotes the maximum eigenvalue of $\Upsilon$. We define $W^{(n)}$ as the sample of $W$'s, i.e., $W^{(n)} = \{W_i\}_{i \in [n]}$, where $[n] = 1,\cdots,n$. We write $U \stackrel{d}{=} V$ for two random variables $U$ and $V$ if they share the same distribution.

Setup

Potential outcomes for treated and control groups are denoted by $Y(1)$ and $Y(0)$, respectively. Treatment status is denoted by $A$, with $A=1$ indicating treated and $A=0$ untreated. The stratum indicator is denoted by $S$, based on which the researcher implements the covariate-adaptive randomization. The support of $S$ is denoted by $\mathcal{S}$, a finite set. After randomization, the researcher can observe the data $\{Y_i,S_i,A_i,X_i\}_{i \in [n]}$ where $Y_i = Y_i(1)A_i + Y_i(0)(1-A_i)$ is the observed outcome, and $X_i$ contains covariates besides $S_i$ in the dataset. The support of $X$ is denoted as $\text{Supp}(X)$. In this paper, we allow $X_i$ and $S_i$ to be dependent. For $s \in S$, let $n_{s} = \sum_{i \in [n]}1\{S_i = s\}$, $n_{1,s} = \sum_{i \in [n]}A_i1\{S_i=s\}$, and $n_{0,s} = n_{s} - n_{1,s}$. Let $\aleph_{a,s} = \{i \in [n]: A_i = a,S_i=s\}$ denote the set of individuals in stratum $s$ with treatment status $a$ and $\aleph_{s} = \aleph_{1,s} \cup \aleph_{0,s} $. Our parameter of interest is the average treatment effect defined as $\tau = \mathbb{E}(Y(1)-Y(0))$.

We make the following assumptions on the data generating process (DGP) and the treatment assignment rule.

Efficiency Trade-off in RAs with Many Regressors

This section derives the joint asymptotic distribution of the fully saturated adjusted and unadjusted estimators under the assumption that the number of regressors grows no faster than the sample size. We then identify the efficiency trade-off in RAs based on the comparison of two asymptotic variances. The variance-covariance matrix will also be used to construct the optimal linear combination estimator in Section (ref).

The unadjusted estimator, denoted by $\hat \tau^{unadj}$, is the fully saturated estimator proposed by BCS18 and defined as

where for $s \in S$ and $a = 0,1$, $\hat \pi_s = n_{1,s}/n_s$ and $\hat \tau_{a,s}^{unadj} = \frac{1}{n_{a,s}} \sum_{i \in \aleph_{a,s}}Y_i$.

Following their lead, the fully saturated regression adjusted estimator considered in this paper is denoted by $\hat \tau^{adj}$ and defined as

where $\hat \tau_{a,s}$ is computed as the intercept in the OLS regression of $Y_i$ on $(1,\breve X_i)$ using all observations in $\aleph_{a,s}$, $\hat p_s = n_{s}/n$, and $\breve X_i = X_i - \overline{X}_{S_i}$ with $\overline{X}_s = \frac{1}{n_{s}}\sum_{i \in \aleph_s} X_i$. We further denote $\hat \beta_{a,s}$ as the OLS coefficient of $\breve X_i$ in the OLS regression. We note that $Y_i$ should not be demeaned in the above regression because our key parameter of interest is the intercept rather than the slope coefficient of $\breve X_i$. In fact, by the Frisch-Waugh-Lovell Theorem, the regression coefficient $\hat \beta_{a,s}$ is the same as if $X_i$ is demeaned by the cluster-treatment level mean ($\frac{1}{n_{a,s}}\sum_{i \in \aleph_{a,s}} X_i$) or not demeaned at all. Here, by demeaning $X_i$ with the cluster level mean, we effectively just change the intercept estimator $\hat \tau_{a,s}$ in the linear regression.

Asymptotic Properties

Denote the dimension of $X$ as $k$. In the following, we derive the asymptotic properties of $\hat \tau^{adj}$ and $\hat \tau^{unadj}$ jointly in the case that $k=k_n$ increases no faster than the sample size $n$.

To clearly state our assumptions, we need to introduce more notation. Let $\breve X_{\aleph_{a,s}}$ be an $n_{a,s} \times k_n$ matrix which is constructed by stacking $\breve X_i^\top$ for $i \in \aleph_{a,s}$. We define $A_{\aleph_{a,s}} \in \Re^{n_{a,s}}$ similarly. Then, let $P_{a,s} = \breve X_{\aleph_{a,s}} (\breve X_{\aleph_{a,s}}^\top\breve X_{\aleph_{a,s}})^{-1}\breve X_{\aleph_{a,s}}^\top$ be the projection matrix of $\breve X_{\aleph_{a,s}}$, and $M_{a,s} = I_{n_{a,s}} - P_{a,s}$. Let $ M_{a,s, i,j}$ and $ P_{a,s, i,j}$ be the $(i,j)$th entry of $M_{a,s}$ and $ P_{a,s}$, respectively. Define

Efficiency Trade-off

We note that the $(1,1)$ and $(2,2)$ elements of $\Sigma_{\mathcal U}$ represent the asymptotic variances of $\mathcal U^{adj}$ and $\mathcal U^{unadj}$, respectively, which we denote as $Var(\mathcal U^{adj})$ and $Var(\mathcal U^{unadj})$. We denote $Var(\mathcal V^{adj})$ and $Var(\mathcal V^{unadj})$ in the same manner. Following the argument in Remark (ref), we can show that $Var(\mathcal U^{adj}) = Var(\mathcal U^{unadj})$ if $k_n \log k_n = o(n)$. Moreover, it can be shown that $Var(\mathcal V^{adj}) \leq Var(\mathcal V^{unadj})$, which indicates the benefit of RA from incorporating information of $X_i$. Therefore, when $k_n \log k_n = o(n)$, $\hat \tau^{adj}$ is weakly more efficient than $\hat \tau^{unadj}$, a result that the previous literature has established under more stringent conditions.

However, when $k_n$ is of the same order of $n$, $Var(\mathcal U^{adj})$ can be larger than $Var(\mathcal U^{unadj})$. To illustrate this, let's continue with the Gaussian covariates example mentioned in Remark (ref). Suppose for $a=0,1$ and $s\in \mathcal S$, $\pi_s = 1/2$, $\kappa_{a,s} = \kappa$, and $\varepsilon_i(a)$ is homoskedastic in the sense that $\mathbb E (\varepsilon_i^2(a)|X_i,S_i=s) = \mathbb E (\varepsilon_i^2(a)|S_i=s)$, then we have

where $\zeta = \int_{\lambda_-}^{\lambda_+} \frac{\sqrt{(\lambda_+ - \lambda)(\lambda - \lambda_-)}}{2\pi\lambda^2} d\lambda$ and $\lambda_{\pm} = (1 \pm \sqrt{\kappa})^2$. Figure (ref) plots the values of the variance inflation factor (VIF) as a function of $\kappa$. We can see that the VIF is more than $12.5\%$, $25\%$, $50\%$, and $100\%$ if $\kappa$ is greater than $1/5$, $1/3$, $1/2$, $2/3$, respectively.

Consequently, an efficiency trade-off in RAs between $Var(\mathcal U^{adj})$ and $Var(\mathcal V^{adj})$ can be identified in Theorem (ref). This trade-off can be elucidated by decomposing the difference between the asymptotic variances of $\hat \tau^{adj}$ and $\hat \tau^{unadj}$ as

where the second term is always non-positive but the first term can be positive when $k_n$ is of the same order of the sample size. Under some circumstances, the first term may dominate the second, resulting in the estimator with the theoretically “no-harm” adjustment being even less efficient than the unadjusted estimator. For example, see Model 1 of Panel (b) in Figures (ref) and (ref). This defeats the purpose of using RAs.

This also implies that the asymptotic variance estimator of $\hat \tau^{adj}$ in the literature, e.g., the one proposed by YYS22, tends to under-estimate the true asymptotic variance when there are many regressors, as it neglects the VIF induced by the dimensionality. Therefore, when there are many regressors, the causal inference using such a variance estimator may over-reject under the null hypothesis, as illustrated in Section (ref).

Optimal Linear Combination Estimator

As shown above, ignoring the VIF can result in size distortion in causal inference. When we account for the VIF, however, the adjusted estimator may become even less efficient than the unadjusted one due to the estimation errors incurred in RAs with many regressors. In this section, we propose an optimal linear combination estimator, denoted by $\hat \tau^*$, designed to address both issues. We show that $\hat \tau^*$ is consistent, asymptotically normal, and weakly more efficient than both $\hat \tau^{adj}$ and $\hat \tau^{unadj}$. We also present the local asymptotic power analysis for the standard Wald test based on $\hat \tau^*$. Finally, we provide a consistent estimator $\hat \Sigma$ for $\Sigma = \Sigma_{\mathcal U} + \Sigma_{\mathcal V} + \Sigma_{\mathcal W}1_21_2^\top$ in Section (ref) below.

Asymptotic Properties of $\hat \tau^*$

We construct the optimal linear combination estimator as

which is a weighted average of $\hat \tau^{adj}$ and $\hat \tau^{unadj}$ with the weight $\hat \omega$ minimizing the asymptotic variance.

Denote $\hat \Sigma =

$ as a consistent estimator for $\Sigma \equiv \Sigma_{\mathcal U} + \Sigma_{\mathcal V} + \Sigma_{\mathcal W}1_21_2^\top =

$ and suppose $\hat \tau^{adj}$ and $\hat \tau^{unadj}$ are not asymptotically equivalent, i.e., $ \Sigma_{1,1} + \Sigma_{2,2} - 2 \Sigma_{1,2} > 0$. Then, the optimal weight is

and the corresponding estimator of the asymptotic variance for $\hat \tau^*$ is $\hat \sigma^2_* = (\hat w, 1- \hat w) \hat \Sigma (\hat w, 1- \hat w)^\top$.

Local Power

Suppose we want to test the null hypothesis $H_{0}:\tau = \tau_0$ versus the two-sided alternative $H_{1}:\tau \neq \tau_0$ with level $\alpha \in (0,1)$, and we are under the local alternative in the sense that $\tau = \tau_0 + \Delta/\sqrt{n}$. Then, the hypothesis can be rewritten as $H_0: \Delta = 0$ versus $H_1: \Delta \neq 0$. In addition, by Theorem (ref), we have

where $

= \Sigma^{-1/2}

$. Let us consider a limit experiment in which a researcher observes

knows the values of $(a,b)$, and wants to test $\Delta = 0$ versus the two-sided alternative. Then, by the factorization criterion, it is easy to see that $N^* = \frac{aN_1+bN_2}{\sqrt{a^2+b^2}} \sim N(\sqrt{a^2+b^2}\Delta,1)$ is a sufficient statistic for $\Delta$. Define the usual two-sided Wald test based on $\hat \tau^*$ as

where $\mathcal{C}_{\alpha}$ is the $(1-\alpha)$ quantile of the chi-squared distribution with one degree of freedom. Then, the following corollary implies $\mathbb W_n \rightsquigarrow 1\{ (N^*)^2 \geq \mathcal C_\alpha\}$. In addition, $1\{ (N^*)^2 \geq \mathcal C_\alpha\}$ is known as the uniformly most powerful (UMP) unbiased test (see, for example, LR06) in the sense that it is more powerful than any level $\alpha$ unbiased test that is a (potentially nonlinear) combination of $(N_1,N_2)$. Consequently, it means $\mathbb W_n$ is asymptotically the UMP test over a class of level $\alpha$ unbiased tests that only depend on

This implies $\hat \tau^*$ is optimal among not only linear but also nonlinear combination estimators.

To formalize the statement, we need to introduce the following notation. Suppose $(N_1,N_2)$ follow the joint distribution in (ref) and let

be a class of unbiased tests with level $\alpha$.

Covariance Matrix Estimator

We now propose a consistent estimator $\hat \Sigma$ for $\Sigma =\Sigma_{\mathcal U} + \Sigma_{\mathcal V} + \Sigma_{\mathcal W}1_21_2^\top$. We first note that $$\Sigma_{2,2} = \mathbb E\left[\frac{var(Y_i(1)|S_i)}{\pi_{S_i}} + \frac{var(Y_i(0)|S_i)}{1-\pi_{S_i}}\right] + \Sigma_{\mathcal W}$$ is the asymptotic variance of the unadjusted estimator and does not depend on $X_i$. Now let

This estimator is consistent as long as Assumption (ref) holds. See, for example, BCS18 for details. Therefore, it suffices to construct consistent estimators for the $(1,1)$ and $(1,2)$ elements of $ \Sigma_{\mathcal U}$ and $ \Sigma_{\mathcal V}$, and $ \Sigma_{\mathcal W}$.

Denote the residual of the regression and its leave-one-out version are denoted as $\hat \varepsilon_{a,s,i} = Y_i - \hat \tau_{a,s} - \breve{X}_i' \hat \beta_{a,s}$ and $\acute \varepsilon_{a,s,i} = \hat \varepsilon_{a,s,i}/M_{a,s,i,i}$ for $i \in \aleph_{a,s}$, respectively. Then, we construct our estimator $\hat \Sigma_{\mathcal U}$ following the cross-fit approach developed by J22. It is also possible to construct the estimator based on the sample splitting approach proposed by CJN18. The cross-fit estimator, as shown below, is asymptotically valid as long as $\kappa_{a,s}<1$, but is not guaranteed to be positive definite in finite sample. In contrast, the sample splitting estimator is guaranteed to be positive semidefinite but requires $\kappa_{a,s}<1/2$. When both of them are valid, they are asymptotically equivalent. Interested readers can refer to J22 for more discussions. For the rest of the paper, we focus on cross-fit estimators of the $(1,1)$ and $(1,2)$ elements of $\hat \Sigma_{\mathcal U}$, which is defined as

where

To define $\hat \Sigma_{\mathcal V}$, we note that $\mathbb{E}(Y_i(a)|X_i,S_i)-\mathbb{E}(Y_i(a)|S_i)$ can be approximated by $\breve{X}_i^\top \beta_{a,s}$ for $i \in \aleph_{a,s}$, where $\beta_{a,s}$ is defined in Assumption (ref)(iv). KSS2020 considered the estimation and inference for quadratic functions of the coefficient in one linear regression with many regressors. However, in our setting, $\Sigma_{\mathcal V}$ depends on multiple linear coefficients $\beta_{a,s}$ which are estimated from different strata of observations, and each regression model may only be approximately linear. Therefore, the estimator $\hat \Sigma_{\mathcal V}$ includes both the quadratic and cross-product terms of the estimated coefficients $\hat \beta_{a,s}$. Denote $\Gamma_{a,s} = \sum_{i \in \aleph_{a,s}}\breve{X}_i\breve{X}_i^\top$ and $\Gamma_{s} = \sum_{i \in \aleph_{s}}\breve{X}_i\breve{X}_i^\top$. Following the spirit of KSS2020, our estimator of the $(1,1)$ and $(1,2)$ elements of $\hat \Sigma_{\mathcal V}$ is defined as

where

The estimator of $\Sigma_{\mathcal W}$ is standard:

Indeed, stratum-specific RAs reduce the effective sample size, exacerbating the dimensionality issue. When $k_n$ exceeds the effective sample size $n_{a,s}$, researchers can employ Ridge-regularized linear regression for adjustment. Studying the statistical properties of the corresponding regression-adjusted ATE estimator is an interesting avenue for future research.

Fixed or Moderate Number of Regressors

In this section, we consider the CARs specified in Assumption (ref) with a fixed or moderately diverging dimension of $X_i$, i.e., $k_n$. We then introduce another set of regularity conditions to replace Assumptions (ref)--(ref). Under these new conditions, the Wald test based on the same estimator $\hat \tau^*$ - constructed using the same $\hat \tau^{adj}$, $\hat \tau^{unadj}$, and $\hat \Sigma$ as defined above - maintains exact asymptotic size under the null and is weakly more powerful than the Wald test based on the unadjusted estimator $\hat \tau^{unadj}$.

Particularly, in this dimension regime, there's no need to assume that the RAs are approximately correctly specified, as in Assumption (ref)(iv). So the results here are not merely special cases of those in Section (ref). Additionally, they are also novel to the literature of RAs with fixed number of regressors (e.g., YYS22). This is because the cross-fit estimator of the covariance matrix and the optimal linear combination estimator, motivated by our high-dimensional analysis, are new.

Simulations

We conduct a simulation study to evaluate the finite sample performance of our proposed inference method. The nominal level is $\alpha=5\%$ throughout the section.

Data Generating Processes

For $a\in\{0,1\}$, we generate potential outcomes according to the equation

where $\mu_{a}$, and $m_{a}\left(Z_{i}\right)$ and $\sigma_a(Z_i)$ are specified as follows. In each of the following specifications, $\{Z_i, \varepsilon_{i}(1), \varepsilon_{i}(0)\}$ are i.i.d. We have considered four models following the designs by CJN18. For the first three models, the dimension of $Z_i$ is denoted as $d_n$, and we let $X_i$ be either the full set or the first $k_n$ elements of $Z_i$. For the last model, $Z_i$ is of fixed dimension, and we let $X_i$ be either the full set or the first $k_n$ elements of polynomial series of $Z_i$. The polynomial series will be specified below. By varying the choice of $k_n$, we can illustrate the uniformity of our inference method.

For each model, strata are determined by dividing the support of $Z_{1i}$ into $|\mathcal{S}|$ intervals of equal length. Specifically, we have $S_i = \sum_{j = 1}^{|\mathcal{S}|} 1\{Z_{1i} \leq g_j\}, $ where $g_j=2j/|\mathcal{S}|-1$. For all strata, we set $\pi_s=1/2$. The treatment status is determined according to one of the following four CAR schemes:

Simulation Results

Rejection probabilities are computed using $10,000$ replications. We compare our methods with three existing methods recently introduced by BCS17, YYS22, and LTM20. We list these methods below.

Tables (ref) reports the size results for $n=400$ and $800$, under $H_{0}:\mu_1=\mu_0=0$ and $S=2$. In the first three models, we include all elements of $Z_i$ as regressors, i.e., $X_i=Z_i$. For the last model, we include all polynomial series of $Z_i$ as regressors. Therefore, the number of covariates used in each regression is approximately $\mathbf{40\%}$ of the effective sample size (i.e., $\kappa_{a,s} = 0.4$). We find that our methods $\hat{\tau}^{adj}$ and $\hat{\tau}^{\ast}$ have good size control under four models, improving as $n$ increases. In contrast, YYS has high rejection rates, over 10%, in all cases, even when $n$ is as large as 800. Our theory in Section (ref) shows this size distortion is caused by ignoring the estimation errors from many regressors. Furthermore, LTM has rejection probabilities over 7% in many cases. In particular, in Model 3 with correlated covariates, LTM has rejection rates over 9% in all randomization schemes, even when $n$ is 800. This is because in this model, the sparsity condition does not hold. Finally, BCS does not use any regressors, and thus, is immune to the size distortion due to the many regressors.

Table (ref) shows the power results under $H_{1}:\mu_1-\mu_0=0.2$. We find that, first, $\hat{\tau}^{\ast}$ always has higher power than $\hat{\tau}^{adj}$ and BCS, as our theory predicts. Second, YYS ignores the cost of RAs, resulting in a smaller variance estimator, and thus, higher rejection rates under both the null and alternative. Third, LTM has lower power than $\hat{\tau}^{\ast}$ in Models 1-3, where the sparsity condition for the regression coefficients fails. In fact, in this case, Lasso excludes many informative covariates in the adjustments. In model 4, LTM has better power than $\hat{\tau}^{\ast}$ since the sparsity condition is satisfied. By exploiting the sparsity, it is possible to show that the ATE estimator with the $\ell_1$ regularized adjustments achieves the semiparametric efficiency bound. However, even in this case that favors the Lasso regression adjustments, the power of our method based on $\hat \tau^*$ is very close to LTM.

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

Figures (ref) and (ref) show how the rejection rates of $\hat{\tau}^{adj}$, $\hat{\tau}^{\ast}$, YYS and LTM vary with the number of regressors under the null and alternative hypotheses. We consider Models 1-4 with $n=800$ and $n=400$, respectively. We set the number of regressors $k_n=0,2,\cdots,80$ ($k_n=0,1,\cdots,40)$ for $n=800$ ($n=400$), corresponding to $\kappa_{a,s} = 0,0.01,\cdots,0.4$. The x-axis is $\kappa_{a,s}$. The y-axis is the rejection probabilities. We display only the figures under the SBR scheme. Similar patterns are found under the other randomization schemes, which are shown in the Online Supplement.

Panels (a) of the figures show that our methods (i.e., $\hat \tau^{adj}$ and $\hat \tau^*$) have uniform size control over $\kappa_{a,s}$. Even when the dimensions of regressors vary from 0 to 80 (40) for $n=800$ ($n=400$), our methods have close-to-nominal rejection rates under $H_{0}$, consistent with our theories. In contrast, YYS has rejection probabilities that increase linearly with $\kappa_{a.s}$ in all models. Similarly, LTM displays this size distortion in Models 2 and 3, in which the sparsity condition fails.

Panels (b) of the figures show the power of the five methods under $H_{1}:\mu_1-\mu_0=0.2$. Several observations emerge. First, BCS's power, our benchmark, does not vary with $\kappa_{a,s}$ because it uses no regressors. Second, for all models, $\hat{\tau}^{\ast}$ is always more powerful than BCS and usually more powerful than $\hat \tau^{adj}$, as our theory predicts. Third, as the number of regressors increases, the misspecification error decreases, and $\hat{\tau}^{\ast}$ becomes more powerful. However, $\hat{\tau}^{adj}$'s power drops even below that of the unadjusted estimator as $k_n$ increases for Model 1. This reflects that many covariates do not always improve estimation accuracy, and the “no-harm” regression adjustments can in fact harm the estimation precision. Fourth, for Models 1-3, $\hat{\tau}^{\ast}$ is generally more powerful than LTM as the dimension of regressors increases. For Model 4, LTM slightly outperforms $\hat{\tau}^{\ast}$ when $\kappa_{a,s}$ is large for $n=800$ due to the sparsity of the model.

Practical Implication of Dimensionality

When the dimension of covariates, $k_n$, is relatively small compared to $n^{1/2}$, our inference method remains valid and "no-harm," even when the linear regression adjustment is misspecified. This also holds for conventional inference methods like that of YYS22.

When $k_n$ exceeds $n^{1/2}$ and the linear regression is approximately correctly specified (i.e., Assumption (ref)(iii)), our inference method remains valid and “no-harm”, even when $k_n$ is smaller but comparable to $n$, whereas the conventional method does not. Approximately correct specification is satisfied if the baseline covariates $Z_i$ are fixed-dimensional and continuous, with $X_i$ containing sieve basis functions of $Z_i$, or if the baseline covariates are categorical (discrete), with $X_i$ including fully saturated dummies for all categories.

When $k_n$ exceeds $n^{1/2}$ and the linear regression adjustment is not approximately correctly specified, the adjusted and linear combination estimators, $\hat{\tau}^{adj}$ and $\hat{\tau}^*$, may exhibit non-negligible biases. These biases can potentially be corrected using the methods developed by LD21 and CMO23, provided $k_n = o(n^{2/3})$. However, the statistical properties of these bias-corrected ATE estimators, and whether they remain 'no-harm' under CARs, have not yet been studied.

When $k_n$ exceeds $n^{2/3}$ and the linear regression adjustment is not approximately correctly specified, there is currently no known method for valid causal inference using regression-adjusted estimators under any randomization scheme, which is left for future research.

Empirical Application

Studying the effects of child health and nutrition on educational outcomes helps us understand the connection between health and economic development in less developed countries (glewwe2007, glewwe2007; dupas2017, dupas2017). To investigate the impact of an iron supplement program on students' schooling attainment, CCFNT16 conducted a randomized experiment in Peru with a CAR framework.

Their experiment involved 215 students from a rural secondary school in Peru during the 2009 academic year. Students were stratified by five grade levels, with each grade having students randomly assigned to one of three groups: two treatment groups and one control group. In the treatments, one group viewed a video where a well-known soccer player endorsed iron supplements, and the other group saw a video with a doctor's endorsement. The control group watched a video that did not mention iron supplements.

In this section, we focus on the impact of iron supplement program on the number of iron pills taken, whether a student is anemic at endline, and the student's cognitive ability. Throughout the analysis, students in the two treatment arms are grouped into the treatment group as in CCFNT16.

We use a fully saturated model, including two dummy covariates (gender and electricity at home) from the regressors in CCFNT16 and their interaction as the third covariate. This ensures that the model is correctly specified. Although the number of regressors is not large, it is still a significant proportion of the sample size given that the total sample size is not large either. $\kappa_{1,s}$ range from approximately 8% to 15%, and $\kappa_{0,s}$ range from about 16% to 30%.\footnote{We left out the strata where the number of observations in either treatment or control group is less than or equal to 4.}

Table (ref) presents the ATEs and the standard errors (in parentheses) estimated by different methods.\footnote{The description of these methods is similar to that in Section (ref).} We make three observations from the results. First, the variance estimates of the YYS method are smaller than those of $\hat{\tau}^{adj}$ in general. For example, the standard error for pills taken is approximately about 17.6% lower with YYS compared to $\hat{\tau}^{adj}$, despite identical ATE estimates. This indicates that YYS method may fail to fully capture estimation errors that arise when using many regressors.

Second, the standard errors of $\hat{\tau}^{\ast}$ are consistently lower than those of both $\hat{\tau}^{adj}$ and BCS for all outcomes, with some differences being substantial. For example, for cognitive ability, the standard error of $\hat{\tau}^{\ast}$ is about 8.2% smaller than that of BCS. This result is expected because $\hat{\tau}^{\ast}$ is a weighted average of $\hat \tau^{adj}$ and the unadjusted estimator (BCS), designed to minimize asymptotic variance. The standard errors of $\hat{\tau}^{\ast}$ are also generally lower than those of LTM, though by a smaller margin.

Third, consistent with CCFNT16, our results indicate that the iron supplement program increases the number of iron pills taken but has no significant effect on anemia or cognitive ability.

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

Conclusion

This paper shows that the estimation error of the “no-harm” regression adjustments can contaminate the ATE estimator and degrade estimation efficiency when the number of regressors is of the same order as the sample size. We then propose a new ATE estimator which is guaranteed to be weakly more efficient than both the estimators with and without RAs. We also propose a consistent estimator of its asymptotic variance and construct the corresponding Wald-statistic.