A Ranked Sparsity Extension to the Bayesian Information Criterion: A Tool for Selecting Variables from Multiple Data Modalities

1. Introduction

1.1. Background – Sparsity and Skepticism

When selecting a model, statisticians must traditionally balance many trade-offs. We are constantly trying to balance simplicity with predictive accuracy, parsimony with complexity, sparsity with saturation, and bias with variance. While we have a multitude of tools at our disposal, in most realistic situations, the reliability of these instruments can change given information we cannot access. Consider sparsity, for instance. A model can have many parameters, but a sparse model can only have a handful by definition. Does nature produce mostly sparse models? How could we possibly know this? Often times, we can understand more about a mechanism, or predict more accurately, by collecting new features. Perhaps, if we collect enough features about a generating mechanism, the remaining stochasticity might completely vanish, allowing us to predict newly generated data perfectly.

Putting the merits of determinism aside, we often take for granted that scientists know when to stop considering new features. With the advent of new computational methods and theory that help sift through arbitrarily large sets of features, scientists are understandably becoming increasingly ambitious in this regard. However, this Icarean ambition to include as many features as possible, while well-intentioned, tends to drive down the proportion of variables that contain real signal, ultimately leaving the scientific endeavor dead in the water1.

We define the “sparsity level” as the proportion of candidate variables (a.k.a. features, covariates, or predictors) that are inactive in the true generating model. Conversely, the “saturation level” is defined as the proportion of candidate variables that are active. The sparsity level is governed by a mix of what cannot be known about nature’s true generating model and what can (sometimes) be known about the ambition of a particular scientific project.

It is certainly easier to describe sparse models; with a postulated sparse (or parsimonious) model, we can be more precise, collect less data, and improve the interpretability of our final model. And we do have tools to help us decide between sparse or saturated models: model selection criteria. Here, we define such criteria broadly as measures that are used to evaluate whether a certain predictor should be included in the model or not. These criteria typically have a penalty term, which acts as a “cost-of-admission” for each predictor (i.e., predictors are required to add a certain level of improvement in fit in order to be admitted into the model). These popular tools come with many trade-offs we will describe. Bootstrapping and cross-validation (CV) are increasingly being used to overcome the limitations of traditional model selection criteria (most importantly the asymptotic assumptions). However, these methods can be cumbersome, both computationally and in terms of their interpretation.

In practice, the goal of an analysis often dictates which method/criterion is used to select a final model. If a sparse model is desired, the Bayesian Information Criterion (BIC) [1] can be used. If prediction is the only concern, the Akaike Information Criterion (AIC) [2] or CV can be used. Once it is recognized that we cannot realistically know whether the true model is sparse or not, this real trade-off (that depends on the true generating mechanism) is relegated to the background of practicality, where a glossed-over presumption is made about the generating mechanism: a prior belief that all of the covariates are equally worthy of entering into the model. We refer to this presumption as “covariate equipoise.” In situations where multiple levels of sparsity can be expected among covariates, we will show how the presumption of covariate equipoise can have detrimental effects on every one of these model selection tools. With this impetus, we have developed several simple solutions to this problem for many settings – settings which we refer to as “ranked sparsity”. A related technique for this problem, the Sparsity-Ranked Lasso (SRL), is a useful regularization-based approach for selecting interactions and polynomials. However, sometimes this regularization has unwanted side effects, especially in the case of interactions as the fitting procedure is not robust to affine transformations of the main predictors (see the `Centering to different values’ section of Peterson [3]). In this work, we present a separate, information-criterion-based ranked-sparsity methodology: the Ranked-sparsity extension to BIC (RBIC).

1.2. A Motivating Example

The intuition behind ranked sparsity is that in many cases, covariates should not be presumed to be created equal. Consider the following simplistic example. A research team wishes to predict a response y using 20 covariates and 500 observations. Carefully investigating each covariate’s functional form and examining residual diagnostics, their statistician selects several of the 20 covariates as important predictors in a model that optimally fits the data interpretably. However, wondering if better predictions might be possible, the team also tries some state-of-the-art technology to train an army of machine learning (ML) models that they use for an ostensibly unbeatable (though completely opaque) predictive model. Sure enough, on new data, the predictions from the ensemble are more accurate than those from the statistical model. Let’s consider why this can happen.

One key factor is that the army of ML models is not exactly concerned about treating each covariate equally. In trying to predict y, each member of this ML army (each weak learner) is free to grab from an urn of practically infinite predictors (e.g. random splits of each covariate interacted with other random splits of each covariate). The ensemble (the super learner) then observes which ones work best, and uses weights, tuning parameters, and extra-sample validation to cleverly aggregate the results and optimize them for prediction. With 20 predictors in her original bin of candidates (not including additional polynomial and spline terms), the statistician did not consider any interactions between covariates. After all, in order to check all pairwise interactions, she would have had to sift through 190 more predictors, each of which would have made the model’s interpretation more nuanced. To consider three-way interactions renders the selection process increasingly Icarean; there are 1140 of those! If she were to consider all possible interactions of any order, she would have to consider ${\sum }_{i=1}^{20}\left({\displaystyle \genfrac{}{}{0pt}{}{20}{i}}\right)=1,048,575$ candidate variables!

Given this explosion of derived interactions, it is no wonder that the prevailing guidance is to consider an interaction in a model-selection framework only when there is substantive evidence to justify it beforehand. Harrell, for instance, notes that higher-order interactions are typically ignored unless specified a priori based on knowledge of the subject matter, and recommends pre-specifying the small set of plausible interactions rather than testing all of them [4]. In other words, absent such justification, we assume a priori that nature’s generating mechanism is completely sparse in terms of interactions. However, when we zoom out and look at science more broadly, how do scientists ever identify a clinically relevant interaction, if the convention is never to model one without prior justification?

This heuristic is a tacit acknowledgment of ranked sparsity: most of the signal in a linear regression problem can be found without considering interactions, and so the full set of possible interactions tends to add little to no signal (and a great deal of noise). This presumption nonetheless sacrifices a non-negligible amount of signal that flexible machine learning methods can often recover, albeit at the cost of interpretability. Statisticians are hamstrung because model selection tools do not work well when there is a large differential in the sparsity level among candidate predictors; the performance of model selection criteria is directly tied to the “true” sparsity level. In Section 1.5, we illustrate how model selection criteria with higher penalties work better in situations with a high sparsity level, whereas ones with low penalties work better in situations with a high saturation level, demonstrating that the optimal penalty depends greatly on the true sparsity level in the covariate space. Therefore, when there is an expected disparity in sparsity among candidate predictors (such as interactions vs. main effects), model selection tools can and should account for this by using different penalties.

Our goals in this paper are as follows. After describing our notation and providing necessary background on existing tools, we motivate ranked sparsity methodology for covariate groups of different sizes and interactions, first by tying the performance of penalization in model selection criteria to the true saturation level, then by arguing why the expected saturation level for interactions and large feature sets should generally be presumed to be lower. Next, we identify RBIC as a novel extension of BIC that can be used to select from multiple covariate groups (and interactions) with different expected sparsity levels. We compare its performance to other methods in a simulation study and in two real data applications, and describe its user-friendly implementation via the ICtoolkit R package. We conclude with a discussion.

1.3. Traditional Model Selection Criteria

For a data generating model $G:\mathit{x}\to y$, where $\mathit{x}$ is a vector of p covariates, the goal of model selection is to build a linear model to predict future y observations via $f(\mathit{\beta },\mathit{x})={\mathit{x}}^T\mathit{\beta }$. Many model selection criteria have been explored over the course of the last century as solutions to this problem, and most can be broadly classified as either consistent or efficient [5]. Consistent model selection criteria, including, for example, BIC [1], the minimum description length criterion (MDL) [6], Extended BIC [7], aim to select the most parsimonious model with a requisite structure; as $n\to \infty$, a consistent criterion is guaranteed to select the correct model (with probability one), under the strong assumption that this model is included in the candidate collection [8]. In contrast, efficient model selection criteria are not, in general, consistent. Rather than recover a true structure, they aim to minimize the error of the model when predicting the response for new data; an efficient criterion will select a model whose predictive risk converges to that of the best candidate model available, in the sense that the ratio of the two tends to 1 as $n\to \infty$ [9]. Crucially, this benchmark is the best available model, so unlike consistency, efficiency does not presume that the true model lies within the candidate set [10]. Efficient model selection criteria include, for example, AIC [2], small-sample-corrected AIC (AICc) [11], and Mallows’ Cp [12].

When should we use an efficient criterion and when should we use a consistent one? When should we abandon them both? Besides the goals of a particular project, other major factors that come into play for this decision include (but are not limited to) the following: the sample size n (how much data we collected), the number of candidate predictors p (how ambitious we are), and the number of true predictors s (the number of truly active candidate predictors). With this notation, the sparsity level is $1-{\displaystyle \frac{s}{p}}$, and the saturation level is ${\displaystyle \frac{s}{p}}$. As $p\to n$, the asymptotic assumptions on which many model selection criteria are based become increasingly violated and their performance tends to suffer as a result. As ${\displaystyle \frac{s}{p}}\downarrow$ (i.e., the saturation level decreases), the performance will also vary by penalization (efficient vs. consistent criteria) in ways we explore in Section 1.5, but first we introduce model search frameworks.

1.4. Searching for the Correct Model with Information Criteria

A model selection criterion can only be calculated for a model that has been fit, which is not a problem for small p where all of the possible combinations of the p predictors can be quickly fit and optimized over (a process called “best-subsets” selection). However, in some cases (when p is large and/or when n is massive), it becomes impractical to fit all candidate models, the number of which increases on the order of $2^p$. In lieu of fitting every possible model, one can use a stepwise approach. In forward stepwise (greedy) selection, the starting point is a null model. In the first step ($k=1$), univariate models are fit for each feature and selected from using an information criterion (IC). Then at each subsequent step k, models with the remaining $p-k+1$ features are fit, their IC values computed, and the best one selected. The process stops when all of the newly fit models increase the IC. Backwards stepwise selection is similar, but its starting point is the full model and covariates are taken out one-at-a-time until the IC stops improving. Finally, a forwards and backwards approach consists of checking not only for the addition of new features, but also the removal of features that had been included in a previous step.

Stepwise approaches are considerably faster than best-subsets, but they still take a lot of computation time in high-dimensional or massive n settings (this will become evident in our second application). In these settings, the Least Absolute Shrinkage and Selection Operator (the lasso) [13] becomes especially useful. The lasso simultaneously estimates coefficients for each of the p features and selects from them, such that they are either “active” (i.e., ${\widehat{\mathit{\beta }}}_j\ne 0$), or inactive ${\widehat{\mathit{\beta }}}_j=0$. The estimated nonzero coefficients do suffer from a bias that is introduced by the lasso’s penalty term $\lambda$, but this bias is often warranted as it significantly cuts down on the variance of having too saturated of a model.

With features $\left[{\mathit{x}}_1,{\mathit{x}}_2,...,{\mathit{x}}_p\right]={\mathit{X}}_{n\times p}$, and a centered response variable $\mathit{y}$, where some (but not all) features are related to $\mathit{y}$, the lasso solution can be obtained by minimizing the following expression with respect to $\mathit{\beta }$:

$${\mathit{y}-\mathit{X}\mathit{\beta }}^2+\lambda {\displaystyle \sum _{j=1}^p}\left|{\mathit{\beta }}_j\right|.$$

The lasso is not a final solution to choosing an optimal model; a model must be selected along the coefficient path: a sequence of $\lambda$ values that fall upon a one-dimensional grid. In other words, the tuning parameter that dictates how much to shrink the nonzero coefficients, $\lambda$, must be selected somehow. This is where AIC, BIC, and CV come back into play; before this point in the process, the lasso search has been completely informed by goodness-of-fit (subject to shrinkage). In this context, all of the same properties related to model saturation hold for these criteria as was discussed in the previous section.

It is important to distinguish the difference between model search methods (best-subsets, stepwise approaches, and the lasso), and model selection criteria. The former is used to search the model space, while the latter is used to choose the optimal model from a set of candidates. Both are necessary, although currently, the only available tools for ranked sparsity focus on the search process. Earlier work [14] introduced the Sparsity Ranked Lasso, but this only addresses the search aspect of the model selection process. Here we take up the question of incorporating ranked parsimony into the selection process as well.

1.5. Simulation Illustration – Model Selection Criteria and Sparsity Levels

1.5.1. Simulation Setup

In this section, we illustrate the finite-sample behaviors of AIC (as an exemplary efficient criterion) and BIC (as an exemplary consistent criterion) for various sparsity/saturation levels in a simulation. Specifically, for a fixed n:

1)

A higher “true” saturation level requires a criterion with a lower penalty, with

a)

BIC outperforming AIC in a setting with a low saturation level, and

b)

AIC outperforming BIC in a setting with a high saturation level.

2)

As the number of candidate variables p increases, the performance of AIC deteriorates.

3)

Cross-validation (CV) can circumvent this issue by selecting the optimal penalty based on extra-sample performance.

We set $n=250$, $p\in$ {125, 175, 225, 245}, and the number of signal variables $s|p=p\times (.01,.02,.04,.1,.25,.5,.75,.95,.99)$ (i.e., we vary the saturation level $s/p$). For each combination of these, we generate p covariates independently from the standard normal distribution, then generate $\mathit{y}=\mathit{X}\mathit{\beta }+N(\mathbf{0},{\sigma }^2=5^2)$, where ${\mathit{\beta }}^T=({\mathbf{1}}_s^T,{\mathbf{0}}_{p-s}^T)$. There are s signal variables and $p-s$ null variables for each simulation. Since the number of non-zero coefficients is changing across levels of s, the signal-to-noise ratio ($\mathrm{SNR}={\mathit{\beta }}^{T}Var\left(\mathit{X}\right)\mathit{\beta }/{\sigma }^2$) changes as well for different sparsity levels (Table 1). Sparse settings are typically settings where the SNR is low, whereas truly saturated settings typically indicate a higher SNR, so we believe this reflects a useful shift. However, we also investigate what happens when the SNR is fixed. We henceforth denote the saturation level as $\omega$, where $\omega =s/p$.

We use the lasso search method and select the tuning parameter $\lambda$ according to AIC and BIC to compare their performance, relative to CV as a baseline. We utilize the lasso as a model selection tool (as opposed to a stepwise procedure) for its computational efficiency. We also generate $10*n$ new observations and responses, calculating each model’s root-mean-squared prediction error (RMSPE), defined as $\sqrt{N^{-1}{\sum }_{i=1}^N{(y_i-{\widehat{y}}_i)}^2}$ over the $N=10n$ held-out observations, as a gold-standard for comparison. While it does not directly consider whether the model has the right covariates or not, model misspecification (either under- or over-specification) will still inflate the RMSPE in cases where the model is badly specified. To showcase our points (1), (2), and (3) above, we plot the mean RMSPE across 1000 simulations by the saturation level in a series of plots.

1.5.2. Results

In Figure 1, we observe that as expected, there is a large discrepancy in the predictive quality of the model fits selected by AIC and BIC depending on the saturation level, but that this manifests differently for different n and p combinations. When $p=125$ (top-left), AIC and BIC performed relatively similarly for sparse settings, but as the saturation level increased, AIC outperformed BIC. As $p\to n$ however, the performance of AIC deteriorates to the point where BIC greatly outperforms AIC in sparse settings, and there is little, if any, benefit in using AIC even in saturated settings. In all cases, (10-fold) cross-validation performed optimally at selecting the optimal $\lambda$ value. Thus, we have shown in this example that the optimal penalty depends, to an extent at least, on the sparsity level (which we typically cannot know ahead of time).

Table 1. Design of the Section 1.5 saturation-level simulation. Rows index the true saturation level $\omega =s/p$ and columns the number of candidate predictors p, with the sample size fixed at $n=250$. The upper block reports the number of true signal variables s, and the lower block the resulting effective signal-to-noise ratio (SNR); sparser settings (smaller $\omega$) correspond to lower SNR. These configurations are used to compare AIC, BIC, and cross-validation for tuning lasso models.

	$p=125$	$p=175$	$p=225$	$p=245$
Number of Signal Variables
$\omega =0.01$	1	1	2	2
$\omega =0.02$	2	3	4	4
$\omega =0.04$	5	7	9	9
$\omega =0.1$	12	17	22	24
$\omega =0.25$	31	43	56	61
$\omega =0.5$	62	87	112	122
$\omega =0.75$	93	131	168	183
$\omega =0.95$	118	166	213	232
$\omega =0.99$	123	173	222	242
Effective SNR
$\omega =0.01$	0.04	0.04	0.08	0.08
$\omega =0.02$	0.08	0.12	0.16	0.16
$\omega =0.04$	0.20	0.28	0.36	0.36
$\omega =0.1$	0.48	0.68	0.88	0.96
$\omega =0.25$	1.24	1.72	2.24	2.44
$\omega =0.5$	2.48	3.48	4.48	4.88
$\omega =0.75$	3.72	5.24	6.72	7.32
$\omega =0.95$	4.72	6.64	8.52	9.28
$\omega =0.99$	4.92	6.92	8.88	9.68

Table 1. Design of the Section 1.5 saturation-level simulation. Rows index the true saturation level $\omega =s/p$ and columns the number of candidate predictors p, with the sample size fixed at $n=250$. The upper block reports the number of true signal variables s, and the lower block the resulting effective signal-to-noise ratio (SNR); sparser settings (smaller $\omega$) correspond to lower SNR. These configurations are used to compare AIC, BIC, and cross-validation for tuning lasso models.

	$p=125$	$p=175$	$p=225$	$p=245$
Number of Signal Variables
$\omega =0.01$	1	1	2	2
$\omega =0.02$	2	3	4	4
$\omega =0.04$	5	7	9	9
$\omega =0.1$	12	17	22	24
$\omega =0.25$	31	43	56	61
$\omega =0.5$	62	87	112	122
$\omega =0.75$	93	131	168	183
$\omega =0.95$	118	166	213	232
$\omega =0.99$	123	173	222	242
Effective SNR
$\omega =0.01$	0.04	0.04	0.08	0.08
$\omega =0.02$	0.08	0.12	0.16	0.16
$\omega =0.04$	0.20	0.28	0.36	0.36
$\omega =0.1$	0.48	0.68	0.88	0.96
$\omega =0.25$	1.24	1.72	2.24	2.44
$\omega =0.5$	2.48	3.48	4.48	4.88
$\omega =0.75$	3.72	5.24	6.72	7.32
$\omega =0.95$	4.72	6.64	8.52	9.28
$\omega =0.99$	4.92	6.92	8.88	9.68

A similar simulation was performed that holds the effective SNR fixed at 0.5 across all saturation levels. In this setting, very low saturation corresponds to a small number of large, easily detected coefficients, while higher saturation distributes the same total signal across many small coefficients near the noise floor (i.e., small enough that the signal is difficult to distinguish from sampling variability). These supplemental simulations again showed AIC performing increasingly poorly as p approached n. However, once saturation exceeded roughly 0.1, BIC remained surprisingly competitive with AIC: with individual coefficients near the noise floor, AIC’s looser penalty admitted weak signals that contributed little to prediction, while BIC’s stronger penalty trimmed to the detectable subset. As a result, the criterion-to-criterion contrast was less pronounced than in the main simulation.

2. Ranked Sparsity — Motivation

In the previous section, we illustrated how AIC and BIC perform very differently for different sparsity levels when the sample size is held constant, and that this difference increases with the number of candidate predictors. The performance of a penalty is highly dependent on the sparsity level. Now, we turn our attention to cases where we anticipate there being a different sparsity level across all of the covariates.

2.1. For Covariate Groups of Different Sizes

Say we have K covariate sets of varying sizes $p_k$. These covariate groups might include a small group of clinical features, a large group of genomic features, another even larger group of methylation markers, as was the case in VanHawkins et al. [15].

When covariates come from different groups or modalities, we typically expect different levels of sparsity for each group. In such cases, it is advisable to use a technique like the IPF-lasso [16] to perform cross-validation separately across all groups. However, this process is computationally expensive and requires $K-1$ additional tuning parameters to properly penalize each group according to its saturation level. An alternative approach is to use ranked sparsity methods such as the sparsity-ranked lasso [14]. The SRL uses a Bayesian motivation to naturally penalize groups of covariates according to their dimension size; in VanHawkins et al. [15] and Peterson [17], the SRL performed as well as IPF-lasso without the need for additional group-specific tuning.

Specifically, Bayesian interpretation of the lasso solution involves characterizing each ${\mathit{\beta }}_j$ as having independent and identical prior double exponential (Laplace) distributions, such that their joint prior density is

$$\pi \left(\mathit{\beta }\right)={\displaystyle \prod _{j=1}^p}{\displaystyle \frac{\lambda }{2}}{e}^{-\lambda |{\mathit{\beta }}_j|}.$$

Peterson and Cavanaugh [14] showed that when the covariates are part of distinct groups of disparate sizes $p_k$, we can modify this prior to have group-specific ${\lambda }_k$ values, such that the joint prior density becomes

$$\pi \left(\mathit{\beta }\right)={\displaystyle \prod _{k=1}^K}{\displaystyle \prod _{j=1}^{p_k}}{\displaystyle \frac{{\lambda }_k}{2}}{e}^{-{\lambda }_k\left|{\mathit{\beta }}_j\right|}.$$

When ${\lambda }_k=\lambda \sqrt{p_k}$ (which can be readily set up using penalty weights), larger groups are naturally penalized more heavily than smaller groups so that the prior information contributed to the posterior by each group is held constant across groups. This property is desirable when we expect larger groups to be more sparse, as it ties the penalization directly to the group dimension. However, this method is only applicable in penalized regression settings; additional tools are required for this purpose in other model selection settings (e.g. best-subsets or stepwise selection). We will return to this in Section 3.2.

2.2. For Interactions

Take now the context of selecting from all possible pairwise interactions. We define “main effects” as linear coefficients in a model, and “interaction effects” as coefficients on the product of covariates that correspond to the main effects. In this section, we will show that the maximum saturation level attainable for the interaction effects is necessarily limited by hierarchy assumptions about the true generating model. Model hierarchy is typically broken down into “strong”, e.g. $y_i={\beta }_0+{\beta }_1{x}_{1i}+{\beta }_2{x}_{2i}+{\beta }_3{x}_{1i}{x}_{2i}$, “weak”, e.g. $y_i={\beta }_0+{\beta }_1{x}_{1i}+{\beta }_3{x}_{1i}{x}_{2i}$, and “anti” hierarchical models, e.g. $y_i={\beta }_0+{\beta }_3{x}_{1i}{x}_{2i}$.

Let $s_M$ again refer to the number of signals in the main effects of $p_M$ possible main effect predictors, and let $s_I$ refer to the number of signals in the interactions of $p_I=\left({\displaystyle \genfrac{}{}{0pt}{}{p_M}{2}}\right)$ possible interactions. Also, let ${\omega }_M$ and ${\omega }_I$ refer to the saturation levels in the main effects and the interactions, respectively. Assume that $p_M\ge 2$ and that $p_M{\omega }_M\in \mathbb{Z}$. By definition, we have that $s_M=p_M{\omega }_M$ and $s_I=\left({\displaystyle \genfrac{}{}{0pt}{}{p_M}{2}}\right){\omega }_I$, and therefore

$${\omega }_I={\displaystyle \frac{s_I}{\left(\genfrac{}{}{0pt}{}{p_M}{2}\right)}}\left(1\right).$$

Under strong hierarchy, the signals in the interactions must be combinations of the signals in the main effects, which implies

$$\begin{array}{c}\hfill s_I\le \left({\displaystyle \genfrac{}{}{0pt}{}{p_M{\omega }_M}{2}}\right)\left(2\right).\end{array}$$

Together, (1) and (2) imply

$$max{\omega }_I={\displaystyle \frac{\left(\genfrac{}{}{0pt}{}{p_M{\omega }_M}{2}\right)}{\left(\genfrac{}{}{0pt}{}{p_M}{2}\right)}}={\displaystyle \frac{{\omega }_M(p_M{\omega }_M-1)}{p_M-1}}.$$

This upper bound is plotted as a function of ${\omega }_M$ in Figure 2 (left), and it is evident from the figure that $max{\omega }_I<{\omega }_M\forall {\omega }_M\in (0,1)$.2 Under weak hierarchy, instead of the expression in (2), we have

$$s_I\le {\displaystyle \sum _{j=1}^{s_M}}{p}_M-j={\displaystyle \frac{1}{2}}{p}_M{\omega }_M(2p_M-p_M{\omega }_M-1).$$

The corresponding bound on the saturation level of the interactions then becomes

$$\begin{array}{cc}\hfill max{\omega }_I& ={\displaystyle \frac{\frac{1}{2}{p}_M{\omega }_M(2p_M-p_M{\omega }_M-1)}{\left(\genfrac{}{}{0pt}{}{p_M}{2}\right)}}\hfill \\ & ={\displaystyle \frac{{\omega }_M(2p_M-p_M{\omega }_M-1)}{p_M-1}}.\hfill \end{array}$$

We see that in the weakly hierarchical case, the maximum saturation level in the interactions is still bounded by the signals in the main effects, though it is not restricted to be less than the saturation level in the main effects (Figure 2, right). Under the heredity principle of Chipman [18], strongly hierarchical structures are the most probable a priori, weakly hierarchical ones less so, and anti-hierarchical ones the least. The bounds above trace this same ordering: under strong hierarchy the interaction saturation is provably below the main-effect saturation; under weak hierarchy it remains tied to the main-effect signals, though no longer strictly below them; and under anti-hierarchy the main-effect constraint disappears altogether. These derivations therefore suggest that one should not apply the same penalty to both main effects and interactions and expect optimal performance: because ranked sparsity concerns the expected, a priori interaction saturation, and the prior mass arguably concentrates on the more hierarchical regimes, we should still expect interactions to be sparser than main effects even though any particular generating model need not obey this.

Crucially, a heavier penalty like those we propose for ranked sparsity raises the evidentiary bar for interactions and higher-order terms, but does not forbid them. A weakly or anti-hierarchical truth can therefore still be recovered when the data warrant it.

2.3. Conclusions

We have shown that the optimal penalty for any model selection criterion is highly dependent on the true sparsity level. We have also shown that when there are different groups of covariates with different expected sparsity levels, as is especially true when main effects are considered alongside interactions, it is important to account for this in the model selection process. In the next section, we explore how to do this with a novel ranked sparsity method, RBIC.

3. Ranked Sparsity via RBIC

3.1. The Extended Bayesian Information Criterion

3.1.1. Description & Motivation

The Extended Bayesian Information Criterion (EBIC) of Chen and Chen [7] extends BIC to cases in which the candidate pool p is large relative to the sample size. One of the standard conditions for BIC’s consistency is that p is fixed; EBIC adds a term that grows with p to restore consistency when p increases with n, and in particular when $p>\sqrt{n}$, where BIC is otherwise inconsistent. By placing a uniform prior on all candidate models, BIC induces a strongly nonuniform prior on model size, since there are $\left({\displaystyle \genfrac{}{}{0pt}{}{p}{m}}\right)$ candidate models of size m; this combinatorial factor concentrates sharply near $m=p/2$ as p grows. BIC therefore encodes an implicit prior that the saturation level $\omega =m/p$ is near 0.5, which is rarely a reasonable prior belief. We refer the reader to Peterson and Cavanaugh [19] for visualizations of this prior and additional background.

EBIC corrects this prior imbalance by penalizing models in proportion to the size of their model class. Letting m denote the dimension of a candidate model with parameter vector $\mathit{\beta }$,

$$\begin{array}{cc}\hfill \mathrm{EBIC}& =-2logl\left(\widehat{\mathit{\beta }}\right)+mlogn+2\gamma log\left({\displaystyle \genfrac{}{}{0pt}{}{p}{m}}\right)\hfill \\ & =\mathrm{BIC}+2\gamma log\left({\displaystyle \genfrac{}{}{0pt}{}{p}{m}}\right).\hfill \end{array}$$

The tuning parameter $\gamma$ controls the strength of EBIC’s correction; Peterson and Cavanaugh [19] visualize how $\gamma$ shapes the implied prior on model size. Three cases are useful. If $\gamma =0$, EBIC reduces to the original BIC. If $\gamma =1$, EBIC places a uniform prior on the marginal distribution of model size and is consistent whenever p grows polynomially with n ($p=O\left(n^{\kappa }\right)$ for any $\kappa \ge 0$ not depending on n). Intermediate values trade between these extremes; for example, $\gamma =0.5$ retains consistency for $\kappa <1$. Setting $\gamma =1$ can be too stringent in practice, so Chen and Chen [7] propose a data-adaptive choice that targets consistency without being unnecessarily aggressive: solving $\gamma =1-1/\left(2\kappa \right)$ from $p=n^{\kappa }$ and clamping at zero,

$${\gamma }_{\mathrm{EBIC}}=max\left({\displaystyle \frac{log(p/\sqrt{n})}{logp}},0\right),$$

which equals zero when $p\le \sqrt{n}$ (the regime in which BIC is already consistent) and rises smoothly with p. The corresponding additional penalty over BIC is $2{\gamma }_{\mathrm{EBIC}}log\left({\displaystyle \genfrac{}{}{0pt}{}{p}{m}}\right)$. Figure 3 shows how this term varies with p, n, and m: it grows with p, decays with n (eventually to zero), and is largest near $m=p/2$, which means that for high-p scenarios the penalty is effectively increasing with m for sparse models.

An intuitive framing of EBIC is to consider each candidate covariate as a ball in a single urn of size p, where any candidate model is a random draw of m balls, giving $\left({\displaystyle \genfrac{}{}{0pt}{}{p}{m}}\right)$ possible models of that size. In the next section, we extend this one step further to the setting of K groups of covariates.

3.2. Ranked Sparsity Extension for BIC

3.2.1. Description & Motivation

We use this same “balls and urns” metaphor to intuitively motivate a ranked sparsity extension of BIC (RBIC). Instead of one urn with p covariates, we postulate that we have multiple (K) urns, each with $p_k$ covariates. Then, for any candidate model with $m_k$ parameters of $p_k$ possible for all $k\in \{1,2,...,K\}$,

$$\begin{array}{cc}\hfill \mathrm{RBIC}& =-2logl\left(\widehat{\mathit{\beta }}\right)+{\displaystyle \sum _{k=1}^K}{m}_{k}logn+{\displaystyle \sum _{k=1}^K}2{\gamma }_{k}log\left({\displaystyle \genfrac{}{}{0pt}{}{p_k}{m_k}}\right)\hfill \\ & =\mathrm{BIC}+{\displaystyle \sum _{k=1}^K}2{\gamma }_{k}log\left({\displaystyle \genfrac{}{}{0pt}{}{p_k}{m_k}}\right),0\le {\gamma }_k\le 1\forall k.\hfill \end{array}$$

This formulation of RBIC effectively penalizes a covariate’s coefficient differently depending on which covariate subspace (or urn) it comes from. In the case of pairwise interactions, the penalty (when $p>2$) is higher for interaction effects than for main effects.

By default, we suggest setting ${\gamma }_k={\gamma }_{\mathrm{EBIC}}$ to benefit from its consistency property for high-dimensional data:

$$\mathrm{RBIC}=\mathrm{BIC}+2{\gamma }_{\mathrm{EBIC}}{\displaystyle \sum _{k=1}^K}log\left({\displaystyle \genfrac{}{}{0pt}{}{p_k}{m_k}}\right).$$

Setting ${\gamma }_k=0\forall k$ yields the original BIC, and if $K=1$, RBIC is identical to EBIC. In Application 1, we will sweep ${\gamma }_k$ from 0 to ${\gamma }_{\mathrm{EBIC}}$ to illustrate how its performance changes.

3.2.2. Implementation

For any likelihood-based model with well-specified covariate groups, RBIC is amenable to model selection via best-subsets or any stepwise selection framework. However, several advancements in the computational efficiency of information-criterion-based model searching, such as those implemented by the R packages leaps [20] and bestglm [21], are unavailable. These packages take advantage of the fact that all models with m parameters have the exact same penalty ($mlogn$ for BIC, or $2m$ for AIC). With RBIC, each covariate may be penalized differently, so not all models with m parameters have the same penalty.

We have developed the R package ICtoolkit [22] to provide a unified interface for computing AIC, AICc, BIC, EBIC, and RBIC. The package supports models of class lm, glm, glmnet [23], and ncvreg [24]. For lm and glm objects, each criterion returns a scalar value. For penalized regression objects (glmnet and ncvreg), the criteria are computed across the full regularization path, returning a vector of values (one per $\lambda$), enabling information-criterion-based tuning parameter selection as an alternative to cross-validation.

The RBIC implementation in ICtoolkit accepts a user-specified grouping structure via a named list P_index, where each element is a character vector of variable names belonging to that group. For example, if main effects belong to one group and pairwise interactions to another, the user specifies this as P_index = list(main = c("x1", "x2", ...), interactions = c("x1:x2", ...)); this can be automated with the expand_scope function. The package then computes the group-specific combinatorial penalties $log\left({\displaystyle \genfrac{}{}{0pt}{}{p_k}{m_k}}\right)$ for each group k and applies the appropriate $\gamma$ value. The default $\gamma$ is the data-adaptive rule of Chen and Chen [7], though the package also supports per-group $\gamma$ values, fixed numeric values, and user-supplied functions.

For stepwise model selection, ICtoolkit provides an ic_step() function that extends the approach of MASS::stepAIC [25] to support all five criteria. The function supports forward, backward, and bidirectional search for lm and glm objects. When the criterion is set to AIC or BIC, ic_step() produces identical model selections to MASS::stepAIC. When set to EBIC or RBIC, the function applies the appropriate extended penalties at each candidate evaluation step.

As described in in Section 1.4, we can also use information criteria to select the overall tuning parameter $\lambda$ in a lasso-based selection framework as long as covariate groups are specified (in this case we will use BIC, EBIC, and RBIC). For penalized models, ICtoolkit computes the criterion at every point along the regularization path, and the optimal $\lambda$ is selected as the value that minimizes the chosen criterion.

3.2.3. Simulation Study

In order to study how RBIC compares to its competitors, we simulate data under the Friedman generating model [26], shown below.

$$y=10sin\left(\pi x_1{x}_2\right)+20{(x_3-0.5)}^2+10x_4+5x_5+N(0,3^2)$$

In addition to these 5 “signal” terms, 5 noise variables were generated as well (all independent uniform random variables). For 1500 Monte Carlo simulations, models were fit to $n=80$ observations using forwards-and-backwards stepwise regression using BIC, EBIC, and RBIC.

We investigate how RBIC performs in the linear regression setting (no candidate polynomials or interactions), the quadratic regression setting (squared terms and pairwise interactions), and the cubic regression setting (cubed terms, three-way interactions, and combinations of quadratic and main effect terms). Note that in the linear regression setting, $K=1$, so RBIC and EBIC are identical. To quantify model performance, we calculate the RMSPE on 10,000 newly generated observations, as well as the number of false negatives (NFN), the number of false positives (NFP), and the coefficient-wise selection percentage. Based on the form of the Friedman model, the covariates with legitimate nonzero coefficients are $x_1$, $x_2$, $x_4$, $x_5$, $x_1^2$, $x_2^2$, $x_3^2$, and the pairwise interaction $x_1\ast x_2$; all other covariates, if selected, were considered false positives.

Table 2 reports the means and standard deviations of these performance metrics across simulations, and Table 3 shows the per-variable selection percentages. As expected, BIC works well only in the linear regression setting, but fails badly in the quadratic or cubic regression setting (in terms of prediction and the number of false positives). BIC’s poor performance is due to the violation in asymptotic assumptions as the number of candidate predictors is higher than the sample size for these settings. EBIC and RBIC both control the number of false positives and enjoy similarly better prediction performance, though RBIC saw a small but consistent gain in predictive efficacy relative to EBIC for both quadratic and cubic regression. This improvement is accompanied by a slightly lower number of false negatives, but with little change in the number of false positives between EBIC and RBIC. All methods suffered from higher RMSPE when considering spurious cubic terms, with RBIC suffering comparably less than EBIC and BIC.

In summary, RBIC’s ranked sparsity penalty provides an improvement in predictive performance relative to EBIC when applied to interactions, especially as the candidate pool becomes more complex and higher-dimensional (i.e., with more interactions). In the next section, we apply RBIC to a real data set to illustrate how it can be used for interaction selection in practice.

4. Application

4.1. Application 1: Interaction Selection in COVID-Era ICU Caregiver Data

4.1.1. Description

We illustrate RBIC’s utility for interaction selection using data from Amass et al. [27], who surveyed family members of patients admitted to the ICU during the COVID-19 pandemic ($n=298$). The outcome is the Impact of Event Scale (IES) total score, a continuous measure of post-traumatic stress. We consider 10 candidate main-effect predictors: four continuous variables (age, SOFA score, Charlson comorbidity index, ICU length of stay) and six categorical variables (sex, ethnicity, race, education, patient death, and psychotropic medication use). Race and education are multi-level factors contributing 3 degrees of freedom each; all other categorical predictors are binary.

To test whether RBIC appropriately controls the inclusion of spurious complex terms, we compare forward stepwise selection under two candidate scopes:

• 2-way scope: all main effects, quadratic polynomials of the continuous variables, and all $\left({\displaystyle \genfrac{}{}{0pt}{}{10}{2}}\right)=45$ pairwise interactions (59 terms, 103 df total).
• 3-way scope: the 2-way scope augmented by all $\left({\displaystyle \genfrac{}{}{0pt}{}{10}{3}}\right)=120$ three-way interactions (179 terms, 399 df total).

The key question is whether expanding the candidate pool from 2-way to 3-way interactions causes complex terms to proliferate without improving model performance.

4.1.2. Group Structure for RBIC

Candidate terms are grouped by their structural complexity, with group sizes reflecting degrees of freedom rather than term counts (since multi-level factor terms consume more than one degree of freedom each). Under the 2-way scope, $K=2$ groups are defined: main effects ($p_1=14$ df) and level-1 terms comprising pairwise interactions and quadratic polynomials ($p_2=89$ df). Under the 3-way scope, a third group is added for three-way interactions ($p_3=296$ df). Since the level-1 and level-2 pools are substantially larger than the main-effects pool, RBIC’s combinatorial penalty $log\left({\displaystyle \genfrac{}{}{0pt}{}{p_k}{m_k}}\right)$ is correspondingly steeper for interactions than for main effects, encoding the ranked sparsity prior that complex terms should face a higher bar for inclusion.

Because $\gamma$ controls the strength of RBIC’s penalty relative to BIC, with $\gamma =0$ recovering BIC and ${\gamma }_{\mathrm{EBIC}}\approx 0.385$, we report across $\gamma \in \{0,0.10,0.25,{\gamma }_{\mathrm{EBIC}}\}$, with the same $\gamma$ applied uniformly across groups. The $\gamma =0$ case (denoted RBIC.00) serves as a reminder that RBIC reduces to BIC; the intermediate values RBIC.10 and RBIC.25 interpolate between BIC and the default RBIC, allowing us to characterize how RBIC’s performance changes as the analyst dials the strength of the ranked-sparsity penalty.

4.1.3. Results

Table 4 presents the average model characteristics and held-out predictive accuracy across 50 random 2/3-1/3 train/test splits, for each criterion under both candidate scopes. First, we notice that AIC overfits substantially in both scopes, most catastrophically under the 3-way scope where it selects a three-way interaction on every split and incurs a mean test RMSPE of 67.9 (compared to roughly 5.8 for the BIC-family criteria). Second, among the BIC-family criteria, as $\gamma$ increases from 0 (BIC) through .10, .25, and to the default of ≈ 0.385, the mean degrees of freedom, the percentage of splits selecting a three-way interaction, and the model’s permissiveness toward complex terms all decrease monotonically. RBIC at $\gamma =0.10$ achieves the lowest mean test RMSPE in both scopes (5.802 in the 2-way scope; 5.822 in the 3-way scope) while selecting fewer parameters than BIC and roughly halving the spurious three-way selection rate (8% vs. 16% under the 3-way scope). Finally, under the 2-way scope, where there is no spurious three-way risk to protect against, RBIC.25, the default RBIC, and EBIC all over-shrink and incur higher mean test RMSPE than BIC. Under the 3-way scope, those same heavier penalties suppress three-way selection almost entirely (2%, 0%, and 0% respectively, vs. BIC’s 16%).

The $\gamma$ parameter can therefore function as an analyst-tunable knob, with $\gamma \approx 0.10$ in this case retaining most of BIC’s permissiveness on main effects while halving the spurious three-way rate, and the default ${\gamma }_{\mathrm{EBIC}}$ providing the strongest protection against over-fitting to complex higher-order interactions (at the cost of marginally worse prediction).

4.2. Application 2: Multimodal Survival Modeling in Lung Adenocarcinoma

4.2.1. Description

Whereas Application 1 examined ranked sparsity along the complexity axis (main effects vs. interactions), this application demonstrates RBIC on the modality axis: two covariate groups of vastly different dimension that arise from distinct data sources. We use the lung adenocarcinoma data of Shedden et al. [28], a multi-site study of $n=442$ patients with gene expression profiling and clinical annotation. The outcome is overall survival (time to death or censoring; 236 events). The candidate predictors comprise two groups:

• Clinical ($p_1=7$ df): age (centered and scaled), sex, adjuvant chemotherapy status, smoking history (never vs. ever), tumor grade (3 levels, 2 df), and race (White vs. non-White). Missing values are imputed via k-nearest neighbors; categorical predictors are dummy-coded.
• Genetic ($p_2=22,283$): Affymetrix HG-U133A microarray probe-set expression levels, centered and scaled to unit variance.

The total candidate pool is $P=p_1+p_2=22,290$. The extreme asymmetry ($p_2/p_1\approx 3,183$) makes this a prototypical multimodal setting in which standard model selection criteria, which assume covariate equipoise, will be overwhelmed by the genetic group. We compare model selection via two search strategies (lasso and Sparsity-Ranked Lasso) crossed with four selectors (BIC, EBIC, RBIC, and cross-validation), and via forward stepwise selection under BIC, EBIC, and RBIC.

4.2.2. Group Structure Construction and the Sparsity-Ranked Lasso

For RBIC, we define $K=2$ groups whose sizes match the number of candidate predictors in each modality: $p_1=7$ (clinical) and $p_2=22,283$ (genetic). Because $p_2\gg p_1$, the combinatorial penalty $log\left({\displaystyle \genfrac{}{}{0pt}{}{p_k}{m_k}}\right)$ is substantially steeper for genetic predictors than for clinical predictors at every model size, encoding the ranked sparsity prior that clinical covariates (which are fewer and typically well-motivated) should face a lower bar for inclusion.

For the regularization-path approaches, we fit penalized Cox proportional hazards models using the lasso penalty via ncvsurv. The Sparsity-Ranked Lasso (SRL) modifies the standard lasso by assigning group-specific penalty factors ${\lambda }_g=\lambda ·\sqrt{p_g}$, where g denotes the group to which predictor j belongs. Under this scheme, the $p_1=7$ clinical predictors each receive a penalty factor of $\sqrt{7}\approx 2.6$, while the $p_2=22,283$ genetic predictors receive a factor of $\sqrt{22,283}\approx 149.3$, giving clinical variables roughly a 56-fold advantage in the penalization scale. Combining the SRL search with RBIC selection along the regularization path brings ranked sparsity into both stages of the model selection process: the search and the selection.

4.2.3. Results: Regularization Paths

Table 5 presents the (Search × Selector) grid: Lasso and SRL each crossed with BIC, EBIC, RBIC, and 10-fold CV. When using the standard lasso penalty, the information criteria-based selectors produced highly sparse models. BIC selected a one-parameter model containing a single genetic covariate, while EBIC and RBIC both produced null models. CV-based selection, which does not penalize for model complexity, produced a much larger model with 29 covariates (28 genetic, 1 other).

The use of SRL, which applies differential penalty factors according to data modality (genetic vs. clinical), produces a slightly different pattern of selection. Both BIC and RBIC yield the same small $(p=3)$, interpretable model consisting entirely of clinical covariates. EBIC produces a null model, indicating that it remains overly conservative in its penalization relative to BIC and RBIC, even when paired with SRL. CV-based selection under SRL yields a more moderate model than its lasso counterpart, with a dramatic reduction in the total number of selected covariates alongside an increase in clinical representation (all 7 clinical covariates selected compared to 1 under the standard lasso).

These results illustrate a key limitation of the standard lasso in this high-dimensional multimodal setting: when the clinical and genetic predictors are penalized equally, the sheer scale of the genetic predictor space renders individual clinical covariates nearly impossible to justify for inclusion under any information criterion. SRL addresses this directly by incorporating modality-specific penalty factors, allowing clinically relevant predictors to compete for inclusion without being crowded out by the much larger genetic predictor pool.

4.2.4. Results: IC Path Figure

Figure 4 displays BIC, EBIC, and RBIC along the regularization path for both SRL and standard lasso, starting on the left from the null model. These plots show that these methods are functioning appropriately in that they are selecting covariates to be added into models as penalization decreases, which steadily increases the IC values presented on the y-axis.

Of note, BIC values appear to reach a peak and then steadily decrease as more parameters are introduced into the model, indicating that BIC does not sufficiently penalize model complexity in this high-dimensional setting. This is a well-documented limitation of standard BIC when the number of candidate predictors is large relative to sample size. In contrast, EBIC and RBIC maintain a more stable plateau or continue to increase across the regularization path, suggesting they are better suited for variable selection in this context. EBIC and RBIC exhibit qualitatively similar IC path behavior and all three criteria achieve their minima in a relatively sparse region of the path before complexity penalties begin to dominate.

4.2.5. Results: Forward Stepwise Selection

Table 6 summarizes the size and composition of each criterion’s stepwise-selected model. Forward stepwise selection under BIC resulted in a notably oversized model. By step 143, BIC had selected 141 genetic predictors and 2 clinical predictors before failing at step 144, likely due to near-multicollinearity among the large set of included genetic variables. This outcome is not unexpected. BIC does not account for the size of the search space and proceeds in a greedy, sequential manner. As the model grows large relative to the sample size, the asymptotic penalty underlying BIC becomes increasingly unreliable. The computational cost reflects this as BIC required approximately 14 hours on 63 cores to reach this point. In contrast, EBIC and RBIC each completed selection in approximately 6 minutes, terminating after just 4 and 5 steps, respectively.

EBIC selected 4 variables in total: 1 clinical (age) and 3 genetic. RBIC selected 5 variables: 3 clinical (age, sex, and adjuvant chemotherapy) and 2 genetic. The key difference between the two criteria lies in the group-specific combinatorial penalty applied by RBIC, which imposes a heavier penalty on genetic predictors than on clinical ones, reflecting the prior assumption that genetic signal is sparse. As a result, RBIC is more permissive toward clinically meaningful predictors. Its selected clinical variables are well-established prognostic factors in this setting. Notably, EBIC and RBIC share two overlapping genetic markers, suggesting consistency in identifying the strongest genetic signals despite differences in penalization structure.

Comparing the forward stepwise results to those from the regularization path reveals consistency across methods, as presented in Table 7. Under BIC, the sparse regularization path selected age, sex, and adjuvant chemotherapy, while lasso selected a single genetic marker, FAM117A. Forward stepwise selection under BIC recovered age and adjuvant chemotherapy, along with a large number of genetic predictors (138 not shown in the table). Under EBIC, both the SRL and lasso regularization paths selected no variables, while forward stepwise selection identified age and three genetic markers (FAM117A, ZC2HC1A, and NDST1), reflecting EBIC’s conservative penalization. RBIC showed the greatest consistency across methods: the SRL path selected age, sex, and adjuvant chemotherapy; lasso selected no covariates; and forward stepwise selection recovered all three clinical variables alongside two genetic markers (FAM117A and ZC2HC1A). The clinical variables selected by BIC and RBIC with SRL are identical to those identified by RBIC forward selection, and FAM117A appears in every model that includes any genetic variables regardless of criterion or method, suggesting that FAM117A carries consistent prognostic signal across selection approaches. This is further supported by previous analysis of the lung adenocarcinoma data, which found that FAM117A had a protective effect on survival [14].

5. Discussion

In this work, we have argued that “covariate equipoise,” the implicit prior that every candidate predictor is equally worthy of inclusion in a model, is built into AIC, BIC, and the lasso, and that this assumption breaks down whenever the analyst already has structural information that some covariates should face a higher bar than others. RBIC operationalizes this structure in the information-criterion framework: by partitioning the candidate pool into K groups and applying a group-specific combinatorial penalty $log\left({\displaystyle \genfrac{}{}{0pt}{}{p_k}{m_k}}\right)$, RBIC requires increasingly stronger evidence for inclusion as a group’s pool grows. It is a natural cousin of the Sparsity-Ranked Lasso [14], which operates in the regularization framework, and of the IPF-lasso [16], which tunes group-specific ${\lambda }_k$ via cross-validation. Setting $K=1$ recovers Chen and Chen [7]’s EBIC; setting $\gamma =0$ recovers BIC. RBIC is therefore not a wholesale replacement for these criteria but a structured generalization that exposes assumptions previously hidden in their construction.

We have built an empirical case for that generalization across two simulation studies and two applications. Section 1.5 grounded the motivation: the optimal IC penalty depends sharply on the true saturation level, and no single criterion dominates. The Friedman simulation in Section 3.2.3 then showed that, when a polynomial/interaction expansion introduces candidates with a natural ranked-sparsity structure, RBIC consistently outperformed EBIC on prediction in both quadratic and cubic settings while preserving comparable false-positive control. The gain came specifically from the combinatorial term reflecting group structure. Application 1 traced this advantage on the complexity axis, with a $\gamma$ sweep showing that RBIC at $\gamma =0.10$ retained most of BIC’s flexibility on main effects while halving the spurious three-way interaction selection rate and achieving optimal out-of-sample predictions. Application 2 then demonstrated the modality axis, where the number of genetic covariates was over three orders of magnitude larger than the clinical pool. The most striking lessons concerned the interplay between search and selection. Forward stepwise under BIC blew up: 141 genetic predictors had entered the model before it failed at step 144, after roughly 14 hours of computation. RBIC and EBIC stepwise terminated within minutes at parsimonious models, with RBIC admitting three clinical predictors and two genetic predictors while EBIC admitted only one clinical predictor. The regularization path analyses revealed that the standard lasso, which applies the same penalty to all covariates, was overwhelmed by the genetic pool: BIC selected a single genetic predictor and EBIC and RBIC both produced null models. The SRL, which applies group-specific penalty factors throughout the model search path, allowed BIC and RBIC to select the same small model of three clinical predictors, while EBIC remained null.

As a practical recommendation, the default ${\gamma }_{\mathrm{EBIC}}$ is a conservative-but-safe starting point that protects against false discoveries from the largest groups. In settings where prediction matters more than complete suppression of complex terms, exploring other $\gamma$ may expose a useful middle ground.

RBIC has several limitations deserving attention. RBIC requires the group structure to be specified a priori; it operationalizes a prior conceptualization of a ranking, but does not discover one. The combinatorial penalty $log\left({\displaystyle \genfrac{}{}{0pt}{}{p_k}{m_k}}\right)$ is well-defined when groups partition the candidate pool, but overlapping or nested groupings (a polynomial term that could plausibly belong either to a “main effect” or “complexity-1” group) require the analyst to commit to one assignment, and that decision can shape downstream results. The consistency results of Chen and Chen [7] extend to RBIC only via the $K=1$ reduction; a formal proof for $K>1$, particularly when group sizes grow at different rates with n, remains open.

Finally, ranked sparsity must be incorporated into the search process, either via SRL or via RBIC-guided stepwise selection — pairing ordinary lasso with RBIC is ineffective, because the menu of potential models offered to RBIC is already skewed towards terms it finds unappealing. In settings where the number of potential models is very large, SRL is currently the faster option, as it benefits from highly efficient algorithms that have been developed for the lasso. Whether a stepwise RBIC implemetation could be competitive from a computational standpoint remains an open question.

These limitations suggest natural directions for future work. The most immediate is a data-driven choice of $\gamma$, either via cross-validation or an empirical-Bayes calibration; ICtoolkit already supports per-group ${\gamma }_k$, but the paper has not exercised that flexibility. A formal consistency theory for $K>1$ would be valuable, especially in regimes in which different groups grow at different polynomial rates. On the implementation side, the constant-penalty trick that makes leaps [20] and bestglm [21] efficient does not apply to RBIC, but a branch-and-bound with group-specific lower bounds appears tractable and would extend RBIC’s reach into best-subsets selection, building on efficient branch-and-bound best-subset implementations such as BranchGLM [29]. There are also clear conceptual neighbors worth exploring further. The modified BIC (mBIC) of Bogdan et al. [30] adds a multiplicity-aware term to BIC and arrives at a penalty closely related to EBIC’s combinatorial correction; RBIC’s group-specific motivation could be applied as a ranked extension of that same multiplicity adjustment to K pools of different size.

The broader contribution of this work is to make a hidden assumption visible. The implicit prior that every candidate covariate deserves the same admission penalty is rarely a justified scientific belief; more often it is an analytic convenience. When something is already known about how the candidate space is structured (interactions versus main effects, genetic versus clinical, derived versus measured), RBIC offers a modest tool for incorporating that knowledge explicitly. The cost is committing to a ranking; the benefit, as our applications and simulations have shown, is a shield against spurious complexity — one that lets the statistician fly as close to the sun as the data allow without burning their wings.