1 Introduction

Hand, foot, and mouth disease (HFMD) is a viral illness common in infants and children younger than 5 years old. Since the first HFMD case was reported in New Zealand in 1957, many parts of the world have reported epidemics of HFMD. The Asia Pacific Region is one of the major HFMD-endemic regions. Most countries there have reported HFMD epidemics in the past 30 years, including Australia (Burry et al. 1968), Japan (Daisuke and Masahiro 2011), South Korea (Ryu et al. 2010), Malaysia (Chua and Kasri 2011), Singapore (Ang et al. 2009), Thailand (Puenpa et al. 2011), Vietnam (Tu et al. 2007), and China (Xing et al. 2014a, b). Although HFMD is a mild viral infection, children under 5 years of age are particularly susceptible to the neurologic complications of this infection, which can result in encephalitis and death. During 2008–2013, 9,048,244 probable cases of HFMD were reported to the Chinese Center for Disease Control and Prevention (China CDC), of which 355,563 (3.9%) were laboratory confirmed, 92,503 (1.02%) were severe, and 2,707 (0.03%) were fatal (Xing et al. 2014a, b). Pathogen detection results have shown that EV71 causes approximately 50% of cases of HFMD in China and that it is a significant cause of morbidity and mortality (Pallansch and Roos 2007; Lin et al. 2002). Currently, there is no available vaccine against HFMD, even though a newly developed formalin-inactivated EV71 (FI-EV71) vaccine has been tested in clinical trials and has shown efficacy against EV71 (Tsou et al. 2015). Therefore, the early detection of the aberration of infectious disease occurrence and the estimated disease outbreak risk is significant for slowing down the epidemic spread and reducing the morbidity and death caused by the disease.

Some methods for forecasting disease outbreaks of HFMD have been successfully applied in some countries and regions. In Japan, an index calculated from the number of cases per week per sentinel medical institution was used as a disease warning threshold. If the index in an area exceeds the critical value for the onset of an epidemic and continues until the index in that area is lower than the critical value for the end of an epidemic, it can predict that there is an unusual aberration of HFMD cases (Murakami et al. 2004). In China, a nationwide web-based automated system for outbreak detection and rapid response was developed in 2008. In the China Infectious Disease Automated-alert and Response System (CIDARS), the fixed-threshold detection method, the temporal detection method, and the spatial detection method were applied to detect three disease aberrations (Yang et al. 2011). HFMD surveillance has been included in CIDARS. In some regions of China with a high HFMD prevalence, researchers have proposed many statistical models to estimate local development trends of HFMD. These research results give support for local early warnings of HFMD. In Shandong province, a conditional logistic regression was used to explore an early warning index of HFMD death based on clinical treatment, past medical history, clinical symptoms, and signs of ill children (Liu et al. 2013). A hybrid model combining a seasonal autoregressive integrated moving average (ARIMA) model and a nonlinear autoregressive neural network (NARNN) was proposed to predict the expected incidence of cases in Shenzhen city from December 2012 to May 2013 (Yu et al. 2014). In Dalian city, a negative binomial regression model was used to estimate the weekly baseline number of HFMD cases and identified the optimal alerting threshold between different tested threshold values during the epidemic and non-epidemic year (An et al. 2016). Most of these early warning models use historical case data to forecast an imminent outbreak; in fact, these models use only retrospective research (Edmond et al. 2011). The operation of these models also closely depends on the data integrity of the risk factor variables. However, in many cases, the risk factor variables are difficult to collect completely in time.

A Bayesian belief network (BBN) is a graphical model for probabilistic relationships among a set of variables; it can express the mutual dependencies between variables in terms of quality and quantity, and it provides a good explanation for knowledge representation and reasoning (Holický et al. 2013). This method can also produce uncertainty estimates even with missing values. Because of its simplicity and soundness, BBN has been used in many fields, especially in medical and industrial diagnosis (Liao et al. 2010; Liu et al. 2012). So this study establishes an early warning model based on BBN to predict the outbreak risk of severe HFMD and death. Hunan province, one of the regions with the highest prevalence of HFMD in China, was selected as the study area. According to previous study, children under 5 years old accounted for most of the HFMD cases (Deng et al. 2013), and EV71 predominated in laboratory-confirmed cases, accounting for 45% of mild, 80% of severe, and 93% of fatal cases (Pallansch and Roos 2007; Lin et al. 2002). In order to analysis the potential seasonal drivers of HFMD, social-economic factors such as gross domestic product and gross regional product were also put into account (Xing et al. 2014a, b). Meanwhile, a variety of studies identified that meteorological factors including temperature, rainfall, relative humidity, air pressure, and hours of sunshine all have a strong association of the prevalence of HFMD (Deng et al. 2013; Xing et al. 2014a, b; Yu et al. 2014). Therefore, in the proposed model, the monthly pathogen detection ratio, demographic factors, socioeconomic factors, and monthly meteorological factors are all taken into account to estimate the outbreak risk of death and severe HFMD in the next month. Based on the characteristic of BBN, once the structure is established and parameters are learned, the risk of death and sever HFMD can be estimated based on any one node (factor) in the BBN, which can handle the situation where the data of some variables are unavailable.

2 Data and method

2.1 Study area

In this study, Hunan, a province in south central China, was selected as the study area. Hunan is located south of the middle course of the Yangtze River and south of Lake Dongting, situated between the 108°47′–114°16′ east longitudes and the 24°38′–30°08′ north latitudes. It covers an area of 211,800 square kilometers and is divided into 14 cities. Mountains and hills occupy more than 80% of the area, with plains comprising less than 20% of the whole province. Hunan’s climate is subtropical, where sunshine and rainfall are both abundant concentrated. In spring, the temperature would fluctuate in a wide range, whereas in summer and autumn, it is pretty hot and humid, and in winter it’s usually cool and damp. According to the meteorological data, the annual average temperature of Hunan is around 16–19 degree centigrade, with average 38 °C (37–46 °F) in January and average 27–30 °C (81–86 °F) in July. The average annual precipitation is 1200–1700 mm and has an uneven distribution (information from Wikipedia article, “Hunan”) (Fig. 1).

Fig. 1
figure 1

Location of Hunan province in China

Full size image

In terms of the socioeconomic condition, according to the statistical yearbook of Hunan province, during the study period, there were around 66,380,000 people living in Hunan with a GDP of 32903.71 RMB per capita, and the people living in cities accounted for 46.65% of the total population. Since its geographical environment advantages and rich natural resources, the development of Hunan’s economy is good.

During 2009–2013, 547,027 HFMD cases were reported in Hunan. The 5-year prevalence reached 167.83 per 100,000 persons. There were 5,638 severe cases (1.03% of all cases). There were 296 cases of death, a mortality rate of 0.05%. The annual HFMD prevalence varied significantly over these 5 years, ranging from 287.32 (in 2012) to 54.31 per 100,000 persons (in 2009). According to the results of previous studies, there was also a significant difference in HFMD prevalence among the 14 cities in Hunan. The HFMD prevalence in Hunan also showed a seasonal trend in 2009–2013; April to July and September to November are the peak months of HFMD in Hunan (Yang et al. 2015). Hence, a clustering method for panel data was used to classify experimental scenarios based on (1) the monthly rate, severity, and mortality of HFMD in the all cities during the study period; (2) the monthly increase in the rate, severity, and mortality in the all cities during the study period; and (3) the rate of change of the rate, severity, and mortality increment in the all cities during the study period. Four experimental scenarios included a risk assessment of death and severe HFMD in the next month (1) in high-incidence cities during a peak period, (2) in high-incidence cities during a non-peak period, (3) in low-incidence cities during a peak period, and (4) in low-incidence cities during a non-peak period (Fig. 2).

Fig. 2
figure 2

High-incidence and low-incidence cities in Hunan (Xiangxi is not included due to low data quality)

Full size image

2.2 Data sources

The period of the study was from January 1, 2011, to December 31, 2013, and the study unit was a city. The study included the rate, severity, and mortality of HFMD and the rate of pathogen EV71 detection in mild cases. People older than 5 years of age were ignored in the study, because HFMD mainly infects children under the age of five. In addition, considering that prekindergarten children and kindergarten students have different routes of exposure to pathogen EV71, we separately calculated the rates of pathogen EV71 detection in mild cases for children 0–3 years old and children 3–5 years old. All HFMD cases and pathogen detection laboratory data were collected from the local CDC.

The environmental factors in this study were meteorological and socioeconomic. In accordance with previous studies (Yang et al. 2015; Gao et al. 2014), the meteorological data mainly included the average monthly air pressure, maximum monthly wind speed, average monthly wind speed, average monthly relative humidity, minimum monthly relative humidity, and average monthly temperature, monthly temperature difference, monthly rainfall, and monthly sunshine. The meteorological data of all cities were obtained from the China Meteorological Data Sharing Service System and were interpolated by the inverse distance weighted interpolation method, based on the data of the 756 stations nationwide. The socioeconomic factors included the density of the susceptible population (i.e., children under the age of five) and urbanization levels. These data were extracted from the Hunan Statistical Yearbook. ArcGIS 10.0 was used to store the geographic data and display the results of this study.

2.3 Methodology

The framework for establishing an early warning model based on BBN to predict the outbreak risk of severe HFMD and death in Hunan comprises three parts: (1) defining the risk level; (2) building the BBN construct; and (3) calculating the joint probability of the target variable in the BBN and assessing the outbreak risk. In this study, the BBN software tool used was BN PowerSoft (Cheng and Greiner 1999).

2.3.1 Classification of risk level

In this study, we introduced the dynamic control chart percentile method (Yang et al. 2014; Xiong et al. 2013) to determine the risk level of HFMD outbreak; this model had been successfully used in CIDARS. In the dynamic control chart percentile method, to account for the stability of the data, the most recent one-month period is used as the current observation period. The corresponding historical period included the same one-month period, the preceding one-month period, and the following 1-month period over the preceding 3 years of historical data. The current observation period and historical data block were dynamically moved forward month by month.

The reported death and severe incidence in the nine blocks of historical data were first arranged in ascending order, \( X_{1} ,X_{2} , \ldots ,X_{d} , \ldots ,X_{n} \left( {n = 9} \right) \). Then the percentile of the nine blocks of historical data was set as the indicator of potential aberration; \( {\text{p}} \) is a specified percentile of historical data, and \( X_{P} \) can be calculated as follows:

$$ {\text{d}} = 1 + (n - 1)p $$
(1)
$$ X_{P} = X\left( {\left[ d \right]} \right) + \left( {X\left( {\left[ {d + 1} \right]} \right) - X\left( {\left[ d \right]} \right)} \right)\left( {d - \left[ d \right]} \right) $$
(2)

([] means take the integer of the number.)

Finally, the method defined the outbreak risk level of severe HFMD and death by comparing the reported death and severe incidence in the current 1-month period to different control limits of the corresponding historical period at the city level. That is, the risk was regarded as high and the risk level could be set at 2 if the actual incidence exceeded the upper control limit (\( P_{80} \)). If the actual incidence was greater than \( P_{50} \) and less than \( P_{80} \), the risk was moderate, and the risk level was set as 1. If the actual incidence was less than \( P_{50} \), the risk level was set at 0. The selection of these threshold values takes the ones used in CIDARS as the standard.

2.3.2 Building the BBN construct

A Bayesian network is represented by \( {\text{BN}} = \left\langle {N,A,\varTheta } \right\rangle \), where \( \left\langle {{\text{N}},{\text{A}}} \right\rangle \)is a directed acyclic graph (DAG)—each node \( {\text{n}} \in {\text{N}} \) represents a domain variable (corresponding perhaps to a database attribute), and each arc \( {\text{a}} \in {\text{A}} \) between nodes represents a probabilistic dependency between the associated nodes. Associated with each node \( n_{i} \in N \) is a conditional probability distribution (CPtable), collectively represented by \( \varTheta = \left\{ {\theta_{i} } \right\} \), which quantifies how much a node depends on its parents (see Pearl 1988). In this study, the BBN is presented composed of a set of nodes, representing the risk variables of severe HFMD and death, and a set of edges, representing the probabilistic causal dependence among the variables. Before inputting the variables into the BBN, continuous variables were discretized by the entropy method. In the entropy method, the minimum entropy value in a numerical attribute is chosen as a split point, and then the classified results of the intervals recursively and eventually gain discretized results.

In the BBN, the nodes with edges directing into them are called ‘‘child’’ nodes, and the nodes from which the edges depart are called ‘‘parent’’ nodes. That is, if there is an edge from node N1 to another node N2, N1 is called the parent of N2. Nodes without arches directing into them are called ‘‘root’’ nodes. Therefore, a proper BBN network construct is very significant for assessing a future increase in the outbreak risk of severe HFMD and death, based on the BBN.

Generally, three methods—the expert knowledge method, the network structure learning algorithm, and an integration of the above two methods—are applied to build the Bayesian belief network. The expert knowledge method builds a network according to experts’ comments and the results of literature reviews, while the network structure learning algorithm is based on the independence test proposed by Cheng and Greiner (1999). The network structure learning algorithm has three phases: drafting, thickening, and thinning (Liao et al. 2017). In the first phase, the algorithm computes the mutual information of each pair of nodes as a measure of closeness, and it creates a draft based on this information. In the second phase, the algorithm adds arcs when the pairs of nodes cannot be d-separated. The node X is said to be d-separated from node Y if there is no active undirected path between X and Y. The result of the second phase is an independence map (I-map) of the underlying dependency model. In the third phase, if an open path separate from the current arc remains, it will be brought from the set of directed arcs joining vertices. Then the independence of the corresponding block set that blocks each open path between these two nodes by a set of minimum number of nodes is tested. The arc will be removed from the I-map if the two nodes of the arc can be d-separated. The result of the third phase is the minimal I-map. In the study, we applied an integration of the expert knowledge method and the network structure learning algorithm to construct the network.

2.3.3 Assessing the outbreak risk of severe HFMD and death

After the structure of Bayesian network was constructed (i.e., determining what depends on what), the next step is learning the parameters (i.e., the strength of these dependencies, as encoded by the entries in CPtables). For each BBN node, a conditional probability formula or table, which represents the probabilities of each value of the node given the conditions of its parents, is supplied. Meanwhile, the distributions of the parents can be calculated given the values of their children (Cooper and Herskovits 1992). After the conditional probability table is obtained, the joint probability distributions can be decomposed into local probability distributions (Liu et al. 2012). Therefore, in the study, the assessment of the outbreak risk of severe HFMD and death was performed simply by calculating the conditional probabilities among the variables, based on their cause and consequence relationships. The conditional probability distribution in the BBN can be expressed as follows:

$$ {\text{P}}\left( {X_{1} = x_{1} , \ldots ,X_{n} = x_{n} } \right) = {\text{P}}\left( {x_{1} , \ldots ,x_{n} } \right) = \mathop \prod \limits_{i = 1}^{n} P\left( {x_{i} |Parents\left( {X_{i} } \right)} \right) $$
(3)

where \( {\text{P}}\left( {x_{1} , \ldots ,x_{n} } \right) \) refers to the probability of a specific combination of values \( x_{1} , \ldots ,x_{n} \) from the set of variables \( X_{1} , \ldots ,X_{n} \), where Parent (\( X_{i} \)) refers to the set of \( X_{i}^{ '} s \) immediate parent nodes. Thus, \( P\left( {x_{i} |Parents\left( {X_{i} } \right)} \right) \) reflects the conditional probability, which is related to the node \( X_{i} \) based on its parent nodes (Liao et al. 2010).

3 Results

3.1 BBN structure for assessing the outbreak risk of severe HFMD and death

As described in Sect. 2.1, four BBN structures were applied to assess the risk of death and severe HFMD in the next month in different experimental scenarios. Before entering the data into the models, all the data used in this research should be discretized. Table 1 shows the results of each variable classification attribute in the BBN. The experimental scenarios included (1) in high-incidence cities during a peak period, (2) in high-incidence cities during a non-peak period, (3) in low-incidence cities during a peak period, and (4) in low-incidence cities during a non-peak period. In each experimental scenario, 70% of the surveillance data was used to build the early warning model for the outbreak of severe HFMD and death in the next month. Thirty percent of the surveillance data was applied to evaluate the accuracy of the model.

Table 1 The discretization results of each classification attribute in the BBN
Full size table

Figure 3 illustrates the network structures of the BBN used to estimate the outbreak risk of severe HFMD and death in the next month in the four experimental scenarios. In each network structure, we selected the variable ‘‘Risklevel’’ as the parent node of the various factors. All independent variables were represented by nodes; the arcs between the nodes indicated the relationships of the variables. As seen in Fig. 3, the effect on the outbreak risk of severe HFMD and death in the next month varied among the selected factors. However, in all the experimental scenarios, there were some common risk factors for this outbreak risk. The common risk factors included the total current monthly rates of pathogen EV71 detection in mild cases (the variable “EV71_rate”), the current monthly rates of pathogen EV71 detection for children aged 0–3 years (the variable “EV71_rate_1”) and aged 3–5 years (also the variable “EV71_rate_1”). In addition to the influence of the spread of pathogen EV71, some meteorological variables, including the current maximum monthly wind speed (the variable “MaxWinSpeed”), the average monthly temperature (the variable “AveTemp”), and the average monthly relative humidity (the variable “AveReHumidy”), affected the outbreak risk of severe HFMD and death in the next month in the high-incidence cities. During the peak period, the outbreak risk of severe HFMD and death in the next month in all cities was also related to the current monthly temperature and the maximum monthly wind speed. Nevertheless, the meteorological factors that had a direct impact on the outbreak risk of severe HFMD and death in the next month during the non-peak time were changed into the current monthly rainfall (the variable “Rainfall”) and the average monthly relative humidity. The influence of all socioeconomic factors on the outbreak risk of severe HFMD and death in the next month also varied in the different experimental scenarios.

Fig. 3
figure 3

BBN structures a in high-incidence cities during a peak period, b in high-incidence cities during a non-peak period, c in low-incidence cities during a peak period, and d in low-incidence cities during a non-peak period

Full size image

In our BBN structures, the total current monthly rates of pathogen EV71 detection in mild cases was one of risk factors with the most direct influence on the outbreak risk of severe HFMD and death in the next month. Different from how they influenced this outbreak risk, socioeconomic factors mainly affected the total current monthly rates of pathogen EV71 detection in mild cases in all the experimental scenarios. However, meteorological factors also had an impact on the spread of pathogen EV71 in some experimental scenarios. For example, the total current monthly rates of pathogen EV71 detection in mild cases in the high-incidence cities were directly related to the current monthly rainfall. During the peak period, the current maximum monthly wind speed also influenced the total current monthly rates of pathogen EV71 detection in mild cases.

3.2 Outbreak risk assessment of severe HFMD and death in Hunan

As we constructed four BBN network structures for these four study areas, a probabilistic inference reflecting the dependence relationships between the risk level of severe HFMD and death and the relevant risk factors was essential. Table 2 pinpoints the conditional probability distribution of the four BBNs for estimating the HFMD outbreak risk in the four study areas. It is obvious that there are two dominant aspects.

Table 2 Conditional probability distribution table of the BBN for estimating HFMD outbreak risk in Hunan different areas (Density of the susceptible population and Minimum Relative Humidity has no relationship with the outbreak risk assessment of HFMD severe and death in all areas, so they have been ignored in this table.)
Full size table

First, the primary variables that impact the risk levels differ to a great extent in these four different areas, and various probabilities are inferred by the Bayesian posterior probability estimation. In other words, there are various probabilities of the different values of the variables in estimating the risk levels (i.e., the three risk levels 0, 1, 2). In high-incidence cities during the peak period, as in the BBN structure described above, the variables that most notably produce an effect on the risk level of severe HFMD and death are the rate of pathogen EV71 detection in mild cases for children 0–5 and 3–5 years old, the urbanization level, the average wind speed, the maximum wind speed, the average relative humidity, the temperature difference, the average temperature, the average air pressure, and sunshine hours. The probabilistic inference results of all these relative variables linked in the above structure vary because of their different influences on the HFMD risk level. For instance, when the value of the rate of pathogen EV71 detection in mild cases for children aged 3–5 is in class one (< 0.05), the risk levels for 0, 1, and 2 are 0.422, 0.305, and 0.388, respectively, so we can conclude that the probability of no risk in this area is greater. However, when the values of the rate of pathogen EV71 detection for children 3–5 years old increase, the probabilistic inference results of all risk level are reduced, which means that the higher the rate of pathogen EV71 detection for children 3–5 years old, the lower the risk forecast for any risk level.

A similar situation occurs with other influencing variables, such as urbanization level, while the average temperature variable is just the opposite. When the value of this variable is in class 5 (> 27.321), the probability experiences a peak at the risk level of 0, with 0.462, which means that the areas with these temperatures are safe. In high-incidence cities during a non-peak period, the variables that most notably produce an effect on the risk level of severe HFMD and death are the rate of pathogen EV71 detection in mild cases for children aged 0–3 and 3–5, the urbanization level, the average relative humidity, and rainfall. Similarly, the probabilistic inference results of these relative variables may vary. If the rainfall is in class 4 (3.157–5.044), the probability will be 0.545 at a risk level of 2, so these places should be aware of the corresponding protective measures taken by CDC or governments. In low-incidence cities during a peak period, the variables are the rate of pathogen EV71 detection in mild cases for children 0–5 and 3–5 years old, the average temperature, and the average air pressure. In low-incidence cities during a non-peak period, the variables are the rate of pathogen EV71 detection in mild cases for children aged 0–5, the rate of pathogen EV71 detection in mild cases for children aged 3–5 years, the urbanization level, the average wind speed, the average relative humidity, the temperature difference, the average air pressure, and rainfall.

Secondly, the same variable can have a different impact on different areas. In Table 2, we can see that the rate of pathogen EV71 detection in mild cases for children aged 0–5 has an impact in all four study areas, while other variables impact just some areas. For example, rainfall has effects on both high-incidence and low-incidence cities during a non-peak period. In addition, in all areas, the density of the susceptible population and the minimum relative humidity have no relationship with the outbreak risk assessment of severe HFMD and death, so they have been ignored in Table 2. Otherwise, the probabilities of the same variable will vary in different areas.

Take a factor, Average Relative Humidity (ARH), as an example: For high-incidence cities during peak period area, when the value of ARH falls in “70.850–73.994”, it belongs to the “class 2” according to the discretization results in Table 1. Then, we can utilize Table 2, which shows that the probability of risk level 0 is 0.333, the probability of risk level 1 is 0.176, and the probability of risk level 2 is 0.175, to predict the most probable risk level of HFMD outbreak in high-incidence cities during peak is 0. Meanwhile, for high-incidence cities during non-peak period area, when the ARH value is in class 2 (47.044–53.535), the probability of risk level 2, 0.455, is the highest. Other variables experience the same situation.

3.3 Accuracy analysis

The model performance was evaluated with receiver operating characteristic (ROC) analysis, which has been widely used in evaluating the predictive accuracy of a binary classifier. A typical ROC curve plots the true positive rate (sensitivity) against the false positive rate (1-specificity) for the entire range of possible thresholds, thus providing a unified representation for assessing the overall model performance. The area under the curve (AUC) can be used as a single performance measure to decide whether the model prediction is better than random (0.5). A perfect model will yield an AUC value of 1. Training AUC values is fairly high across models and is higher than the testing AUC values, as anticipated. The testing AUC values, which demonstrate the model’s actual predictive powers, suggest whether the model predictions are better than random.

Because the ROC curve applies only to binary classifiers, we need to reclassify the HFMD outbreak risk level (0 or 1 or 2) into two classes. Situation 1: we define level = 0 as the class that means there is no HFMD outbreak, and we merge level = 1 and level = 2 as the other class to represent the outbreak. Situation 2: we merge level = 0 and level = 1 to signify that the epidemic is not serious, and level = 2 is a representation that the HFMD outbreak risk is very high.

As can be seen from Fig. 4, which clearly illustrates the performance of the BBN models applied in the four experimental scenarios, in terms of 10 times the average testing AUC value, no matter how the risk levels were reclassified, the models achieved the best predictive accuracy in the first scenario (high-incidence cities during a peak period), where the average testing AUC value was 0.9636 for reclassification 1 and 0.9139 for reclassification 2. This was followed by the third scenario (low-incidence cities during a peak period) and the fourth scenario (low-incidence cities during a non-peak period). The second scenario (high-incidence cities during a non-peak period) produced a comparatively low AUC value of 0.7285, which was also actually a good performance value. In terms of the standard deviation, which illustrates the stability of the predictive precision, the same order can be found, in which the most stable predictive accuracy applied to the first scenario, followed by the third and fourth scenarios; the second scenario produced the highest standard deviation value, and its predictive accuracy was the most unstable.

Fig. 4
figure 4

ROC curve for belief network models applied in four experimental scenarios

Full size image

A conclusion that can be drawn from the analysis above is that whether the situation involves a peak or non-peak period has a crucial influence on the performance of the BBN Model. In a peak period, the model’s AUC value is high and the SD value is low, which means an outstanding predictive performance. Otherwise, in a non-peak period, the model achieves a comparatively low AUC value and a high SD value, which means that the performance of the model is not as good, and the predictive result is very unstable.

Meanwhile, as a comparison with other predictive models, we also applied the logistic regression model, which is a strong scientific method utilizing one or more variables to forecast another dependent variable value. We also applied the rough set method, proposed by Pawlak in 1982, which has been used successfully in diagnosis and outcome prediction, for the same data of high-incidence cities during a peak period. In order to verify the accuracy and reliability of the different models, we process bootstrap resampling for 30 times, producing 30 sets of experimental data, and input them to each model and evaluated the predictive results through the ROC curve. The ROC curves for each model’s average prediction of 30 replicate runs, with a mean test AUC value and corresponding standard deviation, are provided in Fig. 5.

Fig. 5
figure 5

ROC curve for Bayesian belief network model, logistic regression model, and rough set method prediction averaged on each of 30 replicate runs. (The 1:1 line indicates the condition if the prediction is completely out of random chance. (AUC = 0.5))

Full size image

In the lateral view, no matter how the risk level was reclassified, the Bayesian network model achieved a very high performance evaluation (the average AUC test was 0.9628 for situation 1 and was 0.9234 for situation 2) and was better than the logistic regression model (where the average AUC test was 0.8679 for situation 1 and 0.8342 for situation 2) and the rough set method (where the average AUC test was 0.606 for situation 1 and 0.5841 for situation 2). In the longitudinal view, the stability of the Bayesian network was good (the standard deviation was 0.0382 for situation 1 and was 0.0656 for situation 2), regardless of the classification of the disease, and the logistic regression model (where the standard deviation was 0.1166 for situation 1 and was 0.1019 for situation 2) and the rough set method (where the standard deviation was 0.1079 for situation 1 and was 0.0879 for situation 2) did not have this advantage.

4 Discussion and conclusions

In this study, an early warning model was developed to predict the outbreak risk of severe HFMD and death in Hunan Province, China. In accordance with the temporal spatial epidemic trend of HFMD in Hunan between January 2010 and December 2013, the study was divided into four experimental scenarios. In each scenario, the proposed model applied probabilistic inference based on a BBN structure to assess the outbreak risk of severe HFMD and death in the next month. In four BBN structures, the relationships among outbreak risk, the pathogen detection rate, demographics, and socioeconomic and meteorological factors varied. ROC analysis was applied to evaluate the model’s performance. The results showed that the performance of the proposed method was better than those of the rough set algorithm and the logistic regression.

Different from previous early warning models for HFMD, the proposed model used a network to express the mutual dependencies between the outbreak risk of severe HFMD and death and the risk factors. This clear representation of variable relationships will help us to understand and explore the causes of outbreaks of severe HFMD and death. More importantly, this model can make uncertainty estimates even with missing values. Although the proposed method was tested only in Hunan, it can be used as an early warning tool to solve HFMD outbreak prevention and control problems in any region, since BBN can provide a simple and unified method to address different problems in different areas.

However, there are also some limitations. The BBN in the proposed model can only deal with categorical data, so all continuous data needs to be discretized before being input to the model. This data discretization may lead to information loss. Moreover, while an accurate network is significant to the accuracy of the model, it is difficult to choose suitable network construction algorithms. Because the transmission path of HFMD in different areas varied, it was also hard to include all the risk factors in the study. Therefore, in a later study, we should take more risk factors of HFMD outbreaks into account.