Variable Selection Using Bayesian Additive Regression Trees

Motivated by the boosting algorithm and built on the Bayesian classification and regression tree algorithm ([8]), BART ([9]) is a Bayesian approach to nonparametric function estimation and inference using a sum of trees. Consider the problem of making inference about an unknown function $f_{0}$ that predicts a continuous response $y$ using a $p$-dimensional vector of predictors $\boldsymbol{x}=(x_{1},\cdots,x_{p})^{\intercal}$ when

BART models $f_{0}$ by a sum of $M$ Bayesian regression trees, i.e.,

where $g(\boldsymbol{x};T_{m},\boldsymbol{\mu}_{m})$ is the output of $\boldsymbol{x}$ from a single regression tree. Each $g(\boldsymbol{x};T_{m},\boldsymbol{\mu}_{m})$ is determined by the binary tree structure $T_{m}$ consisting of a set of splitting rules and a set of terminal nodes, and the vector of parameters $\boldsymbol{\mu}_{m}=(\mu_{m,1},\cdots,\mu_{m,b_{m}})$ associated with the $b_{m}$ terminal nodes of $T_{m}$, such that $g(\boldsymbol{x};T_{m},\boldsymbol{\mu}_{m})=\mu_{m,l}$ if $\boldsymbol{x}$ is associated with the $l^{th}$ terminal node of the tree $T_{m}$.

A regularization prior consisting of three components, (1) the $M$ independent trees structures $\{T_{m}\}_{m=1}^{M}$, (2) the parameters $\{\boldsymbol{\mu}_{m}\}_{m=1}^{M}$ associated with the terminal nodes given the trees $\{T_{m}\}_{m=1}^{M}$ and (3) the error variance $\sigma^{2}$ which is assumed to be independent of the former two, is specified for BART in a hierarchical manner. The posterior distribution of BART is sampled through a Metropolis-within-Gibbs algorithm ([17, 14, 18]). The Gibbs sampler involves $M$ successive draws of $(T_{m},\boldsymbol{\mu}_{m})$ followed by a draw of $\sigma^{2}$. A key feature of this sampler is that it employs Bayesian backfitting ([16]) to sample each pair of $(T_{m},\boldsymbol{\mu}_{m})$.

BART estimates $f_{0}(\boldsymbol{x})=E(Y\mid\boldsymbol{x})$ by taking the average of the posterior samples of $f(\boldsymbol{x})$ after a burn-in period. As the trees structures are being updated through the MCMC chain, the model space is being searched by BART. In light of this, the estimation of $f_{0}(\boldsymbol{x})$ can be considered as a model averaging estimation with each posterior sample of the trees structures treated as a model and therefore BART can be regarded as a Bayesian model selection approach.

The continuous BART model (1–2) can be extend to a binary regression model for a binary response $Y$ by specifying

where $f$ is the sum-of-trees function in (2) and $\Phi$ is the cumulative distribution function of the standard normal distribution. Posterior sampling can be done by applying the Metropolis-within-Gibbs sampler of continuous BART to the latent variables which are obtained via the data augmentation approach of [1].

In this section, we review the BART variable inclusion proportion, a variable importance measure produced by BART, and three existing variable selection methods based on BART.

In Section 2.2.2, we show an example where the approach of [5] fails to include relevant binary predictors. The intuition of this approach is that a relevant predictor is expected to appear more often in the model built on the original dataset than that built on a null dataset. In this section, we revisit Example 2 and show that the intuition does hold for both relevant continuous and relevant binary predictors, but that the increase in usage of a relevant binary predictor is often offset by the increase of relevant continuous predictors, when the number of predictors is not large, thereby making the true VIP of a relevant binary predictor insignificant with respect to (w.r.t.) the corresponding null VIPs. In light of this observation, we modify BART VIP with the type information of predictors to help identify relevant binary predictors.

Direct studying BART VIP is challenging because each $v_{j}$ consists of $K$ different components $\{c_{jk}/c_{\cdot k}\}_{k=1}^{K}$. To make the analysis more interpretable, we introduce an approximation of BART VIP, which only contains one component.

The following lemma shows that $\tilde{v}_{j}$ defined above can be used to approximate BART VIP $v_{j}$ under mild conditions.

Consider the approach of [5]. For every $j=1,\cdots,p$, $k=1,\cdots,K$, $r=1,\cdots,L_{\text{rep}}$ and $l=1,\cdots,L$, let $c_{l,jk}^{*}$ be the number of splitting nodes using the predictor $x_{j}$ in the $k^{\text{th}}$ posterior sample of the BART model built on the $l^{\text{th}}$ null dataset and $c_{r,jk}$ be that of the $r^{\text{th}}$ repeated BART model built on the original dataset. Following the notation in Definition 3.1, write $c_{l,j\cdot}^{*}=\sum_{k=1}^{K}c_{l,jk}^{*}$, $c_{l,\cdot\cdot}^{*}=\sum_{j,k=1}^{p,K}c_{l,jk}^{*}$, $c_{r,j\cdot}=\sum_{k=1}^{K}c_{r,jk}$, and $c_{r,\cdot\cdot}=\sum_{j,k=1}^{p,K}c_{r,jk}$. By Lemma 3.2, for each predictor $x_{j}$, we can use $\tilde{v}_{lj}^{*}=c_{l,j\cdot}^{*}/c_{l,\cdot\cdot}^{*}$ to approximate $v_{lj}^{*}$, the VIP obtained from the BART model on the $l^{\text{th}}$ null dataset, and $\tilde{v}_{rj}=c_{r,j\cdot}/c_{r,\cdot\cdot}$ to approximate $v_{rj}$, the VIP obtained from the $r^{\text{th}}$ repeated BART model on the original dataset. Denote by $\bar{\tilde{v}}_{j}=(\sum_{r=1}^{L_{\text{rep}}}\tilde{v}_{rj})/L_{\text{rep}}$ the averaged approximate VIP for the predictor $x_{j}$ across $L_{\text{rep}}$ repetitions of BART on the original dataset. According to the local threshold of [5], a predictor $x_{j}$ is selected only if the average VIP $\bar{\tilde{v}}_{j}$ exceeds the $1-\alpha$ quantile of the empirical distribution $\sum_{l=1}^{L}\delta_{\tilde{v}_{lj}^{*}}(\cdot)$, i.e.,

Similar to the approximation of BART VIP, we can approximate $(\sum_{r=1}^{L_{\text{rep}}}c_{r,j\cdot}/c_{r,\cdot\cdot})/L_{\text{rep}}$ in (5) by $\bar{c}_{L_{\text{rep}},j\cdot}/\bar{c}_{L_{\text{rep}},\cdot\cdot}$ and therefore the selection criteria (5) can be rewritten as

where $\bar{c}_{L_{\text{rep}},j\cdot}=(\sum_{r=1}^{L_{\text{rep}}}c_{r,j\cdot})/L_{\text{rep}}$ and $\bar{c}_{L_{\text{rep}},\cdot\cdot}=(\sum_{r=1}^{L_{\text{rep}}}c_{r,\cdot\cdot})/L_{\text{rep}}$.

The ratio $\bar{c}_{L_{\text{rep}},j\cdot}/c_{l,j\cdot}^{*}$ on the left hand side (LHS) of the inequality inside the sum of (6) represents the average increment in usage of a predictor $x_{j}$ in a BART model built on the original dataset, compared to a null dataset; the ratio $\bar{c}_{L_{\text{rep}},\cdot\cdot}/c_{l,\cdot\cdot}^{*}$ on the right hand side (RHS) represents the average increment in the total number of the splitting nodes over all the posterior samples of a BART model built on the original dataset, compared to a null dataset. Figure S.3 of the Supplementary Material shows the overall counts ratio $\bar{c}_{L_{\text{rep}},\cdot\cdot}/c_{l,\cdot\cdot}^{*}$ and the counts ratio $\bar{c}_{L_{\text{rep}},j\cdot}/c_{l,j\cdot}^{*}$ for each predictor in Example 2. Both relevant continuous ($x_{11}$, $x_{12}$ and $x_{13}$) and relevant binary ($x_{1}$ and $x_{2}$) predictors are indeed used more frequently in the BART model built on the original dataset, i.e., $\bar{c}_{L_{\text{rep}},j\cdot}/c_{l,j\cdot}^{*}>1$, but the increment in usage of a relevant binary predictor is not always greater than the increment in the total number of splits, i.e., $\bar{c}_{L_{\text{rep}},j\cdot}/c_{l,j\cdot}^{*}\ngtr\bar{c}_{L_{\text{rep}},\cdot\cdot}/c_{l,\cdot\cdot}^{*}$. A binary predictor can only be used once at most in each path of a tree because it only has one split value, and a few splits are sufficient to provide the information about a binary predictor. As a result, $\bar{c}_{L_{\text{rep}},j\cdot}$, the number of splits using a binary predictor, is limited by its low cardinality, even if it is relevant. In addition, the number of splits $c_{l,j\cdot}^{*}$ using a binary predictor in a BART model built on a null dataset is roughly $1/p$ of the total number of splits in the ensemble ([5]). As such, when the number of predictors is small, $c_{l,j\cdot}^{*}$ is close to $\bar{c}_{L_{\text{rep}},j\cdot}$ (see the left subfigure of Figure S.4 of the Supplementary Material), thereby making the counts ratio $\bar{c}_{L_{\text{rep}},j\cdot}/c_{l,j\cdot}^{*}$ for a binary predictor close to $1$. The overall counts ratio $\bar{c}_{L_{\text{rep}},\cdot\cdot}/c_{l,\cdot\cdot}^{*}$ is also close to $1$ because of the BART regularization prior. In fact, given a fixed number of trees in an ensemble, the total number of splits in the ensemble is controlled by the splitting probability, the probability of splitting a node into two child nodes, which is the same for BART models built on the original dataset and null datasets, so the total number of splits for the original dataset and a null dataset, $\bar{c}_{L_{\text{rep}},\cdot\cdot}$ and $c_{l,\cdot\cdot}^{*}$, are of similar magnitude (see the right subfigure of Figure S.4). As such, the increment in usage of a relevant binary predictor $\bar{c}_{L_{\text{rep}},j\cdot}/c_{l,j\cdot}^{*}$ is easily offset by the increment in the total number of splits $\bar{c}_{L_{\text{rep}},\cdot\cdot}/c_{l,\cdot\cdot}^{*}$.

Figure S.4 of the Supplementary Material also shows that the increment in the total number of splits is primarily driven by those relevant continuous predictors, so essentially the increment in usage of a relevant binary predictor is offset by that of relevant continuous predictors. To alleviate this behavior, we modify Definition 2.1 using the type information of predictors.

Define $\delta_{jj^{\prime}}=\mathbb{I}(x_{j}\text{ and }x_{j^{\prime}}\text{ are of the same type})$ for any $j=1,\cdots,p$ and $j^{\prime}=1,\cdots,p$. Write $\bar{c}_{L_{\text{rep}},\text{type of }x_{j},\cdot}=(\sum_{r=1}^{L_{\text{rep}}}\sum_{j^{\prime}=1}^{p}c_{r,j^{\prime}\cdot}\cdot\delta_{jj^{\prime}})/L_{\text{rep}}$ and $c_{l,\text{type of }x_{j},\cdot}^{*}=\sum_{j^{\prime}=1}^{p}c_{l,j^{\prime}\cdot}^{*}\cdot\delta_{jj^{\prime}}$. The permutation-based variable selection approach (Algorithm 1) using BART within-type VIP as the variable importance measure selects a predictor $x_{j}$ only if

For a binary predictor, the selection criteria (7) removes the impact of continuous predictors, thereby not diluting the effect of a relevant binary predictor. Figure 2 shows the variable selection result of this approach. As opposed to the approach of [5], both the two relevant binary predictors $x_{1}$ and $x_{2}$ are clearly identified by this approach.

While BART within-type VIP works well in including relevant predictors in presence of a similar number of binary and continuous predictors, it becomes problematic if the number of predictors of one type is low. For example, when there are $99$ continuous predictors and $1$ binary predictor, the BART within-type VIP of the binary predictor is always $1$ for both the original dataset and a null dataset. As such, the binary predictor is never selected by the permutation-based approach, whether it is relevant or not. Furthermore, when there are more than two types of predictors, the computation of BART within-type VIP becomes complicated. Hence, we propose a new type of variable importance measure based on the Metropolis ratios calculated in the Metropolis-Hastings steps for sampling new trees, which are not affected by the issues above.

As mentioned in Section 2.1, a key feature of the Metropolis-within-Gibbs sampler is the Bayesian backfitting procedure for sampling $(T_{m},\boldsymbol{\mu}_{m})$, $1\leq m\leq M$, where each $T_{m}$ is fit iteratively using the residual response $r_{-m}=y-\sum_{m^{\prime}\neq m}g(\boldsymbol{x};T_{m^{\prime}},\boldsymbol{\mu}_{m^{\prime}})$, with all the other trees $T_{-m}=\{T_{m^{\prime}}\}_{m^{\prime}\neq m}$ held constant. Thus, each $(T_{m},\boldsymbol{\mu}_{m})$ can be obtained in two sequential steps:

where $\boldsymbol{r}_{-m}$ is the vector of residual responses. The distribution of $T_{m}\mid\boldsymbol{r}_{-m},\sigma^{2}$ has a closed form up to a normalizing constant and therefore can be sampled by a Metropolis-Hastings algorithm. [8] develop a Metropolis-Hastings algorithm that proposes a new tree based on the current tree using one of the four moves: BIRTH, DEATH, CHANGE and SWAP. The BIRTH and DEATH proposals are the essential moves for sufficient mixing of the Gibbs sampler, so both the CRAN R package BART ([24]) and our package BartMixVs ([22]) implement these two proposals, each with equal probability. A BIRTH proposal turns a terminal node into a splitting node and a DEATH proposal replaces a splitting node leading to two terminal nodes by a terminal node. Thus, the proposed tree $T^{*}$ is identical to the current tree $T_{m}$ except growing/pruning one splitting node.

We consider using the Metropolis ratio for accepting a BIRTH proposal to construct a variable importance measure. For convenience, we suppress the subscript $m$ in the following. The Metropolis ratio for accepting a BIRTH proposal at a terminal node $\eta$ of the current tree $T$ can be written as

where

The BART prior splits a node of depth $d$ into two child nodes of $d+1$ with probability $\gamma(1+d)^{-\beta}$ for some $\gamma\in(0,1)$ and $\beta\in[0,\infty)$. Given that, the untruncated Metropolis ratio $r(\eta)$ can be explicitly expressed as the product of three ratios:

where $b$ is the number of terminal nodes in the current tree $T$ and $d_{\eta}$ is the depth of the node $\eta$ in the current tree $T$. The derivation can be found in Section S.2 of the Supplementary Material.

Figure S.5 of the Supplementary Material shows the nodes ratios for different numbers of terminal nodes and the depth ratios for different depths when the hyper-parameters $\gamma$ and $\beta$ are set as default, i.e., $\gamma=0.95$ and $\beta=2$. From the figure, we can see that the product of nodes ratio and depth ratio is mostly affected by the depth ratio which decreases as the proposed depth $d_{\eta}$ increases, implying that $r(\eta)$ assigns a greater value to a shallower node which typically contains more samples. Since the likelihood ratio indicates the conditional improvement to fitting brought by the new split, the untruncated Metropolis ratio $r(\eta)$ considers both the proportion of samples affected at the node $\eta$ and the conditional improvement to fitting brought by the new split at the node $\eta$. However, it is not appropriate to directly use $r(\eta)$ to construct a variable importance measure, due to the occurrence of extremely large $r(\eta)$ which dominates the others. Instead, we use the Metropolis ratio in (8) to construct the following BART Metropolis importance.

The intuition of BART MI is that a relevant predictor is expected to be consistently accepted with a high Metropolis ratio. Figure 3 shows that the Metropolis ratios of the splits using relevant predictors ($x_{1}$, $x_{2}$, $x_{11}$, $x_{12}$ and $x_{13}$) as the split variable have a higher median and smaller standard error than those using irrelevant predictors, implying that relevant predictors are consistently accepted with a higher Metropolis ratio. BART MI and BART VIP have similar forms. The key difference is the “kernel” they use. While BART VIP uses the number of splits $c_{jk}$ which tends to be biased against binary predictors, BART MI uses the average Metropolis acceptance ratio $\tilde{u}_{jk}$ which is similar between the relevant binary predictors and relevant continuous predictors, as shown in Figure 3. Hence, the average Metropolis acceptance ratio $\tilde{u}_{jk}$ not only helps BART MI distinguish relevant and irrelevant predictors, but also helps BART MI not be biased against predictors of a certain type.

The vector of average Metropolis acceptance ratios $(\tilde{u}_{1k},\cdots,\tilde{u}_{pk})$ does not sum to $1$ and each of them varies from $0$ to $1$. We normalize the vector $(\tilde{u}_{1k},\cdots,\tilde{u}_{pk})$ in (10) based on the idea that $\tilde{u}_{jk}$’s are correlated in the sense that once some predictors that explain the response more are accepted with higher Metropolis ratios, the rest of predictors will be accepted with lower Metropolis ratios or not be used. In addition, we take the average over $K$ posterior samples in (10) because averaging further helps identify relevant predictors. In fact, an irrelevant predictor can be accepted with a high Metropolis ratio at some MCMC iterations, but it can be quickly removed from the model by the DEATH proposal. As a result, the irrelevant predictor can have a few large $\tilde{u}_{jk}$’s and many small $\tilde{u}_{jk}$’s, thereby resulting in a small BART MI $v_{j}^{\text{MI}}$, as shown in Figure 3.

We further explore BART MI under null settings. We create a null dataset of Example 2 and fit a BART model on it. Figure S.6 of the Supplementary Material shows that not all the MIs, $v_{j}^{\text{MI}}$’s, converge to $1/p=0.05$, implying that $1/p$ may not be a good selection threshold for $v_{j}^{\text{MI}}$’s. We repeat the experiment of [5] for BART MI under null settings to explore the variation among $v_{j}^{\text{MI}}$’s. Specifically, we create $100$ null datasets of the data generated from Example 2, and for each null dataset, we run BART $50$ times with different initial values of the hyper-parameters. Let $v_{ijk}^{\text{MI}}$ be the BART MI for the predictor $x_{j}$ from the $k^{\text{th}}$ BART model on the $i^{\text{th}}$ null dataset. We investigate three nested variances listed in Table 1.

Results of the experiment are shown in Figure S.7 of the Supplementary Material. The non-zero $s_{ij}$’s imply that there is variation in $v_{j}^{\text{MI}}$’s among the repetitions of BART with different initial values, so it is necessary to average MIs over a certain number of repetitions of BART to get stable MIs. In practice, we find that among the repetitions of BART, irrelevant predictors occasionally get outliers of MI. Thus, instead of averaging MIs over BART repetitions, we take the median of them. Similar to [5], we also find that $10$ repetitions of BART is sufficient to get stable MIs (see Figure S.8 of the Supplementary Material). Figure S.7 also shows that the second type of standard deviation $s_{j}$ is significantly greater than the median of $s_{ij}$’s for each predictor, which suggests that BART tends to overfit noise. Since the overall MI $\bar{v}_{\cdot\cdot\cdot}^{\text{MI}}$ is always $1/p$, the relatively large $s$ indicates that not all the $\bar{v}_{\cdot j\cdot}^{\text{MI}}$ are approximately $1/p$, suggesting that it is not possible to determine a single selection threshold for all the MIs. Therefore, we employ the permutation-based approach (Algorithm 1) to select thresholds for each MI. A slight modification is applied to line 4 of Algorithm 1: we take the median, rather than the average, of BART MIs, over the repetitions of BART on the original dataset. The variable selection result of this approach for Example 2 is shown in Figure 4.

In this section, we propose an approach to select a subset of predictors with which the BART model gives the best prediction. In general, to select the best subset of predictors from the model space consisting of $2^{p}$ possible models, one needs a search strategy over the model space and decision-making rules to compare models. Here we use the backward elimination approach to search the model space, mainly because it is a deterministic approach and is efficient with a moderate number of predictors. Typically, backward selection starts with the full model with all the predictors, followed by comparing the deletion of each predictor using a chosen model fit criterion and then deleting the predictor whose loss gives the most statistically insignificant deterioration of the model fit. This process is repeated until no further predictors can be deleted without a statistically insignificant loss of fit. Some popular model fit criterions such as AIC and BIC are not available for BART, because the maximum likelihood estimates are unavailable and the number of parameters in the model is hard to determine. To overcome this issue, we propose a backward selection approach with two easy-to-compute selection criterions as the decision-making rules.

We first split the dataset $\{y_{i},\boldsymbol{x}_{i}\}_{i=1}^{n}$ into a training set and a test set before starting the backward procedure. To measure the model fit, we make use of the mean squared errors (MSE) of the test set

where $\hat{f}$ is the estimated sum-of-trees function. Figure 5 shows that $\text{MSE}_{\text{test}}$ can distinguish between acceptable models and unacceptable models, where acceptable models are defined as those models including all the relevant predictors and unacceptable models are defined as those missing some relevant predictors. At each backward selection step, we run BART models on the training set and compute the MSE of the test set. The model with the smallest $\text{MSE}_{\text{test}}$ is chosen as the “winner” at that step. However, due to the lack of stopping rules under such nonparametric setting, the backward selection has to continue until there is only one predictor in the model, and ultimately returns $p$ “winner” models with different model sizes ranging from $1$ to $p$.

Given the $p$ “winner” models, we compare them using the expected log pointwise predictive density based on leave-one-out (LOO) cross validation:

where $\boldsymbol{y}_{-i}=\{y_{j}\}_{j\neq i}$, $\theta=\{\{T_{m},\boldsymbol{\mu}_{m}\}_{m=1}^{M},\sigma^{2}\}$ represents the set of BART parameters, $f\left(y_{i}\mid\boldsymbol{y}_{-i}\right)$ is the predictive density of $y_{i}$ given $\boldsymbol{y}_{-i}$, $f(y_{i}\mid\theta)$ is the likelihood of $y_{i}$, and $f(\theta\mid\boldsymbol{y}_{-i})$ is the posterior density of BART given the observations $\boldsymbol{y}_{-i}$. The quantity $\text{elpd}_{\text{loo}}$ not only measures the predictive capability of a model, but also penalizes the complexity of the model, i.e., the number of predictors used in the model. Hence, we choose the model with the largest $\text{elpd}_{\text{loo}}$ as the best model.

Direct computing (11) involves fitting BART $n$ times. To avoid this, we use the approach of [27] that estimates each $f(y_{i}\mid\boldsymbol{y}_{-i})$ using importance sampling with the full posterior distribution $f(\theta\mid\boldsymbol{y})$ as the sampling distribution and the smoothed ratios $\{f(\theta^{k}\mid\boldsymbol{y}_{-i})/f(\theta^{k}\mid\boldsymbol{y})\propto 1/f(y_{i}\mid\theta^{k})\}_{k=1}^{K}$ as the importance ratios, where $\{\theta^{k}\}_{k=1}^{K}$ are the $K$ posterior samples from the full posterior distribution $f(\theta\mid\boldsymbol{y})$. The smoothness is achieved by fitting a generalized Pareto distribution to the upper tail of each set of the importance ratios $1/f(y_{i}\mid\theta^{k})\}_{k=1}^{K}$, $i=1,\cdots,n$. Thus, $\text{elpd}_{\text{loo}}$ in (11) can be estimated by using the likelihoods $\{f(y_{i}\mid\theta^{k})=\phi(y_{i}\mid\hat{f}^{k}(\boldsymbol{x}_{i}),\sigma^{k})\}_{i,k=1}^{n,K}$ based on $K$ posterior samples $\{\hat{f}^{k}\}_{k=1}^{K}$ from the BART model fit on the dataset $\{y_{i},\boldsymbol{x}_{i}\}_{i=1}^{n}$, where $\hat{f}^{k}$ is the $k^{\text{th}}$ posterior sample of the sum-of-trees function (2) and $\phi(\cdot\mid\mu,\sigma)$ is the normal density with mean $\mu$ and standard error $\sigma$.

The approach above can also be extended to the probit BART model (3) by changing the $\text{MSE}_{\text{test}}$ to the mean log loss (MLL) of the test set

and replacing the normal likelihoods with the Bernoulli density $\{f(y_{i}\mid\theta^{k})=\hat{p}_{i}^{k}y_{i}+(1-\hat{p}_{i}^{k})(1-y_{i})\}_{i,k=1}^{n,K}$, where $\hat{p}_{i}^{k}=\Phi(\hat{f}^{k}(x_{i}))$ is the estimated probability based on the $k^{\text{th}}$ posterior sample of probit BART for the $i^{\text{th}}$ observation and $\hat{p}_{i}=\Phi(\sum_{k=1}^{K}\hat{f}^{k}(x_{i})/K)$. We summarize this algorithm in Algorithm 2. Although this algorithm requires fitting BART $p(p-1)/2$ times, models at the same backward selection step can be fitted in parallel, which implies that the time complexity $O(p(p-1)/2)$ can be reduced to $O(p)$ if there are sufficient computing resources.

In this setting, we consider two scenarios. Scenario C.C.1 is an example from [11]; Scenario C.C.2, which includes correlation between predictors, is borrowed from [29]. For both of them, we consider $n=500$, $p\in\{50,200\}$ and $\sigma^{2}=1$.