Beyond support in two-stage variable selection

Abstract

Numerous variable selection methods rely on a two-stage procedure, where a sparsity-inducing penalty is used in the first stage to predict the support, which is then conveyed to the second stage for estimation or inference purposes. In this framework, the first stage screens variables to find a set of possibly relevant variables and the second stage operates on this set of candidate variables, to improve estimation accuracy or to assess the uncertainty associated to the selection of variables. We advocate that more information can be conveyed from the first stage to the second one: we use the magnitude of the coefficients estimated in the first stage to define an adaptive penalty that is applied at the second stage. We give the example of an inference procedure that highly benefits from the proposed transfer of information. The procedure is precisely analyzed in a simple setting, and our large-scale experiments empirically demonstrate that actual benefits can be expected in much more general situations, with sensitivity gains ranging from 50 to 100 % compared to state-of-the-art.

Appendices

Variational equivalence

We show below that the quadratic penalty in $\varvec{\beta }$ in Problem (2) acts as the Lasso penalty $\lambda \left\| \varvec{\beta }\right\| _{1}$.

Proof

The Lagrangian of Problem (2) is:

$$\begin{aligned} L(\varvec{\beta }) = J(\varvec{\beta }) + \lambda \sum _{j=1}^p \frac{1}{\tau _{j}} \beta _{j}^{2} + \nu _0 \bigg ( \sum _{j=1}^p \tau _j -\left\| \varvec{\beta }\right\| _{1} \bigg ) - \sum _{j=1}^p \nu _j \tau _j. \end{aligned}$$

Thus, the first order optimality conditions for $\tau _j$ are

$$\begin{aligned} \frac{\partial L}{\partial \tau _j}(\tau _j^\star ) = 0&\Leftrightarrow - \lambda \frac{\beta _{j}^2}{{\tau _j^\star }^2} + \nu _0 - \nu _j = 0 \\&\Leftrightarrow - \lambda \beta _{j}^2 + \nu _0 \,{\tau _j^\star }^2 - \nu _j \,{\tau _j^\star }^2 = 0 \\&\Rightarrow - \lambda \beta _{j}^2 + \nu _0 \,{\tau _j^\star }^2 = 0 , \end{aligned}$$

where the term in $\nu _j$ vanishes due to complementary slackness, which implies here $\nu _j \tau _j^\star = 0$. Together with the constraints of Problem (2), the last equation implies $\tau _j^\star = \left| \beta _{j}\right|$, hence Problem (2) is equivalent to

$$\begin{aligned} \min _{\varvec{\beta }\in \mathbb {R}^p} J (\varvec{\beta }) + \lambda \left\| \varvec{\beta }\right\| _{1} , \end{aligned}$$

which is the original Lasso formulation. $\square$

Efficient implementation

Permutation tests rely on the simulation of numerous data sampled under the null hypothesis distribution. The number of replications must be important to estimate the rather extreme quantiles we are typically interested in. Here, we use $B=1000$ replications for the $q=|\mathscr {S}_{\hat{\lambda }}|$ variables selected in the screening stage. Each replication involving the fitting of a model, the total computational cost for solving these B systems of size q on the q selected variables is $O(Bq(q^{3}+q^{2}n))$. However, block-wise decompositions and inversions can bring computing savings by a factor q.

First, we recall that the adaptive-ridge estimate, computed at the cleaning stage, is computed as

$$\begin{aligned} {\,\hat{\varvec{\beta }}} = \left( {\mathbf {X}}^{\top }\mathbf {X}+ \varvec{{\varLambda }}\right) ^{-1} \mathbf {X}^{\top }\mathbf {y}, \end{aligned}$$

where $\varvec{{\varLambda }}$ is the diagonal adaptive-penalty matrix defined at the screening stage, whose jth diagonal entry is ${\lambda }/{\tau _{j}^\star }$, as defined in (1–3).

In the F-statistic (4), the permutation affects the calculation of the larger model $\hat{\mathbf {y}}_{1}$, which is denoted $\hat{\mathbf {y}}_{1}^{(b)}$ for the bth permutation. Using a similar notation convention for the design matrix and the estimated parameters, we have $\hat{\mathbf {y}}_{1}^{(b)}={\mathbf {X}^{(b)}}\hat{\varvec{\beta }}^{(b)}$. When testing the relevance of variable j, ${\mathbf {X}^{(b)}}$ is defined as the concatenation of the permuted variable $\mathbf {x}_{j}^{(b)}$ and the other original variables: ${\mathbf {X}^{(b)}}= (\mathbf {x}_{j}^{(b)}, \mathbf {X}_{\backslash j}) = (\mathbf {x}_{j}^{(b)}, \mathbf {x}_1, ..., \mathbf {x}_{j-1}, \mathbf {x}_{j+1}, ..., \mathbf {x}_{p})$. Then, $\hat{\varvec{\beta }}^{(b)}$ can be efficiently computed by using $a^{(b)}\in \mathbb {R}$, $\mathbf {v}^{(b)}\in \mathbb {R}^{q-1}$ and $\hat{\varvec{\beta }}_{\backslash j}\in \mathbb {R}^{q-1}$ defined as follows:

$$\begin{aligned} a^{(b)}&= (\Vert \mathbf {x}_{j}^{(b)}\Vert ^2_2 + {\varLambda }_{jj} ) - {\mathbf {x}_{j}^{(b)}}^{\top } \mathbf {X}_{\backslash j}{(\mathbf {X}_{\backslash j}^\top \mathbf {X}_{\backslash j}+ {\varvec{{\varLambda }}_{\backslash j}})}^{-1} \mathbf {X}_{\backslash j}^\top \mathbf {x}_{j}^{(b)}\\ \mathbf {v}^{(b)}&= -{ (\mathbf {X}_{\backslash j}^\top \mathbf {X}_{\backslash j}+ {\varvec{{\varLambda }}_{\backslash j}})}^{-1} \mathbf {X}_{\backslash j}^\top \mathbf {x}_{j}^{(b)}\\ \hat{\varvec{\beta }}_{\backslash j}&= {(\mathbf {X}_{\backslash j}^\top \mathbf {X}_{\backslash j}+ {\varvec{{\varLambda }}_{\backslash j}})}^{-1} \mathbf {X}_{\backslash j}^\top \mathbf {y}. \end{aligned}$$

Indeed, using the Schur complement, one writes $\hat{\varvec{\beta }}^{(b)}$ as follows:

$$\begin{aligned} \hat{\varvec{\beta }}^{(b)}= \frac{1}{a^{(b)}} \begin{pmatrix} 1 \\ \mathbf {v}^{(b)}\end{pmatrix} \begin{pmatrix} 1&{\mathbf {v}^{(b)}}^{\top } \end{pmatrix} \begin{pmatrix} {\mathbf {x}_{j}^{(b)}}^{\top } \mathbf {y}\\ \mathbf {X}_{\backslash j}^{\top } \mathbf {y}\end{pmatrix} + \begin{pmatrix} 0 \\ \hat{\varvec{\beta }}_{\backslash j}\end{pmatrix} . \end{aligned}$$

Hence, $\hat{\varvec{\beta }}^{(b)}$ can be obtained as a correction of the vector of coefficients $\hat{\varvec{\beta }}_{\backslash j}$ obtained under the smaller model. The key observation to be made here is that $\mathbf {x}_{j}^{(b)}$ does not intervene in the expression ${(\mathbf {X}_{\backslash j}^\top \mathbf {X}_{\backslash j}+ {\varvec{{\varLambda }}_{\backslash j}})}^{-1}$, which is the bottleneck in the computation of $a^{(b)}$, $\mathbf {v}^{(b)}$ and $\hat{\varvec{\beta }}_{\backslash j}$. It can therefore be performed once for all permutations. Additionally, ${(\mathbf {X}_{\backslash j}^\top \mathbf {X}_{\backslash j}+ {\varvec{{\varLambda }}_{\backslash j}})}^{-1}$ can be cheaply computed from $\left( {\mathbf {X}}^{\top }\mathbf {X}+ \varvec{{\varLambda }}\right) ^{-1}$ as follows:

$$\begin{aligned}&{(\mathbf {X}_{\backslash j}^\top \mathbf {X}_{\backslash j}+ {\varvec{{\varLambda }}_{\backslash j}})}^{-1} = \left[ \left( {\mathbf {X}}^{\top }\mathbf {X}+ \varvec{{\varLambda }}\right) ^{-1} \right] _{\backslash j \backslash j} \\&\qquad -\left[ \left( {\mathbf {X}}^{\top }\mathbf {X}+ \varvec{{\varLambda }}\right) ^{-1} \right] _{\backslash j j} \left[ \left( {\mathbf {X}}^{\top }\mathbf {X}+ \varvec{{\varLambda }}\right) ^{-1} \right] _{j j}^{-1} \\&\qquad \times \left[ \left( {\mathbf {X}}^{\top }\mathbf {X}+ \varvec{{\varLambda }}\right) ^{-1} \right] _{j{\backslash j}} . \end{aligned}$$

Thus we compute $\left( {\mathbf {X}}^{\top }\mathbf {X}+ \varvec{{\varLambda }}\right) ^{-1}$ once, firstly correct for the removal of variable j, secondly correct for permutation b, thus eventually requiring $O(B(q^{3}+q^{2}n)))$ operations.