Spatiotemporal Changes and Influencing Factors of Hand, Foot, and Mouth Disease in Guangzhou, China, From 2013 to 2022: Retrospective Analysis

Introduction

Hand, foot, and mouth disease (HFMD) is an infectious disease caused by enteroviruses, which mostly occurs in children younger than 5 years of age. The main characteristics are small herpes or small ulcers on the hands, feet, mouth, and other parts. Most children recover spontaneously in about a week, but a few children can face complications, such as myocarditis, pulmonary edema, and aseptic meningoencephalitis. In some severe cases, the illness progresses quickly, resulting in death [1]. Due to the scarcity of effective therapeutic medications, symptomatic care is mainly used. HFMD was first reported in New Zealand in 1957. In the past 30 years, it has shown the characteristics of widespread epidemics and infections in the Asia-Pacific region, causing more than 2 million hospitalizations in Asia every year as well as sporadic outbreaks in Europe and the United States [2-5]. Since May 2008, the disease has been included in the management of Category C notifiable infectious diseases under China’s “Infectious Disease Prevention and Control Law.” The average number of reported cases per year exceeds 1 million, ranking it among the top 3 notifiable infectious diseases in China all year round. Dynamic sequence analysis shows that the reported incidence of HFMD is increasing in China, with the largest increase occurring in the southern region and the southeast coastal region as a hot spot and high-risk area [6,7]. Guangzhou is a regional central city in South China and China’s southern gateway to the world. It has a dense population and frequent mobility and has a typical marine monsoon climate in the south subtropics. Many factors have jointly contributed to the prevalence of HFMD in Guangzhou. The number of HFMD incidences in Guangzhou is among the highest in China all year round [8,9].

Due to differences of socioeconomic factors in different regions, the spread and distribution of the HFMD virus show significant spatial heterogeneity [10]. In this case, spatial epidemiology plays a key role as an important branch of epidemiology, using geographic information system (GIS) technology to accurately describe and analyze the spatial distribution characteristics and development patterns of HFMD in different regions [11]. Then, we can more accurately determine key areas for epidemic prevention and control, guide the deployment and adjustment of prevention and control work, effectively control the spread of the epidemic, and protect public health to the greatest extent.

So far, the spatiotemporal distribution of HFMD and its influencing factors have been extensively studied [12-14], but few articles have combined the understanding of spatiotemporal distribution with the control of potential influencing factors. To control the epidemic of HFMD more effectively, more in-depth research in this area is necessary. Conventional linear regression methods are often used to examine the influence of potential factors, but they have some limitations; they cannot capture the changing trends of data that may exist in different regions or time points, and data between adjacent locations may interact with each other. On the contrary, geographically weighted regression (GWR) is an extension of ordinary least squares regression (OLS) for spatial data analysis. It can consider the impact of spatial distribution on the model, fit a local model for each observation point to realize the spatial change of the model, and reflect the impact of geographical location more accurately, thereby improving the reliability and interpretability of analysis results [15]. Moreover, the geographically and temporally weighted regression (GTWR) model builds upon GWR by incorporating a temporal dimension. This enhancement enables the model to account for both spatial and temporal nonstationarity, providing a comprehensive understanding of the spatiotemporal dynamics of HFMD and its influencing factors [16]. Therefore, this study described the spatiotemporal distribution of HFMD in Guangzhou through spatial autocorrelation analysis and space-time scan statistics analysis. Additionally, it analyzed the impact of population, economy, education, and medical care on the spatiotemporal distribution of HFMD in Guangzhou from 2013 to 2022 through multivariate OLS, GWR, and GTWR to provide a reference for precise prevention and control of HFMD in southern China.

Methods

Study Area

Guangzhou is located in the lower reaches of the Pearl River in South China, bordering on the South China Sea (latitude 22°26′ to 23°56′N; longtitude 112°57′ to 114°3′E). Guangzhou has 11 county-level administrative districts, a total of 136 streets, and 35 towns (Figure 1). Due to its coastal geographical location as well as population diversity and mobility, Guangzhou has become a prosperous economic center, with a permanent population of more than 18.73 million in 2022. The annual per capita GDP exceeds 150,000 Yuan (USD 20,637). According to the Guangzhou Urban Master Plan (2017‐2035) [17], Guangzhou’s 11 districts are divided into 3 levels. The central urban area includes Liwan, Yuexiu, Tianhe, and Haizhu districts; the subcenter includes Baiyun, Huangpu, and Nansha districts; and the peripheral areas include Huadu, Panyu, Conghua, and Zengcheng districts.

‎

Data Sources

Since the Ministry of Health included HFMD in the management of notifiable infectious diseases in 2008, every HFMD case is required to be reported to the notifiable infectious disease information system in real time. Information on HFMD cases in Guangzhou from 2013 to 2022 was obtained through the Infectious Disease Information Management System of the Guangzhou Center for Disease Control and Prevention. Only HFMD cases whose current address was in Guangzhou were included in the study. The population, economy, education, and medical care data of each administrative districts came from the Guangzhou Municipal Bureau of Statistics [18].

Spatial Autocorrelation Analysis

Spatial autocorrelation analysis is a statistical method used to explore the distribution of data in geographical space. It mainly covers the following two aspects: (1) global spatial autocorrelation analysis evaluates the overall degree of spatial aggregation or dispersion trend of a phenomenon in the entire region, and (2) local spatial autocorrelation analysis identifies and quantifies local clustering or dispersion patterns of a phenomenon in geographic space [19]. The global Moran’s I was used to measure the strength of the overall spatial autocorrelation. The values ranged from −1 to 1, with values close to 1 indicating positive spatial autocorrelation, values close to −1 indicating negative spatial autocorrelation, and values close to 0 indicating no spatial autocorrelation. The local Moran’s I (I_i) was used to measure the correlation between a region and its neighboring regions. The values ranged from −1 to 1, with positive values indicating high-high (HH) or low-low aggregation areas and negative values indicating high-low or low-high aggregation areas. HH aggregation area was where an area had a higher incidence than the average incidence and was surrounded by adjacent areas higher than average incidence. The high-low aggregation area was clustered with an area with an incidence above the average incidence, but surrounded by adjacent areas with incidences below the average. After statistically testing the local indicator of spatial association (LISA) through the Z test, the LISA aggregation map was drawn to intuitively show the different aggregation forms in geographical space.

The global Moran’s I and I_i were calculated as follows [20]:

${\displaystyle {\displaystyle I=\frac{n}{W}\frac{\sum _{i=1}^n\sum _{j=1}^n{w}_{ij}\left(x_i-\overline{x}\right)\left(x_j-\overline{x}\right)}{\sum _{i=1}^n{\left(x_i-\overline{x}\right)}^2}\left(1\right)}}$

${\displaystyle {\displaystyle I_i=\frac{\left(x_i-\overline{x}\right)\sum _{j=1}^n{w}_{ij}\left(x_j-\overline{x}\right)}{\sum _{j=1}^n{\left(x_j-\overline{x}\right)}^2}\left(2\right)}}$

where n was the number of geographic cells; W was the weight matrix; w_ij was the spatial weight between geographic cells i and j; x_i and x_j were the values of the variables over geographic cells i and j; and x̄ was the mean value of the variables.

Space-Time Scan Statistics

Space-time scan statistics is a statistical method for detecting spatiotemporal aggregation. Its basic principle is to divide the study area into different spatiotemporal windows by establishing dynamically adjustable cylinder windows. The base of the cylinder represents the geographical space, and the height of the cylinder represents the time range [21]. Considering the epidemiological characteristics of HFMD [22], this study used the maximum radius not exceeding 30% of the total population as the spatial scale and the maximum aggregation of 3 months as the time scale. Contrast distributions were generated through stochastic simulation; the expected number of events per spatiotemporal window was calculated and compared with the number of observed events; thus, the test statistic log-likelihood ratio (LLR) was calculated. The P values were calculated by Monte Carlo simulation. When P<.05, a greater LLR value indicated a greater probability that the scanning window area was an aggregation area. The area with the largest LLR value was selected as a class I aggregation area, and other areas were selected as the class II aggregation areas.

Influencing Factors Analysis

The incidence of HFMD is related to many social factors. Based on previous studies and combined with local characteristics [23-25], four aspects including population, economy, education, and medical care were selected to analyze the impact of the spatiotemporal distribution of HFMD cases. Specific indicators are available in Multimedia Appendix 1.

GWR and GTWR are spatial statistical methods used to explore the spatial relationship between dependent variables and independent variables. Considering that factors such as population, economy, education, and medical care are unevenly distributed in space, GWR takes geographical coordinates as part of the regression model. When fitting the model, it will adapt to local changes for each geographical location and fit different parameters accordingly.

GWR model and GTWR model were calculated as follows [26,27]:

${\displaystyle {\displaystyle y_i={\beta }_0\left(u_i,v_i\right)+\sum _j^k{\beta }_j\left(u_i,v_i\right)x_{ij}+{\epsilon }_i\left(3\right)}}$

${\displaystyle {\displaystyle y_i^t={\beta }_0^t\left(u_i,v_i\right)+\sum _j^k{\beta }_j^t\left(u_i,v_i\right)x_{ij}^t+{\epsilon }_i^t\left(4\right)}}$

where y_i is the observed value of the dependent variable at the spatial location; β₀(u_i, v_i) is the intercept term; β_j(u_i, v_i) is the coefficient of the independent variable j at the spatial location i; x_ij is the observed value of the independent variable j at the spatial location i; and ε_i is the error term representing the unexplained part of the model and is the observed value of the dependent variable at time t and spatial position i; the other indicators are the same.

The corrected Akaike information criterion (AIC) and R² values were used to evaluate the OLS, GWR, and GTWR models. A lower AIC value indicates that the model has a good fit to the data while considering complexity, and the closer the R² value is to 1, the stronger the ability of the model to explain the variables.

Statistical Analysis

Multifactorial OLS and GWR analyses were performed using R (version 4.3.2; R Foundation for Statistical Computing). Spatial autocorrelation analysis was performed using ArcGIS (version 10.8; Environmental Systems Research Institute). GTWR analysis was performed using the GTWR plugin, developed by Professor Bo Huang's team. Space-time scanning analysis was performed using SaTScan 10.1 (developed by Martin Kulldorff together with Information Management Services Inc) [28].

Ethics Approval

The study was approved by the ethics review committee of Guangzhou Centers for Disease Control and Prevention (GZCDC-ECHR-2023P0072). The incidence data used in this study were derived from the reported cases of HFND and had been anonymized and deidentified during the data export process.

Results

HFMD Prevalence

From 2013 to 2022, a total of 599,353 cases of HFMD were reported in Guangzhou, Guangdong Province. The ratio of male to female incidences was 1.51:1. The median age of patients was 2.42 (IQR 1.33‐1.75) years, and 93.64% (561,218/599,353) were children aged 5 years and younger (Multimedia Appendix 2). The ratio of scattered children (aged< 3 years) to nursery children (aged 3‐6 years) was 1.87:1. Cases aged ≥10 years accounted for only 1.55% (9304/599,353) of the total number of cases. The average annual reported incidence rate was 403.51/100,000, and the average annual change percentage was −12.94%, showing an overall downward trend. HFMD cases were characterized by a bimodal distribution, with the peak period from May to July and the second peak period from August to October (Figure 2). The average annual incidence of HFMD in Huadu District was the highest (614.90/100,000), which was 2.51 times higher than that in Yuexiu District (244.81/100,000). The top 3 distrcits with cumulative cases were Baiyun District (118,005 cases), Panyu District (79,379 cases), and Zengcheng District (68,637 cases). Multimedia Appendix 3 and Figure 3 show the changes of HFMD incidence at the district and street (town) levels.

‎

‎

Spatial Autocorrelation Analysis of HFMD Cases

Table 1 shows that the global Moran’s I ranges from 0.048 to 0.29, except for the Moran’s I of 0.015 (P=.45) in 2020, and the 10 years global Moran’s I value is 0.059 (P=.02), indicating that the HFMD cases in Guangzhou are not randomly distributed and have a high degree of positive spatial autocorrelation and a spatial aggregation trend. In Figure 3, the LISA aggregation map shows that the HH aggregation area centered on Shiling Town and Xinhua Street in the western region of Guangzhou continued to exist from 2013 to 2017, but the number of streets (towns) involved decreased. The HH aggregation area was centered on Yongning Street and Licheng Street in the eastern region from 2018 to 2022, and the number of streets (towns) involved has increased. From 2013 to 2016, the LL aggregation area was mainly distributed in the southern region, and from 2017 to 2022, it was mainly distributed in the southwest region of Guangzhou (Figure 4).

‎

Space-Time Scanning Analysis of HFMD Cases

As shown in Table 2 and Figure 5, space-time scanning analysis detected 1 to 3 cluster areas every year, with a total of 23 clusters, all of which were statistically significant (P<.001). Class I aggregation area involved an average of 41 streets and towns per year.

From a temporal perspective, clusters were observed in all months except from December to March, and class I clusters mainly occurred from May to July. From a spatial perspective, the centers of class I aggregation areas from 2013 to 2016 were all located in Huadu District, including Xinya Street, Huashan Town, Xiuquan Street, and Huacheng Street, with an average radiation radius of 27.74 km, covering Huadu District, Baiyun District, and some streets (towns) in Tianhe District. From 2017 to 2020, the class I aggregation area moved and expanded to the northeast of Guangzhou, with Lutian Town in Conghua District as the center, with an average radiation radius of 64.1 km, covering Conghua District, Zengcheng District, Huadu District, Baiyun District, and some streets (towns) in Huangpu District. From 2021 to 2022, the class I aggregation area moved and decreased to the southeast of Guangzhou, with Yongning Street and Xintang Town in Zengcheng District as the center, with an average radiation radius of 28.06 km, covering Zengcheng District, Huangpu District, Panyu District, and some streets (towns) in Tianhe District.

‎

Influencing Factors of HFMD Incidence

Before modeling, correlation analysis was conducted with HFMD incidence as the dependent variable and all potential explanatory variables as independent variables. The collinearity diagnosis results indicated that variables in the same group had multicollinearity problems. The performance of the model was improved by deleting or merging variables with high collinearity, reducing the variance inflation factor of the remaining variables was reduced to less than 10.

Multivariate OLS was used to analyze the relationship between demographic, economic, educational, and medical factors and HFMD incidence (Table 3). Population density was found to be a protective factor for HFMD. In terms of economic factors, average annual income of employees was a protective factor. The impact of medical factors on HFMD incidence was not significant. In terms of educational factors, the average number of students in each primary school and kindergarten were risk factors. The fitting effects of different models are shown in Table 3. The GTWR model had the highest R² value and the lowest AIC value, indicating that the GTWR model predicted the relationship between the incidence of HFMD and the above factors better than the OLS and GWR models. In the scatter plot (Multimedia Appendix 4) with the predicted values from the GTWR model on the x-axis and the standardized residuals on the y-axis, the residuals were randomly distributed around 0, with no apparent patterns or trends.

Some results of the GTWR model are shown in Figure 6 and Figure 7. The impact of population density on the incidence of HFMD in the west and north of Guangzhou gradually increases. Additionally, the high-impact area of employee annual wages on HFMD gradually shifts from the northern to the southern regions of Guangzhou. The average number of students per primary school migrates from the west to the south, while the average number of students per kindergarten migrates from the west and south to the east.

‎

‎

Discussion

Principal Results

This study used geospatial analysis methods and GWR models to comprehensively describe the epidemic trend of HFMD in Guangzhou and explore the spatiotemporal distribution patterns and its influencing factors. From 2013 to 2022, a total of 599,353 patients with HFMD were reported in Guangzhou, and children aged 5 years and younger accounted for more than 90% of the total cases. The incidence showed a bimodal distribution, mainly concentrated from May to July, with the subpeak period from August to November, which is consistent with the existing studies in various countries [2,29,30]. There are fewer cases in those aged 10 years and older because almost all HFMD cases in Guangzhou are mild, meaning that even if adults are infected, the symptoms may not be enough to significantly affect their daily lives, resulting in lower reporting rates among adults. The incidence rate fluctuated from 558.40 cases to 403.51 cases per 100,000 people, with an average annual decrease of 12.94%, which is similar to the incidence trend of HFMD in China in recent years [6,31]. However, HFMD is still one of the most serious infectious diseases [32]. During the COVID-19 epidemic, due to lockdown measures and social distancing [33], infants and young children were exposed to fewer viruses and bacteria, resulting in their immune systems not being adequately exercised. Therefore, the prevention and control of HFMD should continue to focus on young children [34], and attention should be paid to home and living environment hygiene, and good hygiene habits should be cultivated in children from an early age.

We conducted spatial autocorrelation analysis and space-time scan statistics analysis of HFMD at street and township scales in Guangzhou, China. The results showed that high-risk areas were mainly distributed in peripheral urban areas (eg, Huadu, Panyu, Conghua, and Zengcheng). The reason may be that the peripheral urban areas are far away from the economic and cultural center of Guangzhou. With the increasing urbanization rate of the central urban areas, the young and middle-aged people in the peripheral urban areas gradually transfered to the central urban areas. The older people and the infants left behind were relatively weak in medical awareness and knowledge of disease prevention. At the same time, the living environment and hygiene habits in the outer urban areas were lower compard to those in central urban areas. In addition, some studies have shown that the change and transformation of industrial structures have an impact on the high-risk areas of HFMD, and the tertiary industry has a greater impact on HFMD than the primary industry [35-37]. Since 2017, Conghua and Zengcheng districts have gradually transformed into commercial districts, industrial parks, tourism, and other service industries [38,39]. The rapid development of the service industry, especially the increase in catering and amusement places, has increased the frequency and scope of children’s gatherings, providing more opportunities for the spread of HFMD. Therefore, the high-risk areas of HFMD were gradually shifted from Huadu District to Conghua and Zengcheng Districts.

Existing studies widely believe that socioeconomic factors are connected to the incidence of HFMD, with population density as a potential determinant in most areas of China [24,40]. A study conducted in China found that the higher population density in developed areas made it easy for the virus to spread [10]. A spatiotemporal analysis of southern China also indicated that HFMD outbreaks occurred earlier as population density gradually increased. Moreover, the areas most prone to clustering were typically the most prosperous [23]. However, this study showed the opposite. The population density in Guangzhou increased by 820 people/km² from 2013 to 2022, and the average annual wage increased from 69,692 to 152,324 Yuan (USD 9,588 to 20,956), while HFMD incidence decreased. Another example is that the population density in the central urban area of Guangzhou exceeds 10,000 people/km², of which Yuexiu District is more than 30,000 people/km², but the incidence rate in the central urban area is lower than that in outer urban areas all year round. This may be due to Guangzhou being different from most regions in China in terms of population, economy, and culture. First, as a very developed city, although the population of Guangzhou is relatively concentrated in a small geographical area, the population structure is dominated by young and middle-aged individuals; specifically, the proportion of people older than 5 years in the central urban area is more than 90%, and the proportion of infants and young children in the key population of HFMD is not high [18]. Second, a cross-sectional survey in Guangzhou [41] showed that just 3 months after the EV71 vaccine was approved, nearly half of parents were willing to pay for their children’s vaccines, which was higher than the willingness and vaccination rate in other cities during the same period [42,43]. A study of the effectiveness of the vaccine conducted by Guangzhou [44] showed that the 2-dose EV71 vaccination was an effective measures to prevent 3-year-old children from being infected with HFMD. A questionnaire survey of 592 parents or legal guardians in 6 kindergartens in Guangzhou showed that the level of hand hygiene and infectious disease awareness of parents or legal guardians was commendable. All these factors show that Guangzhou residents attach great importance to HFMD [45]. In addition, the fecal-oral transmission is one of the routes for intestinal infectious diseases, such as HFMD [2]. Guangzhou’s dietary habits predominantly favor healthy cooking techniques like steaming and stewing. Therefore, the ingredients are of high quality and hygienic standards, and heating them for a long time helps destroy bacteria and viruses in food, thereby reducing the risk of infectious diseases.

We also found that the average number of students per primary school or kindergarten increased the risk of HFMD. Many children in kindergartens congregate in relatively small spaces, and the stagnant indoor air and close contact between children help spread HFMD pathogens. At the same time, inadequate cleaning and disinfection of daily necessities, such as toys, tableware, tables, chairs, and floors, will cause the breeding of pathogens; hand hygiene compliance among children around the world is suboptimal, regardless of whether they are in relatively rich or poor countries [34], which also increases the risk of infection. For example, from 2013 to 2017, the average number of preschool education places in Huadu District continued to grow and ranked in the forefront of Guangzhou, and the number of primary school students also increased from 109,375 to 141,029. Therefore, during this period, the high-risk areas of HFMD in Guangzhou were mainly distributed in Huadu District, while the low-risk areas were distributed in central urban areas, such as Liwan District, where the average number of children per kindergarten was the lowest all year round.

In summary, targeted prevention and control measures can be implemented to further reduce the HFMD incidence. First, while steadily advancing prevention and control measures in the central areas, special attention needs to be paid to prevention and control work in peripheral urban areas. This includes strengthening the management and supervision of places where infants and young children are concentrated, such as kindergartens and schools in peripheral urban areas; building and improving the epidemic monitoring system in peripheral urban areas; and strengthening information sharing and communication among medical institutions, schools, health departments and other units. Furthermore, it is suggested that health authorities implement a free HFMD vaccination program for children younger than 5 years in peripheral urban areas to increase the vaccination rate of infants and young children. Finally, the education department should implement measures to foster a safer environment in preschool education institutions and reduce transmission risks for infectious diseases like HFMD. This includes limiting the number of students in each class to reduce the concentration of infants and young children and carrying out special training on common infectious diseases for teachers in institutions.

Limitations

This study has some limitations. First, the Guangzhou Statistical Yearbook only collects population data of district-level administrative units, while the population data of streets and towns are estimated based on the data in the 2020 Guangzhou Census Yearbook; therefore, the incidence rate of streets and townships was not accurate enough. Second, the influencing factors of HFMD are complex and diverse. In addition to socioeconomic factors, there are also meteorological and other factors [25,46]. Although there is not much difference in temperature and humidity between districts in Guangzhou, they may also have a certain impact. Third, the symptoms of most HFMD cases are mild and similar to clinical symptoms of other childhood diseases, such as chickenpox [47]. Some parents may not be aware that their infants and young children are sick or may not take them to see a doctor. These cases may affect the results of the regression analysis.

Conclusions

From 2013 to 2022, there were obvious differences in HFMD incidence among different populations, seasons, and regions in Guangzhou. Children aged 5 years and younger are the high-risk group for HFMD. The hot and humid summer and autumn seasons are the peak periods for the disease. The peripheral urban areas are the main aggregation areas of Guangzhou. Research shows that the average number of students per primary school and kindergarten is the main reason for the increase in the HFMD incidence, while population density and per capita income play a protective role. The results of this study contribute to a deeper understanding of the epidemic characteristics of HFMD in southern China, provide an important basis for formulating targeted prevention and control strategies, and can be a reference for other countries and regions facing similar challenges.