2.1 Study design and participants

This study was performed using a prospectively collected registry of ED patients who were diagnosed with sepsis and septic shock according to the Sepsis-3 definitions in a tertiary care teaching hospital in Korea. From January 2016 to February 2020, adults (aged ≥18 years) who visited the ED with an actual or suspected infection and an increase in SOFA score of ≥2 points were enrolled by ED physicians on duty. If the patients had baseline SOFA scores, we used the standard of an increase in the SOFA score of at least two points. If the patients had no previous SOFA score, we reviewed their medical records with laboratory results and determined the association between the present infection and the SOFA score. The exclusion criteria were: (a) age <18 years, (b) death within 12 hours of ED presentation, (c) unknown outcome (30-day mortality), and (d) ED visit for trauma care. Finally, two infectious disease experts carefully reviewed all the patients. Clinical variables were collected from the patients and different ML algorithms were used to predict 30-day mortality. The prognostic values of the different ML algorithms were compared with that of the logistic regression models and established clinical scoring systems. This study followed the transparent reporting of a multivariable prediction model for individual prognosis or diagnosis (TRIPOD) reporting guidelines [16].

2.2 Dataset and definitions

Sepsis is as a life-threatening organ dysfunction caused by a dysregulated host response to infection [4]. Septic shock is a subset of sepsis where profound circulatory, cellular, and metabolic abnormalities pose a greater risk of mortality compared to sepsis alone [4, 17]. The qSOFA score uses three criteria: low blood pressure (systolic blood pressure: ≤100 mmHg), high respiratory rate (≥22 breaths/min), and altered mental status (Glasgow coma score: <15), assigning one point for each criterion with the final score ranging from zero to three points. A positive qSOFA score was defined as the presence of ≥2 qSOFA points near the onset of infection. The diagnostic criterion for sepsis includes an increase in the SOFA score by ≥2 points due to the present infection. The SOFA score is used to track a person’s status to determine the extent of organ function or the rate of failure [18, 19, 20]. The score is based on six different scores, each for the respiratory, cardiovascular, hepatic, coagulation, renal, and neurological systems. When the physiological parameters do not match any row, zero points are given. The criteria for septic shock include vasopressor requirements to maintain a mean arterial pressure of 65 mmHg and serum lactate level >2 mmoL/L despite adequate fluid resuscitation. The National Early Warning Score (NEWS) (2012) and NEWS2 (2017) are tools developed by the Royal College of Physicians to improve the detection and response to clinical deterioration in adult patients and are critical to patient safety and improving patient outcomes [21]. The modified early warning score (MEWS) is a composite score used by physicians to determine the severity of illness in various clinical situations [22]. It evaluates the risk of mortality based on vital signs (systolic blood pressure, heart rate, respiration rate, and body temperature), saturation of percutaneous oxygen (SpO₂), and level of consciousness. Determining a MEWS involves assigning a number between zero and three to each of the six items. A total score equal to or more than five points of MEWS is associated with an increased mortality. In the present study, we defined follow-up (F/U) lactate levels within 12 hours as the maximum values of lactate levels measured within 12 hours of ED presentation, except initial lactate levels.

2.3 Data splitting and preprocessing

Fig. 1 shows this study’s ML flowchart. Variables with more than 40% missing data were excluded [23], and 25% of the data were randomly separated with stratifications of sepsis-related death. The hold-out data were used only in the final model evaluation as a test set, and the remaining 75% of the data were used as a training set with a leave-one-out cross-validation strategy. Multivariate imputation by chained equations [24] and isolation forests [25] were used for imputation and outlier detections, respectively. Min-max scaling and one-hot encoding were applied, and Supplementary material provides the details of the preprocessing procedures. Included variables for model development with information of their missing rate and variable type are presented in Supplementary Table 1.

Fig. 1.Machine learning flowchart. MICE: multiple imputation by chained equations; XGBoost: extreme gradient boosting; SVM: support vector machine; LightGBM: light gradient boosting machine; MLP: multilayer perceptron.

2.4 Feature selection and feature importance analyses

Recursive feature elimination [26] was utilized to obtain the best feature set maximizing the model performance of the area under the receiver operating characteristic curve (AUROC) in cross-validation. Recursive feature elimination works with the following iterative procedures: (1) training a classifier, (2) computing the ranking criterion of all features, and (3) removing the feature with the lowest rank. A light gradient boosting machine (LightGBM) [28], which can handle categorical variables without data conversion, was used for the feature importance analyses. Contributions of each feature to the model prediction were measured with Shapley Additive exPlanations (SHAP) values, where positive and negative values indicated a positive and negative effect on the prediction score, respectively [27]. The relative importance of the features was evaluated and ranked using the mean absolute SHAP value. Hierarchical clustering was performed to minimize the effects of multicollinearity, causing the underestimation of the relative importance of the features [29]. Supplementary material describes the details.

2.5 Model construction

A conventional statistical model-logistic regression—was built as a baseline comparator. Four popular and promising models were constructed: extreme gradient boosting (XGBoost) [30], support vector machine (SVM) [31], LightGBM, and multilayer perceptron (MLP) [32]. XGBoost (version 1.3.2, The XGBoost Contributors, New York, NY, USA) and LightGBM (version 3.2.0, Microsoft, Wilmington, Del, USA) are gradient-boosted tree-based ensemble models using a depth-wise algorithm and a leaf-wise growth algorithm, respectively. SVM is a model for classification, regression, and outlier detection using an optimal hyperplane in a multi-dimensional space. MLP is a feed-forward neural network model with a basic architecture consisted of fully connected layers. Bayesian optimization [33] was used for hyperparameter settings to maximize the AUROC in cross-validation. Supplementary material provides the details of the constructed models. All processes including development of SVM and MLP models were implemented in Python version 3.8.2 (Python Software Foundation, Wilmington, Del, USA) with scikit-lean (version 0.24.1) and tensorflow (version 2.0.0) libraries.

2.6 Statistical analysis

Descriptive statistics were represented as the number (percentage), mean (SD), or median (Q1 and Q3). We used the Shapiro-Wilk test for normality and Levene’s test for homoscedasticity. The Chi-squared test, independent t-test, or Mann-Whitney U test were used for comparison analysis.

The AUROC was chosen as the primary evaluation metric with a threshold of 0.50. The AUROC is classified into excellent (0.9–1.0), good (0.8–0.9), fair (0.7–0.8), poor (0.6–0.7), and fail (0.5–0.6). The AUROC with its 95% confidence interval (CI) was calculated and compared between constructed models using Delong’s method [34]. p-values were adjusted using the Benjamini-Hochberg method for multiple comparisons.

During cross validation for developing the models, each fold was calibrated to validation data after fitting at training data using isotonic regression. The isotonic regression fits a non-decreasing real function to one dimensional data. Prediction scores from the calibrated models allow reflection for predictive probability. Error after the post-hoc calibration was evaluated using calibration curves and the Brier score which was the mean squared error of the predicted probability [35]. The AUPRC, balanced accuracy, and F1 score were calculated as the secondary outcome metrics for model performance. The significance level was set at p < 0.05.

3.1 Baseline characteristics

During the study period, 1075 consecutive patients with sepsis were screened. Fig. 2 shows a flowchart of the study population. Patients aged <18 years (n = 65), death within 12 hours of ED presentation (n = 41), unknown 30-day mortality (n = 129), and ED visits for trauma care (n = 30) were excluded. After excluding 265 patients, a final total of 810 patients were finally enrolled in the analysis.

Fig. 2.Flowchart of the study population.

Table 1 shows the baseline characteristics of the patients for 30-day mortality outcomes (Supplementary Tables 2 and 3 for 7-day and 14-day mortalities, respectively). The median age of the enrolled patients was 75 years (interquartile range (IQR): 65 to 82 years), and 337 (41.6%) were women. Among the patients, 259 (32.0%) were non-survivors. There were significant differences in age, vital signs except for heart rate, comorbidities (diabetes mellitus, malignancy, chronic liver disease, and cerebrovascular disease), source of infection (respiratory and genitourinary infections), laboratory findings (platelet, bilirubin, lactate, and arterial blood gas analysis), septic shock, and clinical severity scores (NEWS, NEWS2, MEWS, and SOFA scores) between non-survivors and survivors. The missing data rates were less than 3.1% except for three variables: F/U lactate levels within 12 hours (missing rate, 35.1%), lactate clearance (missing rate, 20.2%), and procalcitonin (missing rate, 16.7%) (Supplementary Table 1).

3.2 Performance of the models and scoring systems

LightGBM achieved the highest AUROC values in prediction for 7-day (0.891 (0.841–0.941)), 14-day (0.893 (0.844–0.941)), and 30-day (0.871 (0.823–0.919)) mortalities among the ML models (Table 2 and Fig. 3(A)), but there were no significant differences between LightGBM and MLP in predictions for 7-day, 14-day, and 30-day mortalities. All the ML models except XGBoost in predictions for 7-day and 14-day mortalities exhibited significantly higher AUROC values compared to the logistic regression.

SOFA achieved the highest AUROC values among the scoring systems in predicting 7-day (0.680 (0.594–0.766)), 14-day (0.647 (0.565–0.730)), and 30-day (0.658 (0.579–0.736)) mortalities. MEWS had no discrimination ability in predicting 30-day mortality, of which the confidence interval of AUROC spanned 0.50. All the ML models and logistic regression had significantly higher AUROC values compared to SOFA scores.

ML models exhibited higher performance in the AUPRC (Table 2 and Fig. 3(B)) and Brier score (Fig. 4 for 30-day mortality and Supplementary Figs. 1,2 for 7-day and 14-day mortalities) compared to the scoring systems and logistic regression. In particular, LightGBM had the highest AUPRC values in prediction for 7-day (0.702 (0.547–0.847)) and 14-day (0.786 (0.670–0.884)) mortalities, while MLP had the highest AUPRC values in predicting 30-day mortality (0.774 (0.662–0.865)). MLP in the prediction for 7-day (0.789) and 14-day (0.774) mortalities and LightGBM in predicting 30-day mortality (0.762) exhibited the highest balanced accuracy. Supplementary Table 4 presents the cross-validation performance.

Fig. 3.Test set performance in prediction for 30-day mortality. A, Solid lines and shades representing receiver operating characteristics curves and its 95% confidence intervals. An asterisk (*) indicates a significant difference (p < 0.05) in comparison with logistic regression. B, Solid lines and shades representing precision-recall curves and its 95% confidence intervals. Only the confidence intervals of the baseline model (logistic regression (LogReg)) are represented with a polka dot pattern in both plots. AUC, area under the receiver operating characteristic curve; CI, confidence interval; XGBoost, extreme gradient boosting; SVM, support vector machine; LightGBM, light gradient boosting machine; MLP, multilayer perceptron.

Fig. 4.Post-hoc calibration of models in the prediction for 30-day mortality. Orange and blue lines represent the calibration curves of raw and calibrated models (isotonic regression), respectively. The legend displays the mean (SD) value of brier scores for model prediction, while solid lines and shades represent the mean and ± SD of the calibration curves. (A) Calibration curves of logistic regression; (B) Calibration curves of XGBoost; (C) Calibration curves of support vector machine; (D) Calibration curves of light gradient boosting machine; (E) Calibration curves of multilayer perceptron. XGBoost, extreme gradient boosting; SVM, support vector machine; LightGBM, light gradient boosting machine; MLP, multilayer perceptron.

3.3 Importance of variables

Thirty-three variables were selected as the best feature set by the recursive feature elimination in predicting 30-day mortality, and all selected variables in the predictions for 7-day and 14-day mortality were included. Fig. 5 shows the overall contribution of the variables (left bar plot) and the impacts of individual values for the model prediction (the right violin plot). Fig. 6 shows the partial SHAP dependence plots for the top six important variables, while Supplementary Fig. 3 presents the plots for the other variables.

The most important variable was septic shock, which had a positive effect on the model prediction, followed by F/U lactate levels within 12 hours followed by the initial lactate levels. When the F/U lactate level within 12 hours was less than 3.7 mmoL/L, it consistently exhibited a negative effect on the model prediction. However, it reversed at 3.7 mmoL/L, and SHAP values increased sharply with positive effects for the model prediction. Initial lactate levels exhibited a similar tendency with a cut-off of 5.3 mmoL/L. The following three variables were considered: malignancy, age, and arterial oxygen saturation (SaO₂). Malignancy had a positive effect on model prediction. Age exhibited a positive correlation with model prediction, resulting in the shape of the sigmoid curve and the cut-off value of 75.1 years. SaO₂ demonstrated a negative correlation with model prediction and a cut-off value of 87.5%.

Fig. 5.Importance of the selected variables. The left bar plot shows the overall contribution of the variables to model prediction. The right violin plot shows the impacts of individual values of the variables and directionality for model prediction. The red color represents the large value in continuous variables or the affirmative response in categorical variables. Supplementary Table 3 presents the variable types. F/U, follow up; HCO₃, bicarbonate; pH, potential of hydrogen; SaO₂, Arterial oxygen saturation; SpO₂, Saturation of percutaneous oxygen; PaO₂, partial pressure of oxygen; PaCO₂, partial pressure of carbon dioxide; SHAP, Shapley Additive exPlanations; WBC, white blood cell.

Fig. 6.Partial SHAP dependence plots for the representative top six important variables. (A), (D) Categorical variables are plotted with a scatter plot and a box plot with whiskers of 1.5 times the interquartile ranges. (B), (C), (E), (F) Continuous variables are plotted with a scatter plot and a regression line represented with the orange line of mean and shade of SD. A red diamond represents a cut-off value of the variable. Histograms on the right and top of each plot are distributions of the SHAP and values of variables. SHAP, Shapley Additive exPlanations; F/U, follow up; SaO₂, Arterial oxygen saturation.