#### Abstract

To investigate the variability of HFMD in each county of Wenzhou, a spatial-temporal ARMA model is presented, and a general Bayesian framework is given for parameter estimation. The proposed model has two advantages: (i) allowing time series to be correlated, thus it can describe the series both spatially and temporally; (ii) implementing forecast easily. Based on the HFMD data in Wenzhou, we find that HFMD had positive spatial autocorrelation and the incidence seasonal peak was between May and July. In the county-level analysis, we find that after first-order difference the spatial-temporal ARMA model provides an adequate fit to the data.

#### 1. Introduction

Hand, foot, and mouth disease (HFMD) is a common infectious disease which usually affects children, particularly those less than 5 years old and younger. It is characterized by a distinct clinical presentation of fever or vesicular exanthema on their hands, feet, mouths, or buttocks [1–5]. The transmission of HFMD occurs from person to person through direct contact with saliva, faeces, vesicular fluid, or respiratory droplets of an infected person and indirectly by contaminated articles [1]. After a susceptible individual is infected he firstly enters the incubation period of HFMD, which is about 3 to 7 days. After the incubation period, the infected will show some clinical symptoms, such as having a fever, poor appetite, malaise, and sore throat, and few people may develop dehydration, febrile seizures, encephalitis, meningitis, cardiomyopathy, and so forth. And the infected people will fully recover after 7 to 10 days [1]. At present, there are still no available effective vaccines or drugs against HFMD human use, but such vaccines are being developed [6]. In 2012, for instance, an epidemic in mainland China involved 2,168,737 cases and 567 deaths [2]. Particularly, in 2012, there were 147,941 HFMD cases and 17 deaths in Zhejiang province, and it ranks the first in the “Ten legal infectious disease” [7]. HFMD has become an emerging public health concern in the affected countries and a focus of increasing amounts of research [4]. Therefore, it is important to use mathematical models to improve our understanding of infectious disease dynamics of HFMD and to help us develop preventive measures to control infection spread qualitatively and quantitatively.

There are several types of analytical models that are valuable to understand and predict the transmission of HFMD. One is compartmental differential equation model [8–15], which is important to understand the spread dynamics of HFMD among the susceptible populations and to enable policy makers to take effective measures to curb the disease spread and reduce the adverse impact of the disease [9, 10]. The other is statistical models which can help us find novel information concerning pathogen detection and some probable coinfection factors in HFMD and have been applied to understand HFMD’s spatiotemporal transmission and discover the relationship between HFMD occurrence and climate [3, 16–29]. Of them, Hu et al. [4] explored the spatial association of HFMD incidence with several potential determinants and found that child population density and climate factors are potential determinants of the HFMD incidence in most areas in China. Liu et al. [5] conducted a spatial and space-time scan statistics analysis in Shandong province to explore the distribution characteristics and detect spatial and spatial-temporal clusters (hotspots) of HFMD cases and found that the incidence seasonal peak was between April and July. Deng et al. [28] analyzed the epidemic characteristic of HFMD in Guangdong province and detected spatial-temporal clusters and found that climate factors and demographic changes might be contributors affecting the epidemic situation of HFMD. All these studies increase our understanding of the distribution and severity of the disease. However, potential factors influencing the incidence of HFMD remain little understood.

Wenzhou is a prefecture-level city in southeastern Zhejiang province in China. At the time of the 2010 Chinese census, 9,122,100 people lived in Wenzhou [30]. It includes 3 municipal districts (Lucheng, Longwan, and Ouhai) and 8 counties (Cangnan, Dongtou, Pingyang, Ruian, Taishun, Wencheng, Yongjia, and Yueqing) with a total land area of 11,784 square kilometers (Figure 1). Since Wenzhou has a humid subtropical climate zone with an annual average C (F), it is of particular public health significance to update molecular epidemiology of HFMD in Wenzhou. It is reported that, in 2012, there were 41,438 HFMD cases and 17 deaths in Wenzhou [7].

In [15], the authors established an SEIQRS epidemic model with periodic transmission rate to investigate the spread of seasonal HFMD in Wenzhou and found that the HFMD becomes an endemic disease in Wenzhou and for controlling the HFMD spread, it is beneficial to increase the quarantined rate or decrease the treatment cycle. However, spatial effect has not been considered in [15]. And in [4, 5, 28], the authors utilized spatial models to determine the risk factors of the incidence of HFMD and found that monthly average temperature, relative humidity, and total sunshine were important factors to affect the incidence of HFMD. However, their models have low power of prediction because future incidence rate is dependent on future risk factors.

The goal of this paper is to explore a model that could describe spatial-temporal effects and that has the ability to forecast. Unlike [4, 5, 28], we do not take the risk factors into account. Based on the monthly observed data of HFMD in each county of Wenzhou [7], we present a spatial-temporal autoregressive moving average (ARMA) model and give a general Bayesian framework for parameter estimation.

The paper is organized as follows. In Section 2, we give a method to preprocess the original data and propose a spatial-temporal ARMA model. In Section 3, the HFMD data from the Wenzhou Center for Disease Control and Prevention is analyzed based on descriptive statistics and spatial-temporal ARMA model. Finally, we give a brief discussion.

#### 2. Model Derivation

##### 2.1. Data Collection and Preprocessing

Data of HFMD in Wenzhou during 2010 to 2012 were collected by Wenzhou Center for Disease Control and Prevention, including the basic demographic characteristics of HFMD cases, the pathogen type (EV71, CoxA16, and other EV) of some HFMD cases, and incident child cases in each county per month. According to the records obtained, all the cases were children aged between 0 years and 14 years, with county-specific case numbers varying from 0 to 1852 incident child cases.

The HFMD dataset lists the number of incident cases in each county per month. It reflects the occurrence of the disease in different regions. However, it cannot reflect the risk of infecting the disease because of the different population sizes among counties or municipal districts. To reduce the influence of population size, cumulative incidence (CI) is utilized to reflect the risks of infecting HFMD in each county. It measures the disease frequency during a period of time [31]. We denote , , and as the CI of HFMD, the number of cases, and the child population size at the time in the th county, respectively. Then . In our case, . However, cannot be used to compare the disease risk between different counties directly because of random effects. Some counties reported 0 HFMD cases in some months; however, it could not say that there is no risk of HFMD in these counties. In fact, they had HFMD cases in other months. Thus, we use a hierarchical Bayes model to adjust CI [32]. The spirit of the idea to improve accuracy in the model is smoothing strength among counties. The model is as follow: where is the overall level of the disease risk at time , is the uncorrelated heterogeneity following normal distribution with mean 0 and variance , and is the component to describe spatial correlation, which is defined by the intrinsic Gaussian autoregression distribution [33, 34]. That is, where is the spatial adjacent matrix defining the connectivity between counties. if the th county and the th county are adjacent; otherwise . To perform Bayesian analysis, we assign gamma distributions with a large variance as the priors for the parameters and , and the trick is mostly used in Bayesian spatial analysis [34]. Thus, the priors for and are given by The estimation of the parameters can be obtained by MCMC, and we use as the adjusted CI, where and are the posterior mean of the parameters and . Since , modeling such data is unstable and the fitted values may exceed the interval (0,1). We make the transformation ; then . The data will be utilized to model the risk of HFMD in Wenzhou.

##### 2.2. Spatial-Temporal ARMA Model

Let , where is an known function, and we assume is a stationary time series. For fixed , we assume follows a multiplicative seasonal model with seasonal period as a model with autoregressive characteristic polynomial and moving average characteristic polynomial , where For better expressing the model, the backshift operator, denoted , is needed. operates on the time index of a series and shifts time back one time unit to form a new series. In particular, . Thus, the model can be formulated as where is a white noise process with variance . However, model (5) does not introduce spatial correlation among different regions. Thus, we modify the model as follows: where is defined in (2). Model (6) is much different from [35, 36]. We introduce to describe spatial variability rather than correlated . The advantage of this is that we can estimate the parameters in the model easily. Let , , , , and . Then the parameters in the model (6) are . Denote , , , and . We denote that the likelihood function of given and is . can be obtained by several methods. One of the most used methods is Kalman filter. The overall strategy for computing the likelihood for a given set of parameter values is to use the Kalman filter equations to generate recursively the prediction errors and their variances and then use the prediction error decomposition of the likelihood function [37].

Prior distribution of can be elicited from the experts’ opinions or historical data. When the historical data is available, power prior is a good choice, which has emerged as a useful class of informative priors for historical data [38]. Assume that the joint prior distribution of is
is following the intrinsic Gaussian autoregression distribution. Then the density function of is defined as
where is a neighborhood of the th county. Thus the joint density of is
Based on the likelihood function , (7), and (9), the posterior distribution of is
With specific priors of , **,** and , we can use the corresponding method to generate the posterior samples of . In the case study, we use the normal distributions as the priors for and , and the Gamma distributions as the prior for . Then the full conditional posterior distributions of , , and are all log-concave densities. Thus adaptive rejection sampling can be used for posterior sampling [39]. In general, Metropolis-Hastings algorithm can be used to implement the Markov chain sampling. After running Markov chain sampling procedure times and discarding the initial burn-in iterations, then we have iterations kept. Since the generated sample is not independent, we need to monitor the autocorrelations of the generated values and select a sampling lag after which the corresponding autocorrelation is low; that is, the length of the thinning interval is . Considering the length of the thinning interval, the final number of iterations kept is , and these independent samples will be used for posterior analysis. Assume the generated sample is
Then for any function of the parameters of interest we can(1)obtain a sample of the desired parameters by simply considering
(2)estimate the posterior mean of by . And the posterior standard deviation is estimated by . Other measures of interest might be the posterior median or quantiles (e.g., 2.5% and 97.5% percentiles will provide a 95% credible interval).

#### 3. Results and Interpretation

##### 3.1. Prevalence of HFMD

Table 1 lists the demographic characteristics of HFMD cases and the pathogen types of some cases from 2010 to 2012. There are a total of 103,671 HFMD cases reported in Wenzhou. The incidence rates are 2,111.1 per 100,000 (18 years old or younger) in 2010, 1,872.0 per 100,000 (18 years old or younger) in 2011, and 2,702.3 per 100,000 (18 years old or younger) in 2012. The majority of the reported HFMD cases are less than 5 years old (about 89.4%), which means that younger children are more vulnerable to HFMD.

Most of HFMD cases (approximately 62.3%) were boys and the male-to-female incidence ratio was 1.638 : 1 in 2010, 1.691 : 1 in 2011, and 1.633 : 1 in 2012, respectively. The values of testing the difference between incidence in male and that in female were smaller than 0.001. Thus, the incidence in male was higher than that in female. Among 2,401 pathogen study between 2010 and 2012, CoxA16, EV71, and other EV accounted for 10.08%, 73.26%, and 16.66%, respectively. Clearly, the majority of HFMD pathogens were EV71.

Figure 2 shows monthly reported cases of HFMD in Wenzhou from 2010 to 2012. From Figure 2, we can see that incidence of HFMD is periodic and that a peak attains between May and July, because the temperature will affect the incidence of HFMD as studied in [4, 5, 28].

##### 3.2. Spatial Autocorrelation of HFMD Cases

Figure 3 shows the annual incidence rate per county accounting for the spatial variability of population size. It clearly indicates that the distribution of HFMD is heterogeneous at county level, and Ouhai municipal district is the most severe region. Table 2 lists the results of the spatial autocorrelation test, which demonstrates that high global spatial autocorrelation of HFMD was detected at county level in Wenzhou from 2010 to 2012 (Moran’s , values ).

**(a)**

**(b)**

**(c)**

##### 3.3. Spatial-Temporal ARMA Model of HFMD

For each county, we test the stationarity of and find that all the series are nonstationary. Then first-order difference is taken and the adjusted Dickey-Fuller test shows that the time series is stationary when the significant level is 0.1 (Table 3). Thus , where . As shown in Figure 2, we can specify the seasonal period . Then autocorrelation function and partial autocorrelation function are used to determine the values of , and initially. That is, , and . We compared all the submodels by AIC criteria [40] and chose the best model with the smallest AIC. The optimal model with the smallest AIC is , and . Thus the model used in the data analysis is where .

Table 4 lists the estimates of and in each county, where the numbers in the parentheses are the posterior standard deviations of the parameters. Figure 4 shows the comparison between the data and the fitted values, where “O” means the observed data and “F” is the fitted values.

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

**(f)**

**(g)**

**(h)**

**(i)**

**(j)**

**(k)**

#### 4. Conclusions and Discussions

In our study, we observed that children ≤5 years old accounted for most of the HFMD cases during the study period in Wenzhou, with an average male-to-female sex ratio 1.65 and the incidence single peak season was between May and July (see Table 1 and Figure 2), which was similar to other studies [3–5, 26, 28]. The dominant pathogen was EV71 (73.26%), which was significantly higher than other districts of China. The difference might be partly attributed to climatic, geographic, social factors, and so forth [24, 41, 42].

In Wenzhou, HFMD had positive spatial autocorrelation at county level with high Morans *I* with value . The highest incidence area is Ouhai, which is located in the middle of Wenzhou, and the second highest group areas are Lucheng, Longwan, and Ruian that are adjacent to Ouhai (see Figure 3). To investigate the spatial variability, we present a spatial-temporal ARMA model, in which we describe the spatial effects by introducing correlated random effects. The model provides an adequate fit to the data of HFMD in Wenzhou. Other models can also be considered, for example,
where is intrinsic Gaussian autoregression distribution. Model (14) introduces the spatial effects by directly. However, it will make the parameter estimation much more difficult, since the number of parameters increases. How to analyze model (14) needs further study.

#### Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

#### Acknowledgments

This research is supported by the National Science Foundation of China (61373005 and 11201345) and Zhejiang Provincial Natural Science Foundation (LY12A01014).