Computational and Mathematical Methods in MedicineVolume 2012 (2012), Article ID 758674, 6 pagesdoi:10.1155/2012/758674

Research Article

Comparing Statistical Models to Predict Dengue Fever Notifications

Arul Earnest,1,2 Say Beng Tan,¹ Annelies Wilder-Smith,3,4 and David Machin5,6

¹Centre for Quantitative Medicine, Office of Clinical Sciences, Duke-NUS Graduate Medical School Singapore, Singapore, 169857²Tan Tock Seng Hospital, Singapore, 308433³Institute of Public Health, University of Heidelberg, Germany⁴National University of Singapore, Singapore, 119077⁵University of Leicester, UK⁶University of Sheffield, UK

Received 16 September 2011; Revised 7 December 2011; Accepted 8 December 2011

Academic Editor: Chris Bauch

Copyright © 2012 Arul Earnest et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

Dengue fever (DF) is a serious public health problem in many parts of the world, and, in the absence of a vaccine, disease surveillance and mosquito vector eradication are important in controlling the spread of the disease. DF is primarily transmitted by the female Aedes aegypti mosquito. We compared two statistical models that can be used in the surveillance and forecast of notifiable infectious diseases, namely, the Autoregressive Integrated Moving Average (ARIMA) model and the Knorr-Held two-component (K-H) model. The Mean Absolute Percentage Error (MAPE) was used to compare models. We developed the models using used data on DF notifications in Singapore from January 2001 till December 2006 and then validated the models with data from January 2007 till June 2008. The K-H model resulted in a slightly lower MAPE value of 17.21 as compared to the ARIMA model. We conclude that the models' performances are similar, but we found that the K-H model was relatively more difficult to fit in terms of the specification of the prior parameters and the relatively longer time taken to run the models.

1. Introduction

The incidence of dengue fever (DF) has grown dramatically around the world in recent decades, with some 2.5 billion people now at risk of the disease [1]. Dengue haemorrhagic fever (DHF) is a potentially lethal complication, with an estimated 500 000 people requiring hospitalization each year, a very large proportion of whom are children. About 2.5% of those affected die [1].

DF is a viral vector-borne disease that is common in the tropics and subtropics and is primarily spread by the female Aedes aegypti mosquito. Mosquito vector control is important in restricting its spread. It has been found that controlling the vector population before disease is detected reducing transmission with a reduction of the Aedes aegypti population in a 3-month period, from 16% to 2%, as measured by the premises index [2]. However, predicting the incidence of vector-borne diseases like DF remains difficult, as DF shows strong variations over time [3–5]. In Singapore, seasonal trends are seen with peaks occurring generally in June or September. DF is characterized by both epidemic peaks that appear every 3–5 years, as well as seasonal oscillations within a year. Possible reasons for changes in outbreak patterns include change in number of infections due to interventions to eradicate the mosquitoes, as well as change in the number of people who are susceptible to the disease through prior infections [6]. Seasonal trends in DF can be caused by several factors, including climatic variables such as temperature and precipitation [7–10].

Autoregressive Integrated Moving Average (ARIMA) models have been used in applications such as the assessment of seasonal variation in selected medical conditions [11], and as a surveillance tool for outbreak detection [12]. ARIMA (AR, D, MA) models make use of previous observations to make predictions of future values using lag parameter values. Lags of the differenced series appearing in the forecasting equation are termed Auto Regressive (AR), those of the forecast errors, Moving Average (MA), and a time series that needs to be differenced to achieve stationarity, Differenced (D). The prediction process uses constantly updated information (in our example weekly DF cases) to predict the course of dengue in subsequent weeks.

Time series analysis of infectious diseases within the Bayesian framework has been considered in some studies [13–16]. One such example demonstrated that Klebsiella pneumoniae is related to the quantity of a third-generation antibiotic use (cephalosporin) in a hospital, with a lag of three months [17]. Others included a Knorr-Held (K-H) two-component model to incorporate both seasonal and epidemic characteristics of notifiable infectious diseases [15], as well as a Bayesian hierarchical time series model to detect outbreaks of Rubella and Salmonella infections [14].

Studies have compared ARIMA models with dynamic models for infectious diseases (fitted via maximum likelihood methods) [18, 19]. However, to the best of our knowledge, a direct comparison between the single-component (ARIMA) and two-component (K-H) models has not been undertaken.

2. Methods

The purpose of this paper is to compare the two-component K-H with the single-component ARIMA model in predicting weekly DF notifications. Different formulations of models within each type are compared, together with a sensitivity analysis of the K-H model, fitted within a Bayesian framework.

2.1. Data

The Singapore Infectious Diseases Act (1977) requires medical practitioners to notify all cases of DF to the Ministry of Health (MoH) within 24 hours. We obtained data from the published “Weekly Infectious Disease Bulletin”, available from the MoH website which uses the World Health Organization 2009 criteria for DF which is also detailed there [20]. All notified and registered DF cases were laboratory confirmed, with laboratory assays from Polymerase Chain Reaction (PCR) and/or NS1 antigen (in the first 5 days of illness) and/or a positive Dengue Immunoglobulin M after day 5 of illness.

We studied weekly DF notifications in Singapore till June 2008. Data from January 2001 to December 2006 was used to estimate the model parameters. Thereafter, we performed external validation of the models using data from January 2007 to June 2008.

2.2. ARIMA Model

If represents the number of cases of DF in week , then AR relates this observation to an earlier , where . MA relates the error (defined as the difference between the observed, , and that predicted, , notifications) at week to week (), where . D allows the differenced series, , to be modelled in the event of nonstationarity in the time series, where . Here , , and are the “orders” of the respective ARIMA components. Partial autocorrelation (PAC) and autocorrelation (AC) plots are used to determine and , respectively.

We describe the ARIMA (3,1,1) model equation used in our analysis. The number of cases of DF at week is denoted as , where is the first week for which DF is to be predicted where is the predicted number of DF cases for week , and , , and are the DF counts in the three immediate preceding weeks, termed lag 1, lag 2, and lag 3, respectively, and , are the error term at time and , respectively. In essence, we used observed values up till time to predict for dengue fever cases at week . is a constant and , , and are the coefficients for the three autoregressive terms in the model, is the first order moving average parameter and these are estimated within Stata V11.0 [21] via full or unconditional maximum likelihood estimates. For the ARIMA models, we used the Mean Absolute Percentage Error (MAPE) described below to compare predictive accuracy of the models.

2.3. Two-Component K-H Model

The K-H model distinguishes between the endemic, , and epidemic, , components of DF such that the number of cases observed and the corresponding prediction model is formulated as Here and have independent Poisson distributions with a composite parameters () and (), in which handles over dispersion, hence is also Poisson with parameter . This in turn corresponds to a negative binomial distribution with dispersion parameter . The mixing parameter, , is assumed to have a Gamma distribution with parameters () and ().

The endemic parameter, , is modelled as a harmonic wave (to handle strong seasonality inherent in infectious disease surveillance data) with see [15], where is the base frequency of the curve, which is suitable for weekly data and is a constant. The logarithmic transformation is necessary to ensure stationarity in the variance of the series.

The epidemic component is derived from the parameter sequence , which is assumed to be a piecewise constant [15] with unknown number of location and unknown location of the changepoints , that is, where are the unknown changepoints, such that for all . For , there is no changepoint and for all [15]. The piecewise function is needed to provide flexibility in the model in terms of modelling the outbreaks of dengue fever in addition to possible seasonal trends that we observe.

The two-component model formulation is completed by specifying prior distributions for the parameters in the model as follows: denotes a normal and Ga a Gamma distribution. was set to 10⁶, representing highly dispersed independent normal priors for each coefficient. is an identity matrix. For , , independent exponential distributions with mean and variance were specified. was then assigned a gamma hyperprior . The marginal prior distribution for is then a gamma-gamma distribution [22]. In our study, we used and , which corresponded to the gamma-gamma marginal of turning out to be an -distribution with degrees of freedom, which then indicates that the marginal prior probability of an outbreak occurring (i.e., ) is 0.39, while always favouring smaller values of , with the density function monotonically decreasing. The dispersion parameter for the negative binomial distribution, , which was designed to handle extra-Poisson variation in the data, was assigned a gamma hyperprior as well, with the following parameter, and were assigned values of 1 and 0.1, respectively in the original analysis corresponding to a prior mean and standard deviation of 10.

The K-H models were fitted using the customised Bayesian software Twins V1.0 [15]. Markov Chain Monte Carlo (MCMC) methods, in particular the Metropolis-Hastings algorithm, were used to estimate the parameters. For each model, we ran 200 iterations as burn-ins. These burn-in samples were discarded and not used in the analysis. We ran a further 60,000 iterations, but only saved every 20th observation, resulting in a final 3000 sample size. This was to circumvent the problem of autocorrelated samples.

2.4. Model Comparison

We compared the ARIMA model with the K-H model and as well conducted a sensitivity analysis on the K-H model using the MAPE: where is the total number of weeks of data.

The Bayesian analyses were based on several assumptions regarding the prior distributions, and we assessed the robustness of our results in a sensitivity analysis. For the sensitivity analyses, we considered 4 different scenarios which involved varying values of , or , while keeping the other variables at their original values: Model 1: and , Model 2: and , Model 3: and , and Model 4: and . The prior values for the sensitivity analysis were selected to represent a range of realistic scenarios where the probabilities of an outbreak were expected to be different. In particular, we selected priors where the probability of observing an outbreak ranged from 0.001 (for and ) to 0.5 ( and ).

3. Results

Figure 1 highlights the weekly distribution of DF notifications in Singapore from January 2001 to June 2008. It is evident that DF notification exhibits both seasonal trends (e.g., regular peaks around June or September and troughs seen in the first 4 months of the year) and epidemic trends (most markedly shown during the 2005 epidemic, when average weekly counts exceeded 600 cases).

Figure 1: Weekly cases of dengue fever (DF) in Singapore.

The autocorrelation plots for DF (Figure 2(a)) indicated that correlations gradually declined over the weeks to insignificant values after 12 weeks. The partial autocorrelations plots (Figure 2(b)) showed a spike at week 1 and week 4 indicating possible inclusion of AR terms of the order of up to four in the ARIMA model. We evaluated the various combinations, including autocorrelation terms 3 and 4 in our analysis.

(a) Correlogram

(b) Partial correlogram

Figure 2: Plots of autocorrelation and partial correlation for dengue fever (DF).

We explored various formulations of the ARIMA model, and we summarise some of the more important ones in Table 1. As can be seen, ARIMA (3,1,0) provided the lowest MAPE value of 19.86. Including a moving average term did not improve the fit of the model, as with adding an autocorrelation term of four. Adding a 12-month seasonal component (not shown) also did not lower the MAPE. The parameters for the final ARIMA model are shown in Table 2. We found all three autoregressive terms (), (), and () to be statistically significant. The parameters for the K-H model are also provided in Table 2.

Table 1: Comparison of MAPE values across various ARIMA models.

Model
Model specification
MAPE

1
ARIMA (1,0,0)
23.61

2
ARIMA (2,0,0)
23.09

3
ARIMA (3,0,0)
23.20

4
ARIMA (4,0,0)
23.23

5
ARIMA (3,1,0)
19.86

6
ARIMA (3,1,1)
19.96

Table 2: Parameters for the final models.

ARIMA model
Coefficient
95% confidence interval
P value

Constant ()
0.28
−3.86
4.41
0.896

AR 1 ()
−0.10
−0.16
−0.04
0.001

AR 2 ()
0.10
0.04
0.17
0.002

AR 3 ()
0.23
0.17
0.29
<0.001

K-H model
Coefficient
95% credible interval

25.1
18.4
32.3

3.3
1.9
3.6

−0.2
−0.3
0.3

−0.5
−0.7
−0.4

1.0
1.0
1.0

7.6
0.01
15.0

1.3
0.7
2.0

The comparison between the ARIMA and K-H model is shown in Figure 3. Table 3 shows the results from comparing the two models. Overall, the K-H model performed marginally better than the ARIMA model (MAPE of 17.21 and 17.54 resp.). In particular, the model predicted well (out-of-sample) for certain periods, including the early endemic periods between weeks 1 to 12. Fine-tuning the parameters for the K-H model allowed us to make better predictions for the epidemic periods, as we show in the sensitivity analysis (Table 4). For instance, the model predicted well for the epidemic periods within the weeks 17 to 24 (sensitivity analysis 4).

Table 3: Comparison of out-of-sample predictions (external validation) between ARIMA and K-H models.

MAPE
ARIMA
K-H

Overall
17.54
17.21

Stratified (in 4 week intervals)

Year 2007

    Weeks 1 to 4
17.07
14.27

    5 to 8
28.60
25.62

    9 to 12
33.41
30.63

    13 to 16
32.52
33.09

    17 to 20
21.83
20.53

    21 to 24
20.64
19.76

    25 to 28
12.86
13.22

    29 to 32
11.53
14.40

    33 to 36
8.54
10.26

    37 to 40
5.07
6.50

    41 to 44
18.49
17.42

    45 to 48
8.54
10.35

    49 to 52
15.70
12.44

Year 2008

    Weeks 1 to 4
11.13
11.16

    5 to 8
29.09
25.63

    9 to 12
16.39
19.41

    13 to 16
15.51
10.77

    17 to 20
19.21
18.38

    21 to 24
9.83
10.07

Table 4: Sensitivity analysis on K-H model parameters.

MAPE
Initial K-H model
Sensitivity analysis

1
2
3
4

Overall
17.21
17.71
17.71
17.50
16.54

Stratified (in 4 week intervals)

Year 2007

    Weeks 1 to 4
14.27
20.12
22.33
20.30
20.03

    5 to 8
25.62
25.41
25.66
25.12
23.39

    9 to 12
30.63
31.06
31.30
30.95
31.07

    13 to 16
33.09
33.20
32.31
32.49
27.89

    17 to 20
20.53
20.41
20.40
20.92
21.82

    21 to 24
19.76
21.14
20.90
21.25
21.58

    25 to 28
13.22
13.45
14.18
13.28
12.91

    29 to 32
14.40
13.98
13.04
13.39
10.65

    33 to 36
10.26
10.76
9.97
10.55
6.11

    37 to 40
6.50
6.69
6.39
5.54
3.30

    41 to 44
17.42
17.62
17.54
16.91
15.94

    45 to 48
10.35
11.30
10.59
11.09
10.67

    49 to 52
12.44
12.99
12.37
12.58
13.12

Year 2008

    Weeks 1 to 4
11.16
11.09
10.85
11.13
11.31

    5 to 8
25.63
25.53
26.01
25.49
25.83

    9 to 12
19.41
20.03
20.25
19.31
16.28

    13 to 16
10.77
10.72
10.47
10.72
11.25

    17 to 20
18.38
19.20
17.97
18.59
19.27

    21 to 24
10.07
9.33
10.98
9.95
9.92

A description of the parameters used in the sensitivity analysis is provided in the 4th page of the manuscript.

Figure 3: Comparison of out-of-sample forecasts of dengue fever (DF) between ARIMA and two-component K-H model (January 2007 to June 2008).

In terms of forecasting one-week ahead DF notifications, both methods performed well (Figure 3). For instance, the K-H model forecasted 53 (observed 58) for 2007 week 1, 356 (observed 371) for 2007 week 26, 112 (observed 115) for 2008 week 1, and 171 (observed 132) for 2008 week 26. It is worth noting that these results are for out-of-sample predictions.

The Bayesian analysis is influenced by the prior specification. As such, we investigated the robustness of our results to different formulation of the priors. These priors represented a wide range of realistic scenarios where the probability of an outbreak is expected to differ. As can be seen from Table 4, it appears that the models have generally similar MAPE values except for sensitivity analysis 4, where the MAPE is actually the lowest at 16.54. In our local setting, we found that specifying a small prior probability of 0.001 for an outbreak to occur provided a better fit of the data.

4. Discussion

We found that the K-H model performed better than the conventional ARIMA time series model; however, this was only marginal. Forecasting weekly cases of DF has immense implication for hospital resources planning. For an infectious disease ward, knowing the normal trend of DF, along with predictions of the following week’s DF can allow hospital planners to better plan for and allocate their manpower and other resources. Intensive media campaigns (e.g., television advertisements) in the weeks prior to a projected increase in DF notifications may prove to reduce the number of new cases.

Though we used the MAPE index to compare the models, other indices are also available. The Mean Squared Error, for instance, is calculated from the sum of the squared error values. Compared to MAPE, the values are not relative to the magnitude of the observation, and the values are not intuitively easy to interpret.

There were several limitations in our study. Firstly, our analysis was dependent on notifiable data. While clinicians are required to report all cases of DF and DHF to the MOH, there is a possibility that the cases could be underreported, especially since mild asymptomatic cases of DF may have not been diagnosed. While this may have led to an under-estimate in the forecasts, the comparisons across the models are still valid, as they make use of the same number of weekly cases.

In our analysis, we compared the predictive capability of the models using one-week ahead forecast of dengue fever notification. It is possible to forecast for periods longer than that, of course the predictions may inherently not be as accurate as a one-week forecast.

In conclusion, we found that both the final models chosen for the ARIMA and K-H models predict the future course of DF in Singapore reliably well, while the former performed marginally better. The ARIMA models were relatively faster to implement and run, while the K-H model was sensitive to the choice of priors, which needs to be carefully made before the study is conducted.

Acknowledgments

Funding for this study was received from the Duke-NUS SRP block grant as well as the Merit Award from the Yong Loo Lin School of Medicine, Singapore.

References

World Health Organisation, “Dengue and dengue haemorrhagic fever,” Revised March 117, World Health Organisation, Geneva, Switzerland, 2009.

E. E. Ooi, K. T. Goh, and D. J. Gubler, “Dengue prevention and 35 years of vector control in Singapore,” Emerging Infectious Diseases, vol. 12, no. 6, pp. 887–893, 2006.

G. Chowell, C. A. Torre, C. Munayco-Escate et al., “Spatial and temporal dynamics of dengue fever in Peru: 1994–2006,” Epidemiology and Infection, vol. 136, no. 12, pp. 1667–1677, 2008.

M. P. Mammen, C. Pimgate, C. J. M. Koenraadt et al., “Spatial and temporal clustering of dengue virus transmission in Thai villages,” PLoS Medicine, vol. 5, no. 11, article no. e205, pp. 1605–1616, 2008.

A. Mondini and F. Chiaravalloti-Neto, “Spatial correlation of incidence of dengue with socioeconomic, demographic and environmental variables in a Brazilian city,” Science of the Total Environment, vol. 393, no. 2-3, pp. 241–248, 2008.

B. K. W. Koh, C. N. Lee, Y. Kita et al., “The 2005 dengue epidemic in Singapore: epidemiology, prevention and control,” Annals of the Academy of Medicine Singapore, vol. 37, no. 7, pp. 538–545, 2008.

J. M. Brunkard, E. Cifuentes, and S. J. Rothenberg, “Assessing the roles of temperature, precipitation, and ENSO in dengue re-emergence on the Texas-Mexico border region,” Salud Publica de Mexico, vol. 50, no. 3, pp. 227–234, 2008.

D. D. Chadee, B. Shivnauth, S. C. Rawlins, and A. A. Chen, “Climate, mosquito indices and the epidemiology of dengue fever in Trinidad (2002-2004),” Annals of Tropical Medicine and Parasitology, vol. 101, no. 1, pp. 69–77, 2007.

C. A. Depradine and E. H. Lovell, “Climatological variables and the incidence of Dengue fever in Barbados,” International Journal of Environmental Health Research, vol. 14, no. 6, pp. 429–441, 2004.

A. Eearnest, S. B. Tan, and A. Wilder-Smith, “Meteorological factors and El Nino Southern Oscillation are independently associated with dengue infections,” Epidemiology and Infection, vol. 12, pp. 1–8, 2011.

R. Moineddin, R. E. G. Upshur, E. Crighton, and M. Mamdani, “Autoregression as a means of assessing the strenght of seasonality in a time series,” Population Health Metrics, vol. 1, article no. 10, 2003.

B. Y. Reis and K. D. Mandl, “Time series modeling for syndromic surveillance,” BMC Medical Informatics and Decision Making, vol. 3, no. 1, p. 2, 2003.

S. I. Hay, M. F. Myers, D. S. Burke et al., “Etiology of interepidemic periods of mosquito-borne disease,” Proceedings of the National Academy of Sciences of the United States of America, vol. 97, no. 16, pp. 9335–9339, 2000.

S. H. Heisterkamp, A. L. M. Dekkers, and J. C. M. Heijne, “Automated detection of infectious disease outbreaks: hierarchical time series models,” Statistics in Medicine, vol. 25, no. 24, pp. 4179–4196, 2006.

L. Held, M. Hofmann, M. Höhle, and V. Schmid, “A two-component model for counts of infectious diseases,” Biostatistics, vol. 7, no. 3, pp. 422–437, 2006.

P. D. O'Neill, “A tutorial introduction to Bayesian inference for stochastic epidemic models using Markov chain Monte Carlo methods,” Mathematical Biosciences, vol. 180, pp. 103–114, 2002.

J. L. Hay and A. N. Pettitt, “Bayesian analysis of a time series of counts with covariates: an application to the control of an infectious disease,” Biostatistics, vol. 2, no. 4, pp. 433–444, 2001.

D. He, E. L. Ionides, and A. A. King, “Plug-and-play inference for disease dynamics: measles in large and small populations as a case study,” Journal of the Royal Society Interface, vol. 7, no. 43, pp. 271–283, 2010.

A. A. King, E. L. Ionides, M. Pascual, and M. J. Bouma, “Inapparent infections and cholera dynamics,” Nature, vol. 454, no. 7206, pp. 877–880, 2008.

Infectious Diseases Guidelines and Ministry of Health Singapore, http://www.moh.gov.sg/content/moh_web/home/Publications/guidelines/infectious_diseases_guidelines/2011/a_guide_on_infectious_diseases_of_public_health_importance_in_singapore.html, 2011

Stata V11.0. Stata Corp College Station, Texas, USA

J. M. Bernado and A. F. M. Smith, Bayesian Theory, Chapman & Hall, London, UK, 1994.