Contribution of Satellite-Derived Aerosol Optical Depth PM2.5 Bayesian Concentration Surfaces to Respiratory-Cardiovascular Chronic Disease Hospitalizations in Baltimore, Maryland

The fine particulate matter baseline (PMB), which includes PM2.5 monitor readings fused with Community Multiscale Air Quality (CMAQ) model predictions, using the Hierarchical Bayesian Model (HBM), is less accurate in rural areas without monitors. To address this issue, an upgraded HBM was used to form four experimental aerosol optical depth (AOD)-PM2.5 concentration surfaces. A case-crossover design and conditional logistic regression evaluated the contribution of the AOD-PM2.5 surfaces and PMB to four respiratory-cardiovascular hospital events in all 99 12 km2 CMAQ grids, and in grids with and without ambient air monitors. For all four health outcomes, only two AOD-PM2.5 surfaces, one not kriged (PMC) and the other kriged (PMCK), had significantly higher Odds Ratios (ORs) on lag days 0, 1, and 01 than PMB in all grids, and in grids without monitors. In grids with monitors, emergency department (ED) asthma PMCK on lag days 0, 1 and 01 and inpatient (IP) heart failure (HF) PMCK ORs on lag days 01 were significantly higher than PMB ORs. Warm season ORs were significantly higher than cold season ORs. Independent confirmation of these results should include AOD-PM2.5 concentration surfaces with greater temporal-spatial resolution, now easily available from geostationary satellites, such as GOES-16 and GOES-17.

In 2004, the U.S. Centers for Disease Control and Prevention (CDC) and EPA established and logistically supported the CDC Public Health Air Surveillance Evaluation (PHASE) project [18,[20][21][22]. One important PHASE project outcome was the development of the first-generation HBM that statistically fused PM 2.5 monitor concentration readings with CMAQ PM 2.5 model predictions [23][24][25]. In urban areas, PMB gives more "weight" to PM 2.5 monitor readings than CMAQ PM 2.5 model predictions. In rural areas, CMAQ PM 2.5 model predictions exert more influence than PM 2.5 monitor readings on PMB, since there are fewer monitors or no monitors. Ambient air monitors are usually found in urban areas. In the last 15 years, PMB has turned out to be a more representative PM 2.5 concentration surface, compared to the interpolation of PM 2.5 monitor data, as a method to resolve spatial gaps between ambient air monitors [16,18,22]. CDC subsequently incorporated PMB into its Environmental Public Health Tracking (EPHT) network of state and New York City partners By incorporating AOD-PM 2.5 concentration values into the currently utilized PMB, we hypothesized there would be a further improvement in the fused AOD-PM 2.5 -PMB surface [16]. Our intention, in this preliminary work, was to test this hypothesis by using these four experimental AOD-PM 2.5 and PMB fused concentration surfaces with linked health outcome data from Baltimore, Maryland, and New York City, New York, in a case-crossover epidemiologic study design data files analyzed by using conditional logistic regression [16]. From the beginning, our expectation was to complete the Baltimore and New York City epidemiologic studies at the same time [16]. Unexpected circumstances related to restricted access to the Maryland emergency department (ED) and inpatient hospitalization IP confidential hospital data files, delayed the completion of the Baltimore AOD-PM 2.5 study component, and this delay made it possible for the Baltimore investigators to examine additional questions, such as completing a fine-grain evaluation of areas with and without ambient PM 2.5 monitors. The New York City study, published in 2016, did not find differences between the four experimental AOD-PM 2.5 concentration surfaces and PMB [16].
These were the objectives of the Baltimore study: firstly, replicate the New York City study by using all CMAQ 12 km 2 grids, and completing further analyses in grids with and without PM 2.5 ambient air monitors. Secondly, determine if the four experimental AOD-PM 2.5 concentration surfaces differed from PMB in grids with and without monitors. Thirdly, evaluate warm season versus cold season differences.

Description of Study Participants
This study used the same experimental methods described in the New York City publication [16]. The same procedural steps used in the two study sites are summarized in this Methods section, with more detail included on differences unique to the Baltimore study. Figure 1 is a map of the 11 (south-north) by 9 (west-east) 12 km 2 CMAQ grid cells, shown in blue circles, and the l7 Federal Reference Method (FRM) PM 2.5 ambient air monitors, shown in red triangles. Baltimore City had 6 PM 2.5 ambient air monitors, more than in Anne Arundel (4), Prince George's (3), Baltimore (2), Montgomery (1) and Harford (1) Counties. Other Counties in the study area did not have a single ambient air monitor.  2.5 concentration data files were downloaded from the EPA Air Quality System (AQS) database [17]. The locations of the Maryland study area FRM PM 2.5 monitors were mapped to the 12 km 2 CMAQ grid system prior to their incorporation into the upgraded HBM. If a grid cell contained two or more monitor concentration readings on any given day, a single daily mean was computed and retained. Grid cells without monitor measurements on a given day were set to the "missing" data condition.

CMAQ PM 2.5 -The 2004-2006
Baltimore study area CMAQ PM 2.5 model predictions were obtained from EPA's Community Modeling and Analysis System website [24,25]. CMAQ PM 2.5 model predictions did not have missing values.

AOD-PM 2.5 -Moderate
Resolution Imaging Spectroradiometer (MODIS) Collection 5 AOD 2004 data files, at 10 km 2 resolution, were downloaded from the National Aeronautics and Space Administration (NASA) Level 1 and Atmosphere Archive and Distribution System (LAADS) [45]. AOD data were available for Terra and Aqua satellites. There were two MODIS AOD observations per day, one for each satellite. Terra satellite had a late morning observation time, while the Aqua satellite had an early afternoon observation time. MODIS AOD data were re-mapped from NASA's 10 km 2 native grid system to CMAQ's native 12 km 2 grid system before the AOD data were uploaded in the upgraded HBM. Additional AOD dataset preparation details, including the technique used to estimate surface PM 2.5 concentrations from the AOD data, have been published [16,46,47].
The HBM runs for New York City and Baltimore produced two AOD-estimated PM 2.5 datasets: one AOD dataset with missing observations, and another dataset where the missing observations were replaced with kriged values (ordinary kriging). The original HBM required separating and sequencing the Bayesian computations into discrete hierarchical steps and then modeling at each level, before combining PM 2.5 monitor measurements with CMAQ model predictions [48,49]. The original HBM was modified to process input surfaces with missing data and simultaneously combine multiple datasets [46]. This latter HBM upgraded feature was "critical" to permit the incorporation of AOD-PM 2.5 input surfaces, in addition to PM 2.5 monitor concentration readings and CMAQ PM 2.5 model predictions.

AOD-PM 2.5 Fused Surfaces
The upgraded HBM was used to produce PMB and the four experimental AOD-PM 2.5 concentration surfaces, at the same time and utilizing the same statistical procedures, for Braggio et al. Page 4 Atmosphere (Basel). Author manuscript; available in PMC 2021 May 11. both New York City [16] [50]. All hospitals in Maryland must report, by State statute, their ED visits/IP hospital events to the HSCRC quarterly. Available temporal variables only included years (2004)(2005)(2006), quarters (winter, spring, summer, fall), and days of the week (Sunday through Saturday), but excluded month. Patient's residential information was limited to 5-digit United States Postal Service Zone Improvement Plan (ZIP) codes (Maryland Department of Planning; MDP) [51]. Each patient record contained demographic variables, primary and secondary diagnoses, entered as International Classification of Diseases, Ninth Revision, Clinical Modification (ICD-9-CM) billing codes [52]. The ICD-9-CM codes were used to select asthma (493), myocardial infarction (MI, 410) and heart failure (HF, 428), as primary discharge diagnoses, and the co-morbid conditions of diabetes mellitus (250), hypertension (401), and atherosclerosis (414, 440), as secondary diagnoses. The Maryland Department of Health (MDH) Institutional Review Board approved this Baltimore AOD-PM 2.5 data linkage and analysis study.

Cases-Controls-
The case-crossover design was used to form three different controls for each case [53][54][55][56][57]. Controls differed from cases only in the assigned exposure period. A case had one quarterly exposure value. Each control had one monthly mean exposure value for each of the four quarters: the first control had monthly mean exposure data for January, April, July, and October; the second control had monthly mean exposure data for February, May, August, and November; the third control combined monthly mean exposure data for March, June, August, and December. Each monthly mean exposure data represented a different calendar quarter (January-March, winter; April-June, spring; July-September, summer; October-December, fall). For each of the three controls, the selected four monthly mean exposures values were averaged, and the overall mean was used as the annual background exposure estimate.
Maclure's [56] case-crossover design paper proposed using cases as their own controls but assigning a different exposure period to the same cases that were also used as controls. Since the same cases are used as controls, the cases and controls should not differ from each other on patient attributes such as gender, age, race, and insurance coverage. This is what we did in this study, because the month temporal variable was not available in the confidential hospital files obtained from the Maryland HSCRC, but was available in the New York City study [16].
The overall average for four months, with each month serving as a proxy for each of the four quarters in one year, provided an estimate of annual (background) exposure level. Since a different month was used to represent each quarter in each of the three controls, the exposure assignment algorithm used in this study is "functionally" equivalent to the case-crossover bidirectional design used in the New York City study [16,53]. In this study, cases with 1st and 2nd quarter exposure values always preceded all three controls because their mean monthly ranks (Q 1 = 2.0, Q 2 = 5.0) were always numerically smaller than the mean monthly ranks of all three controls (1st = 5.5, 2nd = 6.5, 3rd = 7.5). However, when cases had 3rd and 4th quarter exposure values they always came after all three controls, because their mean monthly ranks (Q 3 = 8.0, Q 4 = 11.0) were always numerically larger than the mean monthly ranks of all three controls.

Co-Morbid Conditions (Diabetes, Hypertension, Atherosclerosis)-
Diabetes, hypertension, and atherosclerosis, when present in a patient's hospital record, can synergistically contribute to the occurrence of a patient's ED visit or IP hospitalization with a discharge diagnosis of asthma, MI, or HF. Diabetes mellitus, hypertension, and atherosclerosis by themselves have been shown to lead to an ED visit or IP hospitalization [13,16,[58][59][60][61][62].

Apparent
Temperature (AT and AT 2 )-Ambient temperature, relative humidity, and wind speed can contribute to the occurrence of respiratory-cardiovascular chronic disease ED visits or IP hospitalizations. AT is one summary variable that includes ambient temperature, relative humidity, and wind speed. All three weather parameters were obtained from the CMAQ model [25] and made available to the Baltimore investigators by one of the co-authors (ESH). AT was computed using the formula reported on the National Oceanic and Atmospheric Administration website (NOAA) [16,63,64]. AT 2 was computed as the product of AT*AT, once AT was available. Both AT and AT 2 have been shown to influence respiratory-cardiovascular ED visits and IP hospitalizations, even in the absence of elevated ambient PM 2.5 [13]. AT and AT 2 values were computed for each CMAQ grid cell, which were stratified by year, month, and day. [65][66][67]. The single pollen counting station in Baltimore County, Maryland, provided the multi-year pollen readings. They were used as proxy measures for ambient pollen levels in the Baltimore study area [68].

2.5.4.
Holidays-Each major holiday included the day after each holiday. The dates of recognized holidays were obtained from the U.S. Office of Personnel Management website (OPM) [69]. Each holiday was combined in a single annual dummy variable and coded as 1 (holiday and day after) or 0 (no holiday). The holiday dummy variables were entered in each of the three annual data files for 2004-2006.

Snowstorms-Each
snowstorm was coded as a dummy variable, 1 (snowstorm) or 0 (no snowstorm), for each of the three years separately (National Centers for Environmental Information (NCEI) [70]. The snowstorm variable was a proxy for physical exertion during winter snow removal and its precipitation of a cardiovascular event, such as MI taking place, followed by an ED visit or IP hospitalization.

Poverty and Population
Density-Population density and poverty co-occur more often in urban areas compared to rural locations. Brochu and colleagues [71] found a significant inverse association between economic resources and ambient PM 2.5 concentration levels. Bell and Ebisu [72] found a higher poverty percent in census tracts with ambient air monitors compared to census tracts without monitors. Both are associated with barriers to healthcare access, e.g., fewer hospitals in rural than urban areas, and more frequent use of ED medical services by persons with fewer economic resources and no health insurance. Maryland Zip Code Tabulation Area (ZCTA) poverty percent and population density (subsequently converted to the log 10 scale) data were obtained from the US Census Bureau website (USCB) [73].

File Linkage
Given that the ED visit and IP hospitalization files included temporal variables for year (2004)(2005)(2006), quarter (winter, spring, summer, fall), day of week (Sunday through Saturday) and spatial variables for residential five-digit ZIP codes, it was necessary to map these ZIP codes to CMAQ grids. It was also necessary to map ZCTA polygons for poverty percent and population density log 10 measures to CMAQ grids.
ZIP code and ZCTA latitude-longitude centroid coordinates were entered in a Geographic Information System (GIS), which included a multi-layered map of Baltimore City/Maryland Counties and a 11 (south to north) by 9 (west to east) 12 km 2 CMAQ grid map layer of the Baltimore study area, to develop a ZIP code-ZCTA-CMAQ 12 km 2 grid polygon correspondence file. This assignment was done for each year separately [51,73,80].
Latitude-longitude centroid coordinates of 5-digit ZIP codes and ZCTAs were mapped to the interior area of a unique CMAQ 12 km 2 grid cell. This ZIP code/ZCTA polygon assignment to CMAQ km 2 grid cells was done for each of the three years separately, 2004-2006. Base Statistical Analysis System (SAS) software, version 9.4, was used to link PMB and the four experimental AOD-PM 2.5 concentration surfaces, de-identified ED/IP hospitalization records, confounders and effect modifiers by using the same assigned CMAQ grid identifier  and the temporal variables of year (3), quarter (4) and day of the week (7) [81].

Case-Crossover
Analyses-Each stratum included one case and three controls. Based on the time-space grouping variables of year (three), quarter (4), day of week (7), and spatial CMAQ 12 km 2 grid cells (99), there were 8316 possible combinations of space-time locations that a case paired with three controls could be assigned. SAS/STAT Proportional Hazards Regression (PHREG) Procedure (Proc) was used to perform a regression analysis of stratified subpopulation's "survival time" data. These conditional logistic regression analyses, based on the SAS/STAT PHREG Proc, use the Cox proportional hazards model, and the results explain the effect of time-dependent explanatory variables on survival times [82]. SAS/STAT PHREG Proc, version 14.3, and Base SAS, version 9.4, software programs were used to complete all conditional logistic regression analyses [57,81,82]. In these regression analyses, ties (for [censored] survival/failure times or [uncensored] event times under the Cox proportional hazards model) were set using the Breslow methodology [82,83].

Statistical Analyses-The
Chi Square test, in Base SAS, version 9.4, was used to evaluate case-control count data differences in age categories, gender, race, co-morbid conditions and health insurance [81]. By comparing the mean of one group with the lower and upper values for the 95% Confidence Interval (CI) of the reference group mean, it was possible to determine if the two means were significantly different from each other at p ≤ 0.05. A group mean was significantly different from the reference group mean if the value for the first comparison mean was either below the lower limit or above the upper limit of the 95% CI of the reference group mean (p ≤ 0.05). The Means Procedure in Base SAS, version 9.4, was used to compute means and 95% CIs for poverty percent and population density [81].

Variable Selection-
A multi-step variable evaluation procedure was used for all conditional logistic regression runs in all CMAQ grid cells, with and without ambient air monitors [57,84]. The starting statistical run only contained PMB, or one experimental AOD-PM 2.5 concentration surface, and confounders (holidays, snowstorms, pollen, and AT and/or AT 2 ). The index ED visit/IP hospitalization day was identified as lag day 0. The day before was lag day 1, and so on, through to lag day 4. Summary lag day measures for individual lag days of 0 through 4 were also computed. Summary lag days were obtained by taking the average of the included individual lag days. To illustrate, summary lag days 0 and 1, displayed as lag days 01. Lag days 2-4 were referenced as lag days 24. Lag days of 0 through 4 were named lag days 04. The AT 2 was also entered on lag days 0, 1, 01 and 04. Effect modifiers were evaluated in separate conditional logistic regression runs (diabetes mellitus, hypertension, atherosclerosis; gender, age, race; health insurance, poverty, population density; and, season). An effect modifier was retained and included in the subsequent variable assessment runs if it had a computed probability value of p ≤ 0.20.

Variable Evaluation of Effect Modifiers-This variable evaluation phase
involved utilizing PMB or one of the four experimental AOD-PM 2.5 concentration surfaces, and retained confounders on lag days of 0-4, 01, 24, and 04, in all grids and in grids with and without monitors. AT 2 was added on lag days of 0, 1, 01 and 04. Each retained effect modifier was evaluated in a separate run. Each effect modifier with p ≤ 0.20 was evaluated in the final conditional logistic regression runs.

Final
Models-A stepwise procedure, with a variable entry criterion of p ≤ 0.20 and variable stay criterion of p ≤ 0.09, was used to identify the most parsimonious combination of confounders, effect modifiers for each PMB and the four experimental AOD-PM 2.5 concentration surfaces in all grids, and in grids with and without monitors. The Akaike Information Criterion (AIC) statistic was also used to confirm the selection of the "best" conditional logistic regression runs, with lower AIC values representing a better model fit [84]. The null hypothesis was rejected when p ≤ 0.05 [85].

Cases and Controls
Case-control group characteristics are in Table 1 (ED and IP asthma) and Table 2 (IP MI and  IP HF). As stated previously, there were three controls for each case in each of the four groups. There were more ED asthma cases (11,723) than IP asthma cases (3376), with IP MI (4745) and IP HF (6919) between two asthma case groups. ED asthma cases were significantly younger than IP asthma cases (p ≤ 0.05). The 0-14 age category included 43.8% ED asthma cases but only 27.6% IP asthma cases. The 35+ age category contained 61.8% of the IP asthma cases but only 28.7% of the ED asthma cases. There was also a significant age category difference between IP HF cases and IP HF controls (p ≤ 0.05), with lower percentages for cases (18.5%) and controls (19.2%) in the 35-59 age category and higher percentages for cases (44.9%) and controls (46.0%) in the 76+ age category. Other significant differences were found between grids without monitors and grids with monitors for poverty percent and population density for all case and control groups (all p's ≤ 0.05).  CIs) were computed for all lag days, surfaces, and monitor grid conditions. All comparisons were between each AOD-PM 2.5 concentration surface and PMB. PMCK PM 2.5 concentration means were significantly higher than PMB PM 2.5 concentration means in all grids and in grids without monitors for all lag days (all p's ≤ 0.05). The monitor PMCK PM 2.5 concentration surface was significantly lower than the monitor PMB PM 2.5 concentration surface at all lag days (p ≤ 0.05). All other comparisons between PMC, PMCQ and PMCKQ PM 2.5 concentration surfaces were also significantly lower than PMB for all lag days of 0-4, 01, 24 and 04 in all grids, in grids with monitors and in grids without monitors (all p's ≤ 0.05). Supplementary Materials includes three-year PM 2.5 concentration mean values for PMB and each of the four AOD-PM 2.5 concentration surfaces on lag days 0-4, 01, 24 and 04 in all CMAQ grids (Tables S5 and S6), in grids with monitors (Tables S7  and S8) and in grids without monitors (Tables S9 and S10).

Conditional Logistic Regression
Significant conditional logistic regression analyses for the four health outcomes, three CMAQ grid conditions, and the four experimental AOD-PM 2.5 concentration surfaces and PMB only occurred at lag days of 0, 1 and 01 (all p's ≤ 0.01). Significant, but protective, population density (M) effect modifiers were found for ED asthma in all grids for PMB (lags of 01), PMCQ (lags of 0, 1 and 01) and PMCKQ (lags of 01) (all p's ≤ 0.01 Odds Ratios (ORs) and 95% Confidence Intervals (CIs) for the four health outcomes, under the three grid conditions, PMB and the four experimental AOD-PM 2.5 concentration surfaces, are displayed graphically below, and separated into successive graphs by lag days of 0, 1 and 01. To determine if the AOD-PM 2.5 concentration surface ORs differed from the PMB OR, each AOD-PM 2.5 OR was compared to the PMB's OR (95% CI). If the AOD-PM 2.5 OR was either below or above the PMB's 95% CI lower or upper limit, the outcome was significant (p ≤ 0.05). (Figure 2A-D)-For ED asthma (Figure 2A) in all grids (Both) and in grids without monitors (No) each of the four AOD-PM 2.5 ORs were significantly higher than the PMB OR (all p's ≤ 0.05). However, in grids with monitors (Yes), only the PMCK OR was significantly higher than the PMB OR (p ≤ 0.05). For IP asthma ( Figure 2B), IP MI ( Figure 2C) and IP HF ( Figure 2D) in all grids and in grids without monitors, only PMC, PMCK and PMCKQ ORs were significantly higher than PMB ORs (all p's ≤ 0.05). For all four health outcomes, only PMC and PMCK had significantly higher ORs in the no monitor condition than the PMC and PMCK ORs in the monitor condition (all p's ≤ 0.05). Figure 3A-D)-The ED asthma ( Figure 3A), IP asthma ( Figure 3B), IP MI ( Figure 3C) and IP HF ( Figure 3D) OR results for all comparisons between each of the four AOD-PM 2.5 concentration surfaces with the PMB ORs at lag 1 in all grids, and in grids without and with monitors were identical to the comparisons made on lag day 0 and described above (all p's ≤ 0.05). No monitor versus monitor OR comparisons for all four health outcomes were also the same as those previously described above on lag day 0 (p ≤ 0.05). Figure 4A-D)-ED asthma ( Figure 4A), IP asthma ( Figure 4B), IP MI ( Figure 4C) and IP HF ( Figure 4D) comparisons between AOD-PM 2.5 ORs and PMB ORs at lag days 01 were the same as the comparisons at lag day 0 and lag day 1 and described above (all p's ≤ 0.05). For IP asthma, IP MI, and IP HF in the no monitor grid condition PMCQ ORs were significantly higher than the PMB ORs (all p's ≤ 0.05). In addition, in grids with monitors for ED asthma and IP HF PMCK ORs were significantly higher than the PMB ORs (both p's ≤ 0.05). The no monitor OR versus the monitor OR comparisons for PMC and PMCK were the same as those previously described for lag days 0 and 1 (all p's ≤ 0.05). In addition, for ED asthma the no monitor PMCKQ OR was significantly higher than the monitor PMCKQ OR (p ≤ 0.05).

Lag Day and Monitor
Further evidence supporting the robustness of PMC and PMCK concentration surfaces, especially in grids without monitors, is summarized in Figure 5A Atmosphere (Basel). Author manuscript; available in PMC 2021 May 11. to 0% change for IP MI (C) and ED asthma at lag days of 0 and 1 (A). (4) In all four panels, the PMCKQ surface had positive percent change values, thereby suggesting that kriging partly reversed the percent decrease due to the inclusion of CMAQ.

Lag Day and Season
Results of follow-up conditional logistic regression analyses evaluating warm-cold season (S) differences at lag days of 0, 1 and 01 on PMB and the four AOD PM 2.5 concentration surfaces support these conclusions: (1) only the cold season ORs were protective, e.g., below 1.000. (2) All warm season ORs were significantly higher than cold season ORs (all p's ≤ 0.05). (3) During the warm season, only ED asthma, IP MI, and IP HF PMCK ORs were significantly higher than PMB ORs (all p's ≤ 0.05). Figure 6 below shows the percent change between warm season and cold season ORs for the four AOD-PM 2.5 concentration surfaces and PMB on lag days of 0, 1 and 01 for the four hospital-based health outcomes. For all four health outcomes the largest percent increase occurred at lag day of 01. Figures S1-S3 (Supplementary Materials) display warm-cold season ORs for the four AOD-PM 2.5 experimental surfaces and PMB for ED asthma, IP asthma, IP MI, and IP HF at lag days of 0, 1 and 01, respectively.

Discussion
The shared objectives of the Baltimore and New York City [16] studies were to evaluate the differential contribution of the four experimental AOD-PM 2.5 concentration surfaces relative to the current baseline, PMB, on ED asthma, IP asthma, IP MI and HF hospitalizations. Our expectation was to find the most improvement in the new experimental AOD-PM 2.5 concentration surface, which included PMC (not kriged or kriged) fused with PMB (monitor PM 2.5 and CMAQ PM 2.5 model estimates). The New York City study did not find differences between PMB and the four AOD-PM 2.5 concentration surfaces. Contrary to our expectation, these Baltimore study results suggest that PMC and PMCK, and not PMCQ or PMCKQ, are better estimates of ambient PM 2.5 , especially in grids without monitors.
Although the same upgraded HBM was used to generate the PMB and the four AOD-PM 2.5 concentration surfaces for Baltimore and New York City, all five Baltimore surfaces had significantly higher three-year mean PM 2.5 concentration values than the New York City surfaces. Since both the Baltimore and New York City study sites analyzed the same asthma, MI and HF chronic diseases, used the case-crossover design to create three controls for each case, and used the same SAS conditional logistic regression procedure to analyze the linked exposure-health outcome files, the only remaining major difference was the significantly elevated ambient PM 2.5 levels in Baltimore compared to New York City.
The three methodological differences between the two study sites included: (1) more asthma, MI, and HF cases (and associated controls) in the New York City study than in the Baltimore study.
(2) Completion of separate ED asthma and IP asthma conditional logistic regression analyses in the Baltimore study than in the New York City study, because ED asthma cases (and associated controls) were significantly younger than the IP asthma cases (and associated controls). (3) In the Baltimore study, the three controls were formed by assigning a different exposure period to the cases that were also used as controls [56], while Braggio et al. Page 12 Atmosphere (Basel). Author manuscript; available in PMC 2021 May 11. the New York City study selected the three controls within a 28 day strata. Both test site replications attained the same case-crossover endpoint, since each case was preceded or followed by at least one control-the definition of a bidirectional case-crossover design [53,55]. It is unlikely that these methodological differences would be more important than the fact that the Baltimore study three-year mean PMB and four AOD-PM 2.5 concentration surfaces had significantly higher concentration values than the New York City study.
This Baltimore study also found differences in the contribution of PMB and the four experimental AOD-PM 2.5 concentration surfaces to asthma ED visits and IP asthma, IP MI, and IP HF hospitalizations in grids without monitors. In grids without monitors, PMC, PMCK and PMCKQ ORs were significantly higher than PMB ORs at lag days of 0, 1 and 01, for all four ED/IP respiratory-cardiovascular chronic diseases. In grids with monitors, only ED asthma PMCK ORs were significantly higher than PMB ORs for all three lag days of 0, 1 and 01. A similar outcome occurred for IP HF, but only on lag days of 01. Additional no monitor versus monitor analyses showed that only PMC and PMCK had significantly higher ORs in grids without monitors than in grids with monitors at all three lag days of 0, 1 and 01 for all four respiratory-cardiovascular chronic disease hospital events. The only exception was for IP asthma at lag day of 1, where the PMC no monitor versus monitor comparison was not significant.
Since two AOD-PM 2.5 surfaces included CMAQ PM 2.5 model estimates, PMCQ (PMC not kriged) and PMCKQ (PMCK kriged), we now know that a bias is introduced by the addition of CMAQ PM 2.5 model estimates to the experimental AOD-PM 2.5 concentration surfaces. This bias was expressed by shifting the PM 2.5 concentration values of this surface, PMCQ, toward the PMB PM 2.5 concentration values. PMCKQ resembled PMB less, a reversal that can be attributed to kriging.
Results from this Baltimore study support these conclusions: Firstly, only PMCQ ORs in all grids at lag days of 0, 1 and 01 and in grids without monitors at lag days of 0 and 1, were not significantly different from PMC ORs for IP asthma, IP MI and IP HF. Secondly, no monitor PMCQ, PMCKQ and PMB ORs were not significantly different from monitor PMCQ, PMCKQ and PMB ORs at lag days 0, 1 and 01 for ED asthma, IP asthma, IP MI, and IP HF. Thirdly, the difference in the r 2 % statistic between PMCQ and PMB in grids with monitors (97.4%) and grids without monitors (94.3%) was negatively smaller (−3.1%) than it was for PMCKQ (−15.9%), PMC (−32.4%) and PMCK (−35.6%). These results suggest that PMC and PMCK could be used as a replacement for CMAQ PM 2.5 model estimates in an upgraded HBM-generated PMB.
We also found significant season effects in the Baltimore study, which were not found in the New York City study [16]. In this study only warm season ORs were greater than 1.000, while cold season ORs were below 1.000, and therefore protective for those affected cases.
Warm season PMCK ORs at lag days of 01 were significantly higher than cold season PMCK ORs for ED asthma, IP MI, and IP HF. Some published studies have reported significant warm season effects [12,76,78], while other publications have reported significant cold season effects [74,75,77,79]. Kuo and colleagues [78] found an association between PM 2.5 concentration level and asthma hospitalization in children in the fall season. Braggio et al. Page 13 Atmosphere (Basel). Author manuscript; available in PMC 2021 May 11.
One interpretation of warm season effects, based on the positive results of this study and the non-significant outcome in the New York City study [16], is higher ambient PM 2.5 concentration levels during the warm season compared to the cold season [8].
Surprisingly, there are only two publications on PM 2.5 concentration levels and respiratorycardiovascular hospitalizations in urban Baltimore [14,15], and no published papers in rural areas. This study demonstrated, for the first time, that PM 2.5 concentration level does contribute to respiratory-cardiovascular hospital events, ED visits and inpatient stays, both in Baltimore City, with its higher population density and poverty, and in other study area grids that lacked ambient air monitors, and had lower population density and poverty. By analyzing grids without monitors, we were able to demonstrate that in these rural areas PMC and PMCK were associated with ED visits and IP hospitalizations of respiratorycardiovascular chronic conditions. Others have also found associations between ambient PM 2.5 concentration levels and respiratory-cardiovascular hospital events in rural areas [6,71], but due to course PM [86]. Strickland and associates [10] reported an association between PM 2.5 concentration and ED asthma visits, but no change due to levels of urbanicity.
There are methodological limitations and strengths to the Baltimore study. First, these results were based on 2004-2006 ED visits and IP hospitalizations that, in 2020, are more than a decade old. Baltimore ambient PM 2.5 levels were higher in 2004-2006 than in 2020.
Replicating the Baltimore study with more current HBM-generated PMB and experimental AOD-PM 2.5 concentration surfaces and ED/IP respiratory-cardiovascular chronic diseases would go a long way to confirm the generalizability of these exploratory results. A second limitation is the smaller Baltimore study area that only included 99 CMAQ 12 km 2 grids. A replication of this study should enlarge the CMAQ grid study area to contain the entire state of Maryland.
Methodological strengths include the inclusion of confounders in all conditional logistic regression analyses (apparent temperature, snow storms, major holidays, pollen apparent temperature), and the evaluation of effect modifiers (poverty, population density, season). A correspondence file was used to assign ZIP code and ZCTA polygons to a single CMAQ grid. This procedure was used to minimize differences in polygon shapes found between ZIP codes (United States Postal Service) and ZCTAs (United States Census Bureau), and to uniformly complete all analyses by only using CMAQ grids. Evaluation of conditional logistic regression runs in all three CMAQ grid conditions, all grids, and in grids with and without monitors was another methodological strength of this study.
Within this decade, AOD readings with a grid resolution of 10, 5, 2, 1 and <1 km 2 have become easier to obtain and use [27][28][29]31,[33][34][35][36]. Kumar and colleagues evaluated NASA MODIS Level 1 and 2 in grid dimensions of 10, 5, and 2 km 2 [33,34]. They concluded that AOD in 2 km 2 grids provided greater resolution than AOD in either 10 km 2 or 5 km 2 grids [34]. The major advantages of AOD data in 2 km 2 grids are: (1) higher correlation with onthe-ground ambient PM 2.5 monitors; and, (2) a decrease in lost data due to cloud interference. Another advantage is related to "scale" effects. The higher resolution found in Braggio et al. Page 14 Atmosphere (Basel). Author manuscript; available in PMC 2021 May 11. smaller grids permits more precise AOD PM 2.5 measurements of ambient PM 2.5 concentrations than in larger grids [87].
Another benefit of using AOD in <10 km 2 grids is the possibility of replacing ZIP-codeaggregated health data with address-level health data [23,30,[88][89][90][91]. Patient privacy issues require obtaining the patient's written consent before using a person's ED visit or IP hospitalization record in an environmental health epidemiologic study. The dual benefits of having available AOD data in <10 km 2 grids with address-level data include the possibility of evaluating the contribution of ambient PM 2.5 to ED visits and IP hospitalizations for respiratory-chronic diseases in both urban and rural areas, and utilizing remote sensing data to better understand shared pathophysiological mechanisms responsible for the occurrence of asthma, MI and HF events, and the development and testing of newer population-based intervention efforts.

Conclusions
If these Baltimore experimental AOD-PM 2.5 concentration surface results are confirmed, then a shift in using AOD-PM 2.5 concentrations readings in place of CMAQ PM 2.5 model predictions in the PMB may be warranted, particularly in light of the greater temporal (every 5 min during daylight hours) and spatial (2 km 2 ) data now available from geostationary satellites such as GOES-16 and GOES-17.

Supplementary Material
Refer to Web version on PubMed Central for supplementary material. Map shows Maryland's Counties and Baltimore City in the study area. The extent of the study area is defined by the 1-11 (south-north row) by 1-9 (west-east column) Community Multiscale Air Quality (CMAQ) 12 km 2 grids (blue circles      Atmosphere (Basel). Author manuscript; available in PMC 2021 May 11.