1 Introduction↩︎

Modern statistical problems are widely recognized as having become increasingly high-dimensional [1]. Examples such as genomics [2], neuroimaging [3], asset pricing [4], signal processing [5] and collaborative filtering¹[7] have led to settings in which the number of covariates exceeds the number of observations (\(p > n\)). LASSO regression [8] has become a widely applicable and successful technique for addressing such problems. Estimators are derived by optimizing a regularized likelihood function, typically incorporating an \(L_1\) penalty to induce sparsity. Because they are defined as the solution to a constrained optimization problem, LASSO estimators are inherently non-linear and do not admit an explicit, closed-form representation.

The non-linear and non-explicit nature of the estimators produced by LASSO estimates results in being unable to characterize the distribution the estimators it produces [9]. This creates a challenge inference as one is unable to preform accepted procedures for characterizing uncertainty via confidence intervals and p-values [10], [11]. As of 2014² approaches to quantify statistical significance in high-dimensional parameter estimation settings has been less understood. [12] and [13] introduced high-dimensional testing frameworks that utilize compatibility conditions to provide deterministic guarantees for significance \(\alpha\) and power \(1-\beta\) under restricted eigenvalue or compatibility conditions [14]. These results were later strengthened by [15], who established more refined lower bounds for the signal detection threshold. [16] and [17] consider inference for regression models with high-dimensional data. [18] and [19] studied resampling methods for hypotheses testing for these methods.

This report reviews and discusses the theoretical framework introduced by [1] in their Journal of Machine Learning Research paper titled “Confidence Intervals and Hypothesis Testing for High-Dimensional Regression”. Specifically, it examines the key lemmas and theorems underlying their approach to debiased inference for LASSO regression and presents a partial replication³ of the original simulation study and an analysis of the data used in the original data analysis. The analysis compares their proposed method with the multisample splitting method [20], [21] and the ridge-type projection estimator [12], as well as with alternative approaches implemented in the hdi R package [22]. Furthermore, the comparison is extended to include the desparsified LASSO estimator (also known as the LASSO projection estimator, used interchangably) [13], [23], [24], a closely related approach to debiased inference that was not formally evaluated (but referenced) in the original paper.

In Section 2, we outline the model setting and the de-biasing algorithm proposed by [1], followed by an overview of technical prerequisites defined in Section 2.1. Section 2.2 presents the key theoretical results required for formal inference with the debiased estimator. Section 2.3 provides an overview of the primary applications derived from the theorems, including the construction of confidence intervals, hypothesis tests, and extensions to simultaneous inference and settings involving non-Gaussian noise. Section 3 assesses the empirical performance of [1]’s debiased LASSO estimator through a partial replication of their simulation studies (Section 3.1) and a analysis of the riboflavin dataset [25]. By examining coverage probabilities, interval lengths, and Type I error rates, we compare the proposed method against competing approaches noted, but extend the analyses to include the LASSO projection estimator in order to highlight practical trade-offs. Section 4 concludes by combining the theoretical and empirical findings, discussing the practical trade-offs of debiased inference and its implications for high-dimensional statistical analysis in comparison with competing methods such as the LASSO projection estimator.

For reproducibility, the R code used to perform the simulations and data analysis is available at:
https://github.com/benyamindsmith/javanmard_montanari_2014_replication

2 Methodology↩︎

We work in a linear regression model setting where \(n\) i.i.d. pairs \((Y_1, X_1), (Y_2,X_2), \dots, (Y_n, X_n)\) are given with vectors \(X_i \in \mathbb{R}^p\) and response variables defined as \[Y_i= \langle\theta_0, X_i\rangle + W_i, \quad W_i \sim N(0, \sigma^2) \label{eq:mod951dim}\tag{1}\] Where \(\theta_0\in \mathbb{R}^p\) and \(\langle\cdot, \cdot\rangle\) is the standard scalar product in \(\mathbb{R}^p\). In matrix form, denote \(Y=(Y_1, \dots, Y_n)^\top\) and \(\mathbf{X} \in \mathbb{R}^{n\times p}\). This gives \[Y = \mathbf{X}\theta_0 + W, \quad W \sim N(0, \sigma^2 I_{n\times n}) \label{eq:mod95mat95form}\tag{2}\] The goal is to estimate the unknown (but fixed) parameter vector \(\theta_0\in \mathbb{R}^p\). In a classical regression setting where \(n \gg p\), the estimation method of choice is usually ordinary least squares, i.e. \(\hat{\theta}^{\text{OLS}}= (\mathbf{X}^\top\mathbf{X})^{-1}\mathbf{X}^\top Y\). By definition, \(\hat{\theta}^{\text{OLS}}\) is Gaussian with mean \(\theta_0\) and and covariance \(\sigma^2(\mathbf{X}^\top \mathbf{X})^{-1}\). This allows for typical inference techniques such as confidence intervals, hypothesis tests and construction of P-values.

In a high dimensional setting where the number of parameters exceeds the number of observations (\(p > n\)), the matrix \((X^\top X)\) is rank deficient and one is required to resort to biased estimators [14]. The LASSO [8], [26] has become a widely used method that produces sparse estimates by incorporating an \(\ell_1\) penalty into the objective function: \[\hat{\theta}^n(Y, \mathbf{X};\lambda) = \arg \underset{\theta \in \mathbb{R}^p}{\min} \left\{\frac{1}{2n}\|Y - \mathbf{X}\theta\|^2_2 + \lambda \|\theta\|_1\right\} \label{eq:lasso95eq}\tag{3}\]

The \(\ell_1\) penalty used in LASSO shrinks coefficients toward zero, thus introducing bias into the estimators. While this promotes sparsity, it also complicates statistical inference for the selected model as the resulting sampling distribution is difficult to characterize. Algorithm 1 shows the algorithm proposed by [1] to construct the debiased estimator from the LASSO, given by Equation ([eq:lasso95debiased]) as \(\hat{\theta}^{u}=\hat{\theta}^{\,n}(\lambda)+\frac{1}{n}M\mathbf{X}^{\top}\bigl(Y-\mathbf{X}\hat{\theta}^{\,n}(\lambda)\bigr)\). The basic intuition is that the subgradient of the \(\ell_1\) norm at the LASSO solution \(\hat{\theta}^n\) is \(\mathbf{X}^\top(Y - \mathbf{X}\hat{\theta}^n)/(n\lambda)\). By adding a term which is proportional to the subgradient, ([eq:lasso95debiased]) compensates the bias introduced by the \(\ell_1\) penalty in the LASSO.

It is claimed by [1] that the estimator ([eq:lasso95debiased]) guarantees near optimal confidence interval sizes and testing power. Section 2.1 presents preliminary definitions necessary for understanding the core theory presented in [1]. Section 2.2 presents the theorems which demonstrate that \(\hat{\theta}^{u}\) is Gaussian and is an asymptotically unbiased estimator for \(\theta_0\). Section 2.3 shows that the results from theorems in section 2.2 can be applied directly to derive confidence intervals (section 2.3.1), show near optimality of hypotheses tests generated (section 2.3.2),along with being generalized to derivation of simultaneous confidence intervals (Section 2.3.3) and the validity of the estimator \(\hat{\theta}^u\) in non-Gaussian noise settings (Section 2.3.4).

2.1 Prerequisites↩︎

In this section, we will present basic definitions which will be used later in this paper.

In words, \(\text{supp}(\theta_0)\) is the set of non-zero sparse parameters in the model (1 ). We further denote \(s_0 \equiv |S|\) to be the number of elements in \(\text{supp}(\theta_0)\).

The next two definitions provide the necessary foundation for the theorems in Section 2.2, which establish the unbiasedness of \(\hat{\theta}^{u}\).

Note that the minimum coherence can be computed efficiently since \(M \mapsto \mu_*(\mathbf{X}; M)\) is a convex function.

Assuming for the sake of simplicity that the columns of \(\mathbf{X}\) are normalized so that \(\|\mathbf{X}e_i\|_2 = \sqrt{n}\), a small value of the coherence parameter \(\mu_*(\mathbf{X}; I)\) means that the columns of \(\mathbf{X}\) are roughly orthogonal. So \(\mu_*(\mathbf{X}; I) =0\) if and only if \(\mathbf{X}/\sqrt{n}\) is column orthogonal. However, the generalized coherence parameter \(\mu_*(\mathbf{X}; M)\) can be much smaller (i.e. it can remain small, potentially as small as order \(O(\sqrt{(\log p)/n})\), even when the columns of \(\mathbf{X}\) are significantly correlated).

Armed with these definitions, we now present the central theorems that establish the asymptotic normality of the debiased estimator \(\hat{\theta}^u\).

2.2 A Debiased Estimator for \(\theta_0\)↩︎

Having established the generalized coherence \(\mu_*\), we now characterize the distributional properties of the debiased estimator. The following theorems provide the conditions under which \(\hat{\theta}^u\) is asymptotically normal, which will allow for classical inferential procedures such as confidence intervals and hypotheses testing. In order to clarify the distributional properties of \(\hat{\theta}^u\), [1] consider a general debiasing procedure that makes use of an arbitrary debiasing matrix \(M \in \mathbb{R}^{p\times p}\) and define the generalized debiased estimator \(\hat{\theta}^*\) as: \[\hat{\theta}^*(Y, \mathbf{X}; M, \lambda) = \hat{\theta}^n(\lambda) + \frac{1}{n}M\mathbf{X}^\top (Y-\mathbf{X}\hat{\theta}^n(\lambda)) \label{eq:gen95debiased95theta}\tag{4}\]

Theorem 2 decomposes the estimation error of the generalized debiased estimator \((\hat{\theta}^*- \theta_0)\) into a zero mean Gaussian term \(Z/\sqrt{n}\) and a bias term \(\Delta/\sqrt{n}\) whose maximum entry is bounded in probability per Equation ?? .

The next theorem establishes that for a natural probabilistic model of the design matrix \(\mathbf{X}\), the (constant) bound of the compatibility constant \(\phi_0\) and the generalized coherence parameter \(\mu_*(\mathbf{X}; M)\) can be bounded with probability converging rapidly to one as \(n, p \to \infty\).

Theorem 3 shows for a natural probabilistic model of the design matrix \(\mathbf{X}\in \mathbb{R}^{n\times p}\), both \(\phi_0\) and \(\mu_*(\mathbf{X}; M)\) can be bounded with probability that both converge rapidly to 1 as \(n,p \to \infty\). Specifically, it guarantees that for sufficiently large \(n\), the design matrix remains numerically stable enough to recover the true parameters while maintaining the low coherence required to filter out the shrinkage bias inherent in the LASSO. The proof for this theorem is omitted as it is a more nuanced technicality. For the proof, see the appendix of [1].

Theorems 2 and 3 establish the necessary probabilistic guarantees for recovering asymptotically unbiased and Gaussian estimates of \(\theta_0\). The following theorem builds on these results to characterize the limiting distribution of the debiased estimator \(\hat{\theta}^u\) produced by Algorithm 1.

Assuming \(\sigma\), \(C_{\min}\) of order 1, Theorem 4 establishes that for random designs, the maximum size of the bias term \(\Delta_i\) over \(i \in [p]\) is \[\|\Delta\|_{\infty} = O\left(\frac{s_0\log p}{\sqrt{n}}\right)\] On the other hand, the ‘noise term’ \(Z_i\) is roughly of order \(\sqrt{[M \hat{\Sigma}M^\top]_{ii}}\)

We will now present the proof for Theorem 4:

Theorem 4 shows that asymptotic normality and error bounds of the generalized debiased estimator \(\hat{\theta}^*\) can be applied to the debiased estimtor proposed by [1] in Algorithm 1 \(\hat{\theta}^u\) (See Equation (4 )).

2.3 Statistical Inference↩︎

A direct application of Theorem 4 is to derive confidence intervals and statistical hypothesis tests for high dimensional models. The next Lemma will show how this is possible. Note, the sparisty assumption is assumed to be \(s_0 = o(\sqrt{n}/\log p)\).

While the standard LASSO is an excellent tool for variable selection, its inherent bias makes it unsuitable for traditional statistical inference. Lemma 1 provides the necessary theoretical bridge, transforming \(\hat{\theta}^u\) into a standardized statistic that converges to a standard normal distribution. This in turn shows that it is permissible to perform ‘OLS-style’ inference in settings where \(p \gg n\)

Lemma 1 is the bridge which facilitates standard inferential techniques such as confidence intervals and hypothesis testing for \(\hat{\theta}^u\). Furthermore, this result is generalizable to simultaneous inference and remains robust under non-Gaussian noise distributions. These results are presented in subsections 2.3.1,2.3.2,2.3.3 and 2.3.4 without proof for brevity.

2.3.1 Confidence Intervals↩︎

From Lemma 1 it is quite straightforward to construct asymptotically valid confidence intervals. For \(i \in [p]\) and significance level \(\alpha \in (0,1)\), the confidence interval of the point estimate of the debaised LASSO estimator \(\hat{\theta}^u\) is defined as: \[\label{eq:ci95formula} \begin{align} J_i(\alpha) &\equiv \left[ \hat{\theta}_i^u - \delta(\alpha, n), \hat{\theta}_i^u + \delta(\alpha, n) \right], \\ \delta(\alpha, n) &\equiv \Phi^{-1}(1 - \alpha/2) \frac{\hat{\sigma}}{\sqrt{n}} [M \hat{\Sigma} M^\top]_{i,i}^{1/2} \end{align}\tag{6}\] It can also be shown that the confidence interval is asymptotically valid, namely \[\lim_{n \to \infty} \mathbb{P}\left(\theta_{0,i} \in J_i(\alpha) \right) = 1-\alpha\]

2.3.2 Hypothesis Testing↩︎

An important advantage of sparse linear regression methods such as the LASSO is their ability to provide a parsimonious explanation of the data, relying on only a small number of covariates. This can be simply understood as selecting the number of non-zero covariates \(\hat{\theta}^n_i \neq 0\) for all \(i \subseteq [p]\). However, the covariates estimated by the LASSO \(\hat{\theta}^n_i\) do not provide a measure of the statistical significance for the finding that the coefficient is non-zero. In particular, we are interested in testing an individual null hypotheses \(H_{0,i}: \theta_{0,i} = 0\) verses the alternative \(H_{A,i}: \theta_{\theta,i} \neq 0\). Using the debiased estimator \(\hat{\theta}^u\) derived in Algorithm 1, we construct the p-value \(P_i\) for the test \(H_{0,i}\) to be: \[P_{i} = 2 \left(1 - \Phi\left(\frac{\sqrt{n}|\hat{\theta}^{u}_i |}{\hat{\sigma}[M \hat{\Sigma M^\top}]^{1/2}_{i,i}}\right) \right)\] The decision rule is based on the p-value \(P_i\): \[\widehat{T}_{i, \mathbf{X}}(y) = \begin{cases} 1 & \text{if } P_i \le \alpha \quad (\text{reject } H_{0,i}), \\ 0 & \text{otherwise} \quad (\text{accept } H_{0,i}) \end{cases}\] Where \(\alpha\) is the fixed target Type I error probability. We measure the quality of the test \(\widehat{T}_{i, \mathbf{X}}(y)\) in terms of its significance level \(\alpha_i\) and statistical power \(1- \beta_i\) where \(\beta_i\) is the probability of type II error. It can be shown that this test achieves power \(1 - \beta =\alpha\). From a formal minimax perspective, given a family of tests \(T_{i, \mathbf{X}}:\mathbb{R}^n \to \{ 0, 1\}\) indexed by \(i \in [p]\), \(\mathbf{X} \in \mathbb{R}^{n\times p}\), we defined for \(\gamma >0\) a lower bound on the non-zero entries of \(\theta_{0,i}\) \[\begin{align} \alpha_{i, n}(T) &\equiv \sup \left\{\mathbb{P}_{\theta_0}(T_{i,\mathbf{X}}(y)= 1):\theta_0 \in \mathbb{R}^p, \|\theta_0\|_0 \le s_0(n), \theta_{0,i} = 0 \right\}\\ \beta_{i, n}(T; \gamma) & \equiv\sup \left\{\mathbb{P}_{\theta_0}(T_{i,\mathbf{X}}(y)= 0):\theta_0 \in \mathbb{R}^p, \|\theta_0\|_0 \le s_0(n), \theta_{0,i} \ge \gamma \right\} \end{align}\]

2.3.3 Generalizing to Simultaneous Confidence Intervals↩︎

In many situations, it is necessary to perform statistical inference on more than one of the parameters simultaneously. As an example, we might be interested in performing inference about \(\theta_{0,R}\equiv (\theta_{0,i})_{i \in R}\) for some set \(R \subseteq [p]\). The simplest generalization of the method method is when the number of indices \(|R|\) stay finite as \(n, p \to \infty\). The following Lemma is a generalization of Lemma 1.

Mirroring the univariate inferential capabilities provided by Lemma 1, Lemma 2 provides the necessary theoretical results to move from marginal to simultaneous confidence regions. We can then define the valid confidence region \(J_R(\alpha) \subseteq \mathbb{R}^k\) as \[J_{R}(\alpha) \equiv \hat{\theta}_R^u + (Q^{(n)}_{R,R})^{1/2}\mathcal{C}_{k,\alpha}\] Such that \[\lim_{n \to \infty} \mathbb{P}\left(\theta_{0, R} \in J_R(\alpha)\right)= 1-\alpha\] Where \(\mathcal{C}_{k, \alpha}\subseteq \mathbb{R}^k\) is any Borel set such that \(\int_{\mathcal{C}_{k, \alpha}} \phi(x) dx \ge 1-\alpha\) and \(\phi(x) = \frac{1}{(2\pi)^{k/2}}\exp \left(- \frac{\|x\|^2}{2}\right)\). We can further construct a testing procedure that controls the family wise error rate (FWER) based on Bonferroni. \[\widehat{T}^F_{i, \mathbf{X}}(y) = \begin{cases} 1 \quad \text{if }P_i \le \alpha/p \quad (\text{reject } H_{0,i})\\ 0 \quad \text{otherwise} \quad (\text{accept } H_{0,i}) \end{cases}\] It then follows from the Bonferroni inequality that the asymptotic family-wise error rate is bounded above by \(\alpha\): \[\limsup_{n\to\infty} \text{FWER}(\widehat{T}^F,n) \le \alpha\]

2.3.4 Non-Gaussian Noise↩︎

The inferential guarantees for the simultaneous confidence region \(J_R(\alpha)\), the individual hypothesis tests, and the FWER-controlling procedure \(\widehat{T}^F\) naturally extend to non-Gaussian noise, provided the independent errors \(W_i\) have zero mean, variance \(\sigma^2\), and a bounded moment \(\mathbb{E}[|W_i|^{2+a}] \le C\sigma^{2+a}\) for some \(a > 0\). To ensure the asymptotic normality required for valid confidence intervals and \(p\)-values, the scaled summands of the error, denoted as \(\xi_j=m_i^\top X_j W_j/(\sigma[m_i^\top \hat{\Sigma}m_i]^{1/2})\), must satisfy the Lindeberg condition. Namely, for every \(\varepsilon > 0\), almost surely: \[\lim_{n \to \infty} \frac{1}{n} \sum_{j=1}^n \mathbb{E}\left(\xi_j^2 \mathbb{I}_{\{|\xi_j| > \varepsilon\sqrt{n}\}} \mid X\right) = 0\] Then \(\sum_{j=1}^n \xi_k/\sqrt{n} \overset{d}\to N(0,1)\). This condition is satisfied by introducing an \(\ell_\infty\) constraint to the debiasing optimization problem: \(\|Xm_i\|_\infty \le n^\beta\) for a fixed \(\beta \in (1/4, 1/2)\). By bounding the maximum influence of any single observation, the Lindeberg condition holds, preserving the asymptotic Gaussianity that guarantees the validity of the intervals, the tests, and the FWER upper bound \(\limsup_{n\to\infty} \text{FWER}(\widehat{T}^F,n) \le \alpha\).

To construct the debiasing matrix \(M = (m_1, \dots, m_p)^T\) under the non-Gaussian regime, each row vector \(m_i \in \mathbb{R}^p\) is obtained by solving the following convex optimization problem: \[\begin{align} & \underset{m \in \mathbb{R}^p}{\text{minimize}} & & m^T \hat{\Sigma} m \\ & \text{subject to} & & \|\hat{\Sigma}m - e_i\|_\infty \le \gamma, \\ & & & \|Xm\|_\infty \le n^\beta \quad \text{ for arbitrary fixed } 1/4 <\beta<1/2 \end{align}\] Consequently, under the assumptions of Theorem 4, the resulting p-values remain asymptotically valid in non-Gaussian settings; a formal statement and proof are omitted for brevity.

3 Empirical Analysis↩︎

To evaluate the theoretical results of [1], we replicate six of the original simulation scenarios (selected due to computational constraints) and the real-data application on the riboflavin (vitamin \(B_2\)) dataset. Building on this baseline, we compare their method with the Lasso projection estimator [13], [23], [24] (also known as the LASSO projection estimator), which was cited but not implemented in the original study. As the original work only benchmarked against multisample splitting and a ridge projection estimator, this additional comparison provides a more comprehensive assessment of its practical performance.

All computations were performed on a first-generation Intel NUC laptop equipped with a 12th Gen Intel Core i7-1260P processor. The replication study was conducted using R version 4.2.2, leveraging the parallel, doParallel, and foreach packages to enable parallel processing [27]–[29]. The implementation of the debiased LASSO proposed by [1] was sourced from the authors’ repository (https://web.stanford.edu/ montanar/sslasso/). Furthermore, routines for multisample splitting, the ridge projection estimator, and the LASSO projection estimator, alongside the riboflavin dataset, were accessed via the hdi R package [22].

3.1 Simulation Study↩︎

We consider the linear model (2 ) where rows of the design matrix \(\mathbf{X}\) are fixed i.i.d. realizations from the \(N(0, \Sigma)\), where \(\Sigma \in \mathbb{R}^{p\times p}\) is a circulant symmetric matrix with entries \(\Sigma_{jk}\) given as for \(j \le k\) \[\label{eq:circulant95mat} \Sigma_{jk} = \begin{cases} 1 &\text{if } k= j\\ 0.1 & \text{if } k \in\{j+1, \dots, j+5 \}\\ & \text{or } k \in \{j+p-5, \dots, j+p-1\}\\ 0 & \text{for all other } j\le k. \end{cases}\tag{7}\]

Figure 2 illustrates an example of this matrix structure for \(p=600\). To generate the true regression coefficients, we select a uniformly random support set \(S \subseteq [p]\) of size \(|S|=s_0\). The non-zero coefficients are set to \(\theta_{0,i}=b\) for \(i \in S\), and \(\theta_{0,i}=0\) otherwise. Measurement errors are drawn as \(W_i \sim N(0,1)\) for \(i \in [n]\). We evaluate several configurations of \((n, p, s_0, b)\), and for each configuration, the results are averaged over 20 independent realizations of the model with fixed design and regression coefficients. In contrast to [1], who explicitly tune a regularization parameter and model optimization constraint \(\mu\), we rely on the default configurations for all methods to better reflect how these implementations are typically used in practice.

For each configuration, 20 independent realizations of measurement noise and for each parameter \(\theta_{0,i}\) the average length of the corresponding confidence interval is computed. Denoted by \(\text{Avglength}(J_i(\alpha))\) where \(J_i(\alpha)\) is the estimated confidence interval for \(\hat{\theta}_i^u\) defined in Equation (6 ) and the average is taken over all realizations. The following metric considered is \[\ell \equiv p^{-1} \sum_{i \in [p]} \text{Avglength}(J_i(\alpha))\] Along with the average length of intervals for the active and inactive parameters as: \[\ell_S \equiv s_0^{-1} \sum_{i \in S}\text{Avglength}(J_i(\alpha)), \quad \ell_S^c \equiv (p-s_0)^{-1}\sum_{i \in S^c}\text{Avglength}(J_i(\alpha))\] Similarly, the average coverage for the individual parameters is considered. i.e. \[\widehat{\text{Cov}} = p^{-1}\sum_{i \in [p]}\widehat{\mathbb{P}}[ \theta_{0,i}\in J_i(\alpha)],\] \[\widehat{\text{Cov}}_S = s_0^{-1}\sum_{i \in S}\widehat{\mathbb{P}}[ \theta_{0,i}\in J_i(\alpha)],\] \[\widehat{\text{Cov}}_{S^c} = s_0^{-1}\sum_{i \in S^c}\widehat{\mathbb{P}}[ \theta_{0,i}\in J_i(\alpha)]\] Where probabilities (\(\widehat{\mathbb{P}}\)) are calculated from 20 realizations per configuration. Tables 1–[tab:tab4] compare interval coverage, false positive (FP) rates, and statistical power (TP) between this replication and [1]. Tables 2 and [tab:tab4] show the relevant results listed in [1].

As detailed in Tables 1 and 2, the replication produces tighter confidence intervals (e.g., \(\ell \approx 0.15\) versus \(\ell \approx 0.18\)) and achieves a more precise alignment with the 95% nominal coverage compared to the original results, which tended toward over-coverage. When evaluating error rates (Tables [tab:tab3] and [tab:tab4]), both studies successfully control FP rates near the nominal \(\alpha = 0.05\) level. However, the replication achieves substantially higher statistical power (TP) across multiple methods in moderately sparse settings, most notably maintaining full power for multisample-splitting where the original study observed a significant drop.

Crucially, Table [tab:tab3] highlights the performance of the LASSO projection estimator, which was not considered in the original study despite being cited. In direct comparison to the replicated Javanmard and Montanari method, the LASSO projection demonstrates superior sensitivity in low-signal regimes. For instance, under the \((1000, 600, 30, 0.1)\) configuration, the LASSO projection yields a TP of 0.8217 compared to the Javanmard and Montanari method’s 0.7017, while still rigorously controlling the FP rate (0.0479 versus 0.0363). These findings indicate that the LASSO projecton estimator provides a more powerful framework for high-dimensional inference than the original Javanmard and Montanari procedure, improving signal detection without inflating Type I error.

3.2 Data Analysis↩︎

To demonstrate our approach on empirical data, we analyze the riboflavin production dataset previously examined by [1]. Originally introduced by [25], the dataset comprises \(n=71\) samples and \(p=4,088\) covariates representing individual gene expression levels. The continuous response variable measures the log-transformed riboflavin production rate, which is modelled linearly against the log-transformed gene expression covariates.

Figure 3 shows that Multisample-splitting proved overly conservative with estimated confidence intervals which covered the whole range of the covariate distribution, whereas the debiased projection estimators and the [1] approach successfully identified globally significant genes such as ARGF_at and LYSC_at. The Javanmard and Montanari method exhibited the highest precision , consistently yielding the narrowest confidence intervals and isolating signal where the Ridge and Lasso projection estimators were less sensitive. These results highlight the substantial precision gains achievable through the optimized de-biasing procedure proposed by [1].

4 Discussion & Conclusion↩︎

The transition of high-dimensional regression from pure variable selection to formal statistical inference marks a critical methodological advancement. Because the non-linear \(\ell_1\) penalty of the standard LASSO introduces inherent shrinkage bias, classical uncertainty quantification via confidence intervals and p-values was precluded. The framework proposed by [1] successfully bridges this gap by constructing an optimal de-biasing matrix that corrects the LASSO estimate while explicitly minimizing asymptotic variance. The resulting theoretical guarantees for the estimator \(\hat{\theta}^u\) related to its asymptotic normality, resilience to non-Gaussian noise, and adaptability for simultaneous testing provide a rigorous mathematical foundation for classical statistical inference in high-dimensional setting when \(p \gg n\), under suitable sparsity assumptions.

Empirical evaluations across replications of the simulation study and the riboflavin dataset reveal a nuanced narrative regarding the practical performance of these estimators, highlighting a critical divergence between idealized settings and severely sample-starved, real-world applications. In the simulation study, multisample-splitting [20] maintained strong statistical power, and the LASSO projection estimator [13] outperformed the [1] method by demonstrating superior sensitivity in low-signal regimes. However, this dynamic reversed in a real high-dimensional genomic data setting. Furthermore, the Javanmard and Montanari estimator achieved the highest precision and successfully isolated empirical signals where the LASSO projection lacked sensitivity. This discrepancy suggests that while standard projection estimators excel under well-behaved, circulant covariance structures, the explicit convex optimization of the de-biasing matrix \(M\) employed by [1] adapts far more robustly to the complex correlation networks inherent in real biological data. While modern projection methods offer robust power in idealized scenarios, optimized de-biasing matrices produced by [1] prove to be dominant for maximizing signal detection and enabling inference in noisy, high-dimensional settings.

To provide a more equitable evaluation of these estimators, it is recommended that future simulation studies employ design matrices that intentionally weaken or violate the restricted eigenvalue and compatibility conditions. The circulant covariance structure used in the simulation study features minimal off-diagonal correlation (as shown in Figure 2), creating an idealized setting where the LASSO projection estimator excels. Furthermore, these settings fail to capture the complex correlation structures inherent in real-world biological data. By benchmarking against highly correlated data generating processes, future research can more accurately demonstrate the practical necessity of the debaised estimator proposed by [1].

References↩︎