The study was based on a previously reported cohort that recruited patients with pneumonia who were hospitalized in six secondary and tertiary academic hospitals in China between January 1, 2013 and December 31, 2017 (17). The diagnosis of pneumonia was based on the American Thoracic Society guidelines and guidelines for infectious diseases (6, 18). The study inclusion criteria were as follows: (1) oral or intravenous glucocorticoid therapy, (2) diagnosis of pneumonia during admission and hospitalization, and (3) age 18 years or older. Exclusion criteria were as follows: (1) diagnosis of non-infectious lung disease, including lung cancer, interstitial lung disease without infection, pulmonary embolism, or heart failure, and (2) inability to provide informed consent. The Ethics Committee of China–Japan Friendship Hospital approved the retrospective study and organized centralized collaboration and approval of all participating institutions. The flow chart is shown in Figure 1.

The original data were downloaded from the DRYAD database (https://datadryad.org/stash). The dataset contained 716 samples and 127 observations. First, patients without CTD were removed; second, clinical observations with missing values greater than 5% were removed (Supplementary Table 1), and then the missing values were filled by multiple interpolations using the R package “mice”(v3.1.4) (19). We screened for clinical variables including general condition, comorbidities, biochemical indicators, and medication use. Finally, the samples were randomly divided into two parts, with 70% of the samples as the training cohort and 30% as the validation cohort.

A univariate Cox regression analysis was used to screen for prognostic variables by using the R package “survival”(v3.5–0) (20). Prognostic variables were entered into the least absolute shrinkage and selection operator (Lasso) by using the R package “glmnet”(v4.1–3) (21). We used the “CV.glmnet” function (parameter default) for ten-fold cross-validation, and selected lambda.1se screened variables as the Lasso prognostic variables. Meanwhile, prognostic variables were also entered into a random survival forest (RSF) analysis to screen important prognostic variables based on the “rfsrc” function (parameters were set as mtry = 6, ntree = 1000) in R packet “randomForestSRC”(v3.1.0) (22). The top 10 variables with the greatest importance were selected as the prognostic variables of the RSF analysis. Overlapping prognostic variables were obtained by intersecting Lasso prognostic variables with RSF prognostic variables. The overlapping prognostic variables were entered into stepwise Cox regression (direction = both) to screen the hub prognostic variables, and then prognostic models were constructed based on the hub prognostic variables. The following risk score formula was obtained from the model (23):

where n, Expi, and βi indicate the number of hub prognostic variables, variable values, and the stepwise Cox regression coefficient, respectively. In the training cohort, patients were divided into the high- and low-risk groups based on the median risk score, and the difference in prognosis between the two groups was assessed using the Kaplan–Meier analysis.

We constructed prognostic nomogram using the “nomogram” function (default parameters) by package “rms”(v6.2.0) (24). Calibration curves were plotted and the C-index was calculated to assess the predictive efficacy of the model. To assess the clinical benefit of the model, we plotted decision curve analysis (DCA) curves using the “dca” function in the R package “ggDCA”(v1.2) (25). In addition, to test the predictive power of the model in specific patients, we performed subgroup analyses of clinical variables (age, gender, diabetes mellitus, and interstitial lung disease). In the validation cohort, a Kaplan–Meier analysis, an evaluation of C-index, and the plotting of calibration curve were used to test the predictive power of the model.

The study was statistically analyzed using R language (version 4.1.3). The mean ± standard deviation was used to represent continuous normally distributed variables, and the interquartile range was used to represent non-normally distributed information. A student’s t-test was used for the comparison of normally distributed random variables. The Wilcoxon rank test was used to compare non-normally distributed variables. In addition, categorical data were compared using the chi-square test or Fisher’s exact test. Stepwise Cox regression analyses and prognostic model constructions were applied using the R package “survival”. A survival analysis was performed using a Kaplan–Meier analysis and log-rank test.

Data for a total of 712 pneumonia patients treated with glucocorticoids, encompassing 127 variables, were downloaded from DRYAD. A total of 368 patients with CTD pneumonia were selected, and 84 variables had missing values of less than 5%. After multiple interpolations to supplement missing values, 37 variables were selected for analysis. These included age, gender, fever, cough, expectoration, chest pain, dyspnea, heart rate, respiratory rate, systolic pressure, diastolic pressure, cyanosis, CAP, bronchiectasia, interstitial lung disease, coronary heart disease, congestive heart failure, diabetes mellitus, anemia, white blood cells, neutrophils, lymphocyte, hemoglobin, platelets, aspartate aminotransferase, alanine aminotransferase, blood urea nitrogen, serum creatinine, K, Na, oxygen inhalation, high-dose glucocorticoid, immunosuppressant, ganciclovir, sulfanilamide, anti-aspergillus, and anti-pseudomonas. Patients were randomized into two parts, with 247 patients in the training cohort and 121 patients in the validation cohort. The groups were comparable at baseline except for the need for sulfanilamide treatment (Table 1).

A univariate Cox regression analysis obtained 19 prognostic variables (Table 2). We included 19 prognostic variables in the Lasso Cox regression to reduce dimensionality. Ten-fold cross-validation was performed by the minimum criterion to the preferred parameters (Figure 2A). The trajectory of each prognostic index coefficient with the log-transformed lambda in the Lasso algorithm was observed (Figure 2B). We selected 11 variables corresponding to lambda.1se=0.066 as Lasso prognostic variables. Nineteen prognostic variables were also entered into the random prognostic forest (Figures 2C, D), which was eventually filtered to 10 significant prognostic variables as RSF prognostic variables based on variable importance scores. Eight overlapping prognostic variables (fever, respiratory rate, cyanosis, white blood cells, blood urea nitrogen, ganciclovir, anti-aspergillus, and anti-pseudomonas) were obtained by the two algorithms. The overlapping prognostic variables were put into a stepwise Cox regression analysis to obtain five variables (fever, cyanosis, blood urea nitrogen, ganciclovir, anti-pseudomonas), and a prognostic model was constructed based on these five variables. The model equation was risk score=1.389* Fever+0.990* Cyanosis +0.065* BUN +0.568* Ganciclovir +2.206*Anti-pseudomonas. The model coefficients and risk values for the fever variables can be seen in Table 3.

The training cohort was included in the risk model and risk scores were calculated; patients were then divided into the high- and low-risk groups based on the median risk score, with significantly lower survival rates in the high-risk group than in the low-risk group (Figure 3A). Five independent variables were included in the nomogram to assess prognosis according to the results of the multivariable Cox regression analysis (Figure 3B). Predictive models performed well in predicting mortality: the C-index score was 0.808. Calibration curves for nomogram prediction of 30-day, 60-day, and 90-day overall survival were in good agreement with actual observations (Figures 3C–E). DCA results showed that the net benefit level of 30 days (risk threshold: 0~0.85), 60 days (risk threshold: 0~0.85), and 90 days (risk threshold: 0~0.85) predicted using the line chart was significantly higher than that of “no intervention” and “total intervention”, suggesting that the model had better clinical predictive efficacy (Figures 3F–H). In the subgroup analysis of age, gender, and comorbidities (interstitial lung disease, diabetes mellitus), survival time in the high-risk group was significantly lower than that in the low-risk group (P<0.001, Figures 4A–H). In conclusion, the above results indicate that our model has good predictive power.

To verify the stability of the model, this study applied the model to the validation cohort. The high- and low-risk groups were divided according to the median value of the risk score, and the survival of the high-risk group was significantly lower than that of the low-risk group (Figure 5A). The prediction model performed well in predicting mortality: the C-index score was 0.762. The calibration plots of the model predicting overall survival at 30, 60, and 90 days agreed well with the actual observations (Figures 5B–D).