arXiv:1412.4605v3 [] 27 Jan 2017

Valid confidence intervals for post-model-selection predictors
François Bachoc*, Hannes Leeb**, and Benedikt M. Pötscher**
*Department of Mathematics, University Paul Sabatier
**Department of Statistics, University of Vienna
First version: December 2014
This version: December 2016

Abstract
We consider inference post-model-selection in linear regression. In this setting, Berk
et al. (2013a) recently introduced a class of confidence sets, the so-called PoSI intervals,
that cover a certain non-standard quantity of interest with a user-specified minimal coverage probability, irrespective of the model selection procedure that is being used. In this
paper, we generalize the PoSI intervals to confidence intervals for post-model-selection
predictors.
AMS Mathematics Subject Classification 2010: 62F25, 62J05.
Keywords: Inference post-model-selection, confidence intervals, optimal post-model-selection
predictors, non-standard targets, linear regression.

1

Introduction and overview

In statistical practice, the model used for analysis is very often chosen after the data have
been observed, either by ad-hoc methods or by more sophisticated model selection procedures.
Inference following such a model selection step (inference post-model-selection) has proven to
be a challenging problem. ‘Naive’ procedures, which ignore the presence of model selection, are
typically invalid (e.g., in the sense that the actual coverage probability of ‘naive’ confidence
sets for the true parameter can be dramatically smaller than the nominal one), and the
construction of valid procedures is often non-trivial; see Leeb and Pötscher (2005, 2006, 2017),
Kabaila and Leeb (2006), Pötscher (2009) and references therein for an introduction to the
issues involved here. In these references, inference is focused on the true parameter of the datagenerating model (or on components thereof). Shifting the focus away from the true parameter
as the target of inference, Berk et al. (2013a) recently introduced a class of confidence sets,
the so-called PoSI intervals, that guarantee a user-specified minimal coverage probability after
model selection in linear regression, irrespective of the model selector that is being used; see
also Berk et al. (2013b) and Leeb et al. (2015). In this paper, we generalize the PoSI intervals
to intervals for post-model-selection predictors.

1

Prediction following model selection is obviously also of great importance. In the case
where the selected model is misspecified, parameter estimates are typically biased or at least
difficult to interpret; cf. Remark 2.7. But even a misspecified model may perform well for
prediction. In particular, Greenshtein and Ritov (2004) derive, under appropriate sparsity
assumptions, feasible predictors that asymptotically perform as well as the (infeasible) best
candidate predictor even if the available number of explanatory variables by far exceeds the
sample size. These feasible predictors are also covered by the results in the present paper,
among others. Like Greenshtein and Ritov (2004), our analysis does not rely on the assumption that the true data generating model is among the candidates for model selection. We
develop confidence intervals for such predictors, that are easy to interpret and that are optimal in an appropriate sense; cf. Remarks 2.7(ii), 2.8, and 3.2, as well as Greenshtein and
Ritov (2004). A further rationale for extending the PoSI-approach of Berk et al. (2013a) to
problems related to prediction is that this framework seems to provide a more natural habitat
for considering non-standard targets; see the discussion in Remark 2.1 of Leeb et al. (2015)
as well as in Remarks 2.7 and 3.1 given further below.
The crucial feature of the approach of Berk et al. (2013a) is that the coverage target, i.e.,
the quantity for which a confidence set is desired, is not the standard target, i.e., the parameter
in an overall model (or components thereof), but a non-standard quantity of interest that
depends on the selected model and thus on the data. This non-standard quantity of interest
(n)
is denoted by β
throughout the paper (cf. Section 2 for details). Here M̂ stands for the
M̂
(data-dependent) model chosen by the model selector and n stands for sample size. The non(n)
standard target β
provides a certain vector of regression coefficients for those explanatory
M̂

(n)

variables that are ‘active’ in the model M̂ (more precisely, β
represents the coefficients of
M̂
the projection of the expected value vector of the dependent variable on the space spanned
by the regressors included in M̂ ); for a precise definition see eqs. (3) and (4) in Section 2.
For a new set of explanatory variables x0 , we first extend the PoSI-approach to obtain
(n)
confidence intervals for the predictor x00 [M̂ ]β . Here, x0 [M̂ ] denotes the set of explanatory
M̂

(n)

variables from x0 that correspond to the ‘active’ regressors in the model M̂ . We call x00 [M̂ ]β
M̂
the design-dependent (non-standard) coverage target, because different design matrices in
(n)
the training data typically result in different values of x00 [M̂ ]β
even if both training data
M̂

sets lead to selection of the same model M̂ . We construct PoSI confidence intervals for
(n)
x00 [M̂ ]β
that guarantee a user-specified minimal coverage probability, irrespective of the
M̂
model selector that is being used. The design-dependent coverage target minimizes a certain
‘in-sample’ prediction error; cf. Remark 2.8. However, when the goal is to predict a new
response corresponding to a new vector x0 of explanatory variables, this ‘in-sample’ optimality
property may have little relevance and thus the focus on covering the design-dependent target
(n)
x00 [M̂ ]β
may be debatable.
M̂
In view of this, we next consider an alternative coverage target that depends on the selected
(?)
model but not on the training data otherwise, and that we denote by x00 [M̂ ]β . We call
M̂

(?)

x00 [M̂ ]β
the design-independent (non-standard) coverage target. The design-independent
M̂
coverage target minimizes a certain ‘out-of-sample’ prediction error, namely the mean-squared

2

prediction error, over all (infeasible) predictors of a future response y0 that are of the form
x00 [M̂ ]γ(M̂ ), when x0 and the row-vectors of X are sampled from the same distribution; cf.
Remark 3.2. In particular, this target does not suffer from the issues that plague the designdependent coverage target, as discussed at the end of the preceding paragraph. Certain
(?)
optimality properties of a feasible counterpart of x00 [M̂ ]β are derived in Greenshtein and
M̂

Ritov (2004), for a particular model selector M̂ and under appropriate sparsity assumptions;
(?)
a target closely related to x00 [M̂ ]β is also studied in Leeb (2009). For a large class of model
M̂
selectors, we show that the PoSI confidence intervals constructed earlier also cover the designindependent coverage target with minimal coverage probability not below the user-specified
nominal level asymptotically. In that sense, the PoSI confidence intervals are approximately
(?)
valid for the target x00 [M̂ ]β , irrespective of the model selector M̂ in that class. In simulations
M̂
we find that our asymptotic result is representative of the finite-sample situation even for
moderate sample sizes.
When extending the PoSI-approach to confidence intervals for both the design-dependent
(n)
(?)
and the design-independent coverage target, i.e., for both x00 [M̂ ]β
and x00 [M̂ ]β , we find
M̂

M̂

that the resulting intervals necessarily depend not only on x0 [M̂ ] but also on those components
of x0 that are ‘in-active’ in the model M̂ . This may appear surprising at first sight but turns
out to be inherent to the PoSI-approach (because of the need to take the maximum over all
models M in (10)). In any case, this is problematic in situations when, after having selected
a given model, only the ‘active’ components of x0 are observed, e.g., in situations where
observations are costly and model selection is carried out also with the goal of reducing cost
by not having to observe irrelevant components of x0 . To resolve this, we also develop PoSI
confidence intervals that depend on the ‘active’ variables x0 [M̂ ] only. These intervals are
obtained by maximizing over all inactive variables and are hence larger than the intervals for
the case where x0 is known entirely. In simulations, we find that the excess width of these
intervals is moderate. We also provide analytic results regarding the excess width of these
intervals in an asymptotic setting where the number of regressors goes to infinity, see Section
2.4.
Inference post-model-selection is currently a very active area of research and we can only
give a selection of work relevant for, or related to, this paper. Contemporary analyses of confidence sets for (components of) the true parameter of the underlying model include Andrews
and Guggenberger (2009), Kabaila and Leeb (2006), Leeb and Pötscher (2005), Pötscher
(2009), Pötscher and Schneider (2010), and Schneider (2016). These references also point to
numerous earlier results. Also, the work of Lockhart et al. (2014), Wasserman and Roeder
(2009), and Wasserman (2014) should be mentioned here. For the LASSO, in particular, a
de-sparsifying method has recently been developed by Belloni et al. (2011, 2014), van de Geer
et al. (2014), and Zhang and Zhang (2014). Another strand of literature that, like the PoSI
(n)
approach, also focuses on β as the quantity of interest, is developed in Fithian et al. (2015),
M̂
Lee et al. (2016), Lee and Taylor (2014), Tian and Taylor (2015), Tibshirani et al. (2015), and
(n)
are considered that have a
Tibshirani et al. (2016): In these papers, confidence sets for β
M̂
guaranteed coverage probability conditionally on the event that a particular model has been
selected by the model selection procedure. In contrast to PoSI procedures, the confidence

3

intervals obtained in these papers are specific to the model selection procedure used (the
LASSO, in particular, being considered in these references) and generally rely on certain geometric properties of the specific model selection procedure under consideration. In simulation
experiments, we compare the confidence intervals proposed in these references with the intervals developed here and observe some interesting phenomena, see Section 4.3. As prompted
by a referee, we point out here that in the presence of a large number of regressors PoSI
intervals (including intervals considered in the present paper) typically are computationally
more burdensome than the confidence intervals proposed in Lee et al. (2016) for the LASSO
with a fixed value for the tuning parameter; see, however, also the discussion towards the end
of Section 4.3.
The rest of the paper is organized as follows. In Section 2, we introduce the models, the
model-selection procedures, the design-dependent target, and the PoSI confidence intervals for
both the case where all explanatory variables in x0 are observed and the case where only the
components of x0 corresponding to the ‘active’ explanatory variables are available; moreover,
we analyze properties of these intervals in an asymptotic framework where the model dimension increases; cf. Section 2.4. In Section 3, we present the design-independent target and
show that the PoSI confidence intervals introduced earlier also cover the design-independent
target, with minimal coverage probability not below the nominal one asymptotically when
sample size increases. The results of a numerical study are reported in Section 4. Conclusions
are drawn in Section 5. Appendix A contains some comments on the assumptions made on
the errors variance. The proofs of the results in Sections 2 and 3 are given in Appendices B
and C. Appendix D contains some comments on and extensions of the results in Section 3.
In Appendix E we describe algorithms for computing the PoSI confidence intervals, that are
comparable with those proposed by Berk et al. (2013a) in terms of computational complexity.
Finally, Appendix F contains details concerning the numerical calculations used for the results
in Section 4.

2
2.1

Confidence intervals for the design-dependent non-standard
target
The framework

Consider the model
Y =µ+U

(1)

where µ ∈ Rn is unknown and U follows an N 0, σ 2 In -distribution; here σ 2 , 0 < σ < ∞, is
the unknown error variance and In is the identity matrix of size n ≥ 1. An important instance
of this model arises when µ is known to reside in a lower dimensional linear subspace of Rn ,
but we do not make such an assumption at this point. Apart from the data Y , we are given
a (real) n × p matrix X, not necessarily of full column rank, the columns of which represent
potential regressors. This setup allows for p > n as well as for 1 ≤ p ≤ n. The rank of X will
be denoted by d. The design matrix X is treated as fixed throughout Section 2.


4

We consider fitting (potentially misspecified) linear models with design matrices that are
obtained by deleting columns from X. Such a model will be represented by M , a subset of
{1, ..., p}, where the elements of M index the columns of X that are retained. We use the
following notation: For M ⊆ {1, ..., p}, we write M c for the complement of M in {1, ..., p}. It
proves useful to allow M to be the empty set. We write |M | for the cardinality of M . With
m = |M |, let us write M = {j1 , ..., jm } in case m ≥ 1. For M 6= ∅ and for an l × p matrix T ,
l ≥ 1, let T [M ] be the matrix of dimension l×m obtained from T by retaining only the columns
of T with indices j ∈ M and deleting all others; if M = ∅ we set T [M ] = 0 ∈ Rl . In abuse of
notation we shall, for a p × 1 vector v, write v[M ] for (v 0 [M ])0 , i.e., v[M ] = (vj1 , ..., vjm )0 for
m ≥ 1 and v[M ] = 0 ∈ R in case M = ∅. For a given model M , we denote the corresponding
least squares estimator by β̂ M , i.e.,
−1
X[M ]0 Y,
β̂ M = X[M ]0 X[M ]

(2)

where the inverse is to be interpreted as the Moore-Penrose inverse in case X[M ] does not
have full column rank. For any given model M the corresponding least squares estimator β̂ M
is obviously an unbiased estimator of
(n)

β M = X[M ]0 X[M ]
(n)

−1

X[M ]0 µ.

(3)

Note that β̂ M as well as β M reduce to 0 in case M = ∅.
As in Berk et al. (2013a) we further assume that, as an estimator for σ 2 , we have available
an (observable) random variable σ̂ 2 that is independent of PX Y and that is distributed as
σ 2 /r times a chi-square distributed random variable with r degrees of freedom (1 ≤ r < ∞),
with PX denoting orthogonal projection on the column space of X. This assumption is always
satisfied in the important special case where one assumes that d < n and µ ∈ span(X) hold,
upon choosing for σ̂ 2 the standard residual variance estimator obtained from regressing Y on
X and upon setting r = n−d. However, otherwise it is not an innocuous assumption at all and
this is further discussed in Appendix A. Observe that our assumption allows for estimators
σ̂ 2 that not only depend on Y and X, but possibly also on other observable random variables
(e.g., additional data). The joint distribution of Y and σ̂ 2 depends on µ and σ as well as on
sample size n and will be denoted by Pn,µ,σ (see also Appendix D.4).
We are furthermore given a (non-empty) collection M of admissible models M ⊆ {1, ..., p},
the ‘universe’ of models considered by the researcher. Without loss of generality we will
assume that
S any column of X appears as a regressor in at least one of the models M in M,
i.e., that {M : M ∈ M} = {1, ..., p} holds (otherwise we can just redefine X by discarding
all columns that do not appear in any of the models in M); of course, we have excluded here
the trivial and uninteresting case M = {∅}. For such a collection M it is easy to see that the
assumed independence of σ̂ 2 and PX Y is in fact equivalent to independence of σ̂ 2 from the
collection {β̂ M : M ∈ M} of least squares estimators. While not really affecting the results,
it proves useful to assume, throughout the following, that the empty model belongs to M.
We shall furthermore always assume that any non-empty M ∈ M is of full-rank in the sense
that rank X[M ] = |M |. We point out here that our assumptions on M imply that X can not
have a zero column, and hence d ≥ 1 must hold. An important instance of a collection M
5

satisfying our assumptions is the collection of all full-rank submodels of {1, ..., p} (enlarged
by the empty model) provided that no column of X is zero; of course, there are many other
examples, see, e.g., the list in Section 4.5 of Berk et al. (2013a).
A model selection procedure M̂ is now a (measurable) rule that associates with every
(X, Y, σ̂ 2 ) a (possibly empty) model M̂ (X, Y, σ̂ 2 ) ∈ M. In the following we shall, in abuse
of notation, often write M̂ for M̂ (X, Y, σ̂ 2 ). Allowing explicitly dependence of M̂ on σ̂ 2 is
only relevant in case σ̂ 2 depends on extraneous data beyond (X, Y ) and the model selection
procedure actually makes use of σ̂ 2 . [We note here that in principle we could have allowed
M̂ to depend on further extraneous data, in which case Pn,µ,σ would have to be redefined as
the joint distribution of Y , σ̂ 2 , and this further extraneous data.] The post-model-selection
estimator β̂ M̂ corresponding to the model selection procedure is now given by (2) with M
replaced by M̂ .
The non-standard quantity of interest studied in Berk et al. (2013a) is the random vector
(n)
(with random dimension) β
obtained by replacing M by M̂ in (3). The situation we shall
M̂
consider in the present paper is related to Berk et al. (2013a), but is different in several aspects:
Consider a fixed (real) p × 1 vector x0 and suppose we want to predict y0 which is distributed
as N ν, σ 2 , independently of Y . If one is forced to use a fixed model M for prediction,
i.e., to use predictors of the form x00 [M ]γ, the predictor that would then typically be used
(n)
is x00 [M ]β̂ M , which can be viewed as an estimator of the infeasible predictor x00 [M ]β M . Of
course, for this predictor to be reasonable there must be some relation between the training
data (X, Y ) and (x0 , y0 ). This is further discussed in Remark 2.8. In the presence of model
selection the predictor x00 [M ]β̂ M will then typically be replaced by the post-model-selection
predictor x00 [M̂ ]β̂ M̂ which can in turn be seen as a feasible counterpart to the infeasible
predictor
(n)
x00 [M̂ ]β .
(4)
M̂

The quantity in (4) will be our target for inference throughout Section 2 and will be called the
design-dependent (non-standard) target (to emphasize that it depends on the design matrix
X apart from its dependence on M̂ , cf. (3)). A discussion of the merits of this target and its
interpretation is postponed to Remarks 2.7 and 2.8 given below.
Let now 1 − α ∈ (0, 1) be a nominal confidence level. Throughout Section 2 we are
(n)
interested in confidence intervals for the design-dependent target x00 [M̂ ]β
that are of the
M̂
form
CI(x0 ) = x00 [M̂ ]β̂ M̂ ± K(x0 , M̂ )||sM̂ ||σ̂,
(5)
where k·k denotes the Euclidean norm (σ̂ of course representing the nonnegative square root
of σ̂ 2 ), where
−1
s0M = x00 [M ] X[M ]0 X[M ]
X[M ]0 ,
(6)
where sM = 0 ∈ Rn for M = ∅ by our conventions, and where K(x0 , M ) = K(x0 , M, r) =
K(x0 , M, r, X, α, M) denotes a non-negative constant which may depend on x0 , M , r, X, α,
and M, but does not depend on the observations on Y and σ̂ 2 . Here we have used the notation
a±b for the interval [a − b, a + b] (a ∈ R, b ≥ 0). The motivation for the form of the confidence
interval stems from the observation that for fixed M the interval x00 [M ]β̂ M ± qr,1−α/2 ||sM ||σ̂
6

(n)

is a valid 1 − α confidence interval for x00 [M ]β M , where qr,1−α/2 is the (1 − α/2)-quantile of
Student’s t-distribution with r degrees of freedom. Furthermore note that on the event M̂ =∅
the target is equal to zero and the confidence interval reduces to {0}, thus always containing
the target on this event. Finally note that CI(x0 ) constitutes a confidence interval for the
(n)
predictor x00 [M̂ ]β , and should not be mistaken for a prediction interval for a new response
M̂
y0 .
We aim at finding quantities K(x0 , M ) such that the confidence intervals CI(x0 ) satisfy


(n)
0
inf
P
x
[
M̂
]β
∈
CI(x
)
≥ 1 − α.
(7)
n,µ,σ
0
0
n
M̂

µ∈R ,σ>0

Note that if one replaces K(x0 , M̂ ) in (5) by Knaive = qr,1−α/2 , then the confidence interval
(5) reduces to the so-called ‘naive’ confidence interval which is constructed as if M̂ were fixed
a priori (thus ignoring the presence of model selection). It does not fulfill (7) as can be seen
from the numerical results in Section 4, which is in line with the related results in Leeb et al.
(2015).

2.2

The various confidence intervals

For the construction of the quantities K(x0 , M ) we distinguish two cases regarding the observation on x0 : (i) The vector x0 is observed in its entirety (regardless of which model M̂ is
selected), or (ii) only the subvector x0 [M̂ ] of x0 is observed (note that only this subvector is
needed for the computation of the post-model-selection predictor x00 [M̂ ]β̂ M̂ ). The former case
will arise if measuring all the components of x0 is not too costly, whereas the latter case will
be relevant in practical situations where the selected model is determined first and then only
observations for x0 [M̂ ] (and not for the other components of x0 ) are collected, e.g., out of
cost considerations. For example, in a medical application one may want to avoid measuring
prognostic variables that require invasive procedures or that incur high monetary costs, see,
e.g., Castera et al. (2015). Cost considerations in the context of model selection or prediction
are also common in fields such as industrial process control or engineering (Jaupi (2014),
Souders and Stenbakken (1991)).
For the case (i), where x0 is entirely observed, the following straightforward adaptation of
the approach in Berk et al. (2013a) yields a constant K1 (x0 ) = K1 (x0 , r) = K1 (x0 , r, X, α, M)
(not depending on M ) such that the resulting confidence interval (5) satisfies (7): Observe
that
(n)
x00 [M̂ ]β̂ M̂ − x00 [M̂ ]β = s0M̂ (Y − µ) ,
(8)
M̂

define s̄M = sM / ksM k if sM 6= 0, and set s̄M = 0 ∈ Rn if sM = 0. Then obviously we have
the upper bound
(9)
s̄0M̂ (Y − µ) /σ̂ ≤ max s̄0M (Y − µ) /σ̂.
M ∈M

Define K1 (x0 ) to be the smallest constant satisfying


Pn,µ,σ max s̄0M (Y − µ) /σ̂ ≤ K1 (x0 ) ≥ 1 − α.
M ∈M

7

(10)

It is important to note that the probability on the left-hand side of the preceding display
neither depends on µ nor on σ; it also depends on the estimator σ̂ only through the ‘degrees
of freedom’ parameter r: To see this note that s̄0M (Y − µ) = s̄0M PX (Y − µ), since s̄M belongs
to the column space of X. Consequently, the collection of all the quantities s̄0M (Y − µ) is
jointly distributed as N (0, σ 2 C), independently of σ̂ 2 ∼ σ 2 /r χ2 (r), where the covariance
matrix C depends only on x0 and X. Hence the joint distribution of the collection of ratios
|s̄0M (Y − µ)| /σ̂ does neither depend on µ nor σ, and depends on the estimator σ̂ only through
r. It is now plain that K1 (x0 ) only depends on x0 , r, X, α, and M. Furthermore note that
K1 (x0 ) = 0 in case x0 = 0; otherwise, K1 (x0 ) is positive, equality holds in (10), and K1 (x0 )
is the unique (1 − α)-quantile of the distribution of the upper bound in (9). [This follows
from Lemma B.1 in Appendix B and from the observation that, in view of our assumptions
on M, s̄0M = 0 for all M ∈ M holds if and only if x0 = 0.] Furthermore, observe that K1 (x0 )
coincides with a PoSI1 constant of Berk et al. (2013a) in case x0 is one of the standard basis
vectors ei . [This can be seen by comparison with (4.14) in Berk et al. (2013a) and noting that
the maximum inside the probability in (10) effectively extends only over models satisfying
i ∈ M , since s̄M = 0 holds for models M with i ∈
/ M if x0 = ei .] Finally, Knaive ≤ K1 (x0 )
clearly holds provided x0 6= 0 (since s̄0M (Y − µ) /σ̂ follows Student’s t-distribution with r
degrees of freedom if sM 6= 0).
As a consequence of (9) and the discussion in the preceding paragraph we thus immediately
obtain the following proposition.
Proposition 2.1. Let M̂ be an arbitrary model selection procedure with values in M, let
x0 ∈ Rp be arbitrary, and let K1 (x0 ) be defined by (10). Then the confidence interval (5) with
K(x0 , M̂ ) replaced by K1 (x0 ) satisfies the coverage property (7).
The coverage in Proposition 2.1 is guaranteed for all model selection procedures with
values in M, and thus leads to ‘universally valid post-selection inference’ in case M is chosen
to be the set of all full-rank submodels obtainable from X (enlarged by the empty set and
provided X does not have a zero column); cf. Berk et al. (2013a), where similar guarantees
(n)
are obtained for the components of β . [In fact, the construction of K1 (x0 ) implies that
M̂

the collection of intervals x00 [M ]β̂ M ± K1 (x0 )||sM ||σ̂ with M ∈ M provides a simultaneous
(n)
confidence band for x00 [M ]β M .]
Consider next case (ii) where only the components of x0 [M̂ ] are observed. In this case, the
confidence interval of Proposition 2.1 is not feasible in that it cannot be computed in general,
because K1 (x0 ) will depend on all components of x0 (and not only on those appearing in
x0 [M̂ ]) due the maximum figuring in (10) and our assumptions on M. A first solution is to
define
K2 (x0 [M ], M ) = sup {K1 (x) : x[M ] =x0 [M ]} ,
(11)
and then to use the confidence interval (5) with K(x0 , M̂ ) replaced by K2 (x0 [M̂ ], M̂ ). Note
that K2 (x0 [M ], M ), and hence the corresponding confidence interval, depends on x0 only via
x0 [M ], and thus can be computed in case (ii). Of course, K2 (x0 [M ], M ) also depends on r,
X, α, and M, and we shall write K2 (x0 [M ], M, r) if we want to stress dependence on r. It
is easy to see that K2 (x0 [M ], M ) is finite (as it is not larger than the Scheffé constant as we
8

shall see below). Because K2 (x0 [M ], M ) is never smaller than K1 (x0 ), we have the following
corollary to Proposition 2.1.
Corollary 2.2. Let M̂ be an arbitrary model selection procedure with values in M, let x0 ∈ Rp
be arbitrary, and let K2 (x0 [M ], M ) be defined by (11). Then the confidence interval (5) with
K(x0 , M̂ ) replaced by K2 (x0 [M̂ ], M̂ ) satisfies the coverage property (7).
The computation of K2 (x0 [M̂ ], M̂ ) is more costly than that of K1 (x0 ). Indeed, it requires to embed the algorithm for computing K1 (x0 ) in an optimization procedure. Thus,
for the cases where the resulting computational cost is prohibitive, we present in the subsequent proposition larger constants K3 (x0 [M̂ ], M̂ ), K4 , and K5 that are simpler to compute.
Algorithms for computing these constants are discussed in Appendix E. The constant K4 is
obtained by applying a union bound to (10), whereas K3 is obtained by applying a more refined ‘partial’ union bound. [More precisely, for M ∈ M the complement of the probability in
(10) (with K1 (x0 ) replaced by a generic variable t) can be expressed as in (23) in Appendix B.
For given M ∈ M, and after conditioning on the variance estimator (represented by G there),
we apply a union bound by decomposing the maximum over M into a maximum over the
submodels of the given M and a maximum over the models not nested in M . A further union
bound is applied to the latter group of models, giving rise to the bound (24) in Appendix B.
Inspection of this bound shows that the probability appearing in (12) below springs from the
submodels of M , whereas the models not nested in M give rise to the term in (12) involving
the Beta-distribution function.]
∗
For x0 ∈ Rp and M ∈ M define now the distribution function FM,x
for t ≥ 0 via
0
 

1, Pr maxM∗ ∈M,M∗ ⊆M s̄0M∗ V > t
∗
FM,x0 (t) = 1 − min
(12)
+c (M, M) 1 − FBeta,1/2,(d−1)/2 t2
∗
and via FM,x
(t) = 0 for t < 0. Here c (M, M) denotes the number of models M∗ ∈ M
0
that satisfy M∗ * M , V is a random vector that is uniformly distributed on the unit sphere
in the column space of X, and FBeta,1/2,(d−1)/2 denotes the Beta (1/2, (d − 1)/2)-distribution
function, with the convention that in case d = 1 we use FBeta,1/2,0 to denote the distribution
function of pointmass at 1. In view of our assumptions on M it follows that c (M, M) ≥ 1
always holds, except in the case where M = {1, ..., p} (and when this set belongs to M). Next
define the distribution function FM,x0 via
∗
FM,x0 (t) = EG FM,x
(t/G) ,
0

(13)

where G denotes a nonnegative random variable such that G2 /d follows an F -distribution with
(d, r)-degrees of freedom and EG represents expectation w.r.t. the distribution of G. We stress
that FM,x0 depends on x0 only through x0 [M ], and hence the same is true for the constant
K3 (x0 [M ], M ) we define next: For any x0 ∈ Rp and any M ∈ M define K3 (x0 [M ], M ) to be
the smallest constant K satisfying
FM,x0 (K) ≥ 1 − α.

(14)

Furthermore, set K4 = K3 (x0 [∅], ∅). Finally, K5 is the Scheffé constant, i.e., the (1 − α)quantile of G (Scheffé (1959)); see the corresponding discussion in Section 4.8 of Berk et al.
(2013a). Recall that 1 − α ∈ (0, 1) has been assumed.
9

Proposition 2.3. Let x0 ∈ Rp be arbitrary. Then we have the following:
(a) K3 (x0 [M ], M ) exists and is well-defined. If M = {1, ..., p} ∈ M and x0 = 0, then
K3 (x0 [M ], M ) = 0 (and FM,x0 is the c.d.f. of pointmass at zero). If M 6= {1, ..., p} or
x0 6= 0, then (i) 0< K3 (x0 [M ], M ) < ∞ holds, and (ii) equality holds in (14) if and only if
K = K3 (x0 [M ], M ).
(b) For every M ∈ M we have
K2 (x0 [M ], M ) ≤ K3 (x0 [M ], M ) ≤ K4 ≤ K5 .

(15)

K2 (x0 [M 2 ], M 2 ) ≤ K2 (x0 [M 1 ], M 1 ),

(16)

K3 (x0 [M 2 ], M 2 ) ≤ K3 (x0 [M 1 ], M 1 )

(17)

Furthermore,

hold whenever M1 ⊆ M2 , Mi ∈ M.
It is obvious that K3 (x0 [M ], M ) depends, besides x0 [M ] and M , only on r, X, α, and
M, whereas K4 only depends on r, d, α, and M, and K5 depends only on r, d, and α. [Like
with K1 (x0 ), also the other constants introduced depend on the estimator σ̂ only through
r.] We shall write K3 (x0 [M ], M, r), K4 (r), and K5 (r) if we want to stress dependence on
r. Note that K1 (x0 ) = K3 (x0 [Mf ull ], Mf ull ) = K3 (x0 , Mf ull ), provided Mf ull := {1, ..., p}
belongs to M, and that K3 (x0 [M ], M ) = K4 holds for any M ∈ M satisfying |M | = 1 and
s̄M 6= 0. [Indeed, in this case, the probability appearing in (12) equals 1 − FBeta,1/2,(d−1)/2 t2
as can be seen from the proof of Proposition 2.3.] Similarly, K3 (x0 [M ], M ) = K4 holds for
any M ∈ M in case d = 1 as is not difficult to see. The proof of the inequalities involving
the constants K3 and K4 in the above proposition is an extension of an argument in Berk
et al. (2013b) (not contained in the published version Berk et al. (2013a)) to find – in the
case p = d – an upper-bound for their PoSI constant that does not depend on X, but only on
d. [Note that K4 is a counterpart to Kuniv in Berk et al. (2013b).] Inequalities (16) and (17)
simply reflect the fact that observing only x0 [M ] implies that fewer information about x0 is
provided for smaller models M . As a consequence of these inequalities it is possible that, on
the event where a small model M1 is selected, the resulting confidence interval is larger than
it is on the event where a larger model M2 is selected. Again, this simply reflects the fact
that less information on x0 is available under the smaller model. Note, however, that the just
discussed phenomenon is counteracted by the fact that the length of the confidence interval
also depends on ||sM || and that we have ||sM1 || ≤ ||sM2 || for M1 ⊆ M2 ; cf. Figure 1 in Section
4.
Proposition 2.3 implies that (15) holds with M̂ replacing M , which together with Corollary
2.2 immediately implies the following result. We stress that the confidence intervals figuring
in the subsequent corollary depend on x0 only through x0 [M̂ ] and thus are feasible in case
(ii) discussed at the beginning of Section 2.2.
Corollary 2.4. Let M̂ be an arbitrary model selection procedure with values in M, and let
x0 ∈ Rp be arbitrary. Then the confidence interval (5) with K(x0 , M̂ ) replaced by K3 (x0 [M̂ ], M̂ )
(K4 , or K5 , respectively) satisfies the coverage property (7).
We conclude this section with a few remarks regarding extensions.
10

Remark 2.5. (Infeasible variance estimators) (i) For later use we note that all results derived
in Section 2 continue to hold if σ̂ 2 is allowed to also depend on σ but otherwise satisfies the
assumptions made earlier (e.g., if σ̂ 2 = σ 2 Z/r where Z is an observable chi-square distributed
random variable with r degrees of freedom that is independent of PX Y ).
(ii) If we set σ̂ 2 = σ 2 and r = ∞, all of the results derived in Section 2 continue to hold
with obvious modifications. In particular, in Proposition 2.3 the random variable G2 then
follows a chi-squared distribution with d degrees of freedom. We shall denote the constants
corresponding to K1 (x0 ), K2 (x0 [M ], M ), K3 (x0 [M ], M ), K4 , and K5 obtained by setting
σ̂ 2 = σ 2 and r = ∞ by K1 (x0 , ∞), etc. We stress that these constants do not depend on σ.
Remark 2.6. (i) All results carry over immediately to the case where µ can vary only in a
subset M of Rn .
(ii) We have assumed that any non-empty M ∈ M is of full-rank. This assumption could
easily be dropped, but this would lead to more unwieldy results.
(iii) Since the development in Section 2 is based on the bound (9), it is obvious that all
results in Section 2 also hold if M̂ = M̂ (X, Y, σ̄ 2 ) for some arbitrary estimator σ̄ 2 , that may
differ from the estimator σ̂ 2 that governs the length of the confidence intervals considered.

2.3

On the merits of the non-standard targets

Remark 2.7. (i) As already noted, the (non-standard) coverage target in Berk et al. (2013a)
(n)
is β
(where these authors choose to represent it in what they call ‘full model indexing’).
M̂

While β

(n)
M̂

has a clear technical meaning as the coefficient vector that provides the best

approximation of µ by elements of the form X[M̂ ]γ w.r.t. the Euclidean distance, adopting
this quantity as the target for inference confronts one with the fact that the target then
depends on the data Y via M̂ (implying that the target as well as its dimension are random);
(n)
furthermore, different model selection procedures give rise to different targets β . Also note
that, e.g., the meaning of the first component of the target β
model M̂ . The target x00 [M̂ ]β

(n)
M̂

(n)
M̂

M̂

depends on the selected

considered in this paper, while again being random and

sharing many of the properties of β

(n)
M̂

just mentioned, seems to be somewhat more amenable
P
(n)
to interpretation since it is simply the random convex combination M x00 [M ]β M 1(M̂ = M )
(n)
of the (infeasible) predictors x00 [M ]β M (which one would typically use if model M is forced
upon one for prediction and which all have one and the same dimension, not depending on
the data).
(ii) In the classical case, i.e., when µ = Xβ and d = p ≤ n, one can justly argue that the
(n)
target for inference should be x00 β rather than x00 [M̂ ]β
because x00 β is a better (infeasible)
M̂

(n)

predictor in the mean-squared error sense than is x00 [M̂ ]β provided y0 is independent of M̂
M̂
(which will certainly be the case if y0 is independent of Y and σ̂ 2 , or if y0 is independent of Y
and M̂ is only a function of X and Y ). [This is so since the mean-squared error of prediction
(n)
of x00 β is not larger than the one of x00 [M ]β M for every M and since M̂ is independent of y0 .]
However, this argument does not apply if x0 is not observed in its entirety, but only x0 [M̂ ] is
11

observed, because then x00 β is not available. In this case we thus indeed have some justification
(n)
for the target x00 [M̂ ]β even in the classical case. This is in contrast with the situation when,
M̂
as in Berk et al. (2013a), one’s interest focusses on parameters rather than predictors: Similar
as before one can argue that in the classical case the true parameter β should be the target
(n)
(n)
rather than β
but there seems now to be little to justify the non-standard target β
(as
M̂

the preceding argument justifying the target x00 [M̂ ]β
not applicable to the target β

(n)
).
M̂

(n)
M̂

M̂

even in the classical case is obviously
(n)

(iii) In view of the preceding discussion it seems that the non-standard target β of Berk
M̂
et al. (2013a) mainly has a justification in a non-classical setting where µ is not assumed
to belong to the column space of X (implying d < n), or where d < p holds (subsuming in
particular the important case p > n = d), because in these cases β is no longer available as
a target (being not defined or not uniquely defined). However, in a setting, where µ is not
assumed to belong to the column space of X or where p > n = d holds, the assumption on the
variance estimator σ̂ 2 made in Berk et al. (2013a) (as well as in the present paper) becomes
problematic and quite restrictive; see Remark 2.1(ii) in Leeb et al. (2015) as well as Appendix
(n)
(n)
A. Hence, there is some advantage in considering the targets x00 [M̂ ]β
rather than β
as
M̂
M̂
the former has a justification in the classical as well as in the non-classical framework.
(iv) We note the obvious fact that if the target of inference is the standard target x00 β
(assuming the classical case) then the reasoning underlying Proposition 2.1 does not apply
since the difference between the post-model-selection predictor and the standard target is not
independent of β. For the same reason the approach in Berk et al. (2013a) cannot provide a
solution to the problem of constructing confidence sets for the standard target β.
Remark 2.8. (On the optimality of the design-dependent target) (i) The infeasible predictor
(n)
x00 [M ]β M (for fixed M ) is the best predictor for y0 in the mean-squared error sense among
all predictors of the form x00 [M ]γ in case y0 | ν, x0 ∼ N (ν, σ 2 ) and (ν, x00 ) is drawn from the
empirical distribution of (µi , x0i ) where x0i denotes the i-th row of X (‘in-sample prediction’).
[More generally, this is so if (ν, x00 ) is drawn from the empirical distribution of (µi + ai , x0i )
where a is a fixed vector orthogonal to the column space of X.] Otherwise, it does in general
not have this optimality property (but nevertheless its feasible counterpart x00 [M ]β̂ M would
typically be used if one is forced to base prediction on model M ).
(n)
(ii) The optimality property in (i) carries over to the design-dependent target x00 [M̂ ]β
M̂

provided (y0 , x00 )0 is independent of M̂ .

2.4

Behavior of the constants Ki as a function of p

In this section we provide some results on the size of the constants Ki that govern the length
of the confidence intervals. In particular, these results help in answering the question how
tight a bound for K1 and K2 is provided by K3 or K4 .

12

2.4.1

Orthogonal designs

Berk et al. (2013a) show that in the case p = d ≤ n their PoSI constant becomes smallest
for the case of orthogonal design (provided √
the model universe M is sufficiently rich, e.g., M
contains all submodels) and then has rate log p as p → ∞, at least in the known-variance
case; cf. Proposition 5.5 in Berk et al. (2013a) (where the error term o(d) given in this result
should read o(1)). In the next proposition we study the order of magnitude of K1 (x0 ), the
analogue of the PoSI constant and of the closely related constant K2 (x0 [M ], M ) in the case
of orthogonal design. Recall that K1 (x0 ) is only feasible if x0 is observed in its entirety, while
K2 (x0 [M ], M ) is the ideal bound for K1 (x0 ) given only knowledge of x0 [M ]. Note that in
the following result some of the objects depend on p, but we do not always show this in the
notation. Furthermore, φ and Φ denote the p.d.f. and c.d.f. of a standard normal variable,
respectively, and kxk0 denotes the l0 -norm.
Proposition 2.9. Consider the known-variance case (i.e., r = ∞ and σ̂ 2 = σ 2 ) and assume
that for every p ≥ 1 the model universe M used is the power set of {1, ..., p}. Let α, 0 < α < 1,
be given, not depending on p.
(a) For any p ≥ 1 let X = X(p) be an n(p)×p matrix with (non-zero) orthogonal columns.
For any such sequence X one can find a corresponding sequence of (non-zero) p × 1 vectors
x0 such that K1 (x0 , ∞) = K1 (x0 , ∞, X, α, M) satisfies
√
lim inf K1 (x0 , ∞)/ p ≥ ξ
p→∞

p
where ξ = supb>0 φ(b)/ 1 − Φ(b) ≈ 0.6363. Furthermore, for any sequence X as above one
can find another sequence of (non-zero) p × 1 vectors x0 such that K1 (x0 , ∞) = O(1) (for
example, any sequence of (non-zero) p × 1 vectors x0 satisfying supp kx0 k0 < ∞ will do).
(b) Let γ ∈ [0, 1) be given. Then K2 (x0 [M ], M, ∞) = K2 (x0 [M ], M, ∞, X, α, M) satisfies
p
√
lim inf inf p inf
inf
K2 (x0 [M ], M, ∞)/ p ≥ ξ 1 − γ,
p→∞ x0 ∈R X∈X(p) M ∈M,|M |≤γp

where X(p) =

S

n≥p {X

: X is n × p with non-zero orthogonal columns}.

The lower bounds given in the preceding proposition clearly also apply to K3 (x0 [M ], M, ∞)
and K4 (∞) a fortiori. Part (a) of the above proposition shows that, even in the orthogonal
√
case, the growth of K1 (x0 , ∞) is – in the worst-case w.r.t. x0 – of the order p. This is in
contrast to the above mentioned result of Berk et al. (2013a) for the PoSI constant. Part (a)
also shows that there are other choices for x0 such that K1 (x0 , ∞) stays bounded. In this
context also recall that K1 (x0 , ∞) with x0 equal to a p×1 standard basis vector coincides with
a PoSI1 constant and thus equals the (1 − α)-quantile of the distribution of the absolute value
of a standard normal variable in the orthogonal case. Part (b) goes on to show that regardless
√
of x0 and X the growth of the constants K2 (x0 [M ], M, ∞) is of the order p (except perhaps
for very large submodels M ).
2.4.2

Order of magnitude of K3 and K4

The next proposition, which exploits results in Zhang (2013), shows that K4 (∞) is a tight
upper bound for K3 (x0 [M ], M, ∞) at least if p is large. It also provides the growth rates for
13

K4 (∞) and K3 (x0 [M ], M, ∞). As before, the dependence of several objects on p (or n) will
not always be shown in the notation. For the following recall the constants c (M, M) defined
after (12).
Proposition 2.10. Consider the known-variance case (i.e., r = ∞ and σ̂ 2 = σ 2 ) and assume
that
S for every p ≥ 1 a (non-empty) model universe M = Mp is given that satisfies (i)
{M : M ∈ M} = {1, ..., p}, (ii) ∅ ∈ M, (iii) c (M, M) ≥ τ |M| for every M ∈ M with
M 6= {1, . . . , p}, where τ > 0 is a given number (neither depending on M , M, nor p), and
(iv) |M| → ∞ as p → ∞. For n ∈ N, the set of positive integers, let Xn,p (M) denote the set
of all n × p matrices of rank min(n, p) with the property that X[M ] has full column-rank for
every ∅ 6= M ∈ M. Furthermore, let α, 0 < α < 1, be given (neither depending on p nor n).
Let n(p) ∈ N be a sequence such that n(p) → ∞ for p → ∞ and such that Xn(p),p (M) 6= ∅ for
every p ≥ 1. Then we have
lim

sup

sup

sup

p→∞ M ∈M,M 6={1,...,p} x ∈Rp X∈X
0
n(p),p (M)

|1 − (K3 (x0 [M ], M, ∞)/K4 (∞))| = 0,

(18)

where K3 (x0 [M ], M, ∞) = K3 (x0 [M ], M, ∞, X, α, M) and K4 (∞) = K4 (∞, min(n(p), p), α, M).
Furthermore,
r


K4 (∞)/ min(n(p), p) 1 − |M|−2/(min(n(p),p)−1) → 1
as p → ∞.
Remark 2.11. (i) Xn(p),p (M) 6= ∅ implies Xn,p (M) 6= ∅ for n ≥ n(p).
(ii) Xn(p),p (M) is certainly non-empty for n(p) ≥ p, but – depending on M – this can
already be true for n(p) much smaller than p.
The assumptions (i)-(iv) on M in the preceding proposition are shown in the next corollary
to be always satisfied in the important case where M is of the form {M ⊆ {1, . . . , p} : |M | ≤ mp }.
Furthermore, in the special case where M is the universe of all submodels, a simple formula
for the growth rate of K4 (∞) is found.
Corollary 2.12. Consider the known-variance case (i.e., r = ∞ and σ̂ 2 = σ 2 ) and let α,
0 < α < 1, be given (neither depending on p nor n). Let mp ∈ N satisfy 1 ≤ mp ≤ p for every
p ≥ 1 and define the set M(mp ) = {M ⊆ {1, . . . , p} : |M | ≤ mp }. Then M(mp ) satisfies
(i)-(iv) in Proposition 2.10 with τ = 1/3. Consequently, for n(p) as in Proposition 2.10, (18)
holds with M replaced by M(mp ) and
v


u
mp  !−2/(min(n(p),p)−1)
u
X
p
u
→1
K4 (∞)/tmin(n(p), p) 1 −
k
k=0

as p → ∞. In particular, if mp = p for all p ≥ 1, we necessarily have n(p) ≥ p and
√
√
K4 (∞)/ p → 3/2
as p → ∞.
14

In the important case, where p = d ≤ n and M is the entire power set of {1, ..., p}, the
preceding corollary shows that K4 (∞) (and hence a fortiori all the constants K1 (x0 , ∞),...,
K3 (x0 [M ], M, ∞)) are ‘bounded away’ from the Scheffé constant K5 which clearly satisfies
√
K5 / p → 1 for p → ∞. This is in line with a similar finding in Berk et al. (2013a), Section
6.3, for their PoSI constant.
Remark 2.13. In the proof of Proposition 2.3 union bounds were used to obtain the results
for K3 (x0 [M ], M ) and K4 . Hence, one might ask whether or not these constants as bounds
for K2 (x0 [M ], M ) are overly conservative. We now collect evidence showing that improving
K3 (x0 [M ], M ) and K4 will not be easy and is sometimes impossible: First, Lemma B.4 in
Appendix B shows that there exist n×p design matrices X with p = d = 2 and vectors x0 such
that K4 = K1 (x0 ) in case M is the universe of all submodels. Hence, in this case the union
bounds used in the proof of Proposition 2.3 are all exact. Furthermore, in the known-variance
case with p = d ≤ n and where M √
again is the universe of all submodels, the propositions
√
√
√
given above entail that K4 (∞) ∼ p 3/2 ≈ 0.866 p while K1 (x0 , ∞)  ξ p with ξ ≈ 0.6363
is possible; e.g., as the worst-case behavior in the orthogonal case, or with x0 = ei and the
design matrices constructed in the proof of Theorem 6.2 in Berk et al. (2013a) (recall that
K1 (ei , ∞) coincides with a PoSI1 constant). This again shows that there is little room for
improving K3 and K4 . [Further evidence in that direction is provided by the observation
√
√
that the proof of Theorem 6.3 in Berk et al. (2013a) implies that K1∗ / p tends to 3/2 in
probability as p → ∞, where K1∗ is an analogue of K1 (x0 , ∞) that is obtained from (10) (with
r = ∞) after replacing the vectors s̄M by 2p independent random vectors, each of which is
uniformly distributed on the unit sphere of the column space of X (and these vectors being
independent of Y ). In other words, if one ignores the particular structure of the vectors s̄M ,
then the bound K4 (∞) is close to being sharp for large values of p.]
Remark 2.14. The results for p → ∞ in this subsection as well as the related results in Berk
et al. (2013a) should be taken with a grain of salt as they obviously are highly non-uniform
w.r.t. α: Note that – for fixed n and p – any one of the constants Ki will vary in the entire
interval (0, ∞) as α varies in (0, 1) (except for degenerate cases), while the limits in the results
in question do not depend on α at all.

3

Confidence intervals for the design-independent non-standard
target

In this section we again consider the model (1), but now assume that µ = Xβ for some
unknown β ∈ Rp holds and that the n × p matrix X is random, with X independent of U ,
where U again follows an N 0, σ 2 In -distribution with 0 < σ < ∞. We also assume that X
has full column rank almost surely (implying p ≤ n) and that each row of X is distributed
according to a common p-dimensional distribution L (not depending on n) with a finite
and positive definite matrix of (uncentered) second moments, which we denote by Σ. [We
shall refer to the preceding assumptions as the maintained model assumptions of this section.]
Furthermore, we assume again that we have available an estimator σ̂ 2 such that, conditionally
on X, σ̂ 2 is independent of PX Y (or, equivalently, of β̂ = (X 0 X)−1 X 0 Y ) and is distributed as
15

σ 2 /r times a chi-squared distributed random variable with r degrees of freedom (1 ≤ r < ∞).
The collection M of admissible models will be assumed to be the power set of {1, . . . , p} in
this section for convenience, but see Remark 3.8 for possible extensions. Observe that all the
results of Section 2 remain valid in the setup of the present section if formulated conditionally
on X (and if x0 is treated as fixed). [Alternatively, if x0 is random but independent of X, U ,
and σ̂ 2 , the same is true if the results in Section 2 are then interpreted conditionally on X
and x0 .] The joint distribution of Y , X, and σ̂ 2 (and of σ̃ appearing below) will be denoted
by Pn,β,σ (see also Appendix D.4).
In this section we shall consider asymptotic results for n → ∞ but where p is held constant
(for an extension to the case where p is allowed to diverge with n see Appendix D.3). It is thus
important to recall that all estimators, estimated models, etc. depend on sample size n. Also
note that r may depend on sample size n. We shall typically suppress these dependencies on
n in the notation. Furthermore, we note that, while not explicitly shown in the notation, the
rows of X and U (and thus of Y ) may depend on n. [As the results in Section 2 are results
for fixed n, this trivially also applies to the results in that section.] However, recall that L,
and hence Σ, are not allowed to depend on n.
If M1 and M2 are subsets of {1, ..., p} and if Q is a p × p matrix we shall denote by
Q[M1 , M2 ] the matrix that is obtained from Q by deleting all rows i with i ∈
/ M1 as well as
all columns j with j ∈
/ M2 ; if M1 is empty but M2 is not, we define Q[M1 , M2 ] to be the
1 × |M2 | zero vector; if M2 is empty but M1 is not, we define Q[M1 , M2 ] to be the |M1 | × 1
zero vector; and if M1 = M2 = ∅ we set Q[M1 , M2 ] = 0 ∈ R.
To motivate the target studied in this section, consider now the problem of predicting
a

0
2
new variable y0 = x0 β + u0 where x0 , u0 , X, and U are independent and u0 ∼ N 0, σ . For
(?)

a given model M ⊆ {1, ..., p} we consider the (infeasible) predictor x00 [M ]β M where
(?)

β M = β[M ] + (Σ[M, M ])−1 Σ[M, M c ]β[M c ],
with the convention that the inverse is to be interpreted as the Moore-Penrose inverse in
(?)
(?)
case M = ∅. Note that x00 [M ]β M = 0 if M = ∅ and that x00 [M ]β M = x00 β if M =
{1, . . . , p}. A justification for considering this infeasible predictor is given in Remark 3.2
below. For purpose of comparison we point out that, under the assumption µ = Xβ
(n)
(n)
maintained in the present section, β M defined in (3) can be rewritten as β M = β[M ] +

−1
X[M ]0 X[M ]
X[M ]0 X[M c ]β[M c ]. Given a model selection procedure M̂ = M̂ (X, Y, σ̂ 2 )
we define now the (infeasible) predictor
x00 [M̂ ]β

(?)
M̂

as our new target for inference. We call this target the design-independent (non-standard)
target as it does not depend on the design matrix X beyond its dependence on M̂ . We discuss
its merits in the subsequent remarks.
Remark 3.1. As in Remark 2.7(ii) one can argue that the target for inference should be
(?)
(?)
x00 β rather than x00 [M̂ ]β because again x00 β is a better (infeasible) predictor than x00 [M̂ ]β
M̂

M̂

provided that (x00 , u0 ) is independent of M̂ (which, in particular, will be the case if (x00 , u0 ) is
16

independent of X, U , and σ̂, or if (x00 , u0 ) is independent of X, U and M̂ is only a function
of X and Y ). But again, this argument does not apply if x0 is not observed in its entirety,
but only x0 [M̂ ] is observed.
Remark 3.2. (On the optimality of the design-independent target) (i) Assume that additionally x00 ∼ L. If we are forced to use the (theoretical) predictors of the form x00 [M ]γ, then
(?)
straightforward computation shows that x00 [M ]β M provides the smallest mean-squared error
of prediction among all the linear predictors x00 [M ]γ. [Note that this result corresponds to
the observation made in Remark 2.8 with L corresponding to the empirical distribution of
the rows of X.] If, furthermore, x0 is normally distributed, then x0 and u0 are jointly normal
(?)
and thus x00 [M ]β M is the conditional expectation of y0 given x0 [M ] and hence is also the best
predictor in the class of all predictors depending only on x0 [M ].
(?)
(ii) Again assume that x00 ∼ L. The discussion in (i) implies that x00 [M̂ ]β has a meanM̂

squared error of prediction not larger than the one of x00 [M̂ ]γ(M̂ ) for any choice of γ(M̂ ),
provided (x00 , u0 ) is independent of M̂ . If, additionally, x0 is normally distributed, then
(?)
x00 [M̂ ]β is also the best predictor in the class of all predictors depending only on x00 [M̂ ] and
M̂

M̂ .
After having motivated the design-independent target, we shall, in the remainder of this
section, treat x0 as fixed (but see Remark D.2 in Appendix D.2 for the case where x0 is
random). We now proceed to show that the confidence intervals constructed in Section 2 are
(?)
also valid as confidence intervals for the design-independent target x00 [M̂ ]β in an asymptotic
M̂
sense under some mild conditions. While the results in Section 2 apply to any model selection
procedure whatsoever (in case that M is the power set of {1, . . . , p} as is the case in the
present section), we need here to make the following mild assumption on the model selection
procedure.
Condition 3.3. The model selection procedure satisfies: For any M ⊆ {1, . . . , p} with |M | < p
and for any δ > 0,
n
o
sup Pn,β,σ (M̂ = M |X) : β ∈ Rp , σ > 0, kβ[M c ]k /σ ≥ δ → 0
in probability as n → ∞.
Condition 3.3 is very mild and typically holds for model selection procedures such as AICand BIC-based procedures as well as Lasso-type procedures. [This can be established along
the lines of the proof of Corollary 5.4(a) in Leeb and Pötscher (2003).] In addition, we assume
the following condition on the behavior of the design matrix.
√
Condition 3.4. The sequence of random matrices n [(X 0 X/n) − Σ] is bounded in probability.
Condition 3.4 holds, for example, when the rows of X are independent, or weakly dependent, and when the distribution L has finite fourth moments for all its components. We also
introduce the following condition.
17

Condition 3.5. The degrees of freedom parameters r of the sequence of estimators σ̂ 2 satisfy
r → ∞ as n → ∞.
Of course, if we choose for σ̂ 2 the usual variance estimator σ̂ 2OLS then this condition is
certainly satisfied with r = n − p. We are now in the position to present the asymptotic
coverage result. Recall that the confidence intervals corresponding to Ki with 2 ≤ i ≤ 5
depend on x0 only through x0 [M̂ ] (or not on x0 at all).
Theorem 3.6. Suppose Conditions 3.3 and 3.4 hold.
(a) Suppose also that Condition 3.5 is satisfied. Let CI(x0 ) be the confidence interval (5)
where the constant K(x0 , M̂ ) is given by the constant K1 (x0 , r) defined in Section 2. Then
the confidence interval CI(x0 ) satisfies


(?)
0
inf
P
x
[
M̂
]β
∈
CI(x
)
X
≥ (1 − α) + op (1),
(19)
0
n,β,σ
0
p
p
M̂

x0 ∈R ,β∈R ,σ>0

where the op (1) term above depends only on X and converges to zero in probability as n → ∞.
Relation (19) a fortiori holds if the confidence interval CI(x0 ) is based on the constants
K2 (x0 [M̂ ], M̂ , r), K3 (x0 [M̂ ], M̂ , r), K4 (r), or K5 (r), respectively.
(b) Let σ̃ be an arbitrary estimator satisfying
sup
β∈Rp ,σ>0

p

Pn,β,σ (|σ̃/σ − 1| ≥ δ |X ) → 0

(20)

for any δ > 0 as n → ∞. Let further r∗ = rn∗ be an arbitrary sequence in N∪ {∞} satisfying
r∗ → ∞ for n → ∞. Let CI ∗ (x0 ) denote the modified confidence interval which is obtained by
replacing σ̂ by σ̃ and K(x0 , M̂ ) by K1 (x0 , r∗ ) (K2 (x0 [M̂ ], M̂ , r∗ ), K3 (x0 [M̂ ], M̂ , r∗ ), K4 (r∗ ),
or K5 (r∗ ), respectively) in (5) (while keeping M̂ unchanged). Then relation (19) holds with
CI(x0 ) replaced by CI ∗ (x0 ).
Theorem 3.6(a) shows that for any x0 ∈ Rp the interval CI(x0 ) is an asymptotically
valid confidence interval for the design-independent target and additionally that the lower
bound (1 − α) + op (1) for the minimal (over β and σ) coverage probability can be chosen
independently of x0 . Theorem 3.6(b) extends this result to a larger class of intervals. [Note
that Part (a) is in fact a special case of Part (b) obtained by setting σ̃ = σ̂ and r∗ = r
and observing that σ̂ clearly satisfies the condition on σ̃ in Part (b) under Condition 3.5.]
We note that applying Theorem 3.6(b) with σ̃ = σ̂ and r∗ = ∞ shows that Theorem 3.6(a)
also continues to hold for the confidence interval that is obtained by replacing the constants
K1 (x0 , r) (K2 (x0 [M̂ ], M̂ , r), K3 (x0 [M̂ ], M̂ , r), K4 (r), or K5 (r), respectively) by the constants
K1 (x0 , ∞) (K2 (x0 [M̂ ], M̂ , ∞), K3 (x0 [M̂ ], M̂ , ∞), K4 (∞), or K5 (∞), respectively). Measurability issues regarding Theorem 3.6 are discussed in Appendix D.1.
Condition (20) is a uniform consistency property. It is clearly satisfied by σ̂ 2OLS (and
more generally by the estimator σ̂ 2 under Condition 3.5 as already noted above), but it is also
satisfied by the post-model-selection estimator σ̂ 2M̂ = ||Y −X[M̂ ]β̂ M̂ ||2 /(n−|M̂ |) provided the
model selection procedure satisfies Condition 3.3, see Lemma C.2 in Appendix C for a precise
result. As a consequence, Theorem 3.6(b) shows that the post-model-selection estimator σ̂ 2M̂
can be used instead of σ̂ 2 in the construction of the confidence interval.
18

Remark 3.7. (Infeasible variance estimators) Theorem 3.6(a) remains valid if σ̂ 2 is allowed
to depend also on σ but otherwise satisfies the assumptions made earlier or if σ̂ 2 = σ 2 and
r = ∞. Similarly, Theorem 3.6(b) remains valid if σ̃ 2 is allowed to be infeasible. Furthermore,
a remark similar to Remark 2.6(iii) also applies here.
Remark 3.8. (Restricted universe of selected models) Theorem 3.6 can easily be generalized
to the case where a universe M different from the power set of {1, . . . , p} is employed, provided
the full model {1, . . . , p} belongs to M (and M satisfies the basic assumptions made in Section
2).

4

Numerical study

We next present a numerical study of the lengths and the minimal coverage probabilities of
various confidence intervals. We begin, in Section 4.1, with an investigation of the length
of the confidence intervals introduced in Section 2, including the ‘naive’ confidence interval
that ignores the model selection step, as a function of the selected model. In Section 4.2 we
then evaluate numerically the minimal coverage probabilities of these confidence intervals. As
model selectors we consider here AIC, BIC, LASSO, SCAD (Fan and Li (2001)), and MCP
(Zhang (2010)). Finally, in Section 4.3 we compare the confidence intervals introduced in
Section 2 with the confidence interval proposed recently in Lee et al. (2016), which is specific
to the LASSO model selector. Code for the computations in this section is available from the
first author.

4.1

Lengths of confidence intervals

We consider the lengths of the confidence intervals obtained from (5) standardized by σ̂,
i.e., we consider 2K(x0 , M̂ )ksM̂ k for the six cases where K(x0 , M̂ ) is replaced by either one
of the five constants K1 (x0 ), K2 (x0 [M̂ ], M̂ ), K3 (x0 [M̂ ], M̂ ), K4 , K5 of Section 2 or by the
constant Knaive = qr,1−α/2 , the (1 − α/2)-quantile of Student’s t-distribution with r degrees of
freedom. We recall that the constant Knaive yields the ‘naive’ confidence interval that ignores
the model selection step and that we have Knaive ≤ K1 (x0 ) ≤ ... ≤ K5 (the first inequality
holding provided x0 6= 0).
For computing the standardized length, we set α = 0.05, n = 29, d = p = 10, r = n − p,
σ = 1, and obtain X and x0 from a data set of Rawlings et al. (1998) concerning the peak
flow rate of watersheds. This data set contains a 30 × 10 design matrix XRaw corresponding
to ten explanatory variables. For a description of these variables see Appendix F. This data
set is also studied in Kabaila and Leeb (2006) and Leeb et al. (2015). We refer to it as the
watershed data set, and x0 and X are chosen such that (x0 , X 0 )0 is equal to the watershed
design matrix XRaw . It is easily checked that the so-obtained matrix X is indeed of full
column rank (and x0 6= 0). Furthermore, the model universe M is chosen to be the power set
of {1, ..., p}.
For the so chosen values of α, n, p, r, σ, X, x0 , and M, we compute the standardized
lengths 2K(x0 , M )ksM k of the confidence intervals obtained by replacing K(x0 , M ) by Knaive ,
K1 (x0 ), K2 (x0 [M ], M ), K3 (x0 [M ], M ), K4 , and K5 , respectively. To ease the computational
19

3

4

K5
K4
K3
K2
K1
Knaive

0

1

2

length

5

6

7

burden and to enable a simple presentation as in Figure 1 below, we compute the standardized
lengths of the confidence intervals only for M belonging to the family {{1}, ..., {1, ..., 10}}
consisting of ten nested submodels. [This does not mean that we compute the constants
Ki under the assumption of a restricted universe of models; recall that we use M equal to
the power set of {1, ..., p}.] The computation of Knaive , K1 (x0 ), K3 (x0 [M ], M ), K4 , and
K5 is either straightforward or is obtained from the algorithms described in Appendix E.
However, computing K2 (x0 [M ], M ) for M 6= {1, ..., 10} necessitates to compute sup{K1 (x) :
x[M ] = x0 [M ]}. We approximate this supremum by using a three-step Monte Carlo procedure
described in Appendix F.

2

4

6

8

10

submodel size

Figure 1: Standardized lengths of various confidence intervals as function of model
size. Dashed lines are added to improve readability.
The standardized lengths of the confidence intervals corresponding to the constants Knaive ,
K1 ,. . . , K5 are reported in Figure 1 for the ten nested submodels mentioned before. We first
see that, for each of the constants Knaive , K1 , K4 , and K5 , the standardized length of the
confidence interval increases with submodel size, which must hold since these constants do
20

not depend on the submodel M and since the term ||sM || increases with submodel size (for
nested submodels as considered in Figure 1). However, as discussed after Proposition 2.3, the
values of K2 and K3 decrease with increasing submodel size for nested submodels. Figure 1
shows that the combined effect of the increase of ||sM || and the decrease of K2 and K3 with
submodel size can be an increase or a decrease of the standardized lengths of the confidence
intervals. Indeed, the standardized lengths increase globally (i.e., from submodel size 1 to
10), but can decrease locally (for example, the standardized length of the confidence interval
obtained from K2 decreases from submodel size 6 to submodel size 8; for the interval obtained
from K3 the standardized length decreases from submodel size 9 to submodel size 10). In
Figure 1 the decreases of the standardized lengths occur only between submodel sizes for
which ||sM || is almost constant with M (which can be seen from the standardized lengths
obtained from, say, K5 , since they are proportional to ||sM ||). We also see from Figure 1 that
the ‘naive’ interval is much shorter than the other intervals (at the price of typically not having
the correct minimal coverage probability). The difference in standardized length between the
intervals based on K1 and K2 , respectively, is noticeable but not dramatic. A larger increase
in standardized length is noted when comparing the interval based on the costly-to-compute
constant K2 with the one obtained from K3 , especially for submodel sizes 6 to 9. Furthermore,
the standardized lengths of the confidence intervals obtained from K3 are very close to those
obtained from K4 for model size 1 to 9; cf. (18). Finally, in Figure 1 we also see that the
confidence intervals obtained from K1 , K2 , and K3 have the same standardized length when
the model size is 10, and that the same is true for the confidence intervals obtained from K3
and K4 when the model size is 1. This, of course, is not a coincidence, but holds necessarily
as has been noted in the discussion of Proposition 2.3.
Additional computations of confidence interval lengths, with X and x0 now randomly generated, yield results very similar to those in Figure 1. For the sake of brevity, these results are
not shown here. We find, in particular, that the standardized length of the confidence interval
obtained from K3 always increases with submodel size when they are averaged with respect
to X and x0 , but, as in Figure 1, can decrease locally when not averaged. [In these additional
numerical studies we did not consider the constant K2 due to the high computational cost
involved in its evaluation.]

4.2

Minimal coverage probabilities

In this section we consider the case where µ = Xβ and d = p < n, i.e., the case where the
given matrix X has full rank less than n and provides a correct linear model for the data Y .
We then investigate the minimal coverage probabilities (the minimum being w.r.t. β ∈ Rp
and σ ∈ (0, ∞)) of the intervals obtained from the constants Knaive , K1 , K3 , and K4 when
(n)
used as confidence intervals for the target x00 [M̂ ]β on the one hand as well as for the target
M̂

(?)

x00 [M̂ ]β on the other hand. The constants K1 , K3 , and K4 are computed based on M equal
M̂
to the power set of {1, ..., p}. We do not report results for confidence intervals obtained from
K2 , since the computation of K2 is too costly for the study we present below. The results for
confidence intervals obtained from K5 would be qualitatively similar to those for confidence
intervals obtained from K4 , so we do not report them for the sake of brevity.

21

We consider minimal coverage probabilities in the setting where α = 0.05, p = 10, n = 20
or n = 100, and the variance parameter is estimated by the standard unbiased estimator using
the full model, so that r = n − p. For model selection we consider AIC-, BIC-procedures,
the LASSO, SCAD (Fan and Li (2001)), and MCP (Zhang (2010)). Tuning parameters
of the latter three procedures are chosen by cross-validation. For all five procedures we
always protect the first explanatory variable (which corresponds to an intercept term) from
selection. However, note that the information that the first variable is protected is not used
in computing the constants Ki , i.e., we do not use a restricted universe of models but use M
equal to the power set of {1, ..., p}. [Additional simulations with no intercept term and no
protected explanatory variable lead to results very similar to the ones given in Table 1 below.]
Computational details regarding these procedures can be found in Appendix F.
The design matrix X and the vector x0 are generated in the following manner: The 10×10
matrix Σ of (uncentered) second moments is chosen to be of the form


1 0···0


Σ =  0..
,
.
0

Σ̃

where we consider three choices for the 9 × 9 matrix Σ̃. For the first case, Σ̃ is obtained by
removing the first row and column of the 10×10 empirical covariance matrix (standardized by
30−1 = 29) of the variables in the 30×10 watershed design matrix XRaw . For the second case,
we set Σ̃ = Ip̃ + (2a + p̃a2 )Ep̃ with p̃ = 9, a = 10, and with Ep̃ the p̃ × p̃ matrix which has all
entries equal to 1. For the third case Σ̃ coincides with the identity matrix Ip̃ , except
p that the
zero elements in the last row and column of Ip̃ are replaced by the constant c = 0.8/(p̃ − 1)
where p̃ = 9. Similar as in Berk et al. (2013a) and Leeb et al. (2015), we refer to the data
set obtained in the second case as the exchangeable data set (as the covariance matrix Σ̃ is
permutation-invariant), and to the one obtained in the third case as the equicorrelated data
set (as Σ̃ is the correlation matrix of a random vector, the last component of which has the
same correlation with all the other components); see Appendix F for more details. For a given
configuration of n and Σ, we then sample independently n+1 vectors of dimension 10×1 such
that for each of these vectors the first component is 1 and the remaining nine components are
jointly normally distributed with mean zero and covariance matrix Σ̃. The transposes of the
first n of theses vectors now form the rows of the n × p design matrix X, while the (n + 1)-th
of these vectors is used for the p-dimensional vector x0 . [It is easy to see that the mechanism
just described generates matrices of full column rank almost surely. The matrices X actually
generated were additionally checked to be of full column rank.]
Consider now a given configuration of n, Σ, the model selection procedure, the target
(either the design-dependent or the design-independent target), as well as of a matrix X and a
vector x0 that have been obtained in the manner just described. Then we estimate the minimal
(over β and σ) coverage probabilities (conditional on X and x0 ) of the confidence intervals
obtained from the constants Knaive , K1 , K3 , and K4 for the given target under investigation.
The minimal coverage probabilities are estimated by a three-step Monte Carlo procedure
similar to that of Leeb et al. (2015), which is described in detail in Appendix F. We stress here
that the minimal coverage probabilities found by this Monte Carlo procedure are simulation22

based results obtained by a stochastic search over a 10-dimensional parameter space, and thus
only provide approximate upper bounds for the true minimal coverage probabilities.
Table 1 summarizes the estimated minimal coverage probabilities for the various confidence
sets and targets, and for the model-selection procedures and data sets considered in the study.
The conclusions are pretty much the same for the three data sets. First, we observe that,
for n = 20, the differences of minimal coverage probabilities between the design-dependent
and independent targets can be significant, especially for the ‘naive’ intervals and for the
other intervals in case the LASSO, SCAD, or MCP model selectors are used. However, for
n = 100, these differences are very small for all the configurations. This is in line with Lemma
C.1 in Appendix C, which entails that for a large family of model selection procedures, the
difference of coverage probabilities between the two targets vanishes, uniformly in β and σ,
when n increases. For n = 100, the results are thus almost identical for the two targets: For
the five model selection procedures, the confidence intervals obtained from the constants K1 ,
K3 , and K4 are valid, while the ‘naive’ confidence intervals are moderately too short, so that
their minimal coverage probabilities are below the nominal level, with a minimum of 0.84.
For n = 20 and when AIC or BIC is used, the ‘naive’ confidence intervals fail to have the
right coverage probabilities to a somewhat larger extent than in case n = 100. Their minimal
coverage probabilities can be as small as 0.81 for the design-dependent target and 0.74 for
the design-independent target. [Note that, for the design-dependent target, for n = 20 and
n = 100, the coverage probabilities of the ‘naive’ confidence interval are generally smaller
for the equicorrelated data set than for the exchangeable data set. This can possibly be
explained by the fact that Theorems 6.1 and 6.2 in Berk et al. (2013a) suggest that K1 should
be larger for the equicorrelated data set than for the exchangeable data set. Hence, for the
equicorrelated data set, larger confidence intervals seem to be needed to have the required
minimal coverage probability for all model selection procedures.] Furthermore, again for
n = 20 and when AIC or BIC is used, the confidence intervals obtained from the constants
K1 , K3 , and K4 remain valid here for both targets.
However, when n = 20 and the LASSO model selector is used, the results for the
design-independent target are drastically different from those obtained with the AIC- or
BIC-procedures: All confidence intervals have minimal coverage probabilities for the designindependent target that are below, and in most cases significantly below, the nominal level.
The failure of all the confidence intervals is here often more pronounced than the failure of the
‘naive’ confidence intervals when other model selectors are used. Especially for the watershed
data set, the estimated minimal coverage probability is 0.18 for the ‘naive’ interval and 0.48
for the confidence interval based on K1 . The reason for this phenomenon can be traced to
the observation that the LASSO model selector, as implemented here and for the parameters
used in the stochastic search for the smallest coverage probability, selects models that are
significantly smaller than those AIC and BIC select. In particular, the LASSO procedure
often excludes regressors for which the corresponding regression coefficients are not small. In
our simulation study, selecting a small model, that excludes regressors with significant coefficients, makes the difference between the design-dependent and design-independent targets
larger. Since the confidence intervals are designed to cover the former target, they hence have
a hard time to cover the latter when the two targets are significantly different. In other words,
for n = 20 the supremum in the display in Condition 3.3 is not small for the LASSO proce23

Data set

Watershed

Exchangeable

Equicorrelated

n

Model
selector

20
20
20
20
20
100
100
100
100
100
20
20
20
20
20
100
100
100
100
100
20
20
20
20
20
100
100
100
100
100

AIC
BIC
LASSO
SCAD
MCP
AIC
BIC
LASSO
SCAD
MCP
AIC
BIC
LASSO
SCAD
MCP
AIC
BIC
LASSO
SCAD
MCP
AIC
BIC
LASSO
SCAD
MCP
AIC
BIC
LASSO
SCAD
MCP

Target
design-dependent
design-independent
(n)
(?)
0
x0 [M̂ ] β
x0 [M̂ ]0 β
M̂
M̂
Knaive K1
K3
K4 Knaive K1
K3
K4
0.84
0.99 1.00 1.00
0.79
0.97 0.99 0.99
0.84
0.99 1.00 1.00
0.74
0.96 0.98 0.98
0.90
1.00 1.00 1.00
0.18
0.48 0.61 0.61
0.90
0.99 1.00 1.00
0.45
0.77 0.84 0.84
0.89
0.99 1.00 1.00
0.47
0.78 0.85 0.85
0.87
0.99 1.00 1.00
0.88
0.99 1.00 1.00
0.88
0.99 1.00 1.00
0.87
0.99 1.00 1.00
0.88
0.99 1.00 1.00
0.87
0.99 1.00 1.00
0.88
0.99 1.00 1.00
0.88
0.99 1.00 1.00
0.88
0.99 1.00 1.00
0.88
0.99 1.00 1.00
0.83
0.99 1.00 1.00
0.80
0.98 0.99 0.99
0.84
0.99 1.00 1.00
0.76
0.97 0.99 0.99
0.90
1.00 1.00 1.00
0.46
0.86 0.93 0.92
0.91
1.00 1.00 1.00
0.55
0.90 0.94 0.94
0.91
1.00 1.00 1.00
0.54
0.89 0.94 0.94
0.89
0.99 1.00 1.00
0.90
0.99 1.00 1.00
0.90
0.99 1.00 1.00
0.90
0.99 1.00 1.00
0.90
0.99 1.00 1.00
0.90
0.99 1.00 1.00
0.90
0.99 1.00 1.00
0.90
0.99 1.00 1.00
0.90
0.99 1.00 1.00
0.90
0.99 1.00 1.00
0.83
0.99 1.00 1.00
0.79
0.98 0.99 0.99
0.81
0.99 1.00 1.00
0.74
0.98 0.99 0.99
0.88
1.00 1.00 1.00
0.39
0.71 0.79 0.79
0.88
0.99 1.00 1.00
0.67
0.92 0.95 0.96
0.86
0.99 1.00 1.00
0.66
0.93 0.96 0.96
0.84
0.99 1.00 1.00
0.84
0.99 1.00 1.00
0.86
0.99 1.00 1.00
0.86
0.99 1.00 1.00
0.88
1.00 1.00 1.00
0.88
1.00 1.00 1.00
0.88
0.99 1.00 1.00
0.89
1.00 1.00 1.00
0.88
0.99 1.00 1.00
0.89
0.99 1.00 1.00

Table 1: Monte Carlo estimates of the minimal coverage probabilities (w.r.t. β and σ) of
various confidence intervals. The nominal coverage probability is 1 − α = 0.95 and p = 10.

24

dure, so that the asymptotics in Theorem 3.6 does not provide a good approximation for the
finite-sample situation. Finally, for n = 20 and for the design-independent target, the results
for the SCAD and MCP model selectors lie somewhere in between those of the AIC and BIC
and those of the LASSO model selectors. Indeed, for SCAD and MCP, the confidence intervals
often fail to have the required minimal coverage probabilities, but less severely than for the
LASSO. We stress that the preceding conclusions hold for the LASSO, SCAD, and MCP procedures as implemented here where tuning parameters are chosen by cross-validation. Other
implementations of these procedures may of course give different results.
The results in Table 1 concern the coverage probabilities conditional on the design matrix
X and on x0 , and thus depend on the values of X and x0 used. In additional (non-exhaustive)
simulations we have repeated the above analysis for other values of X and x0 and have found
similar results.

4.3

Comparison with the confidence interval of Lee et al. (2016)

In this section we now compare the confidence intervals of Section 2 with a confidence interval
recently introduced in Lee et al. (2016). Again, we consider the case where µ = Xβ and
(n)
d = p < n, and we focus on the design-dependent target x00 [M̂ ]β . As in Lee et al. (2016)
M̂
we consider the known-variance case and set σ = 1 in this section. The confidence interval
of Lee et al. (2016) is dedicated to the LASSO model selector and is given in the R package
accompanying that paper for the case where x0 is a standard basis vector. We hence assume
in the following that x0 is equal to the first standard basis vector e1 . The proposed interval
is then conditionally valid for the design-dependent target in the following sense: Consider
the model selector M̂ obtained by selecting those explanatory variables for which the LASSO
estimator has non-zero coefficients, with the penalty parameter λ in (4.1) of Lee et al. (2016)
being fixed, independently of Y . Then the interval proposed by Lee et al. (2016), which we
¯ satisfies, for any fixed X, for x0 = e1 , and for any fixed M ⊆ {1, ..., p} with
denote by CI,
1 ∈ M,


(n)
¯ M̂ = M = 1 − α,
infp Pn,β,1 x00 [M̂ ]β ∈ CI
(21)
β∈R

M̂

with the convention that the probability in the above display is 1 if Pn,β,1 (M̂ = M ) = 0. The
¯ for a given value of M̂ can be carried out without observing x0 [M̂ c ], which
computation of CI
is also the case for the confidence intervals obtained from K2 , ..., K5 , but not for that obtained
¯ (when the conditioning additionally is also on
from K1 . Furthermore, the computation of CI
the signs, see Lee et al. (2016) for details) entails a cost that grows linearly with p. Thus,
¯ can be implemented for significantly larger values of p than the confidence intervals based
CI
on K1 , ..., K3 currently can be. We note for later use that, in the case x0 = e1 considered
¯ as given in Lee et al. (2016) is not defined on the event that a model M̂
here, the interval CI
is selected that does not contain 1. Hence, we can not speak about unconditional coverage
without amending the definition in Lee et al. (2016). [A possible amendment, consistent with
our conventions and maximizing unconditional coverage among all possible amendments, is
(n)
¯ = {0} on this event. With such an
to recall that x00 [M̂ ]β
= 0 if 1 ∈
/ M̂ and to set CI
M̂
¯ then a fortiori has minimal unconditional coverage probability not less than
amendment, CI

25

1 − α.]
Despite being specific to the LASSO model selector with fixed λ, we nevertheless find
below that the confidence interval of Lee et al. (2016) is not shorter than those based on
K1 , K3 , and K4 (presumably due to the fact that (21) imposes a stricter requirement than
requiring only correct unconditional coverage). In addition, we point out that this confidence
interval can be very sensitive to deviations from the specific model selector it is designed for:
In particular, we show that its coverage can break down, when the LASSO model selector
is used but with a data-dependent penalty parameter λ selected by cross-validation. This is
certainly in stark contrast to the confidence intervals obtained from K1 , K3 , and K4 , that are
valid for any model selection procedure whatsoever.
We first present the results on confidence interval lengths. We conduct a Monte-Carlo
study in the case p = 10, n = 100, and α = 0.05. We use the function fixedLassoInf of the
R package selectiveInference to construct the confidence interval of Lee et al. (2016). In
line with the presentation in that paper, this function is designed for the cases where x0 is
one of the standard basis vectors of Rp and thus we set x00 = (1, 0, ..., 0) as already mentioned
above. We consider two different settings for X and β. In the ‘independent’ setting, we sample
independently 1000 values of X, β, and Y in the following way: We first sample the (transposes
of the) rows of X and the vector β as n + 1 independent draws from the standard Gaussian
distribution on Rp . Given X and β, we then sample Y from the N (Xβ, In )-distribution. For
each of the 1000 values of X, β, and Y so obtained, we run the LASSO model selector, with
λ fixed as a function of X as described at the beginning of Section 7 of Lee et al. (2016).
We use the function glmnet of the R package glmnet to compute the LASSO model selector.
Then, if the first variable is included in the selected model, we record the lengths of the
¯ and of the confidence intervals (for the same target) obtained from
confidence interval CI
K1 , K3 , and K4 (where these three constants are computed with r = ∞ and M equal to the
power set of {1, . . . , p}). If not, we discard the realization of X, β, and Y (this is in line with
¯ is only defined on the event where 1 ∈ M̂ as discussed above). In the
the fact that the CI
‘correlated’ setting we proceed as just described, with the only difference that the rows of X
are sampled according to the Gaussian distribution with mean vector 0 ∈ Rp and covariance
matrix [exp(−|i − j|/10)]10
i,j=1 and that then β is sampled from a random vector b, so that Xb
follows the standard Gaussian distribution within the column space of X. [We note that this
mechanism almost surely generates matrices X that have full column rank.]
The medians and empirical 90%-quantiles of the confidence interval lengths’ distributions
obtained that way are reported in Table 2. The conclusion is that there is no unilateral
hierarchy of the two methods for confidence interval construction (that of Lee et al. (2016)
and that using the intervals based on K1 , K3 , and K4 ) in terms of median length. Depending
on the situation, any of the two methods can provide the smallest median length. The 90%¯ of Lee et al.
quantiles, on the other hand, are always larger for the confidence interval CI
¯
(2016) than for those obtained from K1 , K3 , and K4 . [The feature that CI can be very long
(with small but non-negligible probability) has also been noted in Lee et al. (2016), and is
not shared by the confidence intervals obtained from K1 , K3 , and K4 .] Note finally that we
have obtained the same conclusions in other length simulations, which we do not report for
the sake of brevity.
¯ of Lee et al. (2016) can have conWe now demonstrate that the confidence interval CI
26

Setting
‘Independent’
‘Correlated’

Lengths
Median
90%-quantile
Median
90%-quantile

Confidence interval
¯
CI
K1
K3
K4
0.46 0.78 0.78 0.43
0.51 0.85 0.85 1.06
0.56 0.81 0.81 1.42
0.90 1.30 1.30 14.3

¯ of Lee
Table 2: Medians and empirical quantiles of the lengths of the confidence intervals CI
et al. (2016) and of those obtained from K1 , K3 , and K4 . The nominal coverage probability
is 1 − α = 0.95, n = 100, and p = 10.
ditional coverage probability considerably smaller than the nominal one when λ is selected
by cross-validation. Rather than evaluating the minimal conditional coverage probabilities
conditional on M̂ = M separately for every M satisfying 1 ∈ M , which would be quite costly,
we evaluate the minimal conditional coverage probability where conditioning is on the event
that 1 ∈ M̂ . We denote this quantity by Pcond,min . A simple calculation shows that if we find
that this latter minimal conditional coverage probability is smaller than 1 − α, then it follows
that we must have


(n)
¯ M̂ = M < 1 − α
infp Pn,β,1 x00 [M̂ ]β ∈ CI
β∈R

M̂

for at least some M satisfying 1 ∈ M , showing that property (21) is violated. [To see this, note
that Pn,β,1 ( ·| 1 ∈ M̂ ) is a convex combination (over all M with 1 ∈ M ) of the probabilities
Pn,β,1 ( ·| M̂ = M ), with the (nonnegative) weights summing to 1.]
In order to numerically evaluate the minimal conditional coverage probability Pcond,min we
proceed as follows: We consider eight configurations given by all the possible combinations
of n = 20, 100, p = 2, 10 and α = 0.05, 0.2. Recall that x00 = (1, 0, ..., 0). For each of these
eight configurations, the (transpose of the) rows of X are sampled once from the N (0, Σ)distribution and then remain fixed throughout the minimal coverage probability evaluation.
For p = 2, we take Σ to have 1 as the diagonal and 0.8 as the off-diagonal elements. For
p = 10, we take Σ to coincide with the identity matrix Ip , exceptpthat the zero elements in
the last row and column of Ip are replaced by the constant c = 0.8/(p − 1). For each of
these eight configurations, we carry out a three-step minimal conditional coverage probability
evaluation as described in Appendix F.
The so evaluated minimal conditional coverage probabilities P̂cond,min , say, are presented in
Table 3. For comparison, we also provide similar evaluations of minimal conditional coverage
probabilities using the same procedure as described in Appendix F, but now with λ fixed
as a function of X as in the beginning of Section 7 of Lee et al. (2016). When λ is fixed,
these minimal conditional coverage probabilities are approximately equal to the nominal level
1 − α, in agreement with the results of Lee et al. (2016). However, when λ is selected by
cross-validation, the evaluated minimal conditional coverage probabilities can be way below
the nominal level. In particular, these probabilities can be equal to 0.31 for a nominal level
of 0.80 and to 0.86 for a nominal level of 0.95.
27

p

n
20

2
100
20
10
100

1−α
0.80
0.95
0.80
0.95
0.80
0.95
0.80
0.95

Fixed λ
0.80
0.95
0.80
0.95
0.79
0.94
0.79
0.95

LASSO
CV-selected λ
0.43
0.93
0.31
0.86
0.79
0.93
0.70
0.92

Table 3: Monte-Carlo estimates of the minimal conditional coverage probabilities (w.r.t. β
and σ) of the confidence intervals of Lee et al. (2016), with the LASSO model selector where
λ is either fixed or selected by cross validation.
In addition, for λ selected by cross-validation and in all the configurations of n, p, and
α, for the vector β leading to the minimal conditional coverage probability P̂cond,min , we
(n)
¯ by P̂ (1 ∈
∈ CI)
can also estimate the unconditional coverage probability Pn,β,1 (x00 [M̂ ]β
M̂
¯ in order
M̂ )P̂cond,min + P̂ (1 6∈ M̂ ). [Here we make use of the aforementioned amendment to CI
to allow for a well-defined unconditional coverage probability.] In this estimate, P̂ (1 ∈ M̂ ) is
the proportion of times the first regressor belongs to M̂ , over the Monte-Carlo samples in the
third step (as described in Appendix F). The so evaluated unconditional coverage probabilities
are 0.60 for p = 2, n = 100, 1 − α = 0.80 and 0.91 for p = 2, n = 100, 1 − α = 0.95, which
implies that the confidence intervals of Lee et al. (2016) also have minimal unconditional
coverage probabilities below the nominal level when λ is estimated by cross-validation.
The conclusion of this comparison, and particularly of the evaluations of minimal coverage
probabilities, is that, although the confidence intervals of Lee et al. (2016) are conditionally
valid and convenient to compute, their current applicability appears to be restricted to the
case where the tuning parameter λ is fixed. [An extension of the strategy of Lee et al. (2016)
to cross-validated versions of LASSO has recently been studied in Loftus and Taylor (2015)
and Loftus (2015). This extension, however, comes with significantly higher computational
cost.] This also highlights the benefit of the confidence intervals introduced in Section 2,
which are intrinsically designed to be valid for any model selection procedure whatsoever.

5

Conclusion

We have extended the PoSI confidence intervals of Berk et al. (2013a) to PoSI intervals for
(n)
(?)
and x00 [M̂ ]β , minimize
predictors. The coverage targets of our intervals, i.e., x00 [M̂ ]β
M̂
M̂
a certain in-sample prediction error and, under additional assumptions relating the training
period to the prediction period, a certain out-of-sample prediction error, respectively. For
(n)
in-sample prediction, i.e., for the target x00 [M̂ ]β , our intervals are valid, in finite samples,
M̂
irrespective of the model selection procedure that is being used. For out-of-sample prediction,
28

(?)

i.e., for the target x00 [M̂ ]β , the same is true asymptotically under very mild assumptions
M̂
on the underlying model selector. See also Greenshtein and Ritov (2004) for optimality
results related to the latter target and for its feasible counterpart, under appropriate sparsity
conditions.
Two types of confidence intervals were studied here: The first one (corresponding to the
constant K1 (x0 , M̂ )) depends on all components of the vector x0 (even if only a subset of these
components is ‘active’ in the selected model M̂ ) and thus is feasible only if x0 is observed
completely. The intervals of the second type (corresponding to the constants K2 (x0 [M̂ ], M̂ ),
K3 (x0 [M̂ ], M̂ ), and K4 ) depend only on the active components in the selected model, i.e., on
x0 [M̂ ]. The constants K2 , K3 , and K4 correspond to successively larger confidence intervals.
Computing the constant K2 was found to be quite expensive in practice. For computing
the remaining constants, simple algorithms were presented. The computational complexity
of our algorithms for computing K1 and K3 is governed by the number of candidate models
under consideration, limiting computations to a few million candidate models in practice.
Computation of K4 is easy and not limited by complexity constraints (see, however, the
warning about numerical stability in Remark E.5 in Appendix E). Our algorithms are of
similar computational complexity as those proposed in Berk et al. (2013a).
We furthermore have studied the behavior of the constants Ki and of the corresponding
confidence intervals through analytic results in a setting where model dimension is allowed
to grow with sample size, and also through simulations. These results provide evidence that
K4 , which is relatively cheap to compute, is a reasonably tight bound for the computationally
more expensive constants K1 to K3 . Furthermore, these results show that all the constants
K1 to K4 are ‘bounded away’ from the Scheffé constant.
We have also provided simulation results regarding the coverage probabilities of the various
intervals introduced in the paper. We find that the asymptotic results in Section 3 regarding
the design-independent target already ‘kick-in’ at moderate sample sizes, and these results
demonstrate that the PoSI confidence intervals for the predictors maintain the desired minimal
coverage probability. The simulation study also shows that ‘naive’ confidence intervals, which
ignore the data-driven model selection step and which use standard confidence procedures as
if the selected model were correct and given a priori, are invalid also in the setting considered
here (which is in line with earlier findings in Leeb et al. (2015), where inter alia ‘naive’ con(n)
fidence intervals for components of β
were studied). Furthermore, studying the confidence
M̂
intervals developed for model selection with the LASSO by Lee et al. (2016), and others, we
find that these intervals are invalid if the LASSO penalty is chosen by cross-validation. This
contrasts the established fact that these intervals are valid (conditionally on the event that a
given model is selected), if the penalty is fixed in advance.
Acknowledgement. We thank the referees and an Associate Editor for thoughtful feedback and constructive comments. The first author acknowledges constructive discussions with
Lukas Steinberger and Nina Senitschnig on the topic of the paper. The second author acknowledges partial support from the Austrian Science Fund (FWF) projects P 28233-N32 and
P 26354-N26.

29

A

Appendix: On the assumptions on σ̂ 2

In line with Berk et al. (2013a) we have postulated the existence of an estimator σ̂ 2 that is
independent of PX Y and is distributed as σ 2 /r times a chi-square distributed random variable
with r degrees of freedom (1 ≤ r < ∞). As already noted in Section 2, if we assume that
d < n and µ = Xβ hold, such an estimator always exists and is given by the usual residual
variance estimator obtained from the residuals Y − PX Y . However, if d = n holds (which
typically is the case if p > n) or if µ is not known to belong to the column space of X, such
an estimator is much harder to come by.
Consider first the case where d = n holds. Then it is plain that such an estimator does not
exist if it is to be only a function of Y (and X): Since here PX Y = Y holds, such an estimator
would have to be independent of Y and thus constant with probability one, contradicting the
requirement to be distributed as a positive multiple of a chi-square. In order to nevertheless
be able to come up with an estimator σ̂ 2 with the desired properties, one is hence forced to
assume that one has access to additional data beyond Y that are related to Y in an appropriate
way. A prototypical situation where such a construction is possible is as follows: Assume that
one has available additional data Y ∗ distributed as N (µ∗ , σ 2 In∗ ), independently of Y (for
example, Y ∗ might have been obtained from splitting the original larger sample into Y and
Y ∗ ). Assume further that for Y ∗ one has available a (non-trivial) correct regression model
(i.e., µ∗ = X ∗ δ with X ∗ of full column rank less than n∗ ). Obviously, then an estimator σ̂ 2
satisfying all the required properties can be constructed from this correct regression model for
Y ∗ . However, this raises the question why one would be willing to assume a correct regression
model for one part of the data, but would refuse to do so for the other part. [This might be
defended by reference to a structural break in the mean, which then however would beg the
question why the structural break would not also affect the variance σ 2 .] Alternatively to the
assumption µ∗ = X ∗ δ, one could assume some ‘smoothness’ in µ∗ and then use nonparametric
estimators to produce σ̂ 2 . Again the question arises why one would then not make a similar
assumption for µ and use the nonparametric method also for the first (or the entire) sample.
In the quite special situation where one has replicated observations in Y ∗ available, one can
abandon the dependence on a correct model (or on smoothness assumptions) and nevertheless
produce an estimator σ̂ 2 with the desired properties. All this granted, it seems that the desired
assumptions on σ̂ 2 and the desire to treat the case d = n are not completely at ease.
Second, in case d < n , but it is not assumed that µ = Xβ holds, it is not obvious how
an estimator σ̂ 2 with the desired properties can be constructed without further assumptions
(note that the residual variance estimator obtained from Y − PX Y while being independent
from PX Y will in general not be guaranteed to follow the required distribution). One such
assumption could be that we have available a correct model µ = Zγ, where the column space
of Z contains the column space of X with the rank of Z still less than n; we could then compute
σ̂ 2 from this larger model (i.e., from Y − PZ Y ), the resulting estimator having the desired
properties. While this assumption solves the existence problem for σ̂ 2 , it raises the question
why one would then still want to keep the model selection exercise restricted to submodels
defined by the columns of X, when it is known that the correct, larger, model µ = Zγ holds
(and Z is available). Hence, we are led back essentially to the classical case with Z playing
the role of X. Alternatively, the same constructions as in the preceding paragraph relying on
30

an independent sample Y ∗ are available, but they again suffer from the limitations pointed
out before.
The discussion in this section shows that outside of the framework d < n and µ = Xβ
(in which case p = d can be assumed with little loss of generality) the assumptions on σ̂ 2
made in Berk et al. (2013a) as well as in the present paper are less than innocuous and will be
satisfied only in quite special situations; cf. also the discussion in Remark 2.1(ii) in Leeb et al.
(2015). [For this very reason, the first version of this paper was set in the classical framework.]
For ways to work around this assumption on σ̂ 2 by using “conservative” estimators for the
variance see Bachoc et al. (2016).

B

Appendix: Proofs for Section 2

Lemma B.1. Suppose W is a random m × 1 vector that has a density that is positive almost
everywhere. Let a1 , . . . , aL , for some L ∈ N, be elements of Rm , not all of which are zero.
Define h(w) = maxl=1,...,L |a0l w|, and set H (t) = Pr (h (W ) ≤ t) for t ∈ R. Then H is
continuous on R, satisfies H (t) = 0 for t ≤ 0, and is strictly increasing on [0, ∞).
Proof: For t < 0 the event {h (W ) ≤ t} is empty; for t = 0 this event is an intersection of
the sets {a0l W = 0} where at least one of these sets has probability zero because W possesses
a density and not all al are zero. Consequently, H (t) = 0 for t ≤ 0 follows. Because H
is a distribution function, continuity of H on R will follow if we can establish continuity on
(0, ∞). Now, for every t > 0 the event {h (W ) = t} is contained in the union of the events
{|a0l W | = t} for which al 6= 0 holds. Since any of these events has probability zero, it follows
that Pr ({h (W ) = t}) = 0 and consequently H is continuous on (0, ∞). It remains to establish
the claim regarding strict monotonicity: For t > 0 the set A(t) = {w : h (w) ≤ t} contains
a sufficiently small ball centered at the origin because h (0) = 0 and h is continuous, and
consequently H (t) > 0 follows by the assumption on the density of W . It hence suffices to
show that 0 < t1 < t2 implies H (t1 ) < H (t2 ). Because not all al are zero and h is positively
homogeneous of degree one, we can find an element w1 ∈ A(t1 ) such that h (w1 ) = t1 holds.
But then there exists an l1 such that a0l1 w1 = t1 and |a0l w1 | ≤ t1 for all l hold. In fact, we
may assume that a0l1 w1 = t1 holds (otherwise we change the sign of w1 ). Consider the set B
consisting of all w ∈ Rm such that a0l1 (w − w1 ) > 0 and such that |a0l (w − w1 )| < (t2 − t1 ) /2
for every l. Then B ⊆ A(t2 )\A(t1 ) holds, since for w ∈ B
h (w) ≤ max a0l (w − w1 ) + h (w1 ) < (t2 − t1 ) /2 + t1 = (t1 + t2 ) /2 < t2 ,
l=1,...,L

a0l1 w > a0l1 w1 = t1 > 0,
and hence also h (w) ≥ a0l1 w > t1 hold. But B obviously has positive Lebesgue measure,
implying that H (t2 ) − H (t1 ) = Pr (A(t2 )\A(t1 )) > 0. 
Remark B.2. In the special case where W = W1 W2 with W1 a random m × 1 vector
having a density that is positive almost everywhere, with W2 a random variable that is
independent of W1 , is positive almost surely, and has a density that is almost everywhere
positive on (0, ∞), an alternative, and perhaps simpler, proof is as follows: Set H ∗ (t) =
31

Pr (h (W1 ) ≤ t). We conclude that H ∗ is continuous on R and satisfies H ∗ (t) = 0 for t ≤ 0
by repeating the corresponding arguments in the preceding proof. The same properties for
H (t) = EW2 H ∗ (t/W2 ) then follow immediately. To establish strict monotonicity of H on
[0, ∞) consider 0 ≤ t1 < t2 . It is not difficult to see that we can then find w2 > 0 such that
H ∗ (t1 /w2 ) < H ∗ (t2 /w2 ) holds since otherwise H ∗ would have to be constant on [0, ∞) which
is impossible since H ∗ (0) = 0 and H ∗ is a distribution function. By continuity of H ∗ then
also H ∗ (t1 /w20 ) < H ∗ (t2 /w20 ) must hold for every w20 in a sufficiently small neighborhood of
w2 . Since H ∗ is nondecreasing and since the distribution of W2 puts positive mass on the
aforementioned neighborhood, we can conclude that EW2 H ∗ (t1 /W2 ) < EW2 H ∗ (t2 /W2 ), i.e.,
that H (t1 ) < H (t2 ) holds.
The following lemma will be used in the proof of Proposition 2.3 below.
Lemma B.3. Suppose F ∗ is a distribution function on R that is continuous at zero. Let S
be a random variable that is positive with probability one and has a continuous distribution
function. Then F (t) = ES F ∗ (t/S) is continuous on R.
Proof: Let S ∗ be a random variable which is independent of S and which has distribution
function F ∗ . Then F (t) = ES ES ∗ 1 (S ∗ ≤ t/S) = Pr (SS ∗ ≤ t). Because S ∗ 6= 0 holds almost
surely by the assumption on F ∗ , we have Pr (SS ∗ = t) = ES ∗ ES 1 (S = t/S ∗ ). Since S has a
continuous distribution function, we have ES 1 (S = t/S ∗ ) = 0 almost surely, implying that
Pr (SS ∗ = t) = 0. 
Proof of Proposition 2.3: (a) Observe that in case M = {1, . . . , p} we have c(M, M) =
∗
0 and that x0 = 0 implies s̄M∗ = 0 for every M∗ ; thus it is obvious that FM,x
, and hence
0
also FM,x0 , is the indicator function of [0, ∞), which then implies that K3 exists and is welldefined in this case and that K3 (x0 [M ], M ) = 0. If M = {1, . . . , p} but x0 6= 0, then FM,x0 is
continuous on R, satisfies FM,x0 (t) = 0 for t ≤ 0, and is strictly increasing on [0, ∞) in view
of Lemma B.1, since in this case FM,x0 (t) reduces to Pn,µ,σ maxM∗ ∈M s̄0M∗ (Y − µ) /σ̂ ≤ t
(see (22) and (23) below) and since not all s̄M∗ can be zero (in view of our assumptions on
M). Again the claims then clearly follow in this case.
In case M ∈ M is a proper subset of {1, . . . , p} and p > 1 holds we argue as follows: Note
∗
that then FM,x
(0) = 0 holds since now c(M, M) ≥ 1 holds and since FBeta,1/2,(d−1)/2 (0) = 0.
0
∗
Hence FM,x0 is continuous at t = 0. We may apply Lemma B.3 to conclude that FM,x0
is continuous on R and thus satisfies FM,x0 (0) = 0 (since FM,x0 (t) = 0 for t < 0 by its
∗
definition). Next let 0 ≤ t1 < t2 . Because FM,x
(0) = 0 as noted before and because
0

0
∗
FM,x0 (1) = 1 (since Pr maxM∗ ⊆M s̄M∗ V > 1 = 0 and FBeta,1/2,(d−1)/2 (1) = 1) we thus can
∗
∗
∗
find a positive g0 such that FM,x
(t1 /g0 ) < FM,x
(t2 /g0 ) holds (if not, constancy of FM,x
0
0
0
on [0, ∞) would have to follow). Because of continuity from the right at t1 /g0 it follows that
∗
∗
FM,x
(t1 /g) < FM,x
(t2 /g) also holds for all g < g0 in a sufficiently small neighborhood of
0
0
∗
∗
g0 that is contained in (0, ∞). Because FM,x
(t1 /g) ≤ FM,x
(t2 /g) holds for every g > 0
0
0
and because G has a density that is positive everywhere on (0, ∞), the strict inequality
∗
∗
EG FM,x
(t1 /G) < EG FM,x
(t2 /G) follows. This establishes strict monotonicity of FM,x0 on
0
0
[0, ∞) also in this case, which proves the claims.
Finally, if M ∈ M is a proper subset of {1, . . . , p} and p = 1 holds, then M is empty
∗
and d = 1 must hold, and hence FM,x
reduces to the indicator function of [1, ∞). But then
0
32

FM,x0 (t) = Pr (G ≤ t) which obviously is continuous on R, takes the value zero at t = 0, and
is strictly increasing on [0, ∞), again implying the claims.
(b) Observe that s̄M∗ belongs to the column space of X for every M∗ ∈ M and hence we
have


0
Pn,µ,σ max s̄M∗ (Y − µ) /σ̂ > t
M∗ ∈M


0
(22)
= Pn,µ,σ max s̄M∗ PX (Y − µ) / kPX (Y − µ)k > (σ̂/ kPX (Y − µ)k) t ,
M∗ ∈M

where PX (Y − µ) / kPX (Y − µ)k and kPX (Y − µ)k /σ̂ are independent since the random
variables PX (Y − µ) / kPX (Y − µ)k, kPX (Y − µ)k, and σ̂ are mutually independent. [Observe that PX (Y − µ) is nonzero with probability 1 since d ≥ 1 holds.] Consequently, the
probability given above can be represented as


0
(23)
Pr max s̄M∗ V > t/G
M∗ ∈M

∗
where V and G are independent and otherwise are as in the definition of FM,x
and FM,x0 .
0
Now, using first independence of V and G and then a union bound twice we have for M ∈ M
and t ≥ 0

 Z


0
0
Pr max s̄M∗ V > t/G = Pr max s̄M∗ V > t/g dFG (g)
M∗ ∈M
M∗ ∈M





Z
0
0
s̄M∗ V > t/g dFG (g)
max
s̄M∗ V > t/g + Pr
≤
min 1, Pr
max
M∗ ∈M,M∗ ⊆M
M∗ ∈M,M∗ "M




Z
X

≤
min 1, Pr
max
s̄0M∗ V > t/g +
Pr s̄0M∗ V > t/g  dFG (g)
M∗ ∈M,M∗ ⊆M


Z
=

min 1, Pr
Z

≤



M∗ ∈M,M∗ "M



max

M∗ ∈M,M∗ ⊆M

s̄0M∗ V > t/g +


X

Pr



s̄0M∗ V

2



> t2 /g 2  dFG (g)

M∗ ∈M,M∗ "M


∗
∗
1 − FM,x
(t/g) dFG (g) = EG 1 − FM,x
(t/G) = 1 − FM,x0 (t) ,
0
0


(24)

where
 FG here denotes
 the c.d.f. of G. The last inequality follows from the fact that 
2
0
2
2
Pr s̄M∗ V > t /g is either equal to zero (if s̄M∗ = 0) or is equal to 1−FBeta,1/2,(d−1)/2 t2 /g 2
(if s̄M∗ 6= 0) as is easy to see; for
 the case where M is the empty set also observe that
0
Pr maxM∗ ∈M,M∗ ⊆M s̄M∗ V > t/g = 0 for t ≥ 0 because s̄∅ = 0. In view of (10) the
chain of inequalities in (22)-(24) establishes K1 (x0 ) ≤ K3 (x0 [M ], M ). It follows that K1 (x) ≤
K3 (x[M ], M ) = K3 (x0 [M ], M ) for every x satisfying x[M ] = x0 [M ], implying K2 (x0 [M ], M ) ≤
K3 (x0 [M ], M ). The inequality (16) is obvious and inequality (17) follows
 since for t ≥ 0 we
have (again noting that expressions like Pr maxM∗ ∈M,M∗ ⊆M1 s̄0M∗ V > t for t ≥ 0 are equal

33

to zero if M1 is empty)



0
s̄M∗ V > t + c(M2 , M) 1 − FBeta,1/2,(d−1)/2 t2
Pr
max
M∗ ∈M,M∗ ⊆M2




0
0
max
s̄M∗ V > t
s̄M∗ V > t + Pr
≤ Pr
max
M∗ ∈M,M∗ ⊆M1
M∗ ∈M,M∗ ⊆M2 ,M∗ "M1

+c(M2 , M) 1 − FBeta,1/2,(d−1)/2 t2


X

Pr s̄0M∗ V > t
s̄0M∗ V > t +
≤ Pr
max
M∗ ∈M,M∗ ⊆M1

M∗ ∈M,M∗ ⊆M2 ,M∗ "M1


+c(M2 , M) 1 − FBeta,1/2,(d−1)/2 t2



0
s̄M∗ V > t + c(M1 , M) 1 − FBeta,1/2,(d−1)/2 t2 .
≤ Pr
max
M∗ ∈M,M∗ ⊆M1

∗
The relation K3 (x0 [M ], M ) ≤ K4 is now immediate. Finally, 1 − F∅,x
(t) ≤ 1 for all t ∈ R
0
∗
and 1 − F∅,x0 (t) = 0 for t > 1 lead to



∗
∗
1 − F∅,x0 (t) = EG 1 − F∅,x
(t/G) = EG 1 − F∅,x
(t/G) 1 (t ≤ G)
0
0

≤ EG 1 (t ≤ G) = 1 − Pr (G ≤ t) ,
which proves K4 ≤ K5 . 
Lemma B.4. Assume p = 2 and n ≥ 2. Then there exists a design matrix X with full column
rank and a vector x0 such that K4 = K1 (x0 ) for M the power set of {1, 2}.
Proof: Assume first that n = 2. In view of the definition of K4 = K3 (x0 [∅], ∅) it suffices
to exhibit a 2 × 2 matrix X and a 2 × 1 vector x0 such that equality holds between the far
l.h.s. and the far r.h.s. of (24) for M = ∅ and all t ≥ 0. Inspection of (24) shows that for
this it suffices to find X and x0 such that




X

Pr
max
s̄0M∗ V > c = min 1,
Pr s̄0M∗ V > c 
∅6=M∗ ⊆{1,2}

∅6=M∗ ⊆{1,2}

holds for every c ≥ 0 and that s̄0M∗ 6= 0 for every ∅ 6= M∗ ⊆ {1, 2}. This is achieved for


1 cos (2π/3)
(2)
X=X =
0 sin (2π/3)
(2)0

and x00 = x0 = (cos (4π/3) , sin (4π/3)) X (2) : Then s̄0{1} = − (1, 0), s̄0{2} = − (cos (2π/3) , sin (2π/3)),

and s̄0{1,2} = (cos (4π/3) , sin (4π/3)). Consequently, the event max∅6=M∗ ⊆{1,2} s̄0M∗ V > c
is either the entire space or is the disjoint union of the events { s̄0{1} V > c}, { s̄0{2} V > c}
and { s̄0{1,2} V > c}. In the case n > 2 simply set

0
X = X (2)0 , 0, . . . , 0
and x00 = (cos (4π/3) , sin (4π/3) , 0, . . . , 0) X. 
34

Remark B.5. Further examples of pairs X, x0 satisfying the above lemma can be generated
from the matrix X constructed in the proof by premultiplying X by an orthogonal matrix
and leaving x0 unchanged.
Proof of Proposition 2.9: (a) The distribution of ω = maxM ∈M |s̄0M (Y − µ)| /σ clearly
does not change if X is replaced by AX, where A is an orthogonal n × n matrix. Furthermore,
scaling the columns of X and the corresponding columns of x00 by the same (column-specific)
positive constants does not alter maxM ∈M |s̄0M (Y − µ)| /σ. Hence, we may assume w.l.o.g.
that X consist of the first p standard basis vectors of Rn . Then choose x00 as the 1 × p vector
(1, . . . , 1). It follows that ω can be written as
max

M ∈M

X

Zi /

p
|M |

i∈M

where Zi are i.i.d. standard normal and where we use the convention that the expression in the display is zero if |M | = 0. For a positive real b define the random set M̃ =
{i ∈ {1, ..., p} : Zi ≥ b}. Since for any realization of the random variables Zi we have that
M̃ ∈ M, we must have
v
r
u
p
p
X
X
X
u
√
−1/2
−1
ω/ p ≥ p
Zi / M̃ = p
Zi 1(Zi ≥ b) /tp−1
1(Zi ≥ b).
i=1

i∈M̃

i=1

p
By the law of large numbers we obtain that the r.h.s. converges to φ(b)/ 1 − Φ(b) almost
√
√
surely. Because K1 (x0 , ∞)/ p is the (1 − α)-quantile of ω/ p with α independent of p and
since b > 0 was arbitrary, the first claim follows. The second claim follows immediately by
choosing x0 equal to a p × 1 standard basis vector and by noting that then ω is distributed as
the absolute value of a standard normal variable. [More generally, if supp kx0 k0 < ∞ holds,
then ω is a maximum absolute value of at most k standard normal variables, where k is a
finite number not depending on p.]
(b) For the same reasons as given at the beginning of the proof of part (a) we have for
every p that
√
inf p inf
inf
K2 (x0 [M ], M, ∞, X, α, M)/ p
x0 ∈R X∈X(p) M ∈M,|M |≤γp

√
K2 (x0 [M ], M, ∞, (Ip , 0p×(n−p) )0 , α, M)/ p
x0 ∈R n≥p M ∈M,|M |≤γp
√
= inf p
inf
K2 (x0 [M ], M, ∞, Ip , α, M)/ p
=

inf p inf

inf

(25)

x0 ∈R M ∈M,|M |≤γp

where Ip is the identity matrix of dimension p. By the monotonicity property (16) the far
r.h.s. of (25) equals
√
(26)
inf p
inf
K2 (x0 [M ], M, ∞, Ip , α, M)/ p.
x0 ∈R M ∈M,|M |=bγpc

Now fix an arbitrary x0 ∈ Rp and M ∈ M with |M | = bγpc. Define x∗0 via x∗0i = x0i for i ∈ M
and set x∗0i = 1 else. Then
K2 (x0 [M ], M, ∞, Ip , α, M) ≥ K1 (x∗0 , ∞, Ip , α, M)
35

(27)

holds and the latter quantity is the (1 − α)-quantile of
sX
X
∗
∗
ω = max
x0i Zi /
(x∗0i )2 .
∗
M ∈M

i∈M ∗

i∈M ∗

For a positive real b define now the random set M̆ = {i ∈
/ M : Zi ≥ b}. Similar as above we
then conclude that
s
X
X
−1
∗ √
Zi 1(Zi ≥ b) / p−1
1(Zi ≥ b).
ω / p≥ p
i∈M
/

i∈M
/

While the r.h.s. of the above display depends on M , its distribution does not as it coincides
with the distribution of
v
u
p
p
X
u −1 X
−1
t
Ap = p
Zi 1(Zi ≥ b) / p
1(Zi ≥ b).
i=bγpc+1

i=bγpc+1

√
√
Consequently, K1 (x∗0 , ∞, Ip , α, M)/ p as the (1 − α)-quantile of ω ∗ / p is not smaller than
the corresponding quantile of Ap , and this is true independently
of the choice of x0 and of M
p
√
with |M | = bγpc. Since Ap converges to 1 − γφ(b)/ 1 − Φ(b) almost surely and α does not
depend on p, we can conclude that
p
p
√
lim inf inf p
inf
K1 (x∗0 , ∞, Ip , α, M)/ p ≥ 1 − γφ(b)/ 1 − Φ(b).
p→∞ x0 ∈R M ∈M,|M |=bγpc

Since b > 0 was arbitrary, the proof is then complete in view of (25), (26), and (27). 
Lemma B.6. Let α, 0 < α < 1, be a fixed number, let N ∈ N, and d ∈ N. Let K(N, d, α)
denote the (1 − α)-quantile of the distribution function given by


1 − EG min 1, N 1 − FBeta,1/2,(d−1)/2 t2 /G2
for t ≥ 0 and by 0 for t < 0. Here G is a nonnegative random variable such that G2 follows
a chi-square distribution with d degrees of freedom. Then as min(N, d) → ∞
q

K(N, d, α)/ d 1 − N −2/(d−1) → 1.
Proof: The c.d.f. in the lemma is the c.d.f. of GW , where W is independent of G, is
nonnegative, and has distribution function given by


1 − min 1, N 1 − FBeta,1/2,(d−1)/2 t2
√
as d → ∞. To complete the
for t ≥ 0. Observe that G/ d converges
q to 1 in probability

proof it thus suffices to show that W/ 1 − N −2/(d−1) converges to 1 in probability as
min(N, d) → ∞: For t > 1 we have that
 q


 


.
Pr W/ 1 − N −2/(d−1) > t ≤ N 1 − FBeta,1/2,(d−1)/2 t2 1 − N −2/(d−1)
36

But the r.h.s. of the preceding display has been shown in Zhang (2013) to converge to zero
as min(N, d) → ∞, cf. (A.5) and (A.6) in that paper. For t < 1 we have
 q

h

 
i

−2/(d−1)
Pr W/ 1 − N
≤ t = 1−min 1, N 1 − FBeta,1/2,(d−1)/2, t2 1 − N −2/(d−1)
,
and hence it suffices to show that

 

N 1 − FBeta,1/2,(d−1)/2, t2 1 − N −2/(d−1)
→∞
as min(N, d) → ∞. But this has again been established in Zhang (2013), see (C.4) and (C.5)
in that paper. 
Proof of Proposition 2.10: Observe that K4 (∞) is always positive and hence the ratio
in (18) is well-defined. In view of (12) and (13) and the assumptions on M we have for
M ∈ M with M 6= {1, . . . , p}, for X ∈ Xn(p),p (M), and for x0 ∈ Rp that K3 (x0 [M ], M, ∞)
is not less than K(bτ |M|c , min(n(p), p), α) in the notation of Lemma B.6, where we note
that bτ |M|c ≥ 1 holds at least from a certain p onwards. Since K3 (x0 [M ], M, ∞) ≤ K4 (∞)
always holds, and since K4 (∞) = K(|M| − 1, min(n(p), p), α) in the notation of Lemma B.6,
it suffices to show that
K(bτ |M|c , min(n(p), p), α)/K(|M| − 1, min(n(p), p), α) → 1
as p → ∞. Note that d = min(n(p), p) → ∞ as p → ∞ by the assumption on n(p), and that
|M| → ∞ as p → ∞. By Lemma B.6 we thus need to show that
 


Ad (|M|) := 1 − (|M| − 1)−2/(d−1) / 1 − bτ |M|c−2/(d−1) → 1
as p → ∞. Observe that τ < 1 must hold, and thus |M| − 1 ≥ bτ |M|c > 1 holds for large p.
This, in particular, implies Ad (|M|) ≥ 1 for large p. But then for large p

 

1 ≤ Ad (|M|) ≤ 1 − |M|−2/(d−1) / 1 − ((τ /2) |M|)−2/(d−1) =: Bd (|M|)
holds since also bτ |M|c ≥ (τ /2) |M| > 1 is satisfied for large p. It thus suffices to show that
Bd (|M|) → 1 for p → ∞. Let f > 2/τ be a real number. Then |M| ≥ f holds for large
p. Because Bd (x) is monotone decreasing in x for x > 2/τ > 1 and for every d as is easily
checked by inspection of the derivative, we have that Bd (f ) ≥ Bd (|M|) holds for large p. But
now it is easily checked (Hôpital’s rule) that Bd (f ) converges to log f / (log f + log(τ /2)) as
p → ∞ (and thus d → ∞). Making f arbitrarily large, log f / (log f + log(τ /2)) approaches
1. This completes the proof of (18). The second claim follows immediately from K4 =
K(|M| − 1, min(n(p), p), α), the preceding lemma, the observation that Bd (|M|) → 1, and
that

 

1/Bd (|M|) ≤ 1 − (|M| − 1)−2/(d−1) / 1 − |M|−2/(d−1) ≤ 1
holds. 

37

Proof of Corollary 2.12: Properties (i), (ii), and (iv) are obvious. In case mp = p,
we have that M(mp ) is the power set and hence |M(mp )| = 2mp = 2p . But then we have
|M | < mp = p for M 6= {1, . . . , p}. Consequently,
c (M, M(mp )) = |M(mp )| − 2|M | ≥ 2mp − 2mp −1 = (1/2) |M(mp )| ≥ (1/3) |M(mp )| .
Next consider the case where mp < p. Then certainly |M(mp )| ≥ 2mp +1 − 1 holds. Now, for
M ∈ M(mp ) we have |M | ≤ mp (and M 6= {1, . . . , p}) and thus
c (M, M(mp )) = |M(mp )| − 2|M | ≥ |M(mp )| − 2mp = |M(mp )| (1 − 2mp / |M(mp )|)

≥ |M(mp )| 1 − 2mp /(2mp +1 − 1) ≥ (1/3) |M(mp )| .
Thus (iii) with
since
 1/3 holds. The next claim now follows from Proposition 2.10
Pmpτ =
p
p
|M(mp )| = k=0 k . The final claim is then a trivial consequence, since |M(p)| = 2 in this
case. Note that Xn(p),p (M(p)) 6= ∅ by the assumptions on n(p), implying n(p) ≥ p. 

C

Appendix: Proofs for Section 3

In the subsequent lemma we assume that σ̃ 1 and σ̃ 2 are defined on the same probability space
as are Y , X, and σ̂ 2 . In slight abuse of notation, we shall then denote by Pn,β,σ the joint
distribution of Y , X, σ̂ 2 , σ̃ 1 , and σ̃ 2 . We note that an argument corresponding to a special
case of this lemma has been used in Ewald (2012).
Lemma C.1. Suppose that the maintained model assumptions of Section 3 are satisfied.
Assume further that Conditions 3.3 and 3.4 hold. Let W be the set of all measurable nonnegative functions of the form W (x0 , X, M ). Then, for any two sequences of random variables
σ̃ 1 = σ̃ 1,n and σ̃ 2 = σ̃ 2,n (which may be functions of σ) satisfying
sup
β∈Rp ,σ>0

Pn,β,σ ( |(σ̃ i /σ) − 1| > δ| X) → 0

(28)

in probability as n → ∞ for every δ > 0 and for i = 1, 2, we have that


(?)
0
0
sup
Pn,β,σ x0 [M̂ ]β̂ M̂ − x0 [M̂ ]β
≤ W (x0 , X, M̂ )σ̃ 1 X
M̂

x0 ∈Rp ,β∈Rp ,σ>0,W ∈W

− Pn,β,σ



x00 [M̂ ]β̂ M̂ − x00 [M̂ ]β

(n)
M̂

≤ W (x0 , X, M̂ )σ̃ 2 X



converges to 0 in probability as n → ∞.
Proof: Because the number of variables p is fixed, it suffices to show for arbitrary but
fixed M ⊆ {1, . . . , p} that


(?)
Qn =
sup
Pn,β,σ x00 [M ]β̂ M − x00 [M ]β M ≤ WM σ̃ 1 ; M̂ = M X
x0 ∈Rp ,β∈Rp ,σ>0,W ∈W

− Pn,β,σ



(n)

x00 [M ]β̂ M − x00 [M ]β M ≤ WM σ̃ 2 ; M̂ = M X
38



goes to 0 in probability, where we have used the abbreviation WM = W (x0 , X, M ). We may
assume in what follows that M 6= ∅ since otherwise Qn is zero. Furthermore, Qn does not
change its value if the supremum is restricted to those x0 which have kx0 [M ]k = 1 (since the
expression inside the supremum is identically zero if x0 satisfies x0 [M ] = 0 and since otherwise
the norm of x0 [M ] can be absorbed into WM ). Hence we have


Qn =
sup
Pn,β,σ |σe1 + σe2 | ≤ WM σ̃ 1 ; M̂ = M X
x0 ∈S(M ),β∈Rp ,σ>0,W ∈W



− Pn,β,σ |σe2 | ≤ WM σ̃ 2 ; M̂ = M X

(29)

where we have used the abbreviations S(M ) = {x0 ∈ Rp : kx0 [M ]k = 1},


−1
e1 = σ −1 n1/2 x00 [M ] X[M ]0 X[M ]
X[M ]0 X[M c ] − (Σ[M, M ])−1 Σ[M, M c ] β[M c ]
and

−1
e2 = σ −1 n1/2 x00 [M ] X[M ]0 X[M ]
X[M ]0 (Y − Xβ) .

Note that we have also absorbed a factor n1/2 into WM , which is possible because of the supremum operation w.r.t. WM . Using the inequality |Pr (A ∩ C) − Pr (B ∩ C)| ≤ Pr (Ac ∩ B ∩ C)+
Pr (A ∩ B c ∩ C) we can bound the absolute value inside the supremum in (29) by


Pn,β,σ |σe1 + σe2 | > WM σ̃ 1 ; |σe2 | ≤ WM σ̃ 2 ; M̂ = M X


(30)
+Pn,β,σ |σe1 + σe2 | ≤ WM σ̃ 1 ; |σe2 | > WM σ̃ 2 ; M̂ = M X .
Let now δ n,1 be an arbitrary sequence of positive numbers converging to zero. Then we can
further bound the above expression by

Pn,β,σ σ̃ 1 (WM − δ n,1 )+ ≤ |σe2 | ≤ σ̃ 2 WM X
+Pn,β,σ ( σ̃ 2 WM ≤ |σe2 | ≤ σ̃ 1 (WM + δ n,1 )| X)


+2Pn,β,σ |σe1 | ≥ σ̃ 1 δ n,1 ; M̂ = M X .

(31)

By the assumption on the estimators σ̃ 1 and σ̃ 2 we can find a sequence δ n,2 < 1 of positive
numbers converging to zero such that


sup Pn,β,σ max |(σ̃ i /σ) − 1| > δ n,2 X → 0
(32)
β,σ

i=1,2

in probability as n → ∞. This can easily be seen from a diagonal sequence argument. Now,

39

using (29), (30), (31), and (32), we have
Qn ≤

sup
x0 ∈S(M ),β∈Rp ,σ>0,W ∈W

+

Pn,β,σ (1 − δ n,2 ) (WM − δ n,1 )+ ≤ |e2 | ≤ (1 + δ n,2 ) WM X



Pn,β,σ ( (1 − δ n,2 ) WM ≤ |e2 | ≤ (1 + δ n,2 ) (WM + δ n,1 )| X)

sup
x0 ∈S(M ),β∈Rp ,σ>0,W ∈W

+2



Pn,β,σ |e1 | ≥ (1 − δ n,2 ) δ n,1 ; M̂ = M X + op (1)

sup
x0 ∈S(M ),β∈Rp ,σ>0,W ∈W

≤ 2

sup
x0 ∈S(M ),β∈Rp ,σ>0,W ∈W

+2

Pn,β,σ (1 − δ n,2 ) (WM − δ n,1 )+ ≤ |e2 | ≤ (1 + δ n,2 ) (WM + δ n,1 ) X



Pn,β,σ |e1 | ≥ (1 − δ n,2 ) δ n,1 ; M̂ = M X + op (1)

sup
x0 ∈S(M ),β∈Rp ,σ>0

= 2Qn,1 + 2Qn,2 + op (1).
We first bound Qn,1 as follows: Observe that, conditionally on X, the quantity e2 is normally
−1
distributed with mean zero and variance given by cn (x0 , X) = x00 [M ] n−1 X[M ]0 X[M ]
x0 [M ].
−1
0
By Condition 3.4 the variance cn (x0 , X) converges to c(x0 ) = x0 [M ] (Σ[M, M ]) x0 [M ] > 0
in probability, and in fact even uniformly in x0 ∈ S(M ). Since Σ[M, M ] is obviously positive
definite, 0 < c∗ ≤ c(x0 ) ≤ c∗ < ∞ must hold for all x0 ∈ S(M ). Consequently,


sup
cn1/2 (x0 , X) /c1/2 (x0 ) − 1
(33)
x0 ∈S(M )

converges to zero in probability. Therefore we can find a sequence δ n,3 ∈ (0, 1) converging to
zero for n → ∞ such that the event Dn where (33) is less than δ n,3 has probability converging
to 1. On this event inf x0 ∈S(M ) cn (x0 , X) is then positive for sufficiently large n and we have

40



on Dn and for sufficiently large n
n 

Qn,1 = 2
sup
Φ (1 + δ n,2 ) (WM + δ n,1 ) /cn1/2 (x0 , X)
x0 ∈S(M ),W ∈W

≤

≤

≤

≤


o
−Φ (1 − δ n,2 ) (WM − δ n,1 )+ /cn1/2 (x0 , X)

 
1 + δ n,2
1/2
(WM + δ n,1 ) /c (x0 )
2
sup
Φ
1 − δ n,3
x0 ∈S(M ),W ∈W


1 − δ n,2
1/2
−Φ
(WM − δ n,1 )+ /c (x0 )
1 + δ n,3
 

1 + δ n,2
1/2
2
sup
Φ
(WM + δ n,1 ) /c (x0 )
1 − δ n,3
x0 ∈S(M ),W ∈W


1 − δ n,2
1/2
−Φ
(WM − δ n,1 ) /c (x0 )
1 + δ n,3
 

1 + δ n,2 
z + δ n,1 /c1/2 (x0 )
2
sup
Φ
1 − δ n,3
x0 ∈S(M ),z≥0


1 − δ n,2 
1/2
−Φ
z − δ n,1 /c (x0 )
1 + δ n,3
 



1 + δ n,2 
1 − δ n,2 
1/2
1/2
2 sup Φ
z + δ n,1 /c∗
−Φ
z − δ n,1 /c∗
,
1 − δ n,3
1 + δ n,3
z≥0

where Φ denotes the standard normal c.d.f. But the far right-hand side in the above display
obviously converges to zero for n → ∞ since δ n,1 , δ n,2 , as well as δ n,3 converge to zero. We
have thus established that Qn,1 converges to zero in probability as n → ∞.
We next turn to Qn,2 . In case M = {1, . . . , p}, we have that e1 = 0, and hence Qn,2 = 0.
Otherwise, from Condition 3.3 we can conclude (from a diagonal sequence argument) the
existence of a sequence of positive numbers δ n,4 that converge to zero for n → ∞ such that
n
o
sup Pn,β,σ (M̂ = M |X) : β ∈ Rp , σ > 0, kβ[M c ]k /σ ≥ δ n,4 → 0
in probability as n → ∞. Then
Qn,2 ≤



Pn,β,σ |e1 | ≥ (1 − δ n,2 ) δ n,1 ; M̂ = M X

sup
x0 ∈S(M ),kβ[M c ]k/σ≥δ n,4

+



Pn,β,σ |e1 | ≥ (1 − δ n,2 ) δ n,1 ; M̂ = M X

sup
x0 ∈S(M ),kβ[M c ]k/σ<δ n,4

≤



Pn,β,σ M̂ = M X

sup
kβ[M c ]k/σ≥δ n,4

+

Pn,β,σ ( |e1 | ≥ (1 − δ n,2 ) δ n,1 | X)

sup
x0 ∈S(M ),kβ[M c ]k/σ<δ n,4

≤ op (1) +

Pn,β,σ ( |e1 | ≥ (1 − δ n,2 ) δ n,1 | X) .

sup
x0 ∈S(M ),kβ[M c ]k/σ<δ n,4

41

(34)

Using the Cauchy-Schwartz inequality we obtain for x0 ∈ S(M )


−1
|e1 | ≤ x00 [M ] kβ[M c ]/σk n1/2 X[M ]0 X[M ]
X[M ]0 X[M c ] − (Σ[M, M ])−1 Σ[M, M c ]
≤ kβ[M c ]/σk Bn (X)
where Bn (X) ≥ 0 is Op (1), this following from Condition 3.4 and positive definiteness of
Σ[M, M ]. This shows that the second term on the far right-hand side of (34) is bounded by
1 (δ n,4 Bn (X) ≥ (1 − δ n,2 ) δ n,1 ) .
1/2

If we set now, for example, δ n,1 = δ n,4 , we see that the above quantity converges to zero in
probability as n → ∞, implying that Qn,2 converges to zero in probability as n → ∞. This
completes the proof. 
Proof of Theorem 3.6: (a) Use Lemma C.1 with W (x0 , X, M ) equal to K1 (x0 , r) ksM k
(K2 (x0 [M ], M, r) ksM k, K3 (x0 [M ], M, r) ksM k, K4 (r) ksM k, or K5 (r) ksM k, respectively) and
σ̃ 1 = σ̃ 2 = σ̂ and combine this with Proposition 2.1 (Corollaries 2.2, 2.4, respectively). Note
that r = rn → ∞ because of Condition 3.5, and hence σ̂ satisfies (28).
(b) Let σ̃ 2 be a sequence of random variables such that, conditionally on X, σ̃ 22 is independent of β̂ and is distributed as σ 2 /r∗ times a chi-squared distributed random variable with r∗ degrees of freedom with the convention that σ̃ 2 = σ in case r∗ = ∞. [Such
a sequence exists: Possibly after redefining the relevant random variables on a sufficiently
rich probability space we may find a sequence (Zi )i∈N of i.i.d. standard Gaussian random
P∗
variables that is independent of Y and X. Then define σ̃ 22 = σ 2 ri=1 Zi2 /r∗ if r∗ < ∞
and set σ̃ 22 = σ 2 otherwise.] In view of Remarks 2.5 and 2.6(iii) we have that Proposition
2.1 (Corollaries 2.2, 2.4, respectively) also hold if the confidence interval (5) for the target
(n)
x00 [M̂ ]β uses σ̃ 2 instead of σ̂ (but uses M̂ as in part (a)) and uses the constants K1 (x0 , r∗ )
M̂

(K2 (x0 [M̂ ], M̂ , r∗ ), K3 (x0 [M̂ ], M̂ , r∗ ), K4 (r∗ ), or K5 (r∗ ), respectively). Now apply Lemma
C.1 with W (x0 , X, M ) equal to K1 (x0 , r∗ ) ksM k (K2 (x0 [M ], M, r∗ ) ksM k, K3 (x0 [M ], M, r∗ ) ksM k,
K4 (r∗ ) ksM k, or K5 (r∗ ) ksM k, respectively) and with σ̃ 1 = σ̃. Note that σ̃ 1 satisfies (28) by
assumption, while σ̃ 2 satisfies it because r∗ → ∞ has been assumed. 
Lemma C.2. Suppose that the maintained model assumptions of Section 3 are satisfied and
that X 0 X/n → Σ in probability for n → ∞. Assume further that Condition 3.3 holds and
define σ̂ 2M̂ = ||Y − X[M̂ ]β̂ M̂ ||2 /(n − |M̂ |) for n > p. Then σ̂ 2M̂ satisfies condition (20).
Proof: Clearly
Y − X[M̂ ]β̂ M̂ = Y − PX[M̂ ] Y = PX[M̂ ]⊥ U + PX[M̂ ]⊥ X[M̂ c ]β[M̂ c ] = A + B,
where PX[M̂ ]⊥ denotes orthogonal projection on the orthogonal complement of the column
space of X[M̂ ]. By the triangle inequality we hence have

−1/2

−1/2

kB/σk
σ̂ M̂ /σ − 1 ≤
n − |M̂ |
kA/σk − 1 + n − |M̂ |
≤



−1/2

−1/2
n − |M̂ |
kA/σk − 1 + n − |M̂ |
X[M̂ c ]β[M̂ c ]/σ .
42

We now bound the probability in (20) by the sum of the probabilities that the first and second
term on the far r.h.s. of the preceding display, respectively, exceed δ/2. Because p is fixed
there is a fixed finite number of possible models M̂ and thus for δ > 0 we have the bound for
the first term
 

−1/2
sup Pn,β,σ
n − |M̂ |
kA/σk − 1 ≥ δ/2 X
β∈Rp ,σ>0

=

X

sup

β∈Rp ,σ>0 M

≤

X

sup

p
M β∈R ,σ>0

Pn,β,σ



(n − |M |)−1/2 PX[M ]⊥ U/σ − 1 ≥ δ/2, M̂ = M X

Pn,β,σ




(n − |M |)−1/2 PX[M ]⊥ U/σ − 1 ≥ δ/2 X .



Note that the probabilities in the upper bound on the far r.h.s. of the preceding display do
actually neither depend on β nor σ and are each of the form Pr (|W/w − 1| ≥ δ) where W is
distributed as the square root of a chi-squared random variable with w2 degrees of freedom.
Since w2 = n − |M | goes to infinity for n → ∞ and any fixed M , and since the sum has
a fixed finite number of terms, we can conclude that the upper bound converges to zero in
probability as n → ∞.
Turning to the second term we have, letting λmax denote the largest eigenvalue of a
symmetric matrix,


−1/2
c
c
sup Pn,β,σ
n − |M̂ |
X[M̂ ]β[M̂ ]/σ ≥ δ/2 X
β∈Rp ,σ>0

≤

sup
β∈Rp ,σ>0

≤

sup
β∈Rp ,σ>0

≤

Pn,β,σ




−1/2


1/2
c 0
c
c
n − |M̂ |
λmax X[M̂ ] X[M̂ ] β[M̂ ]/σ ≥ δ/2 X




0
c
X
X/
(n
−
p)
β[
M̂
]/σ
≥
δ/2
X
Pn,β,σ λ1/2
max

X
M 6={1,...,p}

sup
β∈Rp ,σ>0




1/2
Pn,β,σ M̂ = M, λmax
X 0 X/ (n − p) kβ[M c ]/σk ≥ δ/2 X .

Now, since X 0 X/ (n − p) converges to the positive definite matrix Σ in probability, we can
find an event Dn , which has probability converging to 1 for n → ∞, such that on this event
λmax (X 0 X/ (n − p)) is not larger than 4λmax (Σ). Hence, on Dn we can bound each supremum
on the far r.h.s. of the preceding display by


−1/2
sup Pn,β,σ M̂ = M, kβ [M c ] /σk ≥ λmax
(Σ) δ/4 X
β∈Rp ,σ>0

n
o
= sup Pn,β,σ (M̂ = M |X) : β ∈ Rp , σ > 0, kβ[M c ]k /σ ≥ λ−1/2
(Σ)
δ/4
,
max
which goes to zero in probability as n → ∞ by Condition 3.3. 

43

D
D.1

Appendix: Comments on and extension of results in Section
3
Measurability issues

Various statements concerning uncountable suprema (infima) of conditional probabilities occur in Section 3, Appendix C, and in Appendix D.3 below, such as statements that these
quantities converge in probability. It is not difficult to see that – in absence of measurability
– all these statements remain valid if they are properly interpreted as statements referring to
outer probability. This thus relieves one from the need to establish measurability. For this
reason we do not explicitly mention the measurability issues in the presentation of the results
in Section 3 as well as Appendices C and D.3.

D.2

Some remarks on Theorem 3.6

Remark D.1. Under the assumptions of Theorem 3.6(b) we further have that


(n)
0
∗
x
[
M̂
]β
inf
P
∈
CI
(x
)
X
≥ (1 − α) + op (1),
0
n,β,σ
0
p
p
M̂

x0 ∈R ,β∈R ,σ>0

holds, where the op (1) term above depends only on X and converges to zero in probability
as n → ∞. This follows easily from a repeated application of Lemma C.1 in Appendix
C. [Regarding Theorem 3.6(a) recall that the finite-sample coverage result for the target
(n)
x00 [M̂ ]β in Section 2 continues to hold in the context of Section 3 if interpreted conditionally
M̂
on X.]
Remark D.2. (Random x0 ) If x0 is random and independent of X, U , and σ̂ 2 , Theorem
3.6 continues to hold if the result is then being interpreted as conditional on X and x0 . A
particular consequence of this result conditional on X and x0 is then that the confidence
interval CI(x0 ) also satisfies


(?)
0
x
[
M̂
]β
∈
CI(x
)
X
≥ (1 − α) + op (1)
inf
P
0
n,β,σ
0
p
β∈R ,σ>0

M̂

where again op (1) is a function of X only (and Pn,β,σ here represents the distribution of Y ,
X, σ̂ 2 , and x0 ). As noted at the beginning of Section 3, the results in Section 2 continue to
hold if interpreted conditionally on X and x0 . As a consequence, we thus also have that


(n)
0
inf
P
x
[
M̂
]β
∈
CI(x
)
X
≥1−α
0
n,β,σ
0
p
β∈R ,σ>0

M̂

holds. See also Leeb (2009), where prediction intervals for y0 are studied in a similar setting.
Remark D.3. (Relaxing the assumptions on X) The assumption that the rows of X follow
a common distribution L has been used only to define the matrix Σ, which in turn is used in
(?)
the definition of β M . If this assumption is dropped, but instead it is assumed that Condition
(?)
3.4 holds for some positive matrix Σ, which is then used to define β M , Theorem 3.6 continues
to hold. Note that this version of Theorem 3.6 also covers the case of nonrandom design
matrices for which n−1 X 0 X converges to a positive definite limit at rate n−1/2 (or faster).
44

D.3

Extension of Theorem 3.6 to the case p → ∞

We now provide a variant of Theorem 3.6 where we consider the same setup as in Section 3 but
now allow p to depend on n such that p = p(n) → ∞ as n → ∞. The positive definite second
moment matrix Σ now depends on n, i.e., Σ = Σ(n), and has dimension p(n) × p(n). Note
that p(n) ≤ n holds. For a symmetric non-negative definite matrix A, we let A1/2 denote the
unique symmetric non-negative definite square root of A. Furthermore, set Σ(n) = X 0 X/n.
The theorem below relies on the following two conditions, the second of which is somewhat
intransparent.
Condition D.4. We have for n → ∞ that
An (X) :=

(Σ[M, M ])−1/2 Σ(n) [M, M ](Σ[M, M ])−1/2 − I|M | = op (1),

max

M ∈M,M 6=∅

where the norm here is the spectral norm.
Condition D.5. The model selection procedure satisfies for all  > 0 that
 √
h
sup
Pn,β,σ
n(Σ[M̂ , M̂ ])1/2 (Σ(n) [M̂ , M̂ ])−1 Σ(n) [M̂ , M̂ c ]
β∈Rp ,σ>0,x0 6=0


i1
− (Σ[M̂ , M̂ ])−1 Σ[M̂ , M̂ c ] β[M̂ c ] ≥ K1 (x0 , ∞) X
σ
converges to 0 in probability as n → ∞.
The confidence interval in the subsequent theorem assumes knowledge of σ for defining
the length of the interval. However, the model selection procedure M̂ may or may not make
use of this knowledge.
Theorem D.6. Suppose that Conditions D.4 and D.5 hold. Let δ > 0 be given. Denote by
CI δ (x0 ) the interval obtained from (5) by replacing K(x0 , M̂ ) by (1 + δ)K1 (x0 , ∞) and σ̂ by
σ. Then the interval CI δ (x0 ) satisfies


(?)
0
δ
inf
P
x
[
M̂
]β
∈
CI
(x
)|X
≥ 1 − α + op (1),
0
n,β,σ
0
p
p
M̂

x0 ∈R ,β∈R ,σ>0

where the op (1) term above depends only on X (and possibly on the sequence of second moment
matrices Σ = Σ(n)) and converges to zero in probability as n → ∞. This result a fortiori
holds if in the definition of CI δ (x0 ) the constant K1 (x0 , ∞) is replaced by any one of the
constants K2 (x0 [M̂ ], M̂ , ∞), K3 (x0 [M̂ ], M̂ , ∞), K4 (∞), or K5 (∞), respectively.
Proof of Theorem D.6: Obviously it suffices to prove the result for the case where
(?)
K1 (x0 , ∞) is used. Since x00 [M̂ ]β ∈ CI δ (x0 ) trivially holds when x0 = 0, it suffices to show
M̂
that


(?)
0
0
[
M̂
]β
−
x
[
M̂
]
β̂
>
(1
+
δ)K
(x
,
∞)
s
inf
P
x
σ
X
≤ α + op (1).
1
0
n,β,σ
0
0
M̂
M̂
p
β∈R ,σ>0,x0 6=0

M̂

45

Now, we obtain from the triangular inequality,


(?)
Pn,β,σ x00 [M̂ ]β − x00 [M̂ ]β̂ M̂ > (1 + δ)K1 (x0 , ∞) sM̂ σ X
M̂


(n)
≤ Pn,β,σ x00 [M̂ ]β − x00 [M̂ ]β̂ M̂ > K1 (x0 , ∞) sM̂ σ X
M̂


(?)
(n)
> δK1 (x0 , ∞) sM̂ σ X .
+ Pn,β,σ x00 [M̂ ]β − x00 [M̂ ]β
M̂

(35)

M̂

The first probability on the r.h.s. in the preceding display is not larger than α as a consequence
of Proposition 2.1 and Remarks 2.5, 2.6(iii). For the second probability we have


(?)
(n)
> δK1 (x0 , ∞) sM̂ σ X
Pn,β,σ x00 [M̂ ]β − x00 [M̂ ]β
M̂

M̂

i β[M̂ c ]
x00 [M̂ ] (Σ(n) [M̂ , M̂ ])−1 Σ(n) [M̂ , M̂ c ] − (Σ[M̂ , M̂ ])−1 Σ[M̂ , M̂ c ]
σ

1
> √ (Σ(n) [M̂ , M̂ ])−1/2 x0 [M̂ ] δK1 (x0 , ∞) X
n

h
√
≤ Pn,β,σ (Σ[M̂ , M̂ ])−1/2 x0 [M̂ ]
n(Σ[M̂ , M̂ ])1/2 (Σ(n) [M̂ , M̂ ])−1 Σ(n) [M̂ , M̂ c ]
= Pn,β,σ

h

−1

− (Σ[M̂ , M̂ ])

(36)

i β[M̂ c ]
Σ[M̂ , M̂ ]
> (Σ(n) [M̂ , M̂ ])−1/2 x0 [M̂ ] δK1 (x0 , ∞) X
σ
c

!
.

With λmin (.) denoting the smallest eigenvalue of a symmetric matrix, we have
2

(Σ(n) [M̂ , M̂ ])−1/2 x0 [M̂ ] = (Σ(n) [M̂ , M̂ ])−1/2 (Σ[M̂ , M̂ ])1/2 (Σ[M̂ , M̂ ])−1/2 x0 [M̂ ]


2
≥ λmin (Σ[M̂ , M̂ ])1/2 (Σ(n) [M̂ , M̂ ])−1 (Σ[M̂ , M̂ ])1/2 (Σ[M̂ , M̂ ])−1/2 x0 [M̂ ] .

2

On the event where M̂ 6= ∅ it holds that


λmin (Σ[M̂ , M̂ ])1/2 (Σ(n) [M̂ , M̂ ])−1 (Σ[M̂ , M̂ ])1/2 ≥ 1 − An (X).
This gives
(Σ(n) [M̂ , M̂ ])−1/2 x0 [M̂ ]

2

≥ (1 − An (X)) (Σ[M̂ , M̂ ])−1/2 x0 [M̂ ]

2

.

Note that on the event inside the probability on the far r.h.s. of (36) M̂ 6= ∅ and x0 [M̂ ] 6= 0
must hold (in view of our conventions). Furthermore, 1 − An (X) > 1/4 holds with probability
going to one as n → ∞ in view of Condition D.4. Hence, with op (1) denoting a term that
goes to zero in probability as n → ∞ and only depends on X, (35) is not larger than
h
 √
n(Σ[M̂ , M̂ ])1/2 (Σ(n) [M̂ , M̂ ])−1 Σ(n) [M̂ , M̂ c ]
α + op (1) + Pn,β,σ

i1
−1
c
c
− (Σ[M̂ , M̂ ]) Σ[M̂ , M̂ ] β[M̂ ] ≥ δK1 (x0 , ∞)/2 X
σ
The result then follows from Condition D.5. 
46

Remark D.7. If C = C(n) is a subset of Rp(n) × (0, ∞) and only a weaker version of
Condition D.5 is assumed to hold, where the supremum in this condition only extends over
(β, σ) ∈ C and x0 6= 0, then the same proof delivers a version of Theorem D.6 where the
infimum in the display in that theorem now extends over (β, σ) ∈ C and all x0 .

D.4

A comment on notation

In Section 2 the probability measure Pn,µ,σ denotes the joint distribution of Y and σ̂ 2 . In
the case where σ̂ 2 depends on extraneous data besides X and Y , the joint distribution of
Y and σ̂ 2 can in principle depend not only on µ and σ, but also on additional parameters
governing the joint distribution of Y and the extraneous data. [Of course, this does not occur
if σ̂ 2 depends only on X and Y .] Such a dependence on additional parameters is, however,
not relevant in the context of Section 2 and thus is not shown in the notation, since all the
results in that section are based on (10) and since the probability in (10) is free from any
such additional parameters, the reason being that the event inside this probability can be
expressed in terms of (PX Y, σ̂) only, with PX Y and σ̂ being independent by assumption and
their marginal distributions only depending on µ and σ (and the degree of freedom parameter
r).
In the context of Section 3 various probability statements such as (19) may in general
depend on additional parameters governing the joint distribution of (X, Y ) and σ̂ 2 , which,
however, are not shown in the notation. The assumptions and results in this section (e.g.,
Theorem 3.6) are hence to be read as holding for given values (or a given sequence of values
depending on n) of such additional parameters. As a consequence, the results in Section 3 do
not provide uniformity guarantees w.r.t. the additional parameters in general. Of course, if
all the conditions can be checked to hold for all possible sequences of values of the additional
parameters in a certain domain, uniformity w.r.t. the additional parameters is obtained.

E

Appendix: Algorithms for computing the confidence intervals

In this appendix we consider the setting of Section 2. In particular, recall that X is a fixed
n × p matrix of rank d ≥ 1. Let Q be a n × d matrix so that the columns of Q form
an orthonormal basis of the column space of X. Following Berk et al. (2013a) we define
Ỹ = Q0 Y and X̃ = Q0 X, the so-called canonical coordinates of Y and X, cf. Section 5.1 in

Berk et al. (2013a). We then have Ỹ = µ̃ + Ũ with µ̃ = Q0 µ and Ũ = Q0 U ∼ N 0, σ 2 Id .
It is now easy to see that PX̃ Ỹ = Q0 PX Y and QPX̃ Ỹ = PX Y hold. In particular, the
independence of the given σ̂ 2 from the projection of the data vector on the space spanned
by the regressor holds whether we work with the original data or with the data in canonical
−1
coordinates. Furthermore, setting s̃0M = x00 [M ](X̃[M ]0 X̃[M ]) X̃[M ]0 for ∅ 6= M ∈ M
with M as in Section 2 and s̃0M = 0 ∈ Rd for M = ∅, it follows that ks0M k = ks̃0M k and
s0M (Y − µ) = s̃0M (Ỹ − µ̃). For later use define s̃M = s̃M / ks̃M k if ks̃M k =
6 0 and define s̃M = 0
if s̃M = 0. Inspection of (10) and of the definition of FM,x0 now shows that all the constants
Ki remain the same whether they are computed from the original problem using the design
47

matrix X or from the transformed problem using the canonical coordinates X̃ (but using
the originally given σ̂ 2 in both cases). Hence, in the algorithms below we shall work with
the canonical coordinates as this facilitates computation. Note that x0 is unaffected by this
transformation. In the important case d = p(≤ n) the matrices Q and X̃ can be obtained,
for example, from a SVD or a QR decomposition of X, cf. Section 5.1 in Berk et al. (2013a).
[In case p ≥ n = d, one can always set Q = In .]
The following algorithm for computing K1 (x0 ) is similar to that of Berk et al. (2013b) for
computing the PoSI constant. We present it here for completeness. From Proposition 2.3 and
from the arguments used to prove (23) in Appendix B we see that, in case x0 6= 0, K1 (x0 ) is
the solution to


0
EG Pr max s̃M V ≤ t/G G = 1 − α,
M ∈M

where V here is uniformly distributed on the unit sphere of Rd , independently of G (and G
is as in Section 2). The algorithm now replaces the (conditional) probability in the preceding
display by a Monte-Carlo estimator, analytically performs the integration w.r.t. G, and then
numerically solves the resulting equation. We note that in this and the other algorithms to
]
the c.d.f.
follow there is no need for Monte-Carlo integration w.r.t. G. We shall denote by Fd,r

]
2
of G; note that then Fd,r (t) = Fd,r t /d , where Fd,r denotes the c.d.f. of an F -distribution
with (d, r)-degrees of freedom.
Algorithm E.1. In case x0 6= 0, choose I ∈ N and generate independent identically distributed random vectors V1 , . . . , VI , where each Vi is uniformly distributed on the unit sphere
0
in Rd . Calculate the quantities ci = maxM ∈M s̃M Vi with s̃M as defined above. A numerical
approximation to K1 (x0 ) is then obtained by searching for that value of K that solves
I

1X ]
Fd,r
I
i=1



K
ci


= 1 − α.

(37)

In case x0 = 0, set K1 (x0 ) = 0.
Note that for x0 6= 0 at least one of the vectors s̃M , M ∈ M, is non-zero,
  implying
]
that the quantities ci are all non-zero with probability 1; hence the terms Fd,r K
are wellci
defined with probability one. It is now obvious that – on the event where all ci are non-zero
– the solution K of (37) exists, is unique and positive. The costly factor in Algorithm E.1
is the maximization involved in the computation of the quantities ci , while searching for the
value of K that solves (37), for example by bisection searches, incurs only negligible cost.
In our simulations, computing K1 (x0 ) for p = d = 10 (with M the power set of {1, ..., p})
and I = 10, 000 takes around one second on a personal computer; and around 10 minutes
for p = d = 20 and I = 1, 000. In case M is the power set of {1, ..., p}, the complexity of
Algorithm E.1 will be exponential in p and thus will be feasible only for moderately large
values of p. In relation to this we mention that Berk et al. (2013b) found their algorithm
(which is similar to Algorithm E.1 as noted above) to be tractable for up to about p = 20 and
I = 1, 000, in which case the elapsed time was around one hour on 2012 desktop computer
48

equipment. [The longer running time is due to the fact that Berk et al. (2013b) have to search
over p2p−1 unit vectors, while we have to search only over 2p unit vectors.]
The algorithm for computing K3 (x0 [M ], M ) is given next. We provide this algorithm only
for non-empty M 6= {1, ..., p} since in case M = {1, ..., p} we have K3 (x0 [M ], M ) = K1 (x0 ),
which can be computed by Algorithm E.1, and in case M is empty we have K3 (x0 [M ], M ) =
K4 , which can be computed by Algorithm E.3 given below. We now search for the solution
of the equation
∗
1 − α = EG F̆M,x
(t/G)
0
∗
∗
where F̆M,x
is a Monte-Carlo estimator of FM,x
obtained by replacing the probability in0
0

volving V by an empirical Monte-Carlo estimator (and where s̃M instead of s̄M is being used).
∗
Observing that we need only to integrate over the range where F̆M,x
is positive (i.e., where
0
t/G > m∗ defined below), the integrand can be additively decomposed into a ‘jump’ part and
a continuous part. The integral over the jump part can be expressed analytically in terms of
]
, whereas the integral over the continuous part is approximated by an integral
the c.d.f. Fd,r
]
. Recall that
over a step function, which again can be expressed in terms of the c.d.f. Fd,r
c (M, M) has been defined subsequent to (12).

Algorithm E.2. Suppose that M ∈ M satisfies ∅ 6= M 6= {1, ..., p}. Choose I ∈ N, generate
independent identically distributed random vectors V1 , . . . , VI , where each Vi is uniformly dis0
tributed on the unit sphere in Rd , and calculate the quantities ci = maxM∗ ∈M,M∗ ⊆M s̃M∗ Vi
with s̃M∗ as defined above. In case d > 1, find m∗ as the smallest value such that
I


1X
1 (ci > t) + c (M, M) 1 − FBeta,1/2,(d−1)/2 t2 < 1
I
i=1

holds for all t > m∗ . Next, choose J ∈ N, J > 1, and find the values m1 , ..., mJ−1 so that, for
j = 1, ..., J − 1
1 − FBeta,1/2,(d−1)/2 m2∗


 j
= 1 − FBeta,1/2,(d−1)/2 m2j
J

(38)

holds. Set mJ = m∗ . A numerical approximation to K3 (x0 [M ], M ) is then obtained by
searching for that value of K that solves




 

1 X
K
K
K
]
]
]
−
Fd,r
− Fd,r
1 − α = Fd,r
mJ
I
mJ
ci
i:ci >mJ

+ c (M, M) 1 −

FBeta,1/2,(d−1)/2 m2J



X ]  K 
 1 J−1
K
]
Fd,r
− Fd,r
.
J
mj
mJ

(39)

j=1

In case d = 1, K3 (x0 [M ], M ) is the (uniquely determined and positive) constant K that solves
]
1 − α = F1,r
(K) .

49

Note that m∗ exists, is uniquely determined, is always positive, and satisfies m∗ ≤ 1. [In
fact, m∗ < 1 holds, except in case ci = 1 for all i, which is a probability zero event.] Provided
m∗ < 1 holds, the values mj for j ≥ 1 are uniquely defined and satisfy m∗ < mJ−1 < . . . <
m1 < 1. [In case m∗ = 1, then any mj ≥ 1 would solve (38). But in this case the r.h.s.
]
of (39) reduces to Fd,r
(K) anyway and hence there is no need for solving equation (38).]
Furthermore, note that the r.h.s. of (39) can be written as
#

"
I
X

K
1
J
−
1
]
Fd,r
1−
1 (ci > mJ ) − c (M, M) 1 − FBeta,1/2,(d−1)/2 m2J
mJ
I
J
i=1
 
X ] K
 1 J−1
1 X
K
]
2
+
+ c (M, M) 1 − FBeta,1/2,(d−1)/2 mJ
.
Fd,r
Fd,r
I
ci
J
mj
i:ci >mJ

j=1

Observing that the expression in brackets is nonnegative (in fact, positive) because of the
definition of mJ , we see that the r.h.s. of (39) is strictly increasing in K. Furthermore,
inspection of the r.h.s. of (39) shows that it is zero for K = 0 and converges to one for
K → ∞. Consequently, equation (39) has a unique solution for K, which necessarily is
positive. We note that in Algorithm E.2 the cost of searching for m∗ , for the mj ’s, and for K,
for example by bisection searches, is negligible compared to that of computing the quantities
ci , which is again the limiting factor.

The above algorithm is based on approximating 1 − FBeta,1/2,(d−1)/2 t2 for t > m∗ by a
step function from below. If we approximate by a step function from above, this results in the
same algorithm except that now the second sum on the r.h.s. of equation (39) runs from j = 0
to j = J − 1 with the convention that m0 = 1. A similar argument as above shows that the
solution to this modification of (39) exists, is unique and is positive. Note that the solutions
obtained from running both versions of the algorithm in parallel provide a lower as well as an
upper bound for the solution one would obtain if the integration of the continuous part could
be performed without error. These lower and upper bounds allow one to gauge whether or
not J has been chosen large enough such that the effect of the numerical integration error on
K is negligible. Note that running the two versions of the algorithm in parallel is not much
more costly than running just one version, as only (bisection) searches are involved once the
ci ’s have been computed.
The following algorithm for computing K4 is similar to the algorithm in Berk et al. (2013b),
Section 7.2, for computing the universal upper-bound for the PoSI constants. The computational cost of this algorithm is negligible compared to those of Algorithms E.1 and E.2.
Algorithm E.3. In case d > 1, choose J ∈ N, J > 1, and find the values m1 , ..., mJ so that,
for j = 1, ..., J,

j
c (∅, M) 1 − FBeta,1/2,(d−1)/2 m2j = .
(40)
J
Then, K4 is numerically approximated by the (uniquely determined and positive) constant K
that solves
 
1 X ]
K
Fd,r
= 1 − α.
(41)
J
mj
mj >0

50

In case d = 1, K4 is the (uniquely determined and positive) constant K that solves
]
F1,r
(K) = 1 − α.

(42)

Note that in case d > 1 the constants mj always exist and are unique; they are all positive
in case M 6= {∅, {1, . . . , p}} (as then c (∅, M) > 1 must hold in view of our assumptions on
M), and they are positive for j = 1, . . . , J − 1 in case M = {∅, {1, . . . , p}}. Consequently,
the solution K of (41) exists, is unique and positive. In case d = 1 the solution of (42) also
exists, is unique and positive. As before, this algorithm relies on approximation by a step
function from below. A version of the algorithm that uses a step function that approximates
from above is obtained if equation (41) is replaced by
 
J−1
1X ]
K
=1−α
Fd,r
J
mj
j=0

with the convention that m0 = 1.
Remark E.4. For the computation of the constants K1 (x0 , ∞), K3 (x0 [M ], M, ∞), and
K4 (∞) (cf. Remark 2.5) one can use the above algorithms with the only modification that
]
the distribution function Fd,r
is replaced by the distribution function of the square root of a
chi-squared-distributed random variable with d degrees of freedom.
Remark E.5. When the collection M becomes large (e.g., if M is the power set of {1, ..., p}
in case p = d ≤ n and d is larger than 20), Algorithms E.1 or E.2 may not be tractable,
but Algorithm E.3 can still be as it does not require the costly step of searching over the
model universe M. However, it is reported in Berk et al. (2013b) that, for about d ≥ 40,
it can be problematic to compute the extreme quantiles in (40) with standard routines. In
this case, one can of course always use the Scheffé constant K5 . In practice, one may also
consider in such cases (since p is large) to use rule-of-thumb constants smaller than K5 that
are based on asymptotic considerations such as Corollary 2.12: For example, if p = d ≤ n,
M is the power set of {1, ..., p}, but p is very large, this corollary could be read as suggesting
to use the constant K6 = 0.866K5 in (5). A similar advice is given in the framework of Berk
et al. (2013b). [In case p > n and M is as in Corollary 2.12, this corollary can be used to
provide appropriate substitutes for K6 .] However, we would like to issue a warning here: The
asymptotic results for p → ∞ like Corollary 2.12 and the related results in Berk et al. (2013b)
and Berk et al. (2013a) are highly non-uniform w.r.t. α (cf. Remark 2.14), showing that
rule-of-thumb approximations such as K6 have to be taken with a grain of salt; see also the
warning expressed at the end of Section 5.2 of Berk et al. (2013b).

F
F.1

Appendix: Details for computations in Section 4
Description of the variables in the watershed data set

The explanatory variables in XRaw are a constant term (to include an intercept in the model),
rainfall (inches), area of watershed (square miles), area impervious to water (square miles),
51

average slope of watershed (percent), longest stream flow in watershed (thousands of feet),
surface absorbency index (0= complete absorbency; 100 = no absorbency), estimated soil
storage capacity (inches of water), infiltration rate of water into soil (inches/hour) and time
period during which rainfall exceeded 1/4 inch/hour. Logarithms are taken of the explanatory
variables except for the intercept. [In Rawlings et al. (1998), the response corresponding to
these explanatory variables is peak flow rate from watersheds.]

F.2

Three-step Monte Carlo procedure for determining K2 in Section 4.1

For the given X, x0 , and M we proceed as follows: First, we randomly sample 100, 000
independent vectors x ∈ R10 , so that x[M ] = x0 [M ] and x[M c ] follows a Gaussian distribution
with mean vector 0 ∈ R10−|M | and covariance matrix (1/n)(X[M c ]0 X[M c ]). For each of these
vectors, we evaluate K1 (x) with Algorithm E.1 in Appendix E, with I1 = 1, 000 Monte
Carlo samples. In the second step, we keep the 1, 000 vectors x corresponding to the largest
evaluations of K1 (x) and we reevaluate K1 (x) for them, with a number of Monte Carlo samples
equal to I2 = 100, 000 in Algorithm E.1. In the third step, we keep the vector x from the
second step corresponding to the largest value of K1 and we reevaluate K1 (x) for this x, but
this time with a number of Monte Carlo samples equal to I3 = 1, 000, 000 in Algorithm E.1.

F.3

Details for computing AIC, BIC, LASSO, MCP, and SCAD in Section
4.2

For the AIC- and BIC-procedures we use the step() function in R, with penalty parameter k
equal to 2 for AIC and log(n) for BIC. The AIC and BIC objective functions are minimized
through a greedy general-to-specific search over the resulting 2p−1 candidate models (recall
that the intercept is protected).
For the LASSO, the selected model corresponds to the explanatory variables for which
the LASSO estimator has non-zero coefficients. More precisely, we use the lars package in R
and follow suggestions outlined in Efron et al. (2004): To protect the first regressor, we first
compute the residual of the orthogonal projection of Y on the first regressor; write Ỹ for this
residual vector, and write X̃ for the design matrix X with the first column removed. We then
compute the LASSO-estimator for a regression of Ỹ on X̃ using the lars() function; the
LASSO-penalty is chosen by 10-fold cross-validation using the cv.lars() function. In both
functions we set the intercept parameter to FALSE, but otherwise use the default settings.
The selected model is comprised of those regressors in X̃ for which the corresponding LASSO
coefficients are non-zero, plus the first column of X.
For SCAD and MCP, we use the ncvreg package in R. With the function cv.ncvreg()
(with parameters SCAD or MCP) the penalty is selected by 10-fold cross-validation and the corresponding estimated regression coefficients are computed. Like for the LASSO, the function
cv.ncvreg() is applied for a regression of Ỹ on X̃, and the selected model is comprised of
those regressors in X̃ for which the SCAD (or MCP) coefficients are non-zero, plus the first
column of X.

52

F.4

On the exchangeable and the equicorrelated data set

The matrix Σ̃ in the exchangeable case is related to the design matrix X(p̃) (a) defined in
Section 6.1 of Berk et al. (2013a) via Σ̃ = (X(p̃) (a))0 (X(p̃) (a)), where p̃ = 9 and a = 10. In the
equicorrelated case Σ̃ is related to the design matrix X(p̃) (c) defined p
in Section 6.2 of Berk
et al. (2013a) again via Σ̃ = (X(p̃) (c))0 (X(p̃) (c)), where p̃ = 9 and c = 0.8/(p̃ − 1).

F.5

Three-step Monte Carlo procedure for estimating minimal coverage
probabilities in Section 4.2

For each configuration of n, Σ, the model selection procedure, the target (either the designdependent or the design-independent target), as well as of a matrix X and a vector x0 as
described in Section 4.2 and for each of the constants Knaive , K1 , K3 , and K4 we estimate the
minimal (over β and σ) coverage probabilities (conditional on X and x0 ) of the confidence
intervals by a three-step Monte Carlo procedure as follows: We first sample independently
m1 = 1, 000 parameters β from a p-dimensional random vector b where Xb follows a standard Gaussian distribution within the column-space of X. Then, for each of these vectors
β, we draw I1 = 1000 Monte Carlo samples from the full model (i.e., from a N (Xβ, σ 2 In )distribution) using β and σ = 1 as the true parameters. [For invariance reasons it suffices
to consider only the case where σ = 1.] For each Monte Carlo sample, we use the standard
unbiased estimator σ̂ 2 of the error variance (under the full linear model), we carry out the
model-selection procedure M̂ , and we record whether or not the target currently under investigation is covered by the confidence interval obtained from (5) with K(x0 , M̂ ) replaced
by the constant K under investigation. [For K = Ki , i = 1, 3, 4, the value of Ki is obtained
from the algorithms described in Appendix E.] For each β, the I1 recorded results are then
averaged, resulting in m1 Monte Carlo estimates of the coverage probabilities depending on
the m1 sampled vectors β. Then for the m2 = 100 vectors β corresponding to the smallest
estimated coverage probabilities from the first step, we repeat the Monte Carlo procedures,
but this time with I2 = 10, 000 Monte Carlo samples, and we record the vector β that yields
the smallest estimate for the coverage probability in this second step. Performing these two
steps for each of the four constants Knaive , K1 , K3 , and K4 results in four vectors β(1), β(2),
β(3), and β(4). In a third step, we now reevaluate the coverage probability of any of the four
confidence intervals at each of the vectors β(j), j = 1, . . . , 4, this time now with I3 = 100, 000
Monte Carlo samples, and record, for each of the confidence intervals, the minimum of these
four estimates of the coverage probabilities. This is then used as the final estimate of the
minimal coverage probability of the confidence interval under consideration.

F.6

Three-step Monte Carlo procedure for estimating minimal conditional
coverage probabilities in Section 4.3

For each of the eight configurations mentioned in Section 4.3, we carry out the three-step
minimal coverage probability evaluation described just above, with the same values of m1 , m2
and I1 , I2 , I3 (with the only difference that in the third step the (now conditional) coverage
probability is reevaluated only for one value β(1), say, where β(1) corresponds to that value of
β that gives the smallest estimate for the coverage probability in the second step). When we
53

evaluate a conditional coverage probability for a given β in this process, we proceed as follows:
We sample I1 (or I2 or I3 ) values of Y from the N (Xβ, In ) distribution. For each value of
Y we run the model selection procedure M̂ , where λ = λ̂ is first selected by cross-validation
with the cv.glmnet function of the R package glmnet, and where the glmnet function is then
used to compute the selected model with the LASSO with penalty parameter λ̂. Then, if the
selected model does not contain the first explanatory variable, we discard the value of Y , and
¯ or not. The conditional
else, we record whether the design-dependent target belongs to CI
coverage probability is then obtained by taking the average number of times this is the case,
over all the recorded events.

References
Andrews, D. W. K. and Guggenberger, P. (2009). Hybrid and size-corrected subsampling methods. Econometrica, 77 721–762.
Bachoc, F., Preinerstorfer, D. and Steinberger, L. (2016). Uniformly valid confidence
intervals post-model-selection. ArXiv:1611.01043.
Belloni, A., Chernozhukov, V. and Hansen, C. (2011). Inference for high-dimensional
sparse econometric models. Advances in Economics and Econometrics. 10th World Congress
of the Econometric Society, Volume III, 245–295.
Belloni, A., Chernozhukov, V. and Hansen, C. (2014). Inference on treatment effects
after selection among high-dimensional controls. Rev. Econom. Stud., 81 608–650.
Berk, R., Brown, L., Buja, A., Zhang, K. and Zhao, L. (2013a). Valid post-selection
inference. Ann. Statist., 41 802–837.
Berk, R., Brown, L., Buja, A., Zhang, K. and Zhao, L. (2013b). Valid post-selection
inference. Unpublished version, URL http://www-stat.wharton.upenn.edu/~lzhao/
papers/MyPublication/24PoSI-submit.pdf.
Castera, L., Chan, H., Arrese, M., Afdhal, N., Bedossa, P., Friedrich-Rust, M.,
Han, K.-H. and Pinzani, M. (2015). EASL-ALEH clinical practice guidelines: noninvasive tests for evaluation of liver disease severity and prognosis. Journal of Hepatology,
63 237–264.
Efron, B., Hastie, T., Johnstone, I. and Tibshirani, R. (2004). Least angle regression.
Ann. Statist., 32 407–499.
Ewald, K. (2012). On the influence of model selection on confidence regions for marginal
associations in the linear model. Master’s thesis, University of Vienna.
Fan, J. and Li, R. (2001). Variable selection via nonconcave penalized likelihood and its
oracle properties. J. Amer. Statist. Assoc., 96 1348–1360.

54

Fithian, W., Sun, D. and Taylor, J. (2015). Optimal inference after model selection.
ArXiv:1410.2597.
Greenshtein, E. and Ritov, Y. (2004). Persistence in high-dimensional linear predictor
selection and the virtue of overparametrization. Bernoulli, 10 971–988.
Jaupi, L. (2014). Variable selection methods for multivariate process monitoring. In Proceedings of the World Congress of Engineering 2014 (S. Ao, L. Gelman, D. Hukins, A. Hunter
and A. Korsunsky, eds.), vol. II. 1116–1120.
Kabaila, P. and Leeb, H. (2006). On the large-sample minimal coverage probability of
confidence intervals after model selection. J. Amer. Statist. Assoc., 101 619–629.
Lee, J. D., Sun, D. L., Sun, Y., and Taylor, J. E. (2016). Exact post-selection inference,
with application to the lasso. Ann. Statist., 44 907–927.
Lee, J. D. and Taylor, J. (2014). Exact post model selection inference for marginal screening. In Advances in Neural Information Processing Systems 27 (Z. Ghahramani, M. Welling,
C. Cortes, N. Lawrence and K. Weinberger, eds.). Curran Associates, Inc., 136–144.
Leeb, H. (2009). Conditional predictive inference post model selection. Ann. Statist., 37
2838–2876.
Leeb, H. and Pötscher, B. M. (2003). The finite-sample distribution of post-modelselection estimators, and uniform versus non-uniform approximations. Econometric Theory,
19 100–142.
Leeb, H. and Pötscher, B. M. (2005). Model selection and inference: Facts and fiction.
Econometric Theory, 21 21–59.
Leeb, H. and Pötscher, B. M. (2006). Can one estimate the conditional distribution of
post-model-selection estimators? Ann. Statist., 34 2554–2591.
Leeb, H. and Pötscher, B. M. (2017). Testing in the presence of nuisance parameters:
Some comments on tests post-model-selection and random critical values. In Big and Complex Data Analysis: Methodology and Applications (S. Ahmed, ed.). Springer.
Leeb, H., Pötscher, B. M. and Ewald, K. (2015). On various confidence intervals postmodel-selection. Statist. Sci., 30 216–227.
Lockhart, R., Taylor, J., Tibshirani, R. J. and Tibshirani, R. (2014). A significance
test for the LASSO. Ann. Statist., 42 413–468.
Loftus, J. R. (2015). Selective inference after cross-validation. ArXiv:1511.08866.
Loftus, J. R. and Taylor, J. (2015). Selective inference in regression models with groups
of variables. ArXiv:1511.01478.
Pötscher, B. M. (2009). Confidence sets based on sparse estimators are necessarily large.
Sankhya, 71 1–18.
55

Pötscher, B. M. and Schneider, U. (2010). Confidence sets based on penalized maximum
likelihood estimators in Gaussian regression. Electron. J. Statist., 4 334–360.
Rawlings, J. O., Pantula, S. G. and Dickey, D. A. (1998). Applied Regression Analysis:
A Research Tool. 2nd ed. Springer Verlag, New York, NY.
Scheffé, H. (1959). The Analysis of Variance. Wiley, New York.
Schneider, U. (2016). Confidence sets based on thresholding estimators in high-dimensional
Gaussian regression models. Econometric Reviews, 35 1412–1455.
Souders, T. and Stenbakken, G. (1991). Cutting the high cost of testing. IEEE Spectrum,
28 48–51.
Tian, X. and Taylor, J. (2015). Asymptotics of selective inferene. ArXiv:1501.03588.
Tibshirani, R. J., Rinaldo, A., Tibshirani, R. and Wasserman, L. (2015). Uniform
asymptotic inference and the bootstrap after model selection. ArXiv:1506.06266.
Tibshirani, R. J., Taylor, J., Lockhart, R. and Tibshirani, R. (2016). Exact postselection inference for sequential regression procedures. Journal of the American Statistical
Association, 111 600–620.
van de Geer, S., Bühlmann, P., Ritov, Y. and Dezeure, R. (2014). On asymptotically
optimal confidence regions and tests for high-dimensional models. Ann. Statist., 42 1166–
1202.
Wasserman, L. (2014). Discussion: “A significance test for the LASSO”. Ann. Statist., 42
501–508.
Wasserman, L. and Roeder, K. (2009). High-dimensional variable selection. Ann. Statist.,
37 2178–2201.
Zhang, C. (2010). Nearly unbiased variable selection under minimax concave penalty. Ann.
Statist., 38 894–942.
Zhang, C.-H. and Zhang, S. (2014). Confidence intervals for low dimensional parameters
in high dimensional linear models. J. Roy. Statist. Soc. Ser. B, 76 217–242.
Zhang, K. (2013). Rank-extreme association of gaussian vectors and low-rank detection.
ArXiv:1306.0623.

56

