Associations between long-term exposure to fine particulate matter and osteoporotic fracture risks in South Korea

The prevalence of osteoporotic fracture is increasing globally due to rapid population growth and aging. Current evidence suggests adverse impacts of air pollution on bone mineral density loss and osteoporosis, but population-based evidence for the associations between fine particulate matter (particulate matter no larger than 2.5 μm in diameter [PM2.5]) and osteoporotic fracture is limited due to the small number of studies. This longitudinal study assessed the associations between long-term exposure to PM2.5 and osteoporotic fracture incidence in adults aged ⩾40 years, who enrolled in the National Health Insurance Service-National Sample Cohort data in 2002–2019 in South Korea. A time-varying moving window of past exposures of PM2.5 up to ten past years was estimated for participants’ residential addresses using modeled PM2.5. We used Cox proportional hazard models to estimate hazard ratios (HRs) of time-variant moving concentrations of PM2.5 exposure and osteoporotic fracture. The Cox models calculated HRs for an interquartile range (IQR) increase in PM2.5 exposure, adjusting for age, sex, body mass index, health behaviors, medications, disease history, income, and urbanicity. We assessed 161 831 participants over 993 104 person-year of follow-up. Results suggested linear and positive exposure-response associations for past PM2.5 exposure in the prior four years or more. The IQR increase in 5-year moving average PM2.5 was significantly associated with increased osteoporotic risk (HR = 1.079, 95% CI: 1.001, 1.164). The HRs were significant in women (1.102, 95% CI: 1.011, 1.200) and the subset of women aged 50–74 years (1.105, 95% CI: 1.005, 1.214) but not in men overall or by age groups. The association was not significantly different by income, physical activities, urbanicity, or diet. Overall, long-term PM2.5 exposure was associated with increased osteoporotic fracture risks in Korean adults, especially women.


Introduction
The world is rapidly aging and 17% of the world's population is estimated to be over 65 years of age by 2050 (UN 2020). The population of South Korea is rapidly aging; the proportion of persons aged 65 or older was about 16% in 2020 and is expected to increase to 44% by the 2100s (Roser and Rodés-Guirao 2013). The majority of South Korean regions will become a super-aged society, defined as a society with 21% or more of the population age 65 or older, by 2029 . As a result, the elderly population's public health burden and healthcare expenditures are anticipated to increase. Common major health concerns of aging include bone-related health, such as osteoporotic fracture. Osteoporotic fracture patients use two times 2. Method

Study participants
We obtained the National Health Insurance Service-National Sample Cohort (NHIS-NSC) data, which is an individual-level nationwide open cohort dataset. About 97% of the entire population in South Korea is enrolled in the national insurance system (Sung et al 2020) and the NHIS-NSC data provide information of 2% of enrollees samples from the total national insurance records. The participants in the data were sampled based on age, sex, region, insurance type (e.g. employment insurance, self-employed insurance, dependents), and income. The variables include age, sex, residential information (i.e. district), insurance type, income, health behaviors, diagnosis, and medications compiled from medical records, treatments, and patient surveys. The spatial unit for residential addresses in the database is district (roughly analogous to the US city). The average size of the study districts was 396.6 km 2 (range 284.2-1816.2 km 2 ). As the residential address at the district level at the time of enrollment is available for each year and each participant in the dataset, a time-varying variable for spatiotemporal changes in air pollution exposure can be generated, accounting for participants' moves. The national insurance service in South Korea offers health examinations for insured persons through employment and self-insured head of households every one or two years based on the type of occupation. Dependents of employees or self-employed household heads can take the health examinations from age 40 years old or older. The NHIS-NSC provided the health examination data including self-reported disease history, smoking status, alcohol consumption, frequency of physical activities, height, weight, and body mass index (BMI). Alcohol consumption was measured as self-reported total daily alcohol intake (based on one unit of ethanol equivalent to a glass of 'Soju'). Given that the recommended amount of alcohol intake by the World Health Organization (WHO) is <5 units per day for men and <2.5 units per day for women , we defined high alcohol consumptions as ⩾5 units per day for men and ⩾3 units for women.
We identified the initial occurrence of osteoporotic bone fracture using the procedure code (e.g. cast) and  (Alver et al 2010, Mazzucchelli et al 2018, Sung et al 2020, Chiu et al 2021, Heo et al 2022. These previous studies did not include age-related osteoporosis (ICD-9 733.0) and pathologic fracture (ICD-9 733.1) in osteoporotic fractures. Our analyses excluded the ICD-10 codes for age-related osteoporosis with current pathological fracture (M80), as this code was newly introduced in October 2015 in the ICD-10 guideline and cannot be precisely compatible with the ICD-9 codes for osteoporotic fractures. Both hospitalizations and out-patient care for the initial occurrence of osteoporotic fracture were included in our analyses.
This study targeted persons ⩾40 years in the first enrollment year in the dataset. We defined the baseline period, where the past PM 2.5 exposure is measured, as the first ten years from the entry in the NHIS-NSC data. The number of participants in the NHIS-NSC data during 2002-2019 was 1137 861. Among these, the number of persons aged ⩾40 years in the first enrollment year was 353 563. Persons who were lost due to death or loss of insurance status during the baseline period were excluded from the analysis (n = 89 262). Then, persons who had a diagnosis of osteoporotic fracture during the baseline period were excluded (n = 17 157). Finally, persons who did not have health examination data for any year during the follow-up were excluded (n = 85 313). As a result, the number of participants included in the analysis was 161 831 (figure S1).

Machine-learning based air pollution modeling data
We considered PM 2.5 as the exposure of interest. We used the AiMS-CREATE prediction data for these air pollutants. The developed ensemble models constructed three machine learning models (random forest, light gradient boosting, and neural network) to generate monthly averages at a 1 km × 1 km resolution for South Korea for the years 2002-2020.
Due to the limited-service period of PM 2.5 monitoring stations in South Korea (2015-2020), the ensemble machine learning model was first constructed for the period from 2015 to 2020 with a total of 112 predictors including (1) atmospheric, meteorologic, and land-use variables collected from satellite-based data provided by the Google Earth Engine; (2) socioeconomic variables from the regional health determinant databased distributed by the Korea Disease Control and Prevention Agency; and (3) observed concentration of air pollutants (PM 10 , NO 2 , O 3 , SO 2 , and CO) from the monitoring station (2002-2020). More detailed information on the predictors and data collection is reported in the Supplementary Materials.
The researchers who developed the AiMS-CREATE dataset compared the modeled annual average 1 km 2 spatial resolution PM 2.5 to the monitoring data for 2015-2020, using cross-validated R 2 and root mean-square error (RMSE). During the entire modeling period (2015-2020), the R 2 value of the ensemble model of PM 2.5 was 0.871 with RMSE (µg m −3 ) of 2.994, and the annual R 2 ranged from 0.66 to 0.91, and annual RMSE ranged from 2.389 to 3.726. The R 2 for the monthly averages of PM 2.5 (January-December) in 2015-2020 ranged from 0.60 (June) to 0.86 (March). The R 2 for PM 2.5 were similar between urban regions (R 2 = 0.87) and rural regions (R 2 = 0.86).
The aforementioned machine learning models were fitted and optimized based on cross-validations. Then, the optimized models were used to predict the monthly PM 2.5 concentrations during 2002-2014.

Exposure assessments
We obtained monthly concentration estimates of ambient levels of PM 2.5 at 1 × 1 km 2 spatial resolution across contiguous South Korea from the AiMS-CREATE team at Pusan National University and Korea University. We used these data to calculate the annual mean concentration of each air pollutant for each district by averaging the predictions at the grid cells with centroid points inside the boundary of that district. For each calendar year, we assigned the annual average concentrations (2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018)(2019) for each district to NHIS-NSC enrollees according to the district of residence and calendar year to calculate long-term exposures to PM 2.5 . The aggregation of exposure was the spatial resolution of the district as this was the finest spatial resolution information for the participant's residential addresses from the health data.

Exposure-response associations
The unit of statistical analysis was person-years of follow-up. The baseline period was defined as the first ten years from the entry. The follow-up start year for the statistical models was defined as the 11th year from the entry to the cohort and long-term past exposure to air pollution up to ten past years was estimated for the first ten years from the entry to the data. To reflect changing levels of air pollution over time, we used a time-varying variable of moving average concentrations of PM 2.5 and O 3 from lag-1 yr to lag-10 yr. For example, for a participant who was enrolled in 2002, the statistical analysis for examining the risks starts from 2012, and the 1 yr moving average PM 2.5 concentration for examining the osteoporotic fracture risk in 2012 is the past 1 yr exposure in 2011. The 10 yr moving average PM 2.5 concentration for the risk in 2012 is the PM 2.5 exposure levels in 2002-2011. Similar approaches for time-variant exposure levels were applied in previous studies (Heo et al 2022). We used the 6 yr exposure window for the main results and sub-group analysis as positive exposure-response relationships were found for exposure in the past 6 years or more (i.e. 6-, 7-, 8-, 9-, and 10 yr moving averages) and the risks were highest for the 6 yr moving averages. We targeted adults aged ⩾40 years at the cohort entry year, such that the age of participants from the analysis start year was ⩾50 years.
Multivariate Cox proportional hazards models stratified by sex and age (5 yr interval) were separately used for different exposure windows to calculate hazard ratios (HRs) of bone fracture and 95% confidence intervals (CIs) in relation to long-term exposure to air pollution for an interquartile range (IQR) increase. The IQR for PM 2.5 was calculated using the annual PM 2.5 concentrations of all districts in South Korea during 2002-2019. This IQR value was applied to estimating HRs of PM 2.5 concentrations of all examined exposure windows. For the participants included in the Cox models, a right censoring occurred for those known to be deceased or lost enrollment for the NHIS-NSC. We had three different Cox models a priori, with an increasing level of adjustment for potential confounders. Model 1 was stratified by sex and age and included only a PM 2.5 exposure variable. Utilizing the rich medical information from the NHIS-NSC data, model 2 used the following potential confounders: exposure to oral glucocorticoids in the past year (yes, no), diagnosis of rheumatoid arthritis (yes, no), diagnosis of secondary causes of osteoporosis (yes, no), income-based insurance fee (<40th percentile, 40th-79th percentiles, ⩾80th percentile; higher percentile indicating higher income), smoking status (current, former, and never), high alcohol consumption (yes, no), BMI, Charlson comorbidity index, frequency of exercise per week (0, 1-2, 3-4, 5-6, 7 times), and study provinces (finest spatial resolution of the data; roughly analogous to US counties). All variables were included as time-varying variables.
In estimating HRs, the Cox models applied linear exposure-response associations. To examine the shape of exposure-response associations, natural spline functions with 4 degrees of freedom (df) were applied for the Cox model. The degree of freedoms of 4 were selected from the Cox model with the lowest Akaike's Information Criterion (AIC) (Akaike 1969).
Cox models can use calendar year as the time scale in an attempt to draw exposure information from the same calendar period within a risk set. Nonetheless, we used the survival time scale (e.g. survival years Y i for i = 1, …, N) as the time scale in our Cox regression models, which was suggested to be less subject to potential biases in assessing exposure-response associations for temporally variable environmental factors (e.g. PM), compared to other time scales (e.g. attained age or calendar time) (Griffin et al 2012). Also, this approach was chosen since the entry year did not vary substantially among the participants in the dataset.
We examined potential effect modification by individual-and regional-level factors. The examined individual-level factors included age (50-69 y, 70+ y), sex (men, women), sex and age (men 50-69 y, men 70+ y, women 50-69 y, women 70+ y), frequency of exercise (low: never, moderate: 1-4 times per week, high: 5-7 times per week), and income (low: 0-39th, moderate: 40-79th, high: 80-100th percentile). Effect modification by urbanicity was examined by defining metropolitan cities as urban regions and other study regions as sub-urban or rural regions. Also, we investigated potential effect modification by the average percentage of persons using vitamin supplements for the study provinces (e.g. <20%, ⩾20%) across 2007-2009, obtained from the Korea National Health and Nutrition Examination Survey by the Korea Disease Control and Prevention Agency. Effect modification was examined by adding a multiplicative interaction term between air pollution and each effect modifier. The significance of the interaction term was examined by the p-value at a significance level of 0.1.
All statistical analyses for associations were estimated using packages 'survival' and 'splines' in R Studio version 3.3.3 (RStudio 2021).

Sensitivity analysis
Studies apply multi-pollutant models including two or more air pollution variables to explore potential synergistic or additive effects of simultaneous exposures (Fan et al 2020). We conducted two-pollutants models including O 3 , NO 2 or CO in addition to PM 2.5 . We also applied a Cox model including a categorical variable instead of a continuous variable of moving average annual average PM 2.5 concentrations. The same covariates in the aforementioned model 2 were included in this model.

Characteristics of study participants
The average PM 2.5 concentrations in the study districts aggregated for 2002-2019 ranged from 23.0 to 34.2 µg m −3 (figure S2). The correlation coefficients of PM 2.5 (2002-2019) with PM 10 , O 3 , NO 2 , and CO were 0.83, −0.26, 0.33, and 0.64, respectively (table S1). Characteristics of the participants with and without osteoporotic fracture in the entry year or during the follow-up years are shown in table 1. We followed 161 831 participants over 993 104 person-years of follow-up. The average follow-up time was 6.1 (SD = 3) years for all participants. Of the total participants (n = 161 831), 8268 participants (5.1%) had osteoporotic fractures during the follow-up. The incidence of osteoporotic fracture was higher in women (63.7%) than men (36.3%). All examined individual-level covariates in table 1, except for the percentiles of household income, had significantly different distributions between the group experiencing outcome (i.e. osteoporotic fracture) and the group not experiencing the outcome (e.g. outcomes were more common in women than men). The PM 2.5 exposure concentrations for the entire follow-up years for each participant were significantly higher in the group experiencing outcome than in the group not experiencing the outcome for all examined exposure windows. The minimum 1 year moving average PM 2.5 concentration was 17.9 µg m −3 , which was higher than the national annual standard PM 2.5 (15 µg m −3 ). The average PM 2.5 concentrations for all examined exposure windows were higher than the national annual standard.

Associations between PM 2.5 exposure and osteoporotic fracture incidence
The estimated HRs from our linear Cox models are shown in table 2. The crude HRs were statistically significant for all examined moving average PM 2.5 concentrations; the HR for an IQR increase in PM 2.5 in the prior 6 years (i.e. lag1-lag6) was highest among the crude HRs (HR = 1.103, 95% confidence interval (CI): 1.042, 1.168). The adjusted HRs for an IQR increase in PM 2.5 in model 2 showed significant risks from the 4 yr moving average to the 10 yr moving average PM 2.5 . The highest and most significant HR was found for the 5 yr moving average PM 2.5 concentrations: HR = 1.079 (95% CI: 1.001, 1.164). Model 3, the two-pollutant models of PM 2.5 and O 3 , showed marginal significance for HRs for the prior 6-10 years of PM 2.5 concentrations. The inclusion of O 3 in the Cox models as a potential confounder slightly decreased the estimated HRs from model 2, but the results were robust.
The results were robust in the two-pollutant model, adjusting for O 3 exposure, but slightly decreased and showed marginal significant HRs (table S2). Table S3 shows the risk of osteoporotic fracture among the quantile groups of PM 2.5 concentrations. Relative to the lowest quartile PM 2.5 concentrations (i.e. Q1, reference), the crude HR was 6.9% higher in the second lowest quartile (Q2). The crude HR increased by 14.1% (95% CI: 1.059, 1.229) in the second highest quartile (Q3) and 15.7% (95% CI: 1.052, 1.274) in the highest quartile, compared to the lowest quartile. The adjusted HRs model 2 showed significantly higher HRs   in the Q3 (HR = 1.099, 95% CI: 1.006, 1.201) and Q4 groups (HR = 1.123, 95% CI: 1.002, 1.260) than the Q1 group. Figure 1 shows the exposure-response relationships between PM 2.5 and osteoporotic fracture. The relationship curves showed nearly linear and positive relationships for exposures of 6 years. Positive relationships were also found for exposures of 7, 8, 9, and 10 years (figure S3). For past exposure between one year and five years (i.e. lag1 to lag5), an inverted U-shaped curve was found, indicating potential non-linear associations for a shorter exposure period (figure S3). The non-linear curves of log HRs of osteoporotic fracture incidence and 6 year moving average annual PM 2.5 concentrations stratified by age and sex showed that the relationship was positive in women aged ⩾50 years, whereas an inverted U-shaped curve was found for men aged ⩾50 years ( figure S4). Potentially non-linear positive relationship curves were found for the groups of men aged ⩾75 years, women aged 50-74 years, and women aged ⩾75 years. Table 3. Effect modification by sex, age, exercise, household income, urbanicity, and regional vitamin intake rate for the hazard ratios (HRs) of osteoporotic fracture associated with an IQR increase (IQR = 4.9 µg m −3 ) in the 6-year moving annual average concentration of PM2.5 (n = 161 831 Note: a Seven metropolitan cities as urban regions and the rest regions as suburban or rural regions. b The cutoff value of 20% was the median percentage of adults age 50⩾ years using vitamin supplements among the study provinces. c P-values for the multiplicative interaction terms between PM2.5 and each effect modifier in the table. * : significant at a significance level of 0.05. †: significant at a significance level of 0.1. Table 3 shows HRs with an IQR increase in the 6 yr moving average PM 2.5 concentrations by age, sex, frequency of exercise, household income, urbanicity, and regional rates of vitamin supplement intake. Model 2 showed that the HR was significant in women (HR = 1.102, 95% CI: 1.011, 1.200), but not in men (HR = 1.042, 95% CI: 0.941, 1.154). In the models further stratified by age, the HR was significant only for women aged 50-74 years (HR = 1.105, 95% CI: 1.005, 1.214). Nonetheless, the interaction term for sex indicated no difference in the risks by sex. The HR (1.142, 95% CI: 1.027, 1.270) was significant for the moderate-income group (i.e. 40-79th percentiles of income), but the interaction term for income was not significant. Overall, there was no significant effect modification by the examined factors. Table 4 shows the characteristics and findings of the previous studies and the current study of osteoporotic fracture and air pollution. Among the listed studies, six studies including the current study explored the long-term effects of exposure using cohort or case-control study designs (

Discussion
The global incidence of bone fracture substantially increased from 1990 to 2019 as a result of aging and population growth (Wu et al 2021). The disease burden of osteoporotic fracture is expected to increase in South Korea along with the estimated increase in the number of fractures , warranting public health policies and strategies to reduce future burdens from risk factors such as air pollution. This nationwide longitudinal cohort study examined the associations between long-term PM 2.5 exposure and osteoporotic fracture incidence in Korean adults aged ⩾40 years. We found significant associations for the PM 2.5 exposure in the prior four to ten years before the osteoporotic fracture occurrence. Similarly, a previous study in the US suggested that one IQR (4.18 µg m −3 ) increase in PM 2.5 was associated with about a 4.1% increase in hospital admission for bone fracture (Prada et al 2017).
A few studies have suggested associations between PM 2.5 exposure and bone-related health outcomes including osteoporosis and osteoporotic fracture, but the results are inconclusive (Alver et al 2010, Mazzucchelli et al 2018, Prada et al 2020, Sung et al 2020, Chiu et al 2021, Mousavibaygei et al 2023. A recent review article about air pollution, bone mineral density loss, osteoporosis, and osteoporotic fracture suggested evidence for the adverse impact of PM on low bone mineral density, but the evidence was weak for osteoporotic fracture due to the small number of studies (Mousavibaygei et al 2023). For example, a study in Norway found a significantly negative association between PM 2.5 and bone mineral density loss in men aged 75-76 years, but there was no significant association between PM 2.5 and self-reported forearm fractures in the same group. Alver et al (Alver et al 2010) and (Mazzucchelli et al 2018) examined the risks using both PM 2.5 and PM 10 data, but some other studies examined the impacts of PM 2.5 only on osteoporotic fractures. The significant differences in the study designs, statistical models, target populations, chemical composition of particulate matter, and types of osteoporotic fracture among previous studies of osteoporotic fracture would contribute to the heterogeneity of the results for the associations (table 4). More future studies are warranted to better understand the impacts of PM exposure on osteoporotic fracture, including for different populations and forms of PM, such as from different sources.
Study results for PM 2.5 exposure and bone-related health outcomes (e.g. bone mineral density, osteoporosis, osteoporotic fracture) are scarce in South Korea due to the lack of PM 2.5 data. A previous Korean study of women aged ⩾50 years in three metropolitan cities found significant osteoporotic fracture risks for PM 2.5 but not for PM 10 (Sung et al 2020). Our previous study of men and women in seven metropolitan cities in South Korea examined the impacts of PM on osteoporotic fractures using PM 10 monitoring data, but PM 2.5 monitoring data were not available for our study period (Heo et al 2022). Given that PM 2.5 can be more harmful than PM 10 due to the smaller size and deeper penetrations into the lungs and organs, examining varying sizes of PM in relation to osteoporotic fracture risks is important in establishing the regulation standards for PM. However, direct comparison of results of our study and the previous Korean studies of PM 10 is challenging as earlier studies used PM 10 monitoring data and the current study used PM 2.5 modeling data. More research is needed to better understand the reasons for the risk differences between PM 10 and PM 2.5 (e.g. size, emission sources, chemical components, exposure levels). Also, further analysis is needed to compare the associations between PM and osteoporotic fractures using both PM 2.5 modeling and monitoring data as the monitoring data become more available.
Our study also examined the risk differences by various individual-and regional-level factors. In our analyses stratified by sex, the PM-associated osteoporotic fracture risk was significant for women but not for men. Also, we found significant and positive associations between PM 2.5 exposure and osteoporotic fractures in women aged 50-74 years but not in other age and sex groups. These results are similar to a previous study in Taiwan, which found that a 10 µg m −3 increase in long-term PM 2.5 exposure was associated with osteoporotic fractures in women but not for men (Chiu et al 2021). However, the findings for the differences in the osteoporotic fracture risks associated with air pollution by sex are not consistent. Although focusing on the short-term effects of PM 2.5 , a time-series study in Spain did not find significant associations between PM 2.5 and hip fractures in both male and female groups (Mazzucchelli et al 2018). Although elderly women are generally more vulnerable to osteoporosis, male patients also account for a large percentage of osteoporotic fractures (Colón-Emeric and Saag 2006). Future research is needed to further examine and compare the adverse impacts of PM on osteoporotic fracture between different sex groups.
Our study did not find significant effect modifications by household income, urbanicity, or regional vitamin intake rates. Only a few studies have examined potential effect modifications by socioeconomic status (e.g. income), exercise, and diet for the associations between PM and osteoporotic fracture (Prada et al 2020, Heo et al 2022. A previous US study found that the highest income group had significantly lower associations between PM 2.5 and osteoporotic fracture than lower income groups (Prada et al 2020). Our results showed that the associations between PM 2.5 and osteoporotic fracture incidence were significant in the moderate-income group, while the risks were not significant in low-or high-income groups. Lower mobility in lower-income persons and better health conditions and access to healthcare services in higher-income persons may contribute to this result, but further analysis is required. Populations in urban and rural regions have different health conditions and pollution exposure, but our results did not find significant effect modifications by urbanicity. Different definitions of the degree of urbanization (e.g. land use, population density) rather than the units of governmental administrations (e.g. metropolitan cities vs. non-metropolitan cities) would benefit future studies. In summary, the lack of research and heterogeneity of study results for effect modifications warrants further examinations of effect modifications.
The inconsistent exposure windows applied among previous studies indicate that there is currently no consensus for the exposure windows of long-term air pollution exposure and osteoporotic fracture risks. For instance, a study in Taiwan (Chiu et al 2021) considered the 1 yr moving average PM 2.5 concentrations for the one year before the incidence. A case-control study of self-reported forearm fracture in Norway considered PM 2.5 exposure for a fixed period from 1992 to 2001 (i.e. exposure in 10 years) for the participants recruited in 2000-2001(Alver et al 2010. A case-control study of hip fracture for Korean adults aged ⩾30 years averaged the concentrations of PM 10 during a fixed six-year period (e.g. 2010-2016) for the residential districts and assigned the averages to each participant (Oh and Song 2020). In our previous research where we examined exposure windows using the 3 yr moving average PM 10 concentration, the exposure-response relationship curves for osteoporotic fracture indicated an inverted U-shaped curve and the estimated HRs of PM 10 were not significant (Heo et al 2022).
In the current study, we examined a maximum exposure window of ten years past exposure, and the exposure-response relationship curves suggested linear positive associations. In this type of study, a baseline exposure period is often set from the time of initial enrollment in the cohort and the risk of outcome is examined after this select baseline period, which is determined as the follow-up period. The length of follow-up periods for each participant in the analysis can be substantially shortened when the baseline period is set to be long, while the length of follow-up in the analysis may be important to observe the outcome of interest (e.g. osteoporotic fracture). We attempted to use a baseline period that would be relevant for the mechanisms in which air pollution can interrupt the bone remodeling cycles. Cortical bone can take from 100 d and up to three years to completely replace its components dependent on health conditions (e.g. hormonal disorder) (Eriksen 2010). Therefore, we attempted to consider past exposure for three years or longer.
Our exposure-response relationship curves and the estimated HRs showed that the osteoporotic fracture risks associated with PM 2.5 exposure can be significant if the number of years applied for exposure assessment (i.e. the length of exposure windows) increased from five years to ten years. Although it is unclear due to the lack of information on the latent period for PM exposure and bone damage or loss, different mechanisms between PM and bone homeostasis may contribute to different exposure windows. Metals in PM may reduce the normal activation of vitamin D, mineral absorption from the intestines, and bone mineralization (Youness et al 2012), and such effects may slowly take place over a long period. Reduced ultraviolet B photons reaching the ground and reduced vitamin D synthesis due to outdoor air pollution (Prada et al 2020) may slowly weaken bone over the years. There is a growing literature showing decreasing trends for air pollution levels and the associations between air pollution exposure and health outcomes in South Korea (Jiang et al 2020), but evidence is limited about if the decreasing air pollution trends are relevant for the significant risks of osteoporotic fractures in relation to longer exposure windows. To our best knowledge, there is no observational research on air pollution and osteoporotic fracture considering the chemical composition and constituents of PM. Future research is highly required for studying the impact of different sources (e.g. vehicles, industry, household) and chemical compositions in addition to the air pollution trends. Future research is needed for examining a longer exposure period (e.g. >decade) and follow-up period as well as study of whether short-term PM 2.5 exposure (e.g. monthly) is be associated with increases in osteoporotic fracture risks.
Our study has several strengths. To our best knowledge, this is the first longitudinal study to examine the PM 2.5 effects on osteoporotic fracture risks across mainland South Korea. Our analysis accounted for residential moves across each study year and thereby more accurately estimated long-term PM 2.5 exposure for participants than some other long-term air pollution studies. Our statistical models adjusted for urbanicity in examining the associations between PM 2.5 and osteoporotic fracture, given that the health conditions, lifestyle, exposure to pollutants, characteristics of particles, and other risk factors for bone health would differ between urban and rural areas. Using the rich information of the national cohort data, our models adjusted for time-variant variables of secondary causes of osteoporosis and the use of glucocorticoids in addition to various risk factors for osteoporosis or low bone mineral density (Chung and Chan 2022). Both ICD-10 disease codes and procedure codes were used to more adequately identify cases of osteoporotic fractures. Various effect modifiers were examined in our models which can provide insights and evidence for risk differences among different groups. Also, we examined various lag periods for long-term past PM 2.5 exposure with a maximum lag period of ten years in our models. We considered a two-pollutant model using the PM 2.5 and O 3 data, and results suggested slightly increased impacts of PM 2.5 when O 3 was included in the models. As a previous review article noted (Mousavibaygei et al 2023), a multi-pollutant approach is important in understanding synergic or additive associations.
This study also has limitations. The spatial resolution of the final models (i.e. district) for exposure-response associations, which is based on the resolution of the available health data, does not account for spatial heterogeneity in pollution at smaller scales. Applying finer spatial resolutions such as postal code area or buffers of geocoded residence may lead to more accurate exposure estimates, and variation in measurement errors between the monitoring data and modeling data warranting future research. Second, we did not have measures of the bone mineral density of the study participants as the national cohort data did not collect such information. However, using the records of all treatments and medications of the study participants, our models adjusted for the use of medications for bone-related outcomes and osteoporosis. This approach adjusts for some degree of the role of bone mineral density for osteoporosis or osteoporotic fracture. Third, while considering the yearly-based incidence of bone fracture in our Cox models based on the person-year follow-up estimations is methodologically reasonable, we did not consider the season of bone fracture within a year. Incidence of some bone fracture in elderly is generally higher in cold seasons and suggested to be associated with cold temperature, ice, and snow in cold seasons. Temperature is often assessed as a confounder or effect modifier for air pollution as they both have seasonality and are highly correlated with each other. If temporally more detailed health data rather than yearly data such as the national cohort data used in this research become available, future research could aid understanding of more vulnerable periods within a year related to the risks of osteoporotic fracture from PM exposure. Lastly, the cohort data did not have information on dietary habits for each participant, so our analyses instead used the percentages of vitamin intake in each study province. This may be less accurate to assess the confounding of individuals' dietary habits on osteoporotic fracture.

Conclusion
This epidemiological study examined the impact of long-term PM 2.5 exposure on osteoporotic fracture risks in South Korea, using machine learning-based modeled PM 2.5 data. Long-term exposure to PM 2.5 measured as 6 yr moving average annual PM 2.5 concentrations was significantly associated with increases in osteoporotic fracture incidence risks in adults (age ⩾50 years). The results were robust in other examined longer exposure windows (e.g. 7-, 8-, 9-, and 10 yr exposure windows). Sub-group analysis by individual factors did not find significant risk differences by age, sex, or income. The risks did not differ by urbanicity of residential area or vitamin intake. Our results suggest that long-term exposure to PM 2.5 is a risk factor for osteoporotic fracture for Korean adults age ⩾50 years and the importance of prevention of osteoporotic fracture risks through air pollution management policies.

Data availability statement
The data cannot be made publicly available upon publication because they are owned by a third party and the terms of use prevent public distribution.