A prediction model for 30-day deaths of cirrhotic patients in intensive care unit hospitalization.

ABSTRACT: The aim of this study was to establish a prediction model for 30-day deaths of cirrhotic patients in intensive care unit.A case-control study involving 1840 patients was conducted in the Medical Information Mart of the Intensive Care Database III version 1.4. The logistic regression with L1 regularization was used to screen out the variables. The 30-day in-hospital death was used as the dependent variable and the selected variables were used as the independent variable to build a random forest model. The performance of the model was validated by the internal validation.The variables screened by logistic regression analysis were the age, heart rate, respiratory rate, systolic blood pressure, diastolic blood pressure, Oxygen saturation, white blood cells, platelets, red cell distribution width, glucose, blood urea nitrogen, bicarbonate, total bilirubin, hematocrit, alanine transaminase, aspartate transaminase, bilirubin, Simplified Acute Physiology Score II and Sequential Organ Failure Assessment. The areas under the curve of the random forest model based on these variables was 0.908, and the performance of this model were internally validated with an areas under the curve of 0.801. The random forest model displayed that Simplified Acute Physiology Score, Sequential Organ Failure Assessment, blood urea nitrogen, total bilirubin and bilirubin were more important predictors for the 30-day death of cirrhotic patients in intensive care unit.A prediction model for death of cirrhotic patients was developed based on a random forest analysis, providing a tool to evaluate the patients with a high risk of 30-day in-hospital deaths to help clinician make preventive intervention to decrease the mortality.

Introduction

Hepatic cirrhosis is the terminal stage of liver disease, and can cause continuous damage to liver cells.[ It is worth noting that the onset and course of this disease are generally slow, which may conceal for 3 to 5 years or more than 10 years.[ In recent years, the morbidity and mortality of hepatic cirrhosis are both increasing, and has ranked the fourteenth of the most frequent cause of adult deaths worldwide.[ Studies in the UK and Sweden report that the annual incidence rate of hepatic cirrhosis is 15.3–132.6 per 100,000 people.[ A screening program in France shows that the prevalence of hepatic cirrhosis is 0.3%.[ However, due to the insidious onset, this disease is often asymptomatic at early stages; therefore, the actual prevalence may be higher than reported.[ Cirrhosis can lead to various fatal complications, such as hepatocellular carcinoma, hepatic encephalopathy, gastrointestinal hemorrhage, infections and so on.[ All of which may increase the mortality of patients and bring financial burdens to patients’ families and the society.[

Critically ill cirrhosis is a type of clinical critical illness with a high death rate and attracts increasing attention.[ Existing study has focused on patients with critically ill cirrhosis. A survey abroad has indicated the 30-day in-hospital mortality of patients with critically ill cirrhosis and bacterial ascites reaches 33%.[ The high death rate has highlighted the importance of early identifying the patients with high risk, which is helpful to earlier make treatments.

Our study aimed to establish a prediction model to predict the 30-day deaths of cirrhosis patients in intensive care unit (ICU) hospitalization using the data from Medical Information Mart of the Intensive Care Database III version 1.4. The prediction model is helpful to provide clinicians with an early prediction of hospital mortality and subsequent hints for adequate treatments of patients with high risk to decrease the 30-day in-hospital mortality.

Material and methods

Study population

The patients’ data were collected from MIMIC-IIIv1.4, which provided access for the public and was free of charge. The database contained the information of over 50,000 ICU patients who visited Beth Israel Deaconess Medical Centre (BIDMC, Boston, MA, USA) from 2001 to 2012.[ Considering that our data were accessed from MIMIC-III database, an openly available dataset, there was no need of ethic approval and informed consent.

Variable collection

All variables were recorded within 24 hours of ICU admission. For variables measured more than once, the result of the first measurement was included in the analysis.

Data of patients were collected, including age, gender, heart rate, respiratory rate, systolic blood pressure (SBP), diastolic blood pressure (DBP), Oxygen saturation (SpO2), white blood cells (WBC), platelets (PLT), red blood cell (RBC), red cell distribution width (RDW), levels of glucose, creatinine, blood urea nitrogen (BUN), bicarbonate, total bilirubin (TBIL), hematocrit, hemoglobin, alanine transaminase, aspartate transaminase (AST), bilirubin, Simplified Acute Physiology Score (SAPS) II, Sequential Organ Failure Assessment (SOFA), chronic obstructive pulmonary disease (COPD), heart failure, diabetes, and septicemia. Missing values were deleted and outliers were replaced with null values. We fitted and imputed the missing values based on random forest. Sensitivity analysis was performed using the baseline data of the training set before and after the imputation, and no significant differences were shown (Supplementary Digital Content Table S1). The outcome variables of this study were 30-day deaths of ICU patients with liver cirrhosis, and the start date of the record was the date that patients admitted to the hospital.

Development of the model

The prediction model was developed using random forest method. Total samples were split into training set (n = 1288) and testing set (n = 552) with the ratio of 7:3. The logistic regression with L1 regularization was conducted to select variables from the training set for the construction of the prediction model, which was subsequently validated by a 6-fold cross-validation. The selected variables were independent variables, and 30-day in-hospital death was dependent variable. The number of decision trees was 800, and the maximum depth of trees was 5. To determine the performance of the model, we calculated the accuracy, sensitivity, specificity, positive predictive value (PPV), and negative predictive value (NPV), and the area under the curve (AUC). The cut-off points were ascertained based on the Youden Index. The model performance was validated using the data of testing set. The learning curve was used to evaluate the stability of this model and the calibration curve was used to assess imitative effect. The feature importance computed by Gini importance was used to assess the important variables. The Gini importance was computed from the Random Forest structure, and the average over all trees in the forest was the measure of the feature importance. The model has been uploaded in GitHub (https://github.com/abcgut/data_model).

Statistical analysis

The statistical analyses were performed using SAS9.4 and Python 3.7 software. Categorical data were presented as the number of cases and the constituent ratio (N (%)), and Chi-Squared test or Fisher exact probability method was adopted for the comparisons. Normally distributed continuous variables were expressed as mean ± standard deviation (Mean ± SD), and t-test was applied for the comparisons between the survival group and the death group. Those quantitative data of skewed distribution were displayed as median and quartiles (M [Q1, Q3]), and the between-group comparisons were tested by Wilcoxon. All statistical tests were two-sided, and P < .05 was considered to be statistically significant.

Results

Baseline characteristics

A total of 1851 ICU patients with cirrhosis were selected from MIMIC database, with 11 patients missing the data of blood routine and liver function index. After processing the abnormal values and missing values, 1840 patients were finally included. These patients were randomly divided into the training set (n = 1288) and testing set (n = 552) (Fig. 1). As shown in Table 1, no significant differences were shown in the baseline characteristics of training set and testing set (P > .05).

Figure 1

The flowchart of research methodology.

Table 1

Baseline characteristics of training set and testing set.

Variables	Total	Training set	Testing set	Statistics	P
Age (year), Mean ± SD	58.70 ± 12.07	58.65 ± 12.01	58.81 ± 12.21	t = 0.250	.803
Gender, n (%)				χ2 = 1.502	.220
Female	601 (32.66)	432 (33.54)	169 (30.62)
Male	1239 (67.34)	856 (66.46)	383 (69.38)
Heart_rate (n/min), Mean ± SD	91.12 ± 19.08	91.15 ± 18.85	91.03 ± 19.64	t = −0.130	.898
Respiratory_rate(n/min), M (Q1, Q3)	18.00 (15.00, 22.00)	18.00 (14.00, 22.00)	18.00 (15.00, 22.00)	Z = 1.079	.281
SBP (mm Hg), Mean ± SD	119.49 ± 23.35	119.42 ± 23.28	119.66 ± 23.52	t = 0.210	.836
DBP (mm Hg), Mean ± SD	63.33 ± 16.10	63.25 ± 16.44	63.54 ± 15.28	t = 0.370	.715
SpO2 (%), Mean ± SD	97.31 ± 3.81	97.28 ± 4.09	97.39 ± 3.05	t = 0.610	.543
WBC (K/μL), M (Q1, Q3)	8.60 (5.70, 13.25)	8.69 (5.60, 13.25)	8.40 (5.70, 13.25)	Z = −0.323	.747
PLT (K/μL), M (Q1, Q3)	113.50 (73.00, 176.50)	113.00 (73.00, 180.00)	114.00 (73.00, 171.00)	Z = −0.455	.649
RBC (m/μL), Mean ± SD	3.24 ± 0.73	3.23 ± 0.73	3.26 ± 0.73	t = 0.960	.338
RDW (%), Mean ± SD	17.01 ± 2.57	16.95 ± 2.48	17.16 ± 2.76	t = 1.520	.128
Glucose (mg/dL), M (Q1, Q3)	126.00 (101.00, 165.00)	125.00 (100.00, 163.50)	128.00 (102.00, 166.50)	Z = 1.068	.286
Creatinine (mg/dL), M (Q1, Q3)	1.10 (0.80, 2.00)	1.10 (0.80, 2.00)	1.20 (0.80, 2.10)	Z = 1.361	.174
BUN (mg/dL), M (Q1, Q3)	25.00 (15.00, 44.00)	25.00 (15.00, 43.00)	27.00 (16.00, 47.00)	Z = 1.954	.051
Bicarbonate (mEq/L), Mean ± SD	22.19 ± 5.18	22.11 ± 5.25	22.38 ± 5.00	t = 1.050	.293
TBIL (mg/dL), M (Q1, Q3)	2.80 (1.50, 5.80)	2.80 (1.50, 5.60)	2.80 (1.55, 6.20)	Z = 0.605	.545
Hematocrit (%), Mean ± SD	30.47 ± 6.42	30.40 ± 6.43	30.64 ± 6.40	t = 0.720	.473
Hemoglobin (g/dL), Mean ± SD	10.30 ± 2.21	10.27 ± 2.21	10.36 ± 2.21	t = 0.770	.443
ALT (U/L), M (Q1, Q3)	40.00 (25.00, 71.00)	40.00 (24.00, 71.00)	41.90 (25.00, 72.00)	Z = 0.574	.566
AST (U/L), M (Q1, Q3)	80.96 (47.00, 139.00)	81.19 (47.00, 139.00)	80.00 (46.00, 137.00)	Z = −0.108	.914
Bilirubin (μmol/L) M (Q1, Q3)	2.80 (1.50, 5.80)	2.80 (1.50, 5.60)	2.80 (1.55, 6.20)	Z = 0.605	.545
Sapsii, M (Q1, Q3)	39.00 (29.00, 49.00)	39.00 (29.00, 50.00)	38.50 (29.00, 48.00)	Z = -0.443	.658
Sofa_score, M (Q1, Q3)	8.00 (5.00, 10.00)	8.00 (5.00, 10.00)	7.00 (5.00, 10.00)	Z = -−0.786	.432
COPD, n (%)				χ2 = 1.581	.209
No	1673 (90.92)	1164 (90.37)	509 (92.21)
Yes	167 (9.08)	124 (9.63)	43 (7.79)
Heart_failure, n (%)				χ2 = 2.593	.107
No	1453 (78.97)	1030 (79.97)	423 (76.63)
Yes	387 (21.03)	258 (20.03)	129 (23.37)
Diabetes, n (%)				χ2 = 0.022	.881
No	1351 (73.42)	947 (73.52)	404 (73.19)
Yes	489 (26.58)	341 (26.48)	148 (26.81)
Septicemia, n (%)				χ2 = 0.605	.437
No	1290 (70.11)	896 (69.57)	394 (71.38)
Yes	550 (29.89)	392 (30.43)	158 (28.62)

ALT = alanine transaminase, AST = aspartate transaminase, BUN = blood urea nitrogen, COPD = chronic obstructive pulmonary disease, DBP = diastolic blood pressure, PLT = platelets, RBC = red blood cell, RDW = red cell distribution width, SAPS = Simplified Acute Physiology Score, SBP = systolic blood pressure, SOFA = Sequential Organ Failure Assessment, SPO2 = Oxygen saturation, TBIL = total bilirubin, WBC = white blood cells.

Difference analysis of training set

Characteristics of death group and survival group in the training set were compared in Table 2, and results displayed the statistically significant differences between 2 groups in age, heart rate, respiratory rate, SBP, DBP, SpO2, WBC, PLT, RBC, RDW, glucose, creatinine, BUN, bicarbonate, TBIL, AST, bilirubin, SAPS II, SOFA and septicemia (P < .05).

Table 2

Baseline characteristics of training set.

Variables	Total	Survival group	Death group	Statistics	P
Age (yr), Mean ± SD	58.65 ± 12.01	57.92 ± 11.87	60.35 ± 12.18	t = −3.360	<.001∗
Gender, n (%)				χ2 = 0.389	.533
Female	432 (33.54)	296 (33.00)	136 (34.78)
Male	856 (66.46)	601 (67.00)	255 (65.22)
Heart_rate (n/min), Mean ± SD	91.15 ± 18.85	90.24 ± 18.46	93.24 ± 19.57	t = −2.630	.009∗
Respiratory_rate (n/min), M (Q1, Q3)	18.00 (14.00, 22.00)	17.00 (14.00, 20.00)	20.00 (17.00, 23.00)	Z = 7.882	<.001∗
SBP (mmHg), Mean ± SD	119.42 ± 23.28	122.19 ± 23.16	113.04 ± 22.32	t = 6.590	<.001∗
DBP (mmHg), Mean ± SD	63.25 ± 16.44	64.78 ± 15.80	59.73 ± 17.34	t = 4.940	<.001∗
SpO2 (%), Mean ± SD	97.28 ± 4.09	97.79 ± 2.65	96.12 ± 6.10	t = 5.190	<.001∗
WBC (K/μL), M (Q1, Q3)	8.69 (5.60, 13.25)	8.10 (5.40, 12.20)	10.40 (6.30, 16.40)	Z = 5.727	<.001∗
PLT (K/μL), M (Q1, Q3)	113.00 (73.00, 180.00)	117.00 (77.00, 185.00)	107.00 (66.00, 163.00)	Z = −3.044	.002∗
RBC (m/μL), Mean ± SD	3.23 ± 0.73	3.26 ± 0.72	3.14 ± 0.74	t = 2.810	.005∗
RDW (%), Mean ± SD	16.95 ± 2.48	16.57 ± 2.33	17.81 ± 2.60	t = −8.100	<.001∗
Glucose (mg/dL), M (Q1, Q3)	125.00 (100.00, 163.50)	127.00 (103.00, 167.00)	122.00 (94.00, 153.00)	Z = −3.103	.002∗
Creatinine (mg/dL), M (Q1, Q3)	1.10 (0.80, 2.00)	1.00 (0.70, 1.60)	1.60 (0.90, 2.90)	Z = 9.184	<.001∗
BUN (mg/dL), M (Q1, Q3)	25.00 (15.00, 43.00)	20.00 (14.00, 36.00)	36.00 (23.00, 58.00)	Z = 10.824	<.001∗
Bicarbonate (mEq/L), Mean ± SD	22.11 ± 5.25	22.95 ± 4.69	20.18 ± 5.92	t = 8.200	<.001∗
TBIL (mg/dL), M (Q1, Q3)	2.80 (1.50, 5.60)	2.80 (1.30, 4.10)	3.70 (2.10, 12.40)	Z = 8.761	<.001∗
Hematocrit (%), Mean ± SD	30.40 ± 6.43	30.45 ± 6.24	30.30 ± 6.87	t = 0.350	.725
Hemoglobin (g/dL), Mean ± SD	10.27 ± 2.21	10.32 ± 2.17	10.15 ± 2.28	t = 1.270	.205
ALT (U/L), M (Q1, Q3)	40.00 (24.00, 71.00)	39.00 (24.00, 70.00)	42.85 (26.00, 75.00)	Z = 0.920	.357
AST (U/L), M (Q1, Q3)	81.19 (47.00, 139.00)	79.00 (45.00, 127.00)	87.62 (52.00, 165.00)	Z = 2.599	.009∗
Bilirubin (μmol/L), M (Q1, Q3)	2.80 (1.50, 5.60)	2.80 (1.30, 4.10)	3.70 (2.10, 12.40)	Z = 8.761	<.001∗
Sapsii, M (Q1, Q3)	39.00 (29.00, 50.00)	35.00 (26.00, 44.00)	50.00 (39.00, 61.00)	Z = 15.611	<.001∗
Sofa_score, M (Q1, Q3)	8.00 (5.00, 10.00)	7.00 (4.00, 9.00)	10.00 (7.00, 13.00)	Z = 12.114	<.001∗
COPD, n (%)				χ2 = 1.211	.271
No	1164 (90.37)	816 (90.97)	348 (89.00)
Yes	124 (9.63)	81 (9.03)	43 (11.00)
Heart_failure, n (%)				χ2 = 0.011	.918
No	1030 (79.97)	718 (80.04)	312 (79.80)
Yes	258 (20.03)	179 (19.96)	79 (20.20)
Diabetes, n (%)				χ2 = 1.369	.242
No	947 (73.52)	651 (72.58)	296 (75.70)
Yes	341 (26.48)	246 (27.42)	95 (24.30)
Septicemia, n (%)				χ2 = 111.011	<.001∗
No	896 (69.57)	704 (78.48)	192 (49.10)
Yes	392 (30.43)	193 (21.52)	199 (50.90)

ALT = alanine transaminase, AST = aspartate transaminase, BUN = blood urea nitrogen, COPD = chronic obstructive pulmonary disease, DBP = diastolic blood pressure, PLT = platelets, RBC = red blood cell, RDW = red cell distribution width, SAPS = Simplified Acute Physiology Score, SBP = systolic blood pressure, SOFA = Sequential Organ Failure Assessment, SPO2 = Oxygen saturation, TBIL = total bilirubin, WBC = white blood cells.

Establishment and performance of random forest model

The data of training set were used to establish the random forest model. The variables, including age, heart rate, respiratory rate, SBP, DBP, SpO2, WBC, PLT, RDW, glucose, BUN, bicarbonate, TBIL, hematocrit, alanine transaminase, AST, bilirubin, SAPS II and SOFA were selected using logistic regression with L1 regularization. Table 3 shows the AUC of this model was 0.908 (95% confidence interval [CI]: 0.907–0.908), with the accuracy of 0.815, sensitivity of 0.836, specificity of 0.806, PPV of 0.653, and NPV of 0.919. According to the feature importance of this model, SAPS II, SOFA, BUN, TBIL, and bilirubin were more significant predictors for 30-day in-hospital deaths in ICU patients with cirrhosis (Fig. 2).

Table 3

Prediction effect of the random forest model.

	Cutoff	Accuracy (95% CI)	Sensitivity (95% CI)	Specificity (95% CI)	PPV (95% CI)	NPV (95% CI)	AUC (95% CI)
Training set	0.300	0.815 (0.794–0.836)	0.836 (0.800–0.873)	0.806 (0.780–0.832)	0.653 (0.611–0.848)	0.919 (0.900–0.938)	0.908 (0.907–0.908)
Testing set	0.300	0.734 (0.697–0.771)	0.710 (0.640–0.780)	0.744 (0.700–0.787)	0.535 (0.468–0.810)	0.861 (0.824–0.898)	0.801 (0.799–0.802)

AUC = areas under the curve, CI = confidence interval, NPV = negative predictive value, PPV = positive predictive value.

Figure 2

The importance of the included variables.

The internal validation was conducted to assess the efficacy of prediction model using the testing set and the results revealed that the AUC was 0.801 (95% CI: 0.799–0.802), the accuracy, sensitivity, specificity, PPV and NPV of prediction model were 0.734, 0.710, 0.744, 0.535 and 0.861, respectively (Table 3). The receiver operating characteristic (ROC) curves of the model and internal verification were shown in Figure 3.

Figure 3

The ROC curves of the training and testing sets.

The learning curve of random forest model was presented in Figure 4, and it showed that the prediction effect of the model on the training set and testing set tended to be steady with the increase of sample size, indicating that the performance of the model was relatively stable. The imitative effect of current model on the training set and testing set was respectively shown in Figure 5a and 5b.

Figure 4

Learning curve.

Figure 5

Calibration curve. a: training set; b: testing set.

Discussion

In this case-control study, we collected clinical data of 1840 cirrhotic patients in MIMIC-III database and established a model to predict 30-day in-hospital deaths using a random forest analysis. The final results displayed that the random forest model was effective in predicting 30-day death of ICU patients with cirrhosis. SAPS II, SOFA, BUN, TBIL and bilirubin were important factors for the 30-day death of cirrhotic patients. These results showed that our model could help clinicians identify the high-risk patients and earlier make intervention to improve patients’ prognosis.

Most of the studies tended to establish predictive models using prognostic scores to explore the 30-day outcomes of patients.[ Jacqueline et al conducted North American Consortium for the Study of End-Stage Liver Disease- Acute-on-Chronic Liver Failure Score to assess mortality risk in hospitalized cirrhotic patients. Multivariable modeling demonstrated that this score was an independently validated tool to predict 30-day survival in cirrhotic patients. The sensitivity and specificity were 84% and 70%, respectively. Huang and Yao[ established a new predictive model with combination of ascites albumin, neutrophil to lymphocyte ratio, and MELD. Through logistic multivariate regression analysis, ascites albumin, neutrophil to lymphocyte ratio, and MELD were identified as the 3 independent risk factors related to the 30-day death of patients with liver cirrhosis and bacterial ascites. The AUC of this new scoring model is 0.874. Logistic regression model has certain requirements for sample size, which theoretically requires a large sample, otherwise the test formula is unreasonable. Furthermore, logistic regression model cannot solve the problem of multicollinearity. As far as we know, there is rarely study using random forest model to predict the death of cirrhotic patients within 30 days of admission up to now. In this study, we established a random forest model to predict 30-day deaths of cirrhotic patients, and SAPS II, SOFA, BUN, TBIL, and bilirubin were found as important predictors. The random forest model could accommodate numerous variables, and had strong predictive power and better tolerance to data outliers. The AUC of our model was as high as 0.908 with the sensitivity of 83.6% and the specificity of 80.6%, indicating the good performance of the random forest model in clinical application.

The role of SAPS II and SOFA in predicting hospital mortality of ICU patients has been reported in numerous studies.[ Dupont et al[ conducted a retrospective study to assess the predictive abilities of different prognostic scores, and results revealed the superiority of SOFA and Model for End-Stage Liver Disease (MELD) score compared to other prognostic scores for mortality prediction in ICU patients hospitalized with a diagnosis of cirrhosis. SOFA was considered as the best prognostic score to evaluate cirrhotic patients in the ICU according to nearly all of the literature.[ Our study also identified SOFA as an important predictor for death, and SAPS II presented better discriminative ability for death of cirrhotic patients within 30-day hospitalization. A prior prospective study reached a conclusion that SAPA II and SOFA showed better prediction performance than MELD in ICU mortality for cirrhotic patients.[ In the future, larger sample sizes are needed to verify the priorities of different prognostic scoring systems in ICU cirrhotic patients.

Additionally, elevated BUN and bilirubin were found to be independently correlated with hospital mortality.[ Ning et al[ discussed the clinical features and prognosis in Chinese cirrhotic patients with ascites, and found the concentration of BUN was an independent risk factor for 30-day hospital mortality. The serum bilirubin level better reflects the liver's synthetic and excretory functions, thus, the mass of prognostic scoring systems included TBIL and bilirubin as ingredients.[ Our study demonstrated that BUN, TBIL and bilirubin were significant predictors for 30-day admission death. Previous studies provided a specific explanation.[ It is reported that intrahepatic cholestasis, portal flow distortion or shunting, and hemolysis caused by splenomegaly may all lead to the increased level of bilirubin. Recent research compared the value of bilirubin and TBIL for predicting prognosis of cirrhotic patients, and results showed bilirubin performed better predictive value.[

Our study used random forest model to predict the 30-day in-hospital deaths of ICU cirrhotic patients, and logistic regression was used to screen out the important variables to establish the prediction model. The internal validation had confirmed that the random forest model could perform well. However, some limitations should be concerned. First, our data were collected from the Medical Information Mart of the Intensive Care Database III version 1.4, which inevitably existed data missing. Data missing may affect the performance of the model, and some potential valuable variables for prediction may exclude due to severe data missing. Second, since our study included the patients with critically ill cirrhosis, the prediction performance in those with mild or moderate cirrhosis was unclear. Moreover, our study was lack of external validation, which needed to perform in the future.

Conclusion

In conclusion, our study developed and validated a random forest model to predict the 30-day in-hospital death for cirrhotic patients in ICU. Our model had showed a good predictive performance, and suggested more attentions on SAPS II, SOFA, BUN, TBIL, and bilirubin, indicating that it may be popularized in clinical practice. Applying this model, clinicians may be provided with an early prediction of hospital mortality and take adequate treatments in time for high-risk patients to decrease the 30-day in-hospital mortality of cirrhosis patients in ICU.

Author contributions

Conceptualization: Dongling Chen, Yuyuan Hu.

Data curation: Dongling Chen, Yuyuan Hu, Qian Li, Guichun Yin, Xianjun Zhang, Yachun Wang.

Formal analysis: Dongling Chen, Yuyuan Hu, Qian Li, Guichun Yin, Xianjun Zhang, Yachun Wang.

Investigation: Yuyuan Hu, Qian Li, Xianjun Zhang, Yachun Wang.

Methodology: Yuyuan Hu, Qian Li, Guichun Yin, Xianjun Zhang, Yachun Wang.

Project administration: Dongling Chen, Yuyuan Hu, Guichun Yin.

Writing – original draft: Dongling Chen, Yuyuan Hu.

Writing – review & editing: Dongling Chen, Yuyuan Hu.