Non-linear association between long-term air pollution exposure and risk of metabolic dysfunction-associated steatotic liver disease

Abstract

Background: Metabolic Dysfunction-associated Steatotic Liver Disease (MASLD) has become a global epidemic, and air pollution has been identified as a potential risk factor. This study aims to investigate the non-linear relationship between ambient air pollution and MASLD prevalence.

Method: In this cross-sectional study, participants undergoing health checkups were assessed for three-year average air pollution exposure. MASLD diagnosis required hepatic steatosis with at least 1 out of 5 cardiometabolic criteria. A stepwise approach combining data visualization and regression modeling was used to determine the most appropriate link function between each of the six air pollutants and MASLD. A covariate-adjusted six-pollutant model was constructed accordingly.

Results: A total of 131,592 participants were included, with 40.6% met the criteria of MASLD. “Threshold link function,” “interaction link function,” and “restricted cubic spline (RCS) link functions” best-fitted associations between MASLD and PM_2.5, PM₁₀/CO, and O₃ /SO₂/NO₂, respectively. In the six-pollutant model, significant positive associations were observed when pollutant concentrations were over: 34.64 µg/m³ for PM_2.5, 57.93 µg/m³ for PM₁₀, 56 µg/m³ for O₃, below 643.6 µg/m³ for CO, and within 33 and 48 µg/m³ for NO₂. The six-pollutant model using these best-fitted link functions demonstrated superior model fitting compared to exposure-categorized model or linear link function model assuming proportionality of odds.

Conclusion: Non-linear associations were found between air pollutants and MASLD prevalence. PM_2.5, PM₁₀, O₃, CO, and NO₂ exhibited positive associations with MASLD in specific concentration ranges, highlighting the need to consider non-linear relationships in assessing the impact of air pollution on MASLD.

Introduction

Non-alcoholic fatty liver disease (NAFLD) is a complex systemic inflammatory and metabolic disease majorly defined by its liver manifestation [1, 2]. The updated definition, metabolic dysfunction-associated fatty liver disease (MAFLD), was established to emphasize on the pathophysiological nature of metabolic dysfunctions [3] and related adverse clinical outcomes [4, 5]. The prevalence of fatty liver/steatotic liver disease has been growing rapidly over time [6–8], making it a major contributing factors to the loss of life-years worldwide [9]. Thus in 2023, the nomenclature and the diagnostic criteria were updated by three large pan-national liver associations. The newly proposed metabolic dysfunction-associated steatotic liver disease (MASLD) had much clear definitions for epidemiological studies [10].

In addition to lifestyle and behavior risk factors such as alcohol, smoking and dietary imbalance, air pollution has recently been recognized as one of the potential risk factors for liver diseases, particularly in complex liver diseases like MASLD [11, 12]. Prior studies have indicated that nitro dioxide (NO₂) is associated with increased oxidative stress and biomarkers of hepatic inflammation [11, 13]. Additionally, particulate matter (PM), especially those with diameter smaller than 2.5 µm (PM_2.5), can penetrate the alveoli and trigger systemic inflammatory responses [14–16]. Long-term exposure to elevated level of PM_2.5 have been found to be associated with insulin resistant, metabolic syndrome, and possibly MASLD [17–20]. However, two recent large-scale studies have suggested a possible non-linear relationship between air pollutant exposure and the risk of MASLD. In a landmark study by Guo and colleagues, a non-linear association was reported between long term exposure to elevated PM_2.5 concentration and the odds ratio (OR) for MAFLD [17]. Similarly, Sun et al, used a piecewise analysis with break points to model the hazard ratio for NAFLD in relation to PM_2.5 [18]. Furthermore, no previous reports have explored the potential risk of MASLD associated with the combined effects of commonly monitored air pollutants. The objective of this study was to investigate the existence and form of non-linear relationships between six regularly monitored ambient air pollutants, individually and collectively, and prevalence of MASLD.

Materials and methods

Study design and covariates

A cohort of paid health checkup participants from the MJ clinic, were included in this cross-sectional study. The study was approved by the Institutional Review Board of National Cheng Kung University Hospital, Tainan, Taiwan (Approval Number: B-ER-110-456) with a waiver of informed consent. Participants provided detailed demographic, lifestyle and medical history information through self-administered questionnaires [21]. Socio-demographic characteristics and lifestyle were obtained from self-administered questionnaires on the day of health checkup, including vegetable, fruit, sugar drink and fried food intake amount and frequency (seldom, moderate and frequent intake), and cigarette smoking status (never, former or current). As for alcohol drinking, those who consumed two or more alcoholic drinks per day on three or more days a week for more than 1 year [22] (ie, evidently took more than 70 to 140 gm alcohol intake per week) were considered as having excessive alcohol drinking. Participants with excessive alcohol consumption were excluded in current study, as hepatic steatosis in this group may be related to alcohol, they should be classified as MASLD with increased alcohol intake (MetALD) or alcohol-related liver disease (ALD) [10]. This distinction is crucial because the disease course in ALD or MetALD may differ significantly from MASLD.

Regular exercise was defined as taking at least half an hour of exercise once per 2–3 days. Anthropometric measurements and blood tests were conducted and reported by standard protocols of MJ Health Management Institution. The fasting glucose and lipid profile were taken under fasting state. Diabetes mellitus (DM) was defined as a fasting serum glucose ≥126 mg/dL, or HbA1c ≥6.5% [23] or self-reported history. Hypertension (HTN) was defined as grade 1: systolic blood pressure (BP) between 130 to 139 mmHg or diastolic BP between 80 to 89 mmHg; grade 2: ≥140/90 mmHg or under specific drug treatment [24, 25]. Dyslipidemia was defined as either one of the following: (1) total cholesterol ≥240 mg/dL; (2) high-density lipoprotein cholesterol <40 mg/dL; (3) triglyceride ≥200 mg/dL; (4) low-density lipoprotein cholesterol ≥160 mg/dL; or (5) under specific drug treatment [26–28]. Abnormal liver function test was defined as elevated serum alanine aminotransferase (ALT) above upper limit of normal range. For missing information in major socioeconomic risk factors for steatotic liver disease (marriage, education status and annual household income), they were coded as “unreported” and analyzed as a unique category. The proportion of missing information in the other variables (eg, smoking, diet, exercise…, etc) was less than 3% and was accounted for by imputing the most prevalent category.

Study population and inclusion/exclusion criteria

Records of all participants attending health checkups from 2010 to 2017 were retrieved (N = 382,914). We excluded 17,642 records (4.6%) from the analysis due to missing information of home address, living at regions without reliable air pollution exposure data (e.g., living in remote islands), or missing sonography data to determine hepatic steatosis status. For the participants with more than one health checkup record during the study period, the earliest record (i.e., index health checkup) was retained, leaving 195,625 participants to the study cohort.

Exclusion criteria also included: age <18 years (n = 2,711), missing key information to define MASLD (n = 18,087) and missing alcohol intake information (n = 12,107). We further excluded participants with self-reported past histories of liver cirrhosis or hepatocellular carcinoma (n = 267), with viral or other hepatitis (n = 22,369), with excessive alcohol intake (n = 3,728, see “Covariates” session) and those with significant sonography abnormalities (n = 1,447.) In the remaining participants, 3,317 categorized as cryptogenic steatotic liver disease were not analyzed in our study. Finally, 131,592 participants were included in the current analysis. Among them, 53,431 (40.6%) met the definition of MASLD. The flow chart of study participants’ enrollment is illustrated in Fig. 1.

Fig. 1

The flow diagram of enrolling study participants of this study.

HCC: hepatocellular carcinoma, MASLD: metabolic dysfunction-associated steatotic liver disease

Assessment of residential air pollution exposure

The monthly average concentration for each of the 6 major air pollutants including PM with aerodynamic diameter ≤2.5 µm (PM_2.5), PM with aerodynamic diameter ≤10 µm (PM₁₀), carbon monoxide (CO), ozone (O₃), sulfur dioxide (SO₂), and NO₂ from all 74 air quality monitoring stations in Taiwan were retrieved from the air pollution database supervised by Taiwan Environment Protection Administration (TWEPA). Information on the measurement instrument, method and sensitivity for air quality monitoring was available from the TWEPA website [29]. The spatial-temporal concentration distribution of each pollutant was estimated by land-use regression in combination with XGBoost machine learning algorithm (LUR-XGBoost) or Hybrid Kirging-LUR with XGboost algorithm with a 50-m × 50-m grid resolution. The performance of these algorithm demonstrated 10-fold cross-validation r² as 0.94 for PM_2.5, 0.84 for CO, 0.90 for NO₂ and 0.88 for O₃ with similarly excellent performance for external validity [30–33]. Exposure levels of various air pollutants were assessed based on the geocode of each participant’s home address. We calculated 1- to 4-year average exposure levels prior to the date of index health checkup for each study participant (Table S1) and chose 3-year average of air pollutant concentration to indicate long-term air pollutant exposures in this study [17]. To account for the observed trend of gradually decreasing air pollutant concentrations over time, we included the year of each participant’s enrollment in our multivariate models for adjustment.

Definition of metabolic dysfunction-associated steatotic liver disease (MASLD)

MASLD was defined as cases who had hepatic steatosis combined with at least one of the 5 cardiometabolic criteria [10]: 1) body mass index (BMI) ≧23 Kg/m², or waist circumference (WC) >94 cm in men or WC >80 cm in women; 2) fasting serum glucose ≥100 mg/dL or HbA1c ≥5.7% or history of type 2 diabetes or under drug treatment; 3) BP ≥130/85 mmHg or under specific antihypertensive drug treatment; 4) plasma triglycerides (TG) ≥150 mg/dL or under lipid lowering treatment; 5) plasma high density lipoprotein-cholesterol (HDL) ≤40 mg/dL for men or 50 mg/dL for women or under lipid lowering treatment. Hepatic steatosis was determined on abdominal ultrasonography performed by well-trained clinicians at MJ clinics [34].

Identification of potential non-linearity relationship of air pollution exposure and MASLD

To determine the non-linear link functions between air pollutant exposure and MASLD [17, 18], we combined data visualization with regression modeling to determine the best-fitted link function as the following:

1. 1)    We initially assessed the proportionality or linearity of air pollution exposure and MASLD odds using the Box-Tidwell test for each specific pollutant by testing for the significance of interaction between continuous air pollutant concentrations and their corresponding natural log [35].
2. 2)    If the proportional odds assumption held, a ‘linear link function’ was identified. For non-proportional relationships, we generated logit probability plots of MASLD by dividing air pollution concentrations into 100 percentile bins.
3. 3)    We used an iterative method minimizing non-linear least squares to identify inflection points and estimated the slopes of linear logistic regression models before and after these points, performed with the SAS procedure NLIN (SAS Institute, Cary, NC). Next, several non-linear link functions were explored and compared based on the Akaike Information Criterion (AIC). These functions included:
4. 3.1)    “Threshold link function”: assigned exposure below the inflection point as uniform baseline hazard (slope = 0) and a linear link function with proportional odds above the inflection point.
5. 3.2)    “Interaction link function”: two different linear link functions for exposure above and below the inflection point. An additional variable was then introduced, where a value of 0 and 1 was assigned to the participants with exposure below and above the inflection point, respectively. The original continuous exposure level, the newly created dummy variable, and their interaction term were then simultaneously included in the model.
6. 3.3)    “Restricted cubic spline (RCS) link function”: piecewise polynomials capturing curving relationships between exposure and logit probability of MASLD. A default RCS function with 4 knots was initially utilized with the location of knots determined according to previous literature [36]. If the 4-knot RCS model proved to be superior to other linear-based models, RCS functions with 3 to 7 knots were then sequentially evaluated to identify the optimal number of knots.

Any pair of the above link functions was considered equally fitted if a difference in AIC was 10 or less between link functions. In such circumstances, the principle of parsimony applied. The best-fitted link function was chosen as the following order: “linear link function”, “threshold link function”, “interaction link function” and “RCS link function”, with fewer knots was considered simpler.

To identify optimal link functions for each air pollutant and MASLD, we compared observed and predicted logit probabilities in both unadjusted and covariate-adjusted one-pollutant logistic regression models. Subsequently, the covariate-adjusted six-pollutant model was formed using these optimal link functions.

Statistical analysis

Descriptive statistics were presented mean ± standard deviation (SD) for continuous variables and number/percentage for categorical variables when comparing characteristics of individuals with and without MASLD. One-pollutant and six-pollutant logistic regression models, with best-fitted link functions, were presented with estimated covariate adjusted odds ratios (aORs) and 95% CIs. Pollutant exposure was analyzed both continuously and discretely, categorized by inflection points from best-fit link functions. For multiple logistic regression models, hypertension, diabetes mellitus and dyslipidemia were considered as part of the outcome variable within MASLD definition and were left unadjusted. Sensitivity analyses were conducted to assess the robustness of results: 1) The 6-pollutant model with linear link function for all 6 pollutants; 2) A 5-pollutant model that excluded PM₁₀, due to its relatively high correlation with PM_2.5 and SO₂. All analyses were performed using SAS version 9.4 (SAS Institute, Cary, NC), with a significance level of α = 0.05.

Results

In the 131,592 health checkup participants included, 53,431 (40.6%) had MASLD (Fig. 1). Socio-demographic characteristics and clinical parameters for participants with and without MASLD were compared and shown in Table 1. Participants with MASLD were older, predominantly male, and had higher percentage of being married. Furthermore, they exhibited a higher prevalence of current or former smoking and alcohol intake. Table 1 also showed the 3-year average exposure to 6 air pollutants prior to the index checkup. Compared to those without MASLD, participants with MASLD tended to have greater exposure levels of all 6 commonly monitored air pollutants.

Table 1 The demographic variables and clinical parameters of people with and without MASLD.

	No MASLD (n = 78,161)	MASLD (n = 53,431)
	Mean ± SD / n (%)	Mean ± SD / n (%)
Age, years	37.9 ± 11.9	44.3 ± 12.5
Sex
Male	27,937 (35.7)	35,259 (66.0)
Marriage
Single	34,003 (43.5)	14,818 (27.7)
Married	40,278 (51.5)	35,823 (67.1)
Unreported	3,880 (5.0)	2,790 (5.2)
Education
Illiterate	574 (0.7)	848 (1.6)
Elementary or middle schools	4,871 (6.2)	6,468 (12.1)
High schools	25,695 (32.9)	19,390 (36.3)
College or higher	45,795 (58.6)	25,720 (48.1)
Unreported	1,226 (1.6)	1,005 (1.9)
Annual household income
<28,000 USD	32,825 (42.0)	20,517 (38.4)
28,000 to 70,000 USD	30,204 (38.7)	22,093 (41.3)
≧70,000 USD	5,182 (6.6)	4,460 (8.4)
Unreported	9,950 (12.7)	6,361 (11.9)
Smoking
Never	64,087 (82.0)	37,668 (70.5)
Former smoker	3,705 (4.7)	4,716 (8.8)
Current smoker	10,369 (13.3)	11,047 (20.7)
Alcohol drinking
Never	68,135 (87.2)	44,101 (82.5)
Former	1,402 (1.8)	1,460 (2.7)
Current	8,624 (11.0)	7,870 (14.7)
Vegetable Intake
Seldom	13,581 (17.4)	9,087 (17.0)
Moderate	46,646 (59.7)	32,452 (60.7)
Frequent	17,934 (22.9)	11,892 (22.3)
Fruit intake
Seldom	26,169 (33.5)	16,531 (31.0)
Moderate	42,420 (54.3)	29,901 (56.0)
Frequent	9,572 (12.3)	6,999 (13.1)
Sugary drink intake
Seldom	27,053 (34.6)	19,414 (36.3)
Moderate	38,167 (48.8)	25,266 (47.3)
Frequent	12,941 (16.6)	8,751 (16.4)
Fried food intake
Seldom	23,240 (29.7)	15,313 (28.7)
Moderate	41,657 (53.3)	27,961 (52.3)
Frequent	13,264 (17.0)	10,157 (19.0)
Habit of regular exercise	25,329 (32.4)	17,805 (33.3)
Body mass index (BMI), kg/m2	21.3 ± 2.6	26.2 ± 3.4
Waist circumference (WC), cm
Male	77.7 ± 6.9	88.0 ± 8.0
Female	68.2 ± 5.9	79.7 ± 8.1
Diabetes Mellitus	1,090 (1.4)	5,268 (9.9)
Blood pressure (BP)
Systolic BP, mmHg	111.0 ± 15.4	123.5 ± 16.9
Diastolic BP, mmHg	69.5 ± 10.1	77.1 ± 11.2
Hypertension
Grade 1	10,363 (13.3)	13,193 (24.7)
Grade 2	5,672 (7.3)	13,936 (26.1)
Dyslipidemia	8,823 (11.3)	21,308 (39.9)
Triglyceride (TG), mg/dL	80.7 ± 44.2	148.9 ± 94.4
High-density lipoprotein-cholesterol (HDL), mg/dL
Male	57.1 ± 12.4	49.2 ± 9.7
Female	68.3 ± 14.7	57.3 ± 12.6
Low-density lipoprotein-cholesterol (LDL), mg/dL	106.8 ± 29.4[3(<0.1)]	124.1 ± 32.5[29(0.1)]
Cholesterol (Chol), mg/dL	187.0 ± 32.8	202.5 ± 36.1
Platelet count, 103/µL	247.4 ± 56.4[1(<0.1)]	254.0 ± 57.2[2(<0.1)]
Alanine aminotransferase (ALT), U/L	19.7 ± 13.5[1(<0.1)]	37.5 ± 28.6[3(<0.1)]
Abnormal liver function test	5,953 (7.7)	22,052 (41.3)
Glutamyl Transferase (GGT), U/L	19.4 ± 20.0[837(1.1)]	36.8 ± 41.1[438(0.8)]
Total bilirubin, mg/dL	1.0 ± 0.4[38(<0.1)]	0.9 ± 0.4[35(<0.1)]
Albumin, g/dL	4.5 ± 0.2[41(<0.1)]	4.5 ± 0.2[32(<0.1)]
3-year Average Exposure of Air Pollutant
PM2.5, µg/m3	30.3 ± 6.2	31.0 ± 7.0
PM10, µg/m3	51.1 ± 12.3	52.6 ± 13.9
CO, µg/m3	701.8 ± 234.6	707.3 ± 235.2
O3, µg/m3	52.7 ± 6.4	53.0 ± 6.6
SO2, µg/m3	10.6 ± 3.7	11.1 ± 4.3
NO2, µg/m3	37.3 ± 9.6	37.7 ± 9.7
Cardiometabolic criteria of MASLDa
Criteria 1	19,557 (25.0)	46,238 (86.5)
Criteria 2	21,175 (27.1)	33,790 (63.2)
Criteria 3	11,938 (15.3)	23,144 (43.3)
Criteria 4	5,385 (6.9)	21,600 (40.4)
Criteria 5	6,867 (8.8)	13,597 (25.5)

CO: carbon monoxide; MASLD: metabolic dysfunction-associated steatotic liver disease; NO₂: nitrogen dioxide; O₃: ozone; PM_2.5: particulate matter with an aerodynamic diameter ≤2.5 µm; PM₁₀: particulate matter with an aerodynamic diameter ≤10 µm; SO₂: sulfur dioxide; USD: United States dollars

^[] Numbers and percentage of cases with missing information for the covariate

^aMASLD criteria and prevalence according to:

Criteria 1: BMI ≥23 kg/m² (in Asians) or WC >94 cm for men or WC >80 cm for women

Criteria 2: fasting serum glucose ≥100 mg/dL or HbA1c ≥5.7% or type 2 diabetes or under treatment for type 2 diabetes

Criteria 3: BP ≥130/85 mmHg or under specific antihypertensive drug treatment

Criteria 4: plasma TG ≥150 mg/dL or under lipid lowering treatment

Criteria 5: plasma HDL ≤40 mg/dL for men or HDL ≤50 mg/dL for women or under lipid lowering treatment

Figure 2 illustrated the process of determining the best fitted link function between PM_2.5 and MASLD in a logistic regression model. In Fig. 2A, the observed logit probability showed a non-linear relationship with a potential inflection point and the Box-Tidwell test confirmed violation against proportional odds. To determine the inflection point, we fitted two linear link functions below and above a specific cutoff level, then estimated the proper slope of the two linear functions, which finally identified the best-fitted inflection point at 34.64 µg/m³. Subsequently, we constructed one-pollutant regression models (PM_2.5) with each of the 4 link functions: linear, threshold, interaction, and RCS. The model with interaction link function yielded the lowest AIC. However, after adjusting for multiple covariates in one-pollutant model, both interaction and threshold link functions resulted in very similar AIC values (Fig. 2B). According to our predefined rule of link function selection, the model with threshold link function was used to address the relationship between exposure to PM_2.5 and odds of MASLD (Fig. 2C). For the PM_2.5 exposure above 34.64 µg/m³, each 1 µg/m³ increase exhibited an aOR of 1.068 (95% CI: 1.064 to 1.072) for MASLD.

Fig. 2

Observed and predicted logit probability of 3-year average PM_2.5 exposure for MASLD.

(A) Scatter plot of the observed logit probability and comparisons of the predicted logit probability by crude models with various link functions: linear (proportional odds) link function, threshold link function, interaction link function, and restricted cubic spline (RCS) link function; (B) Comparisons of the predicted logit probability by covariate adjusted models with various link functions; and (C) The model with “threshold link function” was finally selected to analyze the covariate adjusted logit probability in association with 3-year average exposure of PM_2.5; and categorization of exposure was based on the threshold level and was used to calculate the covariate adjusted odds ratio and 95% confidence interval (CI) of MASLD. The models in Fig. 2B and 2C were adjusted for age, sex, marriage, education, household income, alcohol, smoking, fried food intake, vegetable intake, fruit intake, sugary drink intake, habit of regular exercise and the year of participant enrollment.

MASLD: metabolic dysfunction-associated steatotic liver disease; PM_2.5: particulate matter with an aerodynamic diameter ≤2.5 µm.

Similar processes were repeated to determine the best-fitted link function for the association of MASLD with PM₁₀, O₃, CO, SO₂, and NO₂, respectively; and the results are shown in Fig. S1–S5. The interaction link function was considered as the best-fitted link function for PM₁₀ with a cutoff concentration level of 57.93 µg/m³ (Fig. S1C) and for CO with cutoff of 643.6 µg/m³ (Fig. S3C). Since RCS function with 4 knots showed the better fit for the relationship between O₃ and MASLD, the RCS functions with 3–7 knots were further compared (Fig. S2B and S2C). As such, the RCS link function with 3 knots was finally considered as the best-fitted link function for O₃. Through the same process, an RCS link function with 3 knots (Fig. S4D) and 4 knots (Fig. S5D) was finally determined to be the best-fitted link function for SO₂ and NO₂, respectively.

Both Fig. 3 and Table 2 showed the results from covariate adjusted six-pollutant model. Based on the best-fitted link function, per 1 µg/m³ increase in PM_2.5 exhibited a significantly increased aOR of 1.036 (95% CI: 1.030 to 1.042) for MASLD when exposure level above 34.64 µg/m³. On the other hand, a 1 µg/m³ increase in PM₁₀ was found to be significantly associated with reduced and increased risk of MASLD at different exposure levels. For exposure level below 57.93 µg/m³, the aOR was 0.992 (95% CI: 0.990–0.993), while for exposure level above 57.93 µg/m³, the aOR was 1.019 (95% CI: 1.016–1.021). The model with interaction link function for CO exhibited an increase and plateau pattern. CO increase per 1 µg/m³ was significantly and positively associated with increased aOR of MASLD at CO <643.6 µg/m³.

Fig. 3

Predicted logit probability from six-pollutant model of 3-year average air pollutant exposure for MASLD.

(A) the “threshold link function” for particulate matter with an aerodynamic diameter ≤2.5 µm, PM_2.5; (B) the “interaction link function” for particular matter with an aerodynamic diameter ≤10 µm, PM₁₀; (C) the “restricted cubic spline (RCS) link function” with 3 knots for ozone, O₃; (D) the “interaction link function” for carbon monoxide, CO; (E) the “RCS link function” with 3 knots for sulfur dioxide, SO₂; and (F) the “RCS link function” with 4 knots for nitrogen dioxide, NO₂. The model was adjusted for age, sex, marriage, education, household income, alcohol, smoking, fried food intake, vegetable intake, fruit intake, sugary drink intake, habit of regular exercise and the year of participant enrollment. The plots were generated by adjusting covariates and other pollutants to their mean values. The cut-off values in Fig. 3A, 3B and 3D were the inflection values based on nonlinear regressions, while cut-off values in Fig. 3C, 3E and 3F were based on observation of restricted cubic spline plots.

MASLD: Metabolic dysfunction-associated steatotic liver disease.

Table 2 Covariate adjusted odds ratios of MASLD in relation to 3-year average exposure to various air pollutants estimated from covariate adjusted six-pollutant logistic regression model with various link functions

Pollutant	Exposure rangea	Model 1 Model based on the best-fitted link function for each pollutant AIC: 153,539	Model 2 Model with categorical exposure based on the cut-off points determined by Model 1 AIC: 154,115
		Adjusted OR (95% CI) per 1 µg/m3 increaseb	Adjusted OR (95% CI)bc
PM2.5		Threshold link function
	<34.64 µg/m3	1.000 (1.000–1.000)	Ref
	≥34.64 µg/m3	1.036 (1.030–1.042)	1.235 (1.190–1.281)
PM10		Interaction link function
	<57.93 µg/m3	0.992 (0.990–0.993)	Ref
	≥57.93 µg/m3	1.019 (1.016–1.021)	1.152 (1.110–1.196)
O3		3-knot RCS link function
	<56 µg/m3	0.996 (0.994–0.999)	Ref
	≥56 µg/m3	1.006 (1.002–1.011)	1.074 (1.043–1.107)
CO		Interaction link function
	<643.6 µg/m3	1.0006 (1.0004–1.0007)	Ref
	≥643.6 µg/m3	1.0001 (1.0000–1.0001)	1.127 (1.095–1.161)
SO2		3-knot RCS link function
	<10 µg/m3	1.001 (0.992–1.010)	Ref
	≥10 µg/m3	0.999 (0.994–1.004)	1.061 (1.031–1.091)
NO2		4-knot RCS link function
	<33 µg/m3	0.992 (0.989–0.996)	Ref
	33–48 µg/m3	1.008 (1.008–1.008)	1.043 (1.012–1.074)
	≥48 µg/m3	1.002 (0.998–1.006)	1.142 (1.094–1.192)

AIC: Akaike information criterion; CI: confidence interval; CO: carbon monoxide; MASLD: metabolic dysfunction-associated steatotic liver disease; NO₂: nitrogen dioxide; O₃: ozone; PM_2.5: particulate matter with an aerodynamic diameter ≤2.5 µm; PM₁₀: particulate matter with an aerodynamic diameter ≤10 µm; OR: odds ratio; RCS: restricted cubic spline; SO₂: sulfur dioxide

^a3-year average of exposure of pollutants

^bAdjusted for age, sex, marriage, education, household income, alcohol, smoking, fried food intake, vegetable intake, fruit intake, sugary drink intake, habit of regular exercise and the year of participant enrollment.

^cThe cut-off values for PM_2.5, PM₁₀ and CO were the inflection values based on nonlinear regressions, while cut-off values for O₃, SO₂ and NO₂ were based on observation of restricted cubic spline plots.

The logistic regression model with RCS link function for O₃ showed a tendency of slight fluctuation followed by an increase across the exposure range. One µg/m³ increase in O₃ was significantly associated with an increased aOR of MASLD at 1.006 (95% CI: 1.002–1.011) only when O₃ ≥56 µg/m³. SO₂ was the sole pollutant demonstrating no association with MASLD, as its aOR did not significantly differ from 0 across the entire exposure range. The link function for NO₂ showed a trend of decrease, rising and reaching a plateau in aOR of MASLD over the exposure range. For NO₂ exposure in the range of 33–48 µg/m³, per 1 µg/m³ increase was associated with a significantly elevated aOR of MASLD at 1.008 (95% CI: 1.008–1.008) (Table 2).

The aORs of MASLD estimated from the six-pollutant model in association with air pollutant exposure categorized by the inflection points based on the best-fitted link functions are also shown in the right column of Table 2. In this model, compared to participants with exposure levels below the inflection points, PM_2.5, PM₁₀, O₃, CO and NO₂ all exhibited significantly increased aOR above the inflection points, while SO₂ categorizations did not.

Table S2 displays aORs of MASLD in association with 1 µg/m³ increase in air pollutant estimated from logistic regression model with linear link functions assumed for all air pollutants. In this model, all pollutants, except for PM₁₀, demonstrated significant increase in aORs of MASLD in association with 1 µg/m³ increase in exposure. Comparatively, the best fitted model (Fig. 3 and left column of Table 2) exhibited better model fitting (AIC: 153,539) than the regression model with categorized exposure (right column of Table 2) (AIC: 154,115) and the regression model with linear link functions assumed (Table S2) (AIC: 154,074), respectively. Due to marginally high correlation between PM_2.5 and PM₁₀ (Table S3), A 5-pollutant model without PM₁₀ was built. As shown in Table S4 and Fig. S6, the associations between air pollutants and MASLD remained consistent with the main model, with the exception of SO₂, which displayed a diminished J-shaped relationship.

Discussions

Although the nature of cross-sectional design precludes the causal inference, the current study demonstrated a non-linear relationship between 3-year exposure and MASLD based on a large-scale sample of 131,592 health checkup participants. After adjustment for potential confounders, a “threshold link function”, “interaction link functions” and “restricted cubic spline (RCS) link functions” showed the best-fitted associations of MASLD with PM_2.5, PM₁₀/CO, and O₃ /SO₂/NO₂. The significantly positive associations between air pollution and MASLD were only found above (PM_2.5, PM₁₀, O₃), below (CO), and within (NO₂) certain levels of exposure in a six-pollutant model. While the results from two sensitivity analyses that excluded PM₁₀ from the model to avoid potential problem of co-linearity or assumed a linear link functions between all air pollutants and MASLD yielded the highly comparable results. Besides, the non-linear model also exhibited a superior ability to predict probability (odds) of MASLD (Fig. 3 and Table 2), compared to the models based on categorized exposure or assumed linear link functions.

In this study, 40.6% of the health checkup participants were identified to have MASLD. The prevalence of NAFLD in Taiwan in meta-analysis based on published data was estimated to be 33.29% (26.42%–40.96%) [37] and it might exceed 40% in health checkup participants [38]. In prior studies, VoPham et al. analyzed the US National Inpatient Sample using a spatiotemporal exposure model based on ZIP codes, revealing an OR of 1.24 (95% CI, 1.15–1.33) for NAFLD per 10 µg/m³ increase of PM_2.5 [19]. Guo et al. investigated the Chinese CMEC cohort, utilizing a satellite-based forest approach for exposure assessment, showing associations with 3-year average PM_2.5 and NO₂ exposure, resulting in ORs of 1.29 (1.25–1.34) and 1.15 (1.12–1.17) for MAFLD per 10 µg/m³ increment, respectively [17]. Sun et al. utilized the Taiwan MJ cohort and employed satellite-based atmospheric optical depth measurement for exposure assessment, reporting hazard ratios of 1.06 (1.04–1.09) and 1.05 (1.03–1.07) for NAFLD, when PM_2.5 increased every 1 µg/m³ above 23.5 µg/m³, dependent on the determinant for fatty liver status [18]. Li et al. utilized UK Biobank data and a land-use regression model for exposure assessment, and found hazard ratios of 1.10 (1.05–1.14), 1.14 (1.09–1.20) and 1.19 (1.13–1.24) for PM_2.5, PM₁₀ and NO₂, respectively, for each interquartile range increase [20]. Our current study using Taiwan MJ cohort and land-use regression in combination with machine learning algorithm for exposure assessment. Similar to most of the abovementioned studies, our data exhibited positive associations of 3-year elevated exposure to PM_2.5, PM₁₀ and NO₂ with MASLD. However, our study advanced the knowledge by demonstrating that the relationships between the above air pollutants and MASLD were not linear, which was implicitly suggested in the studies by Guo and Sun [17, 18], who nonetheless did not further elaborate such potential non-linear association between air pollution exposure and fatty liver disease.

A recent study demonstrated that long term exposure to ambient PMs and NO₂ is associated with both fatty liver disease and liver cirrhosis [39]. However, commonly monitored pollutants like SO₂ and CO on hepatic injury or fatty liver diseases has not been explored in previous research. It is worth noting that air pollutants rarely existed in a single-pollutant form [40], and the typical pollutants used for air quality assessment, as the 6 pollutants in our study, were intercorrelated (e.g., NO₂ is a precursor of O₃). Considering multipollutant exposures collectively is crucial for a comprehensive understanding of the impact of air pollution, especially from public health regulation perspectives [41]. The dose-response relationship between a specific air pollutant and health outcomes may change when accounting for other pollutants, as observed in our study. Simultaneous inclusion of the six air pollutants in our regression model adds further reassurance of validity for the associations between PM_2.5, PM₁₀, CO, O₃, NO₂ and MASLD. Expanding the MASLD prediction model to encompass a broader spectrum of pollutants could enhance its accuracy. However, regular monitoring of additional pollutants beyond the six pollutants would substantially elevate costs. While incorporating interaction terms could refine the model’s performance, the methodologies for statistical testing are still evolving, and the resulting models may become overly complex, leading to unstable coefficient estimates and challenging interpretation [40].

The mechanism underlying the association between air pollution and MASLD remains unclear. However, animal studies demonstrated that exposure to PM_2.5 increased hepatic dyslipidemia and increased oxidative stress, resulting in hepatic inflammation similar to non-alcoholic steatohepatitis and fibrosis [42, 43]. On the other hand, inhaled O₃ was linked to glycolysis dysfunction and glucose intolerance, which may be precursors of metabolic dysfunctions [44]. Chronic exposure to ambient air pollution (particularly PM_2.5) was associated with metabolic syndrome, systemic and hepatic inflammation in epidemiological studies [14, 45]. Nitrogen dioxide as part of traffic related air pollution was also found to be related to elevated cytokeratin-18 and may be related to liver injury [13]. Oxidative stress and systemic inflammation induced by ambient air pollution may be the major mechanism leading to the presentation of MASLD [46–48].

In the six-pollutant model, as anticipated, some health hazards observed in one-pollutant models were attenuated by the presence of other pollutants, resulting in lower ORs. Sulfur dioxide was the only air pollutant with non-positive association within the whole exposure range. In one pollutant model, SO₂ exhibited a J-shaped curve (Fig. S4D) but transformed into a null pattern in the six-pollutant model (Fig. 3E). Upon removing PM₁₀ from the model, the curve exhibited a more flattened J-shaped association. (Fig. S6D). These findings indicate that SO₂-MASLD association might be highly influenced by the representativeness of our sample and the SO₂’s association to other pollutants. Firstly, petrochemical parks and factories in Taiwan are predominantly located in rural, lower socioeconomic areas, whereas MJ clinics are in urban areas. This geographic disparity may restrict the inclusion of individuals from rural, lower economic areas near petrochemical parks, who are likely to have high SO₂ exposure, in our study. Given the positive correlation between lower socioeconomic status and fatty liver disease [49], enrolling participants solely from MJ clinics might introduce selection bias, underrepresenting those with higher SO₂ exposure and fatty liver disease. Secondly, since SO₂ emission primarily originates from industrial sources, its effects might be confounded by other pollutants like PM_2.5 or NO₂. This could explain the altered association observed in multipollutant models. Consequently, future research with more precise SO₂ measurement in regions heavily impacted by petrochemical pollution is essential for a clearer understanding of these associations.

Several strengths were involved in this study. Firstly, it revealed a non-linear correlation between air pollution and MASLD, enabling precise risk estimation in distinct exposure concentration ranges, rather than relying solely on single estimates. Secondly, our exposure assessment method offered superior 50-m × 50-m grid resolution for air pollutant exposure assessment, with enhanced accuracy compared to prior studies. Thirdly, we devised a process to assess the non-linear relationship between chronic air-pollution exposure and MASLD, which was considered essential for understanding air pollution’s health impacts [50, 51]. Moreover, our multipollutant model, integrating all six common pollutants, outperformed regression models with categorized or linear exposure assumptions. By combining the multipollutant non-linear model and exposure assessment approaches in this study, the relationship between air pollution and MASLD can be better illustrated. To the best of our knowledge, this is the first study of its kind in the literature. This improved methodology may advance the exploration of mechanisms linking air pollution and MASLD.

Despite the aforementioned strengths, several limitations warrant consideration. Firstly, air pollution exposure assessment based on residential addresses may introduce exposure misclassification. Although we combined land-use regression and machine learning for exposure assessment, there precision of these measures may still be subject to variability, potentially underestimation air pollution-MASLD associations, aligning with previous published studies [17–20]. Secondly, this cross-sectional study with retrospective exposure assessment lacks disease course information before checkup. The cross-sectional design and absence of a temporal component preclude the establishment of causality in this study. Thirdly, comorbidity ascertainment via self-report also carries a risk of misclassification. Incomplete adjustment for potential confounders, like DM and HTN, could lead to residual confounding and overestimation the impact of air pollution on MASLD. Furthermore, despite our rigorous attempts to adjust for confounding variables, unaccounted factors such as genetic predispositions or occupational exposures may persist in this study. Fourthly, the six-pollutant model included both PM_2.5 and PM₁₀, and their potential collinearity might affect this model. The naive Pearson’s correlation analysis showed a coefficient of r = 0.770 (0.768–0.772) between PM_2.5 and PM₁₀, and r = 0.677 (0.674–0.680) between PM₁₀ and SO₂ (Table S3), which is marginally acceptable. The tolerance was ≥0.27 (with PM₁₀ having the lowest tolerance), and the variance inflation factors (VIF) was ≤3.69 (with PM₁₀ having the largest VIF) for all 6 pollutants and covariates indicating that multicollinearity should not be a problematic issue. A sensitivity analysis with 5-pollutant model excluding PM₁₀ yielded a very similar results (Table S4 and Fig. S6), supporting the validity of our findings. Finally, our results, derived from a health-checkup cohort in Taiwan, warrant validation in diverse populations with varying environmental exposures, genetic predispositions, and lifestyle patterns.

Our study identified non-linear thresholds for air pollutant associations with MASLD, with some inflection points exceeding the World Health Organization (WHO) recommended annual air quality levels [52]. For example, 34.64 µg/m³ compared to WHO’s 5 µg/m³ for PM_2.5, or 33 µg/m³ compared to WHO’s 5 µg/m³ for NO₂. Notably, 23.7% of our participant cases fell within the effective exposure range for PM_2.5 and 48.1% for NO₂. These findings highlight the necessity of mitigating air pollution exposure as a complementary strategy to conventional diet, behavior, and metabolic risk factor management in preventing MASLD.

Conclusions

In conclusion, this large-scale cross-sectional study demonstrated non-linear associations between 3-year exposure to six common air pollutants and MASLD at time of health check-up, in which PM_2.5, PM₁₀, O₃, CO and NO₂ were positively associated with MASLD only at certain ranges of concentrations.

Abbreviations

AIC

Akaike information criterion

ALT

alanine aminotransferase

BMI

body mass index

Chol

cholesterol

CIs

confidence intervals

CO

carbon monoxide

DM

Diabetes mellitus

HTN

Hypertension

HDL

high-density lipoprotein-cholesterol

MASLD

metabolic dysfunction-associated fatty liver disease

NAFLD

Non-alcoholic fatty liver disease

NO₂

nitrogen dioxide

OR

odds ratio

O₃

ozone

PM_2.5

particulate matter with an aerodynamic diameter ≤2.5 µm

PM₁₀

particulate matter with an aerodynamic diameter ≤10 µm

RCS

Restricted cubic spline

SO₂

sulfur dioxide

TG

triglyceride

WHO

World Health Organization

Declarations

Ethics approval and patient consent

The study was approved by the Institutional Review Board of National Cheng Kung University Hospital, Tainan, Taiwan (Approval Number: B-ER-110-456) with a waiver of informed consent.

Consent for publication

Not applicable.

Data availability statement

The data utilized for this study are not publicly available. However, de-identified data may be available on request, subject to approval by the internal review board and under a formal data use agreement. Contact the corresponding author for more information.

Competing interests

The authors declare no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Funding

This study was supported by research grants of Tainan Hospital, Ministry of Health and Welfare, Tainan, Taiwan (TNHPA11103) and National Science and Technology Council, Taiwan (NSTC 112-2314-B-006-067-MY3). The funding bodies are not involved in study design, data process (collection, analysis, interpretation), nor in the preparation or submission of the report.

Authors’ contribution

WCC: conceptualization, data curation, formal analysis, funding acquisition, investigations, software, visualization, writing - original draft preparation; PYW: data curation, software, writing – reviewing and editing; CDW: resources, supervision, writing – reviewing and editing; PNC: supervision, validation, writing – reviewing and editing; PCL: conceptualization, data curation, supervision, validation, writing – reviewing and editing; CYL: conceptualization, methodology, project administration, formal analysis, funding acquisition, resources, supervision, validation, visualization, writing – original draft preparation.

Acknowledgements

The authors acknowledged the funding and administrative support of Tainan Hospital, Ministry of Health and Welfare, Tainan, Taiwan.

References

• 1.  Ando Y, Jou JH. Nonalcoholic Fatty Liver Disease and Recent Guideline Updates. Clin Liver Dis. 2021;17:23–8.
• 2.  Rinella ME, Neuschwander-Tetri BA, Siddiqui MS, Abdelmalek MF, Caldwell S, Barb D, et al. AASLD Practice Guidance on the clinical assessment and management of nonalcoholic fatty liver disease. Hepatology. 2023;77:1797–835.
• 3.  Eslam M, Newsome PN, Sarin SK, Anstee QM, Targher G, Romero-Gomez M, et al. A new definition for metabolic dysfunction-associated fatty liver disease: An international expert consensus statement. J Hepatol. 2020;73:202–9.
• 4.  Nguyen VH, Le MH, Cheung RC, Nguyen MH. Differential Clinical Characteristics and Mortality Outcomes in Persons With NAFLD and/or MAFLD. Clin Gastroenterol Hepatol. 2021;19:2172–81.
• 5.  Kim D, Konyn P, Sandhu KK, Dennis BB, Cheung AC, Ahmed A. Metabolic dysfunction-associated fatty liver disease is associated with increased all-cause mortality in the United States. J Hepatol. 2021;75:1284–91.
• 6.  Chalasani N, Younossi Z, Lavine JE, Diehl AM, Brunt EM, Cusi K, et al. The diagnosis and management of non-alcoholic fatty liver disease: Practice guideline by the American Association for the Study of Liver Diseases, American College of Gastroenterology, and the American Gastroenterological Association. Am J Gastroenterol. 2012;107:811–26.
• 7.  Golabi P, Paik JM, AlQahtani S, Younossi Y, Tuncer G, Younossi ZM. Burden of non-alcoholic fatty liver disease in Asia, the Middle East and North Africa: Data from Global Burden of Disease 2009–2019. J Hepatol. 2021;75:795–809.
• 8.  An SJ, Yang EJ, Oh S, Park KJ, Kim T, Hong Y, et al. The association between urinary bisphenol A levels and nonalcoholic fatty liver disease in Korean adults: Korean National Environmental Health Survey (KoNEHS) 2015–2017. Environ Health Prev Med. 2021;26:91.
• 9.  Paik JM, Golabi P, Younossi Y, Srishord M, Mishra A, Younossi ZM. The Growing Burden of Disability Related to Nonalcoholic Fatty Liver Disease: Data From the Global Burden of Disease 2007–2017. Hepatol Commun. 2020;4:1769–80.
• 10.  Rinella ME, Lazarus JV, Ratziu V, Francque SM, Sanyal AJ, Kanwal F, et al. A multi-society Delphi consensus statement on new fatty liver disease nomenclature. Hepatology. 2023;78:1966–86.
• 11.  Barouki R, Samson M, Blanc EB, Colombo M, Zucman-Rossi J, Lazaridis KN, et al. The exposome and liver disease - how environmental factors affect liver health. J Hepatol. 2023;79:492–505.
• 12.  Cheung AC, Walker DI, Juran BD, Miller GW, Lazaridis KN. Studying the Exposome to Understand the Environmental Determinants of Complex Liver Diseases. Hepatology. 2020;71:352–62.
• 13.  Hsieh S, Leaderer BP, Feldstein AE, Santoro N, McKay LA, Caprio S, et al. Traffic-related air pollution associations with cytokeratin-18, a marker of hepatocellular apoptosis, in an overweight and obese paediatric population. Pediatr Obes. 2018;13:342–7.
• 14.  Zhang Z, Chang LY, Lau AKH, Chan TC, Chuang YC, Chan J, et al. Satellite-based estimates of long-term exposure to fine particulate matter are associated with C-reactive protein in 30 034 Taiwanese adults. Int J Epidemiol. 2017;46:1126–36.
• 15.  Arciello M, Gori M, Maggio R, Barbaro B, Tarocchi M, Galli A, et al. Environmental Pollution: A Tangible Risk for NAFLD Pathogenesis. Int J Mol Sci. 2013;14:22052–66.
• 16.  Yuan Q, Zhang H. Identification of differentially expressed genes and pathways in BEAS-2B cells upon long-term exposure to particulate matter (PM2.5) from biomass combustion using bioinformatics analysis. Environ Health Prev Med. 2023;28:51.
• 17.  Guo B, Guo Y, Nima Q, Feng Y, Wang Z, Lu R, et al. Exposure to air pollution is associated with an increased risk of metabolic dysfunction-associated fatty liver disease. J Hepatol. 2022;76:518–25.
• 18.  Sun S, Yang Q, Zhou Q, Cao W, Yu S, Zhan S, et al. Long-term exposure to air pollution, habitual physical activity and risk of non-alcoholic fatty liver disease: A prospective cohort study. Ecotoxicol Environ Saf. 2022;235:113440.
• 19.  VoPham T, Kim NJ, Berry K, Mendoza JA, Kaufman JD, Ioannou GN. PM2.5 air pollution exposure and nonalcoholic fatty liver disease in the Nationwide Inpatient Sample. Environ Res. 2022;213:113611.
• 20.  Li F, Liao J, Zhu B, Li X, Cheng Z, Jin C, et al. Long-term exposure to air pollution and incident non-alcoholic fatty liver disease and cirrhosis: A cohort study. Liver Int. 2023;43:299–307.
• 21.  Wu X, Tsai SP, Tsao CK, Chiu ML, Tsai MK, Lu PJ, et al. Cohort Profile: The Taiwan MJ Cohort: half a million Chinese with repeated health surveillance data. Int J Epidemiol. 2017;46:1744–1744g.
• 22.  Wen CP, Wai JPM, Tsai MK, Yang YC, Cheng TYD, Lee MC, et al. Minimum amount of physical activity for reduced mortality and extended life expectancy: A prospective cohort study. Lancet. 2011;378:1244–53.
• 23.  ElSayed NA, Aleppo G, Aroda VR, Bannuru RR, Brown FM, Bruemmer D, et al. 2. Classification and Diagnosis of Diabetes: Standards of Care in Diabetes—2023. Diabetes Care. 2023;46:S19–40.
• 24.  Wang TD, Chiang CE, Chao TH, Cheng HM, Wu YW, Wu YJ, et al. 2022 Guidelines of the Taiwan Society of Cardiology and the Taiwan Hypertension Society for the Management of Hypertension. Acta Cardiol Sin. 2022;38:225–325.
• 25.  Whelton PK, Carey RM, Aronow WS, Casey DE, Collins KJ, Himmelfarb CD, et al. 2017 ACC/AHA/AAPA/ABC/ACPM/AGS/APhA/ASH/ASPC/NMA/PCNA Guideline for the Prevention, Detection, Evaluation, and Management of High Blood Pressure in Adults: A Report of the American College of Cardiology/American Heart Association Task Force on Clinical Practice Guidelines. Hypertension. 2018;71:E13–115.
• 26.  Lee H, Lee Y, Kim SU, Kim HC. Metabolic Dysfunction-Associated Fatty Liver Disease and Incident Cardiovascular Disease Risk: A Nationwide Cohort Study. Clin Gastroenterol Hepatol. 2021;19:2138–47.
• 27.  Huang PH, Lu YW, Tsai YL, Wu YW, Li HY, Chang HY, et al. 2022 Taiwan lipid guidelines for primary prevention. J Formos Med Assoc. 2022;121:2393–407.
• 28.  Rhee EJ, Kim HC, Kim JH, Lee EY, Kim BJ, Kim EM, et al. 2018 Guidelines for the Management of Dyslipidemia. Korean J Intern Med. 2019;34:723–71.
• 29.  TWEPA. Introduction of Central Monitoring Sensors. 2020. https://airtw.epa.gov.tw/CHT/EnvMonitoring/Central/Tools.aspx.
• 30.  Wong PY, Lee HY, Chen YC, Zeng YT, Chern YR, Chen NT, et al. Using a land use regression model with machine learning to estimate ground level PM2.5. Environ Pollut. 2021;277:116846.
• 31.  Wong PY, Hsu CY, Wu JY, Teo TA, Huang JW, Guo HR, et al. Incorporating land-use regression into machine learning algorithms in estimating the spatial-temporal variation of carbon monoxide in Taiwan. Environ Model Softw. 2021;139:104996.
• 32.  Babaan J, Hsu FT, Wong PY, Chen PC, Guo YL, Lung SCC, et al. A Geo-AI-based ensemble mixed spatial prediction model with fine spatial-temporal resolution for estimating daytime/nighttime/daily average ozone concentrations variations in Taiwan. J Hazard Mater. 2023;446:130749.
• 33.  Wong PY, Su HJ, Lee HY, Chen YC, Hsiao YP, Huang JW, et al. Using land-use machine learning models to estimate daily NO2 concentration variations in Taiwan. J Clean Prod. 2021;317:128411.
• 34.  Tsai YL, Chen SCC. Joint effect of changes in physical activity and weight on incident non-alcoholic fatty liver disease. J Epidemiol Community Health. 2021;75:1215–21.
• 35.  Box GEP, Tidwell PW. Transformation of the Independent Variables. Technometrics. 1962;4:531–50.
• 36.  Harrell FE. Regression Modeling Strategies. Cham: Springer International Publishing, Switzerland, 2015.
• 37.  Li J, Zou B, Yeo YH, Feng Y, Xie X, Lee DH, et al. Prevalence, incidence, and outcome of non-alcoholic fatty liver disease in Asia, 1999–2019: a systematic review and meta-analysis. Lancet Gastroenterol Hepatol. 2019;4:389–98.
• 38.  Huang JF, Tsai PC, Yeh ML, Huang CF, Huang CI, Hsieh MH, et al. Risk stratification of non-alcoholic fatty liver disease across body mass index in a community basis. J Formos Med Assoc. 2020;119:89–96.
• 39.  Orioli R, Solimini AG, Michelozzi P, Forastiere F, Davoli M, Cesaroni G. A cohort study on long-term exposure to air pollution and incidence of liver cirrhosis. Environ Epidemiol. 2020;4:e109.
• 40.  Billionnet C, Sherrill D, Annesi-Maesano I. Estimating the Health Effects of Exposure to Multi-Pollutant Mixture. Ann Epidemiol. 2012;22:126–41.
• 41.  Dominici F, Peng RD, Barr CD, Bell ML. Protecting human health from air pollution: Shifting from a single-pollutant to a multipollutant approach. Epidemiology. 2010;21:187–94.
• 42.  Zheng Z, Zhang X, Wang J, Dandekar A, Kim H, Qiu Y, et al. Exposure to fine airborne particulate matters induces hepatic fibrosis in murine models. J Hepatol. 2015;63:1397–404.
• 43.  Zheng Z, Xu X, Zhang X, Wang A, Zhang C, Hüttemann M, et al. Exposure to ambient particulate matter induces a NASH-like phenotype and impairs hepatic glucose metabolism in an animal model. J Hepatol. 2013;58:148–54.
• 44.  Miller DB, Karoly ED, Jones JC, Ward WO, Vallanat BD, Andrews DL, et al. Inhaled ozone (O3)-induces changes in serum metabolomic and liver transcriptomic profiles in rats. Toxicol Appl Pharmacol. 2015;286:65–79.
• 45.  Clementi E, Talusan A, Vaidyanathan S, Veerappan A, Mikhail M, Ostrofsky D, et al. Metabolic Syndrome and Air Pollution: A Narrative Review of Their Cardiopulmonary Effects. Toxics. 2019;7:6.
• 46.  Mendez R, Zheng Z, Fan Z, Rajagopalan S, Sun Q, Zhang K. Exposure to fine airborne particulate matter induces macrophage infiltration, unfolded protein response, and lipid deposition in white adipose tissue. Am J Transl Res. 2013;5:224–34.
• 47.  Xu X, Liu C, Xu Z, Tzan K, Zhong M, Wang A, et al. Long-term exposure to ambient fine particulate pollution induces insulin resistance and mitochondrial alteration in adipose tissue. Toxicol Sci. 2011;124:88–98.
• 48.  Niu Z, Liu F, Yu H, Wu S, Xiang H. Association between exposure to ambient air pollution and hospital admission, incidence, and mortality of stroke: an updated systematic review and meta-analysis of more than 23 million participants. Environ Health Prev Med. 2021;26:15.
• 49.  Orkin S, Brokamp C, Yodoshi T, Trout AT, Liu C, Meryum S, et al. Community Socioeconomic Deprivation and Nonalcoholic Fatty Liver Disease Severity. J Pediatr Gastroenterol Nutr. 2020;70:364–70.
• 50.  Dominici F, Daniels M, Zeger SL, Samet JM. Air Pollution and Mortality. J Am Stat Assoc. 2002;97:100–11.
• 51.  Leitenstorfer F, Tutz G. Generalized monotonic regression based on B-splines with an application to air pollution data. Biostatistics. 2007;8:654–73.
• 52.  WHO. What are the WHO Air quality guidelines? 2021. https://www.who.int/news-room/feature-stories/detail/what-are-the-who-air-quality-guidelines.