Applying SARIMA, ETS, and hybrid models for prediction of tuberculosis incidence rate in Taiwan

Data source

Data of the TB cases and population were obtained respectively from the websites of the Taiwan Center for Disease Control (CDC) (Taiwan National Infectious Disease Statistics System, https://nidss.cdc.gov.tw/ch/SingleDisease.aspx?dc=1&dt=3&disease=010) and the Taiwan Ministry of the Interior (Ministry of Interior Taiwan. Statistics: annual report of the Ministry of the Interior, Table S2). Then, the TB morbidity was subsequently calculated. All TB cases were initially diagnosed by clinical symptoms as well as confirmed by bacteriological and pathological examinations. In Taiwan, TB is a nationally notifiable disease, and hospital physicians must report cases promptly, i.e., cases of TB must be registered within 7 days, via the Taiwan TB registry system authorized by the Taiwan CDC for national-level surveillance purposes (Taiwan Centers for Disease Control, Department of Health. Taiwan Tuberculosis Control Report, https://www.cdc.gov.tw/En/InfectionReport/List/SOIzsdQ5fRn3xPZOleIb0w).

The study period included cases diagnosed from 2005 January to 2020 December, and TB incidence data set (Annex A: Availability of raw data) was aggregated and analyzed monthly based on the public data from confirmed cases. The data were divided into training data and validation datasets. Then, respectively, generated simulating models based on the data from 2005 January to 2017 June (in-sample, training dataset). And, the forecasting testing model was based on the data from 2017 July to 2020 December (pseudo-out-of-sample, forecast test dataset). Thereafter, we utilized the optimizing models to project the future prediction toward 2025 December.

STL decomposition

The Seasonal and Trend decomposition using Loess analysis was adopted for decomposing time series into trend, seasonal, and remainder components which developed by Cleveland, Cleveland & McRae (1990).

Construct of the SARIMA and SARIMA-ETS hybrid model

Evaluation metrics

To evaluate the performance of the SARIMA, ETS or SARIMA-ETS models, the fitted values were tested. Several performance measures were employed in determining the prediction efficiency of the automatic ARIMA model, namely, the Akaike information criterion (AIC) (Aho, Derryberry & Peterson, 2014), root means square error (RMSE) (Woschnagg, 2004), mean absolute error (MAE) (Woschnagg, 2004), mean absolute percentage error (MAPE) (Woschnagg, 2004) and the mean absolute scaled error (MASE) (Woschnagg, 2004; Hyndman, 2006). These metrics have been used by many researchers to evaluate accuracy. For these metrics, the smallest values correspond to the optimal methods.

$$\mathrm{M}\mathrm{A}\mathrm{P}\mathrm{E}=(1/\mathrm{s}\mathrm{a}\mathrm{m}\mathrm{p}\mathrm{l}\mathrm{e}\mathrm{s}\mathrm{i}\mathrm{z}\mathrm{e})\times \sum [(|\mathrm{a}\mathrm{c}\mathrm{t}\mathrm{u}\mathrm{a}\mathrm{l}-\mathrm{f}\mathrm{o}\mathrm{r}\mathrm{e}\mathrm{c}\mathrm{a}\mathrm{s}\mathrm{t}|)/|\mathrm{a}\mathrm{c}\mathrm{t}\mathrm{u}\mathrm{a}\mathrm{l}|]\times 100\mathrm{\%}$$

$\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}=\sqrt{\sum {\left(\mathrm{a}\mathrm{c}\mathrm{t}\mathrm{u}\mathrm{a}\mathrm{l}-\mathrm{f}\mathrm{o}\mathrm{r}\mathrm{e}\mathrm{c}\mathrm{a}\mathrm{s}\mathrm{t}\right)}^2\times \left(1/\mathrm{s}\mathrm{a}\mathrm{m}\mathrm{p}\mathrm{l}\mathrm{e}\mathrm{s}\mathrm{i}\mathrm{z}\mathrm{e}\right)}$

$\mathrm{M}\mathrm{A}\mathrm{E}=\sum (|\mathrm{a}\mathrm{c}\mathrm{t}\mathrm{u}\mathrm{a}\mathrm{l}-\mathrm{f}\mathrm{o}\mathrm{r}\mathrm{e}\mathrm{c}\mathrm{a}\mathrm{s}\mathrm{t}|)\times (1/\mathrm{s}\mathrm{a}\mathrm{m}\mathrm{p}\mathrm{l}\mathrm{e}\mathrm{s}\mathrm{i}\mathrm{z}\mathrm{e})$

$\mathrm{M}\mathrm{A}\mathrm{S}\mathrm{E}=(1/\mathrm{s}\mathrm{a}\mathrm{m}\mathrm{p}\mathrm{l}\mathrm{e}\mathrm{s}\mathrm{i}\mathrm{z}\mathrm{e})\sum [|\mathrm{e}|/(1/\mathrm{s}\mathrm{a}\mathrm{m}\mathrm{p}\mathrm{l}\mathrm{e}\mathrm{s}\mathrm{i}\mathrm{z}\mathrm{e})\times \sum |\mathrm{a}\mathrm{c}\mathrm{t}\mathrm{u}\mathrm{a}\mathrm{l}-\mathrm{f}\mathrm{o}\mathrm{r}\mathrm{e}\mathrm{c}\mathrm{a}\mathrm{s}\mathrm{t}|]$

Software

We applied the time-series model, the forecast and related packages in R (version 4.0.1). The used functions and hyperparameters were listed in Table 1.

General information and analysis of the time series

There were 191.526 newly confirmed TB cases in Taiwan, with an average monthly incidence rate of 4.392 (range: 2.287–6.937) per 100,000 population per month, during the study period of 192 months, i.e., from January 2005 to December 2020. Through the STL analysis, we isolated seasonality and overall trend components from the monthly TB data series and also eliminated part of the random noise or reminder components. As shown in Fig. 1, the outcome of the seasonal-trend decomposition method, yields an overall trend with several composers including a 12-month stochastic seasonality; a downward overall trend with highest in 2005 (72.5 per 100,000 person-year) and lowest in 2020 (33.2 per 100,000 person-year) as well as a periodically change of disease incidences (Figs. 1A, 1C). Also applied a box plot suggested that there was seasonality in the data. It could be seen obviously that the variation of TB in Taiwan showed the periodicity within 1 year (12 months) being a cycle with obvious seasonality during the variation process in each year with two peaks occurring in March and May as well as a trough from December to February in each cycle (Figs. 1B, 1C).

The incidence data set was used to analyze and formulate the simulating models. We applied the time series including 192 datasets of logarithmic monthly incidence rates which counterbalanced by an order of 12 at the seasonal level. The outcome of an Augmented Dickey-Fuller test (Dickey-Fuller = −7.6315, Lag order = 1, P_Dickey-Full = 0.01 of alternative hypothesis: stationary) and Ljung-Box test (X-squared = 125.87, df = 1, P_Ljung-box = <2.2e−16 of alternative hypothesis: being not independent) which showed that this time series of from January 2005–December 2020 was with stationary and serial correlation format (Fig. 1A).

Construct and compare the accuracy metrics of modeling performances: simulating training by in-sample sets and evaluate or validate forecasting testing by pseudo out-of-sample sets

For constructing and evaluating the accuracy of the modellings, the simulating processing, and forecasting performances were assessed and verified, respectively. Beforehand, the data series samples were divided into the training sets (in samples) i.e., the former 150 datasets from January 2005 to June 2017, and the forecasting test sets (pseudo-out-of-samples) i.e., the following 42 data series from July 2017 to December 2020. Next, three training models were constructed utilized 150 datasets. Firstly, a best-simulating ARIMA in samples identified the seasonal ARIMA(1,0,0)(2,1,0)₁₂ mode was selected as the best performing ARIMA model with the most minimum values of AIC = −335.62, BIC = −320.99. According to its Ljung-Box Q test which showed the fitness of the ARIMA(1,0,0)(2,1,0)₁₂ models with the residual error, the sequence was achieved white noise of Q* = 62.283, p-value = 3.139e−06. Secondly, a ETS(A,A,A) model (AIC = −68.44, BIC = −17.26) was selected as the best performing ETS model with the Ljung-Box Q testing showed that Q* = 46.357, p-value = 2.033e−07. Thirdly, a hybrid model was constructed and comprised of the following models: ARIMA with a weight 0.5 and ETS with a weight 0.5 (Tables 1, 2).

Subsequently, the forecasting tests were performed, i.e., a 6-to-30-months ahead forecasting steps were computed. Thus, the modeling’ accuracy metrics were estimated respectively based on these yield forecasting values from the simulating model constructed by the 150 in-sample training datasets. Overall, several metrics were employed to evaluate the performance of these preferred models including their in-sample simulating and out-of-sample forecasting for comparison. Regarding the performance of the simulated fitting, the RMSE, MAE, MAPE and MASE were 0.0585, 0.0462, 0.4613 and 0.5982, respectively, in the SARIMA-ETS-hybrid model vs 0.0652%, 0.0514%, 0.5127%, and 0.6649%, respectively, in the SARIMA(1, 0, 0)(2, 1, 0)₁₂ model. Regarding the performance of the out-of-sample forecasting performance, the RMSE, MAE, MAPE and MASE were 0.084%, 0.067%, 0.646%, and 0.870%, respectively, in the SARIMA-ETS model vs 0.0924%, 0.0742%, 0.7135% and 0.9604%, respectively, in the SARIMA(1, 0, 0)(2, 1, 0)₁₂ model (Table 2). Among these modeling outcome assessments, the SARIMA-ETS-hybrid model outperformed the SARIMA model in near future forecasting i.e., 3–12 months ahead based on the accuracy metrics shown. Furthermore, both models showed that the short-term predictions were better than the long-term predictions, e.g., the MAPE of the first 3–12 months in both models’ forecasting test was better than the MAPE latter 18–28 months, as shown in Table 2.

Simulating with SARIMA, ETS, and SARIMA-ETS hybrid model for the best-fitting model

Analyzing by an ARIMA programming function utilized the time series of 192 datasets, the optimal modeling was established i.e., a seasonal ARIMA(3,0,0)(2,1,0)₁₂ model which was stepwise screened as the best performing SARIMA model (Table 1). This model was automatically taking the correlations between the ACF and PACF graphs of the residual sequence, AIC, AICc, and BIC into consideration. With the outcome of the ADF test (Dickey-Fuller = −9.1951, Lag order = 1, P_Dickey-Full = 0.01, stationary), Ljung-Box test (P_Ljung-box = 0.003915), it showed that the examined data was approximately independent and normally distributed with zero means and variances, the residual series successfully attained white noise, and the residual series values of the minimum AIC, AICc, and SBC were −404.25, −403.54, and −382.42, respectively (Scheme 1 and Table 1).

Simulating by an ETS programming function for the time series datasets of 192 monthly tuberculosis incidence rates, an ETS(A,A,A) model was appropriated to be captured (AIC = −46.107399, AICc = −42.306157, BIC = 8.078159; Dickey-Fuller = −9.9554, Lag order = 1, p-value = 0.01, stationary) (Dickey-Fuller = −9.9554, Lag order = 1, P_Dickey-Full = 0.01, alternative hypothesis: stationary; Ljung-Box test, data: Residuals from ETS(A,A,A), Q* = 56.673, df = 8, P_Box-Ljung = 2.086e−09, Model df: 16, Total lags used: 24) (Table 1).

While establishing the ARIMA-ETS-hybrid models, a combined model comprised of weight 0.5 of ARIMA and weight 0.5 of ETS was constructed. The outcomes were subsequently calculated.

Projection of TB incidence

Base on the verification of each best-fitting model, we paralleled utilizing the time series of 192 datasets to a long-distance projecting of forecasting for the future. The future TB incidence rates were predicted, which was projecting ahead to the future forecasting more than 3 years over the next several years. After exponentiating the resulting logarithmic TB incidence values, the projection outcomes of the next several-year TB incidence rates based on the outcome of ETS with best performance than SARIMA-ETS and SARIMA was 26.586 (95% CI [21.601–32.721]) per 100,000-year in 2025 and 20.706 (95% CI [16.328–26.258]) per 100,000-year in 2030 (Table 3, Scheme 1). Moreover, the projection outcome of the future-year TB incidences based on the ETS model showed a 41.82% (range: 28.40–52.72%) reduction in 2025 and a 54.69% (range: 42.54–54.77%) reduction in 2030 compared with 2015 levels (Table 3, Scheme 1).

Limitations

There are some limitations to our forecasting models. First, the predictivity of our models decreased with longer time distances, and the confidence interval become wider, suggesting that our models are better at predicting the near epidemics than the long-term trend prediction. By integrating the forecasts and manipulating process by human judgment for future epidemics might restructure the working states. Second, using other new deep learning algorithms, may help to improve the forecasting accuracy.