Variation in Estimated Ozone-Related Health Impacts of Climate Change due to Modeling Choices and Assumptions

Background: Future climate change may cause air quality degradation via climate-induced changes in meteorology, atmospheric chemistry, and emissions into the air. Few studies have explicitly modeled the potential relationships between climate change, air quality, and human health, and fewer still have investigated the sensitivity of estimates to the underlying modeling choices. Objectives: Our goal was to assess the sensitivity of estimated ozone-related human health impacts of climate change to key modeling choices. Methods: Our analysis included seven modeling systems in which a climate change model is linked to an air quality model, five population projections, and multiple concentration–response functions. Using the U.S. Environmental Protection Agency’s (EPA’s) Environmental Benefits Mapping and Analysis Program (BenMAP), we estimated future ozone (O3)-related health effects in the United States attributable to simulated climate change between the years 2000 and approximately 2050, given each combination of modeling choices. Health effects and concentration–response functions were chosen to match those used in the U.S. EPA’s 2008 Regulatory Impact Analysis of the National Ambient Air Quality Standards for O3. Results: Different combinations of methodological choices produced a range of estimates of national O3-related mortality from roughly 600 deaths avoided as a result of climate change to 2,500 deaths attributable to climate change (although the large majority produced increases in mortality). The choice of the climate change and the air quality model reflected the greatest source of uncertainty, with the other modeling choices having lesser but still substantial effects. Conclusions: Our results highlight the need to use an ensemble approach, instead of relying on any one set of modeling choices, to assess the potential risks associated with O3-related human health effects resulting from climate change.


Adjustment of Air Quality Output from Modeling Systems
Benefits analysts who deal with air pollution generally have more confidence in monitored air pollutant concentrations than modeled concentrations, since monitor values are actual measurements. However, unlike modeled values, monitors do not exist in all grid cells of an air quality model grid. Therefore, following EPA's typical procedures for a future-year analysis, we applied a Voronoi Neighbor Averaging (VNA) spatial adjustment to the without-climate-change O 3 metrics, and VNA spatial and temporal adjustments to the with-climate-change O 3 metrics, using both monitor and modeled values in BenMAP. These spatial and temporal adjustment procedures are described in detail in Sections C.3.2 and C.3.3 in Appendix C of Abt Associates Inc. (2010). The change in O 3 due to climate change c. 2050 (Δx in equation (2) in the paper) was then calculated in each of the cells in the 30 km x 30 km grid used for the analysis. Woods & Poole population projections are available only through 2030, however, whereas our analysis is c. 2050. Therefore, it was necessary to extrapolate Woods & Poole population projections to 2050. Given the large number of projected population series, we used automatic forecasting algorithms that have been implemented in the forecast package for R (Hyndman 2009;R Development Core Team 2009).
In order to generate our forecasts, we used a set of models that belong to the class of exponential smoothing (ES) forecasting methods (see Gardner 2006 andHyndman 2009 for the theoretical background of exponential smoothing models.) We evaluated the following three ES models: simple exponential smoothing, linear exponential smoothing, and damped-trend exponential smoothing. These models are categorized by their trend component: none, additive, and damped, respectively. We estimated all three models for each projected population series and then chose the best-fitting model based on the Bayesian Information Criterion (BIC), a standard measure of goodness of fit of a model to the underlying data. The best model was used to forecast each series out to 2050.
These ES forecasting methods try to extrapolate trends seen in a given set of years beyond the final year of the dataset. Thus the set of years on which the extrapolation is based could affect the resulting extrapolation. We applied the method described above to each of the following three series of years: 2000 -2030; 2010 -2030; and 2020 -2030. We then averaged the results. This gives somewhat more weight to the latter years, which is appropriate, since time trends may change over the longer course of years beginning in 2000 or 2010.
The resulting 2050 population forecast was adjusted to match the Census national population projection for 2050 (U.S. Census Bureau 2010a). For each of the 304 population sub-groups we calculated the 2050 national total, as implied by the extrapolated Woods & Poole population projections. We then calculated percent differences between these population totals and the population totals projected by the Census Bureau. Finally, we adjusted each county-and population subgroup-specific extrapolated Woods & Poole projection using corresponding percent differences. This method allowed us to match the Census Bureau national population projection as well as preserve some of the county-specific demographic patterns and trends.

Descriptions of Selected ICLUS population projections
The base case population projection uses the standard Census projection method (U.S. Census Bureau 2010b). A1 represents a world of fast economic development, low population growth, and high global integration. Fertility is assumed to decline and remain low similar to recent and current experience in many European countries (Sardon 2004). The A2 storyline represents a world of continued economic development, but with a more regional focus and slower economic convergence between regions. Fertility is assumed to be higher than in A1. International migration is assumed to be low because a regionallyoriented world would result in more restricted movements across borders. Domestic migration is high because, like in A1, the continued focus on economic development is likely to encourage movement within the United States.

Pooling of Concentration-Response Functions
For several health endpoints, two or more C-R functions were pooled. In particular, for respiratory hospital admissions we undertook the following pooling procedure: 1. Moolgavkar et al. (1997) estimated C-R functions in Minneapolis for hospital admissions (HA), pneumonia (ICD-9 codes 480-487) and HA, COPD (ICD 490-496). We summed the results from these two non-overlapping subcategories. 2. Schwartz (1994b) also estimated C-R functions in Minneapolis for the same two subcategories. However, this study found a significant effect only for HA, pneumonia. So the estimate of "PM-related HA for respiratory illness" in Minneapolis based on Schwartz (1994b) was taken to be just PM-related HA, pneumonia.
3. The estimates of "PM-related HA for respiratory illness" in Minneapolis from (1) and (2) above were pooled using a fixed effects pooling method. (When choosing fixed effects as the pooling method, pooling weights are generated automatically based on the inverse variance of each input result, with the weights normalized to sum to one. Results with a larger absolute variance get smaller weights. For more details, see Section L.2.1.3 in Abt Associates Inc. 2010). 4. Schwartz (1994a) estimated C-R functions for the same two non-overlapping subcategories in Detroit. We similarly summed these results. 5. Finally, Schwartz (1995) estimated C-R functions for "HA, all respiratory" in New Haven, CT and Tacoma, WA. We pooled the HA, All respiratory results from these C-R functions with the results from steps (3) and (4). (For more details, see Section L.2.1.4 in Abt Associates Inc. 2010).
To obtain the asthma ER visits results, we pooled Peel et al. (2005) and Wilson et al. (2005) using the random/fixed effects method (for more detail see Section L.2.1.4 in Abt Associates Inc. 2010). To obtain the results for school absence days, we pooled Gilliland et al. (2001) and Chen et al. (2000) also using the random/fixed effects method.

Calculation of Baseline Incidence Rates
We obtained individual-level mortality data, including residence county FIPS codes, 1 age at death, month of death, and underlying causes (ICD-10 codes), for years 2004-2006 for the entire United States from the Centers for Disease Control (CDC), National Center for Health Statistics (NCHS). The detailed mortality data allowed us to generate causespecific death counts at the county level for selected age groups. The county-level death counts were then divided by the corresponding county-level population to obtain the mortality rates. To provide more stable estimates, we used three years ( where i represents the specific cause of mortality (e.g., non-accidental mortality), j represents a specific county, and k represents a specific age group.
Mortality rates based on 20 or fewer deaths were considered unreliable (see NYSDOH 1999 for an explanation). If the rate for a given cause of death was unreliable in certain counties in a state, we summed up the deaths attributed to that cause in those counties, as well as the populations in those counties and created an aggregate rate for that cause of death in those counties. If that aggregate "state-level" rate was unreliable, we aggregated to the region level (using the four regions defined by the U.S. Bureau of the Census), and if the region-level rate was still unreliable, we aggregated to the national level. At each level of aggregation, only those counties with unreliable rates for the specified cause of death were included. So, for example, if 5 counties in a given state had unreliable rates for a specific cause of death, a "state-level" rate was created by summing the deaths from that cause across those counties and dividing by the sum of the populations in those counties. If this "state-level" rate was still unreliable, we repeated the process at the region level. The aggregate rate estimates were applied only to counties that had "unreliable" data and that estimates for all other counties were based on county-specific estimates.
To project age-and county-specific mortality rates developed using 2004-2006 data to the year 2050, we calculated growth ratios using a series of Census Bureau projected national mortality rates (U.S. Census Bureau 2010a). The procedure we used was as follows: • For each age group, we calculated the ratio of the Census Bureau national mortality rate projection in year 2050 to the national mortality rate in 2005. Note that the Census Bureau projected mortality rates were derived from crude death rates. The following formula, given by Chiang (1967)  • To estimate mortality rates in 2050 that are both age-group-specific and countyspecific, we multiplied age-group-specific mortality rates for [2004][2005][2006]  Note that future mortality rates are projected to decrease over time.

Hospitalizations
Regional hospitalization counts were obtained from the National Center for Health Statistics' (NCHS) National Hospital Discharge Survey (NHDS) (CDC 2008). NHDS is a sample-based survey of non-Federal, short-stay hospitals (<30 days), and is the principal source of nationwide hospitalization data. Note that the following hospital types are excluded from the survey: hospitals with an average patient length of stay of greater than 30 days, federal, military, Department of Veterans Affairs hospitals, institutional hospitals (e.g. prisons), and hospitals with fewer than six beds. The survey collects data on patient characteristics, diagnoses, and medical procedures. Public use data files for the year 1999 survey were downloaded and processed to estimate hospitalization counts by region (CDC 2010a). NCHS groups states into four regions using the following groupings defined by the U.S. Census Bureau (2001) We used the 2000 Census of Population and Housing to obtain more age specificity, and then corrected the 2000 Census figures so that the total population equaled the total for 1999 forecasted by NHDS. In particular, for each type of hospital admission (ICD code or codes) we: (1) calculated the count of hospital admissions by region in 1999 for the age groups of interest, (2) calculated the 2000 regional populations corresponding to these age groups, (3) calculated regional correction factors that equal the regional total population in 1999 divided by the regional total population in 2000, (4) multiplied the 2000 population estimates by these correction factors, (5) divided the 1999 regional count of hospital admissions by the estimated 1999 population, and (6) applied the regional rates to every county in that region.
Like mortality rates, the hospitalization rates are cause-specific and the hospital admissions endpoints are defined by different combinations of ICD codes that are used in the selected epidemiological studies.

Emergency Room Visits for Asthma
Regional counts of asthma-related emergency room visit counts were obtained from the National Hospital Ambulatory Medical Care Survey (NHAMCS) (CDC 2010b). NHAMCS is a sample-based survey, conducted by NCHS. The target universe of the NHAMCS is in-person visits made in the United States to emergency and outpatient departments of non-Federal, short-stay hospitals (hospitals with an average stay of less than 30 days) or those whose specialty is general (medical or surgical) or children's general. Public use data files for the year 2000 survey were downloaded and processed to estimate hospitalization counts by region (CDC 2010c). We obtained population estimates from the 2000 Census of Population and Housing. The NCHS regional groupings described above were used to estimate regional emergency room visit rates.  Illinois-1  633000  925000  743000  880000  659000  Illinois-2  638000  937000  755000  893000  650000  CMU  522000  745000  599000  679000  545000  Harvard  299000  445000  362000  422000  347000  GNM  50000  67000  44000  35000  -29000  NERL -25000 -50000 -50000 -67000 -84000 WSU -134000 -212000 -153000 -197000 -27000 a Rounded to the nearest 1000. b Based on current rather than projected future baseline incidence rates. Supplemental Material, Figure S1. Age distributions of ICLUS_A1 and ICLUS_A2 population projections to c. 2050