Pre-Diagnosis Air Pollution Exposure and Survival in Stage 1A NSCLC: A U.S. Cohort Study

1. Introduction

Several modifiable social determinants of health (SDOH) that improve lung cancer survival exist beyond smoking cessation (Strickland et al., 2015; Baxter et al., 2013; Eckel et al., 2016a). Ambient air pollution is a modifiable determinant of lung cancer survival(Aksoy, 1980; McKeon et al., 2022a; Pyo, Kim, and Kang 2022) and yet research exploring the dose- relationship association of ambient air pollution on lung cancer incidence and mortality in the United States (U.S.) is limited(McKeon et al. 2022b; C. S. Liu et al. 2023; Pyo, Kim, and Kang 2022; Eckel et al. 2016a). Air pollutants affect a specific type of lung cancer histology(Lamichhane et al. 2017; H. C. Lee et al. 2022; Moon et al. 2020; Pyo, Kim, and Kang 2022) therefore, it is essential to focus specifically on histology type and specific clinical stages of lung cancer to determine survival outcomes(McKeon et al. 2022a; Eckel et al. 2016a), but only few studies attempt to take that into account. Only one study to date has established a dose-response relationship between localized lung cancer survival and ambient air pollution exposure(Eckel et al. 2016a), but the study did not account for weather components that might affect exposure levels within the vicinity (Kim, Sheppard, and Kim 2009; Zanobetti & Peters 2015; Zheng et al. 2019; Tian et al. 2021).

Air pollutant levels differ geographically, affecting the level of exposure among patients in a longitudinal study. Changes in weather conditions also facilitate chemical reactions between primary pollutants (NO2, SO2, CO, and PM) and other atmospheric chemicals, resulting in secondary pollutant production. The weather components, such as temperature maximum, are also correlated with air pollutants as the rise in air pollutants aids in the urban heat island phenomenon(Wang, Guo, and Han 2021) . Hence, accounting for secondary pollutants such as ozone and weather components such as temperature maximum might provide biased estimation results in a given study context. Therefore, it is vital to understand the complex interaction of air pollutants in the presence of weather components, such as precipitation, snow, and temperature, which affect specific exposure levels and determine the survival outcomes of stage 1A TN0M0 NSCLC (Ngarambe et al., 2021).

Some studies in the literature that identify the dose-response relationship between ambient air pollution and lung cancer survival utilize interpolation or other data techniques to replace missing pollutant levels (Eckel et al., 2016a; McKeon et al., 2022a). The drawback of interpolating or extrapolating missing pollutant values without taking into account other environmental factors might inherently misclassify exposure assignments providing uncertain estimates due to the absence of relevant information such as natural events and weather components such as snow, precipitation, temperature interaction with other spatially and temporally dependant pollutants (Y. Liu, Zhou, and Lu 2020; Oji and Adamu 2020; De Sario, Katsouyanni, and Michelozzi 2013; Zanobetti and Peters 2015; Zheng et al. 2019; Tian et al. 2021; Kim, Sheppard, and Kim 2009). Moreover, a lack of sufficient variance in values exists when utilizing other methods compared to the methods of the nearest monitoring stations to assign exposure values and determine health outcomes (Kim, Sheppard, and Kim 2009).

Several factors affect standard treatment care receipt, an important confounder in determining survival outcomes for stage 1A NSCLC(Cao et al. 2018; Raman et al. 2022; Baig et al. 2020; Shen et al. 2021). The National Comprehensive Cancer Network (NCCN) treatment guidelines are referred by about 95% of U.S. oncologists to recommend standard treatment care. Changes in these guidelines, resulting from revisions over the past several decades, may affect who receives standard care for early-stage lung cancer. Differences in standard treatment receipt exist for early-stage lung cancer by treatment guideline revision year, race, geography, and insurance status, as established in an extended prior study (Patel et al., 2024). The trends in the type of treatment receipt and air pollution levels that are spatially and temporally dependent in the presence of weather elements for a more extended study period across diverse US state counties help identify close to the true causal relationships in similar survival studies (Eckel et al. 2016b; C. Liu et al. 2023). Hence, it is also crucial to identify whether ambient air pollution has a dose-response effect on lung cancer survival outcomes depending on the type of treatment received that includes the timespan of several US national treatment guideline revision years, pre-diagnosis exposure values to account for the cumulative effect, and state counties with differential time invariable confounders accounted in statistical analysis. To our knowledge, only limited studies have aimed to identify it for the U.S. representative population(McKeon et al. 2022b; Xu et al. 2013a; Eckel et al. 2016a) however, the studies did not account for the dose-response relationship in the presence of weather components in a homogenous sample of stage 1A NSCLC TN0M0. It also did not account for other primary air pollutants, such as SO2 and CO. Finally, the studies assigned exposure values from the month of diagnosis to death rather than considering pre-diagnosis exposures. This could lead to an absence of accounting for the carry-over effect on health outcomes from before diagnosis exposure.

Therefore, we aimed to evaluate whether exposure to specific levels of air pollutants is associated with survival outcomes among patients with stage 1A TN0M0 non-small cell lung cancer (NSCLC) undergoing treatment of choice, utilizing U.S. population-based cancer data and U.S. environmental air pollution data. Does accounting for any key confounders missing in previous similar studies reduce selection bias and provide close-to-true hazard ratios? How does treatment choice affect survival outcomes in the presence of exposure to the identified air pollutants? We hypothesize that there exists a difference in all-cause mortality hazards among treated individuals exposed to high versus low air pollution levels (Eckel et al., 2016a).

2. Methods

2.1. Study Design

This retrospective cohort study compared the survival outcomes between patients exposed to higher versus lower air pollution and those receiving different treatment types (i.e., limited resection with adjuvant radiotherapy and lobectomy) in single- and multi-pollutant models similar to pre-existing limited studies (Eckel et al. 2016b; Xu et al. 2013a). The pollutant model included NO2, SO2, and CO, adjusted for precipitation, snow, and daily minimum temperature values in both the single-pollutant and multi-pollutant models. The multi-pollutant model included NO2, SO2, and CO, along with weather components, whereas single-pollutant models consisted of one primary pollutant and weather components. The pollutant models were analyzed separately for three-time intervals (one, three, and five years) pre-diagnosis exposure model, for one year, three years, and five years of survival outcomes (post-diagnosis exposure) to determine the robustness of the estimates(C. S. Liu et al. 2023; McKeon et al. 2022b).

2.2. Data Sources and Construction of Analysis Data File

The SEER 18 Research Plus, Environment data, and AHRF were used from 1988 to 2015. The SEER 18 Research Plus data access request was approved on April 18, 2022, with reference number SAR0028589, to access the data through the SEER*Stat account. The AHRF collects data from over 50 national sources, aggregated at the county level, and is compiled by the Health Resources Services Administration’s (HRSA) Bureau of Health Professions for each of the nation's counties, using publicly available data (“Area Health Resources Files” n.d.).The Surveillance Research Program (SRP) of the National Cancer Institute (NCI) Division of Cancer Control and Population Sciences (DCCPS) supports SEER. SEER collects and publishes cancer incidence and survival data for every cancer case reported from 22 U.S. geographic areas, covering approximately 48 percent of the U.S. population, through population-based cancer registries. Registries routinely collect data on patient demographics, primary tumor site, tumor morphology, and stage at diagnosis, as well as the first course of treatment and follow-up for vital status (survival) (“Overview of the SEER Program” n.d.).

Agency-pregenerated daily summary air pollutant data files from 1988 to 2015 were downloaded from the following website: aqs.epa.gov/aqsweb/airdata/download_files.html. The air pollution gases raw data downloaded included ground-level Ozone (O3), Sulfur Dioxide (SO2), Carbon Monoxide (CO), and Nitrogen Dioxide (NO2). For particulate pollutants, the raw data downloaded were Particulate Matter (PM2.5 and PM10). We initially investigated the toxic precursor benzene among Hazardous Air Pollutants (HAPs) and Volatile Organic Compounds (VOCs); however, the high rate of missing values made it unfeasible to include them in the final data analysis file. The raw data files for weather were retrieved by using the following link in the computer to access the open ftp files: ftp://ftp.ncdc.noaa.gov/pub/data/ghcn/daily/by_year/. Zip files from 1988 to 2015 were downloaded by year and unzipped to retrieve the raw files.

The construction of the data file for the final analysis is shown in Figure 1, and the sample selection process is presented in Figure 2. AHRF files were converted from software-independent archival files to software-dependent files and subsequently cleaned before being merged with SEER data, using the year and county Federal Information Processing Standard (FIPS) code as the merging criterion. Similarly, after cleaning weather and air pollution data files for the study period, they were assigned to SEER registry patients by the nearest monitoring station method, as explained in the exposure assignment section of this paper.

2.3. Statistical Analysis and Empirical Model

Descriptive statistics, Kaplan–Meier survival graphs, and the Cox regression model (Eckel et al. 2016b; Xu et al. 2013a; McKeon et al. 2022b; C. S. Liu et al. 2023)were used to determine the sample demographics and time to all-cause mortality, with right censoring due to death or study end. The model examined the association between treatment type, air pollutants, weather, and survival, as well as the interactions between treatment types and air pollutants and between weather and treatment types while adjusting for patient demographics, clinical characteristics, and time-invariant unobserved variables, including the year of diagnosis and county FIPS. The duration dependence of hazards due to unobserved heterogeneity was accounted for in the model by including the year of diagnosis and county-specific, time-invariant, unobservable factors. Single-pollutant models and multi-pollutant models were computed, adjusting for the same covariates and dummy variables to determine whether the estimates were biased due to the independent variables omitted in the unadjusted model. The final model was examined for diagnostic criteria and model fit, including testing for multicollinearity between the exposure variables. After the preliminary analysis and diagnostics, the final regression models included NO2, SO2, CO, precipitation, daily minimum temperature, and snow accumulation variables.

Kaplan-Meier survival curves and dose-response relationships between adjusted NO2, SO2, and CO hazards were plotted by pollutant quartile groups (McKeon et al., 2022b) to determine survival probabilities and dose-response relationships. Survivor functions by pollutant groups were plotted for the nearest air pollution monitors up to 30 miles, the weather station at 20 miles with 25% monthly missing values, and up to 40 miles air pollution, 20 miles weather stations, and 50% missing monthly values(McKeon et al. 2022b; Eckel et al. 2016b). STATA 16 and Microsoft Excel were used for the data analysis.

The following empirical model analyzes survival outcomes for patients treated with fixed-effect dummy variables:

H(t)=h₀(t).exp{β₁.Treatment Type_i + β₂.Patient Demographics_i + β₃.Clinical Characteristics_i + β₄.County_i + β₆.Air Pollutants_i + β₇.Weather Components_i + β₈.Air Pollutants_i × Treatment Type_i + β₉.Weather Components_i × Treatment Type_i + β₁₀.Year of Diagnosis_i}

where h0(t) is the baseline hazards, and exp(βs) is the hazard ratio or rate ratio. The variables Countyi and Year of Diagnosisi are county- and year-of-diagnosis time-invariant, unobservable factors. In the model, i indicates an individual patient. “Treatment Type” is a binary variable that takes the value “lobectomy” if the patient underwent a lobectomy and “limited resection with adjuvant radiotherapy” if the patient underwent a limited resection with adjuvant radiotherapy. Other treatment types were excluded because there were fewer observations within the radiotherapy and limited resection categories.

2.4. Sensitivity Analyses

The robustness of the effect was tested by estimating hazards using the average monthly median and maximum exposure values for one, three, and five years before and after diagnosis obtained from the corresponding daily exposure values. The confounding effect due to omitted exposure variables was assessed by running both single-pollutant and multi-pollutant models.

2.5. Ethical considerations

The University of Louisville ethics committee approved this study (IRB number 22.0281). The study is exempt under 45 CFR 46.101(b) in Category 4: Secondary research, for which consent is not required.

2.6. Sampling Strategy, Exposure Assignment, and Study Variables

2.6.1. Population and sample

The SEER 18 research plus cancer registry patients inclusion and exclusion criteria are explained in a prior published work(Patel et al. 2024), while the final included study sample is described in Figure 2 of the current paper. The final sample included patients with monthly exposure averages calculated from daily air pollution values and weather data, as well as the percentage of missing values for non-missing variables in the regression analysis, in addition to AHRF and SEER 18 files. The inclusion and exclusion criteria for weather and air pollution exposure values are described in the exposure assignment section of this paper. After preliminary analysis, patients with exposure up to five years before diagnosis were included in the final analysis and followed until death or the study cutoff from the date of diagnosis to five years after diagnosis. The reason for including these patients was to mitigate the compositional effect and misspecification error resulting from migration during the more extended study periods. Including patients post-five years after diagnosis and prior to five years before diagnosis, the exposure period is too long, a time frame that is more prone to migration chances. According to the U.S. Census Bureau's mobility data from 2017 to 2021, approximately 4% and 2% of people in the age groups 25-64 and 65, respectively, migrate to a different county. The information on the excluded sample of the study is provided in Appendix A Table A1a,b.

2.6.2. Exposure Assignment

Air pollution and weather exposure assignments for each patient are shown in Figure 3. We utilized the nearest monitor station method to assign pollutant concentration exposure values by closest monitor to each study participant’s location from the county centroid and included values of the three nearest neighboring monitors in the event one of the nearest monitors had a missing value for a given day in which case the data from second and third nearest monitor were utilized to assign exposure values (Rivera-González et al. 2015; McKeon et al. 2022b; Eckel et al. 2016b). Each patient in the final sample was assigned exposure from death or study cut-off (ten years after diagnosis) at one, three, five, and ten years before diagnosis. Exposure assignments were excluded when the nearest air pollution monitoring station was more than 40 miles away, the weather station was more than 20 miles away, and the percentage of missing monthly values exceeded 50%. Preliminary sample analysis of exposure assignments for air pollution ≤ 30 miles, weather ≤ 10 miles, and < 33.33% missing values determined a minimal sample size; therefore, the final analysis sample was least restrictive in terms of the distance of air pollution exposure assignments with the nearest monitoring station < 40 miles, weather station ≤ 20 miles, and missing monthly values ≤ 50%. We initially generated monthly values from daily values by keeping only those observations that were 10, 20, 30, and 40 miles away from the nearest monitoring station, with 50%, 33.33%, 25%, and 20% missing values for each mile within a month and calculated monthly mean, median, maximum, and interquartile range exposure (Figure 4) values for the same for up to 10 years before and 10 years after diagnosis month or until death. However, for our final analysis, we included exposure values, as mentioned previously in this paper.

Exposure assignment errors can be categorized as measurement and misspecification errors. A recent study relevant to the current study determined that long-term exposure assignment measurement errors are inevitable in epidemiological studies and are random. Although randomly present, the classical and Berkson measurement errors obtain biased results towards the null. If the studies find a statistically significant association, the estimates are smaller than the true effect size and are less likely to be undermined (Wei et al., 2022). One measure we have taken to control for larger misspecification errors is to restrict the study period to ten years, encompassing five years before and five years after diagnosis.

2.6.3. Independent Variables

Each weather and air pollution component continuous variable included the yearly average of monthly averages before diagnosis exposure of each patient, up to the time of death or study cut-off (whichever occurred first). The categorical treatment type variable included two categories: lobectomy and limited resection with adjuvant radiotherapy. Due to the limited number of radiotherapy observations and the few resection cases, including those two categories, it was not feasible to include them in the analysis. Surgery codes for wedge resection and segmentectomy were not differentiated in the data prior to 1998 (Razi et al., 2016). Hence, we adopted a conservative approach and combined the two types of surgery codes into one category, “Limited Resection,” as informed by the NCCN treatment guidelines and similar studies(Kates, Swanson, and Wisnivesky 2011; Mery et al. 2005). The radiation sequence with a variable of surgery from the data was utilized to aid in creating the treatment category of limited resection with adjuvant radiotherapy.

2.6.4. Outcome Variable

Survival time was calculated as the number of months of survival from the first diagnosis to death from any cause (all-cause mortality).

2.6.5 . Covariates

Tumor size categories were constructed as described by the American Lung Cancer Society (ALCS). Due to the limited number of observations in the category “up to 3 cm” and the absence of specific values, the category was merged with the “unknown tumor size' category (American Cancer Society, n.d.). SEER 18 Research Plus cancer registry data lacks information on tumor size before 2004, so the patients before the 2004 diagnosis had missing tumor size values. A more conservative approach was adopted in the current study, and observations with missing information were categorized into the unknown tumor size category. Likewise, for the insurance status information, no data were available before 2007, so an unknown category was constructed for insurance status information prior to 2007. Dummy variables for the county FIPS and year of diagnosis were constructed to account for time-invariant unobservable variables. The non-metropolitan rural-urban continuum category comprised small metropolitan, micropolitan, and non-core, as these three categories had very few observations, and there was not much demographic difference. Hence, the rural-urban continuum categorical variable comprises four categories: large central metro, large fringe metro, medium metro, and non-metropolitan.

3. Results

Overall, individuals exposed to above-median levels of air pollutants had a lower survival probability than those exposed to below-median levels, as reflected in the Kaplan-Meier survival estimates (Figure 5, Figure 6 and Figure 7). The single-pollutant model graphs did not appear to exhibit striking differences from their multi-pollutant counterparts, indicating the robustness of the results. Similarly, the 30 miles of air pollution nearest station values reflected similar directions of survival probability both in Kaplan-Meier survival estimates and statistical analysis.

Table 1. a Frequency statistics of study sample by above and below pollutant exposure median. b Descriptive statistics of study sample above and below pollutant exposure median

3.1. Hazards of Death Five Years After Diagnosis

The one- and three-year hazard estimates of death after diagnosis are robust, as indicated in Appendix A Tables 4 and 5. The all-cause mortality hazard of death for those exposed to NO2 increased by 4%, 6%, and 9%, with an average monthly increase of 1 ppb for one, three, and five years before diagnosis exposure, respectively (Table 2a). Those exposed to SO2 had an increase in all-cause mortality hazards by 16% and 17%, with an average increase in monthly averages of 1 ppb for one, three, and five years before diagnosis exposure. Those exposed to CO had an increase in all-cause mortality hazards of 53%, 51%, and 42%, with an average increase in monthly averages of 1 ppb for one, three, and five years before diagnosis exposure, respectively. Death hazards for those exposed to precipitation decreased by 2% and 3%, with an average monthly increase of one-tenth of a millimeter for one, three, and five years before diagnosis, respectively. Similarly, the hazards of death for those exposed to snowfall decreased by 10%, with an average monthly increase of one mm for five years before diagnosis exposure. The hazard effect modestly changed the effect size for single-pollutant models; however, the estimates remained significant.

The sensitivity analysis determined a similar effect direction, size, and statistical significance, except for one year after diagnosis. The hazards for the average maximum exposure values for NO2 and daily minimum temperature were no longer significant (Appendix A Table A2 and Table A3).

4. Discussion

The present study found that patients exposed to higher concentrations of NO2, SO2, and CO ambient air pollution before diagnosis had decreased survival after diagnosis. The results from prior similar studies are consistent (Eckel et al. 2016b; Xu et al. 2013b; McKeon et al. 2022c; C. Liu et al. 2023) with the existing study results regarding estimate direction for air pollutant NO2 in presence of weather elements and other prior excluded primary air pollutants i.e. SO2 and CO, eventhough we did not utilize intrapolation or extrapolation techniques supporting classical and berkson exposure error theories explained in exposure assignment section of this paper. It has also been determined that snowfall and precipitation decrease death events after diagnosis, which aligns with the logic that ambient air pollution concentration is lower during precipitation(De Sario, Katsouyanni, and Michelozzi 2013; Oji and Adamu 2020; Y. Liu, Zhou, and Lu 2020; Zanobetti and Peters 2015; Zheng et al. 2019; Tian et al. 2021). Although our study is the first of its kind and no relevant studies exist, other studies examining different health outcomes in the presence of air pollution exposure and survival outcomes in the absence of weather components are present. These findings align with the existing literature(Altorki et al. 2019; Rueth et al. 2012) and claim that lobectomy has increased surgery-associated morbidity post-lobectomy, if exposed to higher levels of air pollutant (C. Liu et al. 2023). Higher ambient air pollutants also affect lung function, as per a recent study which translates into increased death hazards (K. K. Lee et al. 2020).

The present study has several strengths as it utilizes key primary air pollutants such as SO2, CO, and weather components such as precipitation, snowfall, and daily minimum temperature to account for the confounding effects. Ozone and daily temperature maximum pose multicollinearity problems due to their inherent correlation with primary pollutants(Ngarambe et al. 2021) so they were excluded from the analysis. This exclusion aligns with Eckel et al. 2016 study(Eckel et al. 2016a), findings that ozone had a non-significant effect on survival outcomes, possibly because of multicollinearity. In addition, the study evaluated the effects of air pollutants and weather components before and after diagnosis exposure assignments to determine their cumulative effects.

Some of the limitations is the insufficient sample size for radiotherapy and limited resection. Therefore, determining the actual hazard rate using these treatment categories is difficult. In addition, the AHRF had significant missing values for area-level information relevant to the study, which could not be controlled for in the analysis. However, the county level and year of diagnosis dummy variables address these limitations for time-invariant unobserved variables. While differential yearl analysis i.e. one, three, and five years before and after diagnosis might help estimate if there are significant time varying confounders affecting overall estimates. Some of the missing contextual variables that could help reduce estimation bias were comorbidity score, cardiopulmonary function, lung function, hospital region, patient’s overall functional status, occupation, and surgeon expertise. Patient functional status and cardiopulmonary function are variables that seem to be negatively correlated with air pollution and weather exposure. However, they appear to be positively correlated with survival outcomes. In the absence of these variables, the derived biased estimates are underestimated. Hence, it was vital to account for key confounders in the present study. For the same reason, our study only measured associational relationships because we did not account for these identified unobserved confounders in the analysis, nor was the study designed to be a randomized control trial or natural experiment.

Moreover, the results of the current study are only generalizable to the population representative of the sample. As most monitors are present in metropolitan areas, potentially due to higher pollution levels, the results from the present study cannot be generalized to population outcomes in rural areas.

5. Implications for Practice and Policy

The survival of treated patients with stage 1A NSCLC is negatively associated with increased concentrations of ambient air pollutants such as NO2, SO2, CO, and daily minimum temperature. Hence, it is vital to implement environmental policies that control the emission or source of emission to reduce preventable deaths in stage 1A NSCLC patients with adenocarcinoma or squamous cell carcinoma histology types and other cardiopulmonary patients residing in metropolitan areas. It will not only help improve early stage lung cancer survival rates but also help reduce healthcare cost burdens due to increased air pollution exposure levels and associated reduced lung function or other complications.