Henan Province is located in the middle and lower reaches of the Yellow River in the middle east of China and the south of the North China Plain, between 31°23′-36°22′N and 110°21′-116°39′E. Luoyang City is located in the western part of Henan Province, it is situated between longitude 112°16′-112°7′ and latitude 34°2′-34°5′, with a length of approximately 179 km from east to west and a width of approximately 168 km from north to south. The detailed geographical location of Henan Province is shown in Figure 1.

We collected daily concentration data of six air pollutants (PM_2.5, PM₁₀, NO₂, SO₂, O₃, and CO) in Luoyang from January 1, 2014 to July 30, 2022 and data from the Luoyang Air Testing Website (https://citydev.gbqyun.com/index/luoyang). There are ten air quality monitoring stations in Luoyang city, the daily concentration of air pollutants was simply an arithmetic mean measure across all the monitoring stations, as in most time-series studies. Daily mean temperature data in Luoyang were from the National Meteorological Information Center (http://data.cma.cn/).

Daily hospital admissions data were collected from the First Affiliated Hospital of Henan University of Science and Technology and The Third People's Hospital, which the two most representative hospitals in Luoyang. The clinical diagnostic criteria for Thyroid diseases from International Classification of Diseases. The patient's goiter was determined by ultrasound, and iodine deficiency patients were excluded. Patients' basic information included gender, age and residence.

A generalized additive model (GAM) with a Poisson distribution was adopted to analyze the impact of air pollutions on daily hospital admissions of goiter. The effect of different time lags was examined including eight single-day lags: (i) lag 0, the pre sent day; (ii) lag 1, the previous day; (iii) lag 2, the day before lag 1; (iv) lag 3, the day before lag 2; (v) lag 4, the day before lag 3; (vi) lag 5, the day before lag 4, (vii) lag 6, the day before lag 5, (viii) lag 7, the day before lag 6, and seven moving average exposure lags: (i) lag 01, the 2-day moving average of the present and previous day; (ii) lag 02, the 3-day moving average of the present and previous 2 days; (iii) lag 03, the 4-day moving average of the present and previous 3 days. (iv) lag 04, the 5-day moving average of the present and previous 4 days. (v) lag 05, the 6-day moving average of the present and previous 5 days. (vi) lag 06, the 7-day moving average of the present and previous 6 days.

GAM is a flexible regression analysis method that allows for the inclusion of non-linear relationships in the model and describes these relationships through non-parametric smooth functions (18). Firstly, we established a basic model that includes the long-term trend of time, the day of the week effect, and the holiday effect, as follows:

Among them, Y_t is the actual number of people in the hospital on the t-th day (following the Poisson distribution); μ_t = E[Y_t]is the expected number of hospitalizations on the t-th day; β is the regression coefficient; X_t is the atmospheric pollutant element on the t-th day; S represents a non-parametric smoothing function, where df_t represents the degree of freedom of the long-term and seasonal trends of time;

In order to capture short-term fluctuations, the model incorporates holiday variables (hol) and week variables (dow); Holiday variables are simplified into binary categories, where hol = 1 represents holidays and hol = 0 represents non-holidays. The α in the model is the intercept term. In order to better capture the relationship and trend of changes between independent variables, cubic spline smoothing is used to map the discrete values of the independent variables to a continuous function, so that the model can better fit the data, i.e. k = 3; The degree of freedom of the annual non-parametric smoothing function for time is set to 7, a natural spline function with 6 degrees of freedom for daily mean temperature (19).

The relative risks (RR) with 95% confidence interval (CI) in thyroid diseases admissions associated with a 10.0 μg/m³ increase in daily concentration of NO₂, SO₂, and O₃, and a 1.0 mg/m³ increase in daily concentration of CO were estimated, RR was calculated using the following formula (20):

where β is the regression coefficient from the GAM model, SE is the standard error.

Implementing stratified analysis to explore the impact of environmental pollutants on the risk of thyroid tumors in different subgroup variables included gender (male and female), age (<45, 45–65, and >65 years), and season (Spring, Summer, Autumn, and Winter).

By adjusting the degrees of freedom of the time smoothing function, we can monitor the sensitivity of the estimated effect values to changes in degrees of freedom and determine the stability of the model (21). R software (Version 3.2.3) was used to perform analysis, all the statistical tests were two-tailed and a P < 0.05 was considered statistically significant.

There were 37,630 hospital admissions for goiter in Luoyang from January 1, 2014 to July 30, 2022. Among them, there are 29,571 female (78.58%), 8,059 male (21.42%); 14,141 people (37.58%) from rural areas and 23,489 people (62.42%) from urban areas; The average age of the patient is 51 years old, 15,349 (40.79%) people are under 45 years old, 17,498 (46.5%) people are between 45 and 64 years old, and 4,783(12.71%) people are over 65 years old.

Analysis of the air pollutants indicated the daily mean concentrations were 62.92 μg/m³ for PM_2.5, 109.34 μg/m³ for PM₁₀, 39.26 μg/m³ for NO₂, 24.47 μg/m³ for SO₂, 96.97 μg/m³ for O₃, and 13.37 mg/m³ for CO (Table 1). The time series analysis of pollutant concentration suggests a decreasing and Seasonal trend in PM_2.5, PM₁₀, SO₂, and CO from 2014 to 2022 (Figure 2).

We calculated Spearman correlation coefficient (r) to examine the relationships of air pollutions (Figure 3). The results indicated that daily PM_2.5 and PM₁₀ concentrations had positive correlations with NO₂ (PM2.5: r = 0.67, P < 0.001; PM₁₀: r = 0.71, P < 0.001), SO₂ (PM_2.5: r = 0.4, P < 0.001; PM₁₀: r = 0.43, P < 0.001), and CO (PM_2.5: r = 0.65, P < 0.001; PM₁₀: r = 0.61). O₃ and other atmospheric pollutants show a negative correlation (r < 0, P < 0.001).

The impact of six pollutants had different lag effects on goiter (Table 2). The lag effect for PM_2.5 and CO (lag1–3days) was significant and relatively longer for NO₂ (lag0–3) and O₃ (lag2–7). For PM₁₀, the lag effect was significant only at lag1 and lag3. The accumulated lag effect for PM_2.5 (lag01–06days), PM₁₀ (lag01–03days), O₃ (lag04–07), NO₂ (lag00–07) and CO (lag01–05days) was significant and relatively longer. Besides, the lag effect for SO₂ on goiter was not statistically significant.

Specifically, a 10μg/m³ increase in PM_2.5, PM₁₀, NO₂ and a 1μg/m³ increase in CO was associated with 0.92% (RR: 1.0092; 95%CI: 1.0032–1.015),0.44% (RR:1.0044; 95% CI: 1.0008–1.0081), 5.96% (RR: 1.0596; 95% CI: 1.0413, 1.0783) and 62.4% (RR: 1.624; 95% CI: 1.1347, 2.3243) increased risk of goiter on the optimal lag day.

The dose-response relationship between pollutants concentration and goiter were shown in Figure 4. As the concentration of PM_2.5 and PM₁₀ increases, the risk of goiter showed a near linear increase trend, and an exponential increase trend with NO₂ and SO₂ increase. O₃ exhibits a protective effect against goiter, and as O₃ concentration increases, the risk of goiter decreases continuously. As the concentration of CO increases, the risk of disease shows a trend of first increasing and then slowly decreasing.

In gender stratification, PM_2.5 and CO have statistical significance for the female population, while NO₂ has statistical significance for both males and females. Women are more sensitive to NO₂ than men. When NO₂ concentration increases by 10 μg/m3, the RR of goiter in female and male populations are 1.065(1.046,1.084) and 1.042(1.006,1.079), respectively; When the concentration of PM_2.5 and CO increases by 10 μg/m³ and 1 mg/m³ separately, the RR of goiter in the female population are 1.01(1.003,1.016) and 1.937 (1.179,3.183) (Figure 5).

In different age groups, PM_2.5, PM₁₀, and CO have significant effects on the population aged 45–65. When the concentration of pollutants increases by 10 μg/m³ (CO increases by 1mg/m³), the RR of goiter increased by 1.48%, 0.61%, and 3.1609 times, respectively; NO₂ has statistical significance in all age groups. When NO₂ concentration increases by 10 μg/m³, the RR of goiter in <45 years old, 45–64 years old, and >65 years old age groups increased by 3.96%, 7.79%, and 8.8%, respectively (Figure 5).

PM_2.5 had a significant effect on goiter in Summer. When the pollutant concentration increases by 10 μ g/m³, the RR of goiter increased by 4.87%; However, PM₁₀, NO₂, and SO₂ have a more significant impact on goiter in autumn. When the pollutant concentration increases by 10 μ g/m³, the RR of goiter increased by 0.87%, 6.38%, and 12.44%, respectively; In winter, NO₂ has a more significant impact on goiter. When the concentration of pollutants increases by 10 μg/m³, the RR of goiter increased by 5.48% (Figure 5).

There is an interactive effect between pollutants on diseases. Except for O₃, other pollutants have a certain synergistic effect on thyroid goiter. Adding PM₁₀ to PM_2.5 has the significant impact on goiter, with a 1.0111% (95% CI: 1.0051, 1.0171) risk for goiter; Adding NO₂ to PM_2.5 with a 1.0701% (95% CI: 1.0512, 1.0894) risk for goiter; Adding PM_2.5, PM₁₀, O₃, NO₂, SO₂ separately to CO may increase the risk of thyroid goiter by 2.5746 times(RR:2.5746, 95% CI: 1.4894–4.4505), 2.34 times(RR:2.34, 95% CI: 1.3845–3.955), 2.1893 times(RR:2.1893, 95% CI: 1.3383–3.5815),1.6727 times (RR:1.6727, 95% CI: 1.0053–2.7832) and 2.2812 times (RR:2.2812, 95% CI: 1.3908–3.7416). However, O₃ exhibited antagonistic effects against other pollutants.

The GAM model was used to fit each pollutant data separately, and all smoothing terms reached significance at the p < 0.05 level (Table 3). After adjusting the interference factor temperature in the model, the relative risk has decreased slightly, but it is still statistically significant. The adjusted R² values of each model in Table 3 are around 80%, indicating a good fitting effect of the GAM model.

The degrees of freedom for selecting time are 5, 6, 7, 8, and 9, respectively. The effects of PM_2.5, PM₁₀, NO₂, and CO on the risk of goiter are statistically significant, and according to the increase in degrees of freedom, it can be seen that there is little change in the RR with 95% confidence interval, indicating that the model is relatively stable (Table 4).