Skip to main content
Advertisement
  • Loading metrics

Combined impacts of environmental and socioeconomic covariates on HFMD risk in China: A spatiotemporal heterogeneous perspective

  • Chun-Hu Li ,

    Contributed equally to this work with: Chun-Hu Li, Jun-Jie Mao, You-Jia Wu

    Roles Data curation, Formal analysis, Methodology, Visualization, Writing – original draft

    Affiliation Joint Division of Clinical Epidemiology, Affiliated Hospital of Nantong University, School of Public Health of Nantong University, Nantong, China

  • Jun-Jie Mao ,

    Contributed equally to this work with: Chun-Hu Li, Jun-Jie Mao, You-Jia Wu

    Roles Data curation, Formal analysis, Methodology, Validation, Writing – original draft

    Affiliation Joint Division of Clinical Epidemiology, Affiliated Hospital of Nantong University, School of Public Health of Nantong University, Nantong, China

  • You-Jia Wu ,

    Contributed equally to this work with: Chun-Hu Li, Jun-Jie Mao, You-Jia Wu

    Roles Data curation, Formal analysis

    Affiliation Department of Pediatrics, Affiliated Hospital of Nantong University, Nantong, China

  • Bin Zhang,

    Roles Formal analysis

    Affiliation Department of Infectious Diseases, Affiliated Hospital of Nantong University, Nantong, China

  • Xun Zhuang,

    Roles Formal analysis, Funding acquisition

    Affiliation Department of Epidemiology and Biostatistics, School of Public Health of Nantong University, Nantong, China

  • Gang Qin ,

    Roles Conceptualization, Funding acquisition, Writing – review & editing

    tonygqin@ntu.edu.cn (GQ); liu.hm@ntu.edu.cn (HML)

    Affiliations Joint Division of Clinical Epidemiology, Affiliated Hospital of Nantong University, School of Public Health of Nantong University, Nantong, China, Department of Infectious Diseases, Affiliated Hospital of Nantong University, Nantong, China

  • Hong-Mei Liu

    Roles Conceptualization, Funding acquisition, Writing – review & editing

    tonygqin@ntu.edu.cn (GQ); liu.hm@ntu.edu.cn (HML)

    Affiliation School of Transportation and Civil Engineering of Nantong University, Nantong, China

Abstract

Background

Understanding geospatial impacts of multi-sourced influencing factors on the epidemic of hand-foot-and-mouth disease (HFMD) is of great significance for formulating disease control policies tailored to regional-specific needs, yet the knowledge is very limited. We aim to identify and further quantify the spatiotemporal heterogeneous effects of environmental and socioeconomic factors on HFMD dynamics.

Methods

We collected monthly province-level HFMD incidence and related environmental and socioeconomic data in China during 2009–2018. Hierarchical Bayesian models were constructed to investigate the spatiotemporal relationships between regional HFMD and various covariates: linear and nonlinear effects for environmental covariates, and linear effects for socioeconomic covariates.

Results

The spatiotemporal distribution of HFMD cases was highly heterogeneous, indicated by the Lorenz curves and the corresponding Gini indices. The peak time (R2 = 0.65, P = 0.009), annual amplitude (R2 = 0.94, P<0.001), and semi-annual periodicity contribution (R2 = 0.88, P<0.001) displayed marked latitudinal gradients in Central China region. The most likely cluster areas for HFMD were located in south China (Guangdong, Guangxi, Hunan, Hainan) from April 2013 to October 2017. The Bayesian models achieved the best predictive performance (R2 = 0.87, P<0.001). We found significant nonlinear associations between monthly average temperature, relative humidity, normalized difference vegetation index and HFMD transmission. Besides, population density (RR = 1.261; 95%CI, 1.169–1.353), birth rate (RR = 1.058; 95%CI, 1.025–1.090), real GDP per capita (RR = 1.163; 95%CI, 1.033–1.310) and school vacation (RR = 0.507; 95%CI, 0.459–0.559) were identified to have positive or negative effects on HFMD respectively. Our model could successfully predict months with HFMD outbreaks versus non-outbreaks in provinces of China from Jan 2009 to Dec 2018.

Conclusions

Our study highlights the importance of refined spatial and temporal data, as well as environmental and socioeconomic information, on HFMD transmission dynamics. The spatiotemporal analysis framework may provide insights into adjusting regional interventions to local conditions and temporal variations in broader natural and social sciences.

Author summary

Hand-foot-and-mouth disease (HFMD) is one of the most common pediatric infectious diseases worldwide and it has affected up to two million children annually in China since 2008. This study used Bayesian spatiotemporal modeling to analyze provincial-level HFMD cases, environmental and socioeconomic data during 2009–2018. We found that both environmental (temperature, relative humidity and normalized difference vegetation index) and socioeconomic covariates (population density, birth rate, real GDP per capita and school vacation) are associated with HFMD occurrence under spatiotemporal scales. This model may help to (1) identify outbreaks of HFMD, including their timing, location and severity; (2) understand the transmission dynamics by mapping the spread of HFMD over time and space; (3) indicate areas of higher risk or the emergence of new outbreaks by detecting clusters of HFMD cases in specific geographic locations. Overall, this study provides valuable information for public health decision-makers and researchers, allowing them to better understand the spread of the disease and target their efforts more efficiently to control its transmission and prevent outbreaks.

Introduction

Globally, hand-foot-and-mouth disease (HFMD) is one of the most common pediatric infectious diseases, and remains a major public health concern across the Asia-Pacific region pandemic potential [1]. In China, since its first large-scale epidemic in Fuyang City, Anhui Province in 2008 [2], HFMD has affected up to two million children annually [3,4]. Children younger than 5 years are especially prone to HFMD. Although most cases show mild and self-limiting illness, HFMD has the highest mortality rate among class ‘C’ infectious diseases in China, due to severe neurological and pulmonary complications [4]. HFMD is mainly caused by infection of human enterovirus (EV) A species. Since 2013 in China, coxsackievirus-A6 (CV-A6) and CV-A10 have become the predominant human EV serotypes over enterovirus A71 (EV-A71) and CV-A16 [5,6]. It was consistent with the worldwide (including Europe) outbreak of CV-A6 in 2013 [7]. EV-A71 vaccine has been available since 2016 and is expected to reduce EV71-related HFMD cases gradually. Low vaccine acceptance (< 50%) among the parents may be a hurdle for general EV71 vaccination among young children [8]. Moreover, for the full control of HFMD caused by several viral serotypes, development of multivalent vaccines is urgently needed [9].

Significant geographic disparities among China’s Eastern, Central, and Western regions pose challenges to disease control policies [10]. Previous studies suggested the spatial heterogeneity of HFMD and some possible influencing factors such as environmental and socioeconomic variables [11]. Some studies found that climate factors (i.e. temperature, relative humidity and wind speed) may be associated with HFMD transmission and socioeconomic factors (i.e. gross domestic product [GDP], ratio of urban to rural population and birth population) may be associated with HFMD transmission [1216]. However, most studies were performed at a city scale or a provincial scale, and did not include analyses of socioeconomic factors. Moreover, spatiotemporal analyses at the national scale are scarce.

Given the high endemicity and dynamic complexity of HFMD, our study had three main objectives: (1) to detect and quantify the spatial heterogeneity of HFMD based on monthly data for 10 years from 31 provinces in China; (2) to construct the Bayesian spatiotemporal model to estimate the HFMD risk and evaluate model performance; and (3) to identify the driving factors on HFMD and quantify their effects in different regions.

Materials and methods

Data sources

Monthly data on HFMD cases from 31 provinces in China during 2009–2018 were collected from public health science data center website (www.phsciencedata.cn) of Chinese Center for Disease Control and Prevention (CCDC) [17]. The meteorological data were retrieved from China Surface Climatic Data Daily Data Set (V3.0) by China Meteorological Data Network (data.cma.cn/en), consisted of the observation data from 824 meteorological stations nationwide [18]. The socioeconomic data were obtained from the National Bureau of Statistics of China (www.stats.gov.cn) [19]. The normalized difference vegetation index (NDVI, range from -1 to 1) were retrieved from satellite images provided in the Moderate Resolution Imaging Spectroradiometer (MODIS) product MOD13Q1 (https://giovanni.gsfc.nasa.gov/)

Temporal analyses

Seasonal analysis of HFMD incidence used X-12 Autoregressive Integrated Moving Average model (ARIMA), which was the most widely used method for decomposing time series into a trend, seasonal, and residual components [20]. Morlet wavelet analysis was performed to detect and quantify the periodicity of the time series of HFMD incidence [21]. Seasonal index was calculated as ratio of the average incidence rate to the average incidence rate per month. Heat maps of monthly HFMD incidence were created to identify the peak seasons in different provinces. A seasonal multiple linear regression model was adopted, using harmonic terms representing annual and semiannual periodic cycles, to estimate the peak timing and amplitude of the annual and semi-annual periodicities of HFMD activity in each province [3].

Spatial analyses

The spatial heterogeneity of HFMD was judged of by plotting Lorenz curves and calculating the corresponding Gini indices [22]. The cumulative percentage of population and the cumulative percentage of HFMD cases were plotted graphically, reflecting the imbalance degree in the cases distribution. The Gini index (0–1) summarizes the statistics produced by the Lorentz curve. The value of Gini index close to 1 indicates a high degree of spatial heterogeneity.

The global Moran’s I index was used to measure global spatial autocorrelation based on both HFMD locations and incidence values simultaneously [23]. The value of Moran’s I index ranges from -1 to 1. If the index is significantly greater than 0, it indicates spatial agglomeration among provinces. A negative value indicates spatial difference.

We also perform the hot spot analysis, using ArcGIS Pro version 2.5 (ESRI Inc., Redlands, CA, USA), to identify local spatial autocorrelation [24]. Specifically, Getis-Ord Gi* statistic was calculated for local odds ratio (OR) maps (i.e., hot spot maps) of HFMD incidence in this study.

Space-time scan analyses

The Kulldorff’s retrospective space-time scan statistic [25], based on the discrete Poisson probability model, was applied to detect the geographic areas and time periods of potential clusters with a significantly higher incidence of HFMD than those of neighboring areas. The shape of space-time scanning windows (clusters) is cylinder with radius (geographic) and height (time). The relative risk (RR) and log-likelihood ratio (LLR) were computed. The maximum radius of the circular base was set to 50% of the total population at risk and the maximum height of the cylinder was set to 50% of the total study period. The analyses were performed with SaTScan version 10.1 (Kulldorff and Information Management Services Inc, Cambridge, MA, USA). Then, ArcGIS was used to visualize the relative risk (RR) of HFMD in high-risk cluster areas.

Spatiotemporal modeling and analyses

All the variables collected in this study were initially screened by the multicollinearity test using variance inflation factor (VIF) [26]. If VIF value of the variable > 10, it would be rejected. Bayesian model averaging (BMA) approach was applied for further variable selection [27]. It may identify the best-performing model, based on posterior probability as the weight to average all plausible models considered.

To explore the relationship of environmental, socioeconomic covariates and HFMD risk, we considered four types of regression models that responded to monthly HFMD cases for 31 provinces in China over 2009–18. Model 1 was the ordinary multiple linear regression (MLR) [28]. Model 2 was distributed lag nonlinear model (DLNM) [11]. Model 3 was the hierarchical Bayesian spatiotemporal model (HBSTM) adjusted for the temporal and spatial confounding effects [29]. We used a penalty complexity priors approach for all temporal and spatial random effects hyperparameters [30]. Model 4 was a combination of HBSTM and DLNM to making further interpretations of nonlinear effects of environmental covariates [31]. For each model, a negative binomial distribution was used to account for over-dispersion of HFMD case counts [26]. The goodness of fit (GoF) was assessed by the Watanabe-Akaike information criterion (WAIC), and the deviance information criterion (DIC), for which smaller values indicate better model fit [32]. The analyses were performed using R software version 4.2.1 (R Foundation for Statistical Computing, Vienna, Austria), mainly with the packages “bma”, “inla” and “dlnm”.

Cross-validation

An essential step in establishing realistic epidemic models is cross-validating them against historical spread of outbreaks and underlying patterns [33]. The out-of-sample predicted HFMD incidence rates from Jan 2009 to Dec 2018, as the posterior means and upper 95% credible intervals, were calculated from the final HBSTM+DLNM model (fitted 10 × 12 months, excluding the month for which the prediction is valid)). We estimated the moving outbreak thresholds as the 75th percentile of the distribution of HFMD cases per month over the 120-month period, excluding 1 month at a time. The moving thresholds were applied to produce predicted probabilities of exceeding the outbreak threshold for all 10 years of 2009–2018.

Results

Temporal, spatial and population distribution of HFMD

Throughout 2009–2018, a total of 20048244 HFMD cases were reported from 31 provinces in mainland China (Table A in S1 File and Fig 1A). The seasonal index of 12 months ranged from 0.16 to 2.09 and the annual epidemic peak was observed from April to July (Fig 1B). The time-series seasonal decomposition analysis showed a significant seasonal periodicity (Fig A-A in S1 File). The wavelet spectrum power analysis at the national level revealed a pronounced band at 1 y and an inconsistent band at 0.5 y, suggesting that the long-term temporal variation of HFMD is dominated by the annual pattern (Fig A-B in S1 File). When the data are stratified by geographic region, the Western, Central and Eastern regions represented low, moderate, and high level of transmission intensity (Fig 1C). The three-dimensional trend of the average annual incidence showed an ascending arc trend from north to south (Fig 1D). The Lorenz curves showed that Central region had a higher degree of spatial heterogeneity larger than Eastern and Western regions (Fig 1E). Children under 5 years of age had the highest average incidence of HFMD (22.3 cases per 1000 population), while the boys had 1.35 times higher incidence than the girls (25.4 vs 18.8 cases per 1000 population, P < 0.001) (Fig 1F). The case-fatality rate of HFMD (overall 0.019%) was highest among children aged < 12 months (0.034%) (Fig 1G).

thumbnail
Fig 1. Temporal, spatial and population distribution of HFMD in China during 2009–2018.

(A) Time series of monthly HFMD incidence. (B) The seasonality index is calculated as the ratio of the average number of cases in a given month to the average number of cases per month over 10×12 months, and its 95% CI is calculated using SD. A seasonality index greater than 1.0 indicates a seasonal trend. (C) The geographic distribution of HFMD incidence. The Western region includes 11 provinces located in the interior of China: Shaanxi, Gansu, Qinghai, Ningxia, Xinjiang, Sichuan, Chongqing, Yunnan, Guizhou, Guangxi and Tibet. The Central region includes 9 provinces: Shanxi, Jilin, Heilongjiang, Inner Mongolia, Anhui, Jiangxi, Henai, Hubei, and Hunan. The Eastern region includes 11 coastal provinces: Liaoning, Hebei, Beijing, Tianjin, Shangdong, Jiangsu, Shanghai, Zhejiang, Fujian, Guangdong, and Hainan. The image is generated in ArcGIS Pro version 2.5 (ESRI Inc., Redlands, CA, USA), using a freely downloaded shapefile from the National Geomatics Center of China (https://www.ngcc.cn/ngcc/html/1/). (D) Three-dimensional trend of the average annual HFMD incidence. (E) Lorenz curves of the distribution HFMD cases as a function of population size at the province level. The black line (first diagonal) represents a constant distribution of HFMD cases (no heterogeneity). Age-gender distribution of reported HFMD cases (F) and case fatality rate (G).

https://doi.org/10.1371/journal.pntd.0011286.g001

Regional HFMD seasonal patterns

Over the 10-year period, HFMD incidence has remained relatively steady. HFMD seasonality varies across the country, with the peak transmission season occurring earlier in the year in the south (Fig 2). We modeled the seasonal models based on negative binomial (NB) generalized linear models (GLMs) with harmonic terms [34]. The median model fitness (R2) for all provinces reached 0.63 (range 0.28–0.81) (Fig B in S1 File). The estimating seasonality displayed significant latitudinal gradients in peak time, annual amplitude, and contribution of semi-annual periodicity (measured by the ratio of the semi-annual amplitude to the sum of the annual and semi-annual amplitudes) of HFMD in China. The gradients in Central region were stronger than those in the other two regions, suggesting higher geographic dependence (P < 0.05, (Figs 3 and C in S1 File).

thumbnail
Fig 2. Spatial and temporal variation in HFMD incidence in China during 2009–2018.

Monthly HFMD incidence rates (per 100000 people) between January, 2009, and December, 2018 were aggregated at the province level. Provinces were ordered by their geographical location.

https://doi.org/10.1371/journal.pntd.0011286.g002

thumbnail
Fig 3. Latitudinal gradients in seasonality of HFMD in China.

(A) Peak time. (B) Annual amplitude. (C) Semiannual amplitude. (D) Ratio semiannual component (higher ratio suggests a stronger semiannual periodicity). Symbol size is proportional to the number of cases in each province. Black solid lines represent linear regression fit. Colors represent different geographic regions (red = Eastern region, green = Central region, blue = Western region).

https://doi.org/10.1371/journal.pntd.0011286.g003

Spatial and spatiotemporal clustering of HFMD

We identified significant positive spatial autocorrelation for HFMD incidence at the provincial level during 2009–2018. Global Moran’s I index ranged from 0.416 to 0.608 (Table B in S1 File), indicating that HFMD incidence was geographically clustered. Based on local indicators of spatial autocorrelation analysis, the hot spots (high-high cluster) were observed in the north in 2009, covering seven provinces. Changes were observed since 2010 when the high-risk cluster moved toward south, and since 2014 when the cold spots (low-low cluster) were moving eastward (Fig D in S1 File).

Based on Kulldorff’s space-time scan analysis [35], one most likely high-risk spatiotemporal clustering area and 1–2 secondary-risk clustering areas were detected annually in China during 2009–2018. The most likely cluster areas for HFMD were located in south China (Guangdong, Guangxi, Hunan, Hainan) from April 2013 to October 2017. The secondary cluster areas were located in some areas of central and east China (Beijing, Tianjin, Hebei, Shandong, Jiangsu, Shanghai, Zhejiang) from April to July 2014 (Fig 4 and Table C in S1 File).

thumbnail
Fig 4. Spatiotemporal clustering of HFMD incidents in China from 2009 to 2018.

Kulldorff’s space-time scan statistic was used to identify the geographic areas and time periods of potential clusters with significantly higher incidence of HFMD than those of neighboring areas. The window with the maximum log-likelihood ratio (LLR) was defined as the most likely cluster area, and other windows with significant LLR as the secondary potential clusters. A statistically significant assessment of 9999 replicates was performed using Monte Carlo simulation with a significance level of 0.05. For other parameters, the maximum radius of the base of the circle is set to 50% of the total population at risk and the maximum height of the cylinder is set to 50% of the total study period. More detailed information is provided in Table C in S1 File. The image is generated in Arcgis Pro version 2.5 (ESRI Inc., Redlands, CA, USA), using a freely downloaded shapefile from the National Geomatics Center of China (https://www.ngcc.cn/ngcc/html/1/).

https://doi.org/10.1371/journal.pntd.0011286.g004

Model selection and evaluation

Using variance inflation factor (VIF) as a screening tool, four covariates with VIF > 10 were considered to be collinear and were removed (Fig E and Table D in S1 File). Further, the Bayesian model averaging (BMA) approach was applied to the regression models. As shown in Table E in S1 File, compared with other selected models, model 1 (considering nine explanatory variables) had the highest posterior model probability (PMP) (50.0% of the total posterior probability).

Moving from traditional linear and nonlinear model (MLR and DLNM), we noted progressive improvements in model fitting and predictive performance for hierarchical Bayesian spatiotemporal model and its incorporation form with DLNM (Table 1). To select the best fitting models, we included province-level monthly random effects, to account for varying seasonality between different areas (Fig F in S1 File), and year-specific spatial random effects to account for unexplained interannual spatial heterogeneity (Fig G in S1 File).

thumbnail
Table 1. Bayesian modeling evaluations of the alternative regressions for China’s HFMD case account for model fitness complexity and predictive power.

https://doi.org/10.1371/journal.pntd.0011286.t001

Detection of outbreaks

The models explicitly incorporated data on main environmental and socioeconomic covariates, and showed the out-of-sample predicted vs. observed HFMD incidence rates (per 100000 population) at the provincial level for China. We found that our models could correctly capture the general timing and intensity of HFMD outbreaks, except in two provinces of Sichuan and Hainan (median goodness-of-fit R2 = 0.78, range 0.28–0.89, Fig 5). Fig 6 shows the probability of exceeding the moving outbreak threshold in two representative provinces (i.e., Shandong and Guangdong). It well distinguished between outbreak and non-outbreak months throughout the ten years.

thumbnail
Fig 5. Predicted versus observed HFMD incidence rates per province.

Mean observed HFMD incidence rate per 100000 people (black curve), and corresponding mean out-of-sample predicted incidence rate (red curve) and 95% credible interval (shaded area) simulated from the HBSTM+DLNM model (refitted 12 × 10 months from Jan 2009 to Dec 2018, excluding the month for which the prediction was valid). Moving outbreak thresholds (blue curve) were calculated as the 75th percentile of the distribution of HFMD cases per month over the 120-month period, leaving out one month at a time. Provinces are ordered by their geographical location.

https://doi.org/10.1371/journal.pntd.0011286.g005

thumbnail
Fig 6. Predicted probability of exceeding moving outbreak thresholds.

Predicted probability of observing a HFMD outbreak in Shandong(A) and Guangdong (B) provinces in China. The moving outbreak thresholds were calculated as the 75th percentile of observed HFMD cases per month from Jan 2009 to Dec 2018, excluding the month for which the prediction is valid. The out-of-sample predictive distributions were simulated from the HBSTM+DLNM model. The graduated color bars represent the predicted probabilities ranging from 0 (pale colors) to 1 (deep colors). The months in which the moving outbreak threshold was exceeded are marked with a pentagram.

https://doi.org/10.1371/journal.pntd.0011286.g006

Combined impacts of environmental and socioeconomic covariates

Two kinds of overall impacts of environmental and socioeconomic covariates on HFMD dynamic were estimated. One was the linear effects based on HBSTM and the other was nonlinear and delayed effects obtained from HBSTM+DLNM. We summarized the critical parameters of the linear HBSTM and HBSTM+DLNM in Table 2, including the overall coefficients and their relative risks. This HBSTM explained 85.6% of the variance in HFMD transmission during 2009–2018 while HBSTM+DLNM accounted for 86.9%. In terms of the environmental covariates, temperature, relative humidity and normalized difference vegetation index (NDVI) displayed as significant positive stimulants for HFMD transmission. Besides, among the socioeconomic covariates, the population density, birth rate, real GDP per capita demonstrated positive association with HFMD incidence while school vacation served as a restraining force (nearly half reduction).

thumbnail
Table 2. Linear impacts of environmental and socioeconomic covariates on HFMD risk in China during 2009–2018.

https://doi.org/10.1371/journal.pntd.0011286.t002

Furthermore, we found significant nonlinear associations of environmental covariates and HFMD transmission, allowing interaction with underlying socioeconomic conditions (Fig 7). The results of the linear model (Table 2) seemed to obscure the nonlinear trends of such associations. The nonlinear trends varied by covariates. Fig 7A and 7B captured the thresholds of the climate covariates (monthly average temperature of 2°C, monthly average relative humidity of 60%) when their effects started to change markedly. In the linear regression model, higher NDVI was associated with higher risk of HFMD; however the direction of this association reversed in the nonlinear model. Fig 7D shows a negative association between NDVI and HFMD risk and the curve became relatively flat for the value range beyond 0.6. Fig 7E–7H shows the lag-response association for the extreme values (5% and 95% percentiles) of the environmental covariates. Fig 7I–7L were contour plots of the exposure-lag-response association between the environmental covariates and HFMD risk. In this context, population density (RR = 1.261; 95%CI, 1.169–1.353), birth rate (RR = 1.058; 95%CI, 1.025–1.090), real GDP per capita (RR = 1.163; 95%CI, 1.033–1.310) and school vacation (RR = 0.507; 95%CI, 0.459–0.559) were identified to have positive or negative effects on HFMD respectively (Table 2). Finally, we obtained a provincial estimate of the incidence of Chinese mainland HFMD in the 10-year period 2009–2018 through this model. The darker the province, the higher the incidence. The incidence predicted by the model was close to the actual incidence result, and the incidence center gradually shifted from the north to the south (Fig H in S1 File).

thumbnail
Fig 7. Nonlinear and delayed impacts of environmental covariates on HFMD risk.

(A-D) Cumulative risks over the 3-month period. Lines and shaded areas showed the relative risk (RR) and their 95%CI for the association between each environmental covariate and HFMD incidence, relative to that covariate’s median. (E-H) Lag-response association for the extreme values (5% and 95% percentiles) of the environmental covariates. (I-L) Contour plots of the exposure-lag-response association between the environmental covariates and HFMD incidence.

https://doi.org/10.1371/journal.pntd.0011286.g007

In summary, our study expanded the limited knowledge of the associations between various covariates and HFMD in China from a spatiotemporal heterogeneous perspective. It may prove valuable for advancing our understanding of complex HFMD transmission system, improving decision-making and addressing global challenges. The framework could provide implications for research on other disease-environment systems.

Discussion

There are rare studies on the spatiotemporal patterns of HFMD transmission at a national scale, but that could prove valuable in the design of geographically-tailored surveillance and intervention strategies such as vaccination. We used HFMD incidence data of China, a geographically diverse country, to quantify spatiotemporal patterns of incidence and clustering for HFMD during the period 2009–2018. It is not surprising that our results support a diversity of HFMD epidemic patterns which are driven by both environmental and socioeconomic covariates. In this study, we developed a statistical model which could successfully predict months with HFMD outbreaks versus non-outbreaks in provinces of China from Jan 2009 to Dec 2018. In practice, the CCDC may use this model to generate disease forecasts of HFMD using the weather forecast products and socioeconomic parameters. Public health decision-makers could use such forecasts as an early warning tool to plan interventions to reduce risk of HFMD and other infectious diseases.

We quantified spatial heterogeneity across geographic regions using Lorenz curves and their corresponding Gini indices. The underlying distributions of HFMD seasonal peaks and epidemic amplitudes across China are also characterized. The spatial heterogeneity was highest in the Central regions. Meanwhile, HFMD had the highest latitude-dependency of the seasonal trends and characteristics within the same region. Such results highlighted the importance of incorporating geographic dimension in appropriately disentangling and interpreting the relevant seasonal patterns and distribution of infectious diseases. That is, data aggregation may complicate our understanding of the HFMD dynamics in China. The complication may be resulted from epidemiological patterns which differ across spatially heterogeneous areas.

Earlier studies suggested simple linear relationships between temperature, relative humidity and HFMD incidence in Japan [36], Hong Kong [37], Shandong province [38] and Guangzhou city [39], China. Other studies showed that HFMD was associated with lagged climate covariates in Beijing and Changsha of China [40,41]. Consistent with these studies, we found that hot and humid climate may be a preferred environment for HFMD transmission. Possible explanation relies on the virological evidence that both temperature and humidity are crucial for the survival of EV [42]. Low temperature (< 2°C) may not favor the survival of virus and personal close contacts, which may hinder the transmission of HFMD. Higher humidity allows the virus to persist longer on inanimate surfaces [43]. Furthermore, in the temporal dimension, the disease-climate associations between HFMD and temperature, relative humidity were very similar to the temporal trend of HFMD itself, indicating that the two climate covariates be better explanatory factors for HFMD incidence at a localized temporal scale.

Studies on association between the NDVI and HFMD risk are scarce. Two studies reported a negative association, with each 1σ increase in NDVI reducing the incidence of HFMD by 20% and 37% [44,45], while another one suggested a positive association [46]. We found a negative nonlinear relationship between the two. Besides, it is worth noting that the NDVI increased in the spring and summer, coincided with the HFMD season, thus reflecting a positive association (in the linear model). Our study advanced the previous studies by using hierarchical Bayesian models to assess the delayed and nonlinear effects of NDVI. A study found that high NDVI may be associated with enhanced physical, psychological, and social well-being that can improve immunity, such as by providing a visually complex environment that reduces stress and mental fatigue, or by adding the look and feel of a place to provide a pleasant location for social or physical activity [47].

In regards to socioeconomic conditions, most covariates we tested were positively associated with the HFMD incidence in China, i.e., the population density, birth rate and per capita GDP. In developed areas, the higher population density leads to an easy spread of the virus. The school vacation was the only socioeconomic covariate that was markedly negatively correlated with HFMD incidence, possibly because the reduced the close contact of children [48]. Another study found that school holidays did not substantially reduced HFMD transmission in Hong Kong [49]. This might be explained by the fact that Hong Kong has a number of breaks and school holidays which vary every year, different from those in mainland China. The population-level findings in this study suggest that some measures such as the annual school vacation may be critical to reducing the incidence of HFMD.

For effective disease forecasting models, it is vital to identify the key drivers, lag periods, and appropriate model formulation that reflects the local disease transmission ecology. In addition, it takes time for anomalies in the covariates to manifest and contribute to disease risk. Our study selected a model to infer combined impacts of environmental and socioeconomic covariates on HFMD risk and predict the probability of exceeding an outbreak threshold. The translation of probabilities into discrete warnings needs to be carried out with caution and reflecting the capacity of the public health system to respond to imminent outbreaks.

Admittedly, there are several limitations which should be better addressed in future research. First, due to a lack of available data, there might be some important effect modifiers that we failed to take into consideration, such as viral type, reproduction number and public health interventions. For example, the inactivated EV-A71 vaccine has shown certain efficacy in reducing the incidence of EV71-related HFMD among the vaccination population [50]. Due to changes in dominant disease-causing virus types [5] as well as the short application (since 2016) and low vaccination rate [8], the vaccination program seemed not to cause substantial changes in HFMD epidemics and thus were not considered in the current analysis. Second, our analysis quantified the associations between the influential covariates and HFMD transmission, but such associations were not necessarily causal. For example, the GDP per capita may just represent economic activities, which were irrelevant to transmission. Third, HFMD cases and covariates may vary across different cities within the province. City-level comparative studies are warranted to address this issue which is particularly important for a vast province with unbalanced development. Last but not least, the relative risk map for the future trends of HFMD has not been projected. Actually, in the context of ongoing COVID-19 pandemic, the behaviors of many infectious diseases have been changed [51,52].

In summary, our study expanded the limited knowledge of the associations between various covariates and HFMD in China from a spatiotemporal heterogeneous perspective. It may prove valuable for advancing our understanding of complex HFMD transmission system, improving decision-making and addressing global challenges. The framework could provide implications for research on other disease-environment systems.

Supporting information

S1 File.

Table A in S1 File: Characteristics of all HFMD cases in China from 2009 to 2018. Fig A in S1 File: Seasonality and periodicity of HFMD cases in China during 2009–2018. Fig B in S1 File: Fit of seasonal models in 31 provinces in China. Fig C in S1 File: Seasonal estimates of HFMD in China during 2009–2018. Table B in S1 File: Global spatial autocorrelation analysis. Fig D in S1 File: Local hotspot clusters of the annual incidence of HFMD in China during 2009–2018. Table C in S1 File: Spatiotemporal clustering of HFMD incidents in China during 2009–2018. Table D in S1 File: Province-level potential explanatory variables of regional HFMD in China: EC1-EC8 denote environmental factors and SC1-SC10 denote socioeconomic factors. Fig E in S1 File: Variables screening procedure: remove variables with higher multicollinearity (VIF > 10). Table E in S1 File: BMA model selection. Fig F in S1 File: Contribution of the province-specific monthly random effects to HFMD incidence rate (IR) estimates. Fig G in S1 File: Contribution of year-specific spatial random effects to HFMD incidence rate (IR) estimates. Fig H in S1 File: Posterior predictive mean HFMD incidence rate 2009–2018.

https://doi.org/10.1371/journal.pntd.0011286.s001

(DOCX)

References

  1. 1. Koh WM, Bogich T, Siegel K, Jin J, Chong EY, Tan CY, et al. The Epidemiology of Hand, Foot and Mouth Disease in Asia: A Systematic Review and Analysis. Pediatr Infect Dis J. 2016;35(10):e285–300. Epub 2016/06/09. pmid:27273688; PubMed Central PMCID: PMC5130063 Communicable Disease Public Health Research (CDPHRG12NOV021), the Centre for Infectious Disease Epidemiology and Research, the Ministry of Education Tier 1 grant and the President’s Graduate Fellowship to W.M.K. The funders had no role in the decision to publish. T.B. is employed by commercial company, Standard Analytics. The remaining authors have no financial relationships relevant to this article to disclose. The authors have no conflicts of interest to disclose.
  2. 2. Zhang Y, Zhu Z, Yang W, Ren J, Tan X, Wang Y, et al. An emerging recombinant human enterovirus 71 responsible for the 2008 outbreak of hand foot and mouth disease in Fuyang city of China. Virol J. 2010;7:94. Epub 2010/05/13. pmid:20459851; PubMed Central PMCID: PMC2885340.
  3. 3. Xing W, Liao Q, Viboud C, Zhang J, Sun J, Wu JT, et al. Hand, foot, and mouth disease in China, 2008–12: an epidemiological study. Lancet Infect Dis. 2014;14(4):308–18. Epub 2014/02/04. pmid:24485991; PubMed Central PMCID: PMC4035015.
  4. 4. Hong J, Liu F, Qi H, Tu W, Ward MP, Ren M, et al. Changing epidemiology of hand, foot, and mouth disease in China, 2013–2019: a population-based study. Lancet Reg Health West Pac. 2022;20:100370. Epub 2022/01/18. pmid:35036978; PubMed Central PMCID: PMC8743221.
  5. 5. Weng Y, Chen W, He W, Huang M, Zhu Y, Yan Y. Serotyping and Genetic Characterization of Hand, Foot, and Mouth Disease (HFMD)-Associated Enteroviruses of No-EV71 and Non-CVA16 Circulating in Fujian, China, 2011–2015. Med Sci Monit. 2017;23:2508–18. Epub 2017/05/26. pmid:28539579; PubMed Central PMCID: PMC5452872.
  6. 6. Yang Q, Ding J, Cao J, Huang Q, Hong C, Yang B. Epidemiological and etiological characteristics of hand, foot, and mouth disease in Wuhan, China from 2012 to 2013: outbreaks of coxsackieviruses A10. J Med Virol. 2015;87(6):954–60. Epub 2015/03/11. pmid:25754274.
  7. 7. Bian L, Wang Y, Yao X, Mao Q, Xu M, Liang Z. Coxsackievirus A6: a new emerging pathogen causing hand, foot and mouth disease outbreaks worldwide. Expert Rev Anti Infect Ther. 2015;13(9):1061–71. Epub 2015/06/27. pmid:26112307.
  8. 8. Qi L, Su K, Xia Y, Tang W, Shen T, Li Q. Enterovirus 71 vaccine acceptance among parents of children < 5 years old and their knowledge of hand, foot and mouth disease, Chongqing, China, 2017. PLoS One. 2019;14(11):e0225569. Epub 2019/11/28. pmid:31774839; PubMed Central PMCID: PMC6881008.
  9. 9. He X, Zhang M, Zhao C, Zheng P, Zhang X, Xu J. From Monovalent to Multivalent Vaccines, the Exploration for Potential Preventive Strategies Against Hand, Foot, and Mouth Disease (HFMD). Virol Sin. 2021;36(2):167–75. Epub 2020/10/01. pmid:32997323; PubMed Central PMCID: PMC7525078.
  10. 10. Kang L, He C, Miao L, Liang J, Zhu J, Li X, et al. Geographic disparities in pneumonia-specific under-five mortality rates in Mainland China from 1996 to 2015: a population-based study. Int J Infect Dis. 2017;59:7–13. Epub 2017/03/28. pmid:28342801.
  11. 11. Bo Z, Ma Y, Chang Z, Zhang T, Liu F, Zhao X, et al. The spatial heterogeneity of the associations between relative humidity and pediatric hand, foot and mouth disease: Evidence from a nation-wide multicity study from mainland China. Sci Total Environ. 2020;707:136103. Epub 2019/12/25. pmid:31874401.
  12. 12. Xin X, Hu X, Zhai L, Jia J, Pan B, Han Y, et al. The effect of ambient temperature on hand, foot and mouth disease in Qingdao, China, 2014–2018. Int J Environ Health Res. 2022:1–10. Epub 2022/05/06. pmid:35510292.
  13. 13. Luo C, Ma Y, Liu Y, Lv Q, Yin F. The burden of childhood hand-foot-mouth disease morbidity attributable to relative humidity: a multicity study in the Sichuan Basin, China. Sci Rep. 2020;10(1):19394. Epub 2020/11/12. pmid:33173087; PubMed Central PMCID: PMC7656260.
  14. 14. Hao J, Yang Z, Yang W, Huang S, Tian L, Zhu Z, et al. Impact of Ambient Temperature and Relative Humidity on the Incidence of Hand-Foot-Mouth Disease in Wuhan, China. Int J Environ Res Public Health. 2020;17(2). Epub 2020/01/16. pmid:31936369; PubMed Central PMCID: PMC7013846.
  15. 15. Zhang X, Xu C, Xiao G. Space-time heterogeneity of hand, foot and mouth disease in children and its potential driving factors in Henan, China. BMC Infect Dis. 2018;18(1):638. Epub 2018/12/12. pmid:30526525; PubMed Central PMCID: PMC6286567.
  16. 16. Hu Y, Xu L, Pan H, Shi X, Chen Y, Lynn H, et al. Transmission center and driving factors of hand, foot, and mouth disease in China: A combined analysis. PLoS neglected tropical diseases. 2020;14(3):e0008070. Epub 2020/03/10. pmid:32150558; PubMed Central PMCID: PMC7062235.
  17. 17. Sun J, Wu S, Yan Z, Li Y, Yan C, Zhang F, et al. Using Geographically Weighted Regression to Study the Seasonal Influence of Potential Risk Factors on the Incidence of HFMD on the Chinese Mainland. 2021;10(7):448.
  18. 18. Tang W, Liu S, Kang P, Peng X, Li Y, Guo R, et al. Quantifying the lagged effects of climate factors on vegetation growth in 32 major cities of China. Ecological Indicators. 2021;132:108290.
  19. 19. Li L, Liu Y, Wu P, Peng Z, Wang X, Chen T, et al. Influenza-associated excess respiratory mortality in China, 2010–15: a population-based study. Lancet Public Health. 2019;4(9):e473–e81. Epub 2019/09/09. pmid:31493844; PubMed Central PMCID: PMC8736690.
  20. 20. Mao Y, Zhang N, Zhu B, Liu J, He R. A descriptive analysis of the spatio-temporal distribution of intestinal infectious diseases in China. BMC Infect Dis. 2019;19(1):766. Epub 2019/09/04. pmid:31477044; PubMed Central PMCID: PMC6721277.
  21. 21. Pons-Salort M, Oberste MS, Pallansch MA, Abedi GR, Takahashi S, Grenfell BT, et al. The seasonality of nonpolio enteroviruses in the United States: Patterns and drivers. Proc Natl Acad Sci U S A. 2018;115(12):3078–83. Epub 2018/03/07. pmid:29507246; PubMed Central PMCID: PMC5866597.
  22. 22. Chowell G, Munayco CV, Escalante AA, McKenzie FE. The spatial and temporal patterns of falciparum and vivax malaria in Peru: 1994–2006. Malar J. 2009;8:142. Epub 2009/06/30. pmid:19558695; PubMed Central PMCID: PMC2714521.
  23. 23. Al-Ahmadi K, Al-Zahrani A. Spatial autocorrelation of cancer incidence in Saudi Arabia. Int J Environ Res Public Health. 2013;10(12):7207–28. Epub 2013/12/20. pmid:24351742; PubMed Central PMCID: PMC3881162.
  24. 24. Song C, He Y, Bo Y, Wang J, Ren Z, Guo J, et al. Disease relative risk downscaling model to localize spatial epidemiologic indicators for mapping hand, foot, and mouth disease over China. Stochastic Environmental Research and Risk Assessment. 2019;33(10):1815–33.
  25. 25. Huang JX, Wang JF, Li ZJ, Wang Y, Lai SJ, Yang WZ. Visualized Exploratory Spatiotemporal Analysis of Hand-Foot-Mouth Disease in Southern China. PLoS One. 2015;10(11):e0143411. Epub 2015/11/26. pmid:26605919; PubMed Central PMCID: PMC4659604.
  26. 26. Huang R, Wei J, Li Z, Gao Z, Mahe M, Cao W. Spatial-temporal mapping and risk factors for hand foot and mouth disease in northwestern inland China. PLoS Negl Trop Dis. 2021;15(3):e0009210. Epub 2021/03/25. pmid:33760827; PubMed Central PMCID: PMC8021183.
  27. 27. Zou Y, Lin B, Yang X, Wu L, Muneeb Abid M, Tang J. Application of the Bayesian Model Averaging in Analyzing Freeway Traffic Incident Clearance Time for Emergency Management. Journal of Advanced Transportation. 2021;2021:6671983.
  28. 28. Niu Q, Liu J, Zhao Z, Onishi M, Kawaguchi A, Bandara A, et al. Explanation of hand, foot, and mouth disease cases in Japan using Google Trends before and during the COVID-19: infodemiology study. BMC Infect Dis. 2022;22(1):806. Epub 2022/10/31. pmid:36309663; PubMed Central PMCID: PMC9617033.
  29. 29. Du Z, Lawrence WR, Zhang W, Zhang D, Yu S, Hao Y. Bayesian spatiotemporal analysis for association of environmental factors with hand, foot, and mouth disease in Guangdong, China. Sci Rep. 2018;8(1):15147. Epub 2018/10/13. pmid:30310172; PubMed Central PMCID: PMC6181968.
  30. 30. Simpson D, Rue H, Riebler A, Martins TG, Sørbye SH. Penalising Model Component Complexity: A Principled, Practical Approach to Constructing Priors. 2017;32%J Statistical Science(1):1–28,.
  31. 31. Chen Y, Li N, Lourenco J, Wang L, Cazelles B, Dong L, et al. Measuring the effects of COVID-19-related disruption on dengue transmission in southeast Asia and Latin America: a statistical modelling study. Lancet Infect Dis. 2022;22(5):657–67. Epub 2022/03/06. pmid:35247320; PubMed Central PMCID: PMC8890758 declare no competing interests.
  32. 32. Lowe R, Lee SA, O’Reilly KM, Brady OJ, Bastos L, Carrasco-Escobar G, et al. Combined effects of hydrometeorological hazards and urbanisation on dengue risk in Brazil: a spatiotemporal modelling study. Lancet Planet Health. 2021;5(4):e209–e19. Epub 2021/04/12. pmid:33838736.
  33. 33. Lowe R, Gasparrini A, Van Meerbeeck CJ, Lippi CA, Mahon R, Trotman AR, et al. Nonlinear and delayed impacts of climate on dengue risk in Barbados: A modelling study. PLoS Med. 2018;15(7):e1002613. Epub 2018/07/18. pmid:30016319; PubMed Central PMCID: PMC6049902.
  34. 34. Head JR, Collender PA, Lewnard JA, Skaff NK, Li L, Cheng Q, et al. Early Evidence of Inactivated Enterovirus 71 Vaccine Impact Against Hand, Foot, and Mouth Disease in a Major Center of Ongoing Transmission in China, 2011–2018: A Longitudinal Surveillance Study. Clin Infect Dis. 2020;71(12):3088–95. Epub 2019/12/28. pmid:31879754; PubMed Central PMCID: PMC7819528.
  35. 35. Rao H, Shi X, Zhang X. Using the Kulldorff’s scan statistical analysis to detect spatio-temporal clusters of tuberculosis in Qinghai Province, China, 2009–2016. BMC Infect Dis. 2017;17(1):578. Epub 2017/08/23. pmid:28826399; PubMed Central PMCID: PMC5563899.
  36. 36. Onozuka D, Hashizume M. The influence of temperature and humidity on the incidence of hand, foot, and mouth disease in Japan. Sci Total Environ. 2011;410–411:119–25. Epub 2011/10/22. pmid:22014509.
  37. 37. Ma E, Lam T, Wong C, Chuang SK. Is hand, foot and mouth disease associated with meteorological parameters? Epidemiol Infect. 2010;138(12):1779–88. Epub 2010/09/30. pmid:20875200.
  38. 38. Liu Y, Wang X, Pang C, Yuan Z, Li H, Xue F. Spatio-temporal analysis of the relationship between climate and hand, foot, and mouth disease in Shandong province, China, 2008–2012. BMC Infect Dis. 2015;15:146. Epub 2015/04/19. pmid:25887074; PubMed Central PMCID: PMC4374415.
  39. 39. Chen C, Lin H, Li X, Lang L, Xiao X, Ding P, et al. Short-term effects of meteorological factors on children hand, foot and mouth disease in Guangzhou, China. Int J Biometeorol. 2014;58(7):1605–14. Epub 2013/11/22. pmid:24258319.
  40. 40. Xu M, Yu W, Tong S, Jia L, Liang F, Pan X. Non-Linear Association between Exposure to Ambient Temperature and Children’s Hand-Foot-and-Mouth Disease in Beijing, China. PLoS One. 2015;10(5):e0126171. Epub 2015/05/27. pmid:26010147; PubMed Central PMCID: PMC4444089.
  41. 41. Meng L, Zhou C, Xu Y, Liu F, Zhou C, Yao M, et al. The lagged effect and attributable risk of apparent temperature on hand, foot, and mouth disease in Changsha, China: a distributed lag non-linear model. Environ Sci Pollut Res Int. 2022. Epub 2022/09/13. pmid:36094702.
  42. 42. Arita M, Shimizu H, Nagata N, Ami Y, Suzaki Y, Sata T, et al. Temperature-sensitive mutants of enterovirus 71 show attenuation in cynomolgus monkeys. J Gen Virol. 2005;86(Pt 5):1391–401. Epub 2005/04/16. pmid:15831951.
  43. 43. Kramer A, Schwebke I, Kampf G. How long do nosocomial pathogens persist on inanimate surfaces? A systematic review. BMC Infect Dis. 2006;6:130. Epub 2006/08/18. pmid:16914034; PubMed Central PMCID: PMC1564025.
  44. 44. Stanaway JD. Insights from disease ecology: focus on hand, foot and mouth disease in China: University of Washington; 2013.
  45. 45. Cao C, Li G, Zheng S, Cheng J, Lei G, Tian K, et al. Research on the environmental impact factors of hand-foot-mouth disease in Shenzhen, China using RS and GIS technologies. IEEE International Geoscience and Remote Sensing Symposium; Munich 2012. p. 7240–3.
  46. 46. Li L, Qiu W, Xu C, Wang J. A spatiotemporal mixed model to assess the influence of environmental and socioeconomic factors on the incidence of hand, foot and mouth disease. BMC Public Health. 2018;18(1):274. Epub 2018/02/22. pmid:29463224; PubMed Central PMCID: PMC5819665.
  47. 47. Shanahan DF, Bush R, Gaston KJ, Lin BB, Dean J, Barber E, et al. Health Benefits from Nature Experiences Depend on Dose. Sci Rep. 2016;6:28551. Epub 2016/06/24. pmid:27334040; PubMed Central PMCID: PMC4917833.
  48. 48. Zhao J, Hu X. The complex transmission seasonality of hand, foot, and mouth disease and its driving factors. BMC Infect Dis. 2019;19(1):521. Epub 2019/06/15. pmid:31196004; PubMed Central PMCID: PMC6567494.
  49. 49. Yang B, Lau EH, Wu P, Cowling BJ. Transmission of Hand, Foot and Mouth Disease and Its Potential Driving Factors in Hong Kong. Sci Rep. 2016;6:27500. Epub 2016/06/09. pmid:27271966; PubMed Central PMCID: PMC4895171.
  50. 50. Wu H, Xue M, Wu C, Lu Q, Ding Z, Wang X, et al. Trend of hand, foot, and mouth disease from 2010 to 2021 and estimation of the reduction in enterovirus 71 infection after vaccine use in Zhejiang Province, China. PLoS One. 2022;17(9):e0274421. Epub 2022/09/21. pmid:36126038; PubMed Central PMCID: PMC9488823.
  51. 51. Song S, Wang P, Li J, Nie X, Liu L, Liu S, et al. The indirect impact of control measures in COVID-19 pandemic on the incidence of other infectious diseases in China. Public Health Pract (Oxf). 2022;4:100278. Epub 2022/06/21. pmid:35722540; PubMed Central PMCID: PMC9190177.
  52. 52. Chen YQ, Ji YF, Chen JM. Dual impacts of the COVID-19 nonpharmaceutical interventions on other infectious diseases. J Med Virol. 2022;94(10):4588–90. Epub 2022/06/09. pmid:35676236; PubMed Central PMCID: PMC9347815.