Forecasting the incidence of tuberculosis in China using the seasonal auto-regressive integrated moving average (SARIMA) model.

OBJECTIVES: The aims of this study were to develop a forecasting model for the incidence of tuberculosis (TB) and analyze the seasonality of infections in China; and to provide a useful tool for formulating intervention programs and allocating medical resources.
METHODS: Data for the monthly incidence of TB from January 2004 to December 2015 were obtained from the National Scientific Data Sharing Platform for Population and Health (China). The Box-Jenkins method was applied to fit a seasonal auto-regressive integrated moving average (SARIMA) model to forecast the incidence of TB over the subsequent six months.
RESULTS: During the study period of 144 months, 12,321,559 TB cases were reported in China, with an average monthly incidence of 6.4426 per 100,000 of the population. The monthly incidence of TB showed a clear 12-month cycle, and a seasonality with two peaks occurring in January and March and a trough in December. The best-fit model was SARIMA (1,0,0)(0,1,1)12, which demonstrated adequate information extraction (white noise test, p>0.05). Based on the analysis, the incidence of TB from January to June 2016 were 6.6335, 4.7208, 5.8193, 5.5474, 5.2202 and 4.9156 per 100,000 of the population, respectively.
CONCLUSIONS: According to the seasonal pattern of TB incidence in China, the SARIMA model was proposed as a useful tool for monitoring epidemics.

Introduction

Tuberculosis (TB), caused by the Mycobacterium tuberculosis complex, is a chronic respiratory infectious disease that is transmitted most commonly through cough aerosols [1], [2]. Despite heroic efforts to combat this devastating disease, TB remains a major global health issue with a vast health burden due to the high incidence, medical expenses, drug resistance and co-infections. In 2015, it was estimated that there were 10.4 million new infections and 1.8 million deaths caused by TB worldwide [3], [4]. In China, owing to the lack of healthcare resources and absence of essential medical services, TB has spread widely with serious harmful effects on the health of the population [5]. In the 1990s, the World Health Organization (WHO) recommended the Directly Observed Treatment Short course (DOTS) strategy comprising two key factors: directly observed treatments and short course chemotherapy. This strategy was implemented in 13 provinces, inhabited by more than 560 million people [6]. With the arrival of the 21st century, the Chinese government has made great efforts to achieve comprehensive control of TB by devising and implementing the State Council of China National TB Control Program (2001–2010) in 2001, and expanding the DOTS program to the entire country by 2005 [7]. For the sustainable development of TB control, the State Council of China issued a new TB prevention and control program in the “12th Five-Year Plan” [8]. However, China remained one of the three countries with the highest burden of TB cases in 2015 (combined, 45% of the global total), which was only marginally lower than India [3].

In response to the health issues associated with TB infection, we should prioritize optimization of the allocation of health resources and medical services, exploration of the temporal dynamic changes of transmission and prevalence, and prediction of the future incidence. To achieve this, mathematical and statistical models are required to forecast TB occurrence as an early warning system. The auto-regressive integrated moving average (ARIMA) is a frequently used model for forecasting the incidence of epidemic diseases, including TB. This system is used to model the temporal dependence structure of a time series, taking into account changing trends, periodic changes, and random disturbances in the time series [9], [10], [11], [12], [13], [14], [15], [16].

The aim of this study was to develop the seasonal auto-regressive integrated moving average (SARIMA) model for short-term prediction of the incidence of TB using the data from January 2004 to December 2015. This model will contribute to the provision of a scientific reference of risk assessment for TB prevention and control.

Materials and method

Data collection

Data for the monthly incidence of TB were obtained from the National Scientific Data Sharing Platform for Population and Health of China (excluding Taiwan, Hong Kong and Macao). All TB cases were verified by clinical or laboratory diagnosis. TB is a class B notifiable disease in China and must be reported within 12 h of diagnosis. After the severe acute respiratory syndrome (SARS) outbreak in 2003, medical and health institutions at different levels emphasized importance to reporting cases of TB. Thus, it is believed that the incidence of TB from this platform is relatively credible.

SARIMA model

The seasonal auto-regressive integrated moving average (SARIMA) model is a time series forecasting method proposed by Box and Jenkins in 1970s [17]. The general form of the SARIMA model is as follows: (p,d,q)(P,D,Q), where p and q are the orders of auto-regressive (AR) and moving average (MA), respectively, d is the order of the differences, P, D and Q are the corresponding seasonal orders, and S represents the steps of the seasonal differences, which can be calculated according to following the equations [18]:With Where B is the backward shift operator, Z is the observed value at the time t (t  = 1, 2, … k), and ε is the residual error at time t.

The processes of the SARIMA model involve four steps:

Sample pre-treatment: No-white noise, stationarity and seasonality of a time series are a pre-condition for building a SARIMA model. In the case of a no-white noise time series, there is a correlation with observed values with a non-random distribution, and can be used to build a model. As the most common method of unit root testing, an Augmented Dickey–Fuller (ADF) test was adopted to determine whether the time series is stationary. If the time series is non-stationary, differencing could be used effectively to remove this trend. Based on the sequence diagrams and practical experiences, the seasonality of the time series can be explored.

Identification and estimation: Determining the orders of the SARIMA model on the basis of graphs for the auto-correlation function (ACF) and partial auto-correlation function (PACF). The least square method was used to estimate the model parameters.

Diagnosis and optimization: The availability of the model required verification by parameters and white noise tests. The Q-test was employed to diagnose the white noise of the residual errors, which are independent and normally distributed. The significance of a single parameter was evaluated by means of a t-test.

Prediction: The obtained model was applied to fit the incidence of TB from January 2004 to December 2015. Finally, the optimal model was determined with the help of accuracy tests, and used for predicting the incidence of TB.

Accuracy test

To measure the accuracy of the SARIMA models, the fitted values should be tested using the Schwarz criterion (BC), the error range and the deviation degree. The error range indexes include: root mean squared error (RMSE), mean absolute error (MAE), mean absolute percent error (MAPE) and the Theil inequality coefficient (TIC). The deviation degree indexes include: bias proportion (BP) and variance proportion (VP). The results of the accuracy testing indexes are shown in Table 2 .

Table 2

Comparison of the accuracy of the SARIMA models.

Index	Models
	SARIMA (1,0,0)(1,1,1)12	SARIMA (1,0,0)(0,1,1)12	SARIMA (2,0,0)(0,1,1)12
SC	0.6888	1.5896	1.5931
RMSE	0.4979	0.5066	1.1014
MAE	0.2099	0.3055	0.8557
MAPE	3.4705	4.5347	14.1629
TIC	0.0245	0.0384	0.0882
BP	<0.0001	0.0048	0.5039
VP	0.0544	0.0058	0.0225

SC, Schwarz criterion; RMSE, root mean squared error; MAE, mean absolute error; MAPE, mean absolute percent error; TIC, Theil inequality coefficient; BP, bias proportion; VP, variance proportion.

Statistical analysis

Data analysis was performed with SPSS for Windows, version 22.0 (SPSS Inc., Chicago, IL, USA). EViews for Windows, version 8.0 (IHS Global Inc., Irvine, CA, USA) was for SARIMA modeling and predictions. p <  0.05 was considered to indicate statistical significance.

Results

Characteristics of the TB incidence

Up to December 31, 2015, 12,321,559 TB cases were reported in China, with an average monthly incidence of 6.4426 per 100,000 of the population. As shown in Fig. 1 , the incidence of TB from 2004 to 2015 exhibited an uneven distribution at the province level. Overall, the highest incidence was reported in Xinjiang and Guizhou, while the lowest incidence was observed in Beijing, Tianjing, Hebei, Liaoning, Shandong, Jiangsu, Shanghai, Yunnan and Ningxia . Fig. 2 summarizes the distribution of TB cases in terms of age groups and gender in China from 2004 to 2015. The pattern shows that a higher incidence of male cases compared with female cases in the different age groups. Furthermore, the highest incidence of TB cases was reported in the 20–24 years age group, which accounted for 10.21% (1,257,679 cases) over the study period. Of 12,321,559 TB cases, 8,563,602 (69.50%) were males and 3,757,957 (30.50%) were females, with a male: female ratio of 2.28:1. As shown in Fig. 3 , there was an obscure decline in the monthly distribution of TB incidence, with obvious seasonality. Apparently, the highest incidence peak appeared in January followed by smaller peak in March, which accounted for 10.04% (1,236,741 cases) and 9.96% (1,227,076 cases) of all reported cases, respectively. After March, the incidence of TB declined sharply until December, with the proportion falling to 5.36% (660,416) of all reported cases.

Fig. 1

Distribution of the reported tuberculosis incidence at the province level in China (2004–2015).

Fig. 2

Distribution of the tuberculosis cases in terms of age groups and gender in China (2004–2015).

Fig. 3

Monthly incidence of tuberculosis and the declining trend in China between 2004 and 2015.

Model recognition and diagnosis

As shown in Fig. 3, distinctly, this time series exhibited seasonality, suggesting that the SARIMA model should be used to fit the data. Moreover, the series exhibited the same tendency when the Augmented Dickey–Fuller (ADF) test was applied (t  = −0.6880, p  = 0.8451). It was noted that the series could be transformed to exhibit stationarity of seasonal difference. Therefore, the series was differenced by an order of 12 at the seasonal level to induce stationarity. ADF tests (t  = −4.5632, p  = 0.0003) indicated that the series was steady, in accordance with a non-zero mean (t  = 53.3912, p  = 0.0001) white noise sequence (Q-test, p >  0.05; Fig. 4 ).

Fig. 4

ACF and PACF diagram of monthly incidence of tuberculosis in China (2004–2015) after one season of lag 1 difference.

As shown in Fig. 4, the ACF declined gradually, while the PACF approached zero at all lags exceeding two, demonstrating the truncation status of the PACF. Based on the characteristics of the ACF and PACF distribution, we selected four models provisionally; SARIMA (1,0,0)(1,1,1)12, SARIMA (1,0,0)(0,1,1)12, SARIMA (2,0,0)(1,1,1)12 and SARIMA (2,0,0)(0,1,1)12. However, SARIMA (2,0,0)(1,1,1)12, with non-seasonal AR(2) that was not statistically significant (t  = −1.3660, p  = 0.1747), was identified by parametric estimation (Table 1). By comparison of the advantageous and statistic validation (Table 2), the potential models of the SARIMA (1,0,0)(1,1,1)12 and SARIMA (1,0,0)(0,1,1)12 were used for fitting the data.

Table 1

Coefficient, standard error, t-Statistic and p-values of the SARIMA models of the parameters.

Variable	Coefficient	Standard error	t-Statistic	p-value
Parameters of the SARIMA model (1,0,0)(1,1,1)12
C	−0.2553	0.1021	−2.5013	0.0138
Non-seasonal AR(1)	0.6426	0.0725	8.8665	<0.0001
Seasonal AR(1)	−0.4561	0.0467	−9.7579	<0.0001
Seasonal MA(1)	0.9269	0.0150	61.6900	<0.0001
Parameters of the SARIMA model (1,0,0)(0,1,1)12
C	−0.2064	0.0778	−2.6565	0.0089
Non-seasonal AR(1)	0.6767	0.0640	10.5822	<0.0001
Seasonal MA(1)	−0.5011	0.0744	−6.7373	<0.0001
Parameters of the SARIMA model (2,0,0)(1,1,1)12
C	−0.2385	0.1007	−2.3692	0.0195
Non-seasonal AR(1)	0.7527	0.0915	8.2265	<0.0001
Non-seasonal AR(2)	−0.1257	0.0920	−1.3660	0.1747
Seasonal AR(1)	−0.4031	0.0480	−8.3926	<0.0001
Seasonal MA(1)	0.9219	0.0164	56.2240	<0.0001
Parameters of the SARIMA model (2,0,0)(0,1,1)12
C	−0.2208	0.0863	−2.5577	0.0117
Non-seasonal AR(1)	0.5378	0.0890	6.0495	<0.0001
Non-seasonal AR(2)	0.1839	0.0860	2.1372	0.0345
Seasonal MA(1)	−0.5261	0.0740	−7.1113	<0.0001

SARIMA, seasonal auto-regressive integrated moving average; C, constant terms; AR, auto-regressive; MA, moving average.

Prediction of TB incidence over the next six months

SARIMA (1,0,0)(0,1,1)12 was the most suitable model due to the absence of statistical significance of the white noise diagnostic of the residuals (Q-test, p >  0.05). The values predicted by the model and the observed values for the incidence of TB are presented in Fig. 5 . The MAPE value was 4.5347. The monthly incidences were 6.6335, 4.7208, 5.8193, 5.5474, 5.2202 and 4.9156 per 100,000 of the population, respectively.

Fig. 5

SARIMA (1,0,0)(0,1,1)12 model fitting, verification and forecasting of tuberculosis in China from January 2004 to June 2016.

Discussion

This study revealed that the SARIMA (1,0,0)(0,1,1)12 was the best-fit mathematical model for forecasting monthly TB incidence based on the data from January 2004 to December 2015 in China. The forecast results indicated that the TB incidence was likely to increase only slightly over the subsequent six months. The factors that contribute to a rise in the incidence of TB are complex and further research is required to clarify this issue. The SARIMA (1,0,0)(0,1,1)12 clearly represents a useful tool for monitoring TB incidence and forecasting TB epidemics. Thus, it is important that the health departments adopt this model as the basis of an early warning system that will allow the timely application of enhanced surveillance and reallocation of medical resources [19]. On 16th February 2017, the General Office of the State Council of China issued the “13th Five-Year Plan” including the National Tuberculosis Prevention and Control Planning, which proposed the full use of information management [20]. In particular, this measure will be conducive to the implementation of a high efficiency forecasting model for epidemic surveillance and as a treatment reference.

The SARIMA model has been applied in previous studies to predict the incidence of TB in China [21], [22] with two limitations. First, the number of observed values was insufficient for accurate forecasting, and second, the precision of the obtained models was not comparable with that of the proposed models. In the current study, 144 observations were collected and seven indexes were adopted to generate a robust and valid SARIMA model.

In this study, TB incidence was recorded at monthly intervals to allow analysis of the seasonal distribution. From an overall perspective, the monthly incidence of TB from January 2004 to December 2015 peaked in January and March, with a trough in December. Other studies have shown that the seasonal distribution of TB incidence in China is similar to that in other countries. For instance, the peak of TB prevalence in Portugal is March, with a trough is December [23]. In the United States as a whole, the peak prevalence is also in March, while the trough occurs in November [24], while in New York, the incidence of TB peaked in spring or summer, with the trough occurring in the fall [25]. In Iran, the peak TB incidence appeared in spring [26]. The peak incidence of TB is mainly caused by new disease spread and increased incidence of latent infections; therefore, studies of the peak seasonal incidence peak should take into account the delay in diagnosis [27]. A previous study showed that the median TB patient delay is 93 days (range 68–128 days) [28]. Missed the early diagnosis and variable incubation periods of TB, will not only lead to a longer infectious period and increased transmission rates, but will also limit the reliability of the TB surveillance system and the effects of an early warning system [19], [29]. Furthermore, it is important to note that the Spring Festival in China contributes to the increase in the incidence of TB in March by virtue of the increased travel around the Chinese lunar new year and population mobility [27], [30].

Some limitations of this model should also be noted when interpreting the forecasting results. First, the natural environment and socioeconomic factors that influence the incidence of TB; however, due to data availability, and the focus on time series, these factors were not taken into account in this study. Second, the SARIMA model is suitable only for short-term forecasting, and the model requires dynamic updating with new data to ensure the accuracy and stability of the forecasting. Third, the seasonality of TB incidence in smaller areas was not considered owing to the lack of available data and vast territory involved. In addition, the data were obtained from the National Scientific Date Sharing Platform for Population and Health, which is a passive surveillance system that is subject to potential biases associated with under-reporting of TB cases and might affect the precision of the predictions.

Conclusion

In this study, based on the seasonal pattern of TB incidence in China, we proposed the SARIMA model as a useful tool for monitoring epidemics. The results of our study will be beneficial in public health management of strategies implemented for the prevention and control of TB. Furthermore, this information will contribute to the achievement of the goals of the National Tuberculosis Prevention and Control Planning included in the “13th Five-Year Plan”.

Funding

No funding sources.

Competing interests

None declared.

Ethical approval

Not required.