We aggregate the hourly gridded 2-m temperature data from the high–resolution ERA5-Land reanalysis7 into weekly regional averages of daily mean 2-m temperature for the period of 1950–2022. We obtain global mean surface temperature (GMST) changes from HadCRUT59. We use the same mortality data as in ref. 5: We obtained weekly counts of all-cause mortality by sex and age groups from Eurostat19, and missing data was complemented by contacting the corresponding national agencies for statistics. The final dataset includes 45,184,044 counts of death (22,000,519 for women and 21,913,050 for men) between January 2015 and November 2022 from 823 contiguous regions inhabited by over 543 million Europeans in 35 countries. All–age data by sex was not available in the United Kingdom, and only at the country level in Germany. Data by sex and age groups were not available in the United Kingdom, Ireland, and Germany.
Counterfactual temperatures
For each region, we linearly regress the annual time series of regional summer mean temperature onto the annual time series of GMST for the 4-year running mean from HadCRUT5.
The slope parameter for each region is then multiplied by the value of GMST in 2022 relative to the pre-industrial era (i.e. +1.15 °C9), and this regional summer warming is then subtracted from the regional factual temperatures to obtain the regional counterfactual temperatures. We also generate counterfactual scenarios for the upper and lower bound of the 95% confidence intervals of the regional anthropogenic warming.
It has recently been shown that in much of Western Europe, circulation changes have considerably exacerbated the largely thermodynamic warming that primarily occurs in response to anthropogenic greenhouse gas emissions20,21,22. As of now, it remains unclear if these circulation changes are entirely unforced, i.e. a manifestation of internal climate variability, or partly externally forced20.
Bearing all this in mind, we deem our simple approach to derive counterfactual climates reasonable because, in our view and considering recent research on European warming and its drivers, the potential benefits of more complex detrending methods cannot be validated.
Epidemiological analysis
We use the same two-step epidemiological model as in ref. 5. In the first stage, we used quasi-Poisson regression models, which allow for overdispersed counts of deaths, to calculate the location-specific temperature–lag–mortality relation in each European region11,23. The models include (i) an intercept, (ii) a natural cubic spline of time with 8 degrees of freedom per year to control for the seasonal and long-term trends, and (iii) a cross-basis function to estimate the exposure–lag–response association between weekly temperatures (temp) and mortality counts (mort):
$$\begin{array}{ll}{log}(E({{{mort}}}))={{{intercept}}}+{ns}({{{time}}},8{{{df}}\; {{per}}\; {{year}}})\\\qquad\qquad\qquad\quad+\,{{{cross-basis}}}({{{temp}}}{{;}}\,0,1,2,3\,{{{weeks}}})\end{array}$$
The lag–response function of the cross-basis is modelled with integer lag values of 0, 1, 2, and 3 weeks, and the exposure–response function with a natural cubic spline with three internal knots at the 10th, 50th, and 90th percentiles of the location-specific weekly temperature distribution.
In the second stage, we use a multivariate multilevel meta-regression analysis24 to pool the location-specific coefficients obtained in the first step. The meta-regression includes (i) country random effects and the location-specific (ii) temperature average, (iii) temperature interquartile range, and (iv) percentage of people aged 80+ years as meta–predictors13. We derived the best linear unbiased predictions of the temperature–mortality relationship in each region from the meta-regression25 to obtain the location-specific minimum mortality temperature and to transform the regional temperature and mortality time-series from January 2015 to November 2022 into the weekly and summer heat-related mortality numbers over the years 2015–202226. Heat-related mortality is calculated for the weeks with average temperatures above the location-specific minimum mortality temperature27. Regional heat-related mortality numbers are aggregated to obtain the national and European burdens28,29. Similarly, we computed 1000 Monte Carlo simulations of the regional heat-related mortality numbers and separately aggregated the numbers in each simulation to calculate the 95% confidence intervals at the national and continental levels11,13. We then follow the methodology proposed by Ballester et al.10 to bias-correct the weekly estimates by using a dataset of daily mortality data spanning over 16 European countries.