Time Series Analysis and Forecasting of the Hand-Foot-Mouth Disease Morbidity in China Using An Advanced Exponential Smoothing State Space TBATS Model

Introduction

Hand-foot-and-mouth disease (HFMD) is a common viral infectious disease that mainly affects children under seven.¹ This illness is typically a self-limiting condition and most children are characterized by oral pain, anorexia, low fever, and popular vesicles on the palms and soles, whereas a few children can develop myocarditis, pulmonary edema, and aseptic meningoencephalitis, which may lead to death. Coxsackievirus A16 (CVA16) and enterovirus 71 (EV71) are the most common causative agents causing HFMD outbreaks in the past years, even though some other pathogens can be involved.^2–4 Moreover, it was discovered that around half of the children were coinfected with more than one pathogen. The Coxsackie viruses were first described in 1948,⁵ and particular attention was not paid to HFMD until 1998 because this illness caused several outbreaks, associated with millions of children, especially in the Western Pacific Region (eg, Malaysia, China, Japan, Singapore, Mongolia, Vietnam, and Cambodia^6–11). In China, several large outbreaks of HFMD were observed in 2007 and early 2008, leading to major burdens of disease and economy. In response to such a daunting threat, HFMD was thus listed as a notifiable disease in May 2008 in China,¹² and since then its morbidity often ranks the first among the 40 notifiable infectious diseases in China.¹³ Currently, three monovalent EV-A71 vaccines have been introduced against HFMD since 2016 in China, and the protective efficacies of the above-mentioned vaccines were shown to be higher than 90% against EV71-related HFMD.¹⁴ However, it is documented that this illness still has a significant influence on the general susceptible population, with more than two million children developed HFMD per year in China, and more crucially, the HFMD incidence appears to exhibit an increasing tendency in recent years, partly attributable to the emergence of new Coxsackievirus, climate change, and environmental degradation.^14,15 Therefore, to attenuate or contain the spreading of this illness, early warning for the temporal patterns of HFMD epidemic trends in the upcoming years using a suitable forecasting model plays an important role in developing effective preventive measures.¹⁶

In the past, different mathematical simulation models, including autoregressive integrated moving average (ARIMA) method,¹³ artificial neural networks (ANN_S), exponential smoothing (ES) method,¹⁷ support vector machine (SVM),¹⁸ decomposition methods,¹⁸ and grey model¹⁹ have been applied to forecast the epidemics of communicable diseases. The occurrence of communicable diseases is typically limited by different influencing factors (eg, meteorological variables, government interventions, vaccines, and air quality^15,20), leading to showing complex epidemiological features with multiple seasonal periods, high-frequency seasonality, non-integer seasonality, and dual-calendar effects of the time series.²¹ For this reason, the above frequently used forecasting tools may be incapable of mining information in such a complex time series reasonably well. Importantly, the incidence series of HFMD has been shown to exhibit complex seasonal patterns in different regions or countries.^12,22–24 To overcome the weaknesses and inadequacy of the existing time series models in dealing with complex seasonal patterns, an advanced exponential smoothing state space framework by combining Box-Cox transformations, Fourier representations with time-varying coefficients and autoregressive moving average (ARMA) error correction (TBATS) model was customized to accommodate such a time series with complex seasonal behaviors.²¹ The TBATS model is introduced based on a trigonometric formulation capable of handling the seasonal complexities in a time series. Also, this model is designed to overcome the defects of the conventional BATS models that require a large number of parameters to be estimated in the process of BATS-development and are unable to deal with non-integer seasonality.^21,25,26 However, there is a current lack of using the TBATS framework to model and forecast the temporal patterns of HFMD. Considering the advantage of the TBATS framework and the complex seasonal patterns of HFMD, the purpose of this work is to use the TBATS framework to describe and predict the epidemic behaviors of HFMD. In the forecasting field of communicable diseases, the most commonly used mathematical methods is the ARIMA model,^13,27–30 and it has been demonstrated to show a high forecasting accuracy. Therefore, in this study, the ARIMA model was also created to nowcast and forecast the epidemic trends of HFMD, and its forecasting accuracy level was compared with the TBATS method to further investigate the flexibility and usefulness of the TBATS method.

Materials and Methods

Data Collection

In China, the notifiable infectious diseases were asked to report to the Chinese Center for Disease Control and Prevention (CDC). Since HFMD was named as class C notifiable infectious disease in May 2008, and hence the monthly notified cases between January 2009 and December 2019 spanning 11 years were drawn from the National Notifiable Infectious Disease Surveillance System (NNIDSS) operated by the Chinese CDC, and the population data were drawn from the Statistical Yearbook of China during the same period. Suppose that the number of populations remained relatively steady during the year, the monthly HFMD incidence rate was computed by the average populations of each year (Table S1). Generally, at least 50 and preferably 100 observations or more are required to construct a robust and effective forecasting model.³¹ Therefore, the dataset including over 100 observations was treated as the training horizon to build the model, and the rest as the testing horizon to verify its generalization of the model.

The study protocol was approved by the institutional review board of Xinxiang Medical University (No: XYLL-2019072), and the ethical approval failed to be needed since the reported HFMD cases without personal information are publicly available in China.

SARIMA Model

Typically, time series is characterized by noticeable correlations between successive observed values.³² The most classical approach to consider the association patterns of a time series is the ARIMA model.²⁹ Since the incidence series of infectious diseases often shows marked seasonal variation and periodicity, and thus the seasonal ARIMA (SARIMA) is more appropriate for capturing the dynamic dependence structure in the HFMD incidence.^13,33 In this model, the seasonality of the object series was regarded as the independent variable and the morbidity as the dependent variable.³³ SARIMA method can be presented with the standard notation of SARIMA (p,d,q)(P,D,Q)_s, where p stands for the order of autoregression, d signifies differenced times, and q refers to the order of moving-average, (P,D,Q) are their corresponding seasonal counterparts, and S is the cycle length of the object series (S=12 in this study). Often, the SARIMA model is created based on four modelling steps. First, model identification. The SARIMA is defined for stationary time series.³⁰ Therefore, the stationarity of HFMD incidence series was detected using an augmented Dickey-Fuller (ADF) test, if suggesting a nonstationary series, the logarithm or square root transformed method or/and differenced method would need to be used until a stationary series was achieved.²⁹ Second, parameter estimation. By examining the correlogram and partial correlogram of the differenced series to identify the plausible values of the key parameters of the SARIMA model. Among the possible models, the parameters that provided a less value of Akaike’s information criterion (AIC) and Schwarz’s Bayesian criterion (BIC) and a greater value of Log-Likelihood were considered the preferred.³³ Third, model diagnostic checking. Using the Ljung-Box Q test to check whether the residual series behaved like a white noise series and the key parameters were statistically significant.³³ Finally, when the best model passed all the required statistical tests, and then a forecast into the future could be done using this best-fitting model. Otherwise, the above steps should be repeated until an appropriate model was detected.

TBATS Model

In the past, the well-documented seasonal ES methods have widely been adopted to deal with the single seasonal time series.¹⁷ Nevertheless, these traditional ES methods fail to describe the complex seasonal time series comprising multiple seasonal patterns, non-integer seasonal patterns, and dual-calendar effects.^21,34 Although some researchers attempted to extend the classical ES models to accommodate a second seasonal trait so that the time series comprising two seasonal components can be described, this may suffer from a major flaw that needs to compute substantial values for the preceding seasonal patterns, particularly when a time series shows high-frequency seasonal behaviors, resulting in over-parameterization.^21,26,35 To accommodate a wider variety of seasonal components and to resolve the issue encountered in the extended version of the common ES models. The notation BATS (p, q, m₁, m₂, …, m_T) method is thus proposed, where B refers to the Box-Cox transformation, A represents the ARMA model, T is the trend trait in the object time series, and S signifies the seasonality in the object time series.^21,26,34 The key parameters of the BATS model are composed of the ARMA method (p and q) and the seasonal periods (m₁, …, m_T). The BATS model has the potential to handle multiple seasonal periods and to allow for some types of nonlinear characteristics included in the object time series with Box-Cox transformation.³⁶ Despite its advantages of the BATS model, it is unable to deal with non-integer seasonality and also encounters many parameters that require computation.²¹ Hence, researchers proposed a novel trigonometric representation of seasonal traits based on the Fourier series, and this representation can be in the form of TBATS (ω, p, q, φ, {m₁, k₁}, {m₂, k₂}, …,{m_T, k_T}), where the parameters (p, q, and m) of the TBATS model have the same meaning as the BATS model, k is the number of harmonics used for the seasonal trait, ω stands for the Box-Cox transformation, and φ signifies the dampening parameter value.^21,37 Importantly, The TBATS method can not only significantly reduce the number of parameters needed (with an estimation of 2(k₁+k₂+ … +k_T) initial seasonal values) relative to the traditional BATS method, but can be used to estimate the non-integer seasonal frequencies with the trigonometric functions.^21,37 Also, the TBATS model has the potential to decompose a time series into constituent different components. Especially, this model is capable of identifying and extracting one or more seasonal traits that may fail to be presented in the object series graphs themselves.³⁷

Statistical Analysis

The SARIMA and TBATS methods were created using R software (version 3.4.3, R Development Core Team, Vienna, Austria). We set two predicted segments to account for the model uncertainty in the forecasting process (ie, the first 108 (from January 2009 to December 2017) and 120 (from January 2009 to December 2018) observations were used for the training parts, respectively, and the rest of 24 (from January 2018 to December 2019) and 12 (from January to December 2019) observations were reserved for the testing parts, respectively). In this study, a Lagrangian Multiplier (LM) test was employed to check whether there existed conditional heteroskedastic behavior and volatility (ARCH effect) in the residual series produced by these two models.³⁸ The statistical significance level was set at p<0.05.

We used two types of evaluation metrics to investigate the forecasting accuracy level between these two models. One is the scale-dependent metrics including root mean square error (RMSE) and mean absolute deviation (MAD), another is the percentage error metrics including mean error rate (MER) and mean absolute percentage error (MAPE). The model that has lower values of the above measurement metrics was considered higher accuracy in the time series forecasting of HFMD morbidity.

RMSE= (1))

MAD= (2)

MER= (3))

MAPE= (4))

where X_i is the observed values, is the predicted values of the two models used, stands for the mean of the observed values, and N signifies the number of forecasted values of the two models used.

Results

Descriptive Analysis

During the period 2009–2019, the reported number of HFMD totaled 22,227,647 cases, leading to annual and monthly average incidence rates of 147.710 and 12.309 per 100,000 people, respectively. From Figure 1 we can see that HFMD cases were notified across the whole year, and yet the HFMD incidence showed a high incidence characteristic every other year. Based on the 11 years of observations, the maximum number of HFMD cases were noted in May 2014, with 528,777 cases (38.658 per 100,000 people), and the minimum number of HFMD cases were observed in January 2009, with 7672 cases (0.575 per 100,000 people). When a seasonal-trend decomposition procedure based on Loess (STL) method was adopted to characterize the different components in the HFMD incidence, a notable seasonality and dual peak distribution were observed, there were few cases between January and March (especially in February), a sudden escalation between April and July (particularly in June), then a reduction in August, and yet a secondary peak in September per year (Figure 1).

Figure 1 Time series displaying the HFMD incidence from January 2009 to December 2019 and the decomposed trend, seasonal, and random traits using the classical multiplicative decomposition method.

Constructing the SARIMA Model

By a visual inspection of the time series plot showing the original data, the logarithm transformed data, and the square root transformed data of the HFMD incidence (Figure S1), suggesting a similar changing trend over time between these three series. After repeated attempts (Figure S2), we found that the original data of the HFMD incidence with logarithm transformation seemed to be more suitable for the development of the SARIMA model to make a long-term forecast (from January 2018 to December 2019). As a result, we first seasonally differenced the logarithm transformed series from January 2009 to December 2017 to remove its strong seasonal pattern (Figure S3). The differenced series was shown to be stationary (ADF=−5.720, p<0.001), which suits the prerequisite for constructing the SARIMA model. Then, by investigating the autocorrelation function (ACF) and the partial autocorrelation function (PACF) plots of the differenced series, it appeared that the orders of p, P, q, and Q might be between 0 and 2 in that the correlation coefficients showed a rapidly decayed trend before lag 3 (Figure S3). Given that some models with simple parameters are not capable of effectively capturing the dynamic dependence structure of the HFMD incidence series, thereby eight plausible models were selected in our experiments (Table 1). Of them, it seemed that the SARIMA(1,0,1)(0,1,1)₁₂ specification was expected to be the preferred model as it provided the minimum information criteria (AIC=17.147 and BIC=27.876) among the plausible models. Further diagnostic checking for the residual series of the best SARIMA model indicated that the key parameters were statistical significance (Table 2), and the residuals from this best SARIMA model behaved like a white noise series as the Box-Ljung Q test showed a p-value greater than 0.05 (Table 3 and Figure 2). Besides, no ARCH effect was observed in the residuals since the LM-test indicated a p-value greater than 0.05 (Table 3). All the results above mean that the SARIMA model obtained is very adequate for modelling the HFMD incidence. Similarly, following the modelling procedures described above, the observations between January 2009 and December 2018 were used to construct the SARIMA model to validate the robustness and accuracy of the SARIMA model. The SARIMA(1,0,1)(0,1,1)₁₂ specification was still considered the optimal model and the statistical tests for this preferred SARIMA model are given in Tables S2-S4 and Figures S4-S6.

Table 1 The Identified Eight Possible SARIMA Methods and Their Corresponding Information Criteria

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

Table 3 Box-Ljung Q and LM Tests of the Residual Series from the Optimal SARIMA (1,0,1)(0,1,1)12 and TBATS(0.062, {1,3}, 0.86, {<12,4>}) Methods

Figure 2 Diagnostic checking for the residual series from the best SARIMA (1,0,1)(0,1,1)12 method. (A) Correlogram of the sample autocorrelation function (ACF); (B) Correlogram of the sample partial autocorrelation function (PACF); (C) p values for the Ljung-Box test. The plot above showed that almost all the sample autocorrelations for the residual series fail to touch the significance bounds apart from the one at lag 25 (which is also reasonable as higher-order autocorrelation may exceed the 95% significance bounds by chance alone) and p values at different lags are greater than 0.05 under the Ljung-Box statistic, suggesting that there is little evidence of non-zero autocorrelations in the residual series at various lags.

Developing the TBATS Model

The forecasting accuracy level of the TBATS model relies on the number of harmonics k used for the seasonal component.²¹ In choosing an effective TBATS model, given one seasonal trait at a time, we thus developed the TBATS model based on the training data from January 2009 to December 2017 by gradually adjusting the k each time but remaining all other harmonics constant for each seasonal component. Meanwhile, an automatic time series algorithm of the ARIMA model was employed to the residual series produced by the TBATS model without the ARMA component to discover the suitable orders of p and q.³⁹ By doing so, the one that had the minimum AIC value was identified as the optimal model. In an effort by trial and error, it is found that the TBATS (0.062, {1,3}, 0.86, {<12,4>}) specification minimizes the AIC (705.260), and the resulting smoothing parameters and other key parameters are listed in Table S5 and the extracted components with Box-Cox transformation based on the above-identified parameters are visible in Figure S7. Further analysis for the residual series indicated that the errors belonged to a white noise series and the ARCH effect from the original observed series was largely attenuated in that gradually increasing the lag periods generated a p-value greater than 0.05 under the Box-Ljung Q test and the LM test (Table 3 and Figure 3). Thereby, we are confident that this best TBATS model is sufficient for the morbidity series forecasting of HFMD. Likewise, following the model selection steps, the TBATS (0.024, {1,1}, 0.855, {<12,4>}) specification was selected as the best-fitting model based on the data from January 2009 to December 2018 to account for the model uncertainty, and the resulting key parameters and statistical checking results are provided in Tables S4-S5 and Figures S8-S9.

Figure 3 Diagnostic checking for the residual series from the best TBATS (0.062, {1,3}, 0.86, {<12,4>}) method. (A) Correlogram of the sample autocorrelation function (ACF); (B) Correlogram of the sample partial autocorrelation function (PACF); (C) p values for the Ljung-Box test. The plot above showed that almost all the sample autocorrelations for the residual series fail to touch the significance bounds and p values at different lags are greater than 0.05 under the Ljung-Box statistic, suggesting that there is little evidence of non-zero autocorrelations in the residual series at various lags.

Comparing the Forecasting Ability Between Models

The forecasting ability between the TBATS model and the SARIMA model was compared from two aspects including in-sample fitting horizon and out-of-sample forecasting horizon based on the measurement metrics of MAD, RMSE, MAPE, and MER. Table 4 and Figure S10 detail the results on the comparison, and it is apparent from the results that the best TBATS models produced lower error rates in both the training and testing segments compared to the best SARIMA models, be it in 24-step or 12-step ahead projections (Table S6 and S7), suggesting that the TBATS model is more reliable and robust in forecasting the HFMD incidence in China than the most frequent use of SARIMA model (Figure S11). Accordingly, we again built the TBATS model based on the 11 years of data to predict the HFMD incidence into the next 24 months. The model selection procedures resulted in the best TBATS (0.022, {3, 1}, 0.851, {<12, 4>}) method, and the resulting forecasting results are presented in Table 5 and Figure 4, being suggestive of a relatively stable epidemiological trend in the HFMD incidence between January 2020 and December 2021 in China.

Table 4 Comparisons of the Fitted Parts and the Predicted Parts Between SARIMA Methods and TBATS Methods

Table 5 Forecasts and Their 95% Uncertainty Limits of the HFMD Incidence for the Next 24 Months Based on the TBATS (0.022, {3, 1}, 0.851, {<12, 4>}) Method

Figure 4 Time series plot showing the in-sample simulation and out-of-sample forecasting using the best SARIMA and TBATS methods. (A) The in-sample fitting and out-of-sample forecasting results using the best SARIMA method; (B) The in-sample fitting and out-of-sample forecasting results using the best TBATS method; (C) The next 24-month projections using the best TBATS method built with the incidence data from January 2009 to December 2019.

Discussion

Since 2008, HFMD has become a major public health concern because of its high morbidity and adverse health effect on more than two million children per year in China.^12,13,22 The use of statistical models with high accuracy and robustness as an effective and helpful tool for the estimation of the long-term epidemic trends of infectious diseases is essential to inform prevention initiatives.¹⁶ For this reason, prompting different researches using various mathematical techniques with the same purpose of achieving an accurate forecast of the upcoming HFMD epidemics is being done. However, as far as we know that there is no study by far using an advanced exponential smoothing state space TBATS framework to describe and estimate the temporal patterns of HFMD. The present research explores, for the first time, its potential of the TBATS model for the application in describing and analyzing the epidemic trajectories of HFMD, and the most commonly used SARIMA model was also fitted to the HFMD incidence data to compare their predictive ability between the TBATS model and the SARIMA model. The results to emerge from this time series analysis indicated that the best TBATS model shows a performance improvement over the optimal SARIMA model in both in-sample data and out-of-sample data, particularly in the forecasting aspect, be it in the 12-step ahead (short-term) projection or the 24-step ahead (long-term) projection. In the meantime, we also developed the TBATS and SARIMA models based on the data from January 2009 to December 2016 to investigate their forecasting ability into the next 36 months, suggesting that the TBATS model was still superior to the SARIMA model (Tables S8 and S9), but these two models produced a high error rate (MAPE=0.441 in the TBATS model vs MAPE=0.930 in the SARIMA model), particularly for the SARIMA model, it had been inadequate for predicting the temporal patterns into the future three years as the MAPE value was greater than 0.5 on the testing data, which is deemed as an unacceptable level of accuracy. This provides further support for the usefulness of the TBATS model as an adequate and effective tool in nowcasting and forecasting the temporal trends of HFMD until the upcoming 24 months in China. Furthermore, an external validation of the utility of TBATS model was performed based on the HFMD incidence data from 2009 to 2017 in Henan Province (Table S10), and the resulting results are listed in Table S11 and Figure S11, indicating a same finding as our present study that the predictive results under TBATS model were closer to the actual data compared with that under SARIMA model. Also, the TBATS model can act as a serviceable tool in estimating the short-term or long-term effects of new prevention and control measures (eg, the introduction of new vaccines available, the optimization of the current prophylaxis, and increasing the intensity of vaccine advocacy and education among children and their parents at high-risk seasons^2,24). If this model produces a significant overestimation over the real data, suggesting that the present interventions may take effect; otherwise, additional interventions are needed to contain the epidemics of HFMD. Besides, the SARIMA model was also discovered to report an acceptable projected result although it yielded a higher forecasting error rate compared to the TBATS model. This further corroborates the finding of a great deal of previous work that the SARIMA model is flexible and useful in the forecasting domain of time series.^{27,29,30,33,40}

The versatile SARIMA method is deemed as the most frequently adopted tool in the forecasting domain of time series data with remarkable seasonality and cyclicity in that this model has the advantage that there is no need to make a foregoing assumption on the inherent rule of a time series.^30,41 For example, Tian et al built a SARIMA (1,1,2) (0,1,1)₁₂ model to describe and analyze the temporal patterns of HFMD between May 2008 and August 2018.¹³ Shi et al investigated the spreading mode of Hemorrhagic Fever with Renal Syndrome morbidity from 2000 to 2017 in East China by establishing a SARIMA (0,1,1,)(0,1,1)₁₂ method.²⁷ In spite of its satisfactory predictive performance of the SARIMA model, it has a limited ability to solve nonlinear problems and to make a long-term forecast as this model is designed under the assumption of linearity.³³ Also, as evidenced by the findings to emerge from our study and other publications,^13,21,42 the SARIMA model is not well suited to address the complex time series with multiple seasonal periods, high-frequency seasonality, non-integer seasonality, and dual-calendar effects despite its wide application.²¹ Instead, the TBATS model is tailored for not only dealing with the complex time series described above but allowing for some types of non-linearity using Box-Cox transformations,^21,43 which enable this model to have the potential to describe the long-term epidemic patterns of HFMD. Given its attractive applications of the TBATS model,²¹ it appeared therefore that this model can be adopted to analyze the HFMD incidence in other countries or regions, even for all the time series, and yet additional work on the current topic is still required to verify its versatility. Additionally, of note, current studies have demonstrated the potential of some new models to perform the long-term forecasting and to handle the nonlinear issue, such as ETS model,⁴⁴ long short-term memory neural network (LSTM),⁹ and neural network nonlinear autoregression (NNAR) model.⁴⁵ Consequently, further investigations into the predictive performance comparison between the TBATS model and the above models are needed.

Understanding the epidemiological features of infectious diseases is very useful and crucial for their prevention and control.⁹ In this study, there were reported HFMD cases all the year round in China, and characterized by biennial or triennial larger epidemics in the HFMD incidence. Such an epidemiological feature fits well with earlier reports in mainland China,^12,13,46 also coincident in Taiwan of China,⁴⁷ Singapore,⁴⁸ Vietnam,¹¹ and Malaysia.⁴⁹ The possible reasons underlying such an epidemic pattern are first, most people who have been infected with the pathogens of HFMD can develop protective antibodies against the causative agents, leading to a risk reduction in the infection with the same types of pathogens.¹³ Second, recent work found that the duration of cross-protection infected with the most common two pathogens of CVA16 and EV71 was about 9.95 weeks (95% uncertainty limit 3.31 to 23.40) in 68% of the population, suggesting that such a periodic pattern of HFMD may be related to cross-serotype immunity.^12,50 Besides, the additional reasons for the periodic pattern described above are expected to go on. Whilst we also observed that the HFMD incidence exhibited a notable seasonality and dual peak distribution, there were few cases between January and March, a rapid upsurge between April and July, then a reduction in August, and yet a secondary peak in September per year. Such a seasonal behavior was in line with the finding from Singapore,⁴⁸ Malaysia,⁴⁹ and most provinces or cities of China (eg, Taiwan, Hong Kong, Ningbo, Shenzhen, Zunyi, Shandong, and Guangdong).²² For the double peak pattern of HFMD, current evidence indicated that it was predominantly seen in southern China,¹² in good agreement with the data reported in Hong Kong and Taiwan of China,^47,51 and Vietnam,¹¹ whereas inconsistent with the data reported in most regions of northern China,¹² where a single peak was predominately seen, this concurs well with the finding from Japan and Malaysia.^12,49 The peak activity of HFMD varied by geography, which may be ascribed to the different pathogens,^12,49 meteorological factors (such as temperature, atmospheric pressure, precipitation, sunshine hours, humidity, and wind),^12,15 and socio-economic differences.^12,22 Prior published work discovered that although the leading pathogens of EV71 and CVA16 arise per year in China (which is different from that in Japan and Malaysia where the leading pathogens arise every 3 or 4 years^12,49), they are circulated in different peaks.⁵² The EV71 is mainly associated with the first strong peak behavior, whereas the CVA16 is predominantly associated with the second weak peak activity, and yet these two common pathogens are less circulated in January and February annually, this wholly corroborates with our initial finding that a trough in the HFMD incidence seen between January and March. Another contributory factor correlated with the least reported HFMD cases in February is closely related to the “spring festival effect”, several studies have given well-grounded discussions on how this effect influences the epidemics of infectious diseases.^53,54 Further, a projection into the next two years based on the TBATS model showed that albeit there did not seem to have a large amplitude of oscillations compared with the data from 2009 to 2019, the HFMD incidence was still high, and the predicted seasonal patterns match well with our findings described above. Currently, some emerging Enteroviruses, especially Coxsackievirus A6 (CA6),^3,14 have been increasing triggers causing HFMD outbreaks, and thus it is a matter of great urgency that the introduction of a multivalent vaccine fusing EV71, CA16 with CA6 to prevent HFMD in addition to the preventions strongly recommended by Saguilet al.²⁴

This study sets out to investigate the flexibility and usefulness of the TBATS model in fitting and forecasting the epidemic patterns of HFMD in China. Importantly, the results to emerge from a series of experiments are encouraging. However, this study also suffers from some drawbacks. First, since HFMD is a self-limiting disease,² the patients with mild and moderate symptoms may not seek medical assistance, thus leading to under-reporting and under-diagnosis. Second, daily or weekly data may provide a greater investigation into the temporal discrepancies among years. Nonetheless, it is difficult to obtain such data for shorter periods. Third, a high forecasting performance can be achieved when using the TBATS model to nowcast and forecast the temporal levels of the next two years (24 months). Future work should therefore include the new data into the model to ensure its high forecasting accuracy. Lastly, the present results may not be transferable to the conclusion that the TBATS model has superiority in the prediction of HFMD in other regions or other infectious diseases, and additional work is entailed.

Conclusion

The most striking result to emerge from the data is that the TBATS model has the capacity to outperform the most frequently used SARIMA model in decomposing and forecasting the HFMD incidence in China. This model can be recommended as a flexible and useful tool in the decision-making process of HFMD prevention and control. Moreover, the predicted HFMD incidence for the next two years remains still high, considering the new increasing emerging pathogens causing HFMD outbreaks, it is required to curb and harness this urgent issue within the feasible measures taken.