Volume 35 Issue 6
Jun.  2022
Turn off MathJax
Article Contents

WANG Yong Bin, LI Yan Yan, LU Hao, TAO Ying Jun, LI Yu Hong, WANG Lei, LIANG Wen Juan. Now- and Fore-casting the Secular Epidemiological Trends and Seasonality of the Comeback of Scarlet Fever in China: A 16-year Time Series Analysis[J]. Biomedical and Environmental Sciences, 2022, 35(6): 563-567. doi: 10.3967/bes2022.076
Citation: WANG Yong Bin, LI Yan Yan, LU Hao, TAO Ying Jun, LI Yu Hong, WANG Lei, LIANG Wen Juan. Now- and Fore-casting the Secular Epidemiological Trends and Seasonality of the Comeback of Scarlet Fever in China: A 16-year Time Series Analysis[J]. Biomedical and Environmental Sciences, 2022, 35(6): 563-567. doi: 10.3967/bes2022.076

Now- and Fore-casting the Secular Epidemiological Trends and Seasonality of the Comeback of Scarlet Fever in China: A 16-year Time Series Analysis

doi: 10.3967/bes2022.076
Funds:  This work was supported by the Natural Science Foundation of Henan, China, Innovation and Entrepreneurship Training Project for University Students of Henan Province and Xinxiang Medical University [222300420265, S202110472047, S202010472007, and XYXSKYZ201932] and the Key Scientific Research Project of Universities in Henan [21A330004]
More Information
  • Author Bio:

    WANG Yong Bin, male, born in 1989, Professor, MD, majoring in infectious disease epidemiology;

    LI Yan Yan, female, born in 1996, master’s, majoring in infectious disease epidemiology; LU Hao, male, born in 1999, Undergraduate, majoring in preventive medicine

  • Corresponding author: LIANG Wen Juan, MD, Tel: 86-373-3831646, E-mail: wenwen3_1@126.com
  • &These authors contributed equally to this work.
  • Received Date: 2021-11-12
  • Accepted Date: 2022-06-02
  • &These authors contributed equally to this work.
  • 加载中
  • [1] Chen M, Cai J, Davies MR, et al. Increase of emm1 isolates among group A Streptococcus strains causing scarlet fever in Shanghai, China. Int J Infect Dis, 2020; 98, 305−14.
    [2] Li WT, Feng RH, Li T, et al. Spatial-temporal analysis and visualization of scarlet fever in mainland China from 2004 to 2017. Geospat Health, 2020; 15, 831.
    [3] Liu YH, Chan TC, Yap LW, et al. Resurgence of scarlet fever in China: a 13-year population-based surveillance study. Lancet Infect Dis, 2018; 18, 903−12. doi:  10.1016/S1473-3099(18)30231-7
    [4] Wang YB, Xu CJ, Yao SQ, et al. Estimating the prevalence and mortality of coronavirus disease 2019 (COVID-19) in the USA, the UK, Russia, and India. Infect Drug Resist, 2020; 13, 3335−50. doi:  10.2147/IDR.S265292
    [5] De Livera AM, Hyndman RJ, Snyder RD. Forecasting time series with complex seasonal patterns using exponential smoothing. J Am Stat Assoc, 2011; 106, 1513−27. doi:  10.1198/jasa.2011.tm09771
    [6] Armitage P, Berry G, Matthews JNS. Statistical methods in medical research. 4th ed. Blackwell Publishing. 2002, 127.
    [7] Lee CF, Cowling BJ, Lau EHY. Epidemiology of reemerging scarlet fever, Hong Kong, 2005-2015. Emerg Infect Dis, 2017; 23, 1707−10. doi:  10.3201/eid2310.161456
    [8] Lamagni T, Guy R, Chand M, et al. Resurgence of scarlet fever in England, 2014-16: a population-based surveillance study. Lancet Infect Dis, 2018; 18, 180−7. doi:  10.1016/S1473-3099(17)30693-X
    [9] Liu YH, Ding H, Chang ST, et al. Exposure to air pollution and scarlet fever resurgence in China: a six-year surveillance study. Nat Commun, 2020; 11, 4229. doi:  10.1038/s41467-020-17987-8
    [10] Hyndman R, Athanasopoulos G, Bergmeir C, et al. Forecast: forecasting functions for time series and linear models. https://cran.r-project.org/web/packages/forecast/. [2022-05-10].
  • 21432Supplementary Materials.pdf
  • 加载中
通讯作者: 陈斌, bchen63@163.com
  • 1. 

    沈阳化工大学材料科学与工程学院 沈阳 110142

  1. 本站搜索
  2. 百度学术搜索
  3. 万方数据库搜索
  4. CNKI搜索

Figures(17)  / Tables(13)

Article Metrics

Article views(892) PDF downloads(79) Cited by()

Proportional views
Related

Now- and Fore-casting the Secular Epidemiological Trends and Seasonality of the Comeback of Scarlet Fever in China: A 16-year Time Series Analysis

doi: 10.3967/bes2022.076
Funds:  This work was supported by the Natural Science Foundation of Henan, China, Innovation and Entrepreneurship Training Project for University Students of Henan Province and Xinxiang Medical University [222300420265, S202110472047, S202010472007, and XYXSKYZ201932] and the Key Scientific Research Project of Universities in Henan [21A330004]
  • Author Bio:

  • Corresponding author: LIANG Wen Juan, MD, Tel: 86-373-3831646, E-mail: wenwen3_1@126.com
  • &These authors contributed equally to this work.
&These authors contributed equally to this work.
WANG Yong Bin, LI Yan Yan, LU Hao, TAO Ying Jun, LI Yu Hong, WANG Lei, LIANG Wen Juan. Now- and Fore-casting the Secular Epidemiological Trends and Seasonality of the Comeback of Scarlet Fever in China: A 16-year Time Series Analysis[J]. Biomedical and Environmental Sciences, 2022, 35(6): 563-567. doi: 10.3967/bes2022.076
Citation: WANG Yong Bin, LI Yan Yan, LU Hao, TAO Ying Jun, LI Yu Hong, WANG Lei, LIANG Wen Juan. Now- and Fore-casting the Secular Epidemiological Trends and Seasonality of the Comeback of Scarlet Fever in China: A 16-year Time Series Analysis[J]. Biomedical and Environmental Sciences, 2022, 35(6): 563-567. doi: 10.3967/bes2022.076
  • Scarlet fever (SF) is a common communicable disease that results from group A Streptococcus (GAS) infections[1]. SF accounted for the global loss of life among children 5–15 years of age in the 18th and 19th centuries[2]. A rapid reduction in SF morbidity and mortality occurred due to the scale-up of effective antibiotics and improvements in sanitation and nutrition[3]. The unexpected increase in the incidence of SF has attracted a renewed interest in infectious diseases[3]. Because the triggers that cause SF outbreaks are not fully understood and there is a scarcity of available vaccines protecting susceptible populations from GAS infections, effective prevention and control plans are required to stop the continued spread of SF.

    Time series analysis assists in the development of hypotheses to explain the temporal patterns of different diseases and to analyze the spread, therefore, facilitating the creation of a quality forecasting system. The seasonal autoregressive integrated moving average (SARIMA) model has been widely applied to estimate the epidemiological patterns of contagious diseases because this model has a simple structure, fast applicability, and a relatively high forecasting reliability level [4]. It has been shown that the SARIMA model is able to satisfactorily estimate a simple time series [4], but it is difficult to manage complex time series, such as the data with multiple seasonal periods, high-frequency seasonality, non-integer seasonality, and dual-calendar effects. By comparison, the innovation state-space modelling framework that combines Box-Cox transformations, Fourier series with time-varying coefficients, and autoregressive moving average (ARMA) error correction (known as the TBATS method) is customized for use with the patterns included in a complex time series described above [5]. In addition, the TBATS model is used for linearity and some types of non-linearity in a complex time series based on Box-Cox transformations [5], which makes it possible to perform a multistep ahead prediction. Moreover, the TBATS model is able to decompose a complex seasonal time series into its trend, seasonal, and irregular components [5], which is not able to be undertaken by use of the SARIMA model. Importantly, SF morbidity has been shown to display dual seasonal patterns in some countries. Therefore, this study analyzed the long-term epidemic patterns using the TBATS model. The forecasting power under the TBATS model was compared with the SARIMA model.

    We obtained the monthly SF incidence and population data between January 2004 and December 2019 from the Chinese CDC and the Statistical Yearbook of China, respectively. Then, we partitioned the SF morbidity series into two segments comprising a training dataset from January 2004 to December 2017 to construct the SARIMA and TBATS models and a testing dataset from January 2018 to December 2019 to test the generalization of both models. Two additional datasets were provided to test the robustness of both models: the first 180 data sets from January 2004 to December 2018 and 156 data sets from January 2004 to December 2016 were treated as training datasets, respectively; and the remainder were testing datasets. We estimated the changing epidemiologic trends of SF based on annual percentage change (APC) and average APC (AAPC) using the Joinpoint regression program (version 4.8.0.1). We constructed the TBATS model with the “forecast,” “tseries,”, and “FinTS” packages in R software (version 3.4.3). The incidence rate ratio (IRR) with a 95% uncertainty limit (UL) before and after the SF outbreak was calculated using the method proposed by Armitage and Berry [6]. The mean absolute deviation (MAD), root mean square error (RMSE), mean absolute percentage error (MAPE), mean error rate (MER), and root mean square percentage error (RMSPE) were computed to compare the forecasting ability of both models.

    During the study period, a notable increase in SF morbidity was detected, with an AAPC = 8.942 (95% UL: 5.995–11.971; t = 6.697, P < 0.001; Supplementary Figure S1, available in www.besjournal.com), and the highest morbidity (5.930/100,000 persons) occurred in 2019, which increased by a factor of nearly four compared with the lowest level (1.489/100,000 population) in 2004 (Supplementary Figure S2, available in www.besjournal.com). The SF epidemics remained relatively steady from 2004–2010 (average, 1.937/100,000 persons annually), with an AAPC = −0.840 (95% UL: −7.678 to 6.505; t = −0.231, P = 0.817). An unexpected outbreak was witnessed in 2011, and since then a rapidly increasing trend occurred (average, 4.580/100,000 persons annually [excluding 2013]; the exact causes regarding this annual drop are unknown; Supplementary Figure S2), with an AAPC = 5.952 (95% UL: 0.239–11.992; t = 2.466, P = 0.043), which showed good agreement with a resurgence of SF in Hong Kong, China[7], but inconsistent with England where the resurgence occurred in 2014 [8]. There has been a doubling in the SF incidence during the post-resurgence periods compared with the pre-resurgence periods (IRR = 2.364, 95% UL: 2.358–2.370). Nevertheless, the driving force associated with the increased pathogenicity of GAS fails to be elucidated. A possible explanation may be due to the acquisition of novel prophages harboring new hybridizations of toxin genes and antimicrobial resistance genes, which is related to the emergence and expansion of the predominant genotypes of emm12 and emm1 in China [3]. Another explanation may be associated with the natural periodicity of the SF incidence (SF epidemics are characterized by a cyclic change of approximately 6 years [3]). A third explanation may be linked to the relaxation of the 2-child policy in 2011 [3], which led to an increase in the number of susceptible individuals. A fourth reason may be attributed to improvements in the diagnostic capacity and the increased awareness of medical workers in reporting SF [3]. A fifth reason may be a result of the deterioration of air quality in China [9], despite the gradual improvements in the last 2 years. Finally, there are no vaccines available to prevent infections with GAS until now.

    Figure S1.  Joinpoint regression showing the scarlet fever epidemic trends in the period 2004–2019. *Showed that the annual percent change (APC) is significantly different from zero at the significance level. This plot was plotted using Joinpoint regression program (Version 4.8.0.1). Joinpoint is a statistical tool for the description and analysis of trends with the Joinpoint methods. This Joinpoint Regression Program is designed as a Windows-based statistical software package to analyze Joinpoint methods. Users can use this software to investigate whether there is a statistically notable change in trend (Available from https://surveillance.cancer.gov/joinpoint/. Accessed on 16 May 2022).

    Figure S2.  The yearly incidence of SF from 2004 to 2019. It was seen that SF had the lowest incidence in 2004 (1.489 per 100,000 population) and the highest morbidity in 2019 (5.930 per 100,000 persons).

    A marked semi-annual seasonal behavior occurred in the monthly SF incidence, with a strong peak between May and June, and a weak peak between November and December (Supplementary Figures S3S4, available in www.besjournal.com). We surmised that different climatic features and beginning of spring and autumn semesters contributed to this difference in the margin of peak activities. Our seasonal profile correlates well with previous findings from Hong Kong, China [7]; however, discordant with that in England, which peaked between February and March [8]. This inconsistency may be due to the different school breaks, population density, and different GAS emm gene types in east Asia and Europe [1, 8]. In addition, the SF epidemics retain the lowest level in February every year (Supplementary Figure S3), attributable to the winter holidays and the Spring Festival.

    Figure S3.  Time series plot showing the monthly SF incidence rate and the decomposed seasonal, trend, and irregular components based on the STL technique. (A) Original monthly SF morbidity series. (B) Decomposed trend component. (C) Decomposed seasonal component. (D) Decomposed random component. As illustrated, overall the SF morbidity showed a rising tendency, with a sudden increase in 2011, and having notable seasonal patterns.

    Figure S4.  The decomposed seasonal factors for the scarlet fever morbidity series using the multiplicative decomposition method (The classical Seasonal Decomposition technique can be applied to decompose a time series into a seasonal pattern, a trend component, and a random fluctuation component. This technique is implemented based on the ratio-to-moving-average method. The Seasonal Decomposition technique can be divided into two types including a multiplicative and an additive method. For a multiplicative decomposition method, the seasonal component is a factor by which the seasonally adjusted series is multiplied to produce original observed values. For a time series without seasonal variation, the seasonal pattern is deemed as one. For an additive decomposition method, the seasonal adjustments are added to the seasonally adjusted series to yield the original series. For the morbidity series of infectious diseases, the multiplicative decomposition method typically produces a more reliable decomposed result in application. By employing a multiplicative decomposition method to an object series, the seasonal pattern or seasonal factors included in the object series can be obtained. Such a decomposition is very helpful in removing the seasonal effect from a series so that we can visibly look into the other characteristics of interest in the series that may be masked by the seasonal trait[7,8]).

    The forecasts under the TBATS approach rely largely on the number of harmonics ki applied for each seasonal pattern. As a result, in selecting the number of harmonics ki, considering one seasonal component each time, we then fitted the model on the target data repeatedly via gradually increasing the number of harmonics ki, but holding the remaining harmonics constant for each i until the optimal AIC is obtained. In determining the most suitable orders (p and q) of the ARMA model, we used the automatic procedure proposed by Hyndman and colleagues[10] to fit the forecasting residuals. If the selected model with the ARMA (p, q) residual component generates a smaller AIC than the one model without the ARMA (p, q) residual component, this selected specification would be considered as the best possible model; otherwise, the ARMA (p, q) residual component is deleted. After modelling by trial and error, the TBATS (0.04, {4,0}, 0.882, {<12,5>}) specification was selected as the optimal model in that the minimum AIC (−197.965) was detected in this model, and the identified key parameters of this best TBATS model are reported in Supplementary Table S1, available in www.besjournal.com. Additional statistical diagnoses for the forecast errors are provided in Supplementary Table S2 and Supplementary Figures S5S6, available in www.besjournal.com. The Ljung-Box Q statistics of the forecast errors produced a Q(18) = 11.442 with a P-value of 0.875, indicating no serial correlations in this residual series. Moreover, the ARCH effect was largely removed because the LM(18) = 23.808 with a P-value of 0.302. These results confirmed the adequacy of the model specifications. Similarly, based on the modelling steps described above, the TBATS (0.01, {0,0}, 0.898, {<12,5>}) and TBATS (0.048, {0,0}, 0.902, {<12,5>}) specifications tended to be the preferred models for forecasting the 12- and 36-holdout periods (Supplementary Tables S3S4 and Supplementary Figures S7S10, available in www.besjournal.com).

    Parameters24-step ahead forecasting12-step ahead forecasting36-step ahead forecasting
    TBATS (0.04, {4,0}, 0.882,
    {<12,5>}) model
    TBATS (0.01, {0,0}, 0.898,
    {<12,5>}) model
    TBATS (0.048, {0,0}, 0.902,
    {<12,5>}) model
    Lambda0.04010.01050.0480
    Alpha0.16520.64890.6795
    Beta−0.0252−0.0362−0.0432
    Damping Parameter0.88180.89820.9024
    Gamma-1 Values0.0002−0.0046−0.0048
    Gamma-2 Values−1.26820.00500.0072
    AR coefficientsAR1 = 0.4772, AR2 = 0.2318,
    AR3 = −0.1894, AR4 = 0.3222
    Seed states
    1−2.3530−2.4373−2.3215
    20.10890.11300.1015
    3−0.0246−0.0185−0.0251
    4−0.0748−0.0757−0.0731
    5−0.0069−0.0077−0.0065
    60.02130.02660.0207
    70.02960.03120.0279
    8−0.15470.14890.1353
    9−0.5953−0.6307−0.5756
    10−0.0602−0.0649−0.0598
    110.1463−0.1521−0.1404
    12−0.0231−0.0248−0.0218
    130.0000
    140.0000
    150.0000
    160.0000
    Sigma0.16860.19090.1797
    AIC−197.9650−175.3322−198.4532
    Parameters48-step ahead forecasting72-step ahead forecastingFuture 36-step ahead forecasting
    TBATS (0.037, {0,0}, 0.901,{<12,5>})TBATS (0.029, {0,0}, 0.87,{<12,4>})TBATS (0.023, {0,0}, 0.895, {<12,5>})
    Lambda0.03680.02920.0225
    Alpha0.69000.72530.6444
    Beta−0.0410−0.1128−0.0378
    Damping Parameter0.90130.87020.8955
    Gamma-1 Values−0.0069−0.0115−0.0050
    Gamma-2 Values0.00520.00350.0037
    AR coefficients

    Table S1.  The best TBATS models obtained on the basis of different datasets and their key parameters

    LagsSARIMA modelTBATS model
    Box-Ljung QPLM-testPBox-Ljung QPLM-testP
    10.0440.8330.0010.9790.0050.94352.651< 0.001
    31.9440.5842.3640.5000.2350.97222.234< 0.001
    63.7200.7154.2230.6471.2460.97523.5550.001
    96.0930.7315.9810.7424.5380.87323.6240.005
    1213.0210.3689.4270.6667.3200.83623.4160.024
    1513.0940.59512.3790.65010.4560.79023.5330.073
    1813.6630.75119.8720.34011.4420.87524.5500.138
    2117.3130.69220.2620.50519.0310.58323.8080.302
    2423.2370.50629.3520.20723.0610.51624.0150.461
    2726.0600.51525.4920.54723.4660.66023.9240.635
    3026.0820.67126.4220.65325.2300.71424.9180.729
    3329.1330.66029.8630.62429.0040.66724.9410.842
    3641.4170.24637.2300.41230.7070.71824.3910.929
      Note. SARIMA seasonal autoregressive integrated moving average method, TBATS an advanced innovation state-space modelling framework by combining Box-Cox transformations, Fourier series with time-varying coefficients, and autoregressive moving average (ARMA) error correction.

    Table S2.  Box-Ljung Q and LM tests for the residual series of the SARIMA and TBATS models

    Figure S5.  Estimated autocorrelation function (ACF) and partial ACF (PACF) plots to display the statistical checking for the residual series of the SARIMA and TBATS models. (A) ACF plot for the SARIMA model. (B) PACF plot for the SARIMA model. (C) ACF plot for the TBATS model. (D) PACF plot for the TBATS model. We can see that the correlation coefficients at lags 11, 23, and 36 in Figure 2A and 2B, and at lag 19 in Figure 2C and 2D touched the confidence levels (this is reasonable as the correlation coefficients at higher-order may exceed the significance bounds by chance alone), but all other correlation coefficients at different lags failed to touch the confidence levels, seemingly suggesting effective and sufficient models.

    Figure S6.  The extracted level, slope, and seasonal components of the best TBATS method developed on the basis of the data from January 2004 to December 2017.

    LagsSARIMA (0,1,0)(0,1,1)12 modelTBATS(0.01, {0,0}, 0.898, {<12,5>}) model
    Box-Ljung QPLM-testPBox-Ljung QPLM-testP
    10.0240.8780.0000.9930.0180.89244.923< 0.001
    32.3530.5032.1960.5339.4280.02417.323< 0.001
    64.5290.6054.2870.63827.713< 0.00118.7420.005
    97.7200.5636.0630.73435.665< 0.00122.3800.008
    1215.2920.2269.7470.63839.404< 0.00122.7150.030
    1515.3490.42713.6510.55248.268< 0.00122.6950.091
    1815.9800.59421.2900.26553.640< 0.00125.4630.113
    2118.9520.58821.5690.42570.249< 0.00124.9020.251
    2424.7340.42030.2920.17582.013< 0.00124.8630.413
    2727.5250.43626.2180.50789.634< 0.00126.3980.497
    3027.5360.59527.5880.59297.136< 0.00127.2680.609
    3329.7880.62831.1920.557107.120< 0.00129.0800.663
    3642.1280.22337.8600.384113.480< 0.00129.8720.754
      Note. SARIMA seasonal autoregressive integrated moving average method, TBATS an advanced innovation state-space modelling framework by combining Box-Cox transformations, Fourier series with time-varying coefficients and autoregressive moving average (ARMA) error correction.

    Table S3.  The Box-Ljung Q and LM tests for the best SARIMA (0,1,0)(0,1,1)12 model and TBATS (0.01, {0,0}, 0.898, {<12,5>}) model built based on the date from January 2004 to December 2018

    LagsSARIMA (0,1,0)(0,1,1)12 modelTBATS (0.048, {0,0}, 0.902, {<12,5>}) model
    Box-Ljung QPLM-testPBox-Ljung QPLM-testP
    10.0550.8150.0010.9750.0230.87941.160< 0.001
    31.9210.5892.0110.5706.7330.08114.0420.003
    63.8350.6993.9510.68322.4880.00115.2020.019
    96.4820.6915.7270.76727.2330.00118.0450.035
    1213.4340.3388.3050.76129.9880.00318.2080.110
    1513.4780.56511.0830.74735.4130.00218.3020.247
    1814.0640.72517.4750.49138.7120.00321.1000.274
    2117.6790.66917.6100.67451.1740.00220.8630.467
    2422.7460.53525.2360.39360.873< 0.00120.9720.640
    2725.1300.56721.9880.73865.748< 0.00122.4460.714
    3025.2700.71222.7990.82370.290< 0.00123.1780.808
    3328.7590.67825.6910.81478.981< 0.00124.3820.861
    3640.7730.26932.0640.65683.260< 0.00124.5380.926
      Note. SARIMA seasonal autoregressive integrated moving average method, TBATS, an advanced innovation state-space modelling framework by combining Box-Cox transformations, Fourier series with time-varying coefficients and autoregressive moving average (ARMA) error correction.

    Table S4.  The Box-Ljung Q and LM tests for the best SARIMA (0,1,0)(0,1,1)12 model and TBATS (0.048, {0,0}, 0.902, {<12,5>}) model built based on the date from January 2004 to December 2016

    Figure S7.  Statistical checking for the residual series of the TBATS model developed on the basis of the data from January 2004 to December 2018. (A) Time series plot showing the residual series. (B) Autocorrelation function (ACF) plot. (C) Partial autocorrelation function (PACF) plot. The correlogram shows that although there are some residual ACF and PACF for the forecast errors that exceed the significance bounds, this is the best model used to estimate the upcoming 12-month epidemics of scarlet fever.

    Figure S10.  The extracted level, slope, and seasonal components of the best TBATS method developed on the basis of the data from January 2004 to December 2016.

    Similarly, following the SARIMA modelling steps, the optimal SARIMA models on different datasets were identified (Supplementary Tables S2S5, and Supplementary Figures S5, S11S12, available in www.besjournal.com). Subsequently, the best SARIMA and TBATS models could be used to perform multistep ahead predictions (Figure 1 and Supplementary Figures S13S14, available in www.besjournal.com). Table 1 lists the measurement metrics, which indicate the forecasting reliability levels on different time windows under the preferred SARIMA and TBATS models. The optimal TBATS models provided a smaller MAD, MAPE, RMSE, RMSPE, and MER compared with the optimal SARIMA models, with a performance improvement of almost 50% in the forecasting abilities for estimating both short- and long-term epidemiological trends, albeit the predictive potential showed a slight reduction with the increase in prediction time windows. We further compared the forecasting abilities of both methods for 48- and 72-step ahead predictions, and the comparative results are listed in Supplementary Tables S1, S6S7, and Supplementary Figures S15S16 (available in www.besjournal.com), which show a similar finding, but the predictive results deviated from the epidemic trajectories. In addition, we used the SF incidence data in Liaoning, Heilongjiang, and Shandong provinces, and Inner Mongolia (which are the hardest hit areas by SF in China in the last decade [2, 3]) to assess the predictive quality of these two methods. Likewise, the TBATS method produced lower error rates in all the datasets (Supplementary Table S8, available in www.besjournal.com). Our recent study indicated that the Error-Trend-Seasonal (ETS) model also has a powerful potential in estimating the long-term epidemic behaviors of diseases [4]. As a result, we further developed the ETS model based on the SF morbidity to predict the epidemiological trends, and the results also showed similar findings (the computed MAPE values were 16.101% vs. 38.511%, 21.142% vs. 28.273%, and 23.984% vs. 26.735% in the 12-, 24-, and 36-step ahead forecasts, respectively; Supplementary Table S9, available in www.besjournal.com). These findings further substantiated the utility of the TBATS model. The TBATS model was introduced by adding the trigonometric representation of seasonal components based on the Fourier series into the traditional BATS model, which enabled handling of all complex time series, as well as linear and non-linear information [5], thus indicating the suitability and adequacy of this model. Considering the attractive advantages of the TBATS model, this model can be recommended as a flexible and useful long-term predictive tool in assessing the epidemic patterns of SF in other countries or other contagious diseases; however, further work is required for validation. Moreover, with the rapid advances in the forecasting domain of time series, many hybrid prediction models (e.g., SARIMA-BPNN, SARIMA-GRNN, and SARIMA-LSTM) have also been reported to show an attractive advantage in estimating the long-term epidemic trajectories of diseases. Therefore, what is now needed are studies involving comparisons of the predictive reliable level between the TBATS model and the above-mentioned models.

    Figure 1.  Comparative results of the forecasts based on the SARIMA and TBATS models. (A) The comparison between the 24-step ahead forecasts of the SARIMA model and the observed values. (B) The predicted upcoming 36-month values from January 2020 to December 2022 using the SARIMA model. (C) The comparison between the 24-step ahead forecasts of the TBATS model and the observed values. (D) The predicted upcoming 36-month values from January 2020 to December 2022 using the TBATS model.

    ModelsTesting horizons
    MADMAPERMSEMERRMSPE
    24-step ahead predictions
     SARIMA0.29564.3460.3820.6070.551
     TBATS0.11421.1420.1600.2350.061
    Reduced percentage (%)
     SARIMA vs. TBATS61.35667.14358.11561.28588.929
    12-step ahead predictions
     SARIMA0.21242.9510.2570.4290.467
     TBATS0.08716.1010.1130.1770.193
    Reduced percentage (%)
     SARIMA vs. TBATS58.96262.51356.03158.74158.672
    36-step ahead predictions
     SARIMA0.25855.5270.3960.5460.758
     TBATS0.13323.9840.1740.2820.271
    Reduced percentage (%)
     SARIMA vs. TBATS48.45056.80756.06148.35264.248
      Note. SARIMA, seasonal autoregressive integrated moving average method; TBATS, an advanced innovation state-space modelling framework by combining Box-Cox transformations, Fourier series with time-varying coefficients and autoregressive moving average (ARMA) error correction; MAD, mean absolute deviation; MAPE, mean absolute percentage error; RMSE, root mean square error; MER, mean error rate; RMSPE, root mean square percentage error.

    Table 1.  The comparisons of the predicted results between the SARIMA model and the TBATS model on different testing datasets

    VariablesEstimateS.E.tPStationary R2R2AICBICLjung-Box Q (18)
    StatisticsP
    SARIMA (0,1,0)(0,1,1)12method developed on the basis of the data from January 2004 to December 2018
    SMA10.7820.07210.873< 0.0010.3000.944−6.159−6.14615.9800.525
    SARIMA (0,1,0)(0,1,1) 12 method developed on the basis of the data from January 2004 to December 2017
    SMA10.7950.07510.600< 0.0010.2920.935−6.166−6.15213.6630.691
    SARIMA (0,1,0)(0,1,1) 12 method developed on the basis of the data from January 2004 to December 2016
    SMA10.7710.0849.230< 0.0010.2880.922−6.162−6.12714.0640.663
      Note. SARIMA seasonal autoregressive integrated moving average method, S.E. standard error, AIC Akaike's Information Criterion, BIC Schwarz's Bayesian Information Criterion, SMA1 seasonal moving average method at lag 1.

    Table S5.  The obtained best SARIMA models with their key parameters on the basis of different datasets

    Figure S13.  Comparative results of the forecasts based on the SARIMA and TBATS models. (A) The comparison between the 12-step ahead forecasts of the SARIMA model and the observed values. (B) The comparison between the 12-step ahead forecasts of the TBATS model and the observed values.

    Figure S14.  Comparative results of the forecasts based on the SARIMA and TBATS models. (A) The comparison between the 36-step ahead forecasts of the SARIMA model and the observed values. (B) The comparison between the 36-step ahead forecasts of the TBATS model and the observed values.

    Figure S15.  Comparative results of the forecasts based on the SARIMA and TBATS models. (A) The comparison between the 48-step ahead forecasts of the SARIMA model and the observed values. (B) The comparison between the 48-step ahead forecasts of the TBATS model and the observed values.

    Figure S16.  Comparative results of the forecasts based on the SARIMA and TBATS models. (A) The comparison between the 72-step ahead forecasts of the SARIMA model and the observed values. (B) The comparison between the 72-step ahead forecasts of the TBATS model and the observed values.

    SitesModelsTesting horizons
    MADMAPERMSEMER
    Liaoning province12-step-ahead forecast
    SARIMA (0,1,0)(0,1,1)120.6818.9830.8580.697
    TBATS (0.073, {4,0}, 0.888, {<12,4>})0.3585.8300.4320.367
    Percentage reductions (%)
    TBATS vs. SARIMA47.43035.10049.65047.346
    24-step-ahead forecast
    SARIMA (0,1,0)(0,1,1)120.33042.4530.5090.326
    TBATS(0.121, {0,0}, 0.868, {<12,4>})0.26225.1440.3720.259
    Percentage reductions (%)
    TBATS vs. SARIMA20.60640.77226.91620.552
    36-step-ahead forecast
    SARIMA (0,1,0)(0,1,1)122.303275.1813.2072.393
    TBATS (0.052, {0,0}, 0.875, {<12,5>})0.22426.4480.3050.233
    Percentage reductions (%)
    TBATS vs. SARIMA90.27490.38990.49090.263
    Heilongjiang province12-step-ahead forecast
    SARIMA (0,1,0)(0,1,1)120.22751.8080.2620.383
    TBATS (0.151, {0,0}, 0.889, {<12,5>})0.13032.6570.1670.219
    Percentage reductions (%)
    TBATS vs. SARIMA42.73136.96536.26042.820
    24-step-ahead forecast
    SARIMA (0,1,0)(0,1,1)120.40274.1760.5620.655
    TBATS (0.17, {0,0}, 0.889, {<12,5>})0.11925.8190.1620.194
    Percentage reductions (%)
    TBATS vs. SARIMA70.39865.19271.17470.382
    36-step-ahead forecast
    SARIMA (0,1,0)(0,1,1)121.448253.0722.1142.357
    TBATS(0.161, {0,0}, 0.894, {<12,4>})0.15520.2350.2520.446
    Percentage reductions (%)
    TBATS vs. SARIMA89.29692.00488.07981.078
    Shandong province12-step-ahead forecast
    SARIMA (0,1,1)(0,1,1)120.23639.8230.2820.319
    TBATS(0, {0,0}, 0.932, {<12,5>})0.17519.3900.2800.237
    Percentage reductions (%)
    TBATS vs. SARIMA25.84751.3100.70925.705
    24-step-ahead forecast
    SARIMA (0,1,1)(0,1,1)120.50976.9300.6700.717
    TBATS (0, {0,0}, 0.932, {<12,5>})0.16819.2220.2680.237
    Percentage reductions (%)
    TBATS vs. SARIMA66.99475.01460.00066.946
    36-step-ahead forecast
    SARIMA (0,1,1)(0,1,1)121.683246.8652.3762.491
    TBATS (0, {0,0}, 0.926, {<12,5>})0.31243.2770.4010.461
    Percentage reductions (%)
    TBATS vs. SARIMA81.46282.46983.12381.493
    Inner Mongolia12-step-ahead forecast
    SARIMA (1,0,2)(0,1,1)120.40151.0260.4960.397
    TBATS (0.232, {0,0}, -, {<12,4>})0.13032.6570.1670.219
    Percentage reductions (%)
    TBATS vs. SARIMA67.58135.99966.33144.836
    24-step-ahead forecast
    SARIMA (1,0,2)(0,1,1)120.31938.9840.4180.308
    TBATS (0.232, {0,0}, -, {<12,4>})0.16616.5080.2360.161
    Percentage reductions (%)
    TBATS vs. SARIMA47.96257.65443.54147.727
    36-step-ahead forecast
    SARIMA (1,0,2)(0,1,1)120.48456.5020.5730.494
    TBATS (0.21, {0,0}, -, {<12,4>})0.15817.0020.2240.161
    Percentage reductions (%)
    TBATS vs. SARIMA67.35569.90960.90867.409
      Note. SARIMA seasonal autoregressive integrated moving average method, TBATS an advanced innovation state-space modelling framework by combining Box-Cox transformations, Fourier series with time-varying coefficients and autoregressive moving average (ARMA) error correction, MAD mean absolute deviation, MAPE mean absolute percentage error, RMSE root mean square error, MER mean error rate.

    Table S8.  The comparisons of the predicted results between the SARIMA model and the TBATS model on different testing datasets

    ModelsTesting power
    MADMAPERMSEMERRMSPE
    12-data ahead forecasting
     SARIMA144.73425.049161.6710.2030.296
     TBATS91.79914.772123.6530.1290.193
     ETS148.63525.835177.0920.2090.325
    Decreased percentage (%)
     TBATS vs. SARIMA36.57441.02823.51636.45334.797
     TBATS vs. ETS38.23942.82230.17638.27840.615
    24-data ahead forecasting
     SARIMA357.73453.753406.2780.4600.625
     TBATS260.61938.407301.4410.3350.469
     ETS279.01444.604307.7310.3590.536
    Decreased percentage (%)
     TBATS vs. SARIMA27.14728.54925.80427.17424.960
     TBATS vs. ETS6.59313.8932.0446.68512.500
    36-data ahead forecasting
     SARIMA290.78242.156365.5700.3340.560
     TBATS204.15128.685254.4180.2350.387
     ETS297.35944.118347.6890.3420.564
    Decreased percentage (%)
     TBATS vs. SARIMA29.79231.95530.40529.64130.893
     TBATS vs. ETS31.34534.98126.82631.28731.383
    60-data ahead forecasting
     SARIMA224.93128.874310.0430.2580.424
     TBATS183.78021.558235.6650.2110.250
     ETS146.43716.449212.0280.1680.223
    Decreased percentage (%)
     TBATS vs. SARIMA18.29525.33823.99018.21741.038
     TBATS vs. ETS-20.319-23.699-10.030-20.379-10.800
    84-data ahead forecasting
     SARIMA931.622134.9151016.8131.0401.679
     TBATS167.43220.081204.1720.1870.240
     ETS183.77721.010227.1130.2050.250
    Decreased percentage (%)
     TBATS vs. SARIMA82.02885.11679.92082.01985.706
     TBATS vs. ETS8.8944.42210.1018.7804.000
    108-data ahead forecasting
     SARIMA1483.868202.7291588.8761.5602.500
     TBATS208.77522.849287.5360.2200.310
     ETS268.72630.177354.0450.2830.386
    Decreased percentage (%)
     TBATS vs. SARIMA85.93088.72981.90385.89787.600
     TBATS vs. ETS22.30924.28318.78522.26119.689
      Note. SARIMA seasonal autoregressive integrated moving average method, TBATS an advanced innovation state-space modeling framework by combining Box-Cox transformations, Fourier series with time-varying coefficients and autoregressive moving average error correction, ETS Error-Trend-Seasonal framework, MAD mean absolute deviation, MAPE mean absolute percentage error, RMSE root mean square error, MER mean error rate, RMSPE root mean square percentage error.

    Table S9.  Comparisons of the out-of-sample forecasting powers between TBATS models and SARIMA models, along with ETS models

    This study had some limitations. First, SF is a mild illness and has rarely led to death since the 20th century [3]. Therefore, infected individuals with mild clinical manifestations sometimes fail to seek medical aid, resulting in under-reporting and under-diagnosis. Second, to ensure that the model obtained a satisfactory forecasting result, it is important to note that this model should be updated with new incidence data. Third, to investigate whether our TBATS model was adequate for estimating the SF epidemics in other study regions or other infectious diseases, much work is still needed. Finally, integrating the factors influencing the SF epidemics may improve the predictive power. Nevertheless, we are not able to perform such an analysis due to the unavailability of a multivariate TBATS method and SF-related factors.

    In summary, SF had dual seasonal behaviors, peaking in May–June and November–December, with a recurrence in 2011 in China; since then it started to be increasing in the SF incidence. The TBATS method was advantageous in analyzing the long-term epidemiological seasonality and trends of SF, which can be considered a useful and flexible alternative to aid stakeholders to develop practical solutions to stop the ongoing spread of SF in China. In addition, we re-established the preferred TBATS (0.023, {0,0}, 0.895, {<12,5>}) specification based on the 16 years of data to predict the SF incidence into 2022, although the SF incidence was predicted to reach a plateau in the next 3 years [Supplementary Tables S1 and S10 (available in www.besjournal.com), and Figure 1], the SF incidence remained at a high level, suggesting that additional or comprehensive interventions must be developed to manage this evolving scenario.

    TimeForecastsLower 95% CIUpper 95% CITimeForecastsLower 95% CIUpper 95% CI
    Jan-200.4790.3300.695Jul-210.4260.1471.204
    Feb-200.2240.1430.349Aug-210.2200.0730.647
    Mar-200.3800.2290.627Sep-210.2500.0810.751
    Apr-200.6270.3601.085Oct-210.4120.1321.249
    May-200.8430.4631.521Nov-210.6160.1951.890
    Jun-200.9300.4891.750Dec-210.8160.2542.543
    Jul-200.4430.2210.879Jan-220.4340.1301.400
    Aug-200.2280.1080.474Feb-220.2040.0590.686
    Sep-200.2590.1190.557Mar-220.3510.1001.186
    Oct-200.4240.1900.935Apr-220.5840.1661.988
    Nov-200.6320.2751.429May-220.7910.2222.722
    Dec-200.8350.3541.938Jun-220.8780.2423.073
    Jan-210.4430.1801.069Jul-220.4210.1111.533
    Feb-210.2080.0810.524Aug-220.2180.0550.823
    Mar-210.3570.1370.914Sep-220.2480.0620.953
    Apr-210.5930.2231.543Oct-220.4090.1021.574
    May-210.8010.2962.126Nov-220.6120.1512.369
    Jun-210.8890.3202.412Dec-220.8110.1993.172
      Note. TBATS an advanced innovation state-space modelling framework by combining Box-Cox transformations, Fourier series with time-varying coefficients and autoregressive moving average (ARMA) error correction, CI, confidence interval.

    Table S10.  Projection into the future 36 months from January 2020 to December 2022 using the best TBATS (0.023, {0,0}, 0.895, {<12,5>}) model

    Data Availability All the data supporting the findings of the work are contained within the Supplementary Material (Supplementary Table S11, available in www.besjournal.com).

    Figure S8.  The extracted level, slope, and seasonal components of the best TBATS method developed on the basis of the data from January 2004 to December 2018.

    Figure S9.  Statistical checking for the residual series of the TBATS model developed on the basis of the data from January 2004 to December 2016. (A) Time series plot showing the residual series. (B) Autocorrelation function (ACF) plot. (C) Partial autocorrelation function (PACF) plot. The correlogram shows that although there are some residual ACF and PACF for the forecast errors that exceed the significance bounds, this is the best model used to estimate the upcoming 36-month epidemics of scarlet fever.

    Figure S11.  Statistical checking for the residual series of the SSARIMA model developed on the basis of the data from January 2004 to December 2018. (A) Time series plot showing the residual series. (B) Autocorrelation function (ACF) plot. (C) Partial autocorrelation function (PACF) plot. The correlogram shows that the ACF and PACF for the forecast errors do not exceed the significance bounds except for lags 11, 23, and 36, indicating that there is little evidence of non-zero autocorrelations at lags 1–36.

    Figure S12.  Statistical checking for the residual series of the SARIMA model developed on the basis of the data from January 2004 to December 2016. (A) Time series plot showing the residual series. (B) Autocorrelation function (ACF) plot. (C) Partial autocorrelation function (PACF) plot. The correlogram shows that the ACF and PACF for the forecast errors do not exceed the significance bounds except for lags 11 and 36, indicating that there is little evidence of non-zero autocorrelations at lags 1–36.

    Continued
    Parameters48-step ahead forecasting72-step ahead forecastingFuture 36-step ahead forecasting
    TBATS (0.037, {0,0}, 0.901, {<12,5>})TBATS (0.029, {0,0}, 0.87, {<12,4>})TBATS (0.023, {0,0}, 0.895, {<12,5>})
    Seed states
    1−2.3577−2.3794−2.3993
    20.10620.10640.1117
    3−0.0314−0.0546−0.0067
    4−0.0815−0.1004−0.0690
    5−0.0103−0.0196−0.0068
    60.02200.01600.0300
    70.02680.14110.0336
    80.1425−0.59400.1460
    9−0.5902−0.0762−0.6235
    10−0.0647−0.1476−0.0634
    11−0.1437−0.1522
    12−0.0212−0.0232
    Sigma0.18680.19470.1871
    AIC−197.7402−203.9513−160.8382
      Note. TBATS, an advanced innovation state-space modelling framework by combining Box-Cox transformations, Fourier series with time-varying coefficients and autoregressive moving average (ARMA) error correction, AR autoregressive method, AIC Akaike's Information Criterion.
    VariablesEstimateS.E.tPStationary R2R2AICBICLjung-Box Q (18)
    StatisticsP
    SARIMA (0,1,0)(0,1,1)12 method developed on the basis of the data from January 2004 to December 2015
    SMA10.7900.0948.394< 0.0010.2770.915−6.086−6.07212.5640.765
    SARIMA (0,1,0)(0,1,1) 12 method developed on the basis of the data from January 2004 to December 2013
    SMA10.8130.1356.025< 0.0010.2760.883−6.055−6.03212.5100.768
      Note. SARIMA seasonal autoregressive integrated moving average method, S.E. standard error, AIC Akaike's Information Criterion, BIC Schwarz's Bayesian Information Criterion, SMA1 seasonal moving average method at lag 1.

    Table S6.  The obtained best SARIMA methods with their key parameters on the basis of different datasets

    ModelsTesting horizons
    MADMAPERMSEMERRMSPE
    48-step ahead predictions
     SARIMA1.186252.3481.8502.6703.440
     TBATS0.09818.3070.1390.2200.222
    Reduced percentage (%)
     SARIMA vs. TBATS91.73792.74592.48691.76093.547
    72-step ahead predictions
     SARIMA0.800171.3811.5111.9012.810
     TBATS0.20746.0420.2480.4910.477
    Reduced percentage (%)
     SARIMA vs. TBATS74.12573.13583.58774.17183.025
      Note. SARIMA seasonal autoregressive integrated moving average method, TBATS an advanced innovation state-space modelling framework by combining Box-Cox transformations, Fourier series with time-varying coefficients and autoregressive moving average (ARMA) error correction, MAD mean absolute deviation, MAPE mean absolute percentage error, RMSE root mean square error, MER mean error rate, RMSPE root mean square percentage error.

    Table S7.  The comparisons of the predicted results between the SARIMA model and the TBATS model on different testing datasets

    DateIncidenceDateIncidenceDateIncidenceDateIncidenceDateIncidenceDateIncidence
    01-20040.029809-20060.10105-20090.262301-20120.329409-20140.177205-20170.8192
    02-20040.045210-20060.164706-20090.278202-20120.167510-20140.301106-20170.8006
    03-20040.094511-20060.280907-20090.150803-20120.246611-20140.514907-20170.3811
    04-20040.163912-20060.329108-20090.082104-20120.328112-20140.629808-20170.1711
    05-20040.180901-20070.23109-20090.094605-20120.532501-20150.49909-20170.2162
    06-20040.189602-20070.09510-20090.098806-20120.501802-20150.207910-20170.3078
    07-20040.11303-20070.135611-20090.079807-20120.252303-20150.275611-20170.5717
    08-20040.061904-20070.247412-20090.096108-20120.104704-20150.434212-20170.7641
    09-20040.074205-20070.343501-20100.06909-20120.132305-20150.664201-20180.5421
    10-20040.10306-20070.3902-20100.035310-20120.207106-20150.726902-20180.1547
    11-20040.203307-20070.199203-20100.063611-20120.326407-20150.376503-20180.2705
    12-20040.229908-20070.093204-20100.100812-20120.374108-20150.156804-20180.4862
    01-20050.168609-20070.112205-20100.169201-20130.243709-20150.207605-20180.7702
    02-20050.058210-20070.169506-20100.178802-20130.076210-20150.309606-20180.768
    03-20050.118211-20070.254307-20100.113803-20130.153611-20150.500707-20180.3859
    04-20050.218212-20070.315308-20100.071404-20130.225812-20150.645908-20180.1602
    05-20050.250401-20080.19309-20100.087205-20130.338101-20160.444909-20180.1967
    06-20050.270702-20080.072410-20100.114906-20130.310602-20160.1510-20180.3738
    07-20050.134403-20080.165511-20100.239907-20130.175603-20160.284311-20180.7114
    08-20050.080104-20080.228612-20100.305308-20130.085804-20160.358412-20180.9025
    09-20050.089105-20080.36701-20110.233309-20130.107705-20160.577101-20190.6313
    10-20050.129606-20080.281602-20110.074210-20130.15806-20160.592302-20190.1851
    11-20050.212907-20080.134303-20110.207411-20130.268907-20160.307803-20190.3657
    12-20050.236108-20080.068804-20110.372812-20130.392908-20160.139604-20190.4974
    01-20060.123909-20080.08705-20110.690801-20140.256809-20160.189805-20190.649
    02-20060.064510-20080.12906-20110.725302-20140.089910-20160.273406-20190.718
    03-20060.134511-20080.183107-20110.383903-20140.211311-20160.440907-20190.4188
    04-20060.211712-20080.236308-20110.161404-20140.317112-20160.550408-20190.1711
    05-20060.240601-20090.107909-20110.208805-20140.5201-20170.333309-20190.2481
    06-20060.270102-20090.062510-20110.36906-20140.552502-20170.16810-20190.3839
    07-20060.145403-20090.139311-20110.6307-20140.283203-20170.339411-20190.7296
    08-20060.076804-20090.229212-20110.719608-20140.133604-20170.489512-20190.9323

    Table 11.  The SF incidence time series data from January 2004 to December 2019

Reference (10)
Supplements:
21432Supplementary Materials.pdf

Catalog

    /

    DownLoad:  Full-Size Img  PowerPoint
    Return
    Return