Simultaneous inference with SILM

SILM provides simultaneous inference for the high-dimensional linear model \(Y = X\beta + \epsilon\), based on the de-biased (de-sparsified) lasso and the bootstrap. It implements two closely related papers:

Zhang and Cheng (2017, JASA) bootstrap the linear part of the de-biased lasso with Gaussian multipliers: simultaneous confidence intervals (Sim.CI()), support recovery (SR()), testing for sparse signals (ST()) and stepdown multiple testing (Step()).
Dezeure, Bühlmann and Zhang (2017, TEST) bootstrap the entire de-sparsified lasso (boot.lasso.proj()): residual, wild (Gaussian or Mammen multipliers) or xyz-paired bootstrap, robust standard errors, Westfall–Young multiple testing, simultaneous confidence intervals and group tests.

The examples below use small dimensions so that the vignette builds quickly; the methods are intended for large \(p\).

library(SILM)
set.seed(1)
n <- 100
p <- 60
Sigma <- 0.9^abs(outer(1:p, 1:p, "-"))            # Toeplitz design
X <- matrix(rnorm(n * p), n, p) %*% chol(Sigma)
beta <- c(1.5, -1, 2, rep(0, p - 3))
Y <- as.vector(X %*% beta + rt(n, 4) / sqrt(2))

The model has no intercept. The data above have mean zero; for real data use center = TRUE (or centre the data beforehand).

Zhang and Cheng (2017)

Here \(p > n/2\), so the matrix \(\Theta\) (the inverse of the population Gram matrix) is estimated by the nodewise lasso. It depends only on the design and can be computed once:

Theta <- Theta.hat(X)

Simultaneous 95% confidence intervals for \(\beta_1, \dots, \beta_5\) (note that alpha is the confidence level in Sim.CI()):

ci <- Sim.CI(X, Y, set = 1:5, M = 500, alpha = 0.95, Theta = Theta)
ci$band.st
#>             [,1]       [,2]     [,3]       [,4]       [,5]
#> low.st 0.9853691 -1.2745464 1.154379 -0.5730351 -0.2408087
#> up.st  2.0708458 -0.0271207 2.430306  0.6172880  0.9304817

Support recovery and stepdown multiple testing at familywise level 5%:

SR(X, Y, Theta = Theta)
#> $`de-biased Lasso`
#> [1] 1 3
#> 
#> $`scaled Lasso`
#> [1]  1  3  5 23 30 32 45 50
Step(X, Y, M = 500, alpha = 0.05, Theta = Theta)
#> $`non-studentized test`
#> [1] 1 3
#> 
#> $`studentized test`
#> [1] 1 3

Testing whether any of the variables 4 to 60 has a non-zero coefficient, with a screening sub-sample of size \(0.3n\):

ST(X, Y, sub.size = 0.3, test.set = 4:p, M = 500)[c(3, 4)]
#> $`studentized test`
#> [1] 1.994765
#> 
#> $`studentized test`
#> [1] "fail to reject"

Dezeure, Bühlmann and Zhang (2017)

boot.lasso.proj() bootstraps the entire estimator. The paper recommends the robust standard error and, for heteroscedastic errors, the wild bootstrap. boot.shortcut = TRUE refits the bootstrap lasso at the original tuning parameter, which is much faster.

fit <- boot.lasso.proj(X, Y, B = 200, robust = TRUE, wild = TRUE,
                       boot.shortcut = TRUE, return.bootdist = TRUE)
fit
#> Bootstrapped de-sparsified lasso 
#> family gaussian; robust (divisor n - s) standard errors; wild (gaussian multipliers) bootstrap, B = 200 
#> Estimated noise level (sigmahat): 0.9907 
#> Significant at level 0.05 (adjusted p-values): 1, 3
#> Significant at level 0.01 (adjusted p-values): 1, 3
#> 
#> Coefficients with the smallest adjusted p-values (10 of 60):
#>    estimate     se  p.value p.adjusted
#> 1    1.5350 0.2337 0.004975   0.004975
#> 3    1.7520 0.2474 0.004975   0.004975
#> 2   -0.5757 0.2558 0.014930   0.731300
#> 10  -0.4095 0.1947 0.024880   0.845800
#> 50   0.3699 0.1961 0.094530   0.970100
#> 38  -0.3746 0.2067 0.064680   0.975100
#> 11   0.4186 0.2357 0.064680   0.980100
#> 16   0.3328 0.1998 0.004975   0.990000
#> 23  -0.3449 0.2009 0.204000   0.990000
#> 55   0.3913 0.2547 0.184100   0.995000

Individual and simultaneous confidence intervals:

confint(fit, parm = 1:5)
#>        lower      upper
#> 1  1.0921443  2.0412785
#> 2 -1.2644342 -0.1592080
#> 3  1.2597368  2.3796765
#> 4 -0.6406480  0.5256013
#> 5 -0.2916762  0.7422369
confint(fit, type = "simultaneous", group = 1:5)
#> Simultaneous 95% confidence intervals (group of 5 coefficients; max/min statistic; B = 200)
#>        lower     upper
#> 1  0.8243849 2.2048126
#> 2 -1.3529671 0.1576355
#> 3  0.9997451 2.4610950
#> 4 -0.7612134 0.8228169
#> 5 -0.5525840 1.0492636

Group tests of \(H_{0,G}: \beta_j = 0\) for all \(j \in G\):

groupTest(fit, list(signal = 1:3, noise = 4:p))
#>      signal       noise 
#> 0.004975124 0.830845771

Mammen multipliers and the xyz-paired bootstrap are selected with multiplier = "mammen" and boot.type = "xyz".

lasso.proj() computes the de-sparsified lasso with asymptotic Gaussian p-values and intervals (the “original” estimator of the paper):

lp <- lasso.proj(X, Y, robust = TRUE)
confint(lp, parm = 1:3)
#>       lower       upper
#> 1  1.079160  2.01112613
#> 2 -1.074198 -0.04933745
#> 3  1.253999  2.24667762

Two bootstraps

Sim.CI() bootstraps only the linear part \(\hat\Theta X^\top \epsilon / n\) of the de-biased lasso, which is fast and needs no refitting. boot.lasso.proj() refits the lasso for every bootstrap sample. Dezeure, Bühlmann and Zhang (2017) found that the latter often has better finite-sample coverage. It also handles heteroscedastic errors (with the wild bootstrap and robust standard errors). Their default nodewise tuning also differs: Sim.CI() uses the Z&Z rule (nodewise = "ZnZ"), while boot.lasso.proj() keeps hdi’s cross-validated lambda (do.ZnZ = FALSE).

Reproducing results of the archived versions

SILM 2.0.0 reproduces the archived SILM 1.0.0 and hdi 0.1-10 exactly (the same numbers under the same random seed):

In the development version, exact archived replication additionally requires scaled.lasso = "legacy" for the four core inference functions and whenever betainit = "scaled lasso" is used in either projection function. The new default normalizes the response and uses relative convergence to fix response unit sensitivity. The legacy setting deliberately retains that defect.

SR(), Sim.CI(), Step(), ST(): results obtained with SILM 1.0.0 and hdi 0.1-7 to 0.1-10 (March 2019 to July 2026) are the default (nodewise = "ZnZ", chosen by a pre-registered simulation study; see ?SR). Results obtained with hdi 0.1-6 (January to March 2019) need nodewise = "cv", the tuning described in Zhang and Cheng (2017). For ST(), legacy = TRUE also restores three behaviours of SILM 1.0.0 that were corrected (see ?ST).
lasso.proj(), boot.lasso.proj(): identical to hdi 0.1-10, except for robust = TRUE, where hdi’s results need robust.divisor = "n" (the default "n-s" is equation (5) of Dezeure, Bühlmann and Zhang (2017), chosen by a pre-registered simulation study; see ?lasso.proj), and lasso.proj(family = "binomial"), where they need legacy = TRUE.
Results computed with R < 3.6.0 used a different random sampling algorithm: call RNGkind(sample.kind = "Rounding") before set.seed().

The equivalence is verified continuously against the archived CRAN packages (see the validation/ directory of the source repository).

References

Zhang, X. and Cheng, G. (2017). Simultaneous inference for high-dimensional linear models. Journal of the American Statistical Association, 112, 757–768.

Dezeure, R., Bühlmann, P. and Zhang, C.-H. (2017). High-dimensional simultaneous inference with the bootstrap. TEST, 26, 685–719.

Dezeure, R., Bühlmann, P., Meier, L. and Meinshausen, N. (2015). High-dimensional inference: confidence intervals, p-values and R-software hdi. Statistical Science, 30, 533–558.