Ensuring comparisons of health-care providers are fair: development and validation of risk prediction models for critically ill patients

The current Intensive Care National Audit & Research Centre model

Risk predictions in the ICNARC model are calculated for each admission based on the following predictors:

• age in years at admission to the critical care unit

• location prior to admission to the critical care unit (emergency department, ward, theatre, other critical care unit, other acute hospital or not in hospital) and, for admissions directly from theatre, urgency of surgery (either elective/scheduled or emergency/urgent)

• CPR within the 24 hours prior to admission to the critical care unit

• the ICNARC Physiology Score – an integer score between 0 and 100 based on derangement in 12 physiological parameters during the first 24 hours following admission to the critical care unit

• primary reason for admission to the critical care unit

• some interactions between the ICNARC Physiology Score and primary reason for admission.

The ICNARC model is regularly recalibrated using CMP data to ensure ongoing fit.

The Scottish Intensive Care Society Audit Group database

The SICSAG has maintained a national database of patients admitted to adult critical care units in Scotland since 1995. Initially, only adult, general intensive care and combined high-dependency/intensive care units (critical care) units participated in the audit. More recently, specialist critical care units have joined the audit, with the result that, as of 2014, all adult, general and specialist critical care units in Scotland participate voluntarily in the audit. Data are collected prospectively using a dedicated software system. Annual data extracts are pooled centrally onto servers at the Information Services Division and validation queries relating to discharges, outcomes, ages and missing treatment information are then issued and fed back to individual units for checking by local and regional audit coordinators.

Use of the SICSAG database for this study was approved by the Privacy Advisory Committee, NHS National Services Scotland (application number 53/10).

Inclusion and exclusion criteria

Data were extracted from the SICSAG database for all admissions to all 24 adult general critical care units in Scotland between 1 January 2007 and 31 December 2009. During this study period, specialist cardiothoracic critical care units were not participating in the national audit; admissions to one specialist neurocritical care unit were not included in the data extract. The following admissions were excluded from the analysis: admissions flagged in the database as ‘exclude from severity of illness scoring’; readmissions of the same patient within the same acute hospital stay; admissions missing the outcome of status at discharge from acute hospital; admissions missing age, location prior to admission or primary reason for admission to the critical care unit; and admissions for whom the primary reason for admission was unable to be mapped onto the ICNARC coding method (see Primary reason for admission).

Application of the Intensive Care National Audit & Research Centre model to the Scottish Intensive Care Society Audit Group database

The most appropriate recalibration of the ICNARC model was selected based on the time period of the data included in the analysis – this was a recalibration undertaken in 2009 using CMP data from 194,892 admissions to 187 critical care units between 1 January 2006 and 31 December 2008.

In order to apply the ICNARC model to data from the SICSAG database, certain assumptions and recoding were required. These are detailed in the following subsections: Location prior to admission, Systolic blood pressure, Arterial pH, Neurological status and Primary reason for admission. After applying this recoding, the predicted risk of acute hospital mortality from the ICNARC model was calculated for each admission using standard algorithms developed for the CMP.

The Acute Physiology And Chronic Health Evaluation II model

The APACHE II model was selected as a comparator for this study, as it was the model in use in Scotland at that time. The SICSAG database does not include all the requisite fields to enable a head-to-head comparison against other, more recent, risk prediction models. The APACHE II model was originally developed using data from 19 critical care units in 13 US hospitals,9 and has subsequently been validated and recalibrated using UK data. 13,21 Risk predictions are calculated for each admission based on the following predictors:

• the APACHE II Score – an integer score between 0 and 71 comprising an Acute Physiology Score (0–60 points) based on derangement in 12 physiological parameters during the first 24 hours following admission to the critical care unit, age points (0–6) for age categories of ≤ 44, 45–54, 55–64, 65–74 or ≥ 75 years, and chronic health points (0–5) for very severe conditions in the patient’s medical history

• admission to the critical care unit following emergency surgery

• diagnostic categories based on the primary reason for admission to the critical care unit.

Values of predicted acute hospital mortality were supplied by the Information Services Division, calculated from the original published coefficients9 using the standard algorithms applied for routine reporting of the SICSAG audit results at that time.

Statistical methods

The ICNARC model was validated using measures of calibration, discrimination and overall fit, as described below. The validation was conducted in the full 3-year SICSAG database extract and for each year separately.

Discrimination was assessed by the c-index,22 which is equivalent to the area under the receiver operating characteristic curve. 23 Calibration was assessed graphically and tested using the Hosmer–Lemeshow test for perfect calibration in 10 equal-sized groups by predicted probability of survival. 24 As the Hosmer–Lemeshow test does not provide a measure of the magnitude of miscalibration and is very sensitive to sample size,25,26 calibration was also assessed using Cox’s calibration regression, which assesses the degree of linear miscalibration by fitting a logistic regression of observed survival on the predicted log-odds of survival from the risk prediction model. 27 Accuracy was assessed by the Brier score (the mean-squared error between outcome and prediction)28 and Shapiro’s R (the geometric mean of the probability assigned to the event that occurred),29 and the associated approximate R² statistics (termed the ‘sum-of-squares’ R² and the ‘entropy-based’ R², respectively), which are obtained by scaling each measure relative to the value achieved from a null model. 30

The performance of the ICNARC model was compared with that of the APACHE II model. The difference in c-index between the two models was assessed using the method of DeLong et al. 31 Confidence intervals (CIs) for observed acute hospital mortality were calculated using the method of Wilson. 32

Statistical analyses were performed using Stata/SE, version 13.0 (StataCorp LP, College Station, TX, USA).

Available data

Data were extracted from the SICSAG database for 29,626 admissions to 24 adult, general critical care units between 1 January 2007 and 31 December 2009. The following admissions were excluded: 3599 admissions (12.1%) flagged in the database as ‘exclude from severity of illness scoring’ (Table 1 contains the breakdown of reasons for exclusion); 1324 (4.5%) readmissions of the same patient within the same acute hospital stay; 173 (0.6%) admissions missing the outcome of acute hospital mortality; 869 (2.9%) admissions missing location prior to admission (n = 16) or primary reason for admission to the critical care unit (n = 864); and 392 (1.3%) admissions for whom the primary reason for admission was unable to be mapped. No admissions were missing age. This resulted in a cohort of 23,269 (78.5%) admissions for analysis.

Of the admissions flagged as ‘exclude from severity of illness scoring’, acute hospital mortality was reported for 3529 admissions (98.1%); 731 (20.7%) of these patients died before discharge from the acute hospital (see Table 1 for breakdown). It was not possible to include these patients in the analysis, even using statistical imputation methods to account for missing data, as insufficient predictor data were recorded. Owing to the large number of admissions flagged as ‘exclude from severity of illness scoring’, a post hoc analysis was undertaken to investigate the potential impact of such exclusions using CMP data (see Simulation of exclusion criteria).

Table 2 summarises the case mix and outcomes for the included admissions, overall and for each year. The mean age of admitted patients was 57 years, 56% were male, and two-thirds of patients were admitted for non-surgical reasons. These characteristics were relatively stable over the 3-year period. The distribution of predicted risk of acute hospital death from the ICNARC model (2009 recalibration) is shown in Figure 1. The mean predicted risk of death (expected acute hospital mortality) was 30.1%, which was very close to the overall observed acute hospital mortality of 29.7%. Compared with the CMP data set used to produce the 2009 recalibration of the ICNARC model, patients admitted to Scottish critical care units were, on average, slightly younger (median 61 vs. 64 years), less likely to be admitted following elective/scheduled surgery (10.5% vs. 23.3%) and had higher acute severity of illness (mean ICNARC Physiology Score 19.6 vs. 18.0).

FIGURE 1.

Distribution of predicted risk from the ICNARC risk prediction model (2009 recalibration) for admissions to Scottish critical care units, 2007 to 2009.

Model validation

The measures of model performance of the ICNARC model (2009 recalibration) compared to the APACHE II model are shown in Table 3. The ICNARC model outperformed the APACHE II model on all measures of model performance. The ICNARC model had substantially better discrimination (c-index 0.848 vs. 0.806; p < 0.001; Figure 2) and was also much better calibrated (Figure 3). Cox calibration regression showed an intercept and slope for the ICNARC model very close to the ideal values of 0 and 1, respectively. In contrast, the APACHE II model underpredicted both risk (intercept < 0) and variability (slope < 1). Performance of the ICNARC model remained consistent across the 3 years studied.

FIGURE 2.

Receiver operating characteristic curves for the ICNARC (2009 recalibration) and the APACHE II risk prediction models among admissions to Scottish critical care units, 2007 to 2009.

FIGURE 3.

Calibration plots showing observed against expected mortality in 10 equal-sized groups for the ICNARC (2009 recalibration) and APACHE II risk prediction models among admissions to Scottish critical care units, 2007 to 2009.

Inclusion and exclusion criteria

For the development data set, data were extracted from the CMP database for all admissions to cardiothoracic critical care units between 1 January 2010 and 31 December 2012. Patients aged less than 16 years and readmissions to the critical care unit within the same acute hospital stay were excluded.

The validation data set consisted of admissions to cardiothoracic critical care units between 1 January 2013 and 30 June 2014. The same exclusion criteria were applied.

Outcome and candidate predictors

The outcome for the risk prediction model was acute hospital mortality, defined as death before final discharge from acute hospital and including deaths after direct transfer to another acute hospital from the hospital housing the critical care unit.

Candidate predictors were chosen based on expert clinical opinion and availability in the CMP database. The candidate predictors included were as follows: age; sex; severe conditions in the past medical history; dependency prior to admission to acute hospital; CPR within 24 hours prior to admission to the critical care unit; location prior to admission to the critical care unit; highest heart rate; mean arterial pressure (from the blood pressure measurement with the lowest SBP); highest temperature (central measurement or, if none available, non-central + 0.5 °C); ratio of PaO₂ to fraction of inspired oxygen (FiO₂) (from the arterial blood gas with the lowest PaO₂); lowest arterial pH; partial pressure of carbon dioxide in arterial blood (PaCO₂) (from the arterial blood gas with the lowest pH); highest blood lactate concentration; highest urea value; highest creatinine value; lowest sodium value; highest potassium value; lowest haemoglobin value; lowest white blood cell (WBC) count; lowest platelet count; lowest total GCS score; and mechanical ventilation status.

Severe conditions in the past medical history were defined according to the APACHE II method9 and categorised as liver disease, renal disease, cardiovascular disease, respiratory disease, metastatic disease, haematological malignancy and immunocompromisation. Conditions must have been evident in the 6 months prior to admission to the critical care unit. Dependency prior to admission to acute hospital was assessed according to the ability to complete activities of daily living, categorised as the ability to live without assistance in daily activities, with some (minor or major) assistance with daily activities or with total assistance with all daily activities. Location prior to admission was categorised as theatre, ward (including intermediate care areas), other critical care unit or emergency department. Admissions from theatre were further categorised as following elective/scheduled or emergency/urgent surgery. Where indicated among candidate predictors as lowest or highest, physiological predictors were the lowest or highest value from the first 24 hours following admission to the critical care unit.

Handling of missing data

Although moderate (from 0.1% to 3.2%), missing data were imputed to address potential bias and loss of precision. Fully conditional specification (FCS)36 was used as the multiple imputation method. All the candidate predictors (with or without missing values)37 and the outcome,38 as well as auxiliary variables related to missingness,39 were entered into the imputation model. When required, simple or zero-skewness log-transformation for non-normality was used. Unless the rate of missing information is unusually high, there tends to be little or no practical benefit to using more than 10 imputations40 and so, in the following analysis, 10 repeat imputations were performed. The examination of the imputed data showed the distribution to be broadly similar to that of the observed data, indicating no obvious problems with the imputation process.

Model development

The distributions of all candidate predictors were explored in patients with and without the primary outcome.

For modelling continuous predictors, different approaches were considered, including fractional polynomials, restricted cubic splines and generalised additive models. The best functional form for each predictor was selected based on fit, plausibility, accuracy and clinical knowledge.

After appropriate functional forms were decided in the univariable setting, a full multivariable model, with all continuous and non-continuous predictors, was fitted to determine the association between the predictors and the outcome. This model was redefined by removing predictors with no significant global effect. To test predictors’ global significance and individual linearity, Wald tests (based on Wald statistics for pooled estimates) were applied. Predictors that were non-significant at a cut-off p-value of 0.1 were discarded. The model was refitted and the remaining predictors were retested. The process continued until all the predictors in the model were significant. Using the resultant model as a starting point, a parsimonious model was developed using a backward elimination strategy. At each step one predictor was dropped from the model while comparing the c-index22 and Brier score. 28 The 10 performance estimates (from the 10 multiply imputed data sets) were averaged and their variances pooled according to Rubin’s rules. 41 The least significant predictor was removed and the process continued until no predictors remained in the model. The final model was chosen to balance parsimony and model performance.

As the majority of patients were admitted to the cardiothoracic critical care unit directly from theatre following cardiac surgery, two expanded models with additional predictors and interactions were tested to improve model performance. The first expanded model included an additional predictor for the pathological or physiological process of the primary reason for admission to the critical care unit (e.g. congenital or acquired deformity, degeneration, dissection or aneurysm, obstruction) from the hierarchical ICNARC coding method among admissions following cardiac surgery. 20 The second expanded model tested interactions between the physiological predictors and cardiac surgery as the primary reason for admission. A nominal p-value of 0.001 was used to retain interaction terms in the model. The enhanced models were tested for improvements in discrimination and calibration.

In order to further evaluate the expanded models, net reclassification improvement (NRI) was determined. Reclassification has been proposed as a measure of utility or improvement in a risk prediction model. 42 The proportions of patients with and without the outcome reclassified into lower- or higher-risk categories are compared. The NRI is defined as the proportion of non-survivors moving to a higher-risk category minus the proportion moving to a lower-risk category plus the proportion of survivors moving to a lower-risk category minus the proportion moving to a higher-risk category.

The final model coefficients were estimated using Rubin’s rules, to give a single estimate and standard error.

Model validation

The risk prediction model was then further validated in the temporally distinct validation data set. The 10 performance estimates of the final model were averaged and their variances pooled according to Rubin’s rules. The predictive performance of the model was estimated by bootstrapping the c-index and Brier score. 43 Calibration was assessed by Cox calibration regression27 and graphically using calibration plots, with 20 equal-sized risk groups. Using Rubin’s rules, 10 linear predictions were calculated and averaged from the new model equation for each admission. The predicted probability of acute hospital mortality was calculated from this pooled result.

Statistical analyses were performed using Stata/SE, version 13.0.

Available data

Between 1 January 2010 and 31 December 2012 there were 17,002 eligible admissions to five cardiothoracic critical care units participating in the CMP, which formed the development data set, and between 1 January 2013 and 30 June 2014 there were a further 10,238 eligible admissions to six cardiothoracic critical care units (one additional unit having joined the CMP), which formed the validation data set (Table 5).

In the development data set, the majority of admitted patients were male (69%) with a median age of 66 years. Only 14.9% of admitted patients had any previous severe conditions in the past medical history, with the majority of those present being due to severe cardiovascular disease (8.6% of admissions). Most patients were fully functional prior to hospital admission (82.3%) and reported as needing no assistance with their activities of daily living. Over three-quarters of all admissions were surgical, of which 97% followed cardiothoracic surgery. Most surgery was elective or scheduled and most patients were ventilated during the first 24 hours following admission to the unit. The median length of critical care unit stay was 1.2 days, while the median total length of stay in acute hospital was 11 days. Critical care unit mortality was 7.4% and acute hospital mortality was 11.1%.

The validation data set had similar characteristics, although a lower proportion of admissions had one or more severe conditions in the past medical history (9.9%), particularly severe cardiovascular disease (4.8%), and a correspondingly higher proportion were fully functional prior to hospital admission (88.1%). Mortality was also lower, both in the critical care unit (6.4%) and in acute hospital (9.7%).

Model development

The best functional form for the relationship between each of the 15 continuous predictors and the outcome of acute hospital mortality was explored (see Table 6). All 15 predictors showed significant non-linearity (p < 0.001). Restricted cubic splines were chosen for the final modelling as they showed the best combination of flexibility and precision. To avoid overfitting spurious dips and unrealistic features of the curve, four knots were chosen to model the continuous predictors. When functional form was reassessed in the full multivariable model, the evidence for non-linearity in the relationship for blood lactate concentration was weak and so this predictor was finally analysed as linear. The final functional form for each predictor, including the position of the knots for restricted cubic splines, is shown in Table 6.

The full multivariable model, including the 22 predictors, had a c-index and Brier score of 0.902 and 0.0635, respectively (Table 7). Of the initial 22 predictors, 18 were found to be associated with acute hospital mortality on multivariable analysis (p < 0.1): age; one or more severe conditions in the past medical history; dependency; CPR within 24 hours prior to admission; location prior to admission; heart rate; mean arterial pressure; temperature; PaO₂/FiO₂; arterial pH; PaCO₂; blood lactate concentration; urea level; creatinine level; sodium level; WBC count; platelet count; and GCS score. Removal of the non-significant predictors resulted in minimal change to c-index and Brier score (see Table 7).

The 18 significant predictors were entered into a stepwise model selection (see Table 7). The model which best balanced parsimony with precision consisted of 10 predictors: age; dependency; location prior to admission; mean arterial pressure; arterial pH; blood lactate concentration; creatinine level; WBC count; platelet count; and GCS score. The c-index and Brier score were 0.895 and 0.0656 respectively.

The first expanded model, incorporating reason for admission, performed moderately better than the baseline parsimonious model, with a c-index of 0.899 and Brier score of 0.0652 (see Table 7). The second expanded model, incorporating interactions between admission following cardiac surgery and blood lactate concentration, creatinine level and platelet count, demonstrated a c-index of 0.899 and Brier score of 0.0649 (see Table 7).

After comparing the reclassification of the two expanded models using risk categories defined by thresholds of 0%, 2%, 5%, 10%, 20% and 50% (Tables 8 and 9), the model with interaction terms was superior (Table 10). With this model, a total of 3677 (23%) admissions were reclassified and 2382 of those (65%) were placed in more appropriate categories. The total NRI for the expanded model with interaction terms was 11.1% (standard error 1.1%; p < 0.0001) compared with 6.5% (1.0%; p < 0.0001) for the expanded model with reasons for admission. The calibration regression for the expanded model with interaction terms demonstrated a slope of 0.98 and an intercept of −0.07, indicating a well-calibrated model. This was therefore taken as the final model. The coefficients for the final model are shown in Table 11.

Model validation

The performance in the validation data set of 10,238 admissions from January 2013 to June 2014 was excellent: a c-index of 0.904 (95% CI 0.893 to 0.915) and Brier score of 0.055. The calibration of the model was satisfactory (Figure 4), with a calibration slope of 0.961 and a calibration intercept of –0.183.

FIGURE 4.

Calibration in the validation data set of the final risk prediction model for admissions to cardiothoracic critical care units.

Missing data mechanism

Before missing values were imputed, we studied the missing data mechanism40 by creating indicator variables for missing values for each predictor and fitting logistic regression models with the indicator variable as the outcome and the other predictors as covariates.

Findings suggested that it was not plausible to assume that data were missing completely at random (MCAR) so it was assumed that data were missing at random (MAR), conditional on observed data; the plausibility of the MAR assumption holding was improved by conditioning on auxiliary variables that were found to be predictive of missingness in the imputation model. Furthermore, some methodologists have argued that routine departures from MAR may not be large enough to cause serious bias in the resulting estimates. 71

Imputation of missing data

We imputed missing data using FCS. 36 This method provides a practical and flexible approach to generating imputations based on a set of imputation models, one for each variable with missing values, allowing one to specify the regression equation for the imputation; this is usually linear regression for continuous variables and logistic regression (binary, ordinal, or unordered multinomial) for categorical variables. In addition, under logistic imputation, imputed values for categorical variables will also be categorical, so rounding to plausible values is not necessary. In large data sets (such as the CMP database), it is common for missing values to occur in several variables, both continuous and categorical, so this approach was chosen because of its ability to impute both continuous and categorical variables appropriately, permitting a great deal of flexibility.

The imputation model and analysis model should be compatible; that is, any relationship in the analysis model should also be part of the imputation model. All the potential predictors that will be considered in the analysis, with or without missing values, should therefore be included in the imputation model. 37 The response variable38 and auxiliary variables related to missingness39 were entered into the imputation model as well.

When required, transforming for non-normality gave better imputed values,72 so logarithmic transformation was used. Urine output had a particularly unusual distribution, with a proportion of values equal to zero and a continuous, but heavily skewed, distribution among the remaining values. A semi-continuous approach to imputation is recommended for this situation. 73 This involves separately imputing a binary variable that indicates whether a value is zero or positive, and then a continuous variable for the value if positive. However, this approach did not work well in the case of urine output, as it overestimated the proportion of zero values; finally, the shifted log-transformation f(z) = ln(± z − a) was chosen for this predictor. Different approaches were considered for imputing sedated/paralysed/GCS score (ordinal, multinomial or predictive mean matching). All produced similar results, so finally the ordinal approach based on logistic regression was used.

As well as including all variables planned to be included in the analysis model, the imputation model must also include them in an appropriate way, that is, in the correct functional form and with any interactions that are required. However, when imputing data, one does not necessarily know what terms will be required in a sequence of analyses, and allowing for all possible terms might make the imputation model impractically large. On the other hand, the simplest approach of passive imputation (i.e. using simple linear and logistic regression models and ignoring interactions and non-linearity that are in the analysis model) could result in biased terms in the analysis model. 74 To assess the impact of omitting interactions, we produced a provisional and relatively simple imputation model, including non-linear terms and any interactions of key scientific/clinical interest. The imputed data were then used to build and check an analysis model, investigating the need for non-linear terms in the imputation model as well as the best way to impute them. The following approaches for imputing non-linear terms were explored: the traditional, passive (‘impute then transform’) method;75 ‘just another variable’ (JAV; also termed ‘transform then impute’);75,76 and substantive model compatible FCS (SMC-FCS). 77 SMC-FCS is a modification of the FCS approach to multiple imputation which ensures that each of the imputation models is compatible with the assumed substantive model. 77

In our imputation process, most of the relationships of the physiology variables followed the same pattern, so a more general model was specified for the candidate predictors with missing data.

Simulation studies have shown that the required number of repeated imputations (m) can be as low as 3 for data with 20% of entries missing. 41 We used a more conservative choice of m = 5.

Validation of imputed data

Graphical assessment and differences in means and proportions were used to compare the distributions of the imputed and observed data. Examination of the imputed data showed that the distribution was broadly similar to that for the observed data, indicating no obvious problems with the imputation process. A few imputed values were placed outside of the distribution of observed data; however, these were too few to be considered important, and some differences are to be expected if the data are not MCAR.

Simulation study

Imputing missing data, especially multiple imputation, is becoming standard in risk prediction modelling. 78 However, it is often implemented without adequate consideration of whether or not it offers any advantage over complete case analysis for the research question of interest. Multiple imputation is not always better than complete case analysis for missing covariate problems. 79 Recovering information in estimating the coefficient could be a potential gain, but previous studies have demonstrated that, although it may be important to use multiple imputation to recover information when there are missing data in covariates required for adjustment, multiple imputation has substantially less value when there are missing data (even when MAR) in the exposure of interest. 80

We explored the impact of missing data in the planned analysis using a simulation study under two different scenarios by creating data sets with MCAR and MAR missing values from 300,000 complete data values (a random sample drawn from the original complete case data: see Appendix 1). To test the consistency of the finding in these scenarios, logistic regression models were estimated on the original data set with no missing values (used as a reference) and the complete case data set after creation of missing values, as well as on multiply imputed data sets using FCS. Results were compared in terms of coefficients, standard errors and p-values.

If we impute from a model that does not allow for a potential non-linear association with outcome, we would expect to obtain inconsistent parameter estimates;74 therefore, in a secondary analysis we also investigated, under our these two scenarios of MCAR and MAR, the non-linearity effect on the risk prediction model estimation and the reliability of the prediction using the different approaches as mentioned above.

Results and conclusions

It was expected that multiple imputation would raise the risk estimates in comparison with the complete case analysis, but only slight differences in estimates were found, even among non-response mechanisms (Table 17). The increase in precision of risk estimates under MAR is therefore minimal at best, and the risk estimates hardly change.

Under our scenarios of missing data that with either MCAR or MAR, both complete case analysis and multiple imputation had negligible bias compared to the reference results based on the data set prior to the creation of missing values. The results indicated that a multiple imputation-based method (FCS) produced similar estimates to the complete case analysis; therefore, little information was gained regarding the coefficients for the predictors with missing values when we imputed those missing values, regardless of the number of missing data (at the levels observed within our data sets). Similar findings of inconsistent benefits of multiple imputation were observed in the coefficients for the predictors with no missing values. No differences in the inference were found. The different approaches to imputation of non-linear terms gave similar results (Table 18), including the passive (‘impute then transform’) approach, JAV (‘transform then impute’) and SMC-FCS, although SMC-FCS gave slightly smaller standard errors. These findings could indicate that, in our scenarios, including non-linear terms into the imputation model is not necessary.

In conclusion, coefficient estimates appear to be insensitive to the missing data and the various models used to deal with them. The benefits of using multiple imputation in developing our risk prediction model are likely to be minimal.

However, there are some factors contributing to the apparent stability of the estimates, such as the moderate number of missing data, the large sample size and the fact that the differences in mortality between complete and incomplete are too small to exert a serious impact on the estimates. Another factor might be that imputations do not contain important information or that the covariates are missing not at random; in this circumstance multiple imputation may be biased while complete case analysis may not be.

Finally, we decided that the model building and analysis process would be done with non-imputed (complete case) data and that a parallel analysis would be done at the same time on the multiply imputed data set in order to test the consistency of the results.

Functional form

Linearity of the association between continuous predictors and the outcome should not automatically be assumed because this could lead to incorrect interpretation of the effect and inaccurate predictions when the model is applied to new individuals. 81 In addition, the use of the appropriate functional form for a continuous variable is crucial for valid predictions because the expected value of the outcome can be different for the same value of a continuous variable with different functional form. As many of the physiological predictors are known to have non-linear (often U-shaped) relationships with outcome, they have typically previously been modelled using categorical approaches and combined into a severity score. 9–11,14,82 One aim of the current project was to model the physiological predictors separately as continuous variables, so how to deal with non-linearity was carefully considered.

The predictor–outcome relationship was explored by expanding the variable into multiple terms and testing pooled and individual non-linearity. The shape of the relationship of continuous predictors with the outcome was also studied, by plotting with a smoothed curve (running line smoother) as a reference. The traditional use of a categorical approach failed to detect the continued increase or decrease in risk for subjects at higher/lower levels of risk, making the implausible assumption that the risk does not vary at the extremes. Moreover, ignoring intracategory variation means throwing away information, tends to reduce a study’s power to detect an association and may lead to inaccurate estimates. 83,84 After the hypothesis of linearity was rejected, two different approaches for modelling continuous predictors were considered: non-linear functions such as second-order fractional polynomials85 and restricted cubic splines. 84

Spline fits could be sensitive to the number of knots so, to avoid overfitting, spurious dips and inflection points, as well as unrealistic features of the curve, three, four or five knots were considered. With large numbers of observations, three or four knots may work better, and more than five knots are usually not necessary unless the response to the predictor is extremely complicated. Knot positions were selected according to the recommendations of Harrell. 84 Right-restricted cubic splines (i.e. with the linearity restriction applied only at the right-hand end of the curve) were used when appropriate to allow more flexibility.

For optimal power fractional polynomials, we selected the family of second order because this offers considerable flexibility and a very rich set of possible functions, including U- and J-shaped relationships (an order higher than 2 is rarely required in practice86).

In order to judge the plausibility and accuracy of the fitted curves, we plotted observed log-odds against the alternative modelling approaches and used a running line smoother as a reference. The Akaike information criterion (AIC) and Bayesian information criterion (BIC) were calculated to compare the fit of the strategies for modelling continuous variables, taking three-knot-restricted cubic splines as a reference point, to assess whether or not the increased complexity of the model that resulted from including more knots was worthwhile.

The best functional form for each predictor was selected based on fitting, plausibility, accuracy and prior knowledge about the predictor.

Finally, we explored collapsing extreme points to determine whether or not the shape of the curves could be affected by outlying values. Additional analysis using imputed data was done in parallel to provide reassurance about the results.

Modelling of interactions

Improved modelling of interactions between physiological predictors was one of the main project objectives. The assumption that all predictors act independently on outcome is physiologically untrue and, although most of the significant past medical history–physiology, intervention–physiology and physiology–physiology interactions were not expected to make an important contribution to the model, these interactions could reflect important dependencies as well as intrinsic prognostic features of the data. Clinical members of the Expert Advisory Group, together with additional clinicians with relevant experience of risk prediction models in clinical epidemiology, identified and prioritised key, potential interactions between predictors. These pre-defined past medical history–physiology, intervention–physiology and physiology–physiology interactions were tested and checked graphically for any spurious interactions as a result of mis-modelling the main terms in the model. We tested each interaction’s importance in predicting the outcome using the BIC in order to reduce the chance of overfitting. Only significant interactions were considered in the model-building process.

In addition, the relationship between physiology and outcome may depend on primary reason for admission to the critical care unit, and a model incorporating interactions could take into account these differences. In the current ICNARC model, patients have a different coefficient for their ICNARC Physiology Score depending on primary reason for admission but, in the new model, these interactions were addressed for each individual physiological predictor.

Reason for admission

Reasons for admission to critical care are recorded in the CMP using the ICNARC coding method: a five-tiered (type – surgical or non-surgical/body system/anatomical site/physiological or pathological process/condition) coding system specifically developed for this purpose. 20 Currently, coefficients for the ICNARC model are applied at only two levels of the hierarchical code (either at tier 5, the individual condition, or at tier 2, the body system), but we aimed to improve the overall performance by incorporating intermediate levels of information and, in that way, to allow for variation in the effect on mortality between different physiological or pathological processes (e.g. infection or trauma).

The following steps were taken in modelling primary reason for admission:

1. If a process did not have the required sample size (number of events < 20):

   1. it was added to a related process or

   2. it was combined with other similar processes to create a new group.

2. Indicators were created for each combination of process and system.

3. If a process/system combination did not have the required sample size (number of events < 20):

   1. it was added to a related process/system combination or

   2. it was combined with other similar process/system combinations to create a new group or

   3. it was combined with the same process across all systems as a single process category.

4. When a process/system combination either was non-significant after adjusting for process or did not make an important contribution to the model:

   1. it was added to the process category.

5. The resulting set of categories was redefined to make a grouping with clinical sense.

6. Finally, the set of categories was refined by adding specific conditions:

   1. Process/system categories were split into individual conditions which had sufficient sample size (number of events ≥ 20).

   2. Each individual condition was retained as a new category if it was significant after adjusting for process/system and made an important contribution to the model.

   3. Individual conditions were combined with other individual conditions within each process/system category if the difference in risk was not significant.

Model building process

The model building process consisted of a number of stages (Figure 5).

FIGURE 5.

Stages in the model building process for the risk prediction model for adult critical care units.

Results

The predicted risk from continuous approaches to modelling physiology was a much better approximation of the true risk function than the previous categorical approach (Figure 6), especially if the factor distribution was skewed and, in particular, for relationships with a single turning point, such as for SBP, heart rate, arterial pH, WBC count, potassium level and haemoglobin level.

FIGURE 6.

Alternative approaches to modelling continuous physiological predictors. FP, fractional polynomial; RCS, restricted cubic spline; RRCS, right-restricted cubic spline.

In general, both methods for continuous modelling (fractional polynomial models and restricted cubic splines) agreed closely, but variation between the two approaches occurred in the regions with sparsest data. Although the data contain little information on mortality, the fractional polynomial fit was poorer at extreme predictor values where a strong non-monotonic relationship was suggested (e.g. arterial pH and creatinine level), whereas restricted cubic splines accommodated a more realistic final trend at the ends of the curve. In addition, the SBP, heart rate and neutrophil fitted curves were more sensitive to extreme values using fractional polynomials, while the ensured linearity in the tails by restricted cubic splines avoided this unrealistic end effect and seemed not to be affected by the collapsed values. We decided to use restricted cubic splines because they showed the flexibility of fractional polynomials, but with better behaviour in the tails, they captured the most prominent features of the relationship between predictors and outcomes and because the fit was more plausible than that of the previous categorical approach.

The optimal functional form selected to model the continuous physiological predictors (Table 19) was four knots for heart rate, SBP, temperature, PaO₂/FiO₂, arterial pH, blood lactate concentration, urine output, urea level, potassium level, glucose level, haemoglobin level, WBC count, neutrophil count and platelet count because this accommodated the trend of decreasing/increasing mortality in a non-monotonic way, as well as capturing the behaviour at extremes. For PaCO₂ and sodium level, a simplification using three knots was enough to accommodate the non-linear behaviour and had better AIC and BIC than four knots. Respiratory rate and creatinine level were modelled using right-restricted cubic splines. This was necessary to capture the initial decrease in mortality and ‘spoon’ behaviour (Figure 6). This approach had better fit than five knots and was more plausible than four. Finally, the stability of the fitted curves and the knots was assessed in each of the 10 imputed data sets, showing similar results.

Consistency with both optimal functional form and global significance of the predictors was found after fitting a model adjusted for all physiological predictors from the current ICNARC model.

A more parsimonious physiology model was then developed using a backward process in order to balance parsimony and model performance (Table 20). All candidate physiological predictors, including those not incorporated in the current ICNARC model, were entered into a full model, which was then simplified. Three new candidate physiological predictors (PaCO₂, blood lactate concentration and platelet count) were retained, as well as all physiological predictors from the current ICNARC model.

Changes were made to a number of the current non-physiological predictors, in particular CPR within 24 hours prior to admission, which was expanded to three categories (in-hospital CPR, community CPR and no CPR), and source of admission, which was expanded into nine categories combining information from admission type. After exploring the relationship with outcome, deprivation was collapsed into two categories (quintiles 1–3 and quintiles 4 and 5) and dependency prior to hospital admission was collapsed into three categories [no dependency, partial dependency (combining either minor or major assistance with activities of daily living) or total dependency]. After adjusting for current and new potential predictors, deprivation and body mass index (BMI) were not retained in the main-terms model because of their lower contribution to the model fit. The new predictor of dependency prior to hospital admission was retained. Severe conditions in the past medical history were added into the model, but the only conditions ultimately retained were severe liver disease, owing to its interaction with physiology, and metastatic disease and haematological malignancy, because of their significant effect. Ventilation was also retained because it is the main term in several interactions. This main-terms model had a c-index of 0.8779 and a Brier score of 0.1098 (Table 21).

A total of 56 process/system combinations and 16 individual conditions from the ICNARC coding method form the new reason for admission categories (Table 22). These were selected after the modelling process described above and accounted for 93.6% and 6.4% of admissions respectively.

The incorporation of primary reason for admission categories, plus interactions with physiology (n = 19), in the risk prediction model produced better fit (see Table 21) and, therefore, could reduce bias when applied across different settings with differing case mix. 15 One of the pre-defined past medical history–physiology interactions, six of the pre-defined intervention–physiology interactions and three of the pre-defined physiology–physiology interactions were retained after adjusting for the main term model plus significant primary reason for admission–physiology interactions. Apart from ventilation with arterial blood gas results, most of the significant interactions added did not make an important contribution to the model (Table 23); however, these interactions could reflect important dependencies as well as intrinsic prognostic features of the data.

The addition of non-physiological predictors and interactions did not affect the functional form of the physiological predictors.

The estimates for the model parameters obtained using data from the multiply imputed data set were similar to values estimated from the development data set and, therefore, the bias that could arise from using only the available information was considered to be very small.

Following the development process described above, the significance and importance of the predictors in the final model are shown in Table 23. Full coefficients for the final model [i.e. the new ICNARC model for acute hospital mortality (ICNARC_H-2014)] are presented in Appendix 2. The distribution of predicted acute hospital mortality from the new model is shown in Figure 7.

FIGURE 7.

Distribution of predicted acute hospital mortality from the new ICNARC model in the development data set.

Among the subset of patients with pupil reactivity recorded (n = 104,063), a model using pupil reactivity in place of the sedated/paralysed/GCS score predictor had very similar discrimination (c-index 0.8883 vs. 0.8887).

Assessing the predictive performance

Discrimination, calibration and accuracy are key aspects of the predictive performance of prediction models. Discrimination is the ability of a model to distinguish individuals who experienced the outcome from those who remained event free, and calibration is the agreement between the probability of developing the outcome of interest as estimated by the model and the observed outcome frequencies. Accuracy refers to the difference between predictions and observed outcomes at the level of individuals.

The discrimination of the model was estimated by the c-index22 (equivalent to the area under the receiver operating characteristic curve23) and accuracy was assessed by the Brier score (mean-square error between outcome and prediction). 28 We assessed calibration graphically with predicted probability on the x-axis and the observed outcomes on the y-axis in 10 equal-sized risk groups (calibration plot) and by Cox’s calibration regression (linear recalibration of the predicted log-odds). 27

The use of calibration plots and parameters from Cox calibration regression (intercept and calibration slope) to assess the calibration of the model, and not the Hosmer–Lemeshow c-statistic,24 is because the last statistic reflects the statistical significance of miscalibration and not its magnitude and hence, in general, models validated in a large sample and/or with the ability to provide very low mortality risk predictions will tend to have a worse Hosmer–Lemeshow c-statistic than models validated in a small sample size and/or providing a smaller range of mortality risks. 25,26

Some researchers have proposed reclassification level as a measure of added utility or improvement. Cook recommends producing a ‘reclassification table’ to show how many subjects are reclassified by adding a marker to a model. 91 Pencina et al. extended the reclassification idea by conditioning on the outcome; therefore, reclassification of subjects with and without the outcome should be considered separately. 42 Any ‘upwards’ movement in categories for subjects with the outcome implies improved classification, and any ‘downwards’ movement indicates worse classification. The interpretation is opposite for subjects without the outcome. The improvement in classification was quantified as the sum of differences in proportions of individuals moving up minus the proportion moving down for those with the outcome, and the proportion of individuals moving down minus the proportion moving up for those without the outcome (NRI).

Internal validation

Validation of predictive models, either internally to adjust for optimism and overfitting or externally to assess generalisability, is essential in terms of understanding the reliability of both the choice of variables and the values of coefficients for each variable. The performance of a predictive model could be overestimated when simply determined on the sample of subjects that was used to construct the model; several internal validation methods are available that aim to provide a more accurate estimate of model performance in new subjects and to estimate the potential for optimism and overfitting in model performance. Bootstrapping is now widely regarded as a better approach to validation than data splitting18,43,88 and the c-index is often used to indicate optimism if the value decreases substantially in an independent data set. Harrell et al. 88 presented an algorithm for estimating the optimism or overfitting in predictive models. Their method is based on using bootstrapping to derive the expected bias-corrected c-index, which is referred to as the optimism-corrected c-index.

When a predictive model is based on a very large sample size, as in this project, and relevant variables are included in the final model, reasonable estimates of the coefficient values for each variable are likely. 92 Optimism is small and so the apparent estimates of model performance (c-index and Brier score in the development data set) are attractive because of their stability. 43 However, to assess optimistic performance, the percentage of overfitting was estimated by the optimism-corrected c-index.

External validation

External validation is an important issue because the performance of most developed and internally validated prediction models, when applied to new individuals, is poorer than the performance seen in the sample from which it was developed. Internal validation does not make use of data other than the development data and, therefore, will not provide the degree of heterogeneity that will be encountered in real-life application of the model.

External validation of any type consists of taking the original model, with its predictors and assigned weights estimated from the development data set, and applying it to a different data set, obtaining the measured predictor and outcome values in the new individuals and quantifying the model’s predictive performance. In this project, external validation was done for a different period of time and in different specialist units. We used the validation data set (January to September 2013) to compare the predictive accuracy of the new model with that of the current ICNARC model, both overall as well as across different types of admissions (surgical vs. non-surgical and planned vs. unplanned). An objective of this project was to develop a single, general risk prediction model that could be applied across all adult critical care units, including specialist cardiothoracic and neurocritical care units and stand-alone high-dependency units. However, it is possible that a single, general model may under- or overestimate mortality in selected admission subpopulations or different unit types and so could show worse performance than specifically designed or calibrated models. Therefore, in addition, we used the validation data set to evaluate the use of a general model in specialist units by comparing it to versions of the current ICNARC model specifically recalibrated to each unit type.

Results

The final model showed good performance (c-index 0.8909 and Brier score 0.1027) and internal calibration of the model was satisfactory (see Table 21 and Figure 8). Overfitting was of limited relevance because of the very large data set and, as expected, model optimism was negligible (0.16% estimated overfitting).

FIGURE 8.

Calibration of the new ICNARC model in the development data set.

The performance in the validation data set (January to September 2013) is presented in Table 24 and Figure 9. Compared with the current ICNARC model, the new model demonstrated small improvements in discrimination and accuracy (c-index 0.8853 vs. 0.8693, Brier score 0.1076 vs. 0.1146). The current ICNARC model tended to slightly overestimate overall mortality; observed and mean predicted mortality were 22.6% and 23.2%, respectively, for a standardised mortality ratio (SMR, defined as observed divided by predicted mortality) of 0.97 (95% CI 0.96 to 0.99). A total of 30,739 (42.8%) of the admissions were reclassified by the new model (31.4% of survivors and 17.2% of non-survivors; Tables 25 and 26) and 20,252 of those (65.8%) were placed into more appropriate categories. The total NRI for the new model was 19.9 (p < 0.0001) and the new model showed a better risk classification, allocating more admissions in the extremes (Table 27).

FIGURE 9.

Calibration of the new ICNARC model compared with the current ICNARC model (2013 recalibration) in the external validation data set. a, 2013 recalibration; and b, developed in this project.

When compared across different types of admissions (surgical vs. non-surgical and planned vs. unplanned), the new model improved discrimination (Table 28) and calibration (Figure 10) compared with the current ICNARC model. For surgical admissions, the new model showed excellent discrimination and calibration (c-index 0.8905, SMR 1.00) and for all other subgroups, good discrimination and SMRs close to 1, which indicated good calibration among these subgroups.

FIGURE 10.

Calibration of the new ICNARC model compared with the current ICNARC model (2013 recalibration) on patient subgroups. (a) Surgical; (b) non-surgical; (c) planned; and (d) unplanned. a, 2013 recalibration; b, developed in this project.

To evaluate the performance of the new ICNARC model in specialist units, it was tested in each type of unit using the January–September 2013 validation data set and was compared with specific, recalibrated versions of the current ICNARC model for these units and, for cardiothoracic critical care units only, the new specific model developed in Chapter 3 of this report (Table 29). The new model showed, in general, a better performance than the specific, recalibrated models. Overall performance improved: the Brier score decreased and discriminative ability showed an increase compared with the specific risk prediction models, except for in the case of cardiothoracic units, where discrimination was similar but calibration was very good (Figure 11). The observed mortality was uniformly lower than the predicted mortality (SMR < 1) in both new and specific recalibrated risk prediction models, except for the new risk model for cardiothoracic critical care units, which gave slightly lower predicted mortality (SMR 1.04), and the new ICNARC model when applied to neurocritical care units, which showed excellent calibration (SMR 1.00). Across all types of unit, the group of admissions with the highest observed mortality had higher predicted risk according to the new model than according to the specific recalibrated models, corresponding to more extreme predictions and a greater spread of predicted risk (Figure 11).

FIGURE 11.

Calibration of the new ICNARC model compared with recalibrations of the current ICNARC model in different types of critical care unit. (a) Cardiothoracic critical care units; (b) neurocritical care units; and (c) high-dependency units. a, 2013 recalibration (specific to unit type); b, developed in this project (single model for all unit types).

Conclusions

The new ICNARC model showed good discrimination and calibration and had good performance when externally validated using data from a different time period. In addition, it improved the discrimination and reclassification compared with the current model. When validated on specialist units, both discrimination and calibration were close to that of specific recalibrated versions of the current ICNARC model; therefore, use of this single model, which works across general and specialist units, is recommended.

Methods

Results

Influence of imputation on benchmarking

When different approaches (single imputation, multiple imputation and normal/mean imputation) were compared, there was no significant difference in overall SMR between the approaches (0.99, 95% CI 0.97 to 1.00; 0.99, 95% CI 0.98 to 1.00; and 0.99, 95% CI 0.98 to 1.00 respectively) or with complete case analysis (0.99, 95% CI 0.97 to 1.00) compared with the reference value (0.99, 95% CI 0.97 to 1.00).

Table 32 and Figure 12 present the critical care unit performance according to original, complete case and imputation approaches. The proportion of critical care units performing above and below expectation, according to 95% control limits, is higher after the imputation approaches but close to the reference, varying from 13 (6.0%), with single and multiple imputation, to 14 (6.5%) with normal/mean imputation. Distribution of the SMRs under each performance category was similar to the reference for all imputation approaches with a high agreement (kappa 0.927, 0.963 and 0.927 for single imputation, multiple imputation and normal imputation respectively). On the other hand, complete case analysis led to a substantially higher number of critical care units performing ‘in control’ [193 (90%) vs. 186 (86.5%) for the reference]. Most of these critical care units shifted from in control (based on complete case) to positive alert after imputation, although complete case still showed a good agreement with the reference (kappa 0.724). No positive alarms were found in the data.

TABLE 32 - Critical care unit performance status for different approaches to handling missing data in the simulation study for application of the new ICNARC model

Approach to missing data	Location of critical care units relative to funnel plot lines, n (%)
	Negative alarm	Negative alert	In control	Positive alert	Positive alarm
Reference	5 (2.4)	12 (5.6)	186 (86.5)	12 (5.6)	0 (0)
Complete case	7 (3.0)	8 (3.7)	193 (90.0)	7 (3.0)	0 (0)
Single imputation	6 (2.8)	13 (6.0)	183 (85.6)	13 (6.0)	0 (0)
Multiple imputation	5 (2.4)	13 (6.0)	184 (85.6)	13 (6.0)	0 (0)
Normal/mean imputation	4 (2.0)	14 (6.5)	183 (85.0)	14 (6.5)	0 (0)

FIGURE 12.

Funnel plots for alternative approaches to handling missing data in the simulation study for application of the new ICNARC model. (a) Complete case; (b) multiple imputation; (c) single imputation; and (d) normal/mean imputation.

Evaluation with real missing values

Table 33 shows the results when the different approaches (multiple imputation, single imputation, normal/mean imputation and complete case) were applied to the full data set, including missing values. Using the normal/mean imputation strategy, the model showed a c-index of 0.892 (95% CI 0.890 to 0.895) and a calibration slope of 0.980 with a calibration intercept of –0.021. For the single and multiple imputation strategies, the c-index was slightly higher (0.894, 95% CI 0.892 to 0.896, and 0.895, 95% CI 0.892 to 0.898 respectively) with a similar calibration slope (0.979 and 0.988 respectively) and small difference in calibration intercept (−0.061 and −0.051 respectively). Under our ‘real’ scenario, the imputation of missing predictors prevented a drop of 0.01 in the c-index compared to a complete case approach (0.885, 95% CI 0.882 to 0.888). This 0.01 drop in the c-index may be considered to be a relevant decrease in the accuracy of predicted risks of individual patients because the c-index is relatively insensitive.

TABLE 33 - Discrimination and calibration of different approaches to handling real missing data in the application of the new ICNARC model

Approach to missing data	c-index (95% CI)	Brier score	Cox calibration regression		Predicted mortality (%)	Observed mortality (%)	SMR (95% CI)
			β	α
Complete case (n = 71,756)	0.885 (0.882 to 0.887)	0.107	0.949	–0.068	22.94	22.63	0.99 (0.97 to 1.00)
Single imputation (n = 90,107)	0.894 (0.891 to 0.896)	0.100	0.979	–0.061	21.92	21.48	0.98 (0.97 to 0.99)
Multiple imputation (n = 90,107)	0.895 (0.892 to 0.897)	0.100	0.988	–0.051	21.89		0.98 (0.97 to 0.99)
Normal/mean imputation (n = 90,107)	0.893 (0.890 to 0.895)	0.101	0.980	–0.021	21.52		1.00 (0.99 to 1.01)

Observed acute hospital mortality was 22.6% for the complete case admissions, in contrast to the predicted mortality of 22.9%, giving an SMR of 0.99 (95% CI 0.97 to 1.00). With the imputation approaches (single, multiple and normal/mean imputation), observed mortality was 21.48%, with predicted mortality 21.92%, 21.89% and 21.52%, respectively, for SMRs of 0.98 (95% CI 0.97 to 0.99), 0.98 (95% CI 0.97 to 0.99) and 1.00 (95% CI 0.99 to 1.01 respectively).

Table 34 and Figure 13 show the critical care unit performance relative to control limits for the complete case and imputation approaches in order to evaluate their impact on benchmarking when applied to the ‘real’ scenario. The proportions of critical care units with SMRs below 1 according to 95% and 99.8% control limits were similar for all approaches, including complete case. In contrast to what the simulation results showed, the proportion of critical care units with positive alert was higher based on complete case and normal/mean imputation and lower after single and multiple imputation approaches. Most of the shifts were from positive alert (based on complete case) to in control after imputation, although three critical care units remained as positive alert after normal/mean imputation. No positive alarms were found in the data.

TABLE 34 - Critical care unit performance status for different approaches to handling real missing data in the application of the new ICNARC model

Approach to missing data	Location of critical care units relative to funnel plot lines, n (%)
	Negative alarm	Negative alert	In control	Positive alert	Positive alarm
Complete case (n = 71,756)	5 (2.3)	12 (5.6)	186 (86.5)	12 (5.6)	0 (0)
Single imputation (n = 90,107)	6 (2.8)	11 (5.1)	193 (90)	6 (2.8)	0 (0)
Multiple imputation (n = 90,107)	5 (2.3)	12 (5.6)	192 (89.3)	6 (2.8)	0 (0)
Normal/mean imputation (n = 90,107)	5 (2.3)	11 (5.1)	189 (87.9)	10 (4.7)	0 (0)

FIGURE 13.

Funnel plots for alternative approaches to handling real missing data in the application of the new ICNARC model. (a) Complete case; (b) multiple imputation; (c) single imputation; and (d) normal/mean imputation.

Conclusions

The similarity of multiple imputation, single imputation and normal/mean imputation on the global results suggests a low impact of missing data on the overall outcomes. The differences in the means of the predicted probabilities to the reference were also minimal so will not affect the overall SMR in a large data set. The difference in performance between multiple imputation and single imputation suggests that repetition is an important component if the approach is to use imputed values; we therefore reject any further consideration of single imputation.

The multiple imputation approach was usually the most accurate method of imputation. However, imputation of normal/mean values demonstrated comparable accuracy in addressing isolated missing physiological predictors, except when imputing sedated/paralysed/GCS score, for which multiple imputation was the best approach. Regression methods for imputing isolated missing values were not the best approach and actually normal/mean imputation performed better in that situation.

When the number of missing predictors was high, the quality of the imputed values was not as good and the imputed values required careful evaluation. If the missing values were in strong predictors, multiple imputation resulted in a c-index that was close to the reference value, whereas all other methods, including a specific submodel, led to an underestimated c-index. In theory, using submodels derived from the development data set that contain only the available predictors could give better discrimination than imputing missing values; however, this is less likely if a strong predictors is missing, because any submodel without the predictor loses discriminative ability. In addition, if those risk-increasing predictors are ignored, then all the predicted risks are lower.

Differences between imputation methods depend not only on the number of missing values but also on the number of predictors missing and their importance in the model. Apparently the discriminative ability of the model was based on one particularly strong predictor (sedated/paralysed/GCS score) and the physiology and that is why when we have fewer missing values in these factors we recover information and improve the accuracy.

The lack of consistency of the approaches under different situations is an important finding, suggesting that there is no standard missing-value mechanism and thus no simple method for handing missing values when a model is applied.

For missing data at the point of model application, we could conclude that multiple imputation is the best approach for addressing patterns of missing predictors (i.e. when missing values occur simultaneously in multiple predictors) and normal/mean imputation for isolated missing predictors. Regression-based approaches showed the worst results and submodels lost discriminative ability and accuracy because of the importance of the missing factors in the risk prediction model. In order to be able to apply the model to a new single admission (or small number of admissions) with missing data, in which case multiple imputation is not feasible, normal/mean imputation is likely to be the best approach.

In a simulation study, we found that complete case analysis led to a substantially lower number of providers performing statistically above or below the national average. This may indicate that complete case analysis fails to detect statistically significant outliers (negative and positive alerts), increasing the proportion identified as ‘in control’.

Missing physiological values may lead to underestimation of predicted mortality; therefore, the number and type of missing variables should be taken into consideration when assessing the performance of a critical care unit. However, when represented using funnel plots, the impact of missing data was no longer evident. In a large data set, the impact of missing data could be minimal and not affect the SMR.

Inclusion and exclusion criteria

For the NCAA, data are collected for all individuals (excluding neonates) receiving chest compressions and/or defibrillation and attended by a hospital-based resuscitation team (or equivalent) in response to a 2222 call (2222 is the telephone number used to summon a resuscitation team in UK NHS hospitals). 16

For development of the risk prediction models, data were extracted for all individuals meeting the scope of the NCAA with a date of 2222 call between 1 April 2011 and 30 September 2012. Data for individual hospitals were included if the hospital had commenced participation in the NCAA prior to April 2012 and had validated data for at least 6 months. Individual team visit records meeting any of the following criteria were considered ineligible for inclusion in the risk prediction models: the arrest occurred before the patient’s arrival at hospital (even if the patient was subsequently attended by a hospital-based resuscitation team, usually in the emergency department, and, therefore, met the scope of the NCAA); second and subsequent visits to the same patient during the same hospital stay; a ‘do not attempt cardiopulmonary resuscitation’ (DNACPR) was already documented in the patient’s notes. The following exclusion criteria were applied to individual team visit records: the patient’s last known status was still in hospital; either of the outcomes of ROSC > 20 minutes or hospital survival was missing; and data for the candidate predictors were missing.

For validation of the risk prediction models, data were extracted for all individuals in hospitals that contributed data to the development data set with a date of 2222 call between 1 October 2012 and 31 March 2013, and for all individuals in hospitals that commenced participation in the NCAA between April and September 2012 (and therefore did not contribute data to the development data set), with a date of 2222 call between 1 April 2012 and 31 March 2013. The same eligibility and exclusion criteria were applied at the individual team visit level as for the development data set.

Outcomes and candidate predictors

Risk prediction models were developed for two outcomes: ROSC > 20 minutes and hospital survival. Patients were followed up to discharge from the original acute hospital and any patients transferred to another acute hospital were reported as hospital survivors.

A list of candidate predictors was established from the data set developed and collected for the NCAA. A valid predictor was considered to be any variable collected either prior to or at the time of the arrival of the hospital-based resuscitation team and not related to the quality of care (e.g. whether or not the appropriate resuscitation was delivered). If factors related to the quality of care were included within the risk prediction model, then the expected number of events would be adjusted to account for these factors. Consequently, a poorly performing provider would not be identified as an outlier and these discrepancies in the quality of care would not be recognised.

The full list of candidate predictors is presented in Table 35. Location of arrest was not considered to be a predictor for patients with a reason for admission to/attendance at/visit to hospital of ‘staff’ or ‘visitor’. Prior to any modelling, candidate predictors were examined for data completeness and distribution. Where categories with very few patients were identified, combining with other categories eliminated these. Multicollinearity between candidate predictors was assessed using variance inflation factors.

Age was modelled as a continuous, non-linear relationship using restricted cubic splines with five knots. All other candidate predictors were modelled as categorical variables. After examining plots of the distribution and the association with the outcomes, the continuous predictor length of stay in hospital prior to 2222 call was categorised as 0 days (i.e. cardiac arrest on the same calendar day as admission to/attendance at/visit to hospital), 1 day, 2–7 days, 8–30 days and > 30 days.

Model development

An initial, full model for each outcome was fitted, including all candidate predictors, using multilevel logistic regression with random effects of hospital. The basis for using multilevel models is that, as with the majority of health outcomes, there is ‘clustering’ at the level of health-care providers; that is, outcomes for patients within the same hospital will be, on average, more similar than outcomes for patients in different hospitals. If clustering is ignored, then the resulting model estimates will have standard errors that are too small, contributing to the potential for misleading conclusions.

These models were then simplified in three stages: first, by testing for non-linearity in the relationship for age; second, by testing for differences between pre-specified combinations of categories of predictors to reduce the numbers of categories; and, third, by stepwise reduction of the models to reduce the number of predictors. Combining categories of predictors was conducted in such a way as to ensure the same categories were used in the models for both outcomes. This was achieved by combining categories if the difference in outcome between the categories (adjusted for all other predictors) was non-significant (p > 0.1) for both outcomes. The combinations considered were:

1. Prior length of stay: adjacent categories.

2. Location of arrest:

   1. adjacent categories from: emergency department; emergency admissions unit; ward, obstetric area or other inpatient location; intermediate care area; coronary care unit; high-dependency unit or paediatric high-dependency unit; intensive care unit, combined high-dependency/intensive care unit or paediatric intensive care unit

   2. any combination of categories from: specialist treatment area; imaging department; cardiac catheter laboratory; and, if no difference was found between any of the three previous categories, theatre and recovery

   3. categories for clinic and non-clinical area.

Stepwise reduction was conducted separately for the two outcomes, allowing the models for ROSC > 20 minutes and hospital survival to include different combinations of predictors. At each step, the least significant predictor was removed and the reduced model assessed for discrimination (c-index22), calibration (Hosmer–Lemeshow test24), accuracy (Brier score28 and Shapiro’s R29) and model fit (AIC99). Stepwise reduction was continued until all predictors had been removed and a model was selected to balance simplicity against model performance.

Finally, the models were further enhanced by the consideration of interactions between predictors. The interactions considered were pre-specified following presentation of an initial model with no interaction terms and a request for input from the NCAA Steering Group, representatives from hospitals participating in the NCAA, and the Expert Advisory Group. The interactions considered were:

• age with sex

• age with reason for attendance

• age with presenting rhythm (non-shockable rhythms compared with shockable or unknown)

• location of arrest with presenting rhythm.

Interaction terms were added to the full model and retained if significant at p < 0.01. For the interaction of location of arrest with presenting rhythm, in order to reduce the potentially large number of interaction terms, combining interaction terms for similar groups of categories of both presenting rhythm (e.g. all shockable arrests, all non-shockable arrests) and location of arrest (e.g. emergency department and emergency admissions unit, emergency admissions unit and ward, coronary care unit and cardiac catheter laboratory) was considered.

Comparisons of models (for testing linearity, combining categories, stepwise reduction and adding interactions) were performed with likelihood ratio tests.

Model validation

The resulting models were validated for discrimination, calibration and accuracy in (1) the development data set; (2) the full validation data set; and (3) the validation data from hospitals that commenced participation in the NCAA from April 2012 onwards not included in the development data set (providing true external validation in a smaller sample of hospitals). To reduce overfitting, model estimates were shrunk using the uniform (heuristic) shrinkage method of Van Houwelingen and Le Cessie. 100

Discrimination was assessed by the c-index. 22 Calibration was assessed graphically and tested using the Hosmer–Lemeshow test for perfect calibration in 10 equal-sized groups by predicted probability of survival. 24 As the Hosmer–Lemeshow test does not provide a measure of the degree of miscalibration and is very sensitive to sample size,25,26 calibration was also assessed using Cox’s calibration regression. 27 Accuracy was assessed by the Brier score28 and Shapiro’s R,29 and the associated approximate R² statistics. 30 Measures of model performance were calculated using the marginal predicted probabilities from the risk prediction model (i.e. without taking into account hospital-level effects) to represent the predicted probability of survival for a patient with the given characteristics in an ‘average’ hospital.

The final risk prediction models were refitted to all data (development and validation data sets combined) to maximise precision and generalisability, with shrinkage applied to reported coefficients.

Statistical analyses were performed using Stata/SE, version 10.1 (StataCorp LP, College Station, TX, USA).

Available data

Between 1 April 2011 and 31 March 2013, 148 hospitals participated in the NCAA. During this time there were a total of 28,987 resuscitation team visits following 2222 calls for cardiac arrest reported to the NCAA. After excluding data that were still undergoing central validation (at the level of calendar months within hospitals) and hospitals with less than 6 months’ data, 27,998 team visits in 143 hospitals were included (Figure 14). After removing records that were ineligible for risk predictions (pre-hospital arrests, second and subsequent visits to the same patient and patients with a documented DNACPR decision) and those excluded for missing data, a total of 22,479 team visits in 143 hospitals were included: 14,688 (65.3%) in the development data set and 7791 (34.7%) in the validation data set. Rates of missing data were very low, with only 0.1% of patients excluded from the development data set (0.1% from the validation data set) because of missing predictor variables and 0.1% (0.8% from the validation data set) because of missing outcomes, and it was therefore not necessary to consider more complex statistical methods for handling missing data. The breakdown of exclusions in the development data set, the validation data set and the external validation data set (the subset of the validation data set from hospitals not included in the development data set) are shown in Table 36, and characteristics and outcomes are summarised in Table 37.

FIGURE 14.

Participation of hospitals in the NCAA, 1 April 2011 to 31 March 2013.

Model development

Prior to modelling, the following categories of predictors were combined to remove small categories. For reason for admission to/attendance at/visit to hospital, the categories of staff and visitor were combined. For location of arrest, the following categories were combined: ward, obstetric area and other inpatient location (noting that obstetric patients are distinguished by the separate reason for attendance field); high-dependency unit and paediatric high-dependency unit; and intensive care unit or combined high-dependency/intensive care unit and paediatric intensive care unit (noting that paediatric patients are distinguished by age).

The initial, full model, including the main effects of all candidate predictors (see Table 35), had a c-index of 0.727 for a ROSC > 20 minutes and 0.804 for hospital survival (Table 38). There was no evidence of multicollinearity (all variance inflation factors less than 2). Age was significantly non-linear in both models (p < 0.001 for a ROSC > 20 minutes and p = 0.007 for hospital survival).

After following the pre-specified process for combining categories of predictors, the following categories were combined:

1. Prior length of stay: 8–30 days with > 30 days (p = 0.21 for a ROSC > 20 minutes and p = 0.73 for hospital survival).

2. Location of arrest:

   1. ward, obstetric area or other inpatient location with intermediate care area (p = 0.56, p = 0.47 respectively)

   2. high-dependency unit or paediatric high-dependency unit with intensive care unit, combined high-dependency/intensive care unit or paediatric high-dependency unit (p = 0.16, p = 0.40 respectively)

   3. specialist treatment area with imaging department (p = 0.82, p = 0.34 respectively)

   4. clinic with non-clinical area (p = 0.69, p = 0.39 respectively).

Combining categories had a minimal effect on the measures of model performance and resulted in an improvement (decrease) in the AIC (see Table 38).

The stepwise reduction of the models is shown in Table 38. The predictor ‘patient deteriorating (not yet arrested) at team arrival’ was removed from both models and sex was removed from the model for hospital survival. All other predictors were highly statistically significant.

After testing the pre-specified interactions, a significant interaction (p < 0.001) was found between location of arrest and presenting rhythm in both models; therefore, alternative categorisations for interactions between location of arrest and presenting rhythm were considered. All other interaction terms were non-significant.

The non-linear relationships between age and outcome are illustrated in Figure 15. For ROSC > 20 minutes, the relationship with age was flat up to around age 60 years, with a rapid decrease in the odds of ROSC > 20 minutes at older ages. Hospital survival decreased across the full age range, although this relationship was steeper at older ages.

FIGURE 15.

Relationship between age and ROSC > 20 minutes (a) and hospital survival following in-hospital cardiac arrest; and (b) odds ratios presented relative to age 70 years.

Model validation

The results of the model validation, based on models fitted in the development data set, are shown in Table 39. Discrimination and accuracy were better for hospital survival (c-index ≈ 0.81, R² = 0.21–0.24) than for ROSC > 20 minutes (c-index ≈ 0.73, R² = 0.11–0.17). Calibration was generally good, supported visually by calibration plots (Figure 16), although there was some evidence of worse calibration for ROSC > 20 minutes in the validation data set. Model performance was generally well preserved in the validation data sets compared with the development data set, particularly for hospital survival. Model accuracy was also compared across age groups (Figure 17). Although there was some variation in outcomes (consistent with chance) in the age groups with smaller sample sizes, overall the model fit was good across all age groups. Interactions between age and other predictors were considered but were found to be unnecessary.

FIGURE 16.

Calibration plots for ROSC > 20 minutes (a–c) and hospital survival (d–f) following in-hospital cardiac arrest. Observed outcome (with 95% CI) plotted against predicted outcome in 10 equal-sized groups, based on models using coefficients fitted in the development data set.

FIGURE 17.

Numbers of in-hospital cardiac arrests and observed and predicted outcomes by age in 5-year bands. (a) Development data set; and (b) validation data set.

The final models for ROSC > 20 minutes and hospital survival, refitted to the full data set, are shown in Tables 40 and 41 respectively. The shrinkage factors were 0.964 and 0.970, respectively, indicating very little overfitting.

Recommendation 1: further research should be conducted by linking with death registrations to evaluate mortality at fixed time points and using time to event analyses

To date, the main outcome for national clinical audits, including the CMP and NCAA, has been mortality at acute hospital discharge, an event-based outcome. Time-based outcomes, for example, mortality at 30 or 90 days following admission or duration of survival, would be less prone to bias arising from variation in provision of community health- and social care services, which may impact the timing of patients being discharged from acute hospital. Research from the Netherlands has suggested that comparison of risk-adjusted mortality across critical care units using mortality at 30 or 90 days, rather than at hospital discharge, results in less heterogeneity. 109 Acute hospital mortality has predominantly been selected because of its convenience to record and collect, as follow-up of patients beyond acute hospital discharge has not been seen as practicable. However, with the increased availability of electronic data sets and the establishment of the NHS Number as a national unique identifier, the vast majority of patients admitted to UK critical care units will now be able to be followed up for mortality following discharge from acute hospital by using data linkage with death registrations maintained by the Office for National Statistics. Furthermore, to enable reporting in a useful timeframe for quality improvement, the main outcomes for national clinical audit are necessarily relatively short term in nature. However, recovery from critical illness can be a slow process, with studies reporting substantial ongoing burden of mortality several years after discharge from hospital. 110,111 Data linkage with death registrations would also permit follow-up of longer-term mortality, enabling us to better understand the time course of recovery from critical illness and which risk factors impact on longer-term mortality.

Recommendation 2: further research in this field should make better use of data linkage across national clinical audits

Although mortality is clearly an important and patient-centred outcome, the impact and consequences of critical care go beyond mortality. As data on longer-term, health-related quality of life for survivors of critical care are currently not routinely collected, national clinical audits of chronic health conditions provide an ideal opportunity to better understand the impact and consequences of critical illness on these specific chronic health conditions, gaining some insight into the wider impact of critical care on patients’ subsequent health status. Critical illness-related hyperglycaemia, for example, has previously been linked with subsequent development of type 2 diabetes,112,113 and acute kidney injury, common among critically ill patients, has been strongly linked with subsequent end-stage renal disease. 114 Both these conditions have well-established national clinical audits (the National Diabetes Audit and the UK Renal Registry), and data linkage across these audits would permit a better understanding of the specific consequences of critical illness. In addition, the risk prediction models developed in this report were limited to the available predictors within the CMP and NCAA data sets which, in turn, are limited by what it is feasible to expect providers to routinely collect for the purpose of national clinical audit. Data linkage across national clinical audits would also enhance the available pool of candidate predictors on which risk prediction models could be developed, potentially allowing for improved prediction with no additional data collection burden. For example, existing literature on predictors of mortality for post cardiac surgery patients suggests that outcomes are best predicted by a combination of pre-, intra- and postoperative risk factors. 51 The CMP data set is a reliable source of postoperative information, but has limited pre- or intraoperative data. Data linkage with the National Adult Cardiac Surgery Audit would greatly increase the available pre- and intraoperative data to improve risk prediction among this patient group.

Recommendation 3: further research in this field should make better use of other routinely collected data sets

Making better use of other routinely collected data sets, such as Hospital Episode Statistics (HES), would also potentially expand the available predictors and outcomes. For example, the recently developed US risk prediction model for hospital survival following in-hospital cardiac arrest includes information on pre-arrest comorbidities and interventions. Collection of comorbidities was considered for the NCAA, but it was decided that this would add a considerable data collection burden. Data linkage with HES would permit evaluation of comorbidities and interventions through use of diagnostic and procedural codes. If these factors improve the fit of the model, then options to either undertake this data linkage routinely or to incorporate specific additional fields into the NCAA data set could be explored. In terms of expanding available outcomes, survivors of critical care experience significant morbidity with substantial resultant health-care resource use and costs. 115 Data linkage with HES would enable the cost of subsequent hospitalisations, and its association with severity and/or duration of critical illness and other risk factors, to be estimated

Recommendation 4: future research should consider the necessity for specific data collection to support national clinical audit compared with benchmarking providers using routinely collected data alone

National clinical audits rely on a separate, specific collection of clinical data in addition to the routine data, such as HES, that are collected in the course of a patient’s journey through the health-care system. It has generally been held that these detailed clinical data are essential for reliable benchmarking of specific clinical services. However, national clinical audit is relatively expensive and consideration should therefore be given to whether or not clinical services, such as critical care, can be benchmarked using routinely collected data alone. Considerations would include data capture, that is, whether or not the same patients can be identified from both national clinical audit data and routine data (both in terms of the completeness of data capture for national clinical audit and also whether the users of a specific service can reliably be identified from routine data), and also the performance of risk prediction models based on routinely collected data, compared with those based on detailed clinical data, for identifying providers with potentially outlying performance.