Weekly ILI data for Hubei Province and its prefecture-level cities from January 2020 to December 2023 were obtained from the Hubei Provincial Center for Disease Control and Prevention. Meteorological and air quality data for Hubei Province covering 2020‐2023, including mean temperature, relative humidity, air quality index, and mean dew point, were acquired from the National Science Data Center platform.

To incorporate the Baidu Search Index (BSI), a literature search was conducted on databases such as PubMed and CNKI to identify studies related to ILI. Frequently used keywords in published papers were summarized, resulting in an initial list of 53 candidate keywords potentially associated with influenza epidemic trends (Table S1 in Multimedia Appendix 1). These keywords included terms related to early ILI symptoms, medications, treatments, and other closely associated phrases, such as “fever,” “body temperature,” and “cough.” The initially screened keywords were subsequently excluded and refiltered (Multimedia Appendix 1). Nine keywords showing significant correlation with influenza activity were finally selected (Table S2 in Multimedia Appendix 1). To avoid multicollinearity, these keywords were aggregated into a synthesized Baidu Index as an internet-based surveillance factor, with higher weights assigned to keywords exhibiting stronger correlations [14] (Multimedia Appendix 1). The composite index was calculated using the following formula:

In Formula (1), \ MW(t) represents the synthesized Baidu Index, Mi denotes the Baidu Index of the i-th keyword, Ci is the maximum correlation between the Baidu Index of the i-th keyword and influenza cases.

To dynamically capture weekly variations in population mobility, we selected the Baidu Migration Scale (BMS) Index as a real-time proxy. Derived from the Baidu Map Qianxi platform using location-service data [19], the BMS provides continuous, weekly insights into regional population aggregation and movement directions. For this study, Hubei Province’s immigration and emigration indices from January 2020 to December 2023 were used to measure population movement.

NPIs and policy measures were characterized using the Oxford COVID-19 Government Response Tracker (OxCGRT), an evaluation framework developed by the Blavatnik School of Government research team to systematically track government policies worldwide during the COVID-19 pandemic [20]. To quantify government response intensity, we used the 2020‐2023 Oxford SI for Hubei Province. Given its phased nature—dropping to 0 after China’s COVID-19 policy shift on January 8, 2023—the SI was incorporated as a dynamic variable. When SI>0, it acts as an active feature; when SI=0 (throughout 2023), its input remains 0, effectively ignoring its contribution. Processing details are in Multimedia Appendix 1.

To rigorously prevent data leakage, the 209-week internal dataset (2020‐2023) was chronologically split prior to any sliding window processing or standardization. Based on an 85:15 ratio, the data were divided into a training set comprising the first 177 weeks (wks 1‐177) and a testing set comprising the remaining 32 weeks (wks 178‐209). Subsequently, a 6-week sliding window was applied independently to each dataset to construct the supervised learning sequences. Because the prediction for a given week t requires features from the preceding weeks t-6 to t-1 , this independent sliding mechanism inherently consumes the first 6 weeks of each split. Consequently, this process yielded 171 training sequences (dimension: 171×6×6, generating predictions from wk 7 to 177) and 26 testing sequences (dimension: 26×6×6, generating predictions from wk 184 to 209). For normalization, z score standardization was applied, with parameters calculated exclusively from the training set to scale the testing data. Finally, K-fold cross-validation was strictly confined to the training set, ensuring the testing set remained completely unseen during model tuning. An additional 26-week dataset from 2024 was later used strictly as an external validation set to evaluate prospective forecasting capability.

In this study, the SARIMA model was built using the Python statsmodels library, while deep learning and RF models were developed using the PyTorch framework and scikit-learn library. For the SARIMA model, optimal configurations were identified through grid search and evaluated based on the lowest AIC and BIC values. For deep learning and machine learning models, hyperparameter tuning was conducted using a 5-fold cross-validation coupled with a grid search strategy to identify the configuration yielding the minimum average loss (mean squared error [MSE]) [26]. To mitigate overfitting in neural networks, a Dropout mechanism (P=.2) was incorporated. During model training, the Adaptive Moment Estimation (Adam) algorithm was used as the optimizer to ensure stable convergence to optimal values [31]. Detailed validation procedures, hyperparameter selection processes, and residual diagnostics are provided in Multimedia Appendix 1.

Following model construction and validation, the test set was used for retrospective prediction. Overall model performance was evaluated using mean absolute error (MAE), root-mean-squared error (RMSE), mean absolute percentage error (MAPE), coefficient of determination (R²), and MSE. Furthermore, to assess the robustness of the models across different forecast horizons, specific evaluations were conducted for short-term (8 wks), medium-term (12 wks), and long-term (26 wks) predictions. These analyses were performed on the same test set by extracting and comparing the predicted values of the final 8-, 12-, and 26-time steps against the corresponding ground truth data. For each specific forecast horizon, performance was comprehensively evaluated using R², RMSE, MAE, and MAPE. Specifically, R² reflects the overall goodness-of-fit, while MAPE indicates the model’s capacity for relative error control; together, these 2 metrics serve as the primary dimensions for evaluating prediction robustness over varying time lengths.

The ILI surveillance data used in this study were provided by the Hubei Provincial Center for Disease Control and Prevention. These data are routine, aggregated weekly statistics from the statutory infectious disease surveillance system and are fully deidentified, containing no personally identifiable information. Because this study relied exclusively on aggregated data and did not involve the collection of individual human participant data, formal institutional review board ethical approval was waived. Furthermore, all meteorological data and internet-derived indices (BSI, Baidu Migration Index, and Oxford SI) were sourced from publicly accessible platforms; these are nonpersonal datasets that do not implicate personal privacy or confidentiality. Finally, as this research did not involve the direct participation of human participants, informed consent procedures were not applicable, and no participant compensation was provided (Figure 1).

From January 2020 to December 2023, ILI in Hubei Province predominantly occurred during the winter-spring transition, with reported cases peaking each year in December and January. The first epidemic peak was observed in early 2020, after which incidence declined rapidly and remained at low levels from February 2020 through the end of 2021. In early 2022, ILI incidence began to rise again, followed by a distinct second peak in July 2022, before decreasing once more. After the relaxation of COVID-19 containment measures in 2023, ILI incidence increased sharply, with the highest peak during the observation period occurring in February 2023-March 2023, when the weekly incidence surged to 143.97 per 100,000 population. Another epidemic peak reappeared in the winter of 2023, reaching a maximum weekly incidence of 106.71 per 100,000. Given the substantial fluctuations in ILI incidence associated with NPIs during the pandemic, a logarithmic transformation was applied to the incidence data in subsequent deep learning modeling. This preprocessing strategy helped mitigate the influence of extreme values on model training and enhanced model stability and generalizability (Figure S1 in Multimedia Appendix 1).

The overall performance of the 4 time-series forecasting models in predicting ILI counts is summarized in the table. The results demonstrate that the deep learning models substantially outperformed the traditional statistical model. Among them, the CNN-LSTM hybrid model achieved the best performance, with an R² of 0.7623, RMSE of 0.3433, MAE of 0.1623, and MAPE of 9.2142%. The RF model exhibited comparable accuracy to CNN-LSTM, yielding an R² of 0.7624 and RMSE of 0.3433, though its MAPE (10.6332%) was slightly higher. In contrast, the SARIMA model performed poorest, with an R² of –0.2645, RMSE as high as 0.8863, MAE of 0.7837, and MAPE of 62.7447% (Table 1).

To assess the stability and accuracy of the models over varying prediction periods, performance was compared across three forecast horizons: 8 weeks, 12 weeks, and 26 weeks (Table 2). For short-term forecasting (8 wks), all deep learning models performed well. CNN-LSTM (R²=0.87, RMSE=0.07, MAE=0.06, and MAPE=5.61%), LSTM (R²=0.87, RMSE=0.07, MAE=0.06, and MAPE=5.47%), and RF (R²=0.88, RMSE=0.07, MAE=0.06, and MAPE=5.83%) achieved similar accuracy, with all error metrics remaining at low levels. When the forecast horizon was extended to 12 weeks, model performance declined slightly but remained stable: CNN-LSTM (R²=0.86), LSTM (R²=0.86), and RF (R²=0.87) continued to perform closely, with MAPE values ranging between 6% and 6.3%. In long-term forecasting (26 wks), the CNN-LSTM model demonstrated stronger robustness, attaining an R² of 0.89, RMSE of 0.11, MAE of 0.08, and MAPE of 7.03%, which was notably better than LSTM (R²=0.91, but MAPE=6.86%) and RF (R²=0.91 and MAPE=7.29%). The SARIMA model performed poorly across all forecast horizons. For the 8-week horizon, it produced an R² of –14.82 and an MAPE as high as 76.58%; in the 26-week horizon, its performance improved slightly (R²=–0.56 and MAPE=85.12%), yet remained substantially inferior to the deep learning models (Table 2).

This study systematically compared the performance of five models—SARIMA, LSTM, CNN-LSTM, Transformer, and RF—in forecasting ILI incidence. The results indicate that deep learning models significantly outperformed the traditional statistical model, and multisource data models exhibited superior predictive accuracy and precision compared with univariate time-series models. In univariate ILI time-series forecasting, CNN-LSTM achieved the best performance, followed by RF and LSTM, while the SARIMA model performed poorest. This finding aligns with the conclusion of Li et al [26], who reported the distinct advantage of deep learning methods in capturing nonlinear patterns in influenza data. Influenza incidence is influenced by multiple factors, including viral mutation, climatic variation, and human behavior, resulting in highly nonlinear and non-stationary characteristics. As noted by Nsoesie et al [32], traditional statistical models have inherent limitations when handling abrupt outbreaks and irregular fluctuations in epidemiological data. In this study, the negative R² value obtained by the SARIMA model suggests that its predictions were even less accurate than a simple mean forecast. The SARIMA model relies on assumptions of linearity and stationarity, which limit its ability to adapt to the sudden outbreaks and nonstationary nature of influenza data. In contrast, deep learning models demonstrated clear advantages. The LSTM architecture, through its gating mechanisms, selectively retains and forgets sequential information, effectively capturing long-range temporal dependencies and overcoming the vanishing-gradient problem associated with traditional RNNs [24]. The CNN-LSTM hybrid framework further integrates the local feature-extraction capability of convolutional networks, enabling the recognition of multiscale temporal patterns. Li also reported that CNN-LSTM outperforms single-model architectures in influenza trend prediction [33]. The RF model achieved an R² value comparable to that of CNN-LSTM. Its core strength lies in ensemble learning, which reduces variance and enhances prediction stability; under limited sample conditions, ensemble methods remain competitive. Regarding robustness across forecast horizons, model performance followed a predictable pattern as the prediction period extended. In short-term forecasting (8 wks), LSTM, CNN-LSTM, and RF performed similarly, all maintaining high accuracy. When the horizon was extended to 12 and 26 weeks, performance declined slightly but remained generally stable. In contrast, SARIMA performed poorly overall, though its performance improved gradually as the time series lengthened, reflecting its capacity to capture longer-term trends as a classical time-series model.

In multisource data fusion forecasting, this study systematically evaluated the performance of four models: LSTM, CNN-LSTM, Transformer, and RF. The results indicate substantial differences among models in their ability to integrate multiple features. The LSTM model demonstrated the strongest feature-fusion capacity, maintaining R² values above 0.71 across all feature combinations and exhibiting the most stable performance. Under the full-feature set (ILI +  meteorology + BMS +  SBI), LSTM achieved the highest fit (R²=0.8350 and MAPE=63.78%); after introducing the SI, MAPE decreased sharply to 9.65%. The memory-cell mechanism of LSTM enables it to effectively learn complex temporal relationships among different features, an advantage that is fully realized in multisource data integration scenarios [34]. Specifically, among these integrated variables, the core value of incorporating the SI lies in equipping the models with adaptability to major public health policy shifts. As a dynamic variable, the SI dropped to zero in 2023 following the relaxation of COVID-19 restrictions [20]. However, its continued presence in the feature set implicitly signals the conclusion of policy interventions, enabling the models to effectively distinguish between intervention-influenced baselines (2020‐2022) and natural epidemic baselines (2023 onward). This mechanism proved crucial for maintaining robust forecasting performance in the first full postpandemic influenza season [35]. For instance, integrating the SI into the Transformer model under the full feature set dramatically reduced the MAPE from 93.60% to 22.74%, demonstrating its vital role in correcting post-pandemic baseline shifts. However, the successful integration of such high-dimensional data (including the SI) was not universal across all architectures. The CNN-LSTM hybrid model showed a declining trend in fitting performance as more features were added. It performed best when only SBI was included (R²=0.8113), but its performance deteriorated with additional features, dropping to R²=0.5392 in the full-feature combination. Notably, CNN-LSTM achieved the best performance in univariate ILI time-series forecasting. This pattern may suggest limitations of the hybrid architecture in modeling high-dimensional feature interactions. Convolutional operations in the CNN layer can introduce information loss when input dimensions increase, potentially impairing the subsequent LSTM layer’s ability to capture temporal dependencies [36]. Zhao et al [36] also noted that CNN-LSTM requires carefully designed network architectures and feature engineering when handling high-dimensional multivariate time series [37]. The Transformer model exhibited high sensitivity to feature composition. Some 2-factor combinations led to a pronounced performance decline, for example, R²=0.0961 for ILI +  BMS + SBI and R²=0.3739 for ILI +  meteorology + SBI. Nevertheless, under the complete feature set (ILI +  meteorology + BMS +  SBI + SI), the Transformer achieved relatively balanced overall performance. This behavior is related to the self-attention mechanism of the Transformer. As noted by Vaswani et al [28], the Transformer can directly model dependencies between any positions in a sequence, but it demands sufficient data volume and feature completeness to realize its potential. The RF model displayed comparatively weak overall fitting ability, with most feature combinations yielding R² below 0.7. However, its MAPE remained relatively stable across scenarios (23.89%‐39.47%), indicating consistent control over relative prediction error.

Regarding the contribution of individual external features, the SBI produced the most pronounced improvement in predictive performance. In the CNN-LSTM model, the ILI +  SBI combination achieved R²=0.8113, representing the best performance across all feature combinations for this model. This combination also yielded stable results in the LSTM model and outperformed other single-feature additions in the Transformer model. These findings are consistent with numerous previous studies indicating that search-engine data can effectively complement influenza activity forecasting [33]. The predictive value of search data stems from its ability to reflect real-time public health concerns and symptom-related information-seeking behavior. As noted by Su et al [31] in a study of influenza surveillance in China, when influenza begins to circulate, patients and their families often search for symptom and treatment information online before seeking medical care, resulting in a significant correlation between the SBI and influenza incidence [37]. The SI substantially enhanced the predictive accuracy of all models. After incorporating SI, MAPE decreased markedly across models. The SI synthesizes the intensity of government responses such as school closures, workplace restrictions, and movement controls. Research by Liu et al [38] reported that influenza activity in China during 2020 was substantially lower than in previous years, closely associated with stringent containment measures. The period covered in this study included the COVID-19 pandemic, during which NPIs not only curbed SARS-CoV-2 transmission but also significantly altered influenza transmission dynamics [39]. Incorporating SI into forecasting models allowed quantification of external intervention effects, thereby significantly improving prediction accuracy. BMS elicited the strongest response in the RF model. Whether added as a single factor, in 2-factor sets, or in multifactor combinations, migration data consistently improved model performance. Studies such as that of Fu and Zhu [40] have demonstrated a strong association between population movement patterns and infectious disease spread. In the deep learning models, BMS alone produced moderate performance, but when combined with other features, it exhibited synergistic effects, enhancing overall predictive efficacy. This suggests that population mobility data can also serve as a valuable auxiliary predictor.

Based on the optimal feature combinations identified for each model, this study further evaluated the performance of multisource models over short-term (8 wks), medium-term (12 wks), and long-term (26 wks) forecasting horizons to examine their robustness and applicability across different time scales. Among the models, the LSTM and Transformer models demonstrated relatively stable performance from short- to long-term forecasts. The LSTM model exhibited excellent fitting ability in short- and medium-term predictions, achieving R² values of 0.83 and 0.89 for 8-week and 12-week horizons, respectively, and maintained a high R² of 0.82 in the 26-week long-term forecast, indicating good temporal stability. The Transformer model also showed consistent predictive performance across the 3 time scales. In contrast, the performance of the CNN-LSTM and RF models improved substantially as the forecast horizon lengthened. For CNN-LSTM, R² increased sharply from 0.26 at 8 weeks to 0.81 at 26 weeks. This pattern may be attributed to the convolutional operations in the CNN layer, which require a sufficiently long sequence to effectively extract local temporal features. The shorter 8-week window may have limited feature learning, whereas the 26-week horizon—approaching half a seasonal cycle—enabled the model to capture more complete seasonal patterns [27]. Li et al [33] also noted that CNN-LSTM achieves optimal performance on periodic time series only when provided with adequately long input sequences. The RF model failed completely in the 8-week forecast (R²=–4.14) but recovered to R²=0.67 in the 26-week horizon, while its MAPE remained relatively low throughout (24.04%‐31.92%). A plausible explanation is that RF has limited capacity to model lag effects, which may exert a weaker influence in longer-term forecasts, thereby improving performance [29]. These observations suggest that the CNN-LSTM and RF models are better suited for long-term trend prediction.

Overall, deep learning methods demonstrate clear advantages in influenza forecasting. In univariate time-series prediction, the CNN-LSTM model achieved the best performance, effectively capturing the nonlinear and complex temporal dependencies in influenza data and significantly outperforming the traditional SARIMA model. All deep learning models exhibited improved performance when multisource data were integrated. Among them, the LSTM, Transformer, and RF models performed best when environmental factors, population mobility data, and internet-based data were incorporated simultaneously. The LSTM model attained the highest overall performance, with R² reaching 0.8350 in the full-feature combination. In contrast, although CNN-LSTM excelled in univariate forecasting, its performance declined as input dimensions expanded to 6 (dropping from R²=0.8113 with SBI to R²=0.5392 in the full-feature combination). This degradation likely stems from architectural limitations: standard 1D convolutions struggle to capture cross-variable interactions among highly heterogeneous features, and max-pooling may discard fine-grained interdimensional temporal details [41]. Furthermore, the expanded parameter space on our limited training set elevates overfitting risks despite dropout regularization [42]. Consequently, CNN-LSTM appears better suited for low-dimensional scenarios (1‐3 features). The value of multisource data fusion should be considered judiciously. Not all external data sources enhance predictive performance, nor can all models effectively integrate multisource information. The SBI and the Oxford SI contributed most substantially to performance improvement. The former reflects public health information-seeking behavior, while the latter quantifies the intensity of government containment measures. The SBI consistently improved the R² of all models, whereas the introduction of the SI led to a marked reduction in MAPE across every model. This suggests that incorporating data sources that reflect external environmental changes is crucial for maintaining prediction accuracy during major public health events or policy interventions.

The projection indicates that ILI incidence in Hubei Province will follow an overall trajectory of “rapid decline—gradual decline—stabilization,” suggesting that the epidemic peak will gradually subside and return to pre-epidemic baseline levels. The forecast period corresponds to the 2024 winter-spring influenza season, which coincides with the tail phase of the winter influenza peak. As temperatures rise and population gatherings decrease, ILI incidence is expected to decline rapidly to interepidemic baseline levels, marking a transition to a low-activity state. In summary, ILI incidence in Hubei Province during the 2024 winter-spring period is projected to exhibit typical seasonal waning characteristics, shifting gradually from the winter epidemic peak to a low-incidence phase by late spring and early summer. This overall trend aligns with the established seasonal pattern of influenza in Hubei Province, which is characterized by high activity in winter and spring and low activity in summer.

First, the study focused exclusively on Hubei Province. While this provided an ideal natural experiment to validate model robustness under complex post–COVID-19 dynamics—spanning from extreme low incidence during the 2020 pandemic to a historical peak in 2023—generalizability to other climatic zones and demographic populations requires further investigation. Second, the limited 209-week ILI surveillance timeframe may restrict the full exploitation of data-intensive models. Third, the 2020‐2023 training data coincided with the pandemic, exhibiting extreme heterogeneity. While this enabled models to robustly learn policy-driven dynamics (external validation R²=0.7706), their exposure to typical seasonality was limited. Future applications should use transfer learning with post-pandemic data to recalibrate model weights. Finally, as Baidu’s market share declines and platforms like WeChat and Douyin grow for health information seeking, the future representativeness of Baidu-based indices may diminish. Subsequent studies should incorporate multiplatform internet data to ensure long-term model applicability.