Meteorological and Air Pollution Data
Guangdong Province is located in the South of China. Our study was limited to 21 cities of Guangdong Province (179,700 square kilometers) - Zhanjiang, Maoming, Yangjiang, Zhaoqing, Shaoguan, Heyuan, Meizhou, Qingyuan, Yunfu, Shantou, Shanwei, Chaozhou, Jieyang, Shenzhen, Zhuhai, Foshan, Jiangmen, Dongguan, Zhongshan, Huizhou and Guangzhou. Data on daily maximum, minimum and mean temperatures and relative humidity (RH) were collected from the National Meteorological Information Center of China (http://data.cma.cn/). There were 68 local weather stations recorded daily measures across Guangdong Province from January 1, 2013 to December 31, 2017. We calculated city-wide meteorological measures by averaging data from stations located in a specific city. Without local weather stations, data of Chaozhou city and Foshan city were collected from the nearest monitoring sites located in Jiexi districts and Gaoyao districts, respectively. DTR was calculated by subtracting daily minimum temperature from daily maximum temperature. Data on daily city-wide concentrations of particulate matter < 2.5 µm diameter (PM2.5) were obtained from the Guangdong Provincial Environmental Monitoring Center. All data from 102 central monitoring stations in 21 cities were available.
Hospitalization data
In China, only second and tertiary level hospitals that are qualified to provide specialized medical-care for exacerbation of COPD, asthma and bronchiectasis
There were 227 government-tiered second or tertiary hospitals in 21 cities that have uploaded their daily hospitalization records to the electronic medical record system of Guangdong Government Affairs Service Center. International Classification of Diseases 10th (ICD-10) codes including J44, J45-46 and J47 were used to identified hospitalization for exacerbations of COPD, asthma and bronchiectasis, respectively.
Statistical Analysis
Spearman correlation analysis was performed among DTR, daily mean temperature, RH and daily concentration of PM2.5. The association of DTR with daily hospitalization for exacerbation of CRD was estimated by a two-stage analysis using 21-city data.
In the first stage, we adopted a standard generalized additive model (GAM) with a quasi-Poisson distribution [22] to investigate the city-specific relationship of DTR and hospital admissions for exacerbations of CRD. Seven-day moving average (lag 0-6) was used to present the lag effect of DTR. In the framework of distributed lag non-linear model (DLNM) function, we used a cubic spline for DTR and its lag with 5 and 4 degrees of freedom, respectively. The 25th percentile of DTR is regarded as the centering point. We used a natural cubic function with 8 degrees of freedom per year to control the long-term trends of years and seasonality [23, 24]. Day of the week and the official holiday were included as an indicator to remove the effect of short-term fluctuation [24]. Confounding meteorological measures, including daily mean temperature at lag 0-14, RH at lag 0-3 [26, 27] and daily concentration of PM2.5 at lag 0-3 was also included into the GAM as follows:
Log[E(Yt)] = α + βDTR0-6t + ns(time, df = 8/per year) + day of the week + holiday + ns(temperature0-14t, df = 3) + ns(RH0-3t, df = 3)+ ns(PM2.5 0-3t, df=3)
Where E(Yt) presents the daily hospital admissions for exacerbation of CRD on day t; α is the intercept in specific-region; β is the regression coefficient, and its exponent value indicates the relative risk of hospitalization per unit increase of DTR; DTR0-6t is 7-day moving average of DTR; time presents the long term trend (from 1 to 1826); temperature0-14t is 15-day moving average of daily mean temperature; RH0-3t is 4-day moving average of daily relative humidity; PM2.5 0-3t indicates 4-day moving average of daily concentrations of PM2.5; ns() presents natural cubic function.
In the second stage, we pooled the estimates in 21 cities by performing a meta-analysis. Considering the heterogeneity of city-level estimate, we adapted the random-effects model by maximum likelihood (REML) rather than fixed effects to ensure a more robust estimation.
To explore the seasonal pattern of relationship between DTR and hospitalization for exacerbation of CRD, data was stratified by time: hot season (May to October) and cold season (November to April of the next year). We used a natural cubic function with 4 degrees of freedom within each 6-month subperiod to control the variation of long-term trends [28]. We also divided the hospitalization data into diverse CRD groups to investigate the heterogeneity of associations of DTR with COPD, asthma and bronchiectasis, respectively. Subgroup analyses on sex and age (i.e., <65 vs. ≥65 years old) were performed to identify vulnerable population.
Sensitivity analysis
Several sensitivity analyses were applied to evaluate the robustness of main results: (1) altering numbers of moving average for DTR; (2) varying df of the long-term trend and meteorological measures; (3) replacing the natural cubic spline function by penalized splines function; (4) excluding data of cities with ≤ 5 hospitals in records.
Analyses were performed in STATA (version 12, StataCorp, TX) and R (version 3.6.2, R Development Core Team) with “mgcv” and “mvmeta” packages. Statistically significance was determined as a two-side P value < 0.05.