arXiv:1309.4686v3 [] 1 Feb 2018

Robust Inference on Average Treatment Effects with Possibly More
Covariates than Observations∗
Max H. Farrell†
University of Chicago Booth School of Business
June 9, 2015
Updated: February 1, 20181

Abstract
This paper concerns robust inference on average treatment effects following model selection.
Under selection on observables, we construct confidence intervals using a doubly-robust estimator that are robust to model selection errors and prove their uniform validity over a large class of
models that allows for multivalued treatments with heterogeneous effects and selection amongst
(possibly) more covariates than observations. The semiparametric efficiency bound is attained
under appropriate conditions. Precise conditions are given for any model selector to yield these
results, and we specifically propose the group lasso, which is apt for treatment effects, and derive new results for high-dimensional, sparse multinomial logistic regression. Both a simulation
study and revisiting the National Supported Work demonstration show our estimator performs
well in finite samples.

Keywords: High-dimensional sparse model, heterogeneous treatment effects, uniform inference, model selection, doubly-robust estimator, unconfoundedness, group lasso.
JEL Classification: C21, C31, C52.

∗
An online suppplement contains additional proofs and simulations. I am deeply grateful to Matias Cattaneo
for advice and support. I am indebted to Xuming He, Lutz Kilian, and Jeffrey Smith for thoughtful feedback and
discussions. I thank Victor Chernozhukov for pointing to the relevant latest results, obtained in joint work Alexandre
Belloni and Christian Hansen, and the latter two authors for conversations in the early stages of this project. I
benefited from discussions with Rosa Matzkin, Blaise Melly, and Jack Porter. I also thank the co-editor, Han Hong,
and two reviewers for their detailed comments and suggestions that improved the paper.
†
Correspondence to: University of Chicago Booth School of Business, 5807 South Woodlawn Avenue, Chicago,
IL 60637, United States. Tel: +1 773-834-0161; Email: max.farrell@chicagobooth.edu; Web: http://faculty.
chicagobooth.edu/max.farrell/.
1
The published version, Farrell (2015a), contains an error in the proof [present also in Farrell (2015b)] which is
rectified here. Assumption 3(c) is required for the results of Section 5 to be valid; see the author’s website for further
detail. I am grateful to Whitney Newey for alerting me to this error.

1

Introduction

Model selection has always had a place in empirical economics, whether or not it is formally acknowledged. A key problem in modern empirical work is that researchers face datasets with large
numbers of variables, sometimes more than observations. A complementary problem is that economic theory and prior knowledge may mandate controlling for certain variables, but are generally
silent regarding functional form. These two problems force researchers to search for a model that is
simultaneously parsimonious and adequately flexible. Many formal methods are computationally
infeasible with a large number of variables. A typical response to this challenge is to iteratively
search over a small set of alternative specifications, guided only by the researcher’s taste and intuition. But no matter the approach used, subsequent inference almost never takes accounts for
this “specification search” and the resulting confidence intervals are not robust to model selection
mistakes, and hence are unreliable in empirical work.
This problem is particularly important in estimating average treatment effects under selection
on observables, because in this framework using the right covariates is crucial for identification
and correct inference. In this context, we provide an easy-to-implement and objective method for
covariate selection and post-selection inference on average treatment effects.2 We establish four
main results for multivalued treatments effects with arbitrary heterogeneity in observables and
heteroskedasticity. First, we show that a doubly-robust estimator is robust to model selection
errors. These estimators were initially developed for robustness to parametric misspecification, but
are now known to be robust to selection.3 By taking explicit account of the model selection stage
and its inherent selection errors, we derive precise conditions required for any model selector to
deliver confidence intervals for average treatment effects that are uniformly valid over a large class
of data-generating processes. Second, we show that a simple refitting procedure allows researchers
to augment variables chosen according economic theory with data-driven selection to deliver flexible
inference that remains uniformly valid. Third, we prove that our estimator is asymptotically linear,
and standard conditions imposed in the program evaluation literature, semiparametrically efficient
bound. Fourth, we derive new results for multinomial (and binary) logistic regression, the most
widely used model for treatment assignment.
Inference following model selection is notoriously difficult. In a sequence of papers, Leeb and
Pötscher (2005, 2008a, 2008b, 2009) have shown that inference relying too heavily on model selection
can not be made uniformly valid. Loosely speaking, uniform validity of a confidence interval
captures the idea that the interval should have the same quality (coverage) for many data-generating
processes. This theoretical property is practically important because it implies greater reliability
2
Treatment effects, missing data, measurement error, and data combination models are equivalent under selection
on observables. Thus, all our results immediately apply to those contexts. For reviews of these literatures, see Tsiatis
(2006), Heckman and Vytlacil (2007), Imbens and Wooldridge (2009), and Wooldridge (2010).
3
Doubly-robust estimation and its role in program evaluation is discussed by Robins and Rotnitzky (1995), van der
Laan and Robins (2003), Kang and Schafer (2007, with discussion), Tan (2010), and references therein.

1

in applications. Our proposed methods for post model selection inference build upon the pathbreaking recent work of Belloni, Chernozhukov, and Hansen (2014).
The crucial insight that leads to uniform inference is to change the goal of model selection
away from perfect covariate selection (the oracle property) and to high-quality approximation
of the underlying functions. This fundamental shift in focus allows us to circumvent, without
contradicting, the impossibility results of Leeb and Pötscher. Valid post-selection inference has
attracted considerable attention during the preparation of this paper: in contexts and with methods
quite different from ours, contributions have been made by Belloni, Chernozhukov, and Wei (2013),
Berk, Brown, Buja, Zhang, and Zhao (2013), Zhang and Zhang (2014), Efron (2014), van de Geer,
Buhlmann, Ritov, and Dezeure (2014), and Belloni, Chernozhukov, Fernandez-Val, and Hansen
(2014), among others.
Our approach, based on the doubly-robust estimator, has several key features. The name
“doubly-robust” reflects that it is robust to misspecification of either the treatment equation
(propensity score) or the outcome equation, a property obtained by combining inverse probability weighting and regression imputation. First, we show that this robustness extends to model
selection, enabling us to allow for selection errors in both equations without impacting inference.
Second, we capture arbitrary treatment effect heterogeneity (dependence of the effect on an individual’s observed characteristics), which is crucial in empirical work. With such heterogeneity,
the average treatment effect and the treatment on the treated differ, and hence we present results
for both. Third, the doubly-robust estimator also stems from the semiparametric efficient moment
conditions, and hence we obtain the semiparametric efficiency bound, even under heteroskedasticity, under standard additional conditions. Thus, Pötscher’s (2009) result that sparse estimators
have large confidence sets is also circumvented. Taking all these features together enables us to
obtain uniform inference over such a large class of treatment effects models.
In recent independent work, Belloni, Chernozhukov, and Hansen (2014), propose a similar
approach. Their main focus is inference on the linear part of a partially linear model, which
motivates an estimator quite different from ours, but it will recover the average treatment effect in
the special case of a binary treatment where the effect is constant across observables. However, their
Section 5, developed independently from our work, considers heterogeneous effects and proposes an
estimator based on the efficient influence function, similar to ours. There are two broad differences.
First, we allow for multivalued treatments, which offers a larger set of estimands and can thus
enhance the understanding of program impacts.4 In this context we propose a group lasso based
approach that naturally exploits the already-present structure of treatment effects data to improve
model selection by pooling information across treatment levels. This is particularly natural in the
multivalued case, but even in the binary case there is still a grouped structure in the outcome
4

Discussion and applications may be found in, for example Imbens (2000), Lechner (2001), Imai and van Dyk
(2004), Abadie (2005), Cattaneo (2010), and Cattaneo and Farrell (2011).

2

regressions, though not in treatment assignment (i.e., in propensity score estimation). Second,
although in both cases the doubly-robust estimator is used for average treatment effects5 (following
a quite different model selection step), we show that this estimator has two benefits: (i) it may
require weaker conditions on the first stage (see Assumption 3); and (ii) it does not require using
variables selected for the treatment equation in the outcome model, and vice versa (“post double
selection”), and indeed, doing may require additional assumptions (see Assumption 5).
Our analysis is conducted under selection on observables, which has a long tradition and remains
quite popular in empirical economics.6 Covariates play three crucial roles in this framework. First,
using more observed covariates as proxies, and more flexibly, may help account for unobserved
confounding and hence increase the plausibility of unconfoundedness. Second, some observed variables may not be part of the causal mechanism under study, and should be excluded. Third, the
efficient conditioning set are those variables that drive the outcome, not necessarily those important for treatment assignment. This reasoning mandates contradicting goals for practitioners: a
large, rich set of controls on the one hand, and parsimony on the other. Our approach is a formal,
theory-driven attempt to reconcile this contradiction.
A special feature of our analysis is that we match the empirical realities of large data sets by
considering selection from amongst (possibly) more covariates than observations, so-called highdimensional data. The goal of variable selection is to find a small model that is nonetheless
sufficiently flexible to capture unknown features of the data-generating process required for inference. If a small model can perfectly capture the unknown feature it is said to be exactly sparse.
More realistic is approximate sparsity, when the bias from using a small model is well-controlled,
but nonzero. Sparsity is a natural framework for thinking about model selection. Indeed, any time
only a few of the available variables are used, a sparsity assumption has effectively been made. It is
common empirical practice to report results from several small models, but for these results to be
valid one must assume these specifications give high-quality, sparse representations of the unknown
features. The alternative we provide involves selecting a sparse, yet flexible, model from among a
large set of variables. Results may then be compared with more traditional methods.
With the aim of mimicking common empirical practice we estimate the propensity score with
multinomial logistic regression, coupled with group lasso selection (Yuan and Lin 2006). Our
results are stated in the language of treatment effects, but apply to general data structures and are
of independent interest in the high-dimensional literature.7 Much of the literature has focused on
5
They use different asymptotic variance estimators, and for treatment effects on the treated they do not exploit
the simplification discussed in Remark 1.
6
For other approaches and reviews of the literature, see, e.g., Holland (1986), Hahn (1998), Horowitz and Manski
(2000), Chen, Hong, and Tarozzi (2004, 2008), Bang and Robins (2005), Abadie and Imbens (2006), Wooldridge
(2007), and references therein.
7
Our techniques build on prior studies, in particular Bickel, Ritov, and Tsybakov (2009), Lounici, Pontil, van de
Geer, and Tsybakov (2011), Obozinski, Wainwright, and Jordan (2011), Belloni and Chernozhukov (2011), Belloni,
Chen, Chernozhukov, and Hansen (2012), Belloni and Chernozhukov (2013), and Belloni, Chernozhukov, FernandezVal, and Hansen (2014).

3

linear models (see Buhlmann and van de Geer (2011) for a survey), while prior studies of nonlinear
models often assume exact sparsity or present limited results.8 Furthermore, these studies often
use high-level conditions that can be hard to verify. In contrast, we obtain sharp results for logistic
regression under the same simple and intuitive conditions used for linear modeling by exploiting
mathematical techniques of self-concordant functions put forth by Bach (2010). We also provide
extensions to prior work on linear models needed to apply them in treatment effect estimation.
Finally, we offer numerical evidence on the finite sample performance of our procedure. In a
small simulation study we find that our procedure delivers very accurate coverage of confidence
intervals even for models where covariate selection is difficult, either because of a low signal-tonoise ratio or lack of sparsity, thus highlighting the uniform validity of inference. We also apply
our method to the widely-used National Supported Work Demonstration data (LaLonde 1986) and
find very accurate estimates and tight confidence intervals (see Table 1).
The paper proceeds as follows. Section 2 gives short, self-contained overview. Section 2.3
collects notation. Section 3 describes the treatment effect models. Sparse models are discussed in
Section 4, which shows how several commonly used models fit in this framework. Section 5 presents
our estimation method and complete results on treatment effect inference. Theoretical results for
the group lasso are in Section 6. Section 7 presents the numerical evidence and Section 8 concludes.
The main proofs are presented in the Appendix, while the remainder are available in a supplement.

2

Overview of Results and Notation

Here we give an overview of the paper, including treatment effect inference (Section 2.1), our new
results for the group lasso (Section 2.2), and notation used throughout (Section 2.3).

2.1

Treatment Effects and Results on Post-Selection Inference

We consider a multivalued treatment, with status indicated by D ∈ {0, 1, . . . , T}. Interest lies
in mean effects of the treatment on a scalar outcome Y . Let {Y (t)}Tt=0 be the (latent) potential
outcomes: Y (t) is the outcome a unit would have under D = t and is only observed for units with
P
D = t; that is, Y = Tt=0 1{D = t}Y (t). Many interesting parameters combine means of potential
outcomes, and having multivalued treatments allows for a wider range of estimands. Define the
mean of one potential outcome as µt = E[Y (t)]. To fix ideas, µ1 − µ0 is the average treatment effect
in the binary case (D ∈ {0, 1}). Sections 3 and 5 consider more general average effects, including
8

Examples include van de Geer (2008) and Negahban, Ravikumar, Wainwright, and Yu (2012), whose bounds do
not imply our results. Bach (2010) only gives an error bound on coefficients in exactly sparse logistic regression,
which can not yield our results; and does not consider prediction error or post-selection estimation. In independent
work, (Kwemou 2012) and Belloni, Chernozhukov, and Wei (2013) also apply Bach’s (2010) tools, but are focused
on a different goals. Vincent and Hansen (2014) apply the group lasso to multinomial logistic regression, but do not
derive any theoretical results.

4

effects on treated groups. For simplicity, in this section we focus on a single µt .
We use the selection on observables framework to identify µt . For a vector of covariates X,
define the generalized propensity score and conditional outcome regressions as
pt (x) = P[D = t|X = x]

and

µt (x) = E[Y |D = t, X = x].

For identification it is sufficient to assume that E[Y (t)|D, X] = E[Y (t)|X] (mean independence) and
pt (X) is bounded away from zero (overlap) for all treatment levels. Broadly, these two assumptions
imply that units from one treatment group are good proxies for other treatments and that there
are always such proxies available (see Section 3).
For an i.i.d. sample {(yi , di , x0i )}ni=1 and model-selection-based estimators p̂t (xi ) and µ̂t (xi ), we
estimate µt with

n

µ̂t =

1X
n
i=1




1{di = t}(yi − µ̂t (xi ))
+ µ̂t (xi ) .
p̂t (xi )

This doubly-robust estimator combines regression imputation and inverse probability weighting,
and remains consistent if either the model pt (x) or µt (x) is misspecified. Following widespread
empirical practice, we estimate p̂t (xi ) with multinomial logistic regression and µ̂t (xi ) linearly (see
Section 6). The choice of covariates in p̂t (xi ) and µ̂t (xi ) impacts consistency, efficiency, and finite
sample performance. Covariate selection based on ad hoc, iterative searches is common in empirical
work, but is not formal, objective, or replicable. Balancing tests are also common, but have the
additional drawback of assuming the same covariates are important for outcomes and treatment
assignment, and more generally do not weight the covariates by their importance for bias.
On the other hand, our proposed procedure gives practitioners an easy to implement, fully
objective tool to perform data-driven covariate selection and treatment effect inference, with replicable results.9 Importantly, we do not preclude the addition of variables known to be important
from economic theory or prior knowledge. Our procedure is intended to supplement these variables
with a flexible set of controls, guarding against misspecification or overfitting.
The following theorem is an example of the more general results presented in Section 5.2,
wherein we also define Vt and V̂t .
Theorem 1. Consider a sequence {Pn } of data-generating processes that obey, for each n, Assumptions 1 and 2 below. If the first stage obeys
(i)

Pn

(ii)

Pn

i=1 (p̂t (xi )
i=1 1{di

− pt (xi ))2 /n = oPn (1) and

Pn

i=1 (µ̂t (xi )

− µt (xi ))2 /n = oPn (1),

1/2 Pn
1/2
2
= oPn (n−1/2 ),
= t}(p̂t (xi ) − pt (xi ))2 /n
i=1 1{di = t}(µ̂t (xi ) − µt (xi )) /n

and
9

For the final step, the doubly-robust estimator is available in STATA and the package of Cattaneo, Drukker, and
Holland (2013). The covariate selection stage is easily implemented in R.

5

(iii)
then

Pn

i=1 (µ̂t (xi )

√


− µt (xi ))(1 − dti /pt (xi ))/n = oPn (n−1/2 ),

n(µ̂t − µt ) →d N (0, Vt ) and V̂t /Vt →Pn 1. For each n, let Pn be the set of data-generating

processes obeying Assumptions 1 and 2 and (i) and (ii) above. Then for cα = Φ−1 (1 − α/2)

sup PP µt ∈




q
µ̂t ± cα V̂t /n
− (1 − α) → 0.

P ∈Pn

This result establishes the uniform validity of an asymptotic confidence interval for µt , overcoming all the post model selection inference challenges: robustness to model selection errors, selecting
a model that is small but flexible enough to capture the features of the underlying data generating
process, and still retaining efficiency under standard conditions (see Section 5.3). Intuitively, this
is similar to (but distinct from) overcoming pretesting bias in other contexts. Also, although our
discussion is in terms of covariate selection in high-dimensional, sparse models, the inference result
is generic for any first stage estimator.
The two conditions placed on the first stage are analogous to the commonly-used, high-level
requirement in semiparametrics that first stage components converge faster than n−1/4 . However
exploiting features of the doubly-robust estimator yields weaker conditions. The first is a mild
consistency requirement. The second requires a rate on the product of errors and is thus easier to
satisfy if one function is easier to estimate, e.g. more smooth or more sparse. In model selection,
the rates for the first stage depend on the sample size, the number of covariates considered, and
the sparsity level. Importantly, the rate will depend on the total number of covariates only logarithmically, allowing for a large number. We propose to use the group lasso and prove that these
estimators satisfy (i) and (ii).

2.2

Model Selection Stage

We propose refitting following group lasso selection, and show that it meets all requirements on
the model selector. The group lasso is well-suited to program evaluation applications because
covariates are penalized according to their overall contribution in all treatment groups. This has
two consequences. First, information from all treatments is pooled when doing selection, and hence
a weaker signal may be extracted, which improves the selection properties. Second, the selected
variables are common to all treatment levels. From a practical point of view this is desirable, as
interest rarely lies in a single µt , but rather a collection, and substantial commonality is expected
in the variables important for different treatment levels.
We consider high-dimensional, sparse models for pt (x) and µt (x). These are defined by a pdimensional vector X ∗ based on the original variables X. The X ∗ may consist of any combination
of the original variables, interactions, flexible parametric transformations, and/or nonparametric
series terms (such as splines or polynomials). A model is approximately sparse if there are s < n

6

of these terms that yield a good approximation (s → ∞ is allowed). To build intuition, suppose
that µt (x) obeys a p-dimensional linear model. Then the sparsity assumption is that there is an sdimensional submodel with sufficiently small specification bias. In the nonparametric case, sparsity
is weaker than (but analogous to) the familiar assumption that a small set of basis functions can
approximate the unknown objects well. In practice researchers employ a hybrid of these approaches,
which is covered by our results. Section 4 gives more detail and examples.
We form p̂t (x) and µ̂t (x) in two steps (complete details in Section 6). First, the group lasso
is applied separately to multinomial logistic and least squares regression to select covariates from
X ∗ . We then estimate pt (x) and µt (x) by refitting unpenalized models using the selected variables,
possibly augmented with controls suggested by prior work or economic theory. It is not desirable
for a model selector to discard theory and prior work, and our procedure explicitly avoids this. We
also allow for using logistic-selected variables in the linear model refitting and vice versa, but this
is not necessary for uniformity nor efficiency.
Our main results give precise bounds for the number of covariates selected and the estimation
error, both for the penalized and unpenalized estimates. Section 6 results gives nonasymptotic
bounds, with exact constants. Such results are complex and so we give the following intuitive,
asymptotic result (The notation OPn is defined in Section 2.3).
Corollary 1. Suppose the biases from the best sd - and sy -term approximations to pt (x) and µt (x)
p
p
are order sd /n and sy /n, respectively. Then under the assumptions in Section 6, and δ > 0
described therein, with high probability we have:
1.

Pn


− pt (xi ))2 /n = OPn n−1 sd log(p ∨ n)3/2+δ and

2.

Pn


− µt (xi ))2 /n = OPn n−1 sy log(p ∨ n)3/2+δ .

i=1 (p̂t (xi )
i=1 (µ̂t (xi )

These two results for our proposed group lasso estimators can be directly used to verify the
high-level conditions in Theorem 1 above. Specifically, if sd sy log(p)3+2δ = o(n), conditions (i) and
(ii) of Theorem 1 are met (requiring s2 = o(n), up to log factors, as found in other results in the
literature). Further, it is clear how the doubly-robust estimator can help: if one function is more
smooth or more sparse, sd or sy will be lower, easing the restriction. Section 6.3 gives further
results: showing that the number of variables selected is the same order as the sparsity level, and
provides bounds on the logistic and linear coefficients directly. Both these results are important for
certain steps in treatment effect estimation that aren’t reflected in the simple statement of Theorem
1. These results appear to be entirely new for the multinomial logistic regression, for any version
of the lasso. From a practical point of view, these results provide formal justification for using
multinomial logistic regression, coupled with group lasso selection and post-selection refitting.

7

2.3

Notation

We collect here notation to be used for the rest of the paper. The data generating process (DGP)
is denoted by Pn and is defined by the joint law of the random variables (Y, D, X 0 )0 . For a given
n, {(yi , di , x0i )0 }ni=1 constitute draws from Pn . The DGP may vary with n, along with features such
as parameters, distributions, and so forth, as discussed in Section 4.2. This is generally suppressed
for clarity. We adopt the following conventions.
Treatments. Define the treatment sets NT = {0, 1, 2, . . . , T} and NT = {1, 2, . . . , T}. No order is
assumed in the treatments. For each unit i, di indicates treatment assignment, and define
P
dti = 1{di = t}. Let nt = ni=1 dti be the number of individuals with treatment t and define
n = mint∈NT nt and n = maxt∈NT nt . Further define T = T + 1.
Vectors. Define Np = {1, 2, . . . , p}. For a doubly-indexed collection of scalars {δt,j : t ∈ NT , j ∈
Np }, define δ·,j ∈ RT as the vector that collects over all t for fixed j; δt,· ∈ Rp collects over
j ∈ Np for fixed t; and δ·,· ∈ Rp×T the concatenation of all δt,· . For simplicity, we write δt
for δt,· . When considering the multinomial logistic model, t will vary only over NT but the
notation will be maintained. For a set S ⊂ Np , let δt,S ∈ Rcard(S) be the vector of {δt,j : j ∈ S}
for fixed t and similarly let δ·,S ∈ R|S|×T = {δt,j : t ∈ NT , j ∈ S}.
Norms. Single bars will be either absolute value or cardinality of a set, and will be clear from the
context. For a vector v, let kvk1 and kvk2 denote the `1 and `2 norms, respectively. For the
P
group lasso, define the mixed `2 /`1 norm as |||δ·,· |||2,1 = j∈Np kδ·,j k2 . It will always be the
case that the (“outer”) `1 norm is over the covariates and the (“inner”) `2 norm is over the
treatments (in our application). When discussing the multinomial logistic model, treatments
will be restricted to NT with no change in notation.
Data-Generating Processes. The set of all Pn considered is Pn . For sequences, {Pn } = {Pn :
n ≥ 1, Pn ∈ Pn }. Expectations and probabilities are taken against Pn , though notationally
suppressed. For asymptotic arguments dependence on n is explicit, so that OPn (·) and oPn (·)
have their usual meaning with the understanding that the measure Pn is used for each n.
For a set of scalars {mt }Tt=1 , let p̂t ({mt }NT ) = exp(mt )[1 +

exp(mt )]−1 denote the multiP
nomial logit function. Empirical expectation will be denoted En [wi ] = ni=1 wi /n and En,t [wi ] =
P
Pn
t
i∈It wi /nt =
i=1 di wi /nt .

3

P

t∈NT

Treatment Effects Model

In this section we formally define the treatment effects model and the parameters of interest. Recall
that D ∈ {0, 1, . . . , T} indicates treatment status, {Y (t)}t∈NT are the (latent) potential outcomes,
8

and Y (t) is only observed for units with D = t; that is, Y =

P

t∈NT

Y (t). The building blocks of

many general estimands are the averages
µt = E[Y (t)],

t ∈ NT ,

and

µt,t0 = E[Y (t)|D = t0 ],

t, t0 ∈ NT × NT .

(1)

In the binary case, the average treatment effect is µ1 − µ0 and the treatment on the treated is
µ1,1 − µ0,1 . A multivalued treatment allows for a large range of interesting estimands. To fix ideas,
we keep as running examples two leading cases. First, the so-called dose-response function: the
(T + 1)-vector µ = (µ0 , µ1 , . . . , µT )0 . Second, define τ as the T-vector with element t given by
µt,t − µ0,t . This gives the effect of each treatment relative to the baseline t = 0, only for those who
received that treatment. These are by no means the only interesting estimands constructed from
µt and µt,t0 ; many others are given by Lechner (2001), Heckman and Vytlacil (2007), and others.
The following two conditions are sufficient to identify µt and µt,t0 .
Assumption 1 (Identification). For all t ∈ NT and almost surely X, Pn obeys:
(a) (Mean independence) E[Y (t)|D, X = x] = E[Y (t)|X = x], and
(b) (Overlap) P[D = t|X = x]) ≥ pmin > 0 for all t ∈ NT .
This assumption is a form of “ignorability” coined by Rosenbaum and Rubin (1983). This
model allows arbitrary treatment effect heterogeneity in observables, but not unobservables. This
assumption is standard in the program evaluation literature, and its plausibility has been discussed
at length, so we omit a general discussion (see, e.g., Imbens (2004), Wooldridge (2010, Chapter 21),
and references therein). However, in the context of model selection, three remarks are warranted.
First, in place of 1(a), it is more common to instead assume full conditional independence:
Y ⊥
⊥ D|X. However, as observed by Heckman, Ichimura, and Todd (1997), the weaker mean
independence is sufficient. For our purposes, the “gap” between the two assumptions is important.
Suppose full independence holds only conditional on a set of variables strictly larger than the
variables entering the mean functions (e.g. the excess variables affect higher moments). In this
case, because mean independence is still sufficient, we need not aim to select the larger set. Full
independence is important for the efficiency discussed in Section 5.3.
Second, the covariates may, in general, include instruments for treatment status, but they are
not known as such. This is standard, but left implicit, in discussions of ignorability. If instruments
are present, and selected for estimation, efficiency suffers but unbiasedness is not harmed. Efficiency
bounds in this context typically (implicitly) assume there are no instruments in X. Assumption
1(b) rules out perfect predictors. Section 5.3 offers further discussion.
Finally, the main drawback of Assumption 1(a) is that it does not give identification of average
effects on transformations of Y (t). However, we are expressly interested in model selection on the
mean function of the level of Y (t), and hence Assumption 1(a) is more natural. To operationalize
9

model selection, structure must be placed on E[Y (t)|X = x], and hence functional form conditions
tied to mean independence are not limiting per se. If the parameter of interest is changed, say to
E[log(Y (t))], and a sparsity assumption is made for E[log(Y (t))|X = x], then our method applies.
Assumption 1 yields identification of µt and µt,t0 using either inverse weighting or regression,
and double robustness follows from combining the two strategies. Recall the notation pt (x) =
P[D = t|X = x] and µt (x) = E[Y |D = t, X = x]. Applying Assumption 1 we find that




1{D = t}Y
1{D = t}µt (X)
E ψt Y, D, µt (X), pt (X), µt = E
+ µt (X) −
− µt = 0
pt (X)
pt (X)

(2)

and


E ψt,t0 Y, D, µt (X), pt (X), pt0 (X), µt,t0


1{D = t0 }µt (X) pt0 (X) 1{D = t}(Y − µt (X))
=E
+
− µt,t0 = 0, (3)
pt0
pt0
pt (X)
where pt = P[D = t]. The moment condition (2) holds if either pt (x) or µt (x) is misspecified. For
µt,t0 , if µt (x) is misspecified, both pt (X) and pt0 (X) must be correctly specified, while if µt (x) is
correct, both propensity scores may be misspecified. It is important to note that the forms of ψt (·)
and ψt,t0 (·) are fixed, so the function itself does not depend on the sample size even if its arguments
do. Our estimator is a plug-in version of this moment condition.
Remark 1 (Simplifications for µt,t ). Identification of µt,t does not require Assumption 1. Y (t) is
fully observed for the sub-population of interest and so a simple average will deliver µt,t = E[1{D =
t}Y ]/pt . Note that (3) reduces to this when t = t0 . For τ this means we must only estimate the
function µt (xi ) for t = 0. Intuitively, we must use comparison group observations to proxy for
treated units, but not the other way around. Thus, for certain parameters of interest, Assumption
1 can be weakened to hold only for the comparison group. However, we cover generic estimands,
without necessarily specifying a comparison group, and so we maintain Assumption 1 for simplicity,
rather than keeping track of hosts of special cases.



Remark 2 (Efficient Influence Functions). The efficient influence functions in this model are
exactly ψt (·) and ψt,t0 (·), and so our estimators have the interpretation of being plug-in versions of
these, and indeed, will be asymptotically linear with this influence function (see Section 5.3).

4



Approximately Sparse Models

∗ be p-dimensional transformations of the
We now formalize approximate sparsity. Let XY∗ and XD

covariates X, with p > n allowed. These transformations are specific to the outcome and treatment
models, but may overlap. They do not vary with t, nor depend on the DGP. Some examples are
10

given below in Section 4.1. For the multinomial logistic model it is convenient to work with the
P
log-odds ratio. We take p0 (x) = 1 − t∈NT pt (x) and write

log

pt (x)
p0 (x)



= x∗D 0 γt∗ + BtD ,

t ∈ NT .

(4)

Similarly, write the outcome regressions as
µt (x) = x∗Y 0 βt∗ + BtY ,

t ∈ NT ,

(5)

The terms BtD = BtD (x) and BtY = BtY (x) are bias terms arising from the parametric specification. As discussed below, these encompass the usual nonparametric bias as well. Approximate
sparsity requires that only a small number of the X ∗ are needed to make the bias small. Define
S
S
S∗D = NT supp(γt∗ ) and S∗Y = NT supp(βt∗ ), so that these sets capture all variables important for
treatment and outcomes, respectively. We assume that there are some sd < n and sy < n, such
that for |S∗D | = sd and |S∗Y | = sy , the biases BtD and BtY are sufficiently small. This is made precise
by defining the bounds:
En [(p̂t ({x∗i 0 γt∗ }NT ) − pt (xi ))2 ]1/2 ≤ bds

and

En [BtY (xi )2 ]1/2 ∨ En,t [BtY (xi )2 ]1/2 ≤ bys .

(6)

Note that the former bias bound is placed directly on the propensity score because it is the ultimate
object of interest, rather than on the linearization of the log-odds.
While a great deal of overlap is expected, in practice it is likely that a few covariates will be
more or less important for different treatments, and so we do not require that the supports of
γt∗ , t ∈ NT or βt∗ , t ∈ NT are constant over t, nor that S∗D overlaps with S∗Y . Instead, it may be
better to think of Np \ S∗D and Np \ S∗Y as the “common nonsupports” of the treatment and outcome
∗ and X ∗ by X ∗ (and their
equations. When it is clear from the context we will abbreviate both XD
Y

realizations by x∗i ) and refer to them generically as “covariates”, and further write s for either sd
or sy . We assume En [(x∗i,j )2 ] = 1 without loss of generality (see Remark 4).

4.1

Parametric and Nonparametric Examples

To concretize the sparse model idea, we now discuss how several models commonly used in practice
fit into this framework. These include parametric and nonparametric models for pt (x) and µt (x),
and hybrids of these. A common theme to all examples will be comparison to the oracle model:
the model that knows the true support in advance. Our uniform inference results include all these
examples as special cases because, loosely speaking, we obtain uniformity over DGPs where pt (x)
and µt (x) have sparse representations. We aim for an accessible discussion of each model, and
defer technicalities to the literature (Raskutti, Wainwright, and Yu 2010, Rudelson and Zhou 2013,

11

Belloni, Chernozhukov, and Hansen 2014).
Example 1 (Oracle parametric model). Assume models (4) and (5) hold with BtD = BtY = 0
∗ = X ∗ = X. Let p = s = dim(X). All covariates are used in all modeling. If dimension
and XD
Y

is fixed this is the textbook parametric model, see for example Wooldridge (2010). Alternatively,
the dimension can be diverging, but more slowly than n. We are not aware of any work which
covers this case explicitly, though for the first stage, He and Shao (2000) cover linear and logistic
regression, and their results easily extend to multinomial logistic models.
The vast majority of treatment effect studies adopt this model (with dimension fixed), taking
the set of covariates as given. In our framework, this is equivalent to the researcher having prior
knowledge of which covariates are important and which are not. Such knowledge no doubt plays
an important role, but it cannot cover all situations or all variables. Furthermore, as more data
become available, the researcher does not increase the complexity of their model.



Example 2 (Exactly sparse parametric model). Retain the exact parametric structure of the prior
example, but let dim(X) = p be possibly larger than n, and assume that S∗Y and S∗D are unknown
sets of cardinality less than n. Model selection must be performed. Often, researchers (implicitly)
rely on the oracle property, that S∗Y and S∗D can be found with probability approaching one, and
conduct inference conditioning on this event. This approach cannot be made uniformly valid and
has poor finite sample properties, as shown by Leeb and Pötscher (2005, 2008a, 2008b, 2009).



Example 3 (Approximately sparse parametric model). Again suppose a purely parametric model,
∗ = X ∗ = X and dim(X) = p, possibly greater than n. Suppose that there exist
so that XD
Y
0 and β 0 such that log[p (x)/p (x)] = x∗ 0 γ 0 and µ (x) = x0 β 0 exactly, but instead
coefficients γ·,·
t
0
t
·,·
t
D t
0 | ∝ j −αγ and
of any coefficients being precisely zero, suppose they may be ordered such that |γt,j
0 | ∝ j −αβ , with α and α at least one. Then, there exist s and s that are o(n) such that
|βt,j
γ
γ
y
d
∗ = γ 0 for j ≤ s and
Equations (4) and (5), and other conditions needed, are satisfied for γt,j
d
t,j
∗ = β 0 for j ≤ s and the rest truncated to zero. That is S D and S Y collect the largest
βt,j
y
∗
∗
t,j
P
D
0 , and similarly for B Y .
coefficients and Bt = Np \S∗D xj γt,j

t

Example 4 (Semiparametric model). Assume pt (x) and µt (x) are unknown functions that can
be well-approximated by a linear combination of sd and sy basis functions, respectively (e.g. are
∗ and β ∗ are the coefficients of these approximations, while
sufficiently smooth). In (4) and (5), γ·,·
·,·
∗ = R (X) and X ∗ = R (X) are series terms
BtD and BtY are the usual nonparametric biases. XD
D
Y
Y

used in the approximation. Standard semiparametric analyses, such as Hirano, Imbens, and Ridder
(2003), Imbens, Newey, and Ridder (2007), or Cattaneo (2010), can be viewed in this context as
oracle models that know in advance which terms yield the best approximation, typically assumed to
be the first terms. Instead, we only require that some sd (or sy ) of a set of p series terms give good
approximations. This allows for greater flexibility in applications, where there is no knowledge of
which series terms to use, and the researcher may want to mix terms from different bases.
12



Example 5 (Mixed parametric and semiparametric model). Partition X = (X1 , X2 ). Suppose
that the true log-odds function satisfies log[pt (x)/p0 (x)] = x01 γt1 + ht (x2 ) + Bt1 (x), where Bt1 (x) is
a specification bias and ht (·) is a smooth unknown function. For a set of basis functions RD (x2 ),
there will exist coefficients γt2 such that ht (x2 ) = RD (x2 )0 γt2 + Bt2 (x2 ) and so

log

pt (x)
p0 (x)



= x∗D 0 γt∗ + BtD ,

x∗D = (x01 , RD (x2 )0 )0 ,

0

0

γt∗ = (γt1 , γt2 )0 ,

and

BtD = Bt1 + Bt2 .

We require that some collection of variables and series terms give a good, sparse approximation,
without placing explicit conditions on how many of either. Implicitly, one will restrict the other.
For example, if the dimension of the parametric part is large, then we require that ht (·) can be
more easily approximated. We treat µt (x) the same. This example is closest to actual practice,
where some variables (e.g. dummies) enter in a known way and should not be considered part of a
nonparametric object, while other covariates must be considered flexibly.



It is important to note that misspecification of the type guarded against by double robustness
can arise in any type of model. In parametric cases, this is most often functional form misspecification. While this type of misspecification does not occur in nonparametrics, others are possible,
such as shape restrictions or separability assumptions being incorrect, or omitting relevant variables. None of these errors disappear asymptotically, and all of them are guarded against by use
of the doubly-robust estimator.

4.2

Conceptual considerations in n-varying DGPs

Much of the DGP, including parameters and distributions, is allowed to depend on n. Perhaps the
most salient features that do not depend on n are the set of treatments and the functions ψt and
ψt,t0 . It is likely that our results can be extended to accommodate a growing number of treatments,
∗ , and β ∗ must depend
but that is beyond the scope of our study. In the models (4) and (5), X ∗ , γ·,·
·,·

on n by construction. Our results on estimation of these models are nonasymptotic. For treatment
effect inference, we use triangular array asymptotics to retain the dependence on n of the DGP.
The interpretation of the results does, and should, change depending on what is assumed about
the DGP. To illustrate, let us return to Examples 2 and 4.
First, consider the simple parametric models of Example 2. In this case, µt = E[E[Y (t)|X]] =
E[X 0 ]βt∗ , which depends on n by construction, as the dimension is diverging. It may seem unnatural
that the parameter to be estimated depends on n, as we typically think of “true” parameters being
features of a (large) fixed study population. However, with a diverging number of covariates, there
(n1 )

is no fixed DGP. Indeed, if we estimate µt = µt

based upon n1 observations, and then proceed to
(n1 +n2 )

gather n2 more observations, when we re-estimate our target is now µt
resolution is as follows. First, the parameter of interest is

13

(∞)
µt

(n1 )

6= µt

. One possible

= E[Y (t)], which is defined without

(∞)

reference to covariates. We can view each successive n-dependent µt as an approximation of µt

based upon p = pn covariates. Note well that in our thought experiment, pn1 6= pn1 +n2 , and so
additional variables should have been collected for all n1 + n2 samples.
Contrast this with the semiparametric model in Example 4. It is common to assume the population DGP is fixed over n. The treatment effects may be constructed in terms of the underlying
(∞)

variables, e.g. µt

= E[Y (t)] = E[E[Y (t)|X]], with X ∗ serving only the purpose of aiding in ap-

proximating the regression functions. Model selection is performed on series terms, not underlying
∗ and β ∗ . If µ = E[X ∗ 0 ]β ∗ + E[B Y ] does not depend on n,
variables, to estimate the coefficients γ·,·
t
·,·
t
t
Y

the bias term, by definition, exactly compensates for the n-dependence in E[XY∗ 0 ]βt∗ . We emphasize
that our inference results allow for general n-dependence in the DGP, and interpretation by the
econometrician must take careful account of any conceptual assumptions.

5

Main Results on Treatment Effect Estimation and Inference

In this section we present results on uniformly valid treatment effect inference. We first present the
estimators and conditions required for a generic first stage to yield uniform inference. Although
our focus is on model selection and sparsity, our results are more general, showcasing the benefits
of doubly robust estimation for any model in Section 4 where Assumption 3 below (which does not
refer to selection or sparsity) can be satisfied.

5.1

Estimation Procedure with a Generic Model Selector

The moment functions ψt (·) and ψt,t0 (·) of Equations (2) and (3) have fixed and known form, and
so for estimators p̂t (x) and µ̂t (x), we can define
n

µ̂t =

1X
n
i=1

and

n

µ̂t,t0

1X
=
n
i=1

(




dti (yi − µ̂t (xi ))
+ µ̂t (xi )
p̂t (xi )

0

dti µ̂t (xi ) p̂t0 (xi ) dti (yi − µ̂t (xi ))
+
p̂t0
p̂t0
p̂t (xi )

(7)

)
,

(8)

where p̂t = nt /n. By combining these estimators appropriately we can construct estimators µ̂ and
τ̂ for the dose-response function µ and the vector τ , respectively, and any other estimand. Notice
that when t = t0 µ̂t,t is an average over the appropriate subpopulation: µ̂t,t = En,t [yi ].
Although in this section we allow for generic estimates p̂t (x) and µ̂t (x), it is important to
distinguish between estimates based upon selected sets that have no “additional randomness” and
those that do. Model selection based estimation will naturally have two steps: first data-driven
selection and then refitting to ameliorate the shrinkage bias and allow the researcher to augment

14

the selected variables. Let S̃ D and S̃ Y be the selected sets and Ŝ D and Ŝ Y be the final sets of
variables used in the refitting. We will say that these contain no “additional randomness” if the
added variables (i.e. Ŝ \ S̃, for Y or D) are nonrandomly selected, such as from economic theory or
prior knowledge. On the other hand, the added variables may be selected from a random process
beyond that included in S̃. The leading example would be using logistic-selected variables in the
regressions or vice versa. Then the variables used in µ̂t (xi ) depend not only on the randomness
of S̃ Y , but also on that of S̃ D , and hence on {di }ni=1 . Additional conditions are required for the
estimators with additional randomness.
The choice of method is in part dependent on the assumptions of the underlying model. To
∗ = X∗ .
illustrate, first, return to Example 2, where we have a purely parametric model with X = XD
Y

The researcher may want to set Ŝ D ⊃ S̃ D ∪ S̃ Y , in order to have a better chance that S∗Y ⊂ Ŝ D .
The set Ŝ D now contains additional randomness due to S̃ Y . Conversely, consider Example 4.
It is natural to include “low-order” basis functions for each underlying covariate, say linear and
quadratic polynomials. Thus, the researcher may want to include these in Ŝ, whether or not selected
by the group lasso. However, there is no reason that the series terms useful for approximating the
functions µt (x) would be useful for pt (x), or vice versa, and no additional randomness is injected.
We now state the sufficient conditions used for treatment effect estimation and inference. For
exposition, we present these in three groups: those concerning the underlying DGP, requirements
of p̂t (x) and µ̂t (x) in the “no additional randomness” case, and finally the additional conditions
to allow for “additionally random” selected sets. Begin with conditions on the DGP. Let U ≡
Y (t) − µt (X) and impose the following conditions.
Assumption 2 (Data Generating Process). Pn obeys the following, with bounds uniform in n.
(a) {(yi , di , x0i )0 }ni=1 is an i.i.d. sample from (Y, D, X 0 )0 .
(b) The covariates X ∗ have bounded support, with maxj∈Np |Xj∗ | ≤ X < ∞. Transformations may
depend on n but not the underlying data generating process.
(c) E[|U |4 | X] ≤ U4 .
(d) minj∈Np ,

t∈NT

E[Xj∗ 2 U 2 ] ∧ E[Xj∗ 2 (1{D = t} − pt (X))2 ] is bounded away from zero.

(e) For some r > 0: E[|µt (xi )µt0 (xi )|1+r ] and E[|ui |4+r ] are bounded.
These conditions are mild and intuitive, and not unique to high-dimensional models or model
selection. Assumption 2(a) restricts attention to cross-sectional applications. The condition of
bounded covariates is unlikely to be a limitation in practice. Any X ∗ that are underlying variables
will naturally be bounded in applications. This condition is automatically satisfied for most common
choices of basis functions employed in nonparametric estimation. The rest are moment conditions

15

on the potential outcome models, including allowing the errors to be heteroskedastic and nonGaussian. The uniform bounds in n are needed for array asymptotics.
We now give precise conditions on the model selector sufficient for uniformly valid inference.
Assumption 3 (First Stage Restrictions). The estimators p̂t (x) and µ̂t (x) obey the following for
a sequence {Pn }, uniformly in t ∈ NT .


(a) En [(p̂t (xi ) − pt (xi ))2 ] = oPn (1) and En (µ̂t (xi ) − µt (xi ))2 = oPn (1),
(b) En [(µ̂t (xi ) − µt (xi ))2 ]1/2 En [(p̂t (xi ) − pt (xi ))2 ]1/2 = oPn (n−1/2 ).
(c) En [(µ̂t (xi ) − µt (xi ))(1 − dti /pt (xi ))] = oPn (n−1/2 ).
These two collectively play the same role as the commonly-used, high-level requirement in semiparametrics that each first-step component separately converge at n−1/4 at least.10 Indeed, Belloni,
Chernozhukov, and Hansen (2014) employ just such a condition for each component. However, by
making use of the doubly-robust property we have the weaker conditions shown.11 The first is a
mild consistency requirement. The second requires an explicit rate on the product of errors, and
hence if one function is relatively easy to estimate Assumption 3(b) can be satisfied even if the other
does not converge at n−1/4 . This formalizes the benefit of doubly-robust estimation in general. In
high-dimensional, sparse modeling specifically the rates for the first stage depend on the sample
size, the number of covariates considered, and the sparsity level. Thus, if one function requires
fewer covariates to estimate, i.e. smaller p or s, then greater complexity can be allowed for in the
other (capturing, in particular, their relative smoothness).
The so-called “additional-randomness” estimators are more specific to the (approximately)
sparse model context, and so we now codify the sparsity requirements of Section 4 and then give
the additional conditions required for these estimators.
Assumption 4 (Sparsity). For each n, Pn obeys (4), (5), and (6), with |S∗Y | = sd and |S∗D | = sy .
Assumption 5 (Regularity conditions for union estimators).
o(n1/3 )

and the estimators pt (x) and µ̂t (x) obey the following, uniformly t ∈ NT :


max |ui | En [(p̂t (xi ) − pt (xi ))2 ] = oPn (n−1/2 )
i∈It

For a sequence {Pn }, log(p) =

and

kγ̂t − γt∗ k1 ∨kβ̂t −βt∗ k1 = oPn (log(p∨n)−1/2 ).

These conditions are needed to apply bounds for self-normalized sums (de la Peña, Lai, and
Shao 2009). Belloni, Chen, Chernozhukov, and Hansen (2012) were the first to use these techniques
10

See Newey and McFadden (1994) and Chen (2007), and references therein.
Many studies in the semiparametric literature relax or do not rely on the n1/4 condition, allowing the nonparametric portion to converge at a slower rate, at any rate, or in some cases be inconsistent; examples include Powell,
Stock, and Stoker (1989), Newey (1990), Robins, Li, Tchetgen, and van der Vaart (2008), Cattaneo, Jansson, and
Newey (2014a), and Cattaneo, Crump, and Jansson (2013, 2014b), among others.
11

16

in high-dimensional, sparse models. The first condition is high-level, but can be verified with
conditions on the errors and a bound for estimation. For the former, Belloni, Chen, Chernozhukov,
and Hansen (2012) assume that maxi∈Nn |ui | = OPn (n1/q ) for some q > 2. A larger q eases the
restriction in Assumption 5 but at the expense of stronger conditions on the noise distribution. For
example, if ui are assumed Gaussian, q can be taken to be any (large) positive number.
Remark 3 (Linear Probability Models). Our results cover use of a linear probability model for
pt (x), instead of the multinomial logistic form. All we require is a sufficiently high-quality approximation of the unknown function, and hence if Assumptions 3, and 5 if appropriate,12 are met then
uniform inference is possible using a linear probability model. Our group lasso results (Theorems
7 and 8) can be used directly to verify these conditions. In the same vein, multinomial logistic
regression can be used to estimate µt (x) if the outcome Y is discretely valued.

5.2



Theoretical Results

We now come to our main results on inference on average treatment effects. Most of our discussion
will concern µt and µ; similar points apply to results for µt,t0 and τ . Our first result formalizes
consistency of our estimates under misspecification.
Theorem 2 (Double Robustness). Consider a sequence {Pn } of data-generating processes. Suppose
that for some p0t (x) and µ0t (x), En [(p̂t (xi ) − p0t (xi ))2 ] = oPn (1) and En [(µ̂t (xi ) − µ0t (xi ))2 ] = oPn (1).
Let Assumptions 1 and 2 hold for each n, with the regularity conditions also holding for p0t (x) and
µ0t (x). If p0t (x) = pt (x) or µ0t (x) = µt (x), then |µ̂t − µt | = oPn (1).
This theorem formalizes the double-robustness property of our estimators: the propensity score
or regression may be misspecified if the limiting objects are well-behaved. Compare to Assumption
3(a). The nearly identical result for µt,t0 is omitted to save space.
We now turn to our main inference results. First we demonstrate a Bahadur representation of
a generic µ̂t or µ̂t,t0 . These are shown to be equivalent to a sample average of the moment functions
ψt (·) and ψt,t0 (·), respectively, after proper centering and scaling, evaluated at the true pt (xi ) and
µt (xi ). Using these results, asymptotic normality can be obtained for general estimands. We state
explicit results for the leading examples µ and τ .
An asymptotic variance formula is needed to state the results. Define the conditional variance
of the potential outcomes s σt2 (x) = E[U 2 |D = t, X = x] and the T-square matrix Vµ with elements
Vµ [t, t0 ] = 1{t = t0 }E
12




σt2 (X)
+ E [(µt (X) − µt )(µt0 (X) − µt0 )] ≡ VµW (t) + VµB (t, t0 ).
pt (X)

Assumption 5 can be slightly weakened in this case due to the linear link function.

17

Straightforward plug-in estimators for these two components are given by13
V̂µW (t)


= En

dti (yi − µ̂t (xi ))2
p̂t (xi )2



V̂µB (t, t0 ) = En [(µ̂t (xi ) − µ̂t )(µ̂t0 (xi ) − µ̂t0 )] .

and

Our first result gives the asymptotic behavior of µ̂t and µ̂ for a sequence of DGPs.
Theorem 3 (Estimation of Average Treatment Effects). Consider a sequence {Pn } of data-generating
processes that obey Assumptions 1, 2, and 3 for each n. If µ̂t (xi ) and p̂t (xi ) do not have additional
randomness in the estimated supports, we have:
1.

√

n(µ̂t − µt ) =

−1/2 √

2. Vµ

√
t
i=1 ψt (yi , di , µt (xi ), pt (xi ), µt )/ n

Pn

+ oPn (1);

n(µ̂ − µ) →d N(0, IT ); and

3. V̂µW (t) − VµW (t) = oPn (1) and V̂µB (t, t0 ) − VµB (t, t0 ) = oPn (1).
If, in addition, Assumptions 4 and 5 hold, then the same is true when the supports contain additional
randomness.
Theorem 3 itself may appear standard, but what is nonstandard is that the model selection step
of the estimation has been explicitly accounted for. This immediately gives the following uniform
inference results.
Corollary 2 (Uniformly Valid Inference). Let Pn be the set of data-generating processes satisfying
the conditions of Theorem 3 for a given n and G : RT → R be a fixed, twice uniformly continuously
differentiable function with gradient ∇G such that lim inf n→∞ k∇G (µ)k2 is bounded away from zero.
Then for cα = Φ−1 (1 − α/2), we have:

sup PP



q
0
G(µ) ∈ G(µ̂) ± cα ∇G (µ̂) V̂µ ∇G (µ̂)/n
− (1 − α) → 0.

P ∈Pn

Corollary 2 shows that these procedures are uniformly valid over the class of DGPs we consider,
and hence will be reliable in applications. The crucial insight that leads to uniform inference is to
change the goal of model selection away from perfect covariate selection (the oracle property) and
to high-quality approximation of the underlying functions (here pt (·) and µt (·)). This fundamental
shift in focus allows us to avoid the uniformity problems demonstrated by Leeb and Pötscher.
Assumption 3 formalizes exactly the quality of approximation needed. Such an approximation can
be found for any element in Pn , and hence inference is uniformly valid over that class. This method
of proving uniformity follows Belloni, Chernozhukov, and Hansen (2014) and Romano (2004), and
is distinct from the approach of Andrews and Guggenberger (2009).
13

Estimators can also be based on sample averages of outer products of influence functions, which would include
the covariance term that vanishes in expectation.

18

Results for the treatment effects on the treated are similar. The variance formula for τ is
slightly more cumbersome. Define the T-square matrix Vτ with elements


i
pt (X) h 2
pt (X)pt0 (X) 2
2
Vτ [t, t ] = 1{t = t }E
+E
σt (X) + (µt (X) − µ0 (X) − µt,t + µ0,t )
σ (X)
pt pt0 p0 (X) 0
p2t
0

0



≡ VτW (t) + VτB (t, t0 ).
Straightforward plug-in estimators for these two components are given by
V̂τW (t)


= En



i
p̂t (xi )p̂t0 (xi ) 0
dti h
2
B
0
2
and V̂τ (t, t ) = En
d (yi − µ̂0 (xi )) .
(yi − µ̂0 (xi ) − µ̂t,t + µ̂0,t )
p̂t p̂t0 p̂0 (xi )2 i
p̂2t

Note that we needn’t estimate µt (x) and σt2 (x), due to the simplification in Remark 1. With this
notation, we have the following results. Proofs are so similar to those for Theorem 3 and Corollary
2 that we omit them.
Theorem 4 (Estimation of Treatment Effects on Treated Groups). Consider a sequence {Pn } of
data-generating processes that obey Assumptions 1, 2, and 3 for each n. Then under Pn , as n → ∞,
if µ̂t (xi ) and p̂t (xi ) do not have additional randomness in the estimated supports:
1.

√

n(µ̂t,t0 − µt,t0 ) =

−1/2 √

2. Vτ

√
t
i=1 ψt,t0 (yi , di , µt (xi ), pt (xi ), pt0 (xi ), µt,t0 )/ n

Pn

+ oPn (1);

n(τ̂ − τ ) →d N(0, IT ); and

3. V̂τW (t) − VτW (t) = oPn (1) and V̂τB (t, t0 ) − VτB (t, t0 ) = oPn (1).
If, in addition, Assumptions 4 and 5 hold, then the same is true when the supports contain additional
randomness.
Corollary 3 (Uniformly Valid Inference). Let Pn be the set of data-generating processes satisfying
the conditions of Theorem 4 for a given n and G : RT → R be a fixed, twice uniformly continuously
differentiable function with gradient ∇G such that lim inf n→∞ k∇G (τ )k2 is bounded away from zero.
Then for cα = Φ−1 (1 − α/2), we have:



q
sup PP G(τ ) ∈ G(τ̂ ) ± cα ∇G (τ̂ )0 V̂τ ∇G (τ̂ )/n
− (1 − α) → 0.

P ∈Pn

5.3

Efficiency Considerations

The prior theoretical results are aimed at delivering robust inference. In this section, we briefly
discuss the efficiency of our estimator according to two criteria: semiparametric efficiency and
oracle efficiency. To put each on sound conceptual footing we separate discussion and restrict to
an appropriate set of models.
19

For semiparametric efficiency, pt (x) and µt (x) are nonparametric objects, as in Example 4,
X are fixed-dimension variables and the DGP does not vary with n. If we “upgrade” the mean
independence of Assumption 1(a) to full, namely {Y (t)}NT ⊥
⊥ D|X, then Theorems 3 and Theorem
4 immediately yield asymptotic linearity and semiparametric efficiency, attaining Hahn’s (1998) or
Cattaneo’s (2010) bounds. This requires there be no (known) instruments for treatment status in
X, as implicitly assumed in those works, else the bound may change (Hahn 2004).
Turning to oracle efficiency, an alternative to our robust approach is to prove that the true
support can be found with probability approaching one (the oracle property), then conduct inference
conditioning on this event. This approach cannot be made uniformly valid, but may be of interest
in the exactly sparse models of Example 2 (there is no “true” support in approximately sparse
models), because discovering the true support is equivalent to finding the variables in the causal
mechanism (White and Lu 2011), if one exists. This may be interesting in its own right, or for
future applications by way of hypothesis generation. The post oracle selection estimator is made
efficient by using only the variables important for µt (xi ) = E[Y |D = t, xi ]. This amounts to entirely
removing the instrumental variables indexed by S∗D \ S∗Y , whose inclusion would, in general, reduce
efficiency, though not increase bias. Further, S∗Y \S∗D are excluded from propensity score estimation.
Perfect selection requires two strong conditions: (i) an orthogonality condition on the Gram
matrices that restricts the correlation between the variables in and out of the true support (Bach
2008), and (ii) a beta-min condition bounding the nonzero coefficients away from zero. Intuitively,
highly correlated variables cannot be distinguished, nor can coefficients sufficiently close to zero be
found with certainty. Both bounds may depend on n, and in particular the lower bound on the
coefficients may vanish at an appropriate rate. Under such conditions, it is straightforward to show
that S∗Y and S∗D can be found with probability approaching one.

6

Group Lasso Selection and Estimation

We now give details for group lasso model selection and estimation, and make the refitting precise.
Section 6.1 discusses penalty choices and implementation. Restricted and sparse eigenvalues, key
quantities in our bounds, are discussed in Section 6.2. Our main nonasymptotic results are stated
in Section 6.3. These results are of interest more generally in the literature on high-dimensional
sparse models Finally, Section 6.4 gives asymptotic rates and verifies the conditions of Section 5.
We first select covariates by applying the group lasso penalty to the multinomial logistic loss
(for the propensity scores) and to least squares loss (to estimate the outcome regression). The loss
functions are defined as
M(γ·,· ) =

X



En −dti log p̂t ({x∗i 0 γt }NT )

and

t∈NT

E(β·,· ) =

X
t∈NT

20

En,t [(yi − x∗i 0 βt )2 ].

Then, the group lasso estimates for the propensity score coefficients, denoted γ̃·,· , and the regression
coefficients, β̃·,· , respectively solve
n
o
γ̃·,· = arg min M(γ·,· ) + λD |||γ·,· |||2,1

and

γ·,· ∈RpT

n
o
β̃·,· = arg min E(β·,· ) + λY |||β·,· |||2,1 ,

(9)

β·,· ∈RpT

where λD and λY are the penalty parameters discussed below and |||γ·,· |||2,1 is the mixed `2 /`1 norm.
To ameliorate the downward bias induced by the penalty and to allow for researcher-added
variables, we refit unpenalized models.14 Let S̃ D = {j : kγ̃·,j k2 > 0} and S̃ Y = {j : kβ̃·,j k2 > 0} be
the selected covariates and Ŝ D and Ŝ Y those used in refitting.15 We require Ŝ ⊃ S̃ and |Ŝ| ≤ s for
D and Y (we will prove that |S̃| ≤ s in both cases). The refitting estimators solve
γ̂·,· =

arg min

{M(γ·,· )}

and

β̂·,· =

{E(β·,· )} .

arg min

(10)

β·,· , supp(βt )=Ŝ Y

γ·,· , supp(γt )=Ŝ D

Remark 4 (Weighted Penalties). The group lasso penalty can be weighted in two ways. First, one
P
may weight the `2 portion, as in λD j∈Np kXj γ·,j k2 , where Xj is the design matrix for covariate
j, across all the treatments. Other weight matrices are possible, but with this choice, the estimate
is invariant to within group (treatment) reparameterizations, and is thus scale invariant for each
covariate. We therefore assume En [(x∗i,j )2 ] = 1 without loss of generality.
Second, the `1 norm can be weighted to give a penalty of the form λD

P

j∈Np

wj kγ·,j k2 . Two

common choices for wj are the number of variables in group j or an adaptive penalty from a pilot
estimate. Our groups are equally sized, and although adaptive procedures may improve oracle
properties (Zou 2006, Wei and Huang 2010), our goal is not perfect selection.

6.1



Choice of Penalty

We must specify choices of λD and λY for programs (9). From a theoretical point of view, these
must be chosen so that the penalty dominates the noise, which is captured by the magnitude of
the score in the dual of the |||·|||2,1 norm, with high probability. To acheive this, we set
√
2X T
λD = √
n

log(p ∨ n)3/2+δD
√
1+
T

!1/2

√
4XU T
and λY = √
n

log(p ∨ n)3/2+δY
√
1+
T

!1/2
,

(11)

14
∗
∗
The bias is away from the pseudo-true coefficients of the sparse parametric representation, γ·,·
and β·,·
. There is
no relation to specification biases BtD and BtY .
15
When supp(γt∗ ) and supp(βt∗ ) do not vary much over t, the group lasso is known to have better properties than the
ordinary lasso in terms of selection and convergence. Obozinski, Wainwright, and Jordan (2011) give a sharp bound
on the overlap necessary to yield improvements, while Huang and Zhang (2010), Kolar, Lafferty, and Wasserman
(2011), and Lounici, Pontil, van de Geer, and Tsybakov (2011) also demonstrate advantages of the group lasso
approach. These works show, among other things, that the group lasso advantage increases with large T, and with
the group structure, may perform better with smaller samples. We defer to the works cited for a formal discussion.

21

for some δD > 0 and δY > 0. With these choices, λD > 2 maxj∈Np kEn [(pt (xi ) − dti )x∗i,j ]k2 and
λY > 4 maxj∈Np kEn,t [ui x∗i,j ]k2 with probability 1 − P for a small (and shrinking) P. In generic
terms, λ is of the form Λ(1 + rn ), where Λ is an upper bound on the true score and rn is a rate that
depends on n and p.16 The specific rate chosen serves to balances the rate of convergence against
the concentration effect: a smaller rn would increase the rate of convergence, but at the expensive
of lowering the concentration probability 1 − P. In the Appendix we show that (for appropriate δ
and n or n) the concentration probability is given by
p
4 log(2p)(1 + 64 log(12p)2 )
P=
,
log(p ∨ n)3/2+δ

(12)

There are two practical methods to make these choices for feasible for implementation. When
p̂t (x) and µ̂t (x) are used to estimate average treatment effects, the decreased sensitivity of the final
estimate to the first stage, thanks to the doubly-robust estimator, in turn results in less sensitivity
to the choice of penalty (through the sparsity).17 The first option is an iterative procedure to
estimate the unknown X and U in λY and λD , as employed by Belloni, Chen, Chernozhukov,
and Hansen (2012) (validity of this procedure may be established along the same lines as in that
study). We use maxi≤n maxj∈Np |x∗i,j | for X and estimate U by iteration: given an initial estimate
(0)

(k−1)

µ̂t (x), set Û(k) = En [(yi − µ̂t

(k)

(xi ))4 ]1/4 , where µ̂t (xi ), k > 0, is based on Eqn. (10). In

implementation we found 10 iterations more than sufficient, and based the initial estimate on ridge
regression (with penalty chosen by cross validation). A second option is to select λY and λD directly
by cross-validation. This has the appealing feature that the precise forms of Eqn. (11) need not
be characterized and estimated. If interest lies in the underlying functions pt (x) and µt (x), cross
validation is appropriate as it minimizes a relevant loss function. Formal results establishing the
validity of cross-validation are not available, but it performs well in practice.

6.2

Restricted Eigenvalues

The local behavior of optimizations (9) and (10) is captured by their respective Hessians, which
involve the second moment matrix of the covariates. The eigenvalues of such matrices will be
explicit in our bounds. We are interested in finite sample bounds, and so we will only discuss the
empirical Gram matrices (see Remark 5). Define
Q = En [x∗i x∗i 0 ]

and

16

Qt = En,t [x∗i x∗i 0 ].

(13)

The slight differences in the two are as follows. The full sample has information on the logistic coefficients, so n
appears instead of n. No error bound appears in λD because the errors are bounded by one. The multiple 4 for λY ,
instead of 2, can be traced to the quadratic loss. These forms are determined at heart by the maximal inequality of
Lounici, Pontil, van de Geer, and Tsybakov (2011).
17
To our knowledge, no formal results exist on “optimal” penalty parameter choices for inference in high-dimensional
problems nor are any procedures free of user-specified choices.

22

In high-dimensional data, both are singular, and so we use restricted eigenvalues and sparse eigenvalues (Bickel, Ritov, and Tsybakov 2009).
For the multinomial logistic regression, the minimal restricted eigenvalue is defined by
(P
κ2D

0
t∈NT δt Qδt
∗ k2
kδ·,SD
2

≤ min
δ

)
:δ∈R

pT

\ {0}, δ·,{S∗D }c

2,1

≤ 4 δ·,S∗D

2,1

.

(14)

.

(15)

For least squares estimation we instead use
0
t∈NT δt Qt δt
kδ·,SY∗ k22

(P
κ2Y

≤ min
δ

)
pT

:δ∈R

\ {0}, δ·,{S∗Y }c

2,1

≤ 3 δ·,S∗Y

2,1

Note that Q appears for κD , whereas the Qt are used in κY . The restricted set, or cone constraint,
requires the magnitude of δ·,· off the true support be small relative to the true support, measured in
∗ ) and (β̃ −β ∗ ) obey the respective constraints.
the group lasso norm.18 We will show that (γ̃·,· −γ·,·
·,·
·,·
∗ ) and (β̂ − β ∗ ) from (10) may not obey the cone
In contrast, the refitting errors (γ̂·,· − γ·,·
·,·
·,·

constraint, but are sparse by construction. This motivates the use of sparse eigenvalues. For a set
S ⊂ Np and a p × p matrix Q̃, define
φ{Q̃, S}2 =

min

δ∈Rp , supp(δ)=S

δ 0 Q̃δ
kδk22

and

φ{Q̃, S}2 =

max

δ∈Rp , supp(δ)=S

δ 0 Q̃δ
.
kδk22

(16)

Finally, it will be useful to define a bound on φ{Q̃, S} over all subsets of a certain size. To this
end, for any integer m, define φ(Q̃, m) = maxS⊂Np , |S|≤m φ{Q̃, S}.
We take these quantities to be primitive, and defer to the literature. For example, van de Geer
and Buhlmann (2009), Huang and Zhang (2010), Raskutti, Wainwright, and Yu (2010), Rudelson
and Zhou (2013), and Belloni, Chernozhukov, and Hansen (2014). In particular, Huang and Zhang
(2010) show that the group lasso may need fewer observations to satisfy conditions on φ{Q̃, S}.
Remark 5. Often, invertibility of Q and Qt relies on their convergence to nonsingular population
counterparts.19 Some of the papers cited use this approach and our results can be restated in this
way by conditioning on the event that Q and Qt are close to their counterparts in the appropriate
sense, and adjusting the probability with which the conclusions hold. We instead take bounds to
be infinite if the minimum eigenvalues are zero.



18

The multiplier of 4 in the constraint for κD is traceable to the nonlinear model.
This is standard in fixed-dimension models, and has been used for diverging-dimensions parametric models (He
and Shao 2000) and nonparametrics (Newey 1997, Huang 2003, Cattaneo and Farrell 2013, Belloni, Chernozhukov,
Chetverikov, and Kato 2015, Chen and Christensen 2015). The eigenvalue assumptions employed in those works are
conceptually the same as the the restricted eigenvalues used here, only restricted to the p < n case.
19

23

6.3

Finite Sample Theoretical Results

We now have the necessary notation and assumptions to state our theoretical results on group lasso
estimation, beginning with multinomial logistic regression, followed by a terse treatment of linear
models. Corollary 1 is a special case of the results in this section, see Section 6.4.
Our first result is a nonasymptotic bound on the group lasso estimates from (9).
Theorem 5 (Group Lasso Estimation of Multinomial Logistic Models). Suppose Assumptions 1(b),
2(a), 2(b), 2(c), and 4 hold. Define Ap = pmin /(0 ∨ (pmin − bds )) and

√ 
p

T
d
RM = Ap pmin TAK 6λD |S∗ |κ−1
+
8b
s T ,
D
o−1
n
√
√ 
p
. Then with
for AK > 2κ2D κ2D − (2/3)X T 30λD |S∗ | + 100 |S∗ |κD bds T + 80κ2D (bds )2 Tλ−1
D
probability 1 − P, we have
1. max En [(p̂t ({x∗i 0 γ̃t }NT ) − pt (xi ))2 ]1/2 ≤ RM + bds ,
t∈NT

2. max kγ̃t − γt∗ k1 ≤ RM
t∈NT

q

∗ | φ{Q, S̃ D ∪ S ∗ },
|S̃ D ∪ SD
D



3. and |S̃ D | ≤ 8sLn min{φ(Q, m) : m ∈ ND
}
,
Q
o
n

√ 2
φ(Q,
m)
and Ln = T (RM + bds ) (λD s) .
where ND
=
m
∈
{1,
2,
.
.
.
n}
:
m
>
8sL
n
Q
This theorem is new to the literature, to the best of our knowledge. Much of the detail involves
capturing the finite sample behavior of the Hessian and Gram matrices. We discuss the features of
this result in the following remarks.
• The Hessian of M(γ·,· ) is En [Hi ⊗ x∗i x∗i 0 ] for a T-square matrix Hi that depends on the
coefficients and x∗i through the estimated probabilities p̂t ({x∗i 0 γt }NT ). The error RM depends
on how well-controlled is this matrix. The factors pmin , Ap , and AK capture the behavior of
Hi and κ−1
D accounts for the rest. Under overlap, the true probabilities are bounded below
−T
by pmin , and hence pmin
captures the nonsingularity of the population version of Hi . To get

to this point requires two steps. First, the sparse parametric representations p̂t ({x∗i 0 γt∗ }NT )
must also be bounded away from zero, leading to the factor of Ap . This is essentially a bias
condition, which in the asymptotic case holds trivially: Ap may be chosen arbitrarily close to
one as bds → 0. Second, AK controls the neighborhood in which p̂t ({x∗i 0 γ̃t }NT ) is also bounded
away from zero. Intuitively (and asymptotically), the estimate will be in a small (shrinking)
neighborhood of the p̂t ({x∗i 0 γt∗ }NT ). In asymptotics AK may be chosen arbitrarily close to
2, which stems from the factor of 1/2 in a quadratic expansion of M(·). A lower bound on
AK is required in finite samples to ensure that p̂t ({x∗i 0 γ̃t }NT ) is positive, and hence the twoterm expansion is valid. This is analogous to Belloni and Chernozhukov’s (2011) “restricted
24

nonlinear impact coefficient” approach, also used by Belloni, Chernozhukov, Fernandez-Val,
and Hansen (2014) with a central difference that AK is captured in our bound directly.
• The maximal sparse eigenvalues are crucial to the bound on |S̃ D |. In many prior results, the
latter is bounded using the largest eigenvalue of Q itself, i.e. φ(Q, n). Adapting the technique
of Belloni and Chernozhukov (2013) to the present case, we are able to find a tighter bound,
which yields sparsity proportional to s under weaker conditions. This is crucial for refitting.
• For the linear model the constants in the group lasso bounds can offset the (logarithmic)
suboptimality in rate (Huang and Zhang 2010, Lounici, Pontil, van de Geer, and Tsybakov
2011), and this may be true here as well. This is application dependent however.
The error bounds for post-selection estimation are more complex and depend in part on the
good properties of the initial group lasso fit. The following theorem gives our results.
Theorem 6 (Post-Selection Multinomial Logistic Regression). Suppose the conditions of Theorem
∗ and φ = φ{Q, Ŝ ∪ S ∗ }. Then for
5 hold. To save notation, let SD = ŜD ∪ SD
D
D

(
AK > 2

)

φ2

√
√ p
φ2 − X T(λD |SD | + bds φ T |SD |)

(
∨

φ
√ p
φ − 2RM X T |SD |

)

 p
√ 
T
−1
d T and
0 = (A /p
TA
λ
/2
+
b
define RM
)
|S
|φ
p min
K
D
D
s
00
RM


= {RM } ∨

0
RM

i1/2 
h

T
2
0
.
+ RM RM + Ap pmin TAK RM

Then with probability 1 − P, max En [(p̂t ({x∗i 0 γ̂t }NT ) − pt (xi ))2 ]1/2
t∈NT
1/2 00
∗
D
max kγ̂t − γt k1 ≤ |S |/φ
RM .

≤

00
RM
+ bds ,

and

t∈NT

It is not readily discernible if these bounds improve upon the initial fit. This will depend on the
DGP, the selection success of the initial fit, and any added variables. In this result, further lower
bounds on AK are required to handle the sparse eigenvalues, compared to the restricted version in
Theorem 5. The role played by AK is the same in both cases, as with the other factors.
It is worth noting that, despite the complexity of multinomial logistic regression, the conditions
for Theorems 5 and 6 are simple and intuitive, and match those used for linear models.
We now give our results for group lasso estimation of the conditional outcome regressions. In
computing µt (xi ) for dti 6= 1 we are performing out of sample prediction, which slightly complicates
the bounds. Our first result is on the initial group lasso fit.

25

Theorem 7 (Group Lasso Estimation of Linear Models). Suppose Assumptions 1(b), 2(a), 2(b),
2(c), and 4 hold. To save notation, let SY = S̃ Y ∪ SY∗ . Define

RE =


√
3λY s
y
+ 2bs .
κY

Then with probability 1 − P, we have

1/2
1. max En [(x∗i 0 β̃t − µt (xi ))2 ]1/2 ≤ φ{Q, SY } φ{Qt , SY }
RE + bys ,
t∈NT

2. max β̃t − βt∗
t∈NT

1

1/2

1/2

φ{Q, SY } φ{Qt , SY }
RE ,
≤ |SY | φ{Q, SY }

o
n
P
3. and |S̃ Y | ≤ 32sLn minm∈NY t∈NT φ(Qt , m) ,
Q
n
o

P
√ 2
where NYQ = m ∈ {1, 2, . . . , n} : m > 32sLn t∈NT φ(Qt , m) and Ln = (RE + bys ) (λY s) .
This theorem generalizes Lounici, Pontil, van de Geer, and Tsybakov (2011) to the nonparametric, approximately sparse case, improves the sparsity bound, and gives out of sample prediction
(imputation) results. The analogous generalization for within sample prediction loss (e.g. multi-task
learning), En,t [(x∗i 0 β̃t − µt (xi ))2 ]1/2 , may be found in the Supplement.
For refitting, we are predicting for the entire sample and so we utilize the general results given
by Belloni, Chen, Chernozhukov, and Hansen (2012) for post-selection estimation of least squares.
The following result is a direct implication of their Lemma 7 and our Theorem 7.
Theorem 8 (Post-Selection Linear Regression). Suppose log(p) = o(n1/3 ) in addition to the conditions of Theorem 7. Then for constants A1 , A2 , A3 , and A4 not depending on n nor the DGP:
s
En [(x0i β̂t − µt (xi ))2 ]1/2 ≤ A1

s(T ∧ log(sT))
+ A2
nφ{Q, SY∗ }

s

|ŜY \ SY∗ | log(pT)
+ A3
nφ{Q, SYF P }

q

En [(x∗i 0 β̃t − µt (xi ))2 ]


1/2

.
and max kβ̂t − βt∗ k1 ≤ A4 |ŜY ∪ SY∗ |En [(x0i β̂t − µt (xi ))2 ] φ{Q, ŜY ∪ SY∗ }
t∈NT

As above, the performance of the refitting procedure depends in part on the success of the
initial group lasso fit. Indeed, the middle term is dropped if the true support union is found. The
constants Ak , k=1, 2, 3, 4 are not given explicitly but are known to be absolute bounds (de la
Peña, Lai, and Shao 2009) under Assumption 2. This result is less precise than Theorems 5 and 6,
but sufficient to verify Assumptions 3 and 5.

6.4

Asymptotic Analysis and Verification of High-Level Conditions

This section derives rates of convergence for the group lasso estimates and uses these results to
verify Assumptions 3 and 5 in Section 5. For simplicity, we only state results for the post-selection
26

estimators that we recommend in practice. In reducing the finite sample results of Theorems 6 and
8 to rates we retain the dependence on n, p, s, and the bias. Note that the number of treatments
is fixed, and the overlap assumption ensures that all nt ∝ n. Further, the various (restricted and
sparse) eigenvalues are commonly taken to be bounded (or bounded away from zero) in asymptotic
analyses. This accounts for the remaining factors in the bounds. For multinomial logistic regression,
we obtain the following result.
Corollary 4 (Asymptotics for Multinomial Logistic Regression). Suppose the conditions of Theorem 6 hold and further that (i) λD sd = o(1), (ii) κD is bounded away from zero, and (iii)
minS:|S|=O(s) φ{Q, S} is bounded away from zero and φ(Q, ·) is bounded, uniformly in ND
Q . Then
1. |S̃ D | = OPn (sd ),

2. En [(p̂t ({x∗i 0 γ̂t }NT ) − pt (xi ))2 ] = OPn n−1 sd log(p ∨ n)3/2+δD + (bds )2 ,vand
q
√ 
n−1 s2d log(p ∨ n)3/2+δD + bds sd .
3. kγ̂t − γt∗ k1 = OPn
Similarly, we have the following for the linear models.
Corollary 5 (Asymptotics for Linear Regression). Suppose the conditions of Theorem 8 hold and
√
further that (i) λY sy = o(1), (ii) κY is bounded away from zero, and (iii) uniformly in NT ,
minS:|S|=O(s) φ{Qt , S} ∧ φ{Q, S} is bounded away from zero and φ(Q, ·) ∨ φ(Qt , ·) is bounded uniformly in NYQ . Then
1. |S̃ Y | = OPn (sy ),

2. En [(µ̂t (xi ) − µt (xi ))2 ] = OPn n−1 sy log(p ∨ n)3/2+δY + (bys )2 , and
q
√ 
3. kβ̃t − βt∗ k1 = OPn
n−1 s2y log(p ∨ n)3/2+δY + bys s .
It is now straightforward to verify the requirements of Section 5. Assumption 3(b) requires

(n−1 sd log(p ∨ n)3/2+δD + (bds )2 )(n−1 sy log(p ∨ n)3/2+δY + (bys )2 ) = o n−1 .
Under the common assumption that bs = O(

p
s/n), we require sd sy log(p∨n)3+δD +δY = o(n). Both

this, and the display above, clearly show how the sparsity and smoothness of the two functions
interact due to the double robustness. Assumption 5 can be verified similarly.
These rates of convergence (i.e. part 2 of each corollary) are optimal up to factor log(p∨n)1/2+δ .
At heart, this loss appears to stem from the maximal inequality used to establish the concentration
probability of (12). In practice, this is unlikely to be a limitation. As mentioned above, the use of
group lasso can yield improvements in the constants if the data obey a grouped sparsity pattern,
as is expected for treatment effects data, and may even yield improvements in the detection of the
27

sparse signal, further offsetting the suboptimal log factor (see for example Lounici, Pontil, van de
Geer, and Tsybakov (2011) or Obozinski, Wainwright, and Jordan (2011)). Alternative methods
could, in principle, yield a rate improvement. Chief among these would be lasso-penalized linear
probability models (see also Remark 3) or separate logistic regressions. The group lasso approach
adopted here reflects common practice, and so it may be preferred. In any case, the log factors do
not impact the treatment effect inference.

7

Numerical and Empirical Evidence

7.1

Simulation Study

We conducted a Monte Carlo exercise to study how our estimator behaves as the propensity score
and regression functions change, and the model selection problem becomes more or less difficult.20
For simplicity we focus on the average effect of a binary treatment. We generated 1000 observations
(yi , di , x0i )0 from the models in Example 3, using both p = 1000 and p = 1500. The covariates include
an intercept, with the remainder drawn from N (0, Σ), with covariance Σ[j1 , j2 ] = 2−|j1 −j2 | , 2 ≤
j1 , j2 ≤ p. Errors are standard Normal. The crucial aspects of the DGP are the coefficient vectors
β00 , β10 , and γ 0 , which are defined to vary with the positive scalars ρβ , ργ , αβ , and αγ , as follows:
β00 = ρβ (−1, 1, −1, 2−αβ , −3−αβ , . . . , j −αβ , . . . , p−αβ )0 ,
γ 0 = ργ (1, −1, 1, −2−αγ , 3−αγ , . . . , j −αγ , . . . , −p−αγ )0 ,
with β10 = −β00 . The ρ multipliers affect the signal-to-noise ratio, but not the sparsity. For smaller
values distinguishing the large and small coefficients is more difficult for a given sample. The
exponents α control the sparsity, where a sparse representation is not possible for small values.
Figure 1 shows the empirical coverage rates of 95% confidence intervals for µ1 − µ0 for different
DGPs, for p = 1000 and 1500. Panels (a) and (c) show coverage as the multipliers ρβ and ργ
range over 0.01 (weak signal) to 1 (strong), with αβ = αγ = 2. Panels (b) and (d) vary the
sparsity exponents αβ and αγ over 1/8 (not sparse) to 4 (very sparse), with ρβ = ργ = 1. Of 1000
observations total, the (mean) size of the comparison group declines from roughly 500 to 300 as
ργ increases and 450 to 300 as αγ increases, over their given ranges. Coverage is accurate over all
signal strengths, and breaks down only when neither µt (xi ) nor pt (xi ) is sparse, which is exactly
when Assumption 3(b) (or condition (ii) of Theorem 1) cannot be satisfied. Note that coverage
accuracy is retained when only one function is sparse, showcasing the double-robustness property.
The penalty parameters λD and λY are chosen using the iterative procedure described in Section
6.1, with δD = 4.5 and δY = 5 throughout. Different DGPs exhibit different sensitivity to these
20

The supplemental appendix contains the additional results.

28

values. Results using penalties chosen via 10-fold cross-validation appear in Figure 2, which also
exhibits excellent coverage across all sparse designs.21

7.2

Empirical Application

To illustrate the role that model selection can play in a real-world application, we revisit the
National Supported Work (NSW) demonstration. The NSW has been analyzed numerous times
since LaLonde (1986). Our aim is a simple study of model selection, not a comprehensive or
conclusive evaluation of the NSW. We focus on the subsample used by Dehejia and Wahba (1999)
and the Panel Study of Income Dynamics (PSID) comparison sample, taking as given their data
definitions, sample selection, and trimming rules. Detailed discussion of these choices, and the NSW
program may be found in Dehejia and Wahba (1999, 2002) (hereafter DW99 and DW02) and Smith
and Todd (2005), and references therein. Briefly, the outcome of interest is earnings following a job
training program. The dataset includes a treatment indicator, post-treatment earnings (1978), two
years of pre-treatment earnings (197422 and 1975), as well as age, education, a marital status, and
indicators for Black and Hispanic. Thus, X consists of seven variables. We will keep the estimator
fixed: all estimates will be based on the doubly-robust estimator with standard errors from Section
5.2. We will compare the following specifications for X ∗ :
1. No Selection: X, (earn1974)2 , (earn1975)2 , (age)2 , and (educ)2 ;
2. Informally Selected: The above, plus 1{educ<HS}, 1{earn1974=0}, 1{earn1975=0}, and
(1{earn1974=0}×Hispanic). This specification was selected by DW02 using an informal balance test.
3. Group Lasso Selection: X, 1{educ<HS}, 1{earn1974=0}, 1{earn1975=0}, all possible firstorder interactions, and all polynomials up to order five of the continuous covariates (age, educ,
earn1974, earn1975).
For specifications 1 and 2, the same covariates are in the outcome and treatment models. All specifications include an intercept and we include education and pre-treatment income in the refitting
step following model selection. We follow DW99 and DW02 and trim comparisons with estimated
propensity score larger (smaller) than the maximum (minimum) in the treated sample.23
Table 1 presents results from these three specifications, and includes the experimental arm of
the NSW. The group lasso based estimate performs very well: the point estimate is accurate and
the interval is tight. Selecting from 171 possible covariates allows for a great deal of flexibility, but
21

The R routines appear unstable for nonsparse designs, thus the analogues to Panels (b) and (d) of Figure 1 are
omitted. See the supplement for limited versions. This will be explored for future software development.
22
This naming follows DW99, but the variable may be measured outside 1974, see discussion in the works cited.
23
A formal treatment of trimming is beyond the scope of the present study. The goal of our analysis is illustrative,
and hence we take DW99’s trimming as given. This issue is discussed by DW99, DW02, and Smith and Todd (2005).

29

the sparsity of the estimate keeps the variance well-controlled. The no-selection point estimate is
accurate, but fails to yield significance, while the specification of DW02 yields a significant, but
overly high estimate and wide confidence interval. The benefits of explicit model selection are clear.

8

Discussion

This paper proposed a method that achieves uniformly valid inference on mean effects of a multivalued treatment even after model selection among possibly more covariates than observations.
We demonstrated robustness to model selection errors, misspecification, and heterogeneous effects
in observables. To accomplish this, a doubly-robust estimator was employed and shown to have
excellent properties following model selection. We proved new results on group lasso estimation,
which we argue is natural for treatment effects data. Multinomial logistic regression was studied
in some detail. Numerical evidence shows that our method is quite promising for applications.
A key outstanding question in this work and in the high-dimensional, sparse modeling literature
more generally, is penalty parameter choice. Very little work has been done in this area, which
is a crucial gap in implementability of these techniques. We plan to develop a formal choice
for the penalty parameter that is appropriately optimal. Tuning parameter selection in semi- and
nonparametric analysis, and its impact on estimation and inference, is becoming better understood,
and parallel developments must take place in model selection contexts.

Appendix A

Proofs for Treatment Effect Inference

The proofs in this section are asymptotic. Order symbols hold for the sequence being considered,
as a shorthand for the more formal versions given in e.g. Assumption 3. C will denote a generic
positive constant, which may be a matrix. Define the set of indexes It = {i : di = t}. The online
supplement contains much greater detail. We make frequent use of the linearization
1
1 b−a
1 b − a (b − a)2
= +
= + 2 +
.
a
b
ab
b
b
ab2

(A.1)

Proof of Theorem 2. SEE SUPPLEMENTAL APPENDIX.
Proof of Theorem 3.1 without Additional Randomness. With ψt (·) defined in Eqn. (2), we have
√
√
n(µ̂t − µt ) = nEn [ψt (yi , dti , µt (xi ), pt (xi ), µt )] + R1 + R2 , where
n

1 X t
R1 = √
di (yi − µt (xi ))
n
i=1

and



1
1
−
p̂t (xi ) pt (xi )





n
1 X
dti
(µ̂t (xi ) − µt (xi )) 1 −
R2 = √
.
p̂t (xi )
n
i=1

30

The proof proceeds by showing that both R1 and R2 are oPn (1). Applying the first equality in Eqn.
(A.1), we rewrite R1 as


n
1 X t
pt (xi ) − p̂t (xi )
R1 = √
.
di ui
p̂t (xi )pt (xi )
n
i=1

Applying Assumptions 1(b) and 2(c) and the first-stage consistency condition of Assumption 3(a):


 2

dti σt2 (xi )
2
n
(pt (xi ) − p̂t (xi )) ≤ CEn [(pt (xi ) − p̂t (xi ))2 ] = oPn (1).
E R1 |{xi , di }i=1 = En
p̂t (xi )2 pt (xi )2
Next, again using Eqn. (A.1), we have R2 = R21 + R22 , where
n

1 X
R21 = √
(µ̂t (xi ) − µt (xi ))
n
i=1



pt (xi ) − dti
pt (xi )



and



n
dti
1 X
.
(µ̂t (xi ) − µt (xi ))(p̂t (xi ) − pt (xi ))
R22 = √
p̂t (xi )pt (xi )
n
i=1
√
For the first term, R21 = nEn [(µ̂t (xi ) − µt (xi ))(1 − dti /pt (xi ))] = oPn (1) by Assumption 3(c).
Next,


p
√
1
|R22 | ≤ n max
En [(µ̂t (xi ) − µt (xi ))2 ]En [(p̂t (xi ) − pt (xi ))2 ] = oPn (1).
i≤n p̂t (xi )pt (xi )
by Hölder’s inequality, Assumption 1(b) and the rate condition of Assumption 3(b).
Proof of Theorem 3.1 with Additional Randomness. We must reconsider the remainders R1 and
R2 . For the former, applying Eqn. (A.1), we find R1 = R11 + R12 , where
n

1 X dti ui
(pt (xi ) − p̂t (xi ))
R11 = √
pt (xi )2
n

n

1 X
dti ui
R12 = √
(p̂t (xi ) − pt (xi ))2 .
pt (xi )2 p̂t (xi )
n

and

i=1

i=1

For R11 , we first add and subtract the parametric representation to get R11 = R111 + R112 , where,
n

1 X dti ui
∗0 ∗
∗0
R111 = √
p̂
({x
γ
}
)
−
p̂
({x
γ̂
}
)
t
t
t
N
N
i
t
i
T
T
pt (xi )2
n

and

i=1

R112

n

1 X dti ui
∗0 ∗
√
=
p
(x
)
−
p̂
({x
γ
}
)
.
t
i
t
N
i
t
T
pt (xi )2
n
i=1

By a two-term mean-value expansion R111 = R111a + R111b , with
R111a

n

1 X dti ui X 
=√
p̂t ({x∗i 0 γt∗ }NT )(1 − p̂t ({x∗i 0 γt∗ }NT )) x∗i 0 (γ̂t − γt∗ )
2
pt (xi )
n
i=1

t∈NT

31

n

1 X dti ui 0
R111b = √
v H̄vi ,
pt (xi )2 i
2 n

and

i=1

{x∗i 0 (γ̂t

H({x∗i 0 γt∗

where vi =
−
and H =
+ mt x∗i 0 γ̂t }NT ) for appropriate scalars mt and the
0
T-square Hessian matrix H({x∗i γt }NT ) (defined in Appendix B).
P
For R111a , consider each term in the sum over NT one at a time; let R111a = t∈NT R111a (t).
Let t0 denote the original treatment under consideration. Define


Σt,j = E (x∗i,j )2 σt20 (xi )p̂t ({x∗i 0 γt∗ }NT )2 (1 − p̂t ({x∗i 0 γt∗ }NT ))2 /pt0 (xi )3 .
γt∗ )}NT

Then proceed as follows
(
!)
n
t0 u p̂ ({x∗ 0 γ ∗ } )(1 − p̂ ({x∗ 0 γ ∗ } ))
X
d
1 X
i
t
t
1/2
N
N
∗
i t
i t
T
T
√
R111a (t) =
x∗i,j i
Σt,j (γ̂t,j − γt,j
)
1/2
2Σ
n
0
p
(x
)
i
t
i=1
t,j
j∈ŜD
!


n
t0 u p̂ ({x∗ 0 γ ∗ } )(1 − p̂ ({x∗ 0 γ ∗ } ))
X
d
1
i t
t
1/2
i t NT
i t NT
kγ̂t − γt∗ k1
≤ max Σt,j
max √
x∗i,j i
1/2
j∈Np
j∈Np
2
n
pt0 (xi ) Σ
= O(1)OPn (log(p)) kγ̂t −

i=1
γt∗ k1

t,j

= oPn (1).

Convergence follows under Assumption 5. For the penultimate equality, it follows from Assumptions
1(b), 2(b), and 2(c) that maxj∈Np Σt,j = O(1). Finally, the center factor is shown to be OPn (log(p))
by applying the moderate deviation theory for self-normalized sums of de la Peña, Lai, and Shao
(2009, Theorem 7.4) and in particular Belloni, Chen, Chernozhukov, and Hansen (2012, Lemma
5). To apply this lemma, first note that the summand of the center factor has bounded third
moment and second moment bounded away from zero, from Assumptions 1(b), 2(b), 2(c), and the
requirements of Assumptions 3 and 5. Σt,j normalizes the second moment, and the lemma applies
under Assumptions 4 and the first restriction of Assumption 5.
For R111b , the results of Tanabe and Sagae (1992) coupled with Assumption 3 give vi0 H̄vi ≤
Ckvi k22 . Thus, using Assumption 1(b) to bound maxi≤n pt (xi )−2 < C, we find R111b may be
bounded as follows:
X√


|R111b | ≤ C
n(max |ui |)En |x∗i 0 (γ̂t − γt∗ )|2
i∈It

t∈NT

≤ CT max
t∈NT

√



n(max |ui |)En |p̂t ({x∗i 0 γ̂t }NT ) − p̂t ({x∗i 0 γt∗ }NT )|2 = oPn (1),
i∈It

by the union bound and Assumption 5, using the Assumptions 1(b) and 3(a) to apply Eqn. (B.15)
with the inequality reversed.
A variance bound may be applied to R112 as in the previous proof, and we have |R112 | =
OPn (bs ) = oPn (1) by Markov’s inequality.
Next, R12 is simply bounded by


h
i
√
1
2
|R12 | ≤ n(max |ui |) max
E
(p̂
(x
)
−
p
(x
))
n
t
i
t
i
i∈It
i∈It pt (xi )2 p̂t (xi )
h
i
√
≤ OPn (1) n(max |ui |)En (p̂t (xi ) − pt (xi ))2 = oPn (1),
i∈It

32

where the rate follows from Assumptions 1(b), 2, and 3, and this tends to zero by Assumption 5.
As in the prior proof, write R2 = R21 + R22 . The same bound is used for R22 . However, for
R21 , add and subtract the pseudotrue values to get R21 = R211 + R212 , where
n

R211

1 X ∗0
=√
(xi β̂t −x∗i βt∗ )
n



i=1

For the first term, define Σ̃t,j
R211

pt (xi ) − dti
pt (xi )

n


and

R212

1 X ∗ ∗
=√
(xi βt −µt (xi ))
n
i=1



pt (xi ) − dti
pt (xi )

h
i
= E (x∗i,j )2 (dti − pt (xi ))2 /pt (xi )2 and then proceed as follows:


n 
1 X pt (xi ) − dti X ∗
∗
xi,j (β̂t,j − βt,j
)
=√
pt (xi )
n
i=1
j∈ŜY
(
)
n
∗
X
1 X xi,j (pt (xi ) − dti )/pt (xi )
1/2
∗
√
Σ̃t,j (β̂t,j − βt,j
)
=
1/2
n
Σ̃t,j
i=1
j∈ŜY
!


n
1 X x∗i,j (pt (xi ) − dti )/pt (xi )
1/2
max √
≤ max Σ̃t,j
β̂t − βt∗
1/2
j∈Np
j∈Np
n
Σ̃
i=1
t,j
= O(1)OPn (log(p)) β̂t − βt∗

1

1

= oPn (1),

where the final line follows exactly as above. A variance bound may be applied to R212 as in the
previous proof, and we have |R212 | = OPn (bs ) = oPn (1) by Markov’s inequality.
Proof of Theorem 3.2. This follows from the prior result and Assumption 2(e).
Proof of Theorem 3.3. We begin with V̂W (t). Expanding the square and using Eqn. (A.1), rewrite
V̂µW (t) = En [dti u2i pt (xi )−2 ] + RW,1 + RW,2 + RW,3 where

dti u2i
RW,1 = En
(p̂t (xi ) − pt (xi )) (p̂t (xi ) + pt (xi )) ,
p̂t (xi )2 pt (xi )2
 t

 t

d (µt (xi ) − µ̂t (xi ))2
di ui (µt (xi ) − µ̂t (xi ))
= En i
,
and
R
=
2E
.
n
W,3
p̂t (xi )2
p̂t (xi )2


RW,2

Using Hölder’s inequality, Assumptions 1(b), 2(e), and 3(a), we have the following


p̂t (xi ) + pt (xi )
RW,1 ≤ max
En [dti |ui |4 ]1/2 En [dti (p̂t (xi ) − pt (xi ))2 ]1/2 = oPn (1),
i∈It p̂t (xi )2 pt (xi )2


1
RW,2 ≤ max
En [dti (µ̂t (xi ) − µt (xi ))2 ] = oPn (1),
i∈It p̂t (xi )2


1
and,
RW,3 ≤ 2 max
En [dti |ui |2 ]1/2 En [dti (µ̂t (xi ) − µt (xi ))2 ]1/2 = oPn (1),
i∈It p̂t (xi )2
where En [|ui |4 ] = OPn (1) from the inequality of von Bahr and Esseen (1965). From the same
inequality it follows that En [dti u2i pt (xi )−2 ] − VµW (t)| = oPn (1), under Assumptions 1(b) and 2(c).

33



Next consider the “between” variance estimator, V̂µB . For any tNT and t0 ∈ NT , define
RB,1 (t, t0 ) = En [(µ̂t (xi ) − µt (xi ))(µ̂t0 (xi ) − µt0 (xi ))] ,
RB,2 (t, t0 ) = µ̂t En [µ̂t0 (xi ) − µt0 (xi )] ,

and RB,3 (t, t0 ) = En [µt (xi )(µ̂t0 (xi ) − µt0 (xi ))] .

From Hölder’s inequality, Assumption 3(a), Theorem 3.2, the von Bahr and Esseen inequality, and
Assumptions 2(c) and 2(e) it follows that RB,k (t, t0 ) = oPn (1) for k ∈ N3 and all pairs (t, t0 ) ∈ N2t .
With this in mind, we decompose
V̂µB (t, t0 ) = En [µt (xi )µt0 (xi )] − µ̂t En [µt0 (xi )] − µ̂t0 En [µt (xi )] + µ̂t µ̂t0
+ RB,1 (t, t0 ) + RB,2 (t, t0 ) + RB,2 (t0 , t) + RB,3 (t, t0 ) + RB,3 (t0 , t).
Consistency of V̂µB (t, t0 ) now follows from the von Bahr and Esseen inequality and Theorem 3.2.
Proof of Corollary 2. Suppose the result did not hold. Then, there would exist a subsequence
Pm ∈ Pm , for each m, such that



q
0
lim PPm G(µ) ∈ G(µ̂) ± cα ∇G (µ̂)V̂ ∇G (µ̂)/n
− (1 − α) > 0.
m→∞

But this contradicts Theorem 3, under which (∇G (µ̂)V̂ ∇0G (µ̂)/n)−1/2 (G(µ̂) − G(µ)) is asymptotically standard normal under the sequence Pm .

Appendix B

Proofs for Group Lasso Selection and Estimation of
Multinomial Logistic Models

This section is nonasymptotic. We use generic notation X ∗ , δ, etc. The online supplement has
greater detail.

B.1

Lemmas

The following three lemmas are needed for the proofs of Theorems 5 and 6. Due to space considerations, only a short sketch of the proofs will be given, highlight the main ideas in each. Full details
are available in the online supplement.
Lemma B.1 (Score Bound). For λD and P defined in Eqn. (11) and Eqn. (12) we have


λD
t ∗
≤ P.
P max kEn [(pt (xi ) − di )xi,j ]k2 ≥
j∈Np
2
2 |x ] ≤
Proof. The residuals vt,i = pt (xi )−dti are conditionally mean-zero by definition and satisfy E[vt,i
i
h
i
∗
2
2
1. Using this, Assumption 2(a), and the definition of X, we find that E kEn [vt,i xi,j ]k2 ≤ X T/n,

uniformly in j ∈ Np . Define the mean-zero random variables ξt,j = (En [vt,i x∗i,j ])2 − n1 E[Vt2 Xj∗ 2 ] and

34

set rn = T −1/2 log(p ∨ n)3/2+δ . Then






2
X
X
λD
X Trn 
n
P max kEn [(pt (xi ) − dti )x∗i,j ]k2 ≥
≤ P max
ξt,j ≥
≤ E max
ξt,j  2
j∈Np
j∈Np
j∈Np
2
n
X Trn
t∈NT

t∈NT

where final line follows from Markov’s inequality. Next, applying Lemma 9.1 of Lounici, Pontil,
van de Geer, and Tsybakov (2011), Jensen’s inequality, and Assumption 2(c), we find that



 1/2
X 
X X4
4
 ≤ 4 log(2p)1/2 
+
E max En [vt,i x∗i,j ]  .
j∈Np
n2


E max
j∈Np

X
t∈NT

ξt,j



t∈NT

t∈NT

Again using Lemma 9.1 of Lounici, Pontil, van de Geer, and Tsybakov (2011), and Assumptions
2(a) and 2(b), we bound the expectation in the second term above as follows:


64 log(12p)2 X4
4
∗
E max En [vt,i xi,j ]
≤
.
j∈Np
n2
Collecting these results proves the Lemma.
Lemma B.2 (Estimate Sparsity). With probability at least 1 − P
|S̃ D | ≤

X


4
D
φ{Q,
S̃
}
En (p̂t ({x∗i 0 γ̃t }NT ) − pt (xi ))2 .
2
λD
t∈NT

Proof. From the Karush-Kuhn-Tucker conditions for (9), for all t ∈ NT , if γ̃·,j 6= 0 it must satisfy
En [x∗i,j (p̂t ({x∗i 0 γ̃t }NT ) − dti )] = λD

γ̃t,j
.
kγ̃·,j k2

Taking the `2 -norm over t ∈ NT for fixed j ∈ S̃ D , adding and subtracting the true propensity
score, using the triangle inequality, the score bound (B.1), collecting terms, squaring both sides,
and summing over j ∈ S̃ D (i.e. applying k · k22 over j ∈ S̃ D to both sides) yields
X
X
X X


λ2D ≤ 4
En [x∗i,j (p̂t ({x∗i 0 γ̃t }NT )−pt (xi ))]2 ≤ 4φ{Q, S̃ D }
En (p̂t ({x∗i 0 γ̃t }NT ) − pt (xi ))2 .
t∈NT

j∈S̃ D t∈NT

j∈S̃ D

The result now follows, as the left-hand side is equal to |S̃ D |λ2D .
∗ satisfies
Lemma B.3 (Bounds in `2 /`1 norm). With probability 1 − P the vector δ̃·,· = γ̃·,· − γ·,·
n
o
√
p
d
∗0
2 1/2 .
δ̃·,·
≤ an where an := max κ−1
≤ 5an and δ̃·,S∗
|S∗ |, 2λ−1
D
D bs T En [k{xi δ̃t }NT k2 ]
2,1

2,1

Proof. By the Cauchy-Schwarz inequality and Lemma B.1,
s
h
i2 s X
h
i
X
X
X
t ∗0
t
∗
2
En (pt (xi ) − di )xi,j
δ̃t,j
En (pt (xi ) − di )xi δ̃t ≤
t∈NT

j∈Np

t∈NT

t∈NT

35

≤ max
j∈Np

n


 oX
δ̃·,j
En (pt (xi ) − dti )x∗i,j 2
j∈Np

≤

2

λD
δ̃·,·
2

2,1

, (B.1)

with probability at least 1 − P. Applying the Cauchy-Schwarz inequality, the bias condition of
Assumption 4, and Cauchy-Schwarz again yields
X

h
i
h
i1/2
X

1/2
En (p̂t ({x∗i 0 γt∗ }NT ) − pt (xi ))x∗i 0 δ̃t ≤
En (p̂t ({x∗i 0 γt∗ }NT ) − pt (xi ))2
En (x∗i 0 δ̃t )2

t∈NT

t∈NT

h
i1/2
En (x∗i 0 δ̃t )2

X

≤ bds

t∈NT

√
≤ bds TEn [k{x∗i 0 δ̃t }NT k22 ]1/2 .

(B.2)

Combining Equations (B.1) and (B.2), we have, probability at least 1 − P,
h
i
h
i
X
X
En (p̂t ({x∗i 0 γt∗ }NT ) − dti )x∗i 0 δ̃t =
En (pt (xi ) − dti )x∗i 0 δ̃t
t∈NT

t∈NT

+

X

h
i
En (p̂t ({x∗i 0 γt∗ }NT ) − pt (xi ))x∗i 0 δ̃t

t∈NT

≤

λD
δ̃·,·
2

2,1

∗ + δ̃ ) + λ
∗
By the optimality of δ̃·,· , M(γ·,·
·,·
D γ·,· + δ̃·,·


λD

∗
γ·,·

−
2,1

∗
γ·,·


+ δ̃·,·

2,1

√
+ bds TEn [k{x∗i 0 δ̃t }NT k22 ]1/2 .

2,1

∗
∗
≥ M(γ·,·
+ δ̃·,· ) − M(γ·,·
)≥

∗ )+λ
∗
≤ M(γ·,·
D γ·,·

X

2,1

(B.3)

, and so

h
i
En (p̂t ({x∗i 0 γt∗ }NT ) − dti )x∗i 0 δ̃t ,

t∈NT

applying the convexity of M. Using the bound in Eqn. (B.3) and rearranging we find that


√
λD
∗
∗
0 ≤ λD
γ·,· 2,1 − γ·,· + δ̃·,·
+
δ̃·,·
+ bds TEn [k{x∗i 0 δ̃t }NT k22 ]1/2 .
2
2,1
2,1
Dividing through λD and decomposing the supports, we find that
√


bds T
1
∗
∗
+
+
γ·,· 2,1 − γ·,· + δ̃·,·
δ̃·,·
0≤
En [k{x∗i 0 δ̃t }NT k22 ]1/2
2
λD
2,1
2,1
1
δ̃·,S∗
=
2

1
+
δ̃·,S∗c
2
2,1

2,1

+

∗
γ·,S
∗

2,1

−

∗
γ·,S
∗

+ δ̃·,S∗

2,1

−

δ̃·,S∗c

√
bds T
+
En [k{x∗i 0 δ̃t }NT k22 ]1/2 ,
λD
2,1

∗
because γ·,S
c = 0. Collecting terms and applying the triangle inequality yields
∗

1
δ̃·,S∗c
2

2,1

1
δ̃·,S∗
≤
2

2,1

1
δ̃·,S∗
2

2,1

≤

+
+

√
bds T
En [k{x∗i 0 δ̃t }NT k22 ]1/2
−
+ δ̃·,S∗
+
2,1
λD
2,1
√


bds T
∗
En [k{x∗i 0 δ̃t }NT k22 ]1/2
− γ·,S∗ + δ̃·,S∗
+
λD
2,1

∗
γ·,S
∗
∗
γ·,S
∗

∗
γ·,S
∗

36

1
=
δ̃·,S∗
2

2,1

+

√
bds T
+
En [k{x∗i 0 δ̃t }NT k22 ]1/2 .
λD
2,1

δ̃·,S∗

Therefore with probability at least 1 − P
δ̃·,S∗c

2,1

≤ 3 δ̃·,S∗

√
2bds T
+
En [k{x∗i 0 δ̃t }NT k22 ]1/2 .
λD
2,1

(B.4)

Consider two cases based on the upper bound in (B.4). First, suppose that δ̃·,· obeys the cone
constraint of Eqn. (14) in the definition of κ2D . This implies
p
p
5 |S∗ |
En [k{x∗i 0 δ̃t }NT k22 ]1/2 ,
(B.5)
≤ 5 δ̃·,S∗
≤ 5 |S∗ | δ̃·,S∗ ≤
δ̃·,·
κD
2
2,1
2,1
by
inequality, the restricted eigenvalue definition of Eqn. (14), and noting that
P the Cauchy-Schwarz
0 Qδ̃ = E [k{x∗ 0 δ̃ } k2 ]. Collecting across the second and third inequalities yields
δ̃
t
n
i t NT 2
t∈NT t
δ̃·,S∗

2,1

p
|S∗ |
≤
En [k{x∗i 0 δ̃t }NT k22 ]1/2 .
κD

On the other hand, if the cone constraint fails, then

δ̃·,S∗

2,1

<

(B.6)
1
4

δ̃·,S∗c

2,1

. Using this for

the first and third inequalities, and Eqn. (B.4) for the second, we have
√
5 bds T
5
15
δ̃·,S∗c
δ̃·,S∗
+
En [k{x∗i 0 δ̃t }NT k22 ]1/2
δ̃·,·
≤
≤
4
4
2 λD
2,1
2,1
2,1
√
15
5 bds T
≤
δ̃·,S∗c
En [k{x∗i 0 δ̃t }NT k22 ]1/2
+
16
2 λD
2,1
Combining the right hand side of the first line with third lines yields
δ̃·,S∗c
√

8bds TEn [k{x∗i 0 δ̃t }NT k22 ]1/2 λD . Plugging this back into the last line we obtain the bound
δ̃·,·

2,1

√
bds T
≤ 10
En [k{x∗i 0 δ̃t }NT k22 ]1/2 ,
λD

while instead, plugging it into the failure of the cone constraint yields
√
bds T
δ̃·,S∗
≤2
En [k{x∗i 0 δ̃t }NT k22 ]1/2 .
λD
2,1

2,1

≤

(B.7)

(B.8)

Combining Equations (B.5) and (B.7) gives the first claim of the lemma and Equations (B.6) and
(B.8) give the second.

37

B.2

Proof of Theorem 5

∗ . By the optimality of δ̃ , we have
Define δ̃·,· = γ̃·,· − γ·,·
·,·
∗
∗
M(γ·,·
+ δ̃·,· ) + λD γ·,·
+ δ̃·,·

∗
∗
≤ M(γ·,·
) + λD γ·,·

2,1

2,1

.

Rearranging and subtracting the score, we have
X


∗
∗
M(γ·,·
+ δ̃·,· ) − M(γ·,·
)−
En (p̂t ({x∗i 0 γt∗ }NT ) − dti )x∗i 0 δ̃t
t∈NT


≤ λD

∗
γ·,·

−
2,1

∗
γ·,·


+ δ̃·,·

−

2,1

X



En (p̂t ({x∗i 0 γt∗ }NT ) − dti )x∗i 0 δ̃t .

(B.9)

t∈NT

The proof proceeds by deriving a further upper bound to the right and a quadratic lower bound of
the left. The combination of these will yield a bound on En [(x∗i 0 δ̃t )2 ]1/2 .
Begin with the right side of Eqn. (B.9). For the penalized difference of coefficients we have
∗
∗
∗
c
γ·,S∗c
− γ·,S
= δ̃·,S∗c
, because γ·,S
c + δ̃·,S∗
c = 0. Therefore,
∗
∗
2,1

2,1

∗
γ·,·

2,1

−

∗
γ·,·
+ δ̃·,·

2,1

2,1

=

∗
γ·,S
∗

2,1

≤

∗
γ·,S
∗

2,1

≤
≤

∗
γ·,S
∗

−

∗
γ·,S
+ δ̃·,S∗
∗

−

∗
γ·,S
+ δ̃·,S∗
∗

2,1

−


∗
∗
γ·,S
−
γ·,S
∗
∗

2,1

−

δ̃·,S∗c

2,1

∗
γ·,S
+ δ̃·,S∗
∗
2,1

+ δ̃·,S∗
= δ̃·,S∗
2,1

2,1

2,1

,

where the first inequality reflects dropping the nonpositive final term (the norm is nonnegative)
and the third inequality follows from the triangle inequality. Using this result for the first term
and the bound (B.3) for the second, the right side of Eqn. (B.9) is bounded by
λD δ̃·,S∗
≤
≤

√
λD
+
δ̃·,·
+ bds TEn [k{x∗i 0 δ̃t }NT k22 ]1/2
2
2,1
2,1
( p
!
(p
√ )
√ )
√
|S∗ | 2bds T
λD 5 |S∗ | 10bds T
d
λD
∨
+
∨
+ bs T En [k{x∗i 0 δ̃t }NT k22 ]1/2
κD
λD
2
κD
λD
!
p
√
λD |S∗ |
d
+ 8bs T En [k{x∗i 0 δ̃t }NT k22 ]1/2 ,
(B.10)
6
κD

where the second inequality applies Lemma B.3 and the third bounds the maximum by the sum.
Now turn to the left side of Eqn. (B.9). Our goal is to show that this is bounded below by
a quadratic function. We apply the bounds for Bach’s (2010) modified self-concordant functions.
To show that M(·) belongs to this class, we must
the third 
derivative in terms of the
 bound
P
0
0
0
∗
∗
∗
Hessian. Recall that p̂t ({xi γt }NT ) = exp{xi γt }/ 1 + NT exp{xi γt } and the T-square matrix

38

H({x∗i 0 γt }NT ) has (t, t0 ) ∈ N2T entry given by
(
p̂t ({x∗i 0 γt }NT )(1 − p̂t ({x∗i 0 γt }NT ))
H({x∗i 0 γt }NT )[t,t0 ] =
−p̂t ({x∗i 0 γt }NT )p̂t0 ({x∗i 0 γt }NT )

if t = t0
if t =
6 t0

First, note that M(γ·,· ) can be written as
 
 X

X
∗0
t ∗0
exp{xi γt } −
di (xi γt ) .
M(γ·,· ) = En log 1 +
t∈NT



Define F : RT → R as F (w) = log 1 +

t∈NT



h
i
P
tw
exp(w
)
,
so
that
M(γ
)
=
E
F
(w
)
−
d
t
·,·
n
i
t∈NT
t∈NT i i,t ,

P

where wi,t = x∗i 0 γt and wi = {wi,t }NT . Then for any w ∈ RT , v ∈ RT , and scalar α, define g(α) = F (w + αv) : R → R. We verify the conditions of Bach (2010, Lemma 1) for this
g(α) and F (w). This involves finding the third derivative of g(α), and bounding it in terms of
the second (i.e. the Hessian). To this end, note that the multinomial function has the property that ∂ p̂t ({x∗i 0 γt }NT )/∂γt = p̂t ({x∗i 0 γt }NT )(1 − p̂t ({x∗i 0 γt }NT ))x∗i and ∂ p̂t ({x∗i 0 γt }NT )/∂γt0 ,· =
−p̂t ({x∗i 0 γt }NT )p̂t0 ({x∗i 0 γt }NT )x∗i . From these, we find
X
g 0 (α) = v 0 F 0 (w + αv) =
vt p̂t (w + αv)
and
g 00 (α) = v 0 F 00 (w + αv)v = v 0 H(w + αv)v.
t∈NT

To bound g 000 (α), we again use the derivatives of p̂t ({x∗i 0 γt }NT ) to find the derivatives of elements
H(w). Routine calculations give, for any r 6= s 6= t:
∂H(w)t,t /∂wt = p̂t (w)(1 − p̂t (w))(1 − 2p̂t (w)) = H(w)t,t (1 − 2p̂t (w))
∂H(w)t,t /∂wr = −p̂t (w)p̂r (w)(1 − p̂t (w)) + p̂t (w)2 p̂r (w) = H(w)t,t (p̂t (w)p̂r (w)(1 − p̂t (w))−1 − p̂r (w))
∂H(w)t,s /∂wt = −p̂t (w)p̂s (w)(1 − 2p̂t (w)) = H(w)t,s (1 − 2p̂t (w))
∂H(w)t,s /∂wr = −p̂t (w)p̂s (w)(−2p̂r (w)) = H(w)t,s (−2p̂r (w)).
Each derivative returns the same Hessian element multiplied by term bounded by 2 in absolute
value. Let ar represent this factor. Then we bound
g 000 (α) =

X
r∈NT

≤

X

vr

∂v 0 H(w̃)v
∂wr

=
w̃=w+αv

X

vr v 0 H(w + αv)var

r∈NT

v 0 H(w + αv)v|vr ||ar | ≤ 2v 0 H(w + αv)v

r∈NT

X

√
|vr | = 2kvk1 g 00 (α) ≤ 2 Tkvk2 g 00 (α).

r∈NT

Applying Bach’s (2010) Lemma 1 to each observation, as in Belloni, Chernozhukov, and Wei (2013),
with wi = {x∗i 0 γt∗ }NT and vi = {x∗i 0 δ̃t }NT we get the lower bound
X


∗
∗
M (γ·,·
+ δ̃·,· ) − M(γ·,·
)−
En (p̂t ({x∗i 0 γt∗ }NT ) − dti )x∗i 0 δ̃t
t∈NT

39


vi0 H({x∗i 0 γt }NT )vi  −2kvi k2
≥ En
e
+ 2kvi k2 − 1
4Tkvi k22
 0


vi H({x∗i 0 γt }NT )vi
4
2
3
≥ En
2kvi k2 − kvi k2 ,
3
4Tkvi k22


(B.11)

where the second inequality follows from Belloni, Chernozhukov, and Wei (2013, Lemma 9).
Tanabe and Sagae (1992, Theorem 1) give H({x∗i 0 γt∗ }NT ) ≥ φmin {H({x∗i 0 γt∗ }NT )}IT , in the positive definite sense, where φmin (A) denotes the smallest eigenvalue of A and IT is the T × T identity
matrix. Then
Y
 T
φmin {H({x∗i γt∗ }NT )} ≥ det{H({x∗i 0 γt }NT )} =
p̂t ({x∗i 0 γt∗ }NT ) ≥ pmin Ap ,
t∈NT

where p0 ({x∗i 0 γt∗ }NT ) = 1 −

P

t∈NT

p̂t ({x∗i 0 γt∗ }NT ) and the first inequality is also due to Tanabe and

Sagae (1992). These results imply that vi0 H({x∗i 0 γt }NT )vi
and therefore
 0


 T
vi H({x∗i 0 γt }NT )vi
4
2
3
En
2kvi k2 − kvi k2
≥ pmin Ap
2
3
4Tkvi k2
 T
= pmin Ap

≥ (pmin /Ap )T vi0 IT vi = (pmin /Ap )T kvi k22


1
4
3
2
En 2kvi k2 − kvi k2
4T
3


2
1 En [kvi k2 ]
2 En [kvi k32 ]
1−
. (B.12)
T
2
3 En [kvi k22 ]

Recall that vi = {x∗i 0 δ̃t }NT . To prove a quadratic lower bound, consider two cases, depending
on whether
!
2 En [k{x∗i 0 δ̃t }NT k32 ]
1
1−
2
3 En [k{x∗i 0 δ̃t }NT k22 ]
is above or below 1/AK . In the first case, combining Equations (B.11) and (B.12) gives
∗
∗
M(γ·,·
+ δ̃·,· ) − M(γ·,·
)−

X
t∈NT



 T 1 En [k{x∗i 0 δ̃t }NT k22 ]
.
En (p̂t ({x∗i 0 γt∗ }NT ) − dti )x∗i 0 δ̃t ≥ pmin Ap
T
AK

(B.13)
Now consider the second case, where this bound does not hold. By Assumption 2(b), the CauchySchwarz inequality, and the conclusion of Lemma B.3
X X
√
k{x∗i 0 δ̃t }NT k1 =
x∗i,j δ̃t,j ≤ X δ̃·,· ≤ TX δ̃·,·
t∈NT j∈Np

2,1

1

≤

√

( p
√ )
5 |S∗ | 10bds T
TX
∨
En [k{x∗i 0 δ̃t }NT k22 ]1/2 .
κD
λD

Hence, by subadditivity (to bound the `2 norm by the `1 norm),
En [k{x∗i 0 δ̃t }NT k32 ]

≤

En [k{x∗i 0 δ̃t }NT k22 k{x∗i 0 δ̃t }NT k1 ]

≤

40

En [k{x∗i 0 δ̃t }NT k22 ]3/2

√

( p
√ )
5 |S∗ | 10bds T
TX
∨
.
κD
λD

Thus
1
1
>
AK
2

2 En [k{x∗i 0 δ̃t }NT k32 ]
1−
3 En [k{x∗i 0 δ̃t }NT k22 ]

!

1
≥
2

!
√
√ 
p
2 X T 
d
∗0
2 1/2
1−
,
5λD |S∗ | + 10κD bs T En [k{xi δ̃t }NT k2 ]
3 κD λ D

which is equivalent to


√ −1
p
3 κD λ D 
2
√
:= rn .
> 1−
5λD |S∗ | + 10κD bds T
AK 2 X T


P
∗ + δ ) − M(γ ) −
∗0 ∗
t ∗0
Because M(γ·,·
·,·
·,·
t∈NT En (p̂t ({xi γt }NT ) − di )xi δt is convex in δ·,· , and hence
any line segment lies above the function, we know that En [k{x∗i 0 δ̃t }NT k22 ]1/2 > rn , so we have
En [k{x∗i 0 δ̃t }NT k22 ]1/2

∗
M(γ·,·
+ δ̃·,· ) − M(γ·,· ) −

X
t∈NT



En [k{x∗i 0 δ̃t }NT k22 ]1/2
En (p̂t ({x∗i 0 γt∗ }NT ) − dti )x∗i 0 δ̃t ≥ rn2 ≥ rn2
rn
= rn En [k{x∗i 0 δ̃t }NT k22 ]1/2 .

Combining this result with Equations (B.9) and (B.10), we have


√ −1
p
2
3 κD λD 
√
5λD |S∗ | + 10κD bds T
1−
En [k{x∗i 0 δt }NT k22 ]1/2
AK 2 X T
!
p
√
λD |S∗ |
d
+ 8bs T En [k{x∗i 0 δ̃t }NT k22 ]1/2 ,
≤ 6
κD
which is impossible under the restriction on AK . Therefore, Eqn. (B.13) must hold.24 Combining
this with Equations (B.9) and (B.10), we find that
!
p
√
 T 1 En [k{x∗i 0 δ̃t }NT k22 ]
λD |S∗ |
≤ 6
+ 8bds T En [k{x∗i 0 δ̃t }NT k22 ]1/2 .
pmin Ap
T
AK
κD
Thus, dividing through and applying the union bound we find that
max En [(x∗i 0 δ̃t )2 ]1/2
t∈NT

≤

En [k{x∗i 0 δ̃t }NT k22 ]1/2


T
≤ Ap pmin TAK

6

λD

!
p
√
|S∗ |
d
+ 8bs T .
κD

(B.14)

To bound the propensity score error, we apply the mean value theorem and the form of
∂ p̂t ({x∗i 0 γt }NT )/∂γt . We must linearize with respect to t only (recall that p̂t ({x∗i 0 γ̃t }NT ) depends on
all of γ̃·,· ). To this end, define Mt as the T-vector with entry t given by x∗i 0 γt∗ + m̃t x∗i 0 γ̃t for a scalar
m̃t ∈ [0, 1] and entries t0 ∈ NT \ {t} equal to x∗i 0 γt0 . Then we have
p̂t ({x∗i 0 γ̃t }NT ) − p̂t ({x∗i 0 γt∗ }NT ) = p̂t (Mt )[1 − p̂t (Mt )]x∗i 0 δ̃t ≤ x∗i 0 δ̃t .
24

(B.15)

This analysis is conceptually similar to using Belloni and Chernozhukov’s (2011) restricted nonlinearity impact
coefficient, but our characterization is different.

41

Using this result coupled with the triangle inequality, the bias condition, and Eqn. (B.14), we find
En [(p̂t ({x∗i 0 γ̃t }NT ) − pt (xi ))2 ]1/2 ≤ En [(p̂t ({x∗i 0 γ̃t }NT ) − p̂t ({x∗i 0 γt∗ }NT ))2 ]1/2 + En [(p̂t ({x∗i 0 γt∗ }NT ) − pt (xi ))2 ]1/2
h
i1/2
+ bds
≤ En (x∗i 0 δ̃t )2
!
p
√

T
λD |S∗ |
≤ Ap pmin TAK 6
+ 8bds T + bds .
κD
The `1 bound follows from Eqn. (B.14), the Cauchy-Schwarz inequality, and Eqn. (16):
kγ̃t −

γt∗ k1

q
∗ | kγ̃ − γ ∗ k
≤ |S̃ D ∪ SD
t
t 2,p ≤

∗|
|S̃ D ∪ SD
∗}
φ{Q, S̃ D ∪ SD

!1/2
En [(x∗i 0 (γ̃t − γt∗ ))2 ]1/2 .

Finally, we bound the size of the selected set of coefficients. First, note that optimality of
γ̃·,· ensures that |S̃ D | ≤ n. Then, restating the conclusion Lemma B.2 using the notation of the
Theorem and the rate result (B.14), then bounding φ by φ we find that
∗
|S̃ D | ≤ |SD
|4Ln φ{Q, |S̃ D |}.

The argument now parallels that used by Belloni and Chernozhukov (2013), relying on their result
on the sublinearity of sparse eigenvalues. Let dme be the ceiling function and note that dme ≤ 2m.
D
For any m ∈ ND
Q , suppose that |S̃ | > m. Then,
∗
|S̃ D | ≤ |SD
|4Ln φ{Q, m(|S̃ D |/m)}
l
m
∗
≤ |S̃ D |/m |SD
|4Ln φ{Q, m}
∗
≤ (|S̃ D |/m)|SD
|8Ln φ{Q, m}.
∗ |8L φ{Q, m} whence m 6∈ ND . Minimizing over ND gives the result.
Rearranging gives m ≤ |SD
n
Q
Q

B.3

Proof of Theorem 6

∗ . Many of the arguments parallel those for Theorem 5. The key differences
Define δ̂·,· = γ̂·,· − γ·,·


P
∗ + δ̂ ) − M(γ ∗ ) −
∗0 ∗
t ∗0
are that a quadratic lower bound for M(γ·,·
·,·
·,·
t∈NT En (p̂t ({xi γt }NT ) − di )xi δ̂t
may occur, but is not necessary, and δ̂·,· may not belong to the cone of the restricted eigenvalues,
but obeys the sparse eigenvalue constraints.


P
t ∗0
∗ +δ̂ )−M(γ ∗ )−
∗0 ∗
We first give a suitable upper bound for M(γ·,·
·,·
·,·
t∈NT En (p̂t ({xi γt }NT ) − di )xi δ̂t .
By the Cauchy-Schwarz inequality and the definition of the sparse eigenvalues of Eqn. (16),
r
v
uX X
2
∗ u
δ̂·,·
≤
ŜD ∪ SD
δ̂t,j
t
2,1

t∈NT j∈ŜD ∪S ∗
D

r
≤

s
∗
ŜD ∪ SD

X

n
o−2
∗
φ Q, ŜD ∪ SD
δ̂t0 Qδ̂t

t∈NT

42

r

n
o−1
∗ φ Q, Ŝ ∪ S ∗
En [k{x∗i 0 δ̂t }NT k22 ]1/2 .
ŜD ∪ SD
D
D

=

(B.16)

Following identical steps to Equations (B.1), (B.2), and (B.3), but with δ̂·,· in place of δ̃·,· , and then
using the above bound, we have
√


λD
+ bds TEn [k{x∗i 0 δ̂t }NT k22 ]1/2
En (p̂t ({x∗i 0 γt∗ }NT ) − dti )x∗i 0 δ̂t ≤
δ̂·,·
2
2,1
t∈NT
q


∗|
|
Ŝ
∪
S
√
D
λD
D
≤
+ bds T  En [k{x∗i 0 δ̂t }NT k22 ]1/2 .
∗}
2 φ{Q, ŜD ∪ SD
X

(B.17)
∗ + δ̂ ) − M(γ ∗ ). By optimality of the post selection estimator M(γ̂ ) ≤
Next we turn to M(γ·,·
·,·
·,·
·,·
∗ + δ̂ ) − M(γ ∗ ) ≤ M(γ̃ ) − M(γ ∗ ). By the
M(γ̃·,· ), as S̃ D ⊂ ŜD by construction, and hence M(γ·,·
·,·
·,·
·,·
·,·
mean value theorem, for scalars {mt ∈ [0, 1]}NT we have
i
h
X
∗
∗
M(γ·,·
+ δ̃·,· ) − M(γ·,·
)=
En (dti − p̂t ({x∗i 0 γt∗ + mt x∗i 0 δ̃t }))x∗i 0 δ̃t
t∈NT

h
i
En (dti − p̂t ({x∗i 0 γt∗ }NT ))x∗i 0 δ̃t

X

=

t∈NT

+

X

h
i
En (p̂t ({x∗i 0 γt∗ }NT ) − p̂t ({x∗i 0 γt∗ + mt x∗i 0 δ̃t }))x∗i 0 δ̃t ,

t∈NT

i
h
X
√
λD
δ̃·,·
+ bds TEn [k{x∗i 0 δ̃t }NT k22 ]1/2 +
En mt (x∗i 0 δ̃t )2 .
2
2,1
t∈NT
q


∗|
|ŜD ∪ SD
√
λD
d
≤
T  En [k{x∗i 0 δ̃t }NT k22 ]1/2 + En [k{x∗i 0 δ̃t }NT k22 ],
+
b
s
∗}
2 φ{Q, ŜD ∪ SD

≤

(B.18)
where the first inequality follows from Eqn. (B.3) and the same steps as in (B.15) while the second
applies (B.16) with δ̃·,· and mt ≤ 1.25
Collecting the bounds of (B.17) and (B.18), and the definition of RM gives
∗
∗
M(γ·,·
+ δ̂·,· ) − M(γ·,·
)−

X



En (p̂t ({x∗i 0 γt∗ }NT ) − dti )x∗i 0 δ̂t

t∈NT


≤

q

∗|
|ŜD ∪ SD

λD
∗}
2 φ{Q, ŜD ∪ SD




√
2
+ bds T  En [k{x∗i 0 δ̂t }NT k22 ]1/2 + RM + RM
. (B.19)

25

Applying the steps of Eqn. (B.16) to δ̃·,· is preferred to using the results of Lemma B.3 because it leads to the
∗
tidier expression involving φ{Q, ŜD ∪ SD
}, but the latter method could be substituted.

43

Next, we turn to a lower bound. Consider the same two cases as in the proof of Theorem 5. In
the first case, we have the quadratic lower bound:
∗
∗
M (γ·,·
+δ̂·,· )−M(γ·,·
)−

X
t∈NT



 T En [k{x∗i 0 δ̂t }NT k22 ]
. (B.20)
En (p̂t ({x∗i 0 γt∗ }NT ) − dti )x∗i 0 δ̂t ≥ pmin Ap
TAK

In the other case, this bound may not hold. Arguing as in the proof of Theorem 5, but applying
Eqn. (B.16), we get
√ q
∗ |φ{Q, Ŝ ∪ S ∗ }−1 E [k{x∗ 0 δ̃ } k2 ]1/2 .
k{x∗i 0 δ̂t }NT k1 ≤ TX |ŜD ∪ SD
n
D
D
i t NT 2
Therefore, as above, we find
∗
M(γ·,·
+ δ̂·,· ) − M(γ·,· ) −

X



En (p̂t ({x∗i 0 γt∗ }NT ) − dti )x∗i 0 δ̂t ≥ rn En [k{x∗i 0 δ̂t }NT k22 ]1/2 ,

(B.21)

t∈NT

with

3
rn =
2



2
1−
AK



∗}
φ{Q, ŜD ∪ SD
q
.
√
∗|
X T |ŜD ∪ SD

Collecting the upper bound of (B.19) and the lower bounds (B.20) and (B.21) we have
(

)
o
n
 T 1 En [k{x∗i 0 δ̂t }NT k22 ]
pmin Ap
∧ rn En [k{x∗i 0 δ̂t }NT k22 ]1/2
T
AK
q


∗|


|
Ŝ
∪
S
√
D
λ
D
D
d
∗0
2 1/2
2

≤
+
b
T
E
[k{x
δ̂
}
k
]
+
R
n
M + RM . (B.22)
s
i t NT 2
∗}
2 φ{Q, ŜD ∪ SD

Suppose the linear term is the minimum. The restrictions on AK imply, algebraically, that Eqn.
(B.22) yields


2
rn En [k{x∗i 0 δ̂t }NT k22 ]1/2 ≤ (rn /3) En [k{x∗i 0 δ̂t }NT k22 ]1/2 + RM + RM


≤ (rn /3) En [k{x∗i 0 δ̂t }NT k22 ]1/2 + 2RM .
Canceling the rn and solving yields En [k{x∗i 0 δ̂t }NT k22 ]1/2 ≤ RM . On the other hand, if the quadratic
term is the minimum, define


q
√

T

0
−1
∗
d
∗
RM = Ap pmin TAK 2 λD |ŜD ∪ SD | φ{Q, ŜD ∪ SD } + bs T .
With this notation and the quadratic term being the minimum, Eqn. (B.22) becomes

T
0
0
2
En [k{x∗i 0 δ̂t }NT k22 ] ≤ RM
En [k{x∗i 0 δ̂t }NT k22 ]1 /2 + RM
RM + Ap pmin TAK RM
.

44

Then, because a2 ≤ ab + c implies that a ≤ b +

√

c, we have


1/2

T
0
0
2
En [k{x∗i 0 δ̂t }NT k22 ]1/2 ≤ RM
+ RM
RM + Ap pmin TAK RM
.
Combining the bounds on En [k{x∗i 0 δ̂t }NT k22 ]1/2 from the two cases gives


1/2 

T
0
0
2
.
En [k{x∗i 0 δ̂t }NT k22 ]1/2 ≤ {RM } ∨ RM
+ RM
RM + Ap pmin TAK RM
From this bound on the log-odds, we bound the propensity score and the `1 rate:

1/2 


T
2
0
0
+ bds ;
max En [(p̂t ({x∗i 0 γ̂t }NT ) − pt (xi ))2 ]1/2 ≤ {RM } ∨ RM
+ RM
RM + Ap pmin TAK RM
t∈NT

max kγ̂t −
t∈NT

γt∗ k1

≤

∗|
|S̃ D ∪ SD
∗}
φ{Q, S̃ D ∪ SD

!1/2


{RM } ∨

0
RM

+



0
RM
RM

1/2 

T
2
+ Ap pmin TAK RM
,

by arguments parallel to those used in the proof of Theorem 5.

Appendix C

Proofs for Group Lasso Selection and Estimation of
Linear Models

SEE SUPPLEMENTAL APPENDIX.

9

References

Abadie, A., 2005, Semiparametric difference-in-differences estimators. Review of Economic Studies, 72, 1–19.
Abadie, A., and G. W. Imbens, 2006, Large Sample Properties of Matching Estimators for Average Treatment
Effects. Econometrica, 74, 235–267.
Andrews, D. W. K., and P. Guggenberger, 2009, Incorrect asymptotic size of subsampling procedures based
on post-consistent model selection estimators. Journal of Econometrics, 152, 19–27.
Bach, F. R., 2008, Consistency of the Group Lasso and Multiple Kernel Learning. Journal of Machine
Learning Research, 9, 1179–1225.
, 2010, Self-concordant analysis for logistic regression. Electronic Journal of Statistics, 4, 384–414.
Bang, H., and J. M. Robins, 2005, Doubly Robust Estimation in Missing Data and Causal Inference Models.
Biometrics, 61, 962–972.
Belloni, A., D. Chen, V. Chernozhukov, and C. Hansen, 2012, Sparse models and methods for optimal
instruments with an application to eminent domain. Econometrica, 80, 2369–2429.
Belloni, A., and V. Chernozhukov, 2011, `1 -Penalized quantile regression in high-dimensional sparse models.
The Annals of Statistics, 39, 82–130.
, 2013, Least Squares After Model Selection in High-dimensional Sparse Models. Bernoulli, 19, 521–
547.

45

Belloni, A., V. Chernozhukov, D. Chetverikov, and K. Kato, 2015, Some New Asymptotic Theory for Least
Squares Series: Pointwise and Uniform Results. Journal of Econometrics, forthcoming.
Belloni, A., V. Chernozhukov, I. Fernandez-Val, and C. Hansen, 2014, Program Evaluation with HighDimensional Data. Arxiv preprint arXiv:1311:2645.
Belloni, A., V. Chernozhukov, and C. Hansen, 2014, Inference on Treatment Effects after Selection Amongst
High-Dimensional Controls. Review of Economic Studies, 81, 608–650.
Belloni, A., V. Chernozhukov, and Y. Wei, 2013, Honest Confidence Regions for Logistic Regression with a
Large Number of Controls. arXiv:1304.3969.
Berk, R., L. Brown, A. Buja, K. Zhang, and L. Zhao, 2013, Valid Post-Selection Inference. The Annals of
Statistics, 4, 802–837.
Bickel, P. J., Y. Ritov, and A. B. Tsybakov, 2009, Simultaneous Analysis of LASSO and Dantzig Selector.
The Annals of Statistics, 37, 1705–1732.
Buhlmann, P., and S. van de Geer, 2011, Statistics for High-Dimensional Data, Springer Series in Statistics.
Springer-Verlag, Berlin.
Cattaneo, M. D., 2010, Efficient Semiparametric Estimation of Multi-valued Treatment Effects under Ignorability. Journal of Econometrics, 155, 138–154.
Cattaneo, M. D., R. K. Crump, and M. Jansson, 2013, Generalized Jackknife Estimators of Weighted Average
Derivatives. Journal of the American Statistical Association, 108, 1243–1256.
Cattaneo, M. D., D. M. Drukker, and A. D. Holland, 2013, Estimation of multivalued treatment effects
under conditional independence. The Stata Journal, 13, 407–450.
Cattaneo, M. D., and M. H. Farrell, 2011, Efficient Estimation of the Dose Response Function under Ignorability using Subclassification on the Covariates. In: D. Drukker, (Eds.), Advances in Econometrics:
Missing Data Methods, vol. 27A. Emerald Group Publishing Limited, pp. 93–127.
, 2013, Optimal Convergence Rates, Bahadur Representation, and Asymptotic Normality of Partitioning Estimators. Journal of Econometrics, 174, 127–143.
Cattaneo, M. D., M. Jansson, and W. K. Newey, 2014a, Alternative asymptotics and the partially linear
model with many regressors. working paper.
, 2014b, Small Bandwidth Asymptotics for Density-Weighted Average Derivatives. Econometric
Theory, 30, 176–200.
Chen, X., 2007, Large Sample Sieve Estimation of Semi-Nonparametric Models. In: J. Heckman, and
E. Leamer, (Eds.), Handbook of Econometrics, vol. 6B of Handbook of Econometrics. Elsevier, chap. 76.
Chen, X., and T. M. Christensen, 2015, Optimal uniform convergence rates and asymptotic normality for
series estimators under weak dependence and weak conditions. Journal of Econometrics, forthcoming.
Chen, X., H. Hong, and A. Tarozzi, 2004, Semiparametric Efficiency in GMM Models of Nonclassical Measurament Errors, Missing Data and Treatment Effects. Cowles Foundation Discussion Paper No. 1644.
, 2008, Semiparametric Efficiency in GMM Models With Auxiliary Data. The Annals of Statistics,
36, 808–843.
de la Peña, V. H., T. L. Lai, and Q.-M. Shao, 2009, Self-Normalized Processes: Limit Theory and Statistical
Applications, Probability and Its Applications. Springer.
Dehejia, R. H., and S. Wahba, 1999, Causal Effects in Nonexperimental Studies: Reevaluating the Evaluation
of Training Programs. Journal of the American Statistical Association, 94, 1053–1062.

46

, 2002, Propensity Score-Matching Methods for Nonexperimental Causal Studies. The Review of
Economics and Statistics, 84, 151–161.
Efron, B., 2014, Estimation and Accuracy after Model Selection. Journal of the American Statistical Association, 109, 991–1007.
Farrell, M. H., 2015a, Robust Inference on Average Treatment Effects with Possibly More Covariates than
Observations. Journal of Econometrics, 189, 1–23.
, 2015b, Supplement to ‘Robust Inference on Average Treatment Effects with Possibly More Covariates than Observations’. Journal of Econometrics, 189.
Hahn, J., 1998, On the Role of the Propensity Score in Efficient Semiparametric Estimation of Average
Treatment Effects. Econometrica, 66, 315–331.
, 2004, Functional restriction and efficiency in causal inference. Review of Economics and Statistics,
84, 73–76.
He, X., and Q.-M. Shao, 2000, On Parameters of Increasing Dimensions. Journal of Multivariate Analysis,
73, 1201–35.
Heckman, J., H. Ichimura, and P. Todd, 1997, Matching as an Econometric Evaluation Estimator: Evidence
from Evaluating a Job Training Programme. Review of Economic Studies, 64, 605–654.
Heckman, J., and E. J. Vytlacil, 2007, Econometric Evaluation of Social Programs, Part I. In: J. Heckman,
and E. Leamer, (Eds.), Handbook of Econometrics, vol. VIB. Elsevier Science B.V., pp. 4780–4874.
Hirano, K., G. W. Imbens, and G. Ridder, 2003, Efficient Estimation of Average Treatment Effects using
the Estimated Propensity Score. Econometrica, 71, 1161–1189.
Holland, P. W., 1986, Statistics and Causal Inference. Journal of the American Statistical Association, 81,
945–960.
Horowitz, J. L., and C. F. Manski, 2000, Nonparametric Analysis of Randomized Experiments With Missing
Covariate and Outcome Data. Journal of the American Statistical Association, 95, 77–84.
Huang, J., and T. Zhang, 2010, The Benefit of Group Sparsity. The Annals of Statistics, 38, 1978–2004.
Huang, J. Z., 2003, Local asymptotics for polynomial spline regression. The Annals of Statistics, 31, 1600–
1635.
Imai, K., and D. A. van Dyk, 2004, Causal Inference With General Treatment Regimes: Generalizing the
Propensity Score. Journal of the American Statistical Association, 99, 854–866.
Imbens, G. W., 2000, The Role of the Propensity Score in Estimating Dose-Response Functions. Biometrika,
87, 706–710.
, 2004, Nonparametric Estimation of Average Treatment Effects under Exogeneity: A Review. Review
of Economics and Statistics, 86, 4–29.
Imbens, G. W., W. K. Newey, and G. Ridder, 2007, Mean-Squared-Error Calculations for Average Treatment
Effects. working paper.
Imbens, G. W., and J. M. Wooldridge, 2009, Recent Developments in the Econometrics of Program Evaluation. Journal of Economic Literature, 47, 5–86.
Kang, J. D. Y., and J. L. Schafer, 2007, Demystifying Double Robustness: A Comparison of Alternative
Strategies for Estimating a Population Mean from Incomplete Data. Statistical Science, 22, 523–539.
Kolar, M., J. Lafferty, and L. Wasserman, 2011, Union Support Recovery in Multi-task Learning. Journal of
Machine Learning Research, 12, 2415–2435.

47

Kwemou, M., 2012, Non-asymptotic Oracle Inequalities for the Lasso and Group Lasso in high dimensional
logistic model. Arxiv preprint arXiv:1206.0710.
LaLonde, R. J., 1986, Evaluating the Econometric Evaluations of Training Programs with Experimental
Data. American Economic Review, 76, 604–620.
Lechner, M., 2001, Identification and estimation of causal effects of multiple treatments under the conditional
independence assumption. In: M. Lechner, and E. Pfeiffer, (Eds.), Econometric Evaluations of Active
Labor Market Policies. Physica, Heidelberg, pp. 43–58.
Leeb, H., and B. M. Pötscher, 2005, Model Selection and Inference: Facts and Fiction. Econometric Theory,
21, 21–59.
, 2008a, Can One Estimate the Unconditional Distribution of Post-Model-Selection Estimators?.
Econometric Theory, 24, 338–376.
, 2008b, Sparse estimators and the oracle property, or the return of Hodges’ estimator. Journal of
Econometrics, 142, 201–211.
Lounici, K., M. Pontil, S. van de Geer, and A. B. Tsybakov, 2011, Oracle Inequalities and Optimal Inference
under Group Sparsity. The Annals of Statistics, 39, 2164–2204.
Negahban, S. N., P. Ravikumar, M. J. Wainwright, and B. Yu, 2012, A Unified Framework for HighDimensional Analysis of M -Estimators with Decomposable Regularizers. Statistical Science, 27, 538–557.
Newey, W. K., 1990, Efficient Instrumental Variables Estimation of Nonlinear Models. Econometrica, 58,
809–837.
, 1997, Convergence rates and asymptotic normality for series estimators. Journal of Econometrics,
79, 147–168.
Newey, W. K., and D. L. McFadden, 1994, Large sample estimation and hypothesis testing. In: R. F. Engle,
and D. McFadden, (Eds.), Handbook of Econometrics, vol. 4 of Handbook of Econometrics. Elsevier,
chap. 36, pp. 2111–2245.
Obozinski, G., M. J. Wainwright, and M. I. Jordan, 2011, Support Union Recovery in High-Dimensional
Multivariate Regression. The Annals of Statistics, 39, 1–47.
Pötscher, B. M., 2009, Confidence Sets Based on Sparse Estimators Are Necessarily Large. Sankhyā, 71-A,
1–18.
Pötscher, B. M., and H. Leeb, 2009, On the distribution of penalized maximum likelihood estimators: The
LASSO, SCAD, and thresholding. Journal of Multivariate Analysis, 100, 2065–2085.
Powell, J. L., J. H. Stock, and T. M. Stoker, 1989, Semiparametric Estimation of Index Coefficients. Econometrica, 57, 1403–1430.
Raskutti, G., M. J. Wainwright, and B. Yu, 2010, Restricted Eigenvalue Properties for Correlated Gaussian
Designs. Journal of Machine Learning Research, 11, 2241–2259.
Robins, J., L. Li, E. Tchetgen, and A. van der Vaart, 2008, Higher order influence functions and minimax
estimation of nonlinear functionals. In: D. Nolan, and T. Speed, (Eds.), Probability and Statistics: Essays
in Honor of David A. Freedman, vol. 2. Beachwood, Ohio, USA: Institute of Mathematical Statistics.
Robins, J. M., and A. Rotnitzky, 1995, Semiparametric Efficiency in Multivariate Regression Models with
Missing Data. Journal of the American Statistical Association, 90, 122–129.
Romano, J. P., 2004, On non-parametric testing, the uniform behaviour of the t-test, and related problems.
Scandinavian Journal of Statistics, 31, 567–584.

48

Rosenbaum, P. R., and D. B. Rubin, 1983, The Central Role of the Propensity Score in Observational Studies
for Causal Effects. Biometrika, 70, 41–55.
Rudelson, M., and S. Zhou, 2013, Reconstruction from anisotropic random measurements. IEEE Transactions
on Information Theory, 59, 3434–3447.
Smith, J. A., and P. E. Todd, 2005, Does matching overcome LaLonde’s critique of nonexperimental estimators?. Journal of Econometrics, 125, 305–353.
Tan, Z., 2010, Bounded, efficient and doubly robust estimation with inverse weighting. Biometrika, 97,
661–682.
Tanabe, K., and M. Sagae, 1992, An Exact Cholesky Decomposition and the Generalized Inverse of the
Variance-Covariance Matrix of the Multinomial Distribution, with Applications. Journal of the Royal
Statistical Society. Series B (Methodological), 54, 211–219.
Tsiatis, A. A., 2006, Semiparametric Theory and Missing Data. Springer, New York.
van de Geer, S., 2008, High-Dimensional Generalized Linear Models and the Lasso. The Annals of Statistics,
36, 614–645.
van de Geer, S., and P. Buhlmann, 2009, On the conditions used to prove oracle results for the Lasso.
Electronic Journal of Statistics, 3, 1360–1392.
van de Geer, S., P. Buhlmann, Y. Ritov, and R. Dezeure, 2014, On Asymptotically Optimal Confidence
Regions and Tests for High-Dimensional Models. The Annals of Statistics, 42, 1166–1202.
van der Laan, M., and J. M. Robins, 2003, Unified Methods for Censored Longitudinal Data and Causality.
Springer-Verlag.
Vincent, M., and N. R. Hansen, 2014, Sparse group lasso and high dimensional multinomial classification.
Computational Statistics and Data Analysis, 71, 771–786.
von Bahr, B., and C.-G. Esseen, 1965, Inequalities for the rth absolute moment of a sum of random variables,
1 5 r 5 2. Annals of Mathematical Statistics, 36, 299–303.
Wei, F., and J. Huang, 2010, Consistent group selection in high-dimensional linear regression. Bernoulli, 16,
1369–1384.
White, H., and X. Lu, 2011, Causal Diagrams for Treatment Effect Estimation with Application to Efficient
Covariate Selection. Review of Economics and Statistics, 93, 1453–1459.
Wooldridge, J. M., 2007, Inverse probability weighted estimation for general missing data problems. Journal
of Econometrics, 141, 1281–1301.
, 2010, Econometric Analysis of Cross Section and Panel Data. MIT Press, Cambridge, 2 edn.
Yuan, M., and Y. Lin, 2006, Model selection and estimation in regression with grouped variables. Journal of
the Royal Statistical Society. Series B, 68, 46–67.
Zhang, C.-H., and S. S. Zhang, 2014, Confidence intervals for low dimensional parameters in high dimensional
linear models. Journal of the Royal Statistical Society. Series B, 76, 217–242.
Zou, H., 2006, The Adaptive Lasso and Its Oracle Properties. Journal of the American Statistical Association,
101, 1418–1429.

49

N/A
??
171

–

Before
selection(a)

11
15
20/6

–

After
selection(b)

1211
1058
1735

260

Control

185
185
185

185

Treated

Sample Sizes(c)

1664
2528
1737

1794

ATT

[-276, 3604]
[149, 4908]
[33, 3441]

[110, 3479]

95% CI

Notes: All analyses use the DW99 subsample and PSID comparison group. Specifications vary, but all estimates and standard
errors of from the method defined in Section 5 with the exception of the partially linear model.
(a) Not counting the intercept. The total set of variables considered by DW02 is not known.
(b) For the group lasso estimators, the two numbers given are for those used in the outcome regressions and propensity score,
respectively. For other doubly-robust estimators all variables are used in the propensity score and outcome models.
(c) The full sample begins with 2490 comparisons and 185 treated units. Control observations outside the range of estimated
propensity scores in the treated sample are discarded.

Doubly-Robust Estimates
Specification 1 (No Selection)
DW02 (Informal Selection)
Refitting after Group Lasso Selection

Experimental Benchmark

Specifications:

Number of Variables

Table 1: Analysis of NSW Demonstration: Treatment Effects on the Treated and Confidence Intervals for
Various Specifications

Figure 1: Empirical Coverage of 95% Confidence Intervals, Varying Signal Strength and Sparsity
of pt (x) and µt (x)

µ(x) signal

µ(x) sparsity

pt(x) signal

1.0

pt(x) sparsity

1.0

0.5

0.5

0.0

0.0

(a) p = 1000, varying signal strength

µ(x) signal

(b) p = 1000, varying sparsity

µ(x) sparsity

pt(x) signal

1.0

pt(x) sparsity

1.0

0.5

0.5

0.0

0.0

(c) p = 1500, varying signal strength

(d) p = 1500, varying sparsity

Figure 2: Empirical Coverage of 95% Confidence Intervals, Penalty Chosen with Cross-validation,
Varying Signal Strength of pt (x) and µt (x)

µ(x) signal

1.0

µ(x) signal

pt(x) signal
1.0

0.5

0.5

0.0

0.0

(a) 1000 Covariates

(b) 1500 Covariates

pt(x) signal

Supplemental Appendix for “Robust Inference on Average Treatment
Effects with Possibly More Covariates than Observations”
Max H. Farrell
March 27, 2015
Updated: February 1, 201826

Summary
This is a supplemental appendix for “Robust Inference on Average Treatment Effects with Possibly More
Covariates than Observations” containing complete proofs. Notation is kept in line with the main text,
though equation numbers may change. This file may serve as a drop-in replacement for the main appendix.
Additional simulations results are also presented.

26

The published version of the supplement, Farrell (2015b), contains an error in the proof which is rectified here.
See the author’s website for further detail. I am grateful to Whitney Newey for alerting me to this error.

Appendix D

Proofs for Treatment Effect Inference

The proofs in this section are asymptotic in nature, compared to the nonasymptotic bounds of the next
section. It shall be understood that asymptotic order symbols hold for the sequence being considered, as
a shorthand for the more formal versions given in the assumptions (e.g. Assumption 3). C will denote a
generic positive constant, which may be a matrix. Define the set of indexes It = {i : di = t}.
Proof of Theorem 2. First, we have µ̂t − µt = En [ψt (yi , dti , µ0t (xi ), p0t (xi ), µt )] + R1 + R2 + R3 + R4 , where


dti yi
0
R1 = En [µ̂t (xi ) − µt (xi )],
R2 = En
(p
(x
)
−
p̂
(x
))
,
i
t i
p̂t (xi )p0t (xi ) t
 t

 t 0

di
di µt (xi )
0
(µ̂t (xi ) − µ0t (xi )) ,
(p
R3 = −En
and
R4 = −En
(x
)
−
p̂
(x
))
.
i
t i
p̂t (xi )
p̂t (xi )p0t (xi ) t
Under Assumptions 1(b) and 2, if follows from the Cauchy-Schwarz inequality, the consistency condition of
the theorem and the inequality of von Bahr and Esseen (1965), it follows that |Rk | = oPn (1) for k = 1, 2, 3, 4.
Finally, by the same von Bahr and Esseen inequality |En [ψt (yi , dti , µ0t (xi ), p0t (xi ), µt )]| = oPn (1) if either µt (x)
or pt (x) is correctly specified.
We first prove Theorem 3.1 assuming there is no additional randomness injected into the support estimates. Following this, we redo the proof to account for additional randomness. We then turn to the
remaining portions of Theorem 3 and to Corollary 2, which require shorter arguments.
We make frequent use of the linearization
1
1 b−a
1 b − a (b − a)2
= +
= + 2 +
,
a
b
ab
b
b
ab2

(D.1)

where the first inequality is readily verified, and the second re-applies the first.
Proof √
of Theorem 3.1 without Additional Randomness. With ψt (·) defined in Eqn. (2), we have
µt ) = nEn [ψt (yi , dti , µt (xi ), pt (xi ), µt )] + R1 + R2 , where
n

1 X t
R1 = √
d (yi − µt (xi ))
n i=1 i
and



1
1
−
p̂t (xi ) pt (xi )

√

n(µ̂t −





n
dti
1 X
(µ̂t (xi ) − µt (xi )) 1 −
R2 = √
.
p̂t (xi )
n i=1

The proof proceeds by showing that both R1 and R2 are oPn (1).
For R1 , applying the first equality in Eqn. (D.1), we rewrite R1 as
n

1 X t
R1 = √
d ui
n i=1 i



pt (xi ) − p̂t (xi )
p̂t (xi )pt (xi )


.

Applying Assumptions 1(b) and 2(c) and the first-stage consistency condition of Assumption 3(a):




dti σt2 (xi )
2
≤ CEn [(pt (xi ) − p̂t (xi ))2 ] = oPn (1).
E R12 |{xi , di }ni=1 = En
(p
(x
)
−
p̂
(x
))
t i
t i
p̂t (xi )2 pt (xi )2

Next, again using Eqn. (D.1) we have
1−

dti
pt (xi ) − dti
dt (p̂t (xi ) − pt (xi ))
=
+ i
.
p̂t (xi )
pt (xi )
p̂t (xi )pt (xi )

We use this to re-write R2 = R21 + R22 , where
n

R21

1 X
=√
(µ̂t (xi ) − µt (xi ))
n i=1



pt (xi ) − dti
pt (xi )



and



n
1 X
dti
R22 = √
(µ̂t (xi ) − µt (xi ))(p̂t (xi ) − pt (xi ))
.
p̂t (xi )pt (xi )
n i=1
√
For the first term, R21 = nEn [(µ̂t (xi ) − µt (xi ))(1 − dti /pt (xi ))] = oPn (1) by Assumption 3(c). Next, by
Hölder’s inequality, Assumption 1(b) and the rate condition of Assumption 3(b)


p
√
1
|R22 | ≤ n max
En [(µ̂t (xi ) − µt (xi ))2 ]En [(p̂t (xi ) − pt (xi ))2 ]
i≤n p̂t (xi )pt (xi )
√ p
= OPn (1) n En [(µ̂t (xi ) − µt (xi ))2 ]En [(p̂t (xi ) − pt (xi ))2 ] = oPn (1).
This completes the proof, as |R1 + R2 | = oPn (1) by Markov’s inequality and the triangle inequality.

Proof of Theorem 3.1 with Additional Randomness. We must reconsider the remainders R1 and R2 . For the
former, applying Eqn. (D.1), we find R1 = R11 + R12 , where
n

R11

1 X dti ui
(pt (xi ) − p̂t (xi ))
=√
n i=1 pt (xi )2

n

and

R12

1 X
dti ui
2
=√
(p̂t (xi ) − pt (xi )) .
n i=1 pt (xi )2 p̂t (xi )

For R11 , we first add and subtract the parametric representation to get R11 = R111 + R112 , where,
n

R111


1 X dti ui
p̂t ({x∗i 0 γt∗ }NT ) − p̂t ({x∗i 0 γ̂t }NT )
=√
2
n i=1 pt (xi )

n

and

R112


1 X dti ui
=√
pt (xi ) − p̂t ({x∗i 0 γt∗ }NT ) .
2
n i=1 pt (xi )

By a two-term mean-value expansion R111 = R111a + R111b , with
n

1 X dti ui X 
R111a = √
p̂t ({x∗i 0 γt∗ }NT )(1 − p̂t ({x∗i 0 γt∗ }NT )) x∗i 0 (γ̂t − γt∗ )
2
n i=1 pt (xi )
t∈NT

and

n

1 X dti ui 0
v H̄vi ,
R111b = √
2 n i=1 pt (xi )2 i
where vi = {x∗i 0 (γ̂t − γt∗ )}NT and H = H({x∗i 0 γt∗ + mt x∗i 0 γ̂t }NT ) for appropriate scalars mt where the T-square
matrix H({x∗i 0 γt }NT ) is defined as having the (t, t0 ) ∈ N2T entry given by
(
p̂t ({x∗i 0 γt }NT )(1 − p̂t ({x∗i 0 γt }NT )) if t = t0
∗0
H({xi γt }NT )[t,t0 ] =
−p̂t ({x∗i 0 γt }NT )p̂t0 ({x∗i 0 γt }NT )
if t 6= t0
P
For R111a , consider each term in the sum over NT one at a time; let R111a = t∈NT R111a (t). Let t0 denote
 ∗ 2 2

the original treatment under consideration. Define Σt,j = E (xi,j ) σt0 (xi )p̂t ({x∗i 0 γt∗ }NT )2 (1 − p̂t ({x∗i 0 γt∗ }NT ))2 /pt0 (xi )3 .

Then proceed as follows
!
0
n
1 X dti ui p̂t ({x∗i 0 γt∗ }NT )(1 − p̂t ({x∗i 0 γt∗ }NT )) X ∗
xi,j (γ̂t − γt∗ )
R111a (t) = √
pt0 (xi )2
n i=1
j∈ŜD
(
!)
0
n
t
∗0 ∗
∗0 ∗
X
1 X
1/2
∗ di ui p̂t ({xi γt }NT )(1 − p̂t ({xi γt }NT ))
∗
√
=
xi,j
Σt,j (γ̂t,j − γt,j
)
1/2
2
n i=1
0 (xi ) Σ
p
t
t,j
j∈ŜD
!


0
n
1 X ∗ dti ui p̂t ({x∗i 0 γt∗ }NT )(1 − p̂t ({x∗i 0 γt∗ }NT ))
1/2
kγ̂t − γt∗ k1
≤ max Σt,j
max √
x
1/2
j∈Np
j∈Np
n i=1 i,j
pt0 (xi )2 Σ
t,j

= O(1)OPn (log(p)) kγ̂t − γt∗ k1 = oPn (1).
Convergence follows under Assumption 5. For the penultimate equality, it follows from Assumptions 1(b),
2(b), and 2(c) that maxj∈Np Σt,j = O(1). Finally, the center factor is shown to be OPn (log(p)) by applying
the moderate deviation theory for self-normalized sums of de la Peña, Lai, and Shao (2009, Theorem 7.4)
and in particular Belloni, Chen, Chernozhukov, and Hansen (2012, Lemma 5). To apply this lemma, first
note that the summand of the center factor has bounded third moment and second moment bounded away
from zero, from Assumptions 1(b), 2(b), 2(c), and the requirements of Assumptions 3 and 5. Σt,j normalizes
the second moment, and the lemma applies under Assumptions 4 and the first restriction of Assumption 5.
For R111b , the results of Tanabe and Sagae (1992) coupled with Assumption 3 give vi0 H̄vi ≤ Ckvi k22 .
Thus, using Assumption 1(b) to bound maxi≤n pt (xi )−2 < C, we find R111b may be bounded as follows:
X√


|R111b | ≤ C
n(max |ui |)En |x∗i 0 (γ̂t − γt∗ )|2
i∈It

t∈NT

√

≤ CT max
t∈NT



n(max |ui |)En |p̂t ({x∗i 0 γ̂t }NT ) − p̂t ({x∗i 0 γt∗ }NT )|2 = oPn (1),
i∈It

by the union bound and Assumption 5, using the Assumptions 1(b) and 3(a) to apply Eqn. (E.19) with the
inequality reversed.
A variance bound may be applied to R112 as in the previous proof, and we have |R112 | = OPn (bs ) = oPn (1)
by Markov’s inequality.
Next, R12 is simply bounded by


i
h
√
1
2
|R12 | ≤ n(max |ui |) max
E
(p̂
(x
)
−
p
(x
))
n
t
i
t
i
i∈It
i∈It pt (xi )2 p̂t (xi )
h
i
√
2
≤ OPn (1) n(max |ui |)En (p̂t (xi ) − pt (xi )) = oPn (1),
i∈It

where the rate follows from Assumptions 1(b), 2, and 3, and this tends to zero by Assumption 5.
As in the prior proof, write R2 = R21 + R22 . The same bound is used for R22 . However, for R21 , add
and subtract the pseudotrue values to get R21 = R211 + R212 , where
n

R211

1 X ∗0
=√
(x β̂t − x∗i βt∗ )
n i=1 i



pt (xi ) − dti
pt (xi )

n


and

R212

1 X ∗ ∗
=√
(x β − µt (xi ))
n i=1 i t





For the first term, define Σ̃t,j = E (x∗i,j )2 (dti − pt (xi ))2 /pt (xi )2 and then proceed as follows:
n

R211

1 X
=√
n i=1



pt (xi ) − dti
pt (xi )

 X
j∈ŜY

∗
x∗i,j (β̂t,j − βt,j
)

pt (xi ) − dti
pt (xi )



(
=

X
j∈ŜY

n

1 X x∗i,j (pt (xi ) − dti )/pt (xi )
√
1/2
n i=1
Σ̃

)
1/2

∗
Σ̃t,j (β̂t,j − βt,j
)

t,j

!


n
1 X x∗i,j (pt (xi ) − dti )/pt (xi )
1/2
β̂t − βt∗
≤ max Σ̃t,j
max √
1/2
j∈Np
j∈Np
n i=1
Σ̃
t,j

= O(1)OPn (log(p)) β̂t − βt∗

1

1

= oPn (1),

where the final line follows exactly as above.
A variance bound may be applied to R212 as in the previous proof, and we have |R212 | = OPn (bs ) = oPn (1)
by Markov’s inequality.
Proof of Theorem 3.2. This claim follows directly from the prior result under the moment conditions of
Assumption 2(e).
Proof of Theorem 3.3. We begin with V̂W (t). Expanding the square and using Eqn. (D.1), rewrite V̂µW (t) =
En [dti u2i pt (xi )−2 ] + RW,1 + RW,2 + RW,3 where


dti u2i
RW,1 = En
(p̂
(x
)
−
p
(x
))
(p̂
(x
)
+
p
(x
))
,
t i
t i
t i
t i
p̂t (xi )2 pt (xi )2


 t
 t
di (µt (xi ) − µ̂t (xi ))2
di ui (µt (xi ) − µ̂t (xi ))
RW,2 = En
,
and
RW,3 = 2En
.
p̂t (xi )2
p̂t (xi )2
Using Hölder’s inequality, Assumptions 1(b), 2(e), and 3(a), we have the following


p̂t (xi ) + pt (xi )
RW,1 ≤ max
En [dti |ui |4 ]1/2 En [dti (p̂t (xi ) − pt (xi ))2 ]1/2 = oPn (1),
i∈It p̂t (xi )2 pt (xi )2


1
En [dti (µ̂t (xi ) − µt (xi ))2 ] = oPn (1),
RW,2 ≤ max
i∈It p̂t (xi )2
and,

RW,3 ≤ 2 max
i∈It

1
p̂t (xi )2



En [dti |ui |2 ]1/2 En [dti (µ̂t (xi ) − µt (xi ))2 ]1/2 = oPn (1),

where En [|ui |4 ] = OPn (1) from the inequality of von Bahr and Esseen (1965). From the same inequality it
follows that En [dti u2i pt (xi )−2 ] − VµW (t)| = oPn (1), under Assumptions 1(b) and 2(c).
Next consider the “between” variance estimator, V̂µB . For any tNT and t0 ∈ NT , define
RB,1 (t, t0 ) = En [(µ̂t (xi ) − µt (xi ))(µ̂t0 (xi ) − µt0 (xi ))] ,
RB,2 (t, t0 ) = µ̂t En [µ̂t0 (xi ) − µt0 (xi )] ,

and RB,3 (t, t0 ) = En [µt (xi )(µ̂t0 (xi ) − µt0 (xi ))] .

From Hölder’s inequality, Assumption 3(a), Theorem 3.2, the von Bahr and Esseen inequality, and Assumptions 2(c) and 2(e) it follows that RB,k (t, t0 ) = oPn (1) for k ∈ N3 and all pairs (t, t0 ) ∈ N2t . With this in
mind, we decompose
V̂µB (t, t0 ) = En [µt (xi )µt0 (xi )] − µ̂t En [µt0 (xi )] − µ̂t0 En [µt (xi )] + µ̂t µ̂t0
+ RB,1 (t, t0 ) + RB,2 (t, t0 ) + RB,2 (t0 , t) + RB,3 (t, t0 ) + RB,3 (t0 , t).
Consistency of V̂µB (t, t0 ) now follows from the von Bahr and Esseen inequality and Theorem 3.2.
Proof of Corollary 2. Suppose the result did not hold. Then, there would exist a subsequence Pm ∈ Pm , for

each m, such that



q
lim PPm G(µ) ∈ G(µ̂) ± cα ∇G (µ̂)V̂ ∇0G (µ̂)/n
− (1 − α) > 0.

m→∞

But this contradicts Theorem 3, under which (∇G (µ̂)V̂ ∇0G (µ̂)/n)−1/2 (G(µ̂) − G(µ)) is asymptotically standard normal under the sequence Pm .

Appendix E

Proofs for Group Lasso Selection and Estimation of
Multinomial Logistic Models

Unless otherwise noted, all bounds in this section are nonasymptotic. We will use generic notation X ∗ , δ, s,
etc, as this section deals only with multinomial logistic models.

E.1

Lemmas

Lemma E.1 (Score Bound). For λD and P defined respectively in Eqn. (11) and Eqn. (12) we have


λD
P max kEn [(pt (xi ) − dti )x∗i,j ]k2 ≥
≤ P.
j∈Np
2
2
Proof. The residuals vt,i = pt (xi ) − dti are conditionally mean-zero by definition and satisfy E[vt,i
|xi ] ≤ 1.
Using this, Assumption 2(a), and the definition of X, we find that

X 
X 1



X2 T
2
E[vt,i
E kEn [vt,i x∗i,j ]k22 =
E En [vt,i x∗i,j ]2 =
(x∗i,j )2 ] ≤
n
n
t∈NT

t∈NT

uniformly in j ∈ Np . Define the mean-zero random variables ξt,j as:
ξt,j = (En [vt,i x∗i,j ])2 −

1
E[Vt2 Xj∗ 2 ].
n

Using this definition and the above bound after inserting the definition of λD , setting rn = T −1/2 log(p ∨
n)3/2+δ , and squaring both sides, we have
"
#
√


λ
X
T
D
1/2
P max kEn [(pt (xi ) − dti )x∗i,j ]k2 ≥
= P max kEn [vt,i x∗i,j ]k2 ≥ √ (1 + rn )
j∈Np
j∈Np
2
n


X2 T
∗
2
(1 + rn )
= P max kEn [vt,i xi,j ]k2 ≥
j∈Np
n


X2 T
X2 Trn
∗
2
= P max kEn [vt,i xi,j ]k2 −
≥
j∈Np
n
n
#
"
X
X2 Trn
≤ P max
ξt,j ≥
j∈Np
n
t∈NT
"
#
X
n
≤ E max
ξt,j
,
(E.1)
j∈Np
X2 Trn
t∈NT

where final line follows from Markov’s inequality.

Next, applying Lemma 9.1 of Lounici, Pontil, van de Geer, and Tsybakov (2011) (with their m = 1 and
hence c(m) = 2) followed by Jensen’s inequality and Assumption 2(c), we find that

#
"
!1/2 
X
X
2

E max
ξt,j ≤ (8 log(2p))1/2 E 
max ξt,j
j∈Np

t∈NT

t∈NT

j∈Np

#!1/2

"
≤ (8 log(2p))

1/2

X

E

t∈NT
1/2

≤ 4 log(2p)

2
max ξt,j
j∈Np

!1/2
X 
X X4
4
∗
.
+
E max En [vt,i xi,j ]
j∈Np
n2

(E.2)

t∈NT

t∈NT

√ √
√
2
The leading 4 is 8 2, where 2 is a byproduct of applying the inequality (a−b)2 ≤ 2(a2 +b2 ) to ξt,j
. Again
using Lemma 9.1 of Lounici, Pontil, van de Geer, and Tsybakov (2011) (with their m = 4, and c(m) = 12
since c(4) ≥ (e4−1 − 1)/2 + 2 ≈ 11.54), we bound the expectation in the second term above as follows:

!4/2 


n
∗ 2
2 4
X
v
x
t,i i,j
4
 ≤ 64 log(12p) X ,
E max En [vt,i x∗i,j ]
(E.3)
≤ [8 log(12p)]4/2 E 
max
2
j∈Np
j∈Np
n
n
i=1
using Assumptions 2(a) and 2(b).
Now, inserting the results of Eqns. (E.2) and (E.3) into Eqn. (E.1), we have


λD
4n log(2p)1/2
P max kEn [vt,i x∗i,j ]k2 ≥
≤
j∈Np
4
TX2 rn
1/2

≤

4 log(2p)
√
rn T

X X4
X 64 log(12p)2 X4
+
2
n
n2

t∈NT

!1/2

t∈NT

[1 + 64 log(12p)2 ]1/2 = P,

using the choice rn = T −1/2 log(p ∨ n)3/2+δ .
Lemma E.2 (Estimate Sparsity). With probability at least 1 − P
|S̃ D | ≤

X


4
φ{Q, S̃ D }
En (p̂t ({x∗i 0 γ̃t }NT ) − pt (xi ))2 .
2
λD
t∈NT

Proof. First, by Karush-Kuhn-Tucker conditions for (9), for all t ∈ NT , if γ̃·,j 6= 0 it must satisfy
En [x∗i,j (p̂t ({x∗i 0 γ̃t }NT ) − dti )] = λD

γ̃t,j
.
kγ̃·,j k2

(E.4)

Hence, taking the `2 -norm over t ∈ NT for fixed j ∈ S̃ D , adding and subtracting the true propensity score,
using the triangle inequality, and the score bound (E.1), we find that
λD = En [x∗i,j (p̂t ({x∗i 0 γ̃t }NT ) − dti )]
≤ En [x∗i,j (pt (xi ) − dti )]

2

+ En [x∗i,j (p̂t ({x∗i 0 γ̃t }NT ) − pt (xi ))]
2

≤ λD /2 + En [x∗i,j (p̂t ({x∗i 0 γ̃t }NT ) − pt (xi ))]

2

.
2

Let P̃t be the vector of {p̂t ({x∗i 0 γ̃t }NT )}ni=1 and Pt collect {pt (xi )}ni=1 . Collecting terms, then squaring both

sides and summing over j ∈ S̃ D (i.e. applying k · k22 over j ∈ S̃ D to both sides) yields
X
X X
λ2D ≤ 4
En [x∗i,j (p̂t ({x∗i 0 γ̃t }NT ) − pt (xi ))]2
j∈S̃ D t∈NT

j∈S̃ D

=4

i
X 1 h
0
X
(
P̃
−
P
)
t
t
n2
j∈S̃ D

t∈NT

2
2

D

X φ{Q, S̃ }
2
P̃t − Pt
n
2,n
t∈NT
X


≤ 4φ{Q, S̃ D }
En (p̂t ({x∗i 0 γ̃t }NT ) − pt (xi ))2 .
≤4

t∈NT

The result now follows, as the left-hand side is equal to |S̃ D |λ2D .
∗
satisfies
Lemma E.3 (Bounds in `2 /`1 norm). With probability 1 − P the vector δ̃·,· = γ̃·,· − γ·,·

δ̃·,·

2,1

( p
√ )
5 |S∗ | 10bds T
≤
∨
En [k{x∗i 0 δ̃t }NT k22 ]1/2 .
κD
λD

and

√ )
|S∗ | 2bds T
En [k{x∗i 0 δ̃t }NT k22 ]1/2 .
∨
κD
λD

(p
δ̃·,S∗

≤
2,1

Proof. By the Cauchy-Schwarz inequality and Lemma E.1,
h
i
X
X X


En (pt (xi ) − dti )x∗i 0 δ̃t =
En (pt (xi ) − dti )x∗i,j δ̃t,j
t∈NT

j∈Np t∈NT

X sX

≤

t∈NT

j∈Np

≤ max

n

sX

oX
2
j∈Np

λD
δ̃·,·
2

2
δ̃t,j

t∈NT



En (pt (xi ) − dti )x∗i,j

j∈Np

≤


2
En (pt (xi ) − dti )x∗i,j

δ̃·,j

2

,

(E.5)

2,1

with probability at least 1 − P. Applying the Cauchy-Schwarz inequality, the bias condition of Assumption
4, and Cauchy-Schwarz again yield
X

h
i
h
i1/2
X

1/2
En (p̂t ({x∗i 0 γt∗ }NT ) − pt (xi ))x∗i 0 δ̃t ≤
En (p̂t ({x∗i 0 γt∗ }NT ) − pt (xi ))2
En (x∗i 0 δ̃t )2
t∈NT

t∈NT

≤ bds

X

h
i1/2
En (x∗i 0 δ̃t )2

t∈NT

√
≤ bds T
=

bds

√

sX

h
i
En (x∗i 0 δ̃t )2

t∈NT

TEn [k{x∗i 0 δ̃t }NT k22 ]1/2 .

(E.6)

Combining Equations (E.5) and (E.6), we have, probability at least 1 − P,
h
i
h
i
X
X
En (p̂t ({x∗i 0 γt∗ }NT ) − dti )x∗i 0 δ̃t =
En (pt (xi ) − dti )x∗i 0 δ̃t
t∈NT

t∈NT

+

X

i
h
En (p̂t ({x∗i 0 γt∗ }NT ) − pt (xi ))x∗i 0 δ̃t

t∈NT

λD
δ̃·,·
≤
2

2,1

√
+ bds TEn [k{x∗i 0 δ̃t }NT k22 ]1/2 .

(E.7)

By the optimality of δ̃·,· , we have
∗
∗
+ δ̃·,·
M(γ·,·
+ δ̃·,· ) + λD γ·,·

2,1

∗
∗
≤ M(γ·,·
) + λD γ·,·

2,1

,

implying

λD

∗
γ·,·

2,1

∗
γ·,·

−


+ δ̃·,·

2,1

∗
∗
≥ M(γ·,·
+ δ̃·,· ) − M(γ·,·
)
h
i
X
≥
En (p̂t ({x∗i 0 γt∗ }NT ) − dti )x∗i 0 δ̃t ,
t∈NT

applying the convexity of M. Using the bound in Eqn. (E.7) and rearranging we find that


√
λD
∗
∗
0 ≤ λD
γ·,·
−
γ
+
δ̃
+
δ̃·,·
+ bds TEn [k{x∗i 0 δ̃t }NT k22 ]1/2 .
·,·
·,·
2,1
2
2,1
2,1
Dividing through λD and decomposing the supports, we find that
√


1
bds T
∗
∗
0≤
δ̃·,·
+
γ·,· 2,1 − γ·,· + δ̃·,·
En [k{x∗i 0 δ̃t }NT k22 ]1/2
+
2
λD
2,1
2,1
1
δ̃·,S∗
=
2

1
+
δ̃·,S∗c
2
2,1

∗
γ·,S
∗

+
2,1

2,1

−

∗
γ·,S
∗

+ δ̃·,S∗

−
2,1

δ̃

·,S∗c

√
bds T
+
En [k{x∗i 0 δ̃t }NT k22 ]1/2 ,
λD
2,1

∗
where the second line follows because γ·,S
c = 0. Collecting terms and applying the triangle inequality yields
∗

1
δ̃·,S∗c
2

2,1

1
δ̃·,S∗
≤
2
1
≤
δ̃·,S∗
2
1
δ̃·,S∗
2
1
=
δ̃·,S∗
2

+
2,1

+
2,1

≤

+
2,1

+
2,1

√
bds T
En [k{x∗i 0 δ̃t }NT k22 ]1/2
−
+ δ̃·,S∗
+
2,1
λD
2,1
√
bds T
∗
∗
γ·,S∗ 2,1 − γ·,S∗ + δ̃·,S∗
En [k{x∗i 0 δ̃t }NT k22 ]1/2
+
λD
2,1
√


bds T
∗
∗
En [k{x∗i 0 δ̃t }NT k22 ]1/2
γ·,S∗ − γ·,S∗ + δ̃·,S∗
+
λD
2,1
√
bds T
δ̃·,S∗
+
En [k{x∗i 0 δ̃t }NT k22 ]1/2 .
λD
2,1
∗
γ·,S
∗

∗
γ·,S
∗

Therefore with probability at least 1 − P
δ̃·,S∗c

2,1

≤ 3 δ̃·,S∗

+
2,1

√
2bds T
En [k{x∗i 0 δ̃t }NT k22 ]1/2
λD

(E.8)

Consider two cases based on the upper bound in (E.8). First, suppose that δ̃·,· obeys the cone constraint

of Eqn. (14) in the definition of κ2D , such that
δ̃·,S∗c

2,1

≤ 4 δ̃·,S∗

.
2,1

This implies
δ̃·,·

=
2,1

δ̃·,S∗

δ̃·,S∗c

+
2,1

2,1

≤ 5 δ̃·,S∗
2,1
p
≤ 5 |S∗ | δ̃·,S∗
2
p
5 |S∗ |
≤
En [k{x∗i 0 δ̃t }NT k22 ]1/2 ,
κD

by the Cauchy-Schwarz inequality, the restricted eigenvalue definition of Eqn. (14), and noting that
En [k{x∗i 0 δ̃t }NT k22 ]. Collecting across the second and third inequalities yields
p
|S∗ |
δ̃·,S∗
≤
En [k{x∗i 0 δ̃t }NT k22 ]1/2 .
κD
2,1

(E.9)
P

t∈NT

δ̃t0 Qδ̃t =

(E.10)

On the other hand, if the cone constraint fails, then
δ̃·,S∗

<
2,1

1
δ̃·,S∗c
4

.
2,1

Using this for the first and third inequalities, and Eqn. (E.8) for the second, we have
δ̃·,·

=
2,1

δ̃·,S∗

+
2,1

δ̃·,S∗c

5
δ̃·,S∗c
4
15
≤
δ̃·,S∗
4
15
δ̃·,S∗c
≤
16
≤

2,1

2,1

√
5 bds T
+
En [k{x∗i 0 δ̃t }NT k22 ]1/2
2 λD
2,1
√
5 bds T
+
En [k{x∗i 0 δ̃t }NT k22 ]1/2
2 λD
2,1

Combining the right hand side of the first line with third lines yields
√
bds T
En [k{x∗i 0 δ̃t }NT k22 ]1/2 .
δ̃·,S∗c
≤8
λD
2,1
Plugging this back into the last line we obtain the bound
√
bds T
δ̃·,·
≤ 10
En [k{x∗i 0 δ̃t }NT k22 ]1/2 ,
λD
2,1
while instead, plugging it into the failure of the cone constraint yields
√
bds T
δ̃·,S∗
≤2
En [k{x∗i 0 δ̃t }NT k22 ]1/2 .
λD
2,1

(E.11)

(E.12)

Combining Equations (E.9) and (E.11) gives the first claim of the lemma and Equations (E.10) and (E.12)
give the second.

E.2

Proof of Theorem 5

∗
Define δ̃·,· = γ̃·,· − γ·,·
. By the optimality of δ̃·,· , we have
∗
∗
+ δ̃·,·
M(γ·,·
+ δ̃·,· ) + λD γ·,·

∗
∗
≤ M(γ·,·
) + λD γ·,·

2,1

2,1

.

Rearranging and subtracting the score, we have
X


∗
∗
M(γ·,·
+ δ̃·,· ) − M(γ·,·
)−
En (p̂t ({x∗i 0 γt∗ }NT ) − dti )x∗i 0 δ̃t
t∈NT


≤ λD

∗
γ·,·

−
2,1

∗
γ·,·
+ δ̃·,·


2,1



En (p̂t ({x∗i 0 γt∗ }NT ) − dti )x∗i 0 δ̃t .

X

−

(E.13)

t∈NT

The proof proceeds by deriving a further upper bound to the right and a quadratic lower bound of the left.
The combination of these will yield a bound on En [(x∗i 0 δ̃t )2 ]1/2 .
Let us begin with the right side of Eqn. (E.13). For the penalized difference of coefficients we have
∗
γ·,S
c
∗

∗
γ·,S
c + δ̃·,S c
∗
∗

−
2,1

=
2,1

δ̃·,S∗c

,
2,1

∗
because γ·,S
c = 0. Therefore,
∗
∗
γ·,·

2,1

−

∗
γ·,·
+ δ̃·,·

=
2,1

≤
≤
≤

∗
γ·,S
∗

2,1

∗
γ·,S
∗

2,1

∗
γ·,S
∗

−
−

2,1

−


∗
∗
γ·,S
− γ·,S
∗
∗

∗
γ·,S
+ δ̃·,S∗
∗
∗
γ·,S
+ δ̃·,S∗
∗

−
2,1

δ̃·,S∗c

2,1

∗
γ·,S
+ δ̃·,S∗
∗
2,1

+ δ̃·,S∗
= δ̃·,S∗
2,1

2,1

,
2,1

where the first inequality reflects dropping the nonpositive final term (the norm is nonnegative) and the
third inequality follows from the triangle inequality. Using this result for the first term and the bound (E.7)
for the second, the right side of Eqn. (E.13) is bounded by
λD δ̃·,S∗
≤
≤

√
λD
δ̃·,·
+ bds TEn [k{x∗i 0 δ̃t }NT k22 ]1/2
+
2
2,1
2,1
( p
!
(p
√ )
√ )
d
√
|S∗ | 2bs T
λD 5 |S∗ | 10bds T
d
λD
∨
+
∨
+ bs T En [k{x∗i 0 δ̃t }NT k22 ]1/2
κD
λD
2
κD
λD
!
p
√
λD |S∗ |
d
6
(E.14)
+ 8bs T En [k{x∗i 0 δ̃t }NT k22 ]1/2 ,
κD

where the second inequality applies the results of Lemma E.3 and the third bounds the maximum by the
sum.
Now turn to the left side of Eqn. (E.13). Our goal is to show that this is bounded below by a quadratic
function. We apply the bounds for Bach’s (2010) modified self-concordant functions. To show that M(·)
belongs to this class, we must bound the third derivative in terms of the Hessian. Recall that p̂t ({x∗i 0 γt }NT ) =

exp{x∗i 0 γt }/ 1 +

P
NT


exp{x∗i 0 γt } and the T-square matrix H({x∗i 0 γt }NT ) has (t, t0 ) ∈ N2T entry given by
(

H({x∗i 0 γt }NT )[t,t0 ]

=

p̂t ({x∗i 0 γt }NT )(1 − p̂t ({x∗i 0 γt }NT ))
−p̂t ({x∗i 0 γt }NT )p̂t0 ({x∗i 0 γt }NT )

if t = t0
if t =
6 t0

First, note that M(γ·,· ) can be written as

 X
 
X
dti (x∗i 0 γt ) .
M(γ·,· ) = En log 1 +
exp{x∗i 0 γt } −
t∈NT

t∈NT




P
P
Define F : RT → R as F (w) = log 1 + t∈NT exp(wt ) , so that M(γ·,· ) = En F (wi ) − t∈NT dti wi,t , where
wi,t = x∗i 0 γt and wi = {wi,t }NT . Then for any w ∈ RT , v ∈ RT , and scalar α, define g(α) = F (w + αv) :
R → R. We verify the conditions of Bach (2010, Lemma 1) for this g(α) and F (w). This involves finding the
third derivative of g(α), and bounding it in terms of the second (i.e. the Hessian). To this end, note that the
multinomial function has the property that ∂ p̂t ({x∗i 0 γt }NT )/∂γt = p̂t ({x∗i 0 γt }NT )(1 − p̂t ({x∗i 0 γt }NT ))x∗i and
∂ p̂t ({x∗i 0 γt }NT )/∂γt0 ,· = −p̂t ({x∗i 0 γt }NT )p̂t0 ({x∗i 0 γt }NT )x∗i . From these, we find that
X
g 0 (α) = v 0 F 0 (w + αv) =
vt p̂t (w + αv)
t∈NT

and
g 00 (α) = v 0 F 00 (w + αv)v = v 0 H(w + αv)v.
To bound g 000 (α), we again use the derivatives of p̂t ({x∗i 0 γt }NT ) to find the derivatives of elements H(w).
Routine calculations give, for any r 6= s 6= t:
∂H(w)t,t /∂wt = p̂t (w)(1 − p̂t (w))(1 − 2p̂t (w)) = H(w)t,t (1 − 2p̂t (w))
∂H(w)t,t /∂wr = −p̂t (w)p̂r (w)(1 − p̂t (w)) + p̂t (w)2 p̂r (w) = H(w)t,t (p̂t (w)p̂r (w)(1 − p̂t (w))−1 − p̂r (w))
∂H(w)t,s /∂wt = −p̂t (w)p̂s (w)(1 − 2p̂t (w)) = H(w)t,s (1 − 2p̂t (w))
∂H(w)t,s /∂wr = −p̂t (w)p̂s (w)(−2p̂r (w)) = H(w)t,s (−2p̂r (w)).
Each derivative returns the same Hessian element multiplied by term bounded by 2 in absolute value. Let
ar represent this factor. Then we bound
g 000 (α) =

X
r∈NT

X
∂v 0 H(w̃)v
=
vr v 0 H(w + αv)var
∂wr
w̃=w+αv
r∈NT
X
X
√
≤
v 0 H(w + αv)v|vr ||ar | ≤ 2v 0 H(w + αv)v
|vr | = 2kvk1 g 00 (α) ≤ 2 Tkvk2 g 00 (α).

vr

r∈NT

r∈NT

Applying Bach’s (2010) Lemma 1 to each observation, as in Belloni, Chernozhukov, and Wei (2013)27 , with
wi = {x∗i 0 γt∗ }NT and vi = {x∗i 0 δ̃t }NT we get the lower bound
X


∗
∗
M (γ·,·
+ δ̃·,· ) − M(γ·,·
)−
En (p̂t ({x∗i 0 γt∗ }NT ) − dti )x∗i 0 δ̃t
t∈NT


≥ En
27


vi0 H({x∗i 0 γt }NT )vi  −2kvi k2
e
+
2kv
k
−
1
i
2
4Tkvi k22

Kwemou (2012) also applied Bach’s (2010) to study sparse logistic regression


≥ En

vi0 H({x∗i 0 γt }NT )vi
4Tkvi k22



4
2kvi k22 − kvi k32 ,
3

(E.15)

where the second inequality follows from Belloni, Chernozhukov, and Wei (2013, Lemma 9).
Tanabe and Sagae (1992, Theorem 1) give H({x∗i 0 γt∗ }NT ) ≥ φmin {H({x∗i 0 γt∗ }NT )}IT , in the positive
definite sense, where φmin (A) denotes the smallest eigenvalue of A and IT is the T × T identity matrix. Then
φmin {H({x∗i γt∗ }NT )} ≥ det{H({x∗i 0 γt }NT )} =

Y

p̂t ({x∗i 0 γt∗ }NT ) ≥

t∈NT

where p0 ({x∗i 0 γt∗ }NT ) = 1 −
(1992). These results imply

En



pmin
Ap

T
,

∗0 ∗
t∈NT p̂t ({xi γt }NT ) and the first inequality is also due to Tanabe and Sagae
that vi0 H({x∗i 0 γt }NT )vi ≥ (pmin /Ap )T vi0 IT vi = (pmin /Ap )T kvi k22 and therefore

P

vi0 H({x∗i 0 γt }NT )vi
4Tkvi k22



4
2kvi k22 − kvi k32
3





T



1
4
En 2kvi k22 − kvi k32
4T
3

T


pmin
1 En [kvi k22 ]
2 En [kvi k32 ]
=
1−
.
Ap
T
2
3 En [kvi k22 ]
≥

pmin
Ap

(E.16)

Recall that vi = {x∗i 0 δ̃t }NT . To prove a quadratic lower bound, consider two cases, depending on whether
!
2 En [k{x∗i 0 δ̃t }NT k32 ]
1
1−
2
3 En [k{x∗i 0 δ̃t }NT k22 ]
is above or below 1/AK .
In the first case, combining Equations (E.15) and (E.16) gives
∗
∗
M(γ·,·
+ δ̃·,· ) − M(γ·,·
)−

X



En (p̂t ({x∗i 0 γt∗ }NT ) − dti )x∗i 0 δ̃t ≥

t∈NT



pmin
Ap

T

1 En [k{x∗i 0 δ̃t }NT k22 ]
.
T
AK

(E.17)

Now consider the second case, where this bound does not hold. By Assumption 2(b), the Cauchy-Schwarz
inequality, and the conclusion of Lemma E.3
X X
√
k{x∗i 0 δ̃t }NT k1 =
x∗i,j δ̃t,j ≤ X δ̃·,· ≤ TX δ̃·,·
2,1

1

t∈NT j∈Np

( p
√ )
√
5 |S∗ | 10bds T
≤ TX
∨
En [k{x∗i 0 δ̃t }NT k22 ]1/2 .
κD
λD

Hence, by subadditivity (to bound the `2 norm by the `1 norm),
En [k{x∗i 0 δ̃t }NT k32 ]

≤

En [k{x∗i 0 δ̃t }NT k22 k{x∗i 0 δ̃t }NT k1 ]

≤

En [k{x∗i 0 δ̃t }NT k22 ]3/2

√

( p
√ )
5 |S∗ | 10bds T
TX
∨
.
κD
λD

Thus
1
1
>
AK
2

2 En [k{x∗i 0 δ̃t }NT k32 ]
1−
3 En [k{x∗i 0 δ̃t }NT k22 ]

!

1
≥
2

!
√
√ 
p
2 X T 
d
∗0
2 1/2
1−
5λD |S∗ | + 10κD bs T En [k{xi δ̃t }NT k2 ]
,
3 κD λD

which is equivalent to

√ −1
p
2
3 κD λD 
√
5λD |S∗ | + 10κD bds T
:= rn .
AK 2 X T


P
∗
Because M(γ·,·
+ δ·,· ) − M(γ·,· ) − t∈NT En (p̂t ({x∗i 0 γt∗ }NT ) − dti )x∗i 0 δt is convex in δ·,· , and hence any line
segment lies above the function, we know that En [k{x∗i 0 δ̃t }NT k22 ]1/2 > rn , so we have
En [k{x∗i 0 δ̃t }NT k22 ]1/2 >

∗
M(γ·,·
+δ̃·,· )−M(γ·,· )−

X
t∈NT



1−



En [k{x∗i 0 δ̃t }NT k22 ]1/2
= rn En [k{x∗i 0 δ̃t }NT k22 ]1/2 .
En (p̂t ({x∗i 0 γt∗ }NT ) − dti )x∗i 0 δ̃t ≥ rn2 ≥ rn2
rn

Combining this result with Equations (E.13) and (E.14), we have


2
1−
AK



√ −1
p
3 κD λD 
√
5λD |S∗ | + 10κD bds T
En [k{x∗i 0 δt }NT k22 ]1/2 ≤
2 X T

6

λD

p

|S∗ |

κD

+

8bds

√

!
T En [k{x∗i 0 δ̃t }NT k22 ]1/2 ,

which is impossible under the restriction on AK because it is equivalent to
√

√
p
2 X T 
2
d
2
d 2
2
|S
|κ
b
T
+
80κ
(b
)
T
30λ
|S
|
+
100λ
≤
1−
∗ D s
D
D s
D ∗
AK
3 κ2D λD
√ 

√
p
2X T
=
30λD |S∗ | + 100 |S∗ |κD bds T + 80κ2D (bds )2 Tλ−1
.
D
2
3 κD
Solving this for AK would require that
AK ≤ 2

κ2D

√ 
√
p
κ2D − (2/3)X T 30λD |S∗ | + 100 |S∗ |κD bds T + 80κ2D (bds )2 Tλ−1
D

which contradicts the condition in the Theorem.
Therefore, Eqn. (E.17) must hold.28 Combining this with Equations (E.13) and (E.14), we find that


pmin
Ap

T

1 En [k{x∗i 0 δ̃t }NT k22 ]
≤
T
AK

6

λD

p

|S∗ |

κD

+

8bds

√

!
T En [k{x∗i 0 δ̃t }NT k22 ]1/2 .

Thus, dividing through and applying the union bound we find that
max En [(x∗i 0 δ̃t )2 ]1/2
t∈NT

≤

En [k{x∗i 0 δ̃t }NT k22 ]1/2


≤

Ap
pmin

T

TAK

6

λD

p

|S∗ |

κD

√
d

!

+ 8bs T .

(E.18)

To bound the propensity score error, we apply the mean value theorem and the form of ∂ p̂t ({x∗i 0 γt }NT )/∂γt .
We must linearize with respect to t only (recall that p̂t ({x∗i 0 γ̃t }NT ) depends on all of γ̃·,· ). To this end, define
Mt as the T-vector with entry t given by x∗i 0 γt∗ + m̃t x∗i 0 γ̃t for a scalar m̃t ∈ [0, 1] and entries t0 ∈ NT \ {t}
equal to x∗i 0 γt0 . Then we have
p̂t ({x∗i 0 γ̃t }NT ) − p̂t ({x∗i 0 γt∗ }NT ) = p̂t (Mt )[1 − p̂t (Mt )]x∗i 0 δ̃t ≤ x∗i 0 δ̃t .
28

(E.19)

This analysis is conceptually similar to using Belloni and Chernozhukov’s (2011) restricted nonlinearity impact
coefficient, but our characterization is different.

Using this result coupled with the triangle inequality, the bias condition, and Eqn. (E.18), we find
En [(p̂t ({x∗i 0 γ̃t }NT ) − pt (xi ))2 ]1/2 ≤ En [(p̂t ({x∗i 0 γ̃t }NT ) − p̂t ({x∗i 0 γt∗ }NT ))2 ]1/2 + En [(p̂t ({x∗i 0 γt∗ }NT ) − pt (xi ))2 ]1/2
i1/2
h
≤ En (x∗i 0 δ̃t )2
+ bds
!
p

T
√
λD |S∗ |
Ap
≤
TAK 6
+ 8bds T + bds .
pmin
κD
The `1 bound follows from Eqn. (E.18) by the Cauchy-Schwarz inequality and the definition in Eqn.
(16):
kγ̃t −

γt∗ k1

q
∗ | kγ̃ − γ ∗ k
≤ |S̃ D ∪ SD
t
t 2,p ≤

∗
|S̃ D ∪ SD
|
∗}
D
φ{Q, S̃ ∪ SD

!1/2
En [(x∗i 0 (γ̃t − γt∗ ))2 ]1/2 .

Finally, we bound the size of the selected set of coefficients. First, note that optimality of γ̃·,· ensures
that |S̃ D | ≤ n. Then, restating the conclusion Lemma E.2 using the notation of the Theorem and the rate
result (E.18), then bounding φ by φ we find that
∗
|S̃ D | ≤ |SD
|4Ln φ{Q, |S̃ D |}.

The argument now parallels that used by Belloni and Chernozhukov (2013), relying on their result on the
sublinearity of sparse eigenvalues. Let dme be the ceiling function and note that dme ≤ 2m. For any m ∈ ND
Q,
suppose that |S̃ D | > m. Then,
∗
|S̃ D | ≤ |SD
|4Ln φ{Q, m(|S̃ D |/m)}
l
m
∗
≤ |S̃ D |/m |SD
|4Ln φ{Q, m}
∗
≤ (|S̃ D |/m)|SD
|8Ln φ{Q, m}.

Rearranging gives
∗
m ≤ |SD
|8Ln φ{Q, m}
D
whence m 6∈ ND
Q . Minimizing over NQ gives the result.

E.3

Proof of Theorem 6

∗
Define δ̂·,· = γ̂·,· − γ·,·
. Many of the arguments parallel those for Theorem 5. The key differences are that


P
∗
∗
a quadratic lower bound for M(γ·,·
+ δ̂·,· ) − M(γ·,·
) − t∈NT En (p̂t ({x∗i 0 γt∗ }NT ) − dti )x∗i 0 δ̂t may occur, but
is not necessary, and δ̂·,· may not belong to the cone of the restricted eigenvalues, but obeys the sparse
eigenvalue constraints.


P
∗
∗
We first give a suitable upper bound for M(γ·,·
+ δ̂·,· ) − M(γ·,·
) − t∈NT En (p̂t ({x∗i 0 γt∗ }NT ) − dti )x∗i 0 δ̂t .
By the Cauchy-Schwarz inequality and the definition of the sparse eigenvalues of Eqn. (16),

δ̂·,·

X

=
2,1

δ̂·,j

∗
j∈ŜD ∪SD

r
≤

∗
ŜD ∪ SD

2

sX

X

t∈NT j∈ŜD ∪S ∗
D

2
δ̂t,j

s

r
ŜD ∪

=

∗
SD

X

t∈NT

s

r
∗
ŜD ∪ SD

≤

X

2

δ̂·,j

2

n
o−2
∗
φ Q, ŜD ∪ SD
δ̂t0 Qδ̂t

t∈NT

r
∗
ŜD ∪ SD

=

n
o−1
∗
φ Q, ŜD ∪ SD
En [k{x∗i 0 δ̂t }NT k22 ]1/2 .

(E.20)

Following identical steps to Equations (E.5), (E.6), and (E.7), but with δ̂·,· in place of δ̃·,· , and then using
the above bound, we have
√


λD
+ bds TEn [k{x∗i 0 δ̂t }NT k22 ]1/2
δ̂·,·
En (p̂t ({x∗i 0 γt∗ }NT ) − dti )x∗i 0 δ̂t ≤
2
2,1
t∈NT
q


∗|
|ŜD ∪ SD
√
λ
D
≤
+ bds T  En [k{x∗i 0 δ̂t }NT k22 ]1/2 .
∗}
2 φ{Q, ŜD ∪ SD
X

(E.21)

∗
∗
). By optimality of the post selection estimator M(γ̂·,· ) ≤ M(γ̃·,· ),
+ δ̂·,· ) − M(γ·,·
Next we turn to M(γ·,·
∗
∗
D
∗
). By the mean value
) ≤ M(γ̃·,· ) − M(γ·,·
+ δ̂·,· ) − M(γ·,·
as S̃ ⊂ ŜD by construction, and hence M(γ·,·
theorem, for scalars {mt ∈ [0, 1]}NT we have
h
i
X
∗
∗
M(γ·,·
+ δ̃·,· ) − M(γ·,·
)=
En (dti − p̂t ({x∗i 0 γt∗ + mt x∗i 0 δ̃t }))x∗i 0 δ̃t
t∈NT

=

X

h
i
En (dti − p̂t ({x∗i 0 γt∗ }NT ))x∗i 0 δ̃t

t∈NT

+

X

h
i
En (p̂t ({x∗i 0 γt∗ }NT ) − p̂t ({x∗i 0 γt∗ + mt x∗i 0 δ̃t }))x∗i 0 δ̃t ,

t∈NT

h
i
X
√
λD
δ̃·,·
+ bds TEn [k{x∗i 0 δ̃t }NT k22 ]1/2 +
En mt (x∗i 0 δ̃t )2 .
2
2,1
t∈NT
q


∗|
|ŜD ∪ SD
√
λD
+ bds T  En [k{x∗i 0 δ̃t }NT k22 ]1/2 + En [k{x∗i 0 δ̃t }NT k22 ], (E.22)
≤
∗
2 φ{Q, ŜD ∪ SD }
≤

where the first inequality follows from Eqn. (E.7) and the same steps as in (E.19) while the second applies
(E.20) with δ̃·,· and mt ≤ 1.29
Collecting the bounds of (E.21) and (E.22), and the definition of RM (that is, Eqn. (E.18)) gives
∗
∗
M(γ·,·
+ δ̂·,· ) − M(γ·,·
)−

X



En (p̂t ({x∗i 0 γt∗ }NT ) − dti )x∗i 0 δ̂t

t∈NT


≤

q
∗|
|ŜD ∪ SD

λD
∗}
2 φ{Q, ŜD ∪ SD




√
2
+ bds T  En [k{x∗i 0 δ̂t }NT k22 ]1/2 + RM + RM
. (E.23)

Next, we turn to a lower bound. Consider the same two cases as in the proof of Theorem 5. In the first
29

Applying the steps of Eqn. (E.20) to δ̃·,· is preferred to using the results of Lemma E.3 because it leads to the
∗
tidier expression involving φ{Q, ŜD ∪ SD
}, but the latter method could be substituted.

case, we have the quadratic lower bound:
∗
∗
M (γ·,·
+ δ̂·,· ) − M(γ·,·
)−



En (p̂t ({x∗i 0 γt∗ }NT ) − dti )x∗i 0 δ̂t ≥

X
t∈NT



pmin
Ap

T

1 En [k{x∗i 0 δ̂t }NT k22 ]
.
T
AK

(E.24)

In the other case, this bound may not hold. Arguing as in the proof of Theorem 5, but applying Eqn. (E.20),
we get
√ q
∗0
2 1/2
∗ |φ{Q, Ŝ ∪ S ∗ }−1 E [k{x∗ 0 δ̃ }
.
k{xi δ̂t }NT k1 ≤ TX |ŜD ∪ SD
D
n
D
i t NT k2 ]
Therefore, as above, we find
∗
M(γ·,·
+ δ̂·,· ) − M(γ·,· ) −

X



En (p̂t ({x∗i 0 γt∗ }NT ) − dti )x∗i 0 δ̂t ≥ rn En [k{x∗i 0 δ̂t }NT k22 ]1/2 ,

(E.25)

t∈NT

with
rn =

3
2



∗
}
φ{Q, ŜD ∪ SD
2
q
1−
.
√
AK
∗|
X T |ŜD ∪ SD

Collecting the upper bound of (E.23) and the lower bounds (E.24) and (E.25) we have
(

pmin
Ap

T

)
n
o
1 En [k{x∗i 0 δ̂t }NT k22 ]
∧ rn En [k{x∗i 0 δ̂t }NT k22 ]1/2
T
AK
q


∗|


|ŜD ∪ SD
√
λ
D
d
2
 En [k{x∗i 0 δ̂t }NT k22 ]1/2 + RM + RM
+
b
T
≤
. (E.26)
s
∗}
2 φ{Q, ŜD ∪ SD

For some A1 > 1, replace the restriction on AK in the Theorem with the requirement that

AK > 2






∗ 2
}
φ{Q, ŜD ∪ SD







q
√
√ 


∗
∗
∗
∗
2
d

 φ{Q, ŜD ∪ SD } − (A1 /3)X T λD |ŜD ∪ SD | + bs φ{Q, ŜD ∪ SD } |ŜD ∪ SD | T 




∗
φ{Q, ŜD ∪ SD }
q
∨
.
√

∗ } − (A /3)2R X T |Ŝ ∪ S ∗ | 
φ{Q, ŜD ∪ SD
1
M
D
D


Suppose the linear term is the minimum. The first restriction on AK implies, by simple algebraic manipulations, that
q


∗|
|ŜD ∪ SD
√
λ
rn
D
d

+ bs T  <
,
∗}
2 φ{Q, ŜD ∪ SD
A1
while the second gives RM < (rn /A1 ). Plugging the former into Eqn. (E.26) and then applying the latter
yields


2
rn En [k{x∗i 0 δ̂t }NT k22 ]1/2 ≤ (rn /A1 ) En [k{x∗i 0 δ̂t }NT k22 ]1/2 + RM + RM


≤ (rn /A1 ) En [k{x∗i 0 δ̂t }NT k22 ]1/2 + 2RM .

Canceling the rn and solving yields
En [k{x∗i 0 δ̂t }NT k22 ]1/2 ≤

2RM
.
A1 − 1

On the other hand, if the quadratic term is the minimum, define
q


T

∗|
|
Ŝ
∪
S
√
D
λD
Ap
D
0
TAK 
+ bds T  .
RM
=
∗}
pmin
2 φ{Q, ŜD ∪ SD
With this notation and the quadratic term being the minimum, Eqn. (E.26) becomes
En [k{x∗i 0 δ̂t }NT k22 ]

≤

0
RM
En [k{x∗i 0 δ̂t }NT k22 ]1 /2

Then, because a2 ≤ ab + c implies that a ≤ b +
En [k{x∗i 0 δ̂t }NT k22 ]1/2

√

0
≤ RM
+

+

0
RM
RM


+

Ap
pmin

T

2
TAK RM
.

c, we have
0
RM
RM +



Ap
pmin

T

!1/2
2
TAK RM

.

Taking A1 = 3 and combining the bounds on En [k{x∗i 0 δ̂t }NT k22 ]1/2 from the two cases gives

!1/2 
T



A
p
0
2
0
0
En [k{x∗i δ̂t }NT k22 ]1/2 ≤ {RM } ∨ RM
TAK RM
.
+ RM
RM +


pmin
From this bound on the log-odds estimates, we obtain the bound on the propensity score estimates and
the `1 rate, given by,

!1/2 

T


Ap
0
0
2
max En [(p̂t ({x∗i 0 γ̂t }NT ) − pt (xi ))2 ]1/2 ≤ {RM } ∨ RM
+ RM
RM +
TAK RM
+ bds ,
t∈NT


pmin
and
max kγ̂t − γt∗ k1 ≤

t∈NT

∗
|S̃ D ∪ SD
|
∗}
D
φ{Q, S̃ ∪ SD

!1/2
{RM } ∨




0
RM
+



0
RM
RM +



Ap
pmin

T

2
TAK RM

!1/2 


,



by arguments parallel to those used in the proof of Theorem 5.

Appendix F

Proofs for Group Lasso Selection and Estimation of
Linear Models

Unless otherwise noted, all bounds in this section are nonasymptotic. We will use generic notation X ∗ , δ, s,
etc, as this section deals only with linear models.

F.1

Lemmas

Lemma F.1 (Score Bound). For λY and P defined respectively in Eqn. (11) and Eqn. (12) we have


λY
P max kEn,· [ui x∗i,j ]k2 ≥
≤ P.
j∈Np
4
Proof. The residuals ui are conditionally mean-zero by definition. Using this, Assumption 2(a), the definitions of X and U, and the Cauchy-Schwarz inequality, we find that
X 1
X 
X 1



X2 U2 T
∗ 4 1/2
E En,t [ui x∗i,j ]2 =
E[u2i (x∗i,j )2 ] ≤
E[|Xi,j
E kEn,· [ui x∗i,j ]k22 =
| ] E[|Ui |4 ]1/2 ≤
nt
nt
n
t∈NT

t∈NT

t∈NT

uniformly in j ∈ Np . Define the mean-zero random variables ξt,j as:
ξt,j = (En,t [ui x∗i,j ])2 −
Let rn = T

1
E[U 2 Xj∗ 2 ].
nt

−1/2

log(p ∨ n)3/2+δ . Then, using the definition of λY :




λ2
λY
= P max kEn,· [ui x∗i,j ]k22 ≥ Y
P max kEn,· [ui x∗i,j ]k2 ≥
j∈Np
j∈Np
4
16


X2 U2 T
X2 U2 Trn
= P max kEn,· [ui x∗i,j ]k22 ≥
+
j∈Np
n
n


2 2
2 2
X U T
X U Trn
∗
2
= P max kEn,· [ui xi,j ]k2 −
≥
j∈Np
n
n


X
X2 U2 Trn 
ξt,j ≥
≤ P max
j∈Np
n
t∈NT


X
n
≤ E max
ξt,j 
,
j∈Np
X2 U2 Trn

(F.1)

t∈NT

where final line follows from Markov’s inequality.
Next, applying Lemma 9.1 of Lounici, Pontil, van de Geer, and Tsybakov (2011) (with their m = 1 and
hence c(m) = 2) followed by Jensen’s inequality and Assumption 2(c), we find that



1/2 
X
X


2 
E max
ξt,j  ≤ (8 log(2p))1/2 E 
max ξt,j

j∈Np

j∈Np

t∈NT

t∈NT

1/2

 
≤ (8 log(2p))1/2 E 

X

2 
max ξt,j
j∈Np

t∈NT


 1/2
X X4 U4
X 
4
≤ 4 log(2p)1/2 
+
E max En,t [ui x∗i,j ]  .
j∈Np
n2


t∈NT

(F.2)

t∈NT

√ √
√
2
The leading 4 is 8 2, where 2 is a byproduct of applying the inequality (a−b)2 ≤ 2(a2 +b2 ) to ξt,j
. Again
using Lemma 9.1 of Lounici, Pontil, van de Geer, and Tsybakov (2011) (with their m = 4, and c(m) = 12

since c(4) ≥ (e4−1 − 1)/2 + 2 ≈ 11.54), we bound the expectation in the second term above as follows:

!4/2 


∗ 2
X
u
x
i
4
i,j

E max En,t [ui x∗i,j ]
≤ [8 log(12p)]4/2 E 
max
j∈Np
j∈Np
nt
i∈It


64 log(12p)2 X4 X  4  X X
≤
E |ui | +
E[|ui |2 |uk |2 ]
n4t
i∈It k∈It \{i}
i∈It


64 log(12p)2 X4 X  4  X X
≤
E |ui | +
E[|ui |4 ]1/2 E[|uk |4 ]1/2 
n4t
i∈It k∈It \{i}

i∈It

64 log(12p) X U
(nt + nt (nt − 1))
n4t
64 log(12p)2 X4 U4
,
=
n2t
2

4

4

≤

(F.3)

where the second inequality uses Hölder’s inequality and the final inequality applies Assumptions 2(a) and
2(c).
Now, inserting the results of Eqns. (F.2) and (F.3) into Eqn. (F.1), we have

1/2


X 64 log(12p)2 X4 U4
λY
4n log(2p)1/2  X X4 U4
∗

P max kEn,· [ui xi,j ]k2 ≥
+
≤
j∈Np
4
n2
n2t
TX2 U2 rn
t∈NT

≤
using the choice rn = T

−1/2

4 log(2p)
√
rn T

1/2

t∈NT

[1 + 64 log(12p)2 ]1/2 = P,

log(p ∨ n)3/2+δ .

Lemma F.2 (Estimate Sparsity). With probability 1 − P, as defined in Eqn. (12), the model selected by
solving (9) obeys
16 X
|S̃ Y | ≤ 2
φ{Qt , S̃ Y }En,t [(µt (xi ) − x∗i 0 β̃t )2 ].
λY
t∈NT

Proof. First, by Karush-Kuhn-Tucker conditions for (9), for all t ∈ NT , if β̃·,j 6= 0 it must satisfy
2En,t [x∗i,j (yi − x∗i 0 β̃t )] = λY

β̃t,j
.
kβ̃·,j k2

(F.4)

Hence, taking the `2 -norm over t ∈ NT for fixed j ∈ S̃ Y , using yi = µt (xi ) + ui , the triangle inequality, and
Lemma F.1:
λY = 2 En,t [x∗i,j (yi − x∗i 0 β̃t )]
≤2

En,t [x∗i,j {µt (xi )

−

2

x∗i 0 β̃t }]

≤ 2 En,t [x∗i,j {µt (xi ) − x∗i 0 β̃t }]

2
2

+ 2 En,t [x∗i,j ui ]

2

+ λY /2.

Let Gt be the vector of {µt (xi )}i∈It and G̃t that of {x∗i 0 β̃t }i∈It . Collecting terms, then squaring both sides

and taking

P

j∈S̃ Y

(i.e. applying k · k22 over j ∈ S̃ Y to both sides) yields
X

X X

λ2Y ≤ 16

2
En,t [x∗i,j {µt (xi ) − x∗i 0 β̃t }]

j∈S̃ Y t∈NT

j∈S̃ Y

i
X 1 h
0
X
(G
−
G̃
)
t
t
t
n2t
j∈S̃ Y

= 16

t∈NT

2
2

X φ{Qt , S̃ Y }
2
Gt − G̃t
nt
2
t∈NT
X
≤ 16
φ{Qt , S̃ Y }En,t [(µt (xi ) − x∗i 0 β̃t )2 ].
≤ 16

t∈NT

The claim follows, as the left-hand side is equal to |S̃ Y |λ2Y .

F.2

Proof of Theorem 7

∗
∗
. First, because β·,S
Let δ̃·,· = β̃·,· − β·,·
c = 0
∗
∗
β·,S
c
∗

∗
β·,S
c + δ̃·,S c
∗
∗

−
2,1

=

δ̃·,S∗c

2,1

.
2,1

Therefore:
∗
β·,·

2,1

−

∗
β·,·
+ δ̃·,·

=
2,1

≤
≤
≤

∗
β·,S
∗
∗
β·,S
∗
∗
β·,S
∗

2,1

2,1

−
−

2,1

−


∗
∗
β·,S
− β·,S
∗
∗

∗
β·,S
+ δ̃·,S∗
∗
∗
β·,S
∗

−
2,1

+ δ̃·,S∗

δ̃·,S∗c

2,1

2,1

∗
β·,S
+ δ̃·,S∗
∗
2,1

+ δ̃·,S∗
= δ̃·,S∗
2,1

,

(F.5)

2,1

where the first inequality reflects dropping the nonpositive final term (the norm is nonnegative) and the
third inequality follows from the triangle inequality. Because β̃·,· solves (9)
E(δ̃·,· ) + λY

β̃·,·

2,1

≤ E(δ̃·,· ) + λY

∗
β·,·

2,1

.

Define the ith realization of BtY as byt,i = µt (xi ) − x∗i 0 βt∗ . Inserting yi = x∗i 0 βt∗ + byt,i + ui on each side, we
obtain
X

En,t [(x∗i 0 δ̃t )2 ] − 2

t∈NT

X

En,t [(byt,i + ui )x∗i 0 δ̃t ] +

t∈NT

X

En,t [(byt,i + ui )2 ] + λY

∗
β·,·
+ δ̃·,·

2,1

t∈NT

≤

X

En,t [(byt,i + ui )2 ] + λY

∗
β·,·

2,1

t∈NT

Canceling common factors and rearranging, we find that
X
t∈NT

En,t [(x∗i 0 δ̃t )2 ] ≤ 2

X
t∈NT

En,t [(byt,i + ui )x∗i 0 δ̃t ] + λY



∗
β·,·

−
2,1

∗
β·,·
+ δ̃·,·


(F.6)
2,1

By the Cauchy-Schwarz inequality and Lemma F.1, with probability at least 1 − P
X
X X
2
En,t [ui x∗i 0 δ̃t ] = 2
En,t [ui x∗i,j δ̃t ]
j∈Np t∈NT

t∈NT

sX
X sX
∗
2
2
2
≤
(En,t [ui xi,j ])
δ̃t,j
j∈Np

X

=

t∈NT

t∈NT

2kEn,· [ui x∗i,j ]k2 kδ̃·,j k2

j∈Np

λY
2

≤

.

δ̃·,·

2,1

Next, using the Cauchy-Schwarz inequality, the bias condition, and Jensen’s inequality:
X
X
2
En,t [byt,i x∗i 0 δ̃t ] ≤ 2
En,t [(byt,i )2 ]1/2 En,t [(x∗i 0 δ̃t )2 ]1/2
t∈NT

t∈NT

En,t [(x∗i 0 δ̃t )2 ]1/2

≤ 2bys

X

≤ 2bys

sX

t∈NT

En,t [(x∗i 0 δ̃t )2 ].

t∈NT

Plugging the previous two inequalities into Eqn. (F.6) we find that with probability at least 1 − P:


sX
X
λY
∗
∗
y
∗0
2
δ̃·,·
+ λY
β·,· 2,1 − β·,· + δ̃·,·
+ 2bs
En,t [(xi δ̃t ) ] ≤
En,t [(x∗i 0 δ̃t )2 ]. (F.7)
2
2,1
2,1
t∈NT

t∈NT

Consider two cases, depending on whether
X

En,t [(x∗i 0 δ̃t )2 ] − 2bys

sX

En,t [(x∗i 0 δ̃t )2 ]

t∈NT

t∈NT

is negative or nonnegative. In the first case, rearranging the display above gives
sX
X
En,t [(x∗i 0 δ̃t )2 ] < 2bys
En,t [(x∗i 0 δ̃t )2 ]

(F.8)

t∈NT

t∈NT

For the second case, returning to Eqn. (F.7), rearranging, and discarding positive terms (under the
second case) from the left side, we have


λY
∗
∗
0≤
δ̃·,·
+ λY
β·,· 2,1 − β·,· + δ̃·,·
.
2
2,1
2,1
Canceling λY , decomposing the supports
0≤

1
δ̃·,S∗
2

+
2,1

1
δ̃·,S∗c
2

+
2,1

∗
β·,S
∗

2,1

−

∗
β·,S
+ δ̃·,S∗
∗

−
2,1

δ̃·,S∗c

,
2,1

Collecting terms and applying the final inequality of Eqn. (F.5) yields
1
δ̃·,S∗c
2

≤
2,1

1
δ̃·,S∗
2

+
2,1

∗
β·,S
∗

2,1

−

∗
+ δ̃·,S∗
β·,S
∗

≤
2,1

3
δ̃·,S∗
2

,
2,1

and hence δ̃·,· obeys the cone constraint of Eqn. (15).
Thus, beginning with Eqn. (F.7), decomposing the support of δ̃·,· , using the cone constraint and the
result of Eqn. (F.5), the Cauchy-Schwarz inequality, the definition of κY from Eqn. (15),


sX
X
λY
λY
y
∗
∗
+
2b
β
+
δ̃
+
−
+ λY
β·,·
δ̃·,S∗
δ̃·,S∗c
En,t [(x∗i 0 δ̃t )2 ] ≤
En,t [(x∗i 0 δ̃t )2 ]
·,·
s
·,·
2,1
2
2
2,1
2,1
2,1
t∈NT
t∈NT
sX
λY
λY
+ 2bys
+
+ λY δ̃·,S∗
δ̃·,S∗
≤
En,t [(x∗i 0 δ̃t )2 ]
3 δ̃·,S∗
2
2
2,1
2,1
2,1
t∈NT
sX
≤ 3λY δ̃·,S∗
+ 2bys
En,t [(x∗i 0 δ̃t )2 ]
2,1

t∈NT

≤ 3λY

p

3λY

p

|S∗ | δ̃·,S∗

2

+ 2bys

sX

En,t [(x∗i 0 δ̃t )2 ]

t∈NT

≤

=

κY

|S∗ |

sX

δ̃t0 Qt δ̃t + 2bys

sX

t∈NT

En,t [(x∗i 0 δ̃t )2 ]

t∈NT

!s
p
X
3λY |S∗ |
y
+ 2bs
En,t [(x∗i 0 δ̃t )2 ].
κY

(F.9)

t∈NT

Equations (F.8) and (F.9) show that in both cases defined above, the root left side appears on the right.
Thus, dividing through in both we find that
X

En,t [(x∗i 0 δ̃t )2 ]

3λY

≤

t∈NT

!2
p
|S∗ |
y
+ 2bs
,
κY

because the bound given in Eqn. (F.9) contains that of (F.8). From the union bound we have
!
p
3λY |S∗ |
∗0
2 1/2
y
≤
+ 2bs ,
max En,t [(xi δ̃t ) ]
κY
t∈NT
and therefore, by the triangle inequality
max En,t [(x∗i 0 β̃t
t∈NT

2 1/2

− µt (xi )) ]

≤

3λY

p
κY

|S∗ |

!
+

2bys

+ bys .

(F.10)

The rate above pertains only to the “with-in sample” fit, for those observations with dti = 1. To obtain
a rate on the entire sample, we use the sparse eigenvalues defined in Eqn. (16), as follows:
−1/2

En [(x∗i 0 δ̃t )2 ]1/2 = Q1/2 Qt
≤

1/2

Qt δ̃t
2
!1/2
Y
∗
φ{Q, S̃ ∪ SY }
1/2
Qt δ̃t
φ{Qt , S̃ Y ∪ SY∗ }

2

!1/2

=

φ{Q, S̃ Y ∪ SY∗ }
φ{Qt , S̃ Y ∪ SY∗ }

!1/2

≤

φ{Q, S̃ Y ∪ SY∗ }
φ{Qt , S̃ Y ∪ SY∗ }

En,t [(x∗i 0 δ̃t )2 ]1/2
3λY

p

|S∗ |

κY

!
+

2bys

.

The first conclusion of the Theorem now follows from this rate, the triangle inequality, and the bias
condition, because
En [(x∗i 0 β̃t − µt (xi ))2 ]1/2 ≤ En [(x∗i 0 β̃t − x∗i 0 βt∗ )2 ]1/2 + En [(byt,i )2 ]1/2
!1/2
!
p
3λY |S∗ |
φ{Q, S̃ Y ∪ SY∗ }
y
≤
+ 2bs + bys .
κY
φ{Qt , S̃ Y ∪ SY∗ }
The `1 bound now follows by the Cauchy-Schwarz inequality and the definition in Eqn. (16):
δ̃t

1

q
≤ |S̃ Y ∪ SY∗ | δ̃t

≤
2

|S̃ Y ∪ SY∗ |
φ{Q, S̃ Y ∪ SY∗ }

!1/2
En [(x∗i 0 δ̃t )2 ]1/2 .

Finally, we bound the size of the selected set of coefficients. First, note that optimality of β̃·,· ensures
that |S̃ Y | ≤ n. Then, restating the conclusion Lemma F.2 using the notation of the Theorem and Eqn.
(F.10), then bounding φ by φ we find that
X
X
|S̃ Y | ≤ |SY∗ |16Ln
φ{Qt , S̃ Y } ≤ |SY∗ |16Ln
φ{Qt , |S̃ Y |}.
t∈NT

t∈NT

The argument now parallels that used by Belloni and Chernozhukov (2013), relying on their result on the
sublinearity of sparse eigenvalues. Let dme be the ceiling function and note that dme ≤ 2m. For any m ∈ NYQ ,
suppose that |S̃ Y | > m. Then,
X
|S̃ Y | ≤ |SY∗ |16Ln
φ{Qt , m(|S̃ Y |/m)}
t∈NT

l
m
X
≤ |S̃ Y |/m |SY∗ |16Ln
φ{Qt , m}
t∈NT

≤ (|S̃

Y

|/m)|SY∗

|32Ln

X

φ{Qt , m}.

t∈NT

Rearranging gives
X

m ≤ |SY∗ |32Ln

φ{Qt , m}

t∈NT

whence m 6∈ NYQ . Minimizing over NYQ gives the result.

F.3

Verification of Assumption 3(c) for Group Lasso Estimators

Under conditions imposed therein, Section 6 of the paper verifies that Assumptions 3(a) and 3(b) hold for
the proposed group lasso estimators µ̂t (xi ) and p̂t (xi ). Here we show that 3(c) holds also. No additional
assumptions are required.
Recall that rather than generic µ̂t (x) and µt (xi ), as above, we are now explicitly considering high-

dimensional approximately sparse linear models for µt (xi ). In this context, we add and subtract the pseudotrue values to write
√
nEn [(µ̂t (xi ) − µt (xi ))(1 − dti /pt (xi ))] = A1 + A2 ,
where
n

1 X ∗ ∗
A1 = √
(x β − µt (xi ))
n i=1 i t



pt (xi ) − dti
pt (xi )

n


and

1 X ∗0
A2 = √
(x β̂t − x∗i βt∗ )
n i=1 i



pt (xi ) − dti
pt (xi )


.

For the first term, E[A1 |{xi }ni=1 ] = 0 holds as βt∗ is nonrandom. From Assumption 1(b) and the definition
of the bias term bys we find
"
2 #
n
pt (xi ) − dti
1X ∗ ∗
2
2
n
(x β − µt (xi )) E
≤ CEn (x∗i βt∗ − µt (xi ))2 ≤ C(bys )2 .
E[A1 |{xi }i=1 ] =
n i=1 i t
pt (xi )
Therefore |A1 | = OPn (bys ) = oPn (1), where the second equality is assumed in the bias condition of Assumption
4 and the first equality follows from Markov’s
inequality.


For the second term, define Σ̃t,j = E (x∗i,j )2 (dti − pt (xi ))2 /pt (xi )2 and then proceed as follows:

n 
1 X pt (xi ) − dti X ∗
∗
xi,j (β̂t,j − βt,j
)
A1 = √
pt (xi )
n i=1
j∈ŜY
)
(
n
X
1 X x∗i,j (pt (xi ) − dti )/pt (xi )
1/2
∗
√
Σ̃t,j (β̂t,j − βt,j
)
=
1/2
n i=1
Σ̃
t,j
j∈ŜY
!


n
1 X x∗i,j (pt (xi ) − dti )/pt (xi )
1/2
max √
≤ max Σ̃t,j
β̂t − βt∗
1/2
j∈Np
j∈Np
n i=1
Σ̃
t,j

= O(1)OPn (log(p)) β̂t − βt∗

1

.
1

This quantity is oPn (1) by Corollary 5 in the original paper, which among other results, gives a rate for
the `1 norm of the estimated coefficients. For the final equality, Assumptions 1(b), 2(b), and 2(c) imply
that maxj∈Np Σ̃t,j = O(1), while the center factor is bounded by applying the moderate deviation theory
for self-normalized sums of de la Peña, Lai, and Shao (2009, Theorem 7.4) and in particular Belloni, Chen,
Chernozhukov, and Hansen (2012, Lemma 5). To apply this theory, first note that the summand of the
center factor has bounded third moment and second moment bounded away from zero from Assumptions
1(b) and 2. Σt,j normalizes the second moment, and the theory applies under Assumption 4.

Appendix G

Additional Simulation Results

The DGP is as described in the main text. The mean comparison group sample sizes for the various DGPs
are in Table G.1. Figure G.1 shows the analogue of Figure 1 Panels (a) and (b) with 2000 covariates. The
manually chosen values of δD and δY do not appear well-suited to one particular DGP. Figures G.2, G.3,
and G.4 show the coverage results based using 10-fold cross validation to choose the penalty parameters, for
1000, 1500, and 2000 covariates, respectively. For these three, the exponents αβ and αγ range from one to
four, and hence the functions are always sparse (to a certain degree). For nonsparse functions, the current
R routines are not reliable. This will be explored in future software development. Cross-validation choices
perform very well.

Table G.1: Mean Comparison Group Sample Sizes for Various
Specifications
Multiplier ργ (exponent αγ = 2)
No. of Covariates:
1000
1500
2000

0.01
498.368
497.832
498.368

0.05
488.311
487.987
488.311

0.25
439.159
438.608
439.159

0.50
384.060
383.966
384.060

0.75
338.106
338.408
338.106

1
301.793
301.924
301.793

Exponent αγ (multiplier ργ = 1)
No. of Covariates:
1000
1500
2000

0.125

0.25

0.5

.75

1

2

3

4

456.373
462.646
462.646

420.073
426.711
426.711

341.358
342.527
342.527

312.818
312.544
312.544

305.861
305.965
305.965

301.793
301.924
301.924

302.134
301.943
301.943

302.764
302.229
302.229

Figure G.1: Empirical Coverage of 95% Confidence Intervals, Varying Signal Strength and
Sparsity of pt (x) and µt (x), 2000 Covariates

µ(x) signal

µ(x) sparsity

pt(x) signal

1.0

pt(x) sparsity

1.0

0.5

0.5

0.0

0.0

(a) Varying signal strength

(b) Varying sparsity

Figure G.2: Empirical Coverage of 95% Confidence Intervals, Penalty Chosen with
Cross-Validation, Varying Signal Strength and Sparsity of pt (x) and µt (x), 1000 Covariates

µ(x) signal

µ(x) sparsity

pt(x) signal

1.0

pt(x) sparsity

1.0

0.5

0.5

0.0

0.0

(a) Varying signal strength

(b) Varying sparsity

Figure G.3: Empirical Coverage of 95% Confidence Intervals, Penalty Chosen with
Cross-Validation, Varying Signal Strength and Sparsity of pt (x) and µt (x), 1500 Covariates

µ(x) signal

µ(x) sparsity

pt(x) signal

1.0

pt(x) sparsity

1.0

0.5

0.5

0.0

0.0

(a) Varying signal strength

(b) Varying sparsity

Figure G.4: Empirical Coverage of 95% Confidence Intervals, Penalty Chosen with
Cross-Validation, Varying Signal Strength and Sparsity of pt (x) and µt (x), 2000 Covariates

µ(x) signal

µ(x) sparsity

pt(x) signal

1.0

pt(x) sparsity

1.0

0.5

0.5

0.0

0.0

(a) Varying signal strength

(b) Varying sparsity

