Emulating causal dose-response relations between air pollutants and mortality in the Medicare population

Background Fine particulate matter (PM2.5), ozone (O3), and nitrogen dioxide (NO2) are major air pollutants that pose considerable threats to human health. However, what has been mostly missing in air pollution epidemiology is causal dose-response (D-R) relations between those exposures and mortality. Such causal D-R relations can provide profound implications in predicting health impact at a target level of air pollution concentration. Methods Using national Medicare cohort during 2000–2016, we simultaneously emulated causal D-R relations between chronic exposures to fine particulate matter (PM2.5), ozone (O3), and nitrogen dioxide (NO2) and all-cause mortality. To relax the contentious assumptions of inverse probability weighting for continuous exposures, including distributional form of the exposure and heteroscedasticity, we proposed a decile binning approach which divided each exposure into ten equal-sized groups by deciles, treated the lowest decile group as reference, and estimated the effects for the other groups. Binning continuous exposures also makes the inverse probability weights robust against outliers. Results Assuming the causal framework was valid, we found that higher levels of PM2.5, O3, and NO2 were causally associated with greater risk of mortality and that PM2.5 posed the greatest risk. For PM2.5, the relative risk (RR) of mortality monotonically increased from the 2nd (RR, 1.022; 95% confidence interval [CI], 1.018–1.025) to the 10th decile group (RR, 1.207; 95% CI, 1.203–1.210); for O3, the RR increased from the 2nd (RR, 1.050; 95% CI, 1.047–1.053) to the 9th decile group (RR, 1.107; 95% CI, 1.104–1.110); for NO2, the DR curve wiggled at low levels and started rising from the 6th (RR, 1.005; 95% CI, 1.002–1.018) till the highest decile group (RR, 1.024; 95% CI, 1.021–1.027). Conclusions This study provided more robust evidence of the causal relations between air pollution exposures and mortality. The emulated causal D-R relations provided significant implications for reviewing the national air quality standards, as they inferred the number of potential early deaths prevented if air pollutants were reduced to specific levels; for example, lowering each air pollutant concentration from the 70th to 60th percentiles would prevent 65,935 early deaths per year. Supplementary Information The online version contains supplementary material available at 10.1186/s12940-021-00742-x.


Introduction
Fine particulate matter (PM 2.5 ), ozone (O 3 ), and nitrogen dioxide (NO 2 ) are major air pollutants that pose considerable threats to human health [1,2]. Starting in the 1990s, a large literature of epidemiological research has reported associations between chronic air pollution exposures and mortality, with PM 2.5 and O 3 being the most extensively studied components [3][4][5][6][7][8][9]. Chronic exposure to NO 2 has also been associated with mortality, although the evidence is relatively scarce [10,11]. These findings provide important implications for understanding the health burden attributable to poor air quality. In the United States, it is estimated that each 1 μg·m − 3 increase in PM 2.5 concentration is associated with over 30, 000 deaths each year, equivalent to a loss of 0.13-0.15 years in national life expectancy [12].
The primary objective of epidemiology is to identify a causal connection between exposure and health outcome, thereby informing decisions on policy interventions [13]. For example, the United States Environmental Protection Agency (US EPA) reviews the National Ambient Air Quality Standards (NAAQS) periodically based on the cause-effect relationship that can be inferred from the best available science [14]. However, as observational studies, many air pollution epidemiological investigations, by nature, have been associational rather than causal [15]. Although a growing literature has examined the long-term effect of PM 2.5 on mortality using the formal causal modeling techniques, there is so far little evidence for O 3 and NO 2 [16][17][18]. Indeed, O 3 and NO 2 have received less attention than they deserve; so far there is no standard for long-term O 3 concentrations (only daily) and the standard for annual NO 2 concentrations has remained the same for decades [19].
What has been mostly missing in air pollution epidemiology is the specific shapes of causal dose-response (D-R) relations between air pollution exposures and risk of mortality. Such causal D-R relations can provide profound implications in predicting the health impact at a target level of air pollution concentration [20]. Recently, a study of PM 10 that explicitly used a formal causal modeling approach to estimate the D-R relationship found a higher mortality risk at low to moderate air pollution levels [21]. However, to date no such studies have been done for PM 2.5 , O 3 or NO 2 . Specifying the causal D-R relationship, especially at very low levels, is critically important in measuring the risk of mortality induced directly by the change of air pollution level, thus supporting the potential revision of NAAQS in the US and, globally, the World Health Organization air quality guidelines [19,22].
The present study analyzed 74 million Medicare beneficiaries in the contiguous US with 637 million personyears of follow-up from 2000 to 2016, which covers more than 95% of elders aged 65 years and older in the US who are considered to be most susceptible to air pollution [23]. The Medicare population also accounts for two-thirds of total mortality, allowing us to analyze most deaths induced by air pollution [1,6]. By linking the annual averages of ambient PM 2.5 and NO2 concentrations as well as warm-season (April-September) average of ambient O 3 to the ZIP Codes of beneficiaries' residence, we were able to have proxy measures of chronic exposures for each individual [24]. We proposed a decile binning approach which divided each exposure by deciles and predicted the inverse probability of being assigned to the observed group for each observation, adjusting for the other two concurrent exposures, personal characteristics, meteorological, socioeconomic, behavioral, and medical access variables, and long-term time trend. If propensity score models were correctly specified, we had constructed a valid counterfactual framework and thus estimated the causal D-R relations between chronic exposures to PM 2.5 , O 3 , and NO 2 and the risk of all-cause mortality.

Mortality data
We obtained Medicare enrollment records for beneficiaries aged 65 years and above residing in the contiguous US between 2000 and 2016 from the Centers for Medicare and Medicaid Services, with all-cause mortality as the study outcome. For each beneficiary, we extracted their demographic information (sex, race, age at initial enrollment), Medicaid eligibility, ZIP Code of residence, year of initial enrollment, and year of death if it occurred during the study period. We constructed an open cohort with person-years of follow-up in which each beneficiary was followed each year from the study entry until the end of study, drop out of the cohort, or death, whichever occurred earliest. Note that the same data format has been used to fit time-varying Cox proportional hazard models [6].

Exposure assessment
The daily concentrations of ambient PM 2.5 , O 3 , and NO 2 at 1 km × 1 km grid cells across the contiguous US were predicted and validated using hybrid models that ensembled predictions from random forest, gradient boosting, and neural network. Multiple predictor variables were incorporated in the predictions, including ground monitoring data, satellite data, meteorological conditions, land-use variables, and chemical transport model simulations, etc., with details published elsewhere [25][26][27].
These high-resolution predictions at 1 km × 1 km grid cells allow us to estimated ZIP Code-level exposure levels with a high degree of accuracy, with annual R 2 on held out monitors of 0.89 for PM 2.5 , 0.86 for O 3 , and 0.84 for NO 2 . There are two major types of ZIP Codes in the US: standard ZIP Code and PO Box. Because a standard ZIP Code represents a delivery area, we used the polygon layer generated by Environmental Systems Research Institute (Esri) [28], and estimated the ZIP Code's daily concentrations by averaging the predictions at grid cells whose centroid points were inside the polygon of that ZIP Code. For PO Box, because it is used only for a given facility and therefore can be represented by a single point, we estimated its daily concentrations by linking it to the nearest grid cell.
The exposures of interest were assessed based on the ZIP Code-level estimates. For PM 2.5 and NO 2 , we defined their chronic exposures as annual average concentrations. For O 3 , following previous literature [5,6], we defined the chronic exposure as the average concentration during warm season (April-September) of the year. We assigned the chronic exposures to PM 2.5 , O 3 , and NO 2 to each person-year based on that person's ZIP Code of residence and calendar year.

Covariate information
Meteorological variables including daily air temperature and humidity at 2 m above the ground were extracted from Phase 2 of the North American Land Data Assimilation System, with 12 km × 12 km resolution across the continental US [29]. The average temperatures during warm (April-September) and cold seasons (January-March plus October-December) of each year were calculated from the daily data because both exceedingly low and high temperatures were physically stressful and were also associated with air pollution levels [30,31]. Annual average humidity was calculated from the daily data on a yearly basis. ZIP Code Tabulation Areas (ZCTA)-level socioeconomic variables, including the percentage of Blacks, percentage of Hispanics, median household income, median value of owner occupied housing, percentage of Americans aged 65 and older living below the poverty threshold, percentage of Americans with less than high school education, percentage of owner occupied housing units, and population density, etc., were obtained from 2000 and 2010 US Census and the American Community Survey [32]. These variables were linearly extrapolated by year to account for the time varying nature of socioeconomic status. Countylevel behavioral variables, including body mass index (BMI) and percentage of ever smokers for each year, were obtained from the Behavioral Risk Factor Surveillance System [33]. From the Dartmouth Atlas of Health Care [34], we obtained percentage of Medicare participants who had a hemoglobin A1c test, a low-density lipoprotein cholesterol (LDLC) test, a mammogram, and an eye exam to a primary care physician for each year in each hospital catchment area in the US and assigned it to all ZCTAs in that area. We also computed the distance from each ZIP Code centroid to the nearest hospital. These variables were linked to each personyear by ZIP Code of residence and calendar year. Summary statistics of the covariates are provided in Section 5 of Supplementary Information.

A decile binning approach to emulate causal D-R relations
To emulate the causal D-R relationship, we need a counterfactual framework. For a binary exposure, the causal estimate in a population of interest comes from the difference between the counterfactual outcome under which all the members of the population had been exposed versus the counterfactual outcome had they not been exposed, thus no confounding occurs [35]. In randomized experiments, counterfactuals are constructed by randomly assigning individuals to treatment groups to ensure that exposure is independent of all potential confounders. In observational studies, however, the exposure assignment is not random but instead is considered to be influenced by subject characteristics, and causal methods seek ways to approximate counterfactuals with reference to the observed population [36]. Inverse probability weighting (IPW), for example, is a formal causal modeling technique and is increasingly being used in observational studies [37,38]. For a binary exposure, it uses quasi-experimental design to construct a "pseudo-population" by weighting the population by the inverse probability of the observed exposure given all measured confounders. The "pseudo-population" is then used to estimate the exposure effect. If the systematic difference of characteristics among the exposed and unexposed is adequately adjusted so that the two groups are comparable with respect to any confounders, a causal conclusion is warranted [35,39].
But air pollution exposure is continuous in nature. Estimating the inverse probability weights in the continuous setting is challenging as it needs to 1) correctly specify the distributional form of exposure, 2) deal with non-constant variance (heteroscedasticity), and 3) avoid excessively large or small weights for outliers that are more likely to occur [40,41]. For these reasons, in this section we proposed a decile binning approach to emulate the causal D-R relations between chronic air pollution exposures and mortality by dividing each exposure into ten equal-sized groups by deciles, treating the lowest decile group (i.e., 10% of the study population with the lowest exposure levels) as the reference, and estimating the effects for the other groups compared to the reference. This relaxes the strong assumptions of distribution form and homoscedasticity for continuous exposure by relying solely on deciles. In addition, binning data makes the inverse probability weights robust against outliers [42].
We had a dataset with person-year representations of follow-up which allowed for time-varying exposures and covariates. To reduce the computational burden, first we aggregated the person-years with the same sex, race, age, Medicaid eligibility, living in the same ZIP Code of residence and in the same year. We treated them as a single record because those person-years had identical values for all exposures and covariates and thus could be treated interchangeably in the analysis. As a result, we retained all the information yet significantly reduced the size of the data in which each observation represented a stratum of combination of sex, race, age, and Medicaid eligibility per ZIP Code of residence per year. Numbers of deaths and person-years were cumulated for each stratum.
For each exposure, the analysis under the counterfactual framework was composed of two stages: a design stage where a randomized "pseudo-population" was constructed by weighting the observed population by the inverse probability of exposure given all measured confounders, and an analysis stage where the treatment effect was estimated among the constructed "pseudopopulation" [43]. In the first stage, we binned the exposure into ten equal-sized categories based on deciles. The stabilized inverse probability weight (sw ij ) for stratum j in exposure category i was defined as: where P(X ∈ i) denotes the probability of any observed exposure X being in group i, which equals to 0.1; expit(•) denotes the inverse logistic link function where expitðxÞ ¼ expðxÞ 1þ expðxÞ ; and g(•) the gradient boosting machine (GBM) model with logistic loss function for predicting the probability of the observed categorized exposure X ij given the set of confounders C, weighted by n ij , the number of person-years aggregated in the stratum. The use of GBM for estimating the probability of observed exposure has demonstrated a better predictive accuracy compared to the conventional logistic regression, as it captures nonlinearity and interactions of confounders and is unaffected by the potential autocorrelation [44,45]. The confounder set C includes the other two concurrent exposures, calendar year, the individual-level variables (sex, race, 5-year age group, and Medicaid eligibility), and the area-level meteorological, socioeconomic, behavioral, and medical access variables as detailed in the previous section. The numerator, P(X ∈ i), is used to stabilize the variability of weights to avoid excessively upweighting or downweighting observations [40].
In the second stage of the analysis, for each exposure, we fitted a log linear regression relating the number of deaths and factored exposure category, weighted by the stabilized inverse probabilities estimated from the first stage. We used quasi-Poisson link function to account for overdispersion in the number of death, and included an offset of the number of person-years to account for the different population size at risk in each stratum. As a result, we obtained the marginal effect of each decile group on mortality. If the model for estimating the stabilized inverse probability weight is correctly specified, we have achieved an unbiased estimator of the causal effect for each group [41].
The results are expressed as the relative risks (RR) of mortality for higher decile groups against the lowestdecile group (reference). The number of early deaths avoided by lowering air pollutant concentration of a higher to lower decile group can be calculated as Nα 0 ð where N is the annual averaged number of person-years, α 0 is the baseline annual mortality rate, and RR high and RR low are the relative risks of mortality for the higher decile group and the lower decile group, respectively. More details are provided in Section 4 of Supplementary Information.
We assessed the robustness of the causal doseresponse relations between the chronic air pollution exposures and mortality risk by conducting sensitivity analysis of splitting each exposure into 14 bins.

Results
The demographic characteristics of the national Medicare cohort during 2000-2016 were summarized in Table 1. The cohort included 74,537,533 Medicare beneficiaries with a total of 637,207,589 person-years of follow-up. The average follow-up time was 8.5 years.  Fig. 1. The PM 2.5 concentration was higher in most central and eastern states and the Central Valley of California, and was lower in the northeast US and mountainous region. The warm-season O 3 concentration was highest in the mountainous region and California. The NO 2 concentration was higher in populous cities and along major highways. Over the years 2000-2016, the annual PM 2.5 concentration at ZIP Codes averaged at 9.85 μg·m − 3 , the warm-season O 3 averaged at 39.34 ppb, and the annual NO 2 averaged at 17.30 ppb (Table 2). Figure 2 presents the estimated causal D-R relations between chronic exposures to PM 2.5 , O 3 , NO 2 and the RR of mortality. The exposure concentration corresponding to each effect estimate represents the average concentration within the decile group. The doseresponse relationship between chronic exposure to PM 2.5 and mortality was monotonic and approximately linear, with higher concentration levels associated with greater risk of mortality. Specifically, the RRs of mortality associated with chronic exposure to PM 2.5 ranged from 1.022 [95% confidence interval (CI), 1.018-1.025] at 6.60 μg·m − 3 (the 2nd decile group) to 1.207 (95% CI, 1.203-1.210) at 15.47 μg·m − 3 (the 10th decile group). For O 3 , the risk of mortality monotonically increased from the 2nd (RR, 1.050; 95% CI, 1.047-1.053) to the 9th decile group (RR, 1.107; 95% CI, 1.104-1.110), and dropped at the highest decile group (RR, 1.044; 95% CI, 1.041-1.048). For NO 2 , the dose-response curve wiggled at low levels and started rising from the 6th decile group (RR, 1.005; 95% CI, 1.002-1.018) till the highest decile group (RR, 1.024; 95% CI, 1.021-1.027). Importantly, the risk of mortality associated with chronic PM 2.5 exposure was substantially larger than those with O 3 and NO 2 ; the highest RR for PM 2.5 was greater than those for O 3 and NO 2 . The entire dose-response relationship for NO 2 occurred at concentrations below the national standard of 53 ppb, and most of the PM 2.5 relationship was also below the standard of 12 μg·m − 3 [19]. There is no longterm standard for O 3 . All the numerical results are provided in Section 1 of Supplementary Information.
The dose-response relations remained robust after splitting each exposure into 14 bins, with details provided in Section 2 and Section 3 of Supplementary Information.

Discussion
We proposed a decile binning approach to simultaneously emulate the D-R relations between chronic exposures to major air pollutants and mortality in a general and susceptible older population. Assuming that the IPW models were correctly specified and the counterfactual framework was valid, the D-R curves revealed that in general, higher levels of PM 2.5 , O 3 , and NO 2 were causally associated with a greater risk of mortality. Compared with previous associational D-R curves [3,[5][6][7], the causal D-R curves essentially infer the number of potential lives saved if air pollution concentrations were reduced to targeting levels. For example, lowering each air pollutant concentration from the 70th to 60th percentiles would prevent 65,935 early deaths among elders per year (Section 4 of Supplementary Information), and this is a substantial public health benefit.
A major advance of the present study is that we simultaneously evaluated PM 2.5 , O 3 , and NO 2 , which allowed us to mutually adjust their confounding and also to directly compare their health impacts. We found that PM 2.5 had a substantially larger effect on mortality than O 3 and NO 2 . The finding confirmed previously published results suggesting that PM 2.5 is the most deadly air pollutant and that chronic exposure to PM 2.5 is of greater public health concern [18]. The increasing patterns of the D-R relations for PM 2.5 and NO 2 at levels below the current NAAQS suggest the necessity of more stringent national air quality standards for the protection of public health. Currently the NAAQS lack regulation for longterm O 3 , and clearly the daily standard has not reduced the warm-season average to concentrations with no mortality association [46]. Our results support the argument for establishing a warm-season O 3 standard. The lower risk of mortality for the highest decile group for O 3 may suggest that the O 3 effect was represented by traffic exhausts such as nitrogen oxides and volatile organic compounds, as they play important roles in O 3 actions and are highly reactive at extreme levels. However, further investigations are needed to address this question.
The causal conclusions of this study depend on the key assumption of correct IPW model specification. The validity of this assumption is not testable and relies on outside information. To minimize confounding bias, we adjusted for any known possible confounders such as concurrent air pollutants, Medicaid eligibility (proxy for individual's low socioeconomic status), and seasonal temperatures and humidity (important physical stressors and determinants of air pollution [30]), etc. We also adjusted for area-level confounders of socioeconomic status, ethnicity, smoking status, obesity, population density, access to medical care, and calendar year (to capture unmeasured confounders that had a temporal scale of variation). In predicting the inverse probability of being assigned to the observed exposure decile given the set of confounders, GBM adaptively captured any nonlinearity and interactions and was unaffected by the potential autocorrelation [44]. Although residual confounding can never be ruled out, the consistent doseresponse relationships for PM 2.5 and O 3 obtained across different study designs and populations provide some reassurance that our causal estimates are not substantially biased [3,[5][6][7].
As we have noted, the proposed decile binning approach relaxed assumptions on data distribution and homoscedasticity when constructing inverse probability weights for continuous exposures. Ambient air pollution concentrations usually follow a heteroscedastic distribution possibly with long tails, which results in excessively upweighted observations. To fix this issue, Naimi et al. proposed a quantile binning approach where he estimated weights for binned exposure and then treated those bins as continuous and linear, and found it outperformed other IPW estimators with various parametric forms of the exposure distribution [42]. Adopting his idea, our approach further relaxed the linear assumption by categorizing bins and comparing the effect of each bin to the reference group. If the assumption of correct IPW model specification holds, the estimated effect of each bin is an asymptotically unbiased estimator of the true causal effect [41]. Further, the estimand of our interest, the marginal effect estimates, do not depend upon the distributions of confounders and have arguably greater public health relevance because many confounders might not be measurable at decision time. The marginal effect estimates are also more useful when depicting dose-response relationship for the purpose of understanding the total effect [20].
Assigning ambient air pollution concentrations at ZIP Codes as a marker of individual exposure levels may result in measurement error. Although measuring more personal exposures can overcome the limitation, it also introduces confounding that are difficult to control such as personal behaviors, which may affect personal exposure measurement directly but not affect ambient air pollution estimation. In addition, personal exposure measurements can be compromised by the study outcome and thus is also more vulnerable to reverse causation [24]. For example, patients who die from chronic obstructive pulmonary disease (COPD) generally spent less time outdoors [47]; because ambient air pollutants are filtered by the building envelop and deposit on indoor surfaces, there are lower concentrations of those ambient pollutants indoors [48]. Hence, those patients have lower levels of personal exposures. By contrast, under the null assumption, COPD mortality is not associated with ambient concentration predictions, which are more proxy exposure measurement than the personal measurement. In epidemiology studies, ideally the measure of exposure should be as accurate as possible. In practice, however, this is usually not possible and the issue is to choose an appropriate exposure metric that balances the biological relevance, interpretability, and implications for public health policy. While using a proxy measure for air pollution exposure increases measurement error, it also brings important advantages for causal inference.
Some limitations must be acknowledged. First, we were not able to examine on cause-specific mortality which is not available for the Medicare data. Further studies investigating which major specific causes are driving the death would provide a valuable addition. Second, spatial confounding inherent to proximity-based air pollution measurements could still be present given that ZIP Code was the finest geographical unit we could use to link air pollution levels with each beneficiary. Third, restricted by available data sources, we could not adjust for individual behavior and medical history because such information was not available for the Medicare enrollment data, which may contribute to residual confounding. Fourth, although air pollution levels were estimated from models with excellent out-of-sample prediction ability, they are not perfect and therefore may attenuate effect estimates [24].

Conclusions
In summary, this study simultaneously emulated D-R curves between chronic exposures to PM 2.5 , O 3 , NO 2 and all-cause mortality among the national Medicare cohort during 2000-2016. We proposed a decile binning approach to relax the contentious assumptions of conventional IPW estimators, which yielded more robust causal evidence on adverse effects of air pollution exposure on mortality. Assuming that the IPW models were correctly specified, the estimated D-R curves reveal that in general, higher levels of chronic PM 2.5 , O 3 , and NO 2 exposures were causally associated with a greater risk of mortality, even at levels below the national standards. Among the three pollutants, PM 2.5 posed the greatest public health concern. The estimated D-R relations provide particularly significant implications for US EPA reviewing NAAQS, as the causal D-R curves essentially infer the number of potential lives saved if air pollution concentrations were reduced to specific levels. For example, lowering the air pollutant concentration from the 70th to 60th percentiles would prevent 65,935 early deaths among elders per year.