Risk Prediction Using Bayesian Networks: An Immunotherapy Case Study in Patients With Metastatic Renal Cell Carcinoma.

PURPOSE: To address the need for more accurate risk stratification models for cancer immuno-oncology, this study aimed to develop a machine-learned Bayesian network model (BNM) for predicting outcomes in patients with metastatic renal cell carcinoma (mRCC) being treated with immunotherapy.
METHODS: Patient-level data from the randomized, phase III CheckMate 025 clinical trial comparing nivolumab with everolimus for second-line treatment in patients with mRCC were used to develop the BNM. Outcomes of interest were overall survival (OS), all-cause adverse events, and treatment-related adverse events (TRAE) over 36 months after treatment initiation. External validation of the model's predictions for OS was conducted using data from select centers from the International Metastatic Renal Cell Carcinoma Database Consortium (IMDC).
RESULTS: Areas under the receiver operating characteristic curve (AUCs) for BNM-based classification of OS using baseline data were 0.74, 0.71, and 0.68 over months 12, 24, and 36, respectively. AUC for OS at 12 months increased to 0.86 when treatment response and progression status in year 1 were included as predictors; progression and response at 12 months were highly prognostic of all outcomes over the 36-month period. AUCs for adverse events and treatment-related adverse events were approximately 0.6 at 12 months but increased to approximately 0.7 by 36 months. Sensitivity analysis comparing the BNM with machine learning classifiers showed comparable performance. Test AUC on IMDC data for 12-month OS was 0.71 despite several variable imbalances. Notably, the BNM outperformed the IMDC risk score alone.
CONCLUSION: The validated BNM performed well at prediction using baseline data, particularly with the inclusion of response and progression at 12 months. Additionally, the results suggest that 12 months of follow-up data alone may be sufficient to inform long-term survival projections in patients with mRCC.

INTRODUCTION

Immunotherapy with immune checkpoint inhibitors (ICIs) has shifted the paradigm of cancer therapy. By activating the body's immune response against cancer cells, ICIs have been shown to be particularly effective in advanced malignancies at producing durable responses in patients.[1,2] ICIs are also associated with fewer side effects because of their selectivity for tumor cells rather than indiscriminate cytotoxicity.[3-5] Despite their therapeutic promise across multiple cancer indications, studies show significant variability in individual-level response to immunotherapy drugs, driving a need to accurately identify treatment responders.[1,6] By identifying prognostic variables, markers of survival and adverse events (AE), and/or risk stratification criteria, therapies could be individually tailored to patients to maximize therapeutic benefit and improve cost-effectiveness of innovative therapies.

CONTEXT

Key Objective

To develop an interpretable model for predicting survival and adverse events in patients with metastatic renal cell carcinoma (mRCC) treated with immunotherapy.

Knowledge Generated

Bayesian network–based multivariate individual-level prediction model evaluated on face validity, performance against machine learning classifiers, and real-world generalizability. We identify prognostic variables for survival and safety outcomes and show that tumor response within the first year of initiation of nivolumab for second-line therapy is highly prognostic of long-term outcomes in mRCC.

Relevance

Because of individual-level heterogeneity in response to treatment, many novel cancer therapies are efficacious in only a subset of patients receiving treatment. Our findings indicate that individual-level prediction of survival and adverse events in patients with mRCC treated with immunotherapy may be tractable at the time of treatment initiation or shortly thereafter.

Machine learning methods are increasingly important for prognostic modeling and knowledge discovery in clinical research, offering an opportunity to deliver on the promise of precision medicine.[7] Clinical trials and real-world data sets routinely capture many variables on patient demographics, treatment regimens, and tumor molecular markers among others that could be used to assess optimal targets for treatments. However, traditional methods of analysis typically use only a small set of predefined variables to inform risk and survival predictions. On the other hand, machine learning methods leverage the full suite of available patient data for individual-level risk estimation.[8] Potential applications of this type of analysis include informing adaptive clinical trial design, identifying predictors of surrogate end points, and supporting extrapolation of trial outcomes and clinical decision making. Yet, very few studies have used this more promising approach as a result of concerns about interpretability, limited clinician input, missing observations, and challenges with concurrently modeling several correlated outcomes (ie, multivariate predictions).[10-12]

This proof-of-concept study aimed to address these limitations of traditional modeling approaches by using machine learning to develop a multivariate Bayesian network model (BNM) for predicting survival and safety outcomes in patients with metastatic renal cell carcinoma (mRCC) being treated with ICIs. Currently, two validated prognostic models for targeted therapy in RCC are in widespread use—the Memorial Sloan Kettering Cancer Center (MSKCC) risk score (originally published in 1999[9]) and the International Metastatic Renal Cell Carcinoma Database Consortium (IMDC)/Heng risk score.[10] These models use Karnofsky performance status, time from diagnosis to systemic treatment, and three or four serum markers to stratify patients into favorable, intermediate, or poor prognostic risk groups. Although both models are continually updated to improve their prediction accuracies,[11] they were developed prior to the advent of ICIs for cancer therapy. Thus, they may exclude variables that predict safety, response, and/or survival while on immunotherapy.[16]

This study addresses these gaps by building a transparent risk prediction BNM that also offers insight into prognostically significant variables in cancer immunotherapy. The objectives of this study are as follows:

To develop an interpretable BNM for jointly predicting overall survival (OS), AE, and treatment-related adverse events (TRAE) within 3 years after initiating immunotherapy,

To assess the predictive performance of the BNM,

To validate the BNM to assess real-world performance and generalizability, and

To identify variables that predict OS, AE, and TRAE within 3 years after initiating immunotherapy (ie, prognostically significant variables).

METHODS

This article is concordant with Transparent reporting of a multivariable prediction model for individual prognosis or diagnosis guidelines.[12] Specific study components are detailed below. Additional details can be found in the Appendix.

CheckMate 025 Data Set

CheckMate 025 (CM-025) is an open-label, randomized, phase III trial that compares the safety and efficacy of nivolumab (a PD1 ICI) versus everolimus in patients with mRCC who had received prior antiangiogenic therapy. Details about CM-025 have been previously described.[13] Participants were at least 18 years old with advanced or metastatic clear cell renal cell carcinoma and had received prior antiangiogenic therapy. Key exclusion criteria were Karnofsky performance status of < 70, past or current CNS metastases, and previous treatment with mTOR (mammalian target of rapamycin) inhibitors or glucocorticoids.[3]

For purposes of this study, baseline data (eg, sociodemographic characteristics, tumor characteristics, and prior treatment types) and 3 years of follow-up data for patients who received treatment (eg, treatment response, disease progression, and subsequent therapies) were used to train the BNM (n = 803 with 75 variables). Outcomes of interest were OS at months 12, 24, and 36 (alive, dead, or missing or censored if lost to follow-up); AE in year 1, 2, and 3; and TRAE in year 1, 2, and 3. AE and TRAE were encoded as number of events of CTCAE grade 3 or above binned as one of 0, 1-2, 3-5, 6+. In all cases, we defined AE and TRAE as the number of AE and TRAE observed from the start of a year until death, loss to follow-up, or the end of the year, whichever happened first. A 3-year database lock that included a maximum follow-up of 5 years was used for the analysis.

IMDC Data Set

The BNM was validated with the latest cut of the IMDC data set containing observational data for patients with metastatic RCC who were receiving nivolumab or everolimus (n = 2,152).[10] Only IMDC centers with nivolumab access and sufficiently reliable data were included. Notably, key exclusion criteria from CM-025 did not apply to patients in IMDC data, and several variables showed moderate-to-high imbalance (see Results on external validation), which makes this a sufficiently distinct cohort for assessing broader model generalizability. However, only 26 (35%) of the baseline variables used for the development of the BNM were available in the IMDC data set, and data on AE were not available in IMDC data. Thus, only OS predictions (at 12 months) were used to validate the BNM to assess real-world generalizability.

RESULTS

Model Structure

The structure of the estimated BNM, showing nodes (ie, variables) and directed edges (ie, associations), is presented in Figure 1. To evaluate model interpretability, the face validity of the BNM structure was assessed by comparing the model against expected relationships. Specifically, the level of mutual information within clusters of biologically related variables in the model was examined and was determined to be high. These included clusters for glucose measures (serum glucose, fasting serum glucose, and total glucose), thyroid hormones (T3, T4, and thyroid-stimulating hormone), liver enzymes (AST and ALT), immune cells (lymphocytes, leukocytes, neutrophils, and platelets), and PD-1 and its ligand PD-L1. Similarly, summary metrics such as body mass index, health state utility values, and risk scores (ie, MSKCC and IMDC) were clustered with baseline variables used for their calculation (eg, height, weight, health-related quality-of-life scores, Karnofsky performance status, and hemoglobin, respectively) with a high level of correlation or mutual information. Thus, the model structure learned from CM-025 was consistent with expected relationships from literature. Markov blankets[14] of year-1 outcomes are shown in Figure A1. Variables in the Markov blanket, when nonmissing, are the only variables used for prediction as a result of conditional independence assumptions in BNs.

FIG 1.

Bayesian network model structure learned from CheckMate 025. Edge linewidth represents log-transformed mutual information between variables (thicker linewidth implies stronger statistical dependence). Disconnected nodes are not shown. Adv. events, Adverse events; CRF, case report form; Diag, Diagnosis; Rand. years, randomization (years); Rx, prescription; Sys. Th, Systemic therapy; Tx, therapy.

FIG A1.

Markov blankets for (A) overall survival, (B) adverse events, and (C) treatment-related adverse events in year 1. AE, adverse event; CRF, case report form; MSKCC, Memorial Sloan Kettering Cancer Center; Rx, prescription.

In addition to the graphical structure, we also examined the conditional probability tables underlying the graph. Conditional survival at 24 months, ie, the probability of survival conditional on having survived to month 12, was 65%, which is considerably higher than the marginal probability of survival to 24 months at baseline (46%) (Table 1). Similarly, probability of surviving to 36 months was 33% at baseline, increasing to 46% and 71% conditional on surviving to months 12 and 24, respectively (Table 1). The probability of AE decreased over time, going from 54% to 41% then to 15% conditional on surviving to months 12, 24, and 36, respectively. This indicates a substantial incremental improvement in survival and safety events over time for patients who were able to survive longer.

TABLE 1.

Conditional Probability of Survival at Months 12, 24, and 36

Classification Performance

Observed mean areas under the receiver operating characteristic curve (AUCs) for the BNM were generally either acceptable (0.7 < AUC < 0.8) or fair or poor (0.5 < AUC < 0.7). Specifically, the AUCs were 0.74, 0.71, and 0.68 for years 1, 2, and 3, respectively, for OS with event rates (% death) of 30%, 51%, and 62%; 0.60, 0.68, and 0.71 in years 1, 2, and 3, respectively, for AE; and 0.71, 0.64, and 0.72 in years 1, 2, and 3, respectively, for TRAE. However, the inclusion of response and progression as predictors of 12-month OS increased mean AUC to 0.85—ie, in the excellent category (0.8 < AUC < 0.9). Note that response and progression during years 1, 2, and 3 are a part of the model; the BNM allows flexibility in choosing whether we include these as predictors, outcomes, or neither depending on the analysis of interest.

Notably, similar mean AUCs were observed when alternative machine learning approaches were used to assess the relative performance of the estimated BNM (Table A1). The study found no significant difference on average in model performance across the nine outcomes when comparing the BNM against three machine learning classifiers: Lasso regularized logistic regression (Lasso), support vector machine, and random forest (RF). Comparing the BNM with gold-standard RF classifiers, RF did significantly better at classifying 12-month OS (P value .0029), whereas the BNM performed significantly better at classifying year-1 TRAE (P value < .001).

TABLE A1.

AUCs for Classification of OS, AE, and TRAE From Iterated 10-fold Cross-Validation

External Validation

Notable differences were observed between patients in the CM-025 and IMDC data sets. For instance, patients in the IMDC data set had worse Karnofsky performance status, poorer MSKCC risk scores, higher prevalence of bone metastases, lower weight or BMI, and lower survival rates and were more likely to be receiving everolimus. There were also notable imbalances (ie, standardized differences [d] > 0.25) in similar variables across both data sets (Fig A2). Overall, the IMDC data were found to be appropriate for assessing broader generalizability of the BNM.

FIG A2.

Standardized differences (d) between variables in International Metastatic Renal Cell Carcinoma Database Consortium and CheckMate 025. d < 0.1 implies no imbalance, whereas d > 0.25 implies considerable imbalance. CRF, case report form; MSKCC, Memorial Sloan Kettering Cancer Center.

Results from evaluating real-world performance of the BNM are presented in Figure 2 and Table 2. The estimated AUC for OS at 12 months was 0.71 (compared with 0.74 in CM-025). The study also compared the BNM against the IMDC/Heng risk score commonly used for risk stratification in patients with RCC undergoing systemic therapy. Mean AUCs for OS at 12 months using the BNM and IMDC/Heng risk score (available for 58% of patients in the IMDC data set) were 0.76 and 0.69, respectively. Conversely, in the subset of patients with missing IMDC risk score because of one or more missing risk criteria, externally validated AUC for OS was 0.64 for BNM. Therefore, the BNM outperforms the IMDC risk score for year-1 OS and can handle missing data for risk prediction. It is worth noting that probabilities were recalibrated to account for differences in 12-month mortality between CM-025 and IMDC data, and probability output showed < 10% deviance from ideal calibration in all cases (Fig 3).

FIG 2.

ROC curves for Bayesian network model for 12-month overall survival. CM-025, all available variables in the model; CM-025 (sub), common variables between CM-025 and IMDC from CM-025; IMDC, common variables between CM-025 and IMDC from IMDC. Red curves represent internal cross-validation, and blue curves represent external validation. Faint red lines represent ROC curves from individual runs of 10-fold cross-validation; the bold red line represents their mean. CM-025, CheckMate 025; IMDC, International Metastatic Renal Cell Carcinoma Database Consortium; ROC, receiver operating characteristic.

TABLE 2.

Areas Under the Receiver Operating Characteristic Curve for External Validation of BNM Trained Using CM-025 on IMDC Data

FIG 3.

Probability calibration performance. CM-025, all available variables in the model; CM-025 (sub), common variables between CM-025 and IMDC from CM-025; IMDC, common variables between CM-025 and IMDC from IMDC. CM-025, CheckMate 025; IMDC, International Metastatic Renal Cell Carcinoma Database Consortium.

Prognostic Variables

Prognostic variables were identified using their AUCs and are presented in Table 3. For OS, MSKCC risk score, IMDC/Heng risk score, and Karnofsky performance score were prognostic from months 12, 24, and 36 (univariable AUCs between 0.65 and 0.69; see also Tables A2 and A3). Functional Assessment of Cancer Therapy–Kidney Symptom Index–Disease-Related Symptoms[15] and health-related quality-of-life scores (EuroQol EQ-5D) at baseline were more highly prognostic of OS at months 24 and 36 (univariable AUC > 0.6) than month 12 (AUC < 0.6). Serum hemoglobin was highly prognostic, as expected. Interestingly, only 12-month follow-up related to progression, response, and subsequent radiotherapy and nonimmunotherapy was highly prognostic of OS over 3 years; this is possibly due to the majority of these occurring within the first 12 months after treatment initiation. For AE in year 1, sodium; hemoglobin; calcium; platelet count; and MSKCC, Karnofsky, and IMDC/Heng scores were prognostic. For AE in years 2 and 3, subsequent interventions (nonimmunotherapy, radiotherapy, and surgery) were prognostic. Similar to OS, progression and response in year 1 were prognostic of AE in years 1-3. Prognostic variables for TRAE included actual treatment type (everolimus or nivolumab), tumor size, bone and liver metastasis, and subsequent radiotherapy and immunotherapy. Consistent with the study findings for OS and AE, progression and response in year 1 were also prognostic of TRAE (Appendix).

TABLE 3.

Top 10 Prognostic Variables Ranked by AUC for Outcomes in Year 1

TABLE A2.

Top 10 Prognostic Variables Ranked by AUC for Outcomes in Year 2

TABLE A3.

Top 10 Prognostic Variables Ranked by AUC for Outcomes in Year 3

DISCUSSION

The study focused on developing and validating a proof-of-concept BNM for predicting OS, AE, and TRAE in patients with mRCC.[3] The study also identified prognostic variables that may be useful in clinical practice. The results showed that BNMs could serve as interpretable multivariate predictive models that perform on par with more commonly used black box machine learning approaches. An important differentiation between BNM and the other machine learning models we used is that missing data can be handled intrinsically by the model instead of requiring a prior distinct imputation step, although careful thought needs to be given to whether the mechanism of handling missingness is justified. In our case, CM-025 data contained negligible amounts of total missing data (approximately 3%[16]).

The ability of models to identify prognostic variables over different periods of time may be contingent on changing variability with time. Clinically, being able to predict a patient's 3-month OS is easier than predicting their 24-month OS. AUCs for OS exhibited a logarithmic decline over time (ie, by approximately 26% between years 1 and 2 and by approximately 4% between years 2 and 3). However, adding indicators of treatment response and disease progression in year 1 in the BNM increased the mean AUC for 12-month OS from 0.74 to 0.85. This suggests that the classification of long-term OS may be possible with fair accuracy if baseline data and 12-month follow-up are accounted for in the model. On the other hand, the AUCs for AE increased over time. This is consistent with the general observation that most AE occur within the first 12 months after randomization. It is possible that AE in years 2 and 3 may be present only in patients with poor baseline risk, making them easier to classify.

As previously noted, the model was validated with data from applicable IMDC centers. Comparison of variables between IMDC and CM-025 showed large imbalances in OS, Karnofsky performance status, bone metastases, and BMI, suggestive of a larger fraction of poor-risk individuals in IMDC, and by extension the real world. Notably, prognostic variables such as health-related quality-of-life scores were not available in IMDC. However, the data set was found to be suitable for testing the BNM. The analysis indicated that the overall generalizability of the model was quite high–external validation AUC of 0.71 showed a 4% decrease from cross-validation AUC on CM-025 for OS.

This study has some limitations. First, one primary concern about machine learning approaches is the reliance on black box models that sacrifice transparency and interpretability for greater prediction accuracy.[12] Model transparency and interpretability are essential to clinical research, and domain knowledge is critical for the development and evaluation of patients' risk stratification models. Causal constraints were applied during the final step of machine learning as this had the potential to reduce predictive accuracy in the estimated BNM. Second, although the model was validated across clinically important subgroups, with largely consistent overall validation, patients age at least 75 years or with favorable IMDC risk at baseline had relatively poor external AUCs. The poor generalizability of the model to these two subgroups is possibly due to variable imbalances between CM-025 and IMDC or overfitting and may be explored in future research. Third, we discretized survival outcomes, which may lead to a loss of information. Approaches to BN learning using weighted censored instances may be useful for modeling survival as a censored time-to-event outcome.[17] Our choice to discretize continuous predictors was to avoid distributional assumptions, similar to prior work with BNs.[18-23] Although this may not be ideal because of loss of information, it was negligible for prediction in our case as BNM performance was on par with the other machine learning methods where covariates were not discretized. Last, different types of serious AE were grouped together under two categories (AE and TRAE). However, it is plausible that different sets of prognostic variables apply to different types of AE.

In conclusion, the study developed a proof-of-concept BNM for predicting OS, AE, and TRAE using CM-025 data. The model was then validated using IMDC data, and the performance across clinically important subgroups was highlighted. The practical utility of BNM for interpretable machine learning analysis, multivariate prediction, and handling missing data was also addressed. The results indicate that the model performs on par with black box machine learning models and outperforms IMDC/Heng risk score for classifying OS (with a 10% difference in external AUC for patients with all IMDC risk factors available and a 28% difference in patients with one or more missing risk criteria—assuming that IMDC/Heng risk score cannot be calculated). Future work could extend the BNM to allow for extrapolations beyond the 3-year horizon as the results suggest that long-term predictions with information of 12-month follow-up may be possible. The model could also be used as a clinical risk prediction tool for supporting clinical decision making after prospective validation to determine feasibility for use in the clinic, reproducibility, and accuracy.[19,24] Although combination therapy with nivolumab in the first-line setting is increasingly becoming common, this model would be useful in regions where these combinations are not available or funded and where nivolumab is used for second-line therapy. Nonetheless, we anticipate increased use of transparent machine learning models for leveraging bigger data sets to inform health policy and decision making for innovative therapies in the future.