Evidence of rapid adaptation integrated into projections of temperature-related excess mortality

Data on daily non-external mortality (ICD-9: 0-799; ICD-10: A00-R99) for 11 Spanish capital cities in the period 1978–2017 were provided by the Spanish National Institute of Statistics (see figure S1 for a map indicating city locations available in the supplementary online material at stacks.iop.org/ERL/17/044075/mmedia). Detailed information on the data collection methodology is given in INE (2013). Among the most populated cities in Spain, we selected those with time series of daily mean temperature covering the longest time period. There were no missing values in the mortality time series selected. Data of daily mean temperatures from a representative weather station in each city, corresponding to the average of maximum and minimum daily temperatures, were provided by the Spanish National Meteorology Agency (see table S1 for more information on selected weather stations). None of the temperature series had more than 0.3% of missing values (table S1). Mean summer temperatures (MSTs) were computed from the daily data as June-July-August averages.

Projections of daily mean temperature were derived from global climate model (GCM) simulations carried out in the Scenario Model Intercomparison Project (ScenarioMIP; O'Neill et al 2016) that was part of the 6th phase of the Coupled Model Intercomparison Project (CMIP6; Eyring et al 2016). The data used here comprises future simulations covering 2015–2100 from 5 GCMs (GFDL-ESM4, IPSL-CM6A-LR, MPI-ESM1-2 h, MRI-ESM2-0, UKESM1-0-LL). We considered one low-emission scenario (ssp126) and one high-emission scenario (ssp370). To avoid deviation between observed and simulated historical temperature distributions, the temperature time series were bias-adjusted using the available weather station data, employing a method developed for the third phase of the Inter-Sectoral Impact Model Intercomparison Project (ISIMIP3; Lange 2019, Lange 2021). This method is a parametric quantile-mapping method that preserves the projected temperature trends in all distribution quantiles. For the mapping, the daily mean temperature values were assumed to follow a normal distribution with seasonally dependent mean value and standard deviation.

The statistical analysis consisted of 3 stages. First, we divided the data (1978–2017) into 8 non-overlapping 5-y subperiods. For each city and each period, we fitted a quasi-Poisson regression model using distributed lag non-linear models (DLNMs) to represent the complex temperature-mortality associations. Specifications of the models were adopted from previous studies, showing their suitability for projections of temperature-related mortality (Gasparrini et al 2017, Vicedo-Cabrera et al 2019, Huber et al 2020). Specifically, we used natural cubic splines for both the exposure-response and lag-response curves, with three internal knots placed at the 10th, 75th and 90th percentile and equally spaced in the log-space, respectively. Lags of up to 21 d were accounted for. The model also included an indicator for the day of the week, and we used a natural cubic spline of time with 7 degrees of freedom per year to control for seasonal and long-term trends. We undertook a sensitivity analysis altering the specifications of the models (see table S2 for a description of the alternative modelling choices). We also tested other temporal subdivisions of the data, namely using five 8-y subperiods and four 10-y subperiods, respectively, and repeating the analysis based on 5-y subperiods for 1978–1997 and 1998–2017, separately (table S3).

In the second stage, we combined the reduced model coefficients from the first stage, describing the overall cumulative exposure-response curves in each city and each subperiod, with a multivariate meta-regression model (Sera et al 2019, 2020). This meta-regression model included an intercept only and one random-effect term at city-level. We chose not to include any meta-predictors at this stage to avoid spurious correlation in the subsequent analysis of possible association between climatic indicators and the minimum mortality temperature percentiles (MMTPs) and minimum mortality temperatures (MMTs). MMTPs and MMTs were derived based on best linear unbiased predictors (BLUPs) from the intercept-only meta-regression model, and the city- and subperiod-specific temperature distributions. We chose to use BLUPs to determine MMTPs and MMTs, because estimates based on first-stage model coefficients were unreliable due to the small sample sizes resulting from 5-y subperiods. Following (Tobías et al 2017), we used Monte Carlo simulations (n = 5000) to empirically estimate 95% confidence intervals and variances of the MMTP and MMT estimates. We restricted the searching of the MMTP to the 70th to the 99th percentiles for more robust estimates, given that for some cities and some subperiods the cumulative exposure-response curves showed two distinct minima in the lower and higher temperature range. This restriction was widened to the 1st to 99th percentiles in the sensitivity analysis. Excess mortality estimates (attributable numbers and fractions) by observational subperiod and city were computed based on BLUPs and MMTs using the forward approach described by (Gasparrini and Leone 2014).

The third stage consisted of a mixed-effect meta-regression, pooling MMT estimates for each city i and subperiod p from the previous stage. Given prior evidence on the role of mean summer temperatures (MST) in explaining adaptive processes regarding vulnerability to non-optimal temperatures (Todd and Valleron 2015, Wang et al 2018) we focused our analysis on potential associations of MMT with MST. This choice also makes sense conceptually because the MMT usually falls in the upper range of the temperature distribution (Gasparrini et al 2015), corresponding to the summer season. In order to separate longitudinal (i.e. intra-city) variability from cross-sectional (i.e. inter-city) variability we adopted a Cronbach decomposition (Cronbach and Webb 1975) with the following formula:

$$\begin{align}MM{T_{ip}} &= { }\alpha + {\beta _1}\left( {MS{T_{ip}} - MS{T_{i \cdot }}} \right)\nonumber\\ &\quad + {\beta _2}\left( {MS{T_{i \cdot }} - MS{T_{ \cdot \cdot }}} \right) + {b_i} + {\varepsilon _{ip}},\end{align} \tag{ 1 }$$

where $MS{T_{ip}} - MS{T_{i \cdot }}$ are mean summer temperatures for city i and subperiod p centred on the city mean (i.e. city-specific averages across subperiods) and $MS{T_{i \cdot }} - MS{T_{ \cdot \cdot }}$ the city means centred on the overall mean, which are included as fixed-effect predictors. The random terms ${b_i}$ represent the city effect on the estimates and they are assumed to follow a normal distribution ${b_i}\sim N\left( {0,\theta } \right)$, where $\theta$ is the between cities variance term. The cities random terms allow to consider possible non-independence of the multiple MMTs and MMTPs estimated over the different periods in each city. We assumed ${\varepsilon _{ip}}\sim N\left( {0,{S_{ip}}} \right)$, where ${S_{ip}}$ is the estimation error of MMT for each city and subperiod, computed as the variance of MMTs from the Monte Carlo simulations in the second stage of the analysis. We fitted the same type of meta-regression model to the city- and subperiod-specific MMTP estimates. Heterogeneity was assessed using the I² statistics. In a supplementary analysis, we used centred mean annual temperatures instead of centred MST as meta-predictors.

For future projections of MMTs and excess mortality, we used the last of the observational subperiods as a baseline, i.e. using estimated BLUPs, MMTs, observed temperatures and death counts from 2013 to 2017. First, for each city i, 5-y subperiod p over 2016–2100, GCM g and climate scenario s, we computed the projected minimum mortality temperatures $pMM{T_{ipgs}}$ based on the projected mean summer temperatures $pMS{T_{ipgs}}{\text{ }}$with the following equation:

$$\begin{align}pMM{T_{ipgs}} &= MM{T_{i\left( {2013 - 2017} \right)}}\nonumber\\ &\quad + {\beta _1}\left( {pMS{T_{ipgs}} - MS{T_{i\left( {2013 - 2017} \right)}}} \right),\end{align} \tag{ 2 }$$

where ${\beta _1}$ is taken from equation (1). We repeated these computations for the lower and upper value defining the 95% CI of ${\beta _1}$ to account for uncertainty in the observed MMT-MST association. The rationale of choosing ${\beta _1}$ here was that we meant to capture short-term adaptative processes, rather than the long-term adaptation reflected by the cross-sectional MMT-MST association (i.e. ${\beta _2}$ in equation (1).

In the next step, we computed daily attributable numbers based on projected temperatures in subsequent 5-y periods over 2016–2100, without adaptation (i.e. constant MMTs) and with adaptation (i.e. shifting MMTs). All of these computations were based on the overall cumulative exposure-response functions of the last observational subperiod (2013–2017), and used the average of daily death counts per day of the year from the same period to construct 5-y series of daily projected mortality. The approach taken to account for changing future MMTs in the computation of attributable mortality was to shift projected daily mean temperatures, ${T_{itgs}}{\text{ }}$, for each GCM g and each climate scenario s, by the difference between projected and observed MMTs in each future 5-y period p, according to:

$$\begin{equation}{\tilde T_{itgs}} = {T_{itgs}} - \left( {pMM{T_{ipgs}} - \;MM{T_{i\left( {2013 - 2017} \right)}}} \right).\end{equation} \tag{ 3 }$$

This is equivalent to re-centring the exposure-response functions on the projected MMTs, but easier to implement computationally. Under the assumption of no adaptation, we set $pMM{T_{ipgs}} =$ $MM{T_{i\left( {2013 - 2017} \right)}}$, thus simply using the projected temperature series.

The computations of daily attributable numbers, based on either ${\tilde T_{itgs}}$ (with adaptation) or ${T_{itgs}}$ (without adaptation) then followed the usual approach for the projection of temperature-related excess mortality (Vicedo-Cabrera et al 2019). We summed daily attributable numbers in defined 5-y periods, and specified 10-y time slices (2030s, 2050s, 2090s), separating components due to heat and cold by summing the subsets corresponding to days with ${\tilde T_{itgs}}$ (${T_{itgs}}$) higher or lower than $MM{T_{i\left( {2013 - 2017} \right)}}$. Attributable fractions were computed by dividing the resulting sums by the total number of projected deaths in the respective periods.

Uncertainties in attributable mortality was assessed by drawing Monte Carlo samples (n = 1000) from the multivariate normal distribution defined by the BLUPs and the corresponding co-variance matrices. In projections, this procedure was repeated for each scenario, GCM and, under the assumption of adaptation, each estimate of projected MMTs (central, low, high). Uncertainties are reported in terms of empirical 95% confidence intervals, which were estimated based on the 2.5th and 97.5th percentile of the resulting samples.

All computations were done using R (version 4.0.0) with packages dlnm and mixmeta. The code is available on request from the first author.

Our dataset from 11 cities in Spain included more than 3 million deaths registered over 1978–2017 (table 1). It spanned a range of summer climates with average MST of 19.6 °C (range: 18.9, 20.2) in the coolest city (Valladolid) and 27.1 °C (range: 25.6, 28.3) in the hottest city (Sevilla). On average across cities summers warmed by approximately 2 °C between the first and the last of the considered 5-y subperiods, exceeding the warming observed in annual mean temperatures (table 1).

The intercept-only meta-regression model used to derive BLUPs of the overall cumulative temperature-mortality associations (figure 1) suggested some heterogeneity across cities and subperiods (I² = 52.4%). Yet, apart from changes in the MMT, the shape of the city-specific associations remained quite similar between different subperiods (figure 1). In particular, we did not find any temporal trends in the RRs at low and high temperatures in the relative percentile scale (figure S2). The average attributable fraction aggregated across cities during 1978–2017 was 6.7% (95%CI: 4.9,8.2), with 5.8% (95%CI: 3.9,7.5) attributable to low temperatures, and 0.9% (95%CI: 0.6,1.0) attributable to high temperatures (table 1). City-specific attributable fractions (figure S3) as well as pooled attributable fractions (figure 3) varied by subperiods but did not show any consistent temporal trends over the historical study period.

Zoom In Zoom Out Reset image size

City-specific MMTs increased over time in parallel with MSTs (figure 2 solid and dashed black lines, table S4). The corresponding MMTPs ranged between the 89th and 97th percentile of the subperiod-specific daily mean temperature distributions (figure S4 and table S4). Widening the percentile range for searching of the MMTP only affected the estimates in one of the cities (Sevilla, figure S4).

Zoom In Zoom Out Reset image size

The Cronbach-type meta-regression models suggested that MMTs increased by 0.73 °C (95% CI: 0.65, 0.80) per 1 °C longitudinal (i.e. temporal) increase in the MST, and 0.84 °C (95% CI: 0.76, 0.92) per 1 °C cross-sectional (i.e. spatial) increase in the MST (table 2). For MMTP, the model indicated small decreases on the order of <1 percentile point per 1 °C increase in the MST (table 2).

Table 2. Longitudinal and cross-sectional associations between the minimum mortality temperature (MMT) and the MMT percentile (MMTP) with the mean summer temperature (MST), from Cronbach-type mixed-effect meta-regression models. For comparison, the I² statistics of the intercept-only model are reported.

	$MM{T_{ip}}$				$MMT{P_{ip}}$
	βx	95% CI	p-value	I2	βx	95% CI	p-value	I2
Intercept-only				97.6%				67.5%
Cronbach-type				43.5%				53.1%
$MS{T_{ip}} - MS{T_{i \cdot }}$	0.73	(0.65, 0.80)	<0.001		−0.75	(−1.11, −0.40)	<0.001
$MS{T_{i \cdot }} - MS{T_{ \cdot \cdot }}$	0.84	(0.76, 0.92)	<0.001		−0.55	(−0.96, −0.15)	<0.01

The sensitivity analysis showed that period-specific estimates of MMTs and MMTPs changed slightly according to specifications of the model for temperature-mortality associations (figures S5 and S6). However, coefficients of the meta-regression models describing the longitudinal and cross-sectional associations of the MMT with the MST, which we used to define adaptation in future projections, were relatively robust against changes in these model specifications (table S2). When testing alternative temporal subdivisions of the data, MMT and MMTP estimates generally fell within the 95% CIs of the default estimates (figures S7 and S8). Results from the Cronbach-type meta-regression model for MMT-MST associations were very similar to the default results, except when basing the analysis on the last 20 years rather than the full 40-y period (table S3). For MMTP, results were more strongly dependent on the choice of subperiods. Yet, with only a few exceptions, coefficients suggested a negative association of the MMTP with MST similar to our main analysis (table S3).

Increases in projected decadal averages of the MST above the last observational decade (2008–2017) did not exceed 1 °C in the low-emission scenario ssp126 (figure 2 and table 3). By contrast, in the high-emission scenario ssp370, the MST averaged across cities was projected to increase by 1.6 °C (range: 0.9, 3.9) in the 2050s and by 3.9 °C (range: 2.6, 6.0) in the 2090s relative to 2008–2017 (figure 2 and table 3). The annual average temperature increase was projected to be only slightly lower (table 3).

Without adaptation, future heat-related excess mortality for all cities combined was projected to increase, reaching 1.9% (95%CI: 1.1, 3.5) by the 2050s and 4.8% (95%CI: 2.4, 9.0) by the 2090s under the high-emission scenario ssp370 (figure 3 and table 3). Under the same scenario, cold-related excess mortality was projected to decrease to 4.6% (95%CI: 2.7, 6.4) by the 2050s and 3.4% (95%CI: 1.5, 5.3) by the 2090s (figure 3 and table 3). Total temperature-related excess mortality was estimated to change little until mid-century and then to increase to 8.2% (95%CI: 5.6, 11.6) by the 2090s. Under the low-emission scenario ssp126, without adaptation, attributable fractions were estimated to stay within the ranges observed over 1978–2017, except for a relatively small projected increase in heat-related excess mortality by the 2050s (figure 3 and table 3).

Zoom In Zoom Out Reset image size

Assuming adaptation according to past MMT shifts (figure 2) strongly impacted the projected aggregated excess mortality under ssp370. Heat attributable fractions were 36.9% (95%CI: 21.3, 64.1) lower by the 2050s and 63.5% (95%CI: 50.0, 81.2) lower by the 2090s with adaptation relative to the estimates without adaptation (figure 3 and table 3). Cold attributable fractions were 20.7% (95%CI: 7.6, 60.3) higher by the 2050s and 63.7% (95%CI: 30.2, 166.7) higher by the 2090s. The relative percent differences with and without adaptation for total temperature-related excess mortality were comparatively small, with the confidence interval even for the 2090s estimate (−11.2%; 95%CI: −39.5, 5.5) including zero (figure 3 and table 3). Under ssp126, assuming adaptation in projections changed excess mortality estimates compared to projections without adaptation in the same direction as for ssp370 (figure 3). However, the magnitudes of these changes were small, and all confidence intervals included zero (table 3). For the 2030s, estimates of attributable fractions were very similar between low and high-emission scenarios, and adaptation assumptions had minor effects in both scenarios (figure 3 and table 3). City-specific results were qualitatively the same as results for all cities combined (figures S9–S11).