The Shape of the Concentration-Response Association between Fine Particulate Matter Pollution and Human Mortality in Beijing, China, and Its Implications for Health Impact Assessment.

BACKGROUND: Studies found approximately linear short-term associations between particulate matter (PM) and mortality in Western communities. However, in China, where the urban PM levels are typically considerably higher than in Western communities, some studies suggest nonlinearity in this association. Health impact assessments (HIA) of PM in China have generally not incorporated nonlinearity in the concentration-response (C-R) association, which could result in large discrepancies in estimates of excess deaths if the true association is nonlinear.
OBJECTIVES: We investigated nonlinearity in the C-R associations between with PM with aerodynamic diameter [Formula: see text] ([Formula: see text]) and mortality in Beijing, China, and the sensitivity of HIA to linearity assumptions.
METHODS: We modeled the C-R association between [Formula: see text] and cause-specific mortality in Beijing, China (2009-2012), using generalized linear models (GLM). [Formula: see text] was included through either linear, piecewise-linear, or spline functions to investigate evidence of nonlinearity. To determine the sensitivity of HIA to linearity assumptions, we estimated [Formula: see text]-attributable deaths using both linear- and nonlinear-based C-R associations between [Formula: see text] and mortality.
RESULTS: We found some evidence that, for nonaccidental and circulatory mortality, the shape of the C-R association was relatively flat at lower concentrations of [Formula: see text], but then had a positive slope at higher concentrations, indicating nonlinearity. Conversely, the shape for respiratory mortality was positive and linear at lower concentrations of [Formula: see text], but then leveled off at the higher concentrations. Estimates of excess deaths attributable to short-term [Formula: see text] exposure were, in some cases, very sensitive to the linearity assumption in the association, but in other cases robust to this assumption.
CONCLUSIONS: Our results demonstrate some evidence of nonlinearity in [Formula: see text]-mortality associations and that an assumption of linearity in this association can influence HIAs, highlighting the importance of understanding potential nonlinearity in the [Formula: see text]-mortality association at the high concentrations of [Formula: see text] in developing megacities like Beijing. https://doi.org/10.1289/EHP4464.

Introduction

Chinese cities face some of the worst urban air pollution in the world, in part because of China’s rapid economic expansion and industrial development (Chan and Yao 2008). This growth has been accompanied by an expansion of urban population (Zhang et al. 2010), with several megacities () emerging since the 1990s (Chan and Yao 2008), increasing exposure to dangerous air pollution. This severe air pollution has been associated with substantial risk of all-cause mortality, cardiorespiratory mortality, and cardiorespiratory hospital admission in Chinese cities (Zhang et al. 2010; Kan and Chen 2004; Wang and Mauzerall 2006; Voorhees et al. 2014), including in Beijing (Zhang et al. 2007; Guo et al. 2013).

In particulate matter (PM) epidemiologic studies, researchers often investigate the health effects associated with ambient PM concentrations by estimating and examining a concentration–response (C-R) association. Direct evidence of both the shape and magnitude of this C-R association is critical in conducting public health impact assessments (HIAs) (Ostro 2004). Because estimated C-R associations from time-series studies only capture the effects of short-term, not long-term, exposure of PM on mortality risk, they cannot, by themselves, be used to estimate years of life lost (YLL) in adults without substantial assumptions (Ostro 2004; Cohen et al. 2004). However, they still represent a critical tool in HIA; for example, some researchers have used time-series results to estimate YLL among children (Cohen et al. 2004), and others have incorporated estimates of premature deaths attributable to acute PM exposure as one element of HIA (Voorhees et al. 2014; Matus et al. 2012) or as a basis for comparing short-term and long-term impacts of PM (Ostro 2004). In China-based studies, such time series–based approaches are often the only way of including evidence from China-based epidemiologic studies in HIAs, given the paucity of cohort studies of PM and health in the country. Further, estimates of short-term C-R associations are key in conducting short-term HIAs of the immediate health benefits of Blue Sky days—short-term air pollution control periods like the Asia-Pacific Economic Cooperation Summit in 2014—as in Lin et al. (2017), and in interpreting results from studies that investigate evidence of such benefits (e.g., Lin et al. 2014). Estimates of short-term C-R associations will likely also help inform future epidemiologic assessments of the immediate impacts of the expanding measures taken to reduce PM with aerodynamic diameter () on days when air pollution warnings are at high-alert levels (e.g., limiting vehicle use; Tambo et al. 2016).

Extensive studies have investigated the shape of the association between ambient PM concentration and the logarithm of daily counts or rates of nonaccidental or cardiorespiratory mortality in U.S. and European communities, in particular to look for evidence of a threshold below which PM concentration and mortality were unassociated (Daniels et al. 2000; Schwartz and Zanobetti 2000; Schwartz et al. 2001). In these Western studies, the C-R association was found to be approximately linear, without strong evidence of any such threshold (Daniels et al. 2000; Schwartz and Zanobetti 2000; Schwartz et al. 2001). However, the PM levels in the Western communities used in these studies are considerably lower than those typical of Chinese megacities like Beijing. For example, in a study of the shape of the C-R association between fine PM () and nonaccidental mortality in six U.S. communities, the C-R association was assessed up to (Schwartz et al. 2002). By comparison, Beijing and other large Chinese cities often have average daily concentrations well in excess of (Chen et al. 2011).

Some evidence, indeed, suggests the shape of the C-R association might be nonlinear across the large range of ambient PM concentrations experienced in Chinese cities. For example, one multicity study of PM with aerodynamic diameter () and nonaccidental mortality in 16 Chinese cities found that estimated health effect sizes tended to be higher in cities with lower average pollution concentrations, consistent with a nonlinear C-R association that levels off at very high PM concentrations (Chen et al. 2011). Further, one national study found a leveling off in the C-R association for and total mortality at high concentrations in most Chinese regions (Chen et al. 2017).

The shape or magnitude of this C-R association might also differ between Beijing’s warm and cool seasons. This is in part because pollution sources differ between these two seasons. For example, between mid-November and mid-March, homes in Beijing are heated using coal-based heating, which adds a source of PM pollution, while spring wind patterns often bring desert dust to Beijing (Chan and Yao 2008). Further, air-conditioning and ventilation practices differ between seasons, so the indoor-to-outdoor exposure ratio likely shifts over the seasons, affecting personal exposure for a given outdoor concentration (Janssen et al. 2002). Some studies have found evidence of larger effects on nonaccidental, respiratory, and circulatory mortality in the cool season vs. warm season in Beijing (Li et al. 2013).

Here, we used a number of techniques to investigate both the shape and magnitude of the C-R associations between and several mortality outcomes in Beijing, China. Further, to determine how assumptions of linearity in the C-R association might influence HIA, we estimated excess deaths attributable to short-term exposure to in Beijing during our study period using C-R associations that are both nonlinear and season specific, compared with impact estimates based on C-R associations estimated assuming linearity.

Methods

Data

We obtained daily mortality counts from the Chinese Center for Disease Control and Prevention for all Beijing residents from 1 January 2009 to 31 December 2012. Cause of death was coded according to the International Classification of Disease, Revision 10 (ICD-10) (WHO 2016). We aggregated data to create daily values for three causes of death: nonaccidental death (ICD-10, A00–R99), deaths due to circulatory disease (ICD-10, I00–I99), and deaths due to respiratory disease (ICD-10, J00–J99). While individualized demographic characteristics were not available, Table S1 provides estimated demographic characteristics of the study population as a whole.

We collected daily data on concentrations from the Beijing Meteorological Bureau (http://www.bjmb.gov.cn/) for the Hai Dian Bao Lian station (Figure 1). We also collected hourly data on concentrations from the U.S. Embassy’s monitor in the Chaoyang District (Figure 1). The U.S. Embassy monitor provided hourly concentrations, which we aggregated to daily values by averaging across all hourly values for days with of hourly measurements available; the daily average was set to missing if of hourly measurements were missing for the day. To combine these data to create an exposure dataset for our main analysis, we imputed missing daily concentrations at either monitor, whenever possible, based on a model with data from the other monitor. For all days with observed monitoring data from both stations, we fit the following two linear regression models: where was the monitored concentrations of on day i from Hai Dian Bao Lian station (Monitor A), and was the observed concentrations of on the same day i from the U.S. Embassy’s monitor (Monitor B). For missing values in the dataset from Monitor A, we imputed using the predicted values from Equation 1. The same method was used to impute the missing values from Monitor B based on Equation 2. For each day, we averaged the values from Monitors A and B (either measured or imputed) to generate a daily value for main analysis. As a sensitivity analysis, we re-estimated main results using concentrations observed at each monitor individually. We collected weather data (mean temperature and relative humidity) from the Beijing Meteorological Bureau for a station in the Nanjiao District (Figure 1).

Figure 1.

Locations of two particulate matter with aerodynamic diameter () monitors and weather data monitor in Beijing. Light-gray lines show roads in Beijing.

Statistical Methods

We first explored whether the associations between daily concentration and mortality risk were consistent with a linear model or if this C-R association showed evidence of nonlinearity. For each of the three mortality outcomes considered (all nonaccidental, circulatory, and respiratory), we fit generalized linear models (GLM) with overdispersed Poisson distributions to model the C-R association between the logarithm of expected daily mortality and . To explore evidence of nonlinearity, we investigated three potential functions to include concentrations within the models (more details below). was incorporated as the average of same-day and previous-day concentrations (lag01) based on evidence from previous China-based studies that this lag period is appropriate for capturing acute risks (Kan et al. 2007; Chen et al. 2017). We conducted a sensitivity analysis of our main models using a longer lagged exposure of (lag02).

The general modeling framework followed the form: where is daily mortality count in Beijing on day t; is average of concentrations on days t and t-1 [measured for lag01 given previous findings that most of the association between and mortality in Chinese cities occurs within two days for both circulatory and respiratory outcomes (Chen et al. 2017)], incorporated into the model through one of three functional forms, [see below for the three forms of considered]; and is a matrix incorporating possible confounders, including a natural cubic spline for time with 7 degrees of freedom (df) per year to address potential confounding by seasonal or long-term trends at a time scale longer than approximately 2 months (Dominici et al. 2000), a natural cubic spline for daily mean temperature on day t with 3 df, a natural cubic spline for daily relative humidity on day t with 3 df, an indicator for day of the week, and an indicator for holidays (about 25 holidays each year, including Chinese New Year and National Day). Degrees of freedom for these natural spline functions were determined based on Quasi-Akaike’s Information Criterion (Q-AIC) and on controls used in previous air pollution epidemiologic studies (Schwartz et al. 2001).

While relative humidity provides a measure of the degree of saturation of air, other ways of incorporating air moisture into models might be more appropriate in certain scenarios (Davis, McGregor, and Enfield 2016). We therefore conducted a sensitivity analysis of our main models, testing alternatives to relative humidity used in the main analysis. First, we reran the main results using dew point temperature as a metric of air moisture. Since dew point temperature is highly correlated with air temperature, we adjusted dew point temperature for air temperature, using the residuals of a linear regression of dew point temperature on air temperature (Dominici et al. 2006; Samet et al. 2000). Second, instead of adjusting for air temperature and air moisture through separate model terms, we adjusted for heat index, a measure that incorporates air temperature and air moisture, using a natural cubic spline function with 3 df (as with the separate air temperature and moisture terms in the main model). To calculate the heat index, we used an open-source R software (version 3.4.2; R Development Core Team) implementation (Anderson et al. 2013) of the U.S. National Weather Service's heat index algorithm (http://www.wpc.ncep.noaa.gov/html/heatindex.shtml).

As a previous study found that lagged temperature can confound –mortality associations (Kim, Bell, and Lee 2017), we checked that results were robust to including control for longer lags of temperature. For this analysis, we included temperature in the model using a distributed lag function up to lag 21, with the lag dimension modeled by a cubic spline with 2 df and knots at log intervals and the temperature dimension modeled by a cubic spline with 3 df (as for the main model).

To explore evidence of nonlinearity in the shape of the –mortality association, we explored three potential functional forms for including within this modeling framework:

Linear function: , which assumes a linear association between and the logarithm of expected daily mortality (the functional form assumed in most studies of and acute mortality risk);

Piecewise-linear function: , where is a breakpoint allowing a change in slope in the association between and the logarithm of expected daily mortality, set a priori at [approximately the 75th () percentile of the year-round concentrations in Beijing; Figure 2]; and

Figure 2.

Frequency distribution of particulate matter with aerodynamic diameter () levels in all seasons (year round), cold season, and warm season from 1 January 2009 to 31 December 2012. The dotted line shows World Health Organization (WHO) guidelines for daily (), and the solid line shows Chinese 24-h urban standard for ().

Spline function: included as a natural cubic spline function of , with two interior knots at 75 and (Daniels et al. 2000; Samoli et al. 2005). These knot locations were approximately the 50th () and 75th () percentiles of the year-round concentrations in Beijing over the study period (Figure 2).

These functional forms move from a more constrained form (linear) to a less constrained form (spline function), allowing increasing flexibility in the shape of the C-R association between and mortality. For all models, we used convergence criteria that were more stringent than the software defaults to ensure proper model convergence (Dominici et al. 2002). We used Q-AIC to assess the performance of models incorporating these three functional forms of daily . We conducted diagnostic analyses of our main models to examine possible violations of model assumptions. Specifically, we calculated standardized deviance residuals to identify and investigate the potential influence of outliers (observations with standardized residuals ) and also used Cook's distance to assess evidence of influential outliers (Dunn and Smyth 2018).

We also fit the models described above separately using only warm season (April–September) or cold season (October–March) data. For these seasonal models, we used in the time spline included in the term of the models (Equation 3).

Finally, we investigated whether HIAs for short-term exposure to were sensitive to the assumed shape of the C-R association (i.e., which of the three functional forms was used to incorporate within Equation 3). We estimated the yearly attributable number of deaths for mortality outcomes (all nonaccidental, circulatory, and respiratory) as compared with three counterfactual thresholds: a background concentration of (), the World Health Organization (WHO) guideline for daily average () (WHO 2005), and the Chinese Ambient Air Quality Standard (GB 3095-2012, Class II 24-h standard: ), as Class II standards apply to urban areas. In each case, we calculated the number of deaths estimated given the observed daily concentrations of in Beijing over the study period compared with a scenario in which never exceeded the given counterfactual threshold (i.e., on any day with an observed concentration exceeding the threshold, the concentration was reset at that threshold, while observed concentrations were maintained for any day with an observed concentration below the threshold). In this context, a threshold is a level below which no deaths are attributed to short-term exposure to .

We estimated the attributable fraction (AF) of mortality associated with the observed daily concentrations compared with the counterfactual scenarios of , never exceeding a given counterfactual threshold () using the equation: where is the fraction of mortality attributable to the observed daily concentration (compared with the counterfactual scenario) on day t; is a vector of length one, two, or three (depending on the functional form used for ) of coefficients estimated from linear, piecewise-linear, and spline functions of from fitting the model given in Equation 3; is a matrix of one, two, or three columns (depending on the functional form used for ) describing the difference between the observed concentration () on day t and the counterfactual threshold being considered (). This function has a different form depending on the form used for in Equation 3, but is always a zero matrix on days when concentration is equal to or below the counterfactual threshold (i.e., ). Hence, when the observed concentration is equal to or below the counterfactual level, . On days when concentration exceeds the threshold, for the linear model, , while for the piecewise-linear model and spline models, is a matrix with the differences between the model matrix basis values for vs. , based on the basis used for the function fit in Equation 3.

We then estimated the daily attributable number () of deaths on each study day from the estimated AF of deaths as (Gasparrini and Leone 2014):

Here, and are as defined for Equations 3 and 4, respectively. We summed over our study period and averaged these values to obtain an estimated attributable number of deaths by both season and year for each cause of death and under each format of the model fit with Equation 3 (linear, piecewise-linear, and spline functions of daily concentrations). To compute 95% empirical confidence intervals (CIs) for attributable death estimates, we used Monte Carlo simulations with 5,000 samples (Gasparrini and Leone 2014). All analyses were conducted in R (version 3.4.2; R Development Core Team).

Results

Over our study period, the daily concentration of in Beijing ranged from 13 to , with a mean of (Table 1). Average concentrations did not vary substantially by month or season (Figure 2; Figure S1), with concentrations generally lowest around transition dates between warm and cold seasons (March and September). However, the highest-concentration days were most common in the middle of the cold season, coinciding with the period when coal combustion is used for heating in Beijing (generally from 15 November to 15 March). Over 90% (90.5%) of daily observations for in our study period exceeded the WHO’s daily guideline of , and about half of study days (47.9%) also exceeded the Chinese Ambient Air Quality 24-h urban standards of . observations were missing on 142 d (9.7% of study days) at the Hai Dian Bao Lian Station and on 100 d (6.8% of study days) at the U.S. Embassy Monitor. Daily data were correlated between the two monitors (Figure S2; ). Only 5 d had missing data for both monitors, and so, after imputation, concentrations were available for 99.7% of study days. During our study period, daily average death counts were 193 for nonaccidental mortality, 97 for circulatory mortality, and 20 for respiratory mortality. Death counts were typically higher in the cold season compared with the warm season (Table 1).

Table 1

Summary statistics of daily mortality, particulate matter with aerodynamic diameter () concentrations, and meteorological conditions. The gives the number of nonmissing observations in the study dataset (out of the 1,461 d from 1 January 2009 to 31 December 2012) for each measurement.

	Year round		Cold season (October–March)		Warm season (April–September)
	n	Mean (range)	n	Mean (range)	n	Mean (range)
PM_2.5 (μg/m^3)	1,456a	84.9 (13, 388)	727	89.2 (13, 388)	729	81.0 (13, 269)
Temperature (ºC)	1,461	13.1 (−13, 35)	729	3.2 (−13, 23)	732	22.9 (7, 35)
Relative humidity (%)	1,461	50.5 (9, 97)	729	45.1 (12, 91)	732	56.0 (9, 97)
Daily mortality count
All nonaccidentalb	1,461	192.9 (122, 279)	729	210.1 (139, 279)	732	175.8 (122, 253)
Circulatoryb	1,461	96.7 (49, 164)	729	109.0 (69, 164)	732	84.4 (49, 161)
Respiratoryb	1,461	19.7 (6, 45)	729	22.2 (6, 45)	732	17.3 (6, 41)

The 5 d (difference between 1,461 and 1,456 days) with missing data for were days where data was missing for both monitors and so could not be imputed.

Cause of death was coded according to the International Classification of Disease, Revision 10 (ICD-10) (WHO 2016), including all nonaccidental mortality (ICD-10, A00–R99), circulatory mortality (ICD-10, I00–I99), and respiratory mortality (ICD-10, J00–J99).

Adverse associations with were observed for nonaccidental (Figure 3) and circulatory mortality (Figure 4), and the estimated associations suggested some evidence of nonlinearity. When the linear assumption for the C-R association was relaxed through use of the piecewise-linear or spline models, the association between concentration and circulatory mortality was low at levels (although with large uncertainty), then increased linearly across higher concentrations. For respiratory mortality (Figure 5), the C-R association with also showed evidence of nonlinearity when fit using the piecewise-linear or spline model, although with a different shape. When fit with a piecewise-linear model, the observed association between concentration and respiratory mortality was positive at concentrations up to but then flattened out at high concentrations. This pattern was also evident when fit using a spline model, although in this case, there was also evidence that the C-R association then decreased once concentrations exceeded (Figure 5). However, the C-R associations for respiratory mortality had wide CIs, possibly due to the lower rates of this cause of mortality, making it difficult to determine the association’s shape from the most flexible (spline) model, especially at exceptionally high but infrequently observed concentrations of . Based on the Q-AIC for each model, the more flexible models (piecewise-linear and spline models) had slightly better performance than the linear model for all three mortality outcomes (Table S2). There was moderate evidence of improved model performance for the most flexible (spline) model vs. the linear model when considering circulatory and nonaccidental mortality (difference in Q-AIC of 8.74 and 2.94, respectively; Table S2). For respiratory mortality, however, the difference in Q-AIC was lower between the spline and linear models (difference in Q-AIC of 1.72; Table S2), and the preferred model for this less common mortality outcome (compared with circulatory mortality) was the less flexible but still nonlinear piecewise-linear model (difference in Q-AIC of 8.15 compared with the linear model). For all outcomes, C-R associations in the cold season were more similar to year-round associations than those observed in the warm season (Figures 3–5); in fact, we observed little evidence of an association between concentration and mortality risk at any concentration during the warm season.

Figure 3.

Concentration–response curves for nonaccidental mortality fit using a linear model (left), piecewise-linear (middle), and spline model (right) for lag01 days for particulate matter with aerodynamic diameter () using year-round (top), cold-season (middle), and warm-season data (bottom). Dotted lines show 95% confidence intervals for each model fit. Tick marks on the bottoms of the graphs show the distribution of concentrations across study days. Relative risks (RRs) were estimated using generalized linear models with overdispersed Poisson distributions, controlling for natural cubic splines for time trends [], daily mean temperature (3 df), and daily relative humidity (3 df), and indicator terms for day of the week and holidays.

Figure 4.

Concentration–response curves for circulatory mortality fit using a linear model (left), piecewise-linear model (middle), and spline model (right) for lag01 days for particulate matter with aerodynamic diameter () using year-round (top), cold-season (middle), and warm-season data (bottom). Dotted lines show 95% confidence intervals for each model fit. Tick marks on the bottoms of the graphs show the distribution of concentrations across study days. Relative risks were estimated using generalized linear models with overdispersed Poisson distributions, controlling for natural cubic splines for time trends [], daily mean temperature (3 df), and daily relative humidity (3 df), and indicator terms for day of the week and holidays.

Figure 5.

Concentration–response curves for respiratory mortality fit using a linear model (left), piecewise-linear model (middle), and spline model (right) for lag01 days for particulate matter with aerodynamic diameter () using year-round (top), cold-season (middle), and warm-season data (bottom). Dotted lines show 95% confidence intervals for each model fit. Tick marks on the bottoms of the graphs show the distribution of concentrations across study days. Relative risks were estimated using generalized linear models with overdispersed Poisson distributions, controlling for natural cubic splines for time trends [], daily mean temperature (3 df), and daily relative humidity (3 df), and indicator terms for day of the week and holidays.

Figure 6 demonstrates the implications of this observed nonlinearity on estimates of the percent increase in mortality per increase in lag01 at different concentrations of . With the linear model, a single association is estimated for any increase across the range of observed concentrations. Based on this linear model, we estimated an increase of 0.19% (95% CI: 0.07, 0.31) for nonaccidental mortality, 0.29% (95% CI: 0.12, 0.46) for circulatory mortality, and 0.20% (95% CI: , 0.53) for respiratory mortality for a increase in lag01 days of . When the C-R association was fit using the piecewise-linear and spline models, we were able to estimate these effect estimates conditionally on concentration, as the association was allowed to be nonlinear. For nonaccidental and circulatory mortality, adverse associations were observed when concentrations were (using the piecewise-linear model) or within the range of 130 to (using the spline model). While adverse associations between and respiratory mortality were seen at lower concentrations of [e.g., 0.54% (95% CI: 0.01, 1.07) increase of mortality estimated for a increased from 70 to ], there was little evidence of association at higher concentrations [e.g., (95% CI: , 0.37) increase in mortality estimated when increased from 170 to ].

Figure 6.

Percentage difference (mean and 95% confidence intervals) in daily mortality per increase in particulate matter with aerodynamic diameter () at lag01, at different levels of . The shape of the point estimates refers to different models, with estimates from fitting the linear model shown with circles, the piecewise-linear model with triangles, and the spline model with squares. Results were estimated using generalized linear models with overdispersed Poisson distributions, controlling for natural cubic splines for time trends [], daily mean temperature (3 df), and daily relative humidity (3 df), and indicator terms for day of the week and holidays.

For nonaccidental and respiratory mortality, we found the standardized deviance residuals were approximately normally distributed with a mean of zero, with approximately constant variance (Figure S3, left and center panels). For circulatory mortality, although the standardized residuals had a fairly constant variance and mean of 0 (Figure S4, left panel), there were a few observations with standardized residuals greater than 3 and a suggestion of some right skewedness for standardized residuals based on the Q-Q plot (Figure S4, left and middle panels). For each outcome and model, there were standardized residuals greater than 3, out of 4 y of daily data included in the study. According to the Cook’s distance measures (Figure S3–S5, right panels), there was little evidence that any of these outliers were influential observations in our modeling process for the three mortality outcomes, with all Cook’s distance measurements well below 1 (Dunn and Smyth 2018). As a sensitivity analysis, we estimated mortality risk associated with ambient using concentration measurements from each of the two monitors separately (Figure S6). While the specific central effect estimates changed, the observed trends in effect estimates based on data from the two monitors were similar at the lower concentrations, but differed at high concentrations. The monitor at the U.S. Embassy, located closer to the city center (Figure 1), demonstrated trends consistent with our primary analysis across all concentrations. In further sensitivity analysis, we found results were robust to using dewpoint temperature or heat index rather than relative humidity, to extending the lag of to lag02, and to controlling for a distributed lag function of temperature out to 21 d, in all cases showing similar association sizes and, for nonlinear models, C-R association shapes (Figures S7–S10).

Finally, we explored whether estimates of attributable deaths from short-term exposure were sensitive to whether the C-R association was fit assuming a linear association. We estimated yearly excess deaths attributable to short-term exposure in Beijing during our study period compared with scenarios in which never exceeded one of three counterfactual thresholds. For this analysis, we calculated estimated excess deaths using linear-, piecewise-linear-, and spline-based C-R associations fit using year-round data (Table 2).

Table 2

Estimated number and fraction of deaths attributable to particulate matter with aerodynamic diameter () in Beijing on average per year (based on 2009–2012) compared with counterfactual scenarios in which no study day exceeds either a) a background concentration of , (2) the World Health Organization (WHO) daily guideline of , or (3) the Chinese Ambient Air Quality Standard (CAAQS) 24-hour standard of . Estimated excess deaths are calculated using concentration–response functions generated using linear, piecewise-linear, or spline models fit to year-round data with overdispersed Poisson distributions, controlling for natural cubic splines for time trends [], daily mean temperature (3 df), and daily relative humidity (3 df), and indicator terms for day of the week and holidays. Numbers in parentheses show 95% confidence intervals for either attributable number or attributable fraction.

	Counterfactual PM_2.5 threshold
	7.5 μg/m^3		25 μg/m^3		75 μg/m^3
	Attributable number	Attributable fraction (%)	Attributable number	Attributable fraction (%)	Attributable number	Attributable fraction (%)
Nonaccidental mortality
Linear model	1,033 (361; 1,694)	1.47 (0.52, 2.39)	805 (290; 1,317)	1.14 (0.42, 1.87)	322 (110, 520)	0.46 (0.16, 0.75)
Piecewise-linear model	379 (−647, 1391)	0.38 (−1.16, 1.88)	347 (−400, 1073)	0.39 (−0.73, 1.48)	278 (67, 487)	0.42 (0.13, 0.71)
Spline model	259 (−1,244; 1645)	0.37 (−1.61, 2.38)	283 (−635; 1,208)	0.40 (−0.81, 1.70)	272 (26, 522)	0.39 (0.03, 0.74)
Circulatory mortality
Linear model	795 (320, 1237)	2.25 (0.95, 3.52)	621 (262, 973)	1.76 (0.78, 2.75)	251 (108, 394)	0.71 (0.31, 1.11)
Piecewise-linear model	127 (−568, 833)	0.10 (−2.09, 2.23)	152 (−382, 648)	0.27 (−1.33, 1.80)	205 (55, 350)	0.63 (0.22, 1.03)
Spline model	−60 (−1,062; 902)	−0.17 (−3.05, 2.54)	56 (−573, 690)	0.16 (−1.69, 1.88)	212 (39, 387)	0.60 (0.11, 1.10)
Respiratory mortality
Linear model	111 (−72, 291)	1.54 (−1.04, 4.11)	86 (−58, 225)	1.20 (−0.80, 3.15)	34 (−23, 92)	0.48 (−0.32, 1.25)
Piecewise-linear model	240 (−41, 509)	3.54 (−0.65, 7.42)	177 (−32, 379)	2.59 (−0.41, 5.48)	43 (−17, 102)	0.56 (−0.27, 1.34)
Spline model	270 (−98, 639)	3.74 (−1.71, 9.02)	200 (−46, 437)	2.77 (−0.68, 6.02)	54 (−16, 123)	0.75 (−0.20, 1.70)

Note: Since some models estimated relative risks lower than one within some concentration ranges of , some scenarios resulted in an estimated deficit of deaths (negative values in table cells) rather than estimated excess deaths compared with the counterfactual.

Excess death estimates in the HIA were, in some cases, very robust to whether the C-R association was fit using a linear assumption, even though nonlinearity was evident in the observed associations. For example, we estimated similar excess nonaccidental deaths attributable to using a threshold of (i.e., estimating deaths for for the three model forms considered (322, 278, and 272 deaths/year for linear, piecewise-linear, and spline forms, respectively; Table 2); this is because the estimated shapes for the C-R association for nonaccidental mortality were fairly linear with similar slopes at levels for all three model forms (Figure 3).

In other cases, however, estimates of attributable deaths were sensitive to whether the C-R association was fit using a linear assumption. For example, when using a threshold of , the excess death estimates were different across C-R associations estimated using different model forms (1,033; 379; and 259 deaths/year for linear, piecewise-linear, and spline forms, respectively; Table 2) because of the observed relatively low risk estimate at low levels captured by the piecewise-linear and spline models as compared with the linear model. For respiratory mortality, although nonlinear model forms estimated that the association leveled-off at high concentrations of (Figure 5), we found that assuming linearity in the association resulted in lower estimates of attributable risk. This is because levels in about half of study days were (Figure 2), and nonlinear model forms estimated steeper PM–respiratory mortality associations at concentrations than the linear model (Figure 5). In all cases, the CIs for excess death estimates were wider when the C-R association was fit using a nonlinear assumption rather than a linear one.

We found that most of the excess deaths were during the cold season (Table S3). For example, when using the Chinese 24-h standard of as a counterfactual and a C-R association based on linear models fit separately to seasonal data, we estimated 322 excess nonaccidental annual deaths year round (Table 2), combining 214 excess deaths in the cold season and 108 excess deaths in the warm season (Table S3).

Discussion

In this study, we investigated both the shape and magnitude of the C-R associations between and three mortality outcomes in Beijing, China. We found evidence of nonlinearity in these C-R associations, but with differing shapes depending on the mortality outcome. Within the cold season, the shape of the C-R association was consistent with the shape estimated from year-round data, although the effect size was larger. Conversely, in the warm season, we observed almost no association between mortality risk and ambient concentrations. Based on these estimated C-R associations, short-term exposure to in Beijing is associated with a large number of attributable deaths. However, while in some cases, estimates of attributable deaths associated with acute exposures were robust to assumptions of linearity in the C-R association, in other cases, these impact estimates were extremely sensitive to this assumption, which highlights the importance of understanding potential nonlinearity in the C-R association between PM and mortality at the high PM concentrations in developing megacities like Beijing.

We found evidence of nonlinearity in the –mortality association in Beijing and, interestingly, that the nonlinear C-R association shapes were different for respiratory mortality compared with circulatory and all nonaccidental mortality. A number of factors in Chinese cities could lead to nonlinearity in this C-R association. First, people might vary their behaviors [like wearing a face mask on haze days (Langrish et al. 2012)] depending on pollution levels, especially at the very high concentrations sometimes experienced in Chinese cities, and such behavior changes could influence the observed shape of the C-R association at different concentrations. While face masks commonly used in Beijing do not completely filter , they do prevent penetration of some (Cherrie et al. 2018), and studies have found that both simple and high-efficiency face masks can improve cardiovascular markers, especially heart rate variability, in volunteers exposed to PM air pollution in Chinese cities (Langrish et al. 2012; Shi et al. 2017; Langrish et al. 2009). Other behavior changes may also help explain the observed nonlinearity in the C-R associations. For example, if some people avoid outdoor activities and decrease indoor ventilation during extreme air pollution, the association between personal exposure and the air pollution levels observed at ambient monitors may decrease compared with more moderate days. Future research could further investigate the use, effectiveness, and health benefits of face mask use and other behavior changes during high air pollution days in Chinese megacities to help clarify how these factors influence population-level estimates of C-R associations.

Another possible explanation for both the observed nonlinearity in C-R associations and the difference in shapes between associations for circulatory vs. respiratory mortality is that the composition of can vary between higher and lower pollution days, and the different chemical species and mixtures that comprise can have different impacts on mortality risk (Cao et al. 2012; Huang et al. 2012; Franklin et al. 2008). In Beijing, on days of high () to severe () pollution, concentrations are dominated by secondary inorganic species—sulfate, nitrate, ammonium—due to chemical production and meteorological conditions (Ma et al. 2017). Some studies have suggested that these secondary inorganic species of can have much larger effects on cardiovascular mortality than other PM components (Huang et al. 2012; Son et al. 2012), which may help explain why we found evidence of a strong association between and circulatory mortality risk in the range of days with high concentrations and less evidence of an association for days with lower concentrations (Figure 4). In terms of the effects of PM components on respiratory mortality, evidence is limited and inconsistent, although some studies have found substantially increased risks on respiratory mortality from magnesium, chlorine, and sulfur (Dai et al. 2014; Son et al. 2012; Huang et al. 2012). In Beijing, magnesium is mainly attributed to windblown soil dust, and chlorine can be considered as a tracer for coal combustion in noncoastal cities (Yu et al. 2013). By quantifying source contributions of during the winter in Beijing, one study found that the contribution of soil dust was 20% higher on nonhaze days vs. haze days (Shang et al. 2018), while coal combustion was a steady source of wintertime pollution in Beijing given its use for home heating. This could help support the –respiratory mortality association at lower concentrations observed in this study.

When fitting C-R associations using season-specific data, we found further evidence that supports the hypothesis that the nonlinearity we observed in these C-R associations may be driven by differences in composition. While earlier studies found the total mass of varied strongly across seasons in Beijing (Zheng et al. 2005; He et al. 2001), more recent research indicates that current total concentrations in Beijing generally have smaller seasonal variation than earlier findings, and we similarly found that the average concentrations varied little by season during our study period (Figure S3; Table 1). However, the composition of can vary substantially by season as a result of varying emission sources, meteorological conditions, and photochemical activities (Zhang et al. 2013; Yu et al. 2013). During the warm season, windblown soil and resuspended road dust are major contributors to ; for example, Zheng et al. (2005) reported dust alone constituted 36% of mass in Beijing in April, while the annual average was 20%. Although some components of dust (e.g., aluminum, silicon, and elemental carbon) are associated with mortality risks (Dai et al. 2014; Huang et al. 2012), it is likely that in Beijing, the frequent precipitation in summer and strong wind in spring (Yu et al. 2013) could wash out suspended dust, resulting in less toxic PM in the warm season. The cold season, on the other hand, coincides with the harvest season (late September to mid-November) and the heating season (mid-November to mid-March) in the northern areas of China, when the contribution from biomass aerosol and fossil fuel combustion to sharply increases from an annual average of 7 to 24.5% in January (Zheng et al. 2005). A study of another northern city of China found that several species—sulfur, elemental carbon, potassium, and chlorine—are linked with biomass burning and coal combustion (Ryu et al. 2007; Zheng et al. 2005), all species that have been reported to be significantly associated with increased mortality risks (Huang et al. 2012; Li et al. 2015). These differences in PM composition between cold and warm seasons in Beijing may explain why we found different shapes for the –mortality associations in cold and warm seasons, with little evidence of an association throughout the range for warm season days, while the observed nonlinearity in the C-R association that we observed specifically within the cold season (Figures 3–5) may reflect the role of further, within-season variation in the composition of between days with higher vs. lower concentrations.

When the actual C-R function is nonlinear but not extremely different from linear, methods of determining whether it is nonlinear or linear (e.g., AIC) can struggle to detect the nonlinearity (Roberts and Martin 2006). Here, this may be particularly the case for investigating nonlinearity for respiratory mortality, as this outcome was much rarer on average than circulatory mortality (20 deaths/d vs. 97 deaths/d; Table 1). Also, exposure measurement error generally attenuates the observed nonlinearity in exposure–disease associations (Keogh et al. 2012). In communities with very high pollution levels, including Beijing and other Chinese megacities, future research could use a multicity study to improve power in detecting nonlinearity (Daniels et al. 2000; Schwartz et al. 2002; Samoli et al. 2005). However, if the nonlinearity is driven by differences in PM composition across the range of observed PM concentrations, the shapes of C-R associations could vary substantially across different Chinese cities, for which PM composition patterns can differ, preventing the ability to meaningfully pool multicity data to better estimate a nonlinear C-R association that is relevant throughout Chinese megacities. Further, the effect estimates could differ by population. Indeed, while the nonlinear patterns we found here were consistent with nonlinear C-R association shapes found in previous studies of Beijing (Li et al. 2013), noticeable spatial heterogeneity exists regarding the shape of –mortality curves in China. A study of Shanghai, where the levels were lower than Beijing (daily mean is during their study period), found the shape was nonlinear for total and cardiovascular mortality, and almost linear for respiratory mortality (Kan et al. 2007). Cao et al. (2012) reported linear associations between and total, cardiovascular, and respiratory mortality in Xi’an, a heavily polluted city with daily mean concentration of , although their curves show some evidence of nonlinearity for total and cardiovascular mortality, while a study of the association between and mortality in Wuhan, China, found some evidence, from generalized cross validation, of nonlinearity in the association between and respiratory mortality, but little evidence of nonlinearity for nonaccidental or cardiovascular mortality (Qian et al. 2007). A recent national study of 272 Chinese cities also reported significant regional heterogeneity in the shape of C-R associations for and mortality (Chen et al. 2017). This would suggest that -associated health effects may not be sufficiently captured by models based on PM mass concentrations only, especially within Chinese megacities. If differing composition is contributing to the observed nonlinearity in these C-R associations, future studies that estimate the association between specific chemical species and mixtures of and mortality would contribute substantially in helping to understand how PM pollution affects human health in Chinese megacities and beyond.

Compared with studies conducted in Western communities (e.g., Zanobetti and Schwartz 2009; Franklin et al. 2007), previous China-based studies have found lower effect estimates for the association between short-term exposure to and mortality outcomes when assuming a linear C-R association (Aunan and Pan 2004). When similarly assuming a linear C-R association, we also estimated associations between and mortality outcomes that were lower than those found in Western studies. The difference in effect size assuming a linear C-R association between China- and Western-based studies may be related to different population distribution and potential exposure misclassification. The residents of Beijing have a younger age distribution (only 9–10% aged ; Table S1) compared with many Western cities, and the association between and mortality can be much stronger among older than younger adults (Yang et al. 2013). –mortality associations also differ across specific causes of circulatory and respiratory death (Chen et al. 2017), so differences in underlying health and susceptibility of Beijing residents compared with Western populations may also contribute to observed differences in the size of the C-R associations. Also, Beijing residents may be more likely to take personal interventions (e.g., air purifier, face mask) to reduce personal exposure, which could potentially lead to exposure misclassification in population-based time-series studies, which could bias estimates toward the null.

Further, we found suggestive evidence that the lower effect estimates in China-based studies for respiratory mortality may partially result from assuming linearity in the C-R association. When using a nonlinear model, our effect estimates for respiratory mortality at lower concentrations of (consistent with levels in Western communities), while still not as high as estimates from Western studies, were more consistent with them than when assuming a linear C-R association than when incorporating nonlinearity. For example, Franklin et al. (2007) reported a 1.67% increase for respiratory death per increase in lag01 days of in 27 U.S. communities, where the average daily levels ranged from 9.3 to across all communities. In our study, the corresponding estimate was 0.2% when assuming a linear C-R association, but 0.6% at concentrations when using the piecewise-linear model and 0.56% when concentration changed from 50 to (the lowest increment we considered) based on the spline model (Figure 6). To the best of our knowledge, previous China-based studies have only reported effect estimates based on a linear assumption for the C-R association, concluding a lower association between concentration and respiratory mortality risk compared with Western studies (Kan et al. 2007; Chen et al. 2017).

In HIA, use of effect estimates from time-series studies can underestimate the number of excess premature deaths associated with air pollution because, while they capture immediate deaths triggered by air pollution, they do not capture premature deaths in which a person’s frailty was caused by longer-term exposure to air pollution (Künzli et al. 2001). While C-R associations derived from time-series studies, therefore, typically are not used to estimate impacts in terms of YLL for western HIAs, they can be a critical way to add local estimates of C-R associations into China-based HIAs. Some HIAs of air pollution in China only use effect estimates from cohort studies (mainly conducted in the United States) (Zhang et al. 2007), which raise challenges given differences in demographics, population susceptibility, and severity and composition of pollution (Kan and Chen 2004; Aunan and Pan 2004; Zhang et al. 2010; Hubbell et al. 2009; Kan et al. 2009). In particular, studies of the health effects of PM in China have typically resulted in lower estimates than similar studies in Western communities (Aunan and Pan 2004), and so the effect estimates from cohort studies based in Europe or the United States might overestimate expected impacts, especially at higher PM concentrations (Ostro 2004). Because of these limitations in using evidence from Western cohort studies in Chinese HIAs, a number of such assessments have used estimated C-R associations from Chinese time-series studies as a key input. For example, Matus et al. (2012) estimated YLL from acute exposure to PM in China using evidence from time-series studies, assuming 6 months of life lost for each premature death associated with acute exposure to . Others have estimated attributable risk based on time-series studies in China (Kan and Chen 2004; Wang and Mauzerall 2006), given the paucity of cohort studies in China, and to compare mortality risks from short-term and long-term exposure to PM (Evans et al. 2013).

We found that HIA for short-term exposure to that incorporate an estimated C-R association from time-series studies can be sensitive to model specification of the shape of the association; in Beijing, more flexible models estimated fewer attributable deaths for nonaccidental and circulatory mortality, but more attributable deaths for respiratory mortality, compared with the linear model. Based on these results, using a linear C-R for the PM–mortality association when the true association is nonlinear could substantially change estimates of acute health impacts (Roberts and Martin 2006). Thus, uncertainty resulting from the assumed model form could affect previous HIA of acute PM in China, as they generally have not incorporated nonlinearity in the C-R function when estimating attributable risk (Matus et al. 2012; Voorhees et al. 2014; Guo et al. 2013).

Click here for additional data file.

Click here for additional data file.