Modern data collection technology has made it possible for researchers in many different scientific domains to access high-dimensional and ultra-high-dimensional data at relatively low costs. Sparsity is common in high-dimensional data, especially in genomics, proteomics, biomedical imaging, tumor classification, image analysis, signal processing, and finance (Evans, 2006;Manolio and Collins, 2007;Kooperberg and Leblanc, 2008;Cordell et al., 2009). Variable selection plays a critical role in high-dimensional data modeling. Variable selection can increase estima-tion accuracy and model interpretability with sparsity by successfully identifying the subset of significant predictors. Over the last decades, numerous statistical and machine learning-based methods have been proposed for high-dimensional data, especially when the number of predictors (p) surpasses the number of observations (N ) (Fan and Li, 2006;Donoho, 2000;Fan and Lv, 2010).
Conventional methods such as stepwise selection, forward and backward selection, and subset selection combined with AIC and BIC were devised to address variable selection tasks successfully. However, these methods face computational challenges as the number of predictors (p) grows. For example, subset selection becomes exceptionally demanding when p exceeds 20 due to the potential creation of 2 p subset models. Additionally, forward and backward selection criteria encounter difficulties in backtracking variable addition or removal orders. Researchers have proposed different families of variable selection methods in the high-dimensional setting, some of which are penalty-driven algorithms. This type of penalty-driven method simultaneously achieves variable selection and parameter estimation in various types of regression models. For instance, the least absolute shrinkage selector operator (LASSO) performs variable selection via L 1 -norm penalization on the regression coefficients (Tibshirani, 1996). LASSO is inclined to set some coefficients to zero while promoting sparsity, which can lead to the exclusion of relevant variables, potentially overlooking important features. In situations involving correlated predictors or when the number of predictors surpasses the number of observations, LASSO may exhibit instability in variable selection. Other penalty-driven variable selection algorithms are Smoothed Clipped Absolute Deviance (SCAD) (Fan and Li, 2001) and Minimax Concave Penalty (MCP) (Zhang, 2010). Both algorithms incorporate a nonconvex penalty function to alleviate excessive penalization of large coefficient values in regression. The penalty-driven methods depend on selecting appropriate tuning parameter values, which can be non-trivial.
While the aforementioned variable selection methods have been successfully applied in regression modeling, they primarily focus on selecting main effects. When we consider two-way interactions between main effects, the dimension of predictors significantly increases from p to p + p(p -1)/1. Conventional variable selection methods may not be well-suited for this expanded dimensionality. However, penalty-driven variable selection methods may be more suitable for addressing interaction selection in high-dimensional data, e.g., additive LASSO (Hastie and Tibshirani, 1990). Bien et al. (2013) introduced LASSO with hierarchical interaction selection. Zhao et al. (2009) proposed a composite absolute penalty function and showed that it can perform two-way interaction selection. Yuan et al. (2009) proposed a structured two-way interaction selection and estimation procedure for parametric models. Choi et al. (2010) introduced the parametric interaction selection method under the heredity assumption. Hao et al. (2018) pointed out that these approaches perform well and exhibit oracle properties when the number of predictors is relatively small, typically a few hundred or less. However, when p is much greater than the sample size, these techniques become impractical due to the computational challenges associated with handling the entire O(p 2 ) × N design matrix and solving intricate constrained optimization problems. More recently, interaction screening methods have been proposed to deal with interaction selection in high dimensional data (Hao andZhang (2014, 2017); Kong et al. (2017)). While these methods are computationally efficient, their selection performance is slightly subpar due to forward selection and the marginal relationship between predictors and response. Subsequently, Hao et al. (2018) proposed the Regularization path Algorithm under the Marginality Principle (RAMP), designed to select main effects and interactions simultaneously under certain statistical principles defined in Section 2. The RAMP algorithm considers twoway hierarchical multiplicative interactions within penalized linear or logistic regression, which means that the method has some limitations for higher-order, nonlinear interactions.
Alternatively, tree-based algorithms have the potential to detect and select higher-order non-multiplicative interactions in high-dimensional data while sacrificing the selection process's interpretability. This tree-based approach enables the use of several ensemble methods for interaction selection and model construction, such as bagging (Breiman, 1996) and random subspace sampling (Tin Kam Ho, 1998;Breiman, 2001). As outlined by (Kotsiantis and Kanellopoulos, 2012), these ensemble techniques effectively achieve robust modeling in environments with high noise levels. The bagging technique involves training a multitude of models on different segments of the dataset, which significantly reduces variance and boosts the robustness of the overall ensemble. Additionally, random forest (RF) (Breiman, 2001;Liaw and Wiener, 2002), integrating random subspace sampling with bagging, enhances the ensemble's performance. This integration promotes model de-correlation, thereby elevating the ensemble's ability to generalize effectively in the presence of noise and outliers. However, these tree-based ensemble methods were not originally designed to capture interaction selection in ultrahigh dimensional data. The random intersection trees (RIT), proposed by Shah and Meinshausen (2014), is an algorithm built on the modification of RF that aids the search for stable, high-order non-multiplicative interactions. The iterative random forest (iRF), developed by Basu et al. (2018), is an extension of the RIT that can capture non-hierarchical, higher-order interactions.
This study aims to compare the predictive performance of two approaches for interaction selection: penalty-driven algorithms, which mainly capture two-way multiplicative interactions in regression, and tree-based ensemble algorithms, which excel at selecting higher-order nonmultiplicative interactions. While both penalized regression and tree-based ensemble methods are popular supervised learning techniques, there is a research gap comparing their predictive performance in high and ultra-high-dimensional data, especially in the presence of interactions. We perform comparative simulation studies and empirical data analyses to address this gap.
The rest of the paper is organized as follows. Section 2 describes the variable selection algorithms considered in our comparative analyses. Section 3 summarizes the settings and results of extensive simulation experiments comparing the predictive and variable selection performance of penalized regression methods and ensemble tree-based methods. Section 4 compares the two sets of methods on a real high-dimensional data set. We conclude the paper with a discussion of the main results in Section 5.
This section describes selected predictive modeling algorithms that perform interaction selection: LASSO, SCAD, MCP, RAMP, and iRF. Let p denote the number of predictors, and let N denote the sample size of the training data. A linear regression model with p main effects and their two-way multiplicative interactions is defined as follows:
where ε is the independently and identically distributed error term with E(ε) = 0 and E(ε 2 ) = σ 2 . In matrix notation, we can rewrite (1) as follows:
where y = (Y 1 , Y 2 , . . . , Y N ) T represents the response vector, and ε = (ε 1 , ε 2 , . . . , ε N ) T denotes the corresponding error vector. The augmented coefficient vector is denoted as β = (β T 1 , β T 2 ) T where β 1 = (β 1 , β 2 , . . . , β p ) T represents the main effects coefficients vector, and β 2 = (β 1,2 , . . . , β p-1,p ) T = (β p+1 , . . . , β q ) T corresponds to the vector of interaction coefficients with q = p + p(p -1)/2. Each X i is standardized. The augmented design matrix is given by Z = (Z 1 Z 2 ) where Z 1 and Z 2 represent specific components of the matrix.
.
In statistical modeling, interaction terms involve two fundamental principles for interaction selection: the hierarchical principle and the marginality principle. 1. The hierarchical principle states that when interaction effects are selected, we also include the main effects irrespective of their significance (McCullagh, 2002;Yuan et al., 2009;Zhao et al., 2009;Choi et al., 2010). 2. Conversely, the marginality principle states that when main effects are selected, we consider the inclusion of interaction effects among the selected main effects (Nelder, 1977;Chipman et al., 1997;McCullagh, 2002;Zhao et al., 2009;Choi et al., 2010).
• Strong heredity rule: If both predictors X k and X ℓ are selected, then we include the twoway multiplicative interaction of them. • Weak heredity rule: If either X k or X ℓ is selected, then we include the two-way multiplicative interaction of them.
The least absolute shrinkage selector operator (LASSO) is a penalized regression for variable selection proposed by Tibshirani (1996). LASSO solves the following minimization problem:
where (Zβ) i is the i-th row of Zβ. In (3), we aim to estimate β instead of β 1 , which are the same for the SCAD and MCP objective functions in (4) and (5). These three penalized regression methods do not follow the statistical principles regarding interaction selection described above. Under the LASSO penalty, the regression coefficients of unimportant predictors are shrunk towards zero, with some forced to be exactly zero via a sufficiently large λ. When λ = 0, optimum β is equal to the ordinary least squares (OLS) solution. One of the significant limitations of this augmented design LASSO is that it considers interaction terms while ignoring the main effects and vice versa. Further, LASSO can select only a maximum of N variables when p > N.
Moreover, in the presence of strongly correlated variables, LASSO may choose only one of these variables at random and ignore the others. LASSO excels in high-dimensional settings by creating sparse models and reducing coefficients to zero for effective feature selection. However, it often falls short in identifying interaction terms, particularly when main effects are weak, due to its propensity to select only one variable from groups of correlated predictors.
The Smoothly Clipped Absolute Deviation (SCAD) penalty, proposed by Fan and Li (2001), is a nonconvex regularization technique used in regression models to prevent overfitting by penalizing large coefficients. Unlike L 1 regularization, which imposes a fixed penalty regardless of the β size, SCAD has a unique property that can be advantageous in scenarios where overly penalizing large β is a concern. The key feature of SCAD that makes it useful for not overly penalizing large β s is its flattening effect for β s beyond a certain threshold. The SCAD method solves the following optimization problem:
where
where λ is the regularization parameter, and a is a predefined constant (often set to 3.7 for SCAD). For small β, the SCAD penalty maintains LASSO's penalization rate (and bias), but as the absolute value β j increases, the rate of penalization gradually decreases. The SCAD penalty is continuously differentiable on (-∞, 0) ∪ (0, ∞), but it is singular at zero with zero derivatives outside the range [-aλ, aλ]. Small coefficients are set to zero as a result of this, while a few additional coefficients tend to approach zero while large coefficients are retained.
The Minimax Concave Penalty (MCP), proposed by (Zhang, 2010), is another popular nonconvex penalty method to achieve variable selection and parameter estimation in high-dimensional data. MCP, similar to SCAD, starts by penalizing coefficients at the same rate as the LASSO, gradually reducing the rate toward zero as the absolute value of the coefficient increases. However, unlike SCAD, MCP immediately relaxes the penalization rate, while SCAD maintains a stable rate before reducing it. Although the unbiasedness and selection criteria we impose on the penalty function prohibit the use of fully convex penalties, MCP achieves sparse convexity to a great extent by constraining the maximum concavity (Zhang, 2010). The augmented design SCAD and MCP have a similar issue to the LASSO in the interaction selection in the sense that they do not consider the hierarchical and marginality principles. It should be noted that while MCP is not entirely convex, it still provides a substantial level of sparsity. Mathematically, the MCP optimization problem is formulated as follows:
where
SCAD and MCP are known for their oracle properties, meaning that, under certain conditions, they can correctly identify non-zero coefficients that achieve asymptotic normality as sample size increases. On the other hand, LASSO may not possess this property without modifications.
Adaptive LASSO (Zou, 2006) was proposed to overcome this issue in LASSO. In practical terms, the choice between LASSO, SCAD, and MCP often comes down to the specific data scenario, the size of the dataset, the computational resources available, and the degree of bias one is willing to accept in the final model. LASSO is typically favored for its simplicity and computational efficiency, while SCAD and MCP are often chosen for their ability to reduce bias and handle complex data structures more effectively (Breheny and Huang, 2015).
The RAMP algorithm, proposed by Hao et al. (2018), uses the marginality principle to preserve a model's hierarchical structure following variable selection. The RAMP algorithm can be applied to both linear regression and logistic regression. Consider the interaction model in (1) and define the linear main effects index set M = {1, 2, . . . , p} and the second-order interaction index set I = {(j, k) : 1 j < k p}. The objective function of the RAMP algorithm can be written as follows:
where the penalty is imposed on the candidate interaction effects. The set M ℓ-1 is the active main effects set at iteration ℓ-1. The set H ℓ-1 is defined as the parent (main effects) set corresponding to I ℓ-1 , the active interaction set at iteration ℓ -1. More precisely, H ℓ-1 comprises the main effects which have at least one associated interaction effect (referred to as a "child") in I ℓ-1 .
Conversely, H c ℓ-1 representing the complement of H ℓ-1 , consist of those main effects in M that are not linked by the strong heredity constraint to H ℓ-1 , i.e., H c ℓ-1 = M -H ℓ-1 . Throughout the solution path, M ℓ , I ℓ and H ℓ are iteratively updated until all the main effects and the interaction effects are selected and meet the hierarchy and marginality principles. The RAMP algorithm allows for three penalty functions (LASSO, SCAD, and MCP) and two heredity rules (weak and strong). The RAMP algorithm is given in Algorithm 1.
Algorithm 1 RAMP (Hao et al., 2018).
3: Initialize the main effect set M 0 = ∅ and interaction effect set I 0 = ∅. 4: Path-building: Repeat the following steps for ℓ = 1, 2, . . . , L. Given M ℓ-1 , I ℓ-1 , H ℓ-1 , add the possible interactions among main effects in M ℓ-1 to the current model.
) T , we minimize (6).
RF is an ensemble machine learning algorithm that constructs multiple decision trees for classification or regression, enhancing prediction accuracy and reducing overfitting by aggregating the outputs from various trees (Breiman, 2001). Each tree in the RF is built using a random subset of data and features, contributing to algorithm diversity. Predictions are made through majority voting in classification tasks or averaging in regression. Noted for its robustness and effectiveness across different data types, RF can be computationally intensive and less interpretable. The RF algorithm is given in Algorithm 2.
Algorithm 2 Breiman's random forest (Hastie et al., 2004).
, where X = (X 1 , X 2 , . . . , X p ) is p-predictor vector and Y is the response.
1: Generate B different bootstrap samples, each with N observations, from the original training dataset. 2: For each bootstrap sample numbered b (where b ranges from 1 to B), construct a decision tree fb . This is done by iteratively performing the following steps for each terminal node of the tree until the size of the node is reduced to or below a pre-set minimum, n min • Randomly choose m out of the total p feature variables.
• Identify and execute the optimal split from among the selected m features.
• Divide the current node into two child nodes.
Output: the ensemble of trees f1 , f2 , f3 , . . ., fB is used to produce the final prediction for new points as follows:
A crucial aspect of the RF algorithm involves using a subset of features randomly selected from the entire feature space considered for each tree using a bootstrapped sample to train decorrelated trees. This random subset of features mainly differentiates RF from bagging. While selecting a random subset of features is simple, determining the most effective feature for splitting is more complex. This selection uses different metrics depending on whether the RF is used for classification or regression. The Gini impurity is used to measure the likelihood of incorrect classification at a node by considering the proportion of misclassified observations, thereby guiding the algorithm towards the optimal split for classification trees. For regression trees, the best split is determined by minimizing the total sum of squared deviances between the actual and predicted values (residuals). This approach helps in making accurate and efficient splits for regression tasks. Moreover, before each splitting event, m p input variables are randomly chosen as potential candidates. The default setting for m in classification tasks is √ p, with a minimum node size of one. For regression, m is typically set at p/3, and the minimum node size is determined to be five (Hastie et al., 2004).
An alternative to the multiplicative interaction method, such as RAMP, includes a popular treebased interaction method (Breiman et al., 1984;Breiman, 1996Breiman, , 2001;;Meinshausen, 2010). In particular, RF achieves robust and accurate prediction performance while mitigating overfitting and leveraging high-order interactions. The iterative random forest (iRF) is an extension of the traditional RF to model stable, predictive high-order interactions for classification or regression tasks. The iRF algorithm allows for the iterative refinement of feature importance scores by repeatedly building random forests on reduced feature sets. This iterative process helps identify and prioritize the most informative features for the task. It can be particularly useful when dealing with high-dimensional datasets with many features, as it can improve algorithm performance and reduce overfitting by focusing on the most relevant features.
The iRF algorithm is a set of feature-weighted ensemble decision trees created to find stable high-order interaction terms. We summarize important steps of the iRF algorithm for a binary classification problem as in Basu et al. (2018). Let U i ∈ {0, 1} be a binary response, and let I i ⊆ {1, 2, . . . , p} be the feature-index subset for i = 1, 2, . . . , N, which are used for identifying features in leaf nodes of RF trees. Let D = {(X i , Y i )} N i=1 be the training data, where X = (X 1 , X 2 , . . . , X p ) are predictors, and Y ∈ {0, 1} is a binary response. The algorithm detects higher-order interactions via three steps. In the first step, adaptive regularization is applied to the RF fitting through iterative feature re-weighting, which determines the relative importance of features. In the second step, the feature-reweighted RF is utilized to identify prevalent interactions within the RF using a generalized version of the random intersection tree (RIT) method, as proposed by Shah and Meinshausen (2014). More specifically, the computationally efficient RIT algorithm searches order-s interactions that appear with higher frequency in a given class, which leads to the construction of informative interaction set S ⊆ {1, 2, . . . , p} in class C ∈ {0, 1}. The prevalence of an interaction for binary classification is defined as
where S ⊆ {1, 2, . . . , p}, I (•) denotes an indicator function, and 0 sta(S) 1. The bagged stability score is the proportion of times an interaction appears as an output of RIT in B bootstrap samples, which is defined as
In our simulation and empirical studies, we identified the interactions with sta(S) 0.7. Order-s interaction terms can be addressed and identified through Algorithm 3.
Using bootstrap sampling, the RIT algorithm searches the interaction set S satisfying the following conditions (see Algorithm S1 in Supplementary Material):
where 0 π 0 < π 1 1. In the last step, bagging evaluates the stability of recovered interactions from bootstrap perturbations.
The iterative random forests offer a dynamic approach to capturing and adapting to interactions in evolving data. Tracking variable importance, employing dynamic feature selection and engineering, using visualization tools, and considering ensemble learning techniques contribute to the algorithm's ability to identify and respond to interaction effects over multiple iterations.
Our comparative study will evaluate the prediction performance of the selected methods for regression and classification using the following assessment metrics.
Algorithm 3 iterative random forests (Basu et al., 2018).
Input: D, C ∈ {0, 1}, B,K , w (1) ← (1/p, . . . , 1/p)
Fit RF (w (k) ) on D Iterative reweighted RF 3:
Generalized RIT (through RF (w (K) ))
8: (ii) The trained iRF algorithm: RF (w (K) ) .
We used the mean squared prediction error (MSE) as the metric for assessing the predictive performance in the regression case, i.e., continuous response:
where Y i is the actual response value, Ŷi is the predicted value, and N is the size of the training/testing set.
We use the following metrics:
• Sensitivity (or Recall): measures the proportion of actual positives (true positive cases) correctly identified by the model.
• Specificity: measures the proportion of actual negatives (true negative cases) correctly identified by the model.
• Accuracy: is the proportion of true results (both true positives and true negatives) in the total cases as a measure of the overall correctness of the model.
• Balanced Accuracy: is the average of sensitivity and specificity, and often used for imbalanced response.
Balanced Accuracy (BalAcc) = Sensitivity + Specificity 2 • F1-Score: is the harmonic mean of precision and recall, providing a balance between them.
Precision * Recall Precision + Recall ,
where the precision is defined as TP/(TP+FP). Finally, we assess variable selection performance using the proportion of accurately captured interaction terms (coverage rate for interaction terms). It is worth noting that higher sensitivity, specificity, F-1 score, and balanced accuracy indicate better predictive performance, whereas lower MSE indicates improved predictive performance.
In this section, we introduce various simulation experiments and summarize their results to compare the performance of penalty-driven algorithms such as LASSO, SCAD, MCP, and RAMP with two tree-based ensemble algorithms, RF and iRF, under both the regression and classification scenarios. As mentioned in the previous section, while all of the methods are capable of capturing interaction terms in the model, LASSO, SCAD, MCP, and RF were not primarily designed for this purpose, whereas RAMP and iRF were specifically designed for capturing certain types of interactions. We separately consider regression and classification models in the following subsections. The glmnet package was used for LASSO implementation (Friedman et al., 2010), the ncvreg package was used for SCAD and MCP implementation (Breheny and Huang, 2011), the randomForest package was used for RF implementation (Liaw and Wiener, 2002), the RAMP package was used for RAMP implementation (Feng et al., 2020), and the iRF package was used for RF implementation (Basu and Kumbier, 2018). The tuning parameters for the selected algorithms were determined using the default options provided by each package except for iRF, where the tree depth was set to 20 instead of the default value of 5. We measure predictive performance using the mean squared prediction error (MSE) for the regression case and using accuracy, sensitivity, specificity, F-1 score, and balanced accuracy for the classification case.
We used four different simulation scenarios for regression-based response cases. The simulation scenarios include: 1) a perfect hierarchical structured model with a multiplicative two-way interaction term; 2) a non-hierarchical model (no main effect term); 3) a nonlinear model that does not conform to obeying the marginality principle and; 4) a nonlinear model that obeys the marginality principle (hierarchical structure). Specifically, we consider the following 4 models: 1. Linear hierarchical model:
where the parameter vector {β 1 , β 2 , β 3 , β 1,2 , β 3,4 } = {0.2, 0.3, 0.4, 0.3, 0.3}. In all 4 models, the predictors were generated from a multivariate normal distribution with a zero mean vector and a covariance matrix with standard deviation σ = 1, non-zero covariance among {X 1 . . . , X 5 }, and zero covariance among all other predictors as shown below (Jain and Xu, 2021):
.0 0.3 0.3 0.6 0.6 0 . . . 0 0.3 1.0 0.3 0.2 0.1 0 . . . 0 0.3 0.3 1.0 0.2 0.1 0 . . . 0 0.6 0.2 0.2 1.0 0.1 0 . . . 0 0.6 0.1 0.1 0.1 1.0 0 . . . 0 0 0 0 0 0 1.0 . . . 0 . . . . . . . . . . . . . . . . . . . . . . . .
.
Models 1 and 2 appeared in Jain and Xu (2021), whereas models 3 and 4 are similar to the models used by Van der Laan et al. (2007). In all 4 models, the errors, ε, were generated from the normal distribution with mean zero and standard deviation σ = 0.5. I (A) is the usual indicator function which takes 1 if event A holds and zero otherwise. Under each scenario, we varied the number of predictors, p = {25, 100, 500, 1000}. In all scenarios, the number of simulation iterations was fixed at iter = 100, and the full data sample size was n = 500. The full data was partitioned into a 60% training set and a 40% test set in each simulation iteration. Two evaluation metrics were used to compare the various algorithms under study. The first metric focuses on the efficacy of the algorithm in interaction selection and was defined as the rate at which the true interaction terms present in the underlying model are correctly identified (coverage rate of interaction terms). High coverage means that the algorithms capture most of the interaction terms affecting the outcome. The second metric, taken to be the training and test MSE, was used to assess the algorithm's predictive performance. The results of our simulations for the regression case are presented in Tables 12.
In Table 1, RF shows a better interaction selection coverage across the four models except for the nonlinear non-hierarchical case when p = 25, where iRF performs better than RF and other methods. The RAMP algorithm shows better selection coverage than iRF for the linear hierarchical model. Moreover, the table reveals that other non-hierarchical, penalty-based algorithms nearly perfectly identify interaction terms. Although ensemble methods such as iRF are potent for nonlinear models, they can be computationally intensive, particularly when the predictor count p greatly exceeds the number of observations N . Notably, these methods are not limited to two-way interactions; they can identify and select higher-order interactions.
In terms of predictive accuracy, both the weak and strong RAMP rules demonstrate steady performance across various numbers of predictors (p), with Mean Squared Error (MSE) ranging Table 3: Conditional probabilities π := P (Y = 1|X) for classification models.
Linear hierarchical 1 + exp -(0.5 + A)
Nonlinear hierarchical 1 + exp -(0.5 + D)
from 0.244 to 0.264. This consistent performance is noted across all models, with the exception of the linear non-hierarchical model. Additionally, unlike some alternative approaches, the RAMP rules do not exhibit overfitting issues, as indicated in Table 2. In linear hierarchical models, RAMP-W and RAMP-S consistently demonstrate the lowest test MSE across all p values. For linear non-hierarchical models, SCAD and MCP also tend to perform well, but SCAD and MCP show competitive performance on test sets, especially as p increases.
In nonlinear non-hierarchical and hierarchical models, iRF exhibits the best test MSE at higher p values (p = {500, 1000}), while RF shows better test MSE for the small p cases (p = {25, 100}). There is an overall trend where simpler methods like LASSO tend to perform less effectively as the complexity of the model (nonlinearity) and dimensionality (p) increase. The advanced ensemble methods (RF and iRF) demonstrate resilience to the increase in p for the nonlinear (non-hierarchical and hierarchical) model, maintaining low MSE values and suggesting their suitability for complex, high-dimensional data modeling. It is also notable that the hierarchical model types do not consistently outperform non-hierarchical ones, as might be expected, which could be due to the specific nature of the data or the interactions being modeled.
In the classification context, we simulated data considering balanced and imbalanced responses for four models. Set
The binary response variable Y for the classification simulations was generated using the conditional probability defined in Table 3. The predictor variables and the error term for the classification simulations are generated in the same manner as described in the regression case in Section 3.1. Similar to the regression case, the full data in each simulation iteration was split into a 60% training set and a 40% test set. We again use the rate of correctly identifying interaction terms as our interaction selection performance metric and report the results in Table 4. For predictive performance, we calculated five metrics for both the training and testing sets. The five performance metrics comprised accuracy, balanced accuracy, sensitivity, specificity, and F1-score. The predictive performance results for the classification scenarios are summarized in Figures 1234. We first analyze the interaction selection results. Overall, RF demonstrates superior performance over the RAMP and other penalty-based interaction selection methods. iRF also showed good coverage results when it didn't suffer from convergence issues, i.e. when p = {25, 100}. For large values of p (500 and 1000), iRF failed as indicated by the NA's in Table 4. This is likely due to the computational burden associated with such extensive interaction selection for large p cases. Finally, the penalty-based algorithms, LASSO, SCAD, MCP, and RAMP, showed their best performance under the linear hierarchical model, especially for small p = 25, but they suffered under other models, especially the nonlinear non-hierarchical model.
We now turn our attention to the predictive performance results in Figures 1234. Analyzing the training set curves reveals that the RF and iRF algorithms attain robust performance in the metrics, regardless of the number of predictors (p). Among the other five algorithms, LASSO shows slight superiority across the various performance metrics. The four algorithms, SCAD, MCP, RAMP-W, and RAMP-S, have similar performance in terms of accuracy, balanced accuracy, F1 score, and sensitivity. RAMP-W and RAMP-S have lower training specificity than SCAD and MCP under three of the four scenarios considered, especially as p increases.
For the test set results, all algorithms underperform across the various metrics, but they all maintain reasonable performance in terms of sensitivity, F1 score, and accuracy. Despite the big drop in its specificity, iRF maintains the highest specificity rates across all models and dimensions, except for p = 1000, where RAMP-S and RAMP-W show competitively higher specificity rates. Another noteworthy observation is how the RF's testing specificity drops for higher dimensions, p = 500, 1000. The dramatic change in the performance of the two tree-based algorithms (RF and iRF) when transitioning from training to testing data hints at the potential challenge of overfitting for these algorithms. We also examined the case of imbalanced response for the same four classification models. The results were quite similar to the balanced response case. Thus, these results were deferred to the supplementary material section.
Breast cancer is the most frequent cancer among women, accounting for about one-fifth of all malignancies diagnosed in women worldwide. It is the second most common cause of death from cancer among women. In this section, we compare the performance of the predictive modeling algorithms described in earlier sections using the bcTCGA breast cancer gene expression data from The Cancer Genome Atlas Program (TCGA) of the National Cancer Institute's Center for Cancer Genomics: bcTCGA. The bcTCGA data contain information on 17, 323 gene expressions from 536 women with breast cancer. The BRCA1 gene, a key discovery in breast cancer research, along with the BRCA2 gene (identified a year later), are the core causes for 70% of breast cancer cases. Due to BRCA1 interactions with other genes, identifying genes that interact with it is crucial for further research (Deng and Brodie, 2000;Antoniou et al., 2003;Kuchenbaecker et al., 2017). The selected techniques are particularly valuable in this context for reducing the vast number of potential gene interactions to a more manageable set that is significant for cancer development. These interactions can then be further studied in a laboratory setting and targeted for cancer treatment development. In our regression analysis, the BRCA1 gene expression level serves as the response variable, with the other 17322 genes serving as predictors. This scenario is a typical example of the p n scenario. A two-fold cross-validation approach was used where we repeatedly (100 times) split the dataset into training and testing sets at a ratio of 6 : 4. In each iteration, we computed the training and testing MSE for each of the seven predictive algorithms under study. Table 5 displays the average training and testing MSE for the seven algorithms. SCAD outperforms other algorithms with the minimal test MSE observed and higher standard deviation. The RAMP algorithm demonstrates a marginally better test MSE of 0.200 compared to iRF's 0.220. Both RF and iRF exhibit signs of overfitting.
For the classification case, we use the Wisconsin breast cancer dataset available at the UCI Machine Learning Repository (Wolberg et al., 1995). The dataset consists of 569 observations on 30 predictors and one binary response variable representing tumor diagnosis with 62.742% benign (non-cancerous) cases and 37.258% malignant (cancerous) cases. The features, which describe characteristics of the cell nuclei, are computed from a digitized image of a fine needle aspirate of a breast mass. The features included radius, texture, perimeter, area, smoothness, etc. This is a typical classification problem that can be handled by any of the predictive modeling algorithms studied here. Similar to the regression case, we used a two-fold cross-validation approach by repeatedly (100 times) splitting the dataset into training and testing sets at a ratio of 6 : 4. The predictive performance results averaged over 100 iterations are summarized in Table 6. As highlighted in the table, RAMP achieved the highest scores in test accuracy, balanced accuracy, and sensitivity. MCP excelled in specificity, and SCAD led in F1 score performance. Notably, SCAD's test accuracy, balanced accuracy, and F1 score results were closely comparable to those of RAMP. Finally, it is readily seen that the tree-based algorithms (RF and iRF) suffered overfitting issues, with the test set results being significantly lower than the training set ones.
This study compared two predictive modeling approaches that enable interaction selection: the penalty-based approach (LASSO, SCAD, MCP, and RAMP) and the tree-based approach (RF and iRF). Our findings from extensive simulations and real data applications revealed that RAMP had superior predictive performance under linear hierarchical regression scenarios, exhibiting lower MSE across varying dimensions (different values of p). iRF showed signs of overfitting in prediction performance in most regression and classification scenarios, but it was found to be best-suited for nonlinear hierarchical/non-hierarchical large p regression scenarios (p = {500, 1000}). Similarly, RF suffered from overfitting but was best suited for the nonlinear hierarchical and non-hierarchical small p regression scenarios (p = {25, 100}). SCAD and MCP performed superiorly under linear non-hierarchical scenarios.
While all of the six techniques under comparison are capable of capturing interactions, RAMP and iRF are specifically designed to capture interactions and, hence, were expected to have superior interaction selection performance. Overall, our interaction selection results revealed mixed results across most scenarios under both regression and classification tasks. For regression tasks, non-hierarchical penalty-based algorithms such as LASSO, SCAD, and MCP also demonstrated superior interaction selection performance compared to RAMP (in most scenarios) and iRF (except for non-linear non-hierarchical models). iRF outperformed RAMP in selecting interactions for regression tasks except under linear hierarchical regression models where RAMP was superior to iRF. In classification tasks, RAMP had superior interaction selection performance to the other penalty-based algorithms, while iRF did not deliver stable results, especially for the large p scenarios.
Taking the above results collectively, we were not be able to declare one algorithm as a clear winner across all or even the majority of scenarios. There were at least one or more scenarios that favored each of the algorithms included in this comparative study. However, from the general pattern, we can recommend 1) the use of RF and iRF for regression/classification tasks