 Research article
 Open Open Peer Review
 Published:
A spatialtemporal statistical analysis of health seasonality: explaining HFMD infections within a children population along the Vietnamese south central coast
BMC Public Healthvolume 19, Article number: 937 (2019)
Abstract
Background
Various neglected tropical diseases show spatially changing seasonality at small areas. This phenomenon has received little scientific attention so far. Our study contributes to advancing the understanding of its drivers. This study focuses on the effects of the seasonality of increasing social contacts on the incidence proportions at multiple district level of the childhood handfootmouth disease in Da Nang city, Viet Nam from 2012 to 2016.
Methods
We decomposed the nonstationary time series of the incidence proportions for the nine spatialtemporal (ST) strata in the study area, where S indicates the spatial and T the temporal stratum. The longterm trends and the seasonality are presented by the Fourier series. To study the effects of the monthly average ambient temperature and the period of preschooling, we developed a spatialtemporal autoregressive model.
Results
Seasonality of childhood handfootmouth disease incidence proportions shows two peaks in all spatial strata annually: large peaks synchronously in April and small ones asynchronously during the preschooling period. The peaks of the average temperature are asynchronous with the seasonal peaks of the childhood handfootmouth disease incidence proportions in the period between January and May, with the negative values of the regression coefficients for all spatial strata, respectively: \( {\beta}_{{\mathrm{T}}_11}^{S_1}=0.18\pm 0.07;{\beta}_{{\mathrm{T}}_11}^{S_2}=0.25\pm 0.09;{\beta}_{{\mathrm{T}}_11}^{S_3}=0.14\pm 0.05 \). The increasingly cumulative preschooling period and the seasonal component of the incidence proportions are negatively correlated in the period between August and December, with the negative values of the regression coefficients for all temporal strata, respectively: \( {\beta}_{{\mathrm{T}}_32}^{S_1}=0.40\pm 0.01;{\beta}_{{\mathrm{T}}_32}^{S_2}=0.29\pm 0.00;{\beta}_{{\mathrm{T}}_32}^{S_3}=0.25\pm 0.01 \).
Conclusions
The study shows that social contact amongst children under five years of age is the important driving factor of the dynamics of the childhood handfootmouth disease outbreaks in the study area. The preschooling season when children’s contact with each other increases stimulates the geographical variation of the seasonality of childhood handfootmouth disease infections at small areas in the study area.
Background
Neglected tropical diseases (NTDs) have shown evidently seasonal patterns with various timings, amplitudes and types [1,2,3]. NTDs remain socioeconomic burdens in developing countries, especially in lowincome communities in Asia and Africa. The considerable seasonality of these NTDs is affected by seasonal changes of both environmental and social risk factors [4, 5]. Understanding the regularity or irregularity of NTDs’ outbreaks is necessary to obtain optimal control over the disease spread [6], to reduce their burdens and to reach the UN sustainable development goals for health [7].
Amongst others, the infections of Coxsackievirus A16 (CA16) and human Enterovirus 71 (EV71) that mainly cause handfootmouth disease (HFMD) in children under 5 years of age (U5s) have been recognized as an emerging public health problem in North East and South East Asian countries [8, 9]. So far, studies of the seasonal variations of childhood HFMD infections have mainly focused on climatic driving factors such as temperature, humidity and rainfall [9]. The effects of seasonally changing social contacts on the dynamics of childhood HFMD for small areas have been overlooked, even though HFMD viruses are mainly transmitted by directly physical contacts between infected and noninfected children [10]. Additional research has been done for northeastern Asian countries [11,12,13], whereas some studies were also done for southeastern Asian regions [14, 15]. Despite the increasing application of mathematical models to understand the epidemiology of infectious diseases, spatialtemporal statistical methods (STS) have not yet been well recognized and leveraged for studying seasonality of NTDs, i.e. to quantitatively explain the spatialtemporal dynamics of the disease outbreaks [16], especially for small geographical areas.
In this research, we study the seasonality of childhood HFMD in Viet Nam and its geographical variations at small areas as affected by seasonally changing weather and the social contacts. The study of Horby et al. in Viet Nam in 2011 [17] shows that the physical contact amongst U5s themselves is the most intensive, especially during the period at preschool. Hence, the social contact amongst U5s is measured by the annually cumulated period of time they spend at preschool since their first day at preschool. We distinguish between the trend and the seasonality. The trend establishes the longterm changing pattern of the mean level, whereas the seasonality represents yearly periodic variations. We propose a two stage STS analysis to extract the seasonal patterns and to identify the driving factors of the incidences of childhood HFMD infections in Da Nang city, Viet Nam as an application. In so doing, we aimed to identify and understand the effects of seasonally changing social contacts of U5s on the dynamics of HFMD outbreaks at multiple district level.
Methods
Study area and data
Da Nang city (Fig. 1), the biggest city on the south central coast of Viet Nam has an area of 1,285 km^{2} and a population of more than 1 million people in 2016 [18]. Approximately 8% of Da Nang city’s population consists of U5s [18]. According to the guideline from the Vietnamese Ministry of Health on diagnosis and treatment of HFMD, a confirmed HFMD case is defined when a child has a positive RTPCR assay for CA16 or EV71. Weekly counts of new cases, i.e. the incidences from seven mainland districts of the city from 2012 to 2016 were obtained from the published reports of Da Nang city Preventive medicine center (https://ksbtdanang.vn/). Because these weekly counts have many small number of cases (< 5 cases per week), statistics based upon them are not reliable. To avoid this small number problem, the monthly count of new cases that is the sum of the weekly counts of new cases from this month was analyzed in this study. Note that the first and last weeks of the months can include 1 to 3 days of the previous or the next months. The demographical data were obtained from the demographic statistics of Da Nang city. The monthly population of U5s was interpolated from the yearly data with the assumption of a constant monthly growth rate. The monthly meteorological data from 2012 to 2016 were obtained from the Viet Nam institute of meteorology, hydrology and climate change (http://www.imh.ac.vn).
Spatialtemporal heterogeneity analysis
To analyze the spatialtemporal heterogeneity and identify the predominant risk factors of HFMD, we compared the qstatistic of the nine spatialtemporal stratification to that of without stratification: \( q=1\frac{\sum \limits_{h=1}^H{N}_h{\sigma}_h^2}{N{\sigma}^2} \) [19]; where: N is the total number of HFMD cases in the study area and σ^{2} is the corresponding total variance, N_{h} and \( {\sigma}_h^2 \) are those of various spatialtemporal strata that are considered in the statistical test, H is the total number of these strata. The qstatistic for heterogeneity analysis regarding covariates was provided by the geodetector library of R [19, 20].
Spatialtemporal stratification
To evaluate the effects of U5s as the population at risk on the magnitude of the seasonality of the childhood HFMD incidences, we stratified the districts by the density of U5s. This classification resulted into three spatial strata (Fig. 1): S_{1} is an urban stratum with the highest density of U5s (> 500 per km^{2}), S_{3} is the rural stratum having the lowest density (< 100 per km^{2}) and S_{2} is located in between.
We transformed the monthly counts of HFMD incidences to the monthly HFMD incidence proportions in order to correct the effects of the internal heterogeneity of the population at risk within the spatial strata. Let i = 1:12 index the months, j = 1:5 the years, d = 1:7 the districts and s = 1:3 the spatial strata. \( {IR}_{ij}^s \) is used to denote the incidence proportion at spatial stratum s in month i of year j. To account for unequal numbers of days in months, the \( {IR}_{ij}^s \) were adjusted as: \( {IR}_{ij}^s=\frac{365.25\times \sum \limits_s{I}_{ij}^d}{12\times {D}_i\times \sum \limits_s{N}_{ij}^d\times \mathrm{1,000}} \) [21], where \( {I}_{ij}^d \) is the monthly incidence proportion, \( {N}_{ij}^d \) is the population of U5s, D_{i} is the number of days per month, \( \sum \limits_s{I}_{ij}^d \) is the sum of HFMD incidences of the districts and \( \sum \limits_s{N}_{ij}^d \) is the total population of U5s belonging to the stratum s. This adjustment removes the false variations of \( {IR}_{ij}^s \) attributed to the variations of the atrisk populations and the length of time unit, partly due to the use of the weekly reported data. \( {IR}_{ij}^s \) therefore equals the number of infected U5s per 1,000 U5s at risk per month, provided that removing the infected children from the population at risk only has a marginal correction effect.
Modelling seasonality
The seasonality of \( {IR}_{ij}^s \) was diagnosed by using the seasonal index (SI). The original value of \( {IR}_{ij}^s \) for each stratum s in month i in year j was represented as the percentage of the mean value of the \( {IR}_{ij}^s \) for stratum s in year j. SI of stratum s in month i equals the median of the percentages from all 5 years: \( {SI}_i^s= median\left(\frac{IR_{ij}^s}{12^{1}\sum \limits_{i=1}^{12}{IR}_{ij}^s}\times 100,j=1:5\right) \) [22]. The median allows for the inclusion of the effects of extreme values over the years. The standard error of SI was derived from the bootstrap confident interval assuming a normal distribution. A value of \( {SI}_i^s \) > 100% indicates that the incidence proportion is above the yearly average incidence proportion and vice versa.
For the analysis of the dynamics, we applied to \( {IR}_{ij}^s \) an additive conceptual model with three main components: a trend (μ), a seasonality (ξ) and random noises (ε):
We used the Fourier series to represent different oscillations of \( \mu \left({IR}_{ij}^s\right) \) and \( \xi \left({IR}_{ij}^s\right) \) in (1) by a linear combination of sinusoidal functions with different frequencies, magnitudes and phases [23]:
Here P and Q are the number of the sine and cosine functions representing the trend and seasonal components, respectively, \( {\mu}_0^s \) and \( {\xi}_0^s \) are the intercepts defining the baselines of the trend and seasonal components, \( {A}_{\mu p}^s \), \( {B}_{\mu p}^s \), \( {A}_{\xi q}^s \) and \( {B}_{\xi q}^s \) are the Fourier coefficients and t_{ij} indexes time. The amplitude of the trend equals: \( {C}_{\mu p}^s=\operatorname{sign}\left({B}_{\mu p}^s\right){\left({A_{\mu p}^s}^2+{B_{\mu p}^s}^2\right)}^{1/2} \); that of the seasonal component equals: \( {C}_{\xi q}^s=\operatorname{sign}\left({B}_{\xi p}^s\right){\left({A_{\xi q}^s}^2+{B_{\xi q}^s}^2\right)}^{1/2} \), where sign is the function extracting the sign of a real number [24].
Iterative weighted least square estimation was used to fit (2) and (3) by minimizing the generalized crossvalidation value: GCV = RSS + λ × PEN_{L}, where RSS is the weighted residual sum of squares, λ is the smoothing parameter and PEN_{L} is a penalty term as a measure of Fourier series roughness [25].
Identifying driving factors of seasonality
To contrast the effects of preschooling period as a proxy measure of the social contacts and the weather season on the seasonality of the disease, we further stratified the seasonal components \( \xi \left({IR}_{ij}^s\right) \) into three temporal strata: T_{1} is the dry and preschooling season from January to May, T_{2} is the dry and summer holiday season from June to July and T_{3} is the rainy and preschooling season from August to December. Hence, for all j and s, \( \xi \left({IR}_{ij}^s\right) \) belong to T_{1} for i = 1:5, to T_{2} for i = 6:7 and to T_{3} for i = 8:12, respectively.
To relate \( \xi \left({IR}_{ij}^s\right) \) of T_{1}, T_{2} and T_{3} to the potential risk factors, a spatialtemporal regression model was deployed:
where N is the total number of the weather covariates, \( {\beta}_{T0}^s \) and \( {\beta}_{Tn}^s \) are the regression parameters, \( {X}_{Tnij}^s \) are the weather covariates, \( {\upvarepsilon}_{Tij}^s \) are the spatialtemporal autoregressive regression residuals, P_{Ti} are the number of months at preschools, cumulated from August current year to May next year. Notice that P_{Ti} did not vary with s and j as they were obtained by: \( {P}_{Ti}=\left\{\begin{array}{c}i+5\ \mathrm{for}\ i\in {\mathrm{T}}_1\\ {}0\ \mathrm{for}\ i\in {\mathrm{T}}_2\\ {}i7\ \mathrm{for}\ i\in {\mathrm{T}}_3\ \end{array}\right. \).
To account for spatialtemporal autocorrelation remaining in the residual, we applied the SpatialTemporal AutoRegressive model (STAR) [26] for the error term \( {\upvarepsilon}_{Tij}^s \) in (4):
where M is the total number of the firstorder spatial neighbors of S_{s}, \( {\gamma}_{T1}^s \) and \( {\gamma}_{T2}^s \) are the temporal autoregressive coefficients, \( {\rho}_{Tm}^s \) are the spatial autoregressive coefficients, \( {\left({\sigma}_{T\omega}^s\right)}^2 \) is the variance of the i.i.d. random error \( {\upomega}_{Tij}^s \). The model parameters were estimated by maximum likelihood (ML) [27]. Schwarz’s Bayesian criterion (BIC) was used to indicate the fitness of the models. The Shapiro–Wilk test [28] was used for normality test. The analyses and modellings were executed in R. 3.4.4, using mainly stat4 library [29] and fda library [30].
Results
Seasonality at multiple district level
The seasonality of \( {IR}_{ij}^s \) was evidenced by the variations of the corresponding \( {SI}_i^s \) (Table 1). For all spatial strata, the maximum \( {SI}_i^s \) > 100% fell into two periods, between April and May and between September and October. These indicate that, for example for S_{2} in April, the median incidence proportion was 1.66 times higher than the average incidence proportion of the year. The minimum \( {SI}_i^s \) < 60% were in the period between December and January. These show that, for example for S_{2} in January, the median incidence proportion was 2.5 times lower than the average incidence proportion of the year.
As outcomes from the iterative fitting processes, the optimally fitted trends of S_{1} and S_{2} had two sine and two cosine functions (P = 2) and the optimal logλ = 4 corresponding to the minimum values of the GCV equal to 2.28 and 3.33, respectively. For S_{3}, the fitted model included four of each of the trigonometric functions (P = 4) with logλ = 2 and the minimum GCV equal to 6.1. The fitted models for the seasonality (3) have five of each sine and cosine functions (Q = 5) with logλ = 0 and the minimum GCV equal to 0.87, 2.6 and 4.1 for the S_{1}, S_{2} and S_{3}, respectively. The trends show a gradual decrease from 2012 to 2015, whereas from the end of 2015 onwards, they increase (Fig. 2). The magnitude is highest in S_{2}, being approximately double to those of S_{1} and S_{3}.
Fig 3 presents the seasonal patterns of \( {IR}_{ij}^s \). In T_{1}, all the peaks fell within April as indicated by the maximum \( {SI}_i^s \) in Table 1. In T_{3}, the peaks occurred in August (S_{3}), September (S_{2}) and October (S_{1}). All the deepest troughs were in January, whereas the small troughs occurred in June. The maximum amplitude of seasonality was equal to 4.56 infected U5s per 1,000 U5s at risk per month.
Driving factors of seasonality
The monthly relative humidity and sunny hours in Da Nang city in the study period were significantly correlated at lag zero with the average temperature (r_{lag0} = − 0.82 and 0.84, respectively). Meanwhile, \( \xi \left({IR}_{ij}^s\right) \) from all three spatial strata had significantly lower correlation with the rainfall than with the average temperature at lag zero (r_{lag0} = − 0.13, − 0.10, − 0.23 versus 0.74, 0.74, 0.63 for S_{1}, S_{2}, S_{3}, respectively). The average temperature was apparently the sole explanatory weather variable included or N = 1 in (4). The maximum average temperature of over 30 °C occurred in May and June; the minimum of approximately 22 °C occurred in January every year (Fig. 4). The \( {\xi}_T\left({IR}_{ij}^s\right) \) all satisfied the ShapiroWilk normality test (pvalue < 0.05), thus satisfying the specifications of (4) and (5) without any need for transformation.
The results of the qstatistic show that the spatialtemporal stratification into nine strata reflecting the influences of the density of U5s and the monthly average temperature yielded a significant homogeneity within the strata. The qstatistic of the nine spatialtemporal strata was equal to 0.21 (pvalue = 0.00), while that of the combinations of the months versus the districts was not significant, qstatistic = 0.26 (pvalue = 0.24).
To compare the effects of the average temperature and the preschooling period (Fig 4), the parameters in (4) and (5) were consequently estimated for two submodels: 1) with preschooling period or \( {\beta}_{T\left(N+1\right)}^s={\beta}_{T2}^s \) and 2) without preschooling period or \( {\beta}_{T\left(N+1\right)}^s=0 \). Accordingly, every spatialtemporal stratum had two sets of estimated parameters, except for those in T_{2} because children did not attend preschools in this period. The estimated results are presented in Tables 2 and 3. The effects of the spatialtemporal autocorrelations were taken into account in these estimations, resulting in estimated parameters as presented in Tables 4 and 5 in Appendix. The corresponding BIC values of the selected models were all minimized.
T_{2} is the transitional period between the dry and rainy season and also the period of summer holidays. Thus, the interest of studying the driving factor of the seasonality laid in T_{1} and T_{3}. Table 2 shows that in the submodel one, the average temperature was negatively correlated with \( \xi \left({IR}_{ij}^s\right) \) in T_{1}. This negative \( {\beta}_{T1}^s \) indicates a mismatch between the peak of the seasonality of the disease incidences and the average temperature. In submodel two as can be seen in Table 3, the \( {\beta}_{T2}^s \) had negative values for all strata in T_{3}; whereas, the \( {\beta}_{T1}^s \) were minimal. The average temperature therefore had no significant effect in T_{3}.
Discussion
In this section, we discuss our main findings, the consistency of the seasonal patterns extracted from \( {IR}_{ij}^s \) using our methods compared to other methods and we also highlight the limitations of this study.
In this study, we have decomposed the monthly time series of the adjusted incidence proportions of childhood HFMD infections for the spatialtemporal strata of the mainland of Da Nang city to reveal the trend and the seasonality. The yearly larger outbreaks simultaneously happened in all three spatial strata in April. The variation of the seasonality of the second smaller outbreaks in T_{3} is intriguing. This variation suggests the effects of the season of increasing social contacts among U5s when they go to preschool on the onset of HFMD infection.
We have shown statistical evidences of these effects during the preschooling period. The positive values of the spatialtemporal regression coefficients imply the simultaneity of the peaks and vice versa [21]. The mismatch between the maximum incidence proportions and the maximum average temperature indicates that the increase of the average temperature was the driving factor, not its maximum values. In the last two quarters of all years, the average temperature was not the predominant risk factor. The contrary was evident in T_{1}. At small areas of multiple districts in the study area, the increasing social contact has been shown to be the important driving factor of the geographical variation of the seasonality of the childhood HFMD outbreaks.
The estimated spatialtemporal regression parameters of (4) in T_{3} are interestingly more informative than those in T_{1} because all the smaller peaks occurred in T_{3}. As the period at the preschools in T_{3} cumulatively increased from August to December, the negative values of \( {\beta}_{T2}^s \) indicate the occurrence of the smaller peaks at the beginning of the preschooling period. Recall from the study of Horby et al. [17] that the physical contacts of U5s are most intense within this age group in Viet Nam. Given that the physical contact is the prominent means of passing HFMD viruses, the nature, frequency and duration of the contacts of those children during the preschooling period explain the occurrence and the geographical dynamics of these smaller outbreaks.
The trends show the occurrences of the large outbreaks at the beginning and at the end of the observed period. Study of a longer period of time would show a longterm oscillation of the large outbreaks with the frequency of every two to 3 years in this city. The possible reasons can be the required lapse of time to cumulate the critical susceptible population [31] or the existence of multiple viral strains in the study area [32, 33].
The decomposed trends and seasonal components were highly correlated with those fitted by the method of seasonal and trend decomposition using Loess [34]. The correlation coefficients of the trends from both methods at lag zero varied from 0.85 to 0.90. Those of the seasonality varied from 0.83 to 0.97. The similar results from both decomposition methods imply the consistency of the seasonality derived from the Fourier decomposition method. This method provides a semiparametric approach to extract different components of a time series including the longterm and shortterm irregular and regular fluctuations, taking into account effects of possible confounders [35, 36]. Using this approach instead of the common seasonal ARIMA, the drawbacks of stationarization by differencing and statistical transformation can be avoided [37]. Moreover, Fourier decomposition in combination with the spatialtemporal autoregressive model allow both spatial and temporal autocorrelation existing in the data to be included into the model calibration. This reduces the biases in the estimation of the associations between the incidence proportions and the risk factors. Notwithstanding, from a statistical point of view, the small number of spatial strata in the study area places one of the limitations of the study to include the effects of the spatial autocorrelation. In addition, the separation the temporal trend from the time series was based mainly upon expert judgments. In other words, the efficiency of the Fourier decomposition in many cases relies on the expert’s familiarity with the understudied phenomena.
Conclusions
HFMD has been becoming one of the most important pediatric NTDs in the Northeastern and Southeastern Asian countries. By applying spatialtemporal statistical analyses, the results have shown that at multiple district level, the social contact is the important driving factor of the spatial variation of the disease’s outbreaks. This study provides statistical evidences of the effects of the seasonality of increasing social contact amongst U5s on geographical dynamics of HFMD infection outbreaks. Our findings contribute to extend the understanding of the underlying driving factors of the disease dynamics at small areas. This contribution is necessary to inform the next insightful research into the spatialtemporal dynamics of the transmission of the NTDs within the local population.
Availability of data and materials
The dataset of HFMD incidences in Da Nang city analyzed during the current study is publicly available on the website of Da Nang city Preventive medicine center https://ksbtdanang.vn/. The dataset of the monthly average ambient temperature in Da Nang city is publicly available on the website of Viet Nam institute of meteorology, hydrology and climate change http://www.imh.ac.vn. The R code of the analyses and models are available from the corresponding author on reasonable request.
Abbreviations
 BIC:

Schwarz’s Bayesian criterion
 CA16:

Coxsackievirus A16
 EV71:

Human Enterovirus 71
 HFMD:

Handfootmouth disease
 ML:

Maximum likelihood
 NTDs:

Neglected tropical diseases
 S:

Spatial stratum
 SI:

Seasonal index
 STS:

Spatialtemporal statistical methods
 T:

Temporal stratum
 U5s:

Children under 5 years of age
References
 1.
Sultan B, Labadi K, Guégan JF, Janicot S. Climate drives the Meningitis epidemics onset in West Africa. PLoS Med. 2005;2:e6. https://doi.org/10.1371/journal.pmed.0020006.
 2.
Do TTT, Martens P, Luu NH, Wright P, Choisy M. Climaticdriven seasonality of emerging dengue fever in Hanoi, Vietnam. BMC Public Health. 2014;14:1078. https://doi.org/10.1186/14712458141078.
 3.
Kim BI, Ki H, Park S, Cho E, Chun BC. Effect of climatic factors on hand, foot, and mouth disease in South Korea, 2010–2013. PLoS One. 2016;11:e0157500. https://doi.org/10.1371/journal.pone.0157500.
 4.
Pascual M, Dobson A. Seasonal patterns of infectious diseases. PLoS Med. 2005;2:e5. https://doi.org/10.1371/journal.pmed.0020005.
 5.
Altizer S, Dobson A, Hosseini P, Hudson P, Pascual M, Rohani P. Seasonality and the dynamics of infectious diseases. Ecol Lett. 2006;9:467–84.
 6.
Grassly NC, Fraser C. Seasonal infectious disease epidemiology. Proc R Soc B Biol Sci. 2006;273:2541–50.
 7.
United Nations. SDGsGoal 3: Ensure healthy lives and promote wellbeing for all at all ages. http://www.un.org/sustainabledevelopment/health. Accessed 18 Jan 2018.
 8.
Hotez PJ, Bottazzi ME, Strych U, Chang LY, Lim YAL, Goodenow MM, AbuBakar S. Neglected tropical diseases among the association of Southeast Asian Nations (ASEAN): Overview and update. PLoS Negl Trop Dis. 2015;9:e0003575. https://doi.org/10.1371/journal.pntd.0003575.
 9.
Koh WM, Bogich T, Siegel K, Jin J, Chong EY, Tan CY, Chen MIC, Horby P, Cook AR. The epidemiology of hand, foot and mouth disease in Asia: a systematic review and analysis. Pediatr Infect Dis J. 2016;35:e285–300. https://doi.org/10.1097/INF.0000000000001242.
 10.
Frydenberg A, Starr M. Hand, foot and mouth disease. Aust Fam Physician. 2003;32:594–5.
 11.
Liu W, Ji H, Shan J, Bao J, Sun Y, Li J, Bao C, Tang F, Yang K, Bergquist R, Peng Z, Zhu Y. Spatiotemporal dynamics of handfootmouth disease and its relationship with meteorological factors in Jiangsu province, China. PLoS One. 2015;10:e0131311. https://doi.org/10.1371/journal.pone.0131311.
 12.
Wang J, Cao Z, Zeng DD, Wang Q. Assessing local risk factors of Beijing handfoot mouth disease in China. Online J Public Health Inform. 2017;9(1):e009. https://doi.org/10.5210/ojphi.v9i1.7759.
 13.
Xu C, Zhang X, Xiao G. Spatiotemporal decomposition and risk determinants of hand, foot and mouth disease in Henan, China. Sci Total Environ. 2019;657:50916.
 14.
NikNadia NMN, Sam IC, Rampal S, WanNorAmalina WMZ, NurAtifah G, Verasahib K, Ong CC, MohdAdib M, Chan YF. Cyclical patterns of hand, foot and mouth disease caused by Enterovirus A71 in Malaysia. PLoS Negl Trop Dis. 2016;10:e0004562. https://doi.org/10.1371/journal.pntd.0004562.
 15.
Phung D, Nguyen HX, Nguyen HLT, Do CM, Tran QD, Chu C. Spatiotemporal variation of handfootmouth disease in relation to socioecological factors: a multipleprovince analysis in Vietnam. Sci Total Environ. 2018;610:983–91.
 16.
Sofianopoulou E, PlessMulloli T, Rushton S, Diggle PJ. Modeling seasonal and spatiotemporal variation: the example of respiratory prescribing. Am J Epidemiol. 2017;186:101–8.
 17.
Horby P, Thai PQ, Hens N, Yen NTT, Mai LQ, Thoang DD, Linh NM, Huong NT, Alexander N, Edmunds WJ, Duong TN, Fox A, Hien NT. Social contact patterns in Vietnam and implications for the control of infectious diseases. PLoS One. 2011;6:e16965. https://doi.org/10.1371/journal.pone.0016965.
 18.
Da Nang Statistics Office. Da Nang statistical yearbook 2016. Da Nang: Statistical publishing house; 2017.
 19.
Wang JF, Zhang TL, Fu BJ. A measure of spatial stratified heterogeneity. Ecol Indic. 2016;67:250–6.
 20.
Wang JF, Li XH, Christakos G, Liao YL, Zhang T, Gu X, Zheng XY. Geographical detectorsbased health risk assessment and its application in the neural tube defects study of the Heshun region. China Int J Geogr Inf Sci. 2010;24:107–27.
 21.
Barnett AG, Dobson AJ. Analysing seasonal health data. Berlin: Springer; 2010.
 22.
Dodge Y. Seasonal index. In: Dodge Y, editor. The concise encyclopedia of statistics. New York: Springer New York; 2008. p. 476–9.
 23.
Bloomfield P. Fourier analysis of time series: an introduction. 2nd ed. New York: Wiley; 2000.
 24.
Naumova EN, MacNeill IB. Seasonality assessment for biosurveillance systems. In: Auget JL, Balakrishnan N, Mesbah M, Molenberghs G, editors. Advances in statistical methods for the health sciences: applications to cancer and aids studies, genome sequence analysis, and survival analysis. Boston: Birkhäuser Boston; 2007. p. 437–50.
 25.
Ramsay JO, Hooker G, Graves S. Smoothing: computing curves from noisy data. In: Ramsay JO, Hooker G, Graves S, editors. Functional data analysis with R and MATLAB. New York: Springer New York; 2009. p. 59–82.
 26.
Pace RK, Barry R, Clapp J, Rodriquez M. Spatiotemporal autoregressive models of neighborhood effects. J Real Estate Finance Econ. 1998;17:15–33.
 27.
Cole SR, Chu H, Greenland S. Maximum likelihood, profile likelihood, and penalized likelihood: a primer. Am J Epidemiol. 2014;179:252–60.
 28.
Yap BW, Sim CH. Comparisons of various types of normality tests. J Stat Comput Simul. 2011;81:2141–55.
 29.
R Core Team. R: A language and environment for statistical computing. Vienna: R Foundation for Statistical Computing; 2017. https://www.Rproject.org
 30.
Ramsay JO, Wickham H, Graves S, Hooker G. fda: Functional Data Analysis. R package version 2.4.8. 2018. https://cran.rproject.org/web/packages/fda/index.html.
 31.
Anderson RM. Populations and infectious diseases: ecology or epidemiology? J Anim Ecol. 1991;60:1–50.
 32.
Sarma N. Relapse of hand foot and mouth disease: are we at more risk? Indian J Dermatol. 2013;58:78–9.
 33.
Thao NTT, Donato C, Trang VTH, Kien NT, Trang PMT, Khanh TQ, Nguyet DT, Sessions OM, Cuong HQ, Lan PT, Huong VTQ, van Doorn HR, Vijaykrishna D. Evolution and spatiotemporal dynamics of Enterovirus A71 subgenogroups in Vietnam. J Infect Dis. 2017;216:1371–9.
 34.
Cleveland RB, Cleveland WS, McRae JE, Terpenning I. STL: a seasonaltrend decomposition procedure based on loess. J Off Stat. 1990;6:3–73.
 35.
Diggle P. Time series: a biostatistical introduction. Oxford statistical science series. Oxford: Oxford University Press; 1990.
 36.
Stolwijk AM, Straatman H, Zielhuis GA. Studying seasonality by using sine and cosine functions in regression analysis. J Epidemiol Community Health. 1999;53:235–8.
 37.
Makridakis S, Hibon M. ARMA models and the box–Jenkins methodology. J Forecast. 1997;16:147–63.
Acknowledgements
We would like to extend our thanks to the Faculty of Disease control and prevention, Pasteur Institute in Ho Chi Minh city, Viet Nam.
Funding
Not applicable.
Author information
Affiliations
Contributions
PNT contributed to the conception, data analysis and interpretation of the results. AS contributed to the conception and interpretation of the results. TVN and TTTN contributed to the interpretation of the results. PNT drafted the manuscript and all authors were involved in its review and revision. All authors read and approved the final manuscript.
Corresponding author
Correspondence to Phuong N. Truong.
Ethics declarations
Ethics approval and consent to participate
Not applicable.
Consent for publication
Not applicable.
Competing interests
The authors declare that they have no competing interests.
Additional information
Publisher’s Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Appendix
Appendix
Rights and permissions
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
About this article
Received
Accepted
Published
DOI
Keywords
 Health seasonality
 NTDs
 Childhood HFMD
 STAR
 Fourier analysis
 Disease
 dynamics, health geography, small area analysis.