Spatiotemporal Distribution of Hand, Foot, and Mouth Disease in Guangdong Province, China and Potential Predictors, 2009–2012

Background: Hand, foot, and mouth disease (HFMD) is a common infectious disease among children. Guangdong Province is one of the most severely affected provinces in south China. This study aims to identify the spatiotemporal distribution characteristics and potential predictors of HFMD in Guangdong Province and provide a theoretical basis for the disease control and prevention. Methods: Case-based HFMD surveillance data from 2009 to 2012 was obtained from the China Center for Disease Control and Prevention (China CDC). The Bayesian spatiotemporal model was used to evaluate the spatiotemporal variations of HFMD and identify the potential association with meteorological and socioeconomic factors. Results: Spatially, areas with higher relative risk (RR) of HFMD tended to be clustered around the Pearl River Delta region (the mid-east of the province). Temporally, we observed that the risk of HFMD peaked from April to July and October to December each year and detected an upward trend between 2009 and 2012. There was positive nonlinear enhancement between spatial and temporal effects, and the distribution of relative risk in space was not fixed, which had an irregular fluctuating trend in each month. The risk of HFMD was significantly associated with monthly average relative humidity (RR: 1.015, 95% CI: 1.006–1.024), monthly average temperature (RR: 1.045, 95% CI: 1.021–1.069), and monthly average rainfall (RR: 1.004, 95% CI: 1.001–1.008), but not significantly associated with average GDP. Conclusions: The risk of HFMD in Guangdong showed significant spatiotemporal heterogeneity. There was spatiotemporal interaction in the relative risk of HFMD. Adding a spatiotemporal interaction term could well explain the change of spatial effect with time, thus increasing the goodness of fit of the model. Meteorological factors, such as monthly average relative humidity, monthly average temperature, and monthly average rainfall, might be the driving factors of HFMD.

. The four parts of Guangdong Province and cumulative incidence of hand, foot, and mouth disease (HFMD), 2009-2012.

Surveillance Data of HFMD
Case-based HFMD surveillance data from 2009 to 2012 was obtained from the China Center for Disease Control and Prevention (China CDC). It was confidential data; raw data cannot be published due to state law and ethics. The clinical criteria for diagnosis of HFMD cases was provided in a guidebook published by the Ministry of Health of China in 2008 [22], in which patients were defined as HFMD with the occurrence of the following symptoms: fever, papules and herpetic lesions on the hands or feet, rashes on the buttocks or knees, inflammatory flushing around the rashes and little fluid in the blisters, sparse herpetic lesions on oral mucosa.

Potential Predictors Data
Monthly average relative humidity, the monthly average temperature, monthly average rainfall, monthly average wind speed, and monthly sunshine duration data were obtained from the China Meteorological Data Sharing Service System (http://data.cma.cn/). The monthly county-level meteorological variables were estimated using ordinary spatial kriging methods based on 26 meteorological surveillance stations within Guangdong Province. The 26 meteorological surveillance stations are mapped in Figure 2. In addition, we obtained population and socioeconomic data, i.e., average GDP, from the statistical yearbook.

Surveillance Data of HFMD
Case-based HFMD surveillance data from 2009 to 2012 was obtained from the China Center for Disease Control and Prevention (China CDC). It was confidential data; raw data cannot be published due to state law and ethics. The clinical criteria for diagnosis of HFMD cases was provided in a guidebook published by the Ministry of Health of China in 2008 [22], in which patients were defined as HFMD with the occurrence of the following symptoms: fever, papules and herpetic lesions on the hands or feet, rashes on the buttocks or knees, inflammatory flushing around the rashes and little fluid in the blisters, sparse herpetic lesions on oral mucosa.

Potential Predictors Data
Monthly average relative humidity, the monthly average temperature, monthly average rainfall, monthly average wind speed, and monthly sunshine duration data were obtained from the China Meteorological Data Sharing Service System (http://data.cma.cn/). The monthly county-level meteorological variables were estimated using ordinary spatial kriging methods based on 26 meteorological surveillance stations within Guangdong Province. The 26 meteorological surveillance

Bayesian Spatiotemporal Model
The Bayesian spatiotemporal model was used to analyze the relative risk (RR) of HFMD from 2009 to 2012. The form of the model was specified as the following: represented the expected number of HFMD cases, which was calculated as the product of overall incidence of the province and the population for each county during the study period, andθ ij was the RR. 0 b was the intercept, pij X was the th p meteorological or socioeconomic variable in region i , month j , β p was the regression coefficient. We calculated the variance inflation factor (VIF) for all candidate variables to assess the multicollinearity; variables with VIF less than 5 were selected for inclusion in the model [23]. The Besag, York, and Mollie (BYM) model was used to model the common spatial component [24], which consisted of two components: spatially structured random effect μ i and spatially unstructured random effect υ i . μ i was commonly assumed to follow a conditional autoregressive prior structure (CAR) in infectious diseases spatial epidemiology [25][26][27], which showed that county i had a similar pattern of disease incidence with the adjacent counties.υ i was assumed to be a normal distribution, representing that county i had an independent pattern of disease incidence from the adjacent counties. If a disease has a strong spatial autocorrelation, the spatially structured random effect μ i will have more contribution explaining the spatial variance of

Bayesian Spatiotemporal Model
The Bayesian spatiotemporal model was used to analyze the relative risk (RR) of HFMD from 2009 to 2012. The form of the model was specified as the following: where Y ij was the number of reported HFMD cases in region i, month j (j = 2009-01, . . . , 2012-12). The model assumed that Y ij was subject to a Poisson distribution with a mean λ ij = e ij × θ ij , where e ij represented the expected number of HFMD cases, which was calculated as the product of overall incidence of the province and the population for each county during the study period, and θ ij was the RR. b 0 was the intercept, X pij was the p th meteorological or socioeconomic variable in region i, month j, β p was the regression coefficient. We calculated the variance inflation factor (VIF) for all candidate variables to assess the multicollinearity; variables with VIF less than 5 were selected for inclusion in the model [23]. The Besag, York, and Mollie (BYM) model was used to model the common spatial component [24], which consisted of two components: spatially structured random effect µ i and spatially unstructured random effect υ i . µ i was commonly assumed to follow a conditional autoregressive prior structure (CAR) in infectious diseases spatial epidemiology [25][26][27], which showed that county i had a similar pattern of disease incidence with the adjacent counties. υ i was assumed to be a normal distribution, representing that county i had an independent pattern of disease incidence from the adjacent counties. If a disease has a strong spatial autocorrelation, the spatially structured random effect µ i will have more contribution explaining the spatial variance of the disease risk; otherwise, the component of the spatially unstructured random effect υ i will have more explanatory power. The temporal component was also consisted of two components: monthly structured random effect γ j and monthly unstructured random effect ϕ j . γ j and ϕ j were assumed to follow a 1 order random walk and a normal distribution, respectively [21,28,29]. δ ij was the spatiotemporal interaction term. The structure matrix can be written as the Kronecker product of R δ = R υ ⊗ R ϕ , which assumed that the spatially unstructured random effect υ i and the monthly unstructured random effect ϕ j interact with each other [21].
We developed five different models. Model 1 only included potential predictors. Model 2 and model 3 added a spatial component and a temporal component to model 1, respectively. Model 4 was the combination of model 2 and model 3. Finally, model 5 included all the components, including potential predictors, spatial, temporal, and spatiotemporal interaction components.
The models were fitted using the integrated nested Laplace approximation (INLA) method incorporated in the R-INLA package (R version 3.4.3., R Core Team, Vienna, Austria). The model with the smallest deviance information criterion (DIC) was considered the optimal and used for our study.

Geographical Detectors
In this study, the q-statistic was used to quantify the spatiotemporal heterogeneity of HFMD and detect the interaction relationship between spatial effect and temporal effect [30,31]. The form of the q-statistic was specified as the following: where q denotes the level of spatiotemporal heterogeneity; the value of the statistic is required to be within [0, 1]. If the value approaches 1, it indicates a strong heterogeneity and if the value approaches 0, it indicates a random distribution. N was the number of all units, which could be divided into L stratums. Stratum h was composed of N h units. σ 2 and σ 2 h were the variance over all the units and within stratum h (h = 1, . . . , L), respectively.
We can divide the units of observation into two stratums, space L S (L S = 1, . . . , 21) and time L T (L T = 1, . . . , 48). If q(L S ∩ L T ) > q(L S ) + q(L T ), there is a nonlinear enhancement of space and time, if q(L S ∩ L T ) < Min(q(L S ), q(L T )), there is a nonlinear weakening of them [31].

Descriptive Statistics
Overall, 911,640 HFMD cases were reported in Guangdong Province from 2009 to 2012, and the annual average incidences from 2009 to 2012 were 10.36/10,000, 24.52/10,000, 28.21/10,000, and 35.26/10,000, respectively. Spatial distribution of cumulative incidence of HFMD at county level during the study period is shown in Figure 1. Counties with a higher HFMD incidence were clustered in the Pearl River Delta region. A descriptive summary of the meteorological and socioeconomic variables is shown in Table 1.  Table 2 shows the results of the multicollinearity analysis. The VIF values of all variables were less than 5, so it could be considered that there was no multicollinearity between the variables, and all variables were included in the model.  Table 3 summarizes the DIC values for the five types of models. Model 5 was selected as the final model as it had the lowest DIC value. Further results were all based on model 5. Table 3. The results of model selection. DIC: deviance information criterion.

Model
Component DIC

Spatial Distribution
The estimated RR values of the spatial effect (µ i + υ i ) are mapped in Figure 3. The spatial distribution of HFMD risk exhibited explicit spatial heterogeneity, as the q-statistic was 0.669 (p < 0.001).
The RR values were relatively high in the Pearl River Delta region of Guangdong Province. This finding implied that the counties in the Pearl River Delta region had a relatively higher HFMD risk.
If a county had a stable higher or lower risk compared with the overall level in the region, it was classified as a hot spot or a cold spot, respectively [32]. This classification can make the spatial distribution of RR more intuitive. Hot and cold spots can be calculated according to the posterior probability of spatial effect-if the posterior probability P(exp(µ i + υ i ) > 1|data) > 0.8 , then the county i was classified as a hot spot. Similarly, if the posterior probability <0.2, then the county i was classified as a cold spot. Figure 3 shows the distribution of hot and cold spots of HFMD in Guangdong Province. There were 53 counties that accounted for 43% of all counties, and they were classified as hot spots. They were mainly distributed in the Pearl River Delta region. The cold spot regions were predominantly clustered in the east-west coast and northern mountainous counties (57 counties), accounting for 46% of all counties.  If a county had a stable higher or lower risk compared with the overall level in the region, it was classified as a hot spot or a cold spot, respectively [32]. This classification can make the spatial distribution of RR more intuitive.  Figure 3 shows the distribution of hot and cold spots of HFMD in Guangdong Province. There were 53 counties that accounted for 43% of all counties, and they were classified as hot spots. They were mainly distributed in the Pearl River Delta region. The cold spot regions were predominantly clustered in the east-west coast and northern mountainous counties (57 counties), accounting for 46% of all counties.

Temporal Distribution
The yearly and monthly numbers of cases are shown in Figure 4a,b, respectively. It was noteworthy that the incidence of HFMD had an increasing trend. Besides, seasonal peaks were detected around May and October.

Temporal Distribution
The yearly and monthly numbers of cases are shown in Figure 4a,b, respectively. It was noteworthy that the incidence of HFMD had an increasing trend. Besides, seasonal peaks were detected around May and October.
The RR of monthly effect (γ j + ϕ j ) is shown in Figure 4c. We detected an upward trend in RR between 2009 and 2012. The RR of monthly effect presented significant seasonality. Semiannual peaks were observed during the study period. The peaks occurred from April to July, and October to December in each year, with the first peak higher than the second.

Spatiotemporal Interaction Patterns
The results of interaction detection showed that there was positive nonlinear enhancement between spatial and temporal effects. q(L S ∩ L T )= 0.759 was greater than the sum of q(L S ) = 0.373 and q(L T ) = 0.132. However, in the previous studies, spatial effect was assumed to be fixed and did not change over time [13,16], i.e., spatial and temporal effects were independent of each other. Actually, the distribution of RR in space was not the same in each month, and there were certain changes. The spatiotemporal interaction component δ ij could explain the differences in the time trend of RR for different areas. The results are listed in Figure 5. From the figure, we can see the variations of spatial effect during each month. The results of heterogeneity detection are listed in Table 4. There were 33 months that had no spatial heterogeneity in the distribution of the changed RR. Therefore, most of the variations of spatial effect during each month were irregular. The RR of monthly effect ( γ ϕ + j j ) is shown in Figure 4c. We detected an upward trend in RR between 2009 and 2012. The RR of monthly effect presented significant seasonality. Semiannual peaks were observed during the study period. The peaks occurred from April to July, and October to December in each year, with the first peak higher than the second.

Spatiotemporal Interaction Patterns
The results of interaction detection showed that there was positive nonlinear enhancement between spatial and temporal effects. . However, in the previous studies, spatial effect was assumed to be fixed and did not change over time [13,16], i.e., spatial and temporal effects were independent of each other. Actually, the distribution of RR in space was not the same in each month, and there were certain changes. The spatiotemporal interaction component δ ij could explain the differences in the time trend of RR for different areas. The results are listed in Figure 5. From the figure, we can see the variations of spatial effect during each month. The results of heterogeneity detection are listed in Table 4. There were 33 months that had no spatial heterogeneity in the distribution of the changed RR. Therefore, most of the variations of spatial effect during each month were irregular.

Potential Predictors
The estimated results of regression coefficients of meteorological variables and socioeconomic variables ( β p ) are shown in Table 5. There was a positive association between meteorological variables and relative risk. The increasing of one unit in the monthly average relative humidity, monthly average temperature, and monthly average rainfall was associated with the increase of 1.5%

Potential Predictors
The estimated results of regression coefficients of meteorological variables and socioeconomic variables (β p ) are shown in Table 5. There was a positive association between meteorological variables and relative risk. The increasing of one unit in the monthly average relative humidity, monthly average temperature, and monthly average rainfall was associated with the increase of 1.5% (95% CI: 0.6%-2.4%), 4.5% (95% CI: 2.1%-6.9%), and 0.4% (95% CI: 0.1%-0.8%) in the RR of HFMD, respectively. Among these variables, the monthly average temperature had the greatest impact on the HFMD. However, monthly sunshine duration, monthly average wind speed, and average GDP were not significantly associated with the risk of HFMD in this study.  Figure 6 shows the county-level estimated RR values of associated potential predictors, which were the product terms of the regression coefficients of potential predictors and the average of the standardization values of potential predictors in each region. They represented the RR of HFMD caused by potential predictors to each region. The western coast of Guangdong Province had a higher average relative humidity level, resulting in a higher RR, as shown in Figure 6a. The southern regions with higher temperatures and the southwestern regions with higher rainfall also had greater RR, as shown in Figure 6b,c.  Figure 6 shows the county-level estimated RR values of associated potential predictors, which were the product terms of the regression coefficients of potential predictors and the average of the standardization values of potential predictors in each region. They represented the RR of HFMD caused by potential predictors to each region. The western coast of Guangdong Province had a higher average relative humidity level, resulting in a higher RR, as shown in Figure 6a. The southern regions with higher temperatures and the southwestern regions with higher rainfall also had greater RR, as shown in Figure 6 b,c.

Discussion
The present study explored the spatiotemporal distribution of HFMD and its associations with potential predictors. The results revealed that the counties with high relative risk were mainly in the Pearl River Delta, which is the economic center of Guangdong. This finding was consistent with previous studies [3,17]. A study of Beijing-Tianjin-Hebei found that high risk regions were mainly located in large cities, such as Beijing, Tianjin, Shijiazhuang, and their neighboring regions [17]. In addition, a study in Chongqing found that cluster centers were in nine main urban districts [3].
During the four-year study period from 2009 to 2012, an increasing trend of relative risk was detected. HFMD showed semiannual seasonality in Guangdong, which was also commonly observed in other southern provinces and regions in China, such as Jiangsu (peaks period: April to June, and November to December) [33], Chongqing (April to July, and October to December) [3], and Hong Kong (May to July, and October to December) [15]. Meanwhile, a single peak has been shown to appear in northern China, such as Beijing (May to July) [13] and Shandong Province (April and August) [34]. Different climates and living habits can be potential reasons for the different peak

Discussion
The present study explored the spatiotemporal distribution of HFMD and its associations with potential predictors. The results revealed that the counties with high relative risk were mainly in the Pearl River Delta, which is the economic center of Guangdong. This finding was consistent with previous studies [3,17]. A study of Beijing-Tianjin-Hebei found that high risk regions were mainly located in large cities, such as Beijing, Tianjin, Shijiazhuang, and their neighboring regions [17]. In addition, a study in Chongqing found that cluster centers were in nine main urban districts [3].
During the four-year study period from 2009 to 2012, an increasing trend of relative risk was detected. HFMD showed semiannual seasonality in Guangdong, which was also commonly observed in other southern provinces and regions in China, such as Jiangsu (peaks period: April to June, and November to December) [33], Chongqing (April to July, and October to December) [3], and Hong Kong (May to July, and October to December) [15]. Meanwhile, a single peak has been shown to appear in northern China, such as Beijing (May to July) [13] and Shandong Province (April and August) [34]. Different climates and living habits can be potential reasons for the different peak months. For instance, the climate is warm and humid in southern China, while cold and dry in the north.
Temperature and relative humidity were associated with the relative risk of HFMD, with increased relative risk of the disease following higher temperature and relative humidity. This was consistent with the previous findings in Guangdong [10,35], Beijing [13], and Hong Kong [36]. Average rainfall was also an important influence on HFMD transmission. It was generally consistent with other studies conducted in Beijing and Henan [13,27]. According to previous studies, there are two potential ways that meteorological factors influence HFMD-one is by affecting the survival and reproductive capacity of enteroviruses [17], and another is by altering patterns of human behavior [37]. For example, Huang et al. found that in hot and humid summers, outdoor activities were reduced, and people tended to spend more time indoors in air-conditioned houses, thereby providing more opportunities for contact with each other [35].
The confounders were controlled by adding the spatial and temporal components, which had both effects on the number of reported HFMD cases and meteorological or socioeconomic variables. Spatiotemporal effects represented the residual caused by latent variables that were not included in the model, such as medical conditions and the capacity for HFMD prevention and control in different counties and different periods, which are easier to control than meteorological factors. Thus, areas and months with high relative risk should be paid more attention through the strategies including strengthening the education and supervision of kindergartens.
Major prior studies used traditional models such as spatial autocorrelation and spatiotemporal scan statistics [38]. Spatial autocorrelation cannot detect trends in the time, but just performs spatial analysis. Spatiotemporal scan statistics is a simple and easy method, and it has been widely used in public health. However, covariates cannot be included in the model and it is not suitable for a risk region with an irregular shape. A Bayesian spatiotemporal model is more flexible in modelling spatiotemporal components than the above models with some limitations.
This was a comprehensive study of HFMD in terms of spatiotemporal effects and meteorological variables. We incorporated spatiotemporal effects in the regression that controlled uncertainties resulting from residual confounding, which made the estimation of the impact of meteorological variables more accurate. In addition, the hot spots and high-risk time of HFMD were identified. Through the detection of spatiotemporal interaction, we found that the spatial effect of HFMD had positive nonlinear enhancement with the temporal effect. After adding the spatiotemporal interaction component to the model, the goodness of fit of the model was significantly increased. Compared with the model without the interaction term (model 4), the DIC value was reduced by 77.8% to 45,036.5. Therefore, adding the interaction component can improve the accuracy of the model-this has not been considered in many previous related studies.
The present study has some limitations. First, we used the annual average GDP for each county as the monthly data. Second, there were only 26 meteorological surveillance stations in Guangdong Province; ordinary kriging interpolation results might not cover all variations of the meteorological variables at the county level. However, this was more accurate than using the same meteorological data to all counties in each city. Third, we estimated spatiotemporal variations in HFMD at the scale of counties and months; we did not include the factors at an individual and pathogenic level, such as personal hygiene, educational background, incomes of the children's parents, living conditions, and composition of major pathogens. Further studies should consider the potential impacts of these factors.

Conclusions
The risk of HFMD in Guangdong showed significant spatial heterogeneity, with higher relative risk counties mainly gathered in the Pearl River Delta region. The relative risk of HFMD exhibited an increasing trend and semiannual seasonality from 2009 to 2012. High relative risk areas and months should be the focus of greater attention in the prevention and control of HFMD. There was spatiotemporal interaction between HFMD relative risk. Adding a spatiotemporal interaction term could well explain the change of spatial effect with time, thus increasing the goodness of fit of the model. Meteorological factors, such as monthly average relative humidity, monthly average temperature, and monthly average rainfall, might be the driving factors of HFMD.