Spatial analysis of COVID-19 and traffic-related air pollution in Los Angeles.

Introduction

As of February 23rd, 2021, more than 113 million people worldwide have been diagnosed with COVID-19, resulting in more than 2.5 million deaths (World Health Organization, 2021). Extensive investigation has been conducted on the etiology of COVID-19, yet researchers are still determining how exposure risk factors may influence COVID-19 incidence and mortality. Recent evidence from China, Italy, England, and the United States suggest that exposure to air pollution may play a role in COVID-19 incidence and deaths (Brandt et al., 2020, Coker et al., 2020, Li et al., 2020, Lippi et al., 2020, Travaglio et al., 2021, Wang et al., 2020a, Wu et al., 2020b, Zhang et al., 2020, Zhu et al., 2020). These findings are consistent with prior research suggesting that air pollution, including traffic-related air pollution (TRAP), is associated with many respiratory morbidities (e.g., asthma, chronic pulmonary disease, lung cancer, and respiratory tract infections) (Bai et al., 2018, Dales et al., 2008, Franklin et al., 2015, Jerrett et al., 2008, Sydbom et al., 2001), hospitalizations (Neupane et al., 2010), all-cause mortality (Beelen et al., 2008, Jerrett et al., 2005) and increased risk of respiratory viral infection (Ciencewicki and Jaspers, 2007, Wang et al., 2020b). Nitrogen dioxide (NO2), a tracer of TRAP generated from tailpipe emissions (Quiros et al., 2013, Zeldovich, 2015), has been found to impair the function of alveolar macrophages and epithelial cells, thereby increasing the risk of lung infections (Neupane et al., 2010).

Other factors such as age, race/ethnicity, and other sociodemographic characteristics appear to increase risk for COVID-19 infection, severity, and associated death (Brandt et al., 2020). For example, compared to non-Hispanic whites, cumulative COVID-19 hospitalization rates for Black and Latinx populations are approximately 4.7 and 4.6 times higher in the U.S., respectively (Centers for Disease Control and Prevention, 2020). Black and Latinx U.S. populations are disproportionately exposed to SARS-CoV-2, as they are more likely to serve as essential workers (Martinez et al., 2020, Rogers et al., 2020) and to live in crowded conditions (Burr et al., 2010, Memken and Canabal, 1994). A higher prevalence of metabolic disorders (such as hypertension, diabetes, and obesity) in these populations likely contributes to more severe disease (Commodore-Mensah et al., 2018, Divens and Chatmon, 2019) and death from COVID-19 (Du et al., 2020). In addition, minority populations are more likely to live in areas where there is greater air pollution (Ailshire and García, 2018, Collaco et al., 2020, Gaither et al., 2019).

Building on past research demonstrating an association between Severe Acute Respiratory Syndrome (SARS) and air pollution (Cui et al., 2003), Wu et al. (2020a) reported associations between county-level COVID-19 mortality rates in 3089 counties through June 2020 and long-term average (from years 2000 to 2016) PM2.5 concentration across the United States. They reported that each 1 µg/m−3 increase in PM2.5 concentration was associated with an 11% increase in COVID-19 mortality rate (95% CI: 6%, 17%). In another recent study of COVID-19 in 3122 U.S. counties through July 2020, researchers found an increase in interquartile range (IQR) of 4.6 ppb of NO2 to be associated with a 16.2% (95% CI: 8.7%, 24.0%) increase in mortality rate and an 11.3% (95% CI: 4.9%, 18.2%) increase in case-fatality rate (Liang et al., 2020). In both studies, exposure estimates were based on concentrations at the county level and, therefore, could not account for variation in air pollution observable over smaller areas within cities (Kulhánová et al., 2018, Wu et al., 2019). In addition, the quality and comparability of COVID-19 health outcome information may vary considerably across U.S. counties as reporting protocols may differ among jurisdictions at the local, state, and national levels (Bergman et al., 2020, Bialek et al., 2020), which may lead to case ascertainment bias.

Los Angeles became one of the only metropolitan cities globally to publicly report neighborhood-level COVID-19 cases in March 2020 and mortality in June 2020 (LACDPH, 2020). These data afforded the opportunity to conduct spatial modeling for a large population with a smaller geographical area neighborhood unit of analysis. These smaller geographic areas allow for more accurate pollution exposure estimates than the county-level studies above. Los Angeles has a wide range of air pollution exposure levels with which to investigate intra-urban relationships with COVID-19. Furthermore, because the Los Angeles County Department of Public Health (LACDPH) governs all health statistics, Los Angeles County is likely to have consistent health reporting practices. This diminishes the possibility of case ascertainment bias that may have been present in the national studies comparing among more than 3000 counties.

Here we aim to analyze the relationship of air pollution and COVID-19 case incidence, mortality, and case-fatality rates in neighborhoods of Los Angeles County, using high-resolution exposure models. We focus on NO2 because this gaseous pollutant serves as a marker for traffic pollution, which displays substantial intra-urban variation over small areas in Los Angeles and elsewhere (Su et al., 2009, Su et al., 2020, Zeldovich, 2015). In California, 62% of NOx emissions come from mobile sources such as vehicle traffic (Almaraz et al., 2018).

Material and methods

Setting

This study is situated in Los Angeles (LA) County. In 2019, LA County had a population size of 10,039,107 and was diverse in its racial, ethnic, and socioeconomic composition (U.S. Census Bureau, 2020). For example, 51% are White, 48% are Latinx, 15% are Asian, and 8.3% are Black. The median income of LA County is $64,251 USD. LA County is spread across a large geographic area of 4057 square miles or 10,508 square kilometers. The sprawling landscape induces high levels of travel by automobile and TRAP (Su et al., 2009). In addition, the presence of two major seaports and associated goods movement infrastructure creates additional emissions from diesel vehicles (Kozawa et al., 2009, Su et al., 2016, Su et al., 2020). The first case of COVID-19 in California was identified on January 26th, 2020 (LACDPH, 2020, Los Angeles Times, 2020), and the first community-acquired case in the United States was confirmed in California on February 26th, 2020 (CADPH, 2020, Heinzerling et al., 2020).

Data sources

Table 1 summarizes the data sources and variables used. Cumulative COVID-19 case and mortality counts for March 16th to February 23rd, 2021 were accessed from the Los Angeles County Department of Public Health (LACDPH) COVID-19 dashboard website. These outcome data were split into two time periods: a main study period from March 16th to September 8th, 2020; and a secondary period for sensitivity analyses from September 8th, 2020 to February 23rd, 2021. These data are reported at a neighborhood statistical area unit geography. LACDPH reports infectious disease data for ‘Countywide Statistical Areas (CSAs)’, used by many LA County agencies to report data to the County Board of Supervisors, which include mixed areal classifications such as ‘city’, ‘community’, ‘neighborhood’, or ‘unincorporated area’ (Harris, 2020, LACDPH, 2021) This study refers to CSAs as ‘neighborhoods’. Based on prior research (Su et al., 2020), a land-use regression model was used to produce an annual pollution surface of NO2 across California at a spatial resolution of 30 m using data from 2016. This surface was used previously for another recent health study in Los Angeles (Wing et al., 2020). The land-use regression model had an out-of-sample cross-validation R2 of 0.76. This average annual NO2 surface was used to define the main exposure metric by neighborhood.

Table 1

Data sources and spatiotemporal dimensions for model of association between NO2 and COVID-19 case, mortality, and case-fatality rates in Los Angeles County.

Data source	Attribute(s)	Spatial Dimension	Temporal Dimensiont
LACDPH1	COVID-19 cases and COVID-19 deathsPopulation*, smoking, and obesity	Neighborhood statistical areas (polygon)	March 16th – September 8th, 2020September 8th, 2020 –February 23rd, 2021α2019
ACS2	Age, median income, race/ethnicity, owner-occupancy status, and population*	Census tracts (polygon)	2018
500 Cities Project	Diabetes and hypertension	Census tracts (polygon)	2019
LA County GeoHub	Testing locations and hospital locationsBuilding footprints**	Site location (point or polygon)	September 8th, 20202014
LUR3surface (Su et al., 2020)	NO2 (ppb)	California (raster; 30 m)	2016

*LACDPH population used for regression modeling; ACS population used to validate aggregation methods.

**building footprints used for data aggregation.

LACDPH = Los Angeles County Department of Public Health.

ACS = American Community Survey (U.S. Census Bureau).

LUR = land-use regression.

Temporal dimension describes the range in time for which the data was recorded.

Secondary time period used as sensitivity analysis in comparison to main study period of March 16th to September 8th, 2020.

Potential covariates were identified a priori based on existing literature on risk factors for disease or severity of disease (including death) for COVID-19 (Myers et al., 2020), other pneumonic infectious diseases (Neupane et al., 2010), and previous studies on air pollution and COVID-19 (Liang et al., 2020, Wang et al., 2020b). Demographic covariates, including age, race/ethnicity, median household income, and household owner occupancy, were downloaded from the U.S. Census Bureau’s American Community Survey (ACS) 5-year moving estimate for 2018 (U.S. Census Bureau, 2018). Population counts at the neighborhood level, and smoking and obesity prevalence at the census-tract level were downloaded from the LACDPH website (LACDPH, 2018). Population counts from LACDPH were compared to counts from the ACS to assess variable aggregation methods. Hospital and testing facility locations were acquired from the LA County’s GeoHub website for their potential association with case ascertainment (LAC, 2014). Residential building footprints were also downloaded from the LA County Geohub website to facilitate these demographic data aggregations (LAC, 2014).

We also considered hypertension and diabetes as health outcomes potentially associated with COVID-19 severity (Myers et al., 2020). These health outcomes were modeled in the U.S. Centers for Disease Control and Prevention’s (CDC) 500 Cities Project health dataset (Centers for Disease Control and Prevention, 2019). These covariates were included as sensitivity analyses due to incomplete spatial coverage,

Uncertainties in the testing regime raise questions about potential case-ascertainment bias. As testing became more widely available, rates of testing likely changed from testing only suspected cases to people potentially exposed as a result of occupational or social interactions. This could have affected the case rate and subsequently the case-fatality rates. It is also possible that in the earlier stages of the pandemic, there was more undercounting of the deaths, which would have diminished over time as medical professionals learned how to identify more accurately deaths resulting from COVID-19.

To address the potential shifts in case-ascertainment, case-fatality, and mortality rates that could have occurred over time, we conducted further sensitivity analyses. Specifically, we extended our original study period which captured approximately the first six months of the pandemic (March 16th to September 8th, 2020), to the subsequent six months (September 8th, 2020 to February 23rd, 2021). Thus, we replicated the analysis for the subsequent 6-month period, which had nearly four times the incident cases (875,368 cases) as the first period (230,621 cases). In the latter period, the County changed the neighborhood definitions to exclude or combine about 13 neighborhoods, so the count of neighborhoods was less than in the original period (348 vs. 335). Consequently, the two data sets are not uniformly constructed, but they are quite similar.

Quantification of variables

Very few spatial variables were available at the neighborhood geographies, as most environmental, health, and demographic areal data are published by postal ZIP code or census tract. Environmental Systems Research Institute’s (ESRI) ArcGIS 10.7 (ESRI, 2020) was used to summarize NO2 zonal mean by each of the N = 348 neighborhoods (‘neighborhood’ statistical area geographies, as delineated by the LACDPH) in LA County. To account for misalignment in areal boundaries between COVID-19 case/mortality and selected covariates, all areal covariates were first reaggregated to residential building footprints (acquired from LA County’s Geohub website (LAC, 2014)) and then reaggregated to neighborhoods by using counts per area density-based raster surfaces. This intermediate step was taken to minimize the effect of geographies with highly variable population densities (e.g., a large neighborhood with few total residents). Hospital and testing facility areal densities were calculated using a 10-mile radius kernel density process – generating a raster surface describing the number of hospitals or facilities per sq km within 10 miles of each raster grid cell. History of hypertension and diabetes drawn from the 500 Cities Project covered only 61% (212 of 348) of neighborhoods; therefore, we imputed the global mean for the remaining 136 neighborhoods and report these results as sensitivity analyses.

Statistical modeling

Extracted neighborhood NO2 concentrations were modeled in relation to incident case rate (cases/population), mortality rate (deaths/population), and case-fatality rate (deaths/cases). NO2 concentrations from neighborhoods were scaled to the interquartile range to aid in interpretation of model results (Liang et al., 2020, Wu et al., 2020a). We used three different statistical models to assess sensitivity of our results to model specification: (1) zero-inflated Poisson, (2) zero-inflated negative binomial, and (3) Bayesian conditional autoregressive (CAR) zero-inflated Poisson models. To assess incident case rate, all models treated the count of COVID-19 cases in the neighborhoods as the dependent variable and the total population as the offset. These models were also run for mortality counts with the total population as the offset (mortality rate) and for mortality counts with the total number of COVID-19 cases as the offset (case-fatality rate). Zero-inflated models were selected to account for a high number of neighborhoods with zero counts of cases or deaths and often low total populations; without zero-inflation, these low-count areas could disproportionately influence the model results. We also employed a Bayesian zero-inflated Poisson model whose spatial random effects were assigned a CAR prior distribution to account for potential spatial autocorrelation between neighborhoods. This model incorporated a spatial adjacency matrix of first-order neighbors and employed flat priors.

All three models were run with and without adjustment for covariates. The final model included the following covariates: mean percent owner occupancy, mean median income, mean percent above 65 years old, mean percent nonwhite; mean smoking prevalence, mean obesity prevalence, and mean hospital density per square mile within a 10-mile radius. Covariates selected for final model were identified a priori; however, in the event of highly correlated or colinear covariate pairs, the covariate with the highest bivariate association with the outcome was included. In all models, no covariates were found to be significant predictors of the zero-inflation component. We used R version 3.6.3 to run all statistical analyses (R Core Team, 2020).

Results

The average area of the 348 LA County neighborhoods was 44.7 sq km (SD = 171.8 sq km), with the largest being 1144 sq km (Antelope Valley) in the northern exurban areas of the county and the smallest being 0.67 sq km (San Pasqual) in a more densely populated area near Pasadena, north of downtown LA. The annual mean NO2 across the study region was 11.7 ppb (SD = 7.3 ppb; range of 1.6–31.3 ppb). Concentrations of NO2 derived from the 2016 land-use regression surface are depicted in Fig. 1 A. The mean aggregated NO2 across neighborhoods was 15.6 ppb (SD = 6.0 ppb) with an interquartile range of 8.7 ppb. Between March 16th and September 8th, 2020, the LACDPH recorded 230,621 confirmed cases of COVID-19, and 5653 deaths due to COVID-19 were observed. In a population of about 10 million, this translated into a case rate of 2.2% (Fig. 1B), a mortality rate of 0.054% (Fig. 1C), and a case-fatality rate of 2.5% (Fig. 1D). The period from September 8th, 2020 to February 23rd, 2021 included 875,368 cases and 13,344 deaths and was used as sensitivity analysis.

Fig. 1

Maps of Los Angeles County: A. NO2 from land-use regression (LUR) model, 2016. (Map is zoomed in to demonstrate fine-resolution variability); B. COVID-19 Case rate (cases/population); C. Mortality rate (deaths/population); and D. Case-fatality rate (deaths/cases) for the period between March 16th and September 8th, 2020, depicted at the neighborhood level.

Between neighborhoods (N = 348), the mean percent owner occupancy was 54.7% (SD = 22.0%); mean median income was $47,483 (SD = $68,898); mean percent above 65 years old was 13.9% (SD = 0.1%); mean percent nonwhite was 45.2% (SD = 20.1%); mean smoking prevalence was 12.7% (SD = 2.6%); mean obesity prevalence was 23.7% (SD = 7.4%); and the mean hospital density per square mile within a 10-mile radius was 1.2x10-4 (SD = 1.4x10-4). The mean hypertension and diabetes prevalence were 18.7% (SD = 9.6%), and 7.4% (SD = 4.2%), respectively. Although median income was highly associated with the outcome in the crude model, it was highly correlated with owner occupancy (r = 0.89); therefore, owner occupancy was selected for inclusion as it demonstrated a larger bivariate association with all outcomes.

Ordinary residuals from the zero-inflated Poisson and negative binomial models demonstrated significant spatial autocorrelation at the global (Moran’s I p-value < 0.001) and local (Anselin hot-spots) levels; thus, we used a Bayesian zero-inflated Poisson model with a CAR prior on the random effects to account for spatial dependence of the residuals per neighborhood.

Crude and adjusted model results for (1) zero-inflated Poisson, (2) zero-inflated negative binomial, and (3) zero-inflated Poisson spatial models are shown in Table 2 . In the adjusted zero-inflated Poisson model, we found that the incidence rate ratio (IRR) of NO2 was 1.31 (95% CI: 1.29, 1.33) for the case rate. That is, we found that an increase of 8.7 ppb (IQR) in mean annual NO2 (2016) was associated with a 31% increase in COVID-19 incident case rate. The adjusted zero-inflated negative binomial and spatial models demonstrated a smaller effect of 16% (95% CI: 2%, 32%) and 18% (Credible Interval - CrI: 10%, 32%) increase in case rate, respectively (Table 2A). Adjusted models reduced residual uncertainty compared to the crude model estimates for case rate. The adjusted Poisson, negative binomial, and spatial models all demonstrated an increase in COVID-19 mortality of 35% (95% CI: 23%, 48%), 44% (95% CI: 11%, 86%), and 60% (CrI: 37%, 88%), respectively, across the IQR exposure increment (Table 2B). Again, these adjusted models reduced residual uncertainty compared to crude model estimates for mortality rate. Finally, adjusted Poisson and negative binomial models showed positive yet non-significant results for the association between NO2 and COVID-19 case-fatality; however, the spatial CAR model demonstrated that an IQR increase in mean annual NO2 was associated with a 31% (CrI: 10%, 65%) increase in case-fatality. Sensitivity analyses, which included the addition of history of hypertension and diabetes in the models, had comparable results across all three models. In comparing the two time periods, before and after September 8th, 2020, we found that the results were largely consistent, despite very different case numbers, testing regimes, and improvements in classifying deaths. While some differences exist in the size of the effects, overall the conclusions remain the same. The results of this sensitivity analysis give some assurance that the changes in testing, case ascertainment, and mortality classification over time are not having a substantial effect on the key conclusion that long-term air pollution exposure likely increases the risk of Covid-19 infection and death. Full model results, including incidence risk ratios for all covariates, are described in Appendix A for the main model, Appendix B for sensitivity analysis including hypertension and diabetes, and Appendix C for sensitivity analysis utilizing the period between September 8th, 2020 and February 23rd, 2021.

Table 2

Adjusted association of NO2 (scaled by interquartile range: 8.7 ppb) and COVID-19 from three models: A. Case rate (cases/population); B. Mortality rate (deaths/population); and C. Case-fatality rate (deaths/cases) for Los Angeles County neighborhoods (N = 348) for the period between March 16th and September 8th, 2020. A sensitivity analysis was conducted to assess the inclusion of hypertension and diabetes as model covariates. A second sensitivity analysis assessed COVID-19 outcome data from a second time period between September 8th, 2020 and February 23rd, 2021 for N = 335 neighborhoods.

	Zero-inflated Poisson				Zero-inflated negative binomial				CAR zero-inflated Poisson with spatial random effect
	Crude		Adjusted		Crude		Adjusted		Crude		Adjusted
Main model1(n = 348)	IRR	CI	IRR	CI	IRR	CI	IRR	CI	IRR	CrI	IRR	CrI
A. Case rate	1.82	(1.80, 1.84)	1.31	(1.29, 1.33)	1.47	(1.33, 1.62)	1.16	(1.02, 1.32)	1.77	(1.53, 2.09)	1.18	(1.10, 1.32)
B. Mortality rate	1.72	(1.62, 1.83)	1.35	(1.23, 1.48)	1.77	(1.50, 2.08)	1.44	(1.11, 1.86)	1.94	(1.46, 2.58)	1.60	(1.37, 1.88)
C. Case-fatality rate	0.96	(0.91, 1.01)	1.05	(0.96, 1.15)	1.07	(0.91, 1.25)	1.21	(0.97, 1.50)	1.13	(0.87, 1.42)	1.31	(1.10, 1.65)
Sensitivity Analysis: Including hypertension & diabetes2(n = 348)
A. Case rate	1.82	(1.80, 1.84)	1.28	(1.26, 1.30)	1.47	(1.33, 1.62)	1.18	(1.04, 1.33)	1.77	(1.53, 2.09)	1.27	(1.14, 1.34)
B. Mortality rate	1.72	(1.62, 1.83)	1.35	(1.23, 1.49)	1.77	(1.50, 2.08)	1.57	(1.23, 2.01)	1.94	(1.46, 2.58)	1.44	(1.13, 2.06)
C. Case-fatality rate	0.96	(0.91, 1.01)	1.05	(0.96, 1.15)	1.07	(0.91, 1.25)	1.19	(0.96, 1.49)	1.13	(0.87, 1.42)	1.34	(1.13, 1.69)
Sensitivity Analysis: Second time period3(n = 335)
A. Case rate	1.66	(1.65, 1.67)	1.29	(1.28, 1.29)	1.22	(1.04, 1.42)	1.21	(1.01, 1.45)	1.38	(1.51, 1.66)	1.24	(1.10, 1.42)
B. Mortality rate	1.72	(1.63, 1.81)	1.28	(1.21, 1.36)	1.82	(1.53, 2.18)	1.38	(1.14, 1.67)	1.88	(1.53, 2.35)	1.77	(1.34, 2.14)
C. Case-fatality rate	1.04	(0.99, 1.09)	1.04	(0.98, 1.10)	1.10	(0.95, 1.26)	1.09	(0.94, 1.25)	1.10	(0.98, 1.24)	1.15	(1.00, 1.35)

All models controlled for owner-occupancy rate, percent population > 65 years of age, percent nonwhite, percent smokers, percent obese, and hospital density per sq mi within 10 miles. Full model results with covariates can be seen in Appendix A. COVID-19 outcome data was acquired from the Los Angeles County Department of Public Health (LACDPH).

Model conducted for sensitivity analysis includes all covariates from main model plus percent hypertensive and percent diabetic derived from 500 Cities Project data. Hypertension and diabetes data covered 212 of N = 348 neighborhoods (61% coverage). We imputed the global mean for the remaining 136 neighborhoods. Full model results from sensitivity analysis can be seen in Appendix B.

Model conducted for sensitivity analyses includes all covariates from main model. COVID-19 outcome data acquired from the Los Angeles County Department of Public Health (LACDPH).

Discussion and conclusion

We found annual NO2 to be associated with COVID-19 incidence and mortality in Los Angeles County neighborhoods while adjusting for numerous confounders. These findings were consistent across statistical model specification, although risk estimates displayed some variation between models. In addition, we found in the CAR an association between NO2 and COVID-19 case-fatality; other models also showed a positive but insignificant associations. Covariates in the models largely had the expected sign of effect. Furthermore, our sensitivity analyses, which included the addition of hypertension and diabetes prevalence covariates, had a minimal impact on the effect size or interpretation of our model estimates for NO2 and COVID-19 outcomes. Our sensitivity analysis involving a second, approximately 6-month time period, with slightly different outcome reporting (N = 348 vs N = 335) also demonstrated comparable results to the main model and study period.

Our findings are consistent with two previous studies demonstrating a relationship between air pollution and COVID-19 nationally at the county scale in the U.S. One study investigated the association between NO2 and COVID-19, and they observed remarkably similar findings. Specifically, Liang et al. (2020) reported an increase of 4.6 ppb (IQR across all counties) NO2 to be associated with a 16.2% (CI: 8.7%, 24.0%) increase in mortality rate and an 11.3% (CI: 4.9%, 18.2%) increase in case-fatality rate (Liang et al., 2020). When we scaled to Liang et al.’s IQR of NO2, our models demonstrated similar results of 17.1% (CI: 11.3%, 23.2%), 21.3% (CI: 5.9%, 38.9%), and 30.1% (CI: 12.4%, 50.6%) increases in COVID-19 mortality rate for our zero-inflated Poisson, negative binomial, and CAR models, respectively. In comparison to the Wu et al. (2020a) study, which observed that a 1 µg m−3 increase in air pollutant PM2.5 was associated with an 11% increase in mortality rate (Wu et al., 2020a, Wu et al., 2020b), we found an 8.7 ppb increase in another traffic-related air pollutant, NO2, to be associated with a 35–60% (range of three models; Table 2) increase in mortality rate. The Wu et al. study, however, did not report results scaled to the interquartile range of PM2.5, so we scaled their results to the IQR for PM2.5 from the Liang et al. study (2.6 µg m−3), which uses a similar number of U.S. counties. This resulted in a highly comparable 31.2% increase in mortality for a 2.6 µg m−3 increase in PM2.5. Although Wu et al., Liang et al., and the current research demonstrated similar effect sizes, there may be different biological effects of NO2 and PM2.5.

The comparable effect size between our study and the Liang et al. study is notable given that the Liang et al. used large-area county-level geographies (3122 U.S. counties), and we focused on small-area neighborhoods of LA. The two studies also utilized different data sources, covariates, and model types – with Liang et al. also controlling for multiple pollutants. Our confounders were either similarly associated with our outcomes, like those included in the Wu et al. and Liang et al.’s studies, or were found to be null, reinforcing the validity of our results based on a priori expectation. The larger effect size on mortality rate in our study compared to the other two studies could be due to greater spatial variability resulting from using building-footprint covariate aggregation on smaller-area neighborhood geographies rather than using county-level data.

To our knowledge, only one study from England has undertaken small-area neighborhood analysis of the association of COVID-19 and air pollution using spatial modeling techniques (Travaglio et al., 2021). This study also reported associations with nitrogen dioxides at a sub-regional scale. Smaller-area analyses likely reduce potential exposure measurement error and lead to more consistent ascertainment of cases and deaths than those using the larger county units — both of which likely result in more precise and reliable estimates of health effects from air pollution exposures. By utilizing Bayesian models with CAR priors, we also accounted for spatial autocorrelation or clustering between adjacent administratively defined neighborhoods. This addition is important as transmission for COVID-19 and other infectious diseases is likely to be clustered spatially due to respiratory community spread (see map of incidence in Fig. 1).

Our study has several limitations. Most importantly, our study is limited by population-level counts of COVID-19 cases and deaths. These aggregate data, made publicly available by LACDPH, have facilitated this research but have also introduced some uncertainties. For example, it is difficult to determine how testing rates or the prevalence of asymptomatic cases, which may show significant neighborhood-level variation, could impact our results. We are unable to determine data accuracy, specifically for the earlier phases of the pandemic, when data collection protocols were still being defined. The number of neighborhoods reported by LACDPH has fluctuated, from 348 (in September 2020) to 335 (in February 2021) distinct areas. Our sensitivity analysis on a second time period (September 2020 to February 2021) including N = 335 neighborhoods demonstrated similar results to our initial period (March 2020 to September 202), so these potential data quality issues appear to have minimal effect on the interpretation of our results. Ideally, with data access granted, future research would avoid aggregate-level data in favor of individual-level outcomes. Utilizing individual-level home locations of cases and deaths rather than neighborhood-level aggregate counts would greatly improve air pollution exposure attribution and allow for better ascertainment of potential confounders. In addition, we could not include daily or weekly observations to account for changes in case or mortality rate over time, but rather used cumulative counts for the study period. This was due to inconsistent case reporting and incomplete death reporting due to human subject concerns earlier in the pandemic. Deaths may have been undercounted, as death certificates and coroner reports may incorrectly attribute cause of death (Jewell et al., 2020, Quast and Andel, 2020). Future research may benefit from using excess mortality for comparison. (Banerjee et al., 2020) Indicators of symptom severity from hospital, intensive care unit, and emergency room admittance data are more difficult to acquire, as they are not publicly available, but future analyses on outcomes of severity would allow for better understanding of the effect of NO2 exposure on the progression of the disease.

This study is also limited by the use of land-use regression from 2016 to estimate long-term NO2 exposure. Additional years of land-use regression surfaces could be utilized to better describe long-term trends; however, estimates of NO2 in Los Angeles County in the years immediately preceding the pandemic are likely similar in their spatial pattern over such a short time span. Other covariate data included in these analyses were spatially misaligned. A strength of this study is our use of residential building footprints as an intermediate step in aggregating areal covariates. Although using these building footprints better accounts for population density patterns than more straightforward aggregation techniques (e.g., census-tract directly to neighborhood), the method may cause misalignment errors due to differences in building characteristics (e.g., height, unit size, etc.). Finally, our spatial models accounted for unexplained spatial variability in the between-neighborhood random effect, which suggests there may be additional covariates with a similar spatial pattern that we have not included in these analyses, and further investigation is necessary. While the CAR model accounts for this dependence in the statistical inference, we cannot rule out important missing confounders.

In summary, our findings imply a potentially large association between exposure to air pollution and population-level rates of COVID-19 cases and deaths. Our findings demonstrate comparable results to other recent literature, especially concerning the association of long-term NO2 and COVID-19 mortality rate. Our small-area analyses, covariate aggregation methods using building footprints for accounting for population density variability, and utilization of spatial modeling (CAR model with spatial random effect) make novel contributions to the available literature. These findings are especially important for targeting interventions aimed at limiting the impact of COVID-19 in polluted communities.

In the U.S., more polluted communities often have lower incomes and higher proportions of Black and Latinx people. In addition, Black and Latinx people have higher rates of pre-existing conditions, potentially further exacerbating the risk of COVID-19 transmission and death (Clark et al., 2014, O’Neill et al., 2003). The elevated risk of case incidence and mortality observed in these populations might result partly from higher exposure to air pollution. As COVID-19 data reporting improves and data access is given more readily to researchers, we will further refine these analyses to the individual-level in a spatial framework.

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.

A. Case rate	Zero-inflated Poisson		Zero-inflated negative binomial		CAR zero-inflated Poisson with spatial random effect
	IRR	CI	IRR	CI	IRR	CrI
NO2 IQR	1.307	(1.288, 1.327)	1.161	(1.024, 1.316)	1.180	(1.098, 1.316)
Owner occupancy (%)	0.996	(0.996, 0.997)	0.996	(0.992, 0.999)	0.998	(0.996, 0.999)
> 65 years (%)	0.977	(0.976, 0.978)	0.956	(0.942, 0.969)	0.965	(0.951, 0.989)
Nonwhite (%)	0.995	(0.995, 0.996)	1.001	(0.997, 1.005)	0.994	(0.991, 0.997)
Smokers (%)	1.020	(1.018, 10.022)	0.992	(0.966, 1.018)	1.033	(1.009, 1.073)
Obese (%)	1.040	(1.039, 1.041)	1.033	(1.023, 1.042)	1.038	(1.035, 1.045)
Hospital density	1.123	(1.118, 1.129)	1.036	(0.966, 1.111)	1.189	(1.149, 1.232)
B. Mortality rate	Zero-inflated Poisson		Zero-inflated negative binomial		CAR zero-inflated Poisson with spatial random effect
	IRR	CI	IRR	CI	IRR	CrI
NO2 IQR	1.347	(1.225, 1.481)	1.438	(1.114, 1.857)	1.597	(1.367, 1.879)
Owner occupancy (%)	0.989	(0.986, 0.991)	0.986	(0.979, 0.993)	0.987	(0.980, 0.991)
> 65 years (%)	1.035	(1.027, 1.044)	1.040	(1.014, 1.067)	1.034	(1.012, 1.065)
Nonwhite (%)	0.997	(0.995, 0.999)	0.996	(0.991, 1.002)	0.997	(0.992, 1.002)
Smokers (%)	1.037	(1.024, 1.050)	1.049	(1.005, 1.094)	1.033	(0.981, 1.081)
Obese (%)	1.019	(1.014, 1.024)	1.025	(1.010, 1.040)	1.025	(1.009, 1.041)
Hospital density	1.061	(1.027, 1.096)	1.014	(0.902, 1.140)	1.090	(0.961, 1.289)
C. Case-fatality rate	Zero-inflated Poisson		Zero-inflated negative binomial		CAR zero-inflated Poisson with spatial random effect
	IRR	CI	IRR	CI	IRR	CrI
NO2 IQR	1.049	(0.959, 1.148)	1.207	(0.969, 1.504)	1.308	(1.100, 1.650)
Owner occupancy (%)	0.991	(0.989, 0.994)	0.993	(0.987, 0.998)	0.994	(0.989, 1.000)
> 65 years (%)	1.070	(1.061, 1.080)	1.059	(1.038, 1.080)	1.064	(1.040, 1.086)
Nonwhite (%)	1.001	(0.999, 1.003)	0.998	(0.993, 1.003)	0.998	(0.994, 1.002)
Smokers (%)	1.001	(0.997, 1.022)	1.039	(1.001, 1.079)	1.037	(0.994, 1.070)
Obese (%)	0.983	(0.977, 0.988)	0.985	(0.974, 0.997)	0.990	(0.978, 1.001)
Hospital density	0.954	(0.924, 0.985)	0.927	(0.840, 1.023)	0.979	(0.868, 1.100)

A. Case rate	Zero-inflated Poisson		Zero-inflated negative binomial		CAR zero-inflated Poisson with spatial random effect
	IRR	CI	IRR	CI	IRR	CrI
NO2 IQR	1.280	(1.261, 1.299	1.175	(1.037, 1.332)	1.272	(1.143, 1.344)
Owner occupancy (%)	0.999	(0.998, 0.999)	0.996	(0.992, 1.000)	1.000	(0.999, 1.001)
> 65 years (%)	0.974	(0.973, 0.976)	0.955	(0.942, 0.969)	0.962	(0.951, 0.972)
Nonwhite (%)	0.994	(0.994, 0.995)	1.000	(0.996, 1.004)	1.004	(1.003, 1.006)
Smokers (%)	1.014	(1.012, 1.016)	0.992	(0.965, 1.018)	1.005	(0.993, 1.020)
Obese (%)	1.034	(1.033, 1.035)	1.030	(1.021, 1.040)	1.023	(1.014, 1.032)
Hospital density	1.109	(1.103, 1.115)	1.028	(0.957, 1.103)	1.046	(1.005, 1.119)
Hypertensive (%)	0.981	(0.980, 0.983)	0.982	(0.967, 0.998)	0.978	(0.974, 0.983)
Diabetic (%)	1.073	(1.069, 1.077)	1.044	(1.005, 1.084)	1.060	(1.047, 1.077)
B. Mortality rate	Zero-inflated Poisson		Zero-inflated negative binomial		CAR zero-inflated Poisson with spatial random effect
	IRR	CI	IRR	CI	IRR	CrI
NO2 IQR	1.354	(1.233, 1.487)	1.570	(1.226, 2.012)	1.441	(1.133, 2.058)
Owner occupancy (%)	0.991	(0.989, 0.994)	0.988	(0.982, 0.995)	0.989	(0.979, 0.996)
> 65 years (%)	1.032	(1.023, 1.040)	1.034	(1.008, 1.061)	1.027	(0.997, 1.059)
Nonwhite (%)	0.996	(0.994, 0.998)	0.995	(0.989, 1.001)	0.999	(0.993, 1.006)
Smokers (%)	1.033	(1.020, 1.047)	1.054	(1.010, 1.100)	1.027	(0.973, 1.076)
Obese (%)	1.012	(1.007, 1.017)	1.021	(1.006, 1.036)	1.021	(0.999, 1.046)
Hospital density	1.047	(1.012, 1.083)	1.000	(0.887, 1.127)	1.100	(0.897, 1.242)
Hypertensive (%)	0.983	(0.974, 0.991)	0.985	(0.959, 1.012)	0.992	(0.980, 1.019)
Diabetic (%)	1.071	(1.049, 1.094)	1.040	(0.975, 1.108)	1.028	(0.943, 1.070)
C. Case-fatality rate	Zero-inflated Poisson		Zero-inflated negative binomial		CAR zero-inflated Poisson with spatial random effect
	IRR	CI	IRR	CI	IRR	CrI
NO2 IQR	1.050	(0.959, 1.149)	1.194	(0.958, 1.489)	1.336	(1.127, 1.691)
Owner occupancy (%)	0.992	(0.989, 0.994)	0.991	(0.985, 0.997)	0.994	(0.988, 1.001)
> 65 years (%)	1.071	(1.061, 1.081)	1.062	(1.041, 1.084)	1.061	(1.042, 1.084)
Nonwhite (%)	1.001	(0.999, 1.003)	0.999	(0.994, 1.004)	0.998	(0.991, 1.003)
Smokers (%)	1.007	(0.994, 1.020)	1.040	(1.002, 1.080)	1.046	(0.996, 1.105)
Obese (%)	0.983	(0.978, 0.988)	0.987	(0.975, 0.999)	0.992	(0.981, 1.002)
Hospital density	0.957	(0.927, 0.989)	0.935	(0.845, 1.034)	0.957	(0.858, 1.047)
Hypertensive (%)	1.003	(0.994, 1.012)	0.998	(0.975, 1.022)	1.010	(0.998, 1.023)
Diabetic (%)	1.001	(0.980, 1.021)	0.987	(0.932, 1.045)	0.966	(0.939, 0.999)

A. Case rate	Zero-inflated Poisson		Zero-inflated negative binomial		CAR zero-inflated Poisson with spatial random effect
	IRR	CI	IRR	CI	IRR	CrI
NO2 IQR	1.285	(1.276, 1.294)	1.212	(1.017, 1.445)	1.245	(1.101, 1.422)
Owner occupancy (%)	1.000	(0.999, 1.000)	0.996	(0.991, 1.001)	0.996	(0.995, 0.999)
> 65 years (%)	0.985	(0.985, 0.986)	0.940	(0.923, 0.958)	0.990	(0.984, 0.995)
Nonwhite (%)	0.995	(0.995, 0.995)	0.990	(0.985, 0.995)	0.991	(0.990, 0.994)
Smokers (%)	1.042	(1.041, 1.043)	1.028	(0.990, 1.068)	1.040	(1.032, 1.048)
Obese (%)	1.033	(1.033, 1.034)	1.035	(1.022, 1.049)	1.039	(1.035, 1.044)
Hospital density	1.092	(1.089, 1.095)	0.885	(0.803, 0.974)	1.220	(1.176, 1.244)
B. Mortality rate	Zero-inflated Poisson		Zero-inflated negative binomial		CAR zero-inflated Poisson with spatial random effect
	IRR	CI	IRR	CI	IRR	CrI
NO2 IQR	1.280	(1.206, 1.358)	1.377	(1.135, 1.671)	1.767	(1.339, 2.141)
Owner occupancy (%)	0.997	(0.995, 0.998)	0.999	(0.993, 1.004)	1.000	(0.996, 1.008)
> 65 years (%)	1.025	(1.019, 1.030)	0.980	(0.960, 1.002)	1.002	(0.980, 1.023)
Nonwhite (%)	1.000	(0.999, 1.001)	1.002	(0.997, 1.006)	1.003	(0.999, 1.007)
Smokers (%)	1.043	(1.035, 1.052)	1.040	(1.005, 1.076)	1.023	(0.991, 1.054)
Obese (%)	1.020	(1.017, 1.023)	1.015	(1.004, 1.027)	1.020	(1.006, 1.037)
Hospital density	1.100	(1.076, 1.124)	1.018	(0.929, 1.116)	1.074	(0.984, 1.142)
C. Case-fatality rate	Zero-inflated Poisson		Zero-inflated negative binomial		CAR zero-inflated Poisson with spatial random effect
	IRR	CI	IRR	CI	IRR	CrI
NO2 IQR	1.036	(0.977, 1.098)	1.087	(0.944, 1.251)	1.154	(1.004, 1.348)
Owner occupancy (%)	0.996	(0.995, 0.998)	0.999	(0.995, 1.003)	0.999	(0.996, 1.003)
> 65 years (%)	1.053	(1.047, 1.059)	1.037	(1.022, 1.052)	1.035	(1.018, 1.052)
Nonwhite (%)	1.005	(1.003, 1.006)	1.003	(1.000, 1.007)	1.004	(1.002, 1.007)
Smokers (%)	1.003	(0.995, 1.011)	1.008	(0.983, 1.032)	0.999	(0.977, 1.020)
Obese (%)	0.991	(0.988, 0.995)	0.991	(0.984, 0.999)	0.992	(0.985, 1.001)
Hospital density	1.011	(0.990, 1.033)	1.022	(0.960, 1.087)	1.038	(0.967, 1.105)