Uniform Asymptotic Inference and the Bootstrap After Model Selection

Introduction

There has been a recent surge of work on conducting formally valid inference in a regression setting after a model selection event has occurred, see Berk et al. 2013; Lockhart et al. 2014; Tibshirani et al. 2016; Lee et al. 2016; Fithian et al. 2014; Bachoc et al. 2014, just to name a few. Our interest in this paper stems in particular from the work of Tibshirani et al. 2016, who presented a method to produce valid p-values and confidence intervals for adaptively fitted coefficients from any given step of a sequential regression procedure like forward stepwise regression (FS), least angle regression (LAR), or the lasso (the lasso is meant to be thought of as tracing out a sequence of models along its solution path, as the penalty parameter descends from $\lambda=\infty$ to $\lambda=0$). These authors use a statistic that is carefully crafted to be pivotal after conditioning on the model selection event. This idea is not specific to the sequential regression setting, and is an example of a broader framework that we might call selective pivotal inference, applicable in many other settings, as in, e.g., Taylor et al. 2016; Lee et al. 2016; Lee & Taylor 2014; Loftus & Taylor 2014; Reid et al. 2017; Choi et al. 2014; Fithian et al. 2014; Hyun et al. 2016.

A high-level description of the selective pivotal inference framework for sequential regression is as follows (details are provided in Section 2). FS, LAR, or the lasso is run for some number of steps $k$, and a model is selected, call it $M$. For FS and LAR, this model will always have $k$ active variables, and for the lasso, it will have at most $k$, as variables can be added to or deleted from the active set at each step. We specify a linear contrast of the mean $v^{T}\theta$ of interest, e.g., one giving the coefficient of a variable of interest in the model $M$ at step $k$, in the regression of $\theta$ onto the active variables. By assuming normal errors in (1), and examining the distribution of $v^{T}Y$ conditional on having selected model $M$, which we denote by $\widehat{M}(Y)=M$, we can construct a confidence interval $C_{\alpha}$ satisfying

for a given $\alpha\in$. The interpretation: if we were to repeatedly draw $Y$ from (1) and run FS, LAR, or the lasso for $k$ steps, and only pay attention to cases in which we selected model $M$, then among these cases, the constructed intervals $C_{\alpha}=C_{\alpha}(Y;M)$ contain $v^{T}\theta$ with frequency tending to $1-\alpha$.

The above is a conditional perspective of the selective pivotal inference framework for FS, LAR, and lasso. An unconditional or marginal point of view is also possible, which we now describe. For each possible selected model $M$, a constrast vector $v_{M}$ is specified, and the contrast $v_{M}^{T}\theta$ is considered when model $M$ is selected, $\widehat{M}(Y)=M$. To be concrete, we can again think of a setup such that $v_{M}^{T}\theta$ gives the coefficient of a variable in the model $M$ at step $k$, in the projection of $\theta$ onto the active set. Confidence intervals are then constructed in exactly the same manner as above (without change), and conditional coverage over all models $M$ implies the following unconditional property for $C_{\alpha}$,

The interpretation is different: if we were to repeatedly draw $Y$ from (1) and run FS, LAR, or lasso for $k$ steps, and construct confidence intervals $C_{\alpha}=C_{\alpha}(Y;\widehat{M}(Y))$, then these intervals contain their respective targets $v_{\widehat{M}(Y)}^{T}\theta$ with frequency approaching $1-\alpha$. Notice that, by construction, the target itself may change each time we draw $Y$, though it is the same for all $Y$ that give rise to the same selected model. In terms of the setting for regression contrasts described above, each time we draw $Y$ and carry out the inferential procedure, the interval $C_{\alpha}$ covers the coefficient of a possibly different variable in the active model, in the projection of $\theta$ onto the active variables. Figure 1 demonstrates this point.

(The above inequalities, as in $W_{n}\leq x$ and $W\leq x$, are meant to be interpreted componentwise; we are also implicitly assuming that the limiting distribution $G$ is continuous, otherwise the above inner supremum should be restricted to continuity points $x$ of $G$.) This is much stronger than the notion of pointwise convergence in distribution, which only requires that

for a particular sequence of distributions $F_{n}$, $n=1,2,3,\ldots$.

A recent article by Kasy 2015 emphasizes the importance of uniformity in asymptotic approximations. This authors points out that a uniform version of the continuous mapping theorem follows directly from a standard proof of the continuous mapping theorem (e.g., see Theorem 2.3 in van der Vaart 1998).

Kasy 2015 also remarks that the central limit theorem for triangular arrays, specifically the Lindeberg-Feller central limit theorem (e.g., Proposition 2.27 in van der Vaart 1998) naturally extends to the uniform case. The logic is, roughly speaking: uniform convergence in (2) is equivalent to pointwise convergence over all sequences of distributions $F_{n}$, $n=1,2,3,\ldots$, and triangular arrays, by design, can have a different distribution assigned to each row. Therefore if the Lindeberg condition holds for any possible sequence, then so does the convergence to normality.

where $\Sigma$ does not depend on the sequence $F_{n}$, $n=1,2,3,\ldots$. Then $W_{n}=\sum_{i=1}^{n}\xi_{i}$ converges in distribution to $W\sim N(0,\Sigma)$, uniformly with respect to $\pazocal{P}_{n}$.

In our work, a motivating reason for the study of uniform convergence is the associated property of uniform validity of asymptotic confidence intervals. That is, if $W_{n}=W_{n}(\mu)$ depends on a parameter $\mu=\mu(F_{n})$ of the distribution $F_{n}$, but $W$ does not, then we can consider any $(1-\alpha)$ confidence set $C_{n,\alpha}$ built from a $(1-\alpha)$ probability rectangle $R_{\alpha}$ of $W$,

and the uniform convergence of $W_{n}$ to $W$, really just by rearranging its definition in (2), implies

Meanwhile, pointwise convergence as in (3) only implies

for a particular sequence $F_{n}$, $n=1,2,3,\ldots$. For a confidence set satisfying (4), and a given tolerance $\epsilon>0$, there exists a sample size $n(\epsilon)$ such that the coverage is guaranteed to be at least $1-\alpha-\epsilon$, for $n\geq n(\epsilon)$, no matter the underlying distribution (over the class of distributions in question). Note that this is not necessarily true for a pointwise confidence set as in (5), as the required sample size here could depend on the particular distribution under consideration.

2 Summary of main results

An overview of our main contributions is as follows.

We establish that TG statistics for typical inferences along the FS, LAR, and lasso paths only depend on the data $(X,Y)$ through $\frac{1}{n}X^{T}X$ and $\frac{1}{\sqrt{n}}X^{T}Y$ (Lemmas 3, 4, and 5 in Section 3), which is important since these two quantities have asymptotic limits in a standard low-dimensional asymptotic setup.

Placing mild constraints on the mean and error distribution in (1), and treating the dimension $d$ as fixed, we prove that the TG test statistic is asymptotically pivotal, converging to $U(0,1)$ (the standard uniform distribution), when evaluated at the true population value for its pivot argument. We show that this holds uniformly over a wide class of distributions for the errors, without any real restrictions on the predictors $X$ (first part of Theorem 7 in Section 4).

The resulting confidence intervals are therefore asymptotically uniformly valid, over the same class of distributions (second part of Theorem 7 in Section 4).

The above asymptotic results assume that the error variance $\sigma^{2}$ is known, so for $\sigma^{2}$ unknown, we propose a plug-in approach that replaces $\sigma^{2}$ in the TG statistic with a simple estimate, and alternatively, an efficient bootstrap approach. Both allow for conservative asymptotic inference (Theorem 11 in Section 5).

We present detailed numerical experiments that support the asymptotic validity of the TG p-values and confidence intervals for inference in low-dimensional regression problems that have nonnormal errors (Section 6). Our experiments reveal that the plug-in and bootstrap versions also show good performance, and the bootstrap method can often deliver substantially shorter intervals than those based directly on the TG statistic.

Our experiments also also suggest that the TG test statistic (and plug-in, bootstrap variants) may be asymptotically valid in even broader settings not covered by our theory, e.g., problems with heteroskedastic errors and (some) high-dimensional problems.

We prove that TG statistic does not exhibit a general uniform convergence to $U(0,1)$ when the dimension $d$ is allowed to increase (Theorem 12 in Section 7).

3 Related work

A recent paper by Tian & Taylor 2017 is very related to our work here. These authors examine the asymptotic distribution of the TG statistic under nonnormal errors. Their main result proves that the TG statistic is asymptotically pivotal, under some restrictions on the model selection events in question. We view their work as providing a complementary perspective to our own: they consider a setting where the dimension $d$ grows, but place strong regularity conditions on the selected models; we adopt a more basic setting with $d$ fixed, and prove more broad uniformly valid convergence results for the TG pivot, free of regularity conditions.

In a sequence of papers, Leeb & Potscher 2003; Leeb & Potscher 2006; Leeb & Potscher 2008 prove that in a classical regression setting, it is impossible to find the distribution of a post-selection estimator of the underlying coefficients, even asymptotically. Specifically, they prove for an estimate $\widehat{\beta}$ of some underlying coefficient vector $\beta_{0}$, any quantity of the form $Q_{n}=\sqrt{n}A(\widehat{\beta}-\beta_{0})$, for a linear transform $A$, cannot be used for inference after model selection. Though $Q_{n}$ can be made to be pivotal or at least asymptotically pivotal (once $A$ is chosen once appropriately), this is no longer true in the presence of selection, even if the dimension $d$ is fixed and the sample size $n$ approaches $\infty$. Furthermore, they show that there is no uniformly consistent estimate of the distribution of $Q_{n}$ (either conditionally or unconditionally), which makes $Q_{n}$ unsuitable for inference. This fact is essentially a manifestation of the well-known Hodges phenomenon. The selective pivotal inference framework, and hence our paper, circumvents this problem as we do not claim (nor attempt) to estimate the distribution of $Q_{n}$, and instead make inferences using an entirely different pivot that is constructed via a careful conditioning scheme.

4 Notation

Selective inference

In this section, we review the selective pivotal inference framework for sequential regression procedures. We present interpretations for the inferences from both conditional and unconditional perpsectives, in Sections 2.2 and 2.5, respectively. The other subsections provide the necessary details for understanding the framework, beginning with the selection events encountered along the FS, LAR, and lasso paths.

The active sets are nested across steps, $\widehat{A}_{1}(y)\subseteq\widehat{A}_{2}(y)\subseteq\widehat{A}_{3}(y)\subseteq\ldots$, as FS selects one variable to add to the active set at each step. However, the sign vectors $\widehat{s}_{1}(y),\widehat{s}_{2}(y),\widehat{s}_{3}(y),\ldots$ are not, since these are determined by least squares on the active variables at each step. Hence, as defined, the number of possible models $\widehat{M}(y)$ after $k$ steps of FS is

Moreover, the corresponding partition elements $\Pi_{M}$, $M\in\pazocal{M}$ in (6) are all convex cones. The proof of this fact is not difficult, and requires only a slight modification of the arguments in Tibshirani et al. 2016, given in Appendix A.1 for completeness. The result is easily seen for $k=1$: after one step of FS, assuming without a loss of generality that $X_{1},\ldots,X_{d}$ have unit norm, we can express, e.g.,

2 Inference after selection

and therefore $H_{0}:v^{T}\theta=0$ is a test for the significance of the $j$th normalized coefficient in the linear projection of $\theta$ onto $X_{A}$, written as $\beta_{j}(A)$ for short. (Though the normalization in the denominator is irrelevant for this significance test, it acts as a key scaling factor for the asymptotics in Section 4.) The idea of using a projection parameter for inference, $\beta_{j}(A)$, has also appeared in, e.g., Berk et al. 2013; Wasserman 2014; Lee et al. 2016. Here is now a summary of the testing framework.

Assume i.i.d. $N(0,\sigma^{2})$ errors in (1). Under the null hypothesis, the TG statistic has a standard uniform distribution, over draws of $Y$ that land in $\Pi_{M}$. Mathematically, this is the property

for all $t\in$. The probability above is taken over an arbitrary mean parameter $\theta$ for which $v^{T}\theta=\mu$ (in fact, the TG statistic is constructed so that the law of $T(Y;M,v,\mu)\,|\,\widehat{M}(Y)=M$ only depends on $\theta$ through $v^{T}\theta$, so this is unambiguous). In order for (8) to hold, of course, $v$ and $\mu$ cannot be random, i.e., they cannot depend on $Y$, though they can be functions of $M$.

Thus $T(Y;M,v,\mu)$ serves as a valid p-value (with exact finite sample size) for testing the null hypothesis $H_{0}:v^{T}\theta=\mu$, conditional on $\widehat{M}(Y)=M$.

A confidence interval is obtained by inverting the test in (8). Given a desired confidence level $1-\alpha$, we define $C_{\alpha}$ to be the set of all values $\mu$ such that $\alpha/2\leq T(Y;M,v,\mu)\leq 1-\alpha/2$. Then, by construction, the property in (8) (which we reiterate, assumes i.i.d. $N(0,\sigma^{2})$ errors) translates into

The interpretation of the above statement is straightforward: the random interval $C_{\alpha}$ contains the fixed parameter $v^{T}\theta$ with probability $1-\alpha$, conditional on $\widehat{M}(Y)=M$.

3 The truncated Gaussian pivot

We now describe the truncated Gaussian (TG) pivotal quantity in detail. As defined in Section 2.1, if we write $\widehat{M}(y)$ for the selected model from the given algorithm (FS, LAR, or lasso), run for $k$ steps on $y$, then $\Pi_{M}=\{y:\widehat{M}(y)=M\}$ is a convex cone, for any fixed achieveable model $M$. Hence

for a fixed matrix $Q_{M}$ (here the inequality is meant to be interpreted componentwise). Now to define the pivot $T(\,\cdot\,;M,v,\mu)$ for testing $H_{0}:v^{T}\theta=\mu$, several preliminary quantities must be introduced:

This has the following property, as stated in (8): when $Y$ is drawn from (1) with i.i.d. $N(0,\sigma^{2})$ errors, and $v^{T}\theta=\mu$, the pivot $T(Y;M,v,\mu)$ is uniformly distributed conditional on $\widehat{M}(Y)=M$. See Lemmas 1 and 2 in Tibshirani et al. 2016 for a proof of this result.

4 P-values and confidence intervals

A statistic aligned to have power against the two-sided alternative $H_{1}:v^{T}\theta\not=0$ is simply given by $2\min\{T(Y;M,v,0),\,1-T(Y;M,v,0)\}$. For purely testing purposes, we find the one-sided p-values discussed above to be more natural, and hence these will serve as our default. On the other hand, for constructing confidence intervals, we prefer to invert the two-sided statistics, since these lead to two-sided intervals. As

the previously described confidence interval in (9) is just given by inverting the two-sided pivot.

To summarize: the default in this work, as with Tibshirani et al. 2016, is to consider one-sided hypothesis tests, but two-sided intervals. These are just two slightly different uses of the same pivot.

5 Inference after selection, revisited

We have portrayed selective pivotal inference, in sequential regression procedures, as a method for producing conditional p-values and intervals. An unconditional interpretation of this framework is also possible, which we describe here.

where $1_{\Pi_{M}}(\,\cdot\,)$ denotes the indicator function for the partition element $\Pi_{M}$ (and $T(\,\cdot\,;M,v_{M},\mu_{M})$ is as before, defined in (10)). The unconditional statistic can be used as follows: if a response $Y$ is drawn from (1), then we can form $\pazocal{T}(Y;V,U)=T(Y;\widehat{M}(Y),v_{\widehat{M}(y)},\mu_{\widehat{M}(y)})$ to test the hypothesis $H_{0}:v_{\widehat{M}(Y)}^{T}\theta=\mu_{\widehat{M}(Y)}$.

Assume that the errors in (1) are i.i.d. $N(0,\sigma^{2})$. Then under the proper hypothesis, by summing up the conditional property in (8) across partition elements, we have

for all $t\in$. The assertion above holds for a parameter $\theta$ such that $V^{T}\theta=U$, which we use as shorthand for $v_{M}^{T}\theta=\mu_{M}$ for all $M\in\pazocal{M}$. Note that this full specification, across all $M\in\pazocal{M}$, is critical in order to apply the relevant null probability within each partition element (giving rise to the equality in (12)).

Therefore $\pazocal{T}(Y;V,U)$ serves as a valid p-value (with exact finite sample size)—but for testing what null hypothesis? Formally, it is attached to $H_{0}:V^{T}\theta=U$, an exhaustive specification of $v_{M}^{T}\theta=\mu_{M}$, over all $M\in\pazocal{M}$, but in truth, $\pazocal{T}(Y;V,U)$ carries no information about models other than the selected one, $\widehat{M}(Y)$. For this reason, we actually consider $\pazocal{T}(Y;V,U)$ to be a p-value for the random null hypothesis $H_{0}:v_{\widehat{M}(Y)}^{T}\theta=\mu_{\widehat{M}(Y)}$. This is made more precise through confidence intervals.

A confidence interval is obtained by inverting the test in (12). But the TG statistic at $Y$,

only depends on $U$ through $\mu_{\widehat{M}(Y)}$. Thus, given a desired confidence level $1-\alpha$, let us define $D_{\alpha}$ to be the set of $U$ such that $\alpha/2\leq\pazocal{T}(Y;V,U)\leq 1-\alpha/2$, and $C_{\alpha}$ to be the set of $\mu_{\widehat{M}(Y)}$ such that $\alpha/2\leq T(Y;\widehat{M}(Y),v_{\widehat{M}(Y)},\mu_{\widehat{M}(Y)})\leq 1-\alpha/2$. Then we can see that

so the confidence interval is effectively infinite with respect to the values $\mu_{M}$, $M\not=\widehat{M}(Y)$, and inverting the test in (12) yields

The above expression says that the random interval $C_{\alpha}$ traps the random parameter $v_{\widehat{M}(Y)}^{T}\theta$ with probability $1-\alpha$, and thus, this supports the interpretation of $H_{0}:v_{\widehat{M}(Y)}^{T}\theta=\mu_{\widehat{M}(Y)}$ as the null hypothesis underlying the unconditional TG statistic.

The pivotal property in (12) is derived under the distributional assumption that $V^{T}\theta=U$, i.e., $v_{M}^{T}\theta=\mu_{M}$ for all $M\in\pazocal{M}$, which may seem unnatural, as the catalog $U$ of pivot value can be large (e.g., on the order of $d^{k}$ after $k$ steps of FS), and so this is condition on possibly many contrasts of $\theta$. However, it is worth emphasizing that the unconditional testing property in (12) is really only useful in that it allows us to formulate the unconditional confidence interval property in (13), which is a more natural statement about coverage of a single (random) parameter. When viewing selective inference from an unconditional perpsective, we find it more natural to place the focus on confidence intervals rather than hypothesis testing; in many ways, we find the former the more natural of the two perspectives, unconditionally. Tibshirani et al. 2016 in fact suggest separate nomenclature for the unconditional case, referring to the property in (13) as that of a selection interval (rather than confidence interval), to emphasize that this interval covers a moving target.

The master statistic

Given a response $y$ and predictors $X$, our description thus far of the selected model $\widehat{M}(y)$, statistics $T(y;M,v,\mu)$ and $\pazocal{T}(y;V,U)$, etc., has ignored the role of $X$. This was done for simplicity. The theory to come in Section 4 will consider $X$ to be nonrandom, but asymptotically $X$ must (of course) grow with $n$, and so it will help to be precise about the dependence of the selected model and statistics on $X$. We will denote these quantities by $\widehat{M}(X,y)$, $T(X,y;M,v,\mu)$, and $\pazocal{T}(X,y;V,U)$ to emphasize this dependence. We define

a $d(d+3)/2$-dimensional quantity that we will call the master statistic. As its name might suggest, this plays an important role: all normalized coefficients from regressing $y$ onto subsets of the variables $X$ can be written in terms of $\Omega_{n}$. That is, for an arbitrary set $A\subseteq\{1,\ldots,p\}$, the $j$th normalized coefficient from the regression of $y$ onto $X_{A}$ is

which only depends on $(X,y)$ through $\Omega_{n}$. The same dependence is true, it turns out, for the selected models from FS, LAR, and the lasso. We defer the proof of the next lemma, as with all proofs in this paper, until the appendix.

For each of the FS, LAR, and lasso procedures, run for $k$ steps on data $(X,y)$, the selected model $\widehat{M}(X,y)$ only depends on $(X,y)$ through $\Omega_{n}=(\frac{1}{n}X^{T}X,\frac{1}{\sqrt{n}}X^{T}y)$, the master statistic.

In more detail, for any fixed $M\in\pazocal{M}$, the matrix $Q_{M}(X)$ such that $\widehat{M}(X,y)=M\iff Q_{M}(X)\,y\geq 0$ can be written as $Q_{M}(X)=P_{M}(\frac{1}{n}X^{T}X)\,\frac{1}{\sqrt{n}}X^{T}$, where $P_{M}$ depends only on $\frac{1}{n}X^{T}X$. Hence

This lemma asserts that the master statistic governs model selection, as performed by FS, LAR, and the lasso. It is also central to TG pivot for these procedures. Denoting $M=\widehat{M}(X,y)$, the statistic $T(X,y;M,v,\mu)$ in (10) only depends on $(X,y)$ through three quantities:

The third quantity is always a function of $\Omega_{n}$, by Lemma 3. When $v$ is chosen so that $v^{T}y$ is a normalized coefficient in the regression of $y$ onto a subset of the variables in $X$, the first two quantities are also functions of $\Omega_{n}$. Thus, in this case, the TG pivot only depends on $(X,y)$ through the master statistic $\Omega_{n}$; in fact, it is continuous at any point such that $\frac{1}{n}X^{T}X$ is nonsingular and $y$ does not lie on the boundary of a model selection event.

Fix any model $M\in\pazocal{M}$, and suppose that $v$ is chosen so that $v^{T}y$ is a normalized coefficient from projecting $y$ onto a subset of the variables in $X$. Then the TG statistic only depends on $(X,y)$ by means of $\Omega_{n}$, so that we may write

Further, the function $\psi_{M}$ is continous at any point $(S,z)$ such that $S$ is nonsingular and $P_{M}(S)\,z>0$.

Finally, we show that the conditional pivotal property of the TG statistic in (8) can be phrased entirely in terms of the master statistic.

Assume the conditions of Lemma 4, and additionally that $Y$ is drawn from (1). Construct the master statistic $\Omega_{n}=(\frac{1}{n}X^{T}X,\frac{1}{\sqrt{n}}X^{T}Y)$. Then there is a function $g$ such that

Thus if the errors in (1) are i.i.d. $N(0,\sigma^{2})$, then the conditional pivotal property (8) of the TG statistic can be reexpressed as

Equipped with the last two lemmas, asymptotic theory for the TG test, when $d$ is fixed, is not far off. Under weak conditions on the data model in (1), the central limit theorem tells us that $\frac{1}{\sqrt{n}}X^{T}Y$ converges weakly to a normal random variable. With $\frac{1}{n}X^{T}X$ converging to a deterministic matrix, the continuous mapping theorem will then provide the appropriate asymptotic limit for the statistic $T(X,y;M,v,\mu)=\psi_{M}(\frac{1}{n}X^{T}X,\frac{1}{\sqrt{n}}X^{T}Y)$. This is made more precise next.

Asymptotic theory

We specify the class of distributions that we will be working with for $Y$ in (1). Let $\sigma^{2}>0$ be a fixed, known constant. First we define a set of error distributions

As $n$ grows, we allow the underlying mean $\theta$ to change, but we place a restriction on this parameter so that it has an appropriate asymptotic limit. Specifically, we consider a class $\Theta$ of sequences of mean parameters such that $\frac{1}{\sqrt{n}}X^{T}\theta$ has an asymptotic limit lying in some compact set, with uniform convergence to this limit. Formally, write (in a slight abuse of notation) $\theta\in\Theta$ to denote a sequence of mean parameters in $\Theta$, and let $E(\Theta)$ denote the set of limit points of $\{\frac{1}{\sqrt{n}}X^{T}\theta:\theta\in\Theta\}$. Then, for some constant $B>0$, we require of the class $\Theta$,

2 Uniform convergence results

We begin with a result on the uniform convergence of (the random part of) the master statistic to a normal distribution, both marginally and conditionally.

Assume that $X$ has asymptotic covariance matrix $\Sigma$, as in (14), and satisfies the normalization condition in (15). Let $Y\sim F_{n}(\theta)\in\pazocal{P}_{n}(\theta)$, this class as defined in (16), for a sequence of mean parameters $\theta\in\Theta$, as defined in (17). Denote $\frac{1}{\sqrt{n}}X^{T}\theta\to\eta$ as $n\to\infty$. Then $Z_{n}=\frac{1}{\sqrt{n}}X^{T}Y$ converges in distribution to $Z\sim N(\eta,\sigma^{2}\Sigma)$, uniformly over $\pazocal{P}_{n}(\theta)$, and uniformly over all $\theta\in\Theta$. That is,

This lemma, combined with Lemmas 4 and 5 of the last section, leads us to uniform asymptotic theory for the TG test. We remind the reader that $k$, the number of steps, is to be considered fixed in the next result (as it is throughout the paper).

Assume the conditions of Lemma 6. Suppose FS, LAR, or the lasso is run for $k$ steps on $(X,Y)$. Below we describe the conditional and unconditional asymptotic results separately.

(a, Markovic) Fix any model $M\in\pazocal{M}$. Let $v$ be a vector such that $v^{T}\theta$ gives a normalized coefficient in the projection of $\theta$ onto some subset of the variables in $X$, and let $\mu$ be an arbitrary pivot value. Then under $v^{T}\theta=\mu$, the conditional TG statistic $T(X,Y;M,v,\mu)\,|\,\widehat{M}(X,Y)=M$ converges in distribution to $W\sim U(0,1)$, uniformly over $\pazocal{P}_{n}(\theta)$, and over $\theta\in\Theta$. That is,

Moreover, if we define $C_{n,\alpha}$ to be the set of $\mu$ such that $\alpha/2\leq T(X,Y;M,v,\mu)\leq 1-\alpha/2$, then $C_{n,\alpha}$ is an asymptotically uniformly valid confidence interval for $v^{T}\theta$. That is,

(b) Let $V=\{v_{M}:M\in\pazocal{M}\}$ be a catalog of vectors such that each $v_{M}^{T}\theta$ yields a normalized coefficient in the projection of $\theta$ onto a subset of the variables in $X$, for $M\in\pazocal{M}$, and $U=\{\mu_{M}:M\in\pazocal{M}\}$ be a catalog of pivot values. Then under $V^{T}\theta=U$, the same results as in part (a) hold marginally. That is,

and for $C_{n,\alpha}$ defined to be the set of $\mu$ such that $\alpha/2\leq T(X,Y;\widehat{M}(X,Y),v_{\widehat{M}(X,Y)},\mu)\leq 1-\alpha/2$,

An initial version of this work contained only the unconditional result in part (b) of the theorem. Jelena Markovic pointed out that the conditional result in part (a) should also be possible, and thus this conditional result should also be attributed to her. Between the initial and the current version of this paper, in addition to revising Theorem 7, we have also revised Theorems 11 and 12 to include the appropriate conditional results.

Unknown σ2\sigma^{2} and the bootstrap

The results of the previous section assumed that the error variance $\sigma^{2}$ in the model (1) was known. Here we consider two strategies when $\sigma^{2}$ is unknown. The first plugs a (rather naive) estimate of $\sigma^{2}$ into the usual TG statistic. The second is a computationally efficient bootstrap method. Both, as we will show, yield asymptotically conservative p-values. (In practice, the bootstrap often gives shorter confidence intervals than those based on the TG pivot; see Section 6.)

Given a model $M\in\pazocal{M}$, contrast vector $v$, and pivot value $\mu$, consider the TG statistic $T(X,Y;M,v,\mu)$. Let us abbreviate

where the latter two functions are as defined in Section 2.3. In this notation, we can succintly write the TG statistic as

When $\sigma^{2}$ is unknown, we propose a simple plug-in approach that replaces $\sigma$ with $cs_{Y}$, where

the sample variance of $Y$ (here $\overline{Y}=\sum_{i=1}^{n}Y_{i}/n$ denotes the sample mean), and $c>1$ is a fixed constant. To be explicit, we consider the modified TG statistic

The scaling factor $c$ facilitates our theoretical study of the above plug-in statistic, and practically, we have found that ignoring it (i.e., setting $c=1$) works perfectly well, though a choice of, say, $c=1.0001$ seems to have a minor effect anyway.

When the mean $\theta$ of $Y$ is nonzero, the sample variance $s_{Y}^{2}$ is generally too large as an estimate of $\sigma^{2}$. As we will show, the modified statistic in (19) thus yields asymptotically conservative p-values. Residual based estimates of $\sigma^{2}$ are not as useful in our setting because they depend more heavily on the linearity of the underlying regression model, and they suffer practically when $d$ is close to $n$ (see also the discussion at the start of Section 6).

2 An efficient bootstrap approach

As an alternative to the plug-in method of the last subsection, we investigate a highly efficient bootstrap scheme that does not rely on knowledge of $\sigma^{2}$. Our general framework so far treats $X$ as fixed, and for our bootstrap strategy to respect this assumption, we cannot use, say, the pairs bootstrap, and must perform sampling with respect to $Y$ only. The residual bootstrap is ruled out since we do not assume that the mean $\theta$ follows a linear model in $X$. This leaves us to consider simple bootstrap sampling of the components of $Y$. This is somewhat nonstandard, as the components of $Y$ in (1) are not i.i.d., but it provides a mechanism for provably conservative asymptotic inference, and it is what makes our approach so computationally efficient.

where the probability on the right-hand side is taken with $Y$ (and thus $\widehat{a}_{M},\widehat{b}_{M}$) treated as fixed, and with $Z_{\mu,\sigma^{2}}$ denoting a $N(\mu,\sigma^{2})$ random variable. The main idea is now to approximate the truncated normal distribution underlying the TG statistic with an appropriate one from bootstrap samples,

where $c>1$ is a constant as before, and $\delta_{n}=\gamma n^{-1/4}$ for a small constant $\gamma>0$. Again, we have found that ignoring the scaling factor $c$ (i.e., setting $c=1$) works just fine in practice, though a choice like $c=1.0001$ does not cause major differences anyway. On the contrary, a nonzero choice of the padding factor like $\delta_{n}=10^{-4}n^{-1/4}$ does play an important practical role, since the bootstrap probabilities in the numerator and denominator in (20) can sometimes be zero.

Lastly, it is worth emphasizing that practical estimation of the bootstrap probabilities appearing in (20) is quite an easy computational task, because the regression procedure in question, be it FS, LAR, or the lasso, need not be rerun beyond its initial run on the observed $Y$. After this initial run, we can just save the realized quantities $\widehat{a}_{M},\widehat{b}_{M}$, and then draw, say, $B=1000$ bootstrap samples $Y^{*}$ in order to estimate the probabilities in (20). This is not at all computationally expensive. Moreover, to estimate (20) over multiple trial values of $\mu$ (so that we can invert these bootstrap p-values for a bootstrap confidence interval), only a single common set of bootstrap samples is needed, since we can just shift $v^{T}Y^{*}$ appropriately for each bootstrap sample $Y^{*}$.

3 Asymptotic theory for unknown σ2\sigma^{2}

Treating the dimension $d$ as fixed, we will assume the previous limiting conditions (14), (15) on the matrix $X$, and additionally, that

Assume that $X$ satisfies (14), (15), (21). If $v$ is any vector such that $v^{T}\theta$ gives a normalized regression coefficient from projecting $\theta$ onto some subset of the variables in $X$, then

We specify assumptions on the distribution of $Y$ in (1) that are similar to (but slightly stronger than) those in Section 4.1. For constants $\sigma^{2},\tau,\kappa>0$, we define a set of error distributions

where as before, $F_{\mu}$ denotes the distribution of $\mu+\delta$, for $\delta\sim F$. We define a class $\Theta^{\prime}$ of sequences of mean parameters that satisfies, as before,

for a constant $B>0$, where recall $E(\Theta^{\prime})$ denotes the set of limit points in $\Theta^{\prime}$; also, for each $\theta\in\Theta^{\prime}$, at each $n$, we require

for constants $S,R>0$, where $\overline{\theta}=\sum_{i=1}^{n}\theta_{i}/n$. Note that the assumptions $Y\sim F_{n}(\theta)$, with $F_{n}(\theta)\in\pazocal{P}^{\prime}_{n}(\theta)$ and $\theta\in\Theta^{\prime}$, are not much stronger than our assumptions in Section 4.1: we require the existence of two more moments for the error distribution, and place an additional weak condition on the growth of (components of) $\theta$. These conditions are sufficient to prove the following helpful lemma.

Assume that $X$ satisfies (14), (15). Let $Y\sim F_{n}(\theta)\in\pazocal{P}^{\prime}_{n}(\theta)$, where this class is as defined in (22), and let $\theta\in\Theta^{\prime}$, where this class is as in (23), (24). Then for any fixed $M\in\pazocal{M}$, and $c>1$,

In words, the event $\{cs_{Y}\geq\sigma\}$ has probability tending to 1 conditional on $\widehat{M}(X,Y)=M$, uniformly over $\pazocal{P}^{\prime}_{n}(\theta)$, and over $\theta\in\Theta^{\prime}$. Furthermore, denoting the sample third moment of $Y$ as

we have that for any $\delta>0$, there exists $C>0$ such that for sufficiently large $n$,

The last two lemmas allow us to tie the distribution function of our bootstrap contrast to that of a normal random variable.

Assume that $X$ satisfies (14), (15), (21). Let $Y\sim F_{n}(\theta)\in\pazocal{P}^{\prime}_{n}(\theta)$, as defined in (22), and let $\theta\in\Theta^{\prime}$, as defined in (23), (24). Let $M\in\pazocal{M}$, and let $v$ be such that $v^{T}\theta$ gives a normalized regression coefficient from projecting $\theta$ onto a subset of the variables in $X$. Then for any $\delta>0$, there exists $C>0$ such that sufficiently large $n$,

We are now ready to present uniform asymptotic results for the plug-in and bootstrap TG statistics. We remind the reader the number of steps $k$ is treated as fixed below (as it is throughout).

Assume the conditions of Lemma 10. Suppose FS, LAR, or the lasso is run for $k$ steps on $(X,Y)$. Then under $v^{T}\theta=0$, the conditional plug-in TG statistic $\widetilde{T}(X,Y;M,v,0)\,|\,\widehat{M}(X,Y)=M$ and conditional bootstrap TG statistic $T^{*}(X,Y;M,v,0)\,|\,\widehat{M}(X,Y)=M$ are each asymptotically larger than $U(0,1)$ in distribution, uniformly over $\pazocal{P}^{\prime}_{n}(\theta)$, and over $\theta\in\Theta^{\prime}$. That is,

where $x_{+}=\max\{x,0\}$ denotes the positive part of $x$. Further, given any catalog $V=\{\mu_{M}:M\in\pazocal{M}\}$ of vectors such that each $v_{M}^{T}\theta$ yields a normalized coefficient in the projection of $\theta$ onto a subset of the variables in $X$, for $M\in\pazocal{M}$, the same results hold marginally under $V^{T}\theta=0$.

For simplicity, we analyzed the plug-in and bootstrap statistics simultaneously. Consequently, the conditions assumed to prove asymptotic properties of the plug-in approach are stronger than what we would need if we were to study this method on its own, but there are not major differences in these conditions.

Theorem 11 establishes that the plug-in and bootstrap versions of the TG statistic are asymptotically conservative when viewed as p-values under $v^{T}\theta=0$. If we look more broadly at the distribution of these test statistics under $v^{T}\theta=\mu$, for an arbitrary value of $\mu$, then a technical barrier arises. For each statistic, our proof of its asymptotic conservativeness leverages the fact that the truncated Gaussian survival function decreases (in a pointwise sense), as its underlying variance parameter decreases. To extend these results to the case of an arbitrary pivot value $\mu$, we would need the analogous fact to hold when we replace the survival function of the Gaussian variate $cs_{Y}Z+\mu$ truncated to $[\widehat{a}_{M},\widehat{b}_{M}]$, with that of $\sigma Z+\mu$ tuncated to $[\widehat{a}_{M},\widehat{b}_{M}]$, on the event $\{cs_{Y}\geq\sigma\}$. Yet, without the guarantee that $\widehat{a}_{M}\geq\mu$ (which clearly cannot always be true, for an arbitrary value of $\mu$), it is no longer the case that decreasing the variance from $c^{2}s_{Y}^{2}$ to $\sigma^{2}$ always decreases the survival functions of these two truncated Gaussians; see Appendix A.11. This means that confidence intervals given by directly inverting either the plug-in or bootstrap TG statistic do not have provably correct asymptotic coverage properties, under the current analysis.

From the arguments in the proof of Theorem 11, we can construct one-sided confidence intervals with conversative asymptotic coverage, by forcing them to include $\widehat{a}_{M}$. We do not pursue the details here, as we have found that these one-sided intervals are practically too wide to be of interest.

Importantly, the plug-in and bootstrap TG statistics often display excellent empirical properties, as we will show in the next section. A more refined analysis is needed to establish asymptotic uniformity for the distribution of these statistics under $v^{T}\theta=\mu$. Such asymptotic uniformity, for arbitrary $\mu$, would lead to asymptotic coverage guarantees for confidence intervals produced by inverting these statistics, and we leave this extension to future work.

Examples

We present empirical examples that support the theory developed in the previous sections, and also suggest that there is much room to refine and expand our current set of results. The first two subsections examine a low-dimensional problem setting that is covered by our theory. The last two look at substantial departures from this theoretical framework, the heteroskedastic and high-dimensional settings, respectively. In all examples, the LAR algorithm was used for variable selection and associated inferences; results with the FS and lasso paths were roughly similar. Also, in all examples, where not explicitly stated otherwise, the computed p-values are a test of whether the target population value is 0.

It may be worth discussing two potentially common reactions to our experimental setups, especially for the low-dimensional problems described in the next subsections. First, our plug-in statistic uses $s_{Y}^{2}$ as an estimate for $\sigma^{2}$; why not use an estimate from the full least squares model of $Y$ on $X$, since this would be less conservative? While experiments (not shown) confirm that this works in low-dimensional regression problems, such an estimate becomes anti-conservative as the number of variables grows (particularly, irrelevant ones), and is obviously not applicable in high-dimensional problems. Therefore, we stick with the simple estimate $s_{Y}^{2}$, as this is always applicable and always conservative.

Second, to determine variable significance in a low-dimensional problem, one could of course fit a full regression model and inspect the resulting p-values and confidence intervals. These p-values and intervals could even be Bonferonni-adjusted to account for selection. Of course, this strategy would not be possible for a high-dimensional problem, but if the number of predictors is small enough, then it may work perfectly fine. So when should one use more complex tools for post-selection inference? This is an important question, deserving of study, but it is not the topic of this paper. The examples that follow are intended to portray the robustness of the selective pivotal inference method against nonnormal error distributions; they are not meant to represent the ideal statistical practice in any given scenario.

Figure 3(a) displays QQ plots of p-values for testing the significance of the variable entered into the active model, across 3 steps of LAR. (The QQ plots compare the p-values to a standard uniform distribution.) The p-values were computed using the TG statistic with $\sigma^{2}=1$, the plug-in TG statistic with $s_{Y}^{2}$ as its estimate for $\sigma^{2}$, and the bootstrap TG statistic with 50,000 bootstrap samples used to approximate the probabilities in the numerator and denominator of (20), and padding factor $\delta_{n}=10^{-4}n^{-1/4}$. (The scaling factor was ignored, i.e., set to $c=1$, for the plug-in and bootstrap statistics.) In steps 1 and 2, the p-values are restricted to repetitions in which a correct variable selection was made—i.e., variable 1 or 2 was entered into the active LAR model. In step 3, the p-values are from repetitions in which an incorrect variable selection was made—i.e., one of variables 3 through 10 was entered into the active model. Since the underlying signal was fairly strong and the predictors uncorrelated, such selections happened the majority of the time; specifically, the p-values displayed for steps 1, 2, and 3 comprise approximately 95%, 85%, and 87% of the 500 repetitions, respectively. The p-values in steps 1 and 2 show reasonable power, for all 3 statistics (TG, plug-in, and bootstrap types), and all 4 error distributions. Also, the p-values in step 3 are uniform, as desired, again for all statistics and all error distributions. Though the guarantees (for uniform null p-values) are only asymptotic for the Laplace, uniform, and skew normal error distributions, such asymptotic behavior appears to kick in quite early for these distributions, as the sample size here is only $n=50$. Further, the QQ plots reveal that the p-values for the nonnormal error distributions are not really any farther from uniform than they are in the normal case. This is somewhat remarkable, recalling that the p-values are, by construction, exactly uniform under normal errors.

Figure 3(b) inspects the TG statistic and plug-in and boostrap variants, when the pivot value $\mu$ is set to the true population value. That is, we set $\mu=v^{T}\theta$ in computing the statistics in (18), (19), and (20), in each data instance and each step of LAR. The figure collects the p-values across all 3 steps of LAR, for each of the 4 error distribution types. According to our theory, the distribution of the TG pivotal statistics here should be asymptotically uniform. This is clearly supported by the QQ plots. Interestingly, both plug-in and bootstrap pivotal statistics also appear uniform in the QQ plots, and yet, this is not a case handled by our asymptotic theory: recall, Theorem 11 fixes the pivot value $\mu$ to be 0 (as, otherwise, technical difficulties are encountered in its proof). This gives empirical evidence to the idea that a more refined analysis could extend Theorem 11 to the broader setting (of arbitrary pivot values) handled by Theorem 7. Moreover, it suggests that inverting the plug-in and bootstrap TG statistics should yield intervals with proper coverage, which is verified in the next subsection.

2 Confidence interval examples

We stay in same setting as the last subsection, so that $n=50$, $d=10$, and $\theta=X\beta_{0}$ for a coefficient vector $\beta_{0}$ with its first 2 components equal to $-4$ and $4$, and the rest equal to 0. We invert the TG, plug-in TG, and bootstrap TG statistics to obtain 90% confidence intervals at each LAR step. See Table 1 for a numerical summary. “Coverage” refers to the average fraction of intervals that contained their respective targets over the 500 repetitions, “power” is the average fraction of intervals that excluded zero, and “width” is the median interval width. These are all recorded in an unconditional sense, i.e., no screening of repetitions was performed based on the variables that were selected across the 3 steps of LAR (the conditional coverages however, were quite similar). From the table, we can see that all 3 methods lead to accurate coverage (around 90%) in all cases. We can further see that the intervals from the bootstrap TG statistic are shorter than those from the plug-in TG statistic in all cases, and considerably shorter than both the plug-in and original TG statistics in steps 2 and 3. The power from the bootstrap TG intervals is generally better than that from the plug-in TG intervals; also, it is on par with the power from the original TG statistic in step 1, but somewhat worse in step 2. Recall that the original TG statistic uses knowledge of the error variance ($\sigma^{2}=1$) but the bootstrap and plug-in variants do not.

It is a bit surprising that the bootstrap intervals can be shorter but still have worse power than the original TG intervals. This is easier to understand once the intervals are visualized, as done in Figure 4. The figure shows 100 sample intervals from the first LAR step, under normally distributed errors. Sample intervals from the other error models are shown in Appendix A.13. We see that the bootstrap TG intervals are indeed shorter, but compared to the original TG intervals, they are more symmetric around the target population values. The original TG intervals, being more asymmetric, are often shorter on the side (of the target value) facing 0, and this results in better power.

3 Heteroskedastic errors

4 High-dimensional examples

A negative result in high dimensions

We prove that the TG statistic fails to converge to a uniform distribution, under the null hypothesis, in a data model that has nonnormal errors and is high-dimensional, but otherwise represents a fairly standard setting: the “many means” setting. We write the observation model as

where we interpret $i=1,\ldots,m$ as replications, and $j=1,\ldots,d$ as dimensions. In total there are hence $n=md$ observations. Denote

We assume that the errors $\epsilon_{ij}$, $i=1,\ldots,m$, $j=1,\ldots,d$ in (25) are i.i.d. from the following mixture:

The mixing proportion $\pi$ and mean shift $B$ will both scale with $d$. Moreover, they will be chosen so that (for each $d$) the error variance is

As mentioned, we will consider model selection events of the form

We note that this is exactly the same selection event as that from the first step of FS, LAR, or lasso paths, when run on the regression version of this problem with orthogonal design $X$. It is not hard to check that the TG statistic for conditionally testing $\mu_{j}=0$, given that $\widehat{M}(Y)=(j,s)$, is

As per the spirit of our paper, we can also view this statistic unconditionally; for this it is helpful to define $W_{1}=|\overline{Y}_{1}|,\ldots,W_{d}=|\overline{Y}_{d}|$, and denote by $W_{(1)}\geq\ldots\geq W_{(d)}$ the order statistics. Then from (27), we can see that the unconditional TG statistic for testing the selected mean being 0 is

The framework underlying the TG statistic tells us that if the errors in (25) are i.i.d. $N(0,2)$, then for any fixed model $(j,s)$, the pivot $T(Y;j,s,0)$ is uniformly distributed conditional on $\widehat{M}(Y)=(j,s)$. Further, if $W_{(1)}$ and $W_{(2)}$ are the largest and second largest absolute values of centered normal random variables (each with variance $2/m$), then the unconditional pivot $\pazocal{T}(Y;0)$ is again uniform. But when $W_{(1)},W_{(2)}$ are large, and are defined by the order statistics of nonnormal random variates, the statistic $\pazocal{T}(Y;0)$—which in this case is defined by the extreme tail behavior of the normal distribution—could be nonuniform. The next theorem asserts that such nonuniformity does indeed happen asymptotically if we choose the mixture distribution in (26) appropriately.

Assume the observation model (25), where the errors are all drawn i.i.d. from (26). Let $d$ and $m$ scale in such a manner that $(\log{d})/m\to\infty$. Further, let

so that the error variance is fixed at $\sigma^{2}=2$. Then under the global null hypothesis, $\mu=0$, the unconditional TG statistic $\pazocal{T}(Y;0)$ in (28) does not converge in distribution to $U(0,1)$. In particular, on an event whose limiting probability is at least $1/e$, the statistic $\pazocal{T}(Y;0)$ converges to 0.

Further, the same results hold conditionally on any selected model. That is, for any fixed $(j,s)$, the conditional TG statistic $T(Y;j,s,0)\,|\,\widehat{M}(Y)=(j,s)$ does not converge in distribution to $U(0,1)$, and on an event with limiting probability (conditional on $\widehat{M}(Y)=(j,s)$) at least $1/e$, it converges to 0.

The assumed condition $(\log d)/m\to\infty$ requires the dimension $d$ to diverge to $\infty$, but not necessarily the number of replications $m$, though it clearly allows $m$ to diverge at a sufficiently slow rate. On the other hand, if $d$ were fixed and $m$ diverged to $\infty$, then the result of the theorem would no longer be true, and the limiting distribution of the TG p-value would revert to $U(0,1)$. (To be careful, here we would have cap the mixing probability $\pi$ at $1/2$ in order for the mixture to make sense, since the current definition of $\pi$ diverges with $d$ fixed and $m$ tending to $\infty$.) In fact, this is ensured by our low-dimensional result in Theorem 7: after reformulating the many means problem in appropriate regression notation, all of the conditions of Theorem 7 are met by our current setup when $d$ is fixed. This is supported by the simulation in Figure 7.

The precise scaling $(\log d)/m\to\infty$ is chosen since this implies $\pi=(1/d)^{1/m}\to 0$, i.e., the extreme mixture components $N(-B,1)$ and $N(B,1)$ each have probability tending to 0, an intuitively reasonable property for the error distribution. But we note that this scaling is not important for any other reason, and the proof would still remain correct if $d/m\to\infty$.

In Theorem 3 of Tian & Taylor 2017, the authors show that the TG statistic converges in distribution to a standard uniform random variable, in a high-dimensional problem setting, with some restrictions on the sequences of selection events that are allowed. One might ask what part of our high-dimensional setup here violates their conditions, because both results obviously cannot be true simultaneously. As far as we can tell, the issue lies in the role of $\delta_{n}$ in Assumption 1 of Tian & Taylor 2017. Namely, as we have defined the error distribution in (26), the value of $\delta_{n}$ needed to certify the third condition Assumption 1 of their work is too small for the main assumption in their Theorem 3 to hold. Hence Theorem 3 of Tian & Taylor 2017 does not apply to our current setup.

Discussion

We have studied the selective pivotal inference framework, with a focus on forward stepwise regression (FS), least angle regression (LAR), and the lasso, in regression problems with nonnormal errors. We have shown that the truncated Gaussian (TG) pivot is asymptotically robust in low-dimensional settings to departures from normality, in that it converges to a $U(0,1)$ distribution (its pivotal distribution under normality), and does so uniformly over a broad class of nonnormal error distributions. When the error variance $\sigma^{2}$ is unknown, we have proposed plug-in and bootstrap versions of the TG statistic, both of which yield provably conservative asymptotic p-values.

Our numerical experiments revealed that the statistics under theoretical investigation generally display excellent finite-sample performance, for highly nonnormal error distributions. These experiments also revealed findings not predicted by our theory: (i) the bootstrap TG statistic often produces shorter confidence intervals than those based on the plug-in TG statistic, and even the TG statistic that relies on the error variance $\sigma^{2}$; and (ii) all three TG statistics show strong empirical properties well-outside of the classic homoskedastic, fixed $d$ regression setting that we presumed theoretically.

However, as we have clearly demonstrated, one should not hope for a convergence result in high dimensions that is as general as the result obtained in low dimensions. In a relatively simple many means problem, we showed the nonconvergence of the TG statistic to $U(0,1)$ as $d\to\infty$, whereas in the same problem but with $d$ fixed, the TG statistic converges to its usual $U(0,1)$ limit.

There is still much left to do in terms of understanding the behavior of selective pivotal inference tools that are constructed to have exact finite-sample guarantees under normality, like the TG statistic of Tibshirani et al. 2016, when applied in high-dimensional regression settings with nonnormal data. When the pivot, the central cog of this framework, is constructed under the assumption of normality, this creates robustness issues that are especially worrisome in high dimensions. Appendix A.16 provides a high-level discussion of some of these issues; a more detailed study will be the subject of future research.

We thank Jelena Markovic and Jonathan Taylor for many helpful discussions, and for their overall generosity. An initial version of our work contained only unconditional (i.e., marginal) results in the main theorems (Theorems 7, 11, and 12); Jelena Markovic pointed out that Theorem 7 should also hold conditionally, and the current version of this work has been revised accordingly.

Appendix A Appendix

We describe a modification of the conic conditioning set in Tibshirani et al. 2016 for FS. Our version is different in that we additionally condition on the sign of every active coefficient at every step, rather than just the coefficient of the variable to enter the model at each step. The modifications needed for the LAR and lasso conditioning sets, made on top of the sets for LAR and lasso given in Tibshirani et al. 2016, will follow similarly to that described for FS, and hence we omit the details.

where $\widetilde{X}_{j}$ is the residual from regressing $X_{j}$ onto $X_{A_{k-1}}$, and $r$ is the residual from regression $y$ onto $X_{A_{k-1}}$. By expressing $\widetilde{X}_{j}=P_{A_{k-1}}^{\perp}X_{j}$ and $r=P_{A_{k-1}}^{\perp}X_{j}$, where $P_{A_{k-1}}^{\perp}$ projects onto the orthocomplement of the column space of $X_{A_{k-1}}$, we can rewrite the above constraints as

a set of $2(d-k)$ linear inequalities in $y$. Meanwhile, the subevent $\widehat{s}_{k}(y)=[s_{k,1},\ldots,s_{k,k}]$ can be characterized by $k$ inequalities expressed in block form,

A.2 Proof of Lemma 3

We prove the result for FS; the results for the LAR and lasso paths follows similarly, by inpsecting the form of the linear inequalities that determine their selection events.

Consider the first FS step as described in Appendix A.1. Multiplying through by $\sqrt{n}$, we see that an equivalent set of inequalities that characterize the selection event $\widehat{j}_{1}(y)=j_{1}$, $\widehat{s}_{1}(y)=s_{1}$ is

This is clearly of the desired form $P_{1}(\frac{1}{n}X^{T}X)\,\frac{1}{\sqrt{n}}X^{T}y\geq 0$, for a matrix $P_{1}(\frac{1}{n}X^{T}X)$ dependent only on $\frac{1}{n}X^{T}X$. At the $k$th step of FS, there are two sets of inequalities to be examined: one that describes the variable to enter $\widehat{j}_{k}(y)=j_{k}$, and the second that describes the active signs $\widehat{s}_{k}(y)=[s_{k,1},\ldots,s_{k,k}]$. The first set, multiplying through by $\sqrt{n}$, is

while the second set, again multiplying through by $\sqrt{n}$, is

These inequalities are clearly all summarized by $P_{k}(\frac{1}{n}X^{T}X)\,\frac{1}{\sqrt{n}}X^{T}y\geq 0$, where $P_{k}(\frac{1}{n}X^{T}X)$ is a matrix that depends only on $\frac{1}{n}X^{T}X$. This completes the proof.

A.3 Proof of Lemma 4

Under the conditions of the lemma, the TG pivot for fixed $M$ in (8) depends only on $X,y$ through the master statistic, because, as explained above the lemma, the only dependence in the pivot on $X,y$ is through the quantities $v^{T}y/\|v\|_{2}$, $(Q_{M}(X)\,v)/\|v\|_{2}$, $Q_{M}(X)\,y$, and each of these is in turn a function of the master statistic $\Omega_{n}$. Moreover, we may reexpress the TG statistic in (10) as

for some functions $f_{1},f_{2},f_{3}$, or more succinctly, as $T(X,y;M,v,\mu)=\psi_{M}(\frac{1}{n}X^{T}X,\frac{1}{\sqrt{n}}X^{T}y)$, where

Note that the quantities $v^{T}y/\|v\|_{2}$, $(Q_{M}(X)\,v)/\|v\|_{2}$, $Q_{M}(X)\,y$ depend smoothly on the master statistic $\Omega_{n}=(\frac{1}{n}X^{T}X,\frac{1}{\sqrt{n}}X^{T}y)$ at any point such that $\frac{1}{n}X^{T}X$ is nonsingular. This implies $f_{1},f_{2},f_{3}$ are smooth functions of $(S,z)$ at any point such that $S$ is nonsingular. Lastly, for all $S,z$ such that $P_{M}(S)\,z>0$, we have $f_{1}(S,z)>f_{3}(S,z)$, and thus the denominator of $\phi_{M}(S,z)$ is positive. This proves the desired continuity result on $\psi_{M}$.

A.4 Proof of Lemma 5

where we use $S_{A,A}$ to denote the submatrix of $S$ with rows in $A$ and columns in $A$, and $z_{A}$ to denote the subvector of $z$ with entries in $A$.

A.5 Proof of Lemma 6

Define $Z_{0,n}=\sum_{i=1}^{n}\xi_{i}$, where $\xi_{i}=\frac{1}{\sqrt{n}}x_{i}\epsilon_{i}$, $x_{i}$ is the $i$th row of $X$, and $\epsilon_{i}=Y_{i}-\theta_{i}$, for $i=1,\ldots,n$. Note that $(\xi_{1},\ldots,\xi_{n})\sim F_{n}(0)$, with independent, mean zero components. We compute

which converges to $\sigma^{2}\Sigma$ as $n\to\infty$, by assumption. Further, for any $\delta>0$, consider

Now consider $Z_{n}=\frac{1}{\sqrt{n}}X^{T}Y=Z_{0,n}+\frac{1}{\sqrt{n}}X^{T}\theta$. Writing $\Phi$ and $\phi$ for the standard normal CDF and density,

where the second line is due to the triangle inequality, and the third line is due to the simple bound $|\Phi(x-t)-\Phi(x-s)|=|\int_{x-s}^{x-t}\phi(u)\,du|\leq|t-s|\phi(0)$, for any $x,s,t$. Note that $a\to 0$ by the argument at the start of this proof, and $b\to 0$ by assumption in (17). This shows that $Z_{n}$ converges in distribution to $Z\sim N(\eta,\sigma^{2}\Sigma)$, uniformly over $\pazocal{P}_{n}(\theta)$, and over $\theta\in\Theta$.

Lastly, we establish the conditional result. By repeating the same arguments as above, the uniform Lindeberg-Feller central limit theorem and condition (17) imply that $(Z_{n},A_{n}Z_{n})$ converges to $(Z,AZ)$, uniformly over $\pazocal{P}_{n}(\theta)$, and over $\theta\in\Theta$. Thus, along sequence $F_{n}(\theta)\in\pazocal{P}_{n}(\theta)$, $n=1,2,3,\ldots$ with $\theta\in\Theta$, observe

at a rate that does not depend on the sequence in question. This is true because the numerator and denominator each converge to their normal probability counterparts, and the denominator remains bounded away from zero since $\{z:Az\geq 0\}$ has nonempty interior, and the set of limits of $\frac{1}{\sqrt{n}}X^{T}\theta$ was assumed compact, in (17). Since $x$ was arbitrary, and the distribution of $Z\,|\,AZ\geq 0$ is continuous, we have (e.g., Lemma 2.11 in van der Vaart 1998)

And as the sequence $F_{n}(\theta)\in\pazocal{P}_{n}(\theta)$, $n=1,2,3,\ldots$ with $\theta\in\Theta$ was arbitrary, we have shown the desired uniform convergence.

A.6 Proof of Theorem 7

We begin with the proof of part (a). Let $Z_{n}=\frac{1}{\sqrt{n}}X^{T}Y$ and $Z\sim N(\eta,\sigma^{2}\Sigma)$. Also, let $A_{n}=P_{M}(\frac{1}{n}X^{T}X)$ and $A=P_{M}(\Sigma)$. Recall that $A_{n}Z_{n}\geq 0\iff\widehat{M}(X,Y)=M$, by Lemma 3. Also, $Z_{n}\,|\,A_{n}Z_{n}\geq 0$ converges weakly to $Z\,|\,AZ\geq 0$, uniformly over $\pazocal{P}_{n}(\theta)$ and over $\theta\in\Theta$, by Lemma 6. As $\frac{1}{n}X^{T}X\to\Sigma$ deterministically, we also have that $\Omega_{n}=(\frac{1}{n}X^{T}X,Z_{n})$ converges uniformly in distribution to $\Omega=(\Sigma,Z)$.

The choice of $v$ as specified in the theorem is now important for two reasons. First, by Lemma 4, we can express

The proof of part (b) follows from the expansion

As the number possible models $|\pazocal{M}|$ is finite, we can simply apply the asymptotic pivotal result from part (a) to each $M\in\pazocal{M}$ to establish the asymptotic pivotal property of $\pazocal{T}(X,Y;V,U)$. The confidence interval result is again just a rearrangement of this pivotal property.

A.7 Proof of Lemma 8

By assumption, the vector $v$ can be written as

The denominator converges to $|e_{j}^{T}(\Sigma_{A,A})^{-1}e_{j}|^{3/2}$ by (14). The numerator satisfies

where $a$ is bounded by (21) and $b$ converges to $\|(\Sigma_{A,A})^{-1}e_{j}\|_{2}^{3}$ by (14). This completes the proof.

A.8 Proof of Lemma 9

We start by proving the result about the event $\{cs_{Y}\geq\sigma\}$. First let us study its asymptotic probability marginally. Consider

where $\overline{\epsilon}=\sum_{i=1}^{n}\epsilon_{i}/n$. Hence

for a constant $C_{t}>0$ only depending on $t$. Hence, observe that