Association between Short-Term Exposure to Air Pollution and COVID-19 Mortality: A Population-Based Case-Crossover Study Using Individual-Level Mortality Registry Confirmed by Medical Examiners

Background: Studies have suggested links between ambient air pollution and coronavirus 2019 (COVID-19) mortality, yet confirmation by well-designed epidemiological studies with individual data is needed. Objectives: We aimed to examine whether short-term exposure to air pollution is associated with risk of mortality from COVID-19 for those infected with COVID-19. Methods: The Cook County Medical Examiner’s Office reports individual-level data for deaths from COVID-19 that occur in its jurisdiction, which includes all confirmed COVID-19 deaths in Cook County, Illinois. Case-crossover analysis was conducted to estimate the associations of estimated short-term exposures to particulate matter (PM) with aerodynamic diameter ≤2.5μm (PM2.5) and ozone (O3) on the day of death and up to 21 d before death at location of death with COVID-19. A total of 7,462 deaths from COVID-19 that occurred up to 28 February 2021 were included in the final analysis. We adjusted for potential confounders by time-stratified case-crossover design and by covariate adjustments (i.e., time-invariant factors, meteorological factors, viral transmission, seasonality, and time trend). Results: Of the 7,462 case and 25,457 self-control days, almost all were days with exposure levels below the PM2.5 24-h National Ambient Air Quality Standard (NAAQS) (35 μg/m3); 98.9% had O3 levels below the maximum 8-h NAAQS (35.7 μg/m3 or 70 parts per billion). An interquartile range (IQR) increase (5.2 μg/m3) in cumulative 3-wk PM2.5 exposure was associated with a 69.6% [95% confidence interval (CI): 34.6, 113.8] increase in risk of COVID-19 mortality. An IQR increase (8.2 μg/m3) in 3-d O3 exposure was associated with a 29.0% (95% CI: 9.9, 51.5) increase in risk of COVID-19 mortality. The associations differed by demographics or race/ethnicity. There was indication of modification of the associations by some comorbid conditions. Discussion: Short-term exposure to air pollution below the NAAQS may increase the mortality burden from COVID-19. https://doi.org/10.1289/EHP10836


Introduction
Because the coronavirus disease 2019 (COVID-19) pandemic has claimed millions of lives worldwide since December 2019, a better understanding of drivers of risk for severe disease and death is needed to limit this ongoing public health crisis. Accumulating evidence of positive associations between air pollution and adverse health outcomes 1,2 and possibly synergistic modification of these associations by preexisting health conditions, age, gender, and race/ethnicity [3][4][5] suggest that air pollution could contribute to the inequities of the COVID-19 pandemic. 6,7 Researchers have called for scientific rigor in studies of air pollution and COVID-19 outcomes. 8,9 Reflecting available data, most of the studies providing evidence of adverse effect estimates of air pollution on COVID-19 mortality are ecological cross-sectional studies (Table S1). Studies using this design may be affected by bias from the ecological fallacy in spite of confounder adjustment, because adjustment for individuallevel risk factors cannot be made. 7 To the best of our knowledge, as of December 2021, there have been four studies to date using individual-level data, all addressing longer-term exposure. A Mexico City study found that COVID-19 was more frequently recorded among deaths when air pollution level was higher. 10 This result did not show whether air pollution increases the overall mortality risk via COVID-19, because there were no controls in the analysis. Another study, using the United Kingdom Biobank cohort, searched for putative risk factors in terms of sociodemographics, health behaviors, preexisting conditions, and air pollution. 11 Its findings are limited by the use of air pollution estimates derived for the residence location in 2010. Findings of a Chinese retrospective cohort study are based on limited adjustment for individual risk factors and potential exposure misclassification due to unavailability of residential address and inability to consider time-varying exposure in its analysis. 12 A Southern California retrospective study focusing on the early COVID-19 pandemic might have been affected by relatively low COVID-19 ascertainment and did not consider time-varying exposure between the time of COVID-19 diagnosis and the time of mortality. 13 Collectively, these studies suggest an association of air pollution with increased COVID- 19 mortality, yet studies based on clear definition of time-varying exposures and rigorous confounding adjustment are needed to strengthen the evidence for potential links between air pollution and COVID-19 mortality.
To investigate whether air pollution increases the risk of dying from COVID-19 once infected with COVID-19, we conducted a population-based case-crossover study to estimate associations between short-term exposure to particulate matter (PM) with aerodynamic diameter ≤2:5 lm (PM 2:5 ) and O 3 and COVID-19 confirmed mortality using individual-level data.

Cook County Medical Examiner's Office COVID-19 Mortality Data
The Cook County Medical Examiner's Office in Illinois publicly provides individual-level records of all deaths with COVID-19 infection in its jurisdiction. 14 People who died from COVID-19 with confirmation by the medical examiners through 28 February 2021 were eligible study participants. The number of all COVID-19-related deaths in this data set by that time was consistent with the reported number for Cook County (including the city of Chicago) at the Illinois Department of Public Health website. 15 Data include age, sex, race, Latino ethnicity, geocoded location of death, and listed primary causes and secondary causes of death for the COVID-19-confirmed decedents. Information on age, gender, race, and address was obtained from vital records, hospitals, and police departments. Next of kin provided data related to ethnicity. Manner and cause of death (i.e., primary and secondary causes of death) were identified by the medical examiners based on autopsies, postmortem examinations, medical records, and/or information from police departments. For the analysis, we excluded accidental deaths, observations without address of location of death or date of death, and observations with COVID-19 as secondary cause of mortality. A total of 7,462 deaths with COVID-19 as primary cause and geocoded location of death within ZIP codes of Cook County were identified. Details are below.
Manner of death. We excluded 61 accidental deaths, deaths by suicide, and deaths without notation of manner of death from a total of 9,731 COVID-19-related deaths in Cook County from 16 March 2020 to 28 February 2021 based on date the deceased was pronounced dead.
Date of death. The data set of COVID-19 mortality for Cook County includes two columns for date: "Date of Death" and "Date of Incident." "Date of Death" is the date the deceased was pronounced dead. "Date of Incident" is the date of a relevant event that preceded the date that deceased was pronounced dead. For illustration, we provide two hypothetical examples: Suppose that a person started experiencing severe symptoms at home on 16 March and was hospitalized. He/she tested positive on 19 March. His/her heart stopped beating at a hospital on 24 March, and he/ she was eventually pronounced dead on the same day. "Date of Death" for this person was 24 March, and "Date of Incident" was 16 March. In addition, suppose that another person was shot on 30 March and was pronounced dead on 2 April. "Date of Death" for this person was 2 April, and "Date of Incident" was 30 March. In the data set, "Date of Incident" is generally earlier than "Date of Death." The median of the difference between these two dates was 8 d. There were also cases with differences of >365 d. A total of 28 cases had dates of incident that were later than the dates of pronouncements. For these cases, we assumed that both "Date of Incident" and "Date of Death" are the earliest dates.
For the analysis, we assumed date of COVID-19 death in several ways. First, the earliest one among the two (i.e., "Date of Incident" except for 28 cases) was used because actual harmful effects of shortterm exposure to air pollution may have occurred in an event preceding death pronouncement. For the first example above, short-term exposure to ambient air pollution may have resulted in an onset of severe symptoms on 16 March. Ambient air pollution unlikely plays a causal role on 24 March because a person in a hospital is unlikely to be exposed to ambient air pollution on that day. However, the analysis based on this date may include COVID-19 mortality risk directly increased by air pollution as well as COVID-19 mortality risk indirectly increased by air pollution (i.e., air pollution increases COVID-19 infection and thereby COVID-19 mortality risk). Therefore, we also used "Date of Death," and the analysis was focused on only COVID-19 mortality risk directly affected by air pollution after an onset of COVID-19 infection. We also restricted the analysis of COVID-19 deaths to those for which the difference between "Date of Death" and "Date of Incident" was ≤7 d to disentangle the effect of air pollution on COVID-19 mortality risk from its effect on other adverse outcomes once people became infected with COVID-19. We excluded 794 of deaths without "Date of Incident." Geocoding address of location of death. In addition to residential ZIP codes, addresses where incidents occurred (See "Date of death" section above) were recorded. For deaths without coordinates in the provided data set, we geocoded location of incident using addresses and Google Maps. We manually revisited our data set. We spotted that there were errors in geocoding of locations of deaths and mismatches between residential ZIP codes and incident ZIP codes due to coding errors and typos in address. We manually corrected typos or incomplete addresses by identifying and considering type of locations (e.g., road, intersection, airport, jail, house/nursing home/rehabilitation center, hospital), search results from Google (e.g., real estate websites), and other deaths with similar address and date of incident. Locations of incidents were mostly not hospitals: Among the 7,462 deaths, 4,881 deaths were at home or work, 2,559 deaths were at nursing homes or rehabilitation centers, and 22 deaths were at hospitals. We excluded 834 deaths due to incomplete records; disagreement between residential ZIP codes and incident ZIP codes; or locations not related to homes, workplaces, nursing homes, rehabilitation centers, and hospitals (e.g., roads, intersections, airports, jails). We further restricted our study subjects to those who died in the ZIP codes of Cook County. Some ZIP codes lie in both Cook County and neighboring counties. We did not exclude people who died in those ZIP codes. As a result, we further excluded 26 COVID-19 deaths.
Deaths from COVID-19 vs. Deaths with COVID-19. Our interest was whether air pollution plays a role in the causal relationship between COVID-19 and mortality. So, we focused on only deaths from COVID-19 (as primary cause of death). We excluded 365 deaths with COVID-19 as a secondary cause. COVID-19 may be listed as a secondary cause when COVID-19 did not directly lead to death (e.g., those who died from other diseases, but who were also infected with COVID- 19).
Merging COVID-19 deaths with environmental exposures. In the "Methods" section titled "Environmental Exposures," environmental exposures were assigned to COVID-19 deaths. This methodology resulted in the exclusion of 205 COVID-19 deaths based on "Date of Incident" and the exclusion of 250 COVID-19 deaths based on "Date of Deaths," resulting in a total of 7,507 or 7,462 study subjects, respectively. These reductions were due to missing values of the first 3-wk lagged exposures in the study period.

Comorbid Conditions of COVID-19 Deaths
We considered chronic diseases/disorders recorded as primary or secondary causes as comorbid conditions accompanying deaths from COVID-19. The comorbid conditions we considered were cardiovascular diseases (CVD), chronic obstructive pulmonary disease (COPD), asthma, diabetes, obesity, metabolic disorders, chronic kidney disease (CKD), dementia, autoimmune diseases, immunodeficiency diseases, liver diseases, multiple sclerosis, anemia, epilepsy, seizure, neurological diseases, cancers, alcohol-related disorders, drug-related disorders, hereditary diseases, malnutrition, and history of solid organ transplant. 16 For each individual, we summed the number of chronic diseases/disorders listed on by the medical examiner as an overall comorbidity indicator.

Environmental Exposures
We obtained Air Quality Survey (AQS) monitoring data from the U.S. Environmental Protection Agency (U.S. EPA) for Illinois and the neighboring states of Wisconsin and Indiana from 1 January 2020 to 28 February 2021. 18 We estimated daily 24-h average of PM 2:5 and daily 8-h maximum of O 3 at each geocoded location of COVID-19 death using inverse distance-weighting interpolation.
Our prediction method followed that suggested by Ivy et al. 19 We normalized log-transformed measured values of a pollutant at monitor i for day t (C i,t ) where NC i,t is the normalized value of log-transformed C i,t . l i is the mean of log ðC i,t Þ. r i is the standard deviation of log ðC i,t Þ.
We interpolated these normalized values for each individual's address (i.e., coordinates) in Cook County using inverse distance-squared weighting as follows: where INC j,t is the interpolated normalized value at a location j for day t. D i,j is the distance between location of death j and monitor i. INC j,t was then converted back to a concentration at j for day t ðC j,t Þ by using where l j and r j were estimated using inverse distance weighed interpolation.
To validate the exposure estimation, we conducted leave-one monitoring station-out cross validation (LOOCV) with respect to monitoring stations within Cook County and one station in a county near the Cook County border. Thus, prediction was made for each area where a monitor is located using the other monitored values. We regressed predicted values on monitored values. LOOCV showed good agreement between estimated values and monitor values at AQS sites (R 2 for PM 2:5 = 0:65, R 2 for O 3 = 0:83). Rootmean-square error (RMSE), mean bias (MB), mean error (ME), normalized mean bias (NMB), normalized mean error (NME), Spearman correlation (SC), and Pearson correlation (PC) are shown in Tables S2-S3. The map of locations of AQS stations is provided in Figure 1. The average of C j for PM 2:5 and O 3 for a total of 7,462 COVID-19 confirmed deaths is also shown in Figure 1.
We also estimated population-weighted concentrations for ZIP code of location of death, C ZIP code,t as where P k is a 5-y estimate of population from the Census American Community Survey (ACS) 2018 for census block group k. We analyzed whether associations between air pollution and COVID-19 mortality based on C j,t as the main results are consistent with those based on C ZIPcode,t . Because most of the AQS monitors measured O 3 only for March-October, whereas PM 2:5 was measured year-round, we restricted analysis for O 3 and COVID-19 mortality to March-October 2020. Modeled O 3 from January-February and November-December was used only for adjustment in analysis for PM 2:5 and COVID-19 mortality.
We obtained 4-km resolution daily average of temperature and dew-point temperature modeled by the Parameter-elevation Relationships on Independent Slopes Model (PRISM) from the PRISM Climate Group. 20 ZIP code-level temperature and relative humidity (estimated using temperature and dew-point temperature) were calculated.

Statistical Analyses
We estimated the lag structure of associations between PM 2:5 and O 3 and COVID-19 related mortality using constrained distributed lags up to 21 d (i.e., lag0-21) to allow for effects of cumulative exposure, possible lag effects, and short-term mortality displacement. 21 We constrained distributed lags using a natural cubic spline with three internal equally spaced (the logarithm-based) knots. This lag0-21 exposure can capture exposure effects across the previous 3 wk and short-term mortality displacement roughly up to 21 d of life lost. 21 PM 2:5 and O 3 were coadjusted.
Our study design is a time-stratified case-crossover design that is widely used in studies of short-term environmental exposures. [22][23][24] With population-based representative mortality data, this design with conditional logistic regression produces rate ratios (approximately risk ratios) as time-series analysis with Poisson regression does. 22,24,25 Time-stratified case-crossover analysis compares deaths and their self-control days at prespecified time points. We selected self-control days as in the same calendar year, month, and day of the week of death, resulting in three or four self-control days for each death, which is widely used to adjust for seasonality and long-term time trend. 5 Timeinvariant factors between an index day (i.e., day of the occurrence of COVID-19 related death) and a control day, such as age, sex, race/ethnicity, socioeconomic status, occupation, unlikely changing health behaviors in a short period of time (e.g., smoking), and other unmeasured time-invariant factors, were adjusted by the design. 22,26 For noncommunicable health outcomes (e.g., cardiovascular mortality), time-stratified case-crossover designs may make sufficient adjustment for some unmeasured time-varying confounders such as seasonality and long-term time trend. Time-varying confounders may also include meteorological factors such as temperature, which are sometimes adjusted by matching case and their self-control time intervals (e.g., days) with the same value of temperature but which can also be adjusted in regression analysis. In this study, we adjusted for temperature and relative humidity in the regression analysis by including distributed lag nonlinear terms with a lag period of 21 d (lag0-21). 27 To reduce the number of parameters to be estimated, we used moving-average constraints for a lag structure. So, in regression, we added terms as below Environmental Health Perspectives 117006-3 130(11) November 2022 where NS is a natural cubic spline with 3 degrees of freedom (df), T t is temperature at time t, and RH t is relative humidity at time t. For infectious diseases, it is noteworthy that residual autocorrelation in time-series and case-crossover studies may be an indicator of outbreaks (e.g., cluster infections) for which adjustment is a challenge. 8,9,22,28 This autocorrelation may be seen as the effect of an additional time-varying confounder that may not have been adjusted by the time-stratification mentioned above. Mortality risk of one individual is affected by the possibility of SARS-CoV-2 transmission from neighbors through such as a) infection or b) if already infected, increasing the viral load in his/her cardiopulmonary system. The viral transmission may occur regardless of air pollution, but it may be coincidentally correlated with air pollution. In the context of case-crossover design, it can be seen as the likelihood of the transmission is different between control days and an index day, such that the transmission becomes a time-varying confounder like other time-varying confounders. Examples of this confounding may include cluster outbreaks in nursing homes and other congregate living facilities.
To adjust for this additional time-varying confounding, we developed a spatiotemporal metric based on ZIP code-level daily number of COVID-19 cases in Cook County and surrounding counties to measure the risk from the viral transmission by neighbors. The idea is that the spatial distribution of daily number of COVID-19 cases represents the distribution of transmission across space and time. We obtained ZIP code-level daily number of COVID-19 cases in Illinois and daily number of COVID-19 cases in Lake County in Indiana. For ZIP codes where daily number of COVID-19 cases during the early phase of the outbreak were not publicly available, we predicted values using autoregressive integrated moving-average models based on daily timeseries of cumulative COVID-19 cases of the other days. This metric consists of the two components in regression: and where A 1 is the first component of this metric, A 2 is the second component of this metric, s is spatial unit (i.e., ZIP code), S is the total number of ZIP codes in Cook County and surrounding counties, t is time, C s,t is daily number of COVID-19 confirmed cases, and d s,u is the Euclidian distance between ZIP code s and ZIP code u (for Indiana, county u). The second component is inverse-distance (squared) weighted average of C around ZIP code s (not including s). In regression analysis, we added distributed lag nonlinear terms of this metric with a lag period of 21 d (i.e., lag0-21), given the contagion period of COVID-19. We allowed for nonlinear associations between this metric and COVID-19 related mortality because the shape of relationship between COVID-19 cases and COVID-19 mortality is unknown, and nonlinear terms would lead to adequate adjustment if nonlinear associations exist. To reduce the number of parameters to be estimated, we used moving-average constraints for a lag structure. So, in regression analysis, we added terms as below Because COVID-19 is highly contagious, if the viral transmission as a time-varying confounder were poorly controlled, positive residual autocorrelation would arise. 22,28 If it were adequately adjusted, positive residual serial autocorrelation would disappear. Because tests for autocorrelation in the setting of conditional logistic regression in time-stratified case-crossover designs are currently not tangible, we instead fitted a conditional Poisson regression model using aggregated time-series data from our individual-level data set as an alternative of conditional logistic regression 29 and then estimated autocorrelation functions and partial autocorrelation functions. In the null model (i.e., no variables were included, and only an intercept exists), we found strong positive residual autocorrelations from lag1 to lag5 (>0:10). When variables for meteorological factors and air pollutants were included into a model, positive autocorrelations from lag2 to lag5 decreased to <0:10, but positive autocorrelation from lag1 still remained >0:10. When the viral transmission was further included, positive autocorrelations decreased to <0:10, suggesting that confounding by the infection between individuals such as cluster outbreaks was adequately adjusted.
To test nonlinear associations between air pollutants and COVID-19 mortality, we used constrained distributed lag nonlinear terms. We used a moving average of lag0-2 and a moving average of lag3-21 to reduce the number of parameters to be estimated. Then, we added these two moving averages as a natural cubic spline with 3 df for each. We compared values of Akaike Information Criteria of a model with constrained distributed lag linear terms and a model with constrained distributed lag nonlinear terms.
Our main analysis was an analysis for COVID-19 mortality based on "Date of Death" (n = 7,462). Analyses based on "Date of Incident" (n = 7,507) and restricted to COVID-19 deaths with ≤7 d of the difference between "Date of Death" and "Date of Incident" (n = 3,586 for "Date of Incident" and 3,584 for "Date of Death," respectively) were also performed. Analyses restricted to COVID-19 deaths whose incident occurred at home or work were also performed (n = 4,923 for "Date of Incident" and 4,881 for "Date of Death," respectively). All subsequent analyses were performed for only the main analysis.
We used R software (version 3.5.3; R Development Core Team) and fitted conditional logistic regression using the survival package. Constrained distributed lag linear terms of air pollutants were analyzed using the dlnm package.

Results
During the study period, averages of daily pollution levels for the locations of the case periods for the 7,462 deaths from COVID-19 based on "Date of Death" and the 25,457 self-control days were 8:8 lg=m 3 for PM 2:5 and 19:2 lg=m 3 for O 3 . The interquartile range for PM 2:5 and for O 3 was 5:2 lg=m 3 and 8:2 lg=m 3 , respectively. Almost all had exposure levels below the PM 2:5 24-h NAAQS (35 lg=m 3 ); 99.4% had O 3 levels below the maximum 8-h NAAQS [35:7 lg=m 3 or 70 parts per billion (ppb)]. The Pearson correlation between these two pollutants in the matched sample was −0:01. Time-series and box plots for PM 2:5 and O 3 are presented in Figure 2.
Among the total of 7,462 COVID-19 confirmed deaths, 5,712 (76.5%) were ≥65 y old and 4,313 (57.8%) were men ( Table 1). By race/ethnicity, 2,115 (28.4%), 1,568 (21.0%), and 3,174 (42.5%) were Black, Latino White, and non-Latino White, respectively. According to 5-y estimates of the ACS 2018, the population of Cook County is 13.9% age ≥65 y and 48.5% male. The percentages of the Cook County population that are Black, Latino White or non-Latino White are 23.6%, 14.2%, and 42.5%, respectively. Median number of comorbid conditions, CCI, and age-adjusted CCI were 2, 1, and 4, respectively. The number of comorbid conditions was correlated with CCI (Spearman rank correlation: 0.72) and ageadjusted CCI (0.49). CCI was correlated with age-adjusted CCI (0.58). The Spearman correlations between age and the number of comorbid conditions, CCI, and age-adjusted CCI were 0.14, 0.04, and 0.67, respectively. Table 2 presents associations between PM 2:5 and COVID-19 mortality. An IQR increase (5:2 lg=m 3 ) in 3-wk (lag0-21) PM 2:5 exposure was associated with a 69.6% [95% confidence interval (CI): 34.6, 113.8] increase in COVID-19 mortality based on "Date of Death" and all locations of incidents (n = 7,462). Figure  3A shows associations for lag0 PM 2:5 , lag0-1 PM 2:5 , . . ., and lag0-21 PM 2:5 , showing that the association is higher as the period of PM 2:5 exposure is longer. Estimates for COVID-19 mortality based on "Date of Incident" were higher by roughly 1.6 times to 2 times than estimates for COVID-19 mortality based on "Date of Death," whereas 95% CIs overlapped. Estimates based on COVID-19 deaths occurring at all locations were consistent with estimates based on COVID-19 deaths occurring at home or work. Estimates restricted to deaths whose difference between "Date of Death" and "Date of Incident" was ≤7 d were higher than estimates for all deaths, whereas their 95% CIs overlapped. Table 3 presents associations between O 3 and COVID-19 mortality. An IQR increase (8:2 lg=m 3 ) of 3-d (lag0-2) O 3 exposure adjusted for lag3-21 O 3 was associated with 29.0% (9.9 to 51.5) increase in COVID-19 mortality based on "Date of Death." Both lag0-2 and lag0-21 for O 3 are presented in Table 3 because the association for O 3 decreased after lag0-2, unlike that for PM 2:5 ( Figure 3B). Lag0-21 O 3 exhibited an association of −9:0% (−49:6 to 64.3), suggesting potential short-term mortality displacement. Estimates for COVID-19 mortality based on "Date of Incident" were lower than estimates for COVID-19 mortality based on "Date of Death," whereas 95% CIs overlapped. Estimates based on COVID-19 deaths occurring at all locations were consistent with estimates based on COVID-19 deaths occurring at home or work. Estimates restricted to deaths whose difference between "Date of Death" and "Date of Incident" was ≤7 d were higher than estimates for all deaths, whereas their 95% CIs overlapped. In stratified analyses by age, sex, and race/ethnicity, modification by age and race/ethnicity was statistically significant for lag0-21 PM 2:5 ( Figure 4). PM 2:5 exposure was associated with a 149.1% (53.5 to 304.5), 90.1% (19.6 to 202.2), 89.5% (19.1 to 201.5), and 0% (−38:3 to 62.0) increase in COVID-19 mortality for age groups of <65, 65-74, 75-84, and ≥85 y, respectively. PM 2:5 exposure was associated with a 78.6% (21.8 to 161.9), −17:0% (−49:6 to 36.9), and 175.3% (77.7 to 326.3) increase in COVID-19 mortality for non-Latino White, Latino White, and Black, respectively. Black males were the group with the highest estimate of a 353.4% (142.2 to 748.9) increase in mortality risk, followed by Blacks age 65 y or older, who had a 193.0% (75.3 to 389.5) increase in mortality. For lag0-2 O 3 , there was no statistically significant modification, although some groups showed higher and significant associations (i.e., ≥85 y, female, Black, Black age 65 y or older, and Black female). Figures S1 and S2 further present associations for cumulative lag of PM 2:5 and O 3 (lag0 to lag0-21) by age, sex, and race/ethnicity. Figure 5 shows the associations between the air pollutants and COVID-19 mortality by comorbid conditions. We found statistically significant modification by the number of comorbid conditions for lag0-21 PM 2:5 . PM 2:5 exposure was associated with a −0:1% (−35:8 to 55.6), 17.3% (−24:4 to 82.0), 166.8% (60.9 to 342.4), and 227.0% (96.3 to 444.6) increase in COVID-19 mortality for 0-1, 2, 3, and ≥4 comorbidity conditions, respectively. For ≥2 Charlson index, PM 2:5 was associated with a 110.0% (34.9 to 226.9) increase in mortality risk, which was higher than a 42.3% (−3:5 to 109.9) increase in the mortality for 0 of Charlson index, but this was not statistically significant modification. For lag0-2 O 3 , we did not find such modification, but for 0 of Charlson index, O 3 was associated with a 47.1% (12.3 to 92.6) increase in mortality, which was higher than a 1.5% (−25:3 to 38.0) increase in the mortality for ≥2 of Charlson index. Figure 6 displays the associations between air pollutants and COVID-19 mortality by selected comorbid conditions. We found statistically significant modifications by some conditions for lag0-21 PM 2:5 . For those with cardiovascular diseases, PM 2:5 was associated with a 99.5% (53.5 to 159.2) increase in COVID-19 mortality, which was higher than −0:9% (−40:2 to 64.2) increase for those without cardiovascular diseases. For those with hypertension, PM 2:5 was associated with a 99.5% (51.6 to 162.6) increase in the mortality, which was higher than 16.4% (−24:8 to 80.2) increase for those without hypertension. For lag0-2 O 3 , we found statistically significant modifications by IHD and dementia. For those with IHD, O 3 was associated with 81.0% (27.0 to 158.1) increase in the mortality, which was higher than 18.4% (−1:2 to 42.0) increase in the mortality for those without ischemic heart disease. For those with dementia, O 3 was associated with 165.3% (49.2 to 371.7) increase in mortality, which was higher than 19.5% (1.0 to 41.5) increase in the mortality for those without dementia. Figures S3-S6 further present associations for cumulative lag of PM 2:5 and O 3 (lag0 to lag0-21) by comorbid conditions. When we estimated associations between PM 2:5 and O 3 and COVID-19 mortality using nonlinear terms, the models with linear terms were a better fit based on AIC ( Figure S7). Associations between air pollution and COVID-19 mortality based on C j,t were consistent with those based on C ZIP code,t (Tables S8 and S9).

Discussion
Our findings show associations of both particulate matter and ozone concentrations with COVID-19 mortality, implying that reduction of air pollution may have an immediate and beneficial  Table S4 for corresponding numerical data. The estimated associations between PM 2:5 and COVID-19 mortality are strong and suggest that the effects of air pollution on COVID-19 mortality are higher than on noncommunicable diseaserelated mortality. A recent systematic review commissioned by the World Health Organization found that a 10 lg=m 3 increase in daily PM 2:5 was associated with a 0.65% increase in all-cause mortality and a 0.73% increase in respiratory mortality. 31 A 10lg=m 3 increase in daily O 3 was associated with a 0.43% increase in allcause mortality. 31 In our study, a 5:2 lg=m 3 increase in PM 2:5 was associated with a 69.6% increase of COVID-19 mortality (i.e., COVID-19 as primary cause of death). An 8:2 lg=m 3 increase in O 3 was associated with a 29.0% increase of COVID-19 mortality. These estimates were adjusted for individual-level risk factors by the case-crossover design, temporal confounding by COVID-19 transmission, meteorological factors, and the other pollutants (i.e., PM 2:5 or O 3 ). Consequently, residual confounding is unlikely to explain the comparative strength of the estimates.
Several factors may contribute to these strong effect estimates. First, more severe infectious diseases may be strongly related to air pollution exposures. For adults age 65 y or older in the Wasatch Front of Utah, a 10 lg=m 3 increase in 1-d PM 2:5 was associated with a 38% (6 to 80) increase in severe pneumonia risk and a 50% (3 to 116) increase in mortality risk of pneumonia inpatients. 32 For all ages in the same region, a 10 lg=m 3 increase in PM 2:5 for a few weeks was associated with ≥20% increase in influenza-related acute respiratory infection. 33 Second, the effect sizes reported in the literature typically differ by disease and exposure duration. Short-term time-series studies in South Korean cities, European cities, and Sao Paulo, Brazil, found that longer-term exposure to particulate pollution is highly associated with respiratory mortality but not with cardiovascular mortality. 34,35 Findings were similar in a study in 708 counties of the United States for cardiovascular and respiratory hospital admissions. 36 For a few days of exposure, the associations with hospital admissions for pneumonia and respiratory infections were similar to those for other diseases in general populations of the United States and China. 37,38 In addition, integration of evidence from time-series, cohort, and quasi-experimental studies for associations between fine particles and cardiovascular mortality suggests that the association becomes higher with longerterm exposure, 39 which is also supported by considerations of biological plausibility. 40 Although our main outcome is COVID-19 mortality, the strong effect estimates for PM 2:5 in comparison with effect estimates for other noncommunicable diseases may be related  Table S5 for corresponding numerical data. Note: CI, confidence interval.
to the prior evidence and our results regarding positive modification by cardiovascular diseases.
We found higher effect estimates for PM 2:5 in those with higher number of comorbid conditions, cardiovascular diseases, and hypertension. Although existence of modification by comorbidity status regarding associations between PM 2:5 and mortality differed across studies, some studies found that those with cardiovascular diseases including myocardial infarction and hypertension may be at higher risk of PM 2:5 -related mortality. 4,41 Effect estimates for O 3 was higher in those with IHD, which is consistent with higher estimates for all-cause mortality in those with myocardial infarction. 3 Those with dementia were found to have higher risk of O 3 -related COVID-19 mortality. Although there is a paucity of evidence of associations between short-term exposure to air pollution and dementia, 42,43 people with dementia may have limited capabilities to address COVID-19 (e.g., limited access to information and preparedness for emergency situations). 44 To the best of our knowledge, as of December 2021, there are only three published studies on effects of clearly specified shortterm exposure to air pollution on COVID-19 mortality in the general population, considering confounding adjustment (Table S1). They are all ecological studies for large spatial areas, unlike our analysis with its individual-level geocoded data and exposure assessment. One time-series study found that an 18:4 lg=m 3 increase in PM 2:5 was associated with a 5.8% (3.4 to 8.2) increase in COVID-19 related mortality (i.e., laboratory confirmed diagnosis of COVID-19 or clinical or epidemiological diagnosis of COVID-19) in nine sectors of the province of Santiago in Chile. 45 That study did not find a positive association between O 3 and COVID-19 related mortality. Other ecological cross-sectional time-series study found an association between a 10 lg=m 3 increase in PM 2:5 for 28 subsequent days and an 8.4% (2.1 to 15.3) increase in COVID-19 mortality in 92 counties of states of Washington, Oregon, and California in the United States. 46 This study found that the association was heterogenous across counties. For example, San Bernardino, California, had the strongest association, which was a 65.9% (22.8 to 105.3) increase in COVID-19 mortality for a 10-lg=m 3 increase in PM 2:5 . The strength of the association for PM 2:5 in our study was approximately 42 times that estimated in the Santiago study 45 and 2 times that estimated for San Bernardino, California. 46 Our study has several strengths. We report individual-level associations between short-term exposure to PM 2:5 and O 3 and COVID-19 mortality in a general population, thereby substantially improving on evidence from earlier studies that investigated similar issues using spatially aggregated data. We used a casecrossover design that adjusted for measured and unmeasured time-invariant confounders by design (e.g., age, sex, race/ethnicity, socioeconomic status, occupation, smoking, obesity, etc.). By using time-stratification and measured variables, we adjusted for  Table S6 for corresponding numerical data. Note: CI, confidence interval. time-varying confounders (i.e., meteorological factors, the viral transmission between individuals, seasonality, and long-term time trend). We assigned exposure estimates at individual-level geocoded locations of deaths. Thus, our findings are not subject to the ecological fallacy and unlikely suffer from residual confounding. An exposure time window (up to 3-wk) was assigned identically to all cases. Proper temporality of the relationship between air pollutants and COVID-19 mortality is implicit. Furthermore, classification of COVID-19 deaths was confirmed by the Cook County Medical Examiner's Office. All these factors may have contributed to different effect estimates than were obtained in previous studies. Our study also considered shortterm mortality displacement, which can mask the actual relationship between air pollutants and COVID-19 mortality. 21 We found evidence of short-term mortality displacement for O 3 (i.e., significant positive associations disappeared for longer cumulative lags; Figure 3B), but we did not find diminishing associations for PM 2:5 ( Figure 3A). Based on the concept of mortality displacement, 47,48 this finding for O 3 suggests that reduction of O 3 becomes more important when reducing multiple risk factors of COVID-19 prognosis.
We note several limitations of this study. First, additional causes of mortality for deaths from COVID-19 may not fully capture preexisting comorbid conditions. Documentation of additional causes in death certificates may be suboptimal, and completion of death certificates may differ by time and demographics 49 ; such factors cannot be resolved in our study. Second, although our exposure model for air pollution shows good agreement between predicted values and measured monitor values; nevertheless, there is inevitably some exposure misclassification. The impact of exposure misclassification in short-term air pollution study designs might have affected our estimates toward the null. 50 Third, our data did not include any steps for personal protection. Such individual factors were adjusted by our time-stratified case-crossover design, but we were not able to investigate whether the associations between air pollution and COVID-19 mortality are modified by COVID-19 protective measures. Fourth, we were unable to investigate effect modification by socioeconomic status and other health behaviors (e.g., smoking) because our data did not include such information. Again, these time-invariant confounders were adjusted by the case-crossover design. Fifth, we could not consider all possible pathways for increased risk of COVID-19 mortality because the time of COVID-19 infection of each individual was unknown. Short-term exposure to air pollution may increase the risk of mortality from COVID-19 by at least two pathways: increasing risk of COVID-19 infection and exacerbating COVID-19 prognosis once infected with COVID-19. The data set we used does not have information on when people experienced symptoms, when they were infected with COVID-19, when they tested positive, and when/ whether they were hospitalized. By the nature of our data and data analysis, our main results, which were based on "Date of Death," reflect the second pathway. That said, the increased risk by the first pathway might have been captured to some degree if some of individuals were infected with COVID-19 within 21 d prior to a relevant event preceding the pronouncement of death. Our analyses based on "Date of Incident" may include the first pathway, which  Table S7 for corresponding numerical data. Note: CI, confidence interval.
showed higher effect estimates than those based on "Date of Death" for PM 2:5 but not for O 3 . Nevertheless, our results are silent about the degree of the increased risk by the first pathway. Sixth, we did not investigate exposure periods longer than 3 wk because time-stratified case-crossover designs may have a multicollinearity problem, 22,51 yet long-term exposure to air pollution may also be relevant for future studies. 6,7 In conclusion, our results support the hypothesis that shortterm exposures to PM 2:5 and O 3 , even at levels below current air quality standards, increase COVID-19 mortality risk.