Projections of temperature-related excess mortality under climate change scenarios

Summary Background Climate change can directly affect human health by varying exposure to non-optimal outdoor temperature. However, evidence on this direct impact at a global scale is limited, mainly due to issues in modelling and projecting complex and highly heterogeneous epidemiological relationships across different populations and climates. Methods We collected observed daily time series of mean temperature and mortality counts for all causes or non-external causes only, in periods ranging from Jan 1, 1984, to Dec 31, 2015, from various locations across the globe through the Multi-Country Multi-City Collaborative Research Network. We estimated temperature–mortality relationships through a two-stage time series design. We generated current and future daily mean temperature series under four scenarios of climate change, determined by varying trajectories of greenhouse gas emissions, using five general circulation models. We projected excess mortality for cold and heat and their net change in 1990–2099 under each scenario of climate change, assuming no adaptation or population changes. Findings Our dataset comprised 451 locations in 23 countries across nine regions of the world, including 85 879 895 deaths. Results indicate, on average, a net increase in temperature-related excess mortality under high-emission scenarios, although with important geographical differences. In temperate areas such as northern Europe, east Asia, and Australia, the less intense warming and large decrease in cold-related excess would induce a null or marginally negative net effect, with the net change in 2090–99 compared with 2010–19 ranging from −1·2% (empirical 95% CI −3·6 to 1·4) in Australia to −0·1% (−2·1 to 1·6) in east Asia under the highest emission scenario, although the decreasing trends would reverse during the course of the century. Conversely, warmer regions, such as the central and southern parts of America or Europe, and especially southeast Asia, would experience a sharp surge in heat-related impacts and extremely large net increases, with the net change at the end of the century ranging from 3·0% (−3·0 to 9·3) in Central America to 12·7% (−4·7 to 28·1) in southeast Asia under the highest emission scenario. Most of the health effects directly due to temperature increase could be avoided under scenarios involving mitigation strategies to limit emissions and further warming of the planet. Interpretation This study shows the negative health impacts of climate change that, under high-emission scenarios, would disproportionately affect warmer and poorer regions of the world. Comparison with lower emission scenarios emphasises the importance of mitigation policies for limiting global warming and reducing the associated health risks. Funding UK Medical Research Council.


Updated information on the MCC dataset
The dataset has been collected through the Multi-City Multi-Country (MCC) network, an international collaboration of research teams working on a program aiming to produce epidemiological evidence on associations between weather and health (http://mccstudy.lshtm.ac.uk/). The dataset collects time series data on environmental and health variables, together with metadata, from several locations (either cities, metropolitan areas, provinces, prefectures, or regions) in multiple countries. The current analysis includes data from 451 locations from 23 countries. The datasets from 13 countries (Australia, Brazil, Canada, China, Italy, Japan, South Korea, Spain, Sweden, Taiwan, Thailand, UK, and USA) were described in the online appendix of a previous publication. 1 Here we provide details on datasets from ten additional countries (Chile, Czech Republic, Finland, France, Ireland, Mexico, Moldova, Philippines, Switzerland, and Vietnam) and updates on those from three of the countries previously described (Canada, Spain, and UK). Canada. We collected data from 25 census metropolitan areas (CMA) and the city of Hamilton between the 1 st of January 1986 and 31 st of December 2011. A list of 20 CMA is reported in the online appendix of a previous publication. 1 The other five CMA are Niagara, Oakville, Oshawa, Sarnia, and Sault Ste Marie. Daily mortality, obtained from Statistics Canada through access to the Canadian Mortality Database, is represented by counts of deaths for all causes. Mean daily temperature (in ˚C) and relative humidity (in %), computed as the 24-hour average based on hourly measurements, were obtained from Environment Canada. A single weather station was selected for each city using the airport monitoring station located closest to the CMA centre. Measures of ozone (O3, in ppb), nitrogen dioxide (NO2, in ppb) and particles (PM2·5, in ppb) were available in the same period from the National Air Pollution Surveillance (NAPS) network of Environment Canada. Daily level of pollutants was computed as the 24-hour mean based on hourly measurements in different stations, and then averaged across stations with no missing data, with an average of 4 stations per city. In total, missing data amount for 2·00% and 0·97% of the mortality and temperature series, respectively. These data were used and described in previous publications. 2,3 Chile. We collected data from the city of Santiago de Chile between 1 st of January 2008 and 31 st of December 2014. Daily mortality, obtained from Departamento de Estadísticas e Información de Salud (Ministerio de Salud), is represented by counts of deaths for all causes. Mean daily temperature (in ˚C), computed as 24-hour average based on hourly measurements, were obtained from Sistema de Información Nacional de Calidad del Aire (SINCA), Ministerio del Medio Ambiente. A single weather station (Parque Ohiggins) was selected. In total, missing data amount for 0·00% and 1·17% of the mortality and temperature series, respectively. Czech Republic. We collected data from the 3 biggest cities (Prague, Brno, Ostrava) and a rural region (South Bohemia) between 1 st January 1994 and 31 st December 2015. Daily mortality is represented by counts of allcause deaths (obtained from the Czech Statistical Office and the Institute of Health Information and Statistics). Meteorological data (temperature and relative humidity) over the same period were obtained from stations operated by the Czech Hydrometeorological Institute (measurements in standard climatic terms 7:00, 14:00 and 21:00 local time, and daily means). For Prague, Brno and Ostrava, data from airport stations in given cities were used. Climatic variables representative for the rural region of South Bohemia was calculated as average from 3 stations (Ceske Budejovice, Kostelni Myslova, Pribyslav). In total, missing data amount for 0·00% and 0·00% of the mortality and temperature series, respectively. Finland. We collected data from the city of Helsinki between the 1 st of January 1994 and 31 st of December 2011. Daily mortality, obtained from Statistics Finland, is represented by counts of deaths for all causes. A dataset containing minimum, mean, and maximum daily temperatures was obtained from the Finnish Meteorological Institute. In this dataset, point measurements from the weather measuring stations around the country have been interpolated onto a 10×10 km grid covering the whole of Finland, using a Kriging model. The temperature variables in the Helsinki Temperature Time-series have been extracted from the GIS-database for KKJ-coordinates 6675470:2552920 (KKJ, Finnish National Coordinate System based on ED50). These are the coordinates for weather measuring station Kallion urheilukenttä of Helsinki Region Environmental Services Authority HSY. In total, missing data amount for 0·00% and 0·00% of the mortality and temperature series, respectively. France. We collected data from 18 cities (see full list in Table S1 below) between 1 st of January 2000 and 31 st of December 2010. Daily mortality, obtained from the French National Institute of Health and Medical Research (CepiDC), is represented by counts of deaths for all causes. Mean daily temperature (in ˚C), computed as the average between daily minimum and maximum, were obtained from the Meteo France. A single weather station was selected for each city; for 15 out 18 locations the weather station was located at the nearest airport. Two cities (Lille and Lens) have the same temperature series from the same weather station. In total, missing data amount for 0·00% and 0·06% of the mortality and temperature series, respectively. Ireland. We collected data in 6 regions that cover all the island of Ireland; four in the Republic of Ireland (ROI) (North East of the Republic of Ireland, South East of the Republic of Ireland, North West of the Republic of Ireland, South West of the Republic of Ireland), and two in Northern Ireland (NI) (West of Northern Ireland, East of Northern Ireland), between 1 st of January 1984 and 31 st of December 2007. Daily mortality, obtained from the Irish Central Statistics Office, and Northern Ireland Social Research Agency for ROI and NI, is represented by counts of deaths for non-external causes only (ICD-9: 0-799; ICD-10: A00-R99). Mean daily temperature (in ˚C) and relative humidity (in %), computed as the 24-hour average based on hourly measurements, were obtained from the Met Eireann, and the United Kingdom Meteorological Office for ROI and NI. Data from 10 weather station were used to build regional series for ROI, with at least two weather stations for each of the four ROI regions; while data from two weather stations were averaged to build regional series for the two NI regions. In total, missing data amount for 0·01% and 0·00% of the mortality and temperature series, respectively. Mexico. We collected data from 10 metropolitan areas (see full list in Table S1 below) between 1 st of January 1998 and 31 st of December 2014. Daily mortality, obtained from the National Institute of Statistics, Geography and Informatics is represented by counts of deaths for all causes. Temperature and relative humidity data, were obtained from Servicio Meteorológico Nacional (SMN) and the Instituto Nacional de Ecología y Cambio Climático (INECC). Daily mean temperature as well as the maximum or minimum daily temperature was calculated using hourly data, with a minimum of adequacy of information of 50%. In the case of the data from airports, these were daily averages of temperature and relative humidity. We obtained the maximum and minimum temperatures as well as the average relative humidity across all stations that met the sufficiency criteria of having at least 75% data for the day. Measures of ozone (O3, in ppb), particles (PM10, in g/m 3 ), and fine particles measures (PM2·5, in g/m 3 ) were available in the same period. Daily level of pollutants was computed as the 24-hour mean based on hourly measurements. In total, missing data account for 0·00% and 26·67% of the mortality and temperature series, respectively. Moldova. We collected data from the city of Anenii Noi (1 st  computed as the average between daily minimum and maximum, were obtained from State Hydrometeorological Service, Moldova. A single weather station was selected for each city. In total, missing data amount for 0·00% and 0·00% of the mortality and temperature series, respectively. Philippines. We collected data from four cities (Cebu, Davao, Manila, and Quezon) between the 1 st of January 2006 and 31 st of December 2010. Daily mortality, obtained from Philippine Statistics Authority-National Statistics Office, is represented by counts of deaths for all causes. Mean daily temperature (in ˚C), and relative humidity (in %), computed as the 24-hour average based on hourly measurements, were obtained from Philippine Atmospheric Geophysical and Astronomical Services Administration, and National Oceanic and Atmospheric Administration. A single weather station was selected for each city. In particular, for Davao and Cebu we used the airport monitoring station closest to the city centre, while the weather station was located in Manila Port for Manila, and in the Science Garden for Quezon. In total, missing data amount for 0·14% and 0·01% of the mortality and temperature series, respectively. Spain. We collected data from the 52 provincial capital cities (between 1 st of January 1990 and 31 st of December 2014. Daily mortality, obtained from the Spanish National Institute of Statistics (reference PB242/2016), is represented by counts of deaths for non-external causes only (ICD-9: 0-799; ICD-10: A00-R99). Mean daily temperature (in C), computed as the 24-hour average based on hourly measurements from Spain National Meteorology Agency stations were obtained through the European Climate Assessment & Dataset (available at http://www.ecad.eu/). A single weather station, located within the urban area or at the nearest airport, was selected for each city. Single-day missing values were imputed as the average of the days before and after. For periods longer than two days no imputation was done. In total, missing data amount for 0·00% and 2·00% of the mortality and temperature series after imputation, respectively. Switzerland. We collected data from 8 cities (Basel, Bern, Zurich, Geneva Lausanne, Luzern, Lugano, St Gallen) between 1 st January 1995 to 31 st December 2013. Lugano also includes the small municipalities around the main city of Lugano with similar altitude. Daily mortality, provided by the Federal Office of Statistics (Switzerland), is represented by counts of deaths for all causes, Daily data on several meteorological indicators were collected from the IDA web database (Federal Office of Meteorology and Climatology, MeteoSwiss). A single weather station located within or near the urban area was selected for each city. The meteorological indicators were: mean daily temperature (in ˚C) and relative humidity (in %), computed as the 24-hour average based on hourly measurements. Daily measurements of nitrogen dioxide (NO2, in ppb), particles (PM10, in ppb), and ozone (O3, in ppb), were provided by the Immissionsdatenbank Luft (IDB, Federal Office of the Environment, Bern, Switzerland) and computed as 24-h average for the former two and as 1h-maximum for the latter. In total, missing data amount for 0·00% and 0·00% of the mortality and temperature series, respectively.

UK.
We collected data in 9 regions of England (North East, North West, Yorkshire & Humber, East, East Midlands, West Midlands, London, South East, and South West) and Wales between the 1 st of January 1990 and 31 st of August 2012. Daily mortality, obtained from the Office of National Statistics, is represented by counts of deaths for all causes and for non-external causes only (ICD-9: 0-799; ICD-10: A00-R99). Mean daily temperature (in ˚C) computed from the 24-h average of hourly measurements) were obtained from the British Atmospheric Data Centre. An average of 29 stations contributed data to each regional series, from a minimum of 7 in London to a maximum of 44 in Wales. In total, missing data amount for 0·00% and 0·00% of the mortality and temperature series, respectively. These data were used and described in previous publications. 4,5 Vietnam. We collected data from the cities of Ho Chi Minh City (1 st of January 2010 to 31 st of December 2013), and Hue (1 st of January 2009 to 31 st of December 2013). Since 1992, a mortality data-collecting system based on the commune health center has been introduced in an official book named A6 [Vietnam Ministry of Health. Decision No 822/BYT.QD to issue mortality reporting book A6/YTCS. 1992]. Data from the A6 are collected at the commune health centre level and then forwarded to the provincial and central levels. In this study, daily mortality data from the A6 mortality reporting system, is represented by counts of deaths for all causes and for non-external causes only (ICD-9: 0-799; ICD-10: A00-R99). Mean daily temperature (in ˚C), and relative humidity (in %) computed as computed from the 24-h average of hourly measurements, were obtained from National Oceanic and Atmospheric Administration's (NOAA) National Climate Data Center (NCDC). A single weather station was selected for each city. In total, missing data amount for 0·00% and 0·45% of the mortality and temperature series, respectively.

Detailed information on the analytical framework
The data and the analytical framework are described below. Part of the statistical analysis is based on previously published methods, applied in a recent assessment on the mortality burden attributable to outdoor temperature using the same dataset. 1 Here we will focus on illustrating original extensions, while steps in common with the previous analysis are only briefly summarized, with references to published material for details. Observed series. Historical data on weather and mortality have been collected through the Multi-Country Multi-city (MCC) Collaborative Research Network (http://mccstudy.lshtm.ac.uk/). At the time of the analysis, the dataset included information from 451 locations within 23 countries aggregated in 9 regions, in largely overlapping periods ranging from 1 st of January 1984 to 31 st December 2015 ( Table 1 in the main text). The dataset is composed of observed daily time series for daily mean temperature and mortality counts for all causes or non-external causes only (International Classification of Diseases -ICD-9: 0-799; ICD-10: A00-R99), in addition to series of other weather variables. Location-specific meta-variables include weather indices derived from the series of (e.g., average and range of annual, summer and winter temperature), climatological zones based the Köppen-Geiger classification, 6 and country-specific gross domestic product (GDP) per capita. The variable was measured as the average across the 24h or between maximum and minimum daily temperature from a single or multiple monitoring stations within the administrative boundaries of each location, or in its proximity (e.g., at the nearby airport). Projected temperature series. We obtained climate projections in the form of daily temperature series for historical  and future (2006-2100) periods from a database developed in the Inter-Sectoral Impact Model Intercomparison Project (ISI-MIP, https://www.isimip.org/). 7 The ISI-MIP database includes single runs of General Circulation Models (GCMs) developed within the Coupled Model Intercomparison Project Phase 5 (CMIP5). 8 GCMs offer a representation of past, current, and projected climate consistent with scenarios of increases in radiative atmospheric forces, summarized by Representative Concentration Pathways (RCPs). Four different scenarios (RCP2·6, RCP4·5, RCP6·0, and RCP8·5), defined in the Fifth Assessment Report (AR5) of the United Nations Intergovernmental Panel on Climate Change (IPCC), 9 correspond to increasing concentration trajectories of green-house gases (GHGs) and related changes in climate. Specifically, the ISI-MIP Fast Track database provides temperature series for each RCP scenario of five GCMs, namely GFDL-ESM2M, 10,11 HadGEM2-ES, 12 IPSL-CM5A-LR, 13 MIROC-ESM-CHEM, 14 and NorESM1-M. 15,16 These five GCMs were shown to be representative of the range of projections of future climate across the CMIP5 models. 7 The model outputs were bias-corrected and downscaled through bi-linear interpolation at a 0·5°×0·5° spatial resolution and linear interpolated by day of the year. 17 The modelled daily temperature series for each of the 451 locations in the period 1990-2099 were extracted by linking the coordinates with the corresponding cells of the grid. Projected mortality series. We computed the impact of non-optimal temperature under the assumption of stable populations, thus ignoring demographic changes. For this reason, the projected series of mortality counts was computed as the average for each day of the year from , and then repeated along the same projection period of . Calibration. The five climate models (GCMs) described above are expected to reproduce the observed temperature distribution during the observed period with reasonable accuracy. However, differences between the and series can occur for several reasons, including different resolution of the data (point sources versus gridded) and poor performance of climate models in areas with sparse information from meteorological stations. 18,19 Non-negligible deviations between the two sources can produce biased results in the impact projections, when the series are applied to exposure-response relationships estimated using . For this reason, we additionally corrected using the bias-correction method that was developed and applied in ISI-MIP. 17 This approach produced * , a series recalibrated using the monthly mean and the daily variability around the monthly mean of . This calibration method ensures that the trend (i.e., the warming signal) in the original data is preserved. The accuracy of the calibration was assessed by visually comparing the cumulative distribution functions of and * for each of the 451 locations, and it was found to be very high. Estimating exposure-response relationships. We estimated the non-linear and lagged exposure-response between outdoor temperature and mortality in each location using a time series regression model with quasi-Poisson family: 20 log[ ( )] = + ( ; ) + ( ; ) + ( ; ) The bi-dimensional spline function and coefficients define an exposure-lag-response risk surface for temperature accounting for 21 days of lag. 21 These were then transformed to * and * , representing unidimensional overall cumulative (net) exposure-response curves with reduced lag dimension. 22 We controlled for seasonal and long term trends with a natural cubic spline function and for differences due to day of the week with an indicator function . Details of model specification are provided in a previous article. 1 In this analysis, we replaced the quadratic B-spline applied in the bi-dimensional function to model the temperature dimension with a natural cubic B-spline with the same knots. This choice allows a linear extrapolation of the functions beyond the boundaries of , a step needed to project the risk using * in a period characterized by a warming climate. Location-specific overall cumulative exposure-response curves obtained using the two functions were compared and did not reveal important differences. The exposure-response curves were represented as relative risk (RR), interpreted as the relative increase in risk of mortality associated with a certain temperature versus a reference value. The reference, chosen as the temperature corresponding to minimum mortality risk ( ), is interpreted as the optimal temperature for that specific location. Pooling the information across locations. We then performed a meta-regression of the reduced coefficients * , 22,23 following a model previously presented. 1 To partly account for the heterogeneity in the relationship, we included region, climate type (first letter of the Koppen-Geiger classification), 6 , GDP, and average and range of mean daily temperature across the whole period as meta-predictor. The residual heterogeneity was relatively low, with an 2 index of 39·0%, 24 similar to the original analysis including country indicators. 1 This change allows a stable estimation also in countries with a limited number of locations, now included in this extended analysis, and the partial characterization of within-country heterogeneity due to climatic and socio-economic conditions. We used the estimates from the meta-analytical model to derive the best linear unbiased prediction (BLUP) * of the coefficients in each location. This method allows areas with small daily mortality counts or short series to borrow information from larger populations that share similar characteristics 23 Projecting attributable mortality. For each day of the series of each location, we computed the number of deaths attributed to non-optimal temperature as: This daily attributable mortality was then aggregated within periods and geographical regions, and the corresponding attributable fraction was computed as the ratio with the corresponding total number of deaths. We interpret this quantity as the excess mortality related to exposure to non-optimal temperature. We also separated components due to heat and cold by summing the subsets corresponding to days with temperatures higher or lower than . 1,25 It is important to note that the method used here applies exposure-response relationships estimated using historical data to future modelled temperature series, and then computes the number of attributable deaths using a baseline mortality averaged from the observed series. Therefore, the projected temperature-related impact must be interpreted under a scenario where no adaptation (i.e., modification of the exposure response function) or population changes (i.e., modification of the baseline mortality) occur. While previous studies have reported an attenuation in the risk association with heat in several population during the last decades, 26,27 findings on cold-related mortality are less clear, 26 and the current evidence and analytical methods do not provide a well-grounded approach to quantify adaptation to future climates in terms of changes in flexible temperature-mortality relationships. Furthermore, by not considering population changes, in terms of both increase, ageing and variation in mortality rates, we isolate impacts solely related to climate change from other determinants. The figures reported in this article should be therefore interpreted as potential impacts in a well-defined scenario, and not as predictions of future excess mortality. Quantification of uncertainty. The main sources of uncertainty in the attributable mortality are related to the estimate of the exposure-response relationships and the variability in temperature projections. These quantities are represented by the covariance matrix ( * ) of the reduced model coefficients and the variability of the five series * generated in each GCM, respectively. We quantified this uncertainty by generating 1000 samples of the reduced BLUP coefficients through Monte Carlo simulations, 1,25 assuming a multivariate normal distribution for the estimated spline model coefficients, and then generating results for each of the five GCMs. We reported the results as point estimates, using the average across climate models (GCM-ensemble) obtained by the estimated coefficients, and as empirical confidence intervals, defined as the 2·5 th and 97·5 th percentiles of the empirical distribution across coefficients samples and GCMs. These empirical 95% confidence intervals (eCI) account for both sources of uncertainty. A quantitative comparison of the two components is provided by the ratio between the average standard deviation of the empirical distribution within each GCM, and the standard deviation of the average value between GCMs. Figure S1 illustrates the procedure to estimate temperature-related excess mortality in two cities located in temperate and sub-tropical areas, respectively, using a specific GMC/RCP combination. London, UK (left column), experiences an average increase of 3·4°C in * from 2010-19 to 2090-99. This increases the number of days with temperatures higher than the optimal value (minimum mortality temperature , reported as MMT in the main text), and produces temperatures exceeding those experienced in the current period. Heatrelated excess mortality rises accordingly from 0·4% to 2·7% of total deaths. However, the excess mortality due to cold, despite the decrease from 9·0% to 5·8%, remain the larger temperature-related component. The impact of climate change is more intense in the second city, Manila, Philippines (right column). The excess mortality in this location during the current period is 1·3% and 2·4%, for heat and cold, respectively. Although this city is projected to experience a lower warming (2·7°C), the shift moves almost the entire temperature distribution to levels higher than the . This eliminates any impact associated with cold, while it raises substantially the heat-related excess mortality to 14·8%. This comparison emphasises the differential health impact of global warming across locations with different characteristics, and demonstrate the various factors contributing to its quantification. In addition, the graph highlights the issues related with extrapolating the exposure-response curves beyond the observed temperature range. First, the extrapolation generates a large uncertainty, represented by the wide confidence intervals, which explains the imprecision in the measure of net impacts, especially in tropical regions with a substantial shift in temperature. Second, the log-linear assumption applied when extrapolating the spline function can underpredict potential non-linear dependencies with extremely high future temperatures.

Figure S1: Temperature and excess mortality in different climates
The graphs represent the estimated exposure-response curves with 95%CI (top panels), with minimum mortality temperature ( , continuous vertical line) used as reference, and the two portions identifying increased risks for cold and heat (blue and red, respectively). The dashed part of the curve represents the extrapolation beyond the maximum temperature observed in 2010-19 (dashed vertical line). The mid panels display the distribution of * for the current period (2010-19, grey area) and at the end of the century (2090-99, green area), projected using a specific climate model (NorESM1−M) and scenario (RCP8·5). The bottom panels represent the related distribution of excess mortality, expressed as the excess deaths (%) related to non-optimal temperatures compared with .

Figure S2: Decadal temperature trends by region and scenario
The graph shows the projected increase in temperature (°C, GCM-ensemble average), averaged by decade, region, and climate change scenario, compared to the current period .