In Defense of the Indefensible: A Very Naive Approach to High-Dimensional Inference

TL;DR

This paper demonstrates that a simple two-step procedure using lasso followed by least squares can produce valid high-dimensional inference under certain conditions, challenging conventional wisdom about data peeking.

Contribution

It shows that the variable set selected by lasso is deterministic under assumptions, enabling asymptotically valid inference with naive confidence intervals and score tests.

Findings

Lasso-selected variables are often deterministic under assumptions.

Naive confidence intervals are asymptotically valid for selected coefficients.

Naive score tests can reliably test hypotheses in high-dimensional models.

Abstract

A great deal of interest has recently focused on conducting inference on the parameters in a high-dimensional linear model. In this paper, we consider a simple and very na\"{i}ve two-step procedure for this task, in which we (i) fit a lasso model in order to obtain a subset of the variables, and (ii) fit a least squares model on the lasso-selected set. Conventional statistical wisdom tells us that we cannot make use of the standard statistical inference tools for the resulting least squares model (such as confidence intervals and $p$-values), since we peeked at the data twice: once in running the lasso, and again in fitting the least squares model. However, in this paper, we show that under a certain set of assumptions, with high probability, the set of variables selected by the lasso is identical to the one selected by the noiseless lasso and is hence deterministic. Consequently, the na\"{i}ve two-step approach can yield asymptotically valid inference. We utilize this finding to develop the \emph{na\"ive confidence interval}, which can be used to draw inference on the regression coefficients of the model selected by the lasso, as well as the \emph{na\"ive score test}, which can be used to test the hypotheses regarding the full-model regression coefficients.

Figures

Tables

Table 1. Table 1: Coverage proportions (Cov) and average lengths (Len) of 95% naïve confidence intervals with tuning parameters λ sup subscript 𝜆 sup \lambda_{\text{sup}} and λ 1SE subscript 𝜆 1SE \lambda_{\text{1SE}} , and 95% exact post-selection confidence intervals under the scale-free graph setting with partial correlation ρ ∈ { 0.2 , 0.6 } 𝜌 0.2 0.6 \rho\in\{0.2,0.6\} , sample size n ∈ { 100 , 300 , 500 } 𝑛 100 300 500 n\in\{100,300,500\} , dimension p = 100 𝑝 100 p=100 and signal-to-noise ratio SNR ∈ { 0.1 , 0.3 , 0.5 } SNR 0.1 0.3 0.5 \mathrm{SNR}\in\{0.1,0.3,0.5\} . The error variance is estimated either through the scaled lasso (Sun and Zhang,, 2012 ) (SL) or Theorem  2.4 (NL).

	$\rho$	0.2
	$n$	100			300			500
	SNR	$0.1$	$0.3$	$0.5$	$0.1$	$0.3$	$0.5$	$0.1$	$0.3$	$0.5$
Cov	exact ${\lambda }_{\text{sup}}$ SL	0.905	0.942	0.936	0.949	0.953	0.953	0.946	0.949	0.949
	naïve ${\lambda }_{\text{sup}}$ SL	0.683	0.961	0.944	0.951	0.952	0.951	0.971	0.948	0.947
	naïve ${\lambda }_{\text{1SE}}$ SL	0.596	0.876	0.887	0.905	0.917	0.910	0.944	0.935	0.932
	naïve ${\lambda }_{\text{1SE}}$ NL	0.587	0.870	0.883	0.905	0.919	0.911	0.945	0.939	0.933
Len	exact ${\lambda }_{\text{sup}}$ SL	4.260	1.634	0.823	1.907	0.437	0.330	0.814	0.327	0.255
	naïve ${\lambda }_{\text{sup}}$ SL	1.122	0.693	0.552	0.717	0.421	0.326	0.562	0.325	0.253
	naïve ${\lambda }_{\text{1SE}}$ SL	1.199	0.705	0.555	0.718	0.420	0.326	0.560	0.326	0.253
	naïve ${\lambda }_{\text{1SE}}$ NL	1.201	0.711	0.559	0.723	0.424	0.329	0.563	0.328	0.255
	$\rho$	0.6
	$n$	100			300			500
	SNR	$0.1$	$0.3$	$0.5$	$0.1$	$0.3$	$0.5$	$0.1$	$0.3$	$0.5$
Cov	exact ${\lambda }_{\text{sup}}$ SL	0.956	0.938	0.945	0.922	0.944	0.945	0.948	0.945	0.947
	naïve ${\lambda }_{\text{sup}}$ SL	0.625	0.943	0.953	0.942	0.942	0.942	0.964	0.944	0.944
	naïve ${\lambda }_{\text{1SE}}$ SL	0.643	0.873	0.876	0.918	0.934	0.929	0.954	0.941	0.935
	naïve ${\lambda }_{\text{1SE}}$ NL	0.639	0.870	0.873	0.919	0.934	0.929	0.957	0.944	0.938
Len	exact ${\lambda }_{\text{sup}}$ SL	6.433	1.616	0.770	1.892	0.429	0.328	0.787	0.326	0.253
	naïve ${\lambda }_{\text{sup}}$ SL	1.108	0.692	0.553	0.713	0.420	0.326	0.560	0.325	0.252
	naïve ${\lambda }_{\text{1SE}}$ SL	1.200	0.705	0.554	0.720	0.420	0.326	0.561	0.326	0.253
	naïve ${\lambda }_{\text{1SE}}$ NL	1.203	0.709	0.557	0.725	0.423	0.329	0.564	0.327	0.254

Table 2. Table 2: Coverage proportions (Cov) and average lengths (Len) of 95% naïve confidence intervals with tuning parameters λ sup subscript 𝜆 sup \lambda_{\text{sup}} and λ 1SE subscript 𝜆 1SE \lambda_{\text{1SE}} , and 95% exact post-selection confidence intervals under the stochastic block model setting with partial correlation ρ ∈ { 0.2 , 0.6 } 𝜌 0.2 0.6 \rho\in\{0.2,0.6\} , sample size n ∈ { 100 , 300 , 500 } 𝑛 100 300 500 n\in\{100,300,500\} , dimension p = 100 𝑝 100 p=100 and signal-to-noise ratio SNR ∈ { 0.1 , 0.3 , 0.5 } SNR 0.1 0.3 0.5 \mathrm{SNR}\in\{0.1,0.3,0.5\} . The error variance is estimated either through the scaled lasso (Sun and Zhang,, 2012 ) (SL) or Theorem  2.4 (NL).

	$\rho$	0.2
	$n$	100			300			500
	SNR	$0.1$	$0.3$	$0.5$	$0.1$	$0.3$	$0.5$	$0.1$	$0.3$	$0.5$
Cov	exact ${\lambda }_{\text{sup}}$ SL	0.921	0.934	0.948	0.951	0.934	0.933	0.953	0.958	0.958
	naïve ${\lambda }_{\text{sup}}$ SL	0.540	0.955	0.960	0.953	0.931	0.930	0.966	0.957	0.957
	naïve ${\lambda }_{\text{1SE}}$ SL	0.585	0.850	0.858	0.911	0.922	0.917	0.961	0.944	0.940
	naïve ${\lambda }_{\text{1SE}}$ NL	0.586	0.854	0.865	0.913	0.924	0.918	0.961	0.945	0.942
Len	exact ${\lambda }_{\text{sup}}$ SL	4.721	1.692	0.857	1.873	0.425	0.325	0.758	0.325	0.252
	naïve ${\lambda }_{\text{sup}}$ SL	1.111	0.682	0.545	0.709	0.416	0.323	0.559	0.323	0.251
	naïve ${\lambda }_{\text{1SE}}$ SL	1.188	0.703	0.552	0.715	0.416	0.323	0.557	0.323	0.251
	naïve ${\lambda }_{\text{1SE}}$ NL	1.195	0.710	0.558	0.721	0.420	0.326	0.560	0.325	0.253
	$\rho$	0.6
	$n$	100			300			500
	SNR	$0.1$	$0.3$	$0.5$	$0.1$	$0.3$	$0.5$	$0.1$	$0.3$	$0.5$
Cov	exact ${\lambda }_{\text{sup}}$ SL	0.923	0.945	0.949	0.959	0.956	0.958	0.943	0.948	0.948
	naïve ${\lambda }_{\text{sup}}$ SL	0.569	0.959	0.964	0.964	0.955	0.956	0.963	0.946	0.946
	naïve ${\lambda }_{\text{1SE}}$ SL	0.568	0.841	0.852	0.922	0.933	0.927	0.947	0.946	0.941
	naïve ${\lambda }_{\text{1SE}}$ NL	0.561	0.841	0.850	0.922	0.935	0.929	0.950	0.947	0.943
Len	exact ${\lambda }_{\text{sup}}$ SL	4.744	1.662	0.832	1.987	0.424	0.325	0.807	0.324	0.251
	naïve ${\lambda }_{\text{sup}}$ SL	1.119	0.690	0.545	0.712	0.416	0.323	0.557	0.322	0.250
	naïve ${\lambda }_{\text{1SE}}$ SL	1.194	0.704	0.554	0.713	0.416	0.323	0.555	0.322	0.250
	naïve ${\lambda }_{\text{1SE}}$ NL	1.194	0.707	0.557	0.718	0.420	0.326	0.558	0.324	0.252

Table 3. Table 3: Among repetitions that 𝒜 ^ λ b ≠ ∅ superscript subscript ^ 𝒜 𝜆 𝑏 \hat{\mathcal{A}}_{\lambda}^{b}\neq\emptyset , the proportion of 𝒜 ^ λ b superscript subscript ^ 𝒜 𝜆 𝑏 \hat{\mathcal{A}}_{\lambda}^{b} that equals the most common 𝒜 ^ λ b superscript subscript ^ 𝒜 𝜆 𝑏 \hat{\mathcal{A}}_{\lambda}^{b} , b = 1 , … , 1000 𝑏 1 … 1000 b=1,\dots,1000 , under the scale-free graph and stochastic block model settings with tuning parameters λ sup subscript 𝜆 sup \lambda_{\text{sup}} and λ 1SE subscript 𝜆 1SE \lambda_{\text{1SE}} . In the simulation, ρ ∈ { 0.2 , 0.6 } 𝜌 0.2 0.6 \rho\in\{0.2,0.6\} , sample size n ∈ { 100 , 300 , 500 } 𝑛 100 300 500 n\in\{100,300,500\} , dimension p = 100 𝑝 100 p=100 and signal-to-noise ratio SNR ∈ { 0.1 , 0.3 , 0.5 } SNR 0.1 0.3 0.5 \mathrm{SNR}\in\{0.1,0.3,0.5\} .

$\rho$	0.2
$n$	100			300			500
SNR	$0.1$	$0.3$	$0.5$	$0.1$	$0.3$	$0.5$	$0.1$	$0.3$	$0.5$
Scale-free; ${\lambda }_{\text{sup}}$	0.000	0.999	0.993	1.000	0.999	0.998	1.000	0.999	0.996
Scale-free; ${\lambda }_{\text{1SE}}$	0.832	0.958	0.851	0.961	0.948	0.934	0.988	0.980	0.966
Stochastic block; ${\lambda }_{\text{sup}}$	1.000	0.997	0.992	0.998	1.000	0.999	1.000	1.000	1.000
Stochastic block; ${\lambda }_{\text{1SE}}$	0.838	0.839	0.831	0.964	0.957	0.946	0.992	0.979	0.968
$\rho$	0.6
$n$	100			300			500
SNR	$0.1$	$0.3$	$0.5$	$0.1$	$0.3$	$0.5$	$0.1$	$0.3$	$0.5$
Scale-free; ${\lambda }_{\text{sup}}$	0.995	0.996	0.996	0.999	0.998	0.998	0.999	1.000	1.000
Scale-free; ${\lambda }_{\text{1SE}}$	0.854	0.869	0.855	0.964	0.947	0.935	0.994	0.981	0.962
Stochastic block; ${\lambda }_{\text{sup}}$	0.998	0.998	0.998	1.000	1.000	1.000	1.000	0.999	0.999
Stochastic block; ${\lambda }_{\text{1SE}}$	0.834	0.843	0.826	0.969	0.976	0.968	0.988	0.979	0.971

Table 4. Table 4: The proportion of 𝒜 ^ λ b superscript subscript ^ 𝒜 𝜆 𝑏 \hat{\mathcal{A}}_{\lambda}^{b} that equals 𝒜 λ b superscript subscript 𝒜 𝜆 𝑏 \mathcal{A}_{\lambda}^{b} , b = 1 , … , 1000 𝑏 1 … 1000 b=1,\dots,1000 , under the scale-free graph and stochastic block model settings with tuning parameters λ sup subscript 𝜆 sup \lambda_{\text{sup}} and λ 1SE subscript 𝜆 1SE \lambda_{\text{1SE}} . In the simulation, ρ ∈ { 0.2 , 0.6 } 𝜌 0.2 0.6 \rho\in\{0.2,0.6\} , sample size n ∈ { 100 , 300 , 500 } 𝑛 100 300 500 n\in\{100,300,500\} , dimension p = 100 𝑝 100 p=100 and signal-to-noise ratio SNR ∈ { 0.1 , 0.3 , 0.5 } SNR 0.1 0.3 0.5 \mathrm{SNR}\in\{0.1,0.3,0.5\} .

$\rho$	0.2
$n$	100			300			500
SNR	$0.1$	$0.3$	$0.5$	$0.1$	$0.3$	$0.5$	$0.1$	$0.3$	$0.5$
Scale-free ${\lambda }_{\text{sup}}$	0.958	0.725	0.928	0.649	1.000	1.000	0.971	0.999	0.998
Scale-free ${\lambda }_{\text{1SE}}$	0.417	0.665	0.829	0.567	0.955	0.940	0.804	0.976	0.964
Stochastic block ${\lambda }_{\text{sup}}$	0.964	0.736	0.931	0.678	1.000	1.000	0.964	1.000	1.000
Stochastic block ${\lambda }_{\text{1SE}}$	0.394	0.661	0.815	0.582	0.964	0.957	0.582	0.559	0.964
$\rho$	0.6
$n$	100			300			500
SNR	$0.1$	$0.3$	$0.5$	$0.1$	$0.3$	$0.5$	$0.1$	$0.3$	$0.5$
Scale-free ${\lambda }_{\text{sup}}$	0.952	0.713	0.924	0.687	1.000	1.000	0.968	0.999	0.998
Scale-free ${\lambda }_{\text{1SE}}$	0.396	0.657	0.799	0.572	0.950	0.941	0.779	0.976	0.963
Stochastic block ${\lambda }_{\text{sup}}$	0.950	0.750	0.935	0.687	1.000	1.000	0.967	1.000	0.999
Stochastic block ${\lambda }_{\text{1SE}}$	0.445	0.661	0.806	0.559	0.959	0.948	0.805	0.987	0.974

Table 5. Table 5: Average power and type-I error rate for the hypotheses H 0 , j ∗ : β j ∗ = 0 : superscript subscript 𝐻 0 𝑗 ∗ superscript subscript 𝛽 𝑗 ∗ 0 H_{0,j}^{\ast}:\beta_{j}^{\ast}=0 for j = 1 , … , p 𝑗 1 … 𝑝 j=1,\dots,p , as defined in ( 23 )–( 25 ), under the scale-free graph setting with p = 500 𝑝 500 p=500 . Results are shown for various values of ρ 𝜌 \rho , n 𝑛 n , SNR . Methods for comparison include LDPE, SSLasso, dScore, and the naïve score test with tuning parameter λ min subscript 𝜆 min \lambda_{\text{min}} , λ 1SE subscript 𝜆 1SE \lambda_{\text{1SE}} and λ sup subscript 𝜆 sup \lambda_{\text{sup}} . The error variance is estimated either through the scaled lasso (Sun and Zhang,, 2012 ) (S-Z) or Theorem  2.4 (T2.2). Note that as mentioned in Section  3 , we do not combine λ sup subscript 𝜆 sup \lambda_{\text{sup}} with T2.2

	$\rho$	0.2
	$n$	100			200			400
	SNR	$0.1$	$0.3$	$0.5$	$0.1$	$0.3$	$0.5$	$0.1$	$0.3$	$0.5$
${\text{Power}}_{\text{strong}}$	LDPE ${\lambda }_{\text{min}}$ S-Z	0.400	0.773	0.910	0.627	0.973	1.000	0.923	1.000	1.000
	SSLasso ${\lambda }_{\text{min}}$ S-Z	0.410	0.770	0.950	0.650	0.970	1.000	0.910	1.000	1.000
	dScore ${\lambda }_{\text{min}}$ S-Z	0.330	0.643	0.857	0.547	0.957	1.000	0.887	1.000	1.000
	nScore ${\lambda }_{\text{sup}}$ S-Z	0.403	0.847	0.960	0.727	0.990	0.997	0.940	1.000	1.000
	nScore ${\lambda }_{\text{1SE}}$ S-Z	0.427	0.763	0.893	0.677	0.977	1.000	0.957	1.000	1.000
	nScore ${\lambda }_{\text{1SE}}$ T2.2	0.357	0.763	0.890	0.680	0.977	1.000	0.910	1.000	1.000
${\text{Power}}_{\text{weak}}$	LDPE ${\lambda }_{\text{min}}$ S-Z	0.064	0.083	0.056	0.054	0.059	0.079	0.070	0.079	0.113
	SSLasso ${\lambda }_{\text{min}}$ S-Z	0.081	0.087	0.060	0.066	0.061	0.086	0.069	0.086	0.113
	dScore ${\lambda }_{\text{min}}$ S-Z	0.044	0.056	0.036	0.039	0.039	0.060	0.046	0.056	0.093
	nScore ${\lambda }_{\text{sup}}$ S-Z	0.061	0.091	0.109	0.070	0.109	0.107	0.097	0.103	0.114
	nScore ${\lambda }_{\text{1SE}}$ S-Z	0.080	0.077	0.059	0.060	0.061	0.061	0.083	0.076	0.101
	nScore ${\lambda }_{\text{1SE}}$ T2.2	0.067	0.070	0.070	0.054	0.071	0.069	0.079	0.094	0.123
T1 Error	LDPE ${\lambda }_{\text{min}}$ S-Z	0.051	0.052	0.051	0.049	0.051	0.047	0.050	0.051	0.049
	SSLasso ${\lambda }_{\text{min}}$ S-Z	0.056	0.056	0.056	0.054	0.055	0.053	0.053	0.054	0.054
	dScore ${\lambda }_{\text{min}}$ S-Z	0.035	0.040	0.040	0.033	0.036	0.034	0.035	0.037	0.034
	nScore ${\lambda }_{\text{sup}}$ S-Z	0.069	0.082	0.095	0.064	0.083	0.079	0.065	0.068	0.050
	nScore ${\lambda }_{\text{1SE}}$ S-Z	0.061	0.057	0.048	0.056	0.055	0.040	0.060	0.046	0.046
	nScore ${\lambda }_{\text{1SE}}$ T2.2	0.049	0.052	0.049	0.054	0.053	0.050	0.056	0.049	0.050
	$\rho$	0.6
	$n$	100			200			400
	SNR	$0.1$	$0.3$	$0.5$	$0.1$	$0.3$	$0.5$	$0.1$	$0.3$	$0.5$
${\text{Power}}_{\text{strong}}$	LDPE ${\lambda }_{\text{min}}$ S-Z	0.330	0.783	0.947	0.627	0.980	1.000	0.887	1.000	1.000
	SSLasso ${\lambda }_{\text{min}}$ S-Z	0.347	0.790	0.957	0.623	0.987	1.000	0.867	1.000	1.000
	dScore ${\lambda }_{\text{min}}$ S-Z	0.270	0.673	0.883	0.533	0.960	0.993	0.863	1.000	1.000
	nScore ${\lambda }_{\text{sup}}$ S-Z	0.430	0.790	0.933	0.707	0.977	1.000	0.923	1.000	1.000
	nScore ${\lambda }_{\text{1SE}}$ S-Z	0.357	0.767	0.887	0.677	0.980	0.997	0.937	1.000	1.000
	nScore ${\lambda }_{\text{1SE}}$ T2.2	0.340	0.697	0.907	0.637	0.973	1.000	0.950	1.000	1.000
${\text{Power}}_{\text{weak}}$	LDPE ${\lambda }_{\text{min}}$ S-Z	0.031	0.046	0.063	0.064	0.074	0.076	0.054	0.077	0.119
	SSLasso ${\lambda }_{\text{min}}$ S-Z	0.047	0.063	0.076	0.063	0.090	0.099	0.053	0.076	0.121
	dScore ${\lambda }_{\text{min}}$ S-Z	0.021	0.037	0.047	0.036	0.060	0.044	0.034	0.050	0.083
	nScore ${\lambda }_{\text{sup}}$ S-Z	0.071	0.089	0.136	0.081	0.121	0.104	0.114	0.113	0.123
	nScore ${\lambda }_{\text{1SE}}$ S-Z	0.039	0.060	0.050	0.076	0.074	0.066	0.070	0.067	0.104
	nScore ${\lambda }_{\text{1SE}}$T2.2	0.056	0.070	0.073	0.093	0.087	0.080	0.107	0.096	0.111
T1 Error	LDPE ${\lambda }_{\text{min}}$ S-Z	0.050	0.051	0.051	0.050	0.049	0.051	0.051	0.050	0.047
	SSLasso ${\lambda }_{\text{min}}$ S-Z	0.056	0.056	0.056	0.054	0.055	0.053	0.053	0.054	0.054
	dScore ${\lambda }_{\text{min}}$ S-Z	0.033	0.036	0.034	0.031	0.031	0.035	0.036	0.035	0.033
	nScore ${\lambda }_{\text{sup}}$ S-Z	0.065	0.080	0.093	0.064	0.084	0.088	0.070	0.071	0.054
	nScore ${\lambda }_{\text{1SE}}$ S-Z	0.056	0.060	0.045	0.061	0.051	0.040	0.058	0.048	0.047
	nScore ${\lambda }_{\text{1SE}}$ T2.2	0.049	0.053	0.051	0.060	0.061	0.049	0.066	0.049	0.048

Table 6. Table 6: Average power and type-I error rate for the hypotheses H 0 , j ∗ : β j ∗ = 0 : superscript subscript 𝐻 0 𝑗 ∗ superscript subscript 𝛽 𝑗 ∗ 0 H_{0,j}^{\ast}:\beta_{j}^{\ast}=0 for j = 1 , … , p 𝑗 1 … 𝑝 j=1,\dots,p , as defined in ( 23 )–( 25 ), under the stochastic block model setting with p = 500 𝑝 500 p=500 . Details are as in Table  5 .

	$\rho$	0.2
	$n$	100			200			400
	SNR	$0.1$	$0.3$	$0.5$	$0.1$	$0.3$	$0.5$	$0.1$	$0.3$	$0.5$
${\text{Power}}_{\text{strong}}$	LDPE ${\lambda }_{\text{min}}$ S-Z	0.370	0.793	0.937	0.687	0.990	1.000	0.914	1.000	1.000
	SSLasso ${\lambda }_{\text{min}}$ S-Z	0.393	0.803	0.933	0.687	0.990	1.000	0.892	1.000	1.000
	dScore ${\lambda }_{\text{min}}$ S-Z	0.333	0.783	0.917	0.693	0.993	1.000	0.905	1.000	1.000
	nScore ${\lambda }_{\text{sup}}$ S-Z	0.473	0.857	0.953	0.697	0.993	0.993	0.943	1.000	1.000
	nScore ${\lambda }_{\text{1SE}}$ S-Z	0.400	0.797	0.903	0.713	0.997	1.000	0.910	1.000	1.000
	nScore ${\lambda }_{\text{1SE}}$ T2.2	0.437	0.787	0.923	0.720	0.987	1.000	0.940	1.000	1.000
${\text{Power}}_{\text{weak}}$	LDPE ${\lambda }_{\text{min}}$ S-Z	0.041	0.044	0.051	0.057	0.050	0.071	0.050	0.093	0.071
	SSLasso ${\lambda }_{\text{min}}$ S-Z	0.054	0.056	0.074	0.071	0.056	0.089	0.071	0.101	0.101
	dScore ${\lambda }_{\text{min}}$ S-Z	0.037	0.044	0.057	0.060	0.046	0.077	0.056	0.101	0.094
	nScore ${\lambda }_{\text{sup}}$ S-Z	0.059	0.071	0.107	0.043	0.083	0.070	0.069	0.094	0.106
	nScore ${\lambda }_{\text{1SE}}$ S-Z	0.047	0.059	0.060	0.059	0.047	0.059	0.062	0.106	0.105
	nScore ${\lambda }_{\text{1SE}}$ T2.2	0.057	0.064	0.067	0.047	0.069	0.086	0.043	0.079	0.113
T1ER	LDPE ${\lambda }_{\text{min}}$ S-Z	0.051	0.049	0.048	0.050	0.050	0.050	0.051	0.050	0.049
	SSLasso ${\lambda }_{\text{min}}$ S-Z	0.057	0.056	0.058	0.054	0.054	0.054	0.054	0.053	0.054
	dScore ${\lambda }_{\text{min}}$ S-Z	0.043	0.040	0.041	0.041	0.044	0.042	0.042	0.042	0.041
	nScore ${\lambda }_{\text{sup}}$ S-Z	0.064	0.074	0.090	0.059	0.076	0.076	0.060	0.060	0.049
	nScore ${\lambda }_{\text{1SE}}$ S-Z	0.062	0.058	0.048	0.056	0.052	0.040	0.054	0.047	0.046
	nScore ${\lambda }_{\text{1SE}}$ T2.2	0.052	0.050	0.050	0.050	0.050	0.049	0.049	0.047	0.047
	$\rho$	0.6
	$n$	100			200			400
	SNR	$0.1$	$0.3$	$0.5$	$0.1$	$0.3$	$0.5$	$0.1$	$0.3$	$0.5$
${\text{Power}}_{\text{strong}}$	LDPE ${\lambda }_{\text{min}}$ S-Z	0.327	0.827	0.960	0.700	0.983	0.997	0.968	1.000	1.000
	SSLasso ${\lambda }_{\text{min}}$ S-Z	0.350	0.853	0.957	0.687	0.990	0.997	0.945	0.996	1.000
	dScore ${\lambda }_{\text{min}}$ S-Z	0.297	0.787	0.937	0.697	0.987	0.993	0.968	0.996	1.000
	nScore ${\lambda }_{\text{sup}}$ S-Z	0.420	0.870	0.957	0.720	0.987	1.000	0.947	1.000	1.000
	nScore ${\lambda }_{\text{1SE}}$ S-Z	0.350	0.800	0.927	0.717	0.980	1.000	0.968	1.000	1.000
	nScore ${\lambda }_{\text{1SE}}$ T2.2	0.373	0.797	0.890	0.653	0.980	1.000	0.917	1.000	1.000
${\text{Power}}_{\text{weak}}$	LDPE ${\lambda }_{\text{min}}$ S-Z	0.043	0.049	0.046	0.041	0.077	0.063	0.053	0.066	0.083
	SSLasso ${\lambda }_{\text{min}}$ S-Z	0.061	0.054	0.070	0.053	0.086	0.083	0.067	0.099	0.105
	dScore ${\lambda }_{\text{min}}$ S-Z	0.044	0.047	0.046	0.040	0.081	0.069	0.063	0.077	0.098
	nScore ${\lambda }_{\text{sup}}$ S-Z	0.054	0.073	0.093	0.063	0.093	0.094	0.061	0.094	0.096
	nScore ${\lambda }_{\text{1SE}}$ S-Z	0.059	0.056	0.044	0.054	0.087	0.074	0.067	0.086	0.103
	nScore ${\lambda }_{\text{1SE}}$ T2.2	0.046	0.054	0.060	0.059	0.064	0.084	0.051	0.087	0.116
T1 Error	LDPE ${\lambda }_{\text{min}}$ S-Z	0.049	0.049	0.049	0.049	0.050	0.049	0.049	0.047	0.048
	SSLasso ${\lambda }_{\text{min}}$ S-Z	0.057	0.056	0.056	0.053	0.054	0.054	0.053	0.053	0.053
	dScore ${\lambda }_{\text{min}}$ S-Z	0.033	0.039	0.036	0.031	0.033	0.033	0.032	0.030	0.031
	nScore ${\lambda }_{\text{sup}}$ S-Z	0.063	0.079	0.089	0.062	0.077	0.075	0.060	0.062	0.048
	nScore ${\lambda }_{\text{1SE}}$ S-Z	0.057	0.051	0.047	0.056	0.049	0.039	0.055	0.045	0.046
	nScore ${\lambda }_{\text{1SE}}$ T2.2	0.051	0.051	0.051	0.049	0.048	0.046	0.047	0.045	0.045

Equations

$$\mathbf{y}$$

$${\hat{\mathbf{\beta }}}_{\lambda }=\underset{\mathbf{b}\in {\mathbb{R}}^p}{\mathrm{argmin}}\left\{\frac{1}{2n}\Vert \mathbf{y}-\mathbf{X}\mathbf{b}{\Vert }_2^2+\lambda \Vert \mathbf{b}{\Vert }_1\right\},$$

$$\mathbf{y}$$

$$\mathrm{Pr}\left[{\beta }_j^{(\mathcal{M})}\in {\mathrm{CI}}_j^{(\mathcal{M})}\right]\ge 1-\alpha .$$

$${\mathbf{\beta }}^{(\mathcal{M})}\equiv {\left({\mathbf{X}}_{\mathcal{M}}^{\top }{\mathbf{X}}_{\mathcal{M}}\right)}^{-1}{\mathbf{X}}_{\mathcal{M}}^{\top }\mathbf{X}{\mathbf{\beta }}^{\ast }.$$

$$\frac{1}{n}{\mathbf{X}}^{\top }\mathbf{X}=\left[\begin{array}{cccc}1& 0& 0.6& 0\\ 0& 1& 0.6& 0\\ 0.6& 0.6& 1& 0\\ 0& 0& 0& 1\end{array}\right],$$

$$\mathcal{M}={\hat{\mathcal{A}}}_{\lambda }\equiv \mathrm{supp}\left({\hat{\mathbf{\beta }}}_{\lambda }\right)\equiv \left\{j:{\hat{\beta }}_{\lambda ,j}\ne 0\right\},$$

$${\mathbf{\beta }}^{(\mathcal{M})}$$

$$={\left[\begin{array}{cc}1& 0.6\\ 0.6& 1\end{array}\right]}^{-1}\left[\begin{array}{cccc}0& 1& 0.6& 0\\ 0.6& 0.6& 1& 0\end{array}\right]{\left[\begin{array}{cccc}1& 1& 0& 0\end{array}\right]}^{\top }$$

$$=\left[\begin{array}{c}0.4375\\ 0.9375\end{array}\right]\ne \left[\begin{array}{c}1\\ 0\end{array}\right]={\mathbf{\beta }}_{\mathcal{M}}^{\ast }.$$

$${\mathbf{\beta }}_{\lambda }\equiv \underset{\mathbf{b}\in {\mathbb{R}}^p}{\mathrm{argmin}}\left\{\frac{1}{2n}\mathbb{E}\left[{\Vert \mathbf{y}-\mathbf{X}\mathbf{b}\Vert }_2^2\right]+\lambda {\Vert \mathbf{b}\Vert }_1\right\}.$$

$${\tilde{\mathbf{\beta }}}^{({\hat{\mathcal{A}}}_{\lambda })}\equiv {\left({\mathbf{X}}_{{\hat{\mathcal{A}}}_{\lambda }}^{\top }{\mathbf{X}}_{{\hat{\mathcal{A}}}_{\lambda }}\right)}^{-1}{\mathbf{X}}_{{\hat{\mathcal{A}}}_{\lambda }}^{\top }\mathbf{y},$$

$$\mathrm{Pr}\left({\beta }_j^{({\hat{\mathcal{A}}}_{\lambda })}\in {\mathrm{CI}}_j^{({\hat{\mathcal{A}}}_{\lambda })}\right)\ge 1-\alpha .$$

$$\sqrt{\frac{\log (p)}{n}}\frac{1}{\lambda }\to 0.$$

$${\varphi }^{\ast 2}\equiv \underset{n\to \infty }{\mathrm{liminf}}\frac{{\mathbf{a}}^{\top }\hat{\mathbf{\Sigma }}\mathbf{a}}{\Vert {\mathbf{a}}_{\mathcal{B}}{\Vert }_2^2}>0,$$

$${\varphi }^2(q)\equiv \underset{\Vert {\mathbf{b}}_{{\mathcal{A}}^{\ast c}}{\Vert }_0\le q}{\sup }\frac{{\mathbf{b}}^{\top }\hat{\mathbf{\Sigma }}\mathbf{b}}{\Vert \mathbf{b}{\Vert }_2^2},$$

$$\underset{n\to \infty }{\lim }\frac{b_{\lambda \min }}{\lambda }=\xi >0,$$

$${\Vert {\mathbf{X}}_{{\mathcal{A}}_{\lambda }^c}{\mathbf{\beta }}_{{\mathcal{A}}_{\lambda }^c}^{\ast }\Vert }_2=\mathcal{O}\left(\sqrt{\log (p)}\right).$$

$$\underset{n\to \infty }{\lim }{\Vert {\mathbf{X}}_{{\mathcal{A}}_{\lambda }^c}^{\top }{\mathbf{X}}_{{\mathcal{A}}_{\lambda }}{\left({\mathbf{X}}_{{\mathcal{A}}_{\lambda }}^{\top }{\mathbf{X}}_{{\mathcal{A}}_{\lambda }}\right)}^{-1}\mathrm{sign}\left({\mathbf{\beta }}_{\lambda ,{\mathcal{A}}_{\lambda }}\right)\Vert }_{\infty }<1-\delta .$$

$$\underset{n\to \infty }{\lim }\frac{2\sqrt{2q^{\ast }}}{n{\varphi }^{\ast 2}}{\Vert {\mathbf{X}}_{{\mathcal{A}}_{\lambda }^c}^{\top }{\mathbf{X}}_{{\mathcal{A}}_{\lambda }}\Vert }_{\infty }<1-\delta .$$

$${\mathbf{r}}^w\equiv {\mathbf{e}}^j({\mathbf{X}}_{{\mathcal{A}}_{\lambda }}^{\top }{\mathbf{X}}_{{\mathcal{A}}_{\lambda }}{)}^{-1}{\mathbf{X}}_{{\mathcal{A}}_{\lambda }}^{\top },$$

$$\frac{{\tilde{\beta }}_j^{({\hat{\mathcal{A}}}_{\lambda })}-{\beta }_j^{({\hat{\mathcal{A}}}_{\lambda })}}{{\sigma }_{\mathbf{ϵ}}\sqrt{{\left[({\mathbf{X}}_{{\hat{\mathcal{A}}}_{\lambda }}^{\top }{\mathbf{X}}_{{\hat{\mathcal{A}}}_{\lambda }}{)}^{-1}\right]}_{(j,j)}}}{\to }_d\mathcal{N}\left(0,1\right),$$

$$\frac{1}{n-{\hat{q}}_{\lambda }}{\Vert \mathbf{y}-{\mathbf{X}}_{{\hat{\mathcal{A}}}_{\lambda }}{\tilde{\mathbf{\beta }}}^{({\hat{\mathcal{A}}}_{\lambda })}\Vert }_2^2{\to }_p{\sigma }_{\mathbf{ϵ}}^2,$$

$${\mathrm{CI}}_j^{({\hat{\mathcal{A}}}_{\lambda })}\equiv \left({\tilde{\beta }}_j^{({\hat{\mathcal{A}}}_{\lambda })}-1.96\times {\hat{\sigma }}_{\mathbf{ϵ}}\sqrt{{\left[{\mathbf{X}}_{{\hat{\mathcal{A}}}_{\lambda }}^{\top }{\mathbf{X}}_{{\hat{\mathcal{A}}}_{\lambda }}\right]}_{(j,j)}},{\tilde{\beta }}_j^{({\hat{\mathcal{A}}}_{\lambda })}+1.96\times {\hat{\sigma }}_{\mathbf{ϵ}}\sqrt{{\left[{\mathbf{X}}_{{\hat{\mathcal{A}}}_{\lambda }}^{\top }{\mathbf{X}}_{{\hat{\mathcal{A}}}_{\lambda }}\right]}_{(j,j)}}\right).$$

$$A_{(j,k)}=\{\begin{array}{ll}1& \text{mbox}forj=k\\ \rho & \text{mbox}for(j,k)\in \mathcal{E}\\ 0& \text{mbox}otherwise\end{array},$$

$${\beta }_j^{\ast }=\{\begin{array}{ll}1& \text{mbox}forj=1\\ 0.1& \text{mbox}for2\le j\le 5\\ 0& \text{mbox}otherwise\end{array}.$$

$$\text{Coverage Proportion}\equiv\sum_{b=1}^{B}\sum_{j\in\hat{\mathcal{A}}_{\lambda}^{b}}1\left\{\beta_{j}^{(\hat{\mathcal{A}}_{\lambda}^{b})}\in\mathrm{CI}^{(\hat{\mathcal{A}}_{\lambda}^{b}),b}_{j}\right\}\bigg{/}\left|\hat{\mathcal{A}}_{\lambda}^{b}\right|,$$

$$S^j\equiv {\mathbf{x}}_j^{\top }\left(\mathbf{y}-{\tilde{\mathbf{y}}}^{({\hat{\mathcal{A}}}_{\lambda }\backslash \{j\})}\right)\equiv {\mathbf{x}}_j^{\top }\left({\mathbf{I}}_n-{\mathbf{P}}^{-j}\right)\mathbf{y},$$

$${\tilde{\mathbf{y}}}^{({\hat{\mathcal{A}}}_{\lambda }\backslash \{j\})}\equiv {\mathbf{X}}_{{\hat{\mathcal{A}}}_{\lambda }\backslash \{j\}}{\tilde{\mathbf{\beta }}}^{({\hat{\mathcal{A}}}_{\lambda }\backslash \{j\})},$$

$${\mathbf{P}}^{-j}\equiv {\mathbf{X}}_{{\hat{\mathcal{A}}}_{\lambda }\backslash \{j\}}{\left({\mathbf{X}}_{{\hat{\mathcal{A}}}_{\lambda }\backslash \{j\}}^{\top }{\mathbf{X}}_{{\hat{\mathcal{A}}}_{\lambda }\backslash \{j\}}\right)}^{-1}{\mathbf{X}}_{{\hat{\mathcal{A}}}_{\lambda }\backslash \{j\}}^{\top }$$

$$\underset{n\to \infty }{\lim }\frac{b_{\lambda \min }}{\lambda }=\xi >0,$$

Full text

Keywords— Confidence interval; Lasso; $p$-value; Post-selection inference; Significance testing.

1 Introduction

In this paper, we consider the linear model

[TABLE]

where ${\bm{X}}=[{\bm{x}}_{1},\dots,{\bm{x}}_{p}]$ is an $n\times p$ deterministic design matrix, ${\bm{\epsilon}}$ is a vector of independent and identically distributed errors with $\mathbb{E}\left[\epsilon_{i}\right]=0$ and $\mathrm{Var}\left[\epsilon_{i}\right]=\sigma^{2}_{\bm{\epsilon}}$, and ${\bm{\beta}}^{\ast}$ is a $p$-vector of coefficients. Without loss of generality, we assume that the columns of ${\bm{X}}$ are centered and standardized, such that $\sum_{i=1}^{n}X_{(i,k)}=0$ and $\|{\bm{x}}_{k}\|_{2}^{2}=n$ for $k=1,\ldots,p$.

When the number of variables $p$ is much smaller than the sample size $n$, estimation and inference for the vector ${\bm{\beta}}^{\ast}$ are straightforward. For instance, estimation can be performed using ordinary least squares, and inference can be conducted using classical approaches (see, e.g., Gelman and Hill,, 2006; Weisberg,, 2013).

As the scope and scale of data collection have increased across virtually all fields, there is an increase in data sets that are high dimensional, in the sense that the number of variables, $p$, is larger than the number of observations, $n$. In this setting, classical approaches for estimation and inference of ${\bm{\beta}}^{\ast}$ cannot be directly applied. In the past 20 years, a vast statistical literature has focused on estimating ${\bm{\beta}}^{\ast}$ in high dimensions. In particular, penalized regression methods, such as the lasso (Tibshirani,, 1996),

[TABLE]

can be used to estimate ${\bm{\beta}}^{\ast}$. However, the topic of inference in the high-dimensional setting remains relatively less explored, despite promising recent work in this area. Roughly speaking, recent work on inference in the high-dimensional setting falls into two classes: (i) methods that examine the null hypothesis $H_{0,j}^{\ast}:\beta^{\ast}_{j}=0$; and (ii) methods that make inference based on a sub-model. We will review these two classes of methods in turn.

First, we review methods that examine the null hypothesis $H_{0,j}^{\ast}:\beta^{\ast}_{j}=0$, i.e. that the variable ${\bm{x}}_{j}$ is unassociated with the outcome ${\bm{y}}$, conditional on all other variables. It might be tempting to estimate ${\bm{\beta}}^{\ast}$ using the lasso (2), and then (for instance) to construct a confidence interval around $\hat{\beta}_{\lambda,j}$. Unfortunately, such an approach is problematic, because $\hat{\bm{\beta}}_{\lambda}$ is a biased estimate of ${\bm{\beta}}^{\ast}$. To remedy this problem, we can apply a one-step adjustment to $\hat{\bm{\beta}}_{\lambda}$, such that under appropriate assumptions, the resulting debiased estimator is asymptotically unbiased for ${\bm{\beta}}^{\ast}$. This is similar to the idea of two-step estimation in nonparametric and semiparametric inference (see, e.g., Hahn,, 1998; Hirano et al.,, 2003). With the one-step adjustment, $p$-values and confidence intervals can be constructed around this debiased estimator. Such an approach is taken by the low dimensional projection estimator (LDPE; Zhang and Zhang,, 2014; van de Geer et al.,, 2014), the debiased lasso test with unknown population covariance (SSLasso; Javanmard and Montanari,, 2013; Javanmard and Montanari, 2014a, ), the debiased lasso test with known population covariance (SDL; Javanmard and Montanari, 2014b, ), and the decorrelated score test (dScore; Ning and Liu,, 2017). See Dezeure et al., (2015) for a review of such procedures. In what follows, we will refer to these and related approaches for testing $H_{0,j}^{\ast}:\beta^{\ast}_{j}=0$ as debiased lasso tests.

Next, we review recent work that makes statistical inference based on a sub-model. Recall that the challenge in high dimensions stems from the fact that when $p>n$, classical statistical methods cannot be applied; for instance, we cannot even perform ordinary least squares (OLS). This suggests a simple approach: given an index set $\mathcal{M}\subseteq\{1,\dots,p\}$, let ${\bm{X}}_{\mathcal{M}}$ denote the columns of ${\bm{X}}$ indexed by $\mathcal{M}$. Then, we can consider performing inference based on the sub-model composed only of the features in the index set $\mathcal{M}$. That is, rather than considering the model (1), we consider the sub-model

[TABLE]

In (3), the notation ${\bm{\beta}}^{(\mathcal{M})}$ and ${\bm{\epsilon}}^{(\mathcal{M})}$ emphasizes that the true regression coefficients and corresponding noise are functions of the set $\mathcal{M}$.

Now, provided that $|\mathcal{M}|<n$, we can perform estimation and inference on the vector ${\bm{\beta}}^{(\mathcal{M})}$ using classical statistical approaches. For instance, we can consider building confidence intervals $\mathrm{CI}_{j}^{(\mathcal{M})}$ such that for any $j\in\mathcal{M}$,

[TABLE]

At first blush, the problems associated with high dimensionality have been solved!

Of course, there are some problems with the aforementioned approach. The first problem is that the coefficients in the sub-model (3) typically are not the same as the coefficients in the original model (1) (Berk et al.,, 2013). Roughly speaking, the problem is that the coefficients in the model (1) quantify the linear association between a given variable and the response, conditional on the other $p-1$ variables, whereas the coefficients in the sub-model (3) quantify the linear association between a variable and the response, conditional on the other $|\mathcal{M}|-1$ variables in the sub-model. The true regression coefficients in the sub-model are of the form

[TABLE]

Thus, ${\bm{\beta}}^{(\mathcal{M})}\neq{\bm{\beta}}^{\ast}_{\mathcal{M}}$ unless ${\bm{X}}^{\top}_{\mathcal{M}}{\bm{X}}_{\mathcal{M}^{c}}{\bm{\beta}}^{\ast}_{\mathcal{M}^{c}}={\bm{0}}$. To see this more concretely, consider the following example with $p=4$ deterministic variables. Let

[TABLE]

and set ${\bm{\beta}}^{\ast}=(1,1,0,0)^{\top}$. The above design matrix does not satisfy the strong irrepresentable condition needed for selection consistency of lasso (Zhao and Yu,, 2006). Thus, if we take $\mathcal{M}$ to equal the support of the lasso estimate, i.e.,

[TABLE]

then it is easy to verify that for some $\lambda$, $\mathcal{M}=\{2,3\}$, in which case

[TABLE]

The second problem that arises in restricting our attention to the sub-model (3) is that in practice, the index set $\mathcal{M}$ is not pre-specified. Instead, it is typically chosen based on the data. The problem is that if we construct the index set $\mathcal{M}$ based on the data, and then apply classical inference approaches on the vector ${\bm{\beta}}^{(\mathcal{M})}$, the resulting $p$-values and confidence intervals will not be valid (see, e.g., Pötscher,, 1991; Kabaila,, 1998; Leeb and Pötscher,, 2003, 2005; Leeb and Pötscher, 2006a, ; Leeb and Pötscher, 2006b, ; Leeb and Pötscher,, 2008; Kabaila,, 2009; Berk et al.,, 2013). This is because we peeked at the data twice: once to determine which variables to include in $\mathcal{M}$, and then again to test hypotheses associated with those variables. Consequently, an extensive recent body of literature has focused on the task of performing inference on ${\bm{\beta}}^{(\mathcal{M})}$ in (3) given that $\mathcal{M}$ was chosen based on the data. Cox, (1975) proposed the idea of sample-splitting to break up the dependence of variable selection and hypothesis testing. Wasserman and Roeder, (2009) studied sample-splitting in application to the lasso, marginal regression and forward step-wise regression. Meinshausen et al., (2009) extended the single-splitting proposal of Wasserman and Roeder, (2009) to multi-splitting, which improved statistical power and reduced the number of falsely selected variables. Berk et al., (2013) instead considered simultaneous inference, which is universally valid under all possible model selection procedures without sample-splitting. More recently, Lee et al., (2016) and Tibshirani et al., (2016) studied the geometry of the lasso and sequential regression, respectively, and proposed exact post-selection inference methods conditional on the random set of selected variables. See Taylor and Tibshirani, (2015) for a review of post-selection inference procedures.

The procedures outlined above successfully address the second problem arising from restricting the attention to a sub-model, namely the randomness of the set $\mathcal{M}$. However, they do not address the first problem regarding the difference in the target of inference (5) unless ${\bm{X}}^{\top}_{\mathcal{M}}{\bm{X}}_{\mathcal{M}^{c}}{\bm{\beta}}^{\ast}_{\mathcal{M}^{c}}={\bm{0}}$. In the case of lasso, valid inference for ${\bm{\beta}}^{\ast}_{\mathcal{A}^{\ast}}$ can be obtained if $\mathcal{M}=\hat{\mathcal{A}}_{\lambda}=\mathcal{A}^{\ast}$. However, among other conditions, this requires the strong irrepresentable condition, which is known to be a restrictive condition that is not likely to hold in high dimensions (Zhao and Yu,, 2006).

In a recent Statistical Science paper, Leeb et al., (2015) performed a simulation study, in which they obtained a set $\mathcal{M}$ using variable selection, and then calculated “naïve” confidence intervals for ${\bm{\beta}}^{(\mathcal{M})}$ using ordinary least squares, without accounting for the fact that the set $\mathcal{M}$ was chosen based on the data. Of course, conventional wisdom dictates that the resulting confidence intervals will be much too narrow. In fact, this is what Leeb et al., (2015) found, when they used best subset selection to construct the set $\mathcal{M}$. However, surprisingly, when the lasso was used to construct the set $\mathcal{M}$, the confidence intervals induced by (4) had approximately correct coverage. This is in stark contrast to the existing literature!

In this paper, we present a theoretical justification for the empirical finding in Leeb et al., (2015). The main idea of our paper is to establish selection consistency of the lasso estimate with respect to its noiseless counterpart. This result allows us to perform valid inference for the support of the noiseless lasso without needing post-selection or sample splitting strategies. Furthermore, we use our theoretical findings to also develop the naïve score test, a simple procedure for testing the null hypothesis $H^{\ast}_{0,j}:\beta^{\ast}_{j}=0$ for all $j=1,\dots,p$.

The rest of this paper is organized as follows. In Sections 2 and 3, we focus on post-selection inference: we seek to perform inference on ${\bm{\beta}}^{(\mathcal{M})}$ in (3), where $\mathcal{M}$ is selected based on the lasso, i.e., $\mathcal{M}=\hat{\mathcal{A}}_{\lambda}$ (6). In Section 2, we point out a previously overlooked scenario in selection consistency theory: although $\hat{\mathcal{A}}_{\lambda}$ (6) is random, with high probability it is equal to the support of the noiseless lasso under relatively mild regularity conditions. This result implies that we can use classical methods for inference on ${\bm{\beta}}^{(\mathcal{M})}$, when $\mathcal{M}=\hat{\mathcal{A}}_{\lambda}$. In Section 3, we provide empirical evidence in support of these theoretical findings. In Sections 4 and Section 5, we instead focus on the task of performing inference on ${\bm{\beta}}^{\ast}$ in (1). We propose the naïve score test in Section 4, and study its empirical performance in Section 5. We end with a discussion of future research directions in Section 6. Technical proofs are relegated to the online Supplementary Materials.

We now introduce some notation that will be used throughout the paper. We use “$\equiv$” to denote equalities by definition, and “$\asymp$” for the asymptotic order. We use $1\{\cdot\}$ for the indicator function; “$\vee$” and “$\wedge$” denote the maximum and minimum of two real numbers, respectively. For any real number $a\in\mathbb{R}$, $a_{+}\equiv a\vee 0$. Given a set $\mathcal{S}$, $|\mathcal{S}|$ denotes its cardinality and $-\mathcal{S}\equiv\mathcal{S}^{c}$ denotes its complement. We use bold upper case fonts to denote matrices, bold lower case fonts for vectors, and normal fonts for scalars. We use symbols with a superscript “$\ast$”, e.g., ${\bm{\beta}}^{\ast}$ and $\mathcal{A}^{\ast}\equiv\mathrm{supp}({\bm{\beta}}^{\ast})$, to denote the true population parameters associated with the full linear model (1); we use symbols superscripted by a set in the parentheses, e.g., ${\bm{\beta}}^{(\mathcal{M})}$, to denote quantities related to the sub-model (3). Symbols subscripted by “$\lambda$” and with a hat, e.g., $\hat{\bm{\beta}}_{\lambda}$ and $\hat{\mathcal{A}}_{\lambda}$, denote parameter estimates from the lasso estimator (2) with tuning parameter $\lambda>0$; symbols subscripted by “$\lambda$” and without a hat, e.g., ${\bm{\beta}}_{\lambda}$ and $\mathcal{A}_{\lambda}\equiv\mathrm{supp}({\bm{\beta}}_{\lambda})$, are associated with the noiseless lasso estimator (van de Geer and Bühlmann,, 2009; van de Geer,, 2017),

[TABLE]

For any vector ${\bm{b}}$, matrix ${\bm{\Sigma}}$, and index sets $\mathcal{S}_{1}$ and $\mathcal{S}_{2}$, we use ${\bm{b}}_{\mathcal{S}_{1}}$ to denote the sub-vector of ${\bm{b}}$ comprised of elements of $\mathcal{S}_{1}$, and ${\bm{\Sigma}}_{(\mathcal{S}_{1},\mathcal{S}_{2})}$ to denote the sub-matrix of ${\bm{\Sigma}}$ with rows in $\mathcal{S}_{1}$ and columns in $\mathcal{S}_{2}$.

2 Theoretical Justification for Naïve Confidence Intervals

Recall that ${\bm{\beta}}^{(\hat{\mathcal{A}}_{\lambda})}$ was defined in (5). The simulation results of Leeb et al., (2015) suggest that if we perform ordinary least squares using the variables contained in the support set of the lasso, $\hat{\mathcal{A}}_{\lambda}$, then the classical confidence intervals associated with the least squares estimator,

[TABLE]

have approximately correct coverage, where correct coverage means that for all $j\in\hat{\mathcal{A}}_{\lambda}$,

[TABLE]

We reiterate that in (9), $\mathrm{CI}_{j}^{(\hat{\mathcal{A}}_{\lambda})}$ is the confidence interval output by standard least squares software applied to the data $({\bm{y}},{\bm{X}}_{\hat{\mathcal{A}}_{\lambda}})$. This goes against our statistical intuition: it seems that by fitting a lasso model and then performing least squares on the selected set, we are peeking at the data twice, and thus we would expect the confidence interval $\mathrm{CI}_{j}^{(\hat{\mathcal{A}}_{\lambda})}$ to be much too narrow.

In this section, we present a theoretical result that suggests that, in fact, this “double-peeking” might not be so bad. Our key insight is as follows: under certain assumptions, the set of variables selected by the lasso is deterministic and non-data-dependent with high probability. Thus, fitting a least squares model on the variables selected by the lasso does not really constitute peeking at the data twice: effectively, with high probability, we are only peeking at it once. That means that the naïve confidence intervals obtained from ordinary least squares will have approximately correct coverage, in the sense of (9).

We first introduce the required conditions for our theoretical result.

• (M1)

The design matrix ${\bm{X}}$ is deterministic, with columns in general position (Rosset et al.,, 2004; Zhang,, 2010; Dossal,, 2012; Tibshirani,, 2013). Columns of ${\bm{X}}$ are standardized, i.e., for any $j=1,\dots,p$, ${\bm{x}}_{j}^{\top}{\bm{x}}_{j}=n$. The error ${\bm{\epsilon}}$ in (1) has independent entries and sub-Gaussian tails. The sample size $n$, dimension $p$ and tuning parameter $\lambda$ satisfy

[TABLE]

• (E)

Let $\hat{\bm{\Sigma}}={\bm{X}}^{\top}{\bm{X}}/n$. For any index set $\mathcal{I}$ with $|\mathcal{I}|=\mathcal{O}(q^{\ast})$, let $\mathcal{B}$ be any index set such that $\mathcal{B}\supseteq\mathcal{I}$, $|\mathcal{B}\backslash\mathcal{I}|\leq|\mathcal{I}|$. Then for all ${\bm{a}}\in\mathbb{R}^{p}$ that satisfy $\|{\bm{a}}_{\mathcal{I}^{c}}\|_{1}\leq\|{\bm{a}}_{\mathcal{I}}\|_{1}$, and $\|{\bm{a}}_{\mathcal{B}^{c}}\|_{\infty}\leq\min_{j\in(\mathcal{B}\backslash\mathcal{I})}|a_{j}|$,

[TABLE]

In addition, the restricted sparse eigenvalue,

[TABLE]

satisfies $\phi^{2}(q^{\ast})=\mathcal{O}(\sqrt{\log(p)/q^{\ast}})$.

• (M2)

Recall that $\mathcal{A}^{\ast}\equiv\mathrm{supp}({\bm{\beta}}^{\ast})$. Let $\mathcal{A}_{\lambda}\equiv\mathrm{supp}({\bm{\beta}}_{\lambda})$, $b_{\lambda\min}\equiv\min_{j\in\mathcal{A}_{\lambda}}|\beta_{\lambda,j}|$ with ${\bm{\beta}}_{\lambda}$ defined in (7). The signal strength in $\mathcal{A}_{\lambda}$ satisfies

[TABLE]

and the signal strength outside of $\mathcal{A}_{\lambda}$ satisfies

[TABLE]

• (T)

The strong irrepresentable condition with respect to $\mathcal{A}_{\lambda}$ and ${\bm{\beta}}_{\lambda}$ holds, i.e., there exists $\delta>0$ such that

[TABLE]

Condition (M1) is mild and standard in literature. Note that in (M1), we require the lasso tuning parameter $\lambda$ to approach zero at a slightly slower rate than $\sqrt{\log(p)/n}$ to control the randomness of the error ${\bm{\epsilon}}$. Most standard literature requires $\lambda>C\sqrt{\log(p)/n}$ for some constant $C>0$, which our condition also satisfies.

In (M2), the requirement that $\lim_{n\to\infty}b_{\lambda,\min}/\lambda>0$ indicates that the noiseless lasso regression coefficients are either 0, or asymptotically no smaller than $\lambda$, which is larger than $\sqrt{\log(p)/n}$. Note that this assumption concerns variables chosen by the noiseless lasso and not the true model. In the second part of (M2), \big{\|}{\bm{X}}_{\mathcal{A}_{\lambda}^{c}}{\bm{\beta}}^{\ast}_{\mathcal{A}_{\lambda}^{c}}\big{\|}_{2}=\mathcal{O}(\sqrt{\log(p)}) implies that the total signal strength of weak signal variables that are not selected by the noiseless lasso cannot be too large to be distinguishable from the sub-Gaussian noise. As shown in the proof of Proposition 2.2, the condition $b_{\lambda,\min}=\mathcal{O}(\lambda)$ is important in showing that $\lim_{n\to\infty}\Pr\left[\mathcal{A}_{\lambda}\subseteq\hat{\mathcal{A}}_{\lambda}\right]=1$, whereas \big{\|}{\bm{X}}_{\mathcal{A}_{\lambda}^{c}}{\bm{\beta}}^{\ast}_{\mathcal{A}_{\lambda}^{c}}\big{\|}_{2}=\mathcal{O}\left(\sqrt{\log(p)}\right) is instrumental in showing $\lim_{n\to\infty}\Pr\left[\mathcal{A}_{\lambda}\supseteq\hat{\mathcal{A}}_{\lambda}\right]=1$.

Condition (E) manifests the behavior of the eigenvalues of $\hat{\bm{\Sigma}}$. Specifically, the first part of (E) is the restricted eigenvalue condition (Bickel et al.,, 2009; van de Geer and Bühlmann,, 2009), except that instead of requiring $|\mathcal{I}|=q^{\ast}$ as in Bickel et al., (2009) and van de Geer and Bühlmann, (2009), we here require $|\mathcal{I}|=\mathcal{O}(q^{\ast})$. The second part of (E) is the sparse Reisz, or sparse eigenvalue condition (Zhang and Huang,, 2008; Belloni and Chernozhukov,, 2013). Both the restricted eigenvalue and sparse Reisz conditions are standard and mild conditions in the literature. (E) implies that $\log(p)/q^{\ast}\to\psi>1$.

Condition (T) is the strong irrepresentable condition with respect to $\mathcal{A}_{\lambda}$ and ${\bm{\beta}}_{\lambda}$. Condition (T) is likely weaker than the classical irrepresentable condition with respect to ${\bm{\beta}}^{\ast}$ and $\mathcal{A}^{\ast}$, proposed in Zhao and Yu, (2006), because the classical irrepresentable condition implies that $\mathcal{A}_{\lambda}=\mathcal{A}^{\ast}$ with large $n$, in which case (T) holds as it becomes identical to the classical irrepresentable condition. Lemma 2.1 shows a sufficient condition for (T), which is proven in Section B in the Supplementary Materials.

Lemma 2.1**.**

Suppose conditions (M1), (M2) and (E) hold. Condition (T) holds if

[TABLE]

Condition (13) allows $\|{\bm{X}}^{\top}_{\mathcal{A}_{\lambda}^{c}}{\bm{X}}_{\mathcal{A}_{\lambda}}\|_{\infty}$ to diverge to infinity at a slower rate than $n/\sqrt{q^{\ast}}$. A more stringent version of condition (13) is presented in Voorman et al., (2014).

We use (T) to prove $\lim_{n\to\infty}\|{\bm{\tau}}_{\lambda,\mathcal{A}_{\lambda}^{c}}\|_{\infty}\leq 1-\delta$, where $\lambda n{\bm{\tau}}_{\lambda}={\bm{X}}^{\top}\left({\bm{X}}{\bm{\beta}}^{\ast}-{\bm{X}}{\bm{\beta}}_{\lambda}\right)$ is the stationary condition of (7). Note that $\lim_{n\to\infty}\|{\bm{\tau}}_{\lambda,\mathcal{A}_{\lambda}^{c}}\|_{\infty}\leq 1-\delta$ is the required condition for Proposition 2.2, and (T) might be sufficient but unncessary. However, (T) is more standard and understandable than the condition $\lim_{n\to\infty}\|{\bm{\tau}}_{\lambda,\mathcal{A}_{\lambda}^{c}}\|_{\infty}\leq 1-\delta$.

We now present Proposition 2.2, which is proven in Section A of the online Supplementary Materials.

Proposition 2.2**.**

Suppose conditions (M1), (M2), (E) and (T) hold. Then, we have $\lim_{n\to\infty}\Pr[\hat{\mathcal{A}}_{\lambda}=\mathcal{A}_{\lambda}]=1$, where $\mathcal{A}_{\lambda}\equiv\mathrm{supp}({\bm{\beta}}_{\lambda})$, with ${\bm{\beta}}_{\lambda}$ defined in (7).

The proof of Proposition 2.2 is in the same flavor as Meinshausen and Bühlmann, (2006); Zhao and Yu, (2006); Tropp, (2006) and Wainwright, (2009), and is based on absorbing the contribution of the weak signals ${\bm{X}}_{\mathcal{A}_{\lambda}^{c}}{\bm{\beta}}_{\mathcal{A}_{\lambda}^{c}}^{\ast}$ into the noise vector. However, variable selection consistency asserts that $\Pr[\hat{\mathcal{A}}_{\lambda}=\mathcal{A}^{\ast}]\to 1$, whereas Proposition 2.2 states that $\Pr[\hat{\mathcal{A}}_{\lambda}=\mathcal{A}_{\lambda}]\to 1$. Consequently, variable selection consistency requires the irrepresentable and signal strength conditions with respect to ${\bm{\beta}}^{\ast}$, whereas Proposition 2.2 requires the irrepresentable and signal strength conditions with respect to ${\bm{\beta}}_{\lambda}$. These conditions are likely much milder than those for the variable selection consistency of the lasso, as confirmed in our simulations in Section 3. In these simulations, we estimate $\Pr[\hat{\mathcal{A}}_{\lambda}=\mathcal{A}_{\lambda}]$ in 36 settings with two choices of $\lambda$. As shown in Table 4, in most settings, $\hat{\mathcal{A}}_{\lambda}=\mathcal{A}_{\lambda}$ with high probability, especially in cases with large $n$. We also estimated $\Pr[\hat{\mathcal{A}}_{\lambda}=\mathcal{A}^{\ast}]$ in the same 36 settings, with the same two choices of $\lambda$, as well as an additional choice of $\lambda$ that minimizes the cross-validated MSE; we found that $\Pr[\hat{\mathcal{A}}_{\lambda}=\mathcal{A}^{\ast}]=0$ in all $36\times 3$ settings.

Based on Proposition 2.2, we could build asymptotically valid confidence intervals, as shown in Theorem 2.3.

Proposition 2.2 suggests that asymptotically, we “pay no price” for peeking at our data by performing the lasso: we should be able to perform downstream analyses on the subset of variables in $\hat{\mathcal{A}}_{\lambda}$ as though we had obtained that subset without looking at the data. This intuition will be formalized in Theorem 2.3.

Theorem 2.3, which is proven in Section C in the online Supplementary Materials, shows that $\tilde{\bm{\beta}}^{(\hat{\mathcal{A}}_{\lambda})}$ in (8) is asymptotically normal, with mean and variance suggested by classical least squares theory: that is, the fact that $\hat{\mathcal{A}}_{\lambda}$ was selected by peeking at the data has no effect on the asymptotic distribution of $\tilde{\bm{\beta}}^{(\hat{\mathcal{A}}_{\lambda})}$. This result requires that $\lambda$ be chosen in a non-data-adaptive way. Otherwise, $\mathcal{A}_{\lambda}$ will be affected by the random error ${\bm{\epsilon}}$ through $\lambda$, which complicates the distribution of $\tilde{\bm{\beta}}^{(\hat{\mathcal{A}}_{\lambda})}$. Theorem 2.3 also requires Condition (W), which is used to apply the Lindeberg-Feller Central Limit Theorem. This condition can be relaxed if the noise ${\bm{\epsilon}}$ is normally distributed.

• (W)

$\lambda$, ${\bm{\beta}}^{\ast}$ and ${\bm{X}}$ are such that $\lim_{n\to\infty}\|{\bf r}^{w}\|_{\infty}/\|{\bf r}^{w}\|_{2}\rightarrow 0$, where

[TABLE]

and ${\bm{e}}^{j}$ is the row vector of length $|\mathcal{A}_{\lambda}|$ with the entry corresponding to $\beta^{\ast}_{j}$ equal to one, and zero otherwise.

Theorem 2.3**.**

Suppose $\lim_{n\to\infty}\Pr\left[\hat{\mathcal{A}}_{\lambda}=\mathcal{A}_{\lambda}\right]=1$ and (W) holds. Then, for any $j\in\hat{\mathcal{A}}_{\lambda}$,

[TABLE]

where $\tilde{\bm{\beta}}^{(\hat{\mathcal{A}}_{\lambda})}$ is defined in (8) and ${\bm{\beta}}^{(\hat{\mathcal{A}}_{\lambda})}$ in (5), and ${\sigma_{\bm{\epsilon}}}$ is the variance of ${\bm{\epsilon}}$ in (1).

The error standard deviation $\sigma_{\bm{\epsilon}}$ in (14) is usually unknown. It can be estimated using various high-dimensional estimation methods, e.g., the scaled lasso (Sun and Zhang,, 2012), cross-validation (CV) based methods (Fan et al.,, 2012) or method-of-moments based methods (Dicker,, 2014); see a comparison study of high dimensional error variance estimation methods in Reid et al., (2016). Alternatively, Theorem 2.4 shows that we could also consistently estimate the error variance using the post-selection OLS residual sum of square (RSS).

Theorem 2.4**.**

Suppose $\lim_{n\to\infty}\Pr\left[\hat{\mathcal{A}}_{\lambda}=\mathcal{A}_{\lambda}\right]=1$ and $\log(p)/(n-q_{\lambda})\to 0$, where $q_{\lambda}\equiv|\mathcal{A}_{\lambda}|$. Then

[TABLE]

where $\hat{q}_{\lambda}\equiv|\hat{\mathcal{A}}_{\lambda}|$.

Theorem 2.4 is proven in Section D in the online Supplementary Materials. In (15), ${\bm{y}}-{\bm{X}}_{\hat{\mathcal{A}}_{\lambda}}\tilde{\bm{\beta}}^{\hat{\mathcal{A}}_{\lambda}}$ is the fitted OLS residual on the sub-model (3). Also, $\log(p)/(n-q_{\lambda})\to 0$ is a weak condition: since $\log(p)/n\to 0$, $\log(p)/(n-q_{\lambda})\to 0$ is satisfied if $\lim_{n\to\infty}q_{\lambda}/n<1$.

To summarize, in this section, we have provided a theoretical justification for a procedure that seems, intuitively, to be statistically unjustifiable:

Perform the lasso in order to obtain the support set $\hat{\mathcal{A}}_{\lambda}$; 2. 2.

Use least squares to fit the sub-model containing just the features in $\hat{\mathcal{A}}_{\lambda}$; 3. 3.

Use the classical confidence intervals from that least squares model, without accounting for the fact that $\hat{\mathcal{A}}_{\lambda}$ was obtained by peeking at the data.

Theorem 2.3 guarantees that the naïve confidence intervals in Step 3 will indeed have approximately correct coverage, in the sense of (9).

3 Numerical Examination of Naïve Confidence Intervals

In this section, we perform simulation studies to examine the coverage probability (9) of the naïve confidence intervals obtained by applying standard least squares software to the data ${\bm{y}},{\bm{X}}_{\hat{\mathcal{A}}_{\lambda}}$.

Recall from Section 1 that (9) involves the probability that the confidence interval contains the quantity ${\bm{\beta}}^{(\hat{\mathcal{A}}_{\lambda})}$, which in general does not equal the population regression coefficient vector ${\bm{\beta}}^{\ast}_{\hat{\mathcal{A}}_{\lambda}}$. Inference for ${\bm{\beta}}^{\ast}$ is discussed in Sections 4 and 5.

The results in this section complement simulation findings in Leeb et al., (2015).

3.1 Methods for Comparison

Following Theorem 2.3, for $\tilde{\bm{\beta}}^{(\hat{\mathcal{A}}_{\lambda})}$ defined in (8), and for each $j\in\hat{\mathcal{A}}_{\lambda}$, the 95% naïve confidence interval takes the form

[TABLE]

In order to obtain the set ${\hat{\mathcal{A}}_{\lambda}}$, we must apply the lasso using some value of $\lambda$. By (M1) we need $\lambda\succ\sqrt{\log(p)/n}$, which is slightly larger than the prediction optimal rate, $\lambda\asymp\sqrt{\log(p)/n}$ (Bickel et al.,, 2009; van de Geer and Bühlmann,, 2009). A data-driven way to obtain a larger tuning parameter it to use $\lambda_{\text{1SE}}$, which is the largest value of $\lambda$ for which the 10-fold CV prediction mean squared error (PMSE) is within one standard error of the minimum CV PMSE (see Section 7.10.1 in Hastie et al.,, 2009). However, the $\lambda$ selected by cross validation depends on ${\bm{y}}$ and may induce additional randomness in the set of selected coefficients, $\hat{\mathcal{A}}_{\lambda}$. This randomness can also impact the exact post-selection procedure of Lee et al., (2016). To address this issue, the authors proposed $\lambda_{\text{sup}}\equiv 2\mathbb{E}[\|{\bm{X}}^{\top}{\bm{e}}\|_{\infty}]/n$, where ${\bm{e}}\sim\mathcal{N}_{n}({\bm{0}},\hat{\sigma}_{\bm{\epsilon}}^{2}{\bm{I}})$. As an alternative to $\lambda_{\text{1SE}}$, we also evaluate $\lambda_{\text{sup}}$, where we simulate ${\bm{e}}$ and approximate the expectation based on the average of 1000 replicates. We compare the confidence intervals for ${\bm{\beta}}^{(\hat{\mathcal{A}}_{\lambda})}$ with those based on the exact lasso post-selection inference procedure of Lee et al., (2016), which is implemented in the R package selectiveInference.

In both approaches, the standard deviation of errors, $\sigma_{\bm{\epsilon}}$ in (1), is estimated either using the scaled lasso (Sun and Zhang,, 2012) or by applying Theorem 2.4. However, we do not examine the combination of $\lambda_{\text{sup}}$ with error variance estimated based on Theorem 2.4, because $\lambda_{\text{sup}}$ requires an estimate of the error standard deviation based on Theorem 2.4, whereas Theorem 2.4 requires a suitable choice of $\lambda$ to be valid.

3.2 Simulation Set-Up

For the simulations, we consider two partial correlation settings for ${\bm{X}}$, generated based on (i) a scale-free graph and (ii) a stochastic block model (see, e.g., Kolaczyk,, 2009), each containing $p=100$ nodes. These settings are relaxations of the simple orthogonal and block-diagonal settings, and are displayed in Figure 1.

In the scale-free graph setting, we used the igraph package in R to simulate an undirected, scale-free network $\mathcal{G}=(\mathcal{V},\mathcal{E})$ with power-law exponent parameter $\gamma=5$, and edge density 0.05. Here, $\mathcal{V}=\{1,\ldots,p\}$ is the set of nodes in the graph, and $\mathcal{E}$ is the set of edges. This resulted in a total of $|\mathcal{E}|=247$ edges in the graph. We then order the indices of the nodes in the graph so that the first, second, third, fourth, and fifth nodes correspond to the 10th, 20th, 30th, 40th, and 50th least-connected nodes in the graph.

In the stochastic block model setting, we first generate two dense Erdős-Rényi graphs (Erdős and Rényi,, 1959; Gilbert,, 1959) with five and 95 nodes, respectively. In each graph, the edge density is 0.3. We then add edges randomly between these two graphs to achieve an inter-graph edge density of 0.05. The indices of the nodes are ordered so that the nodes in the five-node graph precede the remaining nodes.

Next, for both graph settings, we define the weighted adjacency matrix, ${\bm{A}}$, as follows:

[TABLE]

where $\rho\in\{0.2,0.6\}$. We then set ${\bm{\Sigma}}={\bm{A}}^{-1}$, and standardize ${\bm{\Sigma}}$ so that $\Sigma_{(j,j)}=1$, for all $j=1,\dots,p$. We simulate observations ${\bm{x}}_{1},\ldots,{\bm{x}}_{n}\sim_{i.i.d.}\mathcal{N}_{p}({\bm{0}},{\bm{\Sigma}})$, and generate the outcome ${\bm{y}}\sim\mathcal{N}_{n}({\bm{X}}{\bm{\beta}}^{\ast},\sigma^{2}_{\bm{\epsilon}}{\bf I}_{n})$, $n\in\{100,300,500\}$, where

[TABLE]

A range of error variances $\sigma_{\bm{\epsilon}}^{2}$ are used to produce signal-to-noise ratios, $\textrm{SNR}\equiv({\bm{\beta}}^{\ast\top}{\bm{\Sigma}}{\bm{\beta}}^{\ast})/\sigma^{2}_{\bm{\epsilon}}\in\{0.1,0.3,0.5\}$.

Throughout the simulations, ${\bm{\Sigma}}$ and ${\bm{\beta}}^{\ast}$ are held fixed over $B=1000$ repetitions of the simulation study, while ${\bm{X}}$ and ${\bm{y}}$ vary.

3.3 Simulation Results

We calculate the average length and coverage proportion of the 95% naïve confidence intervals, where the coverage proportion is defined as

[TABLE]

where $\hat{\mathcal{A}}_{\lambda}^{b}$ and $\mathrm{CI}_{j}^{(\hat{\mathcal{A}}_{\lambda}^{b}),b}$ are the set of variables selected by the lasso in the $b$th repetition, and the 95% naïve confidence interval (16) for the $j$th variable in the $b$th repetition, respectively. Recall that $\beta_{j}^{(\hat{\mathcal{A}}_{\lambda}^{b})}$ was defined in (5). In order to calculate the average length and coverage proportion associated with the exact lasso post selection procedure of Lee et al., (2016), we replace $\mathrm{CI}^{(\hat{\mathcal{A}}_{\lambda}^{b}),b}_{j}$ in (18) with the confidence interval output by the selectiveInference R package.

Tables 1 and 2 show the coverage proportion and average length of 95% naïve confidence intervals and 95% exact lasso post-selection confidence intervals under the scale-free graph and stochastic block model settings, respectively. The result shows that the coverage probability of the exact post-selection confidence interval is more correct than that of the naïve confidence interval when the data are small and signal is weak. But when the data are large and/or signal is relatively strong, both confidence intervals have approximately correct coverage. This corroborates the findings in Leeb et al., (2015), in which the authors consider settings with $n=30$ and $p=10$. The coverage probability of the naïve confidence intervals with tuning parameter $\lambda_{\text{1SE}}$ is a bit smaller than the desired level, especially when the signal is weak. This may be due to the randomness in $\lambda_{\text{1SE}}$. The naïve approach of error variance estimation as in Theorem 2.4 works similarly compared to the scaled lasso across all settings. In addition, Tables 1 and 2 also show that naïve confidence intervals are substantially narrower than exact lasso post-selection confidence intervals, especially when the signal is weak.

In addition, to evaluate whether $\hat{\mathcal{A}}_{\lambda}$ is deterministic, among repetitions that $\hat{\mathcal{A}}_{\lambda}^{b}\neq\emptyset$ (there is no confidence interval if $\hat{\mathcal{A}}_{\lambda}^{b}=\emptyset$), we also calculate the proportion of $\hat{\mathcal{A}}_{\lambda}^{b}=\mathcal{D}$, where $\mathcal{D}$ is the most common $\hat{\mathcal{A}}_{\lambda}^{b}$, $b=1,\dots,1000$. The result is summarized in Table 3, which shows that $\hat{\mathcal{A}}_{\lambda}$ is almost deterministic with tuning parameter $\lambda_{\text{sup}}$. With $\lambda_{\text{1SE}}$, due to the randomness in the tuning parameter, $\hat{\mathcal{A}}_{\lambda}$ is less deterministic, which may explain the smaller coverage probability than the desired level in this case. We also examined the proportion of $\hat{\mathcal{A}}_{\lambda}^{b}=\mathcal{A}_{\lambda}^{b}$, which is shown in Table 4. The observation is similar to that in Table 3. As mentioned in Section 1, $\Pr[\hat{\mathcal{A}}_{\lambda}=\mathcal{A}^{\ast}]=0$ in all the settings considered.

4 Inference for $\beta^{\ast}$ With the Naïve Score Test

Sections 2 and 3 focused on the task of developing confidence intervals for ${\bm{\beta}}^{(\mathcal{M})}$ in (3), where $\mathcal{M}=\hat{\mathcal{A}}_{\lambda}$, the set of variables selected by the lasso. However, recall from (5) that typically ${\bm{\beta}}^{(\mathcal{M})}\neq{\bm{\beta}}^{\ast}_{\mathcal{M}}$, where ${\bm{\beta}}^{\ast}$ was introduced in (1).

In this section, we shift our focus to performing inference on ${\bm{\beta}}^{\ast}$. We will exploit Proposition 2.2 to develop a simple approach for testing $H_{0,j}^{\ast}:\beta_{j}^{\ast}=0$, for $j=1,\ldots,p$.

Recall that in the low-dimensional setting, the classical score statistic for the hypothesis $H_{0,j}^{\ast}:\beta_{j}^{\ast}=0$ is proportional to ${\bm{x}}_{j}^{T}({\bm{y}}-\hat{\bm{y}}^{0})$, where $\hat{\bm{y}}^{0}$ is the vector of fitted values that results from least squares regression of ${\bm{y}}$ onto the $p-1$ features ${\bm{x}}_{1},\ldots,{\bm{x}}_{j-1},{\bm{x}}_{j+1},\ldots,{\bm{x}}_{p}$. In order to adapt the classical score test statistic to the high-dimensional setting, we define the naïve score test statistic for testing $H_{0,j}^{\ast}:\beta_{j}^{\ast}=0$ as

[TABLE]

where

[TABLE]

and

[TABLE]

is the orthogonal projection matrix onto the set of variables in ${\hat{\mathcal{A}}_{\lambda}\backslash\{j\}}$. $\tilde{\bm{\beta}}^{(\hat{\mathcal{A}}_{\lambda}\backslash\{j\})}$ is defined in (8). In (19), the notation $\hat{\mathcal{A}}_{\lambda}\backslash\{j\}$ represents the set $\hat{\mathcal{A}}_{\lambda}$ in (6) with $j$ removed, if $j\in\hat{\mathcal{A}}_{\lambda}$. If $j\notin\hat{\mathcal{A}}_{\lambda}$, then $\hat{\mathcal{A}}_{\lambda}\backslash\{j\}=\hat{\mathcal{A}}_{\lambda}$.

In Theorem 4.1, we will derive the asymptotic distribution of $S^{j}$ under $H_{0,j}^{\ast}:\beta_{j}^{\ast}=0$. We first introduce two new conditions.

First, we require that the total signal strength of variables not selected by the noiseless lasso, (7), is small.

• (M2∗)

Recall that $\mathcal{A}^{\ast}\equiv\mathrm{supp}({\bm{\beta}}^{\ast})$. Let $\mathcal{A}_{\lambda}\equiv\mathrm{supp}({\bm{\beta}}_{\lambda})$, $b_{\lambda\min}\equiv\min_{j\in\mathcal{A}_{\lambda}}|\beta_{\lambda,j}|$ with ${\bm{\beta}}_{\lambda}$ defined in (7). The signal strength in $\mathcal{A}_{\lambda}$ satisfies

[TABLE]

and the signal strength outside of $\mathcal{A}_{\lambda}$ satisfies

[TABLE]

Condition (M2∗) closely resembles (M2), which was required for Theorem 2.3 in Section 2. The only difference between the two is that (M2∗) requires \|{\bm{X}}_{\mathcal{A}_{\lambda}^{c}}{\bm{\beta}}^{\ast}_{\mathcal{A}_{\lambda}^{c}}\|_{2}=\mathchoice{{\scriptstyle\mathcal{O}}}{{\scriptstyle\mathcal{O}}}{{\scriptscriptstyle\mathcal{O}}}{\scalebox{0.7}{\scriptscriptstyle\mathcal{O}}}(1), whereas (M2) requires only that $\|{\bm{X}}_{\mathcal{A}_{\lambda}^{c}}{\bm{\beta}}^{\ast}_{\mathcal{A}_{\lambda}^{c}}\|_{2}=\mathcal{O}\left(\sqrt{\log(p)}\right)$. Recall that in Section 2, we consider inference for the parameters in the sub-model (3). In other words, testing the population regression parameter ${\bm{\beta}}^{\ast}$ in (1) requires more stringent assumptions than constructing confidence intervals for the parameters in the sub-model (3).

The following condition, required to apply the Lindeberg-Feller Central Limit Theorem, can be relaxed if the noise ${\bm{\epsilon}}$ in (1) is normally distributed.

• (S)

$\lambda$, ${\bm{\beta}}^{\ast}$ and ${\bm{X}}$ satisfy $\lim_{n\to\infty}\|{\bf r}^{s}\|_{\infty}/\|{\bf r}^{s}\|_{2}=0$, where ${\bf r}^{s}\equiv\left({\bf I}_{n}-{\bf P}^{(\mathcal{A}_{\lambda}\backslash\{j\})}\right){\bm{x}}_{j}$.

We now present Theorem 4.1, which is proven in Section E of the online Supplementary Materials.

Theorem 4.1**.**

Suppose (M2∗) and (S) hold and $\lim_{n\to\infty}\Pr\left[\hat{\mathcal{A}}_{\lambda}=\mathcal{A}_{\lambda}\right]=1$. For any $j=1,\dots,p$, under the null hypothesis $H_{0,j}^{\ast}:\beta^{\ast}_{j}=0$,

[TABLE]

*where $S^{j}$ was defined in (19), and where $\sigma_{\bm{\epsilon}}$ is the variance of ${\bm{\epsilon}}$ in (1). *

Theorem 4.1 states that the distribution of the naïve score test statistic $S^{j}$ is asymptotically the same as if $\hat{\mathcal{A}}_{\lambda}$ were a fixed set, as opposed to being selected by fitting a lasso model on the data. Based on (22), we reject the null hypothesis $H^{\ast}_{0,j}:\beta^{\ast}_{j}=0$ at level $\alpha>0$ if $|T|>\Phi_{\mathcal{N}}^{-1}(1-\alpha/2)$, where $\Phi^{-1}_{\mathcal{N}}(\cdot)$ is the quantile function of the standard normal distribution function.

We emphasize that Theorem 4.1 holds for any variable $j=1,\dots,p$, and thus can be used to test $H_{0,j}^{\ast}:\beta_{j}^{\ast}=0$, for all $j=1,\dots,p$. (This is in contrast to Theorem 2.3, which concerns confidence intervals for the parameters in the sub-model (3) consisting of the variables in $\hat{\mathcal{A}}_{\lambda}$, and hence holds only for $j\in\hat{\mathcal{A}}_{\lambda}$.)

5 Numerical Examination of the Naïve Score Test

In this section, we compare the performance of the naïve score test (19) to three recent proposals from the literature for testing $H_{0,j}^{\ast}:\beta_{j}^{\ast}=0$: namely, LDPE (Zhang and Zhang,, 2014; van de Geer et al.,, 2014), SSLasso (Javanmard and Montanari, 2014a, ), and the decorrelated score test (dScore; Ning and Liu,, 2017). (Since with high probability $\hat{\mathcal{A}}_{\lambda}\neq\mathcal{A}^{\ast}$, we do not include the exact post-selection procedure in this comparison.) R code for SSLasso, and dScore was provided by the authors; LDPE is implemented in the R package hdi. For the naïve score test, we estimate $\sigma_{\bm{\epsilon}}$, the standard deviation of the errors in (1), using the scaled lasso (Sun and Zhang,, 2012) or Theorem 2.4.

All four of these methods require us to select the value of the lasso tuning parameter. For LDPE, SSLasso, and dScore, we use 10-fold cross-validation to select the tuning parameter value that produces the smallest cross-validated mean square error, $\lambda_{\text{min}}$. As in the numerical study of the naïve confidence intervals in Section 3, we implement the naïve score test using the tuning parameter value $\lambda_{\text{1SE}}$ and $\lambda_{\text{sup}}$. Unless otherwise noted, all tests are performed at a significance level of 0.05.

In Section 5.1, we investigate the powers and type-I errors of the above tests in simulation experiments. Section 5.2 contains an analysis of a glioblastoma gene expression dataset.

5.1 Power and Type-I Error

5.1.1 Simulation Set-Up

In this section, we adapt the scale-free graph and the stochastic block model presented in Section 3.2 to have $p=500$.

In the scale-free graph setting, we generate a scale-free graph with $\gamma=5$, edge density 0.05, and $p=500$ nodes. The resulting graph has $|\mathcal{E}|=6237$ edges. We order the nodes in the graph so that $j$th node is the $(30\times j)$th least-connected node in the graph, for $1\leq j\leq 10$. For example, the 4th node is the 120th least-connected node in the graph.

In the stochastic block model setting, we generate two dense Erdős-Rényi graphs with ten nodes and 490 nodes, respectively; each has an intra-graph edge density of 0.3. The node indices are ordered so that the nodes in the smaller graph precede those in the larger graph. We then randomly connect nodes between the two graphs in order to obtain an inter-graph edge density of 0.05.

Next, for both graph settings, we generate ${\bm{A}}$ as in (17), where $\rho\in\{0.2,0.6\}$. We then set ${\bm{\Sigma}}={\bm{A}}^{-1}$, and standardize ${\bm{\Sigma}}$ so that $\Sigma_{(j,j)}=1$, for all $j=1,\dots,p$. We simulate observations ${\bm{x}}_{1},\ldots,{\bm{x}}_{n}\sim_{i.i.d.}\mathcal{N}_{p}({\bm{0}},{\bm{\Sigma}})$, and generate the outcome ${\bm{y}}\sim\mathcal{N}_{n}({\bm{X}}{\bm{\beta}}^{\ast},\sigma^{2}_{\bm{\epsilon}}{\bf I}_{n})$, $n\in\{100,200,400\}$, where

[TABLE]

A range of error variances $\sigma_{\bm{\epsilon}}^{2}$ are used to produce signal-to-noise ratios, $\textrm{SNR}\equiv({\bm{\beta}}^{\ast\top}{\bm{\Sigma}}{\bm{\beta}}^{\ast})/\sigma^{2}_{\bm{\epsilon}}\in\{0.1,0.3,0.5\}$.

We hold ${\bm{\Sigma}}$ and ${\bm{\beta}}^{\ast}$ fixed over $B=100$ repetitions of the simulation, while ${\bm{X}}$ and ${\bm{y}}$ vary.

5.1.2 Simulation Results

For each test, the average power on the strong signal variables, the average power on the weak signal variables, and the average type-I error rate are defined as

[TABLE]

respectively. In (23)–(25), $p_{jb}$ is the $p$-value associated with null hypothesis $H^{\ast}_{0,j}:\beta_{j}^{\ast}=0$ in the $b$th simulated data set. In the simulations, the graphs and ${\bm{\beta}}^{\ast}$ are held fixed over $B=100$ repetitions of the simulation study, while ${\bm{X}}$ and ${\bm{y}}$ vary.

Tables 5 and 6 summarize the results in the two simulation settings. Naïve score test with $\lambda_{\text{sup}}$ has slightly worse control of type-1 error rate and better power than the other four methods, which have approximate control over the type-I error rate and comparable power. The performance with the scaled lasso is similar to the performance of Theorem 2.4.

5.2 Application to Glioblastoma Data

We investigate a glioblastoma gene expression data set previously studied in Horvath et al., (2006). For each of $130$ patients, a survival outcome is available; we removed the twenty patients who were still alive at the end of the study. This resulted in a data set with $n=110$ observations. The gene expression measurements were normalized using the method of Gautier et al., (2004). We limited our analysis to $p=3600$ highly-connected genes (Zhang and Horvath,, 2005; Horvath and Dong,, 2008). The normalized data can be found at the website of Dr. Steve Horvath of UCLA Biostatistics. We log-transformed the survival response and centered it to have mean zero. Furthermore, we log-transformed the expression data, and then standardized each gene to have mean zero and standard deviation one across the $n=110$ observations.

Our goal is to identify individual genes whose expression levels are associated with survival time, after adjusting for the other $3599$ genes in the data set. With family-wise error rate (FWER) controlled at level 0.1 using the Holm procedure (Holm,, 1979), the naïve score test identifies three such genes: CKS2, H2AFZ, and RPA3. You et al., (2015) observed that CKS2 is highly expressed in glioma. Vardabasso et al., (2014) found that histone genes, of which H2AFZ is one, are related to cancer progression. Jin et al., (2015) found that RPA3 is associated with glioma development. As a comparison, SSLasso finds two genes associated with patient survival: PPAP2C and RGS3. LDPE and dScore identify no genes at FWER of $0.1$.

6 Discussion

In this paper, we examined a very naïve two-step approach to high-dimensional inference:

Perform the lasso in order to select a small set of variables, $\hat{\mathcal{A}}_{\lambda}$. 2. 2.

Fit a least squares regression model using just the variables in $\hat{\mathcal{A}}_{\lambda}$, and make use of standard regression inference tools. Make no adjustment for the fact that $\hat{\mathcal{A}}_{\lambda}$ was selected based on the data.

It seems clear that this naïve approach is problematic, since we have peeked at the data twice, but are not accounting for this double-peeking in our analysis.

In this paper, we have shown that under certain assumptions, $\hat{\mathcal{A}}_{\lambda}$ converges with high probability to a deterministic set, $\mathcal{A}_{\lambda}$. A similar result for random design matrix is presented in Zhao et al., (2019). This key insight allows us to establish that the confidence intervals resulting from the aforementioned naïve two-step approach have asymptotically correct coverage, in the sense of (4). This constitutes a theoretical justification for the recent simulation findings of Leeb et al., (2015). Furthermore, we used this key insight in order to establish that the score test that results from the naïve two-step approach has asymptotically the same distribution as though the selected set of variables had been fixed in advance; thus, it can be used to test the null hypothesis $H_{0,j}^{\ast}:\beta^{\ast}_{j}=0$, $j=1,\dots,p$.

Our simulation results corroborate our theoretical findings. In fact, we find essentially no difference between the empirical performance of these naïve proposals, and a host of other recent proposals in the literature for high-dimensional inference (Javanmard and Montanari, 2014a, ; Zhang and Zhang,, 2014; van de Geer et al.,, 2014; Lee et al.,, 2016; Ning and Liu,, 2017).

From a bird’s-eye view, the recent literature on high-dimensional inference falls into two camps. The work of Wasserman and Roeder, (2009); Meinshausen et al., (2009); Berk et al., (2013); Lee et al., (2016); Tibshirani et al., (2016) focuses on performing inference on the sub-model (3), whereas the work of Javanmard and Montanari, (2013); Javanmard and Montanari, 2014a ; Javanmard and Montanari, 2014b ; Zhang and Zhang, (2014); van de Geer et al., (2014); Zhao and Shojaie, (2016); Ning and Liu, (2017) focuses on testing hypotheses associated with (1). In this paper, we have shown that the confidence intervals that result from the naïve approach can be used to perform inference on the sub-model (3), whereas the score test that results from the naïve approach can be used to test hypotheses associated with (1).

In the era of big data, simple analyses that are easy to apply and easy to understand are especially attractive to scientific investigators. Therefore, a careful investigation of such simple approaches is worthwhile, in order to determine which ones have the potential to yield accurate results, and which do not. We do not advocate applying the naïve two-step approach described above in most practical data analysis settings: we are confident that in practice, our intuition is correct, and this approach will perform poorly when the sample size is small or moderate, and/or the assumptions, which are unfortunately unverifiable, are not met. However, in very large data settings, our results suggest that this naïve approach may indeed be viable for high-dimensional inference, or at least warrants further investigation.

When choosing among existing inference procedures based on lasso, the target of inference should be taken into consideration. The target of inference can either be the population parameters, ${\bm{\beta}}^{\ast}$ in (1), or the parameters induced by the sub-model chosen by lasso, ${\bm{\beta}}^{(\mathcal{M})}$ in (5). Sample-splitting (Wasserman and Roeder,, 2009; Meinshausen et al.,, 2009) and exact post selection (Lee et al.,, 2016; Tibshirani et al.,, 2016) methods provide valid inferences for ${\bm{\beta}}^{(\mathcal{M})}$. The latter is a particularly appealing choice for inference on ${\bm{\beta}}^{(\mathcal{M})}$, as it provides non-asymptotic confidence intervals under minimal assumptions. However, as we discussed in Section 1, ${\bm{\beta}}^{(\mathcal{M})}$ is, in general, different from ${\bm{\beta}}^{\ast}_{\mathcal{M}}$. A set of sufficient conditions for ${\bm{\beta}}^{(\mathcal{M})}={\bm{\beta}}^{\ast}_{\mathcal{M}}$ is the irrepresentable condition together with a beta-min condition. Unfortunately, these assumptions are unverifiable and may not hold in practice. Our theoretical analysis and empirical studies suggests that the naïve two-step approach described above facilitates inference for ${\bm{\beta}}^{\ast}_{\mathcal{M}}$ under less stringent assumptions and without any conditioning or sample splitting. However, this method is also asymptotic and relies on unverifiable assumptions. Debiased lasso tests (Zhang and Zhang,, 2014; van de Geer et al.,, 2014; Javanmard and Montanari,, 2013; Javanmard and Montanari, 2014a, ; Ning and Liu,, 2017) provide asymptotically valid inference for entries of ${\bm{\beta}}^{\ast}$, without requiring a beta-min or irrepresentable condition. However, they require more restrictive sparsity of ${\bm{\beta}}^{\ast}$, as well as sparsity of the inverse covariance matrix of covariates, ${\bm{\Sigma}}^{-1}$, which are also unverifiable. These limitations underscore the importance of recent efforts to relax these sparsity assumptions (e.g. Zhu et al.,, 2018; Wang et al.,, 2020).

We close with some suggestions for future research. One reviewer brought up an interesting point: methods with folded-concave penalties (e.g., Fan and Li,, 2001; Zhang,, 2010) require milder conditions to achieve variable selection consistency, i.e., $\Pr[\hat{\mathcal{A}}_{\lambda}=\mathcal{A}^{\ast}]\to 1$, than the lasso. Inspired by this observation, we wonder whether the proposal of Fan and Li, (2001) and Zhang, (2010) also require milder conditions to achieve $\Pr[\hat{\mathcal{A}}_{\lambda}=\mathcal{A}_{\lambda}]\to 1$. If so, then we could replace the lasso with a folded-concave penalty in the variable selection step, and improve the robustness of the naïve approaches. We believe this could be a fruitful area of future research. In addition, extending the proposed theory and methods to generalized linear models and M-estimators may also be promising areas for future research.

Appendix A Proof of Proposition 2.2

We first state and prove Lemmas A.1–A.5, which are required to prove Proposition 2.2.

Lemma A.1**.**

Suppose (M1) holds. Then, ${\bm{\beta}}_{\lambda}$ and $\hat{\bm{\beta}}_{\lambda}$ as defined in (7) and (2), respectively, are unique.

Proof.

The proof follows from Lemma 3 of Tibshirani, (2013) and is thus omitted. ∎

Lemma A.2**.**

Suppose (M1) holds. Then,

[TABLE]

Proof.

It is well understood that with an application of the union bound to the standard exponential inequality, (M1) implies

[TABLE]

for some constant $C$ not depending on $p$ or $n$ (Vershynin,, 2012). ∎

Lemma A.3**.**

Suppose (M1), (M2) and (E) hold. Then $q_{\lambda}=\mathcal{O}(q^{\ast})$.

Proof.

Recall that $b_{\lambda\min}\equiv\min_{j\in\mathcal{A}_{\lambda}}|\beta_{\lambda,j}|$. Hence

[TABLE]

With (E), Lemma 2.1 in van de Geer and Bühlmann, (2009) shows that

[TABLE]

i.e.,

[TABLE]

Since $b_{\lambda\min}\succeq\lambda$, the above equation implies that $|\mathcal{A}_{\lambda}\backslash\mathcal{A}^{\ast}|=\mathcal{O}(q^{\ast})$, which completes the proof. ∎

Lemma A.4**.**

Suppose (M1) and (E) hold. Then, the estimator $\hat{\bm{\beta}}_{\lambda}$ defined in (2) and its population version, ${\bm{\beta}}_{\lambda}$, defined in (7) satisfy

[TABLE]

Proof.

Theorem 2.1 in van de Geer, (2017) shows that under (M1),

[TABLE]

Note there are three differences between (29) and the formulation in Theorem 2.1 in van de Geer, (2017). First, van de Geer, (2017) assumes $\|{\bm{x}}_{j}\|_{2}^{2}=1,j=1,\dots,p$, whereas we assume $\|{\bm{x}}_{j}\|_{2}^{2}=n$. This difference leads to the change of a $\sqrt{n}$ factor in the first term in (29). Second, Theorem 2.1 in van de Geer, (2017) factors out a $\|\hat{\bm{\Sigma}}\|_{2}$ from the square-root, whereas (29) does not have this step, and consequently has a power of two on $\hat{\bm{\Sigma}}$ within the square-root; see the proof of Theorem 15.2 in van de Geer, (2017) for more details. Finally, although Theorem 2.1 in van de Geer, (2017) assumes Gaussian random errors ${\bm{\epsilon}}$, it continues to hold for sub-Gaussian errors. This is because Lemma 15.5 in van de Geer, (2017) can be proven with sub-Gaussian data as shown in Hsu et al., (2012).

Let $q\equiv|\mathcal{A}_{\lambda}\backslash\mathcal{A}^{\ast}|$. Then,

[TABLE]

With (E), Lemma 2.1 in van de Geer and Bühlmann, (2009) shows that

[TABLE]

Thus,

[TABLE]

Now, since $q\asymp q^{\ast}$ by Lemma A.3, $\phi^{2}(q)=\mathcal{O}\left(\sqrt{\log(p)/q^{\ast}}\right)$ by (E). Thus, (32) implies that

[TABLE]

Let $\mathcal{I}=\mathrm{supp}(\hat{\bm{\beta}}_{\lambda}-{\bm{\beta}}_{\lambda})$. Then, $|\mathcal{I}|=\|\hat{\bm{\beta}}_{\lambda}-{\bm{\beta}}_{\lambda}\|_{0}=\mathcal{O}_{p}(q^{\ast})$. Moreover, $\|\hat{\bm{\beta}}_{\lambda,\mathcal{I}^{c}}-{\bm{\beta}}_{\lambda,\mathcal{I}^{c}}\|_{1}=0\leq\|\hat{\bm{\beta}}_{\lambda,\mathcal{I}}-{\bm{\beta}}_{\lambda,\mathcal{I}}\|_{1}$. Thus, by the restricted eigenvalue condition in (E) that $\phi^{\ast 2}$ is a positive constant,

[TABLE]

which completes the proof. ∎

Lemma A.5**.**

Suppose (M1), (M2) and (T) hold. Then

[TABLE]

where

[TABLE]

is the stationary condition of (7).

Proof.

Rearranging terms of (35), we get

[TABLE]

and

[TABLE]

where ${\bm{P}}^{(\mathcal{A}_{\lambda})}\equiv{\bm{X}}_{\mathcal{A}_{\lambda}}({\bm{X}}_{\mathcal{A}_{\lambda}}^{\top}{\bm{X}}_{\mathcal{A}_{\lambda}})^{-1}{\bm{X}}^{\top}_{\mathcal{A}_{\lambda}}$ is the projection matrix onto the column space of ${\bm{X}}_{\mathcal{A}_{\lambda}}$. Therefore,

[TABLE]

where the second inequality is due to the irrepresentable condition in (T), and the last equality is due to the assumptions that ${\bm{x}}_{j}^{\top}{\bm{x}}_{j}=n$ in (M1) and $\|{\bm{X}}_{{\bm{A}}_{\lambda}^{c}}{\bm{\beta}}^{\ast}_{{\bm{A}}_{\lambda}^{c}}\|_{2}=\mathcal{O}(\sqrt{\log(p)})$ in (M2). Given that $\sqrt{\log(p)}/(\sqrt{n}\lambda)\to 0$ by (M1), we have

[TABLE]

∎

Proof of Proposition 2.2.

We first prove

[TABLE]

According to the stationary conditions of (2),

[TABLE]

Together with (35), we get

[TABLE]

We now bound both terms on the right hand side of (41). By Lemma A.2,

[TABLE]

In addition, Lemma A.4 shows that

[TABLE]

Therefore,

[TABLE]

Therefore, it follows from (41) that $\|{\bm{\tau}}_{\lambda}-\hat{\bm{\tau}}_{\lambda}\|_{\infty}=\mathcal{O}_{p}\left(\sqrt{\log(p)/n}/\lambda\right)$. Given that $\sqrt{\log(p)/n}/\lambda\to 0$ by (M1), we have \|{\bm{\tau}}_{\lambda}-\hat{\bm{\tau}}_{\lambda}\|_{\infty}=\mathchoice{{\scriptstyle\mathcal{O}}}{{\scriptstyle\mathcal{O}}}{{\scriptscriptstyle\mathcal{O}}}{\scalebox{0.7}{\scriptscriptstyle\mathcal{O}}}_{p}(1). Since $\lim_{n\to\infty}\|{\bm{\tau}}_{\lambda,\mathcal{A}_{\lambda}^{c}}\|_{\infty}\leq 1-\delta$ by Lemma A.5, we have $\lim_{n\to\infty}\Pr[\|\hat{\bm{\tau}}_{\lambda,\mathcal{A}^{c}_{\lambda}}\|_{\infty}<1]=1$, and

[TABLE]

To prove the other direction, we assume without loss of generality that $\mathcal{A}_{\lambda}\neq\emptyset$; otherwise the result holds trivially. Now, (M2) and the fact that $\sqrt{\log(p)}/(\sqrt{n}\lambda)\to 0$ by (M1) imply that

[TABLE]

But, by Lemma A.4, $\|\hat{\bm{\beta}}_{\lambda}-{\bm{\beta}}_{\lambda}\|_{2}=\mathcal{O}_{p}\left(\sqrt{\log(p)/n}\right)$. Thus, for any $\xi>0$, there exists a constant $C>0$, not depending on $n$, such that for $n$ sufficiently large,

[TABLE]

Moreover, by (43), for $n$ sufficiently large, $b_{\lambda\min}>C\sqrt{\log(p)/n}$, where $b_{\lambda\min}\equiv\min_{j\in\mathcal{A}_{\lambda}}|\beta_{\lambda,j}|$. Combining (43) and (44), for $n$ sufficiently large, whenever $|\beta_{\lambda,j}|>0$, $|\beta_{\lambda,j}|>C\sqrt{\log(p)/n}$ and hence $\Pr[|\hat{\beta}_{\lambda,j}|>0]>1-\xi$. Therefore

[TABLE]

which completes the proof. ∎

Appendix B Proof of Lemma 2.1

Proof of Lemma 2.1.

We rewrite condition (T) as

[TABLE]

With (E), Lemma 2.1 in van de Geer and Bühlmann, (2009) shows that

[TABLE]

Based on (38), (M2) implies that

[TABLE]

Thus,

[TABLE]

Given that $\sqrt{\log(p)}/(\sqrt{n}\lambda)\to 0$ by (M1), condition (T) is satisfied if

[TABLE]

which completes the proof. ∎

Appendix C Proof of Theorem 2.3

Proof of Theorem 2.3.

By Proposition 2.2, $\Pr[\hat{\mathcal{A}}_{\lambda}=\mathcal{A}_{\lambda}]\to 1$. Therefore, with probability tending to one,

[TABLE]

Thus,

[TABLE]

We proceed to prove the asymptotic distribution of $[({\bm{X}}_{\mathcal{A}_{\lambda}}^{\top}{\bm{X}}_{\mathcal{A}_{\lambda}})^{-1}{\bm{X}}_{\mathcal{A}_{\lambda}}^{\top}{\bm{\epsilon}}]_{j}$. Dividing it by its standard deviation, $\sigma_{\bm{\epsilon}}\sqrt{[({\bm{X}}_{\mathcal{A}_{\lambda}}^{\top}{\bm{X}}_{\mathcal{A}_{\lambda}})^{-1}]_{(j,j)}}$, where $\sigma_{\bm{\epsilon}}$ is the error standard deviation, we get

[TABLE]

where ${\bm{r}}^{w}\equiv{\bm{e}}^{j}({\bm{X}}^{\top}_{\mathcal{A}_{\lambda}}{\bm{X}}_{\mathcal{A}_{\lambda}})^{-1}{\bm{X}}^{\top}_{\mathcal{A}_{\lambda}}\in\mathbb{R}^{n}$, and ${\bm{e}}^{j}$ is the row vector of length $|\mathcal{A}_{\lambda}|$ with the entry that corresponds to $\beta^{\ast}_{j}$ equal to one, and zero otherwise. In order to use the Lindeberg-Feller Central Limit Theorem to prove the asymptotic normality of (49), we need to show that the Lindeberg’s condition holds, i.e.,

[TABLE]

Given that $|r^{w}_{i}|\leq\|{\bf r}^{w}\|_{\infty}$, and that the $\epsilon_{i}$’s are identically distributed,

[TABLE]

Since $\|{\bf r}^{w}\|_{\infty}/\|{\bf r}^{w}\|_{2}\rightarrow 0$ by Condition (W), $\epsilon_{1}^{2}\mathbf{1}\left\{|\epsilon_{1}|\|{\bf r}^{w}\|_{\infty}/(\sigma_{\bm{\epsilon}}\|{\bf r}^{w}\|_{2})>\eta\right\}\rightarrow_{p}0$. Thus, because $\epsilon_{1}^{2}\geq\epsilon_{1}^{2}\mathbf{1}\left\{|\epsilon_{1}|\|{\bf r}^{w}\|_{\infty}/(\sigma_{\bm{\epsilon}}\|{\bf r}^{w}\|_{2})>\eta\right\}$ with probability one and $\mathbb{E}[\epsilon_{1}^{2}]=\sigma_{\epsilon}^{2}<\infty$, we use $\epsilon_{1}^{2}$ as the dominant random variable, and apply the Dominated Convergence Theorem,

[TABLE]

which gives the Lindeberg’s condition.

Thus,

[TABLE]

Using, again, the fact that by Proposition 2.2, $\lim_{n\to\infty}\Pr\left[\mathcal{A}_{\lambda}=\hat{\mathcal{A}}_{\lambda}\right]=1$, we can write

[TABLE]

∎

Appendix D Proof of Theorem 2.4

Proof.

Based on Proposition 2.2 that $\Pr[\hat{\mathcal{A}}_{\lambda}=\mathcal{A}_{\lambda}]\to 1$, we have

[TABLE]

where $\hat{q}_{\lambda}\equiv|\hat{\mathcal{A}}_{\lambda}|$, $q_{\lambda}\equiv|\mathcal{A}_{\lambda}|$. Denoting ${\bm{P}}^{(\mathcal{A}_{\lambda})}\equiv{\bm{X}}_{\mathcal{A}_{\lambda}}({\bm{X}}_{\mathcal{A}_{\lambda}}^{\top}{\bm{X}}_{\mathcal{A}_{\lambda}})^{-1}{\bm{X}}^{\top}_{\mathcal{A}_{\lambda}}$,

[TABLE]

where (52) is based on the fact that ${\bm{\beta}}^{\ast}_{-\mathcal{A}^{\ast}}\equiv{\bm{0}}$. To simplify notations, denote ${\bm{\theta}}\equiv{\bm{X}}_{\mathcal{A}^{\ast}\backslash\mathcal{A}_{\lambda}}{\bm{\beta}}^{\ast}_{\mathcal{A}^{\ast}\backslash\mathcal{A}_{\lambda}}$ and ${\bm{Q}}\equiv{\bm{I}}-{\bm{P}}^{(\mathcal{A}_{\lambda})}$. Expanding (52),

[TABLE]

Because ${\bm{Q}}$ is an idempotent matrix, ${\bm{Q}}$ is positive semidefinite, whose eigenvalues are all zeros and ones. Thus,

[TABLE]

where the last equality is based on (M2). Since $\log(p)/(n-q_{\lambda})\to 0$,

[TABLE]

which means that

[TABLE]

Therefore,

[TABLE]

We now derive the expected value and variance of $\|{\bm{y}}-{\bm{X}}_{\mathcal{A}_{\lambda}}\tilde{\bm{\beta}}^{(\mathcal{A}_{\lambda})}\|_{2}^{2}/(n-q_{\lambda})$. First, because $\mathbb{E}[{\bm{\epsilon}}]={\bm{0}}$ and $\mathrm{Cov}[{\bm{\epsilon}}]=\sigma_{\bm{\epsilon}}^{2}{\bm{I}}$, using expectation of quadratic forms,

[TABLE]

Because ${\bm{Q}}$ is an idempotent matrix, we have $\mathrm{tr}[{\bm{Q}}]=n-q_{\lambda}$. Thus

[TABLE]

We now calculate the variance of $\|{\bm{y}}-{\bm{X}}_{\mathcal{A}_{\lambda}}\tilde{\bm{\beta}}^{(\mathcal{A}_{\lambda})}\|_{2}^{2}/(n-q_{\lambda})$. Since

[TABLE]

where the second term is derived in (56). We now derive the expected value of $\|{\bm{y}}-{\bm{X}}_{\mathcal{A}_{\lambda}}\tilde{\bm{\beta}}^{(\mathcal{A}_{\lambda})}\|_{2}^{4}/(n-q_{\lambda})^{2}$. Based on (54),

[TABLE]

We consider each term in the above formulation separately. First,

[TABLE]

where the last equality is based on (55) and (56). For the second term, based on (53),

[TABLE]

For the third term,

[TABLE]

where, based on (53),

[TABLE]

and

[TABLE]

The last equality is based on (55) and (56). Since $-1\leq\mathrm{Cor}[{\bm{\theta}}^{\top}{\bm{Q}}{\bm{\epsilon}},{\bm{\epsilon}}^{\top}{\bm{Q}}{\bm{\epsilon}}]\leq 1$, based on (61) and (62), we have

[TABLE]

Therefore, collecting (59), (60), (61) and (62), we have

[TABLE]

But,

[TABLE]

Because ${\bm{\epsilon}}$ is independent and identically distributed with $\mathbb{E}[\epsilon_{1}^{2}]=\sigma_{\bm{\epsilon}}^{2}$,

[TABLE]

Moreover, because $(\sum_{i=1}^{n}Q_{(i,i)})^{2}=\sum_{i\neq j}^{n}Q_{(i,i)}Q_{(j,j)}+\sum_{i=1}^{n}Q_{(i,i)}^{2}$,

[TABLE]

where the inequality is based on Jensen’s inequality that $\mathbb{E}[\epsilon_{1}^{2}]^{2}\leq\mathbb{E}[\epsilon_{1}^{4}]$. Because $\sum_{i,j}^{n}Q_{(i,j)}^{2}=\mathrm{tr}[{\bm{Q}}^{2}]=\mathrm{tr}[{\bm{Q}}]=n-q_{\lambda}$, we have

[TABLE]

Since $\epsilon_{1}$ has sub-Gaussian tails, we have $\mathbb{E}[\epsilon_{1}^{4}]=\mathcal{O}(1)$ (see, e.g., Lemma 5.5 in Vershynin,, 2012). Thus, based on (63),

[TABLE]

and hence, based on (57),

[TABLE]

Finally, applying Chebyshev’s inequality, we obtain

[TABLE]

∎

Appendix E Proof of Theorem 4.1

Proof of Theorem 4.1.

By Proposition 2.2, $\Pr\left[\hat{\mathcal{A}}_{\lambda}=\mathcal{A}_{\lambda}\right]\to 1$. Therefore, we also have \Pr\left[\big{(}\hat{\mathcal{A}}_{\lambda}\backslash\{j\}\big{)}=\big{(}\mathcal{A}_{\lambda}\backslash\{j\}\big{)}\right]\to 1, and with probability tending to one,

[TABLE]

where ${\bf P}^{(\mathcal{A}_{\lambda}\backslash\{j\})}\equiv{\bm{X}}_{\mathcal{A}_{\lambda}\backslash\{j\}}({\bm{X}}_{\mathcal{A}_{\lambda}\backslash\{j\}}^{\top}{\bm{X}}_{\mathcal{A}_{\lambda}\backslash\{j\}})^{-1}{\bm{X}}^{\top}_{\mathcal{A}_{\lambda}\backslash\{j\}}$.

Thus, under the null hypothesis $H_{0,j}^{\ast}:\beta^{\ast}_{j}=0$, (71) is equal to

[TABLE]

The equality holds because $({\bf I}_{n}-{\bf P}^{(\mathcal{A}_{\lambda}\backslash\{j\})}){\bm{X}}_{\mathcal{A}_{\lambda}\backslash\{j\}}{\bm{\beta}}^{\ast}_{\mathcal{A}_{\lambda}\backslash\{j\}}={\bm{0}}$.

We first show the asymptotic distribution of ${\bm{x}}_{j}^{\top}({\bm{I}}_{n}-{\bm{P}}^{(\mathcal{A}_{\lambda}\backslash\{j\})}){\bm{\epsilon}}$. Dividing it by its standard deviation $\sigma_{\bm{\epsilon}}\sqrt{{\bm{x}}_{j}^{\top}({\bf I}_{n}-{\bf P}^{(\mathcal{A}_{\lambda}\backslash\{j\})}){\bm{x}}_{j}}$, where $\sigma_{\bm{\epsilon}}$ is the error standard deviation,

[TABLE]

where ${\bf r}^{s\top}\equiv{\bm{x}}_{j}^{\top}({\bf I}_{n}-{\bf P}^{(\mathcal{A}_{\lambda}\backslash\{j\})})$. Now, we use the Lindeberg-Feller Central Limit Theorem to prove the asymptotic normality of (73). Similar to the proof of Theorem 2.3, we need to prove that the Lindeberg’s condition holds, i.e.,

[TABLE]

Given that $|r^{s}_{i}|\leq\|{\bf r}^{s}\|_{\infty}$, and that the $\epsilon_{i}$’s are identically distributed,

[TABLE]

Since $\|{\bf r}^{s}\|_{\infty}/\|{\bf r}^{s}\|_{2}\rightarrow 0$ by Condition (S), $\epsilon_{1}^{2}\mathbf{1}\left\{|\epsilon_{1}|\|{\bf r}^{s}\|_{\infty}/(\sigma_{\bm{\epsilon}}\|{\bf r}^{s}\|_{2})>\eta\right\}\rightarrow_{p}0$. Thus, because $\epsilon_{1}^{2}\geq\epsilon_{1}^{2}\mathbf{1}\left\{|\epsilon_{1}|\|{\bf r}^{s}\|_{\infty}/(\sigma_{\bm{\epsilon}}\|{\bf r}^{s}\|_{2})>\eta\right\}$ with probability one and $\mathbb{E}[\epsilon_{1}^{2}]=\sigma_{\epsilon}^{2}<\infty$, we use $\epsilon_{1}^{2}$ as the dominant random variable, and apply the Dominated Convergence Theorem,

[TABLE]

which in turn gives the Lindeberg’s condition.

Thus,

[TABLE]

We now prove the asymptotic unbiasedness of the naïve score test on ${\bm{\beta}}^{\ast}$. Dividing the second term in (72) by $\sigma_{\bm{\epsilon}}\sqrt{{\bm{x}}_{j}^{\top}({\bf I}_{n}-{\bf P}^{(\mathcal{A}_{\lambda}\backslash\{j\})}){\bm{x}}_{j}}$, we get

[TABLE]

By (M2∗),

[TABLE]

where the second equality holds under $H_{0,j}:\beta^{\ast}_{j}=0$.

Using, again, the fact that by Proposition 2.2, $\lim_{n\to\infty}\Pr\left[\mathcal{A}_{\lambda}=\hat{\mathcal{A}}_{\lambda}\right]=1$, we get

[TABLE]

∎

Acknowledgements

We thank the authors of Javanmard and Montanari, 2014a ; Ning and Liu, (2017) for providing code for their proposals. We are grateful to Joshua Loftus, Jonathan Taylor, Robert Tibshirani, and Ryan Tibshirani for helpful responses to our inquiries.