Forecasting the Tuberculosis Incidence Using a Novel Ensemble Empirical Mode Decomposition-Based Data-Driven Hybrid Model in Tibet, China

Abstract

Objective

The purpose of this study is to develop a novel data-driven hybrid model by fusing ensemble empirical mode decomposition (EEMD), seasonal autoregressive integrated moving average (SARIMA), with nonlinear autoregressive artificial neural network (NARNN), called EEMD-ARIMA-NARNN model, to assess and forecast the epidemic patterns of TB in Tibet.

Methods

The TB incidence from January 2006 to December 2017 was obtained, and then the time series was partitioned into training subsamples (from January 2006 to December 2016) and testing subsamples (from January to December 2017). Among them, the training set was used to develop the EEMD-SARIMA-NARNN combined model, whereas the testing set was used to validate the forecasting performance of the model. Whilst the forecasting accuracy level of this novel method was compared with the basic SARIMA model, basic NARNN model, error-trend-seasonal (ETS) model, and traditional SARIMA-NARNN mixture model.

Results

By comparing the accuracy level of the forecasting measurements including root-mean-square error, mean absolute deviation, mean error rate, mean absolute percentage error, and root-mean-square percentage error, it was shown that the EEMD-SARIMA-NARNN combined method produced lower error rates than the others. The descriptive statistics suggested that TB was a seasonal disease, peaking in late winter and early spring and a trough in autumn and early winter, and the TB epidemic indicated a drastic increase by a factor of 1.7 from 2006 to 2017 in Tibet, with average annual percentage change of 5.8 (95% confidence intervals: 3.5–8.1).

Conclusion

This novel data-driven hybrid method can better consider both linear and nonlinear components in the TB incidence than the others used in this study, which is of great help to estimate and forecast the future epidemic trends of TB in Tibet. Besides, under present trends, strict precautionary measures are required to reduce the spread of TB in Tibet.

Keywords:

• tuberculosis
• incidence rate
• time series analysis
• statistical models
• forecasting

Introduction

Tuberculosis (TB) is a chronic communicable disease caused by the bacillus Mycobacterium tuberculosis, it is a major cause of ill-health, ranking among the 10 leading causes of morbidity and mortality globally.¹ The WHO estimated that about 10.0 million people developed TB in 2019 around the globe.² Furthermore, there were an estimated 1.4 million deaths associated with TB in 2019, leading to an annual approximately 1.8% and 2.8% reductions between 2015 and 2019 in the incidence and death rates, respectively, worldwide.² Currently, the global case notifications were mainly reported from 30 high TB burden countries, accounting for around 87% of the worldwide totals.² Of them, ranking the third is China, where the proportion of the reported cases reaches 8.4% in the world.² Moreover, China also has the second share (14%) of the global burden of drug-resistant TB in 2019, despite progress in cutting TB incidence with an average annual percent change of 3.2% between 1990 and 2017.³ However, recent years have been witnessing a risk of resurgence in the TB incidence in some regions in China, attributable to climate change, alcohol use disorders, diabetes, HIV infection, smoking, and undernourishment, which have suppressed the progress of the end TB strategy with a target of 90% and 95% reductions in the TB incidence and death rates, respectively, in the period 2015–2035.^1,2,4,5 To accelerate progress towards such an ambitious goal, accurately estimating and forecasting the past, current, and future epidemic trends of TB play a pivotal role in forming the targeted prevention and control measures and in adjusting the present interventions.

There is a large volume of studies forecasting the epidemiological trends of communicable diseases using different statistical techniques, such as the seasonal autoregressive integrated moving average (SARIMA) method,⁶ exponential smoothing method,⁷ generalized regression neural network (GRNN) method,⁸ nonlinear autoregressive neural network (NARNN) method,⁹ backpropagation neural network (BPNN) method,¹ multivariate linear regression method,¹⁰ and Elman and Jordan recurrent neural networks.¹¹ Among them, the most frequently used linear method is the SARIMA model,^12–15 whereas the nonlinear method is the NARNN model.^5,9,12,16 However, there are various factors affecting and limiting the epidemiological patterns of diseases.^1,17,18 For this reason, the epidemics of communicable diseases are often characterized by linear and nonlinear patterns.¹⁹ However, the above-mentioned statistical techniques are often insufficient to capture either tendencies or randomness in the time series, which may limit the generalization of the findings from the mentioned models. Therefore, to improve the prediction accuracy and reliability, an alternative approach is expected to be customized for use with both the above two traits included in the object series. Inspired by such an idea, many studies have adopted a hybrid technique comprising a linear model and a nonlinear model to yield more accurate prediction by integrating each model’s superiority.^9,20,21

Over the last years, wavelet analysis has become a popular method capable of decomposing the data into different high-scale and low-frequency components (linear trait) and low-scale and high-frequency components (nonlinear trait) particularly when target series shows complex nonstationary and nonlinear characteristics.²² More recently, a new wavelet analysis approach-ensemble empirical mode decomposition (EEMD) was developed by adding noise to overcome the mode mixing disadvantage of the traditional empirical mode decomposition (EMD).²³ Importantly, different from the most commonly used discrete wavelet transform techniques that may generate different decomposition results since they are required to predetermine basis functions,²⁴ EEMD is a self-adaptive, intuitive, direct, and empirical decomposition method for the nonstationary and nonlinear time series.²³ Considering the superiority of the EEMD technique in decomposing the data and the worsening epidemic trends of TB in Tibet with the highest TB morbidity in China in the past years,²⁵ we thus proposed a novel forecasting method (designated as EEMD-SARIMA-NARNN hybrid method) by integrating the SARIMA and NARNN models based on the EEMD technique to capture both the tendencies and randomness in the time series simultaneously. Specifically, the EEMD technique was first applied to decompose the TB incidence series into several Intrinsic Mode Functions (IMFs) components representing the different nonlinear parts of the data and a residue component denoting the tendency of the data. Second, using the NARNN and SARIMA methods to model the IMFs and residue components, respectively. Finally, the forecasting values from the IMFs and residue components were aggregated as the predicted results from the EEMD-SARIMA-NARNN hybrid technique. By doing so, the linear and nonlinear parts contained in the TB incidence series can be considered by this data-driven mixture method simultaneously. Whilst we also compared the fitted and predicted results from the basic SARIMA model, basic NARNN model, and the traditional SARIMA-NARNN mixture model with that from our proposed novel hybrid method to judge the forecasting performance of this novel data-driven method.

Materials and Methods

Study Area and Data Source

Tibet is located in the southwest of the People’s Republic of China. It is one of the five ethnic autonomous regions in China, and its average elevation is above 4000 m, known as the “roof of the world”. Tibet covers 1228.4 km². At the end of 2019, the resident population was 3.506 million.

The monthly incident cases between January 2006 and December 2017 reported to the Center for Disease Control and Prevention (CDC) were taken from the National Notifiable Infectious Disease Surveillance System (NNIDSS), and the population data during the same period from the Bureau of statistics of Tibet Autonomous Region (http://tjj.xizang.gov.cn/). A total of 144 observations across 12 years were aggregated. Generally, it is required to use preferably 100 observations or more to develop object models so that a reliable and stable result is obtained.²⁶ Therefore, the observations used in this research were partitioned into two parts including a training dataset comprising 132 observations (from January 2006 to December 2016) and a testing dataset comprising 12 observations (from January to December 2017).

Ethics Approval and Consent to Participate

The study protocol was approved by the study institutional review board of the Xinxiang Medical University (No: XYLL-2019,072), and consent is not needed for this study since we only collected the reported monthly TB cases publicly available without the detailed individual information.

A New Hybrid Model of EEMD-SARIMA-NARNN

EEMD Technique

In the past, albeit EMD technique shows remarkable superiority in signal analysis, particularly for the data produced in noisy nonlinear and non-stationary processes, there are two major weaknesses including the edge-effects and mode-mixing,^20,23,27 especially for the mode-mixing disadvantage, which not only results in the mixing of different scale vibration modes but can even cause the loss of the physical meaning of each IMFs.²⁸ The EEMD method was thus developed based on the EMD algorithm to overcome the barriers of the original technique.²³ In the EEMD approach, it uses a powerful method of ensemble mean to improve the accuracy of measurements.²³ Namely, it defines the actual individual IMFs components as the average of an ensemble of experiments, each IMF component is composed of the signal and noise of finite-amplitude.²⁰ The decomposition procedure of the EEMD approach can be in the form as below.

Initially, add white noise w(t) to the given time series x(t), and the obtained new series can be expressed as

(1)

Secondly, decompose the produced series into IMFs with the help of the EMD algorithm.

Thirdly, repeat the above two steps using various white noise series, the resulting results were added to the original series each time.

Finally, average the ensemble of the corresponding IMFs generated by the EMD decomposition technique.

In the process of decomposition, the numbers of the ensemble and the amplitude of the added white noise play a pivotal role in the obtained final results.²⁹ Fortunately, a well-demonstrated statistical rule helps to determine these two key parameters²³:

(2)

where N represents the numbers of the ensemble, is the amplitude of the added white noise, andsignifies the standard deviation of the error. Prior studies have found that the parameter settings with the ensemble numbers of 100 and the amplitude of added white noise of 0.2 times standard deviation can generate a satisfactory result.^23,30

EEMD-SARIMA-NARNN Hybrid Model

To this end that makes full use of the possible linear and nonlinear information contained in the target series, inspired by the idea of “decomposition and ensemble” of the EEMD technique and the powerful ability to extract nonlinear information of the NARNN model,³¹ the EEMD-SARIMA-NARNN hybrid model was thus proposed (Figure S1). In this project, the time series of TB morbidity was decomposed into different IMFs subseries and residue terms. Then, each of IMFs subseries was predicted using a suitable NARNN model, while the residue term was forecasted with an appropriate SARIMA model. As such, the final results can be derived by integrating the resulting values of the SARIMA model with the counterparts of the NARNN approaches. By doing so, the hybrid model can then attain a capacity that unearths all possible components included in the data. The specified equation of a SARIMA-NARNN combined model can be written as

(3)

(4)

where represents the estimates using the proposed hybrid method, is the estimates of residue terms using the SARIMA method, refers to the estimates of IMFs subseries using the NARNN method.

Furthermore, the SARIMA method, NARNN method, and traditional SARIMA-NARNN method were created according to the steps described in supplementary material.

Statistical Analysis

The TB epidemic trends and changes were measured with average annual percentage change (AAPC) based on the joinpoint regression program (Version 4.8.0.1). The basic SARIMA method was constructed with R software (version 3.4.3, R Development Core Team, Vienna, Austria), and the EEMD, basic NARNN, SARIMA-NARNN, and EEMD-SARIMA-NARNN methods were developed using MATLAB software (Version R2014a, MathWorks, Natick, MA, USA). The Lagrangian Multiplier (LM) method was used to test the conditional heteroskedastic behavior and volatility (ARCH effect) of the residual series from the selected four models.³² The root-mean-square error (RMSE), mean absolute deviation (MAD), mean error rate (MER), mean absolute percentage error (MAPE), and root-mean-square percentage error (RMSPE) were chosen to assess the predictive accuracy levels among models. The lower the indices’ values are, the better the models are.

(5)

(6)

(7)

(8)

(9)

where X_i signifies the observed TB values,is the estimates using the chosen approaches,refers to the mean of all the observed TB values, and N stands for the number of simulations and forecasts.

Results

Descriptive Analysis

During the period 2006–2017, the reported TB cases totaled 45,025 with notified 3753 cases annually, leading to annualized and monthly incidence rates of 121.9 and 10.2 per 100,000 persons, respectively. As illustrated in Figure S2, when the seasonal-trend decomposition procedure based on Loess (STL) technique was applied to the TB morbidity series, it was found that the TB epidemic showed a drastic increase by a factor of 1.7 from 2006 to 2017, with AAPC=5.8 (95% confidence intervals: 3.5–8.1; Z=5.8, P<0.001) (Figure S3). The highest TB morbidity rate was noted in 2017 (152.0 cases per 100,000 people), it has increased 2.4-fold compared with 2008 when the morbidity rate was lowest (64.7 cases per 100,000 people). Besides, it was observed that the TB morbidity had notable seasonal effects, with a peak occurred in January until August, especially in January and March, and a trough in September until December, particularly in December (Figures S2 and S4).

The Best-Fitting SARIMA Method

The original TB incidence series showed a unit root (ADF=−1.3, P=0.2) and marked seasonal effects (Figure S2), indicating a nonstationary series, and thus seasonally differencing once to remove the seasonal variation of the data. By doing so, the differenced data met the prerequisite of the SARIMA model development, being stationary (ADF=−3.4, P=0.001). Afterward, by plotting the ACF and PACF graphs for the stationary series, the autocorrelations at lags 1, 2, and 3 and partial autocorrelations at lags 1 and 2 were found to touch the significance levels (Figure S5), and we observed that the values of the key parameters (p, q, P, and Q) of SARIMA model exceeded 2, the model’s parameter suggested no statistical significance, and thereby we selected several plausible models with the key parameters were between 0 and 2. Table 1 demonstrates that the SARIMA(2,0,0)(2,1,0)₁₂ model seems to provide adequate modeling for the object series since this model presented the lowest values of performance indices (AIC=409.2, CAIC=409.7, and BIC=423.2) among the possible nine models, and statistical test of the estimated parameters for this method was indicated significantly different (Table 2). Moreover, both the Ljung–Box Q test and LM test showed a P-value greater than 0.05 (Figure 1A, Tables S1 and S2), suggesting that the residuals behaved like white noise and the ARCH effect had been removed from the residuals. The above tests mean that the SARIMA(2,0,0)(2,1,0)₁₂ model is appropriate for describing this time series. Therefore, this model is taken as the best.

Table 1 The Identified Nine Plausible SARIMA Methods and Their Information Criteria

Models	AIC	CAIC	BIC	Log-Likelihood
SARIMA(2,0,0)(2,1,0)12	409.22	409.74	423.16	−199.61
SARIMA(1,0,2)(2,1,0)12	410.01	410.75	426.73	−199.00
SARIMA(1,0,1)(2,1,0)12	411.83	412.36	425.77	−200.92
SARIMA(0,0,2)(2,1,0)12	440.42	440.94	454.35	−215.21
SARIMA(2,0,1)(2,1,0)12	410.70	411.44	427.42	−199.35
SARIMA(2,0,1)(0,1,1)12	415.85	416.37	429.78	−202.92
SARIMA(2,0,1)(0,1,2)12	416.17	416.91	432.90	−202.09
SARIMA(2,0,0)(1,1,1)12	415.06	415.59	429.00	−202.53
SARIMA(2,0,0)(1,1,0)12	432.80	433.14	443.95	−212.40

Abbreviations: SARIMA, seasonal autoregressive integrated moving average method; AIC, Akaike’s Information Criterion; BIC, Bayesian Information Criterion.

Table 2 Statistical Test of the Estimated Parameters for the Optimal SARIMA (2,0,0)(2,1,0)₁₂ Method

Variables	Estimates	Standard Error	t	p
AR1	0.519	0.089	5.831	<0.001
AR2	0.329	0.092	3.576	<0.001
SAR1	−0.736	0.086	−8.558	<0.001
SAR2	−0.526	0.089	−5.910	<0.001

Abbreviations: SARIMA, seasonal autoregressive integrated moving average method; AR1, autoregressive method at lag 1; AR2, autoregressive method at lag 2; SAR1, seasonal autoregressive method at lag 1; SAR2, seasonal autoregressive method at lag 2.

Figure 1 Estimated autocorrelogram and partial autocorrelogram of the residual sequence to forecast the TB epidemic patterns using (A) SARIMA model, (B) basic NARNN model, (C) traditional SARIMA-NARNN hybrid model, and (D) novel EEMD-SARIMA-NARNN hybrid model. As seen, the estimated autocorrelations and partial autocorrelations of the errors from the different models almost fell into the 95% confidence intervals (CIs) with few exceptions, such as the autocorrelations at lags 15 and 19 and partial autocorrelations at lags 15 in (A), and autocorrelations and partial autocorrelations at lag 11 in (D), these are also reasonable because higher-order autocorrelations and partial autocorrelations easily exceed the 95 CIs by chance.

The Best-Fitting NARNN Method

As described in previous studies, the ANNs models with hidden units between 10 and 20 and feedback delays between 3 and 7 can often provide a satisfactory prediction and ease overfitting to a certain degree.^33–37 Therefore, we trained a series of networks with the mentioned-above ranges of hidden units and feedback delays by trial and error so that a preferred model was discovered. By doing so, the NARNN (16, 6) model with 16 hidden units and 6 feedback delays tended to be indicated the optimal in that this network structure produced the lowest MSE values in the training (0.6), validation (4.5), and testing (3.8), and all R values were greater than 0.8 in the training (1.0), validation (0.8), testing (0.8), and the whole data (0.9) (Table 3 and Figure S6). Also, as shown in Figure S7, the modeling results of this model agreed pretty well with the observed values as most of the errors were between −1 and 1. Additionally, the Ljung–Box Q and LM-test statistics for the residual sequence yielded by this model showed a P-value greater than 0.05 (Tables S1 and S2), and the sample autocorrelation and partial autocorrelation for the in-sample forecast errors at different lags failed to exceed the significance bounds (Figure 1B), implying that the in-sample forecast errors showed a white noise and no ARCH effect. Thus, we are confident that the NARNN (16, 6) model was considered optimal in this study.

Table 3 The Estimated Parameters of the Optimal NARNN Models for the Original Series and Decomposed Series Using the EEMD Method

Target Components	Hidden Units	Delays	MSE			R
			Training	Validation	Testing	Training	Validation	Testing	Overall
Original series	16	6	0.648	4.531	3.760	0.973	0.842	0.848	0.944
Residuals of SARIMA model	19	6	0.537	8.110	2.287	0.806	0.714	0.607	0.715
IMF1	17	5	0.124	0.885	1.806	0.967	0.799	0.827	0.911
IMF2	16	5	0.052	0.101	0.226	0.992	0.979	0.971	0.988
IMF3	15	5	0.001	0.001	0.002	1.000	1.000	1.000	1.000
IMF4	16	6	0.000	0.000	0.000	1.000	1.000	1.000	1.000
IMF5	14	4	0.000	0.000	0.000	1.000	1.000	1.000	1.000
IMF6	14	5	0.000	0.000	0.000	1.000	1.000	1.000	1.000

Abbreviations: SARIMA, seasonal autoregressive integrated moving average method; NARNN, nonlinear autoregressive artificial neural network.

The Best-Fitting Traditional SARIMA-NARNN Mixture Method

In the process of constructing the traditional SARIMA-NARNN method, the residual series from the SARIMA model was utilized to train the best NARNN model. Similarly, we found that the NARNN (19, 6) method with 19 hidden units and 6 feedback delays appeared to be abundant for modeling the object series via repeated modeling as this model had a relatively great R-value and a relatively low MSR value in the training, validation, and testing data (Table 3). Further, we investigated the residual series from the SARIMA(2,0,0)(2,1,0)₁₂-NARNN (16, 6) mixture model using the Ljung–Box Q and LM tests, showing that this residual series belonged to a white noise series and also no ARCH effect was noted in the residuals (Tables S1 and S2). Besides, the sample autocorrelations and partial autocorrelations for the in-sample forecast errors inside the 95% uncertainty bounds also indicated that the errors met the assumption of uncorrelatedness (Figure 1C). Thereby, the SARIMA(2,0,0)(2,1,0)₁₂-NARNN (16, 6) mixture model is considered the best-performing in this work.

The Best-Fitting Novel EEMD-SARIMA-NARNN Mixture Method

The EEMD technique was employed to decompose the TB morbidity series. As listed in Figure 2, the original series was decomposed into six different IMFs subseries and a residue term. Subsequently, the residue term was modeled using the SARIMA method and the IMFs terms were modeled with different NARNN methods, and the resulting results showed that the SARIMA (1,2,0)(1,0,1)₁₂ method was the most suitable for modeling the residue term and NARNN(17, 5), NARNN(16, 5), NARNN(15, 5), NARNN(16, 6), NARNN(14, 4), and NARNN(14, 5) models were taken to be the best methods for modeling the IMF1~IMF6 subseries (Tables 3–4), respectively, and the diagnostic test results for the best SARIMA model and different NARNN models are shown in Figure 3 and S8-S14. Further checking results for the residuals from the new EEMD-SARIMA-NARNN mixture method are given in Tables S1-S2 and Figure 1D. The P-values for LM and Ljung–Box Q tests were greater than 0.05 (Tables S1-S2), suggesting that there was no ARCH effect in the errors and there was little evidence of non-zero autocorrelations in the in-sample forecast errors at various lags, which could also be verified by the ACF and PACF plots as the sample autocorrelations and partial autocorrelations at different lags almost fell into the 95% significance bounds (Figure 1D). The above results intimated that this new data-driven mixture model can capture the short-term dependencies of the data and is adequate for forecasting the data.

Table 4 Statistical Test of the Estimated Parameters for the Optimal SARIMA (1,2,0)(1,0,1)₁₂ Method

Variables	Estimates	Standard Error	t	p
AR1	1.000	0.002	578.647	<0.001
SAR1	0.984	0.029	34.493	<0.001
SMA1	0.692	0.152	4.564	<0.001

Abbreviations: SARIMA, seasonal autoregressive integrated moving average method; AR1, autoregressive method at lag 1; SAR1, seasonal autoregressive method at lag 1; SMA1, seasonal moving average method at lag 1.

Figure 2 Intrinsic mode functions (IMFs) and residue components of the original TB incidence series using the EEMD.

Figure 3 Estimated autocorrelogram of the residual sequence to forecast the six intrinsic mode functions (IMFs) decomposed by the EEMD method using the SARIMA-NARNN hybrid model. (A) IMF1, (B) IMF2, (C) IMF3, (D) IMF4, (E) IMF5, and (F) IMF6. For the estimated autocorrelations of the residual series from the different IMFs, all of them are inside of the 95 CIs except for the one at lag zero, meaning that there is no correlated serial in the residuals of the best-fitting SARIMA-NARNN hybrid model. Thereby, there appears to be adequate and reasonable for forecasting upcoming epidemiological trends using these best-fitting models.

Assessing the Predictive Accuracy Levels Among Models

We used the preferred individual models to predict the testing data, and then comparing the predictive accuracy levels among models from two aspects of the training and forecasting horizons based on five measures including RMSE, MAD, MER, MAPE, and RMSPE. The comparative results showed that our proposed EEMD-SARIMA-NARNN mixture model produced lower values of the above five measures compared with other methods used in this study both in training and forecasting horizons (Table 5). Likewise, visible in Figure 4, it makes sense that this new mixture technique can better capture the temporal dependency structure of the TB incidence series than the others.

Table 5 Comparisons of the Fitted and Forecasted Accuracy Levels Among the Optimal Four Methods

Models	Simulated Part					Predicted Part
	MAD	MAPE	RMSE	MER	RMSPE	MAD	MAPE	RMSE	MER	RMSPE
SARIMA model	0.977	11.222	1.216	0.100	2.702	1.444	11.486	1.658	0.114	3.805
NARNN model	0.770	12.301	1.170	0.078	0.258	1.183	9.335	1.309	0.093	0.101
SARIMA-NARNN hybrid model	0.648	8.634	0.856	0.066	0.155	1.443	11.281	1.803	0.114	0.143
EEMD-SARIMA-NARNN hybrid model	0.438	6.814	0.669	0.045	0.141	0.965	7.958	1.133	0.076	0.096
Reduced percentage
D vs A	55.169	39.280	44.984	55.000	94.782	33.172	30.716	31.665	33.333	97.477
D vs B	43.117	44.606	42.821	42.308	45.349	18.428	14.751	13.445	18.280	4.950
D vs C	32.407	21.079	21.846	31.818	9.032	33.125	29.457	37.160	33.333	32.867

Notes: A refers to the SARIMA method, B takes for the NARNN method, C signifies the traditional SARIMA-NARNN mixture method, D is the novel EEMD-SARIMA-NARNN mixture method.

Abbreviations: SARIMA, seasonal autoregressive integrated moving average method; EEMD, ensemble empirical mode decomposition; NARNN, nonlinear autoregressive artificial neural network; MAD, mean absolute deviation; MAPE, mean absolute percentage error; RMSE, root-mean-square error; MER, mean error rate; RMSPE, root-mean-square percentage error.

Figure 4 Time series plot showing the comparative results between original observations and fitted and predicted values using the selected best-fitting (A) SARIMA model, (B) NARNN model, (C) traditional SARIMA-NARNN hybrid model, and (D) novel EEMD-SARIMA-NARNN hybrid model (The curve at the left of the vertical black dotted line represents the fitted values, whereas the curve at the right denotes the forecasted values). As shown, the produced curve from the EEMD-SARIMA-NARNN hybrid model is closer to the actual both in the fitted and predicted aspects compared with the curves from the other models.

Discussion

Accurately estimating the epidemics of infectious diseases acts as a basis for the decision-making process and the rational allocation of limited health resources under dynamic demand.¹ The importance and originality of this work are that it creates a novel and more effective forecasting model for TB incidence series prediction by integrating the linear characteristics of a SARIMA method and nonlinear characteristics of the NARNN model based on the EEMD technique. By comparing the performance indicators of this novel model comprising RMSE, MAD, MER, MAPE, and RMSPE with the basic SARIMA model, basic NARNN model, and the traditional SARIMA-NARNN mixture model, it was observed that this new data-driven hybrid model showed better superiority than the others both in the fitted and predicted aspects. Often, the MAPE value is considered a useful measure to assess the forecasting accuracy of a model, with a value less than 10% indicating a highly accurate forecast, and a value less than 20% suggesting a good approximation to the observed values.³⁸ Exhilaratingly, our proposed novel hybrid technique provided a highly accurate forecast for the TB morbidity series as the MAPE showed a value less than 10% both in the simulated and predicted horizons, which further reinforces the usefulness of our proposed new method as a helpful forecasting tool. Also, we used the TB incidence data in Qinghai province to provide an external validation for the predictive ability of our proposed hybrid method,¹² Tables S3-S5 and Figures S15-S16 exhibit the resulting results, suggesting a same finding as our current work that the EEMD-SARIMA-NARNN mixture method produces more accurate approximation to the TB incidence data in Qinghai. Besides, recent work of our research group indicates that the error-trend-seasonal (ETS) model can also accurately estimate and forecast the epidemic patterns of TB incidence,⁴ and thus we further compared the mimic and predicted abilities between our proposed hybrid technique and ETS model. The modeling results displayed that the ETS(A,AD,A) model is supposed to be considered the preferred (Tables S6-S7). Similarly, visible in Table S8, the EEMD-SARIMA-NARNN hybrid technique also significantly outperforms the optimal ETS model. Therefore, this novel data-driven hybrid technique can be recommended as an effective and valuable tool in analyzing and estimating the temporal trends of the TB incidence in Tibet. Moreover, although the SARIMA model, NARNN model, ETS model, and SARIMA-NARNN mixture model produced a relatively low predictive performance relative to our proposed hybrid technique, they are also shown to have a good potential to model the TB epidemics because they gave a MAPE value less than 20% both in the mimic and predictive aspects, this is in good agreement with prior researches that used the above common models to analyze the temporal patterns of other infectious diseases (eg COVID-19,³⁹ typhoid and paratyphoid fevers,⁴⁰ schistosomiasis,⁴¹ and hand-foot-mouth disease⁴²). Also, our proposed hybrid model can be used to estimate the current intervention effects for TB, if this model estimates a significantly higher incidence than the actual, meaning that the current interventions play an important role. Otherwise, additional or optimized measures are essential, such as optimization of the current tools, the pursuit of universal health coverage and social protection, the introduction of new vaccine and new prophylaxis, and the increase of the funding for TB.²

The SARIMA model has been recommended as a versatile tool for estimating and forecasting time series with marked seasonality and periodicity, and for a non-stationary time series, after differencing or making a transformation with logarithm or square root, this model can also be used and can achieve satisfactory forecasting. For instance, Chen et al developed a SARIMA method to estimate the epidemiological trends of the influenza incidence, they found that the SARIMA(0,1,0) (0,1,2)₁₂ method showed a MAPE value of 10.597% for predicting the influenza incidence in the urban area of Shenyang, and the SARIMA(1,1,1) (1,1,0)₁₂ method had a MAPE value of 19.414% in the rural area.⁴³ Gao et al observed that the SARIMA(0,1,7)(1,0,1)₁₂ method displayed a good performance for forecasting the cumulative incidences of typhoid (MAPE=13.257%) and paratyphoid fevers (MAPE=19.501%).⁴⁰ Although the SARIMA model displays a good forecasting performance, it cannot describe the nonlinear information included in the TB morbidity in that this model assumes that there is a linear link between successive values of the time series.⁴ However, different from the SARIMA model, the dynamic NARNN method has been deemed as a promising alternative to address any nonlinear issue without any constraints because of its short-run memory function besides the common properties of static BPNN and GRNN.^9,44,45 Thereby, researchers developed the traditional SARIMA-NARNN mixture model to take advantage of each method’s potentials, and yet recent studies found that such a conventional combined mode may fail to model seasonal or trend variations effectively based on the unpreprocessed time series as the conventional model may not capture the linear and nonlinear components reasonably well,⁴⁶ namely, the conventional hybrid models do not necessarily overmatch their individual predictions.^19,46 Also, as evidenced by our and other work, the traditional SARIMA-NARNN or other conventional mixture models (such as SARIMA-BPNN and SARIMA-GRNN) cannot offer a performance improvement over the basic NARNN and ETS methods on the testing data.^19,47–49 For this reason, our proposed data-driven hybrid model is more reasonably customized for forecasting the different components via decomposing the target series into different linear and nonlinear traits based on the EEMD method. As expected, this hybrid technique is shown to be mechanically more robust and accurate for forecasting the TB incidence in Tibet compared with the others. Given its superiority and the low-cost data collection, this data-driven EEMD-SARIMA-NARNN combined technique seems to be transferable to analyze the epidemiological trends of the TB incidence in other regions or other types of infectious diseases, even of all-time series data, whereas much work is still required to investigate its efficiency and usefulness. Besides, of note, some new advanced deep mining methods (such as long short-term memory neural network (LSTM) and convolutional neural network (CNN)) have currently been applied to predict the dengue cases,⁵⁰ daily land surface temperature,²⁸ and ICD-9 codes of medical records,⁵¹ and these studies suggested that the LSTM and CNN also showed high accuracy for forecasting the target series. Therefore, further work is supposed to compare the forecasting ability between our proposed combined method and the above deep mining methods.

In this work, we discovered that the TB morbidity exhibited marked seasonal variations, predominantly peaking in late winter and early spring and a trough in autumn and early winter. Such a seasonal pattern is similar to previous findings from across China and most of the southwest regions of China (such as Qinghai, Xizang, and Yunnan),^5,12,52 along with the results obtained in Japan and Spain.^53,54 A preceding review to examine the seasonal variability based on 12 studies also supports further evidence that the peak patterns of TB are predominantly seen during the spring and summer seasons.⁵⁵ And yet in disagreement with that in Eastern Cape, the United States, Korea, Mongolia, northern India, and Kuwait.^9,56–58 In this study, an interesting finding is that the TB incidence series in Tibet displayed a primarily exclusive semi-annual seasonal pattern with longer high-risk seasons from January until August and low-risk seasons from September until December per year, which is inconsistent with those observed in earlier studies.^5,12,52,55 Seemingly, various factors are related to the mechanisms underlying TB high-risk seasonal behavior. Firstly, Tibet, as a famous tourist city, attracts hundreds of thousands of persons every year, with tourists being over 40 million in 2019.⁵⁹ Among them, the persons touring during the TB high-risk months approximately accounting for most of them. Large-scale population movements may drive the transmission of TB,⁶⁰ which is also in keeping with the temporal patterns present in the high-risk months in the spatial-temporal sequence. Secondly, current studies have shown that meteorological parameters may play a pivotal role in the spread of TB.^1,61 A study found that the TB incidence showed a decreased trend with the increases in temperature and precipitation in Qinghai.⁶² Another study suggested that wind speed is positively correlated with the TB incidence in Jinghong city of Yunnan.⁶¹ The climate of Tibet is characterized by low pressure, strong wind, heavy rainfall, and relatively low temperature (annual mean temperature is 8°C), this climate type may increase the likelihood of contracting TB in the high-risk seasons as people tend to move indoors or crowded spaces for carrying out their activities and work.⁹ Thirdly, it is thought that vitamin D levels may be closely associated with the seasonal variation of TB, this is mediated through the elevated risk of activation of latent TB because of late winter nadirs in vitamin D levels.^52,56 The reduced sunshine in winter causes vitamin D deficiency, which further influences the potential of macrophages to kill the bacillus Mycobacterium tuberculosis in vitro.^56,63 Fourthly, Lunar New Year, the foremost festival in China, may also play an important role in the seasonal peak patterns of TB since this festival is followed by the largest scale collective on the Tibet migrates, particularly in January and February per year. Given that there is an incubation of around 4–8 weeks from infection to medical diagnosis,⁶⁴ this is wholly in agreement with the maximum peak in March. Besides, it can be seen that there are lower reported cases in February than in January and March. This is derived from the fact that people do not tend to seek medical advice during the Lunar New Year period. Finally, in recent years, more and more activities involving the knowledge-attitude-behavior and new control modes of TB were implemented on World Tuberculosis Day in Tibet, which can enhance public awareness of TB and they are more willing to access healthcare. Besides, an investigation into other factors driving the seasonal variation should go on.

In contrast with the fact that the TB morbidity is falling per year at the rates of about 1.6% and 3% since 2005 in the world and China, respectively,² the TB incidence was demonstrated to display a rapidly increasing trend in Tibet, with AAPC=5.8 from 2006 to 2017. To achieve the 2035 target of the End TB Strategy for 90% reductions in the TB incidence compared with 2015, the WHO estimates that the downtrend in the TB incidence requires to be accelerated at the rates of around 5% by 2020, 10% by 2025, and 17% between 2025 and 2035 globally.^2,65 The current epidemic patterns of TB in Tibet has been an extremely daunting challenge for reaching the aim of a world free of TB during different periods.⁶⁵ Many reasons may be responsible for such an increasing trend in the TB epidemics in Tibet. Firstly, TB is an infectious disease of poverty; Tibet has poor essential health-care services and low socio-economic conditions, with the GDP ranking the last among 31 provinces and municipalities of China in 2019,⁵⁹ which may be the most important trigger for the rapid rise in the TB incidence. Secondly, local meteorological and environmental factors may also be important contributors. Thirdly, Buddhist monks and nuns that account for a significant proportion of the population in Tibet have poor health-care knowledge and they frequently undertake their religious activities in relatively small temples, irrespective of the time of the year, which facilitates the continued transmission among people.⁶⁶ Fourthly, in recent years, the laboratory equipment and testing ability are continuously improved, and new techniques and methods of molecular biology detection and rapid diagnosis are popularized to improve the level of patient discovery in Tibet. Fifthly, around 80% of TB patients in Tibet had been resident in farming and pastoral areas,⁶⁶ and many TB patients do not know that they can be examined and treated free of charge. Faced with the catastrophic costs for TB patients and their households, they often choose not to access healthcare, which further helps the spreading of TB among people. Finally, a study showed that the public awareness of core information on TB prevention and treatment was only 34%, significantly lower than the national average level (57.0%),⁶⁷ such a low public awareness may be a contributory factor in relation to the elevated TB incidence in Tibet. Besides, some other reasons (such as alcohol use disorders, diabetes, HIV infection, smoking, and undernourishment) may also play a role. Confronted with such a worsening tendency, some additional or comprehensive intervention strategies must be taken to combat TB in Tibet.

The present study was subject to some potential weaknesses. Firstly, the time series of the TB incidence was derived from a passive monitoring system, and therefore the underreporting may be inevitable. Secondly, we only used a dataset spanning 12 years to construct the models and some detailed information (such as sex and age) of the TB cases cannot be extracted, which limits further sensitivity analysis for our model. Thirdly, there is currently still in the absence of a guidance to determine the suitable key parameters of ANNs (such as hidden units and feedback delays). In applications, a series of experiments need to be performed by trial and error until the best performance indices are detected. Fourthly, incorporating potential influencing factors of TB into the EEMD-SARIMA-NARNN method may favor the improvement of the forecasting performance, and yet this failed to be done in our work. Fifthly, whether this novel mixture method can also achieve a high-accuracy prediction for other infectious diseases or time series, further research should be undertaken to explore its potential use. Finally, in applications, the newly reported TB cases should be added into the model in time to make sure of its forecasting accuracy level.

Conclusions

Our proposed novel data-driven EEMD-SARIMA-NARNN hybrid method can better capture the dynamic dependence characteristics of the TB incidence in Tibet compared with the basic SARIMA model, basic NARNN model, ETS model, and traditional SARIMA-NARNN combined method, which can be recommended as a helpful and instrumental policy-supportive tool for planning the health-care resources ahead of time and directing the intensity and type of current interventions required to end the TB epidemic. Besides, considering the ongoing increases in TB cases, some additional or comprehensive intervention strategies must be taken to combat TB in Tibet.

Data Sharing Statement

All the data supporting the findings of the work are contained within the study.

Author Contributions

All authors made substantial contributions to conception and design, acquisition of data, or analysis and interpretation of data; took part in drafting the article or revising it critically for important intellectual content; agreed to submit to the current journal; gave final approval of the version to be published; and agree to be accountable for all aspects of the work.

Acknowledgments

We appreciate the Chinese CDC for providing the monthly number of TB cases. This project was supported by the Innovation and Entrepreneurship Training Project for University Students of Henan Province (S202010472007) and Xinxiang Medical University (XYXSKYZ201932) and the Key Scientific Research Project of Universities in Henan (21A330004).

Disclosure

The authors report no conflicts of interest in this work.