Spatio-temporal analysis of the relationship between climate and hand, foot, and mouth disease in Shandong province, China, 2008–2012

Background Hand, foot, and mouth disease (HFMD) is the most common communicable disease in China. Shandong Province is one of the most seriously affected areas. The distribution of HFMD had spatial heterogeneity and seasonal characteristic in this setting. The aim of this study was to explore the associations between climate and HFMD by a Bayesian approach from spatio-temporal interactions perspective. Methods The HFMD data of Shandong Province during 2008–2012 were derived from the China National Disease Surveillance Reporting and Management System. And six climatic indicators were obtained from the Meteorological Bureau of Shandong Province. The global spatial autocorrelation statistic (Moran’s I) was used to detect the spatial autocorrelation of HFMD cases in each year. The optimal one among four Bayesian models was further adopted to estimate the relative risk of the occurrence of HFMD via Markov chain Monte Carlo. Results The annual average incidence rate of HFMD was 104.40 per 100,000 in Shandong Province. Positive spatial autocorrelation appeared at county level (Moran’s I ≥0.30, P < 0.001). The best fitting Spatio-temporal interactive model showed that annual average temperature, annual average pressure, annual average relative humidity, annual average wind speed and annual sunshine hours were significantly positive related to the occurrence of HFMD. The estimated relative risk of 36, 87, 91, 79, 65 out of 140 counties for 2008–2012 respectively were significantly more than 1. Conclusions There were obvious spatio-temporal heterogeneity of HFMD in Shandong Province, and the climatic indicators were associated with the epidemic of HFMD. Bayesian approach should be recommended to capture the spatial-temporal pattern of HFMD.


Background
Hand, foot, and mouth disease (HFMD) is a common communicable disease usually affecting children, particularly those aged 5 years and younger [1]. It is most frequently caused by Coxsackievirus A16 and Enterovirus 71 [1,2], and is often characterized by a distinct clinical presentation of fever, or vesicular exanthema on their hands, feet, mouths, or buttocks [3,4].
Over the last decades, many large-scale outbreaks of HFMD were reported in East and Southeast Asia and have caused major public health concerns worldwide, especially in the affected countries [5][6][7][8]. In China, several outbreaks have also been reported, such as Linyi in Shandong (2007) [9], Fuyang in Anhui (2008) [10], Shanghai (2009) [11], Nanchang in Jiangxi (2010) [12], etc. Since its beginning in Fuyang, large-scale epidemic have resulted in the deaths of many children, HFMD was classified as a class C notifiable infectious disease by the Ministry of Health in China on 2 May 2008 [13]. To date, HFMD is still one of the leading causes of child death and a serious public health issue in China.
Currently, there are still no available effective vaccines or antiviral treatments specifically for HFMD. Thus, it is quite important to identify the possible risk factors for HFMD. Many epidemiological studies have been conducted to explore individual-level risk factors (e.g., demographic factors, socio-economic determinants and behavioral factors) [14][15][16][17][18][19][20]. Although these individual factors may play a role in HFMD incidence, group-level factors are more important for public health intervention. Considering the spatial distribution characteristics of HFMD, several researches (i.e., ecological studies) have been conducted at group-level (e.g., county, community, or city) based on spatial analysis methods (e.g., spatial regression model and spatial paneled model), which was helpful to understand the emerging trend of HFMD in certain area and explore the ecological causes of disease incidence [21], and can guide us to take the correct and timely public health interventions to prevent the outbreak. Among these, some studies found that the incidence of HFMD was associated with some climatic indicators (e.g., average temperature, relative humidity and annual sunshine hours) [22][23][24][25][26]. Therefore, the effects of climatic indicators on HFMD should be paid more attention on, especially in the context of climate change. In addition, the occurrence of HFMD presents significant seasonality, i.e., temporal characteristic [23,[27][28][29][30][31]. Therefore, spatio-temporal model should be preferred, which can not only help to recognize the spatial and temporal trend of disease, but also make prediction and guide us to formulate and implement appropriate regional public health intervention strategies to prevent and control this disease.
However, few studies from spatio-temporal interactions perspective, have attempted to explore this relationship and to assess how climate affect HFMD incidence. On the other hand, Bayesian methods can take into account possible correlations and covariates' effects and fully utilize the related information of disease and prior knowledge, it has been extensively used in many studies and has been recognized as a powerful means to provide more robust estimates [32][33][34][35].
Thus, the goal of this study is to investigate the relationships between HFMD incidence and climatic indicators based on Bayesian approaches. Using the data of new cases of HFMD reported during 2008-2012 in Shandong Province, four Bayesian models were constructed and compared, the best fitting one was further adopted to estimate the effects of six key climatic indicators, with the aim of increasing our knowledge of the true underlying geographic distribution of HFMD rates.

Study area
We performed this ecological spatial study of HFMD in Shandong Province (Figure 1), a coastal province in Eastern China with a population of approximately 97.33 million in 2013 [36]. Its area of 156,700 square kilometers was divided into 17 municipal districts which including 140 subdistricts (counties). Considering the county administrative level was often used for the HFMD-decision making process in China [27], county was used as the spatial unit of analysis ( Figure 1).

HFMD data
Data concerning to the reported cases of HFMD were provided by the Shandong Center for Disease Control and Prevention collected from the China Information System for Disease Control and Prevention (CISDCP), which included information about sex, age, occupation, and regional distribution. The diagnosis of HFMD was according to the clinical criteria from the HFMD Control and Prevention Guide offered by the Chinese Ministry of Health [37]. The report cards of HFMD cases filled out by professional doctors were collected by trained reporter and then input into the CISDCP within 24 hours based on the P. R. China infectious disease prevention and cure statute [27]. The data collected in this paper contained the basic demographic and incidence information for 497,876 HFMD cases during 2008 to 2012 in Shandong Province.

Demographic and climatic data
The climatic indicators included annual average temperature, annual average pressure, annual average relative humidity, annual average wind speed, annual sunshine hours and annual precipitation, which were obtained from the Meteorological Bureau of Shandong Province. The corresponding demographic data of each county over the study period referred to Shandong Statistical Yearbook [38][39][40][41][42].

Ethics statement
HFMD data were provided by Shandong CDC extracted from annual reports, which provided summarized count data of patients reported by county and month. All data were anonymized.

Statistical analysis
The frequencies of HFMD were summarized monthly and annually by geographic area (i.e., county), and the incidence rate of HFMD (per 100,000 population) in each county was calculated by HFMD counts divided by the corresponding population. We firstly analyzed the monthly temporal distribution characteristic of HFMD over the study period. The autocorrelation statistic (Moran's I) [43] was then used to detect the global spatial autocorrelation of HFMD cases. The software GeoDa™ 0.9.5-i was used to conduct the analysis.
Considering the epidemic characteristics of HFMD, a Bayesian approach [32][33][34][35]44] was utilized to explore the relationship between HFMD and climatic indicators. For each county i(i = 1, 2, ⋯, n) and each year t, the expected number of HFMD cases (ε it ) was estimated as the provincial overall mean rate in the year t multiplied by the population of the county (ε it = meanrate t × pop it ). The relative risk (RR) λ it was calculated as the observed number of HFMD cases divided by the expected. The observed number of cases (y it ) in the n counties could be treated as one realization of Poisson random variables, with means μ it , i.e., y it~P ois(μ it ), where μ it = ε it × λ it . λ it is a function of the effects of k covariates (x kit ) as well as spatial and temporal random effects.
We constructed four Bayesian models with climatic variables. The first model only included non-spatial random effects (Non-spatial model), Where β 0 represented the mean incidence over all counties and time periods with a flat distribution; β 1 , β 2 ,…, β 6 were the coefficients of climatic variables following non-informative Gaussian priors with zero mean and precision equal to 10 −5 . A nonspatially structured random effect ν i was included to account for extra-Poisson variation due to non-measured important covariates, with independent Gaussian distribution having zero mean and variance σ 2 v . We defined neighborhood as adjacent counties with simple binary adjacency weights, i.e., w ij = 1 if area i and j share a common boundary and w ij = 0 otherwise [44]. The second model considered spatial and non-spatial random effects (Spatial model), Where θ i represented the spatially structured random effects, which account for the spatial dependence with the prior distribution taken as a conditional intrinsic Gaussian autoregressive model. In the third model, spatial, non-spatial and temporal random effects were included (Spatio-temporal model), Where ω t was a random term representing between-year variation and was assumed to be an autoregressive process. This model assumed spatial correlation between the counties was independent of the year. The fourth model added spatio-temporal interactive effects on the basis of spatio-temporal model (Spatio-temporal interactive model), λ it = exp(β 0 + β 1 x 1it + β 2 x 2it + β 3 x 3it + β 4 x 4it + β 5 x 5it + β 6 x 6it + v i + θ i + ω i + psi it ). This model included a different set of random spatial effects for each year following a conditional autoregressive model. We adopted Markov chain Monte Carlo method to estimate the parameters by the public domain software package WinBUGS 1.4.3 software (Imperial College and MRC, London, UK). Totally 100,000 samples for each parameter of interested were generated, with a burn-in of 20,000 iterations to avoid the influence of the initial values. The convergence process was evaluated based on whether iteration trace and iteration history becoming stable, autocorrelation function being close to zero quickly. The Deviance Information Criterion (DIC) was used to assess the goodness-of -fit of models, smaller values of DIC indicated a more appropriate model [25].  Figure 2 displayed the monthly distribution of HFMD cases during the study period, which presented significant seasonality. Obviously, the incidence peak appeared between April and August, which accounted for 83.00% of all reported cases.

Spatial autocorrelation of HFMD cases
The results of the spatial autocorrelation test were listed in Table 2, demonstrating that high global spatial autocorrelation of HFMD existed at county level in Shandong Province within each epidemic year during 2008 to 2012 (Moran's I ≥0.30, P < 0.001).

Relationship between HFMD incidence and climatic indicators
On the basis of DIC (Table 3), Spatio-temporal model and Spatial model were quite similar and superior to Non-spatial model, indicating that spatial heterogeneity existed, while temporal effect can be ignored. However, Spatio-temporal interactive model including spatiotemporal interaction had the lowest DIC. Thus, further analysis should be focused on Spatio-temporal interactive model.
In Table 4, the parameter estimates for association and the corresponding confidence interval (CI) were presented. The results indicated that all variables except annual precipitation were significantly positive related to the risk of HFMD. Among them, annual average pressure was most important with parameter estimate 0.1054 (95% CI: 0.1014 to 0.1078), followed by annual average temperature with parameter estimate 0.0956 (95% CI: 0.08272 to 0.1076).

Discussion
Shandong Province had been one of the most serious HFMD epidemic areas in China, with annual average incidence rate 104.40 per 100,000 during 2008 to 2012. Most HFMD cases (95.04%) were aged less than 6 years old, with an average male-to-female sex ratio 1.69 (see Table 1). Single seasonal peaks except double peaks in 2009, could be found between April and August during the study period (see Figure 2), which were different from other districts of China. For example, in Jiangsu Province, double peaks (the highest occurrence between April and June and the second occurring in November) appeared [45]; in Hong Kong, warmer seasonal peak (May-July) and winter peak (October-December) were detected [17]; in Guangdong Province, HFMD incidence peaked in April-May and September-October [19]. These difference might be partly due to some risk factors, such as climatic, geographic and social factors [15,19]. In addition, spatial autocorrelation test indicated that HFMD had positive spatial autocorrelation at county level in Shandong Province (Table 2). This also confirmed that the occurrence of HFMD might be closely related to certain ecological factors in the specific area.   Considering the goodness-of-fit by the DIC, the model including the spatial and temporal interactive effects was best fitted (DIC =6938.09, Table 3), which demonstrated a strong spatio-temporal heterogeneity in HFMD risk at the county scale, with clusters of high risk areas, as previously reported [27]. And the HFMD incidence rate was correlated with annual average temperature, annual average pressure, annual average relative humidity, annual average wind speed, and annual sunshine hours (Table 4), which were partly compatible with previous ecological analyses from the literature [20,22,46], inferring that HFMD continued to be a disease related to climate. In this study, much more precise estimates of RR were obtained from Spatio-temporal interactive model, confirming the hypothesis that the risk of the disease was related to what was occurring in the neighbourhood of each spatial unit. Figure 3 clearly displayed the spatiotemporal characteristic of RR, the epidemic of HFMD extended continuously over time, and the spatial pattern changed from concentrating on minority counties to expanding to more areas (especially in 2009 and 2010). Also, many high risk counties were located in areas for which HFMD clusters have been detected [27]. Certainly, these climatic factors are hard to change, while adaptive and objective prevention and control mechanism can be established according to the active detected possible impacts, e,g., advanced prevention for high risk areas.
Bayesian spatio-temporal interactive model was a valuable tool for the spatial and temporal interactive assessment of disease patterns that could help to identify county differences, and explore possible risk factors simultaneously. This method could solve most of the problems faced by traditional statistical methods, such as the spatial autocorrelation and the potential dependence between the covariates [47]. In previous studies, Scan statistics methods had been used to determine the spatial or spatio-temporal distribution of HFMD [23,[27][28][29][30]. These approaches might absorb the surrounding regions and generate false-positives areas due to a lack of specificity, that is, some non-cluster areas might be determined to be clustered since they encompassed many neighbourhoods and tended to detect larger clusters than expected [48,49]. However, the scan statistics could be complementary to Bayesian method. That was to say, the scan statistic detected general regions in which the risk was significantly high and the Bayesian posterior distribution further helped to identify the neighbourhoods contributing strongly to the scan statistic circle [50]. Thus, results for cluster analysis should be interpreted with knowledge of the spatial rate distribution, such as spatial Bayesian rates in particular [51]. Furthermore, the analysis unit (i.e., county) was also considered adequately given the disease and risk factors information that was available and the spatial level at which policies were taken. However, one must note that ecologic bias was inevitable in any ecological study [52]. In addition, some results might be biased due to the artificial grouping of observations and variables at the county level. Despite these limitations, this approach was very useful to explore spatially aggregated data and to highlight the most risky areas to conduct more accurate analysis.

Conclusions
Using a Bayesian approach to estimate the contribution of climatic indicators on the spatial-temporal pattern of HFMD should be encouraged in epidemiology. Our results confirmed the spatial-temporal heterogeneity of HFMD distribution, with high risk in particular areas observed in Shandong Province, and the importance of climatic covariate. The results may help public health authorities to set up priorities regarding to be targeted for prevention or control measures. Submit your next manuscript to BioMed Central and take full advantage of: