A Novel Method for Identifying a Parsimonious and Accurate Predictive Model for Multiple Clinical Outcomes

Background : Most methods for developing clinical prognostic models focus on identifying parsimonious and accurate models to predict a single outcome; however, patients and providers often want to predict multiple outcomes simultaneously. For example, older adults are often interested in predicting nursing home admission as well as mortality. We propose and evaluate a novel predictor selection method for multiple outcomes. Methods : Our proposed method selected the best subset of common predictors based on the minimum average normalized Bayesian Information Criterion (BIC) across outcomes: the Best Average BIC (baBIC) model. We compared the predictive accuracy (Harrell’s C-statistic) and parsimony (number of predictors) of the baBIC model with a subset of common predictors obtained from the union of optimal models for each outcome (Union model). We used example data from the Health and Retirement Study (HRS) to demonstrate our method and conducted a simulation study to investigate performance considering correlated and uncorrelated outcomes. Results : In the example data, the average Harrell’s C-statistics across outcomes of the baBIC and Union models were comparable (0.657 vs. 0.662 respectively). Despite the similar discrimination, the baBIC model was more parsimonious than the Union model (15 vs. 23 predictors respectively). Likewise, in two simulation scenarios with correlated and uncorrelated outcomes, the mean C-statistic across outcomes of the baBIC and Union models were very similar, and the baBIC model had on average fewer predictors. In the simulations, the baBIC method performed well by identifying the correct predictors most of the time and excluding the incorrect predictors in the majority of the simulations. Conclusions : Our method identified a common subset of variables to predict multiple clinical outcomes with superior parsimony and comparable accuracy to current methods.

Parsimonious models offer the potential to save the time it takes to gather unnecessary predictors, and expense, either in visit time or in money.
Most current model development methods focus on accurate and parsimonious prediction of single outcomes. Popular methodologies that are easy to use and interpret include stepwise methods like backward elimination or criterion-based selection like the Akaike Information Criterion (AIC) (1) or the Bayesian Information Criterion (BIC) (2). However, obtaining the most parsimonious and accurate model is more complex for the simultaneous prediction of multiple outcomes, a common scenario in clinical settings.
Several studies have demonstrated that older adults care not only about mortality, but also about their quality of life, specifically their ability to function independently (3,4). In the realm of anticoagulation for atrial fibrillation, for example, clinicians may want to simultaneously predict risk of stroke and risk of a major gastrointestinal bleed (5,6). In primary care, clinicians may want to balance risk of microvascular complications from diabetes against the risks of hypoglycemia and falls (7,8).
Yet, there is limited research on how best to develop clinical prognostic models that predict multiple outcomes simultaneously with accuracy and parsimony.
Much of the research on variable selection for multiple outcomes has been done in the highdimensional multivariate regression setting, where the number of predictors and outcomes outweighs the number of observations. Under this setting, the implementation of shrinkage or regularization methods is common (9)(10)(11). Other authors have addressed variable selection for multivariate modelling using a Bayesian framework (12)(13)(14). However, in clinical settings, where the sample size is frequently large relative to the number of predictors and outcomes, a simpler and easy-to-implement procedure that does not require complex software solutions could be of great utility.
An obvious approach (which we label Individual Outcome method) to address the multiple outcomes problem is to simply select a different subset of variables to predict each of the outcomes using selection methods for single outcomes. Although straightforward, this method could be timeconsuming, expensive (due to the cost of acquiring multiple predictors), and potentially lead to overfitting and high variability (9,15). 4 A slight modification to this approach is the Union method. In this method, we take the separate models from the Individual Outcome method, and then force the union of the predictors from each model into the predictor set for each outcome. The online compendium of prognostic indices "ePrognosis" -freely available at ePrognosis.org (16, 17)-receives over 3,000 users per week, and the most-used index is a Union model: the Combined Lee Schonberg Index, created from union of predictors from the Lee Index (18,19) and the Schonberg index (20,21). Like the Individual Outcome method, the Union method has the advantage of being a simple approach, and, additionally, it allows patients and clinicians to focus on a common subset of variables that can accurately predict their outcomes of interest simultaneously. Nevertheless, the Union model could lack parsimony as it includes all variables that predict all outcomes well, including those that are only important for some of the outcomes.
In this paper, we propose and evaluate a novel method for predictor selection in prognostic models of multiple clinical outcomes using the minimum average normalized BIC across outcomes, which we call the Best Average BIC (baBIC). To develop the proposed method, we use the Health and Retirement Study (HRS) data and a common set of health-related and demographic variables to predict time to: (1) Activities of Daily Living (ADL) Dependence, (2) Instrumental Activities of Daily Living (IADL) Difficulty, (3) Mobility Dependence, and (4) Death. We compare the parsimony and accuracy of this model with the models obtained using the Individual Outcome and the Union methods.

Case study: Health and Retirement Study data
We created a nationally representative cohort of 5,531 community-dwelling seniors enrolled in the HRS, who were 70 years old or older at the time of their baseline interview in 2000. The HRS is an ongoing longitudinal survey of a representative sample of all persons in the United States over age 50 that examines changes in health and wealth (22). It is sponsored by the National Institute on Aging (grant number NIA U01AG009740) and conducted by the University of Michigan. We used the public HRS data: Cross-Wave Tracker file (23) and RAND (24,25) HRS data file.
The pool of predictors included 39 health-related and demographic variables measured at baseline.

5
All the predictors were categorical variables. We used 4 clinical outcomes encompassing 15 years of follow-up: (1) time to first ADL dependence (including five ADLs: bathing, dressing, toileting, transferring, and eating), (2) time to first IADL difficulty (including two IADLs: managing money and medication), (3) time to first mobility dependence, and (4) time to death.

The Best Average BIC (baBIC) method
Our proposed method, the Best Average BIC (baBIC) method, selects the best subset of common predictors for M outcomes according to the baBIC. We compared our method with: (1) a method that selects individual subsets of predictors for each outcome (Individual Outcome method), and (2) an enhanced method that creates a best subset of common predictors based on the union of individual subsets obtained in the Individual Outcome approach (Union method).
Information criteria like the BIC are useful for selecting the best subset of predictors because they work well for both a fixed number of predictors and across predictor sets of varying sizes. In contrast, statistics like concordance (or C) statistic are not that useful for the selection across sets of different number of predictors since, in general, models with more predictors will tend to have higher Cstatistic than those with fewer predictors. Another advantage of the BIC is that it will tend to select more parsimonious models since it penalizes larger models more heavily compared, for example, with the AIC, or traditional stepwise regression methods based on the cut-off significance level of 0.05.
The BICs were obtained from survival models. For time to death, we fitted Cox proportional hazards regression models (26). For times to first ADL dependence, IADL difficulty, and mobility dependence, we fitted Fine and Gray competing-risk regression models to appropriately account for the risk of death (27).
In the baBIC method, we averaged the normalized BIC (nBIC) across outcomes. Normalization was important to ensure that a change in BIC from a complex to a simpler model meant roughly the same across multiple outcomes; that is, the BICs were in a comparable scale. The nBIC was computed by dividing the absolute difference between the BIC of a particular model for a specific outcome and the BIC for the "best" individual model for that outcome by the difference between the BIC in the full model (i.e. with all candidate predictors) and the BIC in the best individual model: (see Formula 1 in 6 the Supplementary Files) The nBIC thus ranges between 0 (for the best individual model) and 1 (for the model that contains all candidate predictors), with smaller being better. The nBIC can be larger than 1 for models that have a worse BIC value than the full model, but these models are typically not of interest in our setting. This normalization allowed us to average the nBIC across different outcomes and, at the same time, made this metric more interpretable.
Explicitly, we defined the baBIC for a model with k parameters as: (see Formula 2 in the Supplementary Files) The baBIC criterion has the flexibility that it could be incorporated into selection methods already available for single outcomes. Therefore, this method is not intrinsically linked to any particular method of variable selection and can be used to compare arbitrary sets of candidate models. In order to compute the nBIC for specific outcome and then the baBIC across outcomes, we need to obtain the BIC of the full model and the BIC of the best individual model. The BICs of the full and best individual models can be found using stepwise regression methods like backward elimination or more current selection methods like the Least Absolute Shrinkage and Selection Operator (LASSO) (28).
Comparison of various methods for variable selection including best subset, stepwise, and LASSO remains an area of active investigation in the statistical literature , with no one method dominating the others across a variety of settings (29).
In stepwise regression, the BIC is output at each step of the selection process so it is straightforward to find the BIC for the full model as well as the best individual model BIC value (further details in section "Application of the baBIC method to the case study HRS data" below). Similarly, selection based on minimum BIC can be directly incorporated into the LASSO setting (30,31). That is, after doing LASSO selection, we can compute BIC for each possible λ (as shown above), select the one that gives the minimum BIC as the optimal λ, and extract the corresponding BIC as the BIC of the best individual model. The BIC of the full model corresponds to the model with λ=0.

Application of the baBIC method to the case study HRS data
Using the HRS data, we compared the parsimony and predictive accuracy of the best Individual and Union models obtained using BIC backward elimination vs. LASSO selection based on optimal λ at the minimum BIC. We found that backward elimination based on minimum BIC produced more parsimonious Individual and Union models than LASSO selection, while maintaining very similar predictive accuracy (see Additional File 1). Consequently, we chose the BIC backward elimination to further illustrate the baBIC method in this setting based on the following: (1) In our example data, BIC backward elimination produced more parsimonious models with similar predictive accuracy compared to LASSO selection. (2) Stepwise selection methods are easier to implement and explain, so in general they are more widely used by clinical researchers, and they are also easily accessible in most modern statistical packages (in contrast, for example LASSO selection for survival models has not been implemented yet in SAS and Stata statistical software). (3) In the clinical settings where the sample size is usually larger than the number of predictors, issues reported for stepwise methods like instability of the selection, biased estimation of the coefficients, and multicollinearity are less important (32,33). (4) We are using backward elimination based on minimum BIC instead of using the traditional cut-off significance level of 0.05. For our sample size of 5,531 respondents, the models obtained based on minimum BIC will be more parsimonious, since the selection of a predictor in the model is approximately equivalent to a significance level of 0.01 (34).
We implemented the baBIC method using BIC backward elimination as follows. The method started with all 39 (p) predictors and selected the subset of 38 (p-1) predictors with minimum baBIC. To select the subset of predictors with minimum baBIC, we fitted for each outcome all possible combinations of predictors obtained by removing 1 predictor at a time. We then computed the average of the nBICs across the 4 outcomes within each subset of predictors and selected the subset of 38 (p-1) with the minimum baBIC (Fig 1). In the next step of backward elimination, the method started with 38 (p-1) predictors and selected a subset of 37 (p-2) predictors that again rendered the minimum baBIC. The same process continued until there were only 2 variables left (i.e. "Male" and "Age decile groups"), which were forced in. Lastly, the method selected the final subset of predictors that had the minimum baBIC across all subsets of different number of predictors from p-1 to 2 (Fig. 2).
For the comparative methods, Individual Outcome and Union methods, we followed a similar approach as described above. The only difference being that the backward elimination was based on the minimum BIC for each individual outcome instead of the minimum baBIC across the 4 outcomes. We then obtained the Union model that contained all the predictors that were in at least 1 of the 4 best subsets of the Individual Outcome models.
For all final models, we computed the number of variables and measured predictive accuracy using the Harrell's C-statistic (35). For times to first ADL dependence, IADL difficulty, and mobility dependence, we used Wolbers et al. (36) adaptation of Harrell's C-statistic to the competing risks setting, where death status is switched to censored and the time-to-event is equal to the longest possible time-to-event that any respondent was followed up (i.e. 15 years). Of note, we obtained the same final subset of predictors for all selection methods with and without this simplification in the original case-study data.

Simulation study
Aim: To assess the performance of the proposed baBIC method in the selection of a common subset of variables to predict multiple outcomes with accuracy and parsimony.

Data-generating mechanisms:
We considered two data-generating mechanisms or scenarios, and within each scenario we simulated 4 survival times with high and low correlation among the outcomes. The 4 outcomes in the HRS data were highly correlated based on the pairwise Pearson correlations (range: 0.80-0.91). Thus, to test whether the correlation among the outcomes impacted the selection methods, we generated survival times with high and low correlation. For both scenarios, the data were simulated on 5,531 respondents.
To model the relationship between predictors and outcome, in scenario 1 for all 4 outcomes we used the results from the fitted Cox proportional hazards regression models with the common subset of predictors obtained in the baBIC model for the original HRS data. In scenario 2, we used the results from the fitted Cox models using the corresponding set of predictors obtained in each of the Individual Outcome models for the original HRS data. Additional File 2 shows the relationships between predictors and outcomes under both scenarios. Scenario 1 is the most unrealistic of the two since we assumed that only the common set of predictors included in the baBIC model from the original HRS data could be used to simulate the relationship between the predictors and all 4 outcomes. In other words, the baBIC model is the "correct" underlying model for all 4 outcomes. On the other hand, in scenario 2 we allowed an individual set of predictors to be used during the simulation of each outcome. Thus, scenario 2 captured what we feel to be a more typical setting where the outcomes could share some of the predictors in the "correct" model, but not all of them.
In the simulation study, for times to first ADL dependence, IADL difficulty, and mobility dependence, we chose to fit Cox models instead of Competing-risk regression models due to the computation time constraint that represented fitting the latter. Thus, in order to do this, we used a modified version of the HRS data where those who died were treated as being censored at the longest possible time that any respondent was followed (i.e. 15 years) (36). As mentioned above, we obtained the same final subset of predictors for the original dataset in all selection methods with and without this simplification.
The simulated survival times of correlated outcomes were obtained as follows. First, we simulated 4variate normal random variables that had means of zero, standard deviations of 1, and the correlation structure from the HRS data. Next, we inverted the random values to probabilities. For uncorrelated outcomes, we used probabilities simulated from the uniform distribution. These probabilities were then used as look-up values in the observed time-to-event distributions for each of the outcomes. We More specifically, after fitting the Cox models for each scenario using the modified HRS data, we used the survival probabilities of each respondent at each time-S(t)-to create a coarse lookup table per respondent with the time-to-events corresponding to the 99 th , 98 th , and down to the 1 st percentile in the survival curve. If the respondent survival curve stopped at an S(t) greater than the remaining percentiles in the survival curve, the respondent was censored at their last observed time-to-event.
Finally, we used the simulated probabilities from correlated and uncorrelated outcomes to select the matched survival probability S(t) and corresponding time-to-event for each respondent.
We generated 500 simulations for each scenario with correlated and uncorrelated outcomes. This number of simulations gave us a good balance between feasible computing time and acceptably small Monte Carlo Standard Errors (SEs). After obtaining the simulated survival times, each simulated outcome dataset was merged with the 39 predictors from the original HRS dataset.

Methods:
For each simulated data we obtained the BIC, nBIC, stepwise baBIC, and the final baBIC, Individual Outcome, and Union models. Individual Outcome models. By contrast, the baBIC model with 15 predictors was more parsimonious than the Union model, and all the predictors selected in the baBIC model were also present in the Union Model. These results were also confirmed in the simulation study. In both scenarios with correlated and uncorrelated outcomes, the Union model had on average more predictors than the baBIC model (Table 1). The difference in the numbers of predictors between Union and baBIC models was more subtle in scenario 1 where the simulated survival times were generated using the common set of predictors from the baBIC model in the HRS data (correlated outcomes Union model 15 In the HRS data and simulations with correlated and uncorrelated outcomes, the C-statistics of the Individual Outcome, Union, and baBIC models were clinically similar within each outcome. The average C-statistics across outcomes of the Union and the baBIC models were also comparable in the  Table 2).
As shown in Table 1, the average number of predictors of the baBIC models obtained across both simulation scenarios with correlated and uncorrelated outcomes were slightly smaller than the 15 predictors obtained in the baBIC model of the HRS data. However, the average C-statistics of the baBIC models of both simulation scenarios with correlated and uncorrelated outcomes were very similar to the C-statistic of the baBIC model of the HRS data (e.g. scenario 2, correlated outcomes: 0.664 [SE: 0.0002] vs. HRS data: 0.657) ( Table 2).
When using the baBIC method in the simulations, most of the predictors present in the baBIC model of the HRS data were correctly identified 80% of the times or more. On average in scenario 1, this method selected the same (correct) predictor as in the HRS data 85.7% of the times for correlated outcomes and 88.1% of the times for uncorrelated outcomes. In simulation scenario 2, the baBIC method selected on average the correct predictors 79.1% of the times for correlated outcomes and 82.4% of the times for uncorrelated outcomes. In the simulations, the number of correct predictors in the baBIC models ranged from 9 to 15 (15 being the maximum possible), and the percentage of models with 13 or more predictors ranged from 42.8% (scenario 2, correlated outcome) to 89.4% (scenario 1, uncorrelated outcomes). In Scenario 2, the predictors selected in the Union model, but not in the baBIC model of the HRS data, were included in the simulations less than 25% of the times, whereas these predictors were almost never present in simulations of scenario 1. Finally, the percentage of predictors not included in either the baBIC or Union models of HRS but present in the baBIC models of the simulations was less than 6% across scenarios (Table 3).

Discussion
The baBIC selection method produced a model with a good balance between parsimony and predictive accuracy. In both the HRS data and the simulations, this model was more parsimonious than the Union model, and it showed minimal loss of predictive discrimination. A good compromise between parsimony and accuracy is important since models that are simpler to understand and explain and that predict outcomes well are more likely to be implemented. Models with too few predictors cannot adequately describe the relationship between outcomes and predictors, whereas those with too many predictors can cause overfitting problems. Moreover, as the number of predictors in the model increases, the time and cost of collecting them could also increase. From a practical perspective, busy clinicians are unlikely to use a prognostic model with a daunting list of predictors to collect and enter. Although we did not formally incorporate a penalization associated with the cost of the predictors, other authors have explicitly balanced predictive accuracy against cost of the predictors (12).
In scenario 1, where the simulated survival times were generated using only the common set of predictors from the original baBIC model, the simulated baBIC models were still more parsimonious than the simulated Union models (by about 2 predictors on average). The selection method intrinsically favored the predictors that were used during the data-generating mechanism.
Consequently, during the individual selection process, the 4 outcomes ended up having more 13 common predictors, which in turn reduced the overall number of predictors in the Union model. On the other hand, scenario 2 assumed that each outcome had an individual best set of predictors, markedly increasing the number of predictors in the simulated Union models while maintaining comparable parsimony in the simulated baBIC models to those from scenario 1.
In the simulations, we found that the baBIC method performed well by selecting on average a high percentage of the predictors included in the final baBIC model of the HRS data (i.e. high percentage of correct inclusion), while keeping a low percentage of the predictors that were not in the baBIC model (i.e. high percentage of correct exclusion). In simulation scenario 2, the predictors included in the Union model, but not in the baBIC model of the HRS data showed up in the simulations less than 25% of the times. This was expected since these predictors were used during the simulation of the outcomes, and consequently they were partially favored during the selection method.
The average number of predictors of the final models of the simulations was slightly smaller than that of the final models of the HRS data. This could be explained because our attempt to replicate the structure of the HRS data did not entirely capture its complexity. Despite this, the average C-statistic of the final models of the simulations and the HRS data were clinically similar.
Breiman and Friedman (37) considered the relationship between the outcomes to improve predictive accuracy. In this way, several studies have developed methods for variable selection explicitly accounting for the correlation among multiple outcomes (10,11,13,38,39). In our method, we did not include the between-outcomes correlation. However, by averaging the normalized BIC across outcomes and selecting the best subset of predictors based on the minimum average normalized BIC, we pooled evidence across outcomes and implicitly incorporated their relationships. Furthermore, in the simulation study, we found that the percentage of correct inclusion of predictors (compared with the predictors selected in the baBIC model of the HRS data) was similar within the same simulation scenario with correlated and uncorrelated outcomes. These findings suggest that during variable selection a similar subset of predictors can be obtained regardless of the correlation among the outcomes.
As noted, several studies have used penalized regression under the high dimensional multivariate regression setting, where the numbers of predictors and outcomes may be large compared to the sample size. Regularization methods are particularly suitable for the study of genetic pathways or genome-wide association analysis, where high dimension, low sample size settings are very common (10,15,38,39). In clinical settings, researchers are usually interested in interpretable effect estimates in addition to good predictive performance. Regression coefficients estimated by regularization schemes like those that are an extension of LASSO can be biased, making their interpretation more difficult (40). Furthermore, in our case-study data we obtained less parsimonious Individual and Union models with LASSO selection compared to backward elimination, while maintaining very similar predictive accuracy. Likewise, in an extensive simulation study, Hastie et al. (29) have recently noted that neither LASSO nor best subset selection nor stepwise regression was dominant across a variety of problem settings, and that no method had a large difference in variation explained. Thus, they suggested favoring methods that are easy to compute.
Consequently, we believe that in the clinical practice where the sample size is usually large compared with the number of outcomes and predictors, our baBIC method, which extends the use of popular (non-regularized) variable selection methods to the multivariate settings, has the benefit of easier implementation and interpretation as well as good predictive performance and parsimony.
As in our study, other authors have extended stepwise methods and criterion-based selection to multivariate settings (41-43). Our method differs in that we combined backward elimination with a criterion-based selection method like the BIC. By doing this, we improved computational efficiency by not fitting all possible models-as traditional criterion-based selection methods do-while maintaining predictive performance by using BIC instead of statistical significance (i.e. traditional stepwise methods), which does not always indicate predictive value (44).
More recently, a clinical study identified a common set of predictors across several adverse outcomes (45). The authors identified the predictors that were significantly associated with most or all the outcomes, one of them being a composite of the other outcomes. This method is simple to implement and allows optimizing clinical resources by focusing on a single-combined outcome. However, the authors relied on more ad-hoc strategy to identify the common set of predictors, whereas our approach focuses on selecting the best subset of predictors based on minimizing an extension of the BIC.
It is worth mentioning that our method focused on one of the first steps of building a regression model. That is, we aimed to select a common subset of variables from a pool of many available predictors rather than identify a final predictive model. Thus, we assumed that all aspects of model building are fixed, except the selection of the predictors. In the actual application of this method, researchers will need to consider the rest of the aspects involved in model building; for example, possible inclusion of non-linear terms, interaction and multicollinearity between predictors, and for survival models, validity of the proportional hazard assumption. Additionally, it will be important to assess the performance of the final model using both calibration and discrimination techniques, as well as conducting model validation by internal cross validation (bootstrapping) and external validation. In a real life application, our method could be fully incorporated during the process of model development and validation.

Conclusions
Our baBIC method implemented a straightforward approach to obtain a common set of variables for the prediction of several outcomes. By selecting a common set of predictors for multiple clinical outcomes, researchers will be able to build prognostic models that are both accurate and parsimonious, potentially saving the clinical time and expense associated with gathering additional unnecessary predictors. Although the method shown here was developed for survival data using BIC as a convenient statistic for selection, it could easily be extended to generalized linear models or other information criteria such as the AIC. Moreover, this method can potentially be applied to larger data with a greater number of predictors and outcomes. And under this setting, it would be worth exploring further the implementation of the baBIC method into the LASSO framework. As the number of predictors and outcomes increases, there would be some computational challenges particularly for Competing-risk survival models that have longer run time than other regression models.

Availability of data and materials
All data generated or analyzed during this study, and the codes used to generate the results of this article are included within the article and its additional files.

Competing interests
The authors declare that they have no competing interests.

Funding
This work was supported by the National Institute on Aging (grant numbers: R01 AG047897, R01 AG057751). The funding agency did not participate in the design of the study, collection, analysis, interpretation of data, or in writing the manuscript.

Authors' contributions
LGDR drafted the manuscript, had full access to all the data in the study and took responsibility for the integrity of the data and the accuracy of the data analysis. SJL, AKS, and WJB designed the study and directed its implementation. SG helped on the acquisition and analysis of the data. All authors read and approved the final manuscript.     are the BICs of the full and best individual models.

Additional File 9
File format: .txt Description: SAS code to generate simulated datasets with correlated outcomes and predictors from Best Individual models and compute some statistics (Scenario 2, correlated outcomes)

Additional File 10
File format: .txt Description: SAS code to generate simulated datasets with correlated outcomes and predictors from baBIC model and compute some statistics (Scenario 1, correlated outcomes)

Additional File 11
File format: .txt Description: SAS code to generate simulated datasets with correlated outcomes and predictors from Best Individual models and compute some statistics (Scenario 2, uncorrelated outcomes)

Additional File 12
File format: .txt 33 Description: SAS code to generate simulated datasets with correlated outcomes and predictors from baBIC model and compute some statistics (Scenario 1, uncorrelated outcomes)

Additional File 13
File format: .txt Description: SAS code to perform BIC backward elimination for individual outcomes for Scenario 1 simulated correlated data and to obtain Union model.

Additional File 14
File format: .txt Description: SAS code to perform BIC backward elimination using baBIC for Scenario 1 simulated correlated data

Additional File 15
File format: .txt Description: SAS code to perform BIC backward elimination for individual outcomes for Scenario 1 simulated uncorrelated data and to obtain Union model.

Additional File 16
File format: .txt Description: SAS code to perform BIC backward elimination using baBIC for Scenario 1 simulated uncorrelated data

Additional File 17
File format: .txt Description: SAS code to perform BIC backward elimination for individual outcomes for Scenario 2 simulated correlated data and to obtain Union model.