Machine learning to predict early recurrence after oesophageal cancer surgery

Background Early cancer recurrence after oesophagectomy is a common problem, with an incidence of 20–30 per cent despite the widespread use of neoadjuvant treatment. Quantification of this risk is difficult and existing models perform poorly. This study aimed to develop a predictive model for early recurrence after surgery for oesophageal adenocarcinoma using a large multinational cohort and machine learning approaches. Methods Consecutive patients who underwent oesophagectomy for adenocarcinoma and had neoadjuvant treatment in one Dutch and six UK oesophagogastric units were analysed. Using clinical characteristics and postoperative histopathology, models were generated using elastic net regression (ELR) and the machine learning methods random forest (RF) and extreme gradient boosting (XGB). Finally, a combined (ensemble) model of these was generated. The relative importance of factors to outcome was calculated as a percentage contribution to the model. Results A total of 812 patients were included. The recurrence rate at less than 1 year was 29·1 per cent. All of the models demonstrated good discrimination. Internally validated areas under the receiver operating characteristic (ROC) curve (AUCs) were similar, with the ensemble model performing best (AUC 0·791 for ELR, 0·801 for RF, 0·804 for XGB, 0·805 for ensemble). Performance was similar when internal–external validation was used (validation across sites, AUC 0·804 for ensemble). In the final model, the most important variables were number of positive lymph nodes (25·7 per cent) and lymphovascular invasion (16·9 per cent). Conclusion The model derived using machine learning approaches and an international data set provided excellent performance in quantifying the risk of early recurrence after surgery, and will be useful in prognostication for clinicians and patients.


Introduction
Oesophageal adenocarcinoma carries a poor prognosis. Among the less than 40 per cent of patients who are candidates for curative treatment 1 , the 5-year survival rate remains approximately 25-50 per cent in randomized trials 2 -4 and rarely exceeds 50 per cent in case series.
Early recurrence (less than 1 year) after surgery is a feared outcome, with rates of 20-30 per cent frequently reported 3 -5 , despite the increasing uptake of neoadjuvant chemotherapy (NACT) and neoadjuvant chemoradiotherapy (NACRT). This is of particular concern because recovery from oesophagectomy is often long and the risk of major complications (Clavien-Dindo grade III-V) is as high as 30 per cent 6 . Many patients have not recovered from the primary cancer treatment when they experience recurrence.
In an ideal setting, prediction of early recurrence before embarking on a multimodal surgical pathway would provide useful information for patients and clinicians. However, staging information correlates poorly between preoperative and postoperative settings 7 , and genomic information is not yet able to predict outcome. Even the most robust preoperative models for prediction have a modest performance at best 8 . In contrast, postoperative information, although not able to influence surgical treatment decisions, is more prognostic and potentially informative for patients. It may also be helpful in making decisions on the merits of adjuvant therapy, further refining the high-risk group of patients in whom novel adjuvant treatments are currently being considered.
Naive logistic regression has been the dominant approach to binary outcome prediction in clinical medicine for decades. Adoption of modern modified regression and machine learning techniques has been limited, in part owing to concerns over computational complexity and reliability. However, an increasing body of evidence has demonstrated that they outperform traditional techniques in predictive performance 9,10 , although this is debatable 11 . In part, the appeal of these approaches lies in their ability to model complex non-linear relationships that are common in cancer data, and which are challenging to model effectively with logistic/linear approaches. The increasing accessibility of software design now also allows the relatively straightforward deployment of these black-box techniques.
The Oesophageal Cancer Clinical and Molecular Stratification (OCCAMS) Consortium 12 previously published a multicentre UK cohort study that assessed survival according to Mandard Tumour Regression Grade (TRG) 13 . This study included patients who had undergone oesophagectomy for adenocarcinoma of the oesophagus or gastro-oesophageal junction (GOJ) preceded by NACT. A clinically meaningful response to NACT was limited to TRG 1-2 only, which represented approximately 15 per cent of patients. The present study used this database, supplemented with an international cohort from the Netherlands, and machine learning techniques to develop and validate a clinically useful predictive model for early recurrence in oesophageal adenocarcinoma.

Methods
The OCCAMS Consortium is a UK-wide multicentre consortium set up to facilitate clinical and molecular stratification of oesophagogastric cancer. It has ethical approval for biological sample collection and analysis in conjunction with detailed clinical annotation (Research Ethics Committee number 10/H0305/1). Data collection and participation in research were approved by institutional ethics committees at each OCCAMS site and University Medical Centre (UMC) Utrecht.

Source of data
Data were sourced from six tertiary oesophagogastric centres in the UK, as described previously 12 . Briefly, the records of consecutive patients from each centre who underwent a planned curative oesophagectomy for adenocarcinoma between 2000 and 2013, and also received NACT (platinum-based triplet or cisplatin and 5-fluorouracil) were reviewed and collated. Treatment was decided by a multidisciplinary team at individual institutions. Neoadjuvant treatment was considered for patients with locally advanced (cT2+) or node-positive disease according to local and national guidelines. Clinical, pathological, recurrence and survival data were recorded. Data from one of the original centres were incomplete to the extent that modelling could not take place and were excluded a priori. To include NACRT as a factor in the model, further patients were identified from University Hospitals Southampton and UMC Utrecht, where CROSS (Chemoradiation for Oesophageal Cancer Followed by Surgery Study)-type NACRT 4 has been the standard of care for oesophageal adenocarcinoma for a number of years. Patients whose tumours were deemed unresectable at the time of surgery or who had metastatic disease on postoperative histology (pM1) were excluded from the analysis.
The primary outcome measure was early recurrence, defined as confirmed local, regional or distant recurrence at less than 1 year from the date of surgery 5,8,14 . Missing data were treated as being missing completely at random and handled by listwise deletion. Modelling was based on a complete-case analysis.

Predictor characteristics
Univariable statistics were calculated using non-parametric Mann-Whitney U and χ 2 tests. The predictive models were generated on the whole data set. All available variables were included in the analysis. A circumferential resection margin (CRM) of less than 1 mm was considered to be involved (and hence R1), in accordance with Royal College of Pathologists (RCP) guidelines 15 . Tumour grade and TRG 13 were assessed by dedicated gastrointestinal histopathologists who were blinded to the clinical data. TRG was used to distinguish between responders (TRG 1-2) and non-responders (TRG [3][4][5], in line with the previous publication based on this data set 12 . To increase the yield of information from lymph node data, both the number of positive lymph nodes and total lymph node harvest were considered as absolute numbers. For the regression model, linearity was assumed for continuous variables. The variables used to predict outcome were: age, sex, tumour location, type of neoadjuvant therapy, response to neoadjuvant therapy (TRG), ypT category, lymphovascular invasion, completeness of resection, grade of differentiation, number of positive lymph nodes and total number of lymph nodes examined.

Model building and validation
Elastic net regularized logistic regression (ELR) 16 was used along with two machine learning techniques: random forest (RF) 17 and extreme gradient boosting (XG boost, XGB) 18 . ELR applies a combination of the ridge and lasso penalties 19,20 with the benefits of both (partly minimization of overfitting and variable selection). RF combines a specified number of decision trees (typically around 1000) created on random subsets of the data set, and is probably the most widely used machine learning approach in the medical literature. XGB attempts to improve sequentially by generating models to explain where the original model fails and then repeating this process (typically around 1000 times), while simultaneously applying regularization to minimize overfitting. Having generated individual models, these were combined to generate overall predictions 21 , an approach that theoretically is particularly beneficial when using diverse model types (such as those described above) that capture different elements of patients' risk profiles.
For ELR, the optimal α and λ hyperparameters (penalty severities) were selected by grid search using tenfold cross-validation with five repeats during model generation and log loss as the metric for optimization. The RF model was derived from 1000 decision trees and hyperparameter tuning was conducted in a similar fashion (for number of variables per tree, split rule and minimum node size). The XGB model was again derived by cross-validation of hyperparameters (number of optimization rounds, maximum tree depth, minimum weight in each child node, minimum loss reduction (γ), regularization penalty (η) and subsampling for regularization). Full details of hyperparameter tuning are available in Appendix S2 (supporting information). These three models were then combined to generate the final (ensemble) model by generating a linear blend of predicted probabilities using logistic regression.
Discrimination of the models was assessed using the area under the receiver operator characteristic (ROC) curve (AUC). In the context of this paper, if two random patients were selected, one with recurrence of cancer at less than 1 year and one disease-free at 1 year, the AUC is equivalent to the probability the model will score the patient with recurrence higher than the patient without. Internal validation was performed using 0⋅632 bootstrapping, with 1000 resampled data sets. Bootstrapping was preferred for internal validation over splitting the cohort into derivation and validation sets, as this has been shown to reduce bias and improve overall model performance, particularly with moderately sized data sets 22 -24 . Calibration was assessed visually and formally with the Hosmer-Lemeshow test. As the data set contains multiple centres with small numbers of patients, an internal-external validation procedure was opted for, as advocated by Steyerberg and Harrell 25 . This entails generating models on all centres apart from one and validating the model on the remaining centre. This process is then repeated leaving each centre out sequentially, and a mean calculated. This method demonstrates how the model performs in external data while also allowing the whole data set to be used for training.
Unadjusted tree models (such as RF, which is included in the ensemble model) and other maximum margin methods typically calibrate poorly as a consequence of their methodology, with predicted probabilities biased towards the centre. To allow meaningful interpretation of probability, isotonic regression was used to scale probabilities on the final model, as described previously 26,27 .
In contrast to logistic regression, assessing global variable importance is challenging using machine learning techniques and to an extent they are black boxes. As coefficients, as seen in logistic regression, are not used, an alternative method is required. The VarImp function of the caret R package was used, where ROC curves are generated for the outcome for each individual predictor, and the contribution to the global ROC curve calculated as a percentage. Owing to the nature of higher-order interactions present in the model, variable importance in individual predictions must be calculated independently. The mean marginal contribution of each variable was calculated (change from the mean prediction; Shapley value 28 ) for individual predictions. A similar approach was used by Nanayakkara and colleagues 29 for analysing in-hospital mortality following cardiac arrest.
Data analysis was conducted using R version 3.5.3 (R Foundation for Statistical Computing, Vienna, Austria).   Values in parentheses are 95 per cent confidence intervals.  a Unscaled calibration (intercept 0⋅395, slope 1⋅574) and b scaled calibration (intercept 0⋅143, slope 0⋅988). The shaded area represents two standard errors.

Results
A total of 812 patients from seven centres were included in model training (Fig. 1). Median age was 64 years and most patients were men (84⋅6 per cent). The majority of tumours were at the GOJ (55⋅5 per cent), and there were high proportions of locally advanced tumours (66⋅7 per cent ypT3-4) and node-positive disease (61⋅0 per cent). First recurrence of cancer within 1 year of surgery was identified in 236 patients (29⋅1 per cent). Patients in the early recurrence group were significantly less likely to have responded to neoadjuvant treatment (8⋅5 versus 21⋅7 per cent), had worse ypT and ypN categories, R1 resection rate and grade of differentiation, and were more likely to have lymphovascular invasion (all P < 0⋅001) ( Table 1). Clinicopathological data are summarized by centre and type of adjuvant therapy in Tables S1 and S2 respectively (supporting information).

Model performance: calibration
Calibration on the training set was visually best in the ELR, and worst in the RF and ensemble models (Fig. S1, supporting information). This was corroborated by the Hosmer Lemeshow test (P = 0⋅806 for ELR, P < 0⋅001 for RF, P = 0⋅030 for XGB, P < 0⋅001 for ensemble). Probabilities generated by the final model were scaled using isotonic regression. Calibration before and after scaling is shown in Fig. 3. A calibration table can be found in Table S3 (supporting information). The Hosmer-Lemeshow test gave a χ 2 value of 38⋅0 (P < 0⋅001) before and 4⋅5 (P = 0⋅806) after scaling. Similarly, the Brier score, a measure of overall model performance, also improved from 0⋅119 to 0⋅114.

Variable importance
Coefficients and odds ratios cannot be generated for these models. Therefore, variable importance as a percentage contribution to the model was computed (Table 3). Overall, the most influential predictor variable was number of positive lymph nodes (25⋅7 per cent), followed by lymphovascular invasion (16⋅9 per cent). There was considerable variability in importance across models. For example, age contributed 0⋅3 per cent to the ELR model, 18⋅2 per cent to the RF model, 10⋅2 per cent to the XGB model and 9⋅6 per cent to the final model. The percentage contribution of each variable in each example patient is shown. This is represented as an absolute percentage change from the mean predicted value of 27⋅4 per cent. A calculator for this is packaged with the online model. LN, lymph node; NACT, neoadjuvant chemotherapy; NACRT, neoadjuvant chemoradiotherapy.
It is important to restate that relationships between the variables and outcome are non-linear and their importance varies considerably according to other variables owing to higher-order interactions. As an example, even though lymph node status was found to be the most influential marker overall, combinations of other variables would make other variables most important in individual patients. To illustrate this and demonstrate how variables interact, three example patients were considered (Tables 4  and 5). The technique used measures the change in prediction from the mean prediction (27⋅4 per cent) that can be attributed to each predictor variable. This approach (calculation of Shapley value) originates from cooperative game theory.

Discussion
An easy-to-use and robust clinical model for predicting the risk of early recurrence after surgery for oesophageal adenocarcinoma was derived in this study. It uses routinely collected clinical and pathological data that should be available for every patient; together, these allow considerably more precision in risk estimation than would be possible using individual variables that are known to be influential, such as pathological lymph node involvement. The final model demonstrated excellent discrimination, and validation techniques supported the generalizability of the approach.
In addition to prognostication, this model may be useful for planning adjuvant therapy. Early recurrence after oesophagectomy, often before recovery from surgery is complete, is a devastating outcome for patients. Targeting existing and emerging treatment combinations in this patient group to prolong time to recurrence or prevent recurrence is vital, but can only happen with accurate predictions of the likelihood of relapse. The starting point for consideration of treatment escalation or novel combinations (such as immunotherapy) after surgery is the identification of patients who are at high risk of recurrence. The authors have purposefully avoided dichotomization/stratification based on outcome, and presented raw probability in preference to this. This will allow full discussions between surgeons/oncologists and patients regarding the benefits of adjuvant therapy and tailored to the individual patient's postoperative recovery and wishes. It may also allow stratification of adjuvant trials based on layered levels of risk.
This cohort exhibited an early recurrence rate of 29⋅1 per cent, which is similar to that in previous reports 3 -5,8 where this outcome was specified explicitly. There was also an R1 resection rate of 28⋅4 per cent, in line with previously reported data 34,35 based on an RCP definition of CRM positivity (CRM less than 1 mm is positive). In univariable analysis, all factors expected to correlate with worse prognosis (including ypT, ypN, lymphovascular invasion, R1 resection and grade of differentiation) were significantly worse in patients who developed early recurrence. This validates the present cohort as a true representation of contemporary practice and a sensible place to begin building more complex models.
Discrimination of the different models was similar, with minimal variability in AUC values between models on validation. However, the ensemble model consistently performed the best and is a suitable choice for the final model. The decline in performance from the training set to validation, which was particularly marked for the RF and ensemble models, is a consequence of the tuning process, whereby the optimum values are chosen from a grid of thousands after repeated tests (in this case repeated 10-fold cross-validation). In this setting, the apparent performance of the model on the training set is overestimated and should be disregarded.
There was marked heterogeneity in variable importance between models. This is interesting, particularly in the context of the models performing so similarly overall, and supports the idea of combining them to capture different patient information. The most important variables overall were number of positive lymph nodes and lymphovascular invasion, which accounted for 42⋅6 per cent of performance in the final model. This is not only biologically sensible, but the subject of several recent publications 12,36,37 and ongoing translational work. Although not available for this study, more detail regarding lymphadenopathy, such as downstaging and anatomical location, would probably be informative. It is difficult to reach firm conclusions regarding variables considering the nature of the study. However, the authors draw attention to two facets of the model. First, TRG was the least influential variable across the board, with an importance of almost 0 per cent. This suggests that in itself TRG adds no information over the other measured variables in predicting early outcomes. This is in keeping with emerging data regarding the genomic disparity between primary tumours and their metastasis (lymph node or distant) 38 , and a previous report 12 of the importance of lymph node downstaging to clinical outcome. Second, type of treatment was the third most important determinant of outcome, with NACT having an advantage over NACRT. In this cohort, although the postoperative pathology was considerably more favourable after NACRT, the rate of early recurrence was no less, and tended to be higher (NACRT 35⋅5 per cent, NACT 27⋅5 per cent; P = 0⋅061 (Table S2, supporting information)). This suggests that, although postoperative pathology is more favourable with NACRT, this does not translate to better outcome 39 -41 ; hence ypT3 N1 R0 status after NACT does not have the same meaning as a ypT3 N1 R0 result after NACRT, at least in the early phase after treatment. This is important in postoperative discussions with patients. As the machine learning approaches detailed here allow interactions between variables, the model suggests that NACRT confers a greater risk; however, this increased risk is conditional on the other variables being static rather than an overall increase in risk from having NACRT.
To explore this further, details of recurrence location (locoregional versus distant) would be informative. However, owing to the historical nature of the data for the majority of the patients (collected for the first study) it was not possible to ascertain this reliably for most of the cohort.
The concern with NACRT is that improved locoregional control is at the expense of undertreatment of microscopic distant disease, particularly where the radiotherapy field is limited anatomically (for example with GOJ tumours). The expected consequence of this would be fewer locoregional recurrences and more distant recurrences, although this has not been demonstrated in other comparative studies and a recently published RCT 41 .
The present study lacks the number of patients needed to separately analyse the influence of neoadjuvant treatment on oesophageal and GOJ tumours, however, the individual variable importance calculation available in the web app allows some insight to be gained. Here, the relative negative influence of NACRT (increased risk of recurrence compared with NACT) is, on the whole, more pronounced for GOJ tumours compared with oesophageal tumours (an example of a second-order interaction), despite the recurrence rate being higher for oesophageal than GOJ tumours.
Other risk factors for early recurrence, including perioperative blood transfusion 42 , complications of surgery 43 and preoperative staging, were not available for this study, but are less discriminatory. Nor were precise neoadjuvant regimens available for all patients. It is therefore unclear whether these results would be influenced by completion of treatment as prescribed, or indeed by whether any adjuvant therapy was given. The absence of these factors seems to have minimal effect on the model, suggesting a small margin of effect on outcomes. Combining these factors could potentially increase the performance of the present model if incorporated in the future. Ultimately, differential gene expression and mutation 44,45 may well determine prognostication and treatment pathways 46 , but such data are unlikely to be available universally for some years. Until then, clinical and histopathological data remain the standard.
In that context, gains from mathematical and computer-based techniques are key to precision in delivery of cancer care. Here, several modern approaches that produce viable models were demonstrated. This study used a data set that was relatively small and simple in a machine learning context, and the improvement in performance over a standard logistic regression was small (internal validation AUC 0⋅781). This is nonetheless important as such an improvement is in effect 'free'. The strengths of this study lie in its multicentre nature and the heterogeneity of the cohort. This approach should maximize the utility of the model on external populations. All the data points used should be collected routinely at the majority of institutions, which should allow uptake without change in practice. The College of American Pathologists (CAP) definition of CRM positivity (CRM positive if there is tumour at the resection margin) was derivable for centre G, and performance was preserved in this subgroup if that definition was used instead of the RCP definition (AUC 0⋅813 with model generated on centres A-F (650 patients) and validated on centre G (162)) (Fig. S2, supporting information), supporting utility in both settings. The study focused on predictive model study design and reporting as suggested by the AJCC 47 and TRIPOD 48 statements.
The training set was limited to patients undergoing neoadjuvant therapy for adenocarcinoma of the oesophagus. No attempt was made to apply the model to a chemotherapy-naive population, and it is unlikely to calibrate well in this group owing to the differing influence on survival of yp compared with p staging 49 . It is also unclear whether the model would be valid in patients with squamous cell carcinoma; the authors advocate an early external validation exercise using this patient group. A formal prospective validation/recalibration using the CAP definition of CRM positivity would also be beneficial. Simulation studies have suggested that 100-200 cases (positives) are required for accurate validation 50 , which, assuming a stable incidence, would require approximately 380-760 patients. A further limitation was the significant proportion of the original patients with missing data, which will have introduced a degree of selection bias. Multiple imputation is possible as a means of addressing this, but was considered less appropriate in this study because of the high proportion of missing data for the outcome measure and the lack of an external validation set.
A large multicentre cohort of patients who underwent oesophagectomy has been used to derive an accurate prediction model for early cancer recurrence, with excellent performance on validation. Machine learning techniques represent an attractive proposition for maximizing performance of predictive models.