On the Selection Stability of Stability Selection and Its Applications

We consider a dataset $\mathcal{D}=\{(\boldsymbol{x}^{\top}_{i},y_{i})\}_{i=1}^{n}$, where each element consists of a univariate response $y_{i}\in\mathbb{R}$ and a $p$-dimensional vector of fixed covariates $\boldsymbol{x}^{\top}_{i}\in\mathbb{R}^{p}$. In the context of linear regression, it is typically assumed that pairs $(\boldsymbol{x}^{\top}_{i},y_{i})$ are independent and identically distributed (i.i.d.); when these are assumed to be random, but for simplicity, we assume that the covariate vector is fixed.

The linear regression model is formally expressed as $\boldsymbol{Y}=\beta_{0}+X\boldsymbol{\beta}+\boldsymbol{\varepsilon}$, where $X\in\operatorname{mat}(n\times p)$ represents the design matrix, $\boldsymbol{Y}\in\mathbb{R}^{n}$ is the $n$-dimensional vector corresponding to the univariate response variable, and $\beta_{0}\in\mathbb{R}$ denotes the intercept term. The vector of regression coefficients for the $p$ non-constant covariates is denoted by $\boldsymbol{\beta}\in\mathbb{R}^{p}$, while $\boldsymbol{\varepsilon}\in\mathbb{R}^{n}$ represents the $n$-dimensional vector of random errors. It is assumed that the components of $\boldsymbol{\varepsilon}$ are i.i.d. and independent of the covariates.

Variable selection, in accordance with the terminology of Meinshausen and Bühlmann (2010), typically involves categorizing covariates into two distinct groups: the signal group $S=\{k\neq 0\mid\beta_{k}\neq 0\}$ and the noise group $N=\{k\neq 0\mid\beta_{k}=0\}$, where $S\cap N=\emptyset$. The primary objective of variable selection is to accurately identify the signal group $S$.

The stability selection framework requires the implementation of an appropriate selection method. One commonly used technique for variable selection in the linear regression context is the LASSO estimator, formally defined as

where $\lambda\in\mathbb{R}^{+}$ represents the LASSO regularization parameter. The set of non-zero coefficients can be identified as $\hat{S}(\lambda)=\{k\neq 0\mid\hat{\beta}_{k}(\lambda)\neq 0\}$, derived from solving the LASSO equation through convex optimization.

Although alternative methods exist for variable selection in linear regression models, a detailed exploration of techniques for determining the regularization method and penalty term in stability selection lies beyond the scope of this paper.

Meinshausen and Bühlmann (2010) defined the stable set as

where $\Lambda$ represents the set of regularization values, $0<\pi_{\text{thr}}<1$ is the threshold for decision-making in variable selection, and $\hat{\Pi}_{j}^{\lambda}$ denotes the selection frequency of the $j$th variable given the regularization parameter $\lambda$.

Equation 1 identifies variables whose selection frequencies exceed $\pi_{\text{thr}}$ under the regularization value that maximizes their selection frequencies; therefore, it considers best-case of each variable for decision-making.

We now move forward to stable stability selection. To implement this, we first define a grid of regularization values, $\Lambda$, from which the optimal $\lambda$ for LASSO is to be identified. For this purpose, we use as the default choice the $\lambda$ values generated by the cv.glmnet function from the glmnet package in R (Friedman et al., 2010), applying a 10-fold cross-validation to the complete dataset $\mathcal{D}$. Alternative regularization values can also be employed. However, a comparison of different methods for generating the regularization grid falls outside the scope of this paper. The goal is to identify the optimal regularization value $\lambda_{\text{stable}}\in\Lambda$ that yields highly stable outcomes with the least possible loss in terms of predictive ability.

The following procedure is applied to each regularization value $\lambda\in\Lambda$. In each iteration of the stability selection process, a random sub-sample, comprising half the size of the original dataset $\mathcal{D}$, is randomly selected. The LASSO model is then fitted to this sub-sample using the fixed value of $\lambda$, and the binary selection outcomes are recorded as a row in the binary selection matrix $M(\lambda)\in\operatorname{mat}(B\times p)$, where $B$ denotes the number of sub-samples. Therefore, in the end, we obtain $|\Lambda|$ distinct selection matrices $M(\lambda)$.

In this context, $M(\lambda)_{bj}=1$ indicates that the $j$th variable is identified as part of the signal set $\hat{S}(\lambda)$ when LASSO is applied to the $b$th sub-sample given $\lambda$. In contrast, $M(\lambda)_{bj}=0$ signifies that the $j$th variable is classified as belonging to the noise set $\hat{N}(\lambda)$ after applying LASSO to the $b$th sub-sample given $\lambda$.

The stability estimator proposed by Nogueira et al. (2018) is defined as

where $s_{j}^{2}$ denotes the unbiased sample variance of the binary selection statuses of the $j$th variable, while $q(\lambda)$ denotes the average number of variables selected under the regularization parameter $\lambda$. The stability estimator $\hat{\Phi}$ is bounded by $\left[-\frac{1}{B-1},1\right]$ (Nogueira et al., 2018); the larger the $\hat{\Phi}$, the more stable is $M(\lambda)$.

To interpret the stability values of Equation (2), Nogueira et al. (2018) adopted the guidelines proposed by Fleiss et al. (2004). According to Nogueira et al. (2018), stability values that exceed $0.75$ indicate excellent agreement between selections beyond what would be expected by random chance, while values below $0.4$ signify poor agreement between them. Based on Nogueira et al. (2018), stability values that reside between these thresholds are categorized as indicative of intermediate to good stability.

By obtaining all $|\Lambda|$ selection matrices, we can estimate the stability of each using the formulation given in Equation (2). We propose that the optimal regularization value, denoted as $\lambda_{\text{stable}}$, is the smallest regularization value at which the stability measure surpasses 0.75, that is,

The threshold value of $0.75$ is somewhat arbitrary and may not be attainable in certain situations. We address this limitation in the following paragraph.

As we will see in Section 4, $\lambda_{\text{stable}}$ may not exist in certain practical applications due to procedural instability, which prevents the stability values from exceeding $0.75$. In such cases, we propose

where SD represents the standard deviation of the stability values.

Both upper-bounds of the Per-Family Error Rate, mentioned in Section 1, provided by Meinshausen and Bühlmann (2010) and Shah and Samworth (2013), depend on the values of $\lambda$ and $\pi_{\text{thr}}$. The calibration of these two values requires the arbitrary selection of one of these two parameters, which can be challenging to justify. By adopting the optimal regularization value, that is $\lambda=\lambda_{\text{stable}}$ or $\lambda=\lambda_{\text{stable-1sd}}$, the upper-bound can be explicitly tailored to the data owner’s preferences. If a pre-determined $\pi_{\text{thr}}$ is chosen, the Per-Family Error Rate is deterministically obtained; alternatively, if a fixed pre-determined Per-Family Error Rate value is desired, the corresponding $\pi_{\text{thr}}$ is enforced.

The asymptotic upper-bound provided by Meinshausen and Bühlmann (2010) for the Per-Family Error Rate is given by

where $q_{\Lambda}$ represents the average number of selected variables over the regularization grid $\Lambda$. In a manner consistent with the approach of Bodinier et al. (2023), we can derive a point-wise control version of Equation 5 by restricting $\Lambda$ to a single value, $\lambda$, without affecting the validity of the original equation. Consequently, we can replace $q_{\Lambda}$ with $q(\lambda)$, as introduced in Equation 2. By employing $\lambda_{\text{stable}}$ within this formula and utilizing $q(\lambda_{\text{stable}})$, we establish a two-way control: fixing $\pi_{\text{thr}}$ allows us to determine the upper-bound, and vice versa. A similar approach can be applied to the upper-bounds proposed by Shah and Samworth (2013).

Having $\lambda_{\text{stable}}$ or $\lambda_{\text{stable-1sd}}$ allows a fresh perspective on variable selection through stability selection. Meinshausen and Bühlmann (2010) introduces Equation (1) to identify stable variables, that is, those with selection frequencies that exceed $\pi_{\text{thr}}$ in the best-case scenario. Instead, we propose

Equation (6) represents the variables with selection frequencies higher than $\pi_{\text{thr}}$ under the $\lambda_{\text{stable}}$, which corresponds to the highly stable outcomes. We term this selection criterion stable stability selection because it leverages highly stable outcomes to facilitate variable selection. If $\lambda_{\text{stable}}$ does not exist, $\hat{S}^{\text{stable-1sd}}$ can be defined similarly.

As outlined in Section 1, as the number of sub-samples increases, the convergence of $\hat{\Phi}$ to the true stability $\Phi$ is ensured. To illustrate this, we define $M^{(t)}(\lambda)$ as the matrix containing the selection outcomes from the first $t$ sub-samples for a given $\lambda\in\Lambda$, that is, the first $t$ rows of $M(\lambda)$. The objective is to evaluate the stability of $M^{(t)}(\lambda)$ across successive sub-samples in order to determine the appropriate cut-off point for the process, that is, the number of sub-samples required to achieve convergence in stability. Beyond this threshold, additional sub-samples do not significantly change the stability of $M(\lambda)$. We propose plotting the stability values of the selection matrix over the sequential sub-sampling, that is $\hat{\Phi}(M^{(t)}(\lambda));\;t=2,3,\ldots,B$, against $t$ to monitor the convergence status of the stability values. Given that the asymptotic distribution of $\hat{\Phi}$ is established, a confidence bound can be drawn along the curve to reflect the inherent uncertainty of the stability estimator.

In Section 4, we present the results obtained from the application of the proposed approach to synthetic and real datasets, which are introduced in Section 3.

As synthetic data, we consider a dataset with a sample size of $n=50$ that includes $p=500$ predictor variables. In this scenario, the predictor variables $\boldsymbol{x}^{\intercal}_{i}$ are independently drawn from the Normal distribution $\mathcal{N}(0,\Sigma)$, where $\Sigma\in\operatorname{mat}(p,p)$ has a diagonal of ones and $\sigma_{jk}=\rho^{|j-k|}$ for $j\neq k$ where $\rho\in\{0.2,0.5,0.8\}$. The response variable is linearly dependent solely on the first two predictor variables, characterized by the coefficient vector $\boldsymbol{\beta}=(1.5,1.1,0,\ldots,0)^{\top}$, and the error term $\boldsymbol{\varepsilon}$ is an i.i.d. sample from the standard Normal distribution $\mathcal{N}(0,1)$.

For the first real example, we use the well-established ‘riboflavin’ dataset, which focuses on the production of riboflavin (vitamin B2) from various Bacillus subtilis. This dataset, provided by the Dutch State Mines Nutritional Products, is accessible via hdi R package (Dezeure et al., 2015). It comprises a single continuous response variable that represents the logarithm of the riboflavin production rate, alongside $p=4,088$ covariates, which correspond to the logarithm of expression levels for 4,088 bacterial genes. The primary objective of analyzing this dataset is to identify genes that influence riboflavin production, with the ultimate goal of genetically engineering bacteria to enhance riboflavin yield. Data were collected from $n=71$ relatively homogeneous samples, which were repeatedly hybridized during a feed-batch fermentation process involving different engineered strains and varying fermentation conditions. Bühlmann et al. (2014) employed stability selection with LASSO and identified three genes—LYSC_at, YOAB_at, and YXLD_at—as having a significant impact on the response variable.

As an additional real-world example, we investigate ‘Affymetrix Rat Genome 230 2.0 Array’ microarray data introduced by Scheetz et al. (2006). This dataset comprises $n=120$ twelve-week-old male rats, with expression levels recorded for nearly 32,000 gene probes for each rat. The primary objective of this analysis is to identify the probes most strongly associated with the expression level of the TRIM32 probe (1389163_at), which has been linked to the development of Bardet-Biedl syndrome (Chiang et al., 2006). This genetically heterogeneous disorder affects multiple organ systems, including the retina.

In accordance with the preprocessing steps outlined by Huang et al. (2008), we excluded gene probes with a maximum expression level below the $25$th percentile and those exhibiting an expression range smaller than $2$. This filtering process yielded a refined set of $p=3,083$ gene probes that demonstrated sufficient expression and variability for further analysis.

Here, we present the results obtained from the application of stable stability selection to synthetic data. We used the GitHub repository associated with Nogueira et al. (2018) to implement their stability estimator and to obtain confidence intervals¹¹1https://github.com/nogueirs/JMLR2018.

Figure 1 illustrates the selection stability of the stability selection on a grid of regularization values, applied to the synthetic datasets introduced in Section 3, using three different values of $\rho$ and choosing the number of sub-samples $B=500$. The horizontal dashed red lines indicate the stability thresholds of $0.4$ and $0.75$ as discussed in Section 2.

As illustrated in Figure 1, the regularization value that minimizes the cross-validation error, $\lambda_{\text{min}}$, falls within the poorly stable region. In addition, the regularization value that maintains the cross-validation error within one standard error of $\lambda_{\text{min}}$, referred to as $\lambda_{\text{1se}}$ (Hastie et al., 2009), lies within the poor to intermediate stability regions. In contrast, $\lambda_{\text{stable}}$, introduced in Equation (3), is specifically designed to fall within the region of excellent stability, where applicable.

Notably, in all three correlation scenarios, the three regularization selection methods successfully identify relevant variables with high selection frequencies, with a minimum selection frequency of $0.994$. However, since $\lambda_{\text{stable}}$ is larger than the two other regularization values, it applies a stronger shrinkage to variables. To further examine this in terms of model accuracy, for each correlation scenario, we generate $n^{\prime}=25$ additional test samples from the corresponding distribution described in Section 3. At each iteration of stability selection, we predict the response variable for the test samples using the estimated model and, ultimately, we aggregate all the mean squared error (MSE) values obtained over the corresponding $\lambda$ by averaging.

As demonstrated in Figure 2, across all three correlation scenarios, the improvement in stability when transitioning from $\lambda_{\text{1se}}$ to $\lambda_{\text{stable}}$ far outweighs the reduction in accuracy. This outcome aligns closely with the findings of Nogueira et al. (2018), who noted that “All these observations show that stability can potentially be increased without loss of predictive power” and “We also observe that pursuing stability may help identifying the relevant set of features”.

To demonstrate the convergence of the stability estimator $\hat{\Phi}$, for each correlation scenario, we used $\lambda_{\text{stable}}$ and calculated the stability of the selection matrix $M(\lambda_{\text{stable}})$ through sequential sub-sampling, using Equation (2). The stability values across iterations are shown in Figure 3. The blue shading around the line represents the $95\%$ confidence interval for $\hat{\Phi}$. As can be observed, the interval narrows with increasing iterations, indicating a reduction in the estimator’s uncertainty. Figure 3 reveals that, across all three scenarios, after approximately $B^{*}\approx 200$ iterations, the stability estimator converges, with no significant change thereafter in its value. By monitoring the stability trajectory throughout the process, valuable insight can be gained in determining the optimal cut-off point, that is, the optimal number of sub-samples in terms of stability of the process.

Next, we apply stable stability selection to the riboflavin dataset. As described in Section 3, the dataset consists of $p=4,088$ genes, the main objective being to identify the key genes that influence riboflavin production. As before, we choose the number of sub-samples $B=500$. The predictor variables are standardized using the scale function prior to being input into LASSO.

In this problem, $\lambda_{\text{stable}}$ does not exist, so we use $\lambda_{\text{stable-1sd}}$ instead. After $B$ iterations, four genes YXLD_at ($\hat{\Pi}_{\texttt{YXLD\_at}}^{\lambda_{\text{stable-1sd}}}=0.606$), YOAB_at ($\hat{\Pi}_{\texttt{YOAB\_at}}^{\lambda_{\text{stable-1sd}}}=0.558$), LYSC_at ($\hat{\Pi}_{\texttt{LYSC\_at}}^{\lambda_{\text{stable-1sd}}}=0.540$), and YCKE_at ($\hat{\Pi}_{\texttt{YCKE\_at}}^{\lambda_{\text{stable-1sd}}}=0.532$) have selection frequencies greater than $0.5$.

Figure 4 illustrates the stability of stability selection when applied to the riboflavin dataset. As shown in Figure 4(a), $\lambda_{\text{stable}}$ does not exist in this example, as the stability values do not surpass the $0.75$ threshold. Figure 4(b), generated with $\lambda_{\text{stable-1sd}}$, indicates the convergence slightly above $0.2$ after about $200$ iterations. These results imply that stability selection with LASSO exhibits poor selection stability within this dataset. Consequently, the findings presented by Bühlmann et al. (2014) may require careful reconsideration, particularly since the selection frequencies of the three identified genes are relatively low.

Finally, we apply the proposed methodology to the rat microarray data. As mentioned in Section 3, the data consists of $p=3,083$ gene probes. The main aim for this data is to identify probes that influence the TRIM32 probe. Again, we choose the number of sub-samples $B=500$. The predictor variables are standardized using the scale function prior to being input into LASSO.

Due to the instability of the outcomes, $\lambda_{\text{stable}}$ does not exist in this problem. After $B$ iterations, three probes have selection frequencies greater than $0.5$: probe 1390539_at ($\hat{\Pi}_{\texttt{1390539\_at}}^{\lambda_{\text{stable-1sd}}}=0.640$), probe 1389457_at ($\hat{\Pi}_{\texttt{1389457\_at}}^{\lambda_{\text{stable-1sd}}}=0.570$), and probe 1376747_at ($\hat{\Pi}_{\texttt{1376747\_at}}^{\lambda_{\text{stable-1sd}}}=0.564$).

As shown in Figure 5, similar to the results of the riboflavin dataset, the procedure cannot be considered stable. Specifically, Figure 5(a) reveals that $\lambda_{\text{stable}}$ does not exist, as the stability values do not exceed the threshold of $0.75$. For Figure 5(b), we used $\lambda_{\text{stable-1sd}}$. Figure 5(b) shows that the stability values converged to approximately $0.15$.

The results illustrated in figures 4 and 5 indicate that while stability selection is a valuable approach, more can be revealed by focusing on the stability of stability selection. As the complexity of the data (in terms of interdependencies between variables, number of variables subject to selection, etc.) increases, the stability values demonstrate a significant decline. Therefore, it is essential to evaluate the stability of its selections before drawing any conclusions about the importance of variables based on the outcomes of stability selection. These findings emphasize the need to carefully assess the stability of the entire process when using stability selection. Neglecting this evaluation may result in overconfidence in the results obtained, even when they lack stability. Given that such analyses are often underrepresented in the literature on stability selection, domain experts may benefit from carefully considering previously published findings.

The convergence plots for both synthetic and real examples indicate that, regardless of whether the convergence value is high or low, stability values tend to stabilize after $B^{*}\approx 200$ iterations. This observation could provide a useful rule of thumb for determining the required number of sub-samples in the stability selection framework.

Nogueira et al. (2018) suggested that stability and accuracy can be analyzed using the concept of the Pareto front (Pareto, 1896), which identifies regularization values that are not dominated by any other in terms of both criteria. A regularization value is considered Pareto optimal if no other value on the regularization grid offers both higher stability and higher accuracy. We call this approach stability-accuracy selection, which seeks to balance these two metrics. In this paper, we introduced the stable stability selection, which focuses on identifying the regularization value that achieves a high stability of the procedure with the least possible loss in accuracy.

As an experiment, in our synthetic data analysis we consider as a measure of accuracy (-MSE), that is, the negative of the mean squared error. Given that a Pareto optimal solution is not necessarily unique, we select the Pareto solution that maximizes the sum of accuracy and stability. In Figure 6(a), $\lambda_{\text{Pareto}}$ is located close to $\lambda_{\text{stable}}$, while in Figures 6(b) and 6(c), they coincide exactly. Interestingly, $\lambda_{\text{stable}}$ is also a Pareto solution in Figure 6(a). It is now pertinent to discuss the relationship between $\lambda_{\text{stable}}$ and Pareto optimality.

The two assumptions in Corollary 1 are both natural and justifiable. Increasing the regularization value is expected to allow the model to achieve an optimal balance of sparsity by shrinking irrelevant variables, thereby enhancing stability. The regularization value $\lambda_{\text{stable}}$ is intended to signify the point where the model attains high stability with minimal loss of accuracy. Consequently, we anticipate that the stability curve exhibits a non-decreasing trend prior to reaching $\lambda_{\text{stable}}$. However, since excessive regularization can lead to underfitting, we expect that the stability curve will flatten or decrease after reaching one or more peaks beyond $\lambda_{\text{stable}}$.

Regarding the second assumption, it is expected that increasing regularization helps the model eliminate irrelevant variables, improving its ability to predict the response variable. The regularization values $\lambda_{\text{min}}$ and $\lambda_{\text{1se}}$ are designed to capture this behavior, marking the points where regularization is most effective in terms of predictive ability. Beyond these points, we anticipate that the loss curve will flatten or increase. Since $\lambda_{\text{stable}}$ is intended to trade a small amount of accuracy for greater stability, we expect $\lambda_{\text{stable}}$ to occur after $\lambda_{\text{1se}}$, that is, $\lambda_{\text{stable}}>\lambda_{\text{1se}}>\lambda_{\text{min}}$, which implies that the loss curve is expected to be non-decreasing after $\lambda_{\text{stable}}$.

Thus, according to Corollary 1, $\lambda_{\text{stable}}$ represents a stability-accuracy solution for the problem of regularization tuning, ensuring high stability while minimizing loss in accuracy.