Use of machine learning techniques to predict poor survival after hematopoietic cell transplantation for myelofibrosis Open Access

Myelofibrosis (MF) is a chronic myeloproliferative neoplasm that appears de novo (primary MF [PMF]) or after a diagnosis of essential thrombocythemia or polycythemia vera (secondary MF [SMF]). Managing MF is complex because of its diverse clinical manifestations including an inherent risk of progression to acute myeloid leukemia (AML). Although the median overall survival (OS) is ∼6 years, it varies significantly among patients. Medical treatment focuses on symptom control and quality of life but is not curative and does not reduce the risk of AML.^1,2 

Allogeneic hematopoietic cell transplantation (allo-HCT) remains the only curative option for MF.³ However, its significant morbidity and mortality require a careful risk-benefit analysis to identify appropriate candidates. This has become particularly critical recently, as several effective therapies for MF have been incorporated into clinical practice,^4-6 with others showing promising results in clinical trials.⁴ Existing prognostic models for OS after allo-HCT have been instrumental in guiding clinical decisions.^7,8 However, they do not account for key factors such as patient comorbidities⁷ and emerging transplant strategies, such as haploidentical transplants or posttransplant cyclophosphamide use.^9,10 Furthermore, there is significant room for improving their ability to accurately identify patients at high risk of posttransplant mortality, who may benefit more from alternative treatments or clinical trials.

Machine learning (ML) is a field of artificial intelligence in which prediction is based on modeling of outcomes considering the complex interactions among multiple variables, rather than relying on predefined human-made rules. These techniques have demonstrated their utility to provide accurate personalized survival predictions for patients with MF undergoing conventional drug treatment.^11,12 In this study, we aim to assess whether ML techniques can similarly improve prognostication of OS in the setting of allo-HCT for MF. The ultimate goal is to enhance transplant decision-making by providing more precise and individualized survival predictions.

We retrieved data from adult patients with MF (PMF or SMF) who underwent first allo-HCT between 2005 and 2020 in European Society for Blood and Marrow Transplantation (EBMT) centers. Patients who received transplantation from umbilical cord blood and those with AML transformation were excluded. Myeloablative or reduced intensity conditioning regimens were defined by standard EBMT criteria.¹³ Matched unrelated transplants were matched at allele level for HLA-A, -B, -C, -DRB1, and -DQB1. Centers ensured informed consent in compliance with local regulations to report pseudonymized data to the EBMT. The study was approved by the Chronic Malignancies Working Party of EBMT and conducted in accordance with the Declaration of Helsinki. Informed consent for inclusion in the EBMT registry was obtained for all patients. The study database included 52 variables selected for their prognostic significance based on previous studies.^7,14,15 

The primary goal was to develop a prognostic model for OS using ML techniques and compare its performance with that of a Cox regression–based score developed in the same data set, and with the Center for International Blood and Marrow Transplant Research (CIBMTR) model.⁸ The models were also applied for predicting the secondary outcome of nonrelapse mortality (NRM). Progression-free survival (PFS) and cumulative incidence of relapse were estimated for descriptive purposes. Median follow-up was determined using the reverse Kaplan-Meier method. In patients who died after disease relapse, relapse was considered the primary cause of death.¹⁶ 

OS and PFS were estimated by the Kaplan-Meier method. NRM was defined as the time from the date of transplantation to the date of death (uncensored) or to the date of disease relapse (censored). The cumulative incidences of relapse/NRM (as competing risk for each other) were analyzed separately in a competing risks framework.¹⁷ 

Two independent statisticians applied distinct methodologies to evaluate the factors influencing OS. One used a conventional multivariate Cox proportional hazards regression model, whereas the other used a range of ML techniques. Both approaches were based on the same random distribution of the patient cohort into a training set (75% of the cohort, n = 3887) and a test set (25% of the cohort, n = 1296). Each statistician independently determined the optimal cutoff points for the newly derived risk scores to stratify patients into distinct prognostic groups. The resulting risk classifications were subsequently compared and contextualized to assess their clinical relevance for treatment decision-making.

Statistical analyses were performed in R version 4.3,¹⁸ with the “survival,”¹⁹ “prodlim,”²⁰ and “cmprsk”²¹ packages.

Factors potentially associated with OS were entered into a Cox proportional hazards model. Selection of variables in the final model was based on expert clinical advice and data availability, to assess the independent effect of each covariate. Variables with a high degree of missingness (defined as >50%) were not considered for inclusion in the model, whereas a missing category was created for those variables with a low degree of missingness. Hazard ratios (HRs) were provided, and corresponding P values were calculated using the Wald test. A score of 1 or 2 was assigned to each significant variable for OS based on the HRs obtained from multivariable analysis. The cutoff values were arbitrarily defined as follows: HR of <1.25 = 0 points, HR of 1.25 to 1.50 = 1 point, HR of >1.50 = 2 points. A prognostic scoring system was subsequently developed considering the sum of risk points to discriminate 4 patient risk groups with significant differences in OS. All P values were 2-sided and P < .05 was considered significant.

Random survival forests (RSF) were created with 1000 trees. For cross-validation, sampling was performed without replacement, which, by default, takes 0.632 times the sample size. Missing variables were imputed in the training and test cohorts separately using a missing data algorithm developed by Ishwaran et al.²² Predictions were cross-validated in the training set (using the out-of-bag method) and then validated in the test cohort. This was done to rule out overfitting of performance metrics in the training set related to either variable selection or the imputation process. Dimensionality reduction was performed using variable-importance estimation analysis. Redundant and dependent variables were discarded applying variable-importance estimations and clinical knowledge, achieving a minimally dimensioned yet effective model.

Hyperparameter tuning was explored to optimize the performance of the RSF model. Specifically, a grid search was used to explore the impact of key hyperparameters, including “mtry” (the number of variables randomly selected at each split) and “node size” (the minimum number of samples required in terminal nodes). Partial dependence plots were used to evaluate the marginal effect of individual covariates on survival estimates, while accounting for the average influence of all other variables. These plots help visualize nonlinear relationships and potential threshold effects, providing insight into the contribution of specific predictors to survival probability over time. This approach enhances the model’s interpretability by isolating the independent effect of each covariate on the predicted survival outcomes.

In addition to the baseline RSF model, we evaluated the performance of 3 complementary methods using the 10 variables included in the final model: oblique RSF (ORSF), gradient-boosted survival trees using XGBoost, and a deep neural network–based survival model (DeepSurv).^23-25 

A detailed explanation of the different ML techniques is outlined in the supplemental Data, available on the Blood website. All 3 approaches, ORSF, XGBoost-based survival modeling, and DeepSurv, were evaluated using cross-validation for performance estimation. Early stopping was applied when computationally practical to limit overfitting. The performance of each model was assessed using standard survival metrics (eg, C-index) to facilitate a rigorous comparative analysis.

The discriminative capacity of the ML and Cox OS models was compared in the training and test sets using the Harrell C-index.²⁶ Time-dependent receiver operating characteristic (ROC) area under the curves (AUCs) for OS were calculated using the timeROC package.²⁷ Quantitative scores (continuous risk estimates) and categorical risk groups were analyzed separately for each model. For risk groups, categorical labels were converted into numeric values for compatibility. Time-dependent ROC-AUCs assessed model performance at specific time points, allowing a robust evaluation of prognostic accuracy over time of the RSF and Cox models in terms of their quantitative and group-based predictions.

Because the C-index is not an optimal metric for competing risk models, we also included the Akaike information criterion (AIC) scores for both NRM and OS.²⁸ 

A total of 5183 patients from 288 centers fulfilling the selection criteria were included. Baseline characteristics of all patients, along with the training (n = 3887) and validation (n = 1296) cohorts, are presented in Table 1. Median follow-up was 58.2 months (95% confidence interval [CI], 55.6-59.8) in the training set and 60.0 months (95% CI, 55.7-63.2) in the test set. Median OS was 79.4 months (95% CI, 69.2-89.6) in the training set and 73.7 months (95% CI, 54.7-92.7) in the test set. No significant differences in characteristics were observed between the cohorts, apart from a higher platelet count at allo-HCT and less antithymocyte globulin use in the test cohort.

The estimated OS rate at 1, 5, and 10 years was 70% (95% CI, 69-71), 53% (95% CI, 51-54), and 43% (95% CI, 41-45), respectively.

The probability of PFS after 1, 5, and 10 years was 62% (95% CI, 60-63), 44% (95% CI, 43-46), and 35% (95% CI, 33-37), respectively. The estimated NRM rate at 1, 5, and 10 years was 23% (95% CI, 22-24), 32% (95% CI, 31-33), and 36% (95% CI, 35-38), respectively. Cumulative incidence of relapse at 1, 5, and 10 years was 15% (95% CI, 14-16), 24% (95% CI, 23-25), and 29% (95% CI, 27-31), respectively.

The graphical representation of the main study outcomes in the training and test cohorts can be seen in supplemental Figure 1.

Factors associated with OS in the univariable analysis are shown in supplemental Table 1. In the multivariable analysis, 7 independent factors significantly predicted reduced OS: older patient age, HLA-mismatched donor type, lower Karnofsky performance status (KPS), higher HCT-specific comorbidity index (HCT-CI), JAK2/triple-negative genotype, graft from a female donor to a male patient, and graft from a donor who is cytomegalovirus (CMV) positive to a recipient who is CMV negative (Table 2). Based on the HRs, a score of 2 was assigned to patient age of ≥60 years and the JAK2/triple-negative genotype; and a score of 1 to patient age of 50 to 59 years, haploidentical or mismatched unrelated donors, KPS of <90, and HCT-CI of ≥3. Because the HRs for sex match and CMV serostatus were only modestly increased, no score was assigned to these factors. The total score ranged from 0 to 7 points, with 4 risk categories: low risk (0-1 points), intermediate-1 risk (2-3 points), intermediate-2 risk (4-5 points), and high risk (6-7 points). The corresponding 5-year OS of each category in the training and test set were 82% (95% CI, 74-90) and 65% (95% CI, 41-89) for low risk (6% of the cohort); 62% (95% CI, 58-66) and 65% (IC 95%, 57-73) for intermediate-1 (36% of the cohort); 52% (95% CI, 48-56) and 47% (95% CI, 40-54) for intermediate-2 (48% of the cohort); and 39% (95% CI, 30-47) and 30% (95% CI, 13-47) for high risk (10% of the cohort), respectively (Figure 1A-B).

A RSF model was created to predict OS using the 52 initial variables of the data set. This model achieved a C-index of 0.603 in the training set and 0.632 in the test set. The variable-importance metrics for the model in the training set are shown in supplemental Figure 2. After dimensionality reduction, the model was refined to a smaller set of key prognostic variables: patient age, HCT-CI, KPS, blood blasts percentage, hemoglobin level, leukocyte and platelet counts, donor type, conditioning intensity, and graft-versus-host disease (GVHD) prophylaxis. This model achieved a C-index of 0.599 in the training set and 0.623 in the test set. Despite having performed hyperparameter tuning, no improvement in C-index was achieved (supplemental Table 2). To further elucidate these relationships, we have also included partial dependence plots in supplemental Figure 3.

As shown in supplemental Table 3, RSF achieved higher concordance indices for OS and NRM predictions in both training and test sets compared with 3 alternative methods (ORSF, DeepSurv, and XGBoost). The consistent and superior performance of RSF across both data partitions justified its selection as the primary approach for downstream analyses.

This analysis was performed on the subset of patients who had complete information on the variables included in the Cox-derived score to minimize biases (training set: n = 1773; test set: n = 566). In mortality prediction, the ML model demonstrated modestly better performance in the training set, achieving a C-index of 0.603 compared with the Cox-derived score of 0.594. The test set results confirmed the better discriminative capacity of the ML model, with a score of 0.612, surpassing the Cox-derived score C-index of 0.587 (Table 3). These findings were corroborated by the AIC scores, with the ML model showing lower values than the Cox-derived score in the test set, indicating a better overall model fit (Table 3).

In refining our analysis, we segmented the ML score into 4 equal groups within the training set and applied the same classification thresholds to the test set (Figure 1C-D). The ML model maintained similar C-indices after segmentation, indicating that the model's prognostic accuracy is resilient to simplification (Table 3). A substantial reassignment of patients from the intermediate-2 risk group of the Cox score to other risk groups by the ML model was noted (Figure 2A-B). The time-dependent ROC AUCs comparing both models are presented in supplemental Figure 4.

We compared the performance of the ML model with the CIBMTR scoring system⁸ in a comparable subset of patients with complete annotations for the CIBMTR score. By integrating patient age, hemoglobin level at transplant, and donor type, this model defined 3 risk categories in the original series, with a 3-year posttransplant OS of 69%, 51%, and 34% for low, intermediate, and high-risk groups, respectively.

The ML model achieved better performance, with a C-index of 0.608 vs 0.557 in the training set (n = 1925) and 0.654 vs 0.581 in the test set (n = 618; Table 3; supplemental Figure 5). Additionally, the lower AIC scores observed with the ML approach further validated these findings (Table 3). The time-dependent ROC AUCs comparing both models are presented in supplemental Figure 6. The difference between the CIBMTR score and ML method was mostly driven by prognostic refinement within the CIBMTR intermediate risk group, for which the ML algorithm reclassified most patients into different risk categories (Figure 2C-D).

Notably, the ML model exhibited a consistently better discriminative performance than the Cox-derived models for patients with PMF and patients with SMF, with the advantage being more pronounced in the test set (supplemental Table 4).

In predicting NRM, the ML model achieved comparable AIC scores to those of the Cox-based model in the training set but substantially lower AIC scores in the test set, indicating better overall performance (Table 3). Furthermore, when compared with the CIBMTR score, the ML model demonstrated an even more pronounced improvement in overall model fit (Table 3; Figure 3).

The clinical utility of the ML model was evident in its ability to stratify patients into risk groups. Notably, it assigned 25% of patients to the high-risk group, significantly more than the 10.1% in the Cox-derived score and 8.2% in the CIBMTR model (Figure 2). Moreover, the ML model not only identified a larger proportion of high-risk patients but also showed consistent and generalizable results across training and test sets (Figure 1). In the training set, the 12- and 24-month OS rates for the ML high-risk group were 58.9% and 51.5%, respectively, closely aligning with the Cox-derived scores of 58.3% and 52.7%, respectively. In the test set, the ML high-risk group had OS rates of 61.0% at 12 months and 48.1% at 24 months, closely matching the 61.8% and 50.1% of the Cox model high-risk group.

The ML model also identified a larger high-risk population for NRM compared with the Cox-derived score (Figure 3). In the training set, the ML high-risk group had 12- and 24-month NRM rates of 34.9% and 40.7%, respectively, lower than the 46.3% and 48.8% observed in the Cox model high-risk group. However, in the test set, the ML high-risk group showed 12- and 24-month NRM rates of 36.4% and 42.6%, respectively, nearly matching the Cox score rates of 36.0% and 42.8%.

The comparison of patient distribution between the Cox-derived method and the CIBMTR model is elicited in supplemental Figure 7.

To predict OS after allo-HCT in MF, we developed a prognostic tool based on the RSF model, accessible as an interactive web application (https://gemfin.click/ebmt). Figure 4 illustrates the web-based calculator, showing the risk score for a hypothetical transplant candidate.

We compared OS after allo-HCT in patients who received the “optimal” donor type, conditioning intensity, or GVHD prophylaxis with those who did not. Optimal strategies were defined as those that maximize the survival probability according to the model’s predictions for a given patient, conditional on the other individual and disease characteristics. To ensure reliability, these predictions were evaluated exclusively on the test set.

For the optimal donor type, the univariable Cox proportional hazards model indicated that receiving transplantation from a donor type predicted as optimal by the ML model was associated with a HR of 0.76 (95% CI, 0.64-0.89; P = .001) compared with nonoptimal donor type, suggesting a statistically significant survival benefit. However, after adjusting for potential confounding factors in the multivariable analysis using inverse probability weighting (IPW),²⁹ the HR was 0.96 (95% CI, 0.56-1.65; P = .89), indicating no significant survival advantage for ML-predicted optimal donor type.

The univariable Cox model revealed that patients receiving the conditioning regimen intensity predicted as optimal by the ML model had a HR of 0.89 (95% CI, 0.76-1.06; P = .19) compared with those receiving nonoptimal intensities. However, the IPW-adjusted analysis showed a HR of 1.025 (95% CI, 0.07-14.80; P = .99), indicating the lack of predictive value of the ML model to select the optimal conditioning regimen intensity.

Regarding GVHD prophylaxis, neither the unadjusted Cox model nor the IPW-adjusted analysis showed any significant difference in survival between patients receiving the ML-predicted optimal GVHD prophylaxis and those who did not. The unadjusted Cox model yielded a HR of 0.95 (95% CI, 0.79-1.15; P = .62), whereas the IPW-adjusted analysis resulted in a HR of 0.95 (95% CI, 0.43-2.09; P = .90). These results indicate no discernible impact of the ML-predicted optimal GVHD prophylaxis on patient survival.

In this study, we have developed a ML model to enhance risk stratification for patients with MF undergoing allo-HCT, using a large database of 5183 patients with MF from the EBMT registry. Notably, this model is particularly comprehensive, because it considers the broad spectrum of current transplant practices, including diverse conditioning regimens, GVHD prophylaxis approaches, and donor types, such as haploidentical transplants. After dimensionality reduction, the model was simplified to a set of 10 key variables maintaining a notable discriminative capacity for both OS and NRM in both training and test sets.

Comparative analyses demonstrated a better performance of the ML model over a risk score developed within the same cohort using Cox regression methods. It also showed better discriminative capacity than the CIBMTR score,⁸ along with improved generalizability and an enhanced ability to identify a larger group of patients at high risk for posttransplant mortality. Notably, the improved performance of the ML model was evident in both patients with PMF and those with SMF. The ML model’s discriminative capacity remained higher after dividing patients into equally sized risk groups based on individual risk predictions, making the method comparable with traditional risk grouping strategies used in clinical practice. However, the intermediate risk categories identified by the ML model had similar OS and should be consolidated into a single, broader intermediate category, because the model lacks sufficient discrimination within this range. Although the C-indices of the ML model may be deemed moderate in terms of discriminative capacity, our data support its integration into clinical prognostics, offering a more refined and nuanced approach to managing the complexities of patient risk assessment before allo-HCT.

The clinical relevance of our ML model is evident in its ability to stratify allo-HCT candidates into well-defined risk categories. Notably, it identifies 25% of the cohort as high-risk with poor outcome after allo-HCT (∼35% NRM rate and 40% overall mortality at 1 year). Moreover, the web-based calculator permits the identification of a subset of very high-risk patients with a predicted 1-year OS of <50%, allowing for more tailored therapeutic interventions where needed. Accurate risk stratification is essential for optimizing allo-HCT outcomes, enabling physicians to select candidates who are most likely to benefit from transplantation, thereby improving treatment efficacy and patient survival.³ The ML model’s robust capability to classify these patients enhances its utility in clinical decision-making, ultimately fostering more personalized and effective patient management strategies in the complex landscape of allo-HCT.³⁰ 

Although we have proved the feasibility of modeling risk using certain prognostic variables in the context of allo-HCT, the observational nature of the data caution against using these variables to guide the optimal transplantation procedure. The complexity of allo-HCT, characterized by numerous interacting factors that were not all included in our models, limits the practical utility of predictive models for determining the most effective treatment approaches.³¹ This underscores the need for continued research, randomized clinical studies, as well as sophisticated modeling techniques that can account for the dynamic and multifactorial nature of such medical interventions before they can be reliably implemented in clinical decision-making processes.

Our study has several limitations that warrant consideration. Firstly, allo-HCT is a multifaceted procedure influenced by numerous interrelated and independent variables, including differences in patient and disease characteristics and center-specific protocols. This complexity can significantly constrain the power of prognostic models, especially for detecting early posttransplant mortality, which may be influenced by acute and unforeseeable clinical events. We recognize that further investigation into potential, less obvious correlations among pretransplant risk factors could provide additional insights. However, our primary focus was on constructing a clinically actionable prognostic tool rather than conducting an in-depth mechanistic analysis of risk determinants. Additionally, the EBMT data set had a substantial rate of missing data for some variables. The ML method addressed this through data imputation, whereas our Cox model used the missing indicator method for variables with missing values, excluding those with a high degree of missingness (eg, hematologic parameters and spleen size) from score development. Although cross-validation was used to reduce overfitting in the ML model, it was not applied to the risk score developed using Cox regression. The lack of molecular annotation regarding additional somatic mutations, which has been shown to provide prognostic information after allo-HCT in some studies^7,32 but not others,^33,34 prevented a comparison with the Myelofibrosis Transplant Scoring System.⁷ Furthermore, data on the grade of bone marrow fibrosis and the variant allele frequency of driver mutations at transplant were not available. Future research could benefit from enhancing data completeness to potentially refine the model’s prognostic accuracy.

In conclusion, this investigation compared the effectiveness of 2 different Cox regression–derived models with a ML-driven approach for stratifying risk in patients with MF undergoing allo-HCT. The results demonstrate that the ML-driven model outperformed traditional statistical approaches by providing enhanced generalizability and identifying a broader subset of patients at high risk for adverse outcomes. ML methods facilitate the modeling of complex interactions and nonlinear associations more effectively than traditional statistical methods. However, our findings also underscore the challenges in predicting early posttransplant mortality based on conventional baseline characteristics, which remain difficult to anticipate with current prognostic tools. To improve clinical decision-making, we have developed a novel prognostic tool using ML techniques that can identify 25% of patients at high risk for mortality after transplantation. The web-based calculator (https://gemfin.click/ebmt) represents a significant advance toward personalized medicine for patients with MF, enabling better strategic planning and potentially improving outcomes. As we move forward, refining this tool through the integration of more comprehensive data and ongoing validation will be crucial to fully realize its clinical potential.

The authors are grateful to all the centers and patients contributing to the European Society for Blood and Marrow Transplantation database.

Contribution: J.C.H.-B. and A.M.-O. conceived the idea and developed the project proposal; A.M.-O. performed the machine learning analysis, created figures and tables, and cowrote the first draft of the manuscript with J.C.H.-B.; L.G. and J.C.H.-B. developed the Cox regression–based prognostic model for survival; L.K. and J.T. managed the study data; J.R. contributed to the tables and cowrote part of the first draft; C.P.M. and D.C. designed the online calculator; and all other coauthors contributed data to the study, critically revised the paper, and approved the final version.

Conflict-of-interest disclosure: The authors declare no competing financial interests.

A complete list of the centers from the Chronic Malignancies Working Party of the EBMT that participated in this study appears in the supplemental Appendix.

Correspondence: Juan Carlos Hernández-Boluda, Hematology Department, Hospital Clínico Universitario, Avd Blasco Ibáñez 17, 46010 Valencia, Spain; email: hernandez_jca@gva.es; and Donal P. McLornan, Department of Haematology and Stem Cell Transplantation, University College London Hospitals NHS Trust, 3rd Floor W Wing, 250 Euston Rd, London NW1 2PG, United Kingdom; email: donal.mclornan@nhs.net.