Forecasting and analyzing influenza activity in Hebei Province, China, using a CNN-LSTM hybrid model

Abstract

Background

Influenza, an acute infectious respiratory disease, presents a significant global health challenge. Accurate prediction of influenza activity is crucial for reducing its impact. Therefore, this study seeks to develop a hybrid Convolution Neural Network—Long Short Term Memory neural network (CNN-LSTM) model to forecast the percentage of influenza-like-illness (ILI) rate in Hebei Province, China. The aim is to provide more precise guidance for influenza prevention and control measures.

Methods

Using ILI% data from 28 national sentinel hospitals in the Hebei Province, spanning from 2010 to 2022, we employed the Python deep learning framework PyTorch to develop the CNN-LSTM model. Additionally, we utilized R and Python to develop four other models commonly used for predicting infectious diseases. After constructing the models, we employed these models to make retrospective predictions, and compared each model’s prediction performance using mean absolute error (MAE), root mean square error (RMSE), mean absolute percentage error (MAPE), and other evaluation metrics.

Results

Based on historical ILI% data from 28 national sentinel hospitals in Hebei Province, the Seasonal Auto-Regressive Indagate Moving Average (SARIMA), Extreme Gradient Boosting (XGBoost), Convolution Neural Network (CNN), Long Short Term Memory neural network (LSTM) models were constructed. On the testing set, all models effectively predicted the ILI% trends. Subsequently, these models were used to forecast over different time spans. Across various forecasting periods, the CNN-LSTM model demonstrated the best predictive performance, followed by the XGBoost model, LSTM model, CNN model, and SARIMA model, which exhibited the least favorable performance.

Conclusion

The hybrid CNN-LSTM model had better prediction performances than the SARIMA model, CNN model, LSTM model, and XGBoost model. This hybrid model could provide more accurate influenza activity projections in the Hebei Province.

Background

Influenza (also referred to as the flu) is an acute infectious respiratory disease. It is caused by two major influenza virus types, flu A and flu B. It is highly contagious, with a short incubation period and severe symptoms [1], and has a history of causing global pandemics, such as the Spanish flu in 1918, the Asian flu in 1957, and the swine flu in 2009 [2]. According to the World Health Organization (WHO), there are approximately one billion cases of flu reported worldwide each year. Among these cases, 3 to 5 million are classified as severe, with 290,000 to 650,000 resulting in fatalities. Thus, flu seriously endangers individual and public health and is responsible for a great deal of economic burdens [3]. Flu typically manifests as an annual seasonal epidemic or a sporadic pandemic [4, 5], and understanding its activity patterns is crucial for effectively carrying out prevention and control measures. Currently, the primary measure against the flu epidemic is widespread vaccination. However, due to the variability of the flu virus and unexpected outbreaks, vaccines have effectiveness as low as 20% [6, 7]. As a result, analyzing historical surveillance data and making prediction about flu is of great importance.

Flu was the first infectious to be monitored by the WHO. The indicators monitored include influenza-like-illness (ILI), flu specimen positivity rates, and reported flu cases. The sensitivity of various indicators varies across regions. In southern China, the positivity rate of flu specimens serves as a better reflection of the flu epidemic situation, whereas in the north, the ILI proves to be more sensitive [8]. Conversely, the reported number of flu cases often lags significantly behind [9].

Flu surveillance data can be regarded as time series data. Constructing prediction models using previous data can allow for the creation of public health guidance for future epidemics. Traditional prediction models like logistic regression models, Extreme Gradient Boosting (XGBoost) models, and Seasonal Auto-Regressive Indagate Moving Average (SARIMA) models, have been utilized for predicting various contagious diseases such as hand, foot, and mouth disease [10], tuberculosis [11], and COVID-19 [12], and have demonstrated good prediction performance. Additionally, various deep learning methods including the Recursive Neural Network (RNN), the Convolution Neural Network (CNN), and the RNN-based improved Long Short Term Memory neural network (LSTM) have been applied across numerous domains. For instance, CNNs are used to detect anomalies in circuit signals [13], while LSTMs are utilized for predicting stock prices [14], weather patterns [15], and the number of cases [16]. In comparison to traditional prediction models, these methods have stronger nonlinear mapping and sequence information extraction capabilities [17]. In addition, different neural networks possess varying information processing capabilities; CNNs excel in information extraction, while LSTMs are more suitable for processing long sequence prediction tasks. Combining the strengths of these different models may enhance prediction performance.

Hebei Province, which is situated in northern China, typically experiences concentrated influenza activity during the winter and spring seasons [18]. In recent years, the number of reported flu cases during the flu season in Hebei Province has been increasing [19]. It is crucial to utilize surveillance data for predicting peak outbreaks in advance to effectively prevent and control flu in the region. However, there are few studies on flu prediction models, with most focusing on several sentinel hospitals in urban areas. Furthermore, these studies predominantly rely on traditional prediction models such as ARIMA and multiple regression models [20, 21]. However, compared to machine learning and deep learning models in other studies, their predictive results are not very good. Thus, there is a significant gap in the development of accurate and representative flu prediction models specifically tailored for Hebei Province. Here, we gathered ILI data reported by 28 sentinel hospitals across 11 cities in Hebei Province. We then developed a hybrid CNN-LSTM model (which leveraged the distinct advantages of both the CNN and LSTM architectures) on the PyTorch platform. Finally, due to the SARIMA and XGBoost models are frequently employed in time series forecasting and characterized by high accuracy. We chose the SARIMA model and the XGBoost model as the baseline models, and compared the predictive efficacy of our hybrid model to the SARIMA, XGBoost, CNN, and LSTM models using retrospective forecasting.

Methods

Figure 1 is a diagrammatic representation of our research methodology, including the materials and methods used to develop the SARIMA, XGBoost, CNN, LSTM, and CNN-LSTM models Forecasting Model for ILI fluctuation in the Hebei Province of China.

Fig. 1

Study methodology

Data acquisition and splitting

Data were retrieved from the China Influenza Surveillance Information System, which gathers weekly reported ILI cases and related outpatient/emergency department visits from 28 national sentinel hospitals in the Hebei Province. Our data covered the period from the first week of 2010 to the 52nd week of 2022. Due to variations in the sizes of sentinel hospitals and the population bases of different areas, calculating the ILI percentage helps to minimize these factors [22]. ILI% were calculated based on ILI cases and visit data using the following formula:

$$ILI\%=\frac{number of ILIs}{total number of outpatient/emergency visits}\times 100\%.$$

We then split the dataset using an 80–20 ratio, with the initial 80% used as a training set and the remaining 20% used as a testing set for predicted results.

Traditional models

SARIMA model

SARIMA models consist of four components: autoregression (AR), integration (I), moving average (MA), and seasonal effects (S) [23]. SARIMAs can be classified as either multiplicative or additive models, depending on the relationship between seasonal effects and other components. In a multiplicative model, the expression is SARIMA $(p,d,q){(P,D,Q)}_{s}$, where $p$ and $P$ represent the orders of ordinary AR and seasonal AR, respectively; $d$ and $D$ represent the orders and the steps and orders of the difference function; $q$ and $Q$ represent the ordinary orders of MA and seasonal MA individually; and $S$ represents the seasonal period. The SARIMA calculation formula was as follows [24]:

$$\begin{array}{c}\phi_p(B)\Phi_P(B^S){(1-B)}^d{(1-B^S)}^Dy_t=\theta_q(B)\Theta_Q(B^S)u_t,\\\phi(B)=1-\phi_1(B)-\phi_2{(B)}^2-......-\phi_p{(B)}^p,\\\Phi\left(B^S\right)=1-\Phi_1\left(B\right)^S-\Phi_2\left(B\right)^{2S}-......-\Phi_P\left(B\right)^{PS},\\\theta(B)=1-\theta(B)-\theta_2{(B)}^2-......-\theta_q{(B)}^q,\\\Theta(B^S)=1-\Theta{(B)}^S-\Theta_2\left(B\right)^{2S}-......-\Theta_Q\left(B\right)^{QS},\\E\left(u_t\right)=0,Var\left(u_t\right)=\sigma^2,\\E\left({\mathrm u}_{\mathrm t},{\mathrm u}_{\mathrm S}\right)=0,S\neq t,E\left(y_s,u_t\right)=0,\forall_S< t_{\circ} \\\end{array}$$

In this formula, $B$ is the lagging operator, $\phi$ and $\theta$ are the ordinary AR and MA coefficients, $\Phi$ and $\Theta$ are the seasonal AR and MA coefficients, $s$ is the seasonal period, ${y}_{t}$ is a zero-mean steady random time series, and ${u}_{t}$ is a white noise sequence.

The SARIMA modeling process involves several key steps includes time series preprocessing, stationarity testing, parameter identification, and model fitting [25]. The initial data processing step includes conducting an augmented Dickey-Fuller (ADF) test and a white noise sequence test (Box-Ljung test) to determine whether the data represent a stationary and non-white noise sequence. If the data are non-stationary, differencing is necessary to achieve stationarity. The values of $d$ and $D$ are defined by the number and order of the differencing steps.

Next is the model identification phase, which involves examining the autocorrelation function (ACF) and partial autocorrelation function (PACF). Within these plots, the AR can be identified in a significant autocorrelation at lag ($p$) in the PACF, and the MA can be identified in significant autocorrelation at lag ($q$) in the ACF. The parameters $P$ and $Q$ correspond to the determination of $p$ and $q$. Finally, seasonal patterns are recognized as periodic spikes which appear at multiples of the seasonal lag ($S$). These steps help systematically determine the appropriate parameters for constructing SARIMA models.

XGboost model

XGBoost is a boosting algorithm proposed by Chen Tianqi et al. in 2016 [26]. Based on the gradient boosting framework, it sequentially trains multiple weak classifiers, minimizes the loss function through gradient descent, and continuously enhances the model’s predictive ability. Figure 2 illustrates the basic structure of the XGBoost model. Initially, input features are constructed using the hysteresis method, and a simple prediction model is initialized. Subsequently, the model calculates the error between its predictions and the training set, adjusts parameters based on this error, and constructs a new model. This iterative process constructs multiple models sequentially, gradually improving the performance of a single model. Finally, all models are integrated and weighted according to their contribution to the overall prediction effect, resulting in the final XGBoost model.

Fig. 2

Architecture of the extreme gradient boosting

Deep learning models

Data pre-processing

In neural network models, training process, large differences in data scales can lead to excessively large or small gradient updates. Additionally, outliers in the data may cause certain features to dominate others. Normalization ensures that all features are on a similar scale, reduces the impact of outliers, and enhances the model’s convergence speed and generalization ability [27]. Here, original data were normalized with the formula provided below:

$$X=\frac{x-x_{min}}{x_{max}-x_{min}}$$

where the $X$ is the normalized value of the data, the $x$ is the original data, and ${x}_{max}$ and ${x}_{min}$ are their maximum and minimum values, respectively.

Furthermore, the input data for these models needed to contain both inputs (X) and outputs (Y), and needed to be reconstructed into 3D tensors format (samples, timesteps, features). Given that the ILI% data were a one-dimensional dataset of length n, a sliding-window approach was employed (Fig. 3). Observations of the current k time points consisted of the input data X, and observations of the current time points consisted of the output data Y. We then used the Numpy library in Python to convert data into an (n, n-k, 1) format.

Fig. 3

Sliding window method to generate inputs (X) and outputs (Y)

CNN model

CNN models have superior capabilities in extracting local features [28]. They can be employed in computer vision and time series prediction and analysis [29]. CNN models consist of convolutional layers, pooling layers, activation functions, and dense layers (Fig. 4). For our time series data, 1D kernels were used to slide across the data (x) and compute dot products between the kernel and local data at each step, eventually generating feature maps. We then used max pooling on feature maps to reduce dimensions while retaining important features. After adding non-linearity to pooled data using a tanh function, data were input into the dense layer and produced a CNN model using the following formula:

$$\begin{array}{c}X=(x\ast w)+b_i\\X_c=max\left(X_{i:i+k-1}\right)\\y=W\ast X_c+b_j\end{array}$$

where $(x * w)$ denotes the convolution operation, $b$ is a bias term ${X}_{i:i+k-1}$ is a set of consecutive elements in the input feature sequence, $max( )$ is a maximum extraction operation, $W$ is a weight matrix, and $*$ is a matrix multiplication operation.

Fig. 4

Architecture of the convolutional neural network

LSTM model

LSTM models are improved neural networks based on RNN [30]. They have a similar structure to RNN, but enhanced memory units (Figs. 5 and 6) [31]. By introducing gate control mechanisms, LSTM models address the long-period dependency problems and the gradient vanishing and exploding issues that are present in RNN models. Our LSTM memory units consisted of input (I), forget (F), and output (O) gates, where the F gate decided the discording and preserving of information, the I gate decided what new information to store in the cell state, and the O gate decided what the next hidden state should be. The calculated LSTM formula was as follows [32]:

$$\begin{aligned}\left\{\begin{array}{ll}cf_t&=\sigma\left(W_f\cdot\left[h_{t-1},x_t\right]+b_f\right)\\i_t&=\sigma\left(W_i\cdot\left[h_{t-1},x_t\right]+b_i\right)\\{\widetilde C}_t&=tan h\left(W_C\cdot\left[h_{t-1},x_t\right]+b_C\right)\\C_t&=f_t\ast C_{t-1}+i_t\ast{\widetilde C}_t\\o_t&=\sigma\left(W_o\cdot\left[h_{t-1},x_t\right]+b_o\right)\\h_t&=o_t\ast tan h\left(C_t\right)\end{array}\right.\end{aligned}$$

where ${W}_{C}$ is the weight matrix, and $b$ is the bias term.

Fig. 5

Neural units of the recursive neural network

Fig. 6

Memory cells of the long-short term memory neural network

The F gate takes the previous hidden state ${h}_{t-1}$ and current input ${{\varvec{x}}}_{{\varvec{t}}}$ as input, passes a sigmoid activation function, and outputs${f}_{t}$. In the I gate, the σ layer outputs ${i}_{t}$ and the tanh layer creates a vector of new candidate values${\widetilde{{\varvec{C}}}}_{t}$, and the cell state ${C}_{t}$ vector is updated by combining${f}_{t}$,${i}_{t}$, and${\widetilde{{\varvec{C}}}}_{t}$. Finally, in the O gate, It filters the input information by passing it through the σ layer with the output ${o}_{t}$,and${o}_{t}$, along with the ${C}_{t}$ that has been operated on by tanh layer, completes the update of${h}_{t+1}$.

CNN-LSTM model

Our novel hybrid CNN-LSTM model combined the local feature extraction capability of CNN models and the long-term dependency resolution of LSTM models [33, 34]. Figure 7 visualizes the architecture of our CNN-LSTM model, which was composed of convolutional layers, pooling layers, LSTM layers, the tanh activation function, and a dense layer. The convolutional layers were encoders that used 1D kernels to extract features from the original data and convert them into feature maps. The pooling layer selected and reduced the dimensions of the feature maps, which were then decoded by the LSTM layers and fed into the dense layer to get the model output. This process generated the following formula:

Fig. 7

Hybrid model framework

Model construction and validation

The SARIMA model was constructed using the “forecast” package based on R 4.2.3. The XGBoost model was constructed using the “xgboost” and “sklearn” libraries in Python 3.9. Additionally, the CNN, LSTM, and CNN-LSTM models were developed using the PyTorch 2.2.1 deep learning framework based on Python 3.9.

After the models were constructed, they need to be validated and their parameters determined. For the SARIMA model, the prediction performance varies based on different parameter combinations. By analyzing the Autocorrelation Function (ACF) plot and the Partial Autocorrelation Function (PACF) plot, appropriate values of p, q, P, and Q can be identified. Typically, these values do not exceed 2. Therefore, we can construct multiple alternative models using an exhaustive method to find the best fit. Subsequently, the residuals of the model undergo a Box-Ljung test to identify outliers and assess the distribution. Finally, the model with the lowest Akaike Information Criterion (AIC) value is selected as the optimal model.

For models such as XGBoost, CNN, LSTM, and hybrid CNN-LSTM models, their performance is also influenced by hyperparameters. During the model training phase, the first step involves determining the range of hyperparameters using the training set. Next, the training set is randomly divided into k subsets using the K-fold cross-validation method, where k-1 subsets are used for training and the remaining subset is used for validation. The grid search method iterates through various hyperparameter combinations, training a model on each subset, and calculating the average loss across all k subsets. Finally, the hyperparameter set that yields the smallest average loss is chosen to build the best-performing model.

Comparison of model predictive effects

After the models are constructed and validated, the testing set is used to make retrospective predictions. The mean absolute error (MAE), root mean squared error (RMSE), mean absolute percentage error (MAPE), coefficient of determination (${R}^{2}$), correlation coefficient (r), and explained variance (Evar) were calculated to evaluate the predictive performance of each model against the true values. The formulas used for calculation were as follows:

$$\begin{array}{c}MAE=\frac1N\sum_{i=1}^N\left|x_i-{\widehat x}_i\right|\\RMSE=\sqrt{\frac1N\sum_{i=1}^N\left(x_i-{\widehat x}_i\right)^2}\\MAPE=\frac{100\%}N\sum_{i=1}^N\left|\frac{x_i-{\widehat x}_i}{x_i}\right|\\R^2=1-\frac{\sum_{i=1}^N{(x_i-{\widehat x}_i)}^2}{\sum_{i=1}^N{(x_i-\overline x)}^2}\\r=\frac{\sum_{i=1}^N{((x_i-\overline x)-(y_i-\overline y))}^2}{\sqrt{\sum_{i=1}^N{(x_i-\overline x)}^2}\sqrt{\sum_{i=1}^N{(y_i-\overline y)}^2}}\\Evar=1-\frac{\sum_{i=1}^N\left(\left(x_i-{\widehat x}_i\right)-E\left(\overline x-{\widehat x}_i\right)\right)^2}{\sum_{i=1}^N\left(x_i-\overline x\right)^2}\end{array}$$

where ${x}_{i}$ represents the true value at time $i$, ${\widehat{x}}_{i}$ is the predicted value of the model at time $i$, $N$ is the sample size at the time of testing, and $\overline{x }$ is the mean value of the sequence $x$.

Results

Temporal distribution of ILI%

The fluctuation in ILI and ILI% in Hebei Province between 2010 and 2022 are demonstrated in Figs. 8 and 9. The changes in both are essentially identical. From 2010 to 2019, the seasonal distribution of ILI% showed a unimodal pattern, with the peak of the epidemic in the winter and spring of each year, starting from the 40th week of the first year and continuing to the 15th week of the following year. Between 2020 and 2022, in addition to the peak of the epidemic in winter and spring, there was a small ILI% peak in the spring and summer. Between 2017 and 2019, the ILI% prevalence peak increased year over year, returned to fluctuating levels between 2020 and 2021, and reached a new peak at the end of 2022.

Fig. 8

ILI and ILI% activity from 2010 to 2022

Fig. 9

Annual changes in ILI and ILI% from 2010 to 2022

Data partitioning

Our 2010 to 2022 data spanned 676 weeks. We thus segmented it in the 21st week of 2020 to represent an 80–20 split. The initial 80% of the data, or 540 weeks, was used as the training set, and the remaining 20%, or 136 weeks, was used as the test set. Results are shown in Fig. 10.

Fig. 10

Data division into training and testing sets

SARIMA model

Previous research suggests that ILI has a noticeable seasonality Fig. 7. Specifically, ILI exhibits a seasonal cycle of 52. ADF testing of the training dataset revealed a lack of sequence stationarity. However, after applying first-order seasonal differencing, the sequence became stationary and passed both the ADF and white noise tests. Thus, data were successfully transformed into a stationary white noise sequence suitable for the development of the SARIMA model (Table 1).

Table 1 Sequence stationary and white noise test results

After determining the values of $d$ and $D$ through the previously-mentioned differencing operation, we initially determined the values of $p$, $P$, $q$, and $Q$ by examining the ACF and PACF plots of the differenced series (Fig. 11). Various alternative models were constructed by combining different parameters (Table 2). The best model of those that passed the parameter tests, selected following the principle of minimizing AIC, was the SARIMA (2, 0, 2) (1, 1, 1)₅₂ model. Residual testing showed that residuals at different lag orders exhibited white noise sequences, indicating effective use of sequence information (Fig. 12).

Fig. 11

Sequence’ autocorrelation function and partial autocorrelation function

Table 2 Alternative model parameters and test results

Fig. 12

SARIMA model residual test results

XGBoost model

The construction process of an XGBoost model can be divided into three parts: feature engineering, model construction, and model validation. Since the ILI% data is one-dimensional, we use the current value as the output and construct input features by lagging one seasonal period (52 weeks) before building the XGBoost model. During the model building phase, we utilize the XGBRegressor class to establish a gradient boosted tree regression model with mean squared error as the loss function.

The alternative hyperparameters for the XGBoost model include the maximum depth of the tree (3, 5, or 7), the learning rate (0.01, 0.1, or 0.2), the number of decision trees (100, 200, or 300), and the proportion of subsamples (0.7, 0.8, or 0.9). Through fivefold cross-validation using the grid search method provided by the sklearn library, we determined the optimal combination of hyperparameters: maximum tree depth of 7, learning rate of 0.2, 200 trees, and subsample proportion of 0.7, yielding a minimum average error of 0.1173. The results are shown in Fig. 13. Finally, we built an XGBoost model based on the best hyperparameters described above.

Fig. 13

XGBoost’s losses for different combinations of hyperparameters

CNN, LSTM, and CNN-LSTM models

The ILI% period of 52 led us to select a time step of 52 for data reconstruction using the sliding window method. Following reconstruction, the training set had a data shape of (488, 52, 1) and the test set had a data shape of (84, 1).

The development process of neural network models involved into three parts: model definition, validation, and training. During the models definition phase, we defined the model structures and forward propagation modules based on the PyTorch platform. Specifically, the CNN model includes convolutional layer and pooling layer, dropout layer, and fully connected layer with the ReLU activation function. The LSTM model consists of LSTM layer, dropout layer, and fully connected layer. The CNN-LSTM model combined aspects from both the CNN and LSTM models by incorporating convolutional layer, LSTM layer, dropout layer, and fully connected layer. In the forward propagation modules: For the CNN model, the tensor input passes through the convolutional and pooling layers, followed by activation functions and dropout layers, culminating in the fully connected layer to produce the output. The LSTM model generates output through its LSTM layer, dropout layer, and fully connected layer. The CNN-LSTM model processes the feature map through a convolutional layer, LSTM layer, dropout layer, and fully connected layer to obtain the final output.

In the model validation phase, alternative hyperparameters for the CNN model include the number of convolutional layers (1 or 2), convolution kernel size (3 or 5), number of pooling layers (1 or 2), learning rate (0.01 or 0.001), and batch size (16, 32, or 64). For the LSTM model and CNN-LSTM model, hyperparameters include the number of LSTM layers (1, 2, or 3), number of hidden units (16, 32, or 64), learning rate (0.01 or 0.001), and batch size (32, 64, or 128). After fivefold cross-validation using the grid search method, optimal hyperparameter combinations were determined: The CNN model consisted of 2 convolutional layers with a kernel size of 5 and 1 pooling layer, using a learning rate of 0.001 and a batch size of 16, achieving a minimum average loss of 0.0038. The LSTM model included 2 LSTM layers with 64 hidden units, a learning rate of 0.01, and a batch size of 128, achieving a minimum average loss of 0.0031. Additionally, the CNN-LSTM model featured 1 LSTM layer with 16 hidden units, a learning rate of 0.001, and a batch size of 32, achieving a minimum average loss of 0.0021. These results are illustrated in Fig. 14. Finally, three models were constructed using the best combination of hyperparameters.

Fig. 14

Losses of distinct hyperparameter combinations in three models

In the model training phase, which includes backward propagation and iterative training. MSE served as the loss assessment metric and Adam as the optimizer. After 200 training epochs, the losses converged: CNN model loss was 0.0005, LSTM model loss was 0.0010, and CNN-LSTM model loss was 0.0014. The training loss of three models are shown in Fig. 15.

Fig. 15

Change in training loss over 200 epochs of the three models

Comparison and analysis of models

After constructing the models, we utilized the aforementioned five models to make retrospective predictions (Fig. 16) and compared them with the test set to calculate MAE, RMSE, MAPE, R², r, and Evar. The SARIMA model was used to forecast ILI% over a span of 136 weeks starting from week 21 in 2020. The XGBoost model predicts 84 weeks based on the constructed input features (starting from week 21 in 2021). And the neural network models were used to generate forecasts for the subsequent 84 weeks (starting from week 21 in 2021) via test inputs created by the sliding-window method. All the models demonstrated improved prediction accuracy with regard to ILI% tendencies. We evaluated the models’ predictive performance using MAE, RMSE, MAPE, R², r, and Evar, with specific results detailed in Tables 3, 4, 5 and 6. On the testing set, the SARIMA model had the highest MAE, RMSE, and MAPE, indicating relatively poor prediction performance in terms of R², r, and Evar. Conversely, the CNN, LSTM, XGBoost, and CNN-LSTM models showed smaller prediction errors and higher R², r, and Evar. Among these models, prediction errors ranked from largest to smallest as CNN, LSTM, XGBoost, and CNN-LSTM, while R², r, and Evar followed as CNN, XGBoost, LSTM, and CNN-LSTM, respectively. Overall, CNNs demonstrated the weakest forecasting performance, whereas CNN-LSTM showed the most accurate predictions. Across the entire testing set, compared to the SARIMA model, the CNN-LSTM model exhibited reductions in MAE, RMSE, and MAPE by 0.4388, 0.3548, and 27.2063%, respectively, resulting in corresponding accuracy improvements of 49.21%, 25.47%, and 65%. When compared with the XGBoost model, these reductions were 0.0300, 0.0319, and 0.5149%, respectively, with corresponding improvements of 6.21%, 2.98%, and 3%, respectively. The specific results are shown in Table 3.

Fig. 16

Comparison of forecasting results from all four five models

Table 3 Mean Absolute Errors across different models

Table 4 Root Mean Squared Errors across different models

Table 5 Mean Squared Percentage Errors across different models

Table 6 Other metrics across different models

We used five models for predictions over the next 4 weeks, 26 weeks, and 52 weeks, and calculated prediction errors, R², r, and Evar. These specific results of these are shown in Table 7. As prediction length increased, so did prediction errors for all models (Fig. 17). At each prediction interval, the SARIMA model consistently showed the largest prediction error, accompanied by lower R², r, and Evar. In predictions for the next 4 weeks, the XGBoost model exhibited the smallest prediction error, followed by CNN and LSTM, with CNN-LSTM showing the largest error. Over 26 weeks, CNN-LSTM had the smallest error, followed by XGBoost and LSTM, while CNN had the largest error. Looking ahead to 52 weeks, CNN-LSTM maintained the smallest error, followed by XGBoost and LSTM, with CNN demonstrating the largest error. It’s important to note that R², r, and Evar are influenced by sample size, with these metrics generally showing an initial increase followed by a decrease across different prediction periods for the five models. For instance, in predictions over the next 4 weeks, the XGBoost model showed the highest r and Evar values, followed by LSTM, CNN-LSTM, and CNN, while LSTM exhibited the highest R-squared, followed by XGBoost and CNN-LSTM, and CNN showed the lowest values. Over 26 weeks, CNN-LSTM had the highest scores across all three metrics, followed by XGBoost, with CNN and LSTM close behind. Over 52 weeks, LSTM had the highest r, with CNN-LSTM and XGBoost close, while CNN-LSTM showed the highest R² and Evar values, followed by XGBoost and LSTM, with CNN recording the lowest scores.

Fig. 17

Comparison of the evaluation metrics of the five models at different prediction scales

Table 7 Evaluation metrics of five models in different prediction scales

Discussion

ILI surveillance is crucial for flu prevention and control efforts. By tracking ILI cases and using historical data to construct predictive models, deeper understandings of ILI fluctuation patterns can be obtained. These understandings allow for the implementation of proactive measures to reduce the impact of flu. In this study, we utilized historical ILI% data spanning from 2010 to 2023 to develop a CNN-LSTM hybrid model for retrospective forecasting. We examined and compared its predictive with the SARIMA model, XGBoost model, and standalone neural network models. The study revealed that all five models effectively captured the trend in ILI% changes and demonstrated strong predictive capabilities. Upon further analysis across varying forecast horizons, the CNN-LSTM model exhibited the highest forecasting accuracy, followed by the XGBoost model, LSTM model, and CNN model, with the SARIMA model performing the least accurately.

The SARIMA model is a widely used time series forecasting model that extrapolated future outcomes by integrating autocorrelation and moving average error results. Unlike the XGBoost and neural network models, SARIMA models do not rely on input features and may offer superior forecasting ability for stationary time series. However, SARIMA models’ predictive capabilities diminish when forecasting over long-term horizons. In Qiang Mao’s study, SARIMA effectively predicted tuberculosis incidence in China for short-term forecasts but showed declining accuracy with longer prediction horizons [11]. Which have demonstrated SARIMA’s superior accuracy in short-term prediction tasks. In our comparative study, SARIMA demonstrated poorer predictive capability too. In contrast, neural network models have no specific data requirements and have robust nonlinear mapping capabilities, enabling them to leverage more information for prediction [17]. This disparity underscores why the neural network models in our study had superior prediction accuracy compared to our SARIMA model. Research by S. Siami-Namini et al. indicates that SARIMA’s structure, built on data autocorrelation and dependence, suits linear data prediction better than the other models examined [35]. SARIMA’s performance on nonlinear data is inferior to that of machine learning and deep learning models, which aligns with our findings.

The XGBoost model is a highly popular machine learning technique renowned for its exceptional accuracy and robust generalization in sequence prediction. Utilizing an ensemble learning algorithm, XGBoost integrates multiple decision tree models through gradient boosting, continuously optimizing split points and leaf node divisions to minimize the loss function and enhance model performance. In our study, we employed the XGBoost model to forecast the test set across various time intervals, demonstrating strong predictive capabilities. Among the five models evaluated, XGBoost outperformed the LSTM, CNN, and SARIMA models, though it slightly lagged behind the CNN-LSTM hybrid model in accuracy. In research by Zheng-gang Fang et al., the XGBoost model effectively predicted Covid-19 case numbers in the United States and surpassed the ARIMA model in predictive accuracy, corroborating our findings [36]. Conversely, Junling Luo et al. found that the LSTM model exhibited superior accuracy over XGBoost in predicting Covid-19 cases in Canada [37]. In our study, while the LSTM model’s prediction accuracy was slightly lower than XGBoost’s, this difference may stem from our use of historical ILI% data exclusively for prediction. Notably, the CNN-LSTM model’s predictive efficacy surpassed that of XGBoost in our study, suggesting that combining CNN and LSTM techniques holds promise for improving prediction performance.

CNN and LSTM are both deep learning models with distinct applications. CNN is commonly used in image and speech recognition, whereas LSTM excels in analyzing long time series data. Each model has unique advantages: CNN effectively learns sequence dynamics by extracting local patterns, while LSTM addresses long-term dependency issues using gate mechanisms within its memory units. The CNN-LSTM hybrid model combines these strengths and finds extensive use in fields such as finance and yield forecasting. In our study, CNN-LSTM was applied to forecast flu activity, outperforming individual models and ranking highest among the four models compared. In Muhammad, L. J et al.’s research, CNN-LSTM successfully predicted Covid-19 case numbers in South Africa and Botswana, achieving commendable accuracy [33]. Similarly, studies by Zhang, J et al. and Rui Yan et al. utilized CNN-LSTM to forecast air quality in Beijing, China, comparing it with CNN and LSTM models. Their findings consistently demonstrated CNN-LSTM’s superior predictive performance [38, 39]. Notably, Rui Yan et al. also employed CNN-LSTM for multi-period prediction, aligning closely with our study’s results.

limitations

Flu activity is influenced by various factors, including population immunity, conditions of transmission, and viral variability. In the retrospective projections of this study, the prediction accuracy of the five models declined to varying degrees during 2021 and 2022. This decline possibly due to the stringent preventative measures enacted during the COVID-19 pandemic. In Hu CY’s study, many provinces and cities in China adopted preventive measures such as mask-wearing, closure of large public venues, and prohibition of large gatherings to control COVID-19. As COVID-19 spreads similarly to flu, activities related to flu and other respiratory infections were disrupted [40], and COVID-19 has similar clinical symptoms to influenza, and many COVID-19 patients go to the hospital to see an increase in the number of ILIs [41]. These will lead to significant bias in predictions. For future studies, data processing methods like singular spectrum analysis or one-hot encoding should be considered to mitigate the impact of COVID-19 prevention measures [42, 43].

In additional, since meteorological data such as meteorological data in the study area are not publicly available. Our study solely utilized ILI% data from 28 flu sentinel hospitals in Hebei Province. Although optimal hyperparameters were employed in model construction, the dimensions of the input features could potentially impact prediction results. In Athanasiou et al.’s study, multi-source data including monitoring, meteorological data, and Twitter search data were integrated, yielding excellent prediction outcomes [44]. This suggests that incorporating more diverse and relevant data into models can enhance prediction accuracy. Future studies should consider leveraging multi-source data to improve overall model performance.

In infectious disease modeling, both linear and nonlinear relationships often coexist. Prediction methods for infectious diseases include linear models, machine learning models, and deep learning models, the latter of which possess strong nonlinear extraction capabilities. Although our study only utilized the combined advantages of CNN and LSTM to construct the CNN-LSTM model, there may be limitations in its ability to extract linear relationships. In Yiran Wan’s prediction study on hand, foot, and mouth disease in Chongqing, China, the SARIMA-EEMD-LSTM model was employed, demonstrating higher prediction performance compared to single models [45]. This suggests that integrating different types of models could be beneficial for our future work.

Finally, it is crucial to provide a reasonable estimate of uncertainty in forecasting. As mentioned earlier, influenza incidence is influenced by many factors, and the predictive accuracy of the model can vary with changes in these uncertainties during actual predictions. It is essential to assess a reasonable confidence interval to model’s predictive results. In a study by Morris, Michael et al., they utilized various neural network models to predict influenza activity and employed Gaussian distribution to estimate prediction result intervals, enhancing the realism and effectiveness of their predictions [46]. In our study, we did not adequately estimate the prediction uncertainty interval. Future research should address this by employing appropriate methods to estimate uncertainty factors, thereby incorporating confidence intervals into the model’s prediction results to enhance effectiveness.

Conclusions

Seasonal flu epidemics impose a significant disease burden, emphasizing the need for scientifically accurate predictive models to facilitate early intervention. Our study developed a hybrid CNN-LSTM model to forecast the activation of flu in Hebei Province. We compared the predictive performance of the hybrid model against the SARIMA model, XGBoost model, CNN model, and LSTM model. Our analysis revealed that while all models effectively predicted ILI activity trends, the hybrid CNN-LSTM model demonstrated superior prediction accuracy compared to the others. These findings not only enhance the precision of flu prediction in Hebei Province but also demonstrate the efficacy of hybrid models in specific forecasting tasks. By providing timely scientific guidance before the flu season, proactive measures can be planned and implemented for prevention and control in Hebei Province. These findings not only enhance the precision of flu prediction in Hebei Province but also highlight the efficacy o f hybrid model in specific predictive tasks. By providing timely scientific guidance before the flu season, proactive measures can be planned and implemented for prevention and control in Hebei Province.