Statistical analysis

Participants’ characteristics were described by their systemic steroid status as frequencies and proportions for categorical variables and as means with standard deviations or medians and interquartile ranges (IQR) as appropriate for continuous variables. As for the primary outcome, we used a logistic regression model, and for the secondary outcomes, we used a Cox proportional hazard model for time to in-hospital death and a gamma regression model for length of hospital stay [15, 16]. We estimated the average differences in 30-day mortality and length of hospital stay for every single patient based on the presence of systemic steroid therapy. To account for the confounders, we calculated the propensity score, which is the probability of a patient receiving systemic steroid therapy conditioning on confounders, using a logistic regression model with adjustment for known confounders [17]. Each patient was assigned an average treatment-effect weight calculated using the inverse probability of treatment weighting (IPTW) method. We illustrated a love plot to confirm that the absolute standardized mean difference of each covariate was <0.1 in a weighted dataset and a histogram to confirm similar distributions of propensity scores in both treatment groups [18]. We then performed weighted outcome analyses according to the outcome type. Missing data were imputed using multiple imputations by chained equations, which created 100 imputed datasets and combined the estimate within each dataset [19]. Detailed information regarding our statistical analysis is provided in the Supporting information. Finally, we calculated the risk differences in the number of cases showing acidosis and Clostridioides difficile colitis.

For sensitivity analyses, we used the following statistical models for the primary outcome: 1) a different patient selection algorithm that did not exclude patients with other concomitant respiratory diseases or daily steroid users; 2) additionally excluding non-empirical antibiotic users; and 3) a complete case analysis. In Japanese administrative claims data, we can collect the A-DROP score, a modified version of CURB-65, at admission for patients with primary diagnosis of bacterial pneumonia [20]. Thus, we performed the main analysis among the subgroup of patients with the admission-precipitating diagnosis of bacterial pneumonia (ICD-10 codes J12, J13, J14, J15, J16, J18, J69, and P23) comorbid with COPD present at the time of admission (ICD-10 codes J44.1 and J44.9) and a sensitivity analysis additionally incorporating oxyhemoglobin saturation measured by pulse oximetry ≤90% or partial oxygen pressure in arterial blood ≤60 mmHg and systolic blood pressure ≤90 mmHg extracted from the A-DROP score. We used R version 4.2.3 (R Foundation, Vienna, Austria) for the statistical analyses. A two-sided p-value less than 0.05 was considered statistically significance.