Monthly ST cases were collected for the top eight high-risk ST counties of China for the 2012–2020 period from the China Information System for Disease Control and Prevention, including Yingshang County (Anhui province), Guangning County (Guangdong province), Huaiji County (Guangdong province), Yingde City (Guangdong province), Longling County (Yunnan province), Gengma County (Yunnan province), Mang city (Yunnan province), and Yingjiang County (Yunnan province), as shown in Figure 1. The total cases of the eight counties over the time period of analysis were 2,177, 2,832, 3,395, 2,688, 4,052, 1,936, 1,952, and 2,073, respectively. The criteria for a confirmed ST case included epidemiological exposure patient histories (travel to an epidemic area and contact with chiggers or rodents within 3 weeks before the onset of illness), clinical manifestations (for example, skin rash, lymphadenopathy, high fever, and eschar or ulcers), and also positivity for at least one of the laboratory diagnostic criteria/tests: isolation of O. tsutsugamushi from clinical specimens, or detection of O. tsutsugamushi by polymerase chain reaction (PCR) in clinical specimens, or a 4-fold or greater rise in serum IgG antibody titers between acute and convalescent sera by using indirect immunofluorescence antibody assay (IFA) (Zhang et al., 2013; Li et al., 2020).

Environmental data used in this study included raster maps of the normalized difference vegetation index (NDVI) representing the amount of vegetation at specific locations which was collected from the MODIS-Terra products (MOD13A2, https://modis.gsfc.nasa.gov/data/dataprod/mod13.php) with spatial resolution 1 km. We also used meteorological data from weather monitoring stations in China from the China Meteorological Administration (http://www.cma.gov.cn/), including precipitation, pressure, relative humidity, sunshine duration, daily mean temperature, daily minimum temperature, daily maximum temperature, and mean wind speed. The multivariate El Niño/Southern Oscillation Index (MEI), regarded as global climate change proxy, was collected from the physical sciences laboratory of National Oceanic and Atmospheric Administration (https://psl.noaa.gov/enso/mei/).

The inverse distance weighted method was employed for mapping the six meteorological data with the same spatial resolution of the NDVI. Then, the NDVI and the six meteorological data were extracted by the administrative boundaries of the eight counties, and the mean values of each variable were calculated each month during the study period for further analysis.

Based on the standard ARIMA model, the SARIMAX model considers simultaneously the seasonal variation in ST cases and accounts for the effects of exogenous risk factors for better understanding and fitting the considered time series. The basic equations of SARIMAX are as follows: Y t = ∑ β i X i , t + Z t , (1) Φ P ( B S ) ϕ p ( B ) ∇ S D ∇ d Z t = Θ Q ( B S ) θ q ( B ) ε t , (2) where Y t represents the ST cases at time instant t , while X i , t represents the i t h exogenous variables ( i = 1,2 , … , 8 ) at time instant t   , and β i represent the coefficients of the exogenous variables; Z t representing the main trend of the ST time series satisfies the SARIMA equation (Eq. 2) based on the model structure ( p , d , q ) × ( P , D , Q ) S . In Eq. 2, S represents the periodicity, and ε t denotes the white noise; B represents the backshift operator, e.g., B i Z t = Z t − i , while ∇   represents the differencing process, e.g., ∇ d = ( 1 − B ) d , ∇ S D = ( 1 − B S ) D with the non-seasonal and seasonal differencing orders d and D, respectively; ϕ p ( B ) = 1 − ∑ i = 1 p ϕ i B i and Φ P ( B S ) = 1 − ∑ i = 1 P Φ i B i S represent the non-seasonal and seasonal autoregressive process with the orders p and P, respectively; θ q ( B ) = 1 − ∑ i = 1 q θ i B i and Θ Q ( B S ) = 1 − ∑ i = 1 Q Θ i B i S represent the non-seasonal and seasonal moving average processes with the orders q and Q, respectively. In the current study, the periodicity parameter S was set to 12 months. The main procedure of defining the structure of SARIMAX is briefly described as follows: 1) SARIMA models with the parameters p, P, q, and Q ranged from 0, 1, and 2, and the parameters d and D ranged from 0 and 1 were constructed to fit the time series of ST cases in each of the eight studied counties. The model with the lowest Akaike information criterion (AIC) value was regarded as the best SARIMA model. 2) Multicollinearity was investigated before entering the exogenous variables into the SARIMAX model. We found a very high level of correlation (0.9834 and 0.9912, respectively) between mean temperature and daily minimum temperature (or daily maximum temperature), and therefore, the daily mean temperature was selected for modeling in the current study. In our study, eight exogenous variables with a 0–5 month lag were prepared in an exogenous variable pool. Based on the SARIMA model, a forward stepwise variable selection procedure was implemented to build the SARIMAX model. During each of the variable selection loops, multicollinearity was tested before adding the exogenous variable into the model, and only the one with a variance inflation factor (VIF) smaller than 3 was regarded as a candidate. Furthermore, the variable candidates were added into the model separately, and only the one with the lowest AIC in each loop was added into the model. Finally, a t-test was used to exclude the exogenous variables with largest p-values (larger than 0.05). More detailed information can be found in the literature (He et al., 2018b). The first eight-year ST cases were employed to fit the SARIMAX model, and the last year data were used to test the performance of the built model. To evaluate the goodness of the fit of the model, we considered the R², mean absolute error (MAE), and root mean square error (RMSE).

In order to test the possibility of applying LSTM in ST modeling and forecasting, the exogenous variables included in the final SARIMAX model together with the ST time series were chosen as the LSTM input. LSTM is constructed by a number of connected cells (the basic unit of the network), while each cell consists of three gates, i.e., the forget gate, input gate, and output gate; see Figure 2. Specifically, the hidden state from the previous cell h t − 1 and the current series Y t at time instant t are combined and flow through the forget gate, leading to useless information of loss by f t = σ ( W f ⋅ [ h t − 1 , Y t ] + b f ) ; in other words, the forget gate can retain the useful information in series modeling. The input gate gathers the hidden state and current series to update the cell situation C t of LSTM, and two preparation workflows are designed, i.e., i t = σ ( W i ⋅ [ h t − 1 , Y t ] + b i ) and C ˜ t = t a n h ( W C ⋅ [ h t − 1 , Y t ] + b C ) , and the cell situation can be calculated by C t = f t ∗ C t − 1 + i t ∗ C ˜ t ; in other words, the input gate collects and updates the information flow into the model. The output gate uses the hidden state and current series to update the hidden state as the input of the next cell of LSTM, i.e., o t = σ ( W o ⋅ [ h t − 1 , Y t ] + b o ) and h t = o ∗ tanh ( C t ) ; in other words, the output gate integrates the information and generates the output information for the next time instant. In these equations, σ represents the sigmoid activate function, W and b represent the weights and bias in different parts of LSTM, respectively. Other than some other advanced LSTM models, such as forward and backward variate sensitive LSTM, convolutional LSTM, and convolutional neural network LSTM (Kim and Cho, 2019; Wan et al., 2019; Fouladgar and Främling, 2020), the current study used the basic LSTM as described above to forecast the ST series. To define the best structure of the LSTM model, the values of the hidden layers were considered to vary from 2 to 3, while the hidden dimension varies from 54 to 68 with interval 4, the batch size varies from 45 to 55 with interval 5, and the number of previous months used to forecast the current month ranges from 3 to 12. The optimal model was selected with the largest R² and smallest MAE and RMSE.

Various SARIMA and SARIMAX models were established with different model structures, and temporal lagged exogenous variables at the eight considered counties and the optimal model structure with the corresponding performance are presented in Supplementary Appendix Table A1. According to the optimal model structures, all models have the same seasonal difference order, i.e., D = 1 . Specifically, the models of the four counties in the Yunnan province have very similar seasonal characteristics, i.e., P, D, and Q are exactly the same, except the Q value for Yingjiang County. Regarding the SARIMA model, it showed good performance in modeling ST variations in the model fitting stage in Guangning County, Longling County, and Yingjiang County with R² larger than 0.8, while the SARIMA model showed better performance (in terms of R²) at the model forecasting stage in Longling County, Gengma County, and the Mang city of the Yunnan province (Table 1).

Through exogenous variable selection processes, the SARIMAX models absorb the strength of one or several exogenous variables with temporal lags for ST variation modeling. Atmospheric pressure with a 4-month temporal lag was found to be the sole exogenous variable for SARIMAX modeling in Huaiji County. As shown in Supplementary Appendix Table A1, the significant ST-associated exogenous variables vary between the eight counties under investigation. Atmospheric pressure, sunshine duration, wind speed, and MEI with various temporal lags were found to be correlated with the ST temporal variation in Yingyang County and Guangning County; however, precipitation and relative humidity were also included in ST forecasting at Yingyang County, while NDVI was included in the model at Guangning County. At Gengma County and Yingjiang County, relative humidity and sunshine duration were significant variables associated with ST; specifically, mean temperature was also found to be associated with ST at Yingjiang County. In addition, mean temperature and wind speed were two important impact factors of ST at Yingde city and Mang city, but precipitation and MEI were considered as other important factors for the two cities, respectively. Finally, the ST variation at Longling County was closely correlated with precipitation, relative humidity, and mean temperature with various temporal lags. Compared to the SARIMA model, the results of the SARIMAX model showed that more accuracy performance (in terms of R²) can be achieved at the model fitting stage in all of the eight counties; meanwhile, at the model forecasting stage, only the SARIMAX models in Yingshang County and Yingjiang County showed better forecasting accuracy than the SARIMA models in terms of RMSE and MAE.

By inputting the ST series and the corresponding significantly related exogenous variables shown in Supplementary Appendix Table A1 to the LSTM model, the optimal model structure was obtained according to the smallest R² in the model fitting stage. The results showed that better performance can be yielded by setting the length of input series as 12 months at Huaiji County, Yingde city, Longling County, Gengma County, and Yingjiang County, while the optimal length of input series at Yingshang County, Guangning County, and Mang city were 10, 3, and 11, respectively. The hidden dimension, batch size, and number of layers varied between counties. At the model fitting process, the LSTM model showed better performance than SARIMA at Yingshang County and Longling County in terms of R²; at the model forecasting stage, LSTM showed more accuracy in forecasting ST cases than SARIMA and SARIMAX at Yingshang County and Yingde County, while only than SARIMAX at Guangning County and Longling County (Table 1).

The model fitting and forecasting results were separated by a vertical dash line, as shown in Figure 3. The gray line, yellow line, red line, and green line represent the real ST notification series, SARIMA estimation, SARIMAX estimation, and LSTM estimation, respectively. It can be concluded as follows: 1) the LSTM model used a smaller length of input data for ST fitting than SARIMA and SARIMAX because the green line begins much earlier than the yellow and red lines; for example, most of the green lines at the eight counties begins during 2013 except the Guangning County, while the yellow and red lines begin during 2014 for Yingshang County, Huaiji County, Yingde County, and Yingjiang County, and for the other four counties, the model estimations begin during 2015. 2) Compared to the SARIMA and SARIMAX, LSTM showed smoother estimation values and less local fluctuations.3) Through comparing the real ST time series in the eight counties, four kinds of temporal patterns of ST variations can be concluded as follows: the ST series peaks in October and August in Yingshang County and Longling County, respectively; two peaks can be found in June and October in Guangning County, Huaiji County, and Yingde city, while in Gengma County, Mang city, and Yingjiang County, bio-peak phenomena can also be found in summer and winter during multiple years of the study period.

Although the performance of the three time-series forecasting methods varied between counties at both the model fitting stage and the model forecasting stage, the methodological comparisons among SARIMA, SARIMAX, and LSTM should be discussed as follows: 1) continuous series modeling was the core concept of the three methods, i.e., the previous conditions can be used to forecast the current or future conditions; furthermore, the SARIMA and SARIMAX models took a step ahead that included seasonal or cyclic parts in modeling, i.e., the seasonal auto regression part and seasonal moving average part with the cycle of 12 months. Similarly, the similar SARIMA model was employed for modeling the ST series at Laiwu city and Shandong province, China (Ding et al., 2012a; Yang et al., 2015). On the other hand, LSTM also borrows the continuous characteristics of ST series for forecasting. Compared to the standard artificial neural network (ANN), the forget gate was utilized to remove the information with larger lags, which is similar to the memory of humans that one can remember the recent things but will forget some of the things that occurred long time ago (Yu et al., 2019). In the current work, we have demonstrated that the length of the time-series data used to forecast ST varies between high-risk counties being 12 months in five counties, 11, 10, and 3 months in the other three counties. These results indicate that ST has predominantly an annual cycle, which is in line with the literature that 8–12 months is the general cycle of ST in Guangzhou city, China (Wei et al., 2017). 2) Previous studies have used a number of regression techniques to unravel the ST temporal variation in incidence, and its associations with environmental factors using the spatial Poisson regression model, negative binomial regression model, MaxEnt, random forest, and ANN models (Wardrop et al., 2013; Li et al., 2014; Yang et al., 2014; Kwak et al., 2015; Seto et al., 2017; Yu et al., 2018; Acharya et al., 2019). However, these models cannot capture the seasonal characteristics of ST incidence, which can be overcome by the SARIMAX model that can simultaneously depict the seasonal characteristics and the effects of exogenous factors. Given the complex and often nonlinear interplay between ST factors and ST infection (Elliott et al., 2019), the linear regression component of the SARIMAX model may not be sufficient for modeling and forecasting ST incidence. The LSTM model offers a self-organizing, self-learning, and highly nonlinearity system from the ANN, and thus, it is rather suitable for modeling the ecological system of ST. In view of its model structure, the LSTM was initially postulated as the optimal model among the three methods. However, our results indicate that LSTM not always outperforms SARIMA or SARIMAX indicating that the selection of methods for ST series modeling and forecasting should be carefully considered through model comparison and validation to local conditions.

Another outcome of the current work is the determination of environmental factors in each of the considered eight counties associated with ST cases using SARIMAX models. Including these factors (i.e., exogenous variables mentioned above), the SARIMAX showed better performance than SARIMA at the modeling fitting stage (Table 1), demonstrating that the considered factors were significantly associated with ST variations. Temperature with various temporal lags was found to be positively related to ST variations at Yingde city, Longling County, Mang city, and Yingjiang County, which is in line with the ecological niche modeling of ST in the literature that temperature was found to be the key factor in determining ST occurrence (Yu et al., 2018). A previous study found that the sunshine duration was negatively and positively associated with ST with 1–3 month lags and 4–6 month lags, respectively (Yang et al., 2014). Similar findings can be found in the current study, i.e., sunshine duration with 3-, 2-, and 1-month lags were negatively associated with the ST variations at Yingshang County, Guangning County, and Gengma County, respectively, but the sunshine duration with 0-month lag showed positive relationship with ST at Yingjiang County. Warming the environment promotes the growth of vector larvae and rodents and increased exposure opportunities due to people wearing shorter clothes (Yao et al., 2019). The other study found that the ST case was positively correlated with the duration of sunshine, suggesting an occupational exposure where people possibly have longer time for agricultural field work, leading to an increased probability of exposure (Li et al., 2014). Besides, relative humidity and pressure were supposed to be key factors that influenced the regeneration of the rodents, i.e., a 1–2 month lag effects of relative humidity and pressure were closely related to the ST variation (Sun et al., 2017), which is similar to the findings of the results at Yingshang County, Guangning County, Huaiji County, Longling County, and Gengma County. In our study, we considered the wind speed in our SARIMAX models with 1-, 2-, 3-, 5-month lag since it can be a factor associated with the spawning conditions of mites (Kwak et al., 2015). Indeed, we found that wind speed was a contributor to ST variation at four counties of the current study, especially in Yingde city. Finally, our results indicate a significant effect of the multivariate ENSO index on the time series of ST cases in three of the eight counties investigated (i.e., Yingshang County, Guangning County, and Mang city). This is in line with other studies that found an association between ST incidence and the ENSO Index which is regarded as a global pattern of climatic oscillation affecting the local environment and thus human population behaviors (Wei et al., 2017).

Machine learning techniques have become popular in modeling nonlinear systems, including in the field of public health (dos Santos et al., 2019; Panch et al., 2018), by providing precise predictive models (Passos et al., 2016). Introduction of LSTM has enriched the ability of local public health authorities in ST high-risk areas to anticipate and respond to ST outbreaks. The proposed models, such as SARIMAX and LSTM, can be used locally in high-risk ST counties for ST early warning and precise ST forecasting programs; it can also enable local public health managers to monitor the variation of the environmental factors and deploy public health measures, such as health promotion alerts to the communities to prevent large ST outbreaks.

Certain limitations of the current work should be mentioned. First, due to lack of a longer time series of ST series data, this study only utilized a 96-month time series for model fitting, which may have hindered the performance of LSTM models; it may be the reason why LSTM showed low accuracy in ST estimation. Second, this study explored the basic LSTM model in structure and did not explore modeling combinations that considered SARIMAX + LSTM, which can also give accuracy predictions (Sheng and Jia, 2020); hence, future studies should explore the possibility of this combination for ST forecasting. Third, our SARIMAX models did not consider socioeconomic factors (including the gross domestic product, income, urbanization, population density, educational institutions, land use and land change, and medical institution), which may also play a role in the variations of ST (Ranjan and Prakash, 2018). Finally, although it is reported the warming condition might favor the reproduction of mites and increase the probability of human infection, the mechanism between the global climate change and ST outbreaks is still unclear (Jeung et al., 2016; Kuo et al., 2015); given the short-term impact of temperature mentioned in previous subsection, the long-term impact of global climate change (such as the warming condition) on ST outbreaks and its feasibility in forecasting ST outbreaks are worthy to explore in the future (Zhou et al., 2021); furthermore, the hydrology impact of climate change on ST variations will be another topic for consideration (Zhou et al., 2018).