Study area
This study was conducted in Hefei, which is the capital and largest city of Anhui province in Eastern China with a population of 8.09 million inhabitants (from 2018 census data). Hefei has a humid subtropical climate with a mean temperature of 16.8 °C.
Arthritis data
Daily counts of outpatient admission for OA/RA during 2014–2017 were obtained from The First Affiliated Hospital of University of Science and Technology of China (Anhui Provincial Hospital). The patient data included the date of outpatient admission, age, gender, residential address. Diagnosis of OA (ICD-10: M13.9) and RA (ICD-10: M06.9) was coded according to the International Classification of Disease, 10th Revision (ICD-10). Ethical approval was obtained from the Ethics Committee of Anhui Provincial Hospital prior to data collection.
Weather and air pollutants data
Meteorological data on daily mean temperature, relative humidity, rainfall, barometric pressure and wind velocity during the same period were obtained from Hefei Bureau of Meteorology. Air pollution data including the average daily level of sulfur dioxide (SO2), nitrogen dioxide (NO2), carbon monoxide (CO), ozone (O3), particulate matter of less than 10 μm and 2.5 μm (PM10 and PM2.5) were collected from the Environmental Protection Bureau in Hefei. Consistent with previous study [17], we chose the 50th percentile of temperature (P50, 17.8 °C) as the reference in analyses.
Statistical analysis
We first examined the correlations among weather indicators and pollutants with Spearman’s correlation test. Then, we applied a Poisson generalized linear regression combined with distributed lag non-linear model (DLNM) to examine the non-linear and lagged effects of ambient temperature on outpatient admission for OA/RA, after controlling for long-term trend and seasonality, day of week (DOW), public holidays (Holiday), relative humidity, wind velocity, PM2.5, SO2, NO2 and O3. The core model is expressed as follows:
$$ {Y}_t\sim Poisson\left({\mu}_t\right)\ Log\left({\mu}_t\right)=\alpha + cb\left({Temperature}_{t,l}\right)+ cb\left({Humidity}_{t,l}, df=3\right)+ cb\left({Wind}_{t,l}, df=3\right)+ cb\left({PM}_{2.5t,l}, df=3\right)+ cb\left({SO}_{2t,l}, df=3\right)+ cb\left({NO}_{2t,l}, df=3\right)+ cb\left({O}_{3t,l}, df=3\right)+ ns\left({Time}_t, df=8\right)+{\ng\mathit{DOW}}_t+{\gamma Holiday}_t+\mathrm{Lag}\left({Y}_t,1\right)+\mathrm{Lag}\left({Y}_t,2\right) $$
Where Yt is the number of OA/RA admission on day t; α represents the intercept;
cb() is a cross-basis function used to models both the exposure effect and lag effect at the same time; l refers to the lag days; ns() denotes a natural cubic spline function using the dlnm package in R. To control long-term time and seasonality, we used a natural cubic spline function with 8 degrees of freedom (dfs) per year, along with an indicator of the day of the week (DOW) and holiday effect. A natural cubic spline with 3 dfs was used to for exposure dimension (mean temperature, humidity, wind velocity, PM2.5, SO2, NO2 and O3) and lag dimension (lags 0–4). In order to reduce the influence of model autocorrelation, we added autoregressive terms to the model to improve the fit of the model.
On the basis of the lowest Akaike Information Criterion (AIC), we selected the maximum lag of 4 days to capture any single and cumulative effects of temperature. Because the plot of overall exposure-response did not find the significant relationship between temperature and OA admission (Fig. 1), we only quantified the relative risks (RRs) of temperature change on RA admission by single day lags at low temperature (25th percentile, P25) compared to the reference temperature (50th percentile, P50). Furthermore, we examined the specific cumulative effects of temperature decrease on RA admission by gender (male and female) and age (0–17 years, 18–40 years, 41–65 years and ≥ 66 years). The statistically significant differences between effect estimates in subgroups were examined by the following formula:
(\( \hat{Q1} \)- \( \hat{Q2} \)) ± 1.96 \( \sqrt{{SE_1}^2+{SE_2}^{2.}} \)
where \( \hat{Q1} \) , \( \hat{Q2} \) are the estimates for two categories (such as male and female), and SE12, SE22 represent the corresponding standard errors [18]. The effects of temperature were estimated and reported as RR and its 95% confidence interval (CI) associated with low temperature at different lags.
To test the robustness of our results, sensitivity analyses were performed by varying df for time (7–9 dfs/year), humidity (3–5 dfs) and wind velocity (3–5 dfs), respectively. Data manipulation and analyses were conducted using R software (version 3.1.1), with the “dlnm” package to fit the DLNM [16].