Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
121,328 characters · 19 sections · 143 citation commands
Approximate Residual Balancing: De-Biased Inference of Average Treatment Effects in High Dimensions
In order to identify causal effects in observational studies, practitioners may assume treatment assignments to be as good as random (or unconfounded) conditional on observed features of the units; see rosenbaum1983central and imbens2015causal for general discussions. Motivated by this setup, there is a large literature on how to adjust for differences in observed features between the treatment and control groups; some popular methods include regression, matching, propensity score weighting and subclassification, as well as doubly-robust combinations thereof abadie, heckman1998matching, hirano2003efficient, robins1994estimation, robins2, robins2008higher, rosenbaum2002observational, tan2010bounded, tsiatis2007semiparametric, van2006targeted.
In practice, researchers sometimes need to account for a substantial number of features to make this assumption of unconfoundedness plausible. For example, in an observational study of the effect of flu vaccines on hospitalization, we may be concerned that only controlling for differences in the age and sex distribution between controls and treated may not be sufficient to eliminate biases. In contrast, controlling for detailed medical histories and personal characteristics may make unconfoundedness more plausible. But the formal asymptotic theory in the earlier literature only considers the case where the sample size increases while the number of features remains fixed, and so approximations based on those results may not yield valid inferences in settings where the number of features is large, possibly even larger than the sample size.
There has been considerable recent interest in adapting methods from the earlier literature to high-dimensional settings. belloni2014inference,belloni2013program show that attempting to control for high-dimensional confounders using a regularized regression adjustment obtained via, e.g., the lasso, can result in substantial biases. belloni2014inference propose an augmented variable selection scheme to avoid this effect, while belloni2013program, chernozhukov2016double, farrell2015robust, and van2011targeted build on the work of robins1994estimation,robins2 and discuss how a doubly robust approach to average treatment effect estimation in high dimensions can also be used to compensate for the bias of regularized regression adjustments. Despite the breadth of research on the topic, all the above papers rely crucially on the existence of a consistent estimator of the propensity score, i.e., the conditional probability of receiving treatment given the features, in order to yield $\sqrt{n}$-consistent estimates of the average treatment effect in high dimensions.\footnote{Some of the above methods assume that the propensity scores can be consistently estimated using a sparse logistic model, while others allow for the use of more flexible modeling strategies following, e.g., mccaffrey2004propensity, van2007super, or westreich2010propensity.}
In this paper, we show that in settings where we are willing to entertain a sparse, well-specified linear model on the outcomes, efficient inference of average treatment effects in high-dimensions is possible under more general assumptions than suggested by the literature discussed above. Given linearity assumptions, we show that it is not necessary to consistently estimate treatment propensities; rather, it is enough to rely on de-biasing techniques building on recent developments in the high-dimensional inference literature javanmard2014confidence,javanmard2015biasing,van2014asymptotically,zhang2014confidence. In particular, in sparse linear models, we show that $\sqrt{n}$-consistent inference of average treatment effects is possible provided we simply require overlap, i.e., that the propensity score be uniformly bounded away from 0 and 1 for all values in the support of the pretreatment variables. We do not need to assume the existence of a consistent estimator of the propensity scores, or any form of sparsity on the propensity model.
The starting point behind both our method and the doubly robust methods of belloni2013program, chernozhukov2016double, farrell2015robust, van2011targeted, etc., is a recognition that high dimensional regression adjustments (such as the lasso) always shrink estimated effects, and that ignoring this shrinkage may result in prohibitively biased treatment effect estimates. The papers on doubly robust estimation then proceed to show that propensity-based adjustments can be used to compensate for this bias, provided we have a consistent propensity model that converges fast enough to the truth. Conceptually, this work builds on the result of rosenbaum1983central, who showed that controlling for the propensity score is sufficient to remove all biases associated with observed covariates, regardless of their functional form.
If we are willing to focus on high-dimensional linear models, however, it is possible to tighten the connection between the estimation strategy and the objective of estimating the average treatment effect and, in doing so, extend the number of settings where $\sqrt{n}$-consistent inference is possible. The key insight is that, in a linear model, propensity-based methods are attempting to solve a needlessly difficult task when they seek to eliminate biases of any functional form. Rather, in linear models, it is enough to correct for linear biases. In high dimensions, this can still be challenging; however, we find that it is possible to approximately correct for such biases whenever we assume overlap.
Concretely, we study the following two stage {approximate residual balancing} algorithm. First, we fit a regularized linear model for the outcome given the features separately in the two treatment groups. In the current paper we focus on the elastic net zou2005regularization and the lasso chen1998atomic,tibshirani1996regression for this component, and present formal results for the latter. In a second stage, we re-weight the first stage residuals using weights that approximately balance all the features between the treatment and control groups. Here we follow zubizarreta2015stable, and optimize the implied balance and variance provided by the weights, rather than the fit of the propensity score. Approximate balancing on all pretreatment variables (rather than exact balance on a subset of features, as in a regularized regression, or weighting using a regularized propensity model that may not be able to capture all relevant dimensions) allows us to guarantee that the bias arising from a potential failure to adjust for a large number of weak confounders can be bounded. Formally, this second step of re-weighting residuals using the weights proposed by zubizarreta2015stable is closely related to de-biasing corrections studied in the high-dimensional regression literature javanmard2014confidence,javanmard2015biasing,van2014asymptotically,zhang2014confidence; we comment further on this connection in Section (ref).
This approach also bears a close conceptual connection to work by chan2015globally, deville1992calibration, graham1, graham2, hainmueller, hellerstein1999imposing, imai2014covariate and zhao2016covariate, who fit propensity models to the data under a constraint that the resulting inverse-propensity weights exactly balance the covariate distributions between the treatment and control groups, and find that these methods out-perform propensity-based methods that do not impose balance. Such an approach, however, is only possible in low dimensions; in high dimensions where there are more covariates than samples, achieving exact balance is in general impossible. One of the key findings of this paper is that, in high dimensions, it is still often possible to achieve approximate balance under reasonable assumptions and that---when combined with a lasso regression adjustment---approximate balance suffices for eliminating bias due to regularization.
In our simulations, we find that three features of the algorithm are important: $(i)$ the direct covariance adjustment based on the outcome data with regularization to deal with the large number of features, $(ii)$ the weighting using the relation between the treatment and the features, and $(iii)$ the fact that the weights are based on direct measures of imbalance rather than on estimates of the propensity score. The finding that both weighting and regression adjustment are important is similar to conclusions drawn from the earlier literature on doubly robust estimation robins1997toward, where combining both techniques was shown to extend the set of problems where efficient treatment effect estimation is possible. The finding that weights designed to achieve balance perform better than weights based on the propensity score is consistent with findings in chan2015globally, graham1, graham2, hainmueller, imai2014covariate, and zubizarreta2015stable.
Our paper is structured as follows. First, in Section (ref), we motivate our two-stage procedure using a simple bound for its estimation error. Then, in Section (ref), we provide a formal analysis of our procedure under high-dimensional asymptotics, and we identify conditions under which approximate residual balancing is asymptotically Gaussian and allows for practical inference about the average treatment effect with dimension-free rates of convergence. Finally, in Section (ref), we conduct a simulation experiment, and find our method to perform well in a wide variety of settings relative to other proposals in the literature. A software implementation for R is available at \url{https://github.com/swager/balanceHD}.
Our goal is to estimate average treatment effects in the potential outcomes framework, or Rubin Causal Model rubin1974estimating,imbens2015causal. For each unit in a large population there is pair of (scalar) potential outcomes, $(\yin, \, \yie)$. Each unit is assigned to the treatment or not, with the treatment indicator denoted by $W_i\in\{0, \, 1\}$. Each unit is also characterized by a vector of covariates or features $X_i\in\mmr^p$, with $p$ potentially large, possibly larger than the sample size. For a random sample of size $n$ from this population, we observe the triple $(X_i, \, W_i, \, Y^\obs_i)$ for $i=1, \, \ldots, \, n$, where
is the realized outcome, equal to the potential outcome corresponding to the actual treatment received. The total number of treated units is equal to $\nt$ and the number of control units equals $\nc$. We frequently use the short-hand $\bx_\ccc$ and $\bx_\ttt$ for the feature matrices corresponding only to control or treated units respectively. We write the propensity score, i.e., the conditional probability of receiving the treatment given features, as $e(x)=\mathbb{P}[W_i=1|X_i=x]$ rosenbaum1983central. We focus primarily on the conditional average treatment effect for the treated {sample},
We note that the average treatment effect for the controls and the overall average effect can be handled similarly. Throughout the paper we assume {unconfoundedness}, i.e., that conditional on the pretreatment variables, treatment assignment is as good as random rosenbaum1983central; we also assume a linear model for the potential outcomes in both groups.
Here, we will only use the linear model for the control outcome because we focus on the average effect for the treated units, but if we were interested in the overall average effect we would need linearity in both groups. The linearity assumption is strong, but in high dimensions some strong structural assumption is in general needed for inference to be possible. Then, given linearity, we have
Estimating the first term is easy: $\hat{\mu}_{\rm t} =\overline{Y}_{\rm t}= \sum_{\cb{i : W_i = 1}} Y_i^\obs \, /\nt$ is unbiased for $\mu_{\rm t}$. In contrast, estimating $\mu_{\rm c}$ is a major challenge, especially in settings where $p$ is large, and it is the main focus of the paper.
We begin by reviewing two classical approaches to estimating $\mu_{\rm c}$, and thus also $\tau$, in the above linear model. The first is a weighting-based approach, which seeks to re-weight the control sample to make it look more like the treatment sample; the second is a regression-based approach, which seeks to adjust for differences in features between treated and control units by fitting an accurate model to the outcomes. Though neither approach alone performs well in a high-dimensional setting with a generic propensity score, we find that these two approaches can be fruitfully combined to obtain better estimators for $\tau$.
A first approach is to re-weight the control dataset using weights $\gamma_i$ to make the weighted covariate distribution mimic the covariate distribution in the treatment population. Given the weights we estimate $\hat{\mu}_{\rm c}$ as a weighted average \smash{$\hat{\mu}_{\rm c} = \sum_{\cb{i : W_i = 0}} \gamma_i \, Y_i^\obs$}. The standard way of selecting weights $\gamma_i$ uses the propensity score: \smash{$\gamma_i = {e(X_i)}/(1 - e(X_i)) \, / \, (\sum_{\cb{i : W_j = 0}} { e(X_j)}/(1 - e(X_j))).$} To implement these methods researchers typically substitute an estimate of the propensity score into this expression. Such {inverse-propensity weights} with non-parametric propensity score estimates have desirable asymptotic properties in settings with a small number of covariates hirano2003efficient. The finite-sample performance of methods based on inverse-propensity weighting can be poor, however, both in settings with limited overlap in covariate distributions and in settings with many covariates. A key difficulty is that estimating the treatment effect then involves dividing by $1 - \hat{e}(X_i)$, and so small inaccuracies in $\hat{e}(X_i)$ can have large effects, especially when $e(x)$ can be close to one; this problem is often quite severe in high dimensions.
As discussed in the introduction, if the control potential outcomes $Y_i(0)$ have a linear dependence on $X_i$, then using weights $\gamma_i$ that explicitly seek to balance the features $X_i$ is often advantageous deville1992calibration, chan2015globally, graham1, graham2, hainmueller, hellerstein1999imposing, imai2014covariate, zhao2016covariate, zubizarreta2015stable. This is a subtle but important improvement. The motivation behind this approach is that, in a linear model, the bias for estimators based on weighted averaging depends solely on \smash{$\overline{X}_{\rm t} - \sum_{\cb{i : W_i = 0}} \gamma_i \, X_i$}. Therefore getting the propensity model exactly right is less important than accurately matching the moments of \smash{$\overline{X}_{\rm t}$}. In high dimensions, however, exact balancing weights do not in general exist. When $p \gg \nc$, there will in general be no weights $\gamma_i$ for which $\overline{X}_{\rm t} - \sum_{\cb{i : W_i = 0}} \gamma_i \, X_i = 0$, and even in settings where $p <\nc$ but $p$ is large such estimators would not have good properties. zubizarreta2015stable extends the balancing weights approach to allow for weights that achieve approximate balance instead of exact balance; however, directly using his approach does not allow for $\sqrt{n}$-consistent estimation in a regime where $p$ is much larger than $n$.
A second approach is to compute an estimator $\hat{\beta}_{\rm c}$ for $\beta_{\rm c}$ using the $\nc$ control observations, and then estimate $\mu_{\rm c}$ as $\hat{\mu}_{\rm c} = \overline{X}_{\rm t} \cdot \hat{\beta}_{\rm c}$. In a low-dimensional regime with $p \ll \nc$, the ordinary least squares estimator for $\beta_{\rm c}$ is a natural choice, and yields an accurate and unbiased estimate of $\mu_{\rm c}$. In high dimensions, however, the problem is more delicate: Accurate unbiased estimation of the regression adjustment is in general impossible, and methods such as the lasso, ridge regression, or the elastic net may perform poorly when plugged in for $\beta_{\rm c}$, in particular when $\overline{X}_{\rm t}$ is far away from $\overline{X}_{\rm c}$. As stressed by belloni2014inference, the problem with plain lasso regression adjustments is that features with a substantial difference in average values between the two treatment arms can generate large biases even if the coefficients on these features in the outcome regression are small. Thus, a regularized regression that has been tuned to optimize goodness of fit on the outcome model is not appropriate whenever bias in the treatment effect estimate due to failing to control for potential confounders is of concern. To address this problem, belloni2014inference propose running least squares regression on the union of two sets of selected variables, one selected by a lasso regressing the outcome on the covariates, and the other selected by a lasso logistic regression for the treatment assignment. We note that estimating $\mu_{\rm c}$ by a regression adjustment $\hat{\mu}_{\rm c}=\overline{X}_{\rm t} \cdot \hat{\beta}_{\rm c}$, with $\hat{\beta}_{\rm c}$ estimated by ordinary least squares on selected variables, is implicitly equivalent to using a weighted averaging estimator with weights $\gamma$ chosen to balance the selected features robins2007comment. The belloni2014inference approach works well in settings where both the outcome regression and the treatment regression are at least approximately sparse. However, when the propensity is not sparse, we find that the performance of such double-selection methods is often poor.
Here we propose a new method combining weighting and regression adjustments to overcome the limitations of each method. In the first step of our method, we use a regularized linear model, e.g., the lasso or the elastic net, to obtain a pilot estimate of the treatment effect. In the second step, we do “approximate balancing” of the regression residuals to estimate treatment effects: We weight the residuals using weights that achieve approximate balance of the covariate distribution between treatment and control groups. This step compensates for the potential bias of the pilot estimator that arises due to confounders that may be weakly correlated with the outcome but are important due to their correlation with the treatment assignment. We find that the regression adjustment is effective at capturing strong effects; the weighting on the other hand is effective at capturing small effects. The combination leads to an effective and simple-to-implement estimator for average treatment effects with many features.
We focus on a meta-algorithm that first computes an estimate $\hat{\beta}_{\rm c}$ of $\beta_{\rm c}$ using the full sample of control units. This estimator may take a variety of forms, but typically it will involve some form of regularization to deal with the number of features. Second we compute weights $\gamma_i$ that balance the covariatees at least approximately, and apply these weights to the residuals cassel1976some,robins1994estimation:
In other words, we fit a model parametrized by $\beta_{\rm c}$ to capture some of the strong signals, and then use direct numerical re-balancing of the control data on the features to extract left-over signal from the residuals $Y_i^\obs - X_i \cdot \hat{\beta}_{\rm c}$. Ideally, we would hope for the first term to take care of any strong effects, while the re-balancing of the residuals can efficiently take care of the small spread-out effects. Our theory and experiments will verify that this is in fact the case.
A major advantage of the functional form in (ref) is that it yields a simple and powerful theoretical guarantee, as stated below. Recall that ${\bf X}_{\rm c}$ is the feature matrix for the control units. Consider the difference between $\hat{\mu}_{\rm c}$ and $\mu_{\rm c}$ for our proposed approach: $ \hat{\mu}_{\rm c} - \mu_{\rm c}= (\overline{X}_{\rm t} - {\bf X}_{\rm c}^\top \gamma) \cdot (\hat{\beta}_{\rm c} - \beta_{\rm c})+\gamma \cdot \varepsilon$, where $\varepsilon$ is the intrinsic noise $\varepsilon_i = Y_i(0) - X_i \cdot \beta_{\rm c}$. With only the regression adjustment and no weighting, the difference would be $\hat{\mu}_{{\rm c},{\rm reg}} - \mu_{\rm c}= (\overline{X}_{\rm t} - \overline{X}_{\rm c}) \cdot (\hat{\beta}_{\rm c} - \beta_{\rm c}) + \mathbf{1} \cdot \varepsilon/\nc$, and with only the weighting the difference would be $\hat{\mu}_{{\rm c},{\rm weight}} - \mu_{\rm c}=(\overline{X}_{\rm t} - {\bf X}_{\rm c}^\top \gamma) \cdot \beta_{\rm c} + \gamma \cdot \varepsilon$. Without any adjustment, just using the average outcome for the controls as an estimator for $\mu_{\rm c}$, the difference between the estimator for $\mu_{\rm c}$ and its actual value would be $\hat{\mu}_{{\rm c},{\rm no-adj}} - \mu_{\rm c}=(\overline{X}_{\rm t} - \overline{X}_{\rm c}) \cdot \beta_{\rm c}+\mathbf{1} \cdot \varepsilon/\nc$. The regression reduces the bias from $(\overline{X}_{\rm t} - \overline{X}_{\rm c}) \cdot \beta_{\rm c}$ to $(\overline{X}_{\rm t} - \overline{X}_{\rm c}) \cdot (\hat{\beta}_{\rm c} - \beta_{\rm c})$, which will be substantial reduction if the estimation error $(\hat{\beta}_{\rm c} - \beta_{\rm c})$ is small relative to $\beta_{\rm c}$. The weighting further reduces this to $(\overline{X}_{\rm t} - {\bf X}_{\rm c}^\top \gamma) \cdot (\hat{\beta}_{\rm c} - \beta_{\rm c})$, which may be helpful if there is a substantial difference between $\overline{X}_{\rm t}$ and $\overline{X}_{\rm c}$. This result, formalized below, shows the complimentary nature of the regression adjustment and the weighting. All proofs are given in the appendix.
This result decomposes the error of $\hat{\mu}_{\rm c}$ into two parts. The first is a bias term reflecting the dimension $p$ of the covariates; the second term is a {variance term} that does not depend on it. The upshot is that the bias term, which encodes the high-dimensional nature of the problem, involves a product of two factors that can both other be made reasonably small; we will focus on regimes where the first term should be expected to scale as $\oo(\sqrt{\log(p)/n})$, while the second term scales as $\oo(k\sqrt{\log(p)/n})$ where $k$ is the sparsity of the outcome model. If we are in a sparse enough regime (i.e., $k$ is small enough), Proposition (ref) implies that our procedure will be variance dominated; and, under these conditions, we also show that it is $\sqrt{n}$-consistent.
In order to exploit Proposition (ref), we need to make concrete choices for the weights $\gamma$ and the parameter estimates $\hat{\beta}_{\rm c}$. First, just like zubizarreta2015stable, we choose our weights $\gamma$ to directly optimize the bias and variance terms in Proposition (ref); the functional form of $\gamma$ is given in (ref), where $\zeta \in (0, \, 1)$ is a tuning parameter. We refer to them as approximately balancing weights since they seek to make the mean of the re-weighted control sample, namely ${\bf X}_{\rm c}^\top \gamma$, match the treated sample mean $\overline{X}_{\rm t}$ as closely as possible. The positivity constraint on $\gamma_i$ in (ref) aims to prevent the method from extrapolating too aggressively, while the upper bound is added for technical reasons discussed in Section (ref). Meanwhile, for estimating $\hat{\beta}_{\rm c}$, we simply need to use an estimator that achieves good enough risk bounds under $L_1$-risk. In our analysis, we focus on the lasso chen1998atomic,tibshirani1996regression; however, in experiments, we use the elastic net for additional stability zou2005regularization. Our complete algorithm is described in Procedure (ref).
Finally, although we do use this estimator in the present paper, we note that an analogous estimator for the average treatment effect $\EE{Y(1) - Y(0)}$ can also be constructed:
and $\gamma_{\rm c}$ is constructed similarly. This method can be analyzed using the same tools developed in this paper, and is available in our software package balanceHD. The conditions required for $\sqrt{n}$-consistent inference of $\tau_{ATE}$ using (ref) directly mirror the conditions (sparsity, overlap, etc.) listed in Section (ref) for inference about the average treatment effect on the treated.
The idea of combining weighted and regression-based approaches to treatment effect estimation has a long history in the causal inference literature. Given estimated propensity scores \smash{$\hat{e}(X_i)$}, cassel1976some and robins1994estimation propose using an augmented inverse-propensity weighted (AIPW) estimator,
the difference between this estimator and ours is that they obtain their weights $\gamma_i$ for (ref) via propensity-score modeling instead of quadratic programming. Estimators of this type have several desirable aspects: they are “doubly robust” in the sense that they are consistent whenever either the propensity fit $\hat{e}(\cdot)$ or the outcome fit $\hat{\beta}_{\rm c}$ is consistent, and they are asymptotically efficient in a semiparametric specification hahn1998role,hirano2003efficient,robins1,tsiatis2007semiparametric. However, a practical concern with this class of methods is that they may perform less well when \smash{$1 - \hat{e}(X_i)$} is close to 0 hirano2003efficient, schafer. Several higher-order refinements to the simple AIPW estimator (ref) have also been considered in the literature. In particular, schafer use the inverse-propensity weights \smash{$\hat{e}(X_i)/(1 - \hat{e}(X_i))$} as sample weights when estimating $\hat{\beta}_{\rm c}$, while scharfstein1999adjusting and van2006targeted consider adding these weights as features in the outcome model; see also robins2007comment and tan2010bounded for further discussion.
belloni2013program and farrell2015robust study the behavior of AIPW in high dimensions, and establish conditions under which they can reach efficiency when both the propensity function and the outcome model are consistently estimable. Intriguingly, in low dimensions, doubly robust methods are not necessary for achieving semiparametric efficiency. This rate can be achieved by either non-parametric inverse-propensity weighting or non-parametric regression adjustments on their own chen2008semiparametric,hirano2003efficient; doubly robust methods can then be used to relax the regularity conditions needed for efficiency farrell2015robust,robins2008higher. Conversely, in high-dimensions, both weighting and regression adjustments are required for $\sqrt{n}$-consistency belloni2013program,farrell2015robust,robins1997toward.
Although our estimator (ref) is cosmetically quite closely related to the AIPW estimator (ref), the motivation behind it is quite different. A common description of the AIPW estimator is that it tries to estimate two different nuisance components, i.e., the outcome model $\hat{\mu}_{\rm c}$ and the propensity model $\hat{e}$; it then achieves consistency if either of these components is itself estimated consistently, and efficiency if both components are estimated at fast enough rates. In contrast, our approximate residual balancing estimator bets on linearity twice: once in fitting the outcome model via the lasso, and once in de-biasing the lasso via balancing weights (ref).
By relying more heavily on linearity, we can considerably extend the set of problem under which $\sqrt{n}$-consistent is possible (assuming linearity in fact holds). As a concrete example, a simple analysis of AIPW estimation in high-dimensional linear models would start by assuming that the lasso is $o_P(n^{-1/4})$ consistent in root-mean squared error,\footnote{A more careful analysis of AIPW estimators can trade off the accuracy of the propensity and main effect models and, instead of requiring that both the propensity and outcome models can be estimated at $o_P(n^{-1/4})$ rates, only assumes that the product of the two rates be bounded as $o_P(n^{-1/2})$; see, e.g., farrell2015robust. In high dimensions, this amounts to assuming that the outcome and propensity models are both well specified and sparse, with respective sparsity levels $k_\beta$ and $k_e$ satisfying $k_\beta k_e \ll n / \log(p)^2$. AIPW can thus be preferable to ARB given sparse enough and well specified propensity models, with $k_e \ll \sqrt{n}/\log(p)$.} which can be attained via the lasso assuming a $k$-sparse true model with sparsity level $k \ll \sqrt{n}/\log(p)$; and this is, in fact, exactly the same condition we assume in Theorem (ref). Then, in addition to this requirement on the outcome model, AIPW estimators still need to posit the existence of an $o_P(n^{-1/4})$ consistent estimator of the treatment propensities, whereas we do not need to assume anything about the treatment assignment mechanism beyond overlap. The reason for this phenomenon is that the task of balancing (which is all that is needed to correct for the bias of the lasso in a linear model) is different from the task of estimating the propensity score---and is in fact often substantially easier.\footnote{The above distinctions are framed are in a situation where the statistician starts with a set of high-dimensional covariates, and needs to find a way to control for all of them at once. In this setting, linearity is a strong assumption, and so it is not surprising that making this assumption lets us considerably weaken requirements on other aspects of the problem. In other applications, however, the statistician may have started with low-dimensional data, but then created a high-dimensional design by listing series expansions of the original data, interactions, etc. In this setting, linearity is replaced with smoothness assumptions on the outcome model (since any smooth function can be well approximated using a large enough number of terms from an appropriately chosen series expansion). Here, variants of our procedure can be more directly compared with doubly robust methods and, in particular, the $\gamma_i$ in fact converge to the oracle inverse-propensity weights $e(X_i)/(1 - e(X_i))$; see hirshberg2017balancing and wang2017approximate for a discussion and further results.}
Our approximately balancing weights (ref) are inspired by the recent litearture on balancing weights chan2015globally, deville1992calibration, graham1, graham2, hainmueller, hellerstein1999imposing, hirr, imai2014covariate, zhao2016covariate. Most closely related, zubizarreta2015stable proposes estimating $\tau$ using the re-weighting formula as in Section (ref) with weights
where the tuning parameter is $t$; he calls these weights stable balancing weights. These weights are of course equivalent to ours, the only difference being that Zubizarreta bounds imbalance in constraint form whereas we do so in Lagrange form. The main conceptual difference between our setting and that of zubizarreta2015stable is that he considers problem settings where $p < \nc$, and then considers $t$ to be a practically small tuning parameter, e.g., $t = 0.1\sigma$ or $t = 0.001 \sigma$. However, in high dimensions, the optimization problem (ref) will not in general be feasible for small values of $t$; and in fact the bias term \smash{$\Norm{\overline{X}_{\rm t} - {\bf X}_{\rm c}^\top \gamma}_\infty$} becomes the dominant source of error in estimating $\tau$. We call our weights $\gamma$ “approximately” balancing in order to remind the reader of this fact. In settings where it is only possible to achieve approximate balance, weighting alone as considered in zubizarreta2015stable will not yield a $\sqrt{n}$-consistent estimate of the average treatment effect, and it is necessary to use a regularized regression adjustment as in (ref).
Similar estimators have been considered by graham1, graham2 and hainmueller in a setting where exact balancing is possible, with slightly different objection functions. For example, hainmueller uses $\sum_{i} \gamma_i\log(\gamma_i)$ instead of $\sum_{i} \gamma_i^2$, leading to
This estimator has attractive conceptual connections to logistic regression and maximum entropy estimation, and in a low dimensional setting where $W|X$ admits a well-specified logistic model the methods of graham1,graham2 and hainmueller are doubly robust zhao2017entropy; see also hirr, imbensel, and newey2004higher. In terms of our immediate concerns, however, the variance of $\htau$ depends on $\gamma$ through \smash{$\Norm{\gamma}_2^2$} and not \smash{$\sum \gamma_i\log\p{\gamma_i}$}, so our approximately balancing weights are more directly induced by our statistical objective than those defined in (ref).
Finally, in this paper, we have emphasized an asymptotic analysis point of view, where we evaluate estimators via their large sample accuracy. From this perspective, our estimator---which combines weighting with a regression adjustment as in (ref)---appears to largely dominate pure weighting estimators; in particular, in high dimensions, we achieve $\sqrt{n}$-consistency whereas pure weighting estimators do not. On the other hand, stressing practical concerns, rubin2008objective strongly argues that “designed based” inference leads to more credible conclusions in applications by better approximating randomized experiments. In our context, design based inference amounts to using a pure weighting estimator of the form $\sum \gamma_i Y_i$ where the $\gamma_i$ are chosen without looking at the $Y_i$. The methods considered by chan2015globally, graham1, hainmueller, zubizarreta2015stable, etc., all fit within this design-based paradigm, whereas ours does not.
As we have already emphasized, approximate residual balancing is a method that enables us to makes inferences about average treatment effects without needing to estimate treatment propensities as nuisance parameters; rather, we build on recent developments on inference in high-dimensional linear models cai2015confidence,javanmard2014confidence,javanmard2015biasing, ning2014general,van2014asymptotically,zhang2014confidence. Our main goal is to understand the asymptotics our estimates for $\mu_{\rm c} = \overline{X}_{\rm t} \cdot \beta_{\rm c}$. In the interest of generality, however, we begin by considering a broader problem, namely that of estimating generic contrasts $\xi \cdot \beta_{\rm c}$ in high-dimensional linear models. This detour via linear theory will help highlight the statistical phenomena that make approximate residual balancing work, and explain why---unlike the methods of belloni2013program, chernozhukov2016double or farrell2015robust---our method does not require consistent estimability of the treatment propensity function $e(x)$.
The problem of estimating sparse linear contrasts $\xi \cdot \beta_{\rm c}$ in high-dimensional regression problems has received considerable attention, including notable recent contributions by javanmard2014confidence,javanmard2015biasing, van2014asymptotically, and zhang2014confidence. These papers, however, exclusively consider the setting where $\xi$ is a sparse vector, and, in particular, focus on the case where $\xi$ is the $j$-th basis vector $e_j$, i.e., the target estimand is the $j$-th coordinate of $\beta_{\rm c}$. But, in our setting, the contrast vector $\overline{X}_{\rm t}$ defining our estimand $\mu_{\rm c} = \overline{X}_{\rm t} \cdot \beta_{\rm c}$ is random and thus generically dense; moreover, we are interested in applications where $m_{\rm t} = \mathbb{E}[\overline{X}_{\rm t}]$ itself may also be dense. Thus, a direct application of these method is not appropriate in our problem.\footnote{As a concrete example, Theorem 6 of javanmard2014confidence shows that their debiased estimator \smash{$\hbeta_{\text{c}}^{(\text{debiased})}$} satisfies \smash{$ \sqrt{n_{\rm c}}(\hbeta_{\text{c}}^{(\text{debiased})} - \beta_{\rm c}) = Z + \Delta$}, where $Z$ is a Gaussian random variable with desirable properties and \smash{$\lVert\Delta\rVert_\infty = o(1)$}. If we simply consider sparse contrasts of \smash{$\hat{\beta}_{\rm c}$}, then this error term $\Delta$ is negligible; however, in our setting, we would have a prohibitively large error term $\overline{X}_{\rm t} \cdot \Delta$ that may grow polynomially in $p$.}
An extension of this line of work to the problem of estimating dense, generic contrasts $\theta = \xi \cdot \beta_{\rm c}$ turns out to be closely related to our approximate residual balancing method for treatment effect estimation. To make this connection explicit, define the following estimator:
$\hat{\beta}_{\rm c}$ is a properly tuned sparse linear estimator, and $K$ is a tuning parameter discussed below. If we set $\xi \leftarrow \overline{X}_{\rm t}$, then this estimator is nothing but our treatment effect estimator from Procedure (ref).\footnote{Here, we phrased the imbalance constraint in constraint form rather than in Lagrange form; the reason for this is that, although there is a 1:1 mapping between these two settings, we found the former easier to work with formally whereas the latter appears to yield more consistent numerical performance. We also dropped the constraints $\sum \gamma_i = 1$ and $\gamma_i \geq 0$ for now, but will revisit them in Section (ref).} Conversely, in the classical parameter estimation setting with $\xi \leftarrow e_j$, the above procedure is algorithmically equivalent to the one proposed by javanmard2014confidence,javanmard2015biasing. Thus, the estimator (ref) can be thought of as an adaptation of the method of javanmard2014confidence,javanmard2015biasing that debiases \smash{$\hat{\beta}_{\rm c}$} specifically along the direction of interest $\xi$.
We begin our analysis in Section (ref) by considering a general version of (ref) under fairly strong “transformed independence design” generative assumptions on ${\bf X}_{\rm c}$. Although these assumptions may be too strong to be palatable in practical data analysis, this result lets us make a crisp conceptual link between approximate residual balancing and the debiased lasso. In particular, we find that (Theorem (ref)), \smash{$\htheta$} from (ref) is $\sqrt{n}$-consistent for $\theta$ provided \smash{$\xi^\top \Sigma_{\rm c}^{-1} \xi = \oo(1)$}, where $\Sigma_{\rm c}$ is the covariance of ${\bf X}_{\rm c}$. Interestingly, if \smash{$\Sigma_{\rm c} = I_{p \times p}$}, then in general \smash{$\xi^\top \Sigma_{\rm c}^{-1} \xi = \Norm{\xi}_2^2 = \oo(1)$} if and only if $\xi$ is very sparse, and so the classical de-biased lasso theory reviewed above is essentially sharp despite only considering the sparse-$\xi$ case cai2015confidence. On the other hand, whenever $\Sigma_{\rm c}$ has latent correlation structure, it is possible to have \smash{$\xi^\top \Sigma_{\rm c}^{-1} \xi = \oo(1)$} even when $\xi$ is dense and $\Norm{\xi}_2 \gg 1$, provided that $\xi$ is aligned with the large latent components of $\Sigma_{\rm c}$. We also note that, in the application to treatment effect estimation, \smash{$\overline{X}_{\rm t}^\top \Sigma_{\rm c}^{-1} \overline{X}_{\rm t}$} will in general be much larger than 1; however, in Corollary (ref) we show how to overcome this issue.
To our knowledge, this was the first result for $\sqrt{n}$-consistent inference about dense contrasts of $\beta_{\rm c}$ at the time we first circulated our manuscript. We note, however, simultaneous and independent work by zhu2016linear, who developed a promising method for testing hypotheses of the form $\xi \cdot \beta_{\rm c} = 0$ for potentially dense vectors $\xi$; their approach uses an orthogonal moments construction that relies on regressing $\xi \cdot X_i$ against a $p - 1$ dimensional design that captures the components of $X_i$ orthogonal to $\xi$.
Finally, in Section (ref), we revisit the specific problem of high-dimensional treatment effect estimation via approximate residual balancing under substantially weaker assumptions on the design matrix ${\bf X}_{\rm c}$: Rather than assuming a generative “transformed independence design” model, we simply require overlap and standard regularity conditions. The cost of relaxing our assumptions on ${\bf X}_{\rm c}$ is that we now get slightly looser performance guarantees; however, our asymptotic error rates are still in line with those we could get from doubly robust methods. We also discuss practical, heteroskedasticity-robust confidence intervals for $\tau$. Through our analysis, we assume that $\hat{\beta}_{\rm c}$ is obtained via the lasso; however, we could just as well consider, e.g., the square-root lasso belloni2011square, sorted $L_1$-penalized regression bogdan2015slope,su2016slope, or other methods with comparable $L_1$-risk bounds.
As we begin our analysis of $\htheta$ defined in (ref), it is first important to note that the optimization program (ref) is not always feasible. For example suppose that $p = 2\nc$, that ${\bf X}_{\rm c} = (I_{\nc \times \nc} \ I_{\nc \times \nc})$, and that $\xi$ consists of $n$ times “$1$” followed by $n$ times “$-1$”; then $\Norm{\xi - {\bf X}_{\rm c}^\top \gamma}_\infty \geq 1$ for any $\gamma \in \RR^{\nc}$, and the approximation error does not improve as $\nc$ and $p$ both get large. Thus, our first task is to identify a class of problems for which (ref) has a solution with high probability. The following lemma establishes such a result for random designs, in the case of vectors $\xi$ for which $\xi^\top \Sigma_{\rm c}^{-1} \xi$ is bounded; here $\Sigma_{\rm c} = \Var{X_i \cond W_i = 0}$ denotes the population variance of control features. We also rely on the following regularity condition, which will be needed for an application of the Hanson-Wright concentration bound for quadratic forms following rudelson2013hanson.
The above lemma is the key to our analysis of approximate residual balancing. Because, with high probability, the weights $\gamma^*$ from Lemma (ref) provide one feasible solution to the constraint in (ref), we conclude that, again with high probability, the actual weights we use for approximate residual balancing must satisfy $\Norm{\gamma}_2^2 \leq \Norm{\gamma^*}_2^2 \approx \nc^{-1} \xi^\top \Sigma_{\rm c}^{-1} \xi$. In order to turn this insight into a formal result, we need assumptions on both the sparsity of the signal and the covariance matrix $\Sigma_{\rm c}$.
The above sparsity requirement is quite strong. However, many analyses that seek to establish asymptotic normality in high dimensions rely on such an assumption. For example, javanmard2014confidence, van2014asymptotically, and zhang2014confidence all make this assumption when seeking to provide confidence intervals for individual components of $\beta_{\rm c}$; belloni2014inference use a similar assumption where they allow for additional non-zero components, but they assume that beyond the largest $k$ components with $k$ satisfying the same sparsity condition, the remaining non-zero elements of $\beta_{\rm c}$ are sufficiently small that they can be ignored, in what they refer to as approximate sparsity.\footnote{There are, of course, some exceptions to this assumption. In recent work, javanmard2015biasing show that inference of $\beta_{\rm c}$ is possible even when $k \ll n \, / \, \log(p)$ in a setting where $X$ is a random Gaussian matrix with either a known or extremely sparse population precision matrix; wager2016high show that lasso regression adjustments allow for efficient average treatment effect estimation in randomized trials even when $k \ll n \, / \, \log(p)$; while the method of zhu2016linear for estimating dense contrasts $\xi \cdot \beta_{\rm c}$ does not rely on sparsity of $\beta_{\rm c}$, and instead places assumptions on the joint distribution of $\xi \cdot X_i$ and the individual regressors. The point in common between these results is that they let us weaken the sparsity requirements at the expense of strengthening our assumptions about the $X$-distribution.}
Next, our analysis builds on well-known bounds on the estimation error of the lasso bickel2009simultaneous,hastie2015statistical that require ${\bf X}_{\rm c}$ to satisfy a form of the restricted eigenvalue condition. Below, we make a restricted eigenvalue assumption on \smash{$\Sigma_{\rm c}^{1/2}$}; then, we will use results from rudelson2013reconstruction to verify that this also implies a restricted eigenvalue condition on ${\bf X}_{\rm c}$.
The statement of Theorem (ref) highlights a connection between our debiased estimator (ref), and the ordinary least-squares (OLS) estimator. Under classical large-sample asymptotics with $n \gg p$, it is well known that the OLS estimator, \smash{$\htheta^{(OLS)} = \xi^\top ({\bf X}_{\rm c}^\top {\bf X}_{\rm c})^{-1} {\bf X}_{\rm c}^\top Y$}, satisfies
where $\gamma^*_i$ is as defined in Lemma (ref). By comparing this characterization to our result in Theorem (ref), it becomes apparent that our debiased estimator \smash{$\htheta$} has been able to recover the large-sample qualitative behavior of \smash{$\htheta^{(OLS)}$}, despite being in a high-dimensional $p \gg n$ regime. The connection between debiasing and OLS ought not appear too surprising. After all, under classical assumptions, \smash{$\htheta^{(OLS)}$} is known to be the minimum variance unbiased linear estimator for $\theta$; while the weights $\gamma$ in (ref) were explicitly chosen to minimize the variance of \smash{$\htheta$} subject to the estimator being nearly unbiased.
A downside of the above result is that our main goal is to estimate \smash{$\mu_{\rm c} = \overline{X}_{\rm t} \cdot \beta_{\rm c}$}, and this contrast-defining vector \smash{$\overline{X}_{\rm t}$} fails to satisfy the bound on \smash{$\overline{X}_{\rm t}^\top \Sigma^{-1} \overline{X}_{\rm t}$} assumed in Theorem (ref). In fact, because \smash{$\overline{X}_{\rm t}$} is random, this quantity will in general be on the order of $p/n$. In the result below, we show how to get around this problem under the weaker assumption that \smash{$m_{\rm t}^\top \Sigma_{\rm c}^{-1} m_{\rm t}$} is bounded; at a high level, the proof shows that the the stochasticity \smash{$\overline{X}_{\rm t}$} does not invalidate our previous result. We note that, because \smash{$\bY_{\rm t}$} is uncorrelated with \smash{$\hat{\mu}_{\rm c}$} conditionally on \smash{$\overline{X}_{\rm t}$}, the following result also immediately implies a central limit theorem for \smash{$\htau = \bY_{\rm t} - \hat{\mu}_{\rm c}$} where \smash{$\bY_{\rm t}$} is the average of the treated outcomes.
The asymptotic variance bound $m_{\rm t}^\top \Sigma_{\rm c}^{-1} m_{\rm t}$ is exactly the Mahalanobis distance between the mean treated and control subjects with respect to the covariance of the control sample. Thus, our result shows that we can achieve asymptotic inference about $\tau$ with a $1/\sqrt{n}$ rate of convergence, irrespective of the dimension of the features, subject only to a requirement on the Mahalanobis distance between the treated and control classes, and comparable sparsity assumptions on the $Y$-model as used by the rest of the high-dimensional inference literature, including belloni2014inference,belloni2013program, chernozhukov2016double and farrell2015robust. However, unlike this literature, we make no assumptions on the propensity model beyond overlap, and do not require it to be estimated consistently. In other words, by relying more heavily on linearity of the outcome function, we can considerably relax the assumptions required to get $\sqrt{n}$-consistent treatment effect estimation.
Our discussion so far, leading up to Corollary (ref), gives a characterization of when and why we should expect approximate residual balancing to work. However, from a practical perspective, the assumptions used in our derivation---in particular the transformed independence design assumption---were stronger than ones we may feel comfortable making in applications.
In this section, we propose an alternative analysis of approximate residual balancing based on overlap. Informally, overlap requires that each unit have a positive probability of receiving each of the treatment and control conditions, and thus that the treatment and control populations cannot be too dissimilar. Without overlap, estimation of average treatment effects relies fundamentally on extrapolation beyond the support of the features, and thus makes estimation inherently sensitive to functional form assumptions; and, for this reason, overlap has become a common assumption in the literature on causal inference from observational studies crump, imbens2015causal. For estimation of the average effect for the treated we in fact only need the propensity score to be bounded from above by $1-\eta$, but for estimation of the overall average effect we would require both the lower and upper bound on the propensity score. If we are willing to assume overlap, we can relax the transposed independence design assumption into much more routine regularity conditions on the design, as in Assumption (ref).
Following Lemma (ref), our analysis again proceeds by guessing a feasible solution to our optimization problem, and then using it to bound the variance of our estimator. Here, however, we use inverse-propensity weights as our guess: $\gamma^*_i \propto e(X_i)/(1 - e(X_i))$. Our proof hinges on showing that the actual weights we get from the optimization problem are at least as good as these inverse-propensity weights, and thus our method will be at most as variable as one that uses augmented inverse-propensity weighting (ref) with these oracle propensity weights.
The rate of convergence guaranteed by (ref) is the same as what we would get if we actually knew the true propensities and could use them for weighting robins1994estimation,robins2. Here, we achieve this rate although we have no guarantees that the true propensities $e(X_i)$ are consistently estimable. Finally, we note that, when the assumptions to Corollary (ref) hold, the bound (ref) is stronger than (ref); however, there exist designs where the bounds match wang2017approximate.
Finally, in applications, it is often of interest to have confidence intervals for $\mu_{\rm c}$ and $\tau$ rather than just point estimates; below, we propose such a construction. Much like the sandwich variance estimates for ordinary least squares regression, our proposed confidence intervals are heteroskedasticity robust even though the underlying point estimates were motivated using an argument written in terms of a homoskedastic sampling distribution.
In order to provide inference about $\tau$, we also need error bounds for $\hat{\mu}_{\rm t}$. Under sparsity assumptions comparable to those made for $\beta_{\rm c}$ in Theorem (ref), we can verify that
where $\hat{\beta}_{\rm t}$ is obtained using the lasso with $\lambda_n = 5 \nu\upsilon \sqrt{\log\p{p} / \nc}$. Moreover, $\hat{\mu}_{\rm c}$ and $\hat{\mu}_{\rm t}$ are independent conditionally on $X$ and $W$, thus implying that \smash{$ (\htau - \tau) \, /\, (\hV_{\rm c} + \hV_{\rm t})^{1/2} \Rightarrow \nn\p{0, \, 1}$}. This last expression is what we use for building confidence intervals for $\tau$.
Starting in 1986, California implemented the Greater Avenues to Independence (GAIN) program, with an aim to reduce dependence on welfare and promote work among disadvantaged households. The GAIN program provided its participants with a mix of educational resources such as English as a second language courses and vocational training, and job search assistance. This program is described in detail by hotz2006evaluating. In order to evaluate the effect of GAIN, the Manpower Development Research Corporation conducted a randomized study between 1988 and 1993, where a random subset of GAIN registrants were eligible to receive GAIN benefits immediately, whereas others were embargoed from the program until 1993 (after which point they were allowed to participate in the program). All experimental subjects were followed for a 3-year post-randomization period.
The randomization for the GAIN evaluation was conducted separately by county; following hotz2006evaluating, we consider data from Alameda, Los Angeles, Riverside and San Diego counties. As discussed in detail in hotz2006evaluating, the experimental conditions differed noticeably across counties, both in terms of the fraction of registrants eligible for GAIN, i.e., the treatment propensity, and in terms of the subjects participating in the experiment. For example, the GAIN programs in Riverside and San Diego counties sought to register all welfare cases in GAIN, while the programs in Alameda and Los Angeles counties focused on long-term welfare recipients.
The fact that the randomization of the GAIN evaluation was done at the county level rather than at the state level presents us with a natural opportunity to test our method, as follows. We seek to estimate the average treatment effect of GAIN on the treated; however, we hide the county information from our procedure, and instead try to compensate for sampling bias by controlling for a large amount of covariates. We used spline expansions of age and prior income, indicators for race, family status, etc., for a total of $p = 93$ covariates. Meanwhile, we can check our performance against a gold standard estimate of the average treatment effect that is stratified by county and thus guaranteed to be unbiased.\footnote{More formally, in our experiments, we set the gold standard using the county-stratified oracle estimator on bootstrap samples of the full $n = 19,170$ sample. We use bootstrap samples to correct for the correlation of estimators $\htau$ obtained using the full dataset and subsamples of it. We also note that, given this setup, the quantity we are using as our gold standard is not and estimate of $\tau$, i.e., the conditional average treatment effect on the treated sample, and should rather be thought of as an estimate of $\EE{\tau}$, i.e., the average treatment effect on the treated population. Since we are in a setting with a fairly weak signal, this should not make a noticeable difference in practice.}
We compare the behavior of different methods for estimating the average treatment effect on the treated using randomly drawn subsamples of the original data (the full dataset has $n = 19,170$). In addition to approximate residual balancing, we consider augmented inverse-propensity weighting (ref) and double selection following belloni2014inference as our baselines. We also show the behavior of an “oracle” procedure that gets to observe the hidden county information and then simply estimates treatment effects for each county separately, and the “naive” difference-in-means estimator that ignores the features $X$. In very small samples, the oracle procedure is not always well defined because some samples may result in counties where either everyone or no one is treated.
Figure (ref) compares the performance of the different methods. We see that approximate residual balancing and double selection both do well in terms of mean-squared error. Moreover, confidence intervals built via approximate residual balancing achieve effectively nominal coverage; double selection also gets reasonable coverage and improves with $n$. In contrast, augmented inverse-propensity weighting does not perform well here. The problem appears to be that estimating treatment propensities is quite difficult, and a cross-validated logistic elastic net often learns an effectively constant propensity model.
In addition to {\bf approximate residual balancing} as described in Procedure (ref), the methods we use as baselines are as follows: {\bf naive} difference-in-means estimation $\htau = \overline{Y}_{\rm t} - \overline{Y}_{\rm c}$ that ignores the covariate information $X$; the {\bf elastic net} zou2005regularization, or equivalently, Procedure (ref) with trivial weights $\gamma_i = 1/\nc$; {\bf approximate balancing}, or equivalently, Procedure (ref) with trivial parameter estimates \smash{$\hat{\beta}_{\rm c} = 0$} zubizarreta2015stable; {\bf inverse-propensity weighting}, as discussed in Section (ref), with propensity estimates $\hat{e}(X_i)$ obtained by elastic net logistic regression, with the propensity scores trimmed at $0.05$ and $0.95$; {\bf augmented inverse-propensity weighting}, which pairs elastic net regression adjustments with the above inverse-propensity weights (ref); the {\bf weighted elastic net}, motivated by schafer, that uses inverse-propensity weights as sample weights for the elastic net regression; {\bf targeted maximum likelihood estimation} (TMLE), which fine-tunes the elastic net regression estimates along the direction specified by the inverse-propensity weights van2006targeted, and {\bf ordinary least squares after model selection} where, in the spirit of belloni2014inference, we run lasso linear regression for $Y \cond X, \, W = 0$ and lasso logistic regression for $W \cond X$, and then compute the ordinary least squares estimate for $\tau$ on the union of the support of the three lasso problems.
Unless otherwise specified, all outcome and propensity models were fit using a (linear or logistic) elastic net. Whenever there is a “$\lambda$” regularization parameter to be selected, we use cross validation with the lambda.1se rule from the glmnet package friedman2010regularization. In belloni2014inference, the authors recommend selecting $\lambda$ using more sophisticated methods, such as the square-root lasso belloni2011square. However, in our simulations, our implementation of belloni2014inference still attains excellent performance in the regimes the method is designed to work in. Similarly, our confidence intervals for $\tau$ are built using a cross-validated choice of $\lambda$ instead of the fixed choice assumed by Corollary (ref). Our implementation of approximate residual balancing, as well as all the discussed baselines, is available in the R-package balanceHD.
We consider five different simulation settings. Our first setting is a {\bf two-cluster} layout, with data drawn as $Y_i = (C_i + Z_i) \cdot \beta + W_i + \varepsilon_i$. Here, $W_i = \text{Bernoulli}(0.5)$, $Z_i \sim \nn\p{0, \, I_{p \times p}}$, $\varepsilon_i \sim \nn\p{0, \, 1}$, and $C_i \in \RR^p$ is a cluster center that is one of $C_i \in \cb{0, \, \delta}$, such that $\PP{C_i = 0 \cond W_i = 0} = 0.8$ and $\PP{C_i = 0 \cond W_i = 1} = 0.2$. We consider two settings for the between-cluster vector $\delta$: a “dense” setting where $\delta = 4/\sqrt{n} \ \mathbf{1}$, and a “sparse” setting where $\delta_j = 40/\sqrt{n} \ 1\p{\cb{j = 1 \text{ modulo } 10}}$. Our second {\bf many-cluster} layout is closely related to the first, except now we have 20 cluster centers $C_i \in \cb{c_1, \, ..., \, c_{20}}$, where all the cluster centers are independently generated as $c_k \sim \nn\p{0, \, I_{p \times p}}$. To generate the data, we first draw $C_i$ uniformly at random from one of the 20 cluster centers and then set $W_i = 1$ wit probability $\eta$ for the first 10 clusters and $W_i = 1$ with probability $1 - \eta$ for the last 10 clusters; we tried both $\eta = 0.1$ and $\eta = 0.25$. We illustrate this simulation concept in Figure (ref). In both cases, we chose $\beta$ as one of
The signal strength was scaled such that $\Norm{\beta}_2 = 2$ in the two-cluster layout and $\Norm{\beta}_2 = 3$ in the many-cluster layout.
Our next two simulations are built using more traditional structural models. We first consider a {\bf sparse two-stage} setting closely inspired by an experiment of belloni2014inference. Here $X_i \sim \nn\p{0, \, \Sigma}$ with $\Sigma_{ij} = \rho^{\abs{i - j}}$, and $\theta_i = X_i \cdot \beta_W + \varepsilon_{i1}$. Then, $W_i \sim \text{Bernoulli}(1/(1 + e^{\theta_i}))$, and finally $Y_i = X_i \cdot \beta_Y + 0.5 \, W_i + \varepsilon_{i2}$ where $\varepsilon_{i1}$ and $\varepsilon_{i2}$ are independent standard Gaussian. Following belloni2014inference, we set the structure model as $(\beta_Y)_j \propto 1/j^2$ for $j = 1, \, ..., \, p$; for the propensity model, we consider both a “very sparse” propensity model $(\beta_W)_j \propto 1/j^2$, and also a “dense” propensity model $(\beta_W)_j \propto 1/\sqrt{j}$. A potential criticism of this simulation design is that the signal is perhaps unusually sparse (in the 4-th column of Table (ref), adjusting for differences in the two most important covariates removes 93% of the bias associated with all the covariates); moreover, we note that all the important coefficients of both $\beta_Y$ and $\beta_W$ are close to each other in terms of their indices; thus, the effect of using a correlated design may be mitigated. Thus, we also ran a {\bf moderately sparse two-stage} simulation, just like the above one, except we now used choices for $\beta_Y$ as in (ref), the only difference being that we shifted the indices of the betas, multiplying them by 23 mod $p$ (e.g., the harmonic setup now has $(\beta)_j \propto 1/[10 + (23\,(j - 1) \text{ mod } p)]$). Here, we drew the treatment assignments from a well-specified logistic model, $W_i \sim \text{Bernoulli}(1/(1 + \exp(-\sum_{j = 1}^{100} X_{ij} / 40)))$.
To test the robustness of all considered methods, we also ran a {\bf misspecified} simulation. Here, we first drew $X_i \sim \nn\p{0, \, I_{p \times p}}$, and defined latent parameters $\theta_i = \log(1 + \exp(-2 - 2 * (X_i)_1)) / 0.915$. We then drew $W_i \sim \text{Bernoulli}(1 - e^{-\theta_i})$, and finally $Y_i = (X_i)_1 + \cdots + (X_i)_{10} + \theta_i (2W_i - 1)/2 + \varepsilon_i$ with $\varepsilon_i \sim \nn\p{0, \, 1}$. We varied $n$ and $p$. This simulation setting, loosely inspired by the classic program evaluation dataset of lalonde1986evaluating, is illustrated in Figure (ref); note that the average treatment effect on the treated is much greater than the overall average treatment effect here.
In the first two experiments, for which we report results in Tables (ref) and (ref), the outcome model $Y|X$ is reasonably sparse, while the propensity model has overlap but is not in general sparse. In relative terms, this appears to hurt the double-selection method most. Meanwhile, in Table (ref), we find that the method of belloni2014inference has excellent performance---as expected---when both the propensity and outcome models are sparse. However, if we make the problem somewhat more difficult (Table (ref)), its performance decays substantially, and double selection lags both approximate residual balancing and propensity-based methods in its performance.
Generally, we find that the balancing performs substantially better than propensity score weighting, with or without direct covariate adjustment. We also find that combining direct covariate adjustment with weighting does better than weighting on its own, irrespective of whether the weighting is based on balance or on the propensity score. In these experiments, the weighted elastic net and TMLE also somewhat improve over AIPW.
Encouragingly, approximate residual balancing also does a good job in the misspecified setting from Table (ref). It appears that our stipulation that the approximately balancing weights (ref) must be non-negative (i.e., $\gamma_i \geq 0$) helps prevent our method from extrapolating too aggressively. Conversely, least squares with model selection does not perform well despite both the outcome and propensity models being sparse; apparently, it is more sensitive to the misspecification here. Perhaps the reason AIPW and TMLE do not do as well here is that there are very strong linear effects.
We evaluate coverage of confidence intervals in the “many-cluster” setting for different choices of $\beta$, $n$, and $p$; results are given in Table (ref). Coverage is generally better with more overlap ($\eta = 0.25$) rather than less ($\eta = 0.1$), and with sparser choices of $\beta$. Moreover, coverage rates appear to improve as $n$ increases, suggesting that we are in a regime where the asymptotics from Corollary (ref) are beginning to apply.
In this paper, we introduced approximate residual balancing as a method for unconfounded average treatment effect estimation in high-dimensional linear models. Under standard assumptions from the high-dimensional inference literature, our method allows for $\sqrt{n}$-consistent inference of the average treatment effect without any structural assumptions on the treatment assignment mechanism beyond overlap.
Widely used doubly robust methods, pioneered by robins1994estimation and studied further by several authors belloni2013program,farrell2015robust,schafer, scharfstein1999adjusting,robins2007comment,tan2010bounded,van2006targeted, approach this problem by trying to estimate two different nuisance components, the outcome model and the propensity model. These methods then achieve consistency if either nuisance component is itself consistently estimated, and achieve semiparametric efficiency if both components as estimated fast enough. In contrast, our method “bets” on linearity twice, both in fitting the lasso and in attempting to balance away its bias. In well specified linear models, this bet allows us to considerably extend the class of problems for with $\sqrt{n}$-consistent inference of average treatment effects is possible; thus, if a practitioner believes linearity to be a reasonable assumption in a given problem, our estimator may be a promising choice.
We end by mentioning two important questions left open by this paper. First, it would be important to develop a better understanding of how to choose the tuning parameter $\zeta$ in (ref) that trades off bias and variance in our balancing. Results from Theorems (ref) and (ref) provide some guidance on choosing $\zeta$ (via a constraint-form characterization); however, in our experiments, we achieved good performance by simply setting $\zeta = 1/2$ everywhere. The difficulty in choosing $\zeta$ is that we are trying to trade off an observable quantity (sampling variance) against an unobservable one (residual bias), and so cannot rely on simple methods like cross-validation that require unbiased estimates of the loss criterion we are trying to minimize. It would be of considerable interest to either devise a data-adaptive choice for $\zeta$, or understand why a fixed choice $\zeta = 1/2$ appears to achieve systematically good performance.
\sloppy{ It would also be interesting to extend our approach to generalized linear models, where there is a non-linear link function $\psi$ for which \smash{$\mathbb{E}[Y_i(c) \cond X_i = x] = \psi(x \cdot \beta_{\rm c})$}. In causal inference applications, this setting frequently arises when the outcomes $Y_i^\obs$ are binary, and we are willing to work with a logistic regression model. In this case, the first-order error component from using a pilot estimator \smash{$\hat{\beta}_{\rm c}$} for estimating $\mu_{\rm c}$ with a plug-in estimator \smash{$ n_{\rm c}^{-1}\sum_{\{W_i = 1\}}\psi(X_i \cdot \hat{\beta}_{\rm c})$} would be of the form \smash{$n_{\rm c}^{-1} \sum_{\{W_i = 1\}}\psi'(X_i \cdot \hat{\beta}_{\rm c})X_i(\beta_{\rm c} - \hat{\beta}_{\rm c})$}. An analogue to Proposition (ref) then suggests using an estimator}
However, due to space constraints, we leave a study of this estimator to further work.
{ {0.2pt plus 0.3ex} }