Study area
Taizhou City, a coastal city in the central part of Zhejiang Province, belongs to the mid-subtropical monsoon area and experiences four distinct seasons (Supplementary Figure S1). The territory experiences mild summers, cold winters, abundant rain, and a mild, humid climate due to the meteorological effects of nearby ocean waters and mountains in the northwest.
Data collection
According to the Law on Prevention and Treatment of Infectious Diseases, HFRS is classified as a Class B infectious disease in China, and cases must be reported within 24 h of diagnosis [19]. Data on HFRS from 2008 to 2020 in Taizhou City were collected from the Chinese Notifiable Disease Reporting System.
We collected daily meteorological data from the China Meteorological Data Sharing Service System (http://data.cma.cn/). These data, including daily average temperature, (Avetemp), minimum temperature (Mintemp), maximum temperature (Maxtemp), relative humidity, and total precipitation, were used to calculate the weekly average for each value.
Statistical methods
Normality test and descriptive analysis were conducted to summarize characteristics of all variables. Spearman correlation was used to assess the relationship between HFRS incidence and meteorological factors. This study developed a time series model based on the GAM and used the cross-basis process to describe the distribution of changes in the independent variable dimension and the lag dimension simultaneously [8]. Further, DLNM was used to fit the non-linear and lag effects of weekly Avetemp, Maxtemp, Mintemp, average relative humidity, and cumulative rainfall on the risk of HFRS [20]. The incubation period for HFRS is affected by the host animal, vector density, and meteorological factors, and lasts for several weeks. In our study, the maximum lag period was set to 16 weeks [13, 20, 21]. Since HFRS cases in Taizhou City were relatively rare, Quasi-Poisson regression was used in this model to control for overdispersion. We used a two-stage analysis method. First, we used the DLNM to estimate the association of weekly Avetemp, Maxtemp, Mintemp, relative humidity, and weekly total precipitation (WTP) with the number of HFRS case [22]. The general algebraical definition of the model are as follows:
$$Log[E(Yt)] = \beta + cb(Kt,16,\beta 1) + S1(x) + S2(z)+S3(m)+S4(n)+ S5(week)$$
Among them, t is the observation week; [E(Yt)] is HFRS cases observed in month Yt, βis the intercept of the entire model; cb (Kt, 16, and β1) is the cross-basis function of K, and K is one of the meteorological elements. S1–4 are factors such as Avtemp, Maxtemp, Mintemp, RH, and WTP; β1 is the estimated value of the effect of K in a specific lag week t; the maximum lag week is set to 16; Week is the ordinal variable of the week in the year; s() is the penalty spline function. This study uses cubic spline functions, s1–4, to adjust the confounding factors in the model, and s5(week) to adjust the weekly confounding factors.
Second, we analyzed the interaction between weekly Avetemp, Maxtemp, Mintemp, relative humidity, and accumulated rainfall with GAMs, and then analyzed the different effects of high and low values of the meteorological factors on the cases. The basic model as follows:
$$Log[\mathrm{E}(\mathrm{Yt})] =\upbeta 2 +\mathrm{ S}1(\mathrm{k},\mathrm{x}) +\mathrm{ S}2(\mathrm{z}) +\mathrm{ S}3(\mathrm{m}) +\mathrm{ S}4(\mathrm{n}) +\mathrm{ S}5(\mathrm{week})$$
β2 is the intercept; K is one of the meteorological factors (Avtemp, Maxtemp, Mintemp, RH, and WTP), and X, Z, M, and N denote the other factors. s() means the penalty spline function. s1 (K and X) is the spline function of the interaction between variables K and X.
In the model, the number of cases was used as the dependent variable, and a cross-basis function was established for the number of cases and temperature. Spline interpolation was used to control the influence of confounding factors such as relative humidity, rainfall, and long-term trends. The best degree of freedom (df) was selected based on the spline function results through sensitivity testing and generalized cross-validation criteria [23].
DLNM can describe the complex nonlinearity and hysteresis correlation of temperature-HFRS through the cross basis function. It is necessary to scientifically define the reasonable lag time of the model [21]. We chose one ns (natural cubic B-spline, df = 6) as the exposure–response. Two nodes are located at the 2.5th and 97.5th percentiles of the meteorological factor distribution, and the other is for the exposure–response relationship, based on high temperature [13, 15]. The assumption that Maxtemp and Mintemp may affect the incidence of HFRS is fixed (df = 6) and the maximum lag time was set at 16 weeks to capture the delayed effects of extreme temperatures. In this study, the degrees of freedom and maximum lag times for mean temperature, relative humidity, accumulated rainfall, mean maximum temperature and mean minimum temperature were set in order: (df = 6, lag = 16; df = 4, lag = 16; df = 3, lag = 10; df = 6, lag = 21; df = 4, lag = 20).
We performed sensitivity analysis by changing the df of the weather variables and time points. All analyses were performed using ArcGIS 10.2 (ESRI, Redlands, CA, USA) and R software (packages "dlnm" and "mgcv") (R Foundation for Statistical Computing, Vienna, Austria).