# A Spatial-Temporal ARMA Model of the Incidence of Hand, Foot, and Mouth Disease in Wenzhou, China

Authors: Jie Li; Yanjun Fu; Ancha Xu; Zumu Zhou; Weiming Wang
Journal: Abstract and Applied Analysis
 (2014)

---

## 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(
0
,
0
)
×
(
1
,
0
)
12 model provides an
adequate fit to the data.


---

## Body

## 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 18.08
°C (64.5
°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].Figure 1
The location of Wenzhou in China.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 Section2, 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 P
t
j, Y
t
j, and N
t
j as the CI of HFMD, the number of cases, and the child population size at the time t in the jth county, respectively. Then P
t
j
=
Y
t
j
/
N
t
j. In our case, t
=
1,2
,
…
,
36
,
j
=
1,2
,
…
,
11. However, P
t
j 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:
(1)
P
t
j
=
α
t
+
μ
t
j
+
ν
t
j
,
t
=
1,2
,
…
,
36
,
j
=
1,2
,
…
,
11
,
where α
t is the overall level of the disease risk at time t, ν
t
j is the uncorrelated heterogeneity following normal distribution with mean 0 and variance σ
ν
2, and μ
t
j is the component to describe spatial correlation, which is defined by the intrinsic Gaussian autoregression distribution [33, 34]. That is,
(2)
ν
t
j
∣
σ
ν
2
~
N
(
0
,
σ
ν
2
)
,
μ
t
j
∣
μ
t
i
,
i
≠
j
,
σ
μ
2
~
N
(
∑
i
μ
t
i
w
i
j
∑
i
w
i
j
,
σ
μ
2
∑
i
w
i
j
)
,
where W
=
(
w
i
j
) is the spatial adjacent matrix defining the connectivity between counties. w
i
j
=
1 if the ith county and the jth county are adjacent; otherwise w
i
j
=
0. To perform Bayesian analysis, we assign gamma distributions with a large variance as the priors for the parameters σ
ν
2 and σ
μ
2, and the trick is mostly used in Bayesian spatial analysis [34]. Thus, the priors for σ
ν
2 and σ
μ
2 are given by
(3)
1
σ
ν
2
~
Gamma
(
0.01,0.01
)
,
1
σ
μ
2
~
Gamma
(
0.1,0.0001
)
.
The estimation of the parameters can be obtained by MCMC, and we use R
^
t
j
=
α
^
t
+
μ
^
t
j as the adjusted CI, where α
^
t and μ
^
t
j are the posterior mean of the parameters α
t and μ
^
t
j. Since 0
<
R
^
t
j
<
1, modeling such data is unstable and the fitted values may exceed the interval (0,1). We make the transformation S
t
j
=
log
⁡
(
R
^
t
j
/
(
1
-
R
^
t
j
)
); then -
∞
<
S
t
j
<
∞. The data S
t
j will be utilized to model the risk of HFMD in Wenzhou.


### 2.2. Spatial-Temporal ARMA Model

Letf
t
j
=
f
(
S
t
j
), where f
(
·
) is an known function, and we assume {
f
t
j
,
t
=
1,2
,
…
,
36
} is a stationary time series. For fixed j, we assume {
f
t
j
,
t
=
1,2
,
…
,
36
} follows a multiplicative seasonal ARMA
(
p
,
q
)
×
(
P
,
Q
)
s model with seasonal period s as a model with autoregressive characteristic polynomial ϕ
j
(
x
)
Φ
j
(
x
) and moving average characteristic polynomial θ
j
(
x
)
Θ
j
(
x
), where
(4)
ϕ
j
(
x
)
=
1
-
ϕ
1
j
x
-
ϕ
2
j
x
2
-
⋯
-
ϕ
p
j
x
p
,
Φ
j
(
x
)
=
1
-
Φ
1
j
x
s
-
Φ
2
j
x
2
s
-
⋯
-
Φ
P
j
x
P
s
,
θ
j
(
x
)
=
1
-
θ
1
j
x
-
θ
2
j
x
2
-
⋯
-
θ
q
j
x
q
,
Θ
j
(
x
)
=
1
-
Θ
1
j
x
s
-
Θ
2
j
x
2
s
-
⋯
-
Θ
Q
j
x
Q
s
.
For better expressing the model, the backshift operator, denoted B, is needed. B operates on the time index of a series and shifts time back one time unit to form a new series. In particular, B
f
t
j
=
f
(
t
-
1
)
j. Thus, the model can be formulated as
(5)
ϕ
j
(
B
)
Φ
j
(
B
)
f
t
j
=
θ
j
(
B
)
Θ
j
(
B
)
e
t
j
,
j
=
1,2
,
…
,
11
,
where {
e
t
j
,
t
=
1,2
,
…
,
36
} is a white noise process with variance σ
j
2. However, model (5) does not introduce spatial correlation among different regions. Thus, we modify the model as follows:
(6)
ϕ
j
(
B
)
Φ
j
(
B
)
f
t
j
=
θ
j
(
B
)
Θ
j
(
B
)
e
t
j
+
μ
t
j
,
j
=
1,2
,
…
,
11
,
where μ
t
j is defined in (2). Model (6) is much different from [35, 36]. We introduce μ
t
j to describe spatial variability rather than correlated e
t
j. The advantage of this is that we can estimate the parameters in the model easily. Let ϕ
.
j
=
(
ϕ
1
j
,
…
,
ϕ
p
j
,
Φ
1
j
,
…
,
Φ
P
j
), θ
.
j
=
(
θ
1
j
,
…
,
θ
p
j
,
Θ
1
j
,
…
,
Θ
P
j
), ϕ
=
(
ϕ
.
1
,
…
,
ϕ
.
11
), θ
=
(
θ
.
1
,
…
,
θ
.
11
), and σ
2
=
(
σ
1
2
,
…
,
σ
11
2
). Then the parameters in the model (6) are ω
=
(
ϕ
,
θ
,
σ
2
). Denote f
j
=
(
f
1
j
,
…
,
f
36
j
), μ
1
=
(
μ
11
,
μ
12
,
…
,
μ
111
), f
=
(
f
1
,
…
,
f
11
), and μ
=
(
μ
1
,
μ
2
,
…
,
μ
36
). We denote that the likelihood function of f given ω and μ is L
(
f
∣
ω
,
μ
). L
(
f
∣
ω
,
μ
) 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
(7)
π
(
ω
)
∝
π
ϕ
(
ϕ
)
π
θ
(
θ
)
π
σ
(
σ
2
)
.
μ
t is following the intrinsic Gaussian autoregression distribution. Then the density function of μ
t is defined as
(8)
p
(
μ
t
∣
σ
μ
2
)
∝
1
σ
μ
1
1
exp
⁡
{
-
1
σ
μ
2
∑
i
‍
∑
j
∈
δ
i
‍
(
μ
t
i
-
μ
t
j
)
2
}
,
where δ
i is a neighborhood of the ith county. Thus the joint density of μ is
(9)
p
(
μ
∣
σ
μ
2
)
∝
∏
t
=
1
36
‍
p
(
μ
t
∣
σ
μ
2
)
.
Based on the likelihood function L
(
f
∣
ω
,
μ
), (7), and (9), the posterior distribution of ω is
(10)
p
(
ω
∣
f
,
μ
,
σ
μ
2
)
∝
L
(
f
∣
ω
,
μ
)
π
ϕ
(
ϕ
)
π
θ
(
θ
)
π
σ
(
σ
2
)
×
∏
t
=
1
36
‍
p
(
μ
t
∣
σ
μ
2
)
.
With specific priors of ϕ, θ
, and σ
2, 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 σ
2. Then the full conditional posterior distributions of ϕ, θ, and σ
2 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 M times and discarding the initial A burn-in iterations, then we have (
M
-
A
) iterations kept. Since the generated sample is not independent, we need to monitor the autocorrelations of the generated values and select a sampling lag L
>
1 after which the corresponding autocorrelation is low; that is, the length of the thinning interval is L. Considering the length of the thinning interval, the final number of iterations kept is M
′
=
(
M
-
A
)
/
L, and these independent samples will be used for posterior analysis. Assume the generated sample is
(11)
(
ω
(
1
)
,
ω
(
2
)
,
…
,
ω
(
M
′
)
)
.
Then for any function G
(
ω
) of the parameters of interest ω we can(1)
obtain a sample of the desired parametersG
(
ω
) by simply considering
(12)
(
G
(
ω
(
1
)
)
,
G
(
ω
(
2
)
)
,
…
,
G
(
ω
(
M
′
)
)
)
,
(2)
estimate the posterior mean ofG
(
ω
) by G
(
ω
)
¯
=
(
1
/
M
′
)
∑
s
=
1
M
′
G
(
ω
(
t
)
). And the posterior standard deviation is estimated by (
1
/
(
M
′
-
1
)
)
∑
s
=
1
N
′
[
G
(
ω
(
t
)
)
-
G
(
ω
)
¯
]
2. 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).


## 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 P
t
j, Y
t
j, and N
t
j as the CI of HFMD, the number of cases, and the child population size at the time t in the jth county, respectively. Then P
t
j
=
Y
t
j
/
N
t
j. In our case, t
=
1,2
,
…
,
36
,
j
=
1,2
,
…
,
11. However, P
t
j 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:
(1)
P
t
j
=
α
t
+
μ
t
j
+
ν
t
j
,
t
=
1,2
,
…
,
36
,
j
=
1,2
,
…
,
11
,
where α
t is the overall level of the disease risk at time t, ν
t
j is the uncorrelated heterogeneity following normal distribution with mean 0 and variance σ
ν
2, and μ
t
j is the component to describe spatial correlation, which is defined by the intrinsic Gaussian autoregression distribution [33, 34]. That is,
(2)
ν
t
j
∣
σ
ν
2
~
N
(
0
,
σ
ν
2
)
,
μ
t
j
∣
μ
t
i
,
i
≠
j
,
σ
μ
2
~
N
(
∑
i
μ
t
i
w
i
j
∑
i
w
i
j
,
σ
μ
2
∑
i
w
i
j
)
,
where W
=
(
w
i
j
) is the spatial adjacent matrix defining the connectivity between counties. w
i
j
=
1 if the ith county and the jth county are adjacent; otherwise w
i
j
=
0. To perform Bayesian analysis, we assign gamma distributions with a large variance as the priors for the parameters σ
ν
2 and σ
μ
2, and the trick is mostly used in Bayesian spatial analysis [34]. Thus, the priors for σ
ν
2 and σ
μ
2 are given by
(3)
1
σ
ν
2
~
Gamma
(
0.01,0.01
)
,
1
σ
μ
2
~
Gamma
(
0.1,0.0001
)
.
The estimation of the parameters can be obtained by MCMC, and we use R
^
t
j
=
α
^
t
+
μ
^
t
j as the adjusted CI, where α
^
t and μ
^
t
j are the posterior mean of the parameters α
t and μ
^
t
j. Since 0
<
R
^
t
j
<
1, modeling such data is unstable and the fitted values may exceed the interval (0,1). We make the transformation S
t
j
=
log
⁡
(
R
^
t
j
/
(
1
-
R
^
t
j
)
); then -
∞
<
S
t
j
<
∞. The data S
t
j will be utilized to model the risk of HFMD in Wenzhou.


## 2.2. Spatial-Temporal ARMA Model

Letf
t
j
=
f
(
S
t
j
), where f
(
·
) is an known function, and we assume {
f
t
j
,
t
=
1,2
,
…
,
36
} is a stationary time series. For fixed j, we assume {
f
t
j
,
t
=
1,2
,
…
,
36
} follows a multiplicative seasonal ARMA
(
p
,
q
)
×
(
P
,
Q
)
s model with seasonal period s as a model with autoregressive characteristic polynomial ϕ
j
(
x
)
Φ
j
(
x
) and moving average characteristic polynomial θ
j
(
x
)
Θ
j
(
x
), where
(4)
ϕ
j
(
x
)
=
1
-
ϕ
1
j
x
-
ϕ
2
j
x
2
-
⋯
-
ϕ
p
j
x
p
,
Φ
j
(
x
)
=
1
-
Φ
1
j
x
s
-
Φ
2
j
x
2
s
-
⋯
-
Φ
P
j
x
P
s
,
θ
j
(
x
)
=
1
-
θ
1
j
x
-
θ
2
j
x
2
-
⋯
-
θ
q
j
x
q
,
Θ
j
(
x
)
=
1
-
Θ
1
j
x
s
-
Θ
2
j
x
2
s
-
⋯
-
Θ
Q
j
x
Q
s
.
For better expressing the model, the backshift operator, denoted B, is needed. B operates on the time index of a series and shifts time back one time unit to form a new series. In particular, B
f
t
j
=
f
(
t
-
1
)
j. Thus, the model can be formulated as
(5)
ϕ
j
(
B
)
Φ
j
(
B
)
f
t
j
=
θ
j
(
B
)
Θ
j
(
B
)
e
t
j
,
j
=
1,2
,
…
,
11
,
where {
e
t
j
,
t
=
1,2
,
…
,
36
} is a white noise process with variance σ
j
2. However, model (5) does not introduce spatial correlation among different regions. Thus, we modify the model as follows:
(6)
ϕ
j
(
B
)
Φ
j
(
B
)
f
t
j
=
θ
j
(
B
)
Θ
j
(
B
)
e
t
j
+
μ
t
j
,
j
=
1,2
,
…
,
11
,
where μ
t
j is defined in (2). Model (6) is much different from [35, 36]. We introduce μ
t
j to describe spatial variability rather than correlated e
t
j. The advantage of this is that we can estimate the parameters in the model easily. Let ϕ
.
j
=
(
ϕ
1
j
,
…
,
ϕ
p
j
,
Φ
1
j
,
…
,
Φ
P
j
), θ
.
j
=
(
θ
1
j
,
…
,
θ
p
j
,
Θ
1
j
,
…
,
Θ
P
j
), ϕ
=
(
ϕ
.
1
,
…
,
ϕ
.
11
), θ
=
(
θ
.
1
,
…
,
θ
.
11
), and σ
2
=
(
σ
1
2
,
…
,
σ
11
2
). Then the parameters in the model (6) are ω
=
(
ϕ
,
θ
,
σ
2
). Denote f
j
=
(
f
1
j
,
…
,
f
36
j
), μ
1
=
(
μ
11
,
μ
12
,
…
,
μ
111
), f
=
(
f
1
,
…
,
f
11
), and μ
=
(
μ
1
,
μ
2
,
…
,
μ
36
). We denote that the likelihood function of f given ω and μ is L
(
f
∣
ω
,
μ
). L
(
f
∣
ω
,
μ
) 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
(7)
π
(
ω
)
∝
π
ϕ
(
ϕ
)
π
θ
(
θ
)
π
σ
(
σ
2
)
.
μ
t is following the intrinsic Gaussian autoregression distribution. Then the density function of μ
t is defined as
(8)
p
(
μ
t
∣
σ
μ
2
)
∝
1
σ
μ
1
1
exp
⁡
{
-
1
σ
μ
2
∑
i
‍
∑
j
∈
δ
i
‍
(
μ
t
i
-
μ
t
j
)
2
}
,
where δ
i is a neighborhood of the ith county. Thus the joint density of μ is
(9)
p
(
μ
∣
σ
μ
2
)
∝
∏
t
=
1
36
‍
p
(
μ
t
∣
σ
μ
2
)
.
Based on the likelihood function L
(
f
∣
ω
,
μ
), (7), and (9), the posterior distribution of ω is
(10)
p
(
ω
∣
f
,
μ
,
σ
μ
2
)
∝
L
(
f
∣
ω
,
μ
)
π
ϕ
(
ϕ
)
π
θ
(
θ
)
π
σ
(
σ
2
)
×
∏
t
=
1
36
‍
p
(
μ
t
∣
σ
μ
2
)
.
With specific priors of ϕ, θ
, and σ
2, 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 σ
2. Then the full conditional posterior distributions of ϕ, θ, and σ
2 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 M times and discarding the initial A burn-in iterations, then we have (
M
-
A
) iterations kept. Since the generated sample is not independent, we need to monitor the autocorrelations of the generated values and select a sampling lag L
>
1 after which the corresponding autocorrelation is low; that is, the length of the thinning interval is L. Considering the length of the thinning interval, the final number of iterations kept is M
′
=
(
M
-
A
)
/
L, and these independent samples will be used for posterior analysis. Assume the generated sample is
(11)
(
ω
(
1
)
,
ω
(
2
)
,
…
,
ω
(
M
′
)
)
.
Then for any function G
(
ω
) of the parameters of interest ω we can(1)
obtain a sample of the desired parametersG
(
ω
) by simply considering
(12)
(
G
(
ω
(
1
)
)
,
G
(
ω
(
2
)
)
,
…
,
G
(
ω
(
M
′
)
)
)
,
(2)
estimate the posterior mean ofG
(
ω
) by G
(
ω
)
¯
=
(
1
/
M
′
)
∑
s
=
1
M
′
G
(
ω
(
t
)
). And the posterior standard deviation is estimated by (
1
/
(
M
′
-
1
)
)
∑
s
=
1
N
′
[
G
(
ω
(
t
)
)
-
G
(
ω
)
¯
]
2. 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

Table1 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.Table 1
Demographic characteristics of HFMD cases and the pathogen types of some cases in Wenzhou, 2010–2012.
Year
2010
2011
2012
Total
Case number
32608
28980
42083
103671
Age (years old)
0-1
2810
2647
3916
9373
1-2
9267
6374
8904
24545
2-3
7830
8226
11340
27396
3-4
5807
6173
9327
21307
4-5
3180
2666
4227
10073
5–10
3207
2476
3828
9511
>10
507
418
541
1466
Gender
Male
20250
18212
26101
64563
Female
12358
10768
15982
39108
Sex ratio
1.638 : 1
1.691 : 1
1.633 : 1
Pathogen
CoxA16
41
79
122
242
EV71
648
574
537
1759
Other EV
121
67
212
400Most 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. TheP 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.Figure2 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].Figure 2
Monthly reported cases of HFMD in Wenzhou, 2010–2012.


### 3.2. Spatial Autocorrelation of HFMD Cases

Figure3 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 I
>
0.4, P values <
0.001).Table 2
The results of the spatial autocorrelation test on HFMD cases in Wenzhou, 2010–2012.
Year
Moran'sI
P values
2010
0.5182
<0.001
2011
0.4256
<0.001
2012
0.5813
<0.001Annual incidence rates (per 100,000 (18 years old or younger)) of HFMD per county in Wenzhou, 2010–2012.
(a)
(b)
(c)


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

For each county, we test the stationarity of{
S
t
j
,
t
=
1,2
,
…
,
36
} 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 f
t
j
=
S
t
j
-
S
(
t
-
1
)
j, where f
1
j
=
0. As shown in Figure 2, we can specify the seasonal period s
=
12. Then autocorrelation function and partial autocorrelation function are used to determine the values of p
,
q
,
P, and Q initially. That is, p
=
3
,
q
=
2
,
P
=
2, and Q
=
1. 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 p
=
0
,
q
=
0
,
P
=
1, and Q
=
0. Thus the model used in the data analysis is
(13)
f
t
j
=
α
j
f
(
t
-
12
)
j
+
e
t
+
μ
t
j
,
where e
t
~
N
(
0
,
σ
j
2
).Table 3
Adjusted Dickey-Fuller test of the monthly incidence rate of each county in Wenzhou, 2010–2012.
County
Statistics
P value
Dongtou
-
5.213
<0.01
Cangnan
-
4.106
0.017
Ruian
-
4.115
0.016
Pingyang
-
4.580
<0.01
Yongjia
-
3.627
0.045
Wencheng
-
3.500
0.059
Taishun
-
3.433
0.068
Yueqing
-
3.382
0.076
Ouhai
-
3.804
0.031
Longwan
-
3.888
0.025
Lucheng
-
3.852
0.028Table4 lists the estimates of α
j and σ
j 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.Table 4
Estimates of the parameters in the model (13).
County
α
j
σ
j
Dongtou
0.574
0.332
(0.133)
(0.066)
Cangnan
0.582
0.274
(0.138)
(0.058)
Ruian
0.655
0.254
(0.119)
(0.052)
Pingyang
0.447
0.322
(0.154)
(0.054)
Yongjia
0.559
0.306
(0.137)
(0.061)
Wencheng
0.491
0.652
(0.200)
(0.141)
Taishun
0.314
0.395
(0.224)
(0.054)
Yueqing
0.629
0.228
(0.118)
(0.044)
Ouhai
0.517
0.393
(0.145)
(0.075)
Longwan
0.502
0.404
(0.137)
(0.072)
Lucheng
0.785
0.156
(0.085)
(0.032)Comparison between data and the fitted values based on model (13).
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)
(k)


## 3.1. Prevalence of HFMD

Table1 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.Table 1
Demographic characteristics of HFMD cases and the pathogen types of some cases in Wenzhou, 2010–2012.
Year
2010
2011
2012
Total
Case number
32608
28980
42083
103671
Age (years old)
0-1
2810
2647
3916
9373
1-2
9267
6374
8904
24545
2-3
7830
8226
11340
27396
3-4
5807
6173
9327
21307
4-5
3180
2666
4227
10073
5–10
3207
2476
3828
9511
>10
507
418
541
1466
Gender
Male
20250
18212
26101
64563
Female
12358
10768
15982
39108
Sex ratio
1.638 : 1
1.691 : 1
1.633 : 1
Pathogen
CoxA16
41
79
122
242
EV71
648
574
537
1759
Other EV
121
67
212
400Most 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. TheP 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.Figure2 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].Figure 2
Monthly reported cases of HFMD in Wenzhou, 2010–2012.


## 3.2. Spatial Autocorrelation of HFMD Cases

Figure3 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 I
>
0.4, P values <
0.001).Table 2
The results of the spatial autocorrelation test on HFMD cases in Wenzhou, 2010–2012.
Year
Moran'sI
P values
2010
0.5182
<0.001
2011
0.4256
<0.001
2012
0.5813
<0.001Annual incidence rates (per 100,000 (18 years old or younger)) of HFMD per county in Wenzhou, 2010–2012.
(a)
(b)
(c)


## 3.3. Spatial-Temporal ARMA Model of HFMD

For each county, we test the stationarity of{
S
t
j
,
t
=
1,2
,
…
,
36
} 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 f
t
j
=
S
t
j
-
S
(
t
-
1
)
j, where f
1
j
=
0. As shown in Figure 2, we can specify the seasonal period s
=
12. Then autocorrelation function and partial autocorrelation function are used to determine the values of p
,
q
,
P, and Q initially. That is, p
=
3
,
q
=
2
,
P
=
2, and Q
=
1. 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 p
=
0
,
q
=
0
,
P
=
1, and Q
=
0. Thus the model used in the data analysis is
(13)
f
t
j
=
α
j
f
(
t
-
12
)
j
+
e
t
+
μ
t
j
,
where e
t
~
N
(
0
,
σ
j
2
).Table 3
Adjusted Dickey-Fuller test of the monthly incidence rate of each county in Wenzhou, 2010–2012.
County
Statistics
P value
Dongtou
-
5.213
<0.01
Cangnan
-
4.106
0.017
Ruian
-
4.115
0.016
Pingyang
-
4.580
<0.01
Yongjia
-
3.627
0.045
Wencheng
-
3.500
0.059
Taishun
-
3.433
0.068
Yueqing
-
3.382
0.076
Ouhai
-
3.804
0.031
Longwan
-
3.888
0.025
Lucheng
-
3.852
0.028Table4 lists the estimates of α
j and σ
j 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.Table 4
Estimates of the parameters in the model (13).
County
α
j
σ
j
Dongtou
0.574
0.332
(0.133)
(0.066)
Cangnan
0.582
0.274
(0.138)
(0.058)
Ruian
0.655
0.254
(0.119)
(0.052)
Pingyang
0.447
0.322
(0.154)
(0.054)
Yongjia
0.559
0.306
(0.137)
(0.061)
Wencheng
0.491
0.652
(0.200)
(0.141)
Taishun
0.314
0.395
(0.224)
(0.054)
Yueqing
0.629
0.228
(0.118)
(0.044)
Ouhai
0.517
0.393
(0.145)
(0.075)
Longwan
0.502
0.404
(0.137)
(0.072)
Lucheng
0.785
0.156
(0.085)
(0.032)Comparison between data and the fitted values based on model (13).
(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 Table1 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 MoransI with P value <
0.001. 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,
(14)
ϕ
j
(
B
)
Φ
j
(
B
)
f
t
j
=
θ
j
(
B
)
Θ
j
(
B
)
e
t
j
,
j
=
1,2
,
…
,
11
,
where {
e
t
j
,
j
=
1,2
,
…
,
11
} is intrinsic Gaussian autoregression distribution. Model (14) introduces the spatial effects by {
e
t
j
} 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.



---

*Source: 238724-2014-03-03.xml*