More powerful post-selection inference, with application to the Lasso
Keli Liu Jelena Markovic Robert Tibshirani
stat.ME



1	

SDN based Control and Management of WLANs in the 3GPP 5G Network
    
 Pranav Jha, Abhay Karandikar
Department of Electrical Engineering, Indian Institute of Technology Bombay 

Email: {pranavjha, karandi}@ee.iitb.ac.in

    December 30, 2023
===========================================================================================================================================================



 Investigators often use the data to generate interesting hypotheses and then perform inference for the generated hypotheses. P-values and confidence intervals must account for this explorative data analysis. A fruitful method for doing so is to condition any inferences on the components of the data used to generate the hypotheses, thus preventing information in those components from being used again. Some currently popular methods “over-condition,” leading to wide intervals. We show how to perform the minimal conditioning in a computationally tractable way. In high dimensions, even this minimal conditioning can lead to  intervals  that are too wide to be useful,  suggesting that up to now the cost of hypothesis generation has been underestimated. We show how to generate hypotheses in a strategic manner that sharply reduces the cost of data exploration and results in useful confidence intervals. Our discussion focuses on the problem of post-selection inference after fitting a lasso regression model, but  we also outline its extension to a much more general setting.





§ INTRODUCTION



In many applications we want to study how the mean of a response variable
is affected by various predictors. In idealized form, suppose we have the response vector y=(y_1,…,y_n)∈ℝ^n following

    y_i=μ_i+ε_i    ε_i∼𝒩(0,σ^2)     i=1,…,n.

The goal is to study the relation between the mean function μ
and a set of p explanatory variables x_1,…,x_p. The
standard approach to this problem approximates μ=(μ_1,…,μ_n)∈ℝ^n using a linear
function of the explanatory variables: 

    μ_lin(x)=x^⊤β    β=(𝐗^⊤𝐗)^-1𝐗^⊤μ,

where x∈ℝ^p and 𝐗∈ℝ^n× p is a matrix of features. An unbiased estimate of β is given by the least squares estimator
of y on the matrix 𝐗; in addition, standard theory allows us to compute p-values and confidence intervals for β.

The standard approach was developed for the setting where the investigator, drawing on prior knowledge, has already focused in on a small set of explanatory variables. It is wholly inappropriate for tasks where the majority of variables are noise and exploratory analyses must be done to arrive at the “variables of interest.” In the latter case,
a popular approach for exploration is the lasso of <cit.> which solves the ℓ_1-penalized regression 

    β_0,β_1,…,β_p minimize 1/2∑_i=1^n(y_i-β_0-∑_j=1^px_ijβ_j)^2+λ∑_j=1^p|β_j|.

For fixed choice of the tuning parameter λ, the lasso selects a set of interesting variables (the “active set”) while setting β̂_j=0 for uninteresting variables. How can we compute valid p-values and confidence intervals while accounting for exploration using the lasso? That is the problem that we address here. One might think of taking the active set, and computing the standard/“naïve” p-values and confidence intervals based on least squares theory. But the standard theory assumes pre-specified hypotheses; without accounting for the fact that our hypotheses are generated from the data, it is quite common for a naïve 90% confidence interval to cover less than 50% of the time.

Performing inference after selection is notoriously difficult, e.g., <cit.> show the impossibility of estimating the conditional distribution of a test statistic used for inference, implying failure of the standard bootstrap for this purpose. On the other hand, <cit.> and <cit.> propose a conditional approach that leads to a truncated normal reference
distribution, that gives provably exact type I error and coverage
in finite samples. A drawback of this quite remarkable construction
is the fact that the confidence intervals are often very wide. The intervals “over-condition,” that is they overstate the cost of data exploration using the lasso (by quite a bit it turns out!). 

Consider data from <cit.> who study the relation between prostate specific antigen (PSA) and various clinical measures: log cancer volume (), log prostate weight (), age, log of benign prostatic hyperplasia amount (), seminal vesicle invasion (), log of capsular penetration (), the Gleason score (), and percent of Gleason scores 4 or 5 (). The data were collected from n=97 men
about to receive a radial prostatectomy. Running the lasso on the
prostate data (with λ chosen by 10 fold cross-validation)
yields seven active variables (the log of capsular penetration was
excluded from the model). Figure <ref> plots 90% confidence intervals for the coefficients of the seven selected variables. The coefficients can be defined either with respect to a linear model of all the variables (“full” coefficients) or a linear model of only the selected variables (“partial” coefficients). We summarize  the findings from the figure as follows:



  * The “naive” intervals are based on standard normal theory and ignore selection. They are too short and hence will undercover their targets.

  * The various “TZ” methods account for selection:
operationally, this implies that the least squares estimates
for the selected variables have a truncated normal distribution (hence TZ).  TZ_M and TZ_Ms  are the approaches introduced by <cit.> and <cit.>. They have exact finite-sample coverage but often produce very long intervals.  Of the three variables (,  and ) deemed highly significant by a least squares analysis, two now have 90% confidence intervals containing 0, and the third has a very wide and right skewed interval.

  * The  truncated Z-statistic methods TZ_V and TZ_stab-t — introduced in this paper — have exact finite-sample coverage and give much shorter intervals compared to TZ_M and TZ_Ms.





Where the various TZ methods differ is in the amount of conditioning.  This distinction is crucial; in a nutshell:


More Conditioning  less information for inference   wider confidence intervals

See <cit.> for justification of this maxim. Heuristically, not all of the variation in
the Z-statistic can be used for inference because some of that variation has already been used for exploratory data analysis (in this case to screen out unimportant variables). As a result,  the reference distribution for the Z-statistic changes from a normal to a truncated normal. TZ_Ms sets aside information contained in the lasso active set (M̂) as well as information in the signs of the selected variables (ŝ);
TZ_M sets aside only the former leading to less truncation (the
subscripts indicate what information has been excluded for inference
purposes). TZ_V and TZ_stab-t set aside even less
information. But what is the minimum level of truncation needed for
the amount of explorative data analysis we have done? And how should
we change our data analysis strategy if this minimal truncation
still turns out to be too severe? This paper answers these two questions as follows: 

	
  * To arrive at the minimal truncation, we first make clear in Section <ref> that explorative data analysis usually incurs two separate “search costs”: the choice of variables on which
	to do inference (e.g., the active set from the lasso), and the choice of targets for inference (Which model should we define our regression coefficients with respect to?). This has not been recognized previously, as the two sources are usually lumped together conceptually. But it has significant practical implications. 
	
  * In Section <ref>, we consider the case where the targets for inference for the selected variables are the full population regression coefficients. There is no adaptivity in the choice of target. We show how to construct intervals that perform the minimal truncation necessary to adjust for the lasso selection event; this is
	the TZ_V method. 
	
  * In Section <ref>, we consider the case where the targets for inference for the selected variables are the population partial regression coefficients. There is significant adaptivity in the choice of targets—the definition of “coefficient” changes depending on which model is chosen. TZ_M is shown to
	give the minimal necessary truncation, but its intervals are impractically wide. In such instances, blindly using the lasso to form hypotheses is a bad idea! We provide a framework for judiciously generating hypotheses using the lasso. TZ_stab-t gives the truncation appropriate for this judicious approach and produces intervals of reasonable length.

Section <ref> contains simulation results and Section <ref> illustrates the benefits of our
proposals on a larger real-world dataset (n=1057 and p=210).
While we use the example of the lasso to simplify exposition, the
tools we develop can be adapted to treat other polyhedral selection
rules, e.g., forward stepwise regression. We give some details of
the general idea in Section <ref>.




 §.§ Review and Critique of Post-Selection Inference via the Truncated Gaussian
  

We briefly review the work of <cit.> and <cit.>. Let  be our n× p predictor matrix and y our n-dimensional outcome vector. To focus in on interesting variables, suppose we have run the lasso with a fixed λ and let M̂ be the active set the lasso returns. We can now use M̂ to generate interesting questions. For example, we might be interested in β^(M̂)=(𝐗_M̂^⊤𝐗_M̂)^-1𝐗_M̂^⊤μ, the coefficients in the projection of μ onto the selected variables. Hypothesis tests of the form H_0:β_j^(M̂)=0 have the general linear form H_0:η^⊤μ=0 for some vector η∈ℝ^n. Standard theory
tells us what to do when η is pre-specified: the Z-statistic

    Z=η^⊤y-η^⊤μ/σ‖η‖_2∼𝒩(0,1)

is a pivot from which we can derive p-values and confidence intervals. However, this theory does not hold because our η depends on the data through M̂, i.e., η=η(M̂) and upon substituting this into the above expression, the pivotal relationship no longer holds.

When performing inference for such data generated hypotheses, we should
remove whatever information has already been used up during our initial
exploration in order to avoid “using the data twice.” This is
achieved by conditioning on M̂ (note that conditional on M̂=M,
η(M̂) is back to being a fixed vector). Post-selection
p-values and intervals can then be built in the usual manner if we can identify quantities that are pivotal conditional on M̂. <cit.> showed how to build such a conditional pivot. We summarize their construction: 

	
  * The event {M̂=M} is difficult to characterize
	but <cit.> and <cit.> reveal a simple characterization
	of a more refined event 
	
    {M̂=M,ŝ=s}={ y≤ b},

	where ŝ is the vector of signs of the lasso coefficients
	of the active variables. The RHS above is a polyhedral region in outcome space (the exact form for  and b can be derived from the KKT conditions at the solution). <cit.> and <cit.> construct a quantity that is pivotal conditional on M̂ and ŝ (which implies that it is pivotal conditional on M̂). 
	
	
  * The standard pivot (<ref>) simply centers and scales η^⊤y by its unconditional mean and standard deviation.
	The conditional distribution of η^⊤y given 𝐀y≤ b is difficult to work with because it depends on nuisance parameters. The usual way to get rid of nuisance parameters is by conditioning; in this case, we condition on P_η^⊥y, the piece of y orthogonal to η. <cit.> show that 
	
    η^⊤y  | 𝐀y≤ b, P_η^⊤y∼𝒯𝒩(η^⊤μ,σ^2‖η‖_2^2,[𝒱^-,𝒱^+])

	where 𝒯𝒩(μ,σ^2,[c,d]) denotes the
	𝒩(μ,σ^2) distribution truncated to [c,d]. Above, 𝒱^- and 𝒱^+ are computable functions of 𝐀, b and P_η^⊤y (see Appendix <ref>). This implies that the usual Z-statistic (<ref>) has a truncated normal distribution as its conditional!
	
  * The CDF of a truncated Gaussian distributiontruncated normal distribution with support confined to [c,d] is 
	
    F_μ,σ^2^c,d(x)=Φ((x-μ)/σ)-Φ((c-μ)/σ)/Φ((d-μ)/σ)-Φ((c-μ)/σ),

	with Φ the CDFcumulative distribution function of the standard Gaussian. The CDF of a random variable, evaluated at the value of that random variable, has a uniform distribution on (0,1). Hence 
	
    F_η^⊤μ,σ^2‖η‖_2^2^𝒱^-,𝒱^+(η^⊤y) | {M̂=M,ŝ=s}∼𝒰(0,1).


We obtain confidence intervals by inverting the pivotal quantity (<ref>): we denote these as the TZ_Ms intervals, for truncated Z-statistic intervals conditional on the selected model variables M and their signs s. For a 1-α interval, we find the largest and the smallest
η^⊤μ such that the value of pivotal quantity remains
in the interval [α/2,1-α/2].


The technique of <cit.> and <cit.> allows for
exact inferences after using the lasso (or any other procedure whose
selection event is polyhedral) to generate hypotheses. However, it
is not without shortcomings: 

	
  * The TZ_Ms intervals condition on (M̂,ŝ)
	rather than just M̂. This would be appropriate if we also
	used information in the signs when forming our hypotheses. For example,
	we might test against a one sided alternative, H_1:β_j^(M̂)>0,
	if the sign of the jth variable is positive. However, in most applications,
	the signs are not used in hypothesis generation, hence we have thrown
	away unused information by conditioning on ŝ. As a result,
	the intervals tend to be quite long (Figure <ref>).
	
  * <cit.> makes mention of the previous point and shows how
	to condition only on M̂ by an enumeration of all possible
	sign vectors. Their enumeration method is intractable when |M̂|
	is large as there are 2^|M̂| possible sign vectors.
	In Section <ref>, we provide a workable implementation. 
	
  * In some applications, M̂ is used to generate hypotheses in a very restrictive manner. The restrictions allow us to condition on less than M̂, thus leading to more precise intervals; see Section <ref>. 

Appendix <ref> presents a simple example that shows the effects of different conditioning strategies.



The  above criticisms have focused on the “over-conditioning” of
the TZ_Ms intervals. There is however an entirely different (and
practically pressing) concern: What if the intervals are uselessly
wide even when we adopt the “minimal” conditioning? It turns out
that when we are interested in the partial regression coefficients,
conditioning on M̂ is the “minimal” conditioning; yet
as Figure <ref> shows, the resulting TZ_M intervals
are still unacceptably wide. Such imprecision is in fact the normal
state of affairs when p>n (see Section <ref>).
What such an occurrence should tell us is that we have used up too much
information in our initial data exploration—and did not
leave enough information for inference. 

Thus we need to strike a more
balanced approach to how we use the data to generate hypotheses. <cit.>
offer one such approach: perform our exploratory analysis not on the
actual data but a noise-perturbed version of it (similarly, one could
split the data). That is, we would obtain our active set not by running
lasso on the original data but by running lasso on the original data
with noise added on. The resulting intervals are substantially shorter
than the TZ_Ms intervals but do require greater computation
in the form of MCMC sampling. More concerning, however, is that the
variables we select will vary across realizations of noise and will
differ from the variables we select had we used the original data.
This paper takes a different route by giving a deterministic method
for controlling the cost of data exploration: we will continue to
generate hypotheses using M̂, the lasso active set from the
original data, but we will restrict how M̂ can be used (Section
<ref>). The key to our development is distinguishing
between two often confounded ways in which the data can be used to
shape the direction of our investigation.



 §.§ Forming a Data Driven Query: Two Costs
 

In many applications, investigators come
to the data without pre-specified questions. Instead, they want to use the data to generate interesting hypotheses. Often, the hypothesis generation serves to answer two questions. Out of a possibly larger set of explanatory variables, which ones should we focus on? What aspect of the interesting variables should we study? The two processes that answer these questions are: 

	
  * Variable Selection: The data is used to decide which variables are worthy of attention, e.g., running the lasso and focusing on the active set. 
	
  * Target Formation: Having settled on a subset M⊂{ 1,…,p}
	of variables for careful study, we need to decide how to summarize the relation between the selected variables and μ. That is, what should be the target of our estimation? We could try and estimate
	the regression coefficients from projecting μ onto 𝐗_M:
	
    β^(M)=(𝐗_M^⊤𝐗_M)^-1𝐗_M^⊤μ,

	summarizing the effect of variable j through β_j^(M). In this approach, the way of capturing the effect of variable j depends on what other variables are in M. We call this the partial target, since they are the (population) partial regression coefficients.
	
	In contrast, prior to looking at the data, we could decide to always
	summarize the effect of j via its coefficient β_j^F
	in a regression of μ on all the variables 
	
    β^F=(𝐗^⊤𝐗)^-1𝐗^⊤μ.

	In this case, we have a fixed (non data-driven) target, β^F.
	We call this the full target, since they are the (population)
	regression coefficients in the full model. While we may use variable
	selection to decide which components of β^F to estimate,
	we do not use the data to inform how we should summarize the relation
	between selected variables and μ. 

Both variable selection and target formation carry a cost that needs
to be acknowledged when doing inference. When we condition on the
lasso active set M̂ for post-selection inference, our inferences
are valid regardless of how M̂ is used in variable
selection and target formation. However, if we know the precise way
in which M̂ is used for these two processes, we can often
do (much) better than simply conditioning on M̂. For example,
in Section <ref>, we consider a procedure where
only the “best variables” in M̂ are used in target formation
rather than all of M̂—hence we only need to condition
on the identity of the “best variables.” For these reasons, distinguishing
between the costs of variable selection and target formation is key
to improving over the standard truncated-Z inferences.



§ INFERENCE FOR THE FULL REGRESSION TARGETS
 


In this section, we are concerned with the (easier) case where
the data has been used to choose a set of interesting variables but
the data has not been used to decide how to summarize the
relation between the response and the selected variables. In particular,
we assume that the investigator has decided a priori that the relation
between a variable x_j and the response variable should be captured
by its coefficient in a regression of μ on all the variables:


    β^F=(𝐗^⊤𝐗)^-1𝐗^⊤μ.

The data is only used to choose the components of β^F we
are interested in. This is the full target discussed above.

We focus on the case where the lasso is used to choose the set M̂
of variables. The relation between μ and each variable j∈M̂
is summarized by β_j^F (again, we emphasize that this choice
of summary was agreed upon before exploring the data). Now the goal
is to build a confidence interval C_j for β_j^F. We
note that an interval is built for β_j^F whenever j∈M̂⊂{ 1,…,p}.
Hence, it seems natural to require that our interval have good coverage
properties conditional on the interval existing in the first place:

    ℙ(β_j^F∉ C_j | j∈M̂)≤α

for a given level α. In particular, we do not condition on the entire active set:
when hypothesis generation involves only variable selection, it suffices
to condition on variable j being in the active set (when performing
inference for j). In contrast, <cit.> conditions on both
the selected model M̂ as well as the signs ŝ. By
conditioning on less, more variation remains in the data for determining
the precise value of β_j^F. Since the goal of conditioning
is to restrict the sample space to those data where confidence intervals
for β_j^F will actually be constructed, the event { j∈M̂}
is precisely the minimal conditioning event.

To construct intervals conditional on { j∈M̂},
we provide a particularly simple characterization of this set. For
motivation, consider the case where 𝐗 has orthonormal columns in
which case β_j^F=η_j^⊤μ where η_j=x_j.
Then the event { j∈M̂} reduces to a union
of two disjoint intervals

    j∈M̂|η_j^⊤y|>λ.

That is, as long as the inner product between x_j and y exceeds
a cutoff, lasso will select variable j. Since η_j^⊤y∼𝒩(β_j^F,σ^2),
the conditional distribution of η_j^⊤y |  j∈M̂ is
simply this normal distribution truncated to [-∞,-λ]∪[λ,∞].
A pivot can be obtained from this truncated normal distribution
and then inverted to obtain confidence intervals. 

In the general case, we no longer have η_j=x_j. Instead,

    β_j^F=η_j^⊤μ    η_j=𝐗(𝐗^⊤𝐗)^-1e_j,

where e_j∈ℝ^p is a vector of all zeros with one at its jth coordinate.
It turns out that this general case can be reduced to the special
case of orthonormal predictors by conditioning. In particular, let
us decompose variation in y into two components: variation in the
direction of η_j and variation orthogonal to η_j.

    y=z_j· c_j+ν_j     z_j=η_j^⊤y     c_j=η_j/‖η_j‖_2^2.

We will condition on the component orthogonal to η_j, i.e.,
ν_j, and show that the law of η_j^⊤y | j∈M̂,ν_j
is the 𝒩(β_j^F,σ^2‖η_j‖_2^2)
distribution truncated to [-∞,a_j]∪[b_j,∞]
for some constants a_j and b_j. That is to say, the general
case simply replaces λ with a_j and b_j but much
of our intuition in the orthogonal case still transfers over. Such a simple result is hinted at by <cit.> who derived truncation regions of this form in the context of sequential testing. We shall refer to intervals based on this truncated distribution for the Z-statistic as _V intervals (the V here stands for variable inclusion); intervals which condition on M̂ or both M̂ and ŝ are similarly termed the _M and _Ms intervals respectively (see Section <ref>).


	
	The conditional distribution η_j^⊤y | j∈M̂,ν_j
	is 𝒩(β_j^F,σ^2‖η_j‖_2^2)
	truncated to [-∞,a_j]∪[b_j,∞]
	where
	
    a_j=‖η_j‖_2^2·(x_j^⊤r_j-λ)     b_j=‖η_j‖_2^2·(x_j^⊤r_j+λ)

	r_j=𝐗_-jβ̂_-j-ν_j and ν_j=(𝐈_n-1/‖η_j‖_2^2η_jη_j^⊤)y. Here, x_j is the jth column of 𝐗, 𝐗_-j is the submatrix of 𝐗 after removing the jth column of 𝐗 and β̂_-j is the lasso fit of ν_j to 𝐗_-j with penalty λ.





	The argument is simply an examination of the KKT conditions for the
	lasso. Using our decomposition for y and noting 
	
    𝐗^⊤c_j=e_j/‖η_j‖_2^2

	we have 
	
    𝐗^⊤(𝐗β̂-y) + λŝ=𝐗^⊤(𝐗β̂-ν_j)-e_j·z_j/‖η_j‖_2^2+λŝ.

	 The lasso solution satisfies 
	
    𝐗^⊤(𝐗β̂-ν_j)-e_j·z_j/‖η_j‖_2^2+λŝ=0 
    ŝ_j=sign(β̂_j)    β̂_j≠ 0 
    ŝ_j∈[-1,1]    β̂_j=0.

	Since ν_j is fixed, our goal is to find the range of z_j=η_j^⊤y
	values for which β̂_j≠0. When β̂_j=0, we have 
	
    |x_j^⊤(𝐗_-jβ̂_-j-ν_j)-z_j/‖η_j‖_2^2|<λ.

	This shows that β̂_j=0 if and only if z_j∈[a_j,b_j], where 
	
    a_j=‖η_j‖_2^2·(x_j^⊤r_j-λ)         b_j=‖η_j‖_2^2·(x_j^⊤r_j+λ) 
    
    	r_j=   𝐗_-jβ̂_-j-ν_j.

	Note that in the orthonormal case, x_j^⊤r_j=0 and ‖η_j‖_2^2=1
	so that a_j=-λ and b_j=λ.


When p > n, the matrix 𝐗^⊤𝐗 is no longer invertible. However, the full target can be approximated by β^F = 𝐀𝐗^⊤μ, where 𝐀 is some approximation to the inverse of 𝐗^⊤𝐗; see, e.g., <cit.>. However, under this construction, we no longer have the condition that η_j is orthogonal to x_k for all k≠ j, i.e., as we move along the direction η_j we affect the lasso solution not only for variable j but also for the other variables. As a result, the proof of Proposition <ref> no longer holds and we do not have easily computed formulas characterizing {j ∈M̂}. Nonetheless, it is still feasible to compute the truncation set via pathwise algorithms (Section <ref>).




 §.§ Simulations for the full target case


We conduct some simple simulations to demonstrate the benefit of conditioning only on j∈M̂ (the minimal conditioning) rather than the entire active set M̂.

	
  * Since the full regression coefficients only exist when n>p, we consider two cases: (n,p)=(100,50) and (n,p)=(100,90), the latter of which is close to being non-identifiable.
	
  * The variables for each individual are independent standard multivariate normal: x_i∼𝒩(0,𝐈_p). (Note x_i, i=1,…,n, denote the rows of 𝐗 and x_j, j=1,…,p denote the columns of 𝐗.)
	
  * The outcome y was generated as y_i=x_i^⊤β^F+ε_i, where ε_i∼𝒩(0,1), for all i=1,…,n.
	
  * The first k=5 components of β_j^F were set to either 0, δ_low, or δ_high (see following bullet). The remaining components were set to 0.
	
  * We consider three signal settings: no signal, low signal, and high signal. To arrive at a meaningful definition of low and high signal,
	we considered the distribution of max_jn^-1|x_j^⊤y| under the global null (where y∼𝒩(0,𝐈_n)). Asymptotically,
	this noise distribution concentrates around √(2log p/n);
	for our purposes we simulated from this distribution 1000 times and
	set δ_low as the median of this distribution and
	δ_high as the 99th quantile of this distribution
	plus an additional offset of 0.25. When the signal β_j^F=δ_low,
	it is difficult to distinguish variable j from noise; when β_j^F=δ_high,
	the signal has exceeded the noise threshold and is easy to detect.
	
  * For our penalty parameter λ, we considered two choices. For
	the first choice we try to approximate the λ that would have
	been picked by cross validation. To do this, we ran lasso at each
	setting 100 times, computed the λ that was chosen by 10-fold
	cross validation each time, and set λ to be the median of
	the resulting distribution. For our second choice, we chose the universal threshold value: λ=√(2log p/n). In all instances the “cv” λ is quite a bit smaller than the “universal threshold” λ; cross validation prioritizes predictive performance and is not model selection consistent (it lets in too many noise variables) while the universal threshold is model selection consistent.
	
  * We simulated 1000 independent datasets. For each dataset, we constructed a 90% confidence interval of β_j^F for each variable j in the lasso active set. 

Figures <ref>, <ref>, and <ref> display boxplots of interval lengths produced by each method, as well as giving coverage information. When inverting the truncated normal pivot, numerical inaccuracies occasionally result in “infinite” intervals (this occurs if the observed Z-statistic is too close to an endpoint of the truncation region). The median length and empirical coverage can still be computed using the usual definitions; for display purposes, however, our boxplots truncate the length of intervals with infinite length to the maximum finite length observed in a simulation. We summarize the key takeaways:

	
  * Many practitioners incorrectly assume that Bonferroni corrected intervals can protect against the effects of selection. When there is no signal β^F=0, Figure <ref> shows that the Bonferroni intervals can severely under-cover. 
	
  * Conditioning on j∈M̂ rather than M̂ produces intervals that are significantly shorter. The benefits are greater when λ is smaller (chosen by cross validation). This is because M̂ carries less information about y when λ is larger (say, near the universal threshold). In other words, increasing λ is one way to control our exploration cost, a fact that we shall exploit in Section <ref>. For sufficiently large λ,
	we will obtain M̂=∅ with high probability in which case there is no exploration cost because there is no exploration.
	
  * _V produces infinite length intervals much less frequently than _M or _Ms (almost never). Since truncation is less severe for _V, the observed Z-statistic is less likely to fall near the left-most or right-most boundary of the truncation region.
	
  * The simulation results are consistent with the good performance of the _V intervals observed for the prostate data (Figure <ref>) and HIV data (Figure <ref>).











§ INFERENCE FOR PARTIAL REGRESSION TARGETS



 Unlike Section <ref>, here we are concerned
with situations where the data is used both to select variables and
to determine how the relationship between selected variables and the
response is summarized. For example when the number of variables exceeds
the sample size, the full regression coefficient β^F is no
longer well defined so we cannot always summarize the effect of variable
j using β_j^F. We could still decide a priori to study
only the marginal relationship between x_j and y if j is
selected (e.g., the marginal correlation). In this case the data remains
uninvolved in target formation and we can continue to condition only
on j∈M̂. However, often, our goal is to understand how the variables
work in concert, to define a target that can meaningfully capture
key aspects of the joint distribution between the outcome and our
predictors.

In high dimensional settings, this necessarily means that we need
to use the data to inform our choice of target. Usually, this occurs when the lasso is used to select a set of variables M̂=M and then the effect of variable j∈ M is summarized by β_j^(M), where 

    β^(M)=(𝐗_M^⊤𝐗_M)^-1𝐗_M^⊤μ

is the regression of μ on the selected variables. These are the population partial regression coefficients, and we called this the partial target earlier.

Now, when variable j is selected, how we choose to summarize its effect depends on what other variables have been selected. It is insufficient to condition on { j∈M̂}; this event does not uniquely determine the quantity we are trying to infer. We need to know the identity of the other selected variables. In other words, two conditions must be satisfied for us to build a confidence interval for β_j^(M): (i) j∈M̂ and (ii) M̂=M (otherwise we might be making an interval for β_j^(M') with M'≠ M). Therefore, the conditional coverage guarantee we seek is 

    ℙ(β_j^(M)∉ C_j | j∈M̂, M̂=M)≤α.

Compared to (<ref>), the additional conditioning on
M̂=M is the cost of target formation. <cit.> showed how to construct such intervals with the following proposition.

  For η_j=𝐗_M(𝐗_M^⊤𝐗_M)^-1e_j,
	the conditional distribution η_j^⊤y | j∈M̂,M̂=M,ν_j is 𝒩(β_j^(M),σ^2‖η_j‖_2^2)
	truncated to a disjoint union of intervals ∪_s[l_s,u_s] where the union is over s∈{ -1,1} ^|M|.


The fact that the truncation region is a union of intervals is not
surprising (see Section <ref>); what is amazing is that
the interval endpoints l_s and u_s have explicit formulas
(see Appendix <ref>). Proposition <ref> should be compared
to Proposition <ref>. The truncation [-∞,a_j]∪[b_j,∞]
includes values of η_j^⊤y for which the lasso selects
a model N≠ M so long as j∈ N; the truncation to ∪_s[l_s,u_s]
further excludes such values of η_j^⊤y (it is a subset
of [-∞,a_j]∪[b_j,∞]). In
fact, even in moderate dimension problems (say p≈25), the
intervals [l_s,u_s] which fall within a practical
range, say [-10^2,10^2], are very few (intervals
falling outside this practical range contribute negligibly to p-value
calculations and can be numerically ignored). This is because M̂
is quite variable even in moderate dimensions. As we perturb the value
of η_j^⊤y, it becomes quite likely for one of the p
variables to enter or leave the active set.

All this is to say that when we condition on M̂=M, the range
of our test statistic β̂_j^(M)=η_j^⊤y
is quite severely truncated. This truncation leads to a loss in power,
hence longer intervals. In practice, we find that intervals conditioning
on M̂=M are often unacceptably long. This says that the cost
of target formation is particularly high (relative to variable selection
using the lasso) and we must exercise great caution in how we let
the data inform our queries. As evidence of this claim, one can refer
back to our prostate example (Figure <ref>) and the
simulation results for the full target (Figure <ref>).
There, we see that whether we condition on j∈M̂ or M̂=M
makes a huge difference (the discrepancy is in fact the cost of target
formation!); on the other hand, whether one conditions on the signs
ŝ makes much less of a difference.




 §.§ Stable Target Formation
 


The above discussion suggests that we cannot afford to let the data overly influence our choice of target. In particular,
choosing to summarize the effect of variable j using β_j^M̂
where M̂ is chosen by the lasso can lead to uselessly long
intervals. There are two simple remedies for limiting the role of
the data in target formation: 

	
  * Data Splitting: We can split the data and choose our target based on only a part of the data. 
	
  * Randomization: We can add some noise to the data and choose our target based on the corrupted data. 

Both approaches can lead to drastically shorter intervals. But they
also share some downsides: 

	
  * With a different split of the data (or different noise instance), not only will our intervals be different but the very targets of our inference may change. 
	
  * The model selected from a sample split or on corrupted data can be much worse than the model selected on the full data. We may get shorter intervals but we are also more likely to be asking the wrong questions!

For these reasons, we seek a deterministic approach for limiting the
cost of target formation. As discussed, the problem with summarizing
our effects using β_j^M̂ is the instability of M̂.
If the target of our research changes too often, that is probably not
a good thing!

Intuitively, we want to use a regularized/more stable version of M̂
to choose our target. We can think of M̂ as containing two
types of variables: those with very strong signals and those with
marginal signals (whose significance is ambiguous). We will call the
first type of variable high value targets (these are variables
we would not want to miss out on) and the second type of variable
low value targets. When we perturb the data, it is the low
value targets that are susceptible to entering or leaving the lasso
active set. The reason why target formation is so costly is because
we have allowed low value targets to participate equally as high value
targets in determining our target. The solution seems straightforward
then: allow only the high value targets to influence our choice of
target.

We will choose a subset Ĥ⊂M̂ of high value targets (see details below). How we choose to summarize the effect of a variable j∈M̂ depends on whether j is a high value target: 

	
  * High Value: We summarize the effect of j using β_j^Ĥ
	where 
	
    β^Ĥ=(𝐗_Ĥ^⊤𝐗_Ĥ)^-1𝐗_Ĥ^⊤μ.

	So our choice of target is fully adaptive for high value targets. 
	
  * Low Value: If variable j is selected by the lasso but is not deemed a high value target, we summarize its effect via β_j^Ĥ∪{j} where
	
    β^Ĥ∪{j}=(𝐗_Ĥ∪{j}^⊤𝐗_Ĥ∪{j})^-1𝐗_Ĥ∪{j}^⊤μ

	and 𝐗_Ĥ∪{j} is the matrix containing the high value targets as well as variable j. The coefficient β^Ĥ∪{j}_j is the effect of variable j after partialing out the effect of the high value targets, i.e., it allows us to ask the question whether variable j contributes any explanatory power beyond the variables in Ĥ. 

One may be worried whether we have greatly sacrificed interpretability by defining our inferential targets in this way. For a variable j, we seek its regression coefficient after partialing out the effect of variables in Ĥ; in contrast, the usual “partial” target seeks the regression coefficient of j after partialing out the effect of variables in M̂. Note that these regression coefficients should be quite similar in magnitude if (i) our low value targets
in M̂∖Ĥ have small coefficients or (ii) if the active variables in M̂ are not too highly correlated. But of course, if we define “high value” targets appropriately, then (i) should occur naturally; and the fact that the lasso usually selects only one out of a group of highly correlated variables makes (ii) likely as well. Thus, we should expect expect the inferential target we define
above to be quite close to the familiar partial regression coefficient based on M̂; re-examining again the prostate data (Figure <ref>), this turns out to be the case (the MLEs for the targets defined above are very close to the MLEs of the partial regression coefficients).


The chief effect of the above construction is that our choice of target,
denoted θ^Ĥ, depends not on M̂ but on Ĥ
which is much more stable. We construct a confidence interval for
θ_j^H whenever (i) j∈M̂ and (ii) Ĥ=H.
Hence the conditional confidence guarantee that we seek is 

    ℙ(θ_j^H∈ C_j | j∈M̂,Ĥ=H)≤α.

Our notation serves to emphasize that: 

	
  * The less stable M̂ is used for variable selection, i.e., choosing the variables whose effects we want to summarize. 
	
  * The more stable Ĥ is used for target formation, i.e., determining how the effects of the chosen variables are to be summarized. This makes sense since target formation is the more costly component of hypothesis generation. 

We will see that by stabilizing the target formation process, the
resulting intervals are much shorter especially the intervals for
high value targets. 




 §.§ How to Define a High Value Target?



From among the selected variables M̂, how should we choose
the subset of high value targets Ĥ? We introduce two methods:
the stable-ℓ_1 approach and the stable-t approach. In fact,
there are many possible proposals for creating a more stable version
of M̂. We found the stable-ℓ_1 and stable-t approaches
to be a good compromise between interpretability and computational
feasibility and one that we hope will appeal to a wide group of practitioners. However, it should be clear that our techniques can be easily modified to suit other definitions of stability.



Stable-ℓ_1: We saw in Section <ref>
that the _Ms intervals perform better when the lasso penalty
parameter is high (say near the universal threshold); this is because
the active set tends to be more stable for larger penalties. Thus,
a natural approach to constructing Ĥ is simply to let it
be the lasso active set but with a higher penalty parameter, λ_high,
than was used to select M̂. In this case, there is no guarantee
that Ĥ⊂M̂ (though it often is) but the above development does not require this. The proposal is also attractive because as the following proposition shows, it is computationally no more difficult than conditioning on M̂=M. 

 
The conditional distribution of η^⊤y given (i) j∈M̂, (ii) Ĥ=H, and (iii) P_η^⊥y is the 𝒩(η^⊤μ,σ^2‖η‖_2^2)
	distribution truncated to a set of the form 
	
    ([-∞,a]∪[b,∞])∩(∪_s[l_s,u_s]).

	The endpoints a,b are determined by the constraint j∈M̂
	and the union is over all sign vectors s∈{-1, 1} ^|Ĥ| such that (H,s) is compatible with the λ_high lasso solution for some element of { y':y'=η/‖η‖_2^2· z+P_η^⊥y,z∈ℝ}.
 
The proposition is a simple consequence of propositions
<ref> and <ref>. Conditioning on Ĥ=H
truncates η^⊤y to a union of intervals. The larger the
penalty parameter corresponding to Ĥ, the larger this union
of intervals tends to be, i.e., the less severe the truncation. We
usually choose λ_high close to the universal threshold
value of √(2log p/n).



Stable-t: Our second approach is motivated by the fact
that while the Bonferroni correction does not guarantee valid coverage
in all situations, it usually works pretty well. This suggests taking
Ĥ to be those variables in M̂ with t-statistics
surpassing a Bonferroni corrected threshold. We first fit a OLS model
using all the variables in M̂, i.e., 

    β̂^M̂=(𝐗_M̂^⊤𝐗_M̂)^-1𝐗_M̂^⊤y

and allow j to be a high value target if the t-statistic for
β̂_j^M̂ is large, i.e., if 

    |β̂_j^M̂/σ(𝐗_M̂^⊤𝐗_M̂)_jj^-1|>c

for some cutoff c. If we choose c by Bonferroni, it has the
form |Φ^-1(α/2p)|≈√(2log p) for large p; for such a choice, the stable-t method behaves similar to the stable-ℓ_1 method (with λ_high=√(2log p/n)) in many scenarios. The Bonferroni corrected cutoff ensures that we remove most variables whose importance is ambiguous. Once again, this proposal is attractive because it is computationally no more costly
than conditioning on M̂=M. 

 The conditional
	distribution of η^⊤y given (i) j∈M̂, (ii) Ĥ=H, and (iii) P_η^⊥y is the 𝒩(η^⊤μ,σ^2‖η‖_2^2)
	distribution truncated to a set of the form 
	
    ([-∞,a]∪[b,∞])∩(∪_k[l_k,u_k])

	where the endpoints l_k and u_k have explicit formulae.
 

 Below, all statements are implicitly conditioned on P_η^⊥y which restricts the data to a line. Note that Ĥ = Ĥ (M̂) is a function of M̂ and that

    {Ĥ(M̂) = H, j ∈M̂} = ⋃_M: j ∈ M{Ĥ(M) = H, M̂ = M}.

We already know {M̂ = M} is a union of intervals whose endpoints have known formulae. To finish characterizing the LHS of the above expression, it then suffices to show that {Ĥ(M)=H} is a union of intervals and compute the interval endpoints.

Fix M. Let us write 

    y'=cz+ν     c=η/‖η‖_2^2    ν=P_η^⊥y

and regress y' on 𝐗_M . The t-statistic for
the kth regression coefficient has the form 

    t_k=y'^⊤x̃_k/σ    x̃_k=𝐗_M(𝐗_M^⊤𝐗_M)^-1e_k/‖𝐗_M(𝐗_M^⊤𝐗_M)^-1e_k‖_2.

Note that |t_k|<ξσ if |z· c^⊤x̃_k+ν^⊤x̃_k|<ξσ, or equivalently 

    z∈[c_k,d_k]     c_k=-sign(c^⊤x̃_k)·ν^⊤x̃_k-ξσ/|c^⊤x̃_k|     d_k=-sign(c^⊤x̃_k)·ν^⊤x̃_k+ξσ/|c^⊤x̃_k|.

Similarly, |t_k|>ξσ if 

    z∈[-∞,c_k]∪[d_k,∞].

For Ĥ(M)=H, we require |t_k|>ξσ for all
k∈ H and |t_k|<ξσ for all k∈ M∖ H.
Thus we require 

    z∈(k∈ H⋂([-∞,c_k]∪[d_k,∞]))∩(k∈ M∖ H⋂[c_k,d_k])

which simplifies to a union of intervals.  

	


 §.§ Prostate data revisited



The right panel of Figure <ref> shows the result of applying the stable-t approach to the prostate data. Three high-value targets have been identified,  and for all targets, the  _stab-t  intervals are shorter than the _Ms intervals of <cit.>.




§ SIMULATIONS FOR THE  PARTIAL REGRESSION TARGET


We test our proposals across a variety of simulation settings; we
focus on the case where the data informs both variable selection and
target formation (for examples involving only variable selection,
see Section <ref>). 

	
  * We compare six methods: naive intervals, Bonferroni corrected, _Ms, _M,_stab-t and _stab-ℓ_1. Recall that _Ms condition on (M̂,ŝ) while _M conditions only on M̂. 
	
  * We consider four possibilities for (n,p): (100,50), (100,125), (100,250) and (100,1250). 
	
  * We generate the outcome as y_i=x_i^⊤β+ε_i, i=1,…,n. For the moment we let x_i∼𝒩(0,𝐈_p) and ε_i∼𝒩(0,1); however, results for the case of correlated predictors and non-normal errors are reported in Appendices <ref> and <ref>.
	
  * The first k=5 components of β are non-zero, the rest are set to 0. The non-null components are all set to either δ_low or δ_high. For each (n,p) configuration, we chose δ_low and δ_high based on the null distribution of max_jn^-1|x_j^⊤y| following the rules detailed in Section <ref>. 
	
  * Following Section <ref>, we consider two values of λ: one based on CV and one based on the universal threshold value of √(2log p/n). 
	
  * For the stable-ℓ_1 method, we set λ_high=√(2log p/n) when λ is based on CV and λ_high=1.25λ when λ is already chosen as the universal threshold. 
	
  * For the stable-t method, we set c=|Φ^-1(α/2p)| where α=0.1 (chosen because we are building 90% confidence
	intervals) 

The results for the cases n=100,p=250 and n=100,p=1250 are presented in Figures <ref> and <ref> respectively (results for n=100,p=50 and n=100,p=125 are given in Appendix <ref>); Figure <ref> presents the null case (β=0). We summarize the key takeaways: 

	
  * Bonferroni adjusted intervals have poor coverage (only 47% for 90% confidence interval) when the global null holds.
	
  * The _M method improves on the _Ms intervals but as the simulations demonstrate they are often many times wider than the naïve or Bonferroni-corrected intervals. In practice, the inferences they deliver would not be useful. This suggests that we need to go beyond minimal conditioning (which in the case of partial targets is M̂) and further restrict how data is used for target formation. 
	
  * The stable-ℓ_1 and especially the stable-t method produce intervals that are often only slightly longer than the Bonferroni corrected intervals (and in high signal settings can be quite a bit shorter than Bonferroni). In addition, they tend to be more stable numerically, i.e., less chance of an “infinite” length interval.
	
  * As we saw in Section <ref>, a higher value of λ (i.e., using the universal threshold rather than cross validation) helps to make M̂ more stable and hence reduce the length of the _Ms and _M intervals. However, we see that in truly high dimensional problems such as when n=100 and p=1250, the stable-t method is half the length of the _M method
	even when λ is selected to be the universal threshold. In addition, many practitioners prefer choosing λ by cross-validation in which case use of the stable methods is a necessity if we want reasonable interval lengths. 

In addition to the above results for independent predictors, we also
considered the case of correlated predictors. We explored block-equicorrelation and Toeplitz covariance structures for the predictors with the maximum pairwise correlation set at 0.5 (see Appendix <ref> for details and results).
What we observed was that the stable methods outperform the _M and _Ms intervals by an even larger margin in the presence of strongly correlated predictors.  Intuitively as we increase the correlation between predictors, the set of lasso active variables M̂ becomes less stable (as we move in the direction η_j capturing
the effect of x_j, we are likely to flip the active/inactive
status of variables correlated with x_j), hence increasing the
need for stabilization.

It has become dogma that we need to regularize parameter estimates
in high dimensions. Similarly, we hope that these simulations argue
clearly for the need to regularize/stabilize how we choose targets
for inference. It is usually not advisable to form a model using all
variables in the lasso active set even if it is technically possible
to account for such selection when doing inference.










 §.§ Performance When Assumptions Are Violated



In Appendix <ref>, we investigate whether the good performance of the proposed methods are robust to various violations of our assumptions. We consider non-normality of the errors, the effect of having to estimate
σ^2, and the effect of a data driven choice of λ.
We find that the proposed TZ intervals are robust to various model
violations, and the lengths of the stable intervals seem less likely
to be negatively impacted by model violations than the lengths of
the _M and _Ms intervals.




§ PREDICTING HIV DRUG RESISTANCE
  

The target of HIV drugs often becomes resistant through mutation. For six nucleotide reverse transcriptase inhibitors
(NRTIs) used to treat HIV, <cit.> investigate which mutations are most predictive of drug resistance. We focus on the drug Lamivudine (3TC). The primary response is log susceptibility to drug. The explanatory variables are binary indicators for whether a mutation has occurred at a particular site. The data have n=1057 samples and p=210 sites.

We fit a lasso model with λ=0.0124 (chosen via 10-fold cross
validation), yielding 36 active mutations. One of the mutations, p184, had a particularly large effect: its OLS estimate from a regression on all the variables is 0.477 with the accompanying naive 90% confidence interval being (0.471, 0.485). TZ_V, TZ_M and TZ_Ms all produced the same selection adjusted interval (0.437, ∞);
the infinite upper bound here is due to numerical inaccuracy from
inverting a truncated Gaussian pivot. When considering the partial
coefficient for p184, the results were qualitatively identical (the
truncated-Z intervals still had an infinite upper bound). The effect
estimates for the remaining 35 active variables are an order of magnitude smaller and their 90% confidence intervals are plotted in Figure <ref>. Here are the key points:

	
  * The results are consistent with the analysis of the prostate data
	(see Figure <ref>). The TZ_V and TZ_stab-t
	intervals tend to be roughly half as long as the TZ_M intervals
	and only a bit longer than the naive intervals. TZ_M and TZ_Ms
	frequently produce very long intervals.
	
  * Using the stable-t criterion (with threshold set to |Φ^-1(0.1/2p)|=3.5),
	we selected 6 high value targets: p65, p67, p69, p90, p184 (not pictured),
	and p215. With the exception of p184 discussed above, the TZ_V
	and TZ_stab-t methods yielded very stable intervals for
	the effects of these high value targets; despite the strong signal
	from these 6 mutations, TZ_M and TZ_Ms yielded very imprecise
	intervals for p67 (interval contains 0), p215, and p69 (in the case
	of full).
	
  * The target that the TZ_stab-t intervals attempt to cover
	is different from the partial regression coefficients (defined with
	respect to active set); see Section <ref>. But
	as the bottom panel of Figure <ref> shows, they are very
	similar in practice. This is because, as we have discussed above, the
	two targets differ significantly only if the low value targets have
	large coefficients (and are correlated with the high value targets).
	
  * TZ_stab-t will sometimes yield longer intervals than TZ_M
	and TZ_Ms (e.g., p215) but this occurs very rarely.


 




§ A GENERAL RECIPE
  


For concreteness, this paper has focused on variable selection and
target formation using the lasso. However, many of the ideas are applicable
to any model selection rule with polyhedral selection regions or even
non-polyhedral selection regions. What changes in the general case
is that the computation cost of the truncation region for the Z-statistic
can be much greater. Let us begin with the ingredients that do not
change:

	
  * We are interested in inference on linear functions μ^⊤η_j^Ĥ of the mean μ for variables j∈M̂. Here M̂ is the set of variables we are interested in, while Ĥ is (possibly) another set of variables used to define our inferential target (through η^Ĥ). Above, we saw three cases (i) Ĥ={ 1,…,p} and M̂ is the lasso active set (full targets) (ii) Ĥ=M̂ is the lasso active set (partial targets) and (iii) M̂ is the lasso active set and Ĥ is the subset of M̂ with large t-statistics (stable-t).
	
  * For η=η_j^Ĥ, we condition our inferences
	on P_η^⊥y=(𝐈_n-ηη^⊤/‖η‖_2^2)y; see <cit.> for ways to remove this conditioning. That is, we first restrict our data to the line 
	
    ℒ={ y':y'=P_η^⊥y+z·η/‖η‖_2^2, z∈ℝ} .

	This line represents perturbations of the data in the direction η only. The line is parametrized by z and note that for y'=y, we have z=y^⊤η∼𝒩(η^⊤μ,σ^2η_2^2) is the Z-statistic (up to standardization of variance). 
	
  * For η=η_j^Ĥ, we further condition
	on j∈M̂ and Ĥ=H to adjust for the costs variable
	selection and target formation respectively. The effect of this conditioning
	is to restrict the line ℒ to a range of z-values (i.e.,
	we obtain a truncated Z-statistic); we can necessarily express the
	intersection of ℒ with the sets j∈M̂
	and Ĥ=H as a union of intervals (of z-values).

The recipe above is quite general: we made no restrictions on what
kind of procedure can be used to construct M̂ or Ĥ.
Of course, for complicated M̂ and Ĥ, finding the
range of values that the Z-statistic has been truncated to can be
difficult (if not impossible). We have seen two particularly easy
to compute situations:

	
  * TZ_Ms intervals: Letting M̂=Ĥ be the lasso active set and also conditioning on the sign of the active variables, we obtain a polyhedral region, which when intersected with the line ℒ
	results in an interval. The power of using “polyhedral” selection rules is that it greatly simplifies computation for the truncation interval: in particular, once we have expressed the selection region as {𝐀y≤ b}, formulas in Appendix <ref> give the interval endpoints. Another common polyhedral selection rule is forward stepwise regression. If M̂ is constructed using K steps of forward stepwise regression, then the set {M̂=M,ŝ=s} ={𝐀y≥ 0} where ŝ are the signs of the K-selected variables; see <cit.> for the formula for 𝐀. 
	
  * Full target: When M̂ is the lasso active set and H={1,…,p}, then the set of allowable z-values is of the form (-∞,a)∪(b,∞) where a and b have simple formulas (see Proposition <ref>). 

Outside these special cases, simple formulas usually do not exist
for characterizing the intersection of ℒ with { j∈M̂}
and {Ĥ=H}. Nonetheless, one can usually run
a pathwise algorithm to determine the truncation region. Let 𝒵={ z_1,…,z_G}
be a grid of z-values (which map to points on ℒ).
For each z_k, we compute M̂(z_k) and Ĥ(z_k).
This yields a subset 𝒵_trunc⊂𝒵
for which j∈M̂(z_k) and Ĥ(z_k)=H.
We then approximate 𝒵_trunc using a union of
intervals and this becomes our truncation set. We note that:

	
  * For algorithms such as the lasso, we can take advantage of warm starts.
	That is, having obtained the lasso solution for z_k, we can use it as an initial starting point for z_k+1. In fact, for the lasso, warm starts are not really necessary, because one can examine the constraint {𝐀y≤ b} characterizing the selection region for the lasso to determine which variables will enter/leave
	the active set at the endpoints of any particular interval of z-values characterized by a (M̂,ŝ) pair.
	
  * Very positive and very negative z-values do not affect inference. Hence, it is usually reasonable to consider only grid values less than
	say 20σ.

The above recipe in fact allows us to adjust for selection using any blackbox model selection procedure—all that is needed is that we be able to run this blackbox procedure for all z-values in some grid. Again, the limiting factor here is computation: how quickly we can re-evaluate our model selection procedure on all grid values.




§ FURTHER TOPICS



Scientific research in the modern era has become much less “pre-specified” and much more dependent on exploratory data analysis to generate interesting hypotheses and models. Accounting for this exploration is crucial
for producing replicable results. Unfortunately, the simplest selection adjusted intervals, as in <cit.> and <cit.>, are often too long to be useful—they “over-condition,” that is they over-state the amount of information used up for hypothesis generation. A more refined accounting of costs reveals two primary ways through which the data guides many investigations:

	
  * Selecting interesting variables for further investigation.
	
  * Choosing a “parametrization” of the scientific problem so that, for example, the effect of a variable can now be captured by some coefficient in a chosen model (target formation).

Decoupling these two activities is crucial because the inferential
cost incurred by target formation is much larger than that incurred by variable selection. For example, in the prostate and HIV datasets (as well as across our simulations), we saw that the TZ_V intervals (adjusts for cost of variable selection) are much shorter than the TZ_M intervals (adjusts for both variable selection and target formation). 

There is of course an extreme response to this asymmetric cost scheme: fix your targets beforehand and use the data only for variable selection. This, however, is too limiting. The appropriate response is to exert greater control over the cost of target formation. This can be done, for example, by splitting the data in two and choosing our targets
based on only half the data; <cit.> develop more sophisticated/efficient randomization approaches to controlling the cost. These approaches,
however, simultaneously use less information to choose our variables and less information to form our targets. But given that the cost of variable selection is not very high, it surely makes sense to use all the data to choose our important variables (hence arriving at a better set of “important” variables) and decrease only the cost of target formation. Rather than using a single lever to control the
cost of data exploration, an innovation of this paper is to use to
introduce two levers. Concretely, we found that selecting variables using the lasso but forming inferential targets based on a more stable set of variables, yielded good practical results —“good” here means intervals that are slightly longer than the naive intervals
but significantly shorter than the TZ_M intervals (frequently,
less than half the length). 

Our results have relied on several assumptions, e.g., λ is
fixed. In the regime where p is fixed and n is growing, these
assumptions can be dispensed with:

	
  * <cit.> show that the truncated normal pivot
	underlying the TZ_Ms method holds asymptotically for a wide class
	of (non-normal) error distributions. The arguments therein continue
	to hold for the other truncated-Z approaches in this paper.
	
  * <cit.> and <cit.> show that if an estimator
	of σ^2 is consistent unconditionally, then it remains consistent
	conditional on selection (by the lasso). In the p fixed regime,
	this implies that the pairs bootstrap can be used to consistently
	estimate various covariances.
	
  * λ can be chosen via randomized cross validation. The details are given in <cit.>; the key is that we continue to
	enjoy (asymptotically) a joint multivariate normality between (i)
	the test statistic, (ii) the statistics characterizing the lasso active set and (iii) the randomized cross validation errors (over a grid of λ values). The multivariate normality then allows us to derive the conditional distribution of the test statistic.


When p>n, our theoretical understanding is limited, but empirically,
the truncated-Z intervals seem to work reasonably well. 
Generalizations of this approach to (regularized- ) generalized linear models and the Cox PH model would also be interesting and
useful. The new procedures presented in this paper will be implemented in our  package in the CRAN repository. 



Acknowledgements: The authors would like to thank Jonathan Taylor, who has been the intellectual leader in the post-selection inference effort at Stanford. 




§ APPENDICES

tocsectionAppendices




 §.§ Formulae for the  truncation points from the polyhedral lemma
 

Suppose that a selection region can be characterized by {𝐀y ≤ b}. Let us decompose y as:

    y=z· c+ν     z=η^⊤y     c=η/‖η‖_2^2.

For fixed 𝐀,b and ν, z is truncated to

    Ν^-(y) ≤ z ≤Ν^+(y)


    Ν^-(y)   =max_j:(𝐀c)_j<0b_j-(ν)_j/(𝐀c)_j
    Ν^+(y)   =min_j:(𝐀c)_j > 0b_j-(ν)_j/(𝐀c)_j.

In addition, for ν to be compatible with 𝐀 and b, we require Ν^0(y) ≥ 0 where

    Ν^0(y) =min_j:(𝐀c)_j = 0 (b_j-(ν)_j).

The derivation and additional details can be found in Lemma 5.1 of <cit.>.




 §.§ A simple problem illustrating the conditioning events in detail
 


Here we  show the precise effects of different conditioning strategies in a simple example.  
Suppose n=2, p=2, and 𝐗=([        1        ρ;        0 √(1-ρ^2) ]). Note that each column of 𝐗 has unit norm and that x_1^⊤x_2=ρ
where x_j denotes the jth column. The solid pink
lines in Figure <ref> denote the boundaries of the different
lasso selection regions in y-space (i.e., ℝ^2).
For example, if y falls in the central parallelogram region
(whose boundaries are given by x_1^⊤y=±λ
 and x_2^⊤y=±λ), then neither variable will be active in the lasso solution. Consider:

	
  * Suppose y equals the red dot in Figure <ref>
	labeled “Case 1.” Then the lasso active set is M̂={ 1}
	and the sign of the active variable is ŝ=-1. As a result of
	this initial exploration, we settle on testing the hypothesis H_0:β^(M̂)=μ_1=0.
	Here η=(1,0)^⊤ and P_η^⊥y=y_2.
	Assuming σ=1, our Z-statistic is simply y_1. If we condition
	on (i) M̂={ 1}, (ii) ŝ=-1 and (iii) y_2,
	then we know that y_1 must lie in (-∞,-λ). So the
	usual 𝒩(0,1) distribution of the Z-statistic is now
	𝒩(0,1) truncated to (-∞,-λ).
	
  * As discussed above, we conditioned on ŝ to make the
	math simpler. If we condition only on M̂={ 1}
	and y_2, we still get that y_1∈(-∞,-λ).
	However, if y were as in Case 2 of Figure <ref>,
	conditioning on both M̂ and ŝ leads to y_1∈(λ,∞)
	while conditioning only on M̂ yields y_1∈(-∞,-λ)∪(λ,∞),
	a union of intervals. That is, the less we condition on, the less
	severe the truncation. 
	
  * In Section <ref> we found it useful to condition on even less than just M̂.
	In particular, in some situations, we may simply want to condition
	just on the fact that variable 1 was selected, 1∈M̂. In
	such a case M̂ can equal { 1} or {1,2}.
	In Case 1, conditioning on 1∈M̂ and y_2 restricts y_1
	to (-∞,-λ)∪(λ,∞). 







 §.§ Further simulations for the partial target
 


Figures <ref> and <ref> show results for the cases n=100,p=50 and n=100,p=125 (the simulation settings are described in Section <ref>).








 §.§ Performance of the proposed methods for correlated predictors
 


We study whether our proposed methods continue to work well
when the predictors are correlated. As in the original simulations, there
are k=5 non-null variables but now they may be heavily correlated
with other variables. We consider two correlation schemes:



  * Equi-correlation blocks: Divide all variables among five blocks of
size p/5. There is one non-null variable in each block and it has
a 𝒩(0,1) distribution. The other p/5-1 variables
in each block are generated from the non-null variable in that block
as follows:

    x_j=ρ x_non-null+√(1-ρ^2)ε    ε∼𝒩(0,1).

We set ρ=0.5.

  * Toeplitz: The covariance matrix of the predictors has Toeplitz form
with the first row of the covariance matrix having the form

    (1,ρ,ρ^2,…).

Five of the predictors are randomly chosen to be the non-null predictors.
We set ρ=0.5.


All simulations are performed with n=100 and p=250. Figures <ref> and <ref> summarize
the results. The stable methods continue to outperform the _M and
_Ms intervals in the presence of strongly correlated predictors. In fact, we see that the stable-t method shines even more than in the case of independent
predictors—in particularly, it significantly outperforms
the _M and _Ms intervals even when λ is set to be the universal
threshold value. This performance advantage is due to the fact that the active set of the lasso becomes less stable when variables are correlated.








 §.§ Performance of the proposed methods when assumptions are violated
 


We now investigate whether the good performance of the stable methods
are robust to various violations of our assumptions. In particular,
we consider:

	
  * Non-normal errors: we consider the effect of heavy tails (t_3
	distribution) and skewness (skew normal with skewness coefficient
	10) in our error distribution. See Figures <ref> and <ref>.
	
  * Unknown σ^2: Estimating the noise level when p>n is a
	challenging problem. However, <cit.> shows
	that 
	
    σ̂^2=1/n-p_λ̂∑_i=1^n(y_i-x_i^⊤β̂_λ̂)^2

	performs well where β̂_λ̂ is the lasso solution,
	p_λ̂ is the number of active variables, and λ̂
	is chosen by 10-fold cross validation. We then plug-in σ̂^2
	in place of the unknown noise level. See Figure <ref>.

	
  * Data driven choice of λ: the  methods adjust for
	the effect of selection assuming a pre-specified λ; they
	do not account for the fact that λ is often chosen based
	on the data, e.g., via cross-validation. See Figure <ref>.

For all simulations, we use n=100 and p=250. We summarize the
results as follows:

	
  * All the truncated-Z intervals maintained nominal coverage close to 0.9 inspite
	of heavy tails and skewness. Qualitatively, the improvements of the
	stable methods over the _M and _Ms intervals are the
	same as in the case of normal errors.
	
  * When σ^2 was unknown, all truncated-Z methods maintained
	coverage close to the nominal level by plugging in an estimate σ̂^2
	(see Figure <ref>). However, the lengths of
	the _M and _Ms intervals are much longer than in the
	case of known σ^2 (cf. Figure <ref>);
	the lengths of the stable methods were less impacted.
	
  * When λ was chosen by 10-fold cross validation rather than
	set at a fixed value, the truncated-Z methods still maintained coverage
	close to the nominal level of 0.9 despite ignoring the cost of estimating
	λ. However, just as in the case of unknown σ^2,
	the lengths of the _M and _Ms intervals greatly increased
	while the lengths of the stable intervals were relatively immune.

What these simulations suggest is that (1) the truncated-Z interval
seem to possess coverage properties that are robust to various model
violations and (2) the lengths of the stable intervals are less likely
to be negatively impacted by model violations than the lengths of
the _M and _Ms intervals.









agsm

 