# Comparing Statistical Models to Predict Dengue Fever Notifications

Authors: Arul Earnest; Say Beng Tan; Annelies Wilder-Smith; David Machin
Journal: Computational and Mathematical Methods in Medicine
 (2012)

---

## 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 femaleAedes 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.


---

## Body

## 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 femaleAedes 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 aegyptipopulation 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

IffT represents the number of cases of DF in week T, then AR relates this observation to an earlier fT-J, where J=1,2,…,T-1. MA relates the error (defined as the difference between the observed, f, and that predicted, F, notifications) at week T to week (T-K), where K=1,2,…. D allows the differenced series, ΔT=(fT-fT-L), to be modelled in the event of nonstationarity in the time series, where L=0,1,2,…. Here J, K, and L are the “orders” of the respective ARIMA components. Partial autocorrelation (PAC) and autocorrelation (AC) plots are used to determine J and K, respectively.We describe the ARIMA (3,1,1) model equation used in our analysis. The number of cases of DF at weekT is denoted as fT, where T is the first week for which DF is to be predicted(1)FT=μ+fT-1+φ1fT-1+φ2fT-2+φ3fT-3+θεT-1+εT,
where FT is the predicted number of DF cases for week T, and fT-1, fT-2, and fT-3 are the DF counts in the three immediate preceding weeks, termed lag 1, lag 2, and lag 3, respectively, and εT, εT-1 are the error term at time T and T-1, respectively. In essence, we used observed values up till time T-1 to predict for dengue fever cases at week T. μ is a constant and φ1, φ2, and φ3 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,x, and epidemic, y, components of DF such that the number of cases observed fT=xT+yT and the corresponding prediction model is formulated as(2)FT=XT+YT.
Here XT and YT have independent Poisson distributions with a composite parameters (ωTνT) and (ωTλTFT-1), in which ωT handles over dispersion, hence FT is also Poisson with parameter ωT[νT+λTFT-1]. This in turn corresponds to a negative binomial distribution with dispersion parameter ψ. The mixing parameter, ωT, is assumed to have a Gamma distribution with parameters (ψ+FT) and (ψ+νT+λTFT-1).The endemic parameter,νT, is modelled as a harmonic wave (to handle strong seasonality inherent in infectious disease surveillance data) with(3)log⁡⁡νT=γ0+γ1sin⁡⁡(2πT52)+γ2cos⁡⁡(2πT52),
see [15], where 2π/52 is the base frequency of the curve, which is suitable for weekly data and γ0 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λ=(λ1,…,λn), which is assumed to be a piecewise constant [15] with unknown number of location K and unknown location of the changepoints θ1<⋯<θK, that is,
(4)λT={λ(1),ifT=1,2,…,θ1,λ(k),ifT=θk-1+1,…,θk,λ(K+1),ifT=θK+1,…,n,
where θ1<θ2<⋯<θK are the K unknown changepoints, such that θk∈{1,2,…,n-1) for all k∈(1,2,…,K). For K=0, there is no changepoint and λT=λ(1) for all T=1,…,n [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:(5)γ~N(0,σγ2I),λ~Ga(χξ,δξ),ψ~S(αψ,βψ)N denotes a normal and Ga a Gamma distribution. σγ2 was set to 106, representing highly dispersed independent normal priors for each coefficient. I is an identity matrix. For λ(k), k=1,…,K+1, independent exponential distributions with mean 1/ξ and variance 1/ξ2 were specified. ξ was then assigned a gamma hyperprior Ga(χξ,δξ). The marginal prior distribution for λ(k) is then a gamma-gamma distribution [22]. In our study, we used χξ=10 and δξ=10, which corresponded to the gamma-gamma marginal of λ(k) turning out to be an F-distribution with (2,20) degrees of freedom, which then indicates that the marginal prior probability of an outbreak occurring (i.e., λ(k)≥1) is 0.39, while always favouring smaller values of λ(k), 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, Ga(αψ,βψ).αψ 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:(6)MAPE=1n∑T=1n|fT-FTFT|,
where n 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: αψ=0.1 and βψ=0.1, Model 2: αψ=10 and βψ=1, Model 3: χξ=1 and δξ=1, and Model 4: χξ=10 and δξ=1. 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 χξ=10 and δξ=1) to 0.5 (χξ=1 and δξ=1).


## 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

IffT represents the number of cases of DF in week T, then AR relates this observation to an earlier fT-J, where J=1,2,…,T-1. MA relates the error (defined as the difference between the observed, f, and that predicted, F, notifications) at week T to week (T-K), where K=1,2,…. D allows the differenced series, ΔT=(fT-fT-L), to be modelled in the event of nonstationarity in the time series, where L=0,1,2,…. Here J, K, and L are the “orders” of the respective ARIMA components. Partial autocorrelation (PAC) and autocorrelation (AC) plots are used to determine J and K, respectively.We describe the ARIMA (3,1,1) model equation used in our analysis. The number of cases of DF at weekT is denoted as fT, where T is the first week for which DF is to be predicted(1)FT=μ+fT-1+φ1fT-1+φ2fT-2+φ3fT-3+θεT-1+εT,
where FT is the predicted number of DF cases for week T, and fT-1, fT-2, and fT-3 are the DF counts in the three immediate preceding weeks, termed lag 1, lag 2, and lag 3, respectively, and εT, εT-1 are the error term at time T and T-1, respectively. In essence, we used observed values up till time T-1 to predict for dengue fever cases at week T. μ is a constant and φ1, φ2, and φ3 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,x, and epidemic, y, components of DF such that the number of cases observed fT=xT+yT and the corresponding prediction model is formulated as(2)FT=XT+YT.
Here XT and YT have independent Poisson distributions with a composite parameters (ωTνT) and (ωTλTFT-1), in which ωT handles over dispersion, hence FT is also Poisson with parameter ωT[νT+λTFT-1]. This in turn corresponds to a negative binomial distribution with dispersion parameter ψ. The mixing parameter, ωT, is assumed to have a Gamma distribution with parameters (ψ+FT) and (ψ+νT+λTFT-1).The endemic parameter,νT, is modelled as a harmonic wave (to handle strong seasonality inherent in infectious disease surveillance data) with(3)log⁡⁡νT=γ0+γ1sin⁡⁡(2πT52)+γ2cos⁡⁡(2πT52),
see [15], where 2π/52 is the base frequency of the curve, which is suitable for weekly data and γ0 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λ=(λ1,…,λn), which is assumed to be a piecewise constant [15] with unknown number of location K and unknown location of the changepoints θ1<⋯<θK, that is,
(4)λT={λ(1),ifT=1,2,…,θ1,λ(k),ifT=θk-1+1,…,θk,λ(K+1),ifT=θK+1,…,n,
where θ1<θ2<⋯<θK are the K unknown changepoints, such that θk∈{1,2,…,n-1) for all k∈(1,2,…,K). For K=0, there is no changepoint and λT=λ(1) for all T=1,…,n [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:(5)γ~N(0,σγ2I),λ~Ga(χξ,δξ),ψ~S(αψ,βψ)N denotes a normal and Ga a Gamma distribution. σγ2 was set to 106, representing highly dispersed independent normal priors for each coefficient. I is an identity matrix. For λ(k), k=1,…,K+1, independent exponential distributions with mean 1/ξ and variance 1/ξ2 were specified. ξ was then assigned a gamma hyperprior Ga(χξ,δξ). The marginal prior distribution for λ(k) is then a gamma-gamma distribution [22]. In our study, we used χξ=10 and δξ=10, which corresponded to the gamma-gamma marginal of λ(k) turning out to be an F-distribution with (2,20) degrees of freedom, which then indicates that the marginal prior probability of an outbreak occurring (i.e., λ(k)≥1) is 0.39, while always favouring smaller values of λ(k), 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, Ga(αψ,βψ).αψ 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:(6)MAPE=1n∑T=1n|fT-FTFT|,
where n 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: αψ=0.1 and βψ=0.1, Model 2: αψ=10 and βψ=1, Model 3: χξ=1 and δξ=1, and Model 4: χξ=10 and δξ=1. 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 χξ=10 and δξ=1) to 0.5 (χξ=1 and δξ=1).


## 3. Results

Figure1 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 (Figure2(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.Plots of autocorrelation and partial correlation for dengue fever (DF).
(a)
Correlogram(b)
Partial correlogramWe explored various formulations of the ARIMA model, and we summarise some of the more important ones in Table1. 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 AR(1)=-0.10 (P=0.001), AR(2)=0.10 (P=0.002), and AR(3)=0.23 (P<0.001) 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.
ModelModel specificationMAPE1ARIMA (1,0,0)23.612ARIMA (2,0,0)23.093ARIMA (3,0,0)23.204ARIMA (4,0,0)23.235ARIMA (3,1,0)19.866ARIMA (3,1,1)19.96Table 2
Parameters for the final models.
ARIMA modelCoefficient95% confidence intervalP valueConstant (μ)0.28−3.864.410.896AR 1 (φ1)−0.10−0.16−0.040.001AR 2 (φ2)0.100.040.170.002AR 3 (φ3)0.230.170.29<0.001K-H modelCoefficient95% credible intervalψ25.118.432.3γ03.31.93.6γ1−0.2−0.30.3γ2−0.5−0.7−0.4ξ1.01.01.0K7.60.0115.0λ1.30.72.0The comparison between the ARIMA and K-H model is shown in Figure3. 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.
MAPEARIMAK-HOverall17.5417.21Stratified (in 4 week intervals)Year 2007Weeks 1 to 417.0714.275 to 828.6025.629 to 1233.4130.6313 to 1632.5233.0917 to 2021.8320.5321 to 2420.6419.7625 to 2812.8613.2229 to 3211.5314.4033 to 368.5410.2637 to 405.076.5041 to 4418.4917.4245 to 488.5410.3549 to 5215.7012.44Year 2008Weeks 1 to 411.1311.165 to 829.0925.639 to 1216.3919.4113 to 1615.5110.7717 to 2019.2118.3821 to 249.8310.07Table 4
Sensitivity analysis on K-H model parameters.
MAPEInitial K-H modelSensitivity analysis1234Overall17.2117.7117.7117.5016.54Stratified (in 4 week intervals)Year 2007Weeks 1 to 414.2720.1222.3320.3020.035 to 825.6225.4125.6625.1223.399 to 1230.6331.0631.3030.9531.0713 to 1633.0933.2032.3132.4927.8917 to 2020.5320.4120.4020.9221.8221 to 2419.7621.1420.9021.2521.5825 to 2813.2213.4514.1813.2812.9129 to 3214.4013.9813.0413.3910.6533 to 3610.2610.769.9710.556.1137 to 406.506.696.395.543.3041 to 4417.4217.6217.5416.9115.9445 to 4810.3511.3010.5911.0910.6749 to 5212.4412.9912.3712.5813.12Year 2008Weeks 1 to 411.1611.0910.8511.1311.315 to 825.6325.5326.0125.4925.839 to 1219.4120.0320.2519.3116.2813 to 1610.7710.7210.4710.7211.2517 to 2018.3819.2017.9718.5919.2721 to 2410.079.3310.989.959.92A 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 (Figure3). 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 Table4, 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.



---

*Source: 758674-2012-03-08.xml*