Introduction

Acute hemorrhagic conjunctivitis (AHC) is a viral infectious disease that affects the conjunctiva, which is primarily caused by the Coxsackievirus A24 (CVA24) and the Enterovirus 70 (EV70).1,2 The disease is characterized by a short incubation period (12–48 h), rapid onset with severe eye irritation symptoms, marked conjunctival congestion, and commonly accompanied by subconjunctival hemorrhage and punctate epithelial keratitis.2,3 It commonly spreads through direct contact (eg, shaking hands or sharing items with an infected person).2,4 AHC is highly contagious and can affect individuals of all ages, but children and young adults are more susceptible.3 This disease is characterized by a high incidence rate, rapid transmission, and concentrated outbreaks.1,5,6 In past years, AHC has been widespread globally (eg, Asia, Africa, the Americas, and Europe),1 resulting in an estimated 10 million or more individuals being affected.7 Also, outbreaks have occurred in both developed and developing countries such as Japan,8 India,9 Thailand,10 Brazil,5 China,3 and most of Africa countries.6 AHC has so far caused serious harm to people’s lives, work, and social production, leading to school closures, production shutdowns, and market suspensions in some provinces, China (eg, Hubei, Guangxi, Hainan, Anhui).11 Henan Province, in particular, warrants focused attention due to its large population (approximately 100 million) and strategic location as a central transportation hub in China, where dense population movement may facilitate rapid AHC transmission. Moreover, Henan has historically experienced AHC outbreaks with substantial socioeconomic impact, yet systematic forecasting studies for this province remain scarce.12 Regrettably, there is currently no specific effective treatment. In this scenario, monitoring AHC incidence is the cornerstone of prevention and control. Forecasting can flexibly utilize the monitoring data, which is a cornerstone of risk management. Accurate predictions are pivotal for developing proactive public health policies, triggering early warning systems, guiding the timely deployment of interventions, and facilitating the efficient allocation of medical supplies and personnel during an outbreak.13–16 This enables public health systems to proactively respond, minimizing the impact on individuals and communities and preventing further transmission.

Various significant models are employed for forecasting and developing an early warning system using statistical and machine learning techniques to predict the trend and seasonality of infectious diseases, such as the seasonal autoregressive integrated moving average (SARIMA) model,14,15 exponential smoothing model,14 generalized regression neural network,17 support vector machine regression,18 and innovation state-space framework.19 Among them, the SARIMA model is the most widely used and flexible technique for time series forecasting (including financial data, economic indicators, climate variables, and health field) because of its ability to capture dynamics and account for random shocks, relatively easy interpretability, and relatively reliable predictions.20 Scholars have also extensively used the SARIMA model to forecast the incidence of AHC,14,15 hemorrhagic fever with renal syndrome (HFRS),19 tuberculosis (TB),21 and mumps.22 However, the SARIMA is primarily designed for capturing simple patterns such as trend, seasonality, and short-term flotations.23 This limitation can be costly in a public health context. Failure to capture complex, long-term dependencies may lead to inaccurate forecasts, resulting in either a late response (allowing an outbreak to grow) or an overly cautious response (wasting resources). Moreover, an integer difference in the SARIMA model for a series exhibiting long memory can generate excessive differencing and the discarding of valuable features in the series, leading to deviations in parameter estimation and modelling.20,24 Unlike the SARIMA model, the seasonal autoregressive fractionally integrated moving average (SARFIMA) model using a fractional differencing allows for the estimation of series with long-range dependence and complex seasonality.23,25,26 Meanwhile, the SARFIMA model does not involve complex machine learning algorithms,26 which allows the end-users to understand how the forecasting model is established, increasing their reliance on the model for decision-making purposes. AHC incidence is influenced by a confluence of virological, environmental, and demographic factors,2,11,27 and these factors collectively give rise to long‑memory properties in the incidence series. Specifically, the persistence of predominant viral strains across multiple seasons can sustain transmission even after temporary declines, creating autocorrelation that decays hyperbolically rather than exponentially. Similarly, multi‑year climate oscillations can modulate temperature and humidity over extended periods, thereby imposing slowly evolving seasonal baselines. At the population level, gradual demographic shifts, including urbanization and changes in school enrolment, alter contact patterns and susceptibility pools over years rather than months. These mechanisms imply that shocks to AHC transmission do not dissipate rapidly; instead, their effects persist and accumulate, a hallmark of long‑range dependence. Therefore, capturing such long‑memory dynamics is essential for accurate forecasting, and the SARFIMA model offers a principled approach to this challenge. Although the SARFIMA model has demonstrated good forecasting accuracy in the climatic, financial market, and commercial data,23,26,28–30 an investigation of how the SARFIMA model contributes to early warning and control for AHC incidence is still lacking. To address the gaps, there are two primary aims of this study: (1). To investigate whether the SARFIMA model is useful in monitoring and containing AHC epidemics in Henan province, China. (2). To ascertain the predictive accuracy of the SARFIMA with SARIMA models in estimating AHC incidence, with the goal of providing a more robust tool for public health policymakers to implement timely and cost-effective intervention strategies. To further assess the robustness and generalizability of the proposed models beyond the primary study region, we additionally validated the forecasting performance using national‑level AHC incidence data from mainland China.

Materials and Methods Data Extraction

The monthly incident cases of AHC between January 2011 and June 2023 were provided by The Health Commission of Henan Province (https://wsjkw.henan.gov.cn/zfxxgk/yqxx/) and The National Population and Health Science Data Sharing Platform (https://www.phsciencedata.cn/Share/en/index.jsp). The population data were collected from the Henan Statistical Yearbook. All AHC cases were confirmed based on the clinical diagnostic criteria issued by The National Health Commission of the People’s Republic of China (https://www.nhc.gov.cn/wjw/s9491/wsbz.shtml). According to these criteria, a clinical diagnosis of AHC is established when a patient presents with acute onset of severe eye irritation, marked conjunctival congestion, and subconjunctival hemorrhage, in the absence of other identifiable causes. In the context of this national notifiable disease surveillance system, laboratory confirmation of specific viral strains (such as Coxsackievirus A24 or Enterovirus 70) is not routinely performed for all cases; therefore, the case definition for this analysis is predominantly clinical. In China, AHC is a notifiable disease, and confirmed (ie, clinically diagnosed) cases must be registered within 24 hours. Any errors or duplicate records need to be corrected at the end of the same month. The SARFIMA and SARIMA were built using data from January 2011 to June 2022, and then the forecasting performances of both models were evaluated using data from July 2022 to June 2023. Also, the data from mainland China during the same period was provided to corroborate the robustness of both models. These national-level data were obtained from the National Population and Health Science Data Sharing Platform and included monthly reported AHC case counts aggregated from all 31 provinces, autonomous regions, and municipalities of mainland China. The case definition and reporting procedures were identical to those used for the Henan data. The national dataset underwent the same preprocessing steps as the Henan data, including monthly reconciliation of duplicate records and exclusion of suspected cases not meeting the clinical criteria. The training (January 2011–June 2022) and testing (July 2022–June 2023) periods were kept consistent with the primary analysis to ensure comparability.

The SARIMA Model Construction

The ARIMA model is a popular statistical model that combines three components to fit and estimate a time series, including the autoregressive (AR), the differencing (I), and the moving average (MA) components.20 The SARIMA model is an extension of the ARIMA model by incorporating the seasonal component of a time series, which has all the components of the ARIMA (AR, I, and MA) model and adds a seasonal component that consists of additional seasonal autoregressive (SAR), seasonal differencing (SI), and seasonal moving average (SMA) terms.14 The notation of the SARIMA model can be expressed as SARIMA(p, d, q)(P, D, Q)S, where p, d, and q are the orders of the non-seasonal AR, I, and MA components, respectively, P, D, and Q are the orders of the SAR, SI, and SMA components, respectively, and S is the length of the seasonal cycle (12 for monthly data).20 Determining the appropriate values of p, d, q, P, D, and Q for a SARIMA model includes four steps. First, although the SARIMA model can handle non-stationary series, it is crucial to assess stationarity during the modelling process. This assessment helps in determining model parameters, improving model performance, satisfying model assumptions, and simplifying the model structure. Therefore, the SARIMA model requires the AHC incidence data to be stationary in this study.20 If the data is shown to be non-stationary by a KPSS unit root test (P < 0.05),31 it needs to be differenced to achieve stationarity. Second, initial parameter ranges were determined by examining autocorrelation function (ACF) and partial ACF (PACF) plots of the differenced series.20 Trial and error were then used to determine the preferred orders by comparing Akaike’s information criterion (AIC), corrected AIC (CAIC), Bayesian information criterion (BIC), and the log-likelihood (LL). Specifically, the non-seasonal AR and MA orders (p, q) were each searched over the range 0–2, while the seasonal AR and MA orders (P, Q) were searched over 0–2, with the seasonal period fixed at S = 12 for monthly data. The non-seasonal differencing order d and seasonal differencing order D were determined from the KPSS test results (d = 1, D = 1) and kept fixed during the candidate search. All combinations satisfying the stationarity and invertibility conditions were considered, and the one with the least value of AIC, BIC, and CAIC, coupled with the greatest value of LL, was considered the optimal model.20 Third, the resulting residuals were investigated to determine whether they meet the requirements of white noise by analyzing the statistics of the Ljung-Box Q test, ACF, and PACF.20 Lastly, once a suitable SARIMA model was fitted to the training data, it can be used to make forecasts for future time steps.

The SARFIMA Model Construction

In time series analysis, the previous observations often have a persistent impact on subsequent observations, meaning that the rate of decay in the recent past can influence the future decay rate, which is called hyperbolic decay (HD) time series.26 By unraveling the underlying mechanisms and developing accurate models for HD series, researchers can gain valuable insights into complex systems and make informed predictions.26 The SARIMA model assumes an exponential decay in the ACF plot, but the SARFIMA model considers an HD, which characterizes the presence of long-term memory.24,26 For this reason, the informed SARFIMA model has attracted significant attention due to its ability to capture the persistence in such data by adding a fractional differencing.23,24 The fractional differencing parameters df (non-seasonal) and Df (seasonal) quantify the persistence of shocks over time, capturing long-memory behavior in the non-seasonal and seasonal components, respectively. According to standard ARFIMA theory, both stationarity and invertibility require −0.5 < df < 0.5 and −0.5 < Df < 0.5. Within this admissible range, the interpretation is as follows:25 (1) if df= 0 (or Df= 0), the series exhibits short memory and a mean-reverting process; (2) if −0.5 < df < 0 (or −0.5 < Df < 0), the series is anti-persistent, meaning that shocks tend to reverse more quickly than a random walk; (3) if 0 < df< 0.5 (or 0 < Df< 0.5), the series exhibits positive long-memory dependence, where past shocks decay hyperbolically rather than exponentially. Values of |df| or |Df| exceeding 0.5 generally violate the assumptions of stationarity or invertibility and are rarely employed in practical applications without further transformation. Values approaching 1 indicate a unit root process, requiring integer differencing. The notation of the SARFIMA model is often written as SARFIMA (p, d*q)(P, D*Q)S, where d* = d + df and D* = D + Df. Most commonly, df or Df ∈ (−1, 0.5) is the fractional component, and d or D ≥ 0 is always the integer component.26 The Hurst exponent (H) is a measure of long-term memory or persistence in a series. Often, the fractional differencing is represented as df or Df = H - 0.5. H ranges between 0 and 1, with H < 0.5 indicating anti-persistence or mean-reversion, with H = 0.5 indicating random or uncorrelated behavior, and H > 0.5 indicating persistence or trending behavior.24,32 The Hurst exponent is calculated using the concept of corrected rescaled range (R/S) analysis in our study.32 According to the description by Veenstra,26 the multi-modal nature of the likelihood function means it has multiple local optima; therefore, we initiate optimizations from different starting points to increase the probability of converging to the global optimum. Specifically, we employed 20 random starting points drawn from the stationary and invertible region (−0.5 < df < 0.5 and −0.5 < Df < 0.5), in addition to the default starting values provided by the “arfima” package. This approach ensures that the optimization process is less likely to be trapped in a local optimum. In this study, the estimation algorithm for the SARFIMA model produced multiple solutions with nearly identical likelihoods. The solution with the maximum LL, together with the minimum AIC and BIC, was selected as the final mode.26,33 The other processes of establishing the SARFIMA model (eg, parameter estimation and model diagnostics) followed what was described in the SARIMA model.

Statistical Analysis

The classical multiplicative decomposition was used to obtain the trend, seasonal, and irregular components.34 Among these components, seasonal relative (SR) indicated the amount by which the incidence for that specific period tends to be above (or below) the mean level (SR > 1, suggesting a high-risk season; otherwise, the opposite).35 Average annual percentage change (AAPC) and annual percentage change (APC) with a 95% confidence interval (CI) were computed to indicate the changing secular trend using Joinpoint (Version 4.8.0.1).36 The SARIMA and SARFIMA models were established based on the “forecast” and “arfima” packages in R (version 4.2.0, R Development CoreTeam, Vienna, Austria). Prior to model fitting, we assessed the linearity assumption of the time series using the Tsay test for nonlinearity.37 This test evaluates the null hypothesis that the underlying data-generating process follows a purely linear autoregressive structure, against the alternative of nonlinear dynamics. A non-significant result (P > 0.05) would support the use of linear frameworks such as SARIMA and SARFIMA. The autoregressive conditional heteroscedasticity (ARCH) effect means that different observed data in time series have various variances.38 Investigating this effect can help improve the efficiency of parameter estimation and the accuracy of interval forecast.38 Therefore, the Lagrangian Multiplier (LM) was used to examine whether there was an ARCH effect of the forecast errors from both models. Moreover, the COVID-19 pandemic introduced a distinct and abrupt structural break in health-seeking behavior, surveillance intensity, and population mobility, which could affect AHC reporting,39 and hence the binary COVID-19 indicator assigning to “1” from January 2020 to June 2023 and “0” from January 2011 to December 2019 was included as an external regressor in both the SARIMA and SARFIMA models to adjust for the effect of COVID-19 on predictive accuracy (while a single binary indicator cannot fully capture the time-varying intensity of the pandemic and associated non-pharmaceutical interventions, it provides a parsimonious adjustment for the overall interruption period). The accuracy of both models was evaluated by comparing the predicted values with the actual values of the test data using various performance metrics such as mean absolute deviation (MAD), root mean square error (RMSE), mean absolute percentage error (MAPE), and mean error rate (MER). A lower value for these metrics indicates a more accurate forecast, which directly translates to a lower risk of public health policy errors—such as over- or under-estimating the required resources for an impending outbreak.

(1)

(2)

(3)

(4)

where denotes the actual observed value at time I , denotes the corresponding predicted value , represents the mean of the actual observations over the validation period, and n is the number of predicted time points (12 in the holdout period). All four metrics are reported for both the training (fitting) and testing (forecast) horizons to allow for direct comparison of in-sample fit and out-of-sample predictive performance.

Results Data Description

A total of 28,503 AHC cases were reported in Henan Province during the study period (January 2011 – June 2023), corresponding to an annualized incidence rate of 2.322 per 100,000 population and a monthly incidence rate of 0.193 per 100,000 population. Joinpoint regression analysis revealed a significant overall upward trend in AHC incidence, with an AAPC of 6.973% (95% CI: 3.963% to 10.071%; t = 4.628, P < 0.001) (Figure 1). The temporal pattern was characterized by two distinct segments: a sharp significant increase from 2011 to 2017, with an APC of 14.841% (95% CI: 9.516% to 20.426%; t = 6.891, P < 0.001), followed by a non-significant slight decrease from 2017 to 2022, with an APC of −1.760% (95% CI: −6.547% to 3.272%; t = 0.840, P = 0.428) (Figures 1 and 2a andb). The peak annual case count was 3,019 (incidence rate: 3.08 per 100,000) in 2017, which was 1.977 times higher than the lowest count of 1,527 cases (incidence rate: 1.608 per 100,000) in 2012. Classical multiplicative decomposition confirmed a strong and consistent seasonal pattern. The SR was greater than 1 from March to September each year, indicating high-risk seasons, and less than 1 from January to February and October to December, indicating low-risk seasons (Figures 2c and d and S1).

Figure 1 Annual reported AHC cases and incidence rate (per 100,000 population) in Henan Province, 2011–2022. The incidence exhibited an increasing trend from 2011 to 2017, peaking in 2017, followed by a relatively stable period from 2018 to 2022.

Figure 2 Classical multiplicative decomposition of the AHC incidence series in Henan Province, January 2011–June 2023. Panels show: (a) original observed series, (b) trend component, (c) seasonal component, and (d) irregular (remainder) component. The series exhibits a clear upward secular trend and pronounced annual seasonality.

Linearity Assessment

For the AHC incidence series, the Tsay test for nonlinearity produced a statistic of 0.665 (P = 0.793) in Henan and 0.912 (P = 0.640) in mainland China. The observed high P-values indicate that there is no significant evidence of nonlinearity, thus confirming that the temporal dependence structure of the AHC series can be adequately approximated by SARIMA and SARFIMA models.

The Best-Fitting SARIMA Model

The value of KPSS-statistic was 2.977 greater than the critical value for a significance level of 0.574, and thus P < 0.05 suggested the original series displays non-stationarity. The series was then seasonally differenced (KPSS = 0.626 > 0.574) and non-seasonally differenced (KPSS = 0.015 < 0.574) once, which help achieve stationarity. So, we determined d = 1 and D = 1. Subsequently, by depicting the ACF and PACF for the stationary series (Figure S2), the trial and error identified eight candidates with significant parameters (Table 1). It is apparent that the SARIMA(0,1,2)(0,1,1)12 model generated the minimum AIC (1231.50), CAIC (1232), and BIC (1245.64), as well as the maximum LL (−610.75), thus this model was our preference. For this best-fitting model, since the correlogram shows that all the spikes fell within the significance limits (Figures 3a and b) and the P-values for Ljung-Box Q test were greater than 0.05 at different lags (Table 2), showing no remaining correlations in the residuals. Besides, the P-values for LM test were greater than 0.05 at most lags, indicating that the ARCH effect was removed to a great extent. This best-fitting model passes the required checks and makes forecasts for future 12 data (Table 3). Similarly, in the validation analysis, we determined the best-fitting SARIMA(1,0,0)(2,1,1)12 model for the national AHC incidence series based on the modelling steps, and the model checks and forecasts are provided in Tables S1–S2 and Figures 3a and b.

Table 1 The Selected SARIMA Candidates with Their AIC, CAIC, BIC, and LL

Table 2 Ljung-Box Q and LM Tests for the Resulting Forecast Errors Under the SARIMA and SARFIMA Models

Table 3 Forecasts Between July 2022 and June 2023 Under the SARIMA and SARFIMA Models

Figure 3 Autocorrelation function (ACF) and partial autocorrelation function (PACF) plots of residuals from the optimal SARIMA and SARFIMA models for Henan data. (a) ACF and (b) PACF of SARIMA residuals; (c) ACF and (d) PACF of SARFIMA residuals. No spikes exceed the 95% confidence bounds except for a single significant lag at lag 5 in panels (c) and (d), indicating that the residuals are largely white noise with minor residual periodic structure.

The Best-Fitting SARFIMA Model

By calculating the corrected R/S for the training data, the corrected R/S was 0.757 pinpointing strong long memory in the series, and hence the SARIFMA can be used to model the AHC incidence series. Then, according to the modelling steps, we determined the SARFIMA(0,0.410,2)(0,0.413,1)12 structure as our optimal model, which resulted in seven possible modes (Table S3). By synthetically considering the AIC, BIC, and LL, we found that the SARFIMA(0,0.410,2)(0,0.413,1)12 model with mode 1 because this mode gave a smaller AIC (987.594) and BIC (1011.012), together with a greater LL (485.797). The resulting residuals are plotted in Figures 3c and d. No correlations exceeded the significance bounds except for the one at lag 5 and the P-values for Ljung-Box Q test were greater than 0.05 at various lags (Table 2). So, the errors behave like white noise. Also, the P-values for LM test were greater than 0.05 except for the one at lag 5, suggesting the removal of the ARCH effect to a great extent. These results imply that the optimal model passes the required checks (although a single significant spike at lag 5 was observed in the SARFIMA models residuals, the overall Ljung-Box Q test was non-significant, and the model was deemed adequate. This may indicate a minor periodic effect not fully captured by the model) and forecasts for future 12 data are generated based on the estimated parameters (Table 3). Likewise, in the validation analysis (corrected R/S = 0.695), we identified the preferred SARFIMA(1,0,0)(2,0.409,1)12 model with mode 1 for the national AHC incidence series following the modelling processes, and the model checks and forecasts are listed in Tables S1–S4 and Figures 3c and d.

Forecasting Accuracy and Model Comparison

The forecasting performance of both models was rigorously evaluated on the test set (July 2022 to June 2023). The SARFIMA model consistently demonstrated superior accuracy across all error metrics compared to the SARIMA model. Specifically, for the Henan data, SARFIMA yielded lower errors: MAD = 39.375 vs 48.850; MAPE = 0.237 vs 0.280; RMSE = 50.390 vs 57.358; and MER = 0.184 vs 0.229 (Table 4). Besides, a visual comparison of the fitted and forecasted values against the actual observations further confirmed that the SARFIMA model more closely tracked the actual trend and seasonal fluctuations of AHC incidence (Figure 4). To corroborate the robustness of our findings, the same modelling procedure was applied to the national AHC incidence series from mainland China. The optimal models identified were the SARIMA(1,0,0)(2,1,1)12 and SARFIMA(1,0,0)(2,0.409,1)12 structures (Corrected R/S = 0.695). Diagnostic checks confirmed both models were adequate (Figure S3), and for the national validation data, SARFIMA also outperformed SARIMA on MAD (494.680 vs 524.934), MAPE (0.227 vs 0.251), and MER (0.217 vs 0.231), while the RMSE values were nearly identical between the two models (604.94 vs 604.765) (Table 4 and Figure 4), suggesting that both models exhibited comparable performance in terms of this metric at the national scale.

Table 4 Comparison of the Performance Indices in the Fitting and Forecast Horizons Under the SARIMA and SARFIMA Models

Figure 4 Comparison of observed versus fitted and forecasted AHC incidence from the SARIMA and SARFIMA models. (a) Henan Province; (b) mainland China. The grey shaded areas indicate the 12‑month holdout forecast period (July 2022–June 2023). The SARFIMA model more closely tracks the observed trend and seasonal fluctuations compared with the SARIMA model in both datasets.

Besides, the ExponenTial Smoothing (ETS) approach has recently been shown to be effective for simulating infectious‑disease incidence time series (the details can be seen in the Supplementary ETS model and Table S5);40 therefore, as an additional benchmark we fitted ETS models to the AHC series and compared their forecasts with those from SARIMA and SARFIMA. The best ETS fits were ETS(M,MD,M) for Henan and ETS(M,AD,M) for mainland China (the modeling details can see Tables S6–S9). As reported in Table S10, SARFIMA outperformed both SARIMA and ETS over the fitting period, yielding lower MAD, MAPE, RMSE and MER at both scales. For 12‑month forecasts in Henan, SARFIMA achieved MAD = 39.375 (SARIMA 48.850; ETS 47.764) and MAPE = 0.237 (SARIMA 0.280; ETS 0.271). In mainland China the advantage was larger: SARFIMA MAD = 494.680 (SARIMA 524.934; ETS 643.439) and MAPE = 0.227 (SARIMA 0.251; ETS 0.288). These results indicate that explicitly modelling long‑range dependence via SARFIMA improves fit and forecast accuracy relative to both conventional SARIMA and the flexible ETS framework.

Discussion

This study demonstrates that the SARFIMA model, by effectively capturing the long-range dependence inherent in AHC incidence data, provides more accurate forecasts than the widely used SARIMA and ETS models in most metrics for the Henan dataset, a pattern partially replicated in the national validation—though the RMSE values were comparable between the two models at the national scale. From a public health policy and risk management perspective, this is a crucial finding. The improved accuracy offers a tangible pathway to reduce uncertainty in decision-making, allowing health authorities to move from a reactive to a proactive stance in managing AHC incidence. However, it is important to emphasize that these findings reflect retrospective predictive performance; the extent to which improved forecast accuracy translates into better real‑world intervention outcomes, such as earlier outbreak detection or more efficient resource allocation, remains to be demonstrated in prospective operational settings.

The SARIMA model has been the most common tool for analyzing and forecasting time series data in various fields such as economics, finance, meteorology, and healthcare, among others.20 Even for a non-stationary series, the SARIMA model can transform non-stationary data into stationarity by differencing the series.20 By incorporating trend, seasonal, and random components into the model,20 the SARIMA model can capture and account for the complex dynamics of the data, leading to relatively accurate and reliable forecasts. The SARIMA model has demonstrated successful applications in communicable diseases such as AHC,14,15 HFRS,19 TB,21 COVID-19,16 and mumps.22 This was also confirmed by our study although the predictive quality of the SARIMA model underperforms the SARFIMA model. However, the SARIMA’s limitations in capturing complex patterns like long-memory present a significant risk for public health. Over-differencing can lead to a loss of crucial information, potentially causing policymakers to underestimate the size and duration of an outbreak, leading to delayed interventions and inefficient resource use.20,23,25 The SARFIMA model extends the ARIMA model by the inclusion of fractional differencing,23 enabling it to show advantages over the SARIMA model in forecasting AHC incidence and in guiding potential containment efforts:23–26,28–30 (1) The SARFIMA model can effectively capture both short-term and long-term dependencies in AHC incidence series, providing a more accurate representation of real-world AHC incidence. (2) The SARFIMA model provides more accurate parameter estimation by taking into account the fractional differencing operator. (3) The SARFIMA model has inherent flexibility, allowing for handling irregularly spaced time series data with ease. (4) By incorporating seasonal and nonlinear fractional differences into the SARFIMA model, allowing for complex and multiple seasonal patterns in AHC incidence. (5) The SARFIMA model has the potential to produce more accurate forecasts over longer horizons. As such, it seems that the SARFIMA model is more worthy of being promoted, and encouraging its adoption can make informed decisions by enabling more accurate and reliable modelling and forecasting of AHC incidence, even all types of communicable diseases in comparison to the SARIMA model. Noting that this generalization requires further verification. Besides, recent studies have demonstrated satisfactory applications in estimating epidemics of diseases using new models such as Bayesian structural time series16 and innovation state-space framework.19 Therefore, future work into comparing and validating the predictive performance with the models above is recommended.

Previous studies indicated that there was an overall increase in AHC incidence in most provinces of China (eg, Hainan, Hubei, Anhui, Yunnan, among others),3 Brazil,5 Japan,8 and Egypt.6 Also, our study indicated together an upturn in AHC incidence during the study period (AAPC = 6.973%). However, there were two different incidence stages, including a dramatic increase in the period 2011–2017 with the most cases in 2017 (APC = 14.841%) and a relative plateau with a slight reduction in the period 2017–2022 (APC = −1.760%). This fits well with an earlier study indicating AHC showed a significant increase in most provinces of China except for Beijing, Tianjin, Shanghai, and Xinjiang during 2010–2018,3 also in agreement with Brazil where over 200,000 cases were confirmed during the summer of 2017 and 20185 and Thailand where over 100,000 persons were infected in 2014.10 There are several factors that may contribute to the increasing trend in AHC. First, the emergence of CVA24 variants has caused considerable concern in recent years. A recent study found that global CVA24v variants were divided into 8 genotypes,1 and the increased transmission of these viruses in populations contributes to a higher incidence. Second, urbanization and population growth may lead to overcrowding and poor sanitation (urban resident population up from 3,829 × 105 persons in 2011 to 5,579 × 105 persons in 2021 in Henan),41 facilitating the transmission of AHC. Third, the implementation of the national two-child policy in 2011 has led to a larger susceptible population (number of births up from 121 × 105 in 2011 to the most numbers of 140 × 105 in 2017 in Henan),41 thereby contributing to the increase in AHC endemics. Fourth, the globalization of travel and migration has led to easier and faster movement of individuals across borders. Fifth, timely reporting has significantly improved due to advancements in the medical field.11 Lastly, the absence of a widely accessible and effective vaccine leaves individuals susceptible to the infection,14 further fueling the increase in cases. Besides, a slight decline in AHC incidence in recent years may be associated with increased public awareness campaigns, increased budgets, and enhanced prevention and control.14,15

Seasonal patterns of AHC have been reported in many studies.3,14,15 A clear seasonal profile was noted (peak in spring-early autumn and a trough in late autumn-winter). This fits well with that in mainland China,3 Japan,8 India,9 Thailand,10 and Brazil,5 yet inconsistent with that in Egypt,6 which related to the differences in school breaks, population density, socioeconomic conditions, people’s lifestyles, climatic factors, and the predominant genotypes.1,5,6,11,27 Several factors may be associated with the high-risk seasons of AHC. First, summer and autumn are characterized by higher temperatures and humidity levels, which create a favorable environment for the survival and transmission of AHC-related causative agents.11,27 Second, people tend to spend more time outdoors during the spring and early autumn, engaging in recreational activities and traveling. This increased social interaction and proximity to others enhance the chances of transmitting the virus from person to person.11 Lastly, spring and early autumn coincide with the start of the school year and busy agricultural activities in China.11,27 Furthermore, the low AHC incidence in late autumn-winter could be related to the winter holidays and Spring Festival in China.

Limitations should also be acknowledged: First, a key limitation is the potential for under-reporting, as AHC is a self‑limiting disease and many mild cases may not seek medical care or be reported. Under‑reporting is likely to be non‑random: it may vary by season, by year, and by region. Such differential under‑reporting could distort both the secular trend and the seasonal amplitude. For example, if summer cases are more under‑reported than winter cases, the observed seasonal peak would be attenuated, leading to an underestimation of the high‑risk period duration and an overly optimistic forecast during peak months. Therefore, our forecasts should be interpreted as predictions of reported incidence rather than true community incidence. Users should consider adjusting thresholds or combining the model with syndromic surveillance data to account for potential under‑ascertainment when making operational decisions. Second, the analysis relied on clinically diagnosed cases from the national surveillance system. While this reflects the reality of most public health surveillance data—where laboratory testing is not feasible for all cases—it means that the time series captures the occurrence of the clinical syndrome of AHC, not necessarily infections with a specific viral strain (eg, CVA24 or EV70). This distinction should be considered when interpreting the model’s forecasts for specific viral variants. Third, the predictive ability often decreases as the predicted periods increase, and thus it requires regular updates to the model with new data. Fourth, while our training-testing split provides a standard evaluation framework, it represents a single holdout validation rather than a prospective or rolling-origin forecasting exercise. However, the relatively limited number of independent seasonal cycles in the current dataset constrained the implementation of such a scheme without substantial loss of training information. Fifth, parameter estimation for the SARFIMA model becomes more challenging with limited data. The requirement of long series with 100 or more samples should be satisfied. Sixth, although our external validation using national‑level data confirmed the superior performance of the SARFIMA model for AHC in mainland China, its generalizability to other regions and to other demographic groups remains untested. Therefore, we recommend that local health authorities conduct their own model identification and parameter estimation based on regional data before adopting it for operational forecasting. Seventh, the SARFIMA model may further enhance predictive performance by incorporating meteorological variables, school calendars, population mobility, and public health interventions. However, due to the unavailability of high-resolution, spatially matched data on these covariates over the entire study period, we were unable to include them in the current models. Eighth, our model does not incorporate viral genomic data. Integrating genomic information with incidence series could enhance biological interpretability by revealing how strain evolution influences long-memory transmission patterns. However, routine whole‑genome sequencing is not performed for all AHC cases under the current surveillance system, and available genomic data lack the spatiotemporal resolution required for seamless integration. Lastly, the potential influence of co‑infections with other ocular or respiratory pathogens on clinical diagnosis was not assessed. Co‑circulation of these pathogens may cause overlapping symptoms, leading to potential misclassification of reported cases.

The superior performance of the SARFIMA model carries substantial implications for health policy and risk management. First, the model’s ability to forecast AHC incidence with greater accuracy directly supports evidence-based decision-making. For instance, a one-month-ahead forecast during the peak season (March-September) can allow local health departments to preposition essential supplies (eg, antiviral eye drops, personal protective equipment) and deploy public health messaging campaigns to high-risk populations. Second, the long-range dependence captured by the SARFIMA model is critical for strategic planning. By understanding the persistent nature of AHC transmission, policymakers can make more informed decisions about long-term investments in public health surveillance, laboratory capacity, and community health education. This is particularly important for a disease like AHC, which has no specific treatment, where prevention and rapid response are the only defense. Third, the adoption of the SARFIMA model can refine the definition of “high-risk” seasons and trigger early warning systems. Our analysis identified March to September as high-risk months. The SARFIMA model can be integrated into a surveillance dashboard to provide a probabilistic forecast, automatically flagging when the predicted incidence exceeds a predefined threshold. This would allow for the automated, risk-based activation of response protocols, improving the efficiency and timeliness of public health interventions. However, the threshold for triggering alerts would need to be calibrated using region‑specific historical data, as the national validation revealed metric‑specific variability in performance. We recommend that pilot implementations start with retrospective simulation before moving to prospective use. While the SARFIMA model requires more data and computational expertise than simpler models, the investment is justified by the potential for improved health outcomes and more efficient resource use. Promoting the adoption of such sophisticated but interpretable models is a key step in modernizing public health risk management.

Conclusions

This study reveals an overall rising trend in AHC incidence in Henan, with a pronounced seasonal peak from March to September, underscoring the need for seasonally targeted control strategies. The SARFIMA model demonstrated improved forecasting accuracy over the SARIMA and ETS models in most metrics for the Henan dataset, and partially replicated this pattern in the national validation—though the RMSE values were comparable between the two models at the national scale. These findings suggest that explicitly modelling long-range dependence through fractional differencing can enhance the understanding and prediction of AHC incidence dynamics, providing a methodological advancement over conventional seasonal models. However, several important qualifications must be emphasized. First, the improved forecast accuracy observed in this retrospective study does not directly imply improved public health outcomes or more effective epidemic control; the actual benefits of model-guided interventions would require prospective implementation studies to evaluate whether more accurate forecasts translate into timely action and efficient resource use. Second, the generalizability of the SARFIMA model to other regions or to infectious diseases with different transmission mechanisms remains untested; the substantial heterogeneity in climate, population density, and healthcare systems across China suggests that region‑specific model calibration is essential before operational deployment. Third, the potential impact of under‑reporting—particularly differential under‑reporting across seasons—should be considered when interpreting the forecasts, as our models predict reported rather than true community incidence. Fourth, the predictive performance of both models is expected to degrade over extended forecast horizons, necessitating regular model updating as new surveillance data accumulate. We recommend that public health authorities treat these findings as proof‑of‑concept evidence for the potential utility of long‑memory time‑series models in AHC surveillance. Future research should focus on prospective validation using rolling‑origin evaluation frameworks, integration of covariates such as meteorological and mobility data, and multi‑site validation across diverse epidemiological settings. Until such evidence accumulates, local calibration and cautious interpretation remain essential before incorporating these models into routine surveillance and early‑warning workflows.

Data Sharing Statement

All the data supporting the findings of this work are derived from publicly accessible official surveillance sources. The Henan provincial AHC incidence data are available from The Health Commission of Henan Province (https://wsjkw.henan.gov.cn/zfxxgk/yqxx/), and the national data are accessible via the National Population and Health Science Data Sharing Platform (https://www.phsciencedata.cn/Share/en/index.jsp). The R code used for model fitting, diagnostic testing, and forecasting is available from the corresponding author (Prof. Yongbin Wang) upon reasonable request.

Ethics Approval and Consent to Participate

The institutional review board of Henan Medical University approved this study protocol (No: XYLL-2019072). All procedures were conducted in accordance with the Declaration of Helsinki (and its amendments). The need for informed consent was waived by the study Ethics Committee of Henan Medical University because the AHC cases were shared anonymously and we cannot access any identifying information of the patients (available from: https://www.phsciencedata.cn/Share/).

Acknowledgments

We thanked the Chinese CDC for sharing the AHC incidence series data in mainland China.

Author Contributions

All authors made a significant contribution to the work reported, whether that is in the conception, study design, execution, acquisition of data, analysis and interpretation, or in all these areas; took part in drafting, revising or critically reviewing the article; gave final approval of the version to be published; have agreed on the journal to which the article has been submitted; and agree to be accountable for all aspects of the work.

Funding

This work was supported by the Henan Provincial Natural Science Foundation (No.: 262300421599).

Disclosure

The authors report no conflicts of interest in this work.

References

1. Chen P, Lin X-J, Ji F, et al. Evolutionary phylogeography reveals novel genotypes of coxsackievirus A24 variant and updates the spatiotemporal dynamics in the population with acute hemorrhagic conjunctivitis. Int J Infect Dis. 2022;124:227–15. doi:10.1016/j.ijid.2022.10.007

2. Liu J, Zhang H. Epidemiological investigation and risk factor analysis of acute hemorrhagic conjunctivitis in huangshi port district, Huangshi city. Comput Math Methods Med. 2022;2022:3009589.

3. Liu R, Chen Y, Liu H, Huang X, Zhou F. Epidemiological trends and sociodemographic factors associated with acute hemorrhagic conjunctivitis in mainland China from 2004 to 2018. Virol J. 2022;19(1):34.

4. Fonseca MC, Pupo-Meriño M, García-González LA, et al. Molecular characterization of coxsackievirus A24v from feces and conjunctiva reveals epidemiological links. Microorganisms. 2021;9(3):531.

5. Sousa IP Jr, Burlandy FM, Ferreira JL, et al. Re-emergence of a coxsackievirus A24 variant causing acute hemorrhagic conjunctivitis in Brazil from 2017 to 2018. Arch Virol Apr. 2019;164(4):1181–1185.

6. Ayoub EA, Shafik CF, Gaynor AM, et al. A molecular investigative approach to an outbreak of acute hemorrhagic conjunctivitis in Egypt, October 2010. Virol J. 2013;10:96.

7. Baggen J, Hurdiss DL, Zocher G, et al. Role of enhanced receptor engagement in the evolution of a pandemic acute hemorrhagic conjunctivitis virus. Proc Natl Acad Sci U S A. 2018;115(2):397–402.

8. Harada K, Fujimoto T, Asato Y, Uchio E. Virol