Research location
Sichuan Province, an inland region in Southwest China, has exhibited continuous increases in population each year since 2005; currently, it contains more than 80 million inhabitants. This province covers approximately 486,000 km2 and has 21 prefecture-level cities. The region’s complex topography consists of a basin in the eastern region, a plateau in the western region and a mountainous area in the southwestern region (Fig. 1). The Sichuan Basin, one of the four largest basins in China, has a dense population and a relatively high level of economic development, while western Sichuan is mostly mountainous and sparsely populated. Additionally, there are obvious regional differences in climate. The Sichuan Basin is situated in a subtropical humid monsoon climate zone, in which the average annual temperature is 16–18 ℃ and the average annual precipitation is 1,000–1,300 mm. In addition, the province contains a subhumid subtropical region and an alpine region in the northwest plateau, with long sunshine durations and little annual precipitation. The high levels of emissions, special terrain, and unique climate conditions have resulted in severe air pollution in winter, especially over the eastern basin. The Sichuan Basin accounts for less than 2.7% of the area of China. However, the total emissions of SO2, PM2.5, and NOx account for 12.1, 8.3, and 5.8% of China’s total emissions, respectively [26].
Data sources
On May 2, 2008, the Chinese Ministry of Health included HFMD in the catalog of class C notifiable infectious diseases. Information reporting and case management should be carried out within 24 h according to the standard guidelines [27]. The demographic information and illness-related information of each case in 21 prefecture-level cities in Sichuan Province between 2015 and 2017 were collected from the Reporting System of the Chinese Center for Disease Control and Prevention. We extracted the dates that cases presented with symptoms and aggregated them into the daily HFMD counts for the respective 21 cities in Sichuan Province. According to the preliminary analysis, over 99% of reported HFMD patients were less than 15 years old; therefore, this study included only data from patients under the age of 15 years.
We retrieved daily meteorological monitoring data, including ambient temperature (°C) (mean, minimum, and maximum), mean relative humidity (%), mean wind speed (m/s), total sunshine duration (h), and mean atmospheric pressure (hPa) data, from the China Meteorological Data Sharing Service System. Daily surveillance data for air pollution, namely, SO2 (μg/m3), NO2 (μg/m3), PM2.5 (μg/m3), PM10 (μg/m3), the air quality index (AQI), CO (mg/m3), and O3 (μg/m3) data, were retrieved from the Sichuan Environmental Monitoring Center. The spatial distributions of the meteorological monitoring stations and air quality monitoring sites are shown in Fig. 1. Considering the nature of climatic and air pollution time-series data, a very small proportion of missing values (< 0.1%) were replaced with zeroes for sunshine hours and linear interpolation for other variables. Following the basic principle of matching data from sites closest to the city center, daily HFMD counts, daily meteorological monitoring data, and daily air pollution data were then matched by city-specific codes.
Additionally, we obtained city-level economic (gross domestic product (GDP) per capita), demographic (population density, population growth and number of primary school students), traffic (travel passengers) and health resource characteristics (number of health institutions, number of hospital beds and number of registered physicians) from the Sichuan Statistical Yearbooks during 2015–2017.
We calculated the arithmetic average of each air pollution variable and socioeconomic indicator in the prefecture-level cities as the unit for the whole research period. Then, the arithmetic average of the same variable for each city was designated as a new substitutive variable, denoting the differences in air quality or socioeconomic characteristics among cities. These city-specific variables were included as potential modifiers in the second-stage meta-regression models to assess the modification effects of long-term air pollution levels, which is also known as the explanation of the heterogeneity [28].
Statistical analysis
In this study, daily HFMD counts, meteorological and air pollution surveillance data from 21 cities in Sichuan Province between 2015 and 2017 were collected, and a multicity time-series approach with two stages was adopted. Epidemiological studies have demonstrated that meteorological and air pollution variables may exhibit complex nonlinear associations with HFMD incidence. To accurately capture the short-term effects of these variables, a unified DLNM with the same modeling structure and parameter determination was first applied for each prefecture-level city to simultaneously quantify the city-specific effect estimates of the exposure and lag dimensions. Then, multivariate meta-regression was performed to combine city-specific effect estimates, and air pollution variables at the city level were fitted as a meta-predictor into models to examine potential modification effects.
First-stage analysis
A DLNM based on a quasi-Poisson distribution served as the basic model for detecting possible delayed effects and nonlinear associations between exposures and HFMD incidence for each city in the first stage of analysis [29]. We constructed univariate models to screen for significant associations of meteorological variables and air pollutants with HFMD cases. Then, these significant variables were incorporated into the subsequent multivariate model. Parameters and modeling components were determined on the basis of prior knowledge and a systematic sensitivity analysis strategy [10]. We constructed the final model as follows:
$$\begin{array}{l}{Y}_{t}\sim Quasi-Poisson\left({\mu }_{t}\right)\\ \mathrm{log}\left({\mu }_{t}\right)=\alpha +\sum cb\left(Weather/air pollution, lag\right)+ns\left(Day,df\right)+Auto+Dow+Hod\end{array}$$
where Yt denotes the daily number of HFMD cases at time t in a prefecture-level city and cb(Weather / air pollution, lag) stands for a cross-basis function of each of the climatic variables and air pollutants that describes the exposure-lag-response relationship. We defined 3 degrees of freedom (dfs) for the natural cubic splines of temperature, relative humidity, wind speed, NO2, PM10, and ozone in the exposure dimensions, as well as 4 dfs for delayed effects. Considering the incubation period (3–5 days) and the infectious period (nearly 2 weeks) of HFMD [27], the lag range was set to 0–17 days to adequately cover delayed effects. Long-term and seasonal trends in HFMD were eliminated using 8 dfs/year for the natural cubic spline function. Due to the mechanism by which HFMD is transmitted, first- and second-order lag terms on the logarithmic scale of HFMD counts were incorporated to control for residual autocorrelations and are expressed as Auto [30]. The day of the week is represented by Dow. Hod is a binary variable that controls for the effect of national public holidays. The median value of each exposure was assigned as the reference to calculate the relative risks.
Second-stage analysis
The cumulative exposure–response relationship in each city was obtained by combining all lag effects to achieve parametric dimensionality reduction for the second stage of analysis. A multivariate meta-regression model containing only intercepts was constructed to capture the overall pooled exposure–response relationship in Sichuan Province and to examine the heterogeneity of city-specific associations. Then, city-specific long-term air pollutant levels and socioeconomic variables were individually added as meta-predictors to fit the meta-regression models with a single meta-predictor to further estimate the modification effects [31]. Variables with P < 0.2 in the single meta-predictor model were selected as alternative factors for the subsequent multiple meta-predictor analysis. All candidate subsets were considered and the best subset of meta-predictors with the lowest AIC was identified. Quantitative statistical indicators, namely, the results of the likelihood ratio test that assessed whether the meta-predictor had statistical significance and the extent of improvement in the I2 statistic that indicated the explicable proportion of residual heterogeneity [28], were computed to evaluate the modification effects of long-term air pollution levels.
In addition, we also performed a subgroup analysis based on the three defined regions (the Sichuan Basin, Southwest Sichuan Mountain region, and West Sichuan Plateau). We pooled the effect estimates for each city by region using a multivariate meta-regression model and then visualized the overall pooled exposure–response curves for each region. Limited by the number of cities in the Southwest Sichuan Mountain Region and West Sichuan Plateau; thus, we further explored the potential modification effects of long-term air pollution levels on the associations between environmental factors and HFMD in only the Sichuan Basin.
Multiple climatic variables and air pollutants were assessed in this study. Considering the strong positive correlations among the mean, minimum, and maximum temperatures and among the PM10, PM2.5, and AQI as well as the strong negative correlation between the mean temperature and air pressure across all 21 cities, only the mean temperature and PM10 were included in the final analysis to avoid collinearity. In addition, a considerable proportion (34.6-54.2%) of the actual zero value for sunshine duration was present in multiple cities due to the special topography of the eastern basin; in the southern region of the basin, the proportions approached or even exceeded 50%. When we constructed DLNMs for these cities, the cross-basis matrixes produced multiple collinear variables. Some of these collinear variables were discarded, and the related parameters could not be estimated [31]. Thus, sunshine duration was not included as an exposure variable in the first-stage analysis to ensure the accuracy of parameter estimations. We carried out all analyses with R software (version 4.0.3), with the package dlnm to fit all DLNMs and the package mvmeta to conduct all multivariate meta-regression models. We constructed geographic maps using ArcGIS software (version 10.0, authorization number: EFL734321752).