Spatiotemporal evolution and forecasting of public attention to special education: a GIS and SARIMA-LSTM based analysis

Spatial correlation analysis

Moran’s I is a common indicator in spatial statistical analysis (Balash et al., 2020). It is typically categorized into Global and Local (Mathur, 2015). In this study, Moran’s I values are computed and visualized using GeoDa software to conduct a comprehensive analysis of the spatial correlation of public attention at both the overall and local levels. This enables the identification of spatial clustering patterns and regional disparities in public interest. The equations are:

(1) $$I={\displaystyle \frac{n}{{\sum }_{i=1}^n{\sum }_{j=1}^n{W}_{ij}}\times {\displaystyle \frac{{\sum }_{i=1}^n{\sum }_{j=1}^n{W}_{ij}\left(X_i-\overline{X}\right)\left(X_j-\overline{X}\right)}{{\sum }_{i=1}^n{\left(X_i-\overline{X}\right)}^2}}}$$ (2) $$I_L={\displaystyle \frac{n}{{\sum }_{i=1}^n{\sum }_{j=1}^n{W}_{ij}}\times {\displaystyle \frac{\left(X_i-\overline{X}\right){\sum }_{j=1}^n{W}_{ij}\left(X_j-\overline{X}\right)}{{\sum }_{j=1}^n{\left(X_i-\overline{X}\right)}^2}}}$$where I and I_L denote the Global and Local indices, respectively, X represents the observed indicator of interest, and X_i and X_j refer to the values of X in regions i and j. n is the number of regions under investigation, and W_ij is the spatial weight. In this study, W_ij is set to 1 if regions i and j share a geographical boundary, and zero otherwise.

To ensure the validity of the index, a Z-score is typically used to conduct a significance test. The equation is:

(3) $$Z\left(I\right)={\displaystyle \frac{I-E\left(I\right)}{\sqrt{VAR\left(I\right)}}.}$$

Here, Z(I) represents the standardized test statistic for significance, E(I) denotes the expected value of the Moran’s I index, and VAR(I) is its variance. According to the 95% confidence interval criterion, if |Z| exceeds 1.96, the spatial autocorrelation is considered statistically significant; conversely, if |Z| is less than 1.96, the observed values are regarded as randomly distributed (Ohana-Levi et al., 2019).

Local spatial autocorrelation analysis can be used to identify spatial clustering patterns within regions (Appice & Malerba, 2014), which are typically categorized into: High–High (H–H), Low–Low (L–L), High–Low (H–L), and Low–High (L–H). Among these, H–H and L–L indicate spatial association in public attention levels, while H–L and L–H reflect spatial heterogeneity. In this study, visual representations of local clustering patterns are generated using GeoDa software to illustrate the spatial distribution and evolution of public attention towards special education.

Forecasting model

(1) SARIMA model

The Autoregressive Integrated Moving Average (ARIMA) model is applied for forecasting non-stationary data in real-world applications (Box et al., 2015). It has been extensively applied across various domains in both the social and natural sciences. The idea of ARIMA is to eliminate local levels or trends in a non-stationary series through differencing, thereby transforming it into a stationary series that can be modeled using an ARMA structure (Fattah et al., 2018). The equation is:

(4) $$y_t=\mu +\sum _{i=1}^p{\gamma }_i{y}_{t-1}+\sum _{i=1}^q{\theta }_i{y}_{t-1}+e_t.$$

Here, y_t denotes the value of the differenced time series X at time t, μ is the constant term, while p and q indicate the respective orders of the autoregressive (AR) and moving average (MA) components. The term e_t indicates the error term, γ_i refers to the autoregressive coefficients, and θ_i represents the moving average coefficients.

On this basis, the SARIMA model integrates seasonal effects, which is suitable for analyzing and forecasting time series with cyclical temporal behavior (Kumar Dubey et al., 2021). It captures linear and seasonal structures in the data, which are essential for understanding recurring patterns in public attention.

The SARIMA model is usually expressed as SARIMA(p, d, q) × (P, D, Q)_s, where p, d, and q correspond to the orders of the non-seasonal autoregressive term, differencing, and moving average term, respectively. Meanwhile, P, D, and Q refer to the seasonal counterparts of these parameters, and s signifies the seasonal cycle length. The equation is:

(5) $${\mathrm{\Phi }}_P\left(B^s\right)\varphi \left(B\right){\left(1-B^s\right)}^D{\left(1-B\right)}^d{y}_t=c+{\mathrm{\Theta }}_Q\left(B^s\right){\theta }_q\left(B\right){\epsilon }_t.$$

In this equation, y_t represents the observed value of the time series at time t, and ε_t corresponds to the white noise at that point in time. The constant term is denoted by c. The operator B is the backshift (or lag) operator, where applying B^s to y_t yields y_t−s, indicating a shift of s time steps. The terms ϕ_p(B) and θ_q(B) refer to the polynomials for the non-seasonal AR and MA components of orders p and q, respectively. Similarly, Φ_P(B^s) and Θ_Q(B^s) denote the seasonal AR and MA polynomials of orders P and Q, constructed based on the seasonal lag s. The SARIMA framework is defined by seven key parameters: p, d, q, P, D, Q, and s. Selecting suitable values for these parameters is known as the model identification or order selection process.

The construction of a SARIMA model generally follows a sequence of critical steps, including data preprocessing, stationarity assessment, white noise verification, parameter identification, model diagnostics, and ultimately, forecasting. In this study, the Augmented Dickey–Fuller (ADF) test is applied to evaluate the stationarity of the series, while the Ljung–Box Q test is utilized to determine if the model residuals exhibit the characteristics of a white noise process (Paparoditis & Politis, 2018; Hassani & Yeganegi, 2019).

(2) LSTM model

The LSTM network is a type of recurrent neural network (RNN) specifically designed to capture long-term dependencies in sequential data. It incorporates three gating mechanisms—input, forget, and output gates—that regulate information flow and mitigate the vanishing and unstable gradient problems common in conventional RNNs (Gers, Schmidhuber & Cummins, 2000), making it effective for modeling complex, non-linear fluctuations in public attention data.

Its structure processes information from three inputs: the current input (x_t), the previous hidden state (s_t−1), and the previous cell state (C_t−1) (see Fig. 1). Its core components include the input gate (i_t), forget gate (f_t), output gate (o_t), and memory cell (C_t). These gates regulate information flow within the network, enabling it to retain essential data and discard irrelevant details, thereby supporting the modeling of long-term dependencies (Sherstinsky, 2020).

The forget gate controls the proportion of past cell state information that is maintained or eliminated. The equation is defined as:

(6) $$f_t=\sigma \left(W_f\cdot \left[S_{t-1},x_t\right]+b_f\right).$$

In this context, σ(·) represents the sigmoid activation function, W_f denotes the weight matrix corresponding to the forget gate, and b_f is its associated bias vector.

The input gate determines which new information is added to the cell’s memory. The equation is defined as:

(7) $$i_t=\sigma \left(W_i\cdot \left[S_{t-1},x_t\right]+b_i\right)$$

(8) $$\widetilde{C_t}=tanh\left(W_c\cdot \left[S_{t-1},x_t\right]+b_c\right).$$

In this equation, W_i and W_c denote the weight matrices associated with the input gate and the candidate memory vector, respectively, while b_i and b_c are their corresponding bias terms. $\widetilde{C_t}$ denotes the candidate memory cell state.

The updated cell state C_t is obtained by combining the outputs from the forget gate (f_t), input gate (i_t), and the candidate memory $\widetilde{C_t}$. This process performs a gated update, blending previous memory with newly processed information to support long-term dependency learning. The equation is as:

(9) $$C_t=f_t\ast C_{t-1}+i_t\ast \widetilde{C_t}.$$

The output gate controls what information from the memory cell is sent forward to influence subsequent computations. The equation is given as:

(10) $$o_t=\sigma \left(W_o\cdot \left[S_{t-1},x_t\right]+b_0\right).$$

Here, W_o refers to the weight matrix associated with the output gate, while b_o denotes its corresponding bias term.

To produce the hidden state s_t at the current time step, the model integrates o_t with C_t, modulated through a non-linear transformation. This operation determines the information passed on to subsequent time steps or network layers, as described by the following equation:

(11) $$s_t=o_t\ast tanh\left(C_t\right).$$

(3) SARIMA-LSTM model

To improve forecasting accuracy, this study constructs a model that combines the SARIMA with LSTM. This integrated approach is designed to address the limitations of using a single forecasting model (Peirano, Kristjanpoller & Minutolo, 2021). The SARIMA-LSTM model has been shown to produce more robust and accurate results across various application domains (Bilgili, Pinar & Durhasan, 2025).

Figure 2 illustrates the key steps of the workflow, including data preparation, data testing, SARIMA modeling and forecasting, residual extraction, LSTM modeling of the residual series, and the final combination of outputs to produce the complete forecast.

(4) Model evaluation

To provide a clear understanding of each method, a brief comparison is offered: ARIMA is simple and interpretable but performs poorly with seasonality and non-linear patterns; SARIMA captures seasonality and improves accuracy over ARIMA yet remains less effective for abrupt or complex non-linear changes; and SARIMA-LSTM handles both seasonal and non-linear dynamics, commonly achieving the highest accuracy though with greater computational cost.

The reliability of forecasting results is largely evaluated by comparing forecast values with actual observations. To evaluate the effectiveness of the model, a set of commonly used metrics is applied, including mean absolute percentage error (MAPE), root mean squared error (RMSE), directional accuracy (DA), and the coefficient of determination (R²) (Vivas, Allende-Cid & Salas, 2020; Chicco, Warrens & Jurman, 2021). These indicators provide quantitative insight into the precision, consistency, and trend-following ability of the forecasting model. The equations are:

(12) $$MAPE={\displaystyle \frac{1}{n}\sum _{i=1}^n\left|{\displaystyle \frac{y_i-y_i^{\mathrm{\prime }}}{y_i}}\right|}$$

(13) $$RMSE=\sqrt{{\displaystyle \frac{1}{n}\sum _{i=1}^n{\left(y_i-y_i^{\mathrm{\prime }}\right)}^2}}$$

(14) $$DA={\displaystyle \frac{1}{n}\sum _{i=1}^n{a}_i}$$

where,

$$a^i=\{\begin{array}{c}1,\left(y_i-y_{i-1}\right)\left(y_i^{\mathrm{\prime }}-y_i^{\mathrm{\prime }}\right)\ge 0\\ 0,\left(y_i-y_{i-1}\right)\left(y_i^{\mathrm{\prime }}-y_i^{\mathrm{\prime }}\right)<0\end{array}$$

(15) $$R^2=1-{\displaystyle \frac{{\sum }_{i=1}^n{\left(y_i-y_i^{\mathrm{\prime }}\right)}^2}{{\sum }_{i=1}^n{\left(y_i-\overline{y}\right)}^2}.}$$

In these equations, n refers to the total number of observations. The symbols y_i and $y_i^{\mathrm{\prime }}$ represent the actual and forecast values for the i-th data point, respectively, while y_i−1 and $y_{i-1}^{\mathrm{\prime }}$ correspond to the actual and forecast values from the previous time step (i−1). The term ȳ indicates the average of all actual values in the dataset. Lower values of MAPE and RMSE indicate better forecasting performance. A DA value closer to 1 suggests that the model more accurately captures the direction of change in comparison to the actual observations.

In this study, MAPE and RMSE indicate the magnitude of forecast errors, with lower values reflecting higher accuracy. DA measures the proportion of correctly predicted directional changes, which is important for timely policy responses to shifts in public attention. R² reflects how well the model explains variability in the observed data, indicating its explanatory power. Together, these metrics provide a balanced evaluation of forecasting precision, trend detection, and practical applicability for decision-making in educational policy.

Yearly variation analysis

The yearly variation in public attention to special education from 2020 to 2024 shows a steady upward trend overall (see Fig. 3). The annual average value increased from 1,222 in 2020 to 1,277 in 2022, with relatively modest growth rates of 4.42% and 0.78%, respectively, indicating a period of relative stability in public interest. A significant increase was observed in 2023, with the annual average rising to 1,413, representing the highest growth rate over the 5 years at 10.65%, possibly driven by policy initiatives or concentrated media coverage. In 2024, although the growth rate declined modestly to 6.94%, public attention to special education continued to stay at a relatively high level. This trend reflects a sustained increase in societal concern for special education in recent years, with a particularly notable surge beginning in 2023. Such a pattern offers a solid data foundation for developing forecasting models and the formulation of relevant policy interventions.

To examine the temporal pattern of public attention to special education between 2020 and 2024, an ordinary least squares (OLS) regression analysis was performed. The regression yielded a coefficient of 71.5, with an R² of 0.904 and a p-value of 0.0131, demonstrating a statistically meaningful upward trend at the 5% level. These results provide strong evidence of a consistent rise in public concern over the 5 years. Fluctuations in the annual growth rate may reflect the impact of various factors, including policy developments, major social events, and intensified media exposure. Future studies could delve deeper into these dynamics to better understand their effects and improve the ability to anticipate and respond to shifts in public engagement.

Monthly variation analysis

Figure 4 presents the monthly variation in public attention to special education in China during the period from 2020 to 2024. This forms a distinctive trend marked by a rapid increase in attention, immediately succeeded by a sharp drop, an annual rhythm that remains consistent throughout the observed period. A comparison of the monthly curves across different years shows that, despite varying annual contexts, public interest generally increases steadily from January to May, reaches its peak in June, declines during the July to September period, and in some years experiences a slight rebound after October. Notably, while the peak in 2020 occurred in July, the peaks from 2021 to 2024 consistently appeared in June, aligning with the timing of the national college entrance examination (Gaokao) in those years. These peaks also coincide with related events such as national policy announcements and the release of education-related reports, which may temporarily intensify public interest.

This structural fluctuation reflects the concentrated impact of specific events such as policy announcements, school entry assessments, and focused media campaigns. It also highlights the influence of public discourse and communication rhythms in driving attention. These observed trends provide a theoretical foundation for forecasting models such as SARIMA and SARIMA-LSTM, supporting the improvement of timeliness in educational communication and the scientific allocation of resources.

To investigate the trend over time, a regression model was fitted using the OLS approach. The results yielded an R² value of 0.008, a regression coefficient of 2.19, and a p-value of 0.495. These findings reveal that public attention has exhibited only a marginal upward trend over time, indicating a lack of statistical significance (p = 0.495 > 0.05). Public attention appears to vary without following a clear upward or downward linear pattern. Instead, it reflects fluctuations primarily driven by seasonal factors or specific annual events, rather than a monotonic temporal evolution.

Spatial distribution analysis

Based on Baidu Index data from 2020–2024, public attention to special education shows significant spatial disparities across Chinese provinces. Five provinces (Guangdong, Shandong, Zhejiang, Jiangsu, and Sichuan) recorded average values above 300, while regions such as Hainan, Ningxia, Qinghai, and Xizang remained below 100. Using the quartile method (Xu et al., 2024), attention values were classified into four levels: Level I (high) ranges from 262 to 399; Level II (relatively high) from 201 to 261; Level III (relatively low) from 166 to 200; and Level IV (low) from 15 to 165. Spatial distribution maps of attention levels across 31 provinces were generated using ArcGIS and Dycharts (see Fig. 5).

The map highlights pronounced disparities and evolutionary trends in public attention between 2020 and 2024. Higher levels are concentrated in the eastern and central regions, while lower levels are typical of western and peripheral provinces, reflecting the influence of economic development and population distribution.

Level I provinces are mainly found in the eastern coastal and central regions, including Guangdong, Jiangsu, Zhejiang, Shandong, Sichuan, and Henan. By contrast, Level IV regions such as Xizang, Xinjiang, Qinghai, and Ningxia remain at the lowest levels of attention. Provinces in Levels II and III, mostly in central and western China, show moderate levels with gradual improvement in some cases, such as Shaanxi. Sichuan stands out in the west, reflecting strong information sharing and active public engagement, while Henan, with its large population, shows a broad base of awareness. In comparison, Xizang and Xinjiang show little growth, which may be linked to smaller populations, lower internet access, and fewer educational resources.

Overall, the spatial pattern can be described as “higher in the east, lower in the west.” To narrow this gap, more efforts should be directed to western regions by improving internet access, distributing educational resources more fairly, and strengthening advocacy for special education. These measures would help promote balanced regional development of inclusive education.

Spatial correlation analysis

Based on Eqs. (1) and (3), the global spatial autocorrelation of public attention to special education across 31 provinces from 2020 to 2024 was calculated and tested for statistical significance (see Table 1). The results show that spatial correlation first increased and then declined. Moran’s I reached its highest values in 2021 and 2022 (above 0.22, p < 0.05), indicating clear positive spatial autocorrelation, but fell afterwards and was no longer significant in 2024. This points to a shift from regional concentration to a more balanced national distribution, likely supported by wider information channels and growing awareness campaigns.

Local Moran’s I analysis shows further changes in clustering (see Fig. 6). H–H clusters were mainly found in eastern and central provinces such as Shandong, Jiangsu, and Anhui, but by 2024 Jiangsu had dropped out of this group, suggesting weaker clustering as national attention became more balanced. L–L clusters continued in western provinces such as Xinjiang and Gansu, though Qinghai moved out of this category in 2024, showing some improvement in public awareness. Sichuan consistently displayed an H–L pattern, reflecting its strong influence compared with neighboring provinces. In contrast, Shanghai and Hainan at times showed L–H clustering, meaning lower attention than nearby regions.

Overall, both high-value and low-value clusters have contracted over the past 5 years, suggesting a gradual shift toward spatial balance in public attention. Nevertheless, regional gaps remain: western provinces continue to lag, while key provinces such as Sichuan maintain a leading role. Future policies should provide greater support to low-attention regions by expanding awareness programs and ensuring fairer access to resources, thereby promoting more balanced development of special education nationwide.

Stationarity test

The dataset from January 2020 to December 2023 was used for model training, while data from January to December 2024 served as the validation set for assessing forecasting performance. The ADF test was applied to assess the stationarity of the time series reflecting public attention to special education (see Table 2, Fig. 7). The test statistic was below the critical threshold at all conventional significance levels, with a p-value less than 0.05. This result supports the rejection of the unit root hypothesis, confirming that the series is stationary during the training period and does not require differencing.

Figure 7 illustrates the autocorrelation function (ACF) and partial autocorrelation function (PACF) of the time series. The ACF reveals a gradual decline over time, whereas the PACF displays a sharp drop after the initial lags. This implies that the series is influenced by both seasonal cycles and short-term dependencies. Together with the ADF test results, these patterns provide clear guidance for model parameter selection and confirm that the series meets the basic conditions for SARIMA modeling.

SARIMA model construction

An initial examination of the time series data on public attention to special education indicates that the prerequisites for applying the SARIMA model are met. As such, the differencing order d is set to 0, and the seasonal period s is defined as 12. The MAPE is used as the evaluation metric, as it is a measure that intuitively captures relative forecasting errors. A model with a smaller MAPE is considered to have superior forecasting accuracy and is therefore more suitable for future forecasting and analysis.

This study applied a grid search approach to systematically explore all possible parameter configurations, selecting values of P, D, Q, and s based on the minimum MAPE. The results indicate that the SARIMA (2, 0, 1) × (1, 0, 0)₁₂ model yields the lowest MAPE value and is thus chosen as the optimal model for forecasting public attention to special education.

The Ljung Box Q test was applied to the standardized residuals of the fitted model. At a lag of 12, the test produced a p-value of 0.890, which exceeds the 0.05 significance threshold. This outcome indicates that the null hypothesis, stating that the residuals follow a white noise process, cannot be rejected. As a result, there is no evidence of autocorrelation remaining in the residuals, suggesting that the model is correctly specified.

SARIMA-LSTM model construction

Although the SARIMA model effectively captures both trend and seasonal components, its residual series still exhibits certain non-linear fluctuations. To further improve forecasting accuracy, the residuals from the SARIMA model are treated as a new time series and used as input for non-linear modeling through an LSTM neural network. Table 3 summarizes the LSTM model settings.

Upon completion of model training, the LSTM network was used to forecast the monthly residuals for 2024. These forecast residuals were subsequently combined with the SARIMA model’s forecasts to produce the final SARIMA-LSTM forecasting results. The standardized residuals from the hybrid model were tested using the Ljung Box Q test at lag 12, yielding a p-value of 0.831. Given that this value is well above the 0.05 significance threshold, the result does not provide sufficient statistical grounds to refute the null assumption. This suggests that the error terms behave similarly to a white noise process, with no detectable signs of strong autocorrelation. As a result, the model can be considered well-specified, with no apparent structural deficiencies, and is suitable for reliable forecasting.

Figure 8 presents a comparative analysis of the 2024 forecast outcomes generated by the SARIMA-LSTM model and those from the standalone ARIMA and SARIMA models.

The Fig. 8 shows how the forecasted values correspond to the actual public attention levels in 2024, as generated by the ARIMA, SARIMA, and SARIMA-LSTM models. The SARIMA-LSTM model aligns most closely with actual public attention trends across most months, reflecting its strong performance in fitting the data. The SARIMA model also performs well, particularly in capturing seasonal variations, and generally yields better results than the traditional ARIMA approach.

Although the ARIMA model captures the overall direction of the data reasonably well, it tends to diverge significantly during months with sharp fluctuations, especially in June and December. In comparison, the SARIMA-LSTM model achieves higher accuracy by combining the strengths of both linear and non-linear modeling techniques. These findings underscore the value of deep learning in time series forecasting and suggest that SARIMA-LSTM offers a more effective framework for modeling complex temporal patterns in public behavior.

Model evaluation and validation

To evaluate forecasting performance, four indicators including MAPE, RMSE, DA and R² are used. These measures were applied to assess and compare the forecasting effectiveness of the three models (see Table 4).

The findings demonstrate that, among the three models evaluated, the SARIMA-LSTM approach yields the most favorable overall performance. It achieves the lowest forecasting error, with a MAPE of 6.94%, an RMSE of 121.53, a DA of 81.82%, and an R² of 0.79. These values suggest that the model not only provides a strong fit but also accurately tracks the trend in public attention.

In contrast, the ARIMA model demonstrates the lowest performance across all evaluation criteria, with an R² of –0.59, reflecting a limited capacity to explain the data’s variability. Although the SARIMA model performs better than ARIMA, particularly in its ability to capture seasonal structures and reduce forecasting error, its accuracy still falls short when compared to the SARIMA-LSTM model.

Overall, by integrating linear statistical modeling with the non-linear learning capacity of neural networks, the SARIMA-LSTM model proves to be more effective in handling the complexity of time series data like public attention, making it a promising tool for analyzing evolving behavioral patterns.

To further assess the forecasting reliability of the SARIMA-LSTM model, a forecasting interval coverage probability test was conducted, with the results presented in Fig. 9. The findings show that the model’s forecasts for 2024 closely match the actual observed values. Importantly, during periods of greater volatility, particularly in the early months of the year and toward year-end, the forecast values largely remain within the 95% confidence interval.

This outcome reinforces the stability and reliability of the SARIMA-LSTM model in capturing patterns of public attention. It demonstrates not only the model’s strong fit to the time series but also its ability to account for uncertainty and variability in future trends.