Summer Heat and Mortality in New York City: How Hot Is Too Hot?

Background To assess the public health risk of heat waves and to set criteria for alerts for excessive heat, various meteorologic metrics and models are used in different jurisdictions, generally without systematic comparisons of alternatives. We report such an analysis for New York City that compared maximum heat index with alternative metrics in models to predict daily variation in warm-season natural-cause mortality from 1997 through 2006. Materials and methods We used Poisson time-series generalized linear models and generalized additive models to estimate weather–mortality relationships using various metrics, lag and averaging times, and functional forms and compared model fit. Results A model that included cubic functions of maximum heat index on the same and each of the previous 3 days provided the best fit, better than models using maximum, minimum, or average temperature, or spatial synoptic classification (SSC) of weather type. We found that goodness of fit and maximum heat index–mortality functions were similar using parametric and nonparametric models. Same-day maximum heat index was linearly related to mortality risk across its range. The slopes at lags of 1, 2, and 3 days were flat across moderate values but increased sharply between maximum heat index of 95°F and 100°F (35–38°C). SSC or other meteorologic variables added to the maximum heat index model moderately improved goodness of fit, with slightly attenuated maximum heat index–mortality functions. Conclusions In New York City, maximum heat index performed similarly to alternative and more complex metrics in estimating mortality risk during hot weather. The linear relationship supports issuing heat alerts in New York City when the heat index is forecast to exceed approximately 95–100°F. Periodic city-specific analyses using recent data are recommended to evaluate public health risks from extreme heat.

volume 118 | number 1 | January 2010 • Environmental Health Perspectives Research It has long been known that elevated summer time temperatures are associated with increased daily mortality (Basu and Samet 2002). This relationship is most dramatic during heat waves (several days of extreme heat) with the occurrence of heat stroke deaths [(Centers for Disease Control and Prevention (CDC) 2005; New York City Department of Health and Mental Hygiene (NYCDOHMH) 2006; Semenza et al. 1996] and excess natural-cause mortality (Davis et al. 2003). Many communities issue heat advisories before these periods of excessive heat to provide advance warning of dangerously hot weather and to allow for a timely response by local agencies.
Most cities in the United States use National Weather Service (NWS) excessive heat alerts (advisories, watches, and warnings) based on the 24-to 48-hr forecast maximum heat index, a metric intended to reflect perceived temperature based on the measured temperature and relative humidity. An advisory is issued when the heat index is forecast to exceed 100°F (35°C) and 105°F (38°C) for northern and southern locations, respectively (NWS 2008a). Increasingly, other cities in the United States, such as Philadelphia (Kalkstein et al. 1996), and at least two dozen other cities worldwide (Sheridan and Kalkstein 2004), are using the Heat Health Watch and Warning System (HHWWS) as an alternative for heat alerts. The HHWWS uses daily forecasts of the dominant local weather pattern-the spatial synoptic classification (SSC)-to calculate predicted excess deaths based on location-specific historical mortality data (Sheridan 2002).
One criticism of using heat index-based alerts, rather than a HHWWS, is that the criteria are not specific for the health risks of the local population, particularly in urban areas (Tew et al. 2004). The differences in the heatmortality relationship between communities, with northern cities experiencing a greater mortality impact from higher summertime temperatures than southern cities (Curriero et al. 2002;McMichael et al. 2008), may be more nuanced than the difference in latitude. For example, the urban heat-island effect causes urban centers to experience higher daytime and nighttime temperatures than surrounding areas (Gaffin et al. 2008), and microclimate effects within cities (Harlan et al. 2006) may also play a role. These factors make some populations in urban areas more vulnerable to heat-related morbidity and mortality (Buechley et al. 1972;Smoyer 1998). Regardless of the weather metric used to model risk, the use of mortality data from previous decades (Ebi et al. 2004) may not accurately reflect current risk because of declining overall mortality (Kung et al. 2008) and changes in population vulnerability to heat (Carson et al. 2006;McMichael et al. 2006).
As part of an effort to evaluate the New York City criteria for issuing heat advisories (NWS 2008b) and activating emergency responses to extreme heat, we evaluated heat index, temperature, SSC categories, and additional meteorologic factors as predictors of natural-cause mortality during summertime weather conditions. The analyses used 10 years of daily mortality data on New York City residents for the months of May through September from 1997 through 2006 and compared the goodness of fit among models with various parametric and nonparametric functions and lag structures relating weather to mortality. temperature was highly correlated among the three weather stations (R > 0.93).
The weather metrics we evaluated included minimum, maximum, and average (mean of minimum and maximum) temperature (°F), sea level pressure (mbar), wind speed (miles per hour), and precipitation (any vs. none). An exponential function of temperature above a threshold {0.1 × exp [0.2 × (T max -90)]} where T max is maximum temperature was also considered given its previous use in examining the relationship of temperature and mortality in New York City (Buechley et al. 1972;Marmor 1975). Heat index was calculated using ambient temperature (°F) and relative humidity (%) for ambient temperature of ≥ 80°F and relative humidity of ≥ 40% (Steadman 1979). For this analysis, the daily maximum heat index is the larger value of maximum temperature or maximum heat index on each day, which was the best measure of the highest perceived daily temperature. Daily deviation from maximum temperature was calculated using as a reference period the average maximum temperature at LaGuardia Airport from 1971 through 2000 (NWS 2009).
Daily spatial synoptic classification (SSC) data were obtained for LaGuardia Airport during the same time period (Sheridan 2008). SSC is a classification system that uses temperature, dew point, wind direction, wind speed, and cloud cover measured four times daily to categorize surface weather conditions on each day into six main types and two subtypes (Sheridan 2002). Our analyses focused on seven weather types that occurred during the warm season in New York City: dry polar, dry moderate, dry tropical, moist polar, moist moderate, moist tropical, moist tropical plus, and transitions.
Mortality data. Data on every death of a New York City resident that occurred in the five boroughs (The Bronx, Brooklyn, Manhattan, Queens, and Staten Island) for the period 1997 through 2006 were obtained from the NYCDOHMH Office of Vital Statistics. Data included underlying and contributing cause of death, age, sex, race, place of birth, place of death, and census tract of residence. The NYCDOHMH Institutional Review Board approved this study.
To classify cause of death, we used the International Classification of Diseases, 10th revision (ICD-10; WHO 2007) for the 1999 through 2006 data; the ICD-9 (WHO 2009) was used for the years 1997 to 1998. The outcome of interest was natural-cause mortality, which included deaths with underlying cause of death exclusive of external causes (ICD-10, A00-R99; ICD-9, 001-799).
Analytic methods. We developed Poisson time-series regression models to relate daily natural-cause mortality to daily weather conditions during the warm season (May through September). Because the focus of the study was heat-related mortality, the seasonal restriction simplified the control for seasonal trends by eliminating the larger wintertime mortality variation associated, in part, with influenza and other communicable diseases (Hajat et al. 2002).
To assess the robustness of our findings, two types of analytic models-parametric generalized linear models and nonparametric generalized additive models were analyzed in parallel. Using parametric models, we investigated models with deterministic functions (i.e., linear, quadratic, cubic, sine, and cosine) and their interactions with year. The parametric models were developed using PROC GENMOD (version 9.1.3; SAS Institute Inc., Cary, NC). Models with nonparametric smoothing terms were developed to estimate the shape of potential nonlinear relationships between temperature and mortality and to compare them with the functional forms obtained from parametric models. In the nonparametric models, the degree of smoothness of model terms was estimated as part of fitting. Penalized splines were used to assess the optimum effective degrees of freedom for temporal and meteorologic terms (Wood 2004), as implemented in R statistical software (version 2.7.1; R Core Development Group, Vienna, Austria).
Our model-building process involved first adjusting for long-and short-term temporal trends only. In parametric models, to control for long-term and seasonal temporal trends, we included linear and quadratic terms for the day of year and their interaction terms with year, plus sine and cosine terms (1 period per year) to control for seasonality. In nonparametric models, the number of degrees of freedom used to control for long-term and seasonal temporal trends was based on the number that yielded the minimum of the sum of the absolute values of partial autocorrelation. This method allows control for a mixture of unmeasured confounders, such as infectious diseases and influenza epidemics, without penalizing for the amount of "wiggle" in the fitted temporal trend. Such methods have been applied to adjust for temporal trends in estimating air pollution effects (Peng et al. 2006;Touloumi et al. 2006). We tested models using between 2 and 10 degrees of freedom per season and found that 5 degrees of freedom per season gave the optimum result. Indicator variables for day of week, holidays, the World Trade Center attacks (11 September 2001), and the Northeast blackout (15 August 2003) were added to control for short-term temporal trends. The fitted series of the temporal nonparametric model was similar to that the temporal parametric model.
Next, we developed meteorology models, built upon the base temporal models, for each candidate meteorologic metric intended to reflect daily population heat exposure. The following metrics were assessed: minimum, average, and maximum temperature, maximum heat index, the exponential temperature transformation, SSC (categorical), and deviation from normal maximum temperature (categorical). For each metric, lagged functions of up to 1 week were assessed to identify those that improved model fit, with the constraint of including all intermediate lags of the shortest and longest lags deemed important (Armstrong 2006). We also evaluated the average of multiday lags but found that the models with multiple single-day lags often resulted in at least as good a fit as the multiday averages or combinations of single days and the average of multiday lags. For example, a model with lags of 0 through 3 days from the maximum heat index yielded approximately the same deviance explained as a model with the same day and the average of 1-, 2-, and 3-day lags. Thus, models with single-day lags were preferred to allow interpretation of results for each day and the ability to estimate excess mortality risk for consecutive hot days. In parametric models, linear, quadratic, and cubic terms for each continuous meteorologic variable were included together to allow for a nonlinear relationship of temperature and mortality (Armstrong 2006). In nonparametric models, penalized splines were used to simultaneously fit each of the individual lagged meteorologic variables. This resulted in degrees of freedom of between 1 and approximately 8 over the variable range, depending on the lag. Thus, although the nonlinear shape of the temperature-mortality relationships estimated at each lag day could be unstable when multiple-day lags were included simultaneously (as with the instability of individual estimates in an unconstrained distributed lag model), the results suggest consistency of the shape estimated by these two models.
We compared single-metric meteorology models for fit and explanatory ability (Armstrong 2006). Model fit was evaluated using percent deviance explained, first-order residual autocorrelation, and correlation of raw and predicted values on days with heat index ≥ 90°F (to assess the model fit during hotter weather).
In sensitivity analyses, we evaluated how much the fit of each model based on a single heat-related meteorologic metric was improved by adding additional meteorologic variables, including average sea level pressure (linear), average wind speed (linear), and precipitation (any precipitation, no precipitation), as well as SSC categories.
We also assessed the potential impact of seasonal acclimatization on mortality. We directly examined potential effect modification of the relation of temperature and mortality by testing interaction terms with day volume 118 | number 1 | January 2010 • Environmental Health Perspectives of year and month. Indirect assessment was conducted using the model with deviation from normal maximum temperature, because a given maximum temperature represents a greater deviation from normal when it occurs earlier in the season.

Results
On average, 136 New Yorkers died from natural causes each day during the warm season (May through September) from 1997 through 2006 (Table 1). For the 10-year study period, a gradual downward trend was observed in the average number of daily natural-cause deaths (144 in 1997 to 129 in 2006). In the warm season, average daily mortality was highest at the beginning of May and reached its nadir in mid-August.
During the study period, the average daily temperature at LaGuardia Airport was 71°F, the average maximum temperature was 80°F, and the average maximum heat index was 80°F in the warm season (Table 1). Of the 5 months, July (average daily temperature, 77°F; average maximum heat index, 87°F) and August (average daily temperature, 76°F; average maximum heat index, 86°F) were the warmest months. Temperatures peaked at the end of July.
The different meteorologic metrics that represent daily population heat exposure were moderately to strongly correlated (Table 1). Maximum temperature and maximum heat index were very strongly correlated (R = 0.986). Average barometric pressure and average wind speed were negatively correlated with temperature metrics, but the relation was weak.
During the warm season, 7% of days were included in the SSC category dry tropical, 22% were moist tropical, and 8% were moist tropical plus. The proportion of days classified into each of these categories increased as the maximum heat index increased. At heat indices of ≥ 100°F (n = 39), dry tropical and moist tropical plus described 44% and 49% of days, respectively. The mean maximum heat index was 91°F on dry tropical days, 85°F on moist tropical days, and 91°F on moist tropical plus days.
The parametric model with a single meteorologic factor that had the best overall fit considering all three criteria-greatest proportion of deviance explained, smallest first-order autocorrelation, and largest correlation of raw and predicted values on hot days-included linear, quadratic, and cubic terms for maximum heat index on the same day and previous 3 days (lags 0-3; Table 2). This model fit the data moderately better than did the four other models with a single meteorologic metric, including cubic functions of maximum temperature (lags 0-3), minimum temperature (lags 0-3), average temperature, and SSC (lags 0-3; Table 2). The models with the exponential temperature transformation (lags 0-1) and cate gories of deviation of normal maximum temperature (lags 0-1) fit the data least well among these single meteorologic metric models. However, the performance difference between the best-and worst-fitting model was not large as measured by the proportion of daily mortality variation (i.e., deviance) explained: 27.9% for the SSC model and 24.7% for the deviation from maximum temperature model.
Results of nonparametric models were similar ( Table 2). The models with maximum heat  index or with average, minimum, or maximum temperature on the same day and previous 3 days (lags 0-3) performed equally well. The model that used penalized splines to control for maximum heat index suggested 1 degree of freedom (i.e., linear) for the same day and approximately 5, 8, and 7 degrees of freedom for lags of 1, 2, and 3 days, respectively.
The fitted relationships of daily maximum heat index to mortality were similar for the parametric and nonparametric models ( Figure 1). On the same day, the relationship appeared to be linear across the range of observed heat index. The relations of mortality to heat index on the 3 previous days, however, were relatively flat across a range of moderate seasonal temperatures but increased sharply at a heat index between approximately 95°F and 100°F. The estimated risk ratios for same-day SSC categories were consistent with the relationship in mean maximum heat index among SSC categories, whereas the lagged effects of SSC were less consistent ( Figure 2).
In sensitivity analyses, we assessed how much the fit of the maximum heat index model could be improved with the addition of other meteorologic variables. The first sensitivity analysis added SSC categories (lags 0-3) to the final heat index model. The second added cubic terms for minimum temperature (lags 0-3), mean sea level pressure (lags 0-1), mean wind speed (lags 0-1), and precipitation (lags 0-3). Terms were added sequentially and determined to be statistically significant in the model. Each of these alternative models fit the data better than the heat index only model (Table 2). Although the risk ratios for Figure 1. Risk ratios (and 95% CIs) for the association of daily natural-cause mortality by maximum heat index on the same day (lag 0) and 3 previous days (lags 1-3, respectively) using parametric generalized linear models (A) and nonparametric generalized additive models (B), New York City, 1997 through 2006. Each risk ratio function is adjusted for temporal covariates for year, season, and day of week and for the other lagged heat index terms.  Figure 3). We found little evidence for effect modification of temperature by time of year during this time period. Using the maximum heat index parametric model as the base model, all interaction terms between maximum heat index functions and within-season temporal trends, including day of year and month, were not significant. When stratifying the analysis by early warm season (May through July) to late warm season (August through September), the shape of the risk ratio response function did not differ significantly (data not shown). The results were essentially the same when we used data from 1990 through 2006 in a sensitivity analysis capture years with more hot weather in May.

Discussion
In our analysis, several meteorologic metrics used to describe the extent of population heat exposure performed similarly as predictors of increased natural-cause mortality in New York City during the months of May to September from 1997 to 2006. We observed this association with same-day mortality across the range of seasonal temperatures, whereas days with maximum heat index of approximately 95-100°F or greater were associated with higher mortality on the following 3 days.
We found that using maximum heat index provided moderately better model fit in parametric generalized linear models compared with models using average, minimum, or maxi mum temperature (Hajat et al. 2002;Hajat et al. 2006;Huynen et al. 2001). Similar to a study assessing the impact of weather on the association of particulate matter and mortality, we found that models using smoothed functions of continuous temperature metrics have better fit than models using synoptic weather categories only (Samet et al. 1998). Additionally, the smoothed function approach, in both parametric and nonparametric models, described the relationship of mortality with successively higher heat indices more fully than simpler transformations of temperature (McMichael et al. 2008). It should be noted, however, that all of these models used highly correlated measures of temperature and heat index; the fit and predictive ability of the models were similar.
Our finding that the temperaturemortality relationship persisted for several days is consistent with previous research (Hajat et al. 2002;Hajat et al. 2006;Huynen et al. 2001). The magnitude of this lagged relationship increased sharply at heat indices between approximately 95°F and 100°F. Other studies of cities with climates similar to that of New York City have described an increase in mortality risk at temperatures within this range (Michelozzi et al. 2006). The lagged effect of extreme heat on mortality is also consistent with the observed lethality of heat waves with extremely high temperatures lasting several days (Tan et al. 2007).
We observed a nonlinear relationship of heat index and mortality at lagged days, but not on the same day in this data set. In contrast, two recent multicity studies found essentially linear relationships between temperature and mortality during the warm season (Anderson and Bell 2009;Zanobetti and Schwartz 2008). Both studies, however, limited their models to same-day temperatures (Zanobetti and Schwartz 2008) or 2-day averages of the same and previous day (Anderson and Bell 2009), for which we also found linear relationships in both parametric and nonparametric models. Those studies also noted that the temperature-mortality relationship differed by city, dependent on location and other community characteristics. Thus, although functional forms may vary from city to city, determining the shape of the temperature-mortality relationship needs to take into consideration lagged associations. If the temperature-mortality relationship were fully linear, weather criteria for heat alerts would be completely arbitrary. Our analysis of New York City data indicates that there is a linear component (same-day) and a nonlinear component (lagged) in the temperature-mortality relationship supporting issuing heat alerts when the heat index is forecast to exceed approximately 95°-100°F. At the same time, our analysis does not suggest a bright line below which no excess heat-related mortality occurs. Even if there was a statistical threshold for increased risk at the population level, individuals within the population may have different thresholds for experiencing life-threatening heat illness, depending on living conditions and health status.
Other studies have indicated that heat waves earlier in the season are more lethal (Hajat et al. 2002), possibly due to populations being less acclimatized to hot weather. Alternatively, if early heat waves result in "harvesting" or "mortality displacement," the mortality impact of heat waves later in the  summer is reduced among the now smaller pool of susceptible individuals. We found little evidence of a larger effect of early-season heat waves during our study period. There was no interaction between maximum heat index and month, nor did deviation from normal temperature fit the data as well as maximum heat index. This finding, however, could be a result of few unseasonably hot days occurring in May or June during this time period, although our sensitivity analysis including data from 1990-2006 was consistent with our main finding. Because of the infrequency of early-season heat waves, we do not feel our findings provide strong evidence either for or against a greater hazard of such events. Our analysis of heat-mortality relationships used both parametric and nonparametric modeling approaches. The final models identified by each of these approaches were remarkably similar, yielding risk ratio functions of similar form and magnitude of effects. Although the large numbers of degrees of freedom used for smoothing temperature metrics in the nonparametric model seem more than that needed to describe a biologically plausible relationship between heat index and mortality (and likely reflect the year-to-year variation in the relationship), the overall shape of these terms does not deviate substantially from either linear, hockey-stick, or U-shaped curves. Nonparametric modeling has the flexibility to identify and account for trends in continuous variables that are more difficult to do with parametric functions. This flexibility, however, can enable overfitting of observed variation. The consistency of our findings using two different approaches suggests that we have estimated a relatively robust relationship between hot weather and mortality in New York City.
We focused our analyses on the warm season only in order to better describe the heatmortality relationship, as several others have done (Vaneckova et al. 2008). Other studies have examined year-round weather and mortality data (Anderson and Bell 2009;Curriero et al. 2002;Huynen et al. 2001;McMichael et al. 2008). Factors that potentially confound the temperature-mortality relationship differ during the warm and cold seasons (Ito et al. 2007). In the winter, circulating influenza viruses tend to be the strongest correlates with daily mortality, whereas in the summer, the impact of communicable disease on mortality is negligible (CDC 2008;Mounts et al. 2000). Thus, an analysis limited to the warm season only may reduce the impact of seasonal mortality trends and increase the potential for model misspecification to conflate seasonal effects with short-term temperature effects.
The relationship of hot weather and population exposure to heat stress is complex, and other weather factors, such as precipitation and wind speed, likely play a role along with temperature and humidity. Individuals within a large and diverse city may experience different exposures to heat stress depending on their access to air conditioning (Bouchama et al. 2007), the characteristics of their dwelling (Semenza et al. 1996), and local microclimate effects (Harlan et al. 2006). Because of this complexity, a composite weather metric has the conceptual appeal of using all relevant weather data to assess the public health risks of hot weather. Our results showed that the fit of a model based on only maximum heat index was modestly improved by adding other weather metrics such as minimum temperature, barometric pressure, wind speed, and precipitation or the composite metric of SSC. Some of the added weather metrics are highly correlated with maximum heat index, however, and it is unclear that the resulting modest improvement in model fit justifies the use of a more complex model for risk assessment or guiding the issuance of heat advisories. For instance, the use of barometric pressure (either explicitly or through its association with SSC) in a weather-mortality model introduces an effect that may be mediated through mechanisms other than heat stress. For example, barometric pressure variation has been associated with risk of pulmonary embolism (Meral et al. 2005), abdominal aortic aneurysm (Harkin et al. 2005), and acute myocardial infarction (Houck et al. 2005). The increase in same-day mortality risk associated with lower barometric pressure (Schwartz 2000), therefore, may not be mitigated by measures to reduce heat exposure, such as air conditioning. In addition, using more weather variables in estimating the public health risk of hot weather adds more complexity to the task of translating weather forecasts into heat advisories that the public and public officials can understand and respond to.
To reduce the burden on the forecaster, heat health warning systems used in some cities apply a model that uses historic data, SSC, maximum temperature, season, and day in sequence of a heat wave to predict excess mortality (Kalkstein et al. 1996). As noted above, if these predictions are not based on recent data, excess mortality estimates may not be accurate. The substantial unexplained components of variation in daily mortality and sampling error in estimated risk functions add further error to estimates of daily mortality. In addition, local meteorologists still exercise discretion regarding the issuance of heat warnings, apparently considering other criteria, including the forecast heat index (Ebi et al. 2004). Our results show that the forecast maximum heat index and forecast duration of a heat wave are useful measures of the public health risk posed by heat waves in New York City. Heat index may also be simpler to use and apply consistently in issuing heat warnings compared with alternative metrics.

Conclusions
We found that, among various weather indices, maximum heat index is a useful metric for assessing public health risk due to hot weather in New York City. Although there is no clear threshold for the temperaturemortality relationship during the warm season, the nonlinear relationship found in this analysis supports issuing heat alerts when the maximum heat index is forecast to exceed approximately 95°-100°F, with augmented warnings and response as the forecast heat index and duration of a heat wave increases.
Other jurisdictions should consider conducting similar analyses with recent local data before adopting new heat-alert criteria. Whatever models or metrics are used to estimate public health risk from summertime heat, the benefits of advisories and heat emergency response will, of course, depend on the ability of the public, government agencies, and nongovernmental organizations to take actions that actually reduce heat exposure and health risks.

Research
Many studies have linked both short-term and long-term ambient particulate matter (PM) air pollution to increased morbidity and mortality of cardiovascular diseases in the general population (Clancy et al. 2002;Dockery 2001;Goldberg et al. 2006;Larrieu et al. 2007), but the biologic mechanisms of these associations remain unclear. Altered cardiac autonomic function, as reflected by alterations of heart rate variability (HRV), is considered one of the pathophysiologic pathways through which PM air pollution influences the cardiovascular system (Pope et al. 2004a;Rhoden et al. 2005). Increased PM air pollution has been associated with declines in HRV in older people and in patients with current or underlying cardiovascular diseases (Gold et al. 2000;Liao et al. 1999;Park et al. 2005;Pope et al. 2004b;Riojas-Rodriguez et al. 2006;Schwartz et al. 2005), suggesting that they are susceptible to the ambient PM air pollution. On the other hand, results from previous studies on the relationship between PM exposure and HRV in younger individuals have not been consistent. Some studies have observed negative exposure-response associations between PM exposure and HRV (Cavallari et al. 2007;Magari et al. 2001Magari et al. , 2002aVallejo et al. 2006), whereas others have report positive associations (Magari et al. 2002b;Peretz et al. 2008;Riediker et al. 2004). These different findings suggest a complicated mechanism through which the cardiac autonomic system responds to PM exposure in younger individuals.
The ambient PM air pollution in Beijing during the 2008 Olympic Games (from 8 August to 20 September) decreased markedly after a series of air quality control meas ures were implemented by the Beijing municipal government, including banning the use of more than half of the motor vehicles in Beijing every day. The marked changes of PM air pollution before, during, and after the Beijing 2008 Olympic Games provided a unique opportunity to investigate the effect of a wide range of PM exposure on cardiac autonomic function. We conducted an occupational panel study to evaluate the relationship between the PM 2.5 exposure and changes in HRV in a group of young, healthy taxi drivers who were exposed to relatively high traffic-related PM air pollution before, during, and after the Beijing 2008 Olympic Games. The objectives of this study were a) to determine whether reduced PM air pollution in Beijing during the Olympic Games improved the HRV, and b) to quantitatively evaluate the percent changes of HRV corresponding to the increase of PM levels by controlling potentially confounding variables.

Materials and Methods
Study subjects and protocol. This study was designed as a panel study. Forty-four young, nonsmoking taxi drivers were recruited from several taxi companies in Beijing. A selfadministered questionnaire was used to collect personal information, including name, sex, age, years of employment as a taxi driver, smoking status, education, and history of cardiovascular diseases or other diseases. Drivers were then physically examined and tested for seated blood pressure, resting electrocardiogram (ECG), blood cholesterol, triglycerides, and high-and low-density lipoproteins. To increase the homogeneity of subjects and reduce confounding, we used the following inclusion criteria to select subjects for this analysis: nonsmoking; no history of physiciandiagnosed cardiovascular, pulmonary, neurologic, or endocrine diseases; body mass index (BMI) ≤ 30, normal blood pressure; normal resting ECG and normal blood test results; daytime taxicab driving hours; and employment as a taxi driver for at least 1 year. Using these criteria, a total of 14 taxi drivers were initially selected to participate in the study. The study was approved by the Institutional Review Board of Peking University Health Science Center, and informed consent was obtained from each subject before the study began.
Continuous measurements of personal exposure to PM with an aerodynamic diameter ≤ 2.5 µm (PM 2.5 ) and ambulatory ECG monitoring were conducted on each subject during a 12-hr work shift (0900-2100 hour) before (26 May to 19 June), during (11 August to 5 September), and after (27 October to 14 November) the Beijing 2008 Olympic Games. Measurements were Background: Heart rate variability (HRV), a marker of cardiac autonomic function, has been associated with particulate matter (PM) air pollution, especially in older patients and those with cardio vascular diseases. However, the effect of PM exposure on cardiac autonomic function in young, healthy adults has received less attention. oBjectives: We evaluated the relationship between exposure to traffic-related PM with an aerodynamic diameter ≤ 2.5 µm (PM 2.5 ) and HRV in a highly exposed panel of taxi drivers. Methods: Continuous measurements of personal exposure to PM 2.5 and ambulatory electrocardiogram monitoring were conducted on 11 young healthy taxi drivers for a 12-hr work shift during their work time (0900-2100 hr) before, during, and after the Beijing 2008 Olympic Games. Mixedeffects regression models were used to estimate associations between PM 2.5 exposure and percent changes in 5-min HRV indices after combining data from the three time periods and controlling for potentially confounding variables. results: Personal exposures of taxi drivers to PM 2.5 changed markedly across the three time periods. The standard deviation of normal-to-normal (SDNN) intervals decreased by 2.2% [95% confidence interval (CI), -3.8% to -0.6%] with an interquartile range (IQR; 69.5 µg/m 3 ) increase in the 30-min PM 2.5 moving average, whereas the low-frequency and high-frequency powers decreased by 4.2% (95% CI, -9.0% to 0.8%) and 6.2% (95% CI, -10.7% to -1.5%), respectively, in association with an IQR increase in the 2-hr PM 2.5 moving average. conclusions: Marked changes in traffic-related PM 2.5 exposure were associated with altered cardiac autonomic function in young healthy adults. key words: air pollution, cardiac autonomic function, epidemiology, heart rate variability, particulate matter.  only on weekdays (from Monday to Friday) to reduce variation due to changes in traffic flow across the week (weekend traffic was more irregular).
HRV measurements. A three-channel Holter recorder (model MGY-H7; DM Software Inc., Stateline, NV, USA) was used for ambulatory ECG monitoring. Each subject's skin was carefully cleaned with alcohol pads by a trained technician before the electrodes were connected. Modified V 1 and V 5 leads were used to obtain the cardiac electro physiologic signals. Each subject wore the Holter recorder for 12 hr from the beginning (0900 hour) to the end (2100 hour) of the work shift, and each subject was also given a daily work diary to record personal activities such as eating or resting. Each 12-hr ambulatory ECG tape was processed using the software for the Holter recorder (Holter System, version 12.net; DM Software Inc.). Only normal-to-normal (NN) beat intervals from 600 msec to 1,500 msec with NN ratios between 0.8 and 1.2 were included in the data analysis. All abnormities and noises were excluded based on standard criteria (Task Force of the European Society of Cardiology and the North American Society of Pacing and Electrophysiology 1996) to ensure the quality of data. Three indices of HRV were calculated and used in the analysis: a) standard deviation of NN intervals (SDNN), b) low-frequency (LF) power (0.04-0.15 Hz), and c) high-frequency (HF) power (0.15-0.40 Hz). These HRV indices were all calculated in standard 5-min segments throughout the entire recording, but data collected when drivers were not in the taxicab were deleted to reduce confounding by nondriving-related activities. Because the subjects spent most of their work time inside the taxicab and always kept a stationary sitting posture, the monitoring environment was considered relatively stable. Consequently, 5-min HRV indices could be used to analyze the association between personal PM exposure and HRV with little interference due to changes in activity level or position.
Exposure to air pollutants and micro climate variables. Air pollutant concentrations were measured continuously and simultaneously using monitoring equipment installed inside each taxicab. Real-time concentrations of PM 2.5 were measured by a portable aerosol spectrometer (model 1.109; Grimm Technologies Inc., Douglasville, GA, USA), logged in 1-min intervals, and aggregated as 5-min averages. The daily mean concentration of PM 2.5 for each driver was obtained directly from the equipment by averaging all data points obtained during the day. PM 2.5 mass was also collected gravimetrically on polytetrafluoroethylene filters (Whatman Inc., Florham Park, NJ, USA) using an SKC sampling system (2.5-µm personal environmental monitor and AirChek XR5000 sampling pump; SKC Inc., Eighty Four, PA, USA). Real-time measurements of PM 2.5 were influenced by seasonal changes in weather. Therefore, real-time concentrations were calibrated based on the average PM 2.5 mass concentration obtained from the gravimetric analysis each day to ensure a more reliable measure of individual exposure [see Supplemental Material available online (doi:10.1289/ehp.0900818.S1 via http://dx.doi.org/)].
Concentrations of carbon monoxide (CO), nitrogen dioxide (NO 2 ), and nitric oxide (NO) were also measured in order to evaluate confounding effects of these gaseous pollutants on the relationship between PM 2.5 and HRV. CO was measured with a model T15n enhanced CO measurer (Langan Products Inc., San Francisco, CA, USA), and NO 2 and NO were measured using a passive sampler (Ogawa Air Inc., Osaka, Japan) and analyzed according to the manufacturer's specifications. Real-time temperature and relative humidity (RH) data were measured using a HOBO Pro V2 temperature/RH logger (Onset Corp., Pocasset, MA, USA), logged at 1-min intervals and aggregated as 5-min averages.
Traffic routes were recorded using a model M-241 wireless global positioning system (GPS) logger (HOLUX Technology Inc., Hsinchu, Taiwan). Drivers used a work diary to record their activities, car window status (open vs. closed), and air conditioner use (on vs. off) throughout the monitoring periods.
Statistical analysis. Field data were checked for completeness on site, and efforts were made to ascertain missing data as needed. HRV indices were log 10 -transformed to improve normality and stabilize the variance. Data collected when drivers were not in their taxicab were identified based on work diaries and GPS records and were excluded before analysis. All statistical analyses were performed using SAS software for Windows (version 9.1; SAS Institute Inc., Cary, NC, USA).
We conducted measurements on each subject during each time period in order to compare associations between different PM 2.5 exposure levels and HRV. To facilitate logical and accurate comparisons and interpretation, we tested the data structure for each subject before analysis. The number of 5-min HRV indices for each subject in each of the three time periods should have been approximately equal; otherwise, the distribution of raw 5-min HRV indices might have been biased. We excluded 3 subjects from the final analysis who had substantially fewer observations in one of the three time periods relative to the other two because of less time (< 5 hr) spent inside the taxicab during that shift. Therefore, data from the remaining 11 subjects were used for the analysis on the association between PM 2.5 exposure and HRV.
We first compared raw 5-min HRV indices during the three time periods using oneway analysis of variance (ANOVA) and then used mixed-effects regression models to analyze combined data from all three time periods and estimate associations of varying exposure levels of PM 2.5 with HRV. Subject and day of the year were treated as random effects, whereas age, time of day, log 10 -transformed heart rate, temperature, and RH were included as fixed effects. The gaseous pollutants CO, NO 2 , and NO were also individually evaluated as copollutants in two-pollutant models with PM 2.5 . Both linear and quadratic terms were used to model moving averages of temperature and RH to account for nonlinear associations with HRV. Other potential confounders were also investigated: sex, BMI, years as a taxi driver, day of week, status of cab windows, and status of cab air conditioner. These variables did not show significant confounding effects in any of the regression models. Therefore, final models included random effects terms for subject and day of the year and fixed effects terms for age, time of day, log 10 -transformed heart rate, temperature, and RH only. Seven different PM 2.5 moving averages were generated for lagged intervals ranging from 5 min to 4 hr using the calibrated 5-min average real-time PM 2.5 concentrations. Corresponding log 10 -transformed 5-min HRV indices were regressed on the different PM 2.5 moving averages. The 30-min and 2-hr moving averages were chosen as the main exposure metrics because PM 2.5 effects on HRV were mainly observed for these two lagged intervals. Analyses examining the influence of outlying exposure values on regression results were also performed after excluding observations in the top 5% and bottom 5% of the PM 2.5 moving averages, respectively.
Final results are presented as the estimated percent change in each 5-min HRV index associated with an interquartile range (IQR) increase of each air pollutant as [10 (β × IQR) -1] × 100%, with 95% confidence intervals (CIs), where β and SE are the effect estimate and its standard error. Table 1 presents detailed information on the characteristics of the 11 subjects included in the final analysis. Six were females (55%), the mean age of all subjects was 35.5 years (range, 27-41 years), and the mean time of employment as a taxi driver was 6.0 years. All 11 subjects were in good health with blood pressures, resting heart rate, and blood cholesterol, triglycerides, and lipoproteins in normal ranges. None of the subjects had a history of cardiovascular diseases or smoking.

Results
In Table 2, we show the daily averages of air pollutant levels and microclimate variables inside the taxicab during each of the three time periods. According to the diaries filled in by drivers, they spent an average of 87% of their work time inside the taxicab each day; thus, the results of exposure monitoring represent the taxi drivers' typical personal exposure levels during their work time. PM 2.5 exposure levels are presented as daily averaged real-time PM 2.5 concentrations and daily averaged PM 2.5 mass concentrations obtained from gravimetric analysis. Results from the two methods are highly correlated (r = 0.973, p < 0.001); therefore, we present results for PM 2.5 mass concentrations unless otherwise indicated. Average individual air pollutant exposure levels during the Olympic Games were much lower than levels before and after the Olympic Games, especially levels of PM 2.5 , which decreased to less than half the average level recorded before the Olympic Games (45.2 µg/m 3 vs. 105.5 µg/m 3 , respectively). The study lasted for about half a year, from the end of spring to the end of fall, resulting in differences in temperature and RH across the three time periods. Table 3 presents the distribution of raw 5-min HRV indices across the three time periods. It is evident that all three indices of HRV were higher during the Olympic Games than before and after the Olympic Games. Oneway ANOVA indicated significant differences for all three 5-min HRV indices among the three time periods (p < 0.0001).
We estimated associations between the 5-min HRV indices and different PM 2.5 moving averages (from 5 min to 4 hr) after combining data for each subject across the three time periods and adjusting for potentially confounding variables. IQR increases in PM 2.5 mass concentrations were associated with declines in all three 5-min HRV indices (Figure 1). Table 4 shows adjusted associations between IQR increases in 30-min and 2-hr PM 2.5 moving averages and the 5-min HRV indices. An IQR (69.5 µg/m 3 ) increase in the 30-min PM 2.5 moving average was associated with a 2.2% decline (95% CI, -3.8% to -0.6%) in the SDNN. Associations with LF and HF powers were stronger for an IQR increase in the 2-hr PM 2.5 moving average, with decreases of 4.2% (95% CI, -9.0% to 0.8%) and 6.2% (95% CI, -10.7% to -1.5%), respectively. We also performed the same regression models after excluding observations in the top 5% and bottom 5% of the PM 2.5 moving averages, to evaluate the influence of outlying exposure values on associations. Effect estimates became more negative, especially those for LF and HF powers in association with the 30-min PM 2.5 moving average.   Excluding outlying exposure values a 5-min SDNN -2.2 (-3.8 to -0.6)** -3.5 (-5.7 to -1.2)** -1.1 (-3.3 to 1.1) -1.6 (-4.2 to 1.1) 5-min LF power -2.8 (-6.6 to 1.0) -5.8 (-10.7 to -0.5)* -4.2 (-9.0 to 0.8) -4.8 (-10.6 to 1.3) 5-min HF power -4.1 (-7.6 to -0.5)* -7.1 (-11.6 to -2.4)** -6.2 (-10.7 to -1.5)* -7.8 (-13.1 to -2.3)** To further explore potential heterogeneity from person to person, we fitted subjectspecific regression models using data from all three time periods. Figure 2 shows the distribution of the subject-specific effect estimates for an IQR increase in the 30-min PM 2.5 moving average for all three 5-min HRV indices. The relatively wide range of effect estimates suggests that individual heterogeneity was present despite the strict inclusion criteria used to select participants for the study. For example, 8 effect estimates for percent change in SDNN with an IQR increase in PM 2.5 were negative, and 3 were positive. Results for LF and HF powers were similar.
Furthermore, we fitted exposure-response curves for PM 2.5 on 5-min HRV indices based on regression models combining data across the three time periods. Figure 3 shows the estimated relationship between the 30-min PM 2.5 and 5-min SDNN using a loess smoother with a smoothing parameter of 0.6. The relationship between PM 2.5 and SDNN is not linear. Specifically, there is a positive association between PM 2.5 and percent change in SDNN when PM 2.5 concentrations are relatively low (up to 49.3 µg/m 3 ), but the association becomes negative when PM 2.5 concentrations are higher, and eventually levels off for the highest levels of PM 2.5 . Results for the LF and HF powers were similar to those for SDNN [see Supplemental Material, Figure 1 (doi:10.1289/ehp.0900818.S1)]. These findings may provide a reasonable explanation for the more negative estimated percent changes in 5-min HRV indices after excluding observations in the top 5% and bottom 5% of the PM 2.5 moving averages (Table 4).
We investigated potential confounding by the gaseous pollutants CO, NO 2 , and NO by including each of them in a two-pollutant model with PM 2.5 . In general, adjusting for CO weakened effect estimates for PM 2.5 on HRV, whereas adjusting for NO 2 made them somewhat stronger. Adjusting for the copollutants decreased the precision of the estimates, but the overall results were generally consistent with estimates that were not adjusted for the copollutants [see Supplemental Material, Table 1 (doi:10.1289/ehp.0900818.S1)]. Separate models to estimate associations of CO, NO 2 , and NO with HRV generally did not support associations between these pollutants and the outcomes [see Supplemental Material, Table 2 (doi:10.1289/ehp.0900818.S1)].

Discussion
In this study, we assessed real-time personal exposure to traffic-related PM 2.5 across three time periods associated with substantial differences in exposure levels, and evaluated the relationship between PM 2.5 and HRV as meas ured by ambulatory ECG in a selected group of young, healthy, nonsmoking taxi drivers. The levels of traffic-related PM 2.5 to which taxi drivers were exposed in our study (especially the levels before the Olympic Games) were much higher than ambient or personal PM exposure levels reported in previous studies of populations in other countries (Holguín et al. 2003;Park et al. 2005;Riediker et al. 2004;Riojas-Rodriguez et al. 2006). To the best of our knowledge, this is the first study to show effects of varying levels of traffic-related PM exposure on cardiac autonomic function in a highly exposed occupational panel.
A major strength of this study is that we evaluated the same workers during three different time periods that had markedly different PM air pollution levels, which allowed us to compare the corresponding levels of 5-min HRV indices among different time periods. A comparison of raw 5-min HRV indices indicated that the low PM 2.5 exposure period (during the Olympic Games) was associated with relatively high HRV, whereas higher PM 2.5 exposures (before and after the Olympic Games) were associated with relatively low HRV.
Decreased HRV has been regarded as a risk factor for cardiac morbidity and mortality (La Rovere et al. 2003;Lombardi et al. 2001), and previous studies have linked decreased HRV to ambient PM air pollution (Gold et al. 2000;Liao et al. 1999Liao et al. , 2004Magari et al. 2002a;Park et al. 2005;Pope et al. 2004b), occupational PM exposure (Cavallari et al. 2007;Magari et al. 2001Magari et al. , 2002a, and PM air pollution related to traffic (Schwartz Our results suggest that marked changes of PM air pollution may lead to cardiac autonomic imbalance in young healthy individuals, as indicated by declines in several 5-min HRV indices in taxi drivers. Several previous studies also have associated PM exposure with HRV in younger individuals (Cavallari et al. 2007;Magari et al. 2001Magari et al. , 2002aMagari et al. , 2002bPeretz et al. 2008;Riediker et al. 2004;Vallejo et al. 2006), but results reported in these studies were not consistent with each other. A previous study of nine young, healthy, nonsmoking, male patrol troopers with a mean age of 27.3 years used short-term HRV (10 min) to evaluate the effect of traffic-related PM exposure on cardiac autonomic function (Riediker et al. 2004). The PM 2.5 exposure levels in patrol cars were monitored for 4 consecutive days on each trooper during their work shift (from 1500 hr to midnight shift), and several periods of shortterm HRV were assessed. Results showed that each 10-µg/m 3 increase of in-vehicle PM 2.5 (average of 24 µg/m 3 ) was associated with a significant 11.7% increase in SDNN and a significant 14.8% increase in HF power the next morning. These findings suggest that that young, healthy individuals may respond to PM 2.5 exposure differently than do older people and patients who have reduced cardiovascular dynamism and altered vagal (parasympathetic) function. On the other hand, some studies have reported declines of HRV in younger individuals who were exposed to PM occupationally (Cavallari et al. 2007;Magari et al. 2001) or environmentally (Magari et al. 2002a;Vallejo et al. 2006) at higher levels. Magari et al. (2001) observed a significant 2.66% (95% CI, -3.75% to -1.58%) decline in 5-min SDNN for every 1-mg/m 3 increase in the 4-hr PM 2.5 moving average due to occupational, personal PM 2.5 exposure (arithmetic mean ± SD, 0.697 ± 0.665 mg/m 3 ) in 33 boilermakers with a mean age of 38.1 years. Cavallari et al. (2007) found similar results in a group of male boilermaker welders who were exposed to an average PM 2.5 level of 0.73 ± 0.50 mg/m 3 , using several indices of nighttime HRV (0000-0700 hr) as indicators; a 1-mg/m 3 increase of 8-hr time-weighted average concentration of workday PM 2.5 exposure was associated with a significant decline of 8.32 msec (95% CI, -16.29 to -0.35 msec) in nighttime sqare root of the mean squared differences of successive intervals (rMSSD), after adjusting for nonworking nighttime HRV, age, and smoking.
The results of our study are generally consistent with those that have reported declines of HRV with PM exposure in younger individuals. However, results from regression models for each subject showed heterogeneity among subjects, and we found that several subjects in our study had positive associations for the three 5-min HRV indices with traffic-related PM exposure. These subjects might have responded to PM in a similar way as reported by Riediker et al. (2004). On the other hand, smoothed curves showing associations between PM exposure and 5-min HRV indices also indicated that lower PM exposures were associated with increases in HRV whereas higher PM exposures were associated with decreases in HRV. Thus, differences between studies may also be related to differences in exposure levels. Overall, factors affecting heterogeneity of responses to PM exposure remain unclear, and further study is needed to explore the mechanisms behind the different responses to PM exposure in younger individuals.
Traffic density was different across the three time periods because of the Olympic Games (half of the motor vehicles in Beijing were banned from the streets each day during the Olympic Games), and this change in traffic density might have introduced a confounding effect by influencing the stress level of subjects, because previous studies have found associations between HRV and work-related stress (Kang et al. 2004;Saito et al. 2008;Vrijkotte et al. 2000). To partly control variation in stress levels across different time periods, we used a random effect for day of the year in regression models.
In summary, our study revealed an association between PM exposure and HRV in a group of young, healthy taxi drivers who were exposed to much higher traffic-related PM air pollution than the general population. This finding provides further evidence that PM air pollution affects cardiac autonomic function in young, healthy individuals.