← All posts

Confidence intervals in high dimensions, in one R package: SIHR

The R Journal · with Prabrisha Rakshit (equal contribution), T. Tony Cai, Zijian Guo · arXiv:2109.03365 · package

The problem: a point estimate is not a finding

When a regression has more variables than observations, penalized estimators like the Lasso are the workhorse. Consider the penalized maximum likelihood estimator for a generalized linear model, \(\mathbb E(y_i \mid X_{i\cdot}) = f(X_{i\cdot}^\top\beta)\) with \(f(z) = z\) (linear) or \(f(z) = e^z/(1+e^z)\) (logistic):

\[ \widehat\beta \;=\; \arg\min_{\beta}\ \ell(\beta) \;+\; \lambda_0 \sum_{j=2}^{p} \frac{\|X_{\cdot j}\|_2}{\sqrt n}\,|\beta_j|, \qquad \lambda_0 \asymp \sqrt{\log p / n}. \]

This attains optimal estimation rates and selects variables well. What it cannot do is support inference. The penalty shrinks coefficients toward zero, and that shrinkage is a bias whose magnitude is comparable to — or larger than — the estimator's standard error. Report \(\widehat\beta_j\) with a naive standard error and the resulting interval is centred in the wrong place. Its coverage does not converge to the nominal level; it converges to something else entirely. In genomics, economics, and biostatistics — the fields that live in high dimensions — a point estimate without honest uncertainty is often not enough to make a claim.

The fix: debias, then infer

Take a target like \(x_{\mathrm{new}}^\top\beta\), a prediction at a new covariate value (set \(x_{\mathrm{new}} = e_j\) and you get the individual coefficient \(\beta_j\)). The plug-in \(x_{\mathrm{new}}^\top\widehat\beta\) carries the bias; the correction adds back an estimate of it, weighted by a projection direction \(\widehat u\):

\[ \widehat{x_{\mathrm{new}}^\top\beta} \;=\; x_{\mathrm{new}}^\top\widehat\beta \;+\; \widehat u^\top\, \frac1n\sum_{i=1}^{n} X_{i\cdot}\big(y_i - X_{i\cdot}^\top\widehat\beta\big). \]

Why does this help? Because the error of the corrected estimator decomposes exactly into a term we understand and a term we can control:

\[ \widehat{x_{\mathrm{new}}^\top\beta} - x_{\mathrm{new}}^\top\beta \;=\; \underbrace{\widehat u^\top \tfrac1n\textstyle\sum_i X_{i\cdot}\,\epsilon_i}_{\text{asymptotically normal}} \;+\; \underbrace{(\widehat\Sigma\widehat u - x_{\mathrm{new}})^\top(\beta - \widehat\beta)}_{\text{remaining bias}} . \]

The second term is a product of two small things, and it is small provided \(\widehat\Sigma\widehat u\) is close to \(x_{\mathrm{new}}\). That is precisely what \(\widehat u\) is chosen to arrange, while keeping the variance \(u^\top\widehat\Sigma u\) of the first term as small as possible:

\[ \widehat u = \arg\min_{u}\ u^\top\widehat\Sigma u \quad\text{s.t.}\quad \underbrace{\|\widehat\Sigma u - x_{\mathrm{new}}\|_\infty \le \|x_{\mathrm{new}}\|_2\,\lambda}_{\text{controls the bias}}, \quad \underbrace{\big|x_{\mathrm{new}}^\top\widehat\Sigma u - \|x_{\mathrm{new}}\|_2^2\big| \le \|x_{\mathrm{new}}\|_2^2\,\lambda}_{\text{makes the normal term dominate}} . \]

The first constraint is the familiar one. The second is the crucial one: the paper describes it as what ensures asymptotic normality for any loading \(x_{\mathrm{new}}\), rather than only for the sparse loadings (like \(e_j\)) that coordinate-wise debiasing was built for. With both in force, \(\big(\widehat{x_{\mathrm{new}}^\top\beta} \pm z_{\alpha/2}\sqrt{\widehat V}\big)\) is an honest interval, with \(\widehat V = \tfrac{\widehat\sigma^2}{n}\widehat u^\top\widehat\Sigma\widehat u\).

The same machinery covers logistic regression once a weight \(\omega(\cdot)\) is inserted, so that a single implementation serves all supported models:

\[ \widehat{x_{\mathrm{new}}^\top \beta} = x_{\mathrm{new}}^\top \widehat\beta + \widehat u^\top \frac1n \sum_{i=1}^n \omega\big(X_{i\cdot}^\top\widehat\beta\big)\big(y_i - f(X_{i\cdot}^\top\widehat\beta)\big) X_{i\cdot}. \]

The two logistic options implement different theoretical routes: logistic uses the linearization weighting, logistic_alter the link-specific one. For a case probability \(f(x_{\mathrm{new}}^\top\beta)\), apply \(f\) to the endpoints of the interval for \(x_{\mathrm{new}}^\top\beta\).

One engine, five targets

Here is the design idea that makes the package more than a wrapper. Many quantities of scientific interest are not linear in \(\beta\), yet their plug-in bias is. Take a quadratic functional \(Q_A = \beta_G^\top A \beta_G\), the natural measure of how much signal a group \(G\) of variables carries. Its plug-in bias is dominated by \(2\,\widehat\beta_G^\top A(\widehat\beta_G - \beta_G)\) — which is a linear functional with loading \(x_{\mathrm{new}} = (\widehat\beta_G^\top A, \mathbf 0^\top)^\top\). So you debias once more, with the machinery you already have:

\[ \widehat Q_A \;=\; \widehat\beta_G^\top A \widehat\beta_G \;+\; 2\,\widehat u_A^\top \Big[\frac1n\sum_i \omega\big(X_{i\cdot}^\top\widehat\beta\big)\big(y_i - f(X_{i\cdot}^\top\widehat\beta)\big)X_{i\cdot}\Big], \]

truncated at zero. The same trick — recognize the leading bias as a linear functional, then reuse LF() — generates treatment effects, two-sample inner products, and distances between regression vectors. The variance for the quadratic targets gains a term \(\tau/n\) that upper-bounds the discarded second-order piece \((\widehat\beta_G - \beta_G)^\top A(\widehat\beta_G - \beta_G)\); it also guards against super-efficiency when \(Q_A \approx 0\), which is exactly the null you usually want to test.

Companion methods ci() and summary() return the interval, and the plug-in estimate, bias-corrected estimate, standard error, z-value and p-value respectively.

In practice the whole thing is two lines. The columns of loading.mat are the loadings \(x_{\mathrm{new}}\) you care about, evaluated in a single call:

# n = 100, p = 120; two loadings evaluated at once
Est = LF(X, y, loading.mat, model = 'linear')
ci(Est)
summary(Est)

Does it work?

A word on what evidence exists. The paper is a software article: it does not carry a simulation study reporting coverage percentages or timing benchmarks, and its comparisons with hdi, DoubleML and selectiveInference are argued conceptually rather than measured. What it does provide is six worked examples with known ground truth, and two real analyses. The examples make the point cleanly — in every one, the plug-in is visibly shrunk toward zero and the debiased interval covers the truth.

Every interval covers its truth. For QF, InnProd and Dist the CI shown uses \(\tau = 0.25\); larger \(\tau\) widens it (for Dist, \(\tau = 0.5\) already truncates the lower endpoint at 0 and the effect stops being significant at the 5% level — \(p = 0.055\) versus \(p = 0.033\)). Note that Dist's plug-in is biased away from zero, as a squared quantity should be.

The two applications show the interface doing real work. In motif regression — predicting transcription-factor binding intensity from \(p = 666\) candidate motifs over \(n = 2{,}587\) genes — one call to LF() with loading.mat = diag(p) returns all 666 intervals; 25 lie entirely above zero and 23 entirely below, isolating 48 significant motifs. In a fasting-glucose GWAS on 1,269 heterogeneous-stock mice, the response is dichotomized at 11.1 mmol/L and the markers pruned from 10,346 to 2,341 (pairwise \(|\text{correlation}| < 0.75\)), with gender and age as baseline covariates; a logistic LF() call yields 13 intervals entirely above zero and 16 entirely below.

Why it matters

Debiased inference had strong theory and a high implementation barrier — the corrections are target-specific, the projection direction needs its own tuning, and the details live in a decade of separate papers. Most applied users never got past that. SIHR collapses the barrier: a handful of functions, a consistent interface, valid intervals for the functionals people actually want, and both linear and logistic models. The broader lesson is worth stating plainly — the debiasing step is a reusable primitive, not a one-off trick for individual coefficients. Once you see that a quadratic functional's leading bias is itself a linear functional, the same engine covers group signal, treatment effects, and cross-population comparisons.

Caveats to keep in view. Validity rests on the usual high-dimensional assumptions: a sparse \(\beta\), a correctly specified link, and tuning at \(\lambda \asymp \sqrt{\log p/n}\). The intervals are mildly conservative by design — rescale = 1.1 inflates standard errors to absorb finite-sample bias. Quadratic and two-sample targets depend on the tuning parameter \(\tau\), and the package deliberately reports a vector of intervals across \(\tau\) rather than hiding the choice; treat it as a conservative variance inflation, and note that conclusions can change with it, as Dist above shows. Logistic inference relies on prob.filter to drop observations with fitted probabilities near 0 or 1, and QF(split = TRUE) needs a larger \(n\). Finally, the comparisons against other packages are conceptual: the paper reports no head-to-head coverage or timing table.

© 2026 Zhenyu (Zach) Wang 王振宇