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.
70,538 characters · 11 sections · 73 citation commands
An alternative to synthetic control for models with many covariates under sparsity
The synthetic control method AbadieGardeazabal2003,AbadieDiamondHainmueller2010,AbadieDiamondHeinmueller2015 is one of the most recent additions to the empiricist's toolbox, gaining popularity not only in economics, but also in political science, medicine, etc. It provides a sound methodology in many settings where only long aggregate panel data is available to the researcher. The method has been specifically developed in a context where a single sizeable unit such as a country, a state or a city undergoes a large-scale policy change (referred to as the treatment or intervention hereafter), while only a moderate number of control units (the donor pool) is available to construct a counterfactual through a synthetic unit. This unit is defined as a convex combination of units from the donor pool that best resembles the treated unit before the intervention. Then the treatment effect is estimated from the difference in outcomes between the treated unit and its synthetic unit after the intervention takes place. In contexts such as those described above, the synthetic unit possesses several appealing properties abadie2019using. First, it does not lead to extrapolation outside the support of the data: because weights are non-negative, the counterfactual never takes a value outside of the convex hull defined by the donor pool. Second, one can assess simply its fit, making it easy to judge the quality of the counterfactual. Third, the synthetic unit is sparse: the number of control units receiving a non-zero weight is at most equal to the dimension of the matching variable plus one.
The method has still some limitations, in particular when applied to micro data, for which it was not initially intended. In such cases, the number of untreated units $n_0$ is typically greater than the dimension $p$ of variables $X$ used to construct the synthetic units. Then, as soon as the treated unit falls into the convex hull defined by the donor pool, the synthetic control solution is not uniquely defined abadie2017penalized. Second, and still related to the fact that the method was not developed for micro data, there is, to the best of our knowledge, no asymptotic theory available for synthetic control yet. This means in particular, that inference cannot be conducted in a standard way. A third issue is related to variable selection. The standard synthetic control method, as advocated in AbadieDiamondHainmueller2010, not only minimizes the norm $\left\lVert.\right\rVert_V$ -- defined for a vector $a$ of dimension $p$ and diagonal positive-definite matrix $V$, as $\left\lVerta\right\rVert_V = \sqrt{a^T V a}$ -- between the characteristics of the treated and those of its synthetic unit under constraints, but also employs a bi-level optimization program over the weighting matrix $V$ so as to obtain the best possible pre-treatment fit. Diagonal elements of $V$ are interpreted as a measure of the predicting power of each characteristics for the outcome AbadieDiamondHainmueller2010,abadie2019using. This approach has been criticized for being unstable and yielding unreproducible results, see in particular klossner2018comparative.
We consider an alternative to the synthetic control that addresses these issues. Specifically, we consider a parametric form for the synthetic control weights, $W_i = h(X_i^T\beta_0)$, where we estimate the unknown parameter $\beta_0$. This approach warrants the uniqueness of the solution in low-dimensional cases where $p<n_0$. With micro data, it may thus be seen as a particular solution of the synthetic control method. We show that the average treatment on the treated (ATT) parameter can be estimated with a two-step GMM estimator, where $\beta_0$ is computed in a first step so that the reweighted control group matches some features of the treated units. A key result is the double robustness of the estimator, as defined by RobinsBang2005. Specifically, we show that misspecifications in the synthetic control weights do not prevent valid inference if the outcome regression function is linear for the control group.
We then turn to the high-dimensional case where $p$ is large, possibly greater than $n_0$. This case actually corresponds to the initial set-up of the synthetic control method, and is therefore crucial to take into consideration. We depart from the synthetic control method by introducing an $\ell_1$ penalization term in the minimization program used to estimate $\beta_0$. We thus perform variable selection in a similar way as the Lasso, but differently from the synthetic control method, which relies on the aforementioned optimization over $V$ (leading to overweighting the variables that are good predictors of the outcome and underweighting the others).
We also study the asymptotic properties of our estimator. Building on double robustness, we construct an estimator that is immunized against first-step selection mistakes in the sense defined for example by ChernozhukovHansenSpindler2015b or DoubleML2018. This construction requires an extra step, which models the outcome regression function and provides a bias correction, a theme that has also been developed in ben2018augmented, abadie2017penalized and SyntheticDiffinDiff2019. We show that both in the low- and high-dimensional case, the estimator is consistent and asymptotically normal. Consequently, we develop inference based on asymptotic approximation, which can be used in place of permutation tests when randomization of the treatment is not warranted.
Apart from its close connection with the synthetic control method, the present paper is related to the literature on treatment effect evaluation through propensity score weighting and covariate balancing. Several efforts have been made to include balance between covariates as an explicit objective for estimation with or without relation to the propensity score Hainmueller2012, Graham2012. Our paper is in particular related to that of ImaiRatkovic2014, who integrate propensity score estimation and covariate balancing in the same framework. We extend their paper by considering the case of high-dimensional covariates. Note that the covariate balancing idea is related to the calibration estimation in survey sampling, see in particular deville1992calibration.
It also partakes in the econometric literature addressing variable selection, and more generally the use of machine learning tools, when estimating a treatment effect, especially but not exclusively in a high-dimensional framework. The lack of uniformity for inference after a selection step has been raised in a series of papers by LeebPotscher2005,LeebPotscher2008a,LeebPotscher2008b, echoing earlier papers by Leamer1983 who put into question the credibility of many empirical policy evaluation results. One recent solution proposed to circumvent this post-selection conundrum is the use of double-selection procedures BelloniChernozhukov2013,Farrell2015,ChernozhukovHansenSpindler2015,DoubleML2018. For example, BCH2012Treat highlight the dangers of selecting controls exclusively in their relation to the outcome and propose a three-step procedure that helps selecting more controls and guards against omitted variable biases much more than a simple “post-single-selection” estimator, as it is usually done by selecting covariates based on either their relation with the outcome or with the treatment variable, but rarely both. Farrell2015 extends the main approach of BCH2012Treat by allowing for heterogeneous treatment effects, proposing an estimator that is robust to either model selection mistakes in propensity scores or in outcome regression. However, Farrell2015 proves the root-n consistency of his estimator only when both models hold true (see Assumption 3 and Theorem 3 therein). Our paper is also related to the work of athey2018, who consider treatment effect estimation under the assumption of a linear conditional expectation for the outcome equation. As we do, they also estimate balancing weights, to correct for the bias arising in this high-dimensional setting, but because of their linearity assumption, they do not require to estimate a propensity score. Their method is then somewhat simpler than ours, but it does not enjoy the double-robustness property of ours.
Finally, a recent work by Wager2019 parallel to ours\footnote{Our results have been presented as early as 2016 at the North American and European Summer Meetings of the Econometric Society (see \url{ https://www.econometricsociety.org/sites/default/files/regions/program_ESEM2016.pdf}), and again in 2017 during the IAAE Meeting in Sapporo, see IAAE2017Sapporo.} suggests a Lasso-type procedure assuming logistic propensity score and linear specification for both treated and untreated items. Their method is similar but still slightly different from ours. The main focus in Wager2019 is to prove asymptotic normality under possibly weaker conditions on the sparsity levels of the parameters. Namely, they allow for sparsity up to $o(\sqrt{n} /\log(p))$ or, in some cases up to $o(n /\log(p))$ where $n$ is the sample size, when the eigenvalues of the population Gram matrix do not depend on $n$. Similar results can be developed for our method at the expense of additional technical effort. We have opted not to pursue in this direction since it only helps to include relatively non-sparse models that are not of interest for the applications we have in mind. More recently, the papers of NingImaiSida2018 and tan2020 written independently have been brought to our attention. tan2020 studies double-robustness of a similar estimator but assumes a logistic propensity score model while we leave room for other possibilities. The setting in NingImaiSida2018 is not restricted to the logistic model. However, that paper considers propensity score only as a mean to select the variables that should be fully balanced and that will enter in the final propensity score. This leads to a more complex estimation procedure. The authors of these two papers do not seem to be aware of the prior work of ours IAAE2017Sapporo.
The paper is organized as follows. Section (ref) introduces the set-up and the identification strategy behind our estimator. Section (ref) presents the estimator both in the low- and high-dimensional case and studies its asymptotic properties. Section (ref) examines the finite sample properties of our estimator through a Monte Carlo experiment. Section (ref) revisits lal86's dataset to compare our procedure with other high-dimensional econometric tools and the effect of the large-scale tobacco control program of AbadieDiamondHainmueller2010 for a comparison with synthetic control. Section (ref) concludes. All proofs are in the Appendix.
We are interested in the effect of a binary treatment, coded by $D=1$ for the treated and $D=0$ for the non-treated. We let $Y(0)$ and $Y(1)$ denote the potential outcome under no treatment and under the treatment, respectively. The observed outcome is then $Y=D Y(1) + (1-D) Y(0)$. We also observe a random vector $X \in \mathbb{R}^p$ of pre-treatment characteristics. The quantity of interest is the average treatment effect on the treated (ATT), defined as:
Here and in what follows, we assume that the random variables are such that all the considered expectations are finite. Since no individual is observed in both treatment states, identification of the counterfactual $\mathbb{E}[Y(0) \vert D=1]$ is achieved through the following two ubiquitous conditions.
Assumption (ref), a version of the usual common support condition, requires that there exist control units for any possible value of the covariates in the population. Since the ATT is the parameter of interest, we are never reconstructing a counterfactual for control units so $P[D=1 \vert X] > 0$ is not required. Assumption (ref) states that conditional on a set of observed confounding factors, the expected potential outcome under no treatment is the same for treated and control individuals. This assumption is a weaker form of the classical conditional independence assumption $(Y(0), Y(1)) \protect\mathpalette{\protect\independenT}{\perp} D \vert X$.
In policy evaluation settings, the counterfactual is usually identified and estimated as a weighted average of non-treated unit outcomes:
where $W$ is a random variable called the weight. Popular choices for the weight are the following:
In this paper, we propose another choice of weight $W$, which can be seen as a parametric alternative to the synthetic control. An advantage is that it is well-defined whether or not the number of untreated observations $n_0$ is greater than the dimension $p$ of $X$, whereas the synthetic control estimator is not uniquely defined when $n_0>p$. Formally, we look for weights $W$ that:
Satisfying a balancing condition means that
Up to a proportionality constant, this is equivalent to $\mathbb{E}[X \vert D=1] = \mathbb{E}[WX \vert D=0]$. In words, $W$ balances the first moment of the observed covariates between the treated and the control group. The definition of the observable covariates $X$ is left to the econometrician and can include transformation of the original covariates so as to match more features of their distribution. The idea behind such weights relies on the principle of “covariate balancing” as in, e.g., ImaiRatkovic2014. The following lemma shows that under Assumption (ref) weights satisfying the balancing condition always exist.
The proof of this lemma is straightforward by plugging the expression of $W$ in Equation (ref) and using the law of iterated expectations. Note that $W$ is not a unique solution of (ref). The linear regression weight $W= \mathbb{E}[DX^T] \mathbb{E}[(1-D) X X^T]^{-1}X$ also satisfies the balancing condition but it can be negative and its use is problematic in high-dimensional regime.
Lemma (ref) would suggest solving a binary choice model to obtain estimators of $P[D=1 \vert X]$ and of the weight $W$ as a first step, and then plugging $W$ in (ref) to estimate $\theta_0$. However, an inconsistent estimator of the propensity score leads to an inconsistent estimator of $\theta_0$ and does not guarantee that the corresponding weight will achieve covariate balancing. Finally, estimation of the propensity score can be problematic when there are very few treated units. For these reasons, we consider another approach where estimation is based directly on balancing equations:
An important advantage of this approach over the usual one based on the propensity score estimation through maximum likelihood is its double robustness RobinsBang2005. Indeed, let $W_1$ denote the weights identified by (ref) under a misspecified model on the propensity score. It turns out that if the balancing equations $\eqref{eq:calibr}$ hold for $W_1$ the estimated treatment effect will still be consistent provided that $\mathbb{E}[Y(0)|X]$ is linear in $X$. The formal result is given in Theorem (ref) below.
In (ref), the effect of $X$ is taken out from $Y$ in a linear fashion, while the effect of $X$ on $D$ is taken out by re-weighting the control group to obtain the same mean for $X$. Theorem (ref) shows that an estimator based on (ref) enjoys the double robustness property. Theorem (ref) is similar to the result of Kline2011 for the Oaxaca-Blinder estimator, which is obtained under the assumption that the propensity score follows specifically a log-logistic model in the propensity-score-well-specified case. Theorem (ref) is more general. It can be applied under parametric modeling of $W$ as well as in nonparametric settings.
In this paper, we consider a parametric model for $w(X)$. Namely, we assume that $P[D=1|X]=G(X^T\beta_0)$ for some unknown $\beta_0\in \mathbb{R}^p$ and some known strictly increasing cumulative distribution function $G$. Then $w(X)=h(X^T\beta_0)$ with $h = G/(1-G)$ and $\beta_0$ is identified by the balancing condition
Clearly, $h$ is a positive strictly increasing function, which implies that its primitive $H$ is strictly convex. A classical example is to take $G$ as the c.d.f. of the logistic distribution, in which case $h(u)= H(u)=\exp(u)$ for $u\in \mathbb{R}$. The strict convexity of $H$ implies that $\beta_0$ is the unique solution of a strictly convex program:
This program is well-defined, whether or not $P[D=1|X]=G(X^T\beta_0)$. Note also that definitions (ref) and (ref) are equivalent provided that $\mathbb{E}\left[ h(X^T\beta)\| X\|_2\right] <\infty$ for $\beta$ in a vicinity of $\beta_0$. Indeed, it follows from the dominated convergence theorem that, under this assumption and due to the fact that any convex function is locally Lipschitz, differentiation under the expectation sign is legitimate in (ref).
We are now ready to state the main identification theorem justifying the use of ATT estimation methods developed below. It is a straightforward corollary of Theorem (ref).
At this stage, the parameter $\mu$ in Equation (ref) does not play any role and can, for example, be zero. However, we will see below that in the high-dimensional regime, choosing $\mu$ carefully is crucial to obtain an “immunized” estimator of $\theta_0$ that enjoys the desirable asymptotic properties.
We now assume to have a sample $(D_i,X_i,Y_i)_{i=1...n}$ of i.i.d. random variables with the same distribution as $(D,X,Y)$.
Consider first an asymptotic regime where the dimension $p$ of the covariates is fixed, while the sample size $n$ tends to infinity. We call it the low-dimensional regime. Define an estimator of $\beta_0$ via the empirical counterpart of (ref):
Next, we plug $\hat\beta_{\rm ld}$ in the empirical counterpart of (ref) to obtain the following estimator of $\theta_0$:
Note that if $X$ includes the intercept, $\widetilde{\theta}_{\rm ld}$ satisfies the desirable property of location invariance, namely it does not change if we replace all $Y_i$ by $Y_i+c$, for any $c\in\mathbb{R}$.
Set $Z:=(D,X,Y)$, $Z_i := (D_i,X_i,Y_i)$ and introduce the function $$ g(Z,\theta,(\beta,\mu)):= [D - (1-D) h(X^T \beta) ] [Y - X^T \mu] - D \theta. $$ Then the estimator $\widetilde \theta_{\rm ld}$ satisfies
This estimator is a two-step GMM. It is consistent and asymptotically normal under mild regularity conditions, with asymptotic variance $\mathbb{E} \left[ g^2(Z, \theta_0, (\beta_0, \mu_0))\right]/\mathbb{E}(D)^2$, where
This can be shown by standard techniques NeweyMcFadden1994. Notice that since $h'(X^T \beta_0)>0$ the vector $\mu_0$ is the coefficient of the weighted population regression of $Y$ on $X$ for the control group. This observation is useful for the derivation of the “immunized” estimator in the high-dimensional case, to which we now turn.
We now consider that $p$ may grow with $n$, with possibly $p\gg n$. This can be of interest in several situations. First, in macroeconomic problems, $n$ is actually small, and $p$ may easily be of comparable size. For example, in the Tobacco control program application by AbadieDiamondHainmueller2010 the control group size is limited due to the fact that the observational unit is the state but many pre-treatment outcomes are included among the covariates. Section (ref) revisits this example. Second, researcher may want to consider a flexible form for the weights by including transformations of the covariates. For instance, one may want to interact categorical variables with other covariates or consider, e.g., different powers of continuous variables if one wants to allow for flexible non-linear effects. See Section (ref) for an application considering such transformations. Third, one may want not only to balance the first moments of the distribution of the covariates but also the second moments, the covariances, the third moments and so on to make the distribution more similar between the treated and the control group. In this case, high-dimensional settings seem to be of interest as well.
In high-dimensional regime, the GMM estimator in (ref) is, in general, not consistent. We therefore propose an alternative Lasso-type method by adding in (ref) an $\ell_1$ penalization term:
Here, $\lambda > 0$ is an overall penalty parameter set to dominate the noise in the gradient of the objective function and $\left\{ \psi_{j} \right\}_{j=1, \ldots,p}$ are covariate specific penalty loadings set as to grant good asymptotic properties. The penalty loadings can be adjusted using the algorithm presented in Appendix (ref).
This type of penalization offers several advantages. First, the program (ref) has almost surely a unique solution when the entries of $X$ have a continuous distribution, cf. Lemma 5 in Tibshirani2013, which cannot be granted for its non-penalized version (ref). Second, it yields a sparse solution in the sense that some entries of the vector of estimated coefficients are set exactly to zero if the penalty is large enough, which is not the case for estimators based on $\ell_2$ penalization. The $\ell_0$-penalized estimator shares the same sparsity property but is very costly to compute, whereas (ref) can be easily solved by computationally efficient methods, see, e.g., StatLearning.
The use of covariate specific penalty loadings goes back to BRT; the particular choice of penalty loadings that we consider below is inspired by BelloniChenChernozhukovHansen2012. A drawback of penalizing by the $\ell_1$-norm is that it induces a bias in estimation of the coefficients. But this is not an issue here since we are ultimately interested in estimating $\theta_0$ rather than $\beta_0$. The solution $\hat\beta$ of (ref) only plays the role of a pilot estimator.
The estimator $\hat\beta$ is consistent as $n$ tends to infinity under assumptions analogous to those used in BTW, BRT for the Lasso with quadratic loss, see Theorem (ref) below. As in the low-dimensional case (cf. Section (ref)), one is then tempted to consider the plug-in estimator for the ATT based on Equation (ref) with $\mu=0$:
We refer to this estimator as the naive plug-in estimator. However, as mentioned above, the Lasso estimator $\hat{\beta}$ of the nuisance parameter $\beta_0$ is not asymptotically unbiased. In high-dimensional regime where $p$ grows with $n$, naive plug-in estimators suffer from a regularization bias and may not be asymptotically normal with zero mean, as illustrated for example in BCH2012Treat,ChernozhukovHansenSpindler2015b,DoubleML2018. Therefore, following the general approach of ChernozhukovHansenSpindler2015b,DoubleML2018, we develop an immunized estimator that, at the first order, is insensitive to $\hat{\beta}$. We show that this estimator is asymptotically normal with mean zero and an asymptotic variance that does not depend on the properties of the pilot estimator $\hat{\beta}$. The idea is to choose parameter $\mu$ in (ref) such that the expected gradient of the estimating function $g(Z,\theta,(\beta,\mu))$ with respect to $\beta$ is zero when taken at $(\theta_0,\beta_0)$. This holds for $\mu=\mu_0$, where $\mu_0$ satisfies
Notice that if the corresponding matrix is invertible we get the low-dimensional solution (ref). Clearly, $\mu_0$ depends on unknown quantities and we need to estimate it. To this end, observe that Equation (ref) corresponds to the first-order condition of a weighted least-squares problem,\footnote{The assumptions under which we prove the results below guarantee that $\mu_0$ defined here is unique. Extension to the case of multiple solutions can be worked out as well. It is technically more involved but in our opinion does not add much to the understanding of the problem.} namely
Since $X$ is high-dimensional we cannot estimate $\mu_0$ via the empirical counterpart of (ref). Instead, we consider a Lasso-type estimator
Here, similarly to (ref), $\lambda' > 0$ is an overall penalty parameter set to dominate the noise in the gradient of the objective function and $\left\{ \psi'_{j} \right\}_{j=1, \ldots, p}$ are covariate-specific penalty loadings. Importantly, by estimating $\mu_0$ we do not introduce, at least asymptotically, an additional source of variability since by construction, the gradient of the moment condition (ref) (if we consider it as function of $\beta$ rather than of $\beta_0$) with respect to $(\beta,\mu)$ vanishes at point $(\beta_0,\mu_0)$.
Finally, the immunized ATT estimator is defined as
Intuitively, the immunized procedure corrects the naive plug-in estimator in the case where the balancing program has “missed” a covariate that is very important to predict the outcome: $$\hat{\theta}= \tilde \theta - \left[\frac{1}{n_1} \sum_{i: D_i = 1}^n X_i -\frac{1}{n_1} \sum_{i: D_i = 0}^n h(X_i^T\hat \beta)X_i \right]^T \hat \mu, $$ where $n_1$ is the number of treated observations. This has a flavor of Frish-Waugh-Lowell partialling-out procedure for model selection as observed by BCH2012Treat and further developed in ChernozhukovHansenSpindler2015b.
To summarize, the estimation procedure in high-dimensional regime consists of the three following steps. Each step is computationally simple as it needs at most to minimize a convex function:
The current framework poses several challenges to achieving asymptotically valid inference. First, $X$ can be high-dimensional since we allow for $p \gg n$ provided that sparsity conditions are met (see Assumption (ref) below). Second, the ATT estimation is affected by the estimation of the nuisance parameters $(\beta_0,\mu_0)$ and we wish to neutralize their influence. Finally, the $\ell_1$-penalized estimators we use for $\beta_0$ and $\mu_0$ are not conventional. The estimator of $\beta_0$ relies on a non-standard loss function and, to our knowledge, the properties of $\hat\beta$ that we need are not available in the literature, cf., e.g., vandegeer and the references therein. The estimator of $\mu_0$ is close to the usual Lasso except for the weights that depend on $\hat \beta$. In general, discrepancy in the weights can induce an extra bias. Thus, it is not granted that such an estimator achieves properties close to the Lasso. We show below that it holds true under our assumptions. Our proof techniques may be of interest for other problems of similar type.
Let $\eta=(\beta,\mu)$ denote the vector of two nuisance parameters and recall that $Z = (D,X,Y)$. In what follows, we write for brevity $g(Z,\theta,\eta)$ instead of $g(Z,\theta,(\beta,\mu))$. In particular, for the value $\eta_0:=(\beta_0,\mu_0)$ we have
Hereafter, the notation $a \lesssim b$ means that $a \leq cb$ for some constant $c>0$ independent of the sample size $n$. We denote by $\Phi$ and $\Phi^{-1}$ the cumulative distribution function and the quantile function of a standard normal random variable, respectively. We use the symbol $\mathbb{E}_n(\cdot)$ to denote the empirical average, that is $\mathbb{E}_n(a) = n^{-1} \sum_{i=1}^n a_i$ for $a=(a_1,\dots,a_n)$. Finally, for a vector $\delta=(\delta_1,\dots,\delta_p) \in \mathbb{R}^p$ and a subset $S \subseteq \{1,\dots,p\}$ we consider the restricted vector $\delta_S=(\delta_j\mathbb{I}(j\in S))_{j=1}^p$, where $\mathbb{I}(\cdot)$ denotes the indicator function, and we set $\left\lVert\delta\right\rVert_0 := \text{Card}\left\{1\le j \le p: \delta_j \neq 0 \right\}$, $\left\lVert\delta\right\rVert_1 := \sum_{j=1}^p \vert \delta_j \vert$, $\left\lVert\delta\right\rVert_2 := \sqrt{ \sum_{j=1}^p \delta_j^2}$ and $\left\lVert\delta\right\rVert_\infty := \underset{j=1,...,p}{\max} \vert \delta_j \vert$.
We now state the assumptions used to prove the asymptotic results.
We also need some conditions on the distribution of data. The random vectors $Z_i=(D_i,X_i,Y_i)$ are assumed to be i.i.d. copies of $Z=(D,X,Y)$ with $D\in \{0,1\}$, $X\in \mathbb{R}^p$ and $Y\in \mathbb{R}$. Throughout the paper, we assume that $Z$ depends on $n$, so that in fact we deal with a triangular array of random vectors. This dependence on $n$ is needed for rigorously stating asymptotical results since we consider the setting where the dimension $p=p(n)$ is a function of $n$. Thus, in what follows $Z$ is indexed by $n$ but for brevity we typically suppress this dependence in the notation. On the other hand, all constants denoted by $c$ (with various indices) and $K$ appearing below are independent of $n$.
Note that in view of Assumption (ref), $\kappa_{\Sigma}$ is uniformly bounded: $\kappa_{\Sigma}\le e_1^T \Sigma e_1\le K^2$.
Finally, we define the penalty loadings for estimation of nuisance parameters. The gradients of the estimating function with respect to the nuisance parameters are
For each $i=1,\dots,n$, we define the random vector $U_i \in \mathbb{R}^{2p}$ with entries corresponding to these gradients: \[ U_{i,j} := \left\{
\right. \] where $X_{i,j}$ is the $j$th entry of $X_i$.
Here, the upper bounds on the loadings are $$\psi_{j, \max}=\sqrt{\frac{1}{n} \sum_{i=1}^n \max_{|u|\le K}\left[(1-D_i) h(u) - D_i\right]^2 X_{i,j}^2} ,$$ $$ \psi'_{j,\max} = \sqrt{\frac{1}{n} \sum_{i=1}^n (1-D_i)\max_{|u|\le K} \big(|h'(u)|^2 \left[Y_i - u\right]^2\big) X_{i,j}^2} $$ representing feasible majorants for $\sqrt{n^{-1}\sum_{i=1}^n U_{i,j}^2}$ under our assumptions.
The values $\sqrt{n^{-1}\sum_{i=1}^n U_{i,j}^2}$ depend on the unknown parameters. Thus, we cannot choose $\psi_{j}, \psi'_{j}$ equal to these values but we can take them equal to the upper bounds $\psi_{j}=\psi_{j, \max}$ and $\psi'_{j}=\psi'_{j, \max}$. A more flexible iterative approach of choosing feasible loadings is discussed in Appendix (ref).\\
The following theorem constitutes the main asymptotic result of the paper.
The proof is given in Appendix (ref). An important point underlying the root-$n$ convergence and asymptotic normality of $\hat\theta$ is the fact that the expected gradient of $g$ with respect to $\eta$ is zero at $\eta_0$. In granting this property, we follow the general methodology of estimation in the presence of high-dimensional nuisance parameters developed in BelloniChenChernozhukovHansen2012, BCH2012Treat, BelloniChernozhukovFenanderValHansen2014, ChernozhukovHansenSpindler2015b among other papers by the same authors. The second important ingredient of the proof is to ensure that the estimator $\hat\eta$ converges fast enough to the nuisance parameter $\eta_0$. Its rate of convergence is given in the following theorem.
The proof is given in Appendix (ref). Theorem (ref), together with the fact that $\hat\theta$ is immunized, implies that $\hat\beta$ and $\hat\mu$ have no asymptotic effect on $\hat\theta$. In related work, similar conditions to (ref) and (ref) have been imposed rather than established. We refer in particular to (36) of ChernozhukovHansenSpindler2015b and Assumption 3-(ii) in Farrell2015. Note, on the other hand, that contrary to Farrell2015 (2015, see Assumption 3-(i)), we do not require that both $x\mapsto G(x'\hat\beta)$ and $x\mapsto x'\hat\mu$ are consistent for $x\mapsto P[D=1|X=x]$ and $x\mapsto \mathbb{E}[Y(0)|X=x]$. Theorem (ref) only requires (ref) to hold. By Theorem (ref) above, this is the case if either $P[D=1|X]=G(X^T\beta_0)$ or $\mathbb{E}[Y(0)|X]=X^T\mu_0$, in which case either $x\mapsto G(x'\hat\beta)$ or $x\mapsto x'\hat\mu$ is consistent, but not necessarily both.
The aim of this experiment is two-fold: illustrate the better properties of the immunized estimator over the naive plug-in, and compare it with other estimators. In particular, we compare it with a similar estimator proposed by Farrell2015. We consider the following DGP. The $p$ covariates are distributed as $X \sim \mathcal{N}(0,\Sigma)$, where the $(i,j)$th element of the covariance matrix satisfies $\Sigma_{i,j}=0.5^{\vert i-j \vert}$. The treatment equation follows a logit model, $P(D=1|X) =\Lambda\left( X^T \gamma_0 \right)$ with $\Lambda(u) = 1/(1+\exp(-u))$. The potential outcomes satisfy $$Y(0) = \exp(X^T\mu_0) + \varepsilon, \quad Y(1) = Y(0) + \zeta_0 X^T\gamma_0,$$ where $\zeta_0$ is a constant, $\varepsilon \sim \mathcal{N}(0,1)$ and $\varepsilon$ is independent of $(D,X)$. We assume the following form for the $j$th entry of $\gamma_0$ and $\mu_0$:
where $\rho_\gamma$ and $\rho_\mu$ are positive constants. We are thus in an strictly sparse setting for both equations where only ten covariates play a role in the treatment assignment and twenty in the outcome. Figure (ref) depicts the precise pattern of the corresponding coefficients for $p=30$. The constants $\rho_\gamma$ and $\rho_\mu$ fix the signal-to-noise ratio. More precisely, $\rho_\gamma$ is set so that $R^2=0.3$ in the latent model for $D$, and $\rho_\mu$ is set so that $R^2=0.8$ in the model for $Y(0)$. Finally, we let $\zeta_0= [V\left( \exp(X^T\mu_0) \right)/5V\left(Y(0)\right)]^{1/2}$. This impies that the variance of the individual treatment effect $\zeta_0 X'\gamma_0$ is one fifth of the variance of $Y(0)$. In this set-up, the ATT satisfies $E[Y(1) -Y(0) \vert D=1] = \zeta_0E[Z\Lambda(Z)]/E[\Lambda(Z)]$, with $Z \sim \mathcal{N}(0,\gamma_0^T \Sigma \gamma_0)$. We compute it using Monte Carlo simulations.
We consider several estimators of the ATT. The first is the naive plug-in estimator defined in (ref). Next, we consider our proposed estimator defined in (ref), with $H(x)=h(x)=\exp(x)$. We also consider the estimator proposed by Farrell2015. This estimator is also defined by (ref), but $\hat{\beta}$ and $\hat{\mu}$ therein are obtained by a Logit Lasso and an unweighted Lasso regression, respectively. For these three estimators, the penalty loadings for the first-step estimators are set as in Appendix (ref). The last estimator, called the oracle hereafter, is our low-dimensional estimator defined by (ref), where in the first step we only include the ten covariates affecting the treatment. For all estimators, we construct 95% confidence estimator on the ATT using the normal approximation and asymptotic variance estimators. We estimate the asymptotic variance of the naive plug-in estimator making as if we were in a low-dimensional setting. This means that this estimator would have an asymptotic coverage of 0.95 if $p$ remained fixed.
In our DGP, the variables $X_j$ with $j\geq p-9$ matter in the outcome equation but are irrelevant in the treatment assignment rule. Given that the propensity score is correctly specified, the four estimators are consistent. The oracle estimator should be the most accurate since it incorporates the information on which covariates matter in the propensity score. Because the balancing program misses some covariates that are relevant for the outcome variable (namely, the $X_j$ with $j\geq p-9$), the naive plug-in estimator is expected to be asymptotically biased. The immunized procedure should correct for this bias.
Table (ref) displays the results. We consider several values of $n$ and $p$ that approximate a relatively high-dimensional setting. For every couple $(n,p)$, we report the root mean squared error (RMSE), the bias and the coverage rate of the confidence intervals associated to each estimator. Our estimator performs well in all settings, with a correct coverage rate and often the lowest RMSE over all estimators. The oracle has always a coverage rate close to 0.95 and a bias very close to 0, as one could expect, but it does not always exhibit the lowest RMSE. This is because, intuitively, the immunized estimator and that of Farrell2015 trade off variance with some bias in their first steps, sometimes resulting in slightly lower RMSE on the final estimator. Note that the bias of the final estimator, though asymptotically negligible, results in a slight undercoverage of the confidence intervals. Yet, even with $n=500$ and $p=1,000$, the coverage rate is still of 0.84, much higher than that of the naive plug-in estimator (0.485). This estimator exhibits a large bias and a low coverage rate even for $p=50$ and $n=2,000$. This shows the importance of correcting for the bias of the first-step estimator, even when $p/n$ is quite small. Finally, the estimator of Farrell2015 exhibits similar performance as the immunized estimator, though it displays a slightly larger RMSE and smaller coverage rate with this particular DGP.
We revisit lal86, who examines the ability of econometric methods to recover the causal effect of employment programs.\footnote{For more discussion on the NSW program and the controversy regarding econometric estimates of causal effects based on nonexperimental data, see lal86 and the subsequent contributions by dw99,DehejiaWahba2002,st05.} This dataset was first built to assess the impact of the National Supported Work (NSW) program. The NSW is a transitional, subsidized work experience program targeted towards people with longstanding employment problems: ex-offenders, former drug addicts, women who were long-term recipients of welfare benefits and school dropouts. The quantity of interest is the average effect of the program for the participants on 1978 yearly earnings. The treated group gathers people who were randomly assigned to this program from the population at risk (with a sample size of $n_1 =185$). Two control groups are available. The first one comes from the Panel Study of Income Dynamics (PSID) (sample size $n_0 =2,490$). The second one comes from the experiment (sample size $n_0=260$) and is therefore directly comparable to the treated group. It provides us with a benchmark for the ATT. Hereafter, we use the group of participants and the PSID sample to compute our estimator and compare it with other competitors and the experimental benchmark.
To allow for a flexible specification, we follow Farrell2015 by taking the raw covariates of the dataset (age, education, black, hispanic, married, dummy variable of no degree, income in 1974, income in 1975, dummy variables of no earnings in 1974 and in 1975), two-by-two-interactions between the continuous and dummy variables, two-by-two interactions between the dummy variables and powers up to degree 5 of the continuous variables. Continuous variables are linearly rescaled to $[0,1]$. We end up with 172 variables to select from. The experimental benchmark for the ATT estimate is \$1,794, with a standard error of 671. We compare several estimators: the naive plug-in estimator, the immunized plug-in estimator, the doubly-robust estimator of Farrell2015, the double-post-selection linear estimator of BCH2012Treat, and the plain OLS estimator including all the covariates. For the four penalized estimators, the penalty loadings on the first-step estimators are set as in Appendix (ref).
Table (ref) displays the results. Our immunized estimator and that of Farrell2015 give credible values for the ATT with respect to the experimental benchmark, with also similar standard errors.\footnote{Farrell2015's estimate shown in Table (ref) differs from that displayed in Farrell's paper because contrary to him, we have not automatically included education, the dummy of no degree and the 1974 income in the set of theory pre-selected covariates. When doing so, the results are slightly better but not qualitatively different for this estimator. We chose not to do so as it would bias the comparison with other estimators, which do not include a set of pre-selected variables. } Notably, only these two estimators out of the five considered display a significant positive impact as the experimental benchmark. The immunized estimator estimator offers a substantial improvement on bias and standard error over the naive plug-in estimator, in line with the evidence from the Monte Carlo experiment. Finally, the OLS estimator in Column (6) is a benchmark of a very simple procedure that does not use any selection at all. As no surprise, this estimator has a substantially larger standard error than the other considered estimators.
Proposition 99 is one of the first and most ambitious large-scale tobacco control program, implemented in 1989 in California. It includes a vast array of measures, including an increase in cigarette taxation of 25 cents per pack, and a significant effort in prevention and education. In particular, the tax revenues generated by Proposition 99 were used to fund anti-smoking campaigns. AbadieDiamondHainmueller2010 analyzes the impact of the law on tobacco consumption in California. Since this program was only enforced in California it is a nice example where the synthetic control method applies whereas more standard public policy evaluation tools cannot be used. It is possible to reproduce a synthetic California by reweighting other states so as to mimic California's behavior.
For this purpose, AbadieDiamondHainmueller2010 consider the following covariates: retail price of cigarettes, state log income per capita, percentage of the population between 15 and 24, per capita beer consumption (all 1980-1988 averages). Cigarette consumptions for the years between 1970 and 1975, 1980 and 1988 are also included. Using the same variables, we conduct the same analysis with our estimator. Figure (ref) displays the estimated effect of Proposition 99 using the immunized estimator.
We find almost no effect of the policy over the pre-treatment period, giving credibility to the counterfactual employed. A steady decline takes place after 1988, and in the long-run, tobacco consumption is estimated to have decreased by about 30 packs per capita per year in California as a consequence of the policy. The variance is larger towards the end of the period because covariates are measured in the pre-treatment period and they become less relevant as predictors. Note also that by construction, including 1970 to 1975, 1980 and 1988 cigarette consumptions among the covariates yields a very good fit at these dates because of the immunization procedure. The fit is not perfect, however, because of the shrinkage induced by the $\ell_1$-penalization.
Figure (ref) displays a comparison between the immunized estimator and the synthetic control method. The dashed green line is the synthetic control counterfactual. It does not match exactly the plot of AbadieDiamondHainmueller2010, in which the weights given to each predictors are optimized to fit best the outcome over the whole pre-treatment period. Instead, the green curve in Figure (ref) optimizes the predictor weights using only the years 1970 through 1975, 1980 and 1988. This strategy brings a fairer comparison with our estimator that does not use California's per capita tobacco consumption outside those dates to optimize the fit. In such a case, the years from 1976 to 1987, excluding 1980, can be used as sort of placebo tests.
Both our estimator and the synthetic control are credible counterfactuals, as they are able to closely match California pre-treatment tobacco consumption. They offer a sizable improvement over a sample average over the U.S. that did not implement any tobacco control program. Furthermore, even if our estimator gives a result relatively similar to the synthetic control, it displays a smoother pattern especially towards the end of the 1980s. The estimated treatment effect appears to be larger with the immunized estimate than with the synthetic control. However, it is hard to conclude that this difference is significant because in the absense of any asymptotic theory on the synthetic control estimator it is unclear how one could make a test on the difference between the two. In fact, the availability of standard asymptotic approximation for confidence intervals is to the advantage of our method.
In this paper, we propose an estimator that makes a link between the synthetic control method, typically used with aggregated data and $n_0 $ smaller than or of the same order as $p$, and treatment effect methods used for micro data for which $n_0 \gg p$. Our method accommodates both settings. In the low-dimensional regime, it pins down one of the solutions of the synthetic control problem, which admits an infinity of solutions. In the high-dimensional regime, the estimator is a regularized and immunized version of the low-dimensional one and then differs from the synthetic control estimator. The simulations and applications suggest that it works well in practice.
In our study, we have focused on specific procedures based on $\ell_1$-penalization and proved that they achieve good asymptotic behavior in possibly high-dimensional regime under sparsity restrictions. Other types of estimators could be explored using these ideas. For example, in the high-dimensional regime, our strategy can be used with the whole spectrum of sparsity-related penalization techniques, such as group Lasso, fused Lasso, adaptive Lasso, Slope, among other.