Epsilon Conundrum

What happens if we apply the bootstrap directly to the lasso?

A Proposed Fix to the Epsilon Conundrum

The lasso can be formulated as a Bayesian regression model by setting an appropriate prior (Tibshirani 1996; Park and Casella 2008). Specifically for the lasso, a Laplace distribution:

\begin{align*}p(\boldsymbol{\beta}) = \prod_{j = 1}^{p} \frac{\gamma}{2}\exp(-\gamma \left\lvert\beta_j\right\rvert), \gamma > 0.\end{align*}

With this prior, the lasso estimate, \widehat{\boldsymbol{\beta}}(\lambda), is the posterior mode for \boldsymbol{\beta}
when \gamma = n\lambda/\sigma^2.

To address the Epsilon Conundrum, we propose drawing from the full conditional posterior, p(\beta_j | \boldsymbol{\beta}_{-j}), whenever the bootstrap draw is 0.

We call this the Hybrid Bootstrap.

Epsilon Conundrum: The Hybrid Bootstrap

Benefits of Biased Estimators

The addition of a penalty has desirable implications for the resulting lasso estimators:

Sparsity
Can force some coefficients to exactly zero, leading to sparse, more interpretable models

Bias-Variance Trade-off
Controlled increase in bias often results in a significant reduction in variance, improving prediction

Overfitting Prevention
Shrinking coefficient estimates reduces model complexity, mitigating over fitting

Improved Stability in Presence of Multicollinearity
Stabilizes coefficient estimates when predictors are highly correlated

Traditional Frequentist Intervals

The statistical community seems to instinctively answer NO.

Popular methods for constructing intervals using lasso focus on debiasing in one way or another:

• The de-sparsified lasso focuses on debiasing the original point estimates from a lasso fit (Zhang and Zhang 2014; Javanmard and Montanari 2014).
• Selective Inference accounts for the uncertainty in model selection by conditioning on the selected model (Lee et al. 2016).

The goal of debiasing is to obtain nominal coverage for all values of \beta.

Bayesian Intervals - Ridge Regression Example

Average coverage

For a high-dimensional problem, instead of integrating with respect to the prior a more interesting quantity may be the average with respect to the distribution of \theta values present.

Does \tfrac{1}{p} \sum_{j=1}^p \textrm{Cover}(\theta_j) \approx \int \textrm{Cover}(\theta)p(\theta) \, d\theta ?

Ideal Conditions - Hybrid Bootstrap

Adding Autoregressive Correlation

Alternative Distributions of \beta

Overview

• Compare Hybrid, Selective Inference, and Bootstrapped De-sparsified Lasso (BDSL)
• Here, all methods are carried out at \lambda_{CV} from the original data set

Scheetz2006

• Scheetz et al. (2006) measured RNA levels from the eyes of 120 rats
• A gene with known association to the disease of interest (BBS) is used as the outcome
• The goal is to try to determine other genes whose expression is associated with BBS
• 120 observations
• 18975 features

Results

Details

• 5 Hybrid intervals do not contain point estimate, 4035 BDSL intervals do not, and all 66 Selective Inference intervals do not
• Hybrid identifies 3 significant genes, BDSL identifies 1927, Selective Inference “identifies” all 66 as significant
• Selective Inference took about .3 seconds, Hybrid took about 4.5 minutes, and BDSL took about 7 hours

Conclusion

If you chose a biased estimation method, is it not reasonable that the resulting intervals should also reflect that bias?

We introduced,

1. a method for interval construction that fixes the Epsilon Conundrum, and
2. an alternative perspective that allows for biased intervals.

Together these proposals lead to intervals which clearly reflect the assumptions that went into fitting the original lasso model.

The Hybrid bootstrap is available in the R package ncvreg.

Coverage Theorem

Proof. By definition, a 100(1-\alpha)\% credible region for \theta_j is any set \mathcal{A}_j(\mathbf{y}) such that \int I\{\theta_j \in \mathcal{A}_j(\mathbf{y})\} p(\boldsymbol{\theta}|\mathbf{y})\,d\boldsymbol{\theta}= 1 - \alpha. The coverage probability is defined as \int I\{\theta_j \in \mathcal{A}_j(\mathbf{y})\} p(\mathbf{y}| \boldsymbol{\theta})d\mathbf{y}.

Coverage Theorem

The average coverage, integrated with respect to the prior p(\boldsymbol{\theta}), is therefore \begin{aligned} \int \int I\{\theta_j \in \mathcal{A}_j(\mathbf{y})\} p(\mathbf{y}| \boldsymbol{\theta}) p(\boldsymbol{\theta}) d\mathbf{y}d\boldsymbol{\theta} &= \int \int I\{\theta_j \in \mathcal{A}_j(\mathbf{y})\} p(\boldsymbol{\theta}| \mathbf{y}) p(\mathbf{y}) d\mathbf{y}d\boldsymbol{\theta}\\ &= \int \int I\{\theta_j \in \mathcal{A}_j(\mathbf{y})\} p(\boldsymbol{\theta}| \mathbf{y}) d\boldsymbol{\theta}\, p(\mathbf{y}) d\mathbf{y}\\ &= \int (1 - \alpha) p(\mathbf{y}) d\mathbf{y}\\ &= 1 - \alpha. \end{aligned} Thus, the average coverage equals the nominal coverage.

MCP Example

Full Conditional Posterior

• The full conditional posterior for \beta_j is defined as the distribution for \beta_j conditional on \boldsymbol{\beta}_{-j} = \widehat{\boldsymbol{\beta}}_{-j}(\lambda).

• The partial residuals, \mathbf{r}_{j}, are then defined as \mathbf{y}- \mathbf{X}_{-j}\hat{\boldsymbol{\beta}}_{-j}(\lambda).

• A normal likelihood and Laplace prior are not conjugate.

\begin{aligned} p(\beta_j | \boldsymbol{\beta}_{-j}) &\propto \begin{cases} C_{-} \exp\{-\frac{n}{2\sigma^2} (\beta_j - (z_j + \lambda))^2\}, \text{ if } \beta_j < 0, \\ C_{+} \exp\{-\frac{n}{2\sigma^2} (\beta_j - (z_j - \lambda))^2\}, \text{ if } \beta_j \geq 0 \\ \end{cases} \end{aligned} where C_{-} = \exp(z_j \lambda n/\sigma^2), C_{+} = \exp(-z_j \lambda n/\sigma^2), z_{j} = \frac{1}{n} \mathbf{x}_{j}^{T}\mathbf{r}_{j}, and assuming that \mathbf{X} has been standardized s.t. \mathbf{x}_j^T\mathbf{x}_j = n.

• This formulation is attractive because it allows efficient computation of quantiles from the full conditional posterior.