EconBase
← Back to paper

Inference on Treatment Effects After Selection Amongst High-Dimensional Controls

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.

108,926 characters · 14 sections · 8 citation commands

Rendered from LaTeX for readability, not typeset faithfully. Citation keys are highlighted; maths is left as source; figures, tables and equation environments are summarised rather than reproduced; unrecognised commands are greyed out so nothing is silently dropped. Email addresses are removed.

Inference on Treatment Effects After Selection Amongst High-Dimensional Controls

abstractWe propose robust methods for inference on the effect of a treatment variable on a scalar outcome in the presence of very many controls. Our setting is a partially linear model with possibly non-Gaussian and heteroscedastic disturbances where the number of controls may be much larger than the sample size. To make informative inference feasible, we require the model to be approximately sparse; that is, we require that the effect of confounding factors can be controlled for up to a small approximation error by conditioning on a relatively small number of controls whose identities are unknown. The latter condition makes it possible to estimate the treatment effect by selecting approximately the right set of controls. We develop a novel estimation and uniformly valid inference method for the treatment effect in this setting, called the “post-double-selection" method. Our results apply to Lasso-type methods used for covariate selection as well as to any other model selection method that is able to find a sparse model with good approximation properties. The main attractive feature of our method is that it allows for imperfect selection of the controls and provides confidence intervals that are valid uniformly across a large class of models. In contrast, standard post-model selection estimators fail to provide uniform inference even in simple cases with a small, fixed number of controls. Thus our method resolves the problem of uniform inference after model selection for a large, interesting class of models. We illustrate the use of the developed methods with numerical simulations and an application to the effect of abortion on crime rates. \\ Key Words: treatment effects, partially linear model, high-dimensional-sparse regression, inference under imperfect model selection, uniformly valid inference after model selection

Introduction

Many empirical analyses in economics focus on estimating the structural, causal, or treatment effect of some variable on an outcome of interest. For example, we might be interested in estimating the causal effect of some government policy on an economic outcome such as employment. Since economic policies and many other economic variables are not randomly assigned, economists rely on a variety of quasi-experimental approaches based on observational data when trying to estimate such effects. One important method is based on the assumption that the variable of interest can be taken as randomly assigned once a sufficient set of other factors has been controlled for. Economists, for example, might argue that changes in state-level public policies can be taken as randomly assigned relative to unobservable factors that could affect changes in state-level outcomes after controlling for aggregate macroeconomic activity, state-level economic activity, and state-level demographics; see, for example, \citeasnoun{heckman:metricslabormarkets} or \citeasnoun{imbens:review}.

A problem empirical researchers face when relying on an identification strategy for estimating a structural effect that relies on a conditional on observables argument is knowing which controls to include. Typically, economic intuition will suggest a set of variables that might be important but will not identify exactly which variables are important or the functional form with which variables should enter the model. This lack of clear guidance about what variables to use leaves researchers with the problem of selecting a set of controls from a potentially vast set of control variables including raw regressors available in the data as well as interactions and other transformations of these regressors. A typical economic study will rely on an ad hoc sensitivity analysis in which a researcher reports results for several different sets of controls in an attempt to show that the parameter of interest that summarizes the causal effect of the policy variable is insensitive to changes in the set of control variables. See \citeasnoun{levitt:abortion}, which we use as the basis for the empirical study in this paper, or examples in \citeasnoun{AngristBook} among many other references.

We present an approach to estimating and performing inference on structural effects in an environment where the treatment variable may be taken as exogenous conditional on observables that complements existing strategies. We pose the problem in the framework of a partially linear model

equation[equation omitted — 72 chars of source]

where $d_i$ is the treatment/policy variable of interest, $z_i$ is a set of control variables, and $\zeta_i$ is an unobservable that satisfies $\textnormal{E}[\zeta_i\mid d_i,z_i] = 0$.\footnote{ We note that $d_i$ does not need to be binary.} The goal of the econometric analysis is to conduct inference on the treatment effect $\alpha_0$. We examine the problem of selecting a set of variables from among $p$ potential controls $x_i=P(z_i)$, which may consist of $z_i$ and transformations of $z_i$, to adequately approximate $g(z_i)$ allowing for $p > n$. Of course, useful inference about $\alpha_0$ is unavailable in this framework without imposing further structure on the data. We impose such structure by assuming that exogeneity of $d_i$ may be taken as given once one controls linearly for a relatively small number $s < n$ of variables in $x_i$ whose identities are a priori unknown. This assumption implies that a linear combination of these $s$ unknown controls provides an approximation to $g(z_i)$ which produces relatively small approximation errors.\footnote{We carefully define what we mean by small approximation errors in Section 2.} This assumption, which is termed approximate sparsity or simply sparsity, allows us to approach the problem of estimating $\alpha_0$ as a variable selection problem. This framework allows for the realistic scenario in which the researcher is unsure about exactly which variables or transformations are important for approximating $g(z_i)$ and so must search among a broad set of controls.

The assumed sparsity includes as special cases the most common approaches to parametric and nonparametric regression analysis. Sparsity justifies the use of fewer variables than there are observations in the sample. When the initial number of variables is high, the assumption justifies the use of variable selection methods to reduce the number of variables to a manageable size. In many economic applications, formal and informal strategies are often used to select such smaller sets of potential control variables. Most of these standard variable selection strategies are non-robust and may produce poor inference.\footnote{An example of inference going wrong is given in Figure 1 (left panel), presented in the next section, where a standard post-model selection estimator has a bimodal distribution which sharply deviates from the standard normal distribution. More examples are given in Section 6 where we document the poor inferential performance of a standard post-model selection method. } In an effort to demonstrate robustness of their conclusions, researchers often employ ad hoc sensitivity analyses which examine the robustness of inferential conclusions to variations in the set of controls. Such sensitivity analyses are useful but lack rigorous justification. As a complement to these ad hoc approaches, we propose a formal, rigorous approach to inference allowing for selection of controls. Our proposal uses modern variable selection methods in a novel manner which results in robust and valid inference.

The main contributions of this paper are providing a robust estimation and inference method within a partially linear model with potentially very high-dimensional controls and developing the supporting theory. The method relies on the use of Lasso-type or other sparsity-inducing procedures for variable selection. Our approach differs from usual post-model-selection methods that rely on a single selection step. Rather, we use two different variable selection steps followed by a final estimation step as follows:

itemize• In the first step, we select a set of control variables that are useful for predicting the treatment $d_i$. This step helps to insure robustness by finding control variables that are strongly related to the treatment and thus potentially important confounding factors. • In the second step, we select additional variables by selecting control variables that predict $y_{i}$. This step helps to insure that we have captured important elements in the equation of interest, ideally helping keep the residual variance small as well as intuitively providing an additional chance to find important confounds. • In the final step, we estimate the treatment effect $\alpha_0$ of interest by the linear regression of $y_{i}$ on the treatment $d_i$ and the union of the set of variables selected in the two variable selection steps.

We provide theoretical results on the properties of the resulting treatment effect estimator and show that it provides inference that is uniformly valid over large classes of models and also achieves the semi-parametric efficiency bound under some conditions. Importantly, our theoretical results allow for imperfect variable selection in either of the two variable selection steps as well as allowing for non-Gaussianity and heteroscedasticity of the model's errors.\footnote{In a companion paper that presents an overview of results for $\ell_1$-penalized estimators, \citeasnoun{BCH2011:InferenceGauss}, we provide similar results in the idealized Gaussian homoscedastic framework. }

We illustrate the theoretical results through an examination of the effect of abortion on crime rates following \citeasnoun{levitt:abortion}. In this example, we find that the formal variable selection procedure produces a qualitatively different result than that obtained through the ad hoc set of sensitivity results presented in the original paper. By using formal variable selection, we select a small set of between eight and fourteen variables depending on the outcome, compared to the set of eight variables considered by \citeasnoun{levitt:abortion}. Once this set of variables is linearly controlled for, the estimated abortion effect is rendered imprecise. It is interesting that the key variable selected by the variable selection procedure is the initial condition for the abortion rate. The selection of this initial condition and the resulting imprecision of the estimated treatment effect suggest that one cannot determine precisely whether the effect attributed to abortion found when this initial condition is omitted from the model is due to changes in the abortion rate or some other persistent state-level factor that is related to relevant changes in the abortion rate and current changes in the crime rate.\footnote{Note that all models are estimated in first-differences to eliminate any state-specific factors that might be related to both the relevant level of the abortion rate and the level of the crime rate.} It is interesting that \citeasnoun{FooteGoetzAbortion} raise a similar concern based on intuitive grounds and additional data in a comment on \citeasnoun{levitt:abortion}. \citeasnoun{FooteGoetzAbortion} find that a linear trend interacted with crime rates before abortion could have had an effect renders the estimated abortion effects imprecise.\footnote{\citeasnoun{DLAbortionResponse} provide yet more data and a more complicated specification in response to \citeasnoun{FooteGoetzAbortion}. In a supplement available at http://faculty.chicagobooth.edu/christian.hansen/research/, we provide additional results based on \citeasnoun{DLAbortionResponse}. The conclusions are similar to those obtained in this paper in that we find the estimated abortion effect becomes imprecise once one allows for a broad set of controls and selects among them. However, the specification of \citeasnoun{DLAbortionResponse} relies on a large number of district cross time fixed effects and so does not immediately fit into our regularity conditions. We conjecture the methodology continues to work in this case but leave verification to future research.} Overall, finding that a formal, rigorous approach to variable selection produces a qualitatively different result than a more ad hoc approach suggests that these methods might be used to complement economic intuition in selecting control variables for estimating treatment effects in settings where treatment is taken as exogenous conditional on observables.

Relationship to literature. We contribute to several existing literatures. First, we contribute to the literature on series estimation of partially linear models (\citeasnoun{donald:newey:pl}, \citeasnoun{hardle:pl}, \citeasnoun{robinson}, and others). We differ from most of the existing literature which considers $p\ll n$ series terms by allowing $p \gg n$ series terms from which we select $\widehat s \ll n$ terms to construct the regression fits. Considering an initial broad set of terms allows for more refined approximations of regression functions relative to the usual approach that uses only a few low-order terms. See, for example, \citeasnoun{BCH2011:InferenceGauss} for a wage function example and Section 5 for theoretical examples. However, our most important contribution is to allow for data-dependent selection of the appropriate series terms. The previous literature on inference in the partially linear model generally takes the series terms as given without allowing for their data-driven selection. However, selection of series terms is crucial for achieving consistency when $p \gg n$ and is needed for increasing efficiency even when $p =C n$ with $C<1$. That the standard estimator can be be highly inefficient in the latter case follows from results in \citeasnoun{CJN:PLMStandardError}.\footnote{ \citeasnoun{CJN:PLMStandardError} derive properties of series estimator under $p = Cn$, $C<1$, asymptotics. It follows from their results that under homoscedasticity the series estimator achieves the semiparametric efficiency bound only if $C\to 0$. } We focus on Lasso for performing this selection as a theoretically and computationally attractive device but note that any other method, such as selection using the traditional generalized cross-validation criteria, will work as long as the method guarantees sufficient sparsity in its solution. After model selection, one may apply conventional standard errors or the refined standard errors proposed by \citeasnoun{CJN:PLMStandardError}.\footnote{If the selected number of terms $\widehat s$ is a substantial fraction of $n$, we recommend using \citeasnoun{CJN:PLMStandardError} standard errors after applying our model selection procedure. }

Second, we contribute to the literature on the estimation of treatment effects. We note that the policy variable $d_i$ does not have to be binary in our framework. However, our method has a useful interpretation related to the propensity score when $d_i$ is binary. In the first selection step, we select terms from $x_i$ that predict the treatment $d_i$, i.e. terms that explain the propensity score. We also select terms from $x_i$ that predict $y_i$, i.e. terms that explain the outcome regression function. Then we run a final regression of $y_i$ on the treatment $d_i$ and the union of selected terms. Thus, our procedure relies on the selection of variables relevant for both the propensity score and the outcome regression. Relying on selecting variables that are important for both objects allows us to achieve two goals: we obtain uniformly valid confidence sets for $\alpha_0$ despite imperfect model selection and we achieve full efficiency for estimating $\alpha_0$ in the homoscedastic case. The relation of our approach to the propensity score brings about interesting connections to the treatment effects literature. \citeasnoun{hahn:prop}, \citeasnoun{heckman:ichimura:smith:mathching}, and \citeasnoun{abadie:imbens} have constructed efficient regression or matching-based estimates of average treatment effects. \citeasnoun{hahn:prop} also shows that conditioning on the propensity score is unnecessary for efficient estimation of average treatment effects. \citeasnoun{hirano:imbens:ridder} demonstrate that one can efficiently estimate average treatment effects using estimated propensity score weighting alone. \citeasnoun{robins:dr} have shown that using propensity score modeling coupled with a parametric regression model leads to efficient estimates if either the propensity score model or the parametric regression model is correct. While our contribution is quite distinct from these approaches, it also highlights the important robustness role played by the propensity score model in the selection of the right control terms for the final regression.

Third, we contribute to the literature on estimation and inference with high-dimensional data and to the uniformity literature. There has been extensive work on estimation and perfect model selection in both low and high-dimensional contexts,\footnote{For reviews focused on econometric applications, see, e.g., \citeasnoun{Hansen2005} and \citeasnoun{BellChernHans:Gauss}.} but there has been little work on inference after imperfect model selection. Perfect model selection relies on unrealistic assumptions, and model selection mistakes can have serious consequences for inference as has been shown in \citeasnoun{potscher}, \citeasnoun{leeb:potscher:pms}, and others. In work on instrument selection for estimation of a linear instrumental variables model, \citeasnoun{BellChenChernHans:nonGauss} have shown that model selection mistakes do not prevent valid inference about low-dimensional structural parameters due to the inherent adaptivity of the problem: Omission of a relevant instrument does not affect consistency of an IV estimator as long as there is another relevant instrument. The partially linear regression model ((ref)) does not have the same adaptivity structure, and model selection based on the outcome regression alone produces non-robust confidence intervals.\footnote{The poor performance of inference on a treatment effect after model selection on only the outcome equation is shown through simulations in Section 6.} Our post-double selection procedure creates the necessary adaptivity by performing two separate model selection steps, making it possible to perform robust/uniform inference after model selection. The uniformity holds over large, interesting classes of high-dimensional sparse models. In that regard, our contribution is in the spirit and builds upon the classical contribution by \citeasnoun{romano:uniform} on the uniform validity of t-tests for the univariate mean. It also shares the spirit of recent contributions, among others, by \citeasnoun{mikusheva} on uniform inference in autoregressive models, by \citeasnoun{andrews:cheng} on uniform inference in moment condition models that are potentially unidentified, and by \citeasnoun{andrews:cheng:guggen} on a generic framework for uniformity analysis.

Finally, we contribute to the broader literature on high-dimensional estimation. For variable selection we use $\ell_1$-penalization methods, though our method and theory will allow for the use of other methods. $\ell_1$-penalized methods have been proposed for model selection problems in high-dimensional least squares problems, e.g. Lasso in \citeasnoun{FF:1993} and \citeasnoun{T1996}, in part because they are computationally efficient. Many $\ell_1$-penalized methods have been shown to have good estimation properties even when perfect variable selection is not feasible; see, e.g., \citeasnoun{CandesTao2007}, \citeasnoun{MY2007}, \citeasnoun{BickelRitovTsybakov2009}, \citeasnoun{horowitz:lasso}, \citeasnoun{BC-PostLASSO} and the references therein. Such methods have also been shown to extend suitably to nonparametric and non-Gaussian cases as in \citeasnoun{BickelRitovTsybakov2009} and \citeasnoun{BellChenChernHans:nonGauss}. These methods also produce models with a relatively small set of variables. The last property is important in that it leaves the researcher with a set of variables that may be examined further; in addition it corresponds to the usual approach in economics that relies on considering a small number of controls.

Paper Organization. In Section 2, we formally present the modeling environment including the key sparsity condition and develop our advocated estimation and inference method. We establish the consistency and asymptotic normality of our estimator of $\alpha_0$ uniformly over large classes of models in Section 3. In Section 4, we present a generalization of the basic procedure to allow for model selection methods other than Lasso. In Section 5, we present a series of theoretical examples in which we provide primitive condition that imply the higher-level conditions of Section 3. In Section 6, we present a series of numerical examples that verify our theoretical results numerically, and we apply our method to the abortion and crime example of \citeasnoun{levitt:abortion} in Section 7. In appendices, we provide the proofs.

Notation. In what follows, we work with triangular array data $\{\left(\omega_{i,n}, i=1,...,n\right), n=1,2,3,...\}$ defined on probability space $(\Omega, \mathcal{A}, {\mathrm{P}}_n)$, where ${\mathrm{P}} = {\mathrm{P}}_n$ can change with $n$. Each $\omega_{i,n}= (y_{i,n}', z_{i,n}', d_{i,n}')'$ is a vector with components defined below, and these vectors are i.n.i.d. -- independent across $i$, but not necessarily identically distributed. Thus, all parameters that characterize the distribution of $\{\omega_{i,n}, i=1,...,n\}$ are implicitly indexed by ${\mathrm{P}}_n$ and thus by $n$. We omit the dependence on these objects from the notation in what follows for notational simplicity. We use array asymptotics to better capture some finite-sample phenomena and to insure the robustness of conclusions with respect to perturbations of the data-generating process ${\mathrm{P}}$ along various sequences. This robustness, in turn, translates into uniform validity of confidence regions over certain regions of data-generating processes.

We use the following empirical process notation, ${\mathbb{E}_n}[f] := {\mathbb{E}_n}[f(\omega_i)] := \sum_{i=1}^n f(\omega_i)/n,$ and $\mathbb{G}_n(f) := \sum_{i=1}^n ( f(\omega_i) - {\mathrm{E}}[f(\omega_i)] )/\sqrt{n}.$ Since we want to deal with i.n.i.d. data, we also introduce the average expectation operator: $ \bar {\mathrm{E}}[f] := {\mathrm{E}} {\mathbb{E}_n}[f] = {\mathrm{E}} {\mathbb{E}_n}[f(\omega_i)] = \sum_{i=1}^n {\mathrm{E}}[f(\omega_i)]/n. $ The ${l}_2$-norm is denoted by $\|\cdot\|$, and the ${l}_0$-norm, $\|\cdot\|_0$, denotes the number of non-zero components of a vector. We use $\| \cdot \|_{\infty}$ to denote the maximal element of a vector. Given a vector $\delta \in {\Bbb{R}}^p$, and a set of indices $T \subset \{1,\ldots,p\}$, we denote by $\delta_T \in {\Bbb{R}}^p$ the vector in which $\delta_{Tj} = \delta_j$ if $j\in T$, $\delta_{Tj}=0$ if $j \notin T$. We use the notation $(a)_+ = \max\{a,0\}$, $a \vee b = \max\{ a, b\}$, and $a \wedge b = \min\{ a , b \}$. We also use the notation $a \lesssim b$ to denote $a \leqslant c b$ for some constant $c>0$ that does not depend on $n$; and $a\lesssim_P b$ to denote $a=O_P(b)$. For an event $E$, we say that $E$ wp $\to$ 1 when $E$ occurs with probability approaching one as $n$ grows. Given a $p$-vector $b$, we denote $\text{support}(b) = \{ j \in \{1,...,p\}: b_j \neq 0\} $.

Inference on Treatment and Structural Effects Conditional on Observables

Framework

We consider the partially linear model

eqnarray[eqnarray omitted — 198 chars of source]

where $y_{i}$ is the outcome variable, $d_i$ is the policy/treatment variable whose impact $\alpha_0$ we would like to infer, $z_i$ represents confounding factors on which we need to condition, and $\zeta_i$ and $v_i$ are disturbances. The parameter $\alpha_0$ is the average treatment or structural effect under appropriate conditions given, for example, in \citeasnoun{heckman:metricslabormarkets} or \citeasnoun{imbens:review} and is of major interest in many empirical studies.

The confounding factors $z_i$ affect the policy variable via the function $m(z_i)$ and the outcome variable via the function $g(z_i)$. Both of these functions are unknown and potentially complicated. We use linear combinations of control terms $x_i = P(z_i)$ to approximate $g(z_i)$ and $m(z_i)$, writing ((ref)) and ((ref)) as

eqnarray[eqnarray omitted — 203 chars of source]

where $ x_i'\beta_{g0}$ and $x_i'\beta_{m0}$ are approximations to $g(z_i)$ and $m(z_i)$, and $r_{gi}$ and $r_{mi}$ are the corresponding approximation errors. In order to allow for a flexible specification and incorporation of pertinent confounding factors, the vector of controls, $x_i = P(z_i)$, can have a dimension $p=p_n$ which can be large relative to the sample size. Specifically, our results only require $\log p = o(n^{1/3})$ along with other technical conditions. High-dimensional regressors $x_i = P(z_i)$ could arise for different reasons. For instance, the list of available controls could be large, i.e. $x_i=z_i$ as in e.g. \citeasnoun{koenker:jappliedeconometircs}. It could also be that many technical controls are present; i.e. the list $x_i=P(z_i)$ could be composed of a large number of transformations of elementary regressors $z_i$ such as B-splines, dummies, polynomials, and various interactions as in \citeasnoun{newey:series} or \citeasnoun{chen:Chapter}.

Having very many controls creates a challenge for estimation and inference. A key condition that makes it possible to perform constructive estimation and inference in such cases is termed sparsity. Sparsity is the condition that there exist approximations $ x_i'\beta_{g0}$ and $x_i'\beta_{m0}$ to $g(z_i)$ and $m(z_i)$ in ((ref))-((ref)) that require only a small number of non-zero coefficients to render the approximation errors $r_{gi}$ and $r_{mi}$ sufficiently small relative to estimation error. More formally, sparsity relies on two conditions. First, there exist $\beta_{g0}$ and $\beta_{m0}$ such that at most $s=s_n \ll n$ elements of $\beta_{m0}$ and $\beta_{g0}$ are non-zero so that $$ \|\beta_{m0}\|_0 \leqslant s \text{ and } \|\beta_{g0}\|_0 \leqslant s. $$ Second, the sparsity condition requires the size of the resulting approximation errors to be small compared to the conjectured size of the estimation error: $$ \{\bar {\mathrm{E}}[r^{2}_{gi}]\}^{1/2} \lesssim \sqrt{s/n} \text{ and } \{\bar {\mathrm{E}}[r^{2}_{mi}]\}^{1/2} \lesssim \sqrt{s/n}. $$ Note that the size of the approximating model $s=s_n$ can grow with $n$ just as in standard series estimation.

The high-dimensional-sparse-model framework outlined above extends the standard framework in the treatment effect literature which assumes both that the identities of the relevant controls are known and that the number of such controls $s$ is much smaller than the sample size. Instead, we assume that there are many, $p$, potential controls of which at most $s$ controls suffice to achieve a desirable approximation to the unknown functions $g(\cdot)$ and $m(\cdot)$ and allow the identity of these controls to be unknown. Relying on this assumed sparsity, we use selection methods to select approximately the right set of controls and then estimate the treatment effect $\alpha_0$.

The Method: Least Squares after Double Selection

We propose the following method for estimating and performing inference about $\alpha_0$. The most important feature of this method is that it does not rely on the highly unrealistic assumption of perfect model selection which is often invoked to justify inference after model selection. To the best of our knowledge, our result is the first of its kind in this setting. This result extends our previous results on inference under imperfect model selection in the instrumental variables model given in \citeasnoun{BellChenChernHans:nonGauss}. The problem is fundamentally more difficult in the present paper due to lack of adaptivity in estimation which we overcome by introducing additional model selection steps. The construction of our advocated procedure reflects our effort to offer a method that has attractive robustness/uniformity properties for inference. The estimator is $\sqrt{n}$-consistent and asymptotically normal under mild conditions and provides confidence intervals that are robust to various perturbations of the data-generating process that preserve approximate sparsity.

To define the method, we first write the reduced form corresponding to ((ref))-((ref)) as:

eqnarray[eqnarray omitted — 154 chars of source]

where $\bar \beta_0 := \alpha_0 \beta_{m0} + \beta_{g0}, \ \ \bar r_i := \alpha_0 r_{mi} + r_{gi}, \ \ \bar \zeta_i := \alpha_0 v_i + \zeta_i.$

We have two equations and hence can apply model selection methods to each equation to select control terms. The chief method we discuss is the Lasso method described in more detail below. Given the set of selected controls from ((ref)) and ((ref)), we can estimate $\alpha_0$ by a least squares regression of $y_{i}$ on $d_i$ and the union of the selected controls. Inference on $\alpha_0$ may then be performed using conventional methods for inference about parameters estimated by least squares. Intuitively, this procedure works well since we are more likely to recover key controls by considering selection of controls from both equations instead of just considering selection of controls from the single equation ((ref)) or ((ref)). In finite-sample experiments, single-selection methods essentially fail, providing poor inference relative to the double-selection method outlined above. This performance is also supported theoretically by the fact that the double-selection method requires weaker regularity conditions for its validity and for attaining the efficiency bound\footnote{Semi-parametric efficiency is attained in the homoscedastic case.} than the single selection method.

Now we formally define the post-double-selection estimator: Let $\widehat I_1 = {\rm support }(\widehat \beta_1)$ denote the control terms selected by a feasible Lasso estimator $\widehat \beta_1$ computed using data $(\tilde y_i,\tilde x_i) = (d_{i}, x_i), \ i =1,...,n$. Let $\widehat I_2 = {\rm support }(\widehat \beta_2)$ denote the control terms selected by a feasible Lasso estimator $\widehat \beta_2$ computed using data $(\tilde y_i,\tilde x_i) = (y_{i}, x_i), \ i =1,...,n$. The post-double-selection estimator $\check \alpha$ of $\alpha_0$ is defined as the least squares estimator obtained by regressing $y_{i}$ on $d_i$ and the selected control terms $x_{ij}$ with $j \in \widehat I \supseteq \widehat I_1 \cup \widehat I_2$:

equation[equation omitted — 240 chars of source]

The set $\widehat I$ may contain variables that were not selected in the variable selection steps with indices in $\widehat I_3$ that the analyst thinks are important for ensuring robustness. We call $\widehat I_3$ the amelioration set. Thus, $\widehat I = \widehat I_1 \cup \widehat I_2 \cup \widehat I_3$; let $\widehat s = |\widehat I|$ and $\widehat s_j = |\widehat I_j|$ for $j =1,2,3$.

We define feasible Lasso estimators below and note that other selection methods could be used as well. When a feasible Lasso is used to construct $\widehat I_1$ and $\widehat I_2$, we refer to the post-double-selection estimator as the post-double-Lasso estimator. When other model selection devices are used to construct $\widehat I = \widehat I_1 \text{ and } \widehat I_2$, we shall refer the estimator as the generic post-double-selection estimator.

The main theoretical result of the paper shows that the post-double-selection estimator $\check \alpha$ obeys

equation[equation omitted — 195 chars of source]

under approximate sparsity conditions, uniformly within a rich set of data generating processes. We also show that the standard plug-in estimator for standard errors is consistent in these settings. All of these results imply uniform validity of confidence regions over large, interesting classes of models. Figure (ref) (right panel) illustrates the result ((ref)) by showing that the finite-sample distribution of our post-double-selection estimator is very close to the normal distribution. In contrast, Figure (ref) (left panel) illustrates the classical problem with the traditional post-single-selection estimator based on ((ref)), showing that its distribution is bimodal and sharply deviates from the normal distribution. Finally, it is worth noting that the estimator achieves the semi-parametric efficiency bound under homoscedasticity.

figure[figure omitted — 526 chars of source]

Selection of controls via feasible Lasso Methods

Here we describe feasible variable selection via Lasso. Note that each of the regression equations above is of the form $$ \tilde y_i = \underbrace{\tilde x_i'\beta_0 + r_i}_{f(\tilde z_i)} + \epsilon_i, $$ where $f(\tilde z_i)$ is the regression function, $\tilde x_i'\beta_0$ is the approximation based on the dictionary $\tilde x_i=P(\tilde z_i)$, $r_i$ is the approximation error, and $\epsilon_i$ is the error. The Lasso estimator is defined as a solution to

equation[equation omitted — 153 chars of source]

where $\|\beta\|_{1} = \sum_{j=1}^p | \beta_j|$; see FF:1993 and T1996. The kinked nature of the penalty function induces the solution $\widehat \beta$ to have many zeroes, and thus the Lasso solution may be used for model selection. The selected model $\widehat T = \text{support}(\widehat \beta)$ is often used for further refitting by least squares, leading to the so called post-Lasso or Gauss-Lasso estimator, see, e.g., \citeasnoun{BC-PostLASSO}. The Lasso estimator/selector is computationally attractive because it minimizes a convex function. In the homoskedastic Gaussian case, a basic choice for penalty level suggested by \citeasnoun{BickelRitovTsybakov2009} is

equation[equation omitted — 101 chars of source]

where $c>1$, $1-\gamma$ is a confidence level that needs to be set close to 1, and $\sigma$ is the standard deviation of the noise. The formal motivation for this penalty is that it leads to near-optimal rates of convergence of the estimator under approximate sparsity. The good behavior of the estimator of $\beta_0$ in turn implies good approximation properties of the selected model $\widehat T$, as noted in \citeasnoun{BC-PostLASSO}. Unfortunately, even in the homoskedastic case the penalty level specified above is not feasible since it depends on the unknown $\sigma$.

\citeasnoun{BellChenChernHans:nonGauss} formulate a feasible Lasso estimator/selector $\widehat \beta$ geared for heteroscedastic, non-Gaussian cases, which solves

equation[equation omitted — 165 chars of source]

where $\widehat \Psi = {\rm diag}(\widehat l_1,\ldots,\widehat l_p)$ is a diagonal matrix of penalty loadings. The penalty level $\lambda$ and loadings $\widehat l_j$'s are set as

equation[equation omitted — 243 chars of source]

where $c>1$ and $1-\gamma$ is a confidence level.\footnote{Practical recommendations include the choice $c=1.1$ and $\gamma=.05$.} The $\l_j$'s are ideal penalty loadings that are not observed, and we estimate $\l_j$ by $\widehat \l_j$ obtained via an iteration method given in Appendix A. We refer to the resulting feasible Lasso method as the Iterated Lasso. The estimator $\widehat \beta$ has statistical performance that is similar to that of the (infeasible) Lasso described above in Gaussian cases and delivers similar performance in non-Gaussian, heteroscedastic cases; see \citeasnoun{BellChenChernHans:nonGauss}. In this paper, we only use $\widehat \beta$ as a model selection device. Specifically, we only make use of $$ \widehat T = \text{support}(\widehat \beta), $$ the labels of the regressors with non-zero estimated coefficients. We show that the selected model $\widehat T$ has good approximation properties for the regression function $f$ under approximate sparsity in Section 3.

\citeasnoun{BCW-SqLASSO} propose another feasible variant of Lasso called the Square-root Lasso estimator, $\widehat \beta$, defined as a solution to

equation[equation omitted — 177 chars of source]

with the penalty level

equation[equation omitted — 93 chars of source]

where $c>1$, $\gamma \in (0,1)$ is a confidence level, and $\widehat \Psi = {\rm diag}(\widehat l_1,\ldots,\widehat l_p)$ is a diagonal matrix of penalty loadings. The main attractive feature of ((ref)) is that one can set $\widehat \l_j =\{{\mathbb{E}_n}[\tilde x_{ij}^2]\}^{1/2}$ which depends only on observed data in the homoscedastic case.

In the heteroscedastic case, we would like to choose $\widehat l_j$ so that

equation[equation omitted — 238 chars of source]

As a simple bound, we could use $\widehat \l_j = 2 \{{\mathbb{E}_n}[\tilde x_{ij}^4]\}^{1/4}$ since $$\{{\mathbb{E}_n}[\tilde x_{ij}^2\epsilon_i^2]]/{\mathbb{E}_n}[\epsilon_i^2]\}^{1/2} \leqslant \{{\mathbb{E}_n}[\tilde x_{ij}^4]\}^{1/4}\{{\mathbb{E}_n}[\epsilon_i^4]\}^{1/4}/\{{\mathbb{E}_n}[\epsilon_i^2]\}^{1/2}.$$ This bound gives $l_j + o_P(1) \leqslant \widehat \l_j$ if $\{{\mathbb{E}_n}[\epsilon_i^4]\}^{1/4}/\{{\mathbb{E}_n}[\epsilon_i^2]\}^{1/2} \leqslant 2 + o_P(1)$, which covers a wide class of marginal distributions for error $\epsilon_i$. For example, all $t$-distributions with degrees of freedom greater than five satisfy this condition. As in the previous case, we can also iteratively re-estimate the penalty loadings using estimates of the $\epsilon_i$'s to approximate the ideal penalty loadings:

equation[equation omitted — 114 chars of source]

The resulting Square-root Lasso and post-Square-root Lasso estimators based on these penalty loadings achieve near optimal rates of convergence even in non-Gaussian, heteroscedastic cases. This good performance implies good approximation properties for the selected model $\widehat T$.

In what follows, we shall use the term feasible Lasso to refer to either the Iterated Lasso estimator $\widehat \beta $ solving ((ref))-((ref)) or the Square-root Lasso estimator $\widehat \beta$ solving ((ref))-((ref)) with $c > 1$ and $1-\gamma$ set such that

equation[equation omitted — 107 chars of source]

Theory of Estimation and Inference

Regularity Conditions

In this section, we provide regularity conditions that are sufficient for validity of the main estimation and inference result. We begin by stating our main condition, which contains the previously defined approximate sparsity as well as other more technical assumptions. Throughout the paper, we let $c$, $C$, and $q$ be absolute constants, and let $\ell_n \nearrow \infty, \delta_n \searrow 0$, and $\Delta_n \searrow 0$ be sequences of absolute positive constants. By absolute constants, we mean constants that are given, and do not depend the dgp ${\mathrm{P}}= {\mathrm{P}}_n$.

We assume that for each $n$ the following condition holds on dgp ${\mathrm{P}} = {\mathrm{P}}_n$.

Condition ASTE (${\mathrm{P}}$). (i) $\{(y_{i}, d_i, z_i), i = 1,...,n\}$ are i.n.i.d. vectors on $(\Omega, \mathcal{F}, {\mathrm{P}})$ that obey the model ((ref))-((ref)), and the vector $x_i = P(z_i)$ is a dictionary of transformations of $z_i$, which may depend on $n$ but not on ${\mathrm{P}}$. (ii) The true parameter value $\alpha_0$, which may depend on ${\mathrm{P}}$, is bounded, $|\alpha_0| \leqslant C$. (iii) Functions $m$ and $g$ admit an approximately sparse form. Namely there exists $s \geqslant 1$ and $\beta_{m0}$ and $\beta_{g0}$, which depend on $n$ and ${\mathrm{P}}$, such that

eqnarray[eqnarray omitted — 309 chars of source]

(iv) The sparsity index obeys $s^2 \log^2 (p\vee n)/n \leqslant \delta_n$ and the size of the amelioration set obeys $ \widehat s_3 \leqslant C (1\vee \widehat s_1 \vee \widehat s_2)$. (v) For $\tilde v_i = v_i + r_{mi}$ and $\tilde \zeta_i = \zeta_i + r_{gi}$ we have $|\bar {\mathrm{E}}[ \tilde v_i^2\tilde \zeta_i^2 ] - \bar {\mathrm{E}}[ v_i^2\zeta_i^2 ]| \leqslant \delta_n$, and $\bar {\mathrm{E}}[|\tilde v_i|^q+|\tilde \zeta_i|^q] \leqslant C$ for some $q>4$. Moreover, $\max_{i\leqslant n} \| x_{i}\|^2_\infty s n^{-1/2+2/q} \leqslant \delta_n$ wp $1-\Delta_n$.}

remarkThe approximate sparsity (iii) and the growth condition (iv) are the main conditions for establishing the key inferential result. We present a number of primitive examples to show that these conditions contain standard models used in empirical research as well as more flexible models. Condition (iv) requires that the size $\widehat s_3$ of the amelioration set $\widehat I_3$ should not be substantially larger than the size of the set of variables selected by the Lasso method. Simply put, if we decide to include controls in addition to those selected by Lasso, the total number of additions should not dominate the number of controls selected by Lasso. This and other conditions will ensure that the total number $\widehat s$ of controls obeys $ \widehat s \lesssim_P s$, and we also require that $s^2 \log^2 (p\vee n)/n \to 0$. This condition can be relaxed using the sample-splitting method of \citeasnoun{FanGuoHao2011}, which is done in the Supplementary Appendix. Condition (v) is simply a set of sufficient conditions for consistent estimation of the variance of the double selection estimator. If the regressors are uniformly bounded and the approximation errors are going to zero a.s., it is implied by other conditions stated below; and it can also be demonstrated under other sorts of more primitive conditions. \qed

The next condition concerns the behavior of the Gram matrix ${\mathbb{E}_n} [x_ix_i']$. Whenever $p>n$, the empirical Gram matrix ${\mathbb{E}_n}[x_ix_i']$ does not have full rank and in principle is not well-behaved. However, we only need good behavior of smaller submatrices. Define the minimal and maximal $m$-sparse eigenvalue of a semi-definite matrix $M$ as

equation[equation omitted — 288 chars of source]

To assume that $\phi_{{\rm min}}(m)[{\mathbb{E}_n} [x_ix_i']] >0$ requires that all empirical Gram submatrices formed by any $m$ components of $x_i$ are positive definite. We shall employ the following condition as a sufficient condition for our results.

Condition SE (${\mathrm{P}}$). There is an absolute sequence of constants $\ell_n \to \infty$ such that the maximal and minimal $\ell_n s$-sparse eigenvalues are bounded from below and away from zero, namely with probability at least $1-\Delta_n$, $$\kappa' \leqslant \phi_{{\rm min}}(\ell_n s)[{\mathbb{E}_n} [x_ix_i']] \leqslant \phi_{{\rm max}}(\ell_n s)[{\mathbb{E}_n} [x_ix_i']] \leqslant \kappa'',$$ where $0< \kappa' < \kappa'' < \infty$ are absolute constants.

remarkIt is well-known that Condition SE is quite plausible for many designs of interest. For instance, Condition SE holds if \begin{itemize} • $ x_i$, $i = 1,\ldots,n$, are i.i.d. zero-mean sub-Gaussian random vectors that have population Gram matrix ${\mathrm{E}}[ x_i x_i']$ with minimal and maximal $s\log n$-sparse eigenvalues bounded away from zero and from above by absolute constants where $s(\log n )(\log p)/n \leqslant \delta_n \to 0$; • $ x_i$, $i=1,\ldots,n$, are i.i.d. bounded zero-mean random vectors with $\| x_i\|_\infty \leqslant K_n$ a.s. that have population Gram matrix ${\mathrm{E}}[ x_i x_i']$ with minimal and maximal $s\log n$-sparse eigenvalues bounded from above and away from zero by absolute constants where $K_n^2s(\log^3 n)\{\log(p\vee n)\} /n \leqslant \delta_n \to 0$. \end{itemize} The claim (a) holds by Theorem 3.2 in \citeasnoun{RudelsonZhou2011} (see also \citeasnoun{Zhou2009a} and \citeasnoun{Baraniuketal2008}) and claim (b) holds by Lemma 1 in \citeasnoun{BC-PostLASSO} or by Theorem 1.8 \citeasnoun{RudelsonZhou2011}. Recall that a standard assumption in econometric research is to assume that the population Gram matrix ${\mathrm{E}}[x_i x_i']$ has eigenvalues bounded from above and away from zero, see e.g. \citeasnoun{newey:series}. The conditions above allow for this and more general behavior, requiring only that the $s \log n$ sparse eigenvalues of the population Gram matrix ${\mathrm{E}}[x_i x_i']$ are bounded from below and from above. \qed

The next condition imposes moment conditions on the structural errors and regressors.

Condition SM (${\mathrm{P}}$). There are absolute constants $0< c< C < \infty$ and $4< q < \infty$ such that for $(\tilde y_i, \epsilon_i) = (y_i, \zeta_i) $ and $(\tilde y_i, \epsilon_i) = (d_i, v_i)$ the following conditions hold:

itemize$\displaystyle \bar {\mathrm{E}} [|d_i|^q] \leqslant C, \ \ c \leqslant {\mathrm{E}}[\zeta_i^2\mid x_i, v_i]\leqslant C \ \mbox{and} \ c\leqslant {\mathrm{E}}[v_i^2\mid x_i]\leqslant C \ \mbox{a.s. } 1 \leqslant i \leqslant n$, • $\displaystyle \bar {\mathrm{E}}[|\epsilon_i|^q]+\bar {\mathrm{E}} [\tilde y_i^2] + \max_{1\leqslant j\leqslant p} \{ \bar {\mathrm{E}}[x_{ij}^2{\tilde y}_i^2]+\bar {\mathrm{E}}[|x_{ij}^3 \epsilon_i^3|]+ 1/\bar {\mathrm{E}}[ x_{ij}^2] \} \leqslant C$, • $ \displaystyle \log^3 p / n \leqslant \delta_n$, • $ \displaystyle \max_{1\leqslant j\leqslant p} \{ |({\mathbb{E}_n}-\bar {\mathrm{E}})[ x_{ij}^2\epsilon_i^2]|+|({\mathbb{E}_n}-\bar {\mathrm{E}})[x_{ij}^2{\tilde y}_i^2]|\} + \max_{1\leqslant i\leqslant n}\| x_i\|_\infty^2 \frac{s\log(n\vee p)}{n} \leqslant \delta_n \text{ wp } 1-\Delta_n.$

} These conditions, which are rather mild, ensure good model selection performance of feasible Lasso applied to equations ((ref)) and ((ref)). These conditions also allow us to invoke moderate deviation theorems for self-normalized sums from \citeasnoun{jing:etal} to bound some important error components.

The Main Result

The following is the main result of this paper. It shows that the post-double selection estimator is root-$n$ consistent and asymptotically normal. Under homoscedasticity this estimator achieves the semi-parametric efficiency bound. The result also verifies that plug-in estimates of the standard errors are consistent.

theorem[Estimation and Inference on Treatment Effects] Let $\{{\mathrm{P}}_n\}$ be a sequence of data-generating processes. Assume conditions ASTE (${\mathrm{P}}$), SM (${\mathrm{P}}$), and SE (${\mathrm{P}}$) hold for ${\mathrm{P}} = {\mathrm{P}}_n$ for each $n$. Then, the post-double-Lasso estimator $\check \alpha$, constructed in the previous section, obeys as $n \to \infty$ $$ \sigma_n^{-1} \sqrt{n} (\check \alpha - \alpha_0) \rightsquigarrow N(0,1), $$ where $\sigma^2_n= [\bar {\mathrm{E}} v_i^2]^{-1}\bar {\mathrm{E}}[ v_i^2\zeta_i^2] [\bar {\mathrm{E}} v_i^2]^{-1}$. Moreover, the result continues to apply if $\sigma^2_n$ is replaced by $\widehat \sigma^2_n = [{\mathbb{E}_n} \widehat v_i^2]^{-1}{\mathbb{E}_n}[\widehat v_i^2\widehat \zeta_i^2][{\mathbb{E}_n} \widehat v_i^2]^{-1}$, for $\widehat \zeta_i := [y_i - d_i\check \alpha - x_i'\check \beta]\{n/(n - \widehat s-1)\}^{1/2}$ and $\widehat v_i:=d_i - x_i'\widehat\beta$, $i=1,\ldots,n$ where $\widehat \beta \in \arg\min_\beta \{{\mathbb{E}_n}[(d_i-x_i'\beta)^2]:\beta_j=0, \forall j\notin \widehat I\}$.

A consequence of this result is the following corollary.

corollary[Uniformly Valid Confidence Intervals] (i) Let $\mathbf{P}_n$ be the collection of all data-generating processes ${\mathrm{P}}$ for which conditions ASTE(${\mathrm{P}}$), SM (${\mathrm{P}}$), and SE (${\mathrm{P}}$) hold for given $n$. Let $c(1-\xi) = \Phi^{-1} (1-\xi/2)$. Then as $n \to \infty$, uniformly in ${\mathrm{P}} \in \mathbf{P}_n$ $$ {\mathrm{P}} \left ( \alpha_0 \in [ \check \alpha \pm c(1-\xi) \widehat \sigma_n /\sqrt{n}]\right) \to 1- \xi . $$ (ii) Let $\mathbf{P} = \cap_{n \geqslant n_0} \mathbf{P}_n$ be the collection of data-generating processes for which the conditions above hold for all $n \geqslant n_0$ for some $n_0$. Then as $n \to \infty$, uniformly in ${\mathrm{P}} \in \mathbf{P}$ $$ {\mathrm{P}} \left ( \alpha_0 \in [ \check \alpha \pm c(1-\xi) \widehat \sigma_n /\sqrt{n}]\right) \to 1- \xi. $$

By exploiting both equations ((ref)) and ((ref)) for model selection, the post-double-selection method creates the necessary adaptivity that makes it robust to imperfect model selection. Robustness of the post-double selection method is reflected in the fact that Theorem (ref) permits the data-generating process to change with $n$. Thus, the conclusions of the theorem are valid for a wide variety of sequences of data-generating processes which in turn define the regions $\mathbf{P}$ of uniform validity of the resulting confidence sets. These regions appear to be substantial, as we demonstrate via a sequence of theoretical and numerical examples in Section 5 and 6. In contrast, the standard post-selection method based on ((ref)) generates non-robust confidence intervals.

remarkOur approach to uniformity analysis is most similar to that of \citeasnoun{romano:uniform}, Theorem 4. It proceeds under triangular array asymptotics, with the sequence of dgps obeying certain constraints; then these results imply uniformity over sets of dgps that obey the constraints for all sample sizes. This approach is also similar to the classical central limit theorems for sample means under triangular arrays, and does not require the dgps to be parametrically (or otherwise tightly) specified, which then translates into uniformity of confidence regions. This approach is somewhat different in spirit to the generic uniformity analysis suggested by \citeasnoun{andrews:cheng:guggen}. \qed
remarkUniformity holds over a large class of approximately sparse models, which cover conventional models used in series estimation of partially linear models as shown in Section 5. Of course, for every interesting class of models and any inference method, one could find an even bigger class of models where the uniformity does not apply. In particular, our models do not cover models with many small coefficients. In the series case, a model with many small coefficients corresponds to a deviation from smoothness towards highly non-smooth functions, namely functions generated as realized paths of an approximate white noise process. The fact that our results do not cover such models motivates further research work on inference procedures that have robustness properties to deviations from the given class of models that are deemed important. In the simulations in Section 6, we consider incorporating the ridge fit along the other controls to be selected over using lasso to build extra robustness against “many small coefficients" deviations away from approximately sparse models. \qed

Auxiliary Results on Model Selection via Lasso and Post-Lasso

The post-double-selection estimator applies the least squares estimator to the union of variables selected for equations ((ref)) and ((ref)) via feasible Lasso. Therefore, the model selection properties of feasible Lasso as well as properties of least squares estimates for $m$ and $g$ based on the selected model play an important role in the derivation of the main result. The purpose of this section is to describe these properties. The proof of Theorem 1 relies on these properties.

Note that each of the regression models ((ref))-((ref)) obeys the following conditions.

Condition ASM. Let $\{{\mathrm{P}}_n\}$ be a sequence of data-generating processes. For each $n$, we have data $\{(\tilde y_i,\tilde z_i,\tilde x_i=P(\tilde z_i)) : 1 \leqslant i \leqslant n\}$ defined on $(\Omega, \mathcal{A}, {\mathrm{P}}_n) $ consisting of i.n.i.d vectors that obey the following approximately sparse regression model for each $n$:

eqnarray*[eqnarray* omitted — 290 chars of source]

} Let $\widehat T$ denote the model selected by the feasible Lasso estimator $\widehat\beta$: $$\widehat T = {\rm support}( \widehat \beta ) = \{ j \in \{1,\ldots,p\} \ : \ |\widehat\beta_j| > 0\},$$ The Post-Lasso estimator $\widetilde \beta$ is is ordinary least squares applied to the data after removing the regressors that were not selected by the feasible Lasso:

equation[equation omitted — 200 chars of source]

The following regularity conditions are imposed to deal with non-Gaussian, heteroscedastic errors.

Condition RF. In addition to ASTE, we have

itemize$\log^{3} p / n \to 0 \ \text{ and } \ s \log (p\vee n) /n \to 0$, • $ \bar {\mathrm{E}}[\tilde y_i^2] + \max_{1\leqslant j\leqslant p} \{\bar {\mathrm{E}}[\tilde x_{ij}^2\tilde y_i^2]+\bar {\mathrm{E}}[|\tilde x_{ij}^3 \epsilon_i^3|]+ 1/\bar {\mathrm{E}}[\tilde x_{ij}^2\epsilon_i^2]\} \lesssim 1$, • $\displaystyle \max_{1\leqslant j\leqslant p} \{|({\mathbb{E}_n}-\bar {\mathrm{E}})[\tilde x_{ij}^2\epsilon_i^2]|+|({\mathbb{E}_n}-\bar {\mathrm{E}})[\tilde x_{ij}^2\tilde y_i^2]|\} + \max_{1\leqslant i\leqslant n}\|\tilde x_i\|_\infty^2 \frac{s\log(n\vee p)}{n} = o_P(1)$.

}

The main auxiliary result that we use in proving the main result is as follows.

lemma[Model Selection Properties of Lasso and Properties of Post-Lasso] Let $\{{\mathrm{P}}_n\}$ be a sequence of data-generating processes. Suppose that conditions ASM and RF hold, and that Condition SE $({\mathrm{P}}_n)$ holds for ${\mathbb{E}_n}[\tilde x_i \tilde x_i']$. Consider a feasible Lasso estimator with penalty level and loadings specified as in Section 3.3. (i) Then the data-dependent model $\widehat T$ selected by a feasible Lasso estimator satisfies with probability approaching 1: \begin{equation} \widehat s = | \widehat T | \lesssim s \end{equation} and \begin{equation} \min_{\beta \in \Bbb{R}^p: \ \beta_j = 0 \ \forall j \not \in \widehat T} \sqrt{{\mathbb{E}_n}[ f(\tilde z_i) - \tilde x_i'\beta]^2} \lesssim \sigma \sqrt{ \frac{ s \log (p \vee n)}{n} }. \end{equation} (ii) The Post-Lasso estimator obeys $$ \sqrt{{\mathbb{E}_n}[ f (\tilde z_i) - \tilde x_i'\widetilde \beta]^2} \lesssim_P \sigma \sqrt{ \frac{ s \log (p \vee n)}{n} }. $$ and \begin{equation} \|\widetilde \beta - \beta_0\| \lesssim_P \sqrt{{\mathbb{E}_n}[\{\tilde x_i'\widetilde \beta - \tilde x_i'\beta_0\}^2]} \lesssim_P \sigma \sqrt{ \frac{ s \log (p \vee n)}{n} }. \end{equation}

Lemma (ref) was derived in \citeasnoun{BellChenChernHans:nonGauss} for Iterated Lasso and by \citeasnoun{BCW-SqLASSO2} for Square-root Lasso. These analyses build on the rate analysis of infeasible Lasso by \citeasnoun{BickelRitovTsybakov2009} and on sparsity analysis and rate analysis of Post-Lasso by \citeasnoun{BC-PostLASSO}. Lemma (ref) shows that feasible Lasso methods select a model $\widehat T$ that provides a high-quality approximation to the regression function $f(\tilde z_i)$; i.e. they find a sparse model that can approximate the function at the “near-oracle" rate $\sqrt{s/n} \sqrt{\log (p\vee n)}$. If we knew the “best" approximating model $ T= {\rm support}(\beta_0)$, we could achieve the “oracle" rate of $\sqrt{s/n}$. Note that Lasso methods generally will not recover $T$ perfectly. Moreover, no method can recover $T$ perfectly in general, except under the restrictive condition that all non-zero coefficients in $\beta_0$ are bounded away from zero by a factor that exceeds estimation error. We do not require this condition to hold in our results. All that we need is that the selected model $\widehat T$ can approximate the regression function well and that the size of the selected model, $\widehat s = |\widehat T|$, is of the same stochastic order as $s = |T|$. This condition holds in many cases in which some non-zero coefficients are close to zero.

The lemma above also shows that feasible Post-Lasso achieves the same near-oracle rate as feasible Lasso. The coincidence in rates occurs despite the fact that feasible Lasso will in general fail to correctly select the best-approximating model $T$ as a subset of the variables selected; that is, $T \not \subseteq \widehat T$. The intuition for this result is that any components of $T$ that feasible Lasso misses are unlikely to be important; otherwise, ((ref)) would be impossible. This result was first derived in the context of median regression by \citeasnoun{BC-SparseQR} and extended to least squares in reference cited above.

Generalization: Inference after Double Selection by a Generic Selection Method

The conditions provided so far are simply a set sufficient conditions that are tied to the use of Lasso as the model selector. The purpose of this section is to prove that the main results apply to any other model selection method that is able to select a sparse model with good approximation properties. As in the case of Lasso, we allow for imperfect model selection. Next we state a high-level condition that summarizes a sufficient condition on the performance of a model selection method that allows the post-double selection estimator to attain good inferential properties.

Condition HLMS (${\mathrm{P}}$). A model selector provides possibly data-dependent sets $ \widehat I_1 \cup \widehat I_2 \subseteq \widehat I \subset \{1,...,p\}$ of covariate names such that, with probability $1-\Delta_n$, $|\widehat I| \leqslant C s$ and $$\displaystyle \min_{\beta: \beta_j=0, j\not\in \widehat I_1} \sqrt{{\mathbb{E}_n}[ (m(z_i)-x_i'\beta)^2]} \leqslant \delta_n n^{-1/4} \textrm{ and } \min_{\beta: \beta_j=0, j\not\in \widehat I_2} \sqrt{{\mathbb{E}_n}[ (g(z_i)-x_i'\beta)^2]} \leqslant \delta_n n^{-1/4}.$$

Condition HLMS requires that with high probability the selected models are sparse and generates a good approximation for the functions $g$ and $m$. Examples of methods producing such models include the Dantzig selector CandesTao2007, feasible Dantzig selector gautier:tsybakov, Bridge estimator HHS2008, SCAD penalized least squares FanLi2001, and thresholded Lasso BC-PostLASSO, to name a few. We emphasize that, similarly to the previous arguments, we allow for imperfect model selection.

The following result establishes the inferential properties of a generic post-double-selection estimator.

theorem[Estimation and Inference on Treatment Effects under High-Level Model Selection] Let $\{{\mathrm{P}}_n\}$ be a sequence of data-generating processes and the model selection device be such that conditions ASTE (${\mathrm{P}}$), SM (${\mathrm{P}}$), SE (${\mathrm{P}}$), and HLSM(${\mathrm{P}}$) hold for ${\mathrm{P}} = {\mathrm{P}}_n$ for each $n$. Then the generic post-double-selection estimator $\check \alpha$ based on $\widehat I$, as defined in ((ref)), obeys $$ ([\bar {\mathrm{E}} v_i^2]^{-1}\bar {\mathrm{E}}[ v_i^2\zeta_i^2] [\bar {\mathrm{E}} v_i^2]^{-1})^{-1/2} \sqrt{n} (\check \alpha - \alpha_0) \rightsquigarrow N(0,1). $$ Moreover, the result continues to apply if $\bar {\mathrm{E}}[v_i^2]$ and $\bar {\mathrm{E}}[v_i^2\zeta_i^2]$ are replaced by ${\mathbb{E}_n}[\widehat v_i^2]$ and ${\mathbb{E}_n}[\widehat v_i^2\widehat \zeta_i^2]$ for $\widehat \zeta_i := [y_i - d_i\check \alpha - x_i'\check \beta]\{n/(n - \widehat s-1)\}^{1/2}$ and $\widehat v_i:=d_i - x_i'\widehat\beta$, $i=1,\ldots,n$ where $\widehat \beta \in \arg\min_\beta \{{\mathbb{E}_n}[(d_i-x_i'\beta)^2]:\beta_j=0, \forall j\notin \widehat I\}$.

Theorem (ref) can also be used to establish uniformly valid confidence intervals as shown is the following corollary.

corollary[Uniformly Valid Confidence Intervals] (i) Let $\mathbf{P}_n$ be the collection of all data-generating processes ${\mathrm{P}}$ for which conditions ASTE(${\mathrm{P}}$), SM (${\mathrm{P}}$), SE (${\mathrm{P}}$), and HLSM (${\mathrm{P}}$) hold for given $n$. Let $c(1-\xi) = \Phi^{-1} (1-\xi/2)$. Then as $n \to \infty$, uniformly in ${\mathrm{P}} \in \mathbf{P}_n$ $$ {\mathrm{P}} \left ( \alpha_0 \in [ \check \alpha \pm c(1-\xi) \widehat \sigma_n /\sqrt{n}]\right) \to 1- \xi . $$ (ii) Let $\mathbf{P} = \cap_{n \geqslant n_0} \mathbf{P}_n$ be the collection of data-generating processes for which the conditions above hold for all $n \geqslant n_0$ for some $n_0$. Then as $n \to \infty$, uniformly in ${\mathrm{P}} \in \mathbf{P}$ $$ {\mathrm{P}} \left ( \alpha_0 \in [ \check \alpha \pm c(1-\xi) \widehat \sigma_n /\sqrt{n}]\right) \to 1- \xi. $$

Theoretical Examples

The purpose of this section is to give a sequence of examples -- progressing from simple to somewhat involved -- that highlight the range of the applicability and robustness of the proposed method. In these examples, we specify primitive conditions which cover a broad range of applications including nonparametric models and high-dimensional parametric models. We emphasize that our main regularity conditions cover even more general models which combine various features of these examples such as models with both nonparametric and high-dimensional parametric components.

In all examples, the model is

equation[equation omitted — 205 chars of source]

however, the structure for $g$ and $m$ will vary across examples, and so will the assumptions on the error terms $\zeta_i$ and $v_i$.

We start out with a simple example, in which the dimension $p$ of the regressors is fixed. In practical terms this example approximates cases with $p$ small compared to $n$. This simple example is important since standard post-single-selection methods fail even in this simple case. Specifically, they produce confidence intervals that are not valid uniformly in the underlying data-generating process; see \citeasnoun{leeb:potscher:pms}. In contrast, the post-double-selection method produces confidence intervals that are valid uniformly in the underlying data-generating process.

Example 1. (Parametric Model with Fixed $p$.) Consider $(\Omega, \mathcal{A}, {\mathrm{P}})$ as the probability space, on which we have $(y_i, z_i, d_i)$ as i.i.d. vectors for $i=1,...,n$ obeying the model ((ref)) with

equation[equation omitted — 142 chars of source]

For estimation we use $x_i = (z_{ij}, j=1,...,p)'$. We assume that there are some absolute constants $0< b < B<\infty$, $q_x \geqslant q > 4$, with $4/q_x + 4/q < 1$, such that

equation[equation omitted — 428 chars of source]

Let $\mathbf{P}$ be the collection of all regression models ${\mathrm{P}}$ that obey the conditions set forth above for all $n$ for the given constants $(p, b, B, q_x, q)$. Then, as established in Appendix (ref), any ${\mathrm{P}} \in \mathbf{P}$ obeys Conditions ASTE (${\mathrm{P}}$) with $s=p$, SE (${\mathrm{P}}$), and SM (${\mathrm{P}}$) for all $n \geqslant n_0$, with the constants $n_0$ and $( \kappa', \kappa'', c, C)$ and sequences $\Delta_n$ and $\delta_n$ in those conditions depending only on $(p, b, B, q_x, q)$. Therefore, the conclusions of Theorem 1 hold for any sequence ${\mathrm{P}}_n \in \mathbf{P}$, and the conclusions of Corollary 1 on the uniform validity of confidence intervals apply uniformly in ${\mathrm{P}} \in \mathbf{P}$. \qed

The next examples are more substantial and include infinite-dimensional models which we approximate with linear functional forms with potentially very many regressors, $p \gg n$. The key to estimation in these models is a smoothness condition which requires regression coefficients to decay at some rates. In series estimation, this condition is often directly connected to smoothness of the regression function.

Let ${a}$ and $A$ be positive constants. We shall say that a sequence of coefficients $$\theta = \{ \theta_j, j=1,2,...\}$$ is ${a}$-smooth with constant $A$ if $$ | \theta_j| \leqslant A j^{-{a}}, \ j=1,2,... , $$ which will be denoted as $\theta \in S^{{a}}_A$. We shall say that a sequence of coefficients $\theta = \{ \theta_j, j=1,2,...\}$ is ${a}$-smooth with constant $A$ after $p$-rearrangement if $$ | \theta_{(j)}| \leqslant A j^{-{a}}, \ j=1,2,..., p, \ \ | \theta_j| \leqslant A j^{-{a}}, \ j=p+1, p+2,..., $$ which will be denoted as $\theta \in S^{{a}}_A(p)$, where $\{| \theta_{(j)}|, j=1,...,p\}$ denotes the decreasing rearrangement of the numbers $\{ | \theta_j|, j =1,...,p\}$. Since $S^{{a}}_A \subset S^{{a}}_A(p)$, the second kind of smoothness is strictly more general than the first kind.

Here we use the term “smoothness" motivated by Fourier series analysis where smoothness of functions often translates into smoothness of the Fourier coefficients in the sense that is stated above; see, e.g., \citeasnoun{kerk:picard}. For example, if a function $h: [0,1]^d \mapsto \Bbb{R}$ possesses $r>0$ continuous derivatives uniformly bounded by a constant $M$ and the terms $P_j$ are compactly supported Daubechies wavelets, then $h$ can be represented as $h(z) = \sum_{j=1}^\infty P_j(z) \theta_{hj}$, with $|\theta_{hj}| \leqslant A j^{-r/d-1/2}$ for some constant $A$; see \citeasnoun{kerk:picard}. We also note that the second kind of smoothness is considerably more general than the first since it allows relatively large coefficients to appear anywhere in the series of the first $p$ coefficients. In contrast, the first kind of smoothness only allows relatively large coefficients among the early terms in the series. Lasso-type methods are specifically designed to deal with the generalized smoothness of the second kind and perform equally well under both kinds of smoothness. In the context of series applications, smoothness of the second kind allows one to approximate functions that exhibit oscillatory phenomena or spikes, which are associated with “high order" series terms. An example of this is the wage function example given in \citeasnoun{BCH2011:InferenceGauss}.

Before we proceed to other examples we discuss a way to generate sparse approximations in infinite-dimensional examples. Consider, for example, a function $h$ that can be represented a.s. as $h(z_i) = \sum_{j=1}^{\infty}\theta_{hj}P_j(z_i) $ with coefficients $\theta_h \in S_A^{a}(p)$. In this case we can construct sparse approximations by simply thresholding to zero all coefficients smaller than $1/\sqrt{n}$ and with indices $j \geqslant p$. This generates a sparsity index $ s \leqslant A^{\frac{1}{a}} n^{\frac{1}{2a}}$. The non-zero coefficient could be further reoptimized by using the least squares projection. More formally, given a sparsity index $s>0$, a target function $h(z_i)$, and terms $x_i =( P_j(z_i) : j=1,\ldots,p)'\in {\Bbb{R}}^p$, we let

equation[equation omitted — 160 chars of source]

and define $x_i'\beta_{h0}$ as the best $s$-sparse approximation to $h(z_i)$.

Example 2. (Gaussian Model with Very Large $p$.) Consider $(\Omega, \mathcal{A}, {\mathrm{P}})$ as the probability space on which we have $(y_i, z_i, d_i)$ as i.i.d. vectors for $i=1,...,n$ obeying the model ((ref)) with

equation[equation omitted — 158 chars of source]

Assume that the infinite dimensional vector $w_i = (z_i', \zeta_i, v_i)'$ is jointly Gaussian with minimal and maximal eigenvalues of the matrix (operator) ${\mathrm{E}}[w_i w_i']$ bounded below by an absolute constant $\underline{\kappa}>0$ and above by an absolute constant $\overline{\kappa}< \infty$.

The main assumption that guarantees approximate sparsity is the smoothness condition on the coefficients. Let ${a}>1$ and $0<A<\infty$ be some absolute constants. We require that the coefficients of the expansions in ((ref)) are ${a}$-smooth with constant $A$ after $p$-rearrangement, namely $$ \theta_m = (\theta_{mj}, j=1,2,...) \in S^{{a}}_A(p), \ \ \theta_g = (\theta_{gj}, j=1,2,...) \in S^{{a}}_A(p). $$ For estimation purposes we shall use $x_i = (z_{ij}, j=1,...,p)',$ and assume that $\|\alpha_0\|\leqslant B$ and $p = p_n$ obeys $$ n^{[(1-{a})/{a}]+\chi}\log^2(p\vee n) \leqslant\bar\delta_n, \ \ A^{1/{a}}n^{\frac{1}{2{a}}} \leqslant p \bar\delta_n, \ \text{ and } \ \log^3 p /n \leqslant \bar \delta_n, $$ for some absolute sequence $\bar\delta_n \searrow 0$ and absolute constants $B$ and $\chi > 0$.

Let $\mathbf{P}_n$ be the collection of all dgp ${\mathrm{P}}$ that obey the conditions set forth in this example for a given $n$ and for the given constants $(\underline{\kappa}, \overline{\kappa}, {a}, A, B,\chi)$ and sequences $p=p_n$ and $\bar\delta_n$. Then, as established in Appendix (ref), any ${\mathrm{P}} \in \mathbf{P}_n$ obeys Conditions ASTE (${\mathrm{P}}$) with $s=A^{1/{a}} n^{\frac{1}{2{a}}}$, SE (${\mathrm{P}}$), and SM (${\mathrm{P}}$) for all $n \geqslant n_0$, with constants $n_0$ and $( \kappa', \kappa'', c, C)$ and sequences $\Delta_n$ and $\delta_n$ in those conditions depending only on $(\underline{\kappa}, \bar\kappa, {a}, A, B,\chi)$, $p$, and $\bar\delta_n$. Therefore, the conclusions of Theorem 1 hold for any sequence ${\mathrm{P}}_n \in \mathbf{P}_n$, and the conclusions of Corollary 1 on the uniform validity of confidence intervals apply uniformly for any ${\mathrm{P}} \in \mathbf{P}_n$. In particular, these conclusions apply uniformly in $ {\mathrm{P}} \in \mathbf{P} = \cap_{n \geqslant n_0} \mathbf{P}_n$. \qed

Example 3. (Series Model with Very Large $p$.) Consider $(\Omega, \mathcal{A}, {\mathrm{P}})$ as the probability space, on which we have $(y_i, z_i, d_i)$ as i.i.d. vectors for $i=1,..., n$ obeying the model:

equation[equation omitted — 165 chars of source]

where $z_i$ has support $[0,1]^d$ with density bounded from below by constant $\underline{f}>0$ and above by constant $\bar f$, and $\{ P_j, j =1,2,..\}$ is an orthonormal basis on $L^2[0,1]^d$ with bounded elements, i.e. $ \max_{z \in [0,1]^d}|P_j(z)| \leqslant B$ for all $j =1,2,...$. Here all constants are taken to be absolute. Examples of such orthonormal bases include canonical trigonometric bases.

Let ${a}>1$ and $0<A<\infty$ be some absolute constants. We require that the coefficients of the expansions in ((ref)) are ${a}$-smooth with constant $A$ after $p$-rearrangement, namely $$ \theta_m = (\theta_{mj}, j=1,2,...) \in S^{{a}}_A(p), \ \ \theta_g = (\theta_{gj}, j=1,2,...) \in S^{{a}}_A(p). $$

For estimation purposes we shall use $x_i = (P_j(z_{i}), j=1,...,p)',$ and assume that $p = p_n$ obeys $$ n^{(1-{a})/{a}}\log^2(p\vee n) \leqslant\bar\delta_n, \ \ A^{1/{a}}n^{\frac{1}{2{a}}} \leqslant p \bar\delta_n \ \text{ and } \ \log^3 p /n \leqslant \bar\delta_n, $$ for some sequence of absolute constants $\bar\delta_n \searrow 0$. We assume that there are some absolute constants $b>0$, $B<\infty$, $q > 4$, with $(1-{a})/{a}+ 4/q <0$, such that

equation[equation omitted — 283 chars of source]

Let $\mathbf{P}_n$ be the collection of all regression models ${\mathrm{P}}$ that obey the conditions set forth above for a given $n$. Then, as established in Appendix (ref), any ${\mathrm{P}} \in \mathbf{P}_n$ obeys Conditions ASTE (${\mathrm{P}}$) with $s=A^{1/{a}} n^{\frac{1}{2{a}}}$, SE (${\mathrm{P}}$), and SM (${\mathrm{P}}$) for all $n \geqslant n_0$, with absolute constants in those conditions depending only on $(\underline{f}, \bar f, {a}, A, b, B, q)$ and $\bar\delta_n$. Therefore, the conclusions of Theorem 1 hold for any sequence ${\mathrm{P}}_n \in \mathbf{P}_n$, and the conclusions of Corollary 1 on the uniform validity of confidence intervals apply uniformly for any ${\mathrm{P}} \in \mathbf{P}_n$. In particular, as a special case, the same conclusion applies uniformly in $ {\mathrm{P}} \in \mathbf{P} = \cap_{n \geqslant n_0} \mathbf{P}_n$. \qed

Monte-Carlo Examples

In this section, we examine the finite-sample properties of the post- double-selection method through a series of simulation exercises and compare its performance to that the standard post-single-selection method.

All of the simulation results are based on the structural model

equation[equation omitted — 125 chars of source]

where $p = \dim(x_i) = 200$, the covariates $ x_i \sim N(0,\Sigma)$ with $\Sigma_{kj} = (0.5)^{|j-k|}$, $\alpha_0 = .5$, and the sample size $n$ is set to $100$. In each design, we generate

equation[equation omitted — 99 chars of source]

with E[$\zeta_i v_i] = 0 $. Inference results for all designs are based on conventional t-tests with standard errors calculated using the heteroscedasticity consistent jackknife variance estimator discussed in \citeasnoun{mackinnon:white}. Another option would be to use the standard error estimator recently proposed in \citeasnoun{CJN:PLMStandardError}.

We report results from three different dgp's. In the first two dgp's, we set $\theta_{g,j} = c_y\beta_{0,j}$ and $\theta_{m,j} = c_d\beta_{0,j}$ with $\beta_{0,j} = (1/j)^2$ for $j = 1,...,200$. The first dgp, which we label “Design 1,” uses homoscedastic innovations with $\sigma_y = \sigma_d = 1$. The second dgp, “Design 2,” is heteroscedastic with $\sigma_{d,i} = \sqrt{\frac{(1+x_i'\beta_0)^2}{{\mathbb{E}_n}(1+x_i'\beta_0)^2}}$ and $\sigma_{y,i} = \sqrt{\frac{(1+\alpha_0 d_i + x_i'\beta_0)^2}{{\mathbb{E}_n}(1+\alpha_0 d_i+x_i'\beta_0)^2}}$. The constants $c_y$ and $c_d$ are chosen to generate desired population values for the reduced form $R^2$'s, i.e. the $R^2$'s for equations ((ref)) and ((ref)). For each equation, we choose $c_y$ and $c_d$ to generate $R^2 = 0, .2, .4, .6,$ and $.8$. In the heteroscedastic design, we choose $c_y$ and $c_d$ based on $R^2$ as if ((ref)) and ((ref)) held with $v_i$ and $\zeta_i$ homoscedastic and label the results by $R^2$ as in Design 1. In the third design (“Design 3”), we use a combination of deterministic and random coefficients. For the deterministic coefficients, we set $\theta_{g,j} = c_y(1/j)^2$ for $j \le 5$ and $\theta_{m,j} = c_d (1/j)^2$ for $j \le 5$. We then generate the remaining coefficients as iid draws from $(\theta_{g,j},\theta_{m,j})' \sim N(0_{2 \times 1},(1/p) I_2)$. For each equation, we choose $c_y$ and $c_d$ to generate $R^2 = 0, .2, .4, .6,$ and $.8$ in the case that all of the random coefficients were exactly equal to 0 and label the results by $R^2$ as in Design 1. We draw new $x$'s, $\zeta$'s, and $v$'s at every simulation replication, and we also generate new $\theta$'s at every simulation replication in Design 3.

We consider Designs 1 and 2 to be baseline designs. These designs do not have exact sparse representations but have coefficients that decay quickly so that approximately sparse representations are available. Design 3 is meant to introduce a modest deviation from the approximately sparse model towards a model with many small, uncorrelated coefficients. Using this we shall document that our proposed procedure still performs reasonably well, although it could be improved by incorporation of a ridge fit as one of regressors over which selection occurs. In a working paper version of this paper \citeasnoun{BCH2011:SuppMatTE}, we present results for 26 additional designs. The results presented in this section are sufficient to illustrate the general patterns from the larger set of results.\footnote{ In particular, the post-double-Lasso performed very well across all simulations designs where approximate sparsity provides a reasonable description of the dgp. Unsurprisingly, the performance deteriorates as one deviates from the smooth/approximately sparse case. However, in no design was the post-double-Lasso outperformed by other feasible procedures. In extensive initial simulations, we also found that Square-Root Lasso and Iterated Lasso performed very similarly and thus only report Lasso results. }

We report results for five different procedures. Two of the procedures are infeasible benchmarks: Oracle and Double-Selection Oracle estimators, which use of knowledge of the true coefficient structures $\theta_g$ and $\theta_m$ and are thus unavailable in practice. The Oracle estimator is the ordinary least squares of $y_i - x_i'\theta_g$ on $d_i$, and the Double-Selection Oracle is the ordinary least squares of $y - x_i'\theta_g$ on $d_i - x_i'\theta_m$. The other procedures we consider are feasible. In all of them, we rely on Lasso and set $\lambda$ according to the algorithm outlined in Appendix A with $1-\gamma = .95$. One procedure is the standard post-single selection estimator -- the Post-Lasso -- which applies Lasso to equation ((ref)) without penalizing $\alpha$, the coefficient on $d$, to select additional control variables from among $x$. Estimates of $\alpha_0$ are then obtained by OLS regression of $y$ on $d$ and the set of additional controls selected in the Lasso step and inference using the Post-Lasso estimator proceeds using conventional heteroscedasticity robust OLS inference from this regression. Post-Double-Selection or Post-Double-Lasso is the feasible procedure advocated in this paper. We run Lasso of $y$ on $x$ to select a set of predictors for $y$ and run Lasso of $d$ on $x$ to select a set of predictors for $d$. $\alpha_0$ is then estimated by running OLS regression of $y$ on $d$ and the union of the sets of regressors selected in the two Lasso runs, and inference is simply the usual heteroscedasticity robust OLS inference from this regression. Post-Double-Selection $+$ Ridge is an ad hoc variant of Post-Double-Selection in which we add the ridge fit from equation ((ref)) as an additional potential regressor that may be selected by Lasso. The ridge fit is obtained with a single ridge penalty parameter that is chosen using 10-fold cross-validation. This procedure is motivated by a desire to add further robustness in the case that many small coefficients are suspected. Further exploration of procedures that perform well, both theoretically and in simulations, in the presence of many small coefficients is an interesting avenue for additional research.

We start by summarizing results in Table 1 for $(R^2_y,R^2_d) = (0,.2), (0,.8), (.8,.2),$ and $(.8,.8)$ where $R^2_y$ is the population $R^2$ from regressing $y$ on $x$ (Structure $R^2$) and $R^2_d$ is the population $R^2$ from regressing $d$ on $x$ (First Stage $R^2$). We report root-mean-square-error (RMSE) for estimating $\alpha_0$ and size of 5% level tests (Rej. Rate). As should be the case, the Oracle and Double-Selection Oracle, which are reported to provide the performance of an infeasible benchmark, perform well relative to the feasible procedures across the three designs. We do see that the feasible Post-Double-Selection procedures perform similarly to the Double-Selection Oracle without relying on ex ante knowledge of the coefficients that go in to the control functions, $\theta_g$ and $\theta_m$. On the other hand, the Post-Lasso procedure generally does not perform as well as Post-Double-Selection and is very sensitive to the value of $R^2_d$. While Post-Lasso performs adequately when $R^2_d$ is small, its performance deteriorates quickly as $R^2_d$ increases. This lack of robustness of traditional variable selection methods such as Lasso which were designed with forecasting, not inference about treatment effects, in mind is the chief motivation for our advocating the Post-Double-Selection procedure when trying to infer structural or treatment parameters.

We provide further details about the performance of the feasible estimators in Figures 1, 2, and 3 which plot size of 5% level tests, bias, and standard deviation for the Post-Lasso, Double-Selection (DS), and Double-Selection Oracle (DS Oracle) estimators of the treatment effect across the full set of $R^2$ values considered. Figure 1, 2, and 3 respectively report the results from Design 1, 2, and 3. The figures are plotted with the same scale to aid comparability and for readability rejection frequencies for Post-Lasso were censored at .5. Perhaps the most striking feature of the figures is the poor performance of the Post-Lasso estimator. The Post-Lasso estimator performs poorly in terms of size of tests across many different $R^2$ combinations and can have an order of magnitude more bias than the corresponding Post-Double-Selection estimator. The behavior of Post-Lasso is quite non-uniform across $R^2$ combinations, and Post-Lasso does not reliably control size distortions or bias except in the case where the controls are uncorrelated with the treatment (where First-Stage $R^2$ equals 0) and thus ignorable. In contrast, the Post-Double-Selection estimator performs relatively well across the full range of $R^2$ combinations considered. The Post-Double-Selection estimator's performance is also quite similar to that of the infeasible Double-Selection Oracle across the majority of $R^2$ values considered. Comparing across Figures 1 and 2, we see that size distortions for both the Post-Double-Selection estimator and the Double-Selection Oracle are somewhat larger in the presence of heteroscedasticity but that the basic patterns are more-or-less the same across the two figures. Looking at Figure 3, we also see that the addition of small independent random coefficients results in somewhat larger size distortions for the Post-Double-Selection estimator than in the other homoscedastic design, Design 1, though the procedure still performs relatively well.

In the final figure, Figure 4, we compare the performance of the Post-Double-Selection procedure to the ad hoc Post-Double-Selection procedure which selects among the original set of variables augmented with the ridge fit obtained from equation ((ref)). We see that the addition of this variable does add robustness relative to Post-Double-Selection using only the raw controls in the sense of producing tests that tend to have size closer to the nominal level. This additional robustness is a good feature, though it comes at the cost of increased RMSE which is especially prominent for small values of the first-stage $R^2$.

The simulation results are favorable to the Post-Double-Selection estimator. In the simulations, we see that the Post-Double-Selection procedure provides an estimator of a treatment effect in the presence of a large number of potential confounding variables that performs similarly to the infeasible estimator that knows the values of the coefficients on all of the confounding variables. Overall, the simulation evidence supports our theoretical results and suggests that the proposed Post-Double-Selection procedure can be a useful tool to researchers doing structural estimation in the presence of many potential confounding variables. It also shows, as a contrast, that the standard Post-Single-Selection procedure provides poor inference and therefore can not be a reliable tool to these researchers.

Empirical Example: Estimating the Effect of Abortion on Crime

In the preceding sections, we have provided results demonstrating how variable selection methods, focusing on the case of Lasso-based methods, can be used to estimate treatment effects in models in which we believe the variable of interest is exogenous conditional on observables. We further illustrate the use of these methods in this section by reexamining Donohue III and Levitt's levitt:abortion study of the impact of abortion on crime rates. In the following, we briefly review \citeasnoun{levitt:abortion} and then present estimates obtained using the methods developed in this paper.

\citeasnoun{levitt:abortion} discuss two key arguments for a causal channel relating abortion to crime. The first is simply that more abortion among a cohort results in an otherwise smaller cohort and so crime 15 to 25 years later, when this cohort is in the period when its members are most at risk for committing crimes, will be otherwise lower given the smaller cohort size. The second argument is that abortion gives women more control over the timing of their fertility allowing them to more easily assure that childbirth occurs at a time when a more favorable environment is available during a child's life. For example, access to abortion may make it easier to ensure that a child is born at a time when the family environment is stable, the mother is more well-educated, or household income is stable. This second channel would mean that more access to abortion could lead to lower crime rates even if fertility rates remained constant.

The basic problem in estimating the causal impact of abortion on crime is that state-level abortion rates are not randomly assigned, and it seems likely that there will be factors that are associated to both abortion rates and crime rates. It is clear that any association between the current abortion rate and the current crime rate is likely to be spurious. However, even if one looks at say the relationship between the abortion rate 18 years in the past and the crime rate among current 18 year olds, the lack of random assignment makes establishing a causal link difficult without adequate controls. An obvious confounding factor is the existence of persistent state-to-state differences in policies, attitudes, and demographics that are likely related to the overall state level abortion and crime rates. It is also important to control flexibly for aggregate trends. For example, it could be the case that national crime rates were falling over this period while national abortion rates were rising but that these trends were driven by completely different factors. Without controlling for these trends, one would mistakenly associate the reduction in crime to the increase in abortion. In addition to these overall differences across states and times, there are other time varying characteristics such as state-level income, policing, or drug-use to name a few that could be associated with current crime and past abortion.

To address these confounds, \citeasnoun{levitt:abortion} estimate a model for state-level crime rates running from 1985 to 1997 in which they condition on a number of these factors. Their basic specification is

align[align omitted — 112 chars of source]

where $i$ indexes states, $t$ indexes times, $c \in \{\textnormal{violent, property, murder}\}$ indexes type of crime, $\delta_i$ are state-specific effects that control for any time-invariant state-specific characteristics, $\gamma_t$ are time-specific effects that control flexibly for any aggregate trends, $w_{it}$ are a set of control variables to control for time-varying confounding state-level factors, $a_{cit}$ is a measure of the abortion rate relevant for type of crime $c$,\footnote{This variable is constructed as weighted average of abortion rates where weights are determined by the fraction of the type of crime committed by various age groups. For example, if 60% of violent crime were committed by 18 year olds and 40% were committed by 19 year olds in state $i$, the abortion rate for violent crime at time $t$ in state $i$ would be constructed as .6 times the abortion rate in state $i$ at time $t-18$ plus .4 times the abortion rate in state $i$ at time $t-19$. See \citeasnoun{levitt:abortion} for further detail and exact construction methods.} and $y_{cit}$ is the crime-rate for crime type $c$. \citeasnoun{levitt:abortion} use the log of lagged prisoners per capita, the log of lagged police per capita, the unemployment rate, per-capita income, the poverty rate, AFDC generosity at time $t - 15$, a dummy for concealed weapons law, and beer consumption per capita for $w_{it}$, the set of time-varying state-specific controls. Tables IV and V in \citeasnoun{levitt:abortion} present baseline estimation results based on ((ref)) as well as results from different models which vary the sample and set of controls to show that the baseline estimates are robust to small deviations from ((ref)). We refer the reader to the original paper for additional details, data definitions, and institutional background.

For our analysis, we take the argument that the abortion rates defined above may be taken as exogenous relative to crime rates once observables have been conditioned on from \citeasnoun{levitt:abortion} as given. Given the seemingly obvious importance of controlling for state and time effects, we account for these in all models we estimate. We choose to eliminate the state effects via differencing rather than including a full set of state dummies but include a full set of time dummies in every model. Thus, we will estimate models of the form

align[align omitted — 123 chars of source]

We use the same state-level data as \citeasnoun{levitt:abortion} but delete Alaska, Hawaii, and Washington, D.C. which gives a sample with 48 cross-sectional observations and 12 time series observations for a total of 576 observations. With these deletions, our baseline estimates using the same controls as in ((ref)) are quite similar to those reported in \citeasnoun{levitt:abortion}. Baseline estimates from Table IV of \citeasnoun{levitt:abortion} and our baseline estimates based on the differenced version of ((ref)) are given in the first and second row of Table 2 respectively.

Our main point of departure from \citeasnoun{levitt:abortion} is that we allow for a much richer set $z_{it}$ than allowed for in $w_{it}$ in model ((ref)). Our $z_{it}$ includes higher-order terms and interactions of the control variables defined above. In addition, we put initial conditions and initial differences of $w_{it}$ and $a_{it}$ into our vector of controls $z_{it}$. This addition allows for the possibility that there may be some feature of a state that is associated both with its growth rate in abortion and its growth rate in crime. For example, having an initially high-levels of abortion could be associated with having high-growth rates in abortion and low growth rates in crime. Failure to control for this factor could then lead to misattributing the effect of this initial factor, perhaps driven by policy or state-level demographics, to the effect of abortion. Finally, we allow for more general trends by allowing for an aggregate quadratic trend in $z_{it}$ as well as interactions of this quadratic trend with control variables. This gives us a set of 251 control variables to select among in addition to the 12 time effects that we include in every model.\footnote{The exact identities of the 251 potential controls is available upon request. It consists of linear and quadratic terms of each continuous variable in $w_{it}$, interactions of every variable in $w_{it}$, initial levels and initial differences of $w_{it}$ and $a_{it}$, and interactions of these variables with a quadratic trend.}

Note that interpreting estimates of the effect of abortion from model ((ref)) as causal relies on the belief that there are no higher-order terms of the control variables, no interaction terms, and no additional excluded variables that are associated both to crime rates and the associated abortion rate. Thus, controlling for a large set of variables as described above is desirable from the standpoint of making this belief more plausible. At the same time, naively controlling lessens our ability to identify the effect of interest and thus tends to make estimates far less precise. The effect of estimating the abortion effect conditioning on the full set of 251 potential controls described above is given in the third row of Table 2. As expected, all coefficients are estimated very imprecisely. Of course, very few researchers would consider using 251 controls with only 576 observations due to exactly this issue.

We are faced with a tradeoff between controlling for very few variables which may leave us wondering whether we have included sufficient controls for the exogeneity of the treatment and controlling for so many variables that we are essentially mechanically unable to learn about the effect of the treatment. The variable selection methods developed in this paper offer one resolution to this tension. The assumed sparse structure maintains that there is a small enough set of variables that one could potentially learn about the treatment but adds substantial flexibility to the usual case where a researcher considers only a few control variables by allowing this set to be found by the data from among a large set of controls. Thus, the approach should complement the usual careful specification analysis by providing a researcher an efficient, data-driven way to search for a small set of influential confounds from among a sensibly chosen broad set of potential confounding variables.

In the abortion example, we use the post-double-selection estimator defined in Section (ref) for each of our dependent variables. For violent crime, ten variables are selected in the abortion equation,\footnote{The selected variables are AFDC generosity squared, beer consumption squared, the initial poverty change, initial income, initial income squared, the initial change in prisoners per capita squared interacted with the trend, initial income interacted with the trend, the initial change in the abortion rate, the initial change in the abortion rate interacted with the trend, and the initial level of the abortion rate.} and one is selected in the crime equation.\footnote{The initial level of the abortion rate interacted with time is selected.} For property crime, eight variables are selected in the abortion equation,\footnote{The selected variables are income, the initial poverty change, the initial change in prisoners per capita squared, the initial level of prisoners per capita, initial income, the initial change in the abortion rate, the initial change in the abortion rate interacted with the trend, and the initial level of the abortion rate.} and six are selected in the crime equation.\footnote{The six variables are the initial level of AFDF generosity, the initial level of income interacted with the trend and the trend squared, the initial level of income squared interacted with the trend and the trend squared, and the initial level of the abortion rate interacted with the trend.} For murder, eight variables are selected in the abortion equation,\footnote{The selected variables are AFDC generosity, beer consumption squared, the change in beer consumption squared, the change in beer consumption squared times the trend and the trend squared, initial income times the trend, the initial change in the abortion rate interacted with the trend, and the initial level of the abortion rate.} and none were selected in the crime equation.

Estimates of the causal effect of abortion on crime obtained by searching for confounding factors among our set of 251 potential controls are given in the fourth row of Table 2. Each of these estimates is obtained from the least squares regression of the crime rate on the abortion rate and the 11, 14, and eight controls selected by the double-post-Lasso procedure for violent crime, property crime, and murder respectively. The estimates for the effect of abortion on violent crime and the effect of abortion on murder are quite imprecise, producing 95% confidence intervals that encompass large positive and negative values. The estimated effect for property crime is roughly in line with the previous estimates though it is no longer significant at the 5% level but is significant at the 10% level. Note that the double-post-Lasso produces models that are not of vastly different size than the “intuitive” model ((ref)). As a final check, we also report results that include all of the original variables from ((ref)) in the amelioration set in the fifth row of the table. These results show that the conclusions made from using only the variable selection procedure do not qualitatively change when the variables used in the original \citeasnoun{levitt:abortion} are added to the equation. For a quick benchmark relative to the simulation examples, we note that the $R^2$ obtained by regressing the crime rate on the selected variables are .0395, .1185, and .0044 for violent crime, property crime, and the murder rate respectively and that the $R^2$'s from regressing the abortion rate on the selected variables are .9447, .9013, and .9144 for violent crime, property crime, and the murder rate respectively. These values correspond to regions of the $R^2$ space considered in the simulation where the double selection procedure substantially outperformed simple Lasso procedures.

It is very interesting that one would draw qualitatively different conclusions from the estimates obtained using formal variable selection than from the estimates obtained using a small set of intuitively selected controls. Looking at the set of selected control variables, we see that initial conditions and interactions with trends are selected across all dependent variables. The selection of this set of variables suggests that there are initial factors which are associated with the change in the abortion rate. We also see that we cannot precisely determine the effect of the abortion rate on crime rates once one accounts for initial conditions. Of course, this does not mean that the effects of the abortion rate provided in the first two rows of Table 2 are not representative of the true causal effects. It does, however, imply that this conclusion is strongly predicated on the belief that there are not other unobserved state-level factors that are correlated to both initial values of the controls and abortion rates, abortion rate changes, and crime rate changes. Interestingly, a similar conclusion is given in \citeasnoun{FooteGoetzAbortion} based on an intuitive argument.

We believe that the example in this section illustrates how one may use modern variable selection techniques to complement causal analysis in economics. In the abortion example, we are able to search among a large set of controls and transformations of variables when trying to estimate the effect of abortion on crime. Considering a large set of controls makes the underlying assumption of exogeneity of the abortion rate conditional on observables more plausible, while the methods we develop allow us to produce an end-model which is of manageable dimension. Interestingly, we see that one would draw quite different conclusions from the estimates obtained using formal variable selection. Looking at the variables selected, we can also see that this change in interpretation is being driven by the variable selection method's selecting different variables, specifically initial values of the abortion rate and controls, than are usually considered. Thus, it appears that the usual interpretation hinges on the prior belief that initial values should be excluded from the structural equation.

Conclusion

In this paper, we consider estimation of treatment effects or structural parameters in an environment where the treatment is believed to be exogenous conditional on observables. We do not impose the conventional assumption that the identities of the relevant conditioning variables and the functional form with which they enter the model are known. Rather, we assume that the researcher believes there is a relatively small number of important factors whose identities are unknown within a much larger known set of potential variables and transformations. This sparsity assumption allows the researcher to estimate the desired treatment effect and infer a set of important variables upon which one needs to condition by using modern variable selection techniques without ex ante knowledge of which are the important conditioning variables. Since naive application of variable selection methods in this context may result in very poor properties for inferring the treatment effect of interest, we propose a “double-selection” estimator of the treatment effect, provide a formal demonstration of its properties for estimating the treatment effect, and provide its approximate distribution under technical regularity conditions and the assumed sparsity in the model.

In addition to the theoretical development, we illustrate the potential usefulness of our proposal through a number of simulation studies and an empirical example. In Monte Carlo simulations, our procedure outperforms simple variable selection strategies for estimating the treatment effect across the designs considered and does relatively well compared to an infeasible estimator that uses the identities of the relevant conditioning variables. We then apply our estimator to attempt to estimate the causal impact of abortion on crime following \citeasnoun{levitt:abortion}. We find that our procedure selects a small number of conditioning variables. After conditioning on these selected variables, one would draw qualitatively different inference about the effect of abortion on crime than would be drawn if one assumed that the correct set of conditioning variables was known and the same as those variables used in \citeasnoun{levitt:abortion}. Taken together, the empirical and simulation examples demonstrate that the proposed method may provide a useful complement to other sorts of specification analysis done in applied research.