Meteorological Factors

In this study, we applied the time-stratified case-crossover (ts-CCO) design, regarded as a self-matched case-control study, which compares the exposure in the case period when events occurred with exposures in nearby referent periods, to examine the differences in exposure which may contribute to the differences in the daily count of cases [18 (link)]. Therefore, a time-stratified case-crossover design was adopted to regulate potential confounders (e.g., age, gender, etc) using self-control and exclude long-term impact of air pollutants (e.g., secular trend, seasonality, etc.) by stratification of time. We used the calendar month as the time stratum, to control the effects of long-term trend, seasonality, and day of the week [17 (link)]. We conducted a quasi-Poisson regression, controlling over-dispersion problem, combined with distributed lag non-linear model (DLNM) to estimate the non-linear and delayed influence of air pollution on IS onset. DLNM is based on the definition of “cross-basis”, a bi-dimensional space of functions to reflect the non-linear exposure-responses and lag structure of the association [19 (link), 20 (link)]. In this study, consequently, the combination of DLNM with the ts-CCO design was employed, which allows estimating the short-term, non-linear and delayed effect of air pollutant using cross-basis functions for depicting the relationship between air pollutant and IS onset along the dimensions of exposure and lag simultaneously based on removing control confounders and long-term trend by ts-CCO design.
A quasi-Poisson regression model combined with time-stratified case-crossover design and DLNM was built as follows:
$$Y_t~\mathit{\text{Poisson}}\left({\mu }_t\right)$$
$$\mathit{Log}\left({\mu }_t/{\mathit{\text{Population}}}_{\mathit{ye}}\right)=\alpha +\sum _{i=1}^m\mathit{cb}\left({\mathit{\text{Pollutant}}}_{i,t},{\mathit{df}}_{2i-1},{\mathit{\text{maxlag}}}_i,{\mathit{df}}_{2i}\right)+\mathit{cb}\left({\mathit{\text{Temp}}}_t,{\mathit{df}}_{2m+1},{\mathit{\text{maxlag}}}_{m+1},{\mathit{df}}_{2m+2}\right)+\mathit{cb}\left({\mathit{RH}}_t,{\mathit{df}}_{2m+3},{\mathit{\text{maxlag}}}_{m+2},{\mathit{df}}_{2m+4}\right)+\gamma {\mathit{\text{Holiday}}}_t+\lambda \mathit{\text{Stratum}}$$ where t is the day of observation; Y_t is the count of IS cases on t; μ_t is the expectation of Y_t; Population_ye is the year-end population size; α is an intercept; Pollutant_{i, t}, Temp_t, and RH_t are the ith pollutant concentration, temperature and relative humidity on t, respectively; cb() represents the cross-basis function with three pre-specified parameters of maximal lag maxlag_i, degree of freedom for lag-response natural spline df_2i − 1, and degree of freedom for exposure-response natural spline df_2i for pollutant, temperature or relative humidity; Holiday is used to control the effect of public holidays; Stratum is the time stratum in the time-stratified case-crossover design. We defined natural cubic spline function with 3 df for air pollution and meteorological factors to mimic the exposure-response pattern of air pollution-IS onset associations, as well as lag spaces with 3 df to estimate the lag effects. To capture the complete lag-response curve, the maximal lag of air pollutants was set to 14 days; for the sake of simplification and without loss of generalization, meanwhile, this maximal lag was assigned to the length of the case and control periods. In addition, a 3-day duration was specified to be the maximal lag of meteorological factors. The df and maximum lag days for air pollution determination referred to the Akaike information criterion for quasi-Poisson (Q-AIC), which could produce the relatively superior model.
We initially conducted single-pollutant model to evaluate the association between air pollution and IS onset, and then the significant air pollutants were included in multi-pollutant model. Spearman’s correlation tests were used to estimate the associations between air pollution and meteorological factors, and pollutants with correlation coefficient r > 0.60 were not included in multi-pollutant model simultaneously to address the collinearity between air pollutants. In order to identify the high-risk or low-risk air pollution condition, the influence of extreme air pollution was evaluated and presented as relative risk (RR) by comparing the 99th above or 1st below percentiles of air pollution to the median values. We calculated the single day lag influence and the cumulative lag influence (lag0–1, lag0–6, lag0–8, lag0–10, lag0–12, lag0–13, and lag0–14) to effectively depict the characteristics of the association between air pollution and IS onset. In addition, we conducted stratified analysis to investigate the impact of air pollution on subgroups according to gender (male and female) and age groups (adult: 18–64 years; the elderly: ≥ 65 years).
All analyses were conducted using R version 3.5.1 with the dlnm package for fitting the DLNM, the gnm package for conditional quasi-Poisson regression.
Sensitivity analyses were performed to test the robustness of the selected model, which were as following: df [2 (link)–6 (link)] for air pollution, and df [2 (link)–6 (link)] for lag space were changed; the maximum lag days (12–21 days) for air pollution were extended.

We conducted several sensitivity analyses to check the robustness of our results. Firstly, we repeated the analyses using alternative df values (from 3 df to 4, 5, and 6 df) for mean temperature, relative humidity, and precipitation. Secondly, we used alternative moving average lag structures for all meteorological factors: lag 0–1 (current month and preceding 1 month) and lag 0–2 (current month and preceding 2 months). Thirdly, we refitted the CITS models using data from 2017 to 2020 to avoid inconsistencies in trends between different periods. Finally, for diseases with adequate sample sizes with the elderly, we restricted analyses to people aged between 20 and 54 and people aged above 55, respectively. We used the statistical software R 4.0.1 (Lucent Technologies, Jasmine Mountain, USA) to perform all analyses. A two-sided P-value < 0.05 was considered statistically significant.

The Theil-Sen Median (Sen) method was employed to determine the trend of NEP over the study area. This method is a robust nonparametric trend method, which does not require the data to follow a certain distribution (Fensholt et al., 2012 (link); Zhang Z. et al., 2021 (link)). It has been widely used in the trend analysis of long-time data series.
The Sen’s slope value (β) indicates the magnitude of the NEP trend. A positive β value suggests an upward trend and a negative β value suggests a downward trend. The calculation formula for β is as follows:
where β is Sen’s slope value; NEP_i and NEP_j are NEP in year i and j respectively.
The Mann–Kendall test is used to assess the significance of NEP trends. The test statistics S value is calculated as:
The test statistic Z value was obtained by standardizing S:
Where the function var is calculated as:
where n is the number of sequence samples; m is the number of affiliated groups in the data sequence; t_i is the number of input values inside the affiliated group i.
To validate the significance of the trend, the absolute z-score value |Z| is compared with the critical value Z_1-α/2 at a given significance level α. if |Z| is higher than Z_1-α/2, the trend is considered significant. The Z_1-α/2 values were obtained from the standard normal distribution table. For the significance level α of 0.05 and 0.01, the critical Z_1-α/2 values are 1.96 and 2.58, respectively
Partial correlation is used to analyze the relationship between NEP and meteorological factors. Its formula is:
where R_ij,k is the partial correlation coefficient between variable i and j after variable k is fixed. r_ij, r_ih, r_jh are correlation coefficients for the variables i and j, j and k, and k and i, respectively.

Standard descriptive analysis was conducted to display the distribution of air pollutants, meteorological factors, and daily hospitalizations for respiratory diseases, while the temporal distribution of air pollutants was presented in line charts. Spearman's correlation analysis was used to identify the correlation between air pollutants and meteorological factors.
The generalized additive model was applied to quantify the association between daily ambient CO concentration and daily hospitalizations for respiratory diseases. Quasi-Poisson regression was applied in the model, as daily hospitalizations tended to display an over-dispersed Poisson distribution. A dichotomous variable for public holidays and a categorical variable for the day of the week (DOW) were incorporated into the model to adjust the variation of daily hospitalizations within holidays and each week. Moreover, smoothing terms were used to fit daily hospitalizations in the models to control long-term and seasonal trends of daily hospitalizations and meteorological effects (14 (link)). According to previous studies (15 (link)–17 (link)), we applied 6 degrees of freedom (df) per year for long-term and seasonal trends, 3 df each for the same day's temperature (Tem0) and relative humidity (Humid0). In brief, the model can be represented as follows:
where E(Y_t) represents the estimated daily hospitalizations for respiratory diseases at day t. β represents the log-relative risk of hospitalization associated with a 1 mg/m³ increase in ambient CO concentration. s () is the restricted smoothing spline function for variables with the non-linear association, day indicates the variable of the long-term and seasonal trends, and α is the intercept for the model.
Considering the delayed health effects of air pollutants, we estimated the lag effects of different days in both single-day lag from lag0 to lag7 and multi-day lag from lag0 to lag0–7 (moving average from lag0 to lag7). To improve the comparability of the association between CO exposure and risk of hospitalization for total and specific respiratory diseases, we selected the same CO exposure window for different respiratory diseases in exposure–response relationship analysis, stratified analysis, and sensitivity analysis adjusted for co-pollutants.
Based on the same models that estimated the association between CO exposure and risk of hospitalization for respiratory diseases, the smoothing function with 3 df was used to graphically describe the exposure–response relationship between ambient CO concentration and risk of hospitalization for respiratory diseases. Stratified analysis was conducted to assess the potential effect modification by age (minors/adult/elderly), gender (men/women), and season (warm: May–October, cold: November–April) (18 (link), 19 (link)). To further quantify the potential effect modification, we calculated the significant differences between subgroups based on the widely used method (20 (link)). We also investigated whether the association between CO and hospitalizations was still robust to the adjustment for other co-pollutants including PM_2.5, PM₁₀, SO₂, NO₂, and O₃. Dual- and multi-pollutant models were performed in this study.
In this study, two-sided P < 0.05 was considered statistically significant. The generalized additive model was conducted in R 4.0.2 within the “mgcv” package (21 (link)). Effect estimates were presented as percentage changes and 95% confidence intervals (CIs) in daily hospitalizations in relation to each 1 mg/m³ increase in ambient CO concentration.