Long-term air pollution and other risk factors associated with COVID-19 at the census tract level in Colorado.

Previous nationwide studies have reported links between long-term concentrations of fine particulate matter (PM2.5) and COVID-19 infection and mortality rates. In order to translate these results to the state level, we use Bayesian hierarchical models to explore potential links between long-term PM2.5 concentrations and census tract-level rates of COVID-19 outcomes (infections, hospitalizations, and deaths) in Colorado. We explicitly consider how the uncertainty in PM2.5 estimates affects our results by comparing four different PM2.5 surfaces from academic and governmental organizations. After controlling for 20 census tract-level covariates, we find that our results depend heavily on the choice of PM2.5 surface. Using PM2.5 estimates from the United States EPA, we find that a 1 μg/m3 increase in long-term PM2.5 concentrations is associated with a statistically significant 26% increase in the relative risk of hospitalizations and a 34% increase in mortality. Results for all other surfaces and outcomes were not statistically significant. At the same time, we find a clear association between communities of color and COVID-19 outcomes at the Colorado census tract level that is minimally affected by the choice of PM2.5 surface. A per-interquartile range (IQR) increase in the percent of non-African American people of color was associated with a 31%, 43%, and 56% increase in the relative risk of infection, hospitalization, and mortality respectively, while a per-IQR increase in the proportion of non-Hispanic African Americans was associated with a 4% and 7% increase in the relative risk of infections and hospitalizations. The current disagreement among the different PM2.5 estimates is a key factor limiting our ability to link environmental exposures and health outcomes at the census tract level. These results have strong implications for the implementation of an equitable public health response during the crisis and suggest targeted areas for additional air monitoring in Colorado.

Introduction

Although a full accounting of the biological and societal risk factors driving the spread of the SARS-CoV-2 virus will require years of study to be fully understood, several recent epidemiological studies have provided national assessments of the role of important risk factors of COVID-19-related morbidity and mortality. In particular, a growing body of evidence points to fine particulate matter (PM2.5) air pollution exposure as being a possible driver of increased vulnerability to complications related to COVID-19. PM2.5 is estimated to contribute more to premature mortality in the United States than any other air pollutant and has a well-established relationship with respiratory diseases and decreased respiratory function (GBD 2019 Risk Factors Collaborators, 2020). In April 2020, the Harvard School of Public Health reported that an increase of 1 μg per cubic meter (μg/m3) in long-term average PM2.5 concentration was associated with an 11% increase in the COVID-19 mortality rate at the U.S. county level (Wu et al., 2020). However, estimated effects for the same increase in PM2.5 reported elsewhere have ranged from a 4% decrease (Travaglio et al., 2021) to a 17% increase (Tian et al., 2020). Multiple systematic reviews have found that long-term air pollution may enhance COVD-19 severity, but have emphasized the need for studies that include additional covariates and account for spatial autocorrelation (Ali et al., 2021; Bourdrel et al., 2021). Further, there is no universally-accepted standard for measuring community-scale exposure to PM2.5 air pollution, and studies to date have not compared effect estimates reported from different PM2.5 exposure surfaces.

In addition, to overcome incomplete reporting of COVID-19 morbidity and mortality by race/ethnicity, two recent national studies analyzed publicly available data from U.S. counties to describe racial/ethnic disparities for non-Hispanic African American and Hispanic/Latino populations, respectively (Millett et al., 2020; Rodriguez-Diaz et al., 2020). Nationwide, these studies reported that counties with a greater percentage of non-Hispanic African American and Hispanic/Latino residents have statistically higher COVID-19 burdens. However, while these recent U.S.-based studies of COVID-19 risk factors accounted for a variety of confounding variables, they all suffered from relatively coarse spatial resolution by relying on county-level data. This limits the ability of state and local public health officials to translate these findings for decision making within their own jurisdictions. Furthermore, multiple studies have emphasized the need for higher-resolution studies to examine the relationship between COVID-19 and air pollution in more detail (Cole et al., 2020; Konstantinoudis et al., 2021; Travaglio et al., 2020; Villeneuve and Goldberg, 2020).

This analysis builds on these previous studies by constructing a statewide census tract-level analysis in the U.S. state of Colorado. We fit spatial regression models of census tract-level counts of three COVID-19 outcomes (infection, hospitalization, and mortality) using Bayesian hierarchical random effects models with long-term PM2.5 concentrations as the primary variable of interest. For the main analysis, we consider four different measures of PM2.5 air pollution exposure, along with additional census tract-level confounding factors including race/ethnicity, socioeconomic status, multiple measures of social distancing, age demographics, rates of comorbidity prevalence and hospitalization, time elapsed since the first reported case, nearby hospital capacity, and overall testing rate. The full model uses random effects to account for spatially correlated as well as uncorrelated variation across all census tracts, and includes a county-level random effect to account for uncorrelated variation across counties.

Methods

While ecological studies have important limitations, they are a valuable source of information in rapidly evolving areas of research and when individual-level data or covariates are not available (Wu et al., 2020). To overcome some of the limitations identified in previous ecological studies on COVID-19 and air quality (Villeneuve and Goldberg, 2020; Ali et al., 2021), we examine outcomes at the census tract level, rather than at the county level, include a wide range of potential covariates, and use a statistical model designed to account for spatial autocorrelation. Census tracts generally have a population between 2500 and 8000 people and can vary widely in spatial extent depending on the density of settlement (United States Department of Commerce, 2018). In Colorado, the median census tract population was 4229, and the median spatial extent was 3.6 km2 according to the 2014–2018 American Community Survey five-year estimates.

COVID-19 outcome data

We obtained census tract-level counts of COVID-19 cases, hospitalizations, and deaths from the Colorado Electronic Disease Reporting System (CEDRS) as well as tract-level counts of unique individuals tested for COVID-19 from Colorado's electronic lab reporting (ELR) system. All counts were stratified by age group and gender, and prior to aggregation to census tract, individual records were geocoded using the address information provided within the CEDRS and ELR databases. While the Colorado Department of Public Health and Environment (CDPHE) ascertains patient address information for COVID-19 cases within the CEDRS database, missing address information for individuals receiving a COVID-19 test within the ELR database can be filled in either by matching to other surveillance data sources or substituted by the address of the health care provider. In all, approximately 88% of cases and 93% of individuals with a COVID-19 test were geocoded to a census tract. For our main analysis, we computed the total number of all COVID-19 infections, hospitalizations, and deaths by census tract between March 1, 2020 and August 31, 2020. The COVID-19 data were pulled from the CEDRS and ELR databases on October 22, 2020. Figure S1 displays the crude rate of each outcome per 100,000 population by county in Colorado.

PM2.5 pollution data

The most widely accepted measurements of PM2.5 in the United States come from regulatory-grade monitoring measurements reported to the US Environmental Protection Agency (EPA). While these monitors are excellent at capturing regional air quality, these networks are too coarse to accurately characterize PM2.5 gradients within Colorado, which can vary sharply over 10 km spatial scales (CDPHE Technical Services Program, 2020). In addition, many areas of Colorado are 50 km or more from the nearest regulatory monitor. To overcome these limitations, several research groups and government agencies have developed estimated PM2.5 surfaces that combine EPA monitoring data with other sources of information, such as satellite measurements of aerosol optical depth or the output from chemical transport models (Diao et al., 2019). We identified four different PM2.5 surfaces for use in our analysis. Details about the different surfaces are presented in Table 1 .

Table 1

Characteristics of the PM2.5 surfaces used in this analysis.

Abbreviated Name	Creator	Spatial Resolution	Data Sources Included	R2 v. EPA monitors	Reference
ACAG	Atmospheric Composition Analysis Group at Washington University	0.01 ° × 0.01 ° (Approximately 1 km × 1 km)	EPA monitors, computer modeling (GEOS-Chem), satellite measurements	0.57	van Donkelaar et al. (2019)
EPA	Environmental Protection Agency	Census tract	EPA monitors, computer modeling (CMAQ)	0.73	Berrocal et al. (2012)
CSU	Pierce Group at Colorado State University	15 km × 15 km	EPA monitors	0.70	O'Dell et al. (2019)
Schwartz	Schwartz Group at Harvard University	1 km × 1 km	EPA monitors, computer modeling (CMAQ), satellite measurements, land-use	0.78	Di et al. (2019)

All four surfaces considered in this analysis are available over the entire continental United States. Three of them (ACAG, EPA, and CSU) are publicly available. The fourth (Schwartz) was provided upon request. We transformed the gridded products into census tract-level averages by assuming the PM2.5 concentration was constant across each grid cell and using an area-weighted average. Because the spatial scale of PM2.5 variation is large compared to the typical size of census tracts in Colorado, this re-interpolation step is unlikely to affect the results (Gryparis et al., 2009; Peng & Bell, 2010). For all of the surfaces, we averaged the annual PM2.5 concentrations from 2008 to 2017. Fig. 1 shows the long-term average PM2.5 concentrations for each surface over the entire state of Colorado and over the Denver metro area. While all four surfaces identify Denver and the Front Range region as having the highest average concentration of PM2.5, there are significant differences among the four surfaces, particularly over the Rocky Mountains and in the area to the northeast of Denver. While the EPA and Schwartz surfaces show the best agreement with the EPA monitors (Table 1), these monitors are not an independent measurement, since all four surfaces are built in part on the same EPA monitoring data. Comparisons against EPA monitors therefore cannot definitively answer the question of which PM2.5 surface is most accurate in Colorado. Additional information about how each surface was constructed is included in the supplemental information (Section S1), and a visual comparison of the spatial resolution of each surface with the size of census tracts is shown in Figure S2.

Fig. 1

Comparison of the four different PM2.5 surfaces over the state of Colorado and in the Denver metro area.

Potential confounding variables

We identified a total of 42 additional census tract-level factors based on those included as confounders in two previous U.S. county-level analyses of the association between PM2.5 exposure and COVID-19 mortality (Wu et al., 2020; Knittel & Ozaltun, 2020). However, we detected significant problems with multicollinearity among several of these proposed covariates when we translated them to the census tract level in Colorado. Therefore, we selected a reduced set of 20 covariates through stepwise regression and manual examination of the variance inflation factors (VIFs). Descriptions of all 42 potential confounding variables and details of the variable selection procedure are presented in the supplemental information (Tables S1 – S2).

The selected set of covariates includes six measures from the U.S. Census American Community Survey (ACS) 2014–2018 five-year estimates describing the age, race/ethnicity, and income distribution of each census tract. Race/ethnicity was characterized by the proportion of the population that is non-Hispanic African American and the proportion that is non-African American people of color. The latter demographic group is composed mainly of Colorado's Hispanic/Latino community, but also includes the proportion that is American Indian and Alaska Native, Asian, Native Hawaiian and other Pacific islander, or two or more races. We included three measures describing the age distribution in each census tract: the proportion of the population between 18 and 44 years old, the proportion between 45 and 64 years old, and the proportion over 64 years old. We also included the percent of census tract population living in households with family income in the last 12 months below the federal poverty line.

Because COVID-19 transmission is thought to be affected by meteorology (Merow and Urban, 2020; Poirier et al., 2020), we included the 30-year average summer temperature and summer relative humidity at the census tract level as covariates. Although other air pollutants such as ozone also affect respiratory health, these variables were not included as covariates because they are strongly correlated with PM2.5.

We include five measures from the ACS related to social distancing: population density, the proportion of the population living in overcrowded housing, the proportion employed in an essential industry, and the proportion living in correctional, nursing, or mental health facilities. We also considered two measures of social distancing derived from mobile phone data provided by the SafeGraph COVID-19 Data Consortium and compiled for CDPHE by Citizen Software Engineers. The first is a derived measure of virus exposure combining average daily travel to a census tract with the infection rate of the corresponding county of origin. Higher measures correspond to greater volumes of travel to a census tract originating in counties with higher COVID-19 infection rates. The second is an index of the average daily amount of time spent at home per census tract derived from mobile phone location data. For the main analysis, we averaged these variables over the time period spanning March 1, 2020 to the end of the City and County of Denver's stay-at-home order on May 8, 2020.

We also included two measures accounting for underlying health factors in each census tract. These measures were based on estimates of census tract-level adult prevalence of diabetes, heart disease, obesity, and smoking and the age-adjusted rates of hospitalization for asthma, diabetes, heart disease, and influenza. The prevalence estimates were derived from multiple years of Colorado Behavioral Risk Factor Surveillance System data (2014–2017), and the hospitalization data was computed from 5 years (2014–2018) of discharge data from the Colorado Hospital Association. We used principal component analysis to reduce the comorbidity measures to a pair of uncorrelated principal components. The first comorbidity component was strongly correlated with higher prevalence rates while the second comorbidity component was strongly correlated with higher age-adjusted hospitalization rates.

To control for differential testing patterns, we considered the crude testing rate per census tract population from the 2014–2018 ACS estimates. We also considered the time elapsed in days since the first case was reported for each census tract, as well as the number of certified hospital beds per unit population to control for differences due to timing of initial outbreak and overall health care access.

Statistical analysis

We applied a Bayesian hierarchical random effects model, based on a design first proposed by Besag, York and Mollié (BYM) (Besag et al., 1991), to investigate the association between long-term PM2.5 concentrations and counts of COVID-19 outcomes at the census tract level. The BYM model includes both an intrinsic conditional auto-regressive (ICAR) component to account for spatial autocorrelation between census tracts and an exchangeable random effect component to account for uncorrelated variation across all census tracts. The prior distributions controlling the variance of the ICAR and exchangeable components were selected to give approximately equal emphasis to the two components (Bernardinelli et al., 1995; Morris et al., 2019). Our model also includes a county-level exchangeable random effect component to account for uncorrelated variation due to unmeasured county-level factors. Our primary analyses included 12 separate models, one for each combination of three COVID-19 outcomes and four PM2.5 exposure surfaces. In each model, the census tract-level count of each COVID-19 outcome was assumed to follow a Poisson distribution using the indirectly age-sex standardized expected count as an offset term. We generated the expected count using reference rates with a denominator comprising the total number of unique individuals tested for COVID-19 in each census tract. Thus, the expected count has the effect of adjusting for age-sex confounding in the tested population, but it does not account for the age structure of each census tract, hence our inclusion of the age group variables in the fully adjusted models.

We adjusted each of our 12 models by the same set of 20 census tract-level covariates described in the previous section. PM2.5 measurements in this analysis represent long-term annual average census tract concentrations reported in units of μg/m3. Hence, our reported regression coefficients for the PM2.5 variables represent the effect of a 1 μg/m3 increase in long-term annual average PM2.5 concentration at the census tract level. Further, we scaled the two race/ethnicity variables by their respective interquartile ranges (IQR), and thus reported regression coefficients compare a census tract in the 75th-percentile of a given race/ethnicity variable to a census tract in the 25th-percentile. All other variables were converted to z-scores and their coefficients represent the effect of a one-standard-deviation increase above their respective means.

Sensitivity analysis

For each of our 12 models, we conducted four main sensitivity analyses. First, to assess the impact of our choice of offset term, we explored an alternative offset term that used the census tract population per 100,000 from the American Community Survey 2014–2018 5-year estimates. The main distinction between models with different offset terms is that the models using population per 100,000 as the offset contain no inherent adjustment for age-sex confounding in the tested population. Second, to account for the possibility that spatially-correlated random effects may bias the fixed effect estimates (Beale et al., 2010; Hodges & Reich, 2010), we created a separate set of models by replacing the BYM component with a census tract-level random intercept term. Third, to assess the robustness of the results to study area selection, we stratified census tracts into two groups, the first including only census tracts that fall within the 5 most populous counties in the Denver metro area, and the second including all tracts outside of these counties (Figure S3). Finally, to assess the sensitivity of observed associations during different phases of the pandemic, we disaggregated the analysis into 3 time periods. Results of our study area and time period sensitivity analysis are provided in our supplementary materials (Section S3).

Model selection

Considering our main analysis together with our sensitivity analyses of the choice of offset term and method of addressing spatial autocorrelation, our next task was to select the best-fitting model among four candidate study designs. For convenience, we refer to the competing study designs as relative risk-BYM, relative risk-i.i.d., relative rate-BYM, and relative rate-i.i.d. Here, relative risk describes the use of the expected count as the offset term in the regression model whereas relative rate refers to the use of census tract population per 100,000 in the offset term. The i.i.d. model uses a tract-level random intercept term to account for uncorrelated variation across all census tracts, while the BYM models include both a random intercept and an ICAR component at the census tract level. For each study design, we built a series of nested models, starting with a null model that included only the respective offset term and random effects as regressors. We then estimated separate unadjusted risk/rate ratios for the effect of our four different measures of long-term PM2.5 exposure on each COVID-19 outcome. We then estimated fully-adjusted models including all other census tract-level characteristics to control for potential confounding of the PM2.5 effect.

We considered two primary factors to select an optimal model: the deviance information criterion (DIC) and the Moran's I statistic of the model residuals. Briefly, the DIC is a measure of model fit that accounts for both goodness of fit and the effective number of parameters. Lower DIC values indicate better model fit, and differences in DIC greater than 10 are generally considered meaningful (Spiegelhalter et al., 2002). Moran's I is a correlation coefficient that measures the overall spatial autocorrelation of the model residuals. The null hypothesis of the Moran's I test states that the residuals are not spatially correlated, and thus a small p-value yields evidence that significant residual spatial correlation is present.

All analyses were conducted using the R statistical software environment version 3.6.2, and the main analysis was implemented in the R-INLA package version 19.09.03 (R Core Team, 2019; Rue et al., 2009). To protect privacy, COVID-19 outcomes at the census tract level cannot be released publicly; however, code to complete this analysis can be made available upon request.

Results

In this study, we analyzed counts of COVID-19 outcomes from 1237 populated census tracts in the U.S. state of Colorado. Between March 1 and August 31, 2020, the average Colorado census tract had accumulated 42.8 infections (range: 0 to 338), 5.1 hospitalizations (range: 0 to 58), and 1.4 deaths (range: 0 to 40), while 7 (0.56%), 175 (14.15%), and 637 (51.5%) tracts had counts of zero for each outcome, respectively. Meanwhile, pooling the four PM2.5 surfaces together, the average census tract had a long-term annual average PM2.5 exposure of 6.89 μg/m3 with a range between 1.78 and 9.95 μg/m3, depending on the surface source. Unadjusted risk/rate ratios from each model show a positive association between PM2.5 exposure and COVID-19 infection, hospitalization, and mortality, regardless of the provenance of the exposure surface (Fig. 2 ).

Fig. 2

Unadjusted effect of a 1 μg/m3 increase in long-term average PM2.5 concentration on three COVID outcomes, depending on the source of PM2.5 surface and model specification.

After adjusting for potential confounding variables, the analyses using the EPA, CSU, or Schwartz surfaces typically show positive associations between PM2.5 exposure and COVID-19 infection, hospitalization, and mortality, while the analysis using the ACAG surface shows a negative association between PM2.5 exposure and all 3 COVID-19 outcomes (Fig. 3 ).

Fig. 3

Estimated effect of a 1 μg/m3 increase in long-term average PM2.5 concentration on three COVID outcomes, from fully adjusted models.

Table 2 displays DIC value and Moran's I statistics from fully-adjusted models for each of our four study designs using the EPA PM2.5 surface as the key exposure variable. DIC values and Moran's I test results both indicated the relative risk-BYM design produced the best-fitting model. Given the improved fit from the relative risk-BYM model design and the robust results across all four designs, we focus primarily on this model for the remainder of this paper.

Table 2

Deviance Information Criteria (DIC) and Residual Moran's I Statistic (p-value) by study design and outcome.

Study Design	Infections	Hospitalizations	Mortality
Relative risk - BYM	8542−0.05 (p > 0.99)	5178−0.04 (p > 0.99)	3043−0.06 (p > 0.99)
Relative rate - BYM	8655−0.07 (p > 0.99)	5210−0.04 (p > 0.99)	3064−0.05 (p > 0.99)
Relative risk - i.i.d.	85720.09 (p < 0.001)	51830.02 (p = 0.13)	30510.02 (p = 0.14)
Relative rate - i.i.d.	86700.04 (p < 0.001)	52180.02 (p = 0.09)	30670.01 (p = 0.18)

PM2.5 effects

Our results show that the choice of exposure surface has a large impact on the direction and significance of the adjusted association between long-term PM2.5 exposure and COVID-19 outcomes. The EPA exposure surface displayed the most consistent positive adjusted associations with COVID-19 outcomes. After controlling for all other census tract-level factors, long-term PM2.5 exposure measured by the EPA surface was statistically associated with an increase in COVID-19 infections in two out of four study designs, and with increases in hospitalizations and mortality across all four study designs. According to our best-fitting model, a 1 μg/m3 increase in long-term PM2.5 exposure measured by the EPA surface was associated with a positive, but insignificant increase in the relative risk of COVID-19 infections (RR: 1.10, 95% CI: 0.98–1.24). The same increase in long-term PM2.5 exposure was associated with a 26% increase in the relative risk of hospitalizations (RR: 1.26, 95% CI: 1.06–1.48) and a 34% increase in the relative risk of mortality (RR: 1.34, 95% CI: 1.02–1.77). Results for the CSU and Schwartz surfaces were generally consistent with those of the EPA surface; however, none of the associations for these two surfaces were statistically significant in the relative risk-BYM model. At the same time, the results from the ACAG surface point in the opposite direction from the three other exposure surfaces. In our best-fitting model, a 1 μg/m3 increase in long-term PM2.5 measured by the ACAG surface was associated with a statistically insignificant decrease in the relative risk of COVID-19 infections (RR: 0.99, 95% CI: 0.94–1.04), hospitalizations (RR: 0.98, 95% CI: 0.91–1.06), and mortality (RR: 0.91, 95% CI: 0.78–1.07).

Other effects

While not the primary goal of our analysis, we identified significant associations between several other factors and COVID-19 infections, hospitalizations, and mortality at the census tract level in Colorado. Adjusting for all other census tract-level factors, we found statistically significant evidence that census tracts with a larger proportion of non-Hispanic African American population are associated with increased risk of COVID-19 infections (RR: 1.04, 95% CI: 1.02–1.06) and hospitalizations (RR: 1.07, 95% CI: 1.04–1.11). In addition, census tracts with larger proportions of non-African American people of color were positively and significantly associated with the risk of infections (RR: 1.31, 95% CI: 1.25–1.39), hospitalizations (RR: 1.43, 95% CI: 1.31–1.57), and mortality (RR: 1.56, 95% CI: 1.30–1.87). With respect to social distancing-related factors, we found a one standard deviation increase in our mobility risk index was statistically associated with increases in the risk of infections (RR: 1.18, 95% CI: 1.09–1.27) and hospitalizations (RR: 1.14, 95% CI: 1.04–1.26).

In addition, we found evidence that communities with a larger proportion of essential workers were subject to increased risk of COVID-19 infections (RR: 1.05, 95% CI: 1.02–1.09) and hospitalizations (RR: 1.06, 95% CI: 1.00–1.13), whereas those with greater percentages of populations confined to correctional, nursing, or mental health facilities had lower risk of hospitalizations (RR: 0.92, 95% CI: 0.87–0.98), but higher risk of mortality (RR: 1.12, 95% CI:1.01–1.22). We also found statistically significant evidence that mortality risk is higher in both census tracts with greater population percentages of adults age 65 years and older (RR: 1.30, 95% CI: 1.14–1.48) and tracts with greater population percentages of adults age 18–44 years (RR: 1.36, 95% CI: 1.12–1.67). Lastly, the number of days since the first case in a census tract was diagnosed was positively and significantly associated with higher risk of infections (RR: 1.08, 95% CI: 1.04–1.11), hospitalizations (RR: 1.27, 95% CI: 1.16–1.40), and mortality (RR:1.32, 95% CI: 1.11–1.62), while the crude COVID-19 test rate was negatively and significantly associated with the risk of infection (RR: 0.96, 95% CI: 0.94–0.98) and hospitalization (RR: 0.91, 95% CI: 0.88–0.95). The relative risks and 95% CIs for all fixed effects in the fully adjusted relative risk-BYM model are shown in Fig. 4 , for all 3 COVID-19 outcomes and all four PM2.5 surfaces.

Fig. 4

Effect estimates for all covariates included in the model. Points represent the estimated mean and thin lines represent the 95% confidence interval. Shapes indicate the COVID-19 outcome considered and colors the PM2.5 surface used. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)

Maps of the spatially random effects included in the model are shown in Figures S4 – S6. All three spatial random effects are highly consistent between the different PM2.5 surfaces, suggesting that the random effects are accounting for unmeasured differences between counties, rather than responding to differences in the PM2.5 surfaces.

Discussion

Long-term PM2.5 exposure

Previous studies on the effects of long-term PM2.5 exposure on COVID-19 outcomes have reported a wide range of results. At least six studies have reported estimated effects for a 1 μg/m3 increase in PM2.5 concentrations on COVID-19 mortality at the county level in the United States (Knittel and Ozaltun, 2020; Wu et al., 2020; Tian et al., 2020; Liang et al., 2020; Rodriguez-Diaz et al., 2020; Millett et al., 2020). Estimates from these studies range from a non-significant 1% decrease (Knittel and Ozaltun, 2020) to a statistically significant 17% increase (Tian et al., 2020). The Liang et al., 2020 study used the Schwartz surface for their PM2.5 estimates, while the other five studies used the ACAG surface. In the United Kingdom, Travaglio et al., 2021 found that a 1 μg/m3 increase in PM2.5 can be associated with anywhere from a 4% decrease to a non-significant 1% increase, while Konstantinoudis et al. (2021) found a non-significant 1.4% increase in COVID-19 mortality risk. In the Netherlands, Cole et al., 2020 found a 1 μg/m3 increase in long-term PM2.5 is associated with a 14% increase in COVID-19 mortality.

Studies looking at COVID-19 infections have reported even greater variability in the estimated PM2.5 effects. At the county level in the United States, Millett et al., 2020 and Rodriguez-Diaz et al. (2020) both report that a 1 μg/m3 increase in PM2.5 is associated with a non-significant 1% increase in COVID-19 infections. In the United Kingdom, Travaglio et al., 2021 find that estimated effect can range from a 10% decrease to a 12% increase, depending on how the data are aggregated; in the Netherlands, Andrée, 2020 found that an increase in PM2.5 from 10 to 12 μg/m3 is associated with nearly a 100% increase in COVID-19 infections; in China, Zheng et al., 2021 report that a 1 μg/m3 increase in PM2.5 is associated with a statistically significant 3.2% increase in COVID-19 infections.

These studies used different model designs, covered different time periods, and adjusted for different covariates. It is therefore unsurprising that these studies found different estimates of the effect on PM2.5 on COVID-19 outcomes. Our analysis adds another potential source of uncertainty, that differences in how PM2.5 is estimated over small spatial regions can cause noticeable changes in the results. Among the 12 surface-outcome combinations studied here, results from our best-fitting model suggested 2 associations are positive and statistically significant, 7 are positive but not statistically significant, and 3 are negative but not statistically significant. None of our results from our best-fitting model in the statewide analysis or in the time period or study area sensitivity studies show a statistically significant decrease in COVID-19 outcomes with increases in long-term PM2.5 concentrations.

At the moment, while we can compare the performance of the four PM2.5 surfaces against EPA monitoring data (Section S1), we have no independent data source to compare the PM2.5 surfaces against. Without the ability to independently measure the error in the PM2.5 surfaces, it is impossible to identify which surfaces are most accurate in Colorado. Looking forward, better estimates of the health effects of environmental exposures will require resolving the disagreements among the four PM2.5 estimation surfaces.

Figure S7 shows the change in our PM2.5 effect estimate on hospitalizations when the PM2.5 estimate for a single census tract is changed from the EPA surface to the ACAG surface. We used the EPA and ACAG surfaces in this sensitivity analysis as these surfaces showed the highest and lowest PM2.5 effects. The largest changes in PM2.5 coefficients are generally found when the PM2.5 surfaces are changed in census tracts running along the eastern edge of the Front Range Urban Corridor. This region combined high disagreement between PM2.5 surfaces, high numbers of COVID-19 hospitalizations, and relatively high population. The high leverage of this region suggests that improving the accuracy of PM2.5 estimates in this area could significantly improve our ability to link environmental exposures to health outcomes in Colorado.

High risk populations

Early descriptions of outcomes among patients with COVID-19 in the United States indicated disease severity is highest among adults age 65 years and older, people with underlying health conditions, and is strongly linked to race/ethnicity and other socioeconomic factors (CDC COVID-19 Response Team, 2020a, 2020b; Hsu et al., 2020). Health inequities of the COVID-19 pandemic are thought to be at least partially explained by well-documented racial and ethnic disparities in the prevalence of chronic health conditions (Alcendor, 2020); however, several recent studies have suggested disparities among Black and Hispanic/Latino communities may be better explained by overrepresentation among the essential workforce (Oppel et al., 2020; Rogers et al., 2020). While our ecological study designs means we cannot make inferences about individuals, our results offer a proxy for understanding how these patterns increase vulnerability disproportionately for certain populations in Colorado.

Our findings suggest communities of color in Colorado are subject to higher risk of infection as well as of more severe complications of COVID-19. Using the fully-adjusted relative risk-BYM model we found statistically significant associations between larger proportions of non-Hispanic African Americans and increased risk of COVID-19 infections and hospitalizations regardless of the choice of PM2.5 surface, while the association with the risk of COVID-19 mortality was positive but not statistically significant. We also found that larger proportions of non-African American people of color, a group largely composed of Colorado's Hispanic/Latino community, were statistically associated with increased risk of all three COVID-19 outcomes. Further, we found a separate statistically significant association between the percentage of census tract population employed in an essential industry and the risk of COVID-19 infection, which appears to be further supported by our finding that communities with more individuals in prime working age (18–44 years) are subject to greater risk of mortality.

Our findings are consistent with Millett et al., 2020 who in a similar ecological study of U.S. counties found that higher proportions of non-Hispanic African American residents were associated with increased rates of infection and mortality. With respect to non-Hispanic African Americans, the similarity of our findings to those of a nationwide study are particularly noteworthy given Colorado's smaller overall population of non-Hispanic African Americans compared to other regions of the United States. A per-IQR increase in the proportion of non-Hispanic African Americans among Colorado census tracts represents an increase of only 3.3% (0.3%, 3.6%), while an identical increase among all U.S. counties represents a 9.6% increase (0.7%, 10.3%). By contrast, a similar ecological study of U.S. counties by Rodriguez-Diaz et al. (2020) did not find statistically significant associations between predominantly Hispanic/Latino counties and COVID-19 rates in the Western U.S. Our findings suggest sub-county variation in the population of people of color may be crucial to understanding health disparities of the COVID-19 pandemic in Colorado. Together, the results of our study contribute additional support for the need to address the disproportionate impact of COVID-19 on higher-risk populations as the pandemic response and recovery efforts continue in Colorado.

Conclusion

This study quantifies the environmental, demographic, and socioeconomic associations and disparities for three COVID-19 outcomes in Colorado. We show the uncertainty in census tract PM2.5 concentrations is a key factor limiting our ability to link environmental exposures and health outcomes, and highlight specific regions in Colorado where the uncertainty in PM2.5 concentrations is most impactful.

In the majority of our analyses, we find that increases in long-term exposure to fine particle pollution is associated with a positive but non-statistically significant increase in infections, hospitalizations, and deaths due to COVID-19. It is notable that we see these results even at relatively low concentrations of PM2.5. The average concentration of PM2.5 in Colorado ranges from 1.78 to 9.95 μg/m3, below the World Health Organization's guideline value for PM2.5 and meeting the EPA's annual National Ambient Air Quality Standard. Our findings are consistent with studies showing a consistent effect of PM2.5 on all-cause mortality with no evidence of a threshold (Liu et al., 2019).

However, we are limited in the conclusions we can draw by the uncertainty in PM2.5 concentrations at the census tract level. Researchers examining the links between PM2.5 and COVID-19 should carefully consider how uncertainty in PM2.5 concentrations affects their results. This is especially true when attempting to understand risk factors in Colorado over small spatial scales or in neighborhoods with dramatically different disease rates.

In the long-term, improved understanding of the links between air pollution and human health at the neighborhood scale will require the development of improved PM2.5 surfaces that can more accurately capture variations in PM2.5 concentrations on the scale of kilometers. There is also a clear need for a method to independently validate these PM2.5 surfaces after they are generated.

While the uncertainty in PM2.5 concentrations makes it difficult to establish definitive conclusions, this study provides evidence that local rates of COVID-19 infections, hospitalizations, and mortality are influenced by patterns of exposure to air pollution, racial and ethnic composition, local travel patterns, and risk factors such as aging population structure and underlying comorbidities. As such, this study could help provide a scientific basis for precisely targeting COVID-19 response and recovery efforts in Colorado based upon community-specific risk factors. For example, additional public health interventions, such as recommending ways to decrease personal exposure to fine particle pollution and considering ways to effect more social distancing, could be implemented in specific neighborhoods to better protect those most at risk.

Finally, the results of this analysis indicate the need for expanded investigations regarding the link between air pollution exposure and COVID-19 outcomes. In particular, further analysis of the link between short-term air pollution exposure and COVID-19 is urgently needed to establish whether these more modifiable exposures play a role in infection and disease severity.

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.