Ambient fine particulate matter air pollution and the risk of hospitalization among COVID-19 positive individuals: Cohort study.

BACKGROUND: Ecologic analyses suggest that living in areas with higher levels of ambient fine particulate matter air pollution (PM2.5) is associated with higher risk of adverse COVID-19 outcomes. Studies accounting for individual-level health characteristics are lacking.
METHODS: We leveraged the breadth and depth of the US Department of Veterans Affairs national healthcare databases and built a national cohort of 169,102 COVID-19 positive United States Veterans, enrolled between March 2, 2020 and January 31, 2021, and followed them through February 15, 2021. Annual average 2018 PM2.5 exposure, at an approximately 1 km2 resolution, was linked with residential street address at the year prior to COVID-19 positive test. COVID-19 hospitalization was defined as first hospital admission between 7 days prior to, and 15 days after, the first COVID-19 positive date. Adjusted Poisson regression assessed the association of PM2.5 with risk of hospitalization.
RESULTS: There were 25,422 (15.0%) hospitalizations; 5,448 (11.9%), 5,056 (13.0%), 7,159 (16.1%), and 7,759 (19.4%) were in the lowest to highest PM2.5 quartile, respectively. In models adjusted for State, demographic and behavioral factors, contextual characteristics, and characteristics of the pandemic a one interquartile range increase in PM2.5 (1.9 µg/m3) was associated with a 10% (95% CI: 8%-12%) increase in risk of hospitalization. The association of PM2.5 and risk of hospitalization among COVID-19 individuals was present in each wave of the pandemic. Models of non-linear exposure-response suggested increased risk at PM2.5 concentrations below the national standard 12 µg/m3. Formal effect modification analyses suggested higher risk of hospitalization associated with PM2.5 in Black people compared to White people (p = 0.045), and in those living in socioeconomically disadvantaged neighborhoods (p < 0.001).
CONCLUSIONS: Exposure to higher levels of PM2.5 was associated with increased risk of hospitalization among COVID-19 infected individuals. The risk was evident at PM2.5 levels below the regulatory standards. The analysis identified those of Black race and those living in disadvantaged neighborhoods as population groups that may be more susceptible to the untoward effect of PM2.5 on risk of hospitalization in the setting of COVID-19. Published by Elsevier Ltd.

Introduction

The coronavirus disease 2019 (COVID-19) is associated with higher risk of hospitalization, and mortality. The illness is characterized by a vigorous systemic inflammatory response and higher risk of respiratory failure and multiorgan dysfunction (Bowe et al., 2020b, Xie et al., 2021a). Ambient fine particulate matter air pollution (PM2.5) is associated with several disease entities (e.g. cardiovascular disease, diabetes, kidney disease, etc.) now recognized as susceptibility factors for adverse COVID-19 outcomes (Al-Aly and Bowe, 2020, Al-Aly and Bowe, 2021, Bowe et al., 2019b, Chan et al., 2018, Cohen et al., 2017, Huang et al., 2019, Landrigan et al., 2017, Williamson et al., 2020). Chronic exposure to higher levels of PM2.5 is also associated with immune system disturbances (maladapted immune response) — now hypothesized as a leading mechanism to explain the vast heterogeneity in clinical outcomes observed in individuals infected with COVID-19 (Glencross et al., 2020, Lucas et al., 2020, Mateus et al., 2020:eabd3871.). Taken together these observations suggest the hypothesis that COVID-19 infected individuals chronically exposed to higher levels of PM2.5 are at increased risk of hospitalization.

Several ecologic reports suggested that living in areas with high levels of PM2.5 may be associated with higher risk of hospitalization due to COVID-19 (Conticini et al., 2020, Strobl, 2020, Travaglio et al., 2020, Wu et al., 2020;6.). However, as elegantly discussed by Wu and collaborators, ecologic analyses — while informative — are limited and studies accounting for individual-level characteristics are needed but currently lacking (Benmarhnia, 2020, Bhaskar et al., 2008.20175901., Wu et al., 2020;6.). The use of study designs which account for individual-level characteristics to address the question of whether chronic prior exposure to higher levels of PM2.5 air pollution is associated with increased risk of adverse outcomes among COVID-19 infected individuals will inform the global discussion on the role of the environment in shaping COVID-19 health outcomes. Findings from these studies are important to prepare for and mitigate health loss in future epidemics or pandemics.

In this study, we built a national cohort of 169,102 US Veterans with a COVID-19 positive test and aimed to examine the association between chronic exposure to PM2.5 and risk of hospitalization. We built a non-linear exposure–response function to characterize the shape of the association between PM2.5 and risk of hospitalization across the spectrum of PM2.5 concentrations experienced in the US. We then examined the association in prespecified subgroups and conducted formal interaction analyses to identify sensitive populations.

Materials and methods

Data sources

The Department of Veterans Affairs (VA) has developed a comprehensive COVID-19 Shared Data Resource (CSDR) that continuously collects, validates, and updates information on the characteristics of COVID-19 positive patients as well as COVID-19 testing, death, mechanical ventilation use, and other hospital utilization both within and outside of the VA healthcare system (Department of Veterans Affairs, 2020). To supplement the COVID-19 CSDR, we additionally obtained data from the VA Corporate Data Warehouse (CDW) Outpatient and Inpatient Encounters domains to collect diagnosis and procedure codes for identifying outcomes and clinical characteristics. The CDW Laboratory Results domain was used to collect laboratory data for certain clinical characteristics. Vital signs measurements were collected from the CDW Vital Signs domain; smoking status was collected from the CDW Health Factors domain; and demographic information was obtained from the CDW SPatient domain and VA Vital Status Databases (Bowe et al., 2016, Cai et al., 2021, Xie et al., 2016:ASN., The Department of Veterans Affairs, 2014). The Planning Systems Support Group (PSSG) Enrollee File provided geocoded residential addresses with the corresponding latitude and longitude, as well as the Federal Information Processing Standard code of the county of residence (VA Information Resource Center 2015). In VA PSSG data, 90% of participants' geocoded residential locations are based on mapping to the participant’s street address or an interpolation based on the building number on their street and the location of that street, resulting in high resolution residence information. The remainder is geo-coded to the ZIP code tabulation area centroid (VA Information Resource Center 2015). Data on PM2.5 was obtained from satellite-based PM2.5 estimates V4.NA.02.MAPLE (Mortality–Air Pollution Associations in Low-Exposure Environments) for the year 2018 (van Donkelaar et al., 2019). Data are generated by combining Aerosol Optical Depth (AOD) retrievals from the NASA MODIS and MISR instruments with AOD and surface-to-column information from the GEOS-Chem chemical transport model, and subsequently statistically calibrating estimates to regional ground-based observations of total PM2.5 mass, achieving an R2 = 0.82. The 2015 Area Deprivation Index (ADI) contains a measure of a county’s socioeconomic deprivation, constructed from the domains of income, education, employment, and housing quality (Kind and Buckingham, 2018, University of Wisconsin School of Medicine Public Health, 2013). We additionally utilized data from the 2020 County Health Rankings (CHR) (Remington et al., 2015), a curated set of county contextual determinants of health, to obtain county-level population percentage with limited access to healthy food, percentage with adequate access to exercise opportunities, and percentage of adults reporting excessive drinking. Population density at the county-level was obtained from the US Census Bureau’s Small Area Income and Poverty Estimates data. Normalized difference vegetation index (NDVI), an indicator for green space, was obtained from the National Oceanic and Atmospheric Administration (NOAA) Climate Data Record of Advanced Very High-Resolution Radiometer Surface Reflectance (Vermote, 2019).

Identification of COVID-19 affected individuals and cohort construction

A flowchart of cohort construction may be seen in Fig. A.1 in Appendix A. We selected those with a record of testing COVID-19 positive between March 2, 2020 and January 31, 2021, followed until February 15, 2021 (N = 174,661). The CSDR identifies COVID-19 positive participants based on the VA’s National Surveillance Tool, which provides near real time information on COVID-19 cases based on laboratory results and natural language processing applied to clinical notes (which are subsequently vetted by human review). The VA utilizes identification methods, starting March 14, 2020, including laboratory tests that comply with Center for Disease Control and Prevention (CDC) standards, such as the 2019-nCOV RT-PCR Diagnostic Panel and the SARS-CoV-2 Multiplex Assay, or human-confirmed case review. Prior to this time the VA employed interim CDC case definitions, where presumptive positive cases were classified by having a positive laboratory testing for SARS-CoV-2 RNA, identified as positive after undergoing human review, which were then further verified by compliance with workload codes. We then selected those who had a recorded residential location in the year prior to testing positive that was geospatially linkable with exposure (n = 172,981), and data on all covariates (n = 172,786). We then excluded those with PM2.5 exposure levels at the 1% tails of the distribution to reduce influence of very low/high values, resulting in a final analytic cohort of 169,102.

Outcome and exposure assessment

A COVID-19 hospitalization was defined by the first hospital admission between 7 days prior to and 15 days after their first COVID-19 positive date. This timeframe of assessment enabled a period of follow-up post testing positive to allow for capture of potentially related hospitalizations and allow for a pre-test period to capture hospitalizations that happened before the presence of COVID was detected. Estimates of annual average ground level PM2.5 total mass over the contiguous US in 2018 were available at approximately 1 km2 resolution. Gridded PM2.5 exposure data with latitude and longitude coordinates of geographic location provided bounds of each geographic grid with corresponding PM2.5 mass; using the latitude and longitude of each participant’s residential street address in the year prior to the record of positive COVID-19 status, we assigned the Veteran an exposure based on the PM2.5 grid the geocoded street address was located in.

Covariates

Covariates included factors that could potentially confound the association of PM2.5 air pollution and risk of hospitalization based on evidence from prior studies (Crouse et al., 2017, Kloog et al., 2013, Laden et al., 2006, Miller et al., 2007). State of each participant’s residential location was included to adjust for differences in State-level composition of PM2.5, State based differences in COVID-19 procedures and other policies that may affect the outcome, and geospatial trends not captured by other covariates (Bowe et al., 2018a, Chan et al., 2015). We additionally included participant demographics of age, race, and sex, and a behavioral factor of smoking status. Contextual characteristics included county ADI, population density, percentage with limited access to healthy food, percentage with adequate access to exercise opportunities, percentage of adults reporting excessive drinking in the county they live, NDVI (a measure of green space), and election margin of victory (Xie et al., 2021b). NDVI measures surface vegetation by comparing red and near-infrared spectral bands, from satellite grids at 0.05° by 0.05° resolution, and ranges from −1 (water), to 0 (bare ground), and maximally 1 (dense vegetation), where an average of the top three maximal values during summer months (June, July, and August) in 2010 was taken (Crouse et al., 2017, James et al., 2016, Orioli et al., 2019). Election margin of victory represented the percentage difference in all votes in a county that were Republican vs. Democrat in the 2020 presidential election (McGovern et al., 2020). Contextual characteristics were assigned based on participants’ residential location in 2019. More detailed definitions of the contextual characteristics and their data sources are provided in Table A.1 in Appendix A. Characteristics of the pandemic included month of testing positive, COVID-19 testing rates, positivity rates, and hospital bed occupancy rates. COVID-19 testing rates were defined as the number of COVID-19 tests recorded in the health system that the participant tested positive at during the month of testing positive, while positivity rates were the percentage of these tests that were positive. Hospital bed occupancy rates was defined by the average number of hospital beds occupied out of all beds available in the hospital station the participant was admitted to, or otherwise tested postive at, over the week prior to admittance. Race was categorized as Black, Other, and White on the basis of self-report during clinical encounters, supplemented by administrative data (Sohn et al., 2006). Smoking status was identified based on work by McGinnis et al. (McGinnis et al., 2011). In brief, an algorithm validated against self-reported smoking survey data converts text entries from Health Factors data into categories of never, former, and current smoker (for instance text of “Former” or “Previous” would be assigned as “former smoker”). The report prior but closest to baseline was used; missing smoking status was assessed as never smoker. All continuous covariates were treated as restricted cubic splines.

Statistical analyses

Demographic, contextual, and pandemic characteristics overall and by PM2.5 quartile are presented as frequency (percentage), mean (standard deviation), and median (interquartile range) as appropriate. We created geospatially smoothed maps of the PM2.5 exposure in 2018 and unadjusted hospitalization rates in our cohort using Besag-York-Mollié models (Besag et al., 1991). Adjusted Poisson regression models were then constructed to investigate the association of PM2.5 exposure with risk of hospitalization within the 7 days prior to, and 15 days after, a positive COVID-19 test. We present an unadjusted model and a model adjusted for State, demographic and behavioral characteristics, contextual characteristics and characteristics of the pandemic including month of testing positive, health system testing capacity, positivity rates, and hospital bed occupancy rates. Standard errors were generated using a robust sandwich variance estimator (Knol et al., 2012). In additional analyses we conducted pooled Poisson regression, where time periods were defined as: −7 to 0 days before testing positive, 1 to 8 days after testing positive, and 9 to 15 days after testing positive. Participants were censored at death. Standard errors were generated by 1000 iteration bootstrap. In consideration of the dynamicity of the pandemic, we conducted subgroup analyses by time periods of March 2, 2020 through May 31, 2020, June 1, 2020 through October 31, 2020, and November 1, 2020 through January 31, 2021, temporally aligning periods with surges in cases in the United States. To assess potential non-linearity, we then employed the Shape Constrained Health Impact Function (SCHIF) algorithm, adjusting for all covariates (Bowe et al., 2019b, Burnett et al., 2018, Nasari et al., 2016). The SCHIF employs a series of transformations that result in monotonic exposure response curve shapes including linear, log-linear, exponential, sigmoidal, and more through a logistic weighting function. Amongst transformations tested, an optimal model was selected through best model fit (measured by the log-likelihood), and an ensembled model of the top three best fitting models, weighted by the log-likelihood, was additionally constructed.

We assessed if the risk of hospitalization within the 7 days prior to, and 15 days after, a positive COVID-19 test was associated with PM2.5, varied by participant demographic characteristics of age, Black vs. White race, sex, and ADI as a means of identifying potentially vulnerable populations. With adjusted models we conducted subgroups analyses (reported in supplement) and formal effect modification analyses, and report measures of interaction on the multiplicative and additive scale (VanderWeele and Knol, 2014). We report the relative risk from the interaction model as the association of PM2.5, with and without the presence of the interaction variable, with risk of hospitalization, along with the p-value for the interaction, in a forest plot (VanderWeele and Knol, 2014). To exclude the potential influence of outlier values, in formal effect modification analyses, we trimmed ADI tails where there were low event rates; we report the risk of hospitalization associated with PM2.5 above and below the median age and ADI distribution in the overall cohort.

Negative exposure control and sensitivity analyses

We tested ambient sodium concentration as a negative exposure control where there is no prior evidence nor biologic plausibility of an association between atmospheric sodium and hospitalization (Lipsitch et al., 2010). The negative exposure control investigates whether shared biases, in this instance possible trends in exposure misclassification on the basis of air-monitoring station measurements, may result in unexpected associations. Ambient sodium concentration was available from Environmental Protection Agency monitor station data, where measured sodium levels were linked to a participant’s latitude and longitude of residence based on the latitude and longitude of the nearest air monitoring station measuring sodium (restricted to those living within 30 miles or less). We additionally tested a negative outcome control of inpatient diagnoses of a composite outcome of accidental injuries and poisonings.

In order to test the consistency of results we conducted several sensitivity analyses (Bowe et al., 2018b, Bowe et al., 2018c, Miller et al., 2007, Wu et al., 2020;6:eaba5692.). (A) In consideration of applying exposure estimates from 2018, we selected those who had a record of residential location in 2018 and who had no record of moving residence in 2018 and 2019. (B) To test whether potential variation in regional characteristics may have confounded the association of PM2.5 with outcomes, we estimated whether within-city variation in PM2.5 is associated with outcomes. This model constructs a city average PM2.5 (defined as the average PM2.5 level in a metropolitan statistical area), and a within city difference from that average (defined as the within city effect), to assess whether, given a similar city average, within city differences are associated with differences in risk of outcomes (Bowe et al., 2017, Bowe et al., 2018b, Bowe et al., 2018c; Miller et al., 2007). (C) Alternatively, to test consistency with the outcome definition, we additionally assessed the association of PM2.5 with risk of hospitalization within the 7 days prior to, and, separately, 30 days and 60 days after, testing positive for COVID-19. Participants were restricted to those enrolled as of January 15, 2021, and December 15, 2020, respectively. (D) Because hospitalization that occurs before a COVID-19 positive test may not be associated with COVID-19 status, we exclude those that had a hospitalization before their recorded date of COVID-19 positive test. (E) We used an alternate exposure definition of the 3-year average based on a participant’s location from 2016 to 2018. (F) We alternatively accounted for demographic, behavioral, contextual confounders using an inverse probability of treatment weighting approach. Further details, including assessment of normality (Fig. A.2) and SAS syntax are provided in the Appendix. (G) To investigate sensitivity of results to geography, analyses were rerun in separate models that excluded participants living in geographic regions of the Northeast, Midwest, South, and West. (H) We additionally adjusted for baseline clinical characteristics, detailed in the Appendix.

All analyses were done using SAS Enterprise Guide 8.2. Geographic data visualization was done using R 3.5.1. No imputation was performed unless otherwise specified for sensitivity analysis. This study was approved by the Saint Louis Veterans Affairs Health Care System Institutional Review Board. The funders of this study had no role in study design; collection, analysis, and interpretation of data; writing the report; and the decision to submit the report for publication.

Results

There were 169,102 US Veterans with a record of a COVID-19 positive test in the overall cohort. The demographic, behavioral, and contextual characteristics in the overall cohort and by PM2.5 quartile are presented in Table 1 . The median PM2.5 exposure was 7.4 (interquartile range (IQR): 6.6, 8.5 µg/m3). A map of the annual average PM2.5 in the contiguous United States in 2018 is presented in Fig. A.3. Those living in areas in the highest quartile of PM2.5 were more likely to be Black, a current smoker, and live in areas that had a higher population density, more access to exercise opportunities, and a lower NDVI (Table 1). Overall, there were 25,422 (15.0%) hospitalizations within the 7 days prior to, and 15 days after, a positive COVID-19 test; unadjusted rates by PM2.5 quartiles were 5,448 (11.9%), 5,056 (13.0%), 7,159 (16.1%), and 7,759 (19.4%) in quartile 1, 2, 3, and 4, respectively. Substantial geographic variation in rates of hospitalization was observed; a map of the geographic distribution of hospitalization is included in Fig. A.4.

Table 1

Demographic, behavioral, and contextual characteristics of a national US Veterans cohort who tested positive for COVID-19 overall, and by PM2.5 quartile, based on PM2.5 estimates in 2018.

Characteristic	Overall	PM2.5 Quartile 1 (4.1, 6.6 µg/m3)	PM2.5 Quartile 2 (6.7, 7.4 µg/m3)	PM2.5 Quartile 3 (7.5, 8.5 µg/m3)	PM2.5 Quartile 4 (8.6, 14.2 µg/m3)
Individual Characteristics – n (%)
Number of participants	169,102	45,798 (27.08)	38,943 (23.03)	44,399 (26.26)	39,962 (23.63)
Median age (IQR) – year	64.53 (50.81, 73.58)	67.25 (52.54, 74.34)	63.91 (50.24, 73.45)	63.47 (50.04, 73.27)	63.78 (50.34, 72.98)
Race
Black	41,302 (24.42)	5629 (12.29)	8993 (23.09)	12,738 (28.69)	13,942 (34.89)
Other	7659 (4.53)	1742 (3.80)	1396 (3.58)	1916 (4.32)	2605 (6.52)
White	120,141 (71.05)	38,427 (83.91)	28,554 (73.32)	29,745 (66.99)	23,415 (58.59)
Male sex	152,318 (90.07)	41,702 (91.06)	34,815 (89.40)	39,609 (89.21)	36,192 (90.57)
Smoking status
Current	30,382 (17.97)	7832 (17.10)	6840 (17.56)	8014 (18.05)	7696 (19.26)
Former	48,213 (28.51)	13,869 (30.28)	10,341 (26.55)	12,078 (27.20)	11,925 (29.84)
Never	90,507 (53.52)	24,097 (52.62)	21,762 (55.88)	24,307 (54.75)	20,341 (50.90)
Contextual Characteristics – median (IQR)
Population density – n/mile2	333.80 (104.70, 1117.20)	114.1 (45.00, 326.30)	258.20 (111.50, 656.30)	464.20 (194.10, 1157.90)	1189.40 (333.80, 2344.20)
Area Deprivation Index*	54.56 (43.44, 63.39)	57.19 (47.51, 67.58)	55.45 (47.23, 64.23)	54.31 (44.38, 61.36)	51.99 (39.82, 59.82)
Limited access to healthy food – %	6.48 (3.94, 9.15)	6.77 (4.28, 9.57)	6.84 (4.4, 9.66)	6.48 (4.04, 8.85)	5.69 (3.07, 8.39)
Adults reporting excessive drinking – %	18.55 (16.43, 20.24)	18.34 (16.38, 20.2)	17.83 (15.99, 19.82)	18.4 (16.29, 20.18)	19.31 (17.75, 20.46)
Access to exercise opportunities – %	86.26 (73.08, 93.81)	77.9 (65.12, 88.35)	81.2 (69.07, 92.11)	88.92 (76.5, 95.15)	91.99 (84.35, 97.65)
NDVI**	0.62 (0.49, 0.7)	0.68 (0.56, 0.75)	0.64 (0.55, 0.71)	0.61 (0.49, 0.68)	0.51 (0.38, 0.62
Election margin of victory*** - %	−1.13 (−20.76, 28.38)	20.04 (−5.62, 40.08)	11.65 (−14.53, 34.17)	−5.09(−21.10, 21.41)	−18.26 (−38.12, –23.50)
Pandemic Characteristics – median (IQR)
Numerical month of testing positive	10 (4, 12)	10 (5, 12)	9 (4, 11)	9 (4, 11)	9 (4, 12)
Health system testing capacity†	3633 (2153, 7584)	3564 (2032, 7223)	3486 (2153, 7223)	3486 (1995, 6848)	4431 (2300, 9550)
Rate of positivity in the health system† - %	22.88 (14.53, 29.98)	22.88 (14.45, 30.15)	22.74 (14.26, 29.98)	23.71 (14.63, 30.15)	22.88 (14.51, 28.11)
Hospital bed occupancy in the health system †† %	38.21 (29.44, 45.32	38.65 (28.45, 44.54)	37.77 (29.33, 44.49)	37.77 (29.39, 44.80)	38.96 (31.61, 48.92)
Outcome – n (%)
Hospitalization	25,422 (15.03)	5448 (11.90)	5056 (12.98)	7159 (16.12)	7759 (19.42)

IQR, interquartile range; NDVI, normalized difference vegetation index.

SI unit conversion: 1 mile = 1.609 km.

Area deprivation ranges from low to high disadvantage of 0–100.

NDVI ranges from -1 (water) to 0 (ground, no vegetation) to 1 (high vegetation).

Represents the percentage of votes as Republican-Democrat votes out of all votes in the county.

In the calendar month of testing positive for COVID-19.

In the week prior to the date of testing positive for COVID-19.

The association between PM2.5 and risk of hospitalization in COVID-19 positive individuals

We then estimated the association between long-term PM2.5 exposure and risk of hospitalization within the 7 days prior to, and 15 days after, a positive COVID-19 test among COVID-19 positive individuals. In unadjusted models, a one IQR increase in PM2.5 (1.9 µg/m3) was associated with 1.18 times (95% confidence interval (CI): 1.17, 1.19) the risk of hospitalization in COVID-19 (Table 2 ). A model adjusted for geographic region, demographic and behavioral factors, contextual characteristics, and pandemic characteristics, showed an association of PM2.5 with risk of hospitalization (relative risk (RR) = 1.10 (95% CI: 1.08, 1.12)). Estimation of the pooled relative risk (Table 2) also suggested an association between PM2.5 and risk of hospitalization (1.12 (95% CI: 1.09, 1.15)). The association was evident in the three time periods corresponding to the three waves the pandemic (March to May 2020, June to October 2020, and November 2020 to January 2021) (Table 2). Application of non-linear exposure response function modeling suggested that the association of PM2.5 with risk of hospitalization was higher with higher levels of PM2.5, and the rate of risk increase attenuated at higher levels of PM2.5 (Fig. 1 ).

Table 2

Association between PM2.5 (every one IQR (1.9 µg/m3) increase) and hospitalization in COVID-19 affected individuals among a national cohort of US Veterans who tested positive for COVID-19, based on PM2.5 estimates in 2018 and hospitalization from March 2, 2020 until January 31, 2021.

Model	RR (95% CI)
Poisson
Unadjusted	1.18 (1.17, 1.19)
Adjusted	1.10 (1.08, 1.12)
Pooled Poisson
Unadjusted	1.18 (1.17, 1.20)
Adjusted	1.12 (1.09, 1.15)
Time period (waves of the COVID-19 pandemic)
Wave 1: March 2, 2020 to May 31, 2020	1.06 (1.01, 1.13)
Wave 2: June 1, 2020 to October 31, 2020	1.08 (1.03, 1.12)
Wave 3: November 1, 2020 to January 31, 2021	1.12 (1.09, 1.15)
Negative controls
Exposure to ambient sodium IQR*	1.00 (0.99, 1.00)
Outcome of accidental injuries and poisoning	1.01 (0.92, 1.10)

RR, relative risk; CI, confidence interval; ADI, area deprivation index; NDVI, normalized difference vegetation index; IQR, interquartile range.

Models adjusted for State, age, race, sex, smoking status, ADI, population density, percentage with limited access to healthy food, percentage with adequate access to exercise opportunities, percentage of adults reporting excessive drinking in the county they live, NDVI, election margin of victory, month of testing positive, health system testing capacity, COVID-19 test positivity rates, and hospital bed occupancy rates.

Negative exposure control model (fully adjusted): N = 91,262, IQR = 0.064 µg/m3.

Fig. 1

Non-linear exposure response curve of the association between PM2.5 and risk of hospitalization among a national cohort of US Veterans who tested positive for COVID-19. Bands represent the 95% confidence intervals. The optimal and ensembled models are presented. Results were obtained from an adapted Shape Constrained Health Impact Function Model applied to Poisson regression. Model was adjusted for State, age, race, sex, smoking status, ADI, population density, percentage with limited access to healthy food, percentage with adequate access to exercise opportunities, percentage of adults reporting excessive drinking in the county they live, NDVI, election margin of victory, month of testing positive, health system testing capacity, positivity rates, and hospital bed occupancy.

Negative controls and sensitivity analyses

In order to test the consistency of our design to potential latent biases that may result in spurious associations, we first employed a negative exposure control, ambient sodium. There was no evidence of an association (null association) between ambient sodium and risk of hospitalization (RR: 1.00 (95% CI: 0.99, 1.00) for every one IQR increase in sodium, Table 2). We then examined the association between PM2.5 and risk of accidental injuries and poisonings (as negative outcome control); the results suggested no evidence of an association.

We then conducted the following sensitivity analyses to test the robustness of our findings:

(A) Because change in residential address may result in misclassification of exposure and/or other contextual covariates, we examined the association in cohort participants who did not change residential address (N = 144,651). (B) Because variation in geographic characteristics may confound the association between PM2.5 and health outcomes, we estimated the within city effect (N = 137,281). (C) To further examine the temporality of COVID-19 positivity and hospitalization, we restricted the time window of the hospitalization to within the 7 days prior to, and 30 days after, a positive COVID-19 test (N = 157,298). We also examined the association between PM2.5 and risk of hospitalization within 60 days as an alternative outcome definition (N = 114,658). (D) Because hospitalizations that occurred before a positive COVID-19 status may not have been related to COVID-19, we investigated the association in cohort participants who had a COVID-19 positive status before hospitalization (N = 167,330). (E) We used an alternative definition of exposure by defining it as the average PM2.5 from 2016 to 2018 (N = 151,935). (F) We applied inverse probability of treatment weighting by the generalized propensity score (N = 163,090), where covariate balance is presented in Fig. A.5. (G) To investigate sensitivity of our results to geography, we conducted separate analyses where we removed participants from the Northeast, Midwest, South, and West. (H) We additionally adjusted for baseline clinical characteristics (characteristics presented in Table A.2). The association remained consistent in direction and magnitude across all the sensitivity analyses, consistently yielding estimates and confidence limits congruent with the primary analyses (Table 3 ).

Table 3

Sensitivity analyses of the association between PM2.5 (every one IQR (1.9 µg/m3) increase) and hospitalization in COVID-19 affected individuals among a national cohort of US Veterans who tested positive for COVID-19, based on PM2.5 estimates in 2018 and hospitalizations from March 2, 2020 until January 31, 2021.

Sensitivity analyses	RR (95% CI)
Non-movers (n = 144,651)	1.10 (1.08, 1.13)
Within-city effect (n = 137,281)	1.10 (1.07, 1.13)
30-day hospitalization (n = 157,298)	1.10 (1.08, 1.13)
60-day hospitalization (n = 114,658)	1.11 (1.08, 1.13)
Excluding COVID-19 positive status after hospitalization (n = 167,330)	1.08 (1.05, 1.10)
3-year PM2.5 average (n = 151,935)	1.10 (1.08, 1.12)
Inverse probability of treatment weighting (n = 163,090)	1.12 (1.10, 1.14)
Exclusion of participants in Northeast (n excluded = 20,255)	1.07 (1.05, 1.10)
Exclusion of participants in Midwest (n excluded = 42,853)	1.03 (1.01, 1.06)
Exclusion of participants in South (n excluded = 77,658)	1.13 (1.09, 1.16)
Exclusion of participants in West (n excluded = 28,366)	1.12 (1.08, 1.16)
Additionally adjusted for clinical characteristics	1.09 (1.07, 1.11)

Adjusted for State, age, race, sex, smoking status, ADI, population density, percentage with limited access to healthy food, percentage with adequate access to exercise opportunities, percentage of adults reporting excessive drinking in the county they live, NDVI, election margin of victory, month of testing positive, health system testing capacity, positivity rates, and hospital bed occupancy rates. Sample size for analyses not in the full sample are included where relevant. RR, relative risk; CI, confidence interval; ADI, area deprivation index; NDVI, normalized difference vegetation index.

*Regions defined as: Northeast (CT, MA, ME, NH, NJ, NY, PA, RI, VT); Midwest (IA, IL, IN, KS, MI, MN, MO, NE, ND, OH, SD, WI); South (AL, AR, DC, DE, FL, GA, LA, MD, MS, NC, OK, SC, TN, TX, VA, WV); West (AZ, CA, CO, ID, NM, MT, NV, OR, UT, WA, WY).

The association between PM2.5 and risk of hospitalization in COVID-19 positive individuals by participant demographic characteristics

To identify potentially susceptible population subgroups, we conducted subgroup and formal effect modification analyses of the association between PM2.5 and the risk of hospitalization within the 7 days prior to, and 15 days after, a positive COVID-19 test among COVID-19 infected participants. Results of subgroup analyses, including stratum specific estimates, and measures of effect modification, are presented in Table A.3. Results from formal effect modification analyses, which consider differences in baseline risks, suggested that Black people had a higher risk of hospitalization associated with a one IQR increase in PM2.5 than White people (RR = 1.17 vs. 1.12, p = 0.045) (Fig. 2 ). Those who lived in an area with a higher ADI also had higher risk of hospitalization associated with a one IQR increase in PM2.5 than those living in an area with lower ADI (1.15 vs. 1.08 above and below the median, p < 0.001). Increased risk by effect modification on the additive scale was also observed by race and ADI. There was no evidence of interaction by age or sex on the multiplicative scale (Fig. 2).

Fig. 2

Effect modification of the association between PM2.5 and risk of hospitalization by age, race, gender, and ADI among a national cohort of US Veterans who tested positive for COVID-19. Results were obtained from individual models that incorporated interaction terms between PM2.5 and the effect modifier being investigated. Relative risks are presented for every one IQR (1.9 µg/m3) increase in PM2.5. Risk associated with a one IQR PM2.5 above and below the median age and ADI are presented for interactions with continuous effect modifiers.

Discussion

In this cohort of 169,102 US Veterans with a COVID-19 positive test, exposure to higher levels of PM2.5 was associated with higher risk of COVID-19 hospitalization. The association was evident in wave 1, wave 2, and wave 3 of the pandemic. The risk was evident at PM2.5 concentrations below the 2012 National Ambient Air Quality Standards for annual average PM2.5 (NAAQS 2012). The association was stronger (by effect modification in formal interaction analyses) in those of Black race, and those living in areas with high measures of socioeconomic deprivation.

The results show that PM2.5 was associated with higher risk of hospitalization among COVID-19 infected individuals — suggesting that prior exposure to air pollution may have been associated with more severe COVID-19, necessitating hospitalization. Our conceptual framework underlying the putative association between PM2.5 and hospitalization in COVID-19 is provided in Fig. A.6 in the Appendix. We note that further adjustment for clinical comorbid conditions known to be associated with PM2.5, and recently characterized as risk factors for adverse outcomes in COVID-19, did not change the estimates (main models in Table 2, and sensitivity analyses in Table 3 which additionally adjusted for clinical comorbidities yielded similar point estimates). These observations suggest that the association between PM2.5 and risk of hospitalization in COVID-19 may have been independent of these shared characteristics (known to be associated with both PM2.5 and COVID-19 outcomes and accounted for in our analyses), providing plausibility to hypothesize the presence of a putative alternative pathway. A proposed mechanism may involve disturbances in immune system response in the setting of chronic prior exposure to PM2.5 which may then increase the risk of adverse outcomes in people infected with COVID-19. This hypothetical framework is supported both by studies predating the pandemic showing that exposure to PM2.5 may lead to maladaptive immune system and by a recent wave of studies during this pandemic suggesting that a misfiring immune response may explain the vast heterogeneity in clinical outcomes observed in individuals infected with COVID-19 (Glencross et al., 2020, Lucas et al., 2020, Mateus et al., 2020:eabd3871.).

We characterized an exposure–response function which showed that risk of hospitalization increased with increasing concentrations of PM2.5, and that risk was present at concentrations below the 2012 National Ambient Air Quality Standards for annual average PM2.5 (NAAQS 2012). Overall, the shape of the exposure–response function is consistent with functions representing the relation between PM2.5 and other health outcomes (Bowe et al., 2019a, Bowe et al., 2019b, Bowe et al., 2020a, Bowe et al., 2020c).

Our results suggest that Black race and high ADI strengthened (through effect modification) the association between PM2.5 and risk of hospitalization — suggesting that for the same level of PM2.5  exposure, Black individuals and those living in disadvantaged neighborhoods (areas of high ADI) may be more susceptible and exhibit greater vulnerability (higher risk of hospitalization) in the setting of COVID-19 positivity. These observations should be juxtaposed and interpreted in the context of two other realities: (1) a body of evidence suggesting that Black people and disadvantaged communities are chronically exposed to higher level of air pollution in the US (Colmer et al., 2020, Tessum et al., 2019), and (2) an evolving understanding of substantial disparities in the risk of COVID-19 infection by Black race and socioeconomic condition. The confluence of environmental injustice, disparities in risk of COVID-19 infections, and the heightened risk of COVID-19 infection severity (hospitalization) among Black people and people living in socioeconomically disadvantaged neighborhoods — brought to prominence in this global pandemic — illustrates the compound risk of these synergistic forces on health outcomes of the impacted population groups.

This study has several limitations. Although the national VA cohort provides the ability to account for relevant individual-level characteristics, Veterans are generally older White males, and results may not be generalizable to other populations. Observational studies are often subject to the possibility of residual confounding due to unmeasured covariates; although we employed analytic strategies including application of negative exposure control and estimated within-city effects to reduce this concern, we cannot completely rule out the possibility of residual confounding. Although we subjected our analyses to challenge by multiple sensitivity analyses, including restricting the cohort to those who did not change residential address, misclassification of exposure and other relevant covariates may exist. We also note that our conceptual framework underpinning the hypothesis and research question considered chronic prior exposure to PM2.5 as the exposure of interest, and as such we estimated the association between chronic prior exposure to PM2.5 and risk of hospitalization; it is possible that during the course of the pandemic and in response to lockdown policies and stay at home orders, PM2.5 may have declined and then rebounded in the US. These recent acute (and relatively shorter term) fluctuations in PM2.5 may also have implications on COVID-19 outcomes but are not considered in this analysis (Pozzer et al., 2020). Our analyses did not account for potential geographic variability in the composition and toxic content of PM2.5. Interaction analyses of some subgroups may have been limited in statistical power due to lower population size in our cohort.

The study has several strengths. This is a national cohort of 169,102 US Veterans who tested positive for COVID-19. The VA operates the largest nationally integrated healthcare system in the US; VA data is national in scope, and the range of exposure in this study spans the spectrum of PM2.5 concentrations experienced by individuals in the US. Furthermore, VA data contains extensive records of individual-level clinical characteristics and geocoded residential addresses. We complemented this rich dataset though linkage with external data sources to additionally account for contextual characteristic (e.g. residential greenness, neighborhood disadvantage, etc.). The conceptual framework underlying the research question is supported by the most up-to-date studies on COVID-19. We employed a state-of-the-art non-linear modeling approach (SCHIF) to characterize the exposure response function — enhancing the policy relevance of the analyses. We developed and tested a negative exposure control and negative outcome control to investigate concerns about spurious associations and examined the consistency of results in multiple sensitivity analyses.

Conclusions

In sum, our results suggest that exposure to higher levels of PM2.5 was associated with higher risk of hospitalization among individuals with a COVID-19 positive test. The association was present in each wave of the pandemic. The risk was evident at PM2.5 concentrations well below the National Ambient Air Quality Standards for annual average PM2.5. The association was stronger among those of Black race and those living in socioeconomically disadvantaged neighborhoods. Further investigation is needed to gain a better and deeper understanding of how air pollution affects COVID-19 severity. More broadly, the findings underscore the interconnectedness of several processes — an emerging infectious disease which has caused a global pandemic, the environment, racial and non-racial socioeconomic disparities, and human health. A broader understanding of the dynamic interplay of these processes is needed to optimize the health of humans and the planet before the onset of the next pandemic.

Declarations

Ethics approval and consent to participate

This research project was reviewed and approved by the Institutional Review Board of the VA Saint Louis Health Care System.

Availability of data and materials

VA data is not publicly available as it is Health Insurance Portability and Accountability Act protected but is available by request from the US Department of Veterans Affairs. PM2.5 data is available at . County Health Rankings are available at . The Area Deprivation Index is available at . The NDVI is available at . 2020 election data is available at .

Funding

This research was funded by the United States Department of Veterans Affairs and the Institute for Public Health at Washington University in Saint Louis, Missouri, USA (for ZAA), and by the McKelvey School of Engineering at Washington University in St. Louis (for RVM). The funders of this study had no role in study design; collection, analysis, and interpretation of data; writing the report; and the decision to submit the report for publication.

Authors' contributions

BB and ZAA contributed to the development of the study concept and design. BB, MC, AvD and RVM contributed to data acquisition. BB, YX, MC and ZAA contributed to data analysis and interpretation. BB, YX, and MC contributed to statistical analysis. BB, AKG, and ZAA drafted the manuscript. Critical revision of the manuscript was contributed to by BB, YX, MC, AKG, AvD, RVM, RB, and ZAA. Administrative, technical, and material support was provided by ZAA. ZAA contributed supervision and mentorship. Each author contributed important intellectual content during manuscript drafting or revision and accepts accountability for the overall work by ensuring that questions pertaining to the accuracy or integrity of any portion of the work are appropriately investigated and resolved. ZAA takes responsibility that this study has been reported honestly, accurately, and transparently; that no important aspects of the study have been omitted, and that any discrepancies from the study as planned have been explained. All authors approved the final version of the report.

Declaration of Competing Interest

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