Edited by: Edward Wilson Ansah, University of Cape Coast, Ghana
Reviewed by: Qi Song, Carnegie Mellon University, United States; Prashant Kaushik, Yokohama Ueki, Japan
This is an openaccess article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
The use of highdimensional data has expanded in many fields, including in clinical research, thus making variable selection methods increasingly important compared to traditional statistical approaches. The work aims to compare the performance of three supervised Bayesian variable selection methods to detect the most important predictors among a highdimensional set of variables and to provide useful and practical guidelines of their use. We assessed the variable selection ability of: (1) Bayesian Kernel Machine Regression (BKMR), (2) Bayesian Semiparametric Regression (BSR), and (3) Bayesian Least Absolute Shrinkage and Selection Operator (BLASSO) regression on simulated data of different dimensions and under three scenarios with disparate predictorresponse relationships and correlations among predictors. This is the first study describing when one model should be preferred over the others and when methods achieve comparable results. BKMR outperformed all other models with small synthetic datasets. BSR was strongly dependent on the choice of its own intrinsic parameter, but its performance was comparable to BKMR with large datasets. BLASSO should be preferred only when it is reasonable to hypothesise the absence of synergies between predictors and the presence of monotonous predictoroutcome relationships. Finally, we applied the models to a real case study and assessed the relationships among anthropometric, biochemical, metabolic, cardiovascular, and inflammatory variables with weight loss in 755 hospitalised obese women from the Follow Up OBese patients at AUXOlogico institute (FUOBAUXO) cohort.
We compared the variable selection ability of three Bayesian variable selection models: Bayesian Kernel Machine Regression (BKMR), Bayesian Semiparametric Regression (BSR) and Bayesian Least Absolute Shrinkage and Selection Operator (BLASSO) using both simulated and real data.
For large sample size data, BKMR and BSR may be employed indiscriminately, as they were able to capture complex relationships between predictors and the outcome even in presence of highly correlated variables.
For small sample size data: BKMR correctly selected the important variables, while BSR was strongly dependent on the choice of its tuning parameter and additional considerations should be taken into account to ensure an accurate variable selection.
BLASSO should be preferred only when it is reasonable to hypothesise the absence of interactions between predictors and in the presence of monotonous relationships.
In general, the performance of these models strongly depends on the sample size of the dataset, the correlation structure, and the predictoroutcome relationships.
The advent of exposomics and the abundance of novel toxicological information have made the identification of the associations between both exposures and their molecular responses and a health outcome more challenging. Therefore, many variable selection approaches have increased in importance and popularity (
Among these Bayesian methods, those employing spikeandslab priors (
In this work, we compared three supervised Bayesian models: the Bayesian Kernel Machine Regression (BKMR) (
In this section we summarise the methodologies applied throughout the study. We assume a continuous outcome variable Y and several main predictors Z and confounders/covariates X for all models. The BKMR model has the following form:
As a result of variable selection, the BSR provides the PIPs for the main effects and the interactions. Lastly, the LASSO is a technique for improving the ordinary least squares estimates for a linear exposureresponse function h by imposing the L1 penalty on the regression coefficients, shrinking or setting them to zero:
The R programming language was used for all analyses, along with opensource software packages for BKMR (
We compared the BKMR, BSR, and BLASSO models in terms of their ability with feature selection and performance prediction in three simulated scenarios with: (1) linear, (2) quadratic, and (3) logistic predictorresponse relationships. We considered each scenario under two circumstances: with and without interactions among predictors. Changes in correlations structures and sample size (small,
The correlation structures were set as moderate (
High and lowcorrelation structure for the simulated datasets.
We assessed the robustness of each model by changing the underlying priors in two situations: for BKMR and BSR, we assumed that a mixture component z_{m} was included in the model with a prior probability π following a beta(a_{π} = 2, b_{π} = 2) or beta(a_{π} = 2, b_{π} = 6) distribution such that,
For each of these configurations, we generated K = 100 simulated datasets with M = 10 predictors
Predictorresponse relationships for the three scenarios.
Scenario 1: linear predictorresponse associations  
a)No interactions 

b)With interactions 

Scenario 2: quadratic predictorresponse associations  
a)No interactions 

b)With interactions 

Scenario 3: logistic predictorresponse associations  
a)No interactions 

b)With interactions 


Simulation summary.
Training and validation sets were selected with an 80/20 ratio for each dataset. We computed the average PIP of each predictor from the training datasets and then reported mean PIPs and standard errors (SEs) for all variables of interest. We fit BSR for a few d values, which are considered optimal in Antonelli et al. (
We assessed the performance prediction on the validation sets using the mean squared error (MSE), reporting the mean and the standard deviation (SD) for each configuration. To select the important predictors in each model, we used the median probability model (MPM) (
The results in terms of feature selection under each configuration and considering only the first set of priors are shown in
Mean PIPs for each configuration in Scenario 1 (linear predictorresponse associations). Results for prior setting 1 shown [π following a beta (a_{π} = 2, b_{π} = 2) for BKMR/BSR, λ^{2} following a gamma (α = 1, β = 0.5) for BLASSO].
Mean PIPs for each configuration in Scenario 2 (quadratic predictorresponse associations). Results for prior setting 1 shown (π following a beta (a_{π} = 2, b_{π} = 2) for BKMR/BSR, λ^{2} following a gamma (α = 1, β = 0.5) for BLASSO].
Mean PIPs for each configuration in Scenario 3 (logistic predictorresponse associations). Results for prior setting 1 shown (π following a beta (a_{π} = 2, b_{π} = 2) for BKMR/BSR, λ^{2} following a gamma (α = 1, β = 0.5) for BLASSO].
Under the first scenario, all the models were able to select the predictors
In the second scenario, all the models were able to identify the variables
Finally, in the third scenario, the BKMR was able to select the predictors
To test the predictive ability of the models on the validation set, we used MSE. The results in terms of the mean MSE for each configuration and the first prior setting are shown in
Mean MSE (and SD) for each model for the different simulations’ configurations.
Scenario  Interaction  Sample size  Correlation  BKMR  BSR (d = 1)  BSR (d = 2)  BSR (d = 3)  BSR (d = 4)  BLASSO  

1  No  100  Low  7.6 (2.2)  7.5 (2.2)  7.8 (2.1)  8.2 (2.3)  8.4 (2.3)  7.5 (2.1)  High  10.0 (3.3)  9.7 (3.1)  9.9 (3.3)  10.2 (3.4)  10.2 (3.4)  9.7 (3.0)  1,000  Low  7.8 (0.6)  7.6 (1.1)  7.6 (1.1)  7.6 (1.1)  7.7 (1.1)  7.7 (0.8)  High  9.5 (0.8)  9.4 (1.2)  9.4 (1.2)  9.4 (1.2)  9.5 (1.2)  9.4 (1.2)  Yes  100  Low  16.4 (7.1)  19.0 (7.9)  15.7 (6.7)  16.2 (6.5)  16.7 (6.8)  24.8 (11.4)  High  20.5 (9.9)  20.6 (9.1)  19.9 (9.2)  20.3 (9.5)  20.8 (9.9)  30.1 (17.8)  1,000  Low  12.0 (4.4)  16.9 (2.1)  13.2 (1.8)  13.3 (1.8)  13.3 (1.8)  25.3 (3.7)  High  17.4 (1.1)  17.7 (2.3)  17.0 (2.3)  17.1 (2.3)  17.1 (2.3)  33.8 (5.8)  
2  No  100  Low  12.8 (4.8)  12.6 (4.2)  13.2 (4.3)  14.5 (5.1)  14.3 (5.2)  14.9 (5.9)  High  17.4 (7.2)  17.3 (7.3)  17.7 (7.4)  18.2 (7.5)  18.5 (7.6)  22.0 (9.6)  1,000  Low  11.2 (1.0)  11.5 (1.6)  11.7 (1.6)  11.9 (1.7)  12.0 (1.7)  15.6 (2.2)  High  15.6 (1.0)  15.5 (2.1)  15.7 (2.1)  15.8 (2.1)  15.9 (1.8)  21.4 (3.0)  Yes  100  Low  36.7 (23.5)  38.9 (23.0)  37.6 (35.0)  37.9 (28.2)  39.1 (24.9)  54.9 (35.5)  High  53.9 (22.1)  52.3 (21.2)  54.7 (24.3)  54.9 (23.8)  57.5 (26.0)  84.7 (44.3)  1,000  Low  28.1 (6.6)  31.6 (3.9)  28.3 (3.6)  28.7 (3.7)  29.3 (3.8)  54.5 (8.6)  High  44.1 (4.2)  43.8 (5.8)  43.2 (5.8)  43.5 (5.8)  43.6 (5.9)  82.9 (13.3)  
3  No  100  Low  0.133 (0.044)  0.142 (0.047)  0.143 (0.047)  0.136 (0.046)  0.138 (0.046)  0.141 (0.046)  High  0.147 (0.051)  0.149 (0.052)  0.148 (0.049)  0.145 (0.053)  0.144 (0.051)  0.149 (0.050)  1,000  Low  0.133 (0.013)  0.138 (0.018)  0.134 (0.018)  0.132 (0.017)  0.132 (0.017)  0.143 (0.018)  High  0.143 (0.012)  0.146 (0.019)  0.143 (0.019)  0.141 (0.019)  0.141 (0.019)  0.156 (0.020)  Yes  100  Low  0.133 (0.050)  0.141 (0.054)  0.140 (0.052)  0.137 (0.050)  0.139 (0.048)  0.148 (0.055)  High  0.117 (0.040)  0.120 (0.044)  0.123 (0.048)  0.118 (0.042)  0.121 (0.045)  0.135 (0.044)  1,000  Low  0.118 (0.011)  0.127 (0.016)  0.122 (0.017)  0.120 (0.016)  0.120 (0.017)  0.154 (0.019)  High  0.106 (0.014)  0.118 (0.019)  0.111 (0.020)  0.109 (0.015)  0.108 (0.015)  0.150 (0.020) 
Results for prior setting 1 shown [π following a beta (a_{π} = 2, b_{π} = 2) for BKMR/nBSR, λ^{2} following a gamma (α = 1, β = 0.5) for BLASSO].
Here we present a case study to compare the models using data from the Istituto Auxologico Italiano. The FUOBAUXO study, approved by the Ethical Committee of the Istituto Auxologico Italiano (protocol research project code: 18A301) (
In this study, we focused on the association between several predictors—anthropometric, clinical, and biochemical variables—and weight loss percentage at discharge on obesityaffected women. The predictors included in the models were grouped according to a medical experts’ suggestion in: anthropometric, biochemical, cardiovascular, metabolic, and inflammatory (
Demographic characteristics of the population: overall population, training set, and validation set.
Overall ( 
Training ( 
Validation ( 


Age; mean (SD)  58.9 (12.8)  58.9 (13.0)  59.1 (12.4)  0.806^{a} 
Height; mean (SD)  156.7 (6.8)  156.7 (6.7)  156.6 (7.0)  0.799^{a} 
BMI; mean (SD)  43.8 (6.7)  43.9 (6.9)  43.2 (6.1)  0.233^{a} 
Type I (30 ≤ bmi < 35); 
40 (5.3)  32 (5.3)  8 (5.3)  0.987^{b} 
Type II (35 ≤ bmi < 40); 
202 (26.8)  161 (26.7)  41 (27.2)  
Type III (bmi ≥ 40); 
513 (67.9)  411 (68.0)  102 (67.5)  
Weight at admission; mean (SD)  107.7 (18.9)  108.1 (19.3)  106.0 (16.7)  0.182^{a} 
Weight at discharge; mean (SD)  102.9 (17.7)  103.2 (18.1)  101.4 (15.8)  0.217^{a} 
Weight loss; mean (SD)  −4.8 (3.1)  −4.8 (3.2)  −4.6 (2.3)  0.246^{a} 
Weight loss (%); mean (SD)  −4.4 (2.4)  −4.4 (2.5)  −4.3 (2.0)  0.573^{a} 
The
Variables selected, with posterior inclusion probabilities (PIP) ≥ 50%, by the Bayesian kernel machine regression (BKMR), Bayesian semiparametric regression (BSR), and Bayesian Least Absolute Shrinkage and Selection Operator regression (BLASSO) models in explaining weight loss percentage during hospitalization in the FUOBAUXO cohort.
Type  Variable  BKMR  BSR (d = 3)  BLASSO 

A  Lean body mass (LBM) (kg)  99.9  81.7  82.6 
A  Waist circumference (cm)  59.8  
A  Hip circumference (cm)  85.4  
A  Body mass index (BMI) (kg/cm^{2})  98.1  
A  Heart rate (bpm)  
C  Systolic blood pressure (mmHg)  67.4  
M  Thyroidstimulating hormone (TSH) (IU/L)  99.3  100.0  
M  Glycated hemoglobin (mmol/mol)  
M  Creatinine (mg/dl)  55.8  
M  Total cholesterol (mg/dl)  54.6  
M  Calcifediol (IU/L)  66.8  
M  Fasting plasma glucose (FPG) (mg/dl)  71.1  
M  Uric acid (mg/dl)  77.2  
I  Erythrocyte sedimentation rate (ESR) (mm)  65.9  100.0  
I  Creactive protein (CRP) (mg/dl)  78.3 
Type: Anthropometric (A), Biochemical (B), Cardiovascular (C), Metabolic (M), Inflammatory (I).
In the training data among 39 predictors, 13 were selected by the models. All three models agreed in selecting the lean body mass (mean PIP: BKMR = 99.9%, BSR = 81.7%, BLASSO = 82.7%). However, we identified discrepancies among the variable selections performed by the BKMR, BSR, and BLASSO models. The BLASSO selected all anthropometric predictors except for uric acid. Both the BKMR and BSR selected thyroidstimulating hormone (TSH) and erythrocyte sedimentation rate, showing a nonadditive effect between the two (BKMR and BSR results are shown in
In this study, we compared three supervised Bayesian methods (BKMR, BSR, and BLASSO) that were designed to perform variable selection and evaluate the association of multiple predictors on an outcome. We highlighted advantages and limitations of each approach under different data configurations. Both the BKMR and BSR are novel semiparametric approaches that are able to handle nonlinear and nonadditive relationships, while the BLASSO regression accommodates only additive and linear predictoroutcome associations. In the simulations, we observed that the performance of both the BKMR and BSR strongly depends on the sample size and correlation structure.
When the data sample size was large, both the BKMR and BSR correctly selected the predictors associated with the outcome, except when the true relationship was extremely weak. With a small sample size, the predictors more strongly associated with the outcome were always selected by the BKMR. With a small sample size and when the predictoroutcome relationship was weak, only the BKMR correctly selected the predictors, thus making the BKMR the most robust variable selection approach across scenarios. With a small sample size dataset and independently of the strength of predictoroutcome associations, BSR showed high heterogeneity in the variable selection abilities because of strong influence the tuning d parameter. The BLASSO captured well strong linear and additive relationships between predictors and outcome. BLASSO variable selection ability was hampered by strong nonlinear associations or interactions. When there were high correlations between predictors, PIP values were more unstable for all models. In terms of prediction ability, no differences were observed between the BKMR and BSR (for each value of d considered), while the BLASSO seemed to be characterised by higher average MSE values in most configurations. Despite the presence of highly correlated variables that are not associated with the outcome of interest, the models considered do not detect spurious relationships. This confirms their ability to handle situations of high correlation between predictors.
Our case study had a sample size comparable to the largesamplesize synthetic scenarios. When we applied all three models (BKMR, BSR, and BLASSO) to the case study, we observed that, out of all variables, the lean body mass was selected by all models. The selection of at least one anthropometric parameter, such as lean body mass, emphasised the importance of body composition in the process of weight loss. Several weight loss programs have resulted in shortterm success, but many patients fail in longterm weight loss maintenance because of the loss of lean body mass that frequently and unintentionally occurs (
Both the BKMR and BSR identified the outcome “weight loss” as linked to TSH and the erythrocyte sedimentation rate (ESR), suggesting the existence of a nonlinear relationship between those variables and the outcome. Thyroid hormones have an essential and wellknown role in body weight regulation, mainly through energy expenditure modulation. Indeed, hyperthyroidism or hypothyroidism frequently lead to significant changes in body weight (
The BKMR identified five other variables associated with the outcome, but the BSR did not select them. Based on what we have learnt with the simulated scenarios, those variables may be weakly associated with weight loss, and further studies should investigate the effective role of those variables on weight loss.
In the case study, the BSR identified fasting plasma glucose as associated with weight loss. This variable was correlated with other metabolic variables selected by the BKMR, thus suggesting a connection between metabolic parameters and weight loss. This result was in part expected, given that the beneficial effects of decreased energy intake and weight loss in blood glucose control begin to occur rapidly during rehabilitation programs. The use in the case study of BKMR and BSR, models capable of handling situations with high correlation and complex relationships with the outcome, enabled the identification of important predictors belonging to different clinical groups by selecting only one or a few variables from each group. In contrast, BLASSO selected highly correlated variables largely belonging to a single clinical group, the anthropometric group, with the greatest possibility of identifying spurious relationships between predictors and outcomes.
A few limitations of our study should be noted. Neither the synthetic nor real case scenarios considered the highdimensional framework, which can include more than 40 predictors and more than 1,000 participants. Moreover, in the simulations, we did not consider departures from the normality assumption in the linear models. Although these approaches adapt very well to highdimensional biological or clinical contexts with different underlining distribution functions, further studies should assess the performance of those approaches under the more extended framework. Furthermore, considering the role played by the prior choice in the Bayesian variable selection, we selected two prior settings, which guaranteed a low and medium inclusion probability,
This is a first attempt to compare the performance of these innovative methods beyond typically environmental contexts and to provide practical usage information. In summary, the performance of all these methods depended on the sample size, the correlation structure, and the relationship between the predictors and the outcome.
The BKMR showed excellent performance in all scenarios, both with high and low correlation and sample size, identifying predictors strongly associated with the outcome and without false positives. It was also able to identify true weak relationships.
The performance of the BSR was also excellent in all different scenarios, but showed a very strong dependence on the choice of its own parameter, especially with a small sample size or high correlations between predictors, showing a higher variability in selection abilities.
The BLASSO performed poorly under nonlinear and synergistic scenarios but had good performance under the linear and additive scenarios. Our results suggest that both the BKMR and BSR may be employed in all scenarios, with attention in selecting an appropriate tuning parameter for BSR when sample size is small. The BLASSO should be preferred when it is possible to hypothesise the absence of interactions between predictors and in the presence of monotonous relationships and low correlations between predictors.
The simulation data supporting the conclusions of this article will be made available by the authors, without undue reservation. FUOBAUXO data are available upon request from the Auxologico team. Please contact
The studies involving human participants were reviewed and approved by Ethical Committee of the Istituto Auxologico Italiano (protocol research project code: 18A301). The patients/participants provided their written informed consent to participate in this study.
NP participated in the conceptualization and design of the work, managed the acquisition of data, gave substantial contributions to the analysis or interpretation of data, wrote the first draft of the paper, gave final approval of the version published, and ensured that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved. PQ and EC participated in the conceptualization and design of the work, gave substantial contribution to the interpretation of data, critically revised the work and gave final approval of the version published. RC and MS gave substantial contribution to the interpretation of data, revised the work and gave final approval of the version published. AZ participated in the conceptualization and design of the work, gave essential contribution to the interpretation of data, revised the draft giving important intellectual contribution and gave final approval of the version published and ensured that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved. All authors contributed to the article and approved the submitted version.
During the preparation of this manuscript, EC was supported by the National Institute of Environmental Health Science (NIEHS): R01ES032242 and P30ES023515.
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
We are grateful to the patients and families who participated in the study. Special thanks also go to the Auxologico institute staff, especially Raffaella Sabatino and Raffaella Radin, for their kind collaboration and support with data collection.
The Supplementary material for this article can be found online at: