Estimating the Long-Term Epidemiological Trends and Seasonality of Hemorrhagic Fever with Renal Syndrome in China

Introduction

Hemorrhagic fever with renal syndrome (HFRS) is a zoonosis as a consequence of the infection with several different families of Hantaviruses, which are predominantly carried and transmitted by the natural reservoir of rodents and can be exposed to humans via the urine, droppings, or saliva of virus carriers.¹ The clinical manifestations of HFRS are characterized by different degrees of fever, bleeding headache, abdominal pain, and acute kidney damage.² The Hantaviruses are distributed globally yet varied greatly in different regions.³ In Asia, the HFRS-causing pathogenic agents predominantly include Hantaan virus (HTNV), Seoul virus (SEOV), and Amur virus (AMRV); whereas in Europe, mainly comprising Puumala virus (PUUV), Sochi virus (SOCV), and Dobrava-Belgrade virus (DOBV).^3–5 HFRS is currently found in more than 70 countries, but it is mainly endemic in Asian and European continents, especially in China, Russia, and Korea.^1,3,4 Among them, China is always hit the hardest country with HFRS in the past, accounting for around 70–90% documented case notifications throughout the world each year.⁶ HFRS still poses a major public health concern in China because of its ongoing epidemic in 28 of 31 provinces, with approximately 10,000–30,000 case notifications annually,^7,8 resulting in case-fatality rates (CFRs) from 1%-14%, contingent on which virus is causing the HFRS.^1,4 In China, the reported HFRS cases are predominately associated with the causative pathogens of HTNV and SEOV,⁹ and the HFRS-causing HTNV accounts for around 70% of the number of case notifications.⁸ Since the 1990s, a considerable reduction in the HFRS incidence has been witnessed, attributable to the implementation of such comprehensive interventions as effective rodent control, environmental management, improvement of people living standards, and vaccination in China.⁵ However, HFRS seemingly signals a recurring risk in China in recent years owing to its continued rise in case notifications and the extension of geographical distribution to 31 provinces across China.^5,8,10 Basic to any development of prevention and control measures is an accurate quantitative analysis for the secular epidemiological trends and seasonality of diseases. Therefore, it is essential to make a long-term estimation for the timing, extent, and duration of HFRS by use of a forecasting model with good reliability and strong robustness.

There are a large volume of published researches that developed different mathematical techniques to analyze and estimate the upcoming epidemics of infectious diseases for different forecasting aims. The common forecasting tools include autoregressive integrated moving average (ARIMA) model,¹¹ exponential smoothing (ES) methods,¹² generalized regression neural network method,¹³ autoregressive distributed lag method,¹⁴ grey approaches,¹⁵ back-propagation neural network technique,¹⁶ generalized linear regression models,¹⁷ support vector machine regression method,¹⁸ and autoregressive conditional heteroskedasticity (ARCH) method.¹⁹ The above-mentioned methods may provide a satisfactory result in their respective prediction domain. However, on the one hand, these models are often well suited for performing a short-term forecast as they are under either linear or nonlinear assumption,¹¹ which may limit their ability to be generalized in that a long-term ahead forecast can often offer informative clues in the decision-making process in practice. On the other hand, the existing predictive methods typically conduct a successful forecast for a simple time series (eg, single seasonal time series), and yet it is difficult to handle the complex time series with multiple seasonal periods, high-frequency seasonality, non-integer seasonality, and dual-calendar effects.^20,21 Importantly, it has been demonstrated that HFRS shows bimodal seasonal patterns in many countries such as China,^6,10 Russia,⁴ and Korea,²² which may lead to an inaccurate forecasting result for the upcoming epidemics of HFRS using the common predictive methods. At this time, the advanced innovation state-space modeling framework by combining Box-Cox transformations, Fourier series with time-varying coefficients and autoregressive moving average (ARMA) error correction method (called TBATS) is tailored for extraction with the information included in a complex time series described above.²¹ Furthermore, the TBATS approach can not only be used to deal with the linear issue but can handle some types of non-linearity based on Box-Cox transformations,²³ whilst allowing for time-dependent dynamic seasonality,²⁴ which enables this advanced method to have the potential to make a long-term forecast. Moreover, the TBATS approach also has a powerful ability to decompose a complex seasonal time series into constituent latent subseries,²¹ which fails to be done by the aforementioned common methods. In the past, the advanced TBATS model has been used to nowcast and forecast the flow time series in water distribution systems,²⁵ the demand for the electric energy,²⁶ and spatiotemporal variations and characterization of benzene,²³ and the resulting results indicated that this advanced method could produce a relatively high forecasting accuracy in these fields. However, there is no study to use this advanced method to analyze and evaluate the long-term epidemiological trends and seasonality of HFRS. Considering the attractive performance of the TBATS approach and the epidemic patterns of HFRS in China, the primary goal of this study is to describe and forecast the long-term epidemic trends and seasonality of HFRS with the advanced TBATS method. Meanwhile, to demonstrate the potential of this advanced method and its suitability for the application in forecasting the long-term spreading patterns of HFRS, its predictive accuracy level was compared with the ARIMA method that has been deemed as the most common and useful methodology in the different domains of time series forecasting.^6,7,27–29

Materials and Methods

Data Collection

Currently, there are 40 notifiable contagious diseases in China. The confirmed cases of these diseases must be reported to the National Notifiable Infectious Disease Surveillance System (NNIDSS). HFRS is among them. Thus, we retrospectively took the monthly HFRS incidents and population data between January 1995 and December 2020 from the internet-based NNIDSS and the Statistical Yearbook of China, respectively. Finally, a total of 312 observations covering 26 years were gathered. Often, to build an effective and stable method, preferably 100 observations or more are expected to be adopted.³⁰ In this study, we aim to conduct a long-run forecast for the HFRS epidemics. To this end, we split the whole data into six different training datasets and testing datasets (ie, the forecast horizons included 12, 24, 36, 60, 84, and 108 data) so that a series of experiments were performed to validate the generalization and the forecasting robustness of the methods.

The study protocol was approved by the institutional review board of Xinxiang Medical University (No: XYLL-2019072), and the ethical approval is not warranted for this research as these data without personal information are publicly available.

Developing SARIMA Method

Time series often shows a noticeable link between successive data. The ARIMA method is introduced to make a forecast by mining the intrinsical attributes and inherent rules of time series data.³¹ HFRS frequently has notable seasonal effects,^2,10 and hence the seasonal ARIMA (SARIMA) method should be adopted. In this method, the seasonality of HFRS was thought of as the predictors and the monthly HFRS incidence as the response variable. The SARIMA method is composed of three parts:²⁸ the “AR” signifies the autoregressive approach, the “I” denotes integration, and the “MA” refers to the moving average approach. The basic notation of this method can be written in the form of SARIMA(p, d, q)(P, D, Q)_s, where p is the orders of AR, d is the orders of integration, and q is the orders of MA, (P, D, Q) correspond to the seasonal components of the above three parameters.²⁸ In creating the best SARIMA method, the most key step is to search for its suitable parameters. Currently, most software programs allow an automatic selection of the most appropriate parameters. The “Expert Modeler” in SPSS software provides an elegant way of choosing the preferred parameters of the SARIMA method.³² Of all the plausible parameters, those that maximize the R-square (R²) and stationary R², and minimize the Normalized Schwarz’s Bayesian Information Criterion (NBIC) were first selected.³³ Subsequently, the t-test was used to investigate the statistical significance of the identified parameters.²⁹ Next, autocorrelogram, partial autocorrelogram, and Ljung-Box Q test were applied to the residuals to investigate whether the produced forecasting errors behaved like a white noise series.²⁹ Once the method passed all the required tests, it can be said that the selected parameters are suitable and adequate. Finally, a forecast into the next periods can be done with this best approach.

Developing TBATS Method

The traditional single seasonal ES models have a limited capacity to deal with such a complex time series composing of multiple seasonal behaviors, non-integer seasonal behaviors, and dual-calendar effects,²¹ despite their widespread use in practice.^12,34 Some researchers strive for the extension of the single seasonal Holt-Winters’ method by accommodating a second seasonal pattern to analyze the time series including two seasonal behaviors.²⁰ However, this extended model was shown to suffer from over-parameterization as there are a large number of parameters that require to be calculated for the preliminary seasonal behaviors, especially when the target series is composed of high-frequency seasonal behaviors, further heightening the need for the modifications of the traditional ES approaches to deal with a wider variety of seasonal patterns and to resolve the issues raised above. To do this, researchers pay more attention to linear homoscedastic approaches but consider some types of non-linear patterns with Box-Cox transformations, in view of the complexity by extending non-linear ES approaches to deal with the complex time series described above, and the BATS(p, q, m₁, m₂,…, m_T) method was hence introduced, herein B signifies the Box-Cox transformation, A signifies the ARMA model, T and S signify the trend and seasonal patterns in the target series, respectively. The parameters include the orders (p and q) of the ARMA approach and the seasonal periods (m₁,…, m_T). Albeit the BATS method is successful in modeling multiple seasonal periods in the target series, it fails to address the non-integer seasonality issue and needs to estimate a considerably large number of parameters since the initial seasonal patterns separately include m₁ + m₂ + · · · + m_T parameters, resulting in too many values when the target series contains high-frequency seasonal components. To deal with the problems raised above, researchers made a further generalization of the BATS method by introducing a novel trigonometric representation of seasonal behaviors based on the Fourier series. By doing so, this new method can significantly not only cut the number of parameters that need estimation but can model the non-integer seasonality in the time series.²¹ This advanced modified method can be denoted as TBATS(ω, p, q, φ, {m₁, k₁}, {m₂, k₂},…,{m_T, k_T}), herein k_i represents the corresponding number of Fourier terms used for each seasonality, ω signifies the BoxCox parameter, and φ refers to the damping parameter that facilitates the trend extrapolation of the model when the trend pattern is dampened.²³ Additionally, the use of trigonometric functions also helps the TBATS method to decompose a complex time series into secular trend, seasonality, and error components.²⁵

Statistical Analysis

The changing epidemiological trends of HFRS were estimated using average annual percentage change (AAPC) and annual percentage change (APC) based on the joinpoint regression program (Version 4.8.0.1). The SARIMA method was created using the “Expert Modeler” function in SPSS software (version 17.0, IBM Corp, Armonk, NY) and the TBATS method was built with R software (Version 3.4.3, The R Foundation for Statistical Computing, Vienna, Austria). The incidence rate ratio (IRR) with 95% confidence intervals (CIs) between different variation trends of HFRS was computed by use of comparison of two proportions.³⁵

To see the predictive accuracy, two types of measurement indicators were computed, including the scale-dependent measurement indicators comprising the mean absolute deviation (MAD) and root mean square error (RMSE), and the measurement indicators based on percentage errors comprising the mean absolute percentage error (MAPE), mean error rate (MER), and root mean square percentage error (RMSPE). The smaller values of these measures above corresponded to a better model.

(1)

where represents the monthly HFRS incidence data, refers to the forecasting data, stands for the mean of the HFRS incidence data, and is the number of forecasts.

Results

Statistical Description

There were 576,361 notified HFRS cases over the period 1995–2020, resulting in yearly and monthly average morbidity rates of 1.689 and 0.141 per 100,000 people, respectively. According to the resultant results of joinpoint regression (Figure 1), overall a notable reduction in HFRS morbidity was seen, with AAPC=−6.767 (95% CIs: −10.592 to −2.778; t=−3.278, p=0.001), and yet different stages had different variation trends: a decreasing trend during 1995–2004, with APC=−8.964 (95% CIs: −11.667 to −6.179; t=−6.547, p<0.001); during 2004–2008, with APC=−37.153 (95% CIs: −37.153 to 3.279; t=−1.828, p=0.084); during 2008–2020, with APC=−0.349 (95% CIs: −3.992 to 3.431; t=−0.197, p=0.846), and it seemed to show a cyclical fluctuation with 3–5 years since 2008. The average incidence for the period 2008–2020 decreased by a factor of almost four compared to that for the period 1995–2004 (0.805 vs 3.065 per 100,000 persons; IRR=3.807, 95% CIs: 3.786 to 3.828), decreased by 1.6-fold than that for the period 2004–2008 (0.805 vs 1.324 per 100,000 persons; IRR=1.643, 95% CIs: 1.630 to 1.656). Visibly, as depicted in Figure 2, people could be infected with HFRS all through the year. Nonetheless, the HFRS epidemics presented pronounced dual seasonal patterns per year, with a weak peak between May and June, a secondary strong peak between November and December, and a trough between August and September.

Figure 1 Joinpoint regression plot displaying the HFRS epidemiological trends from 1995 to 2020. *Showed that the annual percent change (APC) is statistically significant.

Figure 2 Probability density plot showing the monthly HFRS incidence. As depicted, HFRS showed notable seasonal variation, with a strong peak in November and December, and a weak peak in May and June.

The Optimal SARIMA Method

In constructing the SARIMA method, the “Expert Modeler” function was run on the in-sample data between January 1995 and December 2019, which automatically identified the SARIMA(0,1,(1,3))(0,1,1)₁₂ specification as the preferred model for the 12-data ahead forecast since among all the possible methods, the parameters of this method corresponded to the maximum R² value of 0.933 and stationary R² value of 0.473, along with the minimum NBIC value of 11.702 (Table 1). Further diagnostic checking for this optimal method showed that the resulting key parameters were statistical significance (MA1=0.355, MA3=0.253, SMA1=0.747; all p-values<0.001), and the resulting forecast errors could be thought of as an uncorrelated series in that there were little sample autocorrelation function (ACF) and partial ACF (PACF) outside the 95% CIs and there was no statistical difference under Ljung-Box Q test (Q₍₁₈₎=13.709, p=0.548) (Table 1 and Figure 3). The above diagnostic results meant that the identified SARIMA (0,1,(1,3))(0,1,1)₁₂ specification was capable of tracking the epidemiological trends of HFRS effectively and sufficiently. Likewise, the SARIMA (1,1,(2,3))(0,1,1)₁₂, SARIMA (0,1,(1,3))(0,1,1)₁₂, SARIMA (0,1,1)(0,1,1)₁₂, SARIMA (1,0,0)(0,1,1)₁₂, and SARIMA (1,0,0)(0,1,1)₁₂ specifications were considered as the best model to forecast the 24, 36, 60, 84, and 108 holdout data, respectively, and the resulting diagnostic overviews for the five best SARIMA versions were visible in Table 1 and Figures S1–S5. Afterward, projecting into different future periods could be done by use of these best SARIMA methods (Figure 4 and Tables S1–S6).

Table 1 The Best-Performing SARIMA Models Obtained on the Various Training Sets and Their Parameter Estimations and Goodness of Fit Tests

Figure 3 Diagnostic test plots for the forecasting errors of the SARIMA model created with the data between 1995 and 2019. (A) Autocorrelation function (ACF) plot, (B) partial autocorrelation function (PACF) plot, and (C) p-values for Ljung-Box statistic. We could see from the sample correlogram that there were little sample ACFs and PACFs touching the significance bounds, except for that at lags 9, 23, and 24, and the p-value was greater than 0.05 under Ljung-Box test statistic. These results meant that the selected SARIMA approach provides an adequate predictive model for the HFRS incidence.

Figure 4 Time series plots showing the forecasted results on different datasets between the SARIMA models and the TBATS models. (A) 12-step ahead forecast, (B) 24-step ahead forecast, (C) 36-step ahead forecast, (D) 60-step ahead forecast, (E) 84-step ahead forecast, and (F) 108-step ahead forecast. The out-of-sample predictions are shown as a shaded area in these plots. It was discovered that the out-of-sample predictions under TBATS approaches agreed better with the observed values over the SARIMA approaches, especially for the long-term out-of-sample predictions.

The Optimal TBATS Method

The predictive results under TBATS method depend largely on the number of Fourier terms used for each seasonal component, and it is difficult to consider all potential combinations when detecting the most adequate harmonics k_i in application.²¹ It was demonstrated that the use of Akaike’s Information Criterion (AIC) is of great help for automated method selection.²³ Therefore, in choosing the best TBATS method, taking one seasonal pattern each time into account, we simulated the model to the training series repeatedly by the gradual increase of the numbers of k_i but keeping all other harmonics constant each time, until the lowest AIC value is detected. In searching for the best TBATS method, appropriate parameters of the orders (p and q) of the ARMA method must also be selected. Fortunately, these two values can be determined by use of the automatic ARIMA algorithm to the forecast errors under the model with optimal k_i combination above.³⁶ If the identified method with ARMA component gives lower AIC compared to that without ARMA component, this model would be retained; otherwise, the ARMA component was not included in the model. After a series of attempts, it was discovered that the TBATS (0.27, {0,0}, -, {<12,4>}) specification was expected to be the preferred model for the forecast of 12 holdout periods because there was a lower AIC value of 5012.207 in this TBATS model than in others, and the resulting key parameters, statistical checking and constituent latent subseries are illustrated in Table S7 and Figures S6 and S7. Similarly, according to the modeling processes above, the TBATS (0.285, {3,0}, -, {<12,4>}), TBATS (0.289, {0,3}, -, {<12,4>}), TBATS (0.302, {0,0}, -, {<12,4>}), TBATS (0.34, {0,0}, -, {<12,4>}), and TBATS (0.325, {0,0}, -, {<12,4>}) specifications were selected as the best models, respectively, for the 24-, 36-, 60-, 84-, and 108-step ahead forecasts, and the resulting additional summary statistics for their key parameters, statistical checking and constituent latent subseries are reported in Table S7 and Figures S6, S8–S12. Subsequently, a prediction into the upcoming HFRS epidemics was able to be achieved based on these best TBATS versions (Figure 4 and Tables S1–S6).

Forecasting Performance Evaluation Between Methods

Table 2 summarizes the measurement metrics testing the predictive accuracy on different datasets under the best SARIMA and TBATS approaches, and the resulting data pinpointed that the advanced best-fitting TBATS methods presented notably lower values of MAD, MAPE, RMSE, RMSPE, and MER over the best-fitting SARIMA methods, with improvements of at least 18% and at most 88% in the accuracies for analyzing both the short- and long-term trends, although the forecasting accuracy displayed a slight decrease with the increase of holdout periods (Table 2 and Figure 4). Considering the advantage of the TBATS method in the estimation of the HFRS epidemics, we re-constructed the best-performing TBATS (0.286, {3,0}, -, {<12,4>}) model based on the 26 entire years of data to nowcast and forecast the HFRS morbidity into 2030 accordingly, and seemingly projecting a plateau of around 90,120 (95% CIs: 25,358 to 244,667) cases, with monthly average 751 (95% CIs: 211 to 2039) case notifications in the next ten years (Figure 5).

Table 2 Comparisons of the Out-of-Sample Forecasting Powers Between SARIMA Models and TBATS Models

Figure 5 Estimated epidemiological trends and seasonality of HFRS between January 2021 and December 2030 using the best-performing TBATS (0.286, {3,0}, -, {<12,4>}) model.

Discussion

HFRS poses a serious public health threat in China. Accurately nowcasting and forecasting its epidemic trends and seasonality in the next years are capable of providing the stakeholders some significant hints to develop effective and feasible plans to address the evolving scenarios.³⁷ Although numerous forecasting tools have been applied to do this,^38–44 most existing predictive models may be inadequate in analyzing time series with complex seasonal patterns due to their different assumptions,⁴⁴ enhancing the need for the use of the statistical approaches to handle complex seasonal traits. The advanced TBATS approach is tailored for suiting all the requirements of complex time series.²¹ However, studies on the use of the TBATS approach to analyze the long-term epidemic trajectories of contagious diseases are rare. Therefore, the current study, for the first time, focused on the exploration of the utility and adequacy in nowcasting and forecasting the long-term epidemiological trends and seasonality of HFRS, and its predictive accuracy level was compared with the most versatile SARIMA method in the field of time series modeling.^11,13,16,33 It was discovered that the advanced TBATS methods get a more clear perspective of capturing the dynamic dependency structure in the spreading of HFRS over the SARIMA methods in a series of modeling experiments. What’s more, the long-term forecasting results to emerge from the advanced TBATS method remained reliable and robust with the increase of prediction time windows, despite a slight rise in the values of measurement metrics, including MAD, MAPE, RMSE, RMSPE, and MER. By contrast, the SARIMA model has been unable to analyze and capture the long-term trends of HFRS after 60-step ahead predictions as it produced a MAPE value (which is frequently computed to assess the forecasting accuracy level, a forecast with MAPE value greater than 50% is considered unacceptable⁴⁵) greater than 50% in both 84 and 108 holdout periods, further providing evidence that the SARIMA method is valid for assessing short- or medium-term trends.⁴⁶ Our recent work suggests that the Error-Trend-Seasonal framework (ETS) technique has many attractive applications in analyzing the long-run temporal behaviors.^33,46 Accordingly, we further established the ETS approach based on the HFRS incidence data to nowcast and forecast its epidemics, and the resulting results indicated that the TBATS method is also more useful and robust for describing the long-term trends of HFRS than the ETS technique, except that in 60-step ahead forecast (Table S8). Also, we used the HFRS incidence data in Liaoning province and Shanxi province (which were hit the hardest with HFRS in the past decades in China) to validate the forecasting reliability of TBATS model. The remarkable results to emerge from these data suggested that the advanced TBATS method still generated a higher forecasting accuracy compared to the SARIMA model (Table S9). As a result, the advanced TBATS method is able to be recommended as a promising and powerful alternative to analyze the long-run trends and seasonality of HFRS in China. In the meantime, this advanced TBATS approach also provides a serviceable aid for directing the present intensity and type of public health measures.⁴⁷ For instance, if the temporal patterns of HFRS were to start to recede, while this advanced method assessed a marked upward tendency, intimating that the current interventions would be at work in containing the epidemics; otherwise, the present measures required to be optimized.

Time series forecasting is a useful aid for forming hypotheses to understand the epidemic behaviors of different diseases and to estimate their spreading dynamics, and thus favoring the establishment of a quality control system.³³ The SARIMA method has been the most popular time series forecasting tool to analyze the morbidity, mortality, and prevalence of infectious diseases (eg, HFRS,^40,41 Human Brucellosis,⁴⁸ COVID-19,³³ hand-foot-mouth disease (HFMD),²⁸ Syphilis,⁴⁹ dengue,⁵⁰ Dysentery⁵¹) over the last years because of its simple structure, fast applicability, relatively high prediction accuracy, and potential to interpret the target series. However, as evidenced by the findings from our study and others, the SARIMA method is only appropriate for analyzing short- or medium-term trends and cannot deal with the nonlinear information in a complex time series due to its linear assumption.^13,16,52,53 Instead, the advanced TBATS method was developed by adding the trigonometric representation of seasonal components based on Fourier series to the traditional BATS approach, which allows it to not only deal with all of the complex time series but also to handle both linear and non-linear components,^23,25 whilst accommodating a dynamic seasonal patterns over time,²⁵ and hence leading to its usefulness and flexibility in analyzing the long-term trends and seasonality of time series. Given its attractive applications of the TBATS method and our previous work,⁵³ it appears that the importance of this advanced method as a powerful forecasting tool is expected to be emphasized in analyzing the long-term temporal levels of HFRS in other regions or other communicable diseases, and yet validation needs to be done. Besides, some new advanced statistical techniques (eg, bayesian structural time series technique,³⁷ flexible transmitter network,⁵⁴ the optimized theta method,⁵⁵ and age-structure mathematical model⁵⁶) are also recently demonstrated to have the potential to make a long-term forecast for the time series. As a result, what is now required is work focused on the comparisons of the forecasting accuracy between the TBATS approach and the above-mentioned models.

Different from prior studies that reported a re-emerging risk in the HFRS morbidity in China,^5,6,8,10 our study indicated an evident reduction in the numbers of HFRS cases over the 26 years of data (AAPC=−6.767, 95% CIs: −10.592 to −2.778; t=−3.278, p=0.001), consistent with the epidemiological trends reported in several countries in Asia,⁵⁷ and yet in disagreement with that in Russia.⁴ This discrepancy may be related to the different causative agents, geographical environments, social economy, preventive measures, and public awareness. In China, many enhanced prevention and control measures, including effective environmental management, elimination of rodents, vaccination, publicize education, improvement of the living and working environments, along with rapid urbanization and mechanization of farming, have favored the decrease in the HFRS-related morbidity and mortality in the past.⁵⁸ Even so, it was discovered that HFRS displayed changing epidemiological trends in different stages, a remarkable reduction during 1995–2004 (APC=−8.964, 95% CIs: −11.667 to −6.179; t=−6.547, p<0.001) and during 2004–2008 (APC=−37.153, 95% CIs: −37.153 to 3.279; t=−1.828, p=0.084), and a relatively stable epidemic trend during 2008–2020 (APC=−0.349, 95% CIs: −3.992 to 3.431; t=−0.197, p=0.846). Importantly, such an epidemiological trend would be predicted to continue in the ten years to come under TBATS model. Besides, it was noted that HFRS showed a cyclical fluctuation with 3–5 years since 2008, in agreement with the findings in previous reports,^59,60 which also indicated that periodic outbreak is an important epidemiological feature in the HFRS incidence, albeit these studies reported a longer time window of 7–12 years.^59,60 It seems that several possible reasons contribute to the changing epidemiology in the HFRS incidence in recent years; one may be associated with the effects of climate variability, responsible for many emerging or re-emerging infectious diseases around the world.^8,61 HFRS has been thought of as a climate-sensitive illness since meteorological parameters play an important indirect or direct role in the infection rates, habitats, reproductive rates, and incubation period of rodents.^8,62 Another may be closely correlated with the rapidly growing tourism in China; tourists can be infected with the pathogenic agents of HFRS in the endemic areas and then can be spread to the non-epidemic areas,⁶³ which seems to well explain the geographical expansion of HFRS in recent years.¹⁰ The third may be due to the emergence of new Hantavirus subtypes.⁹ Recent several studies have suggested that the novel identified Seoul viruses may pose a great threat to the control of HFRS.^9,64 Finally, the continued decrease in the HFRS incidence in the past may also have relaxed the interventions and public awareness, leading to the changing epidemiology and even the recurring risk.⁶⁵

Seasonal profile of HFRS has been reported in different countries.^4,22,58 A marked semi-annual seasonal behavior was also seen in our study, with a weak peak between May and June, and a strong one between November and December. Besides, the predicted HFRS morbidity also showed the same seasonal pattern as that during 1995–2020 (Figure S13). Our seasonal profile concurs well with previous findings from Russia⁴ and South Korea,⁶⁶ also in line with that in most provinces or cities from across China (eg, Changsha city,⁶⁰ Heilongjiang province,⁶⁷ Anqiu city,⁷ Zibo city,⁶⁸ Shanxi province,⁶⁹ Zhejiang province⁷⁰). We surmised that different climatic features and etiological factors contributed to this difference in the amplitude of peak activities. In China, it was discovered that the HFRS-causing HTNV can occur throughout the year, and yet the majority of them became infected with HTNV in fall and winter; by contrast, the HFRS-causing SEOV is mainly reported in spring.^6,10 The main natural reservoir hosts of the above two causative factors are predominantly spread by the striped field mouse (A.agrarius) and brown Norway rat (R.norvegicus), respectively.⁸ Weather factors (eg, temperature, rainfall, relative humidity, precipitation) have been shown to have a great influence on the population dynamics of the main natural reservoir hosts of HTNV and SEOV,^8,62,64,71 and the effects of temperature and relative humidity on HFRS was demonstrated to show a U-shaped curve,⁶⁴ this fits well with the bimodal seasonal distribution present in spring and winter in the HFRS case notifications.

Our findings emerged from the present investigation successfully indicate the usefulness and flexibility of the TBATS approach in describing and forecasting the long-term seasonality and trends of HFRS. Nevertheless, this work also suffers from certain shortfalls. First, under-reporting and under-diagnosis of HFRS cases may be inevitable in passive surveillance, notwithstanding the well-monitored data quality under mandatory reporting in China. Second, we noted from Tables S1–S6 and Figure 5 that the resulting predicted 95% CIs visibly broadened with the increase of the forecasted time windows, despite its high forecasting accuracy and robustness of the TBATS method. Thus, the model should be updated in time by including new series to ensure its forecasting reliability level. Third, from Figure S6 we can see that the forecasting potential of the TBATS model may be further improved, despite its high forecasting accuracy. Therefore, further work needs to be performed to establish whether the predictive ability of the TBATS model can be improved by combing this model and neural network techniques. Fourth, the current findings from this study reflect how well the TBATS method estimates the long-run seasonality and trends of HFRS in China, whether it is valid for analyzing the long-term trends and seasonality in other countries or other communicable diseases, additional studies are required to investigate. Finally, the multivariate TBATS method may provide a higher forecasting accuracy. However, we are unable to do this owing to the unavailability of such a TBATS model with regression covariates and HFRS-associated factors in this study.

Conclusion

The interesting findings to emerge from this work reveal that the TBATS method does a better job of describing and forecasting the long-term seasonality and trends of HFRS, which is capable of being deemed as a feasible and serviceable alternative to help stakeholders to inform future preventive policy or practical solutions to control the HFRS epidemics in China. Besides, under current interventions, it is challenging to eliminate the HFRS epidemics in the next ten years considering the resultant forecasts of around 90,120 (95% CIs: 25,358 to 244,667) cases, with monthly average 751 (95% CIs: 211 to 2039) case notifications in the ten years to come using the TBATS method, intimating that additional or comprehensive prevention and control measures must be formulated to tackle this evolving scenario.