Statistical methods

We used BMA methods to identify combinations of biomarkers and CPB time. All biomarkers were log-transformed and CPB time was included as a linear term. Urine biomarkers were not normalized to urine creatinine, though urine creatinine was included as a candidate predictor.

BMA involves assigning each variable a prior probability of being useful for prediction; these prior variable probabilities induce a prior probability for each combination, where the combinations are defined by allowing CPB time and each biomarker to either be included or excluded. The method combines these prior probabilities and the data via Bayes’ theorem to calculate a posterior probability for each combination (“posterior model probability”) and a posterior probability for each variable (“posterior variable probability”) [19–21]. The posterior model probability is a measure of the degree to which the model is supported by the data [22]. Similarly, the posterior variable probability reflects the support in the data for the variable as a predictor of the outcome [23]. The BMA framework can be used for variable selection on the basis of posterior model probabilities or posterior variable probabilities. The BMA approach considers all possible combinations and applies a “leaps and bounds” algorithm to identify the most promising combinations for further consideration; this process provides computational feasibility for searching the large space of candidate models (8,388,608 candidate models given 23 candidate predictors) [19].

In our implementation of BMA, we assigned each biomarker and CPB time a prior probability of ½, meaning that each predictor was a priori as likely to be in the model as not. It is possible to incorporate prior information into these probabilities, but we elected to treat all of the candidate predictors equally. These prior probabilities yield a prior probability for each combination of (0.5)²³ = 1.19 × 10⁻⁷, as there are 23 candidate variables. To account for possible center differences, we considered center-adjusted combinations by forcing center to be included in each combination evaluated by BMA. We pre-specified to select two combinations on the basis of the BMA analysis: (1) the maximum posterior model probability combination (the combination with the highest posterior model probability) and (2) the median probability combination (the combination consisting of all predictors with posterior variable probability exceeding 50%) [21].

We applied BMA to our data to develop combinations for predicting sustained mild AKI. After identifying the maximum posterior model probability combination and the median probability combination, we fit a center-adjusted logistic regression to the biomarkers included in these combinations, with sustained mild AKI as the outcome. Using the estimates from these regressions, we estimated the center-adjusted and optimism-corrected area under the receiver operating characteristic curve (AUC) of each combination for sustained mild AKI and for our secondary outcome, severe AKI. First, we estimated the apparent center-adjusted AUC for each combination and each outcome [24]. Then, we estimated the optimism in the center-adjusted AUC for each combination and each outcome using a bootstrapping procedure with 1000 replications [25]. In each bootstrap sample we repeated the entire model selection process. We subtracted the average optimism across bootstrap datasets from the apparent center-adjusted AUC to estimate the center-adjusted and optimism-corrected AUC. Figure 1 describes the analysis in detail. Importantly, this approach addresses model selection bias, resubstitution bias, and potential bias due to center differences [26]. We emphasize that without optimism correction, estimated AUCs will tend to be overestimated due to both resubstitution bias (i.e., using the same data to develop and evaluate a combination) and model selection bias (i.e., using the data to select the model). By accounting for these sources of optimistic bias, we have a more realistic assessment of how the combinations may perform in independent data. This procedure does not supplant external validation. Rather, this is a form of internal validation where the full dataset is used to fit the combination and estimate its apparent performance, followed by bootstrapping to quantify the optimistic bias in the apparent performance. Confidence intervals (CIs) were estimated for the center-adjusted and optimism-corrected AUC by bootstrapping the BMA procedure and obtaining a 95% CI for the apparent center-adjusted AUC, and then shifting the confidence interval by the average optimism.

Analysis flow. Legend: Abbreviations: AKI = acute kidney injury; BMA = Bayesian model averaging; AUC = area under the receiver operating characteristic curve

Our primary measure of model performance was the AUC, which measures how well a combination discriminates cases from controls. We acknowledge the limitations of the AUC and that it represents an incomplete assessment. Our goal was to propose combinations with high prognostic capacity and the potential to be developed into useful risk prediction models, and we were particularly concerned with avoiding common sources of bias in identifying prognostic combinations, including possible center differences [27]. The adjustment for center does not allow for individual predicted risks. Therefore, we do not assess model calibration in this work, as we do not propose risk prediction models. However, if these combinations are later developed into risk prediction models, an assessment of calibration will be required.

We considered several model diagnostics, including the posterior model probability of the selected combinations across bootstrap samples, the posterior variable probability of each predictor across bootstrap samples, the posterior variable probability of each predictor omitting each observation in turn, and the performance of the estimated selected combinations across bootstrap samples.

In an exploratory analysis, we compared the performance of the BMA procedure to two common variable selection methods: forward selection and univariate selection. The following algorithm was used to compare the three methods. We randomly split the data into training and test datasets of equal size with equal numbers of sustained mild AKI cases. We then applied each of the three model selection methods to the training data. First, we applied BMA and identified the maximum posterior model probability combination and the median probability combination. Second, we applied forward selection with a p-value threshold of 0.1. Third, we applied univariate selection, forming a combination of all variables with a p-value less than 0.1. All methods used center-adjustment. In each iteration we applied the resulting combinations to the test data and estimated the center-adjusted AUC for the combination using the test data only; thus, we performed internal validation whereby the training dataset was used for fitting while the test dataset was held out for evaluation. We repeated this procedure 1000 times, independently randomly splitting the data into training and test datasets each time. We calculated 95% intervals as the 2.5th and 97.5th percentiles of the AUC across these 1000 replications.

As a secondary analysis, we evaluated the association of the biomarker combinations identified by the BMA methods with death at 1 year and 3 years after surgery. For each biomarker combination and each time point (1 year and 3 years), we fit a logistic regression model with the fixed estimated biomarker combination, adjusting for center. We used the full dataset to estimate the odds ratio describing the association between the combination and death. We can consider the two estimated combinations, M₁ and M₂, where M₁ has p variables (denoted by X), combined via the parameters β₁, …, β₂, and M₂ has q variables (denoted by Y), combined via the parameters α₁, …, α_q:

The odds ratio for the association between the combination and death was estimated by fitting two logistic regressions for each time point (1 year and 3 years):

where δ0C and θ0C are center-specific intercepts.

All analyses were completed using R 3.1.2. The BMA package in R was used for the BMA analyses [28]. The R code for the primary analysis is provided in (Additional file 1: Item S1) and at https://github.com/allisonmeisner/BMAbiomarkers.