Data sources. We obtained respiratory hospitalization data from 1991–2004 among NYS residents from the NYS Department of Health Statewide Planning and Research Cooperative System (SPARCS), a legislatively mandated database of hospital discharge data for approximately 95% of all NYS acute care admissions, excluding admissions to psychiatric and federal hospitals (NYS Department of Health 2002). Data included principal diagnosis, admission and discharge dates, sources of payment, total charges, date of birth, sex, race, ethnicity, and residential street address. About 94% of addresses were geocoded to street level and 5% to ZIP code level. Less than 1% of addresses could not be geocoded and were excluded from the analysis. Family income information was obtained from 1990 and 2000 U.S. Census data at the census block level (U.S. Census Bureau 1992, 2002).
Meteorological data for NYS during 1991–2004, including hourly observations for temperature, barometric pressure, and dew point, were provided by the Data Support Section of the Computational and Information Systems Laboratory at the National Center for Atmospheric Research (NCAR; Boulder, CO) from airport weather stations. Hourly ambient ozone data were obtained from the NYS Department of Environmental Conservation (Albany, NY) ambient air monitors. We used an 8-hr maximum average ozone concentration limited to the 1000–1800 hours, which represents the most likely time for outdoor exposure. Both meteorological and ambient air monitoring locations were geocoded to street level. Projections of the extent and geographical climate distribution were obtained from the IPCC (NCAR 2007a, 2007b, 2009), including projected temperature, barometric pressure, and specific humidity data for 2046–2065 and 2080–2099 (i.e., about 50 and 100 years from the baseline period).
Study population and health outcomes.
The study population included all NYS residents. We assessed respiratory admissions, related hospitalization costs, days hospitalized, and lost productivity from days hospitalized. Specifically, we counted all hospital admissions in the summer (June–August) with a principal diagnosis of respiratory disease among NYS residents during 1991–2004. Based on the International Classification of Disease, 9th Revision, Clinical Modification
(ICD-9; Department of Health and Human Services 1997
), respiratory diseases included chronic bronchitis (ICD-9 code 491), emphysema (ICD-9 code 492), asthma (ICD-9 code 493), and chronic obstructive pulmonary disease (COPD; ICD-9 code 496). For children ≤ 4 years of age, we included acute bronchitis and bronchiolitis (ICD-9 code 466) and bronchitis, not specified as acute or chronic (ICD-9 code 490) because these respiratory illnesses are difficult to distinguish from asthma among very young children. The hospitalization charges listed in SPARCS do not reflect actual inpatient hospitalization costs; therefore, we multiplied the hospitalization charges indicated in SPARCS by the average cost-to-charge ratio (0.54) for NYS obtained from the Healthcare Cost and Utilization Project (2011)
. We used the length of stay for each patient to estimate the economic cost of lost productivity from days hospitalized according to market and household productivity estimates for U.S. adults by age and sex (Grosse et al. 2009
Meteorological and exposure indicators
. We identified 14 weather regions with relatively homogeneous weather and ozone exposures by overlaying and merging the 10 NYS climate divisions (National Climate Data Center, Asheville, NC) with the 11 ozone regions developed for NYS by Chinery and Walker (2009)
. Each hospitalization was assigned to a weather region on the basis of geocoded residential address.
Daily mean apparent temperature (AT; an index of human discomfort resulting from the combined effects of heat and humidity) was calculated in degrees Celsius as AT = –2.653 + 0.994T + [0.0153 × (dewpoint)2
] as described previously (Kalkstein and Valimont 1986
, Steadman 1979
), where T represents daily mean temperature. Although the relationship between temperature or AT and respiratory disease is usually U- or V-shaped (Linares and Diaz 2008
; McMichael et al. 2008
), we used a linear-threshold model to quantify the effect of high temperature. The threshold (T0
) was selected for each region after modeling all possible values (70–105°F) and selecting the one with the lowest model deviance for each region (Armstrong 2006
; McMichael et al. 2008
). We also identified two alternate extreme heat indicators: a
) the 90th percentile of AT based on the summer AT distribution from 1991–2004, and b
) daily AT > 90°F.
Climate scenarios. The IPCC has defined a range of possible future trends in greenhouse gas emissions (IPCC 2007). The scenarios presented in the IPCC Special Report on Emissions Scenarios (SRES) (IPCC 2000) are plausible indications of what the future could be like over decades or centuries (IPCC-Task Group on Data and Scenario Support for Impacts and Climate Analysis A 2007). To represent a wide range of possible future climates, we selected three of the SRES scenarios—high (A2), mid (A1B), and low (B1) emissions—based on alternative assumptions about changes in economy, technology, demographics, and energy demand (IPCC 2000, 2007). A2 assumes a very heterogeneous world with continuously increasing population growth, slow and regionally oriented economic development, and slow technological change. A1B assumes a world of very rapid economic growth, a global population that peaks in midcentury and then gradually declines, and rapid introduction of new and more efficient technologies with a balance across all energy sources. B1 assumes a convergent world, with the same population growth as in scenario A1B, but with more rapid changes toward a service and information economy and with a reduction in material intensity and clean and resource-efficient technologies (IPCC 2000, 2007).
Projection of future summer AT distributions
. We estimated future AT by using temperature, barometric pressure, and specific humidity obtained from the IPCC (NCAR 2007a, 2007b, 2009), which applied the Community Climate System model Version 3 (CCSM3), based on the three climate scenarios described above and constructed according to longitude, latitude, and time with grid cells of 155 km × 155 km. We assumed that regional variation in climate across the 14 weather regions at baseline would remain unchanged. We used the change in spatially averaged mean summer daily AT for each region from baseline to midcentury (2046–2065) and the end of the century (2080–2099) under each climate scenario (Bambrick et al. 2008
; McMichael et al. 2004
. To assess public health impacts, we estimated the relationship between daily temperature variation and respiratory admissions using a two-stage Bayesian model that included a regional analysis and a statewide estimate adjusted for regional confounders. In stage 1, we estimated the association between extreme heat and respiratory disease hospitalization for each of the 14 NYS regions using generalized additive models (GAM) (Hastie and Tibshirani 1989
) with Poisson distributional errors and a log link function in SAS (version 9.2; SAS Institute Inc., Cary, NC). We assumed a log-linear increase in health risk above a temperature threshold (T0
), which was determined by comparing the maximum likelihood estimates over all possible threshold values in the range of data and using the value with the lowest deviance. We used a linear association between hospitalization and each 1°F increase in AT > T0
to estimate the extreme heat effect for each region:
log(count) = α0 + s(AT < T0, df) + β0(AT > T0) + s(date, df) + s(O3, df) + β1 +…+ β8 + ε, 
is the threshold value of AT, β0
is the slope parameter for AT > T0
(representing the risk of hospitalization with each 1°F increase in AT > T0
), and α0
and ε are the intercept and error terms, respectively. Spline curves, indicated by s(AT < T0
, df), s(date, df), and s(O3
, df), were used to model the effects of AT < T0
, long-term trends and seasonal variation (date), and ozone (O3
). Degrees of freedom (df) were determined using an automatic procedure based on minimizing the sum of absolute values for the first 30 items of the partial autocorrelation function of the model residuals (Armstrong 2006
). We also controlled the effect of day of the week (with Monday–Saturday represented by β1
), and the Northeast blackout events that occurred on 14 and 15 August 2003 (β7
). Model fit was assessed by Bartlett’s Kolmogorov–Smirnov statistic (Bartlett 1978
). We also checked the model residuals for autocorrelation and partial autocorrelation functions to rule out seasonality or other patterns (Armstrong 2006
; Lin et al. 2009
In stage 2, we pooled region-specific estimates to generate a statewide estimate using a Bayesian hierarchical model (Dominici et al. 2002
). We controlled for region-level covariates using yearly data estimated from 1990 and 2000 U.S. Census data, including population density, health-care access (minimum distance to clinics), race and ethnicity (percent of black and Hispanic residents), percent of residents with ≤ high school education, mean apparent temperature during June–August, percent living below the poverty level, and percent of the regional population that were elderly (age ≥ 75 years) and living alone. Pooling of information across regions can potentially improve statistical power and the generalizability of the results, as well as account for geographic heterogeneity of effects (Dominici 2002
). All second stage analyses were conducted using tlnise in R, version 0.2-7 (Everson and Morris 2000
Daily excess admissions attributable to extreme heat were computed at baseline (1991–2004) and projected for 2046–2065 and 2080–2009 as
A = R × ΔT × M, 
where A represents the estimated number of daily admissions attributable to AT > T0 for each time period; R is the estimated percentage increase in admissions per 1°F increase in AT > T0 at baseline based on Equation 1 [i.e., R = exp(β0) – 1]; M is the mean number of daily respiratory admissions during June–August at baseline; and ΔT is the difference between the mean daily AT for each time period and the baseline T0 on days when AT > T0. We used the thresholds identified from the baseline data (1991–2004) for each region to project future impacts.
Daily excess hospitalization costs and days of hospitalization attributable to extreme heat were computed as
C = A × D, 
represents either daily temperature-attributable hospitalization costs or days hospitalized; A
is as defined in Equation 2; and D
is the average number of days hospitalized per hospitalization at baseline (Campbell-Lendrum and Woodruff 2006
; McMichael et al. 2004
). The average cost per day of hospitalization was adjusted for inflation (by month and year) for each time period and standardized to 2004 U.S. dollars to ensure that costs were comparable across the different years in our study (InflationData.com 2011
). Because a dollar in the future is considered to be of less value than a current dollar, it is common practice in health cost estimation to discount for future costs incurred across different time periods to derive the current worth of all future amounts. We used an annual discount rate of 3% as recommended by the U.S. Panel on Cost-Effectiveness in Health and Medicine to estimate the 2004 dollar value of the future stream of costs (Goodman 2004
; Phillips and Chen 2002
). Excess lost productivity from days hospitalized was computed by multiplying estimated excess days of hospitalization by age-specific daily production values presented by Grosse et al. (2009)
Stratified analyses were conducted based on individual level data such as sex, age, specific disease group, insurance type, and census-block group–level family income. We also conducted a sensitivity analysis using another global climate model, the Centre National de Recherches Meteorologiques Coupled Global Climate Model Version 3 (CM3) (Salas 2005a
), to compare and validate our results. Moreover, we performed sensitivity analyses to examine whether the estimate of excess heat-related health burden was due to population composition changes in age distribution or race/ethnicity. We used population growth percentages on these specific subgroups to roughly estimate excess admissions from the subgroups without population growth and those with population growth according to U.S. population projections from 1998–2000, 2050–2070, and 2075–2100 (U.S. Census Bureau 2000) and the U.S. Hispanic population from 2000 to 2050 (U.S. Census Bureau 2011a).