Valid post-selection inference

Introduction: The problem with statistical inference after model selection.

Classical statistical theory grants validity of statistical tests and confidence intervals assuming a wall of separation between the selection of a model and the analysis of the data being modeled. In practice, this separation rarely exists, and more often a model is “found” by a data-driven selection process. As a consequence inferential guarantees derived from classical theory are invalidated. Among model selection methods that are problematic for classical inference, variable selection stands out because it is regularly taught, commonly practiced and highly researched as a technology. Even though statisticians may have a general awareness that the data-driven selection of variables (predictors, covariates) must somehow affect subsequent classical inference from $F$- and $t$-based tests and confidence intervals, the practice is so pervasive that it appears in classical undergraduate textbooks on statistics such as Moore and McCabe 2003.

The reason for the invalidation of classical inference guarantees is that a data-driven variable selection process produces a model that is itself stochastic, and this stochastic aspect is not accounted for by classical theory. Models become stochastic when the stochastic component of the data is involved in the selection process. (In regression with fixed predictors the stochastic component is the response.) Models are stochastic in a well-defined way when they are the result of formal variable selection procedures such as stepwise or stagewise forward selection or backward elimination or all-subset searches driven by complexity penalties (such as $C_{p}$, AIC, BIC, risk-inflation, LASSO$,\ldots$) or prediction criteria such as cross-validation, or more recent proposals such as LARS and the Dantzig selector; for an overview see, for example, Hastie, Tibshirani and Friedman 2009. Models are also stochastic but in an ill-defined way when they are informally selected through visual inspection of residual plots or normal quantile plots or other regression diagnostics. Finally, models become stochastic in an opaque way when their selection is affected by human intervention based on post-hoc considerations such as “in retrospect only one of these two variables should be in the model” or “it turns out the predictive benefit of this variable is too weak to warrant the cost of collecting it.” In practice, all three modes of variable selection may be exercised in the same data analysis: multiple runs of one or more formal search algorithms may be performed and compared, the parameters of the algorithms may be subjected to experimentation and the results may be critiqued with graphical diagnostics; a round of fine-tuning based on substantive deliberations may finalize the analysis.

Posed so starkly, the problems with statistical inference after variable selection may well seem insurmountable. At a minimum, one would expect technical solutions to be possible only when a formal selection algorithm is (1) well-specified (1a) in advance and (1b) covering all eventualities, (2) strictly adhered to in the course of data analysis and (3) not “improved” on by informal and post-hoc elements. It may, however, be unrealistic to expect this level of rigor in most data analysis contexts, with the exception of well-conducted clinical trials. The real challenge is therefore to devise statistical inference that is valid following any type of variable selection, be it formal, informal, post hoc or a combination thereof. Meeting this challenge with a relatively simple proposal is the goal of this article. This proposal for valid Post-Selection Inference, or “PoSI” for short, consists of a large-scale family-wise error guarantee that can be shown to account for all types of variable selection, including those of the informal and post-hoc varieties. On the other hand, the proposal is no more conservative than necessary to account for selection, and in particular it can be shown to be less conservative than Scheffé’s simultaneous inference.

The framework for our proposal is in outline as follows—details to be elaborated in subsequent sections: we consider linear regression with predictor variables whose values are considered fixed, and with a response variable that has normal and homoscedastic errors. The framework does not require that any of the eligible linear models is correct, not even the full model, as long as a valid error estimate is available. We assume that the selected model is the result of some procedure that makes use of the response, but the procedure does not need to be fully specified. A crucial aspect of the framework concerns the use and interpretation of the selected model: we assume that, after variable selection is completed, the selected predictor variables—and only they—will be relevant; all others will be eliminated from further consideration. This assumption, seemingly innocuous and natural, has critical consequences: it implies that statistical inference will be sought for the coefficients of the selected predictors only and in the context of the selected model only. Thus the appropriate targets of inference are the best linear coefficients within the selected model, where each coefficient is adjusted for the presence of all other included predictors but not those that were eliminated. Therefore the coefficient of an included predictor generally requires inference that is specific to the model in which it appears. Summarizing in a motto, a difference in adjustment implies a difference in parameters and hence in inference. The goal of the present proposal is therefore simultaneous inference for all coefficients within all submodels. Such inference can be shown to be valid following any variable selection procedure, be it formal, informal, post hoc, fully or only partly specified.

Problems associated with post-selection inference were recognized long ago, for example, by Buehler and Feddersen 1963, Brown 1967, Olshen 1973, Sen 1979, Sen and Saleh 1987, Dijkstra and Veldkamp 1988, Pötscher 1991, Kabaila 1998. More recently specific problems have been the subject of incisive analyses and criticisms by the “Vienna School” of Pötscher, Leeb and Schneider; see, for example, Leeb and Pötscher (Leeb and Pötscher 2003; Leeb and Pötscher 2005; Leeb and Pötscher 2006a; Leeb and Pötscher 2006b; Leeb and Pötscher 2008a; Leeb and Pötscher 2008b; Leeb and Pötscher 2008c), Pötscher 2006, Leeb 2006, Pötscher and Leeb 2009, Pötscher and Schneider (Pötscher and Schneider 2009; Pötscher and Schneider 2010; Pötscher and Schneider 2011), as well as Kabaila and Leeb 2006 and Kabaila 2009. Important progress was made by Hjort and Claeskens 2003 and Claeskens and Hjort 2003.

This article proceeds as follows: in Section 2 we first develop the “submodel view” of the targets of inference after model selection and contrast it with the “full model view” (Section 2.1); we then introduce assumptions with a view toward valid inference in “wrong models” (Section 2.2). Section 3 is about estimation and its targets from the submodel point of view. Section 4 develops the methodology for PoSI confidence intervals (CIs) and tests. After some structural results for the PoSI problem in Section 5, we show in Section 6 that with increasing number of predictors $p$ the width of PoSI CIs can range between the asymptotic rates $O(\sqrt{\log{p}})$ and $O(\sqrt{p})$. We give examples for both rates and, inspired by problems in sphere packing and covering, we give upper bounds for the limiting constant in the $O(\sqrt{p})$ case. We conclude with a discussion in Section 7. Some proofs are deferred to the Appendix, and some elaborations to the online Appendix in the supplementary material [Berk et al. 2013].

Computations will be described in a separate article. Simulation-based methods yield satisfactory accuracy specific to a design matrix up to $p\approx 20$, while nonasymptotic universal upper bounds can be computed for larger $p$.

Targets of inference and assumptions.

It is a natural intuition that model selection distorts inference by distorting sampling distributions of parameter estimates: estimates in selected models should tend to generate more type I errors than conventional theory allows because the typical selection procedure favors models with strong, hence highly significant predictors. This intuition correctly points to a multiplicity problem that grows more severe as the number of predictors subject to selection increases. This is the problem we address in this article.

Model selection poses additional problems that are less obvious but no less fundamental: there exists an ambiguity as to the role and meaning of the parameters in submodels. In one view, the relevant parameters are always those of the full model, hence the selection of a submodel is interpreted as estimating the deselected parameters to be zero and estimating the selected parameters under a zero constraint on the deselected parameters. In another view, the submodel has its own parameters, and the deselected parameters are not zero but nonexistent. These distinctions are not academic as they imply fundamentally different ideas regarding the targets of inference, the measurement of statistical performance, and the problem of post-selection inference. The two views derive from different purposes of equations:

Underlying the full model view of parameters is the use of a full equation to describe a “data generating” mechanism for the response; the equation hence has a causal interpretation.

Underlying the submodel view of parameters is the use of any equation to merely describe association between predictor and response variables; no data generating or causal claims are implied.

In this article we address the latter use of equations. Issues relating to the former use are discussed in the online Appendix of the supplementary material [Berk et al. 2013, Section B.1].

In what follows we elaborate three points that set the submodel interpretation of coefficients apart from the full model interpretation, with important consequences for the rest of this article:

The full model has no special status other than being the repository of available predictors.

The coefficients of excluded predictors are not zero; they are not defined and therefore do not exist.

The meaning of a predictor’s coefficient depends on which other predictors are included in the selected model.

(1) The full model available to the statistician often cannot be argued to have special status because of inability to identify and measure all relevant predictors. Additionally, even when a large and potentially complete suite of predictors can be measured, there is generally a question of predictor redundancy that may make it desirable to omit some of the measurable predictors from the final model. It is a common experience in the social sciences that models proposed on theoretical grounds are found on empirical grounds to have their predictors entangled by collinearities that permit little meaningful statistical inference. This situation is not limited to the social sciences: in gene expression studies it may well occur that numerous sites have a tendency to be expressed concurrently, hence as predictors in disease studies they will be strongly confounded. The emphasis on full models may be particularly strong in econometrics where there is a “notion that a longer regression …has a causal interpretation, while a shorter regression does not” [Angrist and Pischke 2009, page 59]. Even in causal models, however, there is a possibility that included adjuster variables will “adjust away” some of the causal variables of interest. Generally, in any creative observational study involving novel predictors, it will be difficult a priori to exclude collinearities that might force a rethinking of the predictors. In conclusion, whenever predictor redundancy is a potential issue, it cannot a priori be claimed that the full model provides the parameters of primary interest.

(2) In the submodel interpretation of parameters, claiming that the coefficients of deselected predictors are zero does not properly describe the role of predictors. Deselected predictors have no role in the submodel equation; they become no different than predictors that had never been considered. The selected submodel becomes the vehicle of substantive research irrespective of what the full model was. As such the submodel stands on its own. This view is especially appropriate if the statistician’s task is to determine which predictors are to be measured in the future.

(3) The submodel interpretation of parameters is deeply seated in how we teach regression. We explain that the meaning of a regression coefficient depends on which of the other predictors are included in the model: “the slope is the average difference in the response for a unit difference in the predictor, at fixed levels of all other predictors in the model.” This “ceteris paribus” clause is essential to the meaning of a slope. That there is a difference in meaning when there is a difference in covariates is most drastically evident when there is a case of Simpson’s paradox. For example, if purchase likelihood of a high-tech gadget is predicted from age, it might be found against expectations that younger people have lower purchase likelihood, whereas a regression on age and income might show that at fixed levels of income younger people have indeed higher purchase likelihood. This case of Simpson’s paradox would be enabled by the expected positive collinearity between age and income. Thus the marginal slope on age is distinct from the income-adjusted slope on age as the two slopes answer different questions, apart from having opposite signs. In summary, different models result in different parameters with different meanings.

Must we use the full model with both predictors? Not if income data is difficult to obtain or if it provides little improvement in $R^{2}$ beyond age. The model based on age alone cannot be said to be a priori “wrong.” If, for example, the predictor and response variables have jointly multivariate normal distributions, then every linear submodel is “correct.” These considerations drive home, once again, that sometimes no model has special status.

In summary, a range of applications call for a framework in which the full model is not the sole provider of parameters, where rather each submodel defines its own. The consequences of this view will be developed in Section 3.

2 Assumptions, models as approximations, and error estimates.

We state assumptions for estimation and for the construction of valid tests and CIs when fitting arbitrary linear equations. The main goal is to prepare the ground for valid statistical inference after model selection—not assuming that selected models are correct.

Due to frequent reference we call $d=p$ ($\leq n$) “the classical case.”

It is common practice to assume the full model ${\mathbf{Y}}\sim\mathcal{N}_{n}({\mathbf{X}}{\bm{\beta}},\sigma^{2}{\mathbf{I}})$ to be correct. In the present framework, however, first-order correctness, ${\mathbf{E}}[{\mathbf{Y}}]={\mathbf{X}}{\bm{\beta}}$, will not be assumed. By implication, first-order correctness of any submodel will not be assumed either. Effectively,

is allowed to be unconstrained and, in particular, need not reside in the column space of ${\mathbf{X}}$. That is, the model given by ${\mathbf{X}}$ is allowed to be “first-order wrong,” and hence we are, in a well-defined sense, serious about G.E.P. Box’s famous quote. What he calls “wrong models,” we prefer to call “approximations”: all predictor matrices ${\mathbf{X}}$ provide approximations to ${\bm{\mu}}$, some better than others, but the degree of approximation plays no role in the clarification of statistical inference. The main reason for elaborating this point is as follows: after model selection, the case for “correct models” is clearly questionable, even for “consistent model selection procedures” [Leeb and Pötscher 2003, page 101]; but if correctness of submodels is not assumed, it is only natural to abandon this assumption for the full model also, in line with the idea that the full model has no special status. As we proceed with estimation and inference guarantees in the absence of first-order correctness we will rely on assumptions as follows:

For estimation (Section 3), we will only need the existence of ${\bm{\mu}}={\mathbf{E}}[{\mathbf{Y}}]$.

For testing and CI guarantees (Section 4), we will make conventional second-order and distributional assumptions,

The assumptions (3) of homoscedasticity and normality are as questionable as first-order correctness, and we will report elsewhere on approaches that avoid them. For now we follow the vast model selection literature that relies on the technical advantages of assuming homoscedastic and normal errors.

Accepting the assumption (3), we address the issue of estimating the error variance $\sigma^{2}$, because the valid tests and CIs we construct require a valid estimate ${\hat{\sigma}}^{2}$ of $\sigma^{2}$ that is independent of LS estimates. In the classical case, the most common way to assert such an estimate is to assume that the full model is first-order correct, ${\bm{\mu}}={\mathbf{X}}{\bm{\beta}}$ in addition to (3), in which case the mean squared residual (MSR) ${\hat{\sigma}}_{F}^{2}=\|{\mathbf{Y}}-{\mathbf{X}}{\hat{\bm{\beta}}}\|^{2}/(n-p)$ of the full model will do. However, other possibilities for producing a valid estimate ${\hat{\sigma}}^{2}$ exist, and they may allow relaxing the assumption of first-order correctness:

Exact replications of the response obtained under identical conditions might be available in sufficient numbers. An estimate ${\hat{\sigma}}^{2}$ can be obtained as the MSR of the one-way ANOVA of the groups of replicates.

a larger linear model than the full model might be considered as correct; hence ${\hat{\sigma}}^{2}$ could be the MSR from this larger model.

A different possibility is to use another dataset, similar to the one currently being analyzed, to produce an independent estimate ${\hat{\sigma}}^{2}$ by whatever valid estimation method.

A special case of the preceding is a random split-sample approach whereby one part of the data is reserved for producing ${\hat{\sigma}}^{2}$ and the other part for estimating coefficients, selecting models and carrying out post-model selection inference.

A different type of estimate, ${\hat{\sigma}}^{2}$, may be based on considerations borrowed from nonparametric function estimation [Hall and Carroll 1989].

The purpose of pointing out these possibilities is to separate, at least in principle, the issue of first-order model incorrectness from the issue of error estimation under assumption (3). This separation puts the case $n<p$ within our framework as the valid and independent estimation of $\sigma^{2}$ is a problem faced by all “$n<p$” approaches.

Estimation and its targets in submodels.

Following Section 2.1, the value and meaning of a regression coefficient depends on what the other predictors in the model are. An exception occurs, of course, when the predictors are perfectly orthogonal, as in some designed experiments or in function fitting with orthogonal basis functions. In this case a coefficient has the same value and meaning across all submodels. This article is hence a story of (partial) collinearity.

2 Interpreting regression coefficients in first-order incorrect models.

Hence regression coefficients are weighted averages of case-wise slopes, and this interpretation holds without first-order assumptions.

Universally valid post-selection confidence intervals.

These problems are avoided by using one valid estimate ${\hat{\sigma}}^{2}$ that is independent of all submodels.

If $K=t_{r,1-\alpha/2}$ is the ${1-\alpha/2}$ quantile of a $t$-distribution with $r$ degrees of freedom, then the interval is marginally valid with a $1-\alpha$ coverage guarantee

2 Model selection and its implications for parameters.

Thus the set of parameters for which inference is sought is random also.

3 Post-selection coverage guarantees for confidence intervals.

With randomness of the selected model and its parameters in mind, what is a desirable form of post-selection coverage guarantee for confidence intervals? A natural requirement would be a $1-\alpha$ confidence guarantee for the coefficients of the predictors that are selected into the model,

respectively, are both well-defined and can be the subject of coverage guarantees; see the online Appendix of the supplementary material [Berk et al. 2013, Section B.4].

4 Universal validity for all selection procedures.

(2) There exists a model selection procedure that requires the full strength of universally valid PoSI, and this procedure may not be entirely unrealistic as an approximation to some types of data analytic activities: “significance hunting,” that is, selecting that model which contains the statistically most significant coefficient; see Section 4.9.

(3) There is a general question about the wisdom of proposing ever tighter confidence and retention intervals for practical use when in fact these intervals are valid only under tightly controlled conditions. It might be realistic to suppose that much applied work involves more data peeking than is reported in published articles. With inference that is universally valid after any model selection procedure, we have a way to establish which rejections are safe, irrespective of unreported data peeking as part of selecting a model.

(4) Related to the previous point is the fact that today there is a realization that a considerable fraction of published empirical work is unreproducible or reports exaggerated effects; well known in this regard is Ioannidis 2005. A factor contributing to this problem might well be liberal handling of variable selection and absent accounting for it in subsequent inference.

5 Restricted model selection.

The following are examples of model universes with practical relevance; see also Leeb and Pötscher 2008a, Section 1.1, Example 1.

Nested models: ${\mathcal{M}}_{4}=\{\{1,\ldots,j\}|j\in\{1,\ldots,p\}\}$. $|{\mathcal{M}}_{4}|=p$. Example: selecting the degree up to $p-1$ in a polynomial regression.

Models dictated by an ANOVA hierarchy of main effects and interactions in a factorial design.

6 Reduction of universally valid post-selection inference to simultaneous inference.

We show that universally valid post-selection inference (14) follows from simultaneous inference in the form of family-wise error control for all parameters in all submodels. The argument depends on the following lemma that may fall into the category of the “trivial but not immediately obvious.”

This value will be called “the PoSI constant.” It does not depend on any model selection procedures, but it does depend on the design matrix ${\mathbf{X}}$, the universe ${\mathcal{M}}$ of models subject to selection, the desired coverage $1-\alpha$, and the degrees of freedom $r$ in ${\hat{\sigma}}$, hence $K=K({\mathbf{X}},{\mathcal{M}},\alpha,r)$.

where $K=K({\mathbf{X}},{\mathcal{M}},\alpha,r)$ is the PoSI constant.

This follows immediately from Lemma 4.1. Although mathematically trivial we give the above the status of a theorem as it is the central statement of the reduction of universal post-selection inference to simultaneous inference. The following is just a repackaging of Theorem 4.1:

where $K=K({\mathbf{X}},{\mathcal{M}},\alpha,r)$ is the PoSI constant.

Simultaneous inference provides strong family-wise error control, which in turn translates to strong error control for tests following model selection.

The proof is standard; see the online Appendix of the supplementary material [Berk et al. 2013, Section B.3]. The corollary states that, with probability $1-\alpha$, in a selected model all PoSI-significant rejections have detected true alternatives.

7 Computation of the POSI constant.

8 Scheffé protection.

It provides an upper bound for all PoSI constants:

by variable selection. The universality of the Scheffé constant is a tip-off that it may be too loose for some predictor matrices ${\mathbf{X}}$, and obtaining the sharper constant $K({\mathbf{X}})$ may be worthwhile. An indication is given by the following comparison as $r\rightarrow\infty$:

9 PoSI-sharp model selection—“SPAR.”

There exists a model selection procedure that requires the full protection of the simultaneous inference procedure (16). It is the “significance hunting” procedure that selects the model containing the most significant “effect”:

10 One primary predictor and controls—“PoSI1.”

Variable selection is limited to models that contain the primary predictor. We therefore define for any model universe ${\mathcal{M}}$ a sub-universe ${\mathcal{M}}_{j\cdot}$ of models that contain the primary predictor ${\mathbf{X}}_{j}$,

We call this the “PoSI1” situation in contrast to the unconstrained PoSI situation. Similar to PoSI, PoSI1 starts with a “significant triviality bound”:

Importantly, this constant is dominated by the general PoSI constant,

for the obvious reason that the present max-$|t|$ is smaller than the general PoSI max-$|t|$ due to ${\mathcal{M}}_{j\cdot}\subset{\mathcal{M}}$ and the restriction of inference to ${\mathbf{X}}_{j}$. The constant $K_{j\cdot}$ provides the following “PoSI1” guarantee shown as the analog of Theorem 4.1 and Corollary 4.1 folded into one:

and accordingly we have the following post-selection confidence guarantee:

Inequality (24) is immediate from Lemma 4.2. The “triviality bound” of the lemma is attained by the following variable selection procedure which we name “SPAR1”:

It is a potentially realistic description of some data analyses when a predictor of interest is determined a priori, and the goal is to optimize this predictor’s “effect.” This procedure requires the full protection of the PoSI1 constant $K_{j\cdot}$.

The structure of the PoSI problem.

We can reduce the dimensionality of the PoSI problem from $n\times p$ to $d\times p$, where $d=\operatorname{rank}(X)\leq\min(n,p)$, by introducing Scheffé’s canonical coordinates. This reduction is important both geometrically and computationally because the PoSI coverage problem really takes place in the column space of ${\mathbf{X}}$.

2 PoSI coefficient vectors in canonical coordinates.

The PoSI problem (16) is equivalent to a $d$-dimensional coverage problem for linear functions of the multivariate $t$-vector ${\mathbf{Z}}/{\hat{\sigma}}$,

3 Orthogonalities of PoSI coefficient vectors.

Orthogonalities in $\mathcal{L}({\mathbf{X}},{\mathcal{M}})$: the following statements hold assuming that the models referred to are in ${\mathcal{M}}$ (hence are of full rank).

The following vectors form an orthonormal “Gram–Schmidt” series:

Other series are obtained using $(j_{1},j_{2},\ldots,j_{d})$ in place of $(1,2,\ldots,d)$.

4 The PoSI polytope.

henceforth called the “PoSI polytope.” The PoSI coverage problem (30) is equivalent to calibrating $K$ such that

The simplest case of a PoSI polytope, for $d=p=2$, is illustrated in Figure 1 in the online Appendix of the supplementary material [Berk et al. 2013, Section B.7]. More general polytopes

Many properties of the polytopes $\bm{\Pi}_{K}$ are not specific to PoSI because they hold for polytopes (31) generated by simultaneous inference problems for linear functions with arbitrary sets $\mathcal{L}$ of unit vectors. These polytopes:

form scale families of geometrically similar bodies: $\bm{\Pi}_{K}=K\bm{\Pi}_{1}$;

are point symmetric about the origin: $\bm{\Pi}_{K}=-\bm{\Pi}_{K}$;

contain the Scheffé ball: ${\mathbf{B}}_{K}\subset\bm{\Pi}_{K}$;

are intersections of “slabs” of width $2K$:

have $2|\mathcal{L}|$ faces (assuming $\mathcal{L}\cap-\mathcal{L}=\varnothing$), and each face is tangent to the Scheffé ball ${\mathbf{B}}_{K}$ with tangency points $\pm K{\bar{\mathbf{l}}}$ $({\bar{\mathbf{l}}}\in\mathcal{L})$.

Specific to PoSI are the orthogonalities described in Proposition 5.3.

5 PoSI optimality of orthogonal designs.

Among predictor matrices with $\operatorname{rank}({\mathbf{X}})=d$ and model universes ${\mathcal{M}}$ that contain at least one maximal nested sequence of submodels, orthogonal designs with $p=d$ columns yield:

the maximal coverage probability ${\mathbf{P}}[{\mathbf{Z}}/{\hat{\sigma}}\in\bm{\Pi}_{K}]$ for fixed $K$ and

The proposition holds not only for multivariate $t$-vectors and their Gaussian limits but for arbitrary spherically symmetric distributions. Optimality of orthogonal designs translates to optimal asymptotic behavior of their constant $K({\mathbf{X}},\alpha)$ for large $d$:

Consider the Gaussian limit $r\rightarrow\infty$. For ${\mathbf{X}}$ and ${\mathcal{M}}$ as in Proposition 5.4, the asymptotic lower bound for the constant $K$ as $d\rightarrow\infty$ is attained for orthogonal designs for which the asymptotic rate is

By Proposition 5.4 the PoSI problem is bounded below by orthogonal designs, and by Proposition 4.1 it is loosely bounded above by the Scheffé ball (both for all $\alpha$, $d$, and $r$). The question of how close to the Scheffé bound PoSI problems can get for $r\rightarrow\infty$ will occupy us in Section 6.2. Unlike the infimum problem, the supremum problem does not appear to have a unique optimizing design ${\mathbf{X}}$ uniformly in $\alpha$, $d$ and $r$.

6 A duality property of PoSI vectors.

$\mathcal{L}({\mathbf{X}}^{*})=\mathcal{L}({\mathbf{X}}),\bm{\Pi}_{K}({\mathbf{X}}^{*})=\bm{\Pi}_{K}({\mathbf{X}}),K({\mathbf{X}}^{*})=K({\mathbf{X}})$.

Recall that $\mathcal{L}({\mathbf{X}})$ and $\mathcal{L}({\mathbf{X}}^{*})$ contain the normalized versions of the respective adjusted predictor vectors. The theorem follows from the following lemma which establishes the identities of vectors between $\mathcal{L}({\mathbf{X}}^{*})$ and $\mathcal{L}({\mathbf{X}})$. We extend obvious notation from ${\mathbf{X}}$ to ${\mathbf{X}}^{*}$ as follows:

A special case arises when the predictor matrix (in canonical coordinates) is chosen to be symmetric according to Proposition 5.1(7): if ${\mathbf{X}}^{T}={\mathbf{X}}$, then ${\mathbf{X}}^{*}={\mathbf{X}}({\mathbf{X}}^{T}{\mathbf{X}})^{-1}={\mathbf{X}}^{-1}$, and hence:

If ${\mathbf{X}}$ is symmetric in canonical coordinates, then

Illustrative examples and asymptotic results.

In exchangeable designs all pairs of predictor vectors enclose the same angle. In canonical coordinates a convenient parametrization of a family of symmetric exchangeable designs is

where $-1/p<a<\infty$, and ${\mathbf{E}}_{p\times p}$ is a matrix with all entries equal to $1$. The range restriction on $a$ assures that ${\mathbf{X}^{(p)}}$ is positive definite. We will write ${\mathbf{X}}={\mathbf{X}^{(p)}}={\mathbf{X}}(a)={\mathbf{X}^{(p)}}(a)$ depending on which parameter matters in a given context. We will make use of the fact that

PoSI constants of exchangeable design matrices ${\mathbf{X}^{(p)}}(a)$ [defined in (32) above] have the following limiting behavior:

The proof can be found in Appendix .3. The theorem shows that for exchangeable designs the PoSI constant remains much closer to the orthogonal case than the Scheffé case. Thus, for this family of designs it is possible to improve on the Scheffé constant by a considerable margin.

The following detail of geometry for exchangeable designs has a bearing on their PoSI constants: the angle between pairs of predictor vectors as a function of $a$ is $\cos({\mathbf{X}}_{j}^{(p)}(a),{\mathbf{X}}_{j^{\prime}}^{(p)}(a))=a(2+pa)/(pa^{2}+4a+2)$. As the vectors fall into the rank-$(p-1)$ collinearity at $a\downarrow-1/p$, the cosine becomes $-1/(2p-3)$, which converges to zero as $p\rightarrow\infty$. Thus, as $p\uparrow\infty$, exchangeable designs approach orthogonal designs even at their most collinear extreme. For further illustrative materials related to exchangeable designs, see Figures 2 and 3 in the online Appendix of the supplementary material [Berk et al. 2013, Section B.7].

2 Example 2: Where K⁡(𝐗)K({\mathbf{X}}) is close to the Scheffé bound.

For $\sigma^{2}$ known, the designs (33) have PoSI1 constants $K_{p\cdot}({\mathbf{X}^{(p)}}(c),\alpha)$ with the following asymptotic rate:

The proof is in Appendix .4. As $K({\mathbf{X}},\alpha)\geq K_{j\cdot}({\mathbf{X}},\alpha)$ the theorem provides a lower bound on the rate of the full PoSI constant. The value 0.6363… is not maximal, and we have indications that the supremum over all designs may exceed 0.78. Together with the upper bound of Corollary 6.1, this would provide a narrow asymptotic range for worst-case PoSI. Most importantly, the example shows that for some designs PoSI constants can be much larger than the $O(1)$ $|t|$-quantiles used in common practice.

3 Bounding away from Scheffé.

Denote by $\mathcal{L}_{d}$ arbitrary finite sets of $d$-dimensional unit vectors, $\mathcal{L}_{d}\subset S^{d-1}$, such that $|\mathcal{L}_{d}|\leq a_{d}$ where $a_{d}^{1/d}\rightarrow a(>1)$. Denote by $K(\mathcal{L}_{d},\alpha)$ the $(1-\alpha)$-quantile of $\sup_{{\bar{\mathbf{l}}}\in\mathcal{L}_{d}}|{\bar{\mathbf{l}}}^{T}{\mathbf{Z}}|$. Then the following describes an asymptotic worst-case bound for $K(\mathcal{L}_{d},\alpha)$ and its attainment:

The corollary shows that the asymptotic rate of the PoSI constant, if it reaches the Scheffé rate, will always have a multiplier that is strictly below that of the Scheffé constant. We do not know whether there exist designs for which the bound of the corollary is attained, but the theorem says the bound is sharp for unstructured sets $\mathcal{L}$.

Summary and discussion.

We carried out post-selection inference in a limited framework. Several problems remain open, and many natural extensions are desirable:

Computations for $p>20$ are a challenge. Straight enumeration of the set of up to $p2^{p-1}$ linear combinations should be replaced with heuristic shortcuts that yield practically useful upper bounds on $K({\mathbf{X}}_{n\times p}$, ${\mathcal{M}}$, $\alpha,r)$ that are specific to ${\mathbf{X}}$ and the set of submodels ${\mathcal{M}}$, unlike the 0.8660 fraction of the Scheffé bound which is universal.

Situations to which the PoSI framework should be extended include generalized linear models, mixed effects models, models with random predictors, as well as prediction problems. Results for the last two situations will be reported elsewhere.

It would be desirable to devise post-selection inference for specific selection procedures for cases in which a strict model selection protocol is being adhered to.

$R$ code for computing the PoSI constant for up to $p=20$ can be obtained from the authors’ web pages (a manuscript describing the computations is available from the authors).

Appendix: Proofs

.2 Proof of duality: Lemma 5.1 and Theorem 5.1.

which in turn are consequences of (36) and ${\mathbf{X}}^{T}{\mathbf{X}}^{*}={\mathbf{I}}_{p}$. We also know from (28) that

Putting together (37), (38) and (28), we obtain

Because the two vectors are scalar multiples of each other, we also know that

Putting together (39) and (40) we conclude

.3 Proof of Theorem 6.1.

The term $da$ is nonnegative, maximal wrt $m$ for $m=2$, and thereafter maximal wrt $a$ for $a=1/\sqrt{p}$, whence $\max_{a\geq 0,m\geq 2}da=1/(2(\sqrt{p}+1))$ and finally

Combining (41) and (.3) we obtain for ${\mathbf{Z}}\sim{\mathcal{N}}({\mathbf{0}},{\mathbf{I}}_{p})$ the following:

It remains to prove that equality holds in (45). To this end let $Z_{(1)}<Z_{(2)}<\cdots<Z_{(p)}$ denote the order statistics of $Z_{1}$, $Z_{2},\ldots,Z_{p}$. Fix $m$. We have in probability

Choosing $m$ arbitrarily large and combining this with (45) yields the desired conclusion.

.4 Proof of Theorem 6.2.

For fixed $m$ we can explicitly maximize the sum on the right-hand side,

where $Z_{(j)}$ is the $j$th order statistic of $Z_{1}$, $Z_{2},\ldots,Z_{p-1}$, omitting $Z_{p}$. We can also explicitly maximize the factor $c/\sqrt{1-(m-1)c^{2}}$ in (46),

and equality is attained as $c^{2}\uparrow 1/(p-1)$. Therefore, for fixed $m$, we can continue from (46) as follows:

The reason for writing the two sums in this manner is that we will interpret them as approximations to Riemann sums. To this end we borrow from Bahadur 1966 the following approximations for $j=1,\ldots,p-1$:

Reparametrizing $m=rp$, the anticipated Riemann approximation is

The function $f(r)=\frac{1}{\sqrt{1-r}}\phi(\Phi^{-1}(r))$ is maximized at $r^{*}\approx 0.72972$ with$f(r^{*})\approx 0.6363277$. Therefore,

The bound is sharp because it is attained by the models that include the first or last $m^{*}=r^{*}p$ order statistics of ${\mathbf{Z}}$ when $p\rightarrow\infty$ and $c^{2}\uparrow\frac{1}{p-1}$. From (47) we conclude that $K_{1\cdot}({\mathbf{X}})\sim 0.6363\sqrt{p}$.

.5 Proof of Theorem 6.3.

We show that if $a_{p}^{1/p}\rightarrow a(>1)$, then:

we have a uniform asymptotic worst-case bound,

which is attained when $|\mathcal{L}_{p}|=a_{p}$ and ${\bar{\mathbf{l}}}\in\mathcal{L}_{p}$ are i.i.d. $\operatorname{Unif}(S^{p-1})$ independent of ${\mathbf{Z}}$,

These facts imply the assertions about $(1-\alpha)$-quantiles $K(\mathcal{L}_{p})$ of${\max_{{\bar{\mathbf{l}}}\in\mathcal{L}_{p}}}|{\bar{\mathbf{l}}}^{T}{\mathbf{Z}}|$ in Theorem 6.3. We decompose ${\mathbf{Z}}=R{\mathbf{U}}$ where $R^{2}=\|{\mathbf{Z}}\|^{2}\sim\chi^{2}_{p}$ and ${\mathbf{U}}={\mathbf{Z}}/\|{\mathbf{Z}}\|\sim\operatorname{Unif}(S^{p-1})$ are independent. Due to $R/\sqrt{p}\stackrel{{\scriptstyle{\mathbf{P}}}}{{\rightarrow}}1$ it is sufficient to show the following:

attainment of the bound when $|\mathcal{L}_{p}|=a_{p}$ and ${\bar{\mathbf{l}}}\in\mathcal{L}_{p}$ are i.i.d. $\operatorname{Unif}(S^{p-1})$ independent of ${\mathbf{U}}$,

To show (48), we upper-bound the noncoverage probability and show that it converges to zero for $K^{\prime}>\sqrt{1-1/a^{2}}$. To this end we start with a Bonferroni-style bound, as in Wyner 1967,

where $U$ is any coordinate of ${\mathbf{U}}$ or projection of ${\mathbf{U}}$ onto a unit vector. We will show that bound (50) converges to zero. We use the fact that $U^{2}\sim\operatorname{Beta}(1/2,(p-1)/2)$, hence

We bound the Beta function and the integral separately,

where we used $\Gamma(x+1/2)/\Gamma(x)<\sqrt{x}$ (a good approximation, really) and $\Gamma(1/2)=\sqrt{\pi}$.

where we used $x^{-1/2}\leq 1/K^{\prime}$ on the integration interval. Continuing with the chain of bounds from (50) we have

Since $|\mathcal{L}_{p}|^{1/(p-1)}\rightarrow a(>0)$, the right-hand side converges to zero at geometric speed if $a\sqrt{1-K^{\prime 2}}<1$, that is, if $K^{\prime}>\sqrt{1-1/a^{2}}$. This proves (48).

show (49), we upper-bound the coverage probability and show that it converges to zero for $K^{\prime}<\sqrt{1-1/a^{2}}$. We make use of independence of ${\bar{\mathbf{l}}}\in\mathcal{L}_{p}$, as in Wyner 1967,

We will lower-bound the probability ${\mathbf{P}}[|U|>K^{\prime}]$ recalling (51) and again deal with the Beta function and the integral separately,

where we used $\Gamma(x+1)/\Gamma(x+1/2)>\sqrt{x+1/4}$ (again, a good approximation).

where we used $x^{-1/2}\geq 1$. Putting it all together we bound the exponent in (52),

Since $|\mathcal{L}_{p}|^{1/(p-1)}\rightarrow a(>0)$, the right-hand side converges to $+\infty$ at nearly geometric speed if $a\sqrt{1-K^{\prime 2}}>1$, that is, $K^{\prime}<\sqrt{1-1/a^{2}}$. This proves (49).

Acknowledgments.

We thank E. Candes, L. Dicker, M. Freiman, E. George, A. Krieger, M. Low, Z. Ma, E. Pitkin, L. Shepp, N. Sloane, P. Shaman and M. Traskin for very helpful discussions. The acronym “SPAR” is due to M. Freiman. We are indebted to an anonymous reviewer for extensive and constructive criticism that influenced the positioning of this article.

Supplement to “Valid post-selection inference” The online supplement contains the following sections:

The Full Model Interpretation of Parameters (as a contrast to the sub-model interpretation adopted in this article).

“Omitted Variables Bias” (which is not bias in the sense of this article).

Proof of Corollary 4.2 (strong error control).

Alternative PoSI Guarantees (conditional on selection).

PoSI P-Value Adjustment for Model Selection.

Figures (illustrating PoSI polytopes and results of a simulation for exchangeable designs).

References