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.
115,148 characters · 22 sections · 19 citation commands
Robust Inference on Average Treatment Effects with Possibly More Covariates than Observations
\thispagestyle{empty} \setcounter{page}{0}
{\bf Keywords:} High-dimensional sparse model, heterogeneous treatment effects, uniform inference, model selection, doubly-robust estimator, unconfoundedness, group lasso.
{\bf JEL Classification:} C21, C31, C52.
\onehalfspacing
Model selection has always had a place in empirical economics, whether or not it is formally acknowledged. A key problem in modern empirical work is that researchers face datasets with large numbers of variables, sometimes more than observations. A complementary problem is that economic theory and prior knowledge may mandate controlling for certain variables, but are generally silent regarding functional form. These two problems force researchers to search for a model that is simultaneously parsimonious and adequately flexible. Many formal methods are computationally infeasible with a large number of variables. A typical response to this challenge is to iteratively search over a small set of alternative specifications, guided only by the researcher's taste and intuition. But no matter the approach used, subsequent inference almost never takes accounts for this “specification search” and the resulting confidence intervals are not robust to model selection mistakes, and hence are unreliable in empirical work.
This problem is particularly important in estimating average treatment effects under selection on observables, because in this framework using the right covariates is crucial for identification and correct inference. In this context, we provide an easy-to-implement and objective method for covariate selection and post-selection inference on average treatment effects.\footnote{Treatment effects, missing data, measurement error, and data combination models are equivalent under selection on observables. Thus, all our results immediately apply to those contexts. For reviews of these literatures, see \citeasnoun{Tsiatis2006_Book}, \citeasnoun{Heckman-Vytlacil_2007a_Handbook}, \citeasnoun{Imbens-Wooldridge2009_JEL}, and \citeasnoun{Wooldridge2010_book}.} We establish four main results for multivalued treatments effects with arbitrary heterogeneity in observables and heteroskedasticity. First, we show that a doubly-robust estimator is robust to model selection errors. These estimators were initially developed for robustness to parametric misspecification, but are now known to be robust to selection.\footnote{Doubly-robust estimation and its role in program evaluation is discussed by \citeasnoun{Robins-Rotnitzky1995_JASA}, \citeasnoun{vanderLaan-Robins2007_book}, \citeasnoun[with discussion]{Kang-Schafer2007_SS}, \citeasnoun{Tan2010_Bmka}, and references therein.} By taking explicit account of the model selection stage and its inherent selection errors, we derive precise conditions required for any model selector to deliver confidence intervals for average treatment effects that are uniformly valid over a large class of data-generating processes. Second, we show that a simple refitting procedure allows researchers to augment variables chosen according economic theory with data-driven selection to deliver flexible inference that remains uniformly valid. Third, we prove that our estimator is asymptotically linear, and standard conditions imposed in the program evaluation literature, semiparametrically efficient bound. Fourth, we derive new results for multinomial (and binary) logistic regression, the most widely used model for treatment assignment.
Inference following model selection is notoriously difficult. In a sequence of papers, Leeb and P\"otscher Leeb-Potscher2005_ET,Leeb-Potscher2008_ET,Leeb-Potscher2008_JoE,Potscher-Leeb2009_JMA have shown that inference relying too heavily on model selection can not be made uniformly valid. Loosely speaking, uniform validity of a confidence interval captures the idea that the interval should have the same quality (coverage) for many data-generating processes. This theoretical property is practically important because it implies greater reliability in applications. Our proposed methods for post model selection inference build upon the path-breaking recent work of \citeasnoun{Belloni-Chernozhukov-Hansen2014_REStud}.
The crucial insight that leads to uniform inference is to change the goal of model selection away from perfect covariate selection (the oracle property) and to high-quality approximation of the underlying functions. This fundamental shift in focus allows us to circumvent, without contradicting, the impossibility results of Leeb and P\"otscher. Valid post-selection inference has attracted considerable attention during the preparation of this paper: in contexts and with methods quite different from ours, contributions have been made by \citeasnoun{Belloni-Chernozhukov-Wei2013_logit}, \citeasnoun{Berk-etal2013_AoS}, \citeasnoun{Zhang-Zhang2014_JRSSB}, \citeasnoun{Efron2014_JASA}, \citeasnoun{vandeGeer-etal2014_AoS}, and \citeasnoun{Belloni-etal2014_WP}, among others.
Our approach, based on the doubly-robust estimator, has several key features. The name “doubly-robust” reflects that it is robust to misspecification of either the treatment equation (propensity score) or the outcome equation, a property obtained by combining inverse probability weighting and regression imputation. First, we show that this robustness extends to model selection, enabling us to allow for selection errors in both equations without impacting inference. Second, we capture arbitrary treatment effect heterogeneity (dependence of the effect on an individual's observed characteristics), which is crucial in empirical work. With such heterogeneity, the average treatment effect and the treatment on the treated differ, and hence we present results for both. Third, the doubly-robust estimator also stems from the semiparametric efficient moment conditions, and hence we obtain the semiparametric efficiency bound, even under heteroskedasticity, under standard additional conditions. Thus, \possessivecite{Potscher2009_Sankhya} result that sparse estimators have large confidence sets is also circumvented. Taking all these features together enables us to obtain uniform inference over such a large class of treatment effects models.
In recent independent work, \citeasnoun{Belloni-Chernozhukov-Hansen2014_REStud}, propose a similar approach. Their main focus is inference on the linear part of a partially linear model, which motivates an estimator quite different from ours, but it will recover the average treatment effect in the special case of a binary treatment where the effect is constant across observables. However, their Section 5, developed independently from our work, considers heterogeneous effects and proposes an estimator based on the efficient influence function, similar to ours. There are two broad differences. First, we allow for multivalued treatments, which offers a larger set of estimands and can thus enhance the understanding of program impacts.\footnote{Discussion and applications may be found in, for example \citeasnoun{Imbens2000_Bmka}, \citeasnoun{Lechner2001_chapter}, \citeasnoun{Imai-vanDyk2004_JASA}, \citeasnoun{Abadie2005_REStud}, \citeasnoun{Cattaneo2010_JoE}, and \citeasnoun{Cattaneo-Farrell2011_chapter}.} In this context we propose a group lasso based approach that naturally exploits the already-present structure of treatment effects data to improve model selection by pooling information across treatment levels. This is particularly natural in the multivalued case, but even in the binary case there is still a grouped structure in the outcome regressions, though not in treatment assignment (i.e., in propensity score estimation). Second, although in both cases the doubly-robust estimator is used for average treatment effects\footnote{They use different asymptotic variance estimators, and for treatment effects on the treated they do not exploit the simplification discussed in Remark (ref).} (following a quite different model selection step), we show that this estimator has two benefits: (i) it may require weaker conditions on the first stage (see Assumption (ref)); and (ii) it does not require using variables selected for the treatment equation in the outcome model, and vice versa (“post double selection”), and indeed, doing may require additional assumptions (see Assumption (ref)).
Our analysis is conducted under selection on observables, which has a long tradition and remains quite popular in empirical economics.\footnote{For other approaches and reviews of the literature, see, e.g., \citeasnoun{Holland1986_JASA}, \citeasnoun{Hahn1998_Ecma}, \citeasnoun{Horowitz-Manski2000_JASA}, Chen, Hong, and Tarozzi Chen-Hong-Tarozzi2004_WP,Chen-Hong-Tarozzi2008_AoS, \citeasnoun{Bang-Robins2005_Biometrics}, \citeasnoun{Abadie-Imbens2006_Ecma}, \citeasnoun{Wooldridge2007_JoE}, and references therein.} Covariates play three crucial roles in this framework. First, using more observed covariates as proxies, and more flexibly, may help account for unobserved confounding and hence increase the plausibility of unconfoundedness. Second, some observed variables may not be part of the causal mechanism under study, and should be excluded. Third, the efficient conditioning set are those variables that drive the outcome, not necessarily those important for treatment assignment. This reasoning mandates contradicting goals for practitioners: a large, rich set of controls on the one hand, and parsimony on the other. Our approach is a formal, theory-driven attempt to reconcile this contradiction.
A special feature of our analysis is that we match the empirical realities of large data sets by considering selection from amongst (possibly) more covariates than observations, so-called high-dimensional data. The goal of variable selection is to find a small model that is nonetheless sufficiently flexible to capture unknown features of the data-generating process required for inference. If a small model can perfectly capture the unknown feature it is said to be exactly sparse. More realistic is approximate sparsity, when the bias from using a small model is well-controlled, but nonzero. Sparsity is a natural framework for thinking about model selection. Indeed, any time only a few of the available variables are used, a sparsity assumption has effectively been made. It is common empirical practice to report results from several small models, but for these results to be valid one must assume these specifications give high-quality, sparse representations of the unknown features. The alternative we provide involves selecting a sparse, yet flexible, model from among a large set of variables. Results may then be compared with more traditional methods.
With the aim of mimicking common empirical practice we estimate the propensity score with multinomial logistic regression, coupled with group lasso selection Yuan-Lin2006_JRSSB. Our results are stated in the language of treatment effects, but apply to general data structures and are of independent interest in the high-dimensional literature.\footnote{Our techniques build on prior studies, in particular \citeasnoun{Bickel-Ritov-Tsybakov2009_AoS}, \citeasnoun{Lounici-etal2011_AoS}, \citeasnoun{Obozinski-Wainwright-Jordan2011_AoS}, \citeasnoun{Belloni-Chernozhukov2011_AoS}, \citeasnoun{BCCH2012_Ecma}, \citeasnoun{Belloni-Chernozhukov2013_Bern}, and \citeasnoun{Belloni-etal2014_WP}.} Much of the literature has focused on linear models (see \citeasnoun{Buhlmann-vandeGeer2011_book} for a survey), while prior studies of nonlinear models often assume exact sparsity or present limited results.\footnote{Examples include \citeasnoun{vandeGeer2008_AoS} and \citeasnoun{Negahban-etal2012_StatSci}, whose bounds do not imply our results. \citeasnoun{Bach2010_EJS} only gives an error bound on coefficients in exactly sparse logistic regression, which can not yield our results; and does not consider prediction error or post-selection estimation. In independent work, Kwemou2012_logit and \citeasnoun{Belloni-Chernozhukov-Wei2013_logit} also apply \possessivecite{Bach2010_EJS} tools, but are focused on a different goals. \citeasnoun{Vincent-Hansen2014_CSDA} apply the group lasso to multinomial logistic regression, but do not derive any theoretical results.} Furthermore, these studies often use high-level conditions that can be hard to verify. In contrast, we obtain sharp results for logistic regression under the same simple and intuitive conditions used for linear modeling by exploiting mathematical techniques of self-concordant functions put forth by \citeasnoun{Bach2010_EJS}. We also provide extensions to prior work on linear models needed to apply them in treatment effect estimation.
Finally, we offer numerical evidence on the finite sample performance of our procedure. In a small simulation study we find that our procedure delivers very accurate coverage of confidence intervals even for models where covariate selection is difficult, either because of a low signal-to-noise ratio or lack of sparsity, thus highlighting the uniform validity of inference. We also apply our method to the widely-used National Supported Work Demonstration data LaLonde1986_AER and find very accurate estimates and tight confidence intervals (see Table (ref)).
The paper proceeds as follows. Section (ref) gives short, self-contained overview. Section (ref) collects notation. Section (ref) describes the treatment effect models. Sparse models are discussed in Section (ref), which shows how several commonly used models fit in this framework. Section (ref) presents our estimation method and complete results on treatment effect inference. Theoretical results for the group lasso are in Section (ref). Section (ref) presents the numerical evidence and Section (ref) concludes. The main proofs are presented in the Appendix, while the remainder are available in a supplement.
Here we give an overview of the paper, including treatment effect inference (Section (ref)), our new results for the group lasso (Section (ref)), and notation used throughout (Section (ref)).
We consider a multivalued treatment, with status indicated by $D \in \{0, 1, \ldots, \mathcal{T}\}$. Interest lies in mean effects of the treatment on a scalar outcome $Y$. Let $\{Y(t)\}_{t = 0}^\mathcal{T}$ be the (latent) potential outcomes: $Y(t)$ is the outcome a unit would have under $D = t$ and is only observed for units with $D = t$; that is, $Y = \sum_{t = 0}^\mathcal{T} \ensuremath{{\rm 1\hspace*{-0.84ex} \rule{0.06ex}{1.37ex}\hspace*{0.9ex}}\!}\{D = t\} Y(t)$. Many interesting parameters combine means of potential outcomes, and having multivalued treatments allows for a wider range of estimands. Define the mean of one potential outcome as $\mu_t = \mathbb{E}[Y(t)]$. To fix ideas, $\mu_1 - \mu_0$ is the average treatment effect in the binary case ($D \in \{0,1\}$). Sections (ref) and (ref) consider more general average effects, including effects on treated groups. For simplicity, in this section we focus on a single $\mu_t$.
We use the selection on observables framework to identify $\mu_t$. For a vector of covariates $X$, define the generalized propensity score and conditional outcome regressions as \[p_t(x) = \mathbb{P}[ D = t \vert X = x] \qquad \text{ and } \qquad \mu_t(x) = \mathbb{E}[ Y \vert D = t, X = x].\] For identification it is sufficient to assume that $\mathbb{E}[Y(t) \vert D, X] = \mathbb{E}[Y(t) \vert X]$ (mean independence) and $p_t(X)$ is bounded away from zero (overlap) for all treatment levels. Broadly, these two assumptions imply that units from one treatment group are good proxies for other treatments and that there are always such proxies available (see Section (ref)).
For an i.i.d. sample $\{(y_i, d_i, x_i')\}_{i=1}^n$ and model-selection-based estimators $\hat{p}_t(x_i)$ and $\hat{\mu}_t(x_i)$, we estimate $\mu_t$ with \[\hat{\mu}_t = \frac{1}{n} \sum_{i=1}^n \left\{ \frac{ \ensuremath{{\rm 1\hspace*{-0.84ex} \rule{0.06ex}{1.37ex}\hspace*{0.9ex}}\!}\{d_i = t\} (y_i - \hat{\mu}_t(x_i) ) }{ \hat{p}_t(x_i) } + \hat{\mu}_t(x_i) \right\}.\] This doubly-robust estimator combines regression imputation and inverse probability weighting, and remains consistent if either the model $p_t(x)$ or $\mu_t(x)$ is misspecified. Following widespread empirical practice, we estimate $\hat{p}_t(x_i)$ with multinomial logistic regression and $\hat{\mu}_t(x_i)$ linearly (see Section (ref)). The choice of covariates in $\hat{p}_t(x_i)$ and $\hat{\mu}_t(x_i)$ impacts consistency, efficiency, and finite sample performance. Covariate selection based on ad hoc, iterative searches is common in empirical work, but is not formal, objective, or replicable. Balancing tests are also common, but have the additional drawback of assuming the same covariates are important for outcomes and treatment assignment, and more generally do not weight the covariates by their importance for bias.
On the other hand, our proposed procedure gives practitioners an easy to implement, fully objective tool to perform data-driven covariate selection and treatment effect inference, with replicable results.\footnote{For the final step, the doubly-robust estimator is available in STATA and the package of \citeasnoun{Cattaneo-Drukker-Holland2013_stata}. The covariate selection stage is easily implemented in {\sf R}.} Importantly, we do not preclude the addition of variables known to be important from economic theory or prior knowledge. Our procedure is intended to supplement these variables with a flexible set of controls, guarding against misspecification or overfitting.
The following theorem is an example of the more general results presented in Section (ref), wherein we also define $V_t$ and $\hat{V}_t$.
This result establishes the uniform validity of an asymptotic confidence interval for $\mu_t$, overcoming all the post model selection inference challenges: robustness to model selection errors, selecting a model that is small but flexible enough to capture the features of the underlying data generating process, and still retaining efficiency under standard conditions (see Section (ref)). Intuitively, this is similar to (but distinct from) overcoming pretesting bias in other contexts. Also, although our discussion is in terms of covariate selection in high-dimensional, sparse models, the inference result is generic for any first stage estimator.
The two conditions placed on the first stage are analogous to the commonly-used, high-level requirement in semiparametrics that first stage components converge faster than $n^{-1/4}$. However exploiting features of the doubly-robust estimator yields weaker conditions. The first is a mild consistency requirement. The second requires a rate on the product of errors and is thus easier to satisfy if one function is easier to estimate, e.g.\ more smooth or more sparse. In model selection, the rates for the first stage depend on the sample size, the number of covariates considered, and the sparsity level. Importantly, the rate will depend on the total number of covariates only logarithmically, allowing for a large number. We propose to use the group lasso and prove that these estimators satisfy (i) and (ii).
We propose refitting following group lasso selection, and show that it meets all requirements on the model selector. The group lasso is well-suited to program evaluation applications because covariates are penalized according to their overall contribution in all treatment groups. This has two consequences. First, information from all treatments is pooled when doing selection, and hence a weaker signal may be extracted, which improves the selection properties. Second, the selected variables are common to all treatment levels. From a practical point of view this is desirable, as interest rarely lies in a single $\mu_t$, but rather a collection, and substantial commonality is expected in the variables important for different treatment levels.
We consider high-dimensional, sparse models for $p_t(x)$ and $\mu_t(x)$. These are defined by a $p$-dimensional vector $X^*$ based on the original variables $X$. The $X^*$ may consist of any combination of the original variables, interactions, flexible parametric transformations, and/or nonparametric series terms (such as splines or polynomials). A model is approximately sparse if there are $s < n$ of these terms that yield a good approximation ($s\to \infty$ is allowed). To build intuition, suppose that $\mu_t(x)$ obeys a $p$-dimensional linear model. Then the sparsity assumption is that there is an $s$-dimensional submodel with sufficiently small specification bias. In the nonparametric case, sparsity is weaker than (but analogous to) the familiar assumption that a small set of basis functions can approximate the unknown objects well. In practice researchers employ a hybrid of these approaches, which is covered by our results. Section (ref) gives more detail and examples.
We form $\hat{p}_t(x)$ and $\hat{\mu}_t(x)$ in two steps (complete details in Section (ref)). First, the group lasso is applied separately to multinomial logistic and least squares regression to select covariates from $X^*$. We then estimate $p_t(x)$ and $\mu_t(x)$ by refitting unpenalized models using the selected variables, possibly augmented with controls suggested by prior work or economic theory. It is not desirable for a model selector to discard theory and prior work, and our procedure explicitly avoids this. We also allow for using logistic-selected variables in the linear model refitting and vice versa, but this is not necessary for uniformity nor efficiency.
Our main results give precise bounds for the number of covariates selected and the estimation error, both for the penalized and unpenalized estimates. Section (ref) results gives nonasymptotic bounds, with exact constants. Such results are complex and so we give the following intuitive, asymptotic result (The notation $O_{P_n}$ is defined in Section (ref)).
These two results for our proposed group lasso estimators can be directly used to verify the high-level conditions in Theorem (ref) above. Specifically, if $s_d s_y \log(p)^{3 + 2\delta} = o(n)$, conditions (i) and (ii) of Theorem (ref) are met (requiring $s^2 = o(n)$, up to $\log$ factors, as found in other results in the literature). Further, it is clear how the doubly-robust estimator can help: if one function is more smooth or more sparse, $s_d$ or $s_y$ will be lower, easing the restriction. Section (ref) gives further results: showing that the number of variables selected is the same order as the sparsity level, and provides bounds on the logistic and linear coefficients directly. Both these results are important for certain steps in treatment effect estimation that aren't reflected in the simple statement of Theorem (ref). These results appear to be entirely new for the multinomial logistic regression, for any version of the lasso. From a practical point of view, these results provide formal justification for using multinomial logistic regression, coupled with group lasso selection and post-selection refitting.
We collect here notation to be used for the rest of the paper. The data generating process (DGP) is denoted by $P_n$ and is defined by the joint law of the random variables $(Y,D, X')'$. For a given $n$, $\{(y_i, d_i, x_i')'\}_{i=1}^n$ constitute draws from $P_n$. The DGP may vary with $n$, along with features such as parameters, distributions, and so forth, as discussed in Section (ref). This is generally suppressed for clarity. We adopt the following conventions.
For a set of scalars $\{m_t\}_{t = 1}^{\mathcal{T}}$, let $\hat{p}_t( \{m_t\}_{\mathbb{N}_\mathcal{T}} ) = \exp(m_t) [1 + \sum_{t \in \mathbb{N}_\mathcal{T}} \exp( m_t) ]^{-1}$ denote the multinomial logit function. Empirical expectation will be denoted $\mathbb{E}_n[w_i] = \sum_{i=1}^n w_i/n$ and $\mathbb{E}_{n,t}[w_i] = \sum_{i \in \mathbb{I}_t} w_i / n_t = \sum_{i=1}^n d_i^{t} w_i / n_t$.
In this section we formally define the treatment effects model and the parameters of interest. Recall that $D \in \{0, 1, \ldots, \mathcal{T}\}$ indicates treatment status, $\{Y(t)\}_{t \in \overline{\mathbb{N}}_\mathcal{T}}$ are the (latent) potential outcomes, and $Y(t)$ is only observed for units with $D = t$; that is, $Y = \sum_{t \in \overline{\mathbb{N}}_\mathcal{T}} Y(t)$. The building blocks of many general estimands are the averages
In the binary case, the average treatment effect is $\mu_1 - \mu_0$ and the treatment on the treated is $\mu_{1,1} - \mu_{0,1}$. A multivalued treatment allows for a large range of interesting estimands. To fix ideas, we keep as running examples two leading cases. First, the so-called dose-response function: the $(\mathcal{T} + 1)$-vector $\bm{\mu} = (\mu_0, \mu_1, \ldots, \mu_\mathcal{T})'$. Second, define $\bm{\tau}$ as the $\mathcal{T}$-vector with element $t$ given by $\mu_{t,t} - \mu_{0,t}$. This gives the effect of each treatment relative to the baseline $t=0$, only for those who received that treatment. These are by no means the only interesting estimands constructed from $\mu_t$ and $\mu_{t,t'}$; many others are given by \citeasnoun{Lechner2001_chapter}, \citeasnoun{Heckman-Vytlacil_2007a_Handbook}, and others.
The following two conditions are sufficient to identify $\mu_t$ and $\mu_{t,t'}$.
This assumption is a form of “ignorability” coined by \citeasnoun{Rosenbaum-Rubin1983_Bmka}. This model allows arbitrary treatment effect heterogeneity in observables, but not unobservables. This assumption is standard in the program evaluation literature, and its plausibility has been discussed at length, so we omit a general discussion (see, e.g., \citeasnoun{Imbens2004_REStat}, \citeasnoun[Chapter 21]{Wooldridge2010_book}, and references therein). However, in the context of model selection, three remarks are warranted.
First, in place of (ref), it is more common to instead assume full conditional independence: $Y \protect\mathpalette{\protect\independenT}{\perp} D \vert X$. However, as observed by \citeasnoun{Heckman-Ichimura-Todd1997_REStud}, the weaker mean independence is sufficient. For our purposes, the “gap” between the two assumptions is important. Suppose full independence holds only conditional on a set of variables strictly larger than the variables entering the mean functions (e.g.\ the excess variables affect higher moments). In this case, because mean independence is still sufficient, we need not aim to select the larger set. Full independence is important for the efficiency discussed in Section (ref).
Second, the covariates may, in general, include instruments for treatment status, but they are not known as such. This is standard, but left implicit, in discussions of ignorability. If instruments are present, and selected for estimation, efficiency suffers but unbiasedness is not harmed. Efficiency bounds in this context typically (implicitly) assume there are no instruments in $X$. Assumption (ref) rules out perfect predictors. Section (ref) offers further discussion.
Finally, the main drawback of Assumption (ref) is that it does not give identification of average effects on transformations of $Y(t)$. However, we are expressly interested in model selection on the mean function of the level of $Y(t)$, and hence Assumption (ref) is more natural. To operationalize model selection, structure must be placed on $\mathbb{E}[Y(t) \vert X = x]$, and hence functional form conditions tied to mean independence are not limiting per se. If the parameter of interest is changed, say to $\mathbb{E}[\log(Y(t))]$, and a sparsity assumption is made for $\mathbb{E}[\log(Y(t)) \vert X = x]$, then our method applies.
Assumption (ref) yields identification of $\mu_t$ and $\mu_{t,t'}$ using either inverse weighting or regression, and double robustness follows from combining the two strategies. Recall the notation $p_t(x) = \mathbb{P}[D = t \vert X = x]$ and $\mu_t(x) = \mathbb{E}[Y \vert D = t, X = x]$. Applying Assumption (ref) we find that
and
where $p_t = \mathbb{P}[D = t]$. The moment condition (ref) holds if either $p_t(x)$ or $\mu_t(x)$ is misspecified. For $\mu_{t,t'}$, if $\mu_t(x)$ is misspecified, both $p_t(X)$ and $p_{t'}(X)$ must be correctly specified, while if $\mu_t(x)$ is correct, both propensity scores may be misspecified. It is important to note that the forms of $\psi_t(\cdot)$ and $\psi_{t, t'}(\cdot)$ are fixed, so the function itself does not depend on the sample size even if its arguments do. Our estimator is a plug-in version of this moment condition.
We now formalize approximate sparsity. Let $X_Y^*$ and $X_D^*$ be $p$-dimensional transformations of the covariates $X$, with $p>n$ allowed. These transformations are specific to the outcome and treatment models, but may overlap. They do not vary with $t$, nor depend on the DGP. Some examples are given below in Section (ref). For the multinomial logistic model it is convenient to work with the log-odds ratio. We take $p_0(x) = 1 - \sum_{t \in \mathbb{N}_\mathcal{T}} p_t(x)$ and write
Similarly, write the outcome regressions as
The terms $B_t^D = B_t^D(x)$ and $B_t^Y = B_t^Y(x)$ are bias terms arising from the parametric specification. As discussed below, these encompass the usual nonparametric bias as well. Approximate sparsity requires that only a small number of the $X^*$ are needed to make the bias small. Define $ S_*^D = \bigcup_{\mathbb{N}_\mathcal{T}} \operatorname*{supp}(\gamma^*_t)$ and $ S_*^Y =\bigcup_{\overline{\mathbb{N}}_\mathcal{T}} \operatorname*{supp}(\beta^*_t) $, so that these sets capture all variables important for treatment and outcomes, respectively. We assume that there are some $s_d < n$ and $s_y < n$, such that for $| S_*^D | = s_d$ and $|S_*^Y | = s_y $, the biases $B_t^D$ and $B_t^Y$ are sufficiently small. This is made precise by defining the bounds:
Note that the former bias bound is placed directly on the propensity score because it is the ultimate object of interest, rather than on the linearization of the log-odds.
While a great deal of overlap is expected, in practice it is likely that a few covariates will be more or less important for different treatments, and so we do not require that the supports of $\gamma^*_t, t \in \mathbb{N}_\mathcal{T}$ or $\beta^*_t, t \in \overline{\mathbb{N}}_\mathcal{T}$ are constant over $t$, nor that $S_*^D$ overlaps with $S_*^Y$. Instead, it may be better to think of $\mathbb{N}_p \setminus S_*^D$ and $\mathbb{N}_p \setminus S_*^Y$ as the “common nonsupports” of the treatment and outcome equations. When it is clear from the context we will abbreviate both $X_D^*$ and $X_Y^*$ by $X^*$ (and their realizations by $x_i^*$) and refer to them generically as “covariates”, and further write $s$ for either $s_d$ or $s_y$. We assume $\mathbb{E}_n[({x_{i,j}^*})^2] = 1$ without loss of generality (see Remark (ref)).
To concretize the sparse model idea, we now discuss how several models commonly used in practice fit into this framework. These include parametric and nonparametric models for $p_t(x)$ and $\mu_t(x)$, and hybrids of these. A common theme to all examples will be comparison to the oracle model: the model that knows the true support in advance. Our uniform inference results include all these examples as special cases because, loosely speaking, we obtain uniformity over DGPs where $p_t(x)$ and $\mu_t(x)$ have sparse representations. We aim for an accessible discussion of each model, and defer technicalities to the literature Raskutti-Wainwright-Yu2010_JMLR,Rudelson-Zhou2013_IEEE,Belloni-Chernozhukov-Hansen2014_REStud.
It is important to note that misspecification of the type guarded against by double robustness can arise in any type of model. In parametric cases, this is most often functional form misspecification. While this type of misspecification does not occur in nonparametrics, others are possible, such as shape restrictions or separability assumptions being incorrect, or omitting relevant variables. None of these errors disappear asymptotically, and all of them are guarded against by use of the doubly-robust estimator.
Much of the DGP, including parameters and distributions, is allowed to depend on $n$. Perhaps the most salient features that do not depend on $n$ are the set of treatments and the functions $\psi_t$ and $\psi_{t, t'}$. It is likely that our results can be extended to accommodate a growing number of treatments, but that is beyond the scope of our study. In the models (ref) and (ref), $X^*$, $\gamma^*_{\bm{\cdot}, \bm{\cdot}}$, and $\beta^*_{\bm{\cdot}, \bm{\cdot}}$ must depend on $n$ by construction. Our results on estimation of these models are nonasymptotic. For treatment effect inference, we use triangular array asymptotics to retain the dependence on $n$ of the DGP. The interpretation of the results does, and should, change depending on what is assumed about the DGP. To illustrate, let us return to Examples (ref) and (ref).
First, consider the simple parametric models of Example (ref). In this case, $\mu_t = \mathbb{E}[ \mathbb{E}[Y(t) \vert X]] = \mathbb{E}[ X']\beta^*_t$, which depends on $n$ by construction, as the dimension is diverging. It may seem unnatural that the parameter to be estimated depends on $n$, as we typically think of “true” parameters being features of a (large) fixed study population. However, with a diverging number of covariates, there is no fixed DGP. Indeed, if we estimate $\mu_t = \mu_t^{(n_1)}$ based upon $n_1$ observations, and then proceed to gather $n_2$ more observations, when we re-estimate our target is now $\mu_t^{(n_1 + n_2)} \neq \mu_t^{(n_1)}$. One possible resolution is as follows. First, the parameter of interest is $\mu_t^{(\infty)} = \mathbb{E}[Y(t)]$, which is defined without reference to covariates. We can view each successive $n$-dependent $\mu_t$ as an approximation of $\mu_t^{(\infty)}$ based upon $p = p_n$ covariates. Note well that in our thought experiment, $p_{n_1} \neq p_{n_1 + n_2}$, and so additional variables should have been collected for all $n_1 + n_2$ samples.
Contrast this with the semiparametric model in Example (ref). It is common to assume the population DGP is fixed over $n$. The treatment effects may be constructed in terms of the underlying variables, e.g.\ $\mu_t^{(\infty)} = \mathbb{E}[Y(t)] = \mathbb{E}[\mathbb{E}[Y(t) \vert X]]$, with $X^*$ serving only the purpose of aiding in approximating the regression functions. Model selection is performed on series terms, not underlying variables, to estimate the coefficients $\gamma^*_{\bm{\cdot}, \bm{\cdot}}$ and $\beta^*_{\bm{\cdot}, \bm{\cdot}}$. If $\mu_t = \mathbb{E}[{X_Y^*}'] \beta^*_t + \mathbb{E}[B_t^Y]$ does not depend on $n$, the bias term, by definition, exactly compensates for the $n$-dependence in $\mathbb{E}[{X_Y^*}'] \beta^*_t$. We emphasize that our inference results allow for general $n$-dependence in the DGP, and interpretation by the econometrician must take careful account of any conceptual assumptions.
In this section we present results on uniformly valid treatment effect inference. We first present the estimators and conditions required for a generic first stage to yield uniform inference. Although our focus is on model selection and sparsity, our results are more general, showcasing the benefits of doubly robust estimation for any model in Section (ref) where Assumption (ref) below (which does not refer to selection or sparsity) can be satisfied.
The moment functions $\psi_t(\cdot)$ and $\psi_{t, t'}(\cdot)$ of Equations (ref) and (ref) have fixed and known form, and so for estimators $\hat{p}_t(x)$ and $\hat{\mu}_t(x)$, we can define
and
where $\hat{p}_t = n_t / n$. By combining these estimators appropriately we can construct estimators $\hat{\bm{\mu}}$ and $\hat{\bm{\tau}}$ for the dose-response function $\bm{\mu}$ and the vector $\bm{\tau}$, respectively, and any other estimand. Notice that when $t=t'$ $\hat{\mu}_{t,t}$ is an average over the appropriate subpopulation: $\hat{\mu}_{t,t} = \mathbb{E}_{n,t}[ y_i]$.
Although in this section we allow for generic estimates $\hat{p}_t(x)$ and $\hat{\mu}_t(x)$, it is important to distinguish between estimates based upon selected sets that have no “additional randomness” and those that do. Model selection based estimation will naturally have two steps: first data-driven selection and then refitting to ameliorate the shrinkage bias and allow the researcher to augment the selected variables. Let $\tilde{S}^D$ and $\tilde{S}^Y$ be the selected sets and $\hat{S}^D$ and $\hat{S}^Y$ be the final sets of variables used in the refitting. We will say that these contain no “additional randomness” if the added variables (i.e. $\hat{S} \setminus \tilde{S}$, for $Y$ or $D$) are nonrandomly selected, such as from economic theory or prior knowledge. On the other hand, the added variables may be selected from a random process beyond that included in $\tilde{S}$. The leading example would be using logistic-selected variables in the regressions or vice versa. Then the variables used in $\hat{\mu}_t(x_i)$ depend not only on the randomness of $\tilde{S}^Y$, but also on that of $\tilde{S}^D$, and hence on $\{ d_i \}_{i = 1}^n$. Additional conditions are required for the estimators with additional randomness.
The choice of method is in part dependent on the assumptions of the underlying model. To illustrate, first, return to Example (ref), where we have a purely parametric model with $X = X_D^* = X_Y^*$. The researcher may want to set $\hat{S}^D \supset \tilde{S}^D \cup \tilde{S}^Y$, in order to have a better chance that $S_*^Y \subset \hat{S}^D$. The set $\hat{S}^D$ now contains additional randomness due to $\tilde{S}^Y$. Conversely, consider Example (ref). It is natural to include “low-order” basis functions for each underlying covariate, say linear and quadratic polynomials. Thus, the researcher may want to include these in $\hat{S}$, whether or not selected by the group lasso. However, there is no reason that the series terms useful for approximating the functions $\mu_t(x)$ would be useful for $p_t(x)$, or vice versa, and no additional randomness is injected.
We now state the sufficient conditions used for treatment effect estimation and inference. For exposition, we present these in three groups: those concerning the underlying DGP, requirements of $\hat{p}_t(x)$ and $\hat{\mu}_t(x)$ in the “no additional randomness” case, and finally the additional conditions to allow for “additionally random” selected sets. Begin with conditions on the DGP. Let $U \equiv Y(t) - \mu_t(X)$ and impose the following conditions.
These conditions are mild and intuitive, and not unique to high-dimensional models or model selection. Assumption (ref) restricts attention to cross-sectional applications. The condition of bounded covariates is unlikely to be a limitation in practice. Any $X^*$ that are underlying variables will naturally be bounded in applications. This condition is automatically satisfied for most common choices of basis functions employed in nonparametric estimation. The rest are moment conditions on the potential outcome models, including allowing the errors to be heteroskedastic and non-Gaussian. The uniform bounds in $n$ are needed for array asymptotics.
We now give precise conditions on the model selector sufficient for uniformly valid inference.
These two collectively play the same role as the commonly-used, high-level requirement in semiparametrics that each first-step component separately converge at $n^{-1/4}$ at least.\footnote{See \citeasnoun{Newey-McFadden1994_handbook} and \citeasnoun{Chen2007_handbook}, and references therein.} Indeed, \citeasnoun{Belloni-Chernozhukov-Hansen2014_REStud} employ just such a condition for each component. However, by making use of the doubly-robust property we have the weaker conditions shown.\footnote{Many studies in the semiparametric literature relax or do not rely on the $n^{1/4}$ condition, allowing the nonparametric portion to converge at a slower rate, at any rate, or in some cases be inconsistent; examples include \citeasnoun{Powell-Stock-Stoker1989_Ecma}, \citeasnoun{Newey1990_Ecma}, \citeasnoun{Robins-etal2008_IMS}, \citeasnoun{Cattaneo-Jansson-Newey2014_alt}, and Cattaneo, Crump, and Jansson Cattaneo-Crump-Jansson2013_JASA,Cattaneo-Crump-Jansson2014_ET, among others.} The first is a mild consistency requirement. The second requires an explicit rate on the product of errors, and hence if one function is relatively easy to estimate Assumption (ref) can be satisfied even if the other does not converge at $n^{-1/4}$. This formalizes the benefit of doubly-robust estimation in general. In high-dimensional, sparse modeling specifically the rates for the first stage depend on the sample size, the number of covariates considered, and the sparsity level. Thus, if one function requires fewer covariates to estimate, i.e. smaller $p$ or $s$, then greater complexity can be allowed for in the other (capturing, in particular, their relative smoothness).
The so-called “additional-randomness” estimators are more specific to the (approximately) sparse model context, and so we now codify the sparsity requirements of Section (ref) and then give the additional conditions required for these estimators.
These conditions are needed to apply bounds for self-normalized sums delaPena-Lai-Shao2009_book. \citeasnoun{BCCH2012_Ecma} were the first to use these techniques in high-dimensional, sparse models. The first condition is high-level, but can be verified with conditions on the errors and a bound for estimation. For the former, \citeasnoun{BCCH2012_Ecma} assume that $ \max_{i \in \mathbb{N}_n} |u_i| = O_{P_n}(n^{1/q})$ for some $q>2$. A larger $q$ eases the restriction in Assumption (ref) but at the expense of stronger conditions on the noise distribution. For example, if $u_i$ are assumed Gaussian, $q$ can be taken to be any (large) positive number.
We now come to our main results on inference on average treatment effects. Most of our discussion will concern $\mu_t$ and $\bm{\mu}$; similar points apply to results for $\mu_{t,t'}$ and $\bm{\tau}$. Our first result formalizes consistency of our estimates under misspecification.
This theorem formalizes the double-robustness property of our estimators: the propensity score or regression may be misspecified if the limiting objects are well-behaved. Compare to Assumption (ref). The nearly identical result for $\mu_{t,t'}$ is omitted to save space.
We now turn to our main inference results. First we demonstrate a Bahadur representation of a generic $\hat{\mu}_t$ or $\hat{\mu}_{t,t'}$. These are shown to be equivalent to a sample average of the moment functions $\psi_t(\cdot)$ and $\psi_{t,t'}(\cdot)$, respectively, after proper centering and scaling, evaluated at the true $p_t(x_i)$ and $\mu_t(x_i)$. Using these results, asymptotic normality can be obtained for general estimands. We state explicit results for the leading examples $\bm{\mu}$ and $\bm{\tau}$.
An asymptotic variance formula is needed to state the results. Define the conditional variance of the potential outcomes s $\sigma_t^2(x) = \mathbb{E}[U^2 \vert D = t, X = x]$ and the $\overline{\mathcal{T}}$-square matrix $V_{\bm{\mu}}$ with elements \[V_{\bm{\mu}}[t, t'] = \ensuremath{{\rm 1\hspace*{-0.84ex} \rule{0.06ex}{1.37ex}\hspace*{0.9ex}}\!}\{t = t'\} \mathbb{E}\left[ \frac{\sigma_t^2(X)}{p_t(X)}\right] + \mathbb{E}\left[(\mu_t(X) - \mu_t)(\mu_{t'}(X) - \mu_{t'})\right] \equiv V_{\bm{\mu}}^W(t) + V_{\bm{\mu}}^B(t,t').\] Straightforward plug-in estimators for these two components are given by\footnote{Estimators can also be based on sample averages of outer products of influence functions, which would include the covariance term that vanishes in expectation.} \[\hat{V}_{\bm{\mu}}^W(t) = \mathbb{E}_n\left[ \frac{ d_i^{t} (y_i - \hat{\mu}_t(x_i) )^2}{\hat{p}_t(x_i)^2}\right] \qquad \text{ and } \qquad \hat{V}_{\bm{\mu}}^B(t,t') = \mathbb{E}_n\left[ (\hat{\mu}_t(x_i) - \hat{\mu}_t) (\hat{\mu}_{t'}(x_i) - \hat{\mu}_{t'}) \right].\]
Our first result gives the asymptotic behavior of $\hat{\mu}_t$ and $\hat{\bm{\mu}}$ for a sequence of DGPs.
Theorem (ref) itself may appear standard, but what is nonstandard is that the model selection step of the estimation has been explicitly accounted for. This immediately gives the following uniform inference results.
Corollary (ref) shows that these procedures are uniformly valid over the class of DGPs we consider, and hence will be reliable in applications. The crucial insight that leads to uniform inference is to change the goal of model selection away from perfect covariate selection (the oracle property) and to high-quality approximation of the underlying functions (here $p_t(\cdot)$ and $\mu_t(\cdot)$). This fundamental shift in focus allows us to avoid the uniformity problems demonstrated by Leeb and P{\"o}tscher. Assumption (ref) formalizes exactly the quality of approximation needed. Such an approximation can be found for any element in $\bm{P}_n$, and hence inference is uniformly valid over that class. This method of proving uniformity follows \citeasnoun{Belloni-Chernozhukov-Hansen2014_REStud} and \citeasnoun{Romano2004_SJS}, and is distinct from the approach of \citeasnoun{Andrews-Guggenberger2009_JoE}.
Results for the treatment effects on the treated are similar. The variance formula for $\bm{\tau}$ is slightly more cumbersome. Define the $\mathcal{T}$-square matrix $V_{\bm{\tau}}$ with elements
Straightforward plug-in estimators for these two components are given by \[\hat{V}_{\bm{\tau}}^W(t) = \mathbb{E}_n\left[ \frac{ d_i^{t} }{\hat{p}_t^2} \left[ \left(y_i - \hat{\mu}_0(x_i) - \hat{\mu}_{t,t} + \hat{\mu}_{0,t}\right)^2 \right] \right] \text{ and } \hat{V}_{\bm{\tau}}^B(t,t') = \mathbb{E}_n\left[ \frac{\hat{p}_t(x_i) \hat{p}_{t'}(x_i) }{\hat{p}_t \hat{p}_{t'} \hat{p}_0(x_i)^2} d_i^0(y_i - \hat{\mu}_0(x_i) )^2 \right].\] Note that we needn't estimate $\mu_t(x)$ and $\sigma_t^2(x)$, due to the simplification in Remark (ref). With this notation, we have the following results. Proofs are so similar to those for Theorem (ref) and Corollary (ref) that we omit them.
The prior theoretical results are aimed at delivering robust inference. In this section, we briefly discuss the efficiency of our estimator according to two criteria: semiparametric efficiency and oracle efficiency. To put each on sound conceptual footing we separate discussion and restrict to an appropriate set of models.
For semiparametric efficiency, $p_t(x)$ and $\mu_t(x)$ are nonparametric objects, as in Example (ref), $X$ are fixed-dimension variables and the DGP does not vary with $n$. If we “upgrade” the mean independence of Assumption (ref) to full, namely $ \{ Y(t) \}_{\overline{\mathbb{N}}_\mathcal{T}} \protect\mathpalette{\protect\independenT}{\perp} D \vert X$, then Theorems (ref) and Theorem (ref) immediately yield asymptotic linearity and semiparametric efficiency, attaining \possessivecite{Hahn1998_Ecma} or \possessivecite{Cattaneo2010_JoE} bounds. This requires there be no (known) instruments for treatment status in $X$, as implicitly assumed in those works, else the bound may change Hahn2004_REStat.
Turning to oracle efficiency, an alternative to our robust approach is to prove that the true support can be found with probability approaching one (the oracle property), then conduct inference conditioning on this event. This approach cannot be made uniformly valid, but may be of interest in the exactly sparse models of Example (ref) (there is no “true” support in approximately sparse models), because discovering the true support is equivalent to finding the variables in the causal mechanism White-Lu2011_REStat, if one exists. This may be interesting in its own right, or for future applications by way of hypothesis generation. The post oracle selection estimator is made efficient by using only the variables important for $\mu_t(x_i) = \mathbb{E}[Y \vert D = t, x_i]$. This amounts to entirely removing the instrumental variables indexed by $S_*^D \setminus S_*^Y$, whose inclusion would, in general, reduce efficiency, though not increase bias. Further, $S_*^Y \setminus S_*^D$ are excluded from propensity score estimation.
Perfect selection requires two strong conditions: (i) an orthogonality condition on the Gram matrices that restricts the correlation between the variables in and out of the true support Bach2008_JMLR, and (ii) a beta-min condition bounding the nonzero coefficients away from zero. Intuitively, highly correlated variables cannot be distinguished, nor can coefficients sufficiently close to zero be found with certainty. Both bounds may depend on $n$, and in particular the lower bound on the coefficients may vanish at an appropriate rate. Under such conditions, it is straightforward to show that $S_*^Y$ and $S_*^D$ can be found with probability approaching one.
We now give details for group lasso model selection and estimation, and make the refitting precise. Section (ref) discusses penalty choices and implementation. Restricted and sparse eigenvalues, key quantities in our bounds, are discussed in Section (ref). Our main nonasymptotic results are stated in Section (ref). These results are of interest more generally in the literature on high-dimensional sparse models Finally, Section (ref) gives asymptotic rates and verifies the conditions of Section (ref).
We first select covariates by applying the group lasso penalty to the multinomial logistic loss (for the propensity scores) and to least squares loss (to estimate the outcome regression). The loss functions are defined as \[\mathcal{M}(\gamma_{\bm{\cdot}, \bm{\cdot}}) = \sum_{t \in \mathbb{N}_\mathcal{T}} \mathbb{E}_n \left[ - d_i^{t} \log\left( \hat{p}_t(\{{x_i^*}'\gamma_t\}_{\mathbb{N}_\mathcal{T}}) \right) \right] \qquad \text{ and } \qquad \mathcal{E}(\beta_{\bm{\cdot}, \bm{\cdot}}) = \sum_{t \in \overline{\mathbb{N}}_\mathcal{T}} \mathbb{E}_{n,t}[ (y_i - {x_i^*}' \beta_t)^2].\] Then, the group lasso estimates for the propensity score coefficients, denoted $\tilde{\gamma}_{\bm{\cdot}, \bm{\cdot}}$, and the regression coefficients, $\tilde{\beta}_{\bm{\cdot}, \bm{\cdot}}$, respectively solve
where $\lambda_D$ and $\lambda_Y$ are the penalty parameters discussed below and ${\left\vert\kern-0.25ex\left\vert\kern-0.25ex\left\vert \gamma_{\bm{\cdot}, \bm{\cdot}} \right\vert\kern-0.25ex\right\vert\kern-0.25ex\right\vert_{2,1}}$ is the mixed $\ell_2$/$\ell_1$ norm.
To ameliorate the downward bias induced by the penalty and to allow for researcher-added variables, we refit unpenalized models.\footnote{The bias is away from the pseudo-true coefficients of the sparse parametric representation, $\gamma^*_{\bm{\cdot}, \bm{\cdot}}$ and $\beta^*_{\bm{\cdot}, \bm{\cdot}}$. There is no relation to specification biases $B_t^D$ and $B_t^Y$.} Let $\tilde{S}^D = \{j : \| \tilde{\gamma}_{\bm{\cdot}, j} \|_2 > 0 \}$ and $\tilde{S}^Y = \{j : \| \tilde{\beta}_{\bm{\cdot}, j} \|_2 > 0 \}$ be the selected covariates and $\hat{S}^D$ and $\hat{S}^Y$ those used in refitting.\footnote{When $\operatorname*{supp}(\gamma^*_t)$ and $\operatorname*{supp}(\beta^*_t)$ do not vary much over $t$, the group lasso is known to have better properties than the ordinary lasso in terms of selection and convergence. \citeasnoun{Obozinski-Wainwright-Jordan2011_AoS} give a sharp bound on the overlap necessary to yield improvements, while \citeasnoun{Huang-Zhang2010_AoS}, \citeasnoun{Kolar-Lafferty-Wasserman2011_JMLR}, and \citeasnoun{Lounici-etal2011_AoS} also demonstrate advantages of the group lasso approach. These works show, among other things, that the group lasso advantage increases with large $\mathcal{T}$, and with the group structure, may perform better with smaller samples. We defer to the works cited for a formal discussion.} We require $\hat{S} \supset \tilde{S}$ and $|\hat{S}| \leq s$ for $D$ and $Y$ (we will prove that $| \tilde{S}| \leq s$ in both cases). The refitting estimators solve
We must specify choices of $\lambda_D$ and $\lambda_Y$ for programs (ref). From a theoretical point of view, these must be chosen so that the penalty dominates the noise, which is captured by the magnitude of the score in the dual of the ${\left\vert\kern-0.25ex\left\vert\kern-0.25ex\left\vert \cdot \right\vert\kern-0.25ex\right\vert\kern-0.25ex\right\vert_{2,1}}$ norm, with high probability. To acheive this, we set
for some $\delta_D > 0$ and $\delta_Y > 0$. With these choices, $\lambda_D > 2 \max_{j \in \mathbb{N}_p} \| \mathbb{E}_n [(p_t(x_i) - d_i^{t}) x_{i,j}^* ] \|_2$ and $\lambda_Y > 4 \max_{j \in \mathbb{N}_p} \| \mathbb{E}_{n,t} [u_i x_{i,j}^* ] \|_2$ with probability $1 - \mathcal{P}$ for a small (and shrinking) $\mathcal{P}$. In generic terms, $\lambda$ is of the form $\Lambda (1 + r_n)$, where $\Lambda$ is an upper bound on the true score and $r_n$ is a rate that depends on $n$ and $p$.\footnote{The slight differences in the two are as follows. The full sample has information on the logistic coefficients, so $n$ appears instead of $\underline{n}$. No error bound appears in $\lambda_D$ because the errors are bounded by one. The multiple $4$ for $\lambda_Y$, instead of 2, can be traced to the quadratic loss. These forms are determined at heart by the maximal inequality of \citeasnoun{Lounici-etal2011_AoS}.} The specific rate chosen serves to balances the rate of convergence against the concentration effect: a smaller $r_n$ would increase the rate of convergence, but at the expensive of lowering the concentration probability $1 - \mathcal{P}$. In the Appendix we show that (for appropriate $\delta$ and $n$ or $\underline{n}$) the concentration probability is given by
There are two practical methods to make these choices for feasible for implementation. When $\hat{p}_t(x)$ and $\hat{\mu}_t(x)$ are used to estimate average treatment effects, the decreased sensitivity of the final estimate to the first stage, thanks to the doubly-robust estimator, in turn results in less sensitivity to the choice of penalty (through the sparsity).\footnote{To our knowledge, no formal results exist on “optimal” penalty parameter choices for inference in high-dimensional problems nor are any procedures free of user-specified choices.} The first option is an iterative procedure to estimate the unknown $\mathcal{X}$ and $\mathcal{U}$ in $\lambda_Y$ and $\lambda_D$, as employed by \citeasnoun{BCCH2012_Ecma} (validity of this procedure may be established along the same lines as in that study). We use $\max_{i \leq n } \max_{j \in \mathbb{N}_p} |x_{i,j}^*|$ for $\mathcal{X}$ and estimate $\mathcal{U}$ by iteration: given an initial estimate $\hat{\mu}_t^{(0)}(x)$, set $\hat{\mathcal{U}}^{(k)} = \mathbb{E}_n[(y_i - \hat{\mu}_t^{(k-1)}(x_i))^4]^{1/4}$, where $\hat{\mu}_t^{(k)}(x_i)$, $k > 0$, is based on Eqn.\ (ref). In implementation we found 10 iterations more than sufficient, and based the initial estimate on ridge regression (with penalty chosen by cross validation). A second option is to select $\lambda_Y$ and $\lambda_D$ directly by cross-validation. This has the appealing feature that the precise forms of Eqn.\ (ref) need not be characterized and estimated. If interest lies in the underlying functions $p_t(x)$ and $\mu_t(x)$, cross validation is appropriate as it minimizes a relevant loss function. Formal results establishing the validity of cross-validation are not available, but it performs well in practice.
The local behavior of optimizations (ref) and (ref) is captured by their respective Hessians, which involve the second moment matrix of the covariates. The eigenvalues of such matrices will be explicit in our bounds. We are interested in finite sample bounds, and so we will only discuss the empirical Gram matrices (see Remark (ref)). Define
In high-dimensional data, both are singular, and so we use restricted eigenvalues and sparse eigenvalues Bickel-Ritov-Tsybakov2009_AoS.
For the multinomial logistic regression, the minimal restricted eigenvalue is defined by
For least squares estimation we instead use
Note that $Q$ appears for $\kappa_D$, whereas the $Q_t$ are used in $\kappa_Y$. The restricted set, or cone constraint, requires the magnitude of $\delta_{\bm{\cdot}, \bm{\cdot}}$ off the true support be small relative to the true support, measured in the group lasso norm.\footnote{The multiplier of 4 in the constraint for $\kappa_D$ is traceable to the nonlinear model.} We will show that $(\tilde{\gamma}_{\bm{\cdot}, \bm{\cdot}} - \gamma^*_{\bm{\cdot}, \bm{\cdot}})$ and $(\tilde{\beta}_{\bm{\cdot}, \bm{\cdot}} - \beta^*_{\bm{\cdot}, \bm{\cdot}})$ obey the respective constraints.
In contrast, the refitting errors $(\hat{\gamma}_{\bm{\cdot}, \bm{\cdot}} - \gamma^*_{\bm{\cdot}, \bm{\cdot}})$ and $(\hat{\beta}_{\bm{\cdot}, \bm{\cdot}} - \beta^*_{\bm{\cdot}, \bm{\cdot}})$ from (ref) may not obey the cone constraint, but are sparse by construction. This motivates the use of sparse eigenvalues. For a set $S \subset \mathbb{N}_p$ and a $p \times p$ matrix $\tilde{Q}$, define
Finally, it will be useful to define a bound on $\overline{\phi}\{\tilde{Q}, S\}$ over all subsets of a certain size. To this end, for any integer $m$, define $\overline{\overline{\phi}}(\tilde{Q},m) = \max_{S \subset \mathbb{N}_p,\, |S| \leq m} \overline{\phi}\{\tilde{Q}, S\}$.
We take these quantities to be primitive, and defer to the literature. For example, \citeasnoun{vandeGeer-Buhlmann2009_EJS}, \citeasnoun{Huang-Zhang2010_AoS}, \citeasnoun{Raskutti-Wainwright-Yu2010_JMLR}, \citeasnoun{Rudelson-Zhou2013_IEEE}, and \citeasnoun{Belloni-Chernozhukov-Hansen2014_REStud}. In particular, \citeasnoun{Huang-Zhang2010_AoS} show that the group lasso may need fewer observations to satisfy conditions on $\underline{\phi}\{\tilde{Q}, S\}$.
We now have the necessary notation and assumptions to state our theoretical results on group lasso estimation, beginning with multinomial logistic regression, followed by a terse treatment of linear models. Corollary (ref) is a special case of the results in this section, see Section (ref).
Our first result is a nonasymptotic bound on the group lasso estimates from (ref).
This theorem is new to the literature, to the best of our knowledge. Much of the detail involves capturing the finite sample behavior of the Hessian and Gram matrices. We discuss the features of this result in the following remarks.
The error bounds for post-selection estimation are more complex and depend in part on the good properties of the initial group lasso fit. The following theorem gives our results.
It is not readily discernible if these bounds improve upon the initial fit. This will depend on the DGP, the selection success of the initial fit, and any added variables. In this result, further lower bounds on $A_K$ are required to handle the sparse eigenvalues, compared to the restricted version in Theorem (ref). The role played by $A_K$ is the same in both cases, as with the other factors.
It is worth noting that, despite the complexity of multinomial logistic regression, the conditions for Theorems (ref) and (ref) are simple and intuitive, and match those used for linear models.
We now give our results for group lasso estimation of the conditional outcome regressions. In computing $\mu_t(x_i)$ for $d_i^{t} \neq 1$ we are performing out of sample prediction, which slightly complicates the bounds. Our first result is on the initial group lasso fit.
This theorem generalizes \citeasnoun{Lounici-etal2011_AoS} to the nonparametric, approximately sparse case, improves the sparsity bound, and gives out of sample prediction (imputation) results. The analogous generalization for within sample prediction loss (e.g.\ multi-task learning), $\mathbb{E}_{n,t}[({x_i^*}'\tilde{\beta}_t - \mu_t(x_i))^2] ^{1/2}$, may be found in the Supplement.
For refitting, we are predicting for the entire sample and so we utilize the general results given by \citeasnoun{BCCH2012_Ecma} for post-selection estimation of least squares. The following result is a direct implication of their Lemma 7 and our Theorem (ref).
As above, the performance of the refitting procedure depends in part on the success of the initial group lasso fit. Indeed, the middle term is dropped if the true support union is found. The constants $A_k$, k=1, 2, 3, 4 are not given explicitly but are known to be absolute bounds delaPena-Lai-Shao2009_book under Assumption (ref). This result is less precise than Theorems (ref) and (ref), but sufficient to verify Assumptions (ref) and (ref).
This section derives rates of convergence for the group lasso estimates and uses these results to verify Assumptions (ref) and (ref) in Section (ref). For simplicity, we only state results for the post-selection estimators that we recommend in practice. In reducing the finite sample results of Theorems (ref) and (ref) to rates we retain the dependence on $n$, $p$, $s$, and the bias. Note that the number of treatments is fixed, and the overlap assumption ensures that all $n_t \propto n$. Further, the various (restricted and sparse) eigenvalues are commonly taken to be bounded (or bounded away from zero) in asymptotic analyses. This accounts for the remaining factors in the bounds. For multinomial logistic regression, we obtain the following result.
Similarly, we have the following for the linear models.
It is now straightforward to verify the requirements of Section (ref). Assumption (ref) requires \[(n^{-1} s_d \log(p \vee n)^{3/2 + \delta_D} + (b_{s}^d)^2 ) ( n^{-1} s_y \log(p \vee n)^{3/2 + \delta_Y} + ( b_{s}^y)^2 ) = o\left(n^{-1}\right).\] Under the common assumption that $b_{s} = O(\sqrt{s/n})$, we require $s_d s_y \log(p \vee n)^{3 + \delta_D + \delta_Y} = o(n)$. Both this, and the display above, clearly show how the sparsity and smoothness of the two functions interact due to the double robustness. Assumption (ref) can be verified similarly.
These rates of convergence (i.e. part 2 of each corollary) are optimal up to factor $\log(p \vee n)^{1/2 + \delta}$. At heart, this loss appears to stem from the maximal inequality used to establish the concentration probability of (ref). In practice, this is unlikely to be a limitation. As mentioned above, the use of group lasso can yield improvements in the constants if the data obey a grouped sparsity pattern, as is expected for treatment effects data, and may even yield improvements in the detection of the sparse signal, further offsetting the suboptimal $\log$ factor (see for example \citeasnoun{Lounici-etal2011_AoS} or \citeasnoun{Obozinski-Wainwright-Jordan2011_AoS}). Alternative methods could, in principle, yield a rate improvement. Chief among these would be lasso-penalized linear probability models (see also Remark (ref)) or separate logistic regressions. The group lasso approach adopted here reflects common practice, and so it may be preferred. In any case, the $\log$ factors do not impact the treatment effect inference.
We conducted a Monte Carlo exercise to study how our estimator behaves as the propensity score and regression functions change, and the model selection problem becomes more or less difficult.\footnote{The supplemental appendix contains the additional results.} For simplicity we focus on the average effect of a binary treatment. We generated 1000 observations $(y_i, d_i, x_i')'$ from the models in Example (ref), using both $p=1000$ and $p=1500$. The covariates include an intercept, with the remainder drawn from $N(0,\Sigma)$, with covariance $\Sigma[j_1,j_2] = 2^{-|j_1 - j_2|}, 2 \leq j_1, j_2 \leq p$. Errors are standard Normal. The crucial aspects of the DGP are the coefficient vectors $\beta_0^0$, $\beta_1^0$, and $\gamma^0$, which are defined to vary with the positive scalars $\rho_\beta$, $\rho_\gamma$, $\alpha_\beta$, and $\alpha_\gamma$, as follows:
with $\beta_1^0 = - \beta_0^0$. The $\rho$ multipliers affect the signal-to-noise ratio, but not the sparsity. For smaller values distinguishing the large and small coefficients is more difficult for a given sample. The exponents $\alpha$ control the sparsity, where a sparse representation is not possible for small values.
Figure (ref) shows the empirical coverage rates of 95% confidence intervals for $\mu_1 - \mu_0$ for different DGPs, for $p=1000$ and $1500$. Panels (a) and (c) show coverage as the multipliers $\rho_\beta$ and $\rho_\gamma$ range over 0.01 (weak signal) to 1 (strong), with $\alpha_\beta = \alpha_\gamma=2$. Panels (b) and (d) vary the sparsity exponents $\alpha_\beta$ and $\alpha_\gamma$ over 1/8 (not sparse) to 4 (very sparse), with $\rho_\beta = \rho_\gamma = 1$. Of 1000 observations total, the (mean) size of the comparison group declines from roughly 500 to 300 as $\rho_\gamma$ increases and 450 to 300 as $\alpha_\gamma$ increases, over their given ranges. Coverage is accurate over all signal strengths, and breaks down only when neither $\mu_t(x_i)$ nor $p_t(x_i)$ is sparse, which is exactly when Assumption (ref) (or condition (ii) of Theorem (ref)) cannot be satisfied. Note that coverage accuracy is retained when only one function is sparse, showcasing the double-robustness property.
The penalty parameters $\lambda_D$ and $\lambda_Y$ are chosen using the iterative procedure described in Section (ref), with $\delta_D = 4.5$ and $\delta_Y = 5$ throughout. Different DGPs exhibit different sensitivity to these values. Results using penalties chosen via 10-fold cross-validation appear in Figure (ref), which also exhibits excellent coverage across all sparse designs.\footnote{The R routines appear unstable for nonsparse designs, thus the analogues to Panels (b) and (d) of Figure (ref) are omitted. See the supplement for limited versions. This will be explored for future software development.}
To illustrate the role that model selection can play in a real-world application, we revisit the National Supported Work (NSW) demonstration. The NSW has been analyzed numerous times since \citeasnoun{LaLonde1986_AER}. Our aim is a simple study of model selection, not a comprehensive or conclusive evaluation of the NSW. We focus on the subsample used by \citeasnoun{Dehejia-Wahba1999_JASA} and the Panel Study of Income Dynamics (PSID) comparison sample, taking as given their data definitions, sample selection, and trimming rules. Detailed discussion of these choices, and the NSW program may be found in Dehejia and Wahba Dehejia-Wahba1999_JASA,Dehejia-Wahba2002_REStat (hereafter DW99 and DW02) and \citeasnoun{Smith-Todd2005_JoE}, and references therein. Briefly, the outcome of interest is earnings following a job training program. The dataset includes a treatment indicator, post-treatment earnings (1978), two years of pre-treatment earnings (1974\footnote{This naming follows DW99, but the variable may be measured outside 1974, see discussion in the works cited.} and 1975), as well as age, education, a marital status, and indicators for Black and Hispanic. Thus, $X$ consists of seven variables. We will keep the estimator fixed: all estimates will be based on the doubly-robust estimator with standard errors from Section (ref). We will compare the following specifications for $X^*$:
For specifications 1 and 2, the same covariates are in the outcome and treatment models. All specifications include an intercept and we include education and pre-treatment income in the refitting step following model selection. We follow DW99 and DW02 and trim comparisons with estimated propensity score larger (smaller) than the maximum (minimum) in the treated sample.\footnote{A formal treatment of trimming is beyond the scope of the present study. The goal of our analysis is illustrative, and hence we take DW99's trimming as given. This issue is discussed by DW99, DW02, and \citeasnoun{Smith-Todd2005_JoE}.}
Table (ref) presents results from these three specifications, and includes the experimental arm of the NSW. The group lasso based estimate performs very well: the point estimate is accurate and the interval is tight. Selecting from 171 possible covariates allows for a great deal of flexibility, but the sparsity of the estimate keeps the variance well-controlled. The no-selection point estimate is accurate, but fails to yield significance, while the specification of DW02 yields a significant, but overly high estimate and wide confidence interval. The benefits of explicit model selection are clear.
This paper proposed a method that achieves uniformly valid inference on mean effects of a multivalued treatment even after model selection among possibly more covariates than observations. We demonstrated robustness to model selection errors, misspecification, and heterogeneous effects in observables. To accomplish this, a doubly-robust estimator was employed and shown to have excellent properties following model selection. We proved new results on group lasso estimation, which we argue is natural for treatment effects data. Multinomial logistic regression was studied in some detail. Numerical evidence shows that our method is quite promising for applications.
A key outstanding question in this work and in the high-dimensional, sparse modeling literature more generally, is penalty parameter choice. Very little work has been done in this area, which is a crucial gap in implementability of these techniques. We plan to develop a formal choice for the penalty parameter that is appropriately optimal. Tuning parameter selection in semi- and nonparametric analysis, and its impact on estimation and inference, is becoming better understood, and parallel developments must take place in model selection contexts.