Intended for healthcare professionals

CCBY Open access

Development and internal-external validation of statistical and machine learning models for breast cancer prognostication: cohort study

BMJ 2023; 381 doi: (Published 10 May 2023) Cite this as: BMJ 2023;381:e073800
  1. Ash Kieran Clift, clinical research fellow1 2,
  2. David Dodwell, senior clinical research fellow3,
  3. Simon Lord, senior clinical researcher4,
  4. Stavros Petrou, professor of health economics2,
  5. Michael Brady, professor of oncological imaging4,
  6. Gary S Collins, professor of medical statistics5,
  7. Julia Hippisley-Cox, professor of clinical epidemiology and general practice2
  1. 1Cancer Research UK Oxford Centre, Oxford, UK
  2. 2Nuffield Department of Primary Care Health Sciences, Radcliffe Primary Care Building, Radcliffe Observatory Quarter, University of Oxford, Oxford OX2 6GG, UK
  3. 3Nuffield Department of Population Health, University of Oxford, Oxford, UK
  4. 4Department of Oncology, University of Oxford, Oxford, UK
  5. 5Centre for Statistics in Medicine, Nuffield Department of Orthopaedics, Rheumatology and Musculoskeletal Sciences, University of Oxford, Oxford, UK
  1. Correspondence to: A K Clift ashley.clift{at} (or @AshClift on Twitter)
  • Accepted 28 March 2023


Objective To develop a clinically useful model that estimates the 10 year risk of breast cancer related mortality in women (self-reported female sex) with breast cancer of any stage, comparing results from regression and machine learning approaches.

Design Population based cohort study.

Setting QResearch primary care database in England, with individual level linkage to the national cancer registry, Hospital Episodes Statistics, and national mortality registers.

Participants 141 765 women aged 20 years and older with a diagnosis of invasive breast cancer between 1 January 2000 and 31 December 2020.

Main outcome measures Four model building strategies comprising two regression (Cox proportional hazards and competing risks regression) and two machine learning (XGBoost and an artificial neural network) approaches. Internal-external cross validation was used for model evaluation. Random effects meta-analysis that pooled estimates of discrimination and calibration metrics, calibration plots, and decision curve analysis were used to assess model performance, transportability, and clinical utility.

Results During a median 4.16 years (interquartile range 1.76-8.26) of follow-up, 21 688 breast cancer related deaths and 11 454 deaths from other causes occurred. Restricting to 10 years maximum follow-up from breast cancer diagnosis, 20 367 breast cancer related deaths occurred during a total of 688 564.81 person years. The crude breast cancer mortality rate was 295.79 per 10 000 person years (95% confidence interval 291.75 to 299.88). Predictors varied for each regression model, but both Cox and competing risks models included age at diagnosis, body mass index, smoking status, route to diagnosis, hormone receptor status, cancer stage, and grade of breast cancer. The Cox model’s random effects meta-analysis pooled estimate for Harrell’s C index was the highest of any model at 0.858 (95% confidence interval 0.853 to 0.864, and 95% prediction interval 0.843 to 0.873). It appeared acceptably calibrated on calibration plots. The competing risks regression model had good discrimination: pooled Harrell’s C index 0.849 (0.839 to 0.859, and 0.821 to 0.876, and evidence of systematic miscalibration on summary metrics was lacking. The machine learning models had acceptable discrimination overall (Harrell’s C index: XGBoost 0.821 (0.813 to 0.828, and 0.805 to 0.837); neural network 0.847 (0.835 to 0.858, and 0.816 to 0.878)), but had more complex patterns of miscalibration and more variable regional and stage specific performance. Decision curve analysis suggested that the Cox and competing risks regression models tested may have higher clinical utility than the two machine learning approaches.

Conclusion In women with breast cancer of any stage, using the predictors available in this dataset, regression based methods had better and more consistent performance compared with machine learning approaches and may be worthy of further evaluation for potential clinical use, such as for stratified follow-up.


Clinical prediction models already support medical decision making in breast cancer by providing individualised estimations of risk. Tools such as PREDICT Breast1 or the Nottingham Prognostic Index23 are used in patients with early stage, surgically treated breast cancer for prognostication and selection of post-surgical treatment. Such tools are, however, inherently limited to treatment specific subgroups of patients. Accurate estimation of mortality risk after diagnosis across all patients with breast cancer of any stage may be clinically useful for stratifying follow-up, counselling patients about their expected prognosis, or identifying high risk individuals suitable for clinical trials.4

The scope for machine learning approaches in clinical prediction modelling has attracted considerable interest.56789 Some have posited that these flexible approaches might be more suitable for capturing non-linear associations, or for handling higher order interactions without explicit programming.10 Others have raised concerns about model transparency,1112 interpretability,13 risk of algorithmic bias exacerbating extant health inequalities,14 quality of evaluation and reporting,15 ability to handle rare events16 or censoring,17 and appropriateness of comparisons11 to regression based methods.18 Indeed, systematic reviews have shown no inherent benefit of machine learning approaches over appropriate statistical models in low dimensional clinical settings.18 As no a priori method exists to predict which modelling approach may yield the most useful clinical prediction model for a given scenario, frameworks that appropriately compare different models can be used.

Owing to the risks of harm from suboptimal medical decision making, clinical prediction models should be comprehensively evaluated for performance and utility,19 and, if widespread clinical use is intended, heterogeneity in model performance across relevant patient groups should be explored.20 Given developments in treatment for breast cancer over time, with associated temporal falls in mortality, another key consideration is the transportability of risk models—not just across regions and subpopulations but also across time periods.21 Although such dataset shift22 is a common issue with any algorithm sought to be deployed prospectively, this is not routinely explored. Robust evaluation is necessary but is non-uniform in the modelling of breast cancer prognostication.23 A systematic review identified 58 papers that assessed prognostic models for breast cancer,24 but only one study assessed clinical effectiveness by means of a simplistic approach measuring the accuracy of classifying patients into high or low risk groups. A more recent systematic review25 appraised 922 breast cancer prediction models using PROBAST (prediction model risk of bias assessment tool)26 and found that most of the clinical prediction models are poorly reported, show methodological flaws, or are at high risk of bias. Of the 27 models deemed to be at low risk of bias, only one was intended to estimate the risks of breast cancer related mortality in women with disease of any stage.27 However, this small study of 287 women using data from a single health department in Spain had methodological limitations, including possibly insufficient data to fit a model (see supplementary table 1) and uncertain transportability to other settings. Therefore, no reliable prediction model exists to provide accurate risk assessment of mortality in women with breast cancer of any stage. Although we refer to women throughout, this is based on self-reported female sex, which may include some individuals who do not identify as female.

We aimed to develop a clinically useful prediction model to reliably estimate the risks of breast cancer specific mortality in any woman with a diagnosis of breast cancer, in line with modern best practice. Utilising data from 141 675 women with invasive breast cancer diagnosed between 2000 and 2020 in England from a population representative, national linked electronic healthcare record database, this study comparatively developed and evaluated clinical prediction models using a combination of analysis methods within an internal-external validation strategy.2829 We sought to identify and compare the best performing methods for model discrimination, calibration, and clinical utility across all stages of breast cancer.


We evaluated four model building approaches: two regression methods (Cox proportional hazards and competing risks regression) and two machine learning methods (XGBoost and neural networks). The prediction horizon was 10 year risk of breast cancer related death from date of diagnosis. The study was conducted in accordance with our protocol30 and is reported consistent with the TRIPOD (transparent reporting of a multivariable prediction model for individual prognosis or diagnosis) guidelines.31

Sample size calculations

Assuming 100 candidate predictor parameters, an annual mortality rate of 0.024 after diagnosis,32 and a conservative 15% of the maximal Cox-Snell R2, we estimated that the minimum sample size for fitting the regression models was 10 080, with 1452 events, and 14.52 events for each predictor parameter.3334 No standard method exists to estimate minimum sample size for our machine learning models of interest—some evidence, albeit on binary outcome data, suggests that some machine learning methods may require much more data.35

Study population and data sources

The QResearch database was used to identify an open cohort of women aged 20 years and older (no upper age limit) at time of diagnosis of any invasive breast cancer between 1 January 2000 and 31 December 2020 in England. QResearch has collected data from more than 1500 general practices in the United Kingdom since 1989 and comprises individual level linkage across general practice data, NHS Digital’s Hospital Episode Statistics, the national cancer registry, and the Office for National Statistics death registry.

Patient and outcome definitions

The outcome for this study was breast cancer related mortality within 10 years from the date of a diagnosis of invasive breast cancer. We defined the diagnosis of invasive breast cancer as the presence of breast cancer related Read/Systemised Nomenclature of Medicine – Clinical Terms (SNOMED) codes in general practice records, breast cancer related ICD-10 (international classification of diseases, 10th revision) codes in Hospital Episode Statistics data, or as a patient with breast cancer in the cancer registry (stage >0; whichever occurred first). The outcome, breast cancer death, was defined as the presence of relevant ICD-10 codes as any cause of death (primary or contributory) on death certificates from the ONS register. We excluded women with recorded carcinoma in situ only diagnoses as these are non-obligate precursor lesions and present distinct clinical considerations.36 Clinical codes used to define predictors and outcomes are available in the QResearch code group library ( Follow-up time was calculated from the first recorded date of breast cancer diagnosis (earliest recorded on any of the linked datasets) to the earliest of breast cancer related death, other cause of death, or censoring (reached end of study period, left the registered general practice, or the practice stopped contributing to QResearch). The status at last follow-up depended on the modelling framework (ie, Cox proportional hazards or competing risks framework). The maximum follow-up was truncated to 10 years, in line with the model prediction horizon. Supplementary table 2 shows ascertainment of breast cancer diagnoses across the linked datasets.

Candidate predictor parameters

Individual participant data were extracted on the candidate predictor parameters listed in Box 1, as well as geographical region, auxiliary variables (breast cancer treatments), and dates of events of interest. Candidate predictors were based on evidence from the clinical, epidemiological, or prediction model literature.12337383940 The most recently recorded values before or at the time of breast cancer diagnosis were used with no time restriction. Data were available from the cancer registry about cancer treatment within one year of diagnosis (eg, chemotherapy) but without any corresponding date. The intended model implementation (prediction time) would be at the breast cancer multidisciplinary team meeting or similar clinical setting, following initial diagnostic investigations and staging. To avoid information leakage, and since we did not seek model treatment selection within a causal framework,41 breast cancer treatment variables were not included as predictors.

Box 1

Candidate predictor parameters for models

Candidate predictor parameters, definitions, and functional forms explored

  • Age at breast cancer diagnosis—continuous or fractional polynomial

  • Townsend deprivation score at cohort entry—continuous or fractional polynomial

  • Body mass index (most recently recorded before breast cancer diagnosis)—continuous or fractional polynomial

  • Self-reported ethnicity

  • Tumour characteristics:

    • Cancer stage at diagnosis (ordinal: I, II, III, IV)

    • Differentiation (categorical: well differentiated, moderately differentiated, poorly or undifferentiated)

    • Oestrogen receptor status (binary: positive or negative)

    • Progesterone receptor status (binary: positive or negative)

    • Human epidermal growth factor receptor 2 (HER2) status (binary: positive or negative)

    • Route to diagnosis (categorical: emergency presentation, inpatient elective, other, screen detected, two week wait)

  • Comorbidities or medical history on general practice or Hospital Episodes Statistics data (recorded before or at entry to cohort; categorical unless stated otherwise):

    • Hypertension

    • Ischaemic heart disease

    • Type 1 diabetes mellitus

    • Type 2 diabetes mellitus

    • Chronic liver disease or cirrhosis

    • Systemic lupus erythematosus

    • Chronic kidney disease (ordinal: none or stage 2, stage 3, stage 4, stage 5)

    • Vasculitis

  • Family history of breast cancer (categorical: recorded in general practice or Hospital Episodes Statistics data, before or at entry to cohort)

  • Drug use (before breast cancer diagnosis):

    • Hormone replacement therapy

    • Antipsychotic

    • Tricyclic antidepressant

    • Selective serotonin reuptake inhibitor

    • Monoamine oxidase inhibitor

    • Oral contraceptive pill

    • Angiotensin converting enzyme inhibitor

    • β blocker

    • Renin-angiotensin aldosterone antagonists

  • Age (fractional polynomial terms)×family history of breast cancer

  • Ethnicity×age (fractional polynomial terms)

  • Interactions and fractional polynomials were not included in the machine learning models.


Fractional polynomial42 terms for the continuous variables age at diagnosis, Townsend deprivation score, and body mass index (BMI) at diagnosis were identified in the complete data. This was done separately for the Cox and competing risks regression models, with a maximum of two powers permitted.

Missing data

Multiple imputation with chained equations was used to impute missing data for BMI, ethnicity, Townsend deprivation score, smoking status, cancer stage at diagnosis, cancer grade at diagnosis, HER2 status, oestrogen receptor status, and progesterone receptor status under the missing at random assumption.4344 The imputation model contained all other candidate predictors, the endpoint indicator, breast cancer treatment variables, the Nelson-Aalen cumulative hazard estimate,45 and the period of cohort entry (period 1=1 January 2000-31 December 2009; period 2=1 January 2010-31 December 2020). The natural logarithm of BMI was used in imputation for normality, with imputed values exponentiated back to the regular scale for modelling. We generated 50 imputations and used these in all model fitting and evaluation steps. Although missing data were observed in the linked datasets used for model development, in the intended use setting (ie, risk estimation at breast cancer multidisciplinary team after a medical history has been taken), the predictors would be expected to be available for all patients.

Modelling strategy

Models were fit to the entire cohort and then evaluated using internal-external cross validation,28 which involved splitting the dataset by geographical region (n=10) and time period (see figure 1 for summary). For the internal-external cross validation, we recalculated follow-up so that those women who entered the study during the first study decade and survived into the second study period had their follow-up truncated (and status assigned accordingly) at 31 December 2009. This was to emulate two wholly temporally distinct datasets, both with maximum follow-up of 10 years, for the purposes of estimating temporal transportability of the models.

Fig 1
Fig 1

Summary of internal-external cross validation framework used to evaluate model performance for several metrics, and transportability

Cox proportional hazards modelling

For the approach using Cox proportional hazards modelling, we treated other (non-breast cancer) deaths as censored. A full Cox model was fitted using all candidate predictor parameters. Model fitting was performed in each imputed dataset and the results combined using Rubin’s rules, and then this pooled model was used as the basis for predictor selection. We selected binary or multilevel categorical predictors associated with exponentiated coefficients >1.1 or <0.9 (at P<0.01) for inclusion, and interactions and continuous variables were selected if associated with P<0.01. Then these were used to refit the final Cox model. The predictor selection approach benefits from starting with a full, plausible, maximally complex model,46 and then considers both the clinical and the statistical magnitude of predictors to select a parsimonious model while making use of multiply imputed data.4748 This approach has been used in previous clinical prediction modelling studies using QResearch.495051 Clustered standard errors were used to account for clustering of participants within individual general practices in the database.

Competing risks models

Deaths from other, non-breast cancer related causes represent a competing risk and in this framework were handled accordingly.30 We repeated the fractional polynomial term selection and predictor selection processes for the competing risks models owing to potential differential associations between predictors and risk or functional forms thereof. A full model was fit with all candidate predictors, with the same magnitude and significance rule used to select the final predictors.

The competing risks model was developed using jack-knife pseudovalues for the Aalen-Johansen cumulative incidence function at 10 years as the outcome variable52—the pseudovalues were calculated for the overall cohort (for fitting the model) and then separately in the data from period 1 and from period 2 for the purposes of internal-external cross validation. These values are a marginal (pseudo) probability that can then be used in a regression model to predict individuals’ probabilities conditional on the observed predictor values. Pseudovalues for the cumulative incidence function at 10 years were regressed on the predictor parameters in a generalised linear model with a complementary log-log link function525354 and robust standard errors to account for the non-independence of pseudovalues. The resultant coefficients are statistically similar to those of the Fine-Gray model5254 but computationally less burdensome to obtain, and permit direct modelling of probabilities.

All fitting and evaluation of the Cox and competing risks regression models occurred in each separate imputed dataset, with Rubin’s rules used to pool coefficients and standard errors across all imputations.55

XGBoost and neural network models

The XGBoost and neural network approaches were adapted to handle right censored data in the setting of competing risks by using the jack-knife pseudovalues for the cumulative incidence function at 10 years as a continuous outcome variable. The same predictor parameters as selected for the competing risks regression model were used for the purposes of benchmarking. The XGBoost model used untransformed values for continuous predictors, but these were minimum-maximum scaled (constrained between 0 and 1) for the neural network. We converted categorical variables with more than two levels to dummy variables for both machine learning approaches.

We fit the XGBoost and neural network models to the entire available cohort and used bayesian optimisation56 with fivefold cross validation to identify the optimal configuration of hyperparameters to minimise the root mean squared error between observed pseudovalues and model predictions. Fifty iterations of bayesian optimisation were used, with the expected improvement acquisition function.

For the XGBoost model, we used bayesian optimisation to tune the number of boosting rounds, learning rate (eta), tree depth, subsample fraction, regularisation parameters (alpha gamma, and lambda), and column sampling fractions (per tree, per level). We used the squared error regression option as the objective, and the root mean squared error as the evaluation metric.

To permit modelling of higher order interactions in this tabular dataset, we used a feed forward artificial neural network approach with fully connected dense layers: the model architecture comprised an input layer of 26 nodes (ie, number of predictor parameters), rectified linear unit activation functions in each hidden layer, and a single linear activation output node to generate predictions for the pseudovalues of the cumulative incidence function. The Adam optimiser was used,57 with the initial learning rate, number of hidden layers, number of nodes in each hidden layer, and number of training epochs tuned using bayesian optimisation. If the loss function had plateaued for three epochs, we halved the learning rate, with early stopping after five epochs if the loss function had not reduced by 0.0001. The loss function was the root mean squared error between observed and predicted pseudovalues due to the continuous nature of the target variable.58

After identification of the optimal hyperparameter configurations, we fit the models accordingly to the entirety of the cohort data. We then assessed the performance of these models using the internal-external cross validation strategy—this resembled that for the regression models but with the addition of a hyperparameter tuning component (fig 1). During each iteration of internal-external cross validation, we used bayesian optimisation with fivefold cross validation to identify the optimal hyperparameters for the model fitted to the development data from period 1, which we then tested on the held-out period 2 data. This therefore constituted a form of nested cross validation.59

As the XGBoost and neural network models do not constitute a linear set of parameters and do not have standard errors (therefore not able to be pooled using Rubin’s rules), we used a stacked imputation strategy. The 50 imputed datasets were stacked to form a single, long dataset, which enabled us to use the same full data as for the regression models, avoiding suboptimal approaches such as complete case analysis or single imputation. For model evaluation after internal-external cross validation, we used approaches based on Rubin’s rules,55 with performance estimates calculated in each separate imputed dataset using the internal-external cross validation generated individual predictions, and then the estimates were pooled.

Performance evaluation

Predicted risks when using the Cox model can be derived by combining the linear predictor with the baseline hazard function using the equation: predicted event probability=1−Stexp(Xβ) where St is the baseline survival function calculated at 10 years, and Xβ is the individual’s linear predictor. For internal-external cross validation, we estimated baseline survival functions separately in each imputation in the period 1 data (continuous predictors centred at the mean, binary predictors set to zero), with results pooled across imputations in accordance with Rubin’s rules.55 We estimated the final model’s baseline function similarly but using the full cohort data.

Probabilistic predictions for the competing risks regression model were directly calculated using the following transformation of the linear predictors (Xβ, which included a constant term): predicted event probability=1−exp(−exp(Xβ)).

As the XGBoost and neural network approaches modelled the pseudovalues directly, we handled the generated predictions as probabilities (conditional on the predictor values). As pseudovalues are not restricted to lie between 0 and 1, we clipped the XGBoost and neural network model predictions to be between 0 and 1 to represent predicted probabilities for model evaluation.

Discrimination was assessed using Harrell’s C index,60 calculated at 10 years and taking censoring into account—this used inverse probability of censoring weights for competing risks regression, XGBoost, and neural networks given their competing risks formulation.61 Calibration was summarised in terms of the calibration slope and calibration-in-the-large.6263 Region level results for these metrics were computed during internal-external cross validation and pooled using random effects meta-analysis20 with the Hartung-Knapp-Sidik-Jonkmann method64 to provide an estimate of each metric with a 95% confidence interval, and with a 95% prediction interval. The prediction interval estimates the range of model performance on application to a distinct dataset.20 We also computed these metrics by ethnicity, 10 year age groups, and cancer stage (I-IV) using the pooled, individual level predictions.

Using the individual level predictions from all models, we generated smoothed calibration plots to assess alignment of observed and predicted risks across the spectrum of predicted risks. We generated these using a running smoother through individual risk predictions, and observed individual pseudovalues65 for the Kaplan-Meier failure function (Cox model) or cumulative incidence function (all other models).

Meta-regression following Hartung-Knapp-Sidik-Jonkmann random effects models were used to calculate measures of I2 and R2 to assess the extent to which inter-regional heterogeneity in discrimination and calibration metrics could be attributable to regional variation in age, BMI (standard deviation thereof), mean deprivation score, and ethnic diversity (percentage of people of non-white ethnicity).20 These region level characteristics were estimated using the data from period 2.

We compared the models for clinical utility using decision curve analysis.66 This analysis assesses the trade-off between the benefits of true positives (breast cancer deaths) and the potential harms that may arise from false positives across a range of threshold probabilities. Each model was compared using the two default scenarios of treat all or treat none, with the mean model prediction used for each individual across all imputations. This approach implicitly takes into account both discrimination and calibration and also extends model evaluation to consider the ramifications on clinical decision making.67 The competing risk of other, non-breast-cancer death was taken into account. Decision curves were plotted overall, and by cancer stage to explore potential utility for all breast cancers.

Predictions generated from the Cox proportional hazards model and other, competing risks approaches have different interpretations, owing to their differential handling of competing events and their modelling of hazard functions with distinct statistical properties.

Software and code

Data processing, multiple imputation, regression modelling, and evaluation of internal-external cross validation results utilised Stata (version 17). Machine learning modelling was performed in R 4.0.1 (xgboost, keras, and ParBayesianOptimization packages), with an NVIDIA Tesla V100 used for graphical processing unit support. Analysis code is available in repository

Patient and public involvement

Two people who survived breast cancer were involved in discussions about the scope of the project, candidate predictors, importance of research questions, and co-creation of lay summaries before submitting the project for approval. This project was also presented at an Oxfordshire based breast cancer support group to obtain qualitative feedback on the study’s aims and face validity or plausibility of candidate predictors, and to discuss the acceptability of clinical risk models to guide stratified breast cancer care.


Study cohort and incidence rates

A total of 141 765 women aged between 20 and 97 years at date of breast cancer diagnosis were included in the study. During the entirety of follow-up (median 4.16 (interquartile range 1.76-8.26) years), there were 21 688 breast cancer related deaths and 11 454 deaths from other causes. Restricting to 10 years maximum follow-up from breast cancer diagnosis, 20 367 breast cancer related deaths occurred during a total of 688 564.81 person years. The crude mortality rate was 295.79 per 10 000 person years (95% confidence interval 291.75 to 299.88). Supplementary figure 1 presents ethnic group specific mortality curves. Table 1 shows the baseline characteristics of the cohort overall and separately by decade defined subcohort.

Table 1

Summary characteristics of final study cohort overall and separated into temporally distinct subcohorts used in internal-external cross validation. Values are number (column percentage) unless stated otherwise

View this table:

After the cohort was split by decade of cohort entry and follow-up was truncated for the purposes of internal-external cross validation, 7551 breast cancer related deaths occurred in period 1 during a total of 211 006.95 person years of follow-up (crude mortality rate 357.96 per 10 000 person years (95% confidence interval 349.87 to 366.02)). In the period 2 data, 8808 breast cancer related deaths occurred during a total of 297 066.74 person years of follow-up, with a lower crude mortality rate of 296.50 per 10 000 person years (290.37 to 302.76) observed.

Cox proportional hazards model

We selected non-linear fractional polynomial terms for age and BMI (see supplementary figure 2). The final Cox model after predictor selection is presented as exponentiated coefficients in figure 2 for transparency, with the full model detailed in supplementary table 3. Model performance across all ethnic groups is summarised in supplementary table 4: discrimination ranged between a Harrell’s C index of 0.794 (95% confidence interval 0.691 to 0.896) in Bangladeshi women to 0.931 (0.839 to 1.000) in Chinese women, but the low numbers of event counts in smaller ethnic groups (eg, Chinese) meant that overall calibration indices were imprecisely estimated for some.

Fig 2
Fig 2

Final Cox proportional hazards model predicting 10 year risk of breast cancer mortality, presented as its exponentiated coefficients (hazard ratios with 95% confidence intervals). Model contains fractional polynomial terms for age (0.5, 2) and body mass index (2, 2), but these are not plotted owing to reasons of scale. Model also includes a baseline survival term (not plotted—the full model as coefficients is presented in the supplementary file). ACE=angiotensin converting enzyme; CI=confidence interval; CKD=chronic kidney disease; ER=oestrogen receptor; GP=general practitioner; HER2= human epidermal growth factor receptor 2; HRT=hormone replacement therapy; PR=progesterone receptor; RAA=renin-angiotensin aldosterone; SSRI=selective serotonin reuptake inhibitor

Overall, the Cox model’s random effects meta-analysis pooled estimate for Harrell’s C index was the highest of any model, at 0.858 (95% confidence interval 0.853 to 0.864, 95% prediction interval 0.843 to 0.873). A small degree of miscalibration occurred on summary metrics, with a meta-analysis pooled estimate for the calibration slope of 1.108 (95% confidence interval 1.079 to 1.138, 95% prediction interval 1.034 to 1.182) (table 2). Figure 3, figure 4, and figure 5 show the meta-analysis pooling of performance metrics across regions. Smoothed calibration plots showed generally good alignment of observed and predicted risks across the entire spectrum of predicted risks, albeit with some minor over-prediction (fig 6).

Table 2

Summary performance metrics for all four models, estimated using random effects meta-analysis after internal-external cross validation.

View this table:
Fig 3
Fig 3

Results from internal-external cross validation of Cox proportional hazards model for Harrell’s C index. Plots display region level performance metric estimates and 95% confidence intervals (diamonds with lines), and an overall pooled estimate obtained using random effects meta-analysis and 95% confidence interval (lowest diamond) and 95% prediction interval (line through lowest diamond). CI=confidence interval

Fig 4
Fig 4

Results from internal-external cross validation of Cox proportional hazards model for calibration slope. Plots display region level performance metric estimates and 95% confidence intervals (diamonds with lines), and an overall pooled estimate obtained using random effects meta-analysis and 95% confidence interval (lowest diamond) and 95% prediction interval (line through lowest diamond). CI=confidence interval

Fig 5
Fig 5

Results from internal-external cross validation of Cox proportional hazards model for calibration-in-the-large. Plots display region level performance metric estimates and 95% confidence intervals (diamonds with lines), and an overall pooled estimate obtained using random effects meta-analysis and 95% confidence interval (lowest diamond) and 95% prediction interval (line through lowest diamond). CI=confidence interval

Fig 6
Fig 6

Calibration of the four models tested. Top row shows the alignment between predicted and observed risks for all models with smoothed calibration plots. Bottom row summarises the distribution of predicted risks from each model as histograms

Regional differences in the Harrell’s C index were relatively slight. None of the inter-region heterogeneity observed for discrimination (I2=53.14%) and calibration (I2=42.35%) appeared to be attributable to regional variation in any of the sociodemographic factors examined (table 3). The model discriminated well across cancer stages, but discriminative capability decreased with increasing stage; moderate variation was observed in calibration across cancer stage groups (supplementary table 9).

Table 3

Random effects meta-regression of relative contributions of regional variation in age, body mass index, deprivation, and non-white ethnicity on inter-regional differences in performance metrics after internal-external cross validation

View this table:

Competing risks regression

Similar fractional polynomial terms were selected for age and BMI in the competing risks regression model (see supplementary figure 2), and predictor selection yielded a model with fewer predictors than the Cox model. The competing risks regression model is presented as exponentiated coefficients in figure 7, with the full model (including constant term) detailed in supplementary table 5. Ethnic group specific discrimination and overall calibration metrics are detailed in supplementary table 4—the model generally performed well across ethnic groups, with similar discrimination, but there was some overt miscalibration on summary metrics—although some metrics were estimated imprecisely owing to small event counts in some ethnic groups.

Fig 7
Fig 7

Final competing risks regression model predicting 10 year risk of breast cancer mortality, presented as its exponentiated coefficients (subdistribution hazard ratios with 95% confidence intervals). Model contains fractional polynomial terms for age (1, 2) and body mass index (2, 2), but these are not plotted owing to reasons of scale. Model also includes an intercept term (not plotted—see supplementary file for full model as coefficients). CI=confidence interval; ER=oestrogen receptor; GP=general practitioner; HER2=human epidermal growth factor receptor 2; HRT=hormone replacement therapy; PR=progesterone receptor

The random effects meta-analysis pooled Harrell’s C index was 0.849 (95% confidence interval 0.839 to 0.859, 95% prediction interval 0.821 to 0.876). Some evidence suggested systematic miscalibration overall—that is, a pooled calibration slope of 1.160 (95% confidence interval 1.064 to 1.255, 95% prediction interval 0.872 to 1.447). Smoothed calibration plots showed underestimation of risk at the highest predicted values (eg, predicted risk >40%, fig 6). Supplementary figure 3 displays regional performance metrics.

An estimated 41.33% of the regional variation in the Harrell’s C index for the competing risks regression model was attributable to inter-regional case mix (table 3); ethnic diversity was the leading sociodemographic factor associated therewith (table 3). For calibration, the I2 from the full meta-regression model was 56.68%, with regional variation in age, deprivation, and ethnic diversity associated therewith. Similar to the Cox model, discrimination tended to decrease with increasing cancer stage (supplementary table 9).


Table 4 summarises the selected hyperparameter configuration for the final XGBoost model. The discrimination of this model appeared acceptable overall,68 albeit lower than for both regression models (table 2; supplementary figure 4), with a meta-analysis pooled Harrell’s C index of 0.821 (95% confidence interval 0.813 to 0.828, 95% prediction interval 0.805 to 0.837). Pooled calibration metrics suggested some mild systemic miscalibration—for example, the meta-analysis pooled calibration slope was 1.084 (95% confidence interval 1.003 to 1.165, 95% prediction interval 0.842 to 1.326). Calibration plots showed miscalibration across much of the predicted risk spectrum (fig 6), with overestimation in those with predicted risks <0.4 (most of the individuals) before mixed underestimation and overestimation in the patients at highest risk. Discrimination and calibration were poor for stage IV tumours (see supplementary table 9). Regarding regional variation in performance metrics as a result of differences between regions, most of the variation in calibration was attributable to ethnic diversity, followed by regional differences in age (table 3).

Table 4

Description of machine learning model architectures and hyperparameters tuning performed

View this table:

Neural network

Table 4 summarises the selected hyperparameter configuration for the final neural network. This model performed better than XGBoost for overall discrimination—the meta-analysis pooled Harrell’s C index was 0.847 (95% confidence interval 0.835 to 0.858, 95% prediction interval 0.816 to 0.878, table 2 and supplementary figure 5). Post-internal-external cross validation pooled estimates of summary calibration metrics suggested no systemic miscalibration overall, such as a calibration slope of 1.037 (95% confidence interval 0.910 to 1.165), but heterogeneity was more noticeable across region, manifesting in the wide 95% prediction interval (slope: 0.624 to 1.451), and smoothed calibration plots showed a complex pattern of miscalibration (fig 6). Meta-regression estimated that the leading factor associated with inter-regional variation in discrimination and calibration metrics was regional differences in ethnic diversity (table 3).

Stage specific performance and decision curve analysis

Both the XGBoost and neural network approaches showed erratic calibration across cancer stage groups, especially major miscalibration in stage III and IV tumours, such as a slope for the neural network of 0.126 (95% confidence interval 0.005 to 0.247) in stage IV tumours (see supplementary table 9). Overall decision curves showed that when accounting for competing risks, net benefit was generally better for the regression models, and the neural network had lowest clinical utility; when not accounting for competing risks, the regression models had higher net benefit across the threshold probabilities examined (fig 8). Lastly, the clinical utility of the machine learning models was variable across tumour stages, such as null or negative net benefit compared with the scenarios of treat all for stage IV tumours (see supplementary figure 6).

Fig 8
Fig 8

Decision curves to assess clinical utility (net benefit) of using each model. Top plot accounts for the competing risk of other cause mortality. Bottom plot does not account for competing risks

Clinical scenarios and risk predictions

Table 5 illustrates the predictions obtained using the Cox and competing risks regression models for different sample scenarios. When relevant, these are compared with predictions for the same clinical scenarios from PREDICT Breast and the Adjutorium model (obtained using their web calculators: and

Table 5

Risk predictions from Cox and competing risks regression models developed in this study for illustrative clinical scenarios, compared where relevant with PREDICT and Adjutorium*

View this table:


This study developed and evaluated four models to estimate 10 year risk of breast cancer death after diagnosis of invasive breast cancer of any stage. Although the regression approaches yielded models that discriminated well and were associated with favourable net benefit overall, the machine learning approaches yielded models that performed less uniformly. For example, the XGBoost and neural network models were associated with negative net benefit at some thresholds in stage I tumours, were miscalibrated in stage III and IV tumours, and exhibited complex miscalibration across the spectrum of predicted risks.

Strengths and limitations of this study

Study strengths include the use of linked primary and secondary healthcare datasets for case ascertainment, identification of clinical diagnoses using accurately coded data, and avoidance of selection and recall biases. Use of centralised national mortality registries was beneficial for ascertainment of the endpoint and competing events. Our methodology enabled the adaptation of machine learning models to handle time-to-event data with competing risks and inclusion of multiple imputation so that all models benefitted from maximal available information, and the internal-external cross validation framework28 permitted robust assessment of model performance and heterogeneity across time, place, and population groups.

Limitations include no consideration of genetic data such as presence of high risk mutations or multigene or multigenomics data, or breast density, which could have offered additional predictive utility.697071 Model development depended on using variables that are routinely collected in primary care, Hospital Episodes Statistics data, and the national cancer registry. Reliance on clinical coding for variables such as family history of breast cancer may be skewed towards those with more notable pedigrees; furthermore, as those without recorded positive family history were assumed to have none, misclassification might have occurred. Misclassification bias could also occur with prescriptions data because not all drugs are dispensed by a pharmacist or taken by the individual. Importantly, no coefficient in any model has a causal interpretation, and further work would be required to assess the relevance of altering factors such the management of type 2 diabetes,72 or accounting for treatment drop-in737475 after diagnosis if a causal prediction component was desired, or a counterfactual treatment selection function.

Another approach to model evaluation is bootstrapping, which allows estimation of optimism during model fitting, and calculation of bias corrected performance metrics. The best approach for combining bootstrapping with multiply imputed data may be to impute each individual bootstrap sample76—this would have been computationally intractable for this study, particularly for the machine learning models, which would have additional overhead of hyperparameter tuning in addition to imputation in each resample.

Comparisons with other studies

In a previous systematic review, the authors identified 58 papers concerning prognostic models for breast cancer.24 While the Nottingham Prognostic Index retained its performance in several external evaluations, some other models have performed less well on application to external datasets, such as in patients at the highest ranges of age and risk, which underscores the need for robust assessment of model performance.24 The PREDICT Breast model is endorsed by the American Joint Committee on Cancer, and the model is widely used around the world in clinical decision making about adjuvant chemotherapy—however, external evaluations suggest that PREDICT performs less well in older women, and in other subgroups such as women with large, oestrogen negative cancers,77 reinforcing the need for consideration of relevant subgroup performance in the prediction of breast cancer outcomes. More relevant to the present study is the lack of a reliable clinical prediction model for our outcome of interest, applicable across all women with breast cancer. The only model27 for this clinical scenario found to be at low risk of bias in a published systematic review25 likely had too small a sample size to fit a prediction model, unnecessarily dichotomised predictor variables, and the final model was selected after developing more than 35 000 other models.

A recent study to develop and externally validate the Adjutorium model based on automated machine learning focused on treatment recommendation in the postsurgical (adjuvant) setting, reporting that the model derived from the AutoPrognosis approach was superior for discrimination than a Cox proportional hazards model fitted to the same data.9 It is, however, notable that the comparison was the full complexity and flexibility of an ensemble machine learning model against a Cox model with no interactions and a simple, single, predetermined (ie, not data driven) non-linearity with a quadratic term for age. Our approach considered up to two fractional polynomial powers, which may be able to capture more complex non-linearities if present, considered regression model interactions, sought to identify optimal models for prognostication more generally in patients with breast cancer, and explored several performance metrics from the perspective of geographical and temporal transportability.

Previous studies discuss the adaptation of machine learning models to handle time-to-event data using jack-knife pseudovalues,587879 but in our study we performed comparative evaluation of the discrimination, calibration, and clinical utility of these models in large datasets. In the current study we also report a variation of an XGBoost algorithm that can handle competing risks. Recent developments in machine learning modelling approaches include DeepSurv or DeepHit as adaptations of time-to-event modelling,80 whereas our approach directly modelled risk probabilities. Extensions include complex model ensembles such as the Survival Quilts approach, where machine learning models are temporally joined to estimate risks over timespans.8 However, we opted for simpler model architectures that are arguably more transparent than meta-model ensembles, and they are conducive to a more robust validation strategy within a set computational budget. Furthermore, the added benefit of using complex models to obtain (at best) modest yields in an overall performance metric such as a C index, as has been shown in recent healthcare machine learning papers,81 is debatable. Although our goal was to develop robust models that effectively prognosticate for all breast cancers, a comparative evaluation of the PREDICT and Adjutorium models could have been an interesting analysis in the early breast cancer group treated with surgery. This was not, however, possible owing to systematically missing covariates in our dataset.

It should be noted that no single approach is always optimal for any modelling task—more flexible methods could have better performance in other scenarios if features and risk associations in the given data are complex. Results from this specific modelling scenario on relative performance of different approaches may not hold across all other prediction studies, mandating careful consideration of methodology should more than one modelling approach be used.

Implications of results and future research

This study shows how comparative evaluation of modelling techniques within an internal-external validation framework in large, clustered healthcare datasets may provide insight into relative strengths of different strategies for clinical prognostication. Regardless of the flexibility of modelling strategy used, all clinical prediction algorithms should be extensively evaluated and stress tested: showing that a model works overall is subservient to understanding if, where, and how a clinical prediction model will break down.


In this low dimensional clinical setting, a Cox model and a competing risks regression model provided accurate estimation of breast cancer mortality risks in the general population of people of female sex with breast cancer. Subject to independent external validation and cost effectiveness and impact assessments, the two models could have clinical utility, such as informing stratified follow-up regimens.4 It is possible that the most robust clinical applications could be attained with future integration of multimodal data, such as genomic markers.8283 Implementation of the models in another clinical dataset, such as one based on electronic healthcare records, may be possible. This should follow local validation and potential recalibration—discrimination could be similar on application to a different system, but calibration may be affected by local variations in rates or treatment practices. The models do not include predictors such as ethnicity or race or deprivation score, which would otherwise need adaptation or scaling to match local metrics. The included predictors would be available to clinical teams caring for women with breast cancer after initial diagnostic investigations and would need to be aligned with local coding systems (eg, SNOMED).

What is already known on this topic

  • Clinical prediction models are used for adjuvant treatment selection in early stage, surgically treated breast cancer, but models that can prognosticate well in breast cancer of any stage could more widely inform risk based follow-up strategies and support prognosis counselling or clinical trial recruitment

  • Most breast cancer prediction models have methodological limitations, are poorly reported, and are at high risk of bias—the only model deemed at low risk of bias in a recent systematic review for the unselected population of women with breast cancer likely had too small a sample size for development, and uncertain transportability

  • Considerable interest in machine learning for clinical prediction exists, but there has been criticism of model explainability, transparency, robustness of evaluation, fairness of comparisons, and risks of algorithmic bias

What this study adds

  • A Cox proportional hazards and a competing risks regression model may have utility for informing risk stratified clinical strategies for unselected women with breast cancer—these approaches were superior to the models developed using machine learning methods

  • An internal-external cross validation framework can be used to identify best performing modelling strategies by assessing model performance, performance heterogeneity, and transportability

  • Adapting machine learning models to handle censored data in the setting of competing risks can be achieved using a pseudovalues based approach

Ethics statements

Ethical approval

This study was approved by the QResearch Scientific Committee, the approval for which is from the East Midlands research ethics committee (reference: 18/EM/0400, project: OX129).

Data availability statement

To guarantee the confidentiality of personal and health information only the authors have had access to the data during the study in accordance with the relevant licence agreements. Access to the QResearch database is managed, overseen, and approved by the QResearch Scientific Committee.


We thank the patient and public involvement volunteers for their input into this project in the earlier stages.

We acknowledge the contribution of EMIS practices who contribute to QResearch®, EMIS Health, University of Nottingham, and the chancellor, masters, and scholars of the University of Oxford for expertise in establishing, developing, and supporting the QResearch database. The Hospital Episode Statistics datasets and civil registration data are used with permission of NHS Digital, which retains the copyright for those data. The Hospital Episode Statistics data used in this analysis are Copyright © (2021), the Health and Social Care Information Centre and re-used with the permission of the Health and Social Care Information Centre (and the University of Oxford); all rights reserved. This project involves data derived from patient level information collected by the National Health Service as part of the care and support of patients with cancer. The data are collated, maintained, and quality assured by the National Cancer Registration and Analysis Service, which is part of Public Health England (PHE). Access to the data was facilitated by the PHE Office for Data Release. The Office for National Statistics, PHE, and NHS Digital bear no responsibility for the analysis or interpretation of the data.


  • Contributors: JH-C and GSC are joint senior authors. AKC, GSC, and JHC led the study conceptualisation and statistical analysis plans. AKC drafted the first version of the study protocol, which all co-authors reviewed, provided input into, and approved. AKC undertook the data specification (reviewed by JHC), data curation and preparation, and analyses, and wrote the first draft of the paper. All other authors reviewed the manuscript, provided critical feedback, and approved the final version. AKC and JHC had full access to the study data. Advanced statistical input was from GSC and JHC. DD and SL provided clinical (oncological) input. AKC has final responsibility for submission and is the guarantor. The corresponding author attests that all listed authors meet authorship criteria and no others meeting the criteria have been omitted.

  • Funding: This project was funded through a Clinical Research Training Fellowship awarded to AKC by Cancer Research UK (C2195/A31310). We also acknowledge personal grant support from Cancer Research UK: C8225/A21133 to DD, C49297/A27294 to GSC, C5255/A18085 to JHC; and the National Institute for Health and Care Research (NIHR) Biomedical Research Centre. SP receives support as an NIHR senior investigator (NF-SI-0616-10103), and from the NIHR Applied Research Collaboration Oxford and Thames Valley. The funders had no role in considering the study design or in the collection, analysis, interpretation of data, writing of the report, or decision to submit the article for publication.

  • Competing interests: All authors have completed the ICMJE uniform disclosure form at and declare: JHC is an unpaid director of QResearch (a not-for-profit organisation that is a partnership between the University of Oxford and EMIS Health that supply the QResearch database) and is a founder and shareholder of ClinRisk and was its medical director until 31 May 2019 (ClinRisk produces open and closed source software to implement clinical risk algorithms including two breast cancer risk models implemented in the National Health Service into clinical computer systems.

  • The lead author (AKC) affirms that the manuscript is an honest, accurate, and transparent account of the study being reported; that no important aspects of the study have been omitted; and that any discrepancies from the study as planned (and, if relevant, registered) have been explained.

  • Dissemination to participants and related patient and public committees: We plan to invite the patient and public involvement volunteers to use the results to co-create lay summaries and other public facing materials for the purposes of dissemination.

  • Provenance and peer review: Not commissioned; externally peer reviewed.

This is an Open Access article distributed in accordance with the terms of the Creative Commons Attribution (CC BY 4.0) license, which permits others to distribute, remix, adapt and build upon this work, for commercial use, provided the original work is properly cited. See: