capbtabboxtable[][] scaletikzpicturetowidth[1]

# The association between environmental variables and short-term mortality:
evidence from Europe

Jens Robben  Faculty of Economics and Business, KU Leuven, Belgium.  LRisk,
Leuven Research Center on Insurance and Financial Risk Analysis, KU Leuven,
Belgium.  Corresponding author<EMAIL_ADDRESS>Katrien Antonio
Faculty of Economics and Business, KU Leuven, Belgium.  Faculty of Economics
and Business, University of Amsterdam, The Netherlands.  LRisk, Leuven
Research Center on Insurance and Financial Risk Analysis, KU Leuven, Belgium.
LStat, Leuven Statistics Research Center, KU Leuven, Belgium.  RCLR, Research
Centre for Longevity Risk, University of Amsterdam, The Netherlands.  Torsten
Kleinow  Faculty of Economics and Business, University of Amsterdam, The
Netherlands.  RCLR, Research Centre for Longevity Risk, University of
Amsterdam, The Netherlands.

(Version: )

###### Abstract

Using fine-grained, publicly available data, this paper studies the
association between environmental factors, i.e., variables capturing weather
and air pollution characteristics, and weekly mortality rates in small
geographical regions in Europe. Hereto, we develop a mortality modelling
framework where a baseline captures a region-specific, seasonal historical
trend observed within the weekly mortality rates. Using a machine learning
algorithm, we then explain deviations from this baseline using anomalies and
extreme indices constructed from the environmental data. We illustrate our
proposed modelling framework through a case study on more than 550 NUTS 3
regions (Nomenclature of Territorial Units for Statistics, level 3) located in
20 different European countries. Through interpretation tools, we unravel
insights into which environmental features are most important when estimating
excess or deficit mortality with respect to the baseline and explore how these
features interact. Moreover, we investigate harvesting effects of the
environmental features through our constructed weekly mortality modelling
framework. Our findings show that temperature-related features exert the most
significant influence in explaining deviations in mortality from the baseline.
Furthermore, we find that environmental features prove particularly beneficial
in southern regions for explaining elevated levels of mortality over short
time periods.

##### Keywords:

weekly mortality modelling; high-resolution gridded datasets; environmental
data

##### Data and code availability statement:

The datasets used in this paper are publicly available. We consulted the
deaths by week, sex, 5-year age group and NUTS 3 region on Eurostat
(https://ec.europa.eu/eurostat/databrowser/view/demo_r_mweek3/, downloaded on
January 20, 2024). We retrieved the weather-related variables from the
Copernicus Climate Data Store (https://cds.climate.copernicus.eu/, downloaded
on January 22, 2024) and the environmental-related variables from the
Copernicus atmospheric monitoring service
(https://ads.atmosphere.copernicus.eu/, downloaded on January 22, 2024). The
code used for the implementation and analysis of the case study presented in
this paper will be made available on the GitHub repository
https://github.com/jensrobben/EnvVar-Mortality-Europe. The code is written in
R and can be accessed and downloaded for further reference and replication of
the obtained results.

##### Funding statement:

Katrien Antonio acknowledges funding from the FWO and F.R.S.-FNRS under the
Excellence of Science (EOS) programme, project ASTeRISK (40007517). The FWO
WOG network W001021N is acknowleged as well. Katrien Antonio acknowledges
support from the Chair of Excellence on Digital Insurance And Long‑term risk
(DIALog) by CNP Assurances. This study is furthermore part of the research
programme at the Research Centre for Longevity Risk. RCLR is a joint
initiative of NN Group and the University of Amsterdam, with additional
funding from the Dutch government’s Public-Private Partnership (PPP)
programme.

##### Conflict of interest disclosure:

The authors declare no conflict of interest.

##  1 Introduction

Extreme environmental conditions increasingly impose challenges to the life
insurance industry, a concern underscored by EIOPA (EIOPA, 2021), highlighting
the importance of including climate change scenarios in an insurer’s own risk
and solvency assessment, both in the short and long-term. This importance is
underscored by the more frequent, intense, and prolonged heat waves that are
expected to occur in the 21st century (Meehl Tebaldi, 2004). In the
epidemiological and medical literature, studies have consistently unveiled a
direct association between temperature and mortality statistics (Basu Samet,
2002). Epidemiological evidence also underscores the adverse effects of air
pollution on health (Brunekreef Holgate, 2002; Rückerl ., 2011). Understanding
the impact of environmental conditions on mortality is particularly difficult
because of the interplay between weather factors and air pollution (Orru .,
2017). In the actuarial literature, incorporating environmental factors in
mortality modelling is still largely unexplored. This paper presents an
approach to fill this gap and aims to assist actuaries and risk managers to
better understand and quantify environmental-related risks for life and health
insurance policies. Additionally, the proposed method facilitates the
assessment of the impact of future environmental scenarios on mortality
outcomes, valuable for insurers conducting their risk and solvency
assessments.

The objective of this paper is to identify the main environmental factors
influencing excess or deficit mortality with respect to a baseline mortality
model as registered on a weekly, regional level within Europe. We also focus
on detecting relevant interaction effects between air pollution and weather
factors and their association with mortality rates. Since actuaries typically
work with annual and country-level mortality estimates, this paper further
quantifies the added contribution of environmental factors when aggregating
the resulting estimated weekly mortality rates on an annual or country-level
basis. Lastly, the paper explores how future paths for weekly mortality rates
can be generated, reflecting the impact of specific environmental scenarios on
mortality outcomes.

The proposed framework for weekly mortality modelling comprises two building
blocks. The first building block is a weekly, region-specific baseline
mortality model, which captures the overall seasonal trend observed in the
historical weekly mortality rates of a given region. To achieve this, we
employ a Serfling-type weekly mortality model, incorporating sine and cosine
Fourier terms (Serfling, 1963). The second building block aims to analyse
mortality deviations from the baseline model as a function of region-specific
environmental anomalies and extreme environmental indices. Environmental
anomalies measure deviations from a time-dependent, baseline level, while
extreme environmental indices measure the frequency of how often the
environmental factor surpasses a specified high quantile or falls below a
certain low quantile within a week. This second step allows us to unravel the
associations between the observed deviations from the mortality baseline and
environmental features capturing events like heat waves, cold spells, or
heightened levels of air pollution. Due to the complex interplay between
weather and air pollution features and their impact on mortality, we adopt a
machine learning approach to learn and estimate these intricate, non-linear
interaction effects in a data-driven way. By using interpretation tools for
machine learning models, we investigate our research objectives.

The investigation of the impact of environmental features on mortality is
still largely unexplored in the actuarial literature on mortality modelling.
Among the scant literature in actuarial journals, Li Tang (2022) delve into
the modelling of joint extremes in monthly temperature and death counts using
a bivariate peaks-over-threshold approach. Through an empirical study using
monthly temperature data and death counts in the United States, they reveal
that joint extremes in cold weather and old-age death counts demonstrate the
strongest level of dependence. Additionally, Dong . (2022) quantify the impact
of air quality on monthly mortality rates. Their findings indicate that
incorporating levels from the air pollutant PM2.5, i.e., fine particulate
matter consisting of particles that are 2.5 microns or less in diameter,
enhances the performance of estimating excess mortality. Dong . (2024) assess
the impact of climate risk stress scenarios on life insurance portfolios.
Lastly, the American Academy of Actuaries, the Casualty Actuarial Society, the
Canadian Institute of Actuaries, and the Society of Actuaries developed the
Actuaries Climate Index (ACI), a quarterly measure derived from changes in
extreme weather events and sea levels, aiming at offering a practical
monitoring tool for tracking climate trends (American Academy of Actuaries .,
2016).

In Europe, the EuroMOMO project has emerged as an initiative for monitoring
weekly excess mortality across 28 participating European countries or sub-
national regions.111The website can be consulted here:
https://www.euromomo.eu/. EuroMOMO, known as the European network for
MOnitoring excess MOrtality, estimates weekly all-cause excess deaths in a
standardized manner. EuroMOMO’s objective is to coordinate timely mortality
monitoring and analysis across European countries using a standardized
approach for swift detection of excess mortality signals, particularly related
to pandemics, influenza, emerging infections, and extreme environmental
conditions. Hereto, it uses the Serfling model as mortality baseline when
measuring excess mortality (Serfling, 1963), which consists of a quasi-Poisson
regression model on weekly death counts with a log-linear long-term trend and
Fourier terms to address the seasonality inherent in the weekly mortality
pattern. Within the EuroMOMO network, Nielsen . (2018) introduced the FluMOMO
model, an extension of the weekly Serfling-type baseline mortality model which
adds influenza activity and extreme temperatures as explanatory variables to
the model. Widely utilized in the literature, this model serves as a tool for
calculating excess deaths related to various health risks such as flu periods,
heat waves, and pandemics, e.g., the COVID-19 pandemic. For instance, Nielsen
. (2019) estimate the weekly excess mortality attributable to all causes and
influenza during the season 2017-2018 in Europe. Vestergaard . (2020) present
pooled European-wide weekly mortality estimates and excess all-cause mortality
during the COVID-19 pandemic using the EuroMOMO modelling framework.
Furthermore, Schöley (2021) investigates and describes various procedures for
computing weekly mortality baseline models in order to evaluate excess deaths
during the COVID-19 pandemic. Notably, he finds that the popular EuroMOMO-
style Serfling model to assess expected death counts reports one of the
highest excess deaths among the models considered in his study.

Various methodologies have been explored in the epidemiological and medical
literature to investigate the association between environmental factors and
(daily) mortality rates or death counts. One approach involves an
observational study by Keatinge . (2000), who investigate heat-related
mortalities across Europe by analysing region-specific temperatures
corresponding to the lowest daily death rates. Another strand of research
relies on time series regression models, assuming an overdispersed Poisson
distribution for the daily death counts with temperature or a particular air
pollutant as an explanatory variable, along with other confounding factors
(Armstrong, 2006). While earlier studies favoured a linear relationship
between daily deaths and the considered environmental level, recent research
suggests non-linear associations (Braga ., 2002). Braga . (2001) and Pattenden
. (2003) highlight the delayed effects of cold spells and air pollution on
mortality, necessitating more advanced modelling techniques like distributed
lag models (DLMs). These models, exemplified by Schwartz (2000), accommodate
delayed effects, crucial in understanding mortality patterns following extreme
weather events or air pollution spikes. These have been further extended to
Distributed Lag Non-linear Models (DLNMs), offering greater flexibility in
capturing non-linear effects across both predictor space and lag dimensions
(Gasparrini ., 2010). For a more detailed and technical overview of these
methodologies, we refer to Appendix A.

Our paper contributes to the literature in four ways. First, rather than
focusing solely on the impact of an individual environmental factor on
mortality, as is typically done in the epidemiological and medical literature,
we investigate the combined effects of several environmental factors on
mortality rates. This provides a more comprehensive understanding of the joint
impact of environmental factors on mortality. Second, we utilize a machine
learning model with environmental anomalies and extreme environmental indices
as inputs to analyse mortality deviations from the baseline. Therefore, we do
not impose any functional form regarding the effect of the covariates on
mortality rates, enabling us to better capture the complex interactions
present in the environmental data. Third, we leverage fine-grained open data
from the Copernicus Climate Data Store (CDS) for weather factors and the
Copernicus Atmospheric Monitoring Service (CAMS) for air pollution levels.
Particularly, we use the E-OBS daily gridded meteorological data for Europe
from the CDS (CIESIN et. al, 2018) and the European air quality reanalysis
from the CAMS (INERIS et. al, 2022). These datasets offer high temporal and
spatial resolutions, providing detailed information on environmental factors
across Europe. Our paper stands out as one of the pioneering works to utilize
these datasets within the actuarial literature. By detailing how this raw,
fine-grained data can be used to create environmental anomalies and extreme
environmental indices, defined on a particular geographic region and temporal
scale, we offer valuable insights for researchers and practitioners aiming to
integrate such data into their mortality modelling framework. Fourth, our
paper extends the EuroMOMO mortality methodology by integrating environmental
factors into their proposed weekly mortality baseline model.

This paper is organized as follows: Section 2 provides an overview of the
notations and data used in this paper. In Section 3, we detail our mortality
modelling framework. Specifically, Section 3.1 introduces a weekly, region-
specific mortality baseline model, while Section 3.2 explains our approach to
analyse mortality deviations from this baseline using environmental features
by means of a machine learning technique. We outline the calibration strategy
for both the weekly baseline mortality model and the machine learning model in
Section 4. To illustrate the practical application of our methodology, we
conduct a case study in Section 5, applying our weekly mortality modelling
framework incorporating environmental factors to the NUTS 3 regions of 20
European countries. Section 5.1 details the feature engineering process to
create environmental anomalies and extreme environmental indices on a weekly
time scale and NUTS 3 geographical level, from the raw, fine-grained data. In
Section 5.2, we discuss the calibration results of our weekly mortality
modelling framework applied to our case study, while Section 5.3 provides
further insights and applications. Section 6 provides concluding remarks and
summarizes the key findings of our study.

##  2 Notations and data

###  2.1 Notations

Let $d_{x,t,w}^{(r)}$ be the observed death count in region $r$ at age group
$x$ during ISO week $w$ of ISO year $t$.222We follow the ISO 8601 standard
maintained by the International Organization for Standardization, see
https://www.iso.org/standard/70907.html. We denote the set of age groups or
buckets under consideration as $\mathcal{X}$, the range of ISO years as
$\mathcal{T}$ and the set of ISO weeks in ISO year $t$ as $\mathcal{W}_{t}$.
Particularly, the number of ISO weeks in ISO year $t$ is 52 or 53.333A year
$t$ consists of 53 ISO weeks if January 1st of that year falls on a Thursday,
or if it is a leap year with January 1st being a Wednesday. Furthermore, we
encode each region by an integer, i.e., we denote the set of regions as
$\mathcal{R}=\\{1,2,\ldots,R\\}$. For the sake of brevity in notation, we omit
explicit reference to gender, although the methodologies described in this
paper apply to male, female, and unisex data. Moreover, we denote
$E_{x,t,w}^{(r)}$ for the exposure-to-risk in region $r\in\mathcal{R}$ at age
group $x\in\mathcal{X}$ during ISO week $w\in\mathcal{W}_{t}$ of ISO year
$t\in\mathcal{T}$.

The region-specific weekly force of mortality, denoted by $\mu_{x,t,w}^{(r)}$,
represents the instantaneous rate of mortality in region $r$ at age group $x$
and in ISO week $w$ of ISO year $t$. We define the observed weekly death rate,
i.e., $m_{x,t,w}^{(r)}$, as:

$\displaystyle m_{x,t,w}^{(r)}=\dfrac{d_{x,t,w}^{(r)}}{E_{x,t,w}^{(r)}}.$

Furthermore, $\smash{q_{x,t,w}^{(r)}}$ refers to the mortality rate and
represents the probability that an individual from region $r$ and age group
$x$ and who is alive at the start of ISO week $w$ in ISO year $t$, will die
within the next week. Under the assumption that the force of mortality
$\smash{\mu_{x,t,w}^{(r)}}$ is constant within a week, we can estimate the
weekly mortality rate as (Pitacco, 2009):

$\displaystyle q_{x,t,w}^{(r)}\approx 1-\exp\left(-\mu_{x,t,w}^{(r)}\right),$
(2.1)

for each $x\in\mathcal{X}$, $t\in\mathcal{T}$, $w\in\mathcal{W}_{t}$ and
$r\in\mathcal{R}$.

###  2.2 Data sources: death counts, exposures and environmental data

##### Death counts.

We extract the deaths by week, sex, 5-year age group and NUTS 3 region
throughout the years 2013-2019 from Eurostat for the following 20 European
countries:444 We concentrate on the time period from 2013 to 2019 due to the
availability of air pollution data, as introduced later on in this section.
Austria (AT), Belgium (BE), Switzerland (CH), Czech Republic (CZ), Denmark
(DK), Estonia (EE), Spain (ES), Finland (FI), France (FR), Hungary (HU), Italy
(IT), Liechtenstein (LI), Lithuania (LT), Luxembourg (LU), Latvia (LV), Norway
(NO), Poland (PL), Portugal (PT), Sweden (SE), and Slovakia (SK).555Countries
such as Germany, the Netherlands, and the United Kingdom are excluded from our
analysis due to the unavailability of weekly death counts at the NUTS 3 level
for the years 2013-2019. The database can be consulted on
https://ec.europa.eu/eurostat/databrowser/view/demo_r_mweek3/default/table?lang=en.
This database uses the ISO-time format for weeks and years. The Nomenclature
of Territorial Units for Statistics (NUTS) serves as a geocode standard used
to designate the administrative divisions of countries for statistical
purposes. The NUTS 3 level is the lowest level in the NUTS classification and
corresponds for example to the departments in France (101 regions) and
provinces and metropolitan cities in Italy (107 regions). We do not consider
the overseas regions in Spain, Italy, and Norway in our analysis. Figure 1
visualizes the NUTS 3 regions in the 20 European countries in our study. In
the remainder of this paper, we will often use the terms ‘regions’ instead of
‘NUTS 3 regions’ and ‘years’ and ‘weeks’ instead of ‘ISO years’ and ‘ISO
weeks’, respectively, for the sake of conciseness.

In this paper, we specifically concentrate on older ages due to their
increased vulnerability to the impacts of environmental factors. This approach
aligns with previous studies, such as Keatinge . (2000) and Sunyer . (1996).
Hereto, we aggregate the death counts $\smash{d_{x,t,w}^{(r)}}$ per gender,
week, year, and region, and construct one age group consisting of the ages 65
and older, denoted in the sequel of this paper as $\smash{d_{t,w}^{(r)}}$ for
notational convenience. Figure 1 shows the death counts at age group 65+
across the NUTS 3 regions in the 15th week of the year 2015. We observe three
regions with significantly higher observed death counts, aligning with the
three most populated regions: Madrid, Barcelona, and Rome. Figure 1
illustrates the weekly death counts spanning the years 2013-2019 across the
NUTS 3 regions of Barcelona (ES511), Milano (ITC4C), and Stockholm (SE110).
The figure reveals a seasonal trend, with more deaths during winter weeks and
relatively lower death counts during the summer weeks.

()

()

()

Figure 1: Panel (a) shows the NUTS 3 regions in 20 European countries,
excluding overseas regions. Panel (b) illustrates the female death counts in
the NUTS 3 regions at the 15-th ISO week of the year 2015. Panel (c)
visualizes the weekly female death counts throughout the years 2013-2019 in
three NUTS 3 regions in Europe: Barcelona (red), Milano (green), and Stockholm
(blue).

##### Exposure-to-risk.

In order to calculate weekly mortality rates, it is necessary to formulate a
weekly exposure measure, denoted as $\smash{E_{t,w}^{(r)}}$, for the age group
65+ in each NUTS 3 region throughout the years 2013-2019. Hereto, we retrieve
the population count of people aged 65+ on January 1 of each year $t$ by sex
and NUTS 3 region from Eurostat.666The database can be accessed on the web
page:
https://ec.europa.eu/eurostat/databrowser/view/demo_r_pjanaggr3/default/table?lang=en.
We construct an estimate for the weekly exposure similar to the approach
proposed in Jdanov . (2021), and obtain:

$\displaystyle E_{t,w}^{(r)}=\dfrac{P_{t,0}^{(r)}+P_{t+1,0}^{(r)}}{2\cdot
52.18},$ (2.2)

for $t\in\\{2013,..,2019\\}$, with $P_{t,0}^{(r)}$ the population count of
people aged 65+ at January, $1$ of year $t$ in NUTS 3 region
$r\in\mathcal{R}$.777The subscript 0 refers to the start, i.e., January 1, of
year $t$. The number $52.18$ represents the average number of weeks per
year.888While the current approach assumes a constant weekly exposure within a
year, a more advanced method could employ linear interpolation for a smoother
transition of weekly exposures between consecutive years.

##### Weather data.

To retrieve weather factors, we consult the Copernicus Climate Data Store
(CDS). The CDS provides access to a wide range of weather-related data,
including observations, reanalysis, and projections with different spatial and
temporal resolutions.999Reanalysis is the process of using climate models to
combine observational data with historical climate information to produce
comprehensive time series of climate variables on a (sub)daily time scale.
These datasets have been shown to be a good alternative to weather station
data in temperature-mortality assessments across different regions (de
Schrijver ., 2021).

Particularly, we use the E-OBS land-only, gridded meteorological data for
Europe from 1950 onwards, derived from in-situ observations, see Copernicus
Climate Change Service, Climate Data Store (2020). The E-OBS dataset is
constructed using observations gathered by meteorological stations across
Europe and provided by the National Meteorological and Hydrological Services
(NMHS) and other affiliated institutes. It is a daily, high-resolution gridded
dataset, where weather factors are defined on grids with a spatial resolution
of $0.10^{\circ}$ ($\approx 11$ km) in both longitude and latitude direction.

Table 1 lists the weather factors that we use in this paper, i.e., we extract
the daily maximum temperature Tmax, the daily average temperature Tavg, the
daily minimum temperature Tmin, the daily average relative humidity Hum, the
total daily precipitation Rain, and the daily average wind speed Wind. We opt
for these weather factors due to their relevance in epidemiological and
medical research, particularly concerning their impact on mortality. For
instance, Barnett (2010) finds that the daily minimum, average, and maximum
temperatures all have a comparable level of predictive power for explaining
daily death counts. Findings from studies such as Braga . (2002) and Alberdi .
(1998) did not find a significant impact of humidity and wind speed,
respectively, on daily mortality. Nevertheless, we explore these weather
factors for verification purposes.

Weather factor |  Explanation
---|---
Tmax |  Daily maximum air temperature, measured in degrees Celsius (∘C) at 2 meters above surface.
Tavg |  Daily average air temperature, measured in ∘C at 2 meters above surface.
Tmin |  Daily minimum air temperature, measured in ∘C at 2 meters above surface.
Hum |  Daily average relative humidity, measured at 2 meters above surface. Relative humidity is defined as the percentage of actual humidity to saturation humidity and falls within the range of $[0,100]$. Saturation humidity occurs when the air contains the maximum amount of moisture possible at a specific temperature.
Rain |  The daily total precipitation expressed in $mm$, including rain, snow, and hail, and measured as the height of equivalent liquid water per square meter.
Wind |  Daily average wind speed in $\nicefrac{{m}}{{s}}$, measured at 10 metres above the surface.
Table 1: Weather variables retrieved from the E-OBS gridded meteorological
dataset on the CDS.

We extract the E-OBS daily weather factors in Europe across a grid that covers
the 20 European countries highlighted in Figure 1, and spanning the years
2013-2019. As an example, Figures 2 and 2 display, respectively, the daily
maximum temperature and the total precipitation amount at two randomly
selected dates within the considered time range 2013-2019. Figure 2
illustrates the maximum temperature on August 2 on a timeline from 2013-2019
to show the fluctuation in the maximum temperature on that particular day
across the historical period under consideration.

()

()

Aug. 2, 2013 Aug. 2, 2014 Aug. 2, 2015 Aug. 2, 2016 Aug. 2, 2017 Aug. 2, 2018
Aug. 2, 2019

()

Figure 2: Panel (a) shows the maximum temperature on August 2, 2015, panel (b)
the total precipitation amount on February 13, 2014, and panel (c) the maximum
temperature on August 2 in the years 2013-2019, in Europe.

##### Air pollution data.

To access air pollution measurements, we consult the Atmosphere Data Store
(ADS), a platform that facilitates the implementation of the Copernicus
Atmosphere Monitoring Service (CAMS). This monitoring service offers
comprehensive, historical, and real-time data on the atmospheric composition,
such as air quality, greenhouse gases, and aerosols.

We rely on the CAMS European air quality reanalyses dataset (INERIS et. al,
2022), which delivers hourly air quality reanalyses for Europe based on
validated observations for the years 2013-2020. This reanalysis integrates
both model data and observations from the European Environment Agency (EEA)
using data assimilation techniques, with the inclusion of satellite data to
enhance in-situ observations. Similar to the E-OBS data, this air quality
dataset is a high-resolution gridded dataset. Measurements of several air
pollutants are available on a fine grid with a spatial resolution of
$0.10^{\circ}$ ($\approx 11\text{ km}$) in both longitude and latitude
directions and a temporal resolution of one hour. Table 2 lists the air
pollutants that we will use in this paper, i.e., O3, NO2, PM10, and PM2.5. We
select these air pollutants because of their demonstrated significant short-
term impact on mortality in the literature. Pascal . (2014) identify short-
term impacts of PM10 and PM2.5 levels on daily mortality in nine French
cities. In addition, Orellano . (2020) find short-term impacts of ozone and
nitrogen dioxide on daily mortality based on a systematic review including 196
studies.

Air pollutants |  Explanation
---|---
O3 |  Hourly ozone levels, measured in micrograms per cubic meter ($\mu g/m^{3}$).
NO2 |  Hourly nitrogen dioxide levels ($\mu g/m^{3}$).
PM10 |  Hourly particular matter levels of particles with a diameter of 10 microns or less ($\mu g/m^{3}$).
PM2.5 |  Hourly fine particular matter levels of particles with a diameter of 2.5 microns or less ($\mu g/m^{3}$).
Table 2: Air pollutants retrieved from the CAMS European air quality
reanalyses dataset.

We gather hourly CAMS air pollution measurements in Europe for the years
2013-2019 and focus on the air pollutants in Table 2. Unlike the E-OBS
dataset, this data covers both land and sea. To illustrate, Figures 3 and 3
present the O3 and PM10 concentrations, respectively, at two randomly selected
times. We highlight the boundaries of the 20 European countries in red to
provide a clear visualization of the land component of interest.

()

()

Figure 3: Panel (a) visualizes the ozone concentration on September 7, 2016 at
8:00 AM and panel (b) the PM10 concentration on February 1, 2019 at 10:00 AM.
We highlight the boundaries of the 20 European countries in red

##  3 Model specification

The proposed weekly mortality model aims to explain excess mortality or
mortality deficit, relative to a baseline mortality model, from environmental
data. Hereto, the baseline mortality model estimates a weekly seasonal
mortality pattern observed in the weekly mortality rates of each considered
NUTS 3 region. As such, we obtain an estimate for the expected number of death
counts in each week of the year and for each region, see Section 3.1.
Mortality excess may arise from abnormal temperatures, elevated air pollution,
or extreme humidity levels, among others. Therefore, we also engineer features
that quantify deviations from typical environmental conditions for each week
of the year, referred to as anomalies. Furthermore, to allow for, e.g., heat
waves or cold spells, we create extreme environmental indices that indicate
how many days within a week a particular environmental feature exceeds a high
quantile or is lower than a low quantile. We then use these environmental
anomalies and extreme indices, along with lagged versions of them, to explain
excess or deficit mortality relative to the mortality baseline, see Section
3.2. Section 5.1 provides a detailed discussion of the feature engineering
step outlined above. To capture the complex interaction effects between the
environmental features and their possibly non-linear impact on mortality, we
employ a machine learning model. This offers the additional advantage of
automating the feature selection process in a dataset with numerous features.
Figure 4 illustrates the proposed methodology.

Feature engineering $t$

Figure 4: Visual representation of the proposed weekly mortality model,
consisting of a baseline mortality model (left) and a machine learning model
(right). The baseline model is region-specific and captures the overall
seasonal pattern in the observed mortality rates. Next, we explain the
mortality deviations from this region-specific baseline model with a machine
learning algorithm using a set of engineered region-specific environmental
features.

###  3.1 The baseline weekly mortality model

The weekly, region-specific baseline mortality model estimates the overall,
seasonal trend observed in the region’s weekly death counts, see Figure 1 for
an example of this seasonality. We adhere to the framework proposed by
Serfling (1963), incorporating seasonality through Fourier terms. We exclude
the dependence on the age group $x$ since our focus will be on the age group
65+, as motivated in Section 2.1. In this paper, we include population-
exposures as an offset into a Poisson regression model. Hereto, we assume that
the observed weekly death counts in region $r$ are realizations from a Poisson
distributed random variable $D_{t,w}^{(r)}$:

$\displaystyle
D_{t,w}^{(r)}\sim\text{Poisson}\left(E_{t,w}^{(r)}\>\mu_{t,w}^{(r)}\right),$
(3.1)

where we consider, for each region $r\in\mathcal{R}$, the following structure
for $\mu_{t,w}^{(r)}$ (Serfling, 1963; Schöley, 2021):

$\displaystyle\log\mu_{t,w}^{(r)}=$
$\displaystyle\beta_{0}^{(r)}+\beta_{1}^{(r)}t+\beta_{2}^{(r)}\sin\left(\frac{2\pi
w}{52.18}\right)+\beta_{3}^{(r)}\cos\left(\frac{2\pi w}{52.18}\right)+$ (3.2)
$\displaystyle\>\>\>\beta_{4}^{(r)}\sin\left(\frac{2\pi
w}{26.09}\right)+\beta_{5}^{(r)}\cos\left(\frac{2\pi w}{26.09}\right),$

where $52.18$ is the average number of weeks per year.101010Alternatively, the
exact number of ISO weeks in the year of consideration can be used. Note that
this baseline structure is region-, week-, and year-specific. Further, the
region-specific parameters $\smash{\beta_{p}^{(r)}}$, for $p=0,1,..,5$,
control for the level, trend, and seasonal variation. In Section 4.1, we show
how to estimate this baseline model. We estimate the expected death counts
from the baseline model in Equation (3.1) as:

$\displaystyle\hat{b}^{(r)}_{t,w}:=\hat{E}\left[D_{t,w}^{(r)}\right]=E_{t,w}^{(r)}\>\hat{\mu}_{t,w}^{(r)},$
(3.3)

for each region $r\in\mathcal{R}$, with $\smash{\hat{\mu}_{t,w}^{(r)}}$
denoting the fitted baseline structure from Equation (3.2). The black solid
line in the left panel of Figure 4 shows an example of this for one selected
region.

###  3.2 Modelling mortality deviations from the baseline model

We aim to explain the deviations from the baseline expected death counts using
a set of features engineered from the available region-specific environmental
factors, as visualized in the right panel of Figure 4. Hereto, we treat the
fitted baseline death counts, $\smash{\hat{b}^{(r)}_{t,w}}$, obtained from
Section 3.1, as fixed and use them as an offset in a predictive modelling
technique. We impose the following distributional assumption:

$\displaystyle
D_{t,w}^{(r)}\sim\text{Poisson}\left(\hat{b}^{(r)}_{t,w}\>\phi_{t,w}^{(r)}\right),$
(3.4)
$\displaystyle\phi_{t,w}^{(r)}=f\Big{(}\texttt{long}^{(r)},\>\texttt{lat}^{(r)},$
$\displaystyle\texttt{season}_{t,w},\>\boldsymbol{c}_{t,w}^{(r)},\>\boldsymbol{e}_{t,w}^{(r)},\>$
$\displaystyle
l^{1}\left(\boldsymbol{c}_{t,w}^{(r)}\right),\>l^{1}\left(\boldsymbol{e}_{t,w}^{(r)}\right),\>\ldots,l^{s}\left(\boldsymbol{c}_{t,w}^{(r)}\right),\>l^{s}\left(\boldsymbol{e}_{t,w}^{(r)}\right)\Big{)}.$

Here, $f(\cdot)$ denotes the outcome of a selected predictive model, such as a
regression tree, random forest, gradient boosting, or neural network. Our
preference for a machine learning model stems from its ability to identify
non-linear relationships of and potential interaction effects among the
environmental features. $\smash{\boldsymbol{c}_{t,w}^{(r)}}$ and
$\smash{\boldsymbol{e}_{t,w}^{(r)}}$ refer to the vectors of, respectively,
the environmental anomalies and extreme indices constructed for week
$w\in\mathcal{W}_{t}$ of year $t\in\mathcal{T}$ in region $r\in\mathcal{R}$.
These anomalies quantify deviations from region-specific baseline conditions,
on which we further elaborate in Section 5.1. To enhance our analysis, we
explicitly incorporate the longitude, i.e., $\smash{\texttt{long}^{(r)}}$, and
latitude coordinate, i.e., $\smash{\texttt{lat}^{(r)}}$, of the centroid of
each region $r$ as well as the season,
$\texttt{season}_{t,w}\in\\{\texttt{Spring},\texttt{Summer},\texttt{Autumn},\texttt{Winter}\\}$,
into the machine learning model.111111We define the seasons as follows: spring
covers March 15 to June 15; summer includes June 15 to September 15; autumn
spans September 15 to December 15; and the winter encompasses December 15 to
March 15. As our analysis focuses on weekly death counts, we define the season
as the one in which the first day of the considered week falls. This inclusion
allows the learning method to take potential interaction effects into account
between the region and the environmental features. For instance, the impact of
temperature on mortality deviations from the baseline may differ between
northern and southern regions. A similar rationale applies to seasonal
effects. Lastly, we allow for lagged effects by including lagged values of the
environmental features in the predictive model $f(\cdot)$. The functions
$l^{u}(\cdot)$, for $u=1,2,..,s$, shift the features $u$ weeks back in time.
Lagged features play a crucial role in capturing various temporal phenomena
and dependencies in data. For instance, they enable modelling the impacts on
mortality of consecutive hot or cold weeks, consecutive elevated air pollution
levels, and potential harvesting effects, among others (Schwartz, 2001).
Harvesting effects occur when, for example, temperature related excess
mortality in the previous week leads to a mortality deficit in the current
week.

We estimate expected death counts from the calibrated predictive model in
Equation (3.4) as:

$\displaystyle\hat{d}^{(r)}_{t,w}:=\hat{b}_{t,w}^{(r)}\>\hat{\phi}_{t,w}^{(r)},$
(3.5)

for each region $r\in\mathcal{R}$, with $\smash{\hat{\phi}_{t,w}^{(r)}}$ the
fitted version of $\smash{\phi_{t,w}^{(r)}}$ in Equation (3.4).

##  4 Calibration strategy and interpretation tools

###  4.1 Calibrating the baseline model

We calibrate the baseline weekly mortality model that captures the region-
specific, seasonal pattern in the observed weekly death counts. Under the
Poisson assumption from Equation (3.1) (Nelder Wedderburn, 1972), we consider
all regions and add a penalty term to obtain smooth variations in the
estimated parameters $\smash{\hat{\beta}_{p}^{(r)}}$ across neighbouring
regions, for $p=0,1,..,5$. Hereto, for notational convenience, we first
denote:

$\displaystyle\boldsymbol{\beta}^{(r)}$
$\displaystyle:=\left(\beta_{0}^{(r)},\beta_{1}^{(r)},\beta_{2}^{(r)},\beta_{3}^{(r)},\beta_{4}^{(r)},\beta_{5}^{(r)}\right)\in\mathbb{R}^{6\times
1}$ $\displaystyle\boldsymbol{z}_{t,w}$
$\displaystyle:=\left(1,t,\sin\left(\frac{2\pi
w}{52.18}\right),\cos\left(\frac{2\pi w}{52.18}\right),\sin\left(\frac{2\pi
w}{26.09}\right),\cos\left(\frac{2\pi
w}{26.09}\right)\right)\in\mathbb{R}^{6\times 1},$

where $\boldsymbol{\beta}^{(r)}$ represents the region-specific parameter
vector, and $\boldsymbol{z}_{t,w}$ the covariate vector for week $w$ in year
$t$ in the Poisson GLM, as detailed in Equations (3.1) and (3.2). Moreover,
for $p=0,1,..,5$, we write $\boldsymbol{\beta}_{p}$ to denote the parameter
vector across the considered regions, i.e.,
$\smash{\boldsymbol{\beta}_{p}:=(\beta_{p}^{(r)})_{r\in\mathcal{R}}}$. We then
calibrate the baseline weekly mortality model by minimising the following
objective function:

$\displaystyle\hat{\boldsymbol{\beta}}=\operatorname*{argmin}_{\boldsymbol{\beta}}\>-l_{P}(\boldsymbol{\beta})+\displaystyle\sum_{p=0}^{5}\lambda_{p}\boldsymbol{\beta}_{p}^{T}\boldsymbol{S}\boldsymbol{\beta}_{p},$
(4.1)

where $\boldsymbol{\beta}:=(\boldsymbol{\beta}^{(r)})_{r\in\mathcal{R}}$ and
$l_{P}(\boldsymbol{\beta})$ is the Poisson log-likelihood, defined, up to a
constant, as:

$\displaystyle
l_{P}(\boldsymbol{\beta})=\displaystyle\sum_{r\in\mathcal{R}}\sum_{t\in\mathcal{T}}\sum_{w\in\mathcal{W}_{t}}\left(d_{t,w}^{(r)}\cdot\left(\log
E_{t,w}^{(r)}+\left(\boldsymbol{\beta}^{(r)}\right)^{T}\boldsymbol{z}_{t,w}\right)-E_{t,w}^{(r)}\cdot\exp\left(\left(\boldsymbol{\beta}^{(r)}\right)^{T}\boldsymbol{z}_{t,w}\right)\right)$
(4.2)

with T the transpose operator, and where we assume that the
$\smash{d_{t,w}^{(r)}}$ are independent realizations from the Poisson
distribution in Equation (3.1). Additionally, the objective function to
minimize in Equation (4.1) contains penalty terms
$\smash{\lambda_{p}\boldsymbol{\beta}_{p}^{T}\boldsymbol{S}\boldsymbol{\beta}_{p}}$,
for $p=0,1,..,5$, to impose smooth variations in the parameter vector
$\smash{\boldsymbol{\beta}_{p}}$ across neighbouring regions. Hereto, we use
the following entries in the penalty matrix $\smash{S:=(s_{ij})_{i,j}}$, with
$i,j\in\mathcal{R}$:

$s_{ij}=\begin{cases}|\mathcal{N}_{i}|&\>\>\>\text{if}\>\>i=j\\\
-1&\>\>\>\text{if $i\neq j$ are neighbouring regions}\\\
0&\>\>\>\text{elsewhere},\end{cases}$ (4.3)

with $\mathcal{N}_{i}$ the set of neighbours of region $i\in\mathcal{R}$. The
smoothness penalty in Equation (4.3) leads to (see Appendix B):

$\displaystyle\boldsymbol{\beta}_{p}^{T}\boldsymbol{S}\boldsymbol{\beta}_{p}=$
$\displaystyle\sum_{i\in\mathcal{N}_{1}}\left(\beta_{p,1}-\beta_{p,i}\right)^{2}\>+\>\sum_{\mathclap{i\in\mathcal{N}_{2}\backslash\\{1\\}}}\left(\beta_{p,2}-\beta_{p,i}\right)^{2}\>+\>\ldots\>+\>\sum_{\mathclap{i\in\mathcal{N}_{R-1}\backslash\\{1,..,R-2\\}}}\left(\beta_{p,R-1}-\beta_{p,i}\right)^{2},$

where, e.g., $\mathcal{N}_{k}\backslash\\{1,..,k-1\\}$ represents the set of
neighbours of region $k$, excluding the first $k-1$ regions, for
$k=1,2,...,R-1$. We hence observe that the imposed penalty matrix
$\boldsymbol{S}$ in Equation (4.3) penalizes the sum of the squared
differences between the parameters of neighbouring regions. Appendix B shows
the smoothness penalty matrix as constructed for a set of five Spanish NUTS 3
regions. The parameter $\lambda_{p}$ in Equation (4.1) is a smoothing
parameter controlling the degree of smoothness in the parameter
$\boldsymbol{\beta}_{p}$ across neighbouring regions. A very large
$\smash{\lambda_{p}}$ results in the same $\smash{\beta_{p}^{(r)}}$ across all
regions, while a $\lambda_{p}$ close to $0$ results in no penalization on the
parameter $\smash{\beta_{p}^{(r)}}$, which corresponds to the parameters of a
traditional, unregularized Poisson GLM.

We estimate the parameters $\smash{\beta_{p}^{(r)}}$ in Equation (4.1) using
the penalized iteratively re-weighted least squares algorithm, where we select
optimal smoothing parameters $\lambda_{p}$ using Un-Biased Risk Estimation
(UBRE) scores (Wood, 2017).121212We use the argument paraPen from the gam
function in the R-package mgcv to implement the proposed spatially smoothed
Poisson GLM, see Wood (2015). We denote the optimal smoothing parameters as
$\hat{\lambda}_{p}$, for $p=0,1,..,5$, and the optimal parameter vector as
$\smash{\hat{\boldsymbol{\beta}}}$. We then estimate the expected baseline
death counts for each region $r\in\mathcal{R}$, year $t\in\mathcal{T}$, and
week $w\in\mathcal{W}_{t}$ as:

$\displaystyle\smash{\hat{b}^{(r)}_{t,w}}=E_{t,w}^{(r)}\cdot\exp\left(\left(\hat{\boldsymbol{\beta}}^{(r)}\right)^{T}\boldsymbol{z}_{t,w}\right).$

###  4.2 Calibrating the mortality deviations model

Fixing the expected weekly baseline death counts, as obtained in Section 4.1,
we now calibrate a predictive model to explain deviations from this baseline
using environmental features. As discussed in Section 3.2, we employ a machine
learning approach due to its capacity to handle high dimensional features, to
capture non-linear relationships and to detect interactions among the
environmental features. Machine learning algorithms typically rely on a set of
parameters of which some are carefully selected through a tuning process
(tuning parameters), while others are set to predetermined values (hyper-
parameters). In this paper, we opt for an extensive grid search that involves
exploring a predefined grid of parameter values to identify optimal tuning
parameter configurations.

##### XGBoost algorithm.

While any predictive modelling technique can be deployed in our framework, we
will specifically focus on the extreme gradient boosting machine (Chen
Guestrin, 2016) in the case study outlined in Section 5. Appendix D.1 explains
the XGBoost algorithm in case of a Poisson distributed outcome, see Equation
(3.4), in Appendix D.1 and details the considered tuning parameters. We hereby
choose the negative Poisson log-likelihood as the loss function:

$\displaystyle\mathcal{L}\left(d_{t,w}^{(r)},f\left(\boldsymbol{x}_{t,w}^{(r)}\right)\right)=\hat{b}_{t,w}^{(r)}\cdot
f\left(\boldsymbol{x}_{t,w}^{(r)}\right)-d_{t,w}^{(r)}\cdot\left(\log
f\left(\boldsymbol{x}_{t,w}^{(r)}\right)+\log\hat{b}_{t,w}^{(r)}\right)+\log
d_{t,w}^{(r)}!,$ (4.4)

for $r\in\mathcal{R}$, $t\in\mathcal{T}$ and $w\in\mathcal{W}_{t}$. Here,
$\smash{\boldsymbol{x}_{t,w}^{(r)}}$ is the input vector of the XGBoost model
consisting of the (lagged) weather anomalies and extreme weather indices
$\smash{\boldsymbol{c}_{t,w}^{(r)}}$, the (lagged) air pollution anomalies and
extreme air pollution indices $\smash{\boldsymbol{e}_{t,w}^{(r)}}$, the season
$\texttt{season}_{t,w}$ and the longitude $\smash{\texttt{long}^{(r)}}$ and
latitude $\smash{\texttt{lat}^{(r)}}$ coordinates, see Section 3.2.

##### Parameter tuning with cross-validation.

We use $K$-fold cross-validation to select optimal values for the six tuning
parameters listed in Table D.1, inspired by Hastie . (2001). For simplicity,
let $\mathcal{T}=\\{0,1,2,...,T-1\\}$ represent the considered time range
consisting of $T$ years. We choose the number of folds ($K$) to be equal to
the number of years in our training data, i.e., $K=T$. Each fold corresponds
to a specific year $t\in\mathcal{T}$ as validation data and the remaining
years as training data. Consequently, for each fold, the validation data
consists of either $53\cdot|\mathcal{R}|$ observations for a leap year or
$52\cdot|\mathcal{R}|$ observations otherwise. Appendix D.2 describes the
tuning strategy in further detail, with a visualisation in Figure D.1.

##### Training

We obtain a carefully tuned XGBoost model, calibrated on the available
historical data covering the entire duration of $T$ years, utilising the
optimally chosen values for the tuning parameters. We denote the outcomes
estimated from the XGBoost model as
$\smash{\hat{f}_{\text{XGBoost}}\left(\boldsymbol{x}_{t,w}^{(r)}\right)}$ for
$r\in\mathcal{R}$, $t\in\mathcal{T}$ and $w\in\mathcal{W}_{t}$. Consequently,
the estimated death counts equal:

$\displaystyle\hat{d}_{t,w}^{(r)}=\hat{b}_{t,w}^{(r)}\cdot\hat{f}_{\text{XGBoost}}\left(\boldsymbol{x}_{t,w}^{(r)}\right),\hskip
14.22636ptr\in\mathcal{R},\>t\in\mathcal{T},\>w\in\mathcal{W}_{t}.$

As such, we can interpret the outcomes from the XGBoost model as a multiplier
that can either augment or diminish the baseline number of death counts.

###  4.3 Interpretation tools

Machine learning algorithms, like the XGBoost model, are often seen as black
boxes. Therefore, the use of interpretation tools becomes crucial to unravel
insights from these black box algorithms (Molnar, 2019). We further assume
that the training data, i.e.,

$\displaystyle\left\\{\left(\boldsymbol{x}_{t,w}^{(r)},\hat{b}_{t,w}^{(r)},d_{t,w}^{(r)}\right)\right\\}_{r\in\mathcal{R},t\in\mathcal{T},w\in\mathcal{W}_{t}},$
(4.5)

consists of feature vectors containing $q$ variables:

$\displaystyle\boldsymbol{x}_{t,w}^{(r)}=\left(x_{1,t,w}^{(r)},x_{2,t,w}^{(r)},\ldots,x_{q,t,w}^{(r)}\right)\in\mathbb{R}^{q}.$

For simplicity, we refer to $X_{l}$ as the $l$-th feature, for $l=1,..,q$, as
random variable and to $f(\cdot)$ as the fitted XGBoost model
$\hat{f}_{\text{XGBoost}}(\cdot)$.

##### Feature importance

In this interpretation tool, our objective is to unveil the features that
significantly contribute to the predictions. We follow the approach proposed
by Breiman Ihaka (1984) and measure the importance of a particular feature
$X_{l}$ within a regression tree by aggregating the reductions in the
considered loss function across all splits associated to this feature. Given
that the XGBoost model consists of multiple regression trees, see Appendix
D.1, we extend this analysis to measure the reduction in loss caused by
feature $X_{l}$ across all trees in the ensemble. Appendix D.3 outlines the
mathematical formulation of the feature importance. The higher the feature
importance, the more important the feature is to the prediction process.

##### Accumulated Local Effects

Apley Zhu (2020) propose Accumulated Local Effects (ALE) plots to visualize
the impact of a particular feature on the predictions generated by a machine
learning model. In contrast to the partial dependence plots introduced by
Friedman (2001), ALE plots prove to be more effective in illustrating feature
effects when dealing with correlated features (Molnar, 2019). Given the strong
correlation present in the environmental features, this paper places specific
emphasis on the relevance of ALE plots for visualising and interpreting the
impact of a feature on the predictions.

The ALE effect of a feature $X_{l}$ at value $x$ is defined as:

$\displaystyle
f_{l,\text{ALE}}(x)=\displaystyle\int_{z_{0,l}}^{x}\mathbb{E}\left[\frac{\partial
f(X_{1},X_{2},...,X_{p})}{\partial
X_{l}}\Big{|}X_{l}=z_{l}\right]dz_{l}-c_{l}$ (4.6)

where $z_{0,l}$ is a lower bound on the range of $X_{l}$, typically chosen as
the minimum value observed in the data. We interpret the integrand in Equation
(4.6) as the expected change in the prediction function $f$ when slightly
changing the feature $X_{l}$ around the value $z_{l}$, also known as the local
effect. Afterwards we accumulate these local effects by integrating them from
$z_{0,l}$ to $x$ to unravel the global effect of the feature $X_{l}$ on the
response. Taking the conditional expectation in Equation (4.6) isolates the
effect of the feature $X_{l}$ on the prediction function $f$ from the effects
of all other, (possibly) correlated features. Further, Apley Zhu (2020)
introduce the constant $c_{l}$ such that the expectation of the random
variable $f_{l,\text{ALE}}(X_{l})$ equals zero.

To investigate interaction effects between two environmental features, say
$X_{k}$ and $X_{l}$, we calculate the ALE interaction effect (Apley Zhu,
2020). Let $z_{0,k}$ and $z_{0,l}$ be the lower bounds on the ranges of
$X_{k}$ and $X_{l}$ respectively, then:

$\displaystyle
f_{k,l,\text{ALE}}(x,y)=\displaystyle\int_{z_{0,k}}^{x}\int_{z_{0,l}}^{y}\mathbb{E}\left[\frac{\partial^{2}f(X_{1},X_{2},...,X_{p})}{\partial
X_{l}\partial
X_{k}}\Big{|}X_{k}=z_{k},X_{l}=z_{l}\right]dz_{l}dz_{k}-c_{k,l},$ (4.7)

where now the conditional expectation is taken with respect to the joint
distribution of the random variables $X_{j}$ for $j\neq k,l$. Equation (4.7)
considers the second order partial effects and accumulates these effects in
two dimensions. As such, the main effects of the features are not considered
in the calculation and the ALE interaction effect should therefore be
interpreted as an additional interaction effect on top of the main effects.
Further, the constant $c_{k,l}$ is determined such that the expectation of
$f_{k,l,\text{ALE}}(X_{k},X_{l})$ equals zero. Appendix D.3 outlines the
estimation procedure for calculating the ALE (interaction) effect on the
training data.

##  5 Case study on the NUTS 3 regions in 20 European countries

###  5.1 Feature engineering

To examine the association between environmental factors and weekly mortality
rates, we use data on weekly death counts within the NUTS 3 geographical
regions for the years 2013-2019, as detailed in Section 2.2. Hereto, we first
aggregate the temporal and spatial dimensions of the gridded, environmental
data into daily features, registered at NUTS 3 level, see Sections 5.1.1 and
5.1.2 respectively. To capture the effects of extreme environmental
conditions, such as heat waves and cold spells, we engineer extreme
environmental indices that indicate how many days within a week the
environmental factor exceeds a certain high quantile or is lower than a
certain low quantile (see Section 5.1.3). Furthermore, since the focus in this
paper is on explaining deviations from a baseline model for weekly death
counts, we also engineer environmental anomalies that quantify deviations from
normal, baseline conditions for each week in the year (see Section 5.1.4).
Lastly, we aggregate the environmental features on a weekly basis and create
lagged versions of it in Section 5.1.5. Figure 5 illustrates the feature
engineering process for the weather related factors listed in Table 1 and the
air pollution factors in Table 2. The flow chart refers to the section that
further details each step at the upper right corner of the corresponding box.

Air pollution factor S2.2 Daily aggregation S5.1.1 Weather factor S2.2 Spatial
aggregation S5.1.2 Environmental anomalies S5.1.3 Extreme environmental
indices S5.1.4 Weekly aggregation \+ lagged weekly features
S5.1.5HourlyDailyWeekly(long,lat)NUTS 3

Figure 5: Flow chart explaining the feature engineering process for the
weather-related factors and the air pollution factors, with the relevant
section that details each step at the upper right corner of the box.

####  5.1.1 Daily aggregation of hourly air pollutant levels

We consider an air pollution factor from Table 2 and denote it as
$x_{t,w,d,h}^{(\text{long,lat})}$. This feature represents the concentration
of a specific air pollutant in year $t$, week $w$, day $d$, and hour $h$,
located at longitude-latitude coordinates (long,lat). We then compute the
daily minimum, average, and maximum concentrations of the air pollutant,
measured at the coordinates (long,lat) as follows:

$\displaystyle\overset{\wedge}{x}_{t,w,d}^{(\text{long,lat})}$
$\displaystyle=\min\left\\{x_{t,w,d,h}^{(\text{long,lat})}\>\big{|}\>h=0,1,...,23\right\\}$
$\displaystyle\overline{x}_{t,w,d}^{(\text{long,lat})}$
$\displaystyle=\text{avg}\left\\{x_{t,w,d,h}^{(\text{long,lat})}\>\big{|}\>h=0,1,...,23\right\\}$
$\displaystyle\overset{\vee}{x}_{t,w,d}^{(\text{long,lat})}$
$\displaystyle=\max\left\\{x_{t,w,d,h}^{(\text{long,lat})}\>\big{|}\>h=0,1,...,23\right\\}.$

For the weather factors listed in Table 1, we do not need to make the
conversion from an hourly to a daily time scale, since these factors are
already defined on a daily time scale at specific longitude-latitude
coordinates.

####  5.1.2 Population-weighted spatial aggregation to NUTS 3 level

Let $\tilde{x}_{t,w,d}^{(\text{long,lat})}$ denote the measurement of a
particular environmental factor at coordinates (long,lat), year $t$, week $w$,
and day $d$. This feature can encompass minimum, maximum, or average air
pollution measurements, as defined in Section 5.1.1, or any weather-related
factor from Table 1.

We consider the grid $\mathcal{G}^{(1)}$ of longitude-latitude coordinates
(long,lat) on which the feature $\tilde{x}_{t,w,d}^{(\text{long,lat})}$ is
registered. We then introduce the function:

$\displaystyle
m_{\mathcal{G}}:\mathcal{G}\rightarrow\mathcal{R}:\>\text{(long,lat)}\mapsto
r,$

which maps longitude-latitude coordinates from a grid $\mathcal{G}$ to the
NUTS 3 region $r$ it belongs. Now, define the set $\mathcal{I}_{1}(r)$ that
consists of all (long,lat) coordinates of the grid $\mathcal{G}^{(1)}$
contained in the NUTS 3 region $r$:

$\displaystyle\mathcal{I}_{1}(r)=\left\\{\text{(long,lat)}\in\mathcal{G}^{(1)}\mid
m_{\mathcal{G}^{(1)}}(\text{long},\text{lat})=r\right\\}.$ (5.1)

We create a weighted aggregated feature at the NUTS 3 geographical level as
follows:

$\displaystyle\tilde{x}_{t,w,d}^{(r)}=\displaystyle\smashoperator[r]{\sum_{\begin{subarray}{c}\text{(long,lat)}\in\mathcal{I}_{1}(r)\end{subarray}}^{}}\omega_{\text{(long,lat)}}\cdot\tilde{x}_{t,w,d}^{(\text{long},\text{lat})},$
(5.2)

where $\omega_{\text{(long,lat)}}$ are weights that sum up to one within each
grid $\mathcal{I}_{1}(r)$. These weights aim to compute a population-weighted
version of the daily environmental factors, a method that has also been
explored in the existing literature (Balakrishnan ., 2019; de Schrijver .,
2021). This weighting strategy is necessary due to the uneven distribution of
a population within a NUTS 3 region, which tends to concentrate around
specific locations. Therefore, we put more weight on longitude-latitude
coordinates where more people are located, relative to the population in that
NUTS 3 region. We further assume the weights $\omega_{\text{(long,lat)}}$ to
remain constant over time. This assumption is consistent with the idea that
the proportion of the population within region $r$ associated with grid point
(long,lat) $\in\mathcal{I}_{1}(r)$ remains relatively stable throughout the
time span under consideration.131313Our current method assigns each grid point
(long,lat) $\in\mathcal{G}^{(1)}$ exclusively to one NUTS 3 region $r$.
Alternatively, an interpolation method could be explored, enabling grid points
near the borders of NUTS 3 regions to also contribute to the calculation of
the feature at NUTS 3 level in Equation (5.2).

We compute the weights using the Socioeconomic Data and Applications Center
(SDAC) in NASA’s Earth Observing System Data and Information System
(EOSDIS).141414This data center is available on
https://sedac.ciesin.columbia.edu/. Particularly, we use the high-dimensional
population count dataset, defined on a grid $\mathcal{G}^{(2)}$ with a spatial
resolution of 2.5 arc-minutes ($\approx 5\text{ km}$) and a temporal dimension
of five years (CIESIN et. al, 2018). We extract the population gridded dataset
at the year $2015$ to construct the population weights in Equation
(5.2).151515We select the year $2015$ since it falls within our considered
time range. Appendix C.1 details the strategy on how to compute the population
weights, with a visualisation in Figure C.1.

####  5.1.3 Daily extreme environmental indices

As discussed in Section 3.1, we develop a weekly baseline mortality model and
aim to explain deviations from this baseline by using region-specific
environmental features. As motivated in Section 3, we hereto construct
environmental anomalies and extreme environmental indices, along with lagged
versions, to explain excess or deficit mortality relative to the baseline.

To capture the effects of extreme environmental conditions on mortality, we
construct daily, region-specific features that indicate whether the daily,
environmental measurement in that region exceeds a certain high quantile or
falls below a certain low quantile. Hereto, we first calculate the region-
specific $5\%$ and $95\%$ quantile of the entire set of historical
observations on the daily minimum temperature Tmin, the average temperature
Tavg, and the maximum temperature Tmax throughout the considered time period
2013-2019. We define the extreme high temperature index for region $r$ on day
$d$ in week $w$ of year $t$ as:

$\displaystyle\small\text{T.ind}_{t,w,d}^{(r,95\%)}={{\mathbbm{1}}}\left\\{\texttt{Tmax}_{t,w,d}^{(r)}\geq
q_{\texttt{Tmax}}^{(r,95\%)}\right\\}+{{\mathbbm{1}}}\left\\{\texttt{Tavg}_{t,w,d}^{(r)}\geq
q_{\texttt{Tavg}}^{(r,95\%)}\right\\}+{{\mathbbm{1}}}\left\\{\texttt{Tmin}_{t,w,d}^{(r)}\geq
q_{\texttt{Tmin}}^{(r,95\%)}\right\\},$ (5.3)

with, e.g., $\mathbbm{1}\\{\texttt{Tmax}_{t,w,d}^{(r)}\geq
q_{\texttt{Tmax}}^{(r,95\%)}\\}$ an indicator function which equals one when
the maximum temperature in region $r$ on day $d$ in week $w$ of year $t$
exceeds its 95$\%$ quantile, and zero otherwise. We refer to this feature as
the hot-day index. The index can take values 0, 1, 2, or 3, where a value of
three corresponds to a very hot day. We define a similar cold-day index using
the $5\%$ quantile, where a value of three corresponds to a very cold day.
Similarly, we create daily indices for surpassing high thresholds of the other
environmental factors and refer to them as extreme environmental indices.

####  5.1.4 Baseline models for daily environmental features

Alongside the extreme environmental indices, constructed in Section 5.1.3, we
engineer environmental anomalies that measure deviations from daily baseline
conditions. Some of the daily environmental features exhibit a clear seasonal
pattern. For these features, we first fit a region-specific baseline model
using robust linear regression with a singular set of sine and cosine Fourier
terms as covariates for each region separately (Gervini Yohai, 2002).161616By
adopting robust linear regression, the fitted baseline model becomes less
susceptible to the impact of extreme levels within the observed environmental
features. We use the lmRob function in the R-package robust. In this way, we
capture the region-specific baseline environmental conditions on each day
within a year. We have:

$\displaystyle\tilde{x}_{t,w,d}^{(r)}=\alpha_{0}^{(r)}+\alpha_{1}^{(r)}\sin\left(\frac{2\pi
w}{365.25}\right)+\alpha_{2}^{(r)}\cos\left(\frac{2\pi
w}{365.25}\right)+\epsilon_{t,w,d}^{(r)},$ (5.4)

where $\smash{\epsilon_{t,w,d}^{(r)}}$ is a normally distributed error term
and $365.25$ represents the average number of days per
year.171717Alternatively, the exact number of days in a year can be used, with
366 days in a leap year and 365 days otherwise. In case a feature does not
exhibit a clear annual, seasonal pattern, we do not consider a baseline model
for that feature.

Figure 6 shows an example of the daily baseline trend for three weather
factors as registered in three different NUTS 3 regions. Figure 6 visualizes
the seasonal baseline trend for the daily maximum temperature levels in
Barcelona (ES511) in red. Figure 6 shows the seasonal baseline trend for the
daily wind speed levels in Milano (ITC4C) in green. These estimated baseline
trends align with our expectations, showing higher maximum temperatures and
lower wind speeds during summer days compared to winter days. Since the total
daily rainfall does not exhibit a clear seasonal pattern, we consider a zero
baseline trend for the variable Rain, visualized by the blue line in Figure 6.
In the remainder of this section, we will work with the excesses or deviations
from the baseline levels, also referred to as anomalies. These anomalies are
the residuals of the region-specific baseline models for the different
environmental features in Equation (5.4).

()

()

()

Figure 6: We visualize the baseline model across the years 2013-2019, fitted
by means of robust linear regression with one sine-cosine Fourier pair for the
daily maximum temperatures observed in Barcelona (a) and for the daily wind
speed levels in Milano (b). In panel (c), we show the daily precipitation
amounts in Stockholm with a zero baseline trend in blue.

####  5.1.5 Weekly averages of daily environmental anomalies and extreme
indices

We align the daily time scale of the environmental anomalies (see Section
5.1.4) and the extreme environmental indices (see Section 5.1.3), with the
weekly time scale of the death counts. Hereto, let
$\smash{\tilde{a}_{t,w,d}^{(r)}}$ be any environmental anomaly or extreme
environmental index, registered in region $r$ at day $d$ of week $w$ in year
$t$. We then compute the weekly averages of this environmental feature as:

$\displaystyle\overline{a}_{t,w}^{(r)}=\text{avg}\left\\{\tilde{a}_{t,w,d}^{(r)}\mid
d=1,2,...,7\ \right\\},$

for all $t\in\mathcal{T}$, $w\in\mathcal{W}_{t}$, and $r\in\mathcal{R}$. This
weekly averaging process captures effects such as the heightened impact of
prolonged periods of environmental stress on mortality rates, e.g., multiple
hot days within a week. Table 3 lists the final features which we use as input
features in the machine learning model. Alongside these features, we also
introduce the lagged feature values, using the suffix _l1, e.g.,
w_avg_Tmax_anom_l1.

Feature 5.5cmWeekly average of the daily w_avg_Tmax_anom maximum temperature
anomalies. w_avg_Tmin_anom minimum temperature anomalies. w_avg_T.ind95 hot-
day index, see Equation (5.3). w_avg_T.ind5 cold-day index. w_avg_Hum_anom
average relative humidity anomalies. w_avg_Hum.ind95 high humidity index.
w_avg_Hum.ind5 low humidity index. w_avg_Rain total precipitation levels.
w_avg_Rain.ind95 high precipitation index. w_avg_Rain.ind5 low precipitation
index. w_avg_Wind_anom average wind speed anomalies. w_avg_Wind.ind95 high
wind speed index. w_avg_Wind.ind5 low wind speed index.

Feature 5.5cmWeekly average of the daily w_avg_O3_anom average ozone
anomalies. w_avg_O3.ind95 high ozone index. w_avg_O3.ind5 low ozone index.
w_avg_NO2_anom average nitrogen dioxide anomalies. w_avg_NO2.ind95 high
nitrogen dioxide index. w_avg_NO2.ind5 low nitrogen dioxide index.
w_avg_PM10_anom average PM10 anomalies. w_avg_PM10.ind95 high PM10 index.
w_avg_PM10.ind5 low PM10 index. w_avg_PM2.5_anom average PM2.5 anomalies.
w_avg_PM2.5.ind95 high PM2.5 index. w_avg_PM2.5.ind5 low PM2.5 index.

Table 3: Weekly region-specific, environmental features available at the NUTS
3 geographical level.

Figure 7 illustrates the correlations among the weekly environmental features,
as detailed in Table 3. The most notable positive correlations appear between
the weekly averages of daily PM2.5 and PM10 anomalies (0.95), and between the
weekly averages of daily maximum and minimum temperature anomalies (0.73).
Additionally, unsurprisingly, strong positive correlations are observed
between the weekly averages of various environmental anomalies and their
respective high quantile indices. Conversely, notable negative correlations
emerge between the weekly averages of daily environmental anomalies and their
corresponding low quantile indices.

Figure 7: Correlations between the weekly environmental features, constructed
through the feature engineering process in Figure 5, and presented in Table 3.

Figure 8 presents the weekly aggregated features at the NUTS 3 level during
the 30th ISO week of 2018. Figure 8 illustrates the weekly average of the
daily maximum temperature anomalies (w_avg_Tmax_anom), while Figure 8 depicts
the weekly average of the daily hot day index (w_avg_Tind95). In Figure 8, we
show the weekly average of the daily average relative humidity anomalies
(w_avg_Hum_anom), and Figure 8 displays the weekly average of the daily
average ozone anomalies (w_avg_O3_anom). We conclude from Figure 8 that the
30th ISO week of 2018 exhibited significantly elevated temperatures in the
Nordic countries and the northern regions of France, surpassing region-
specific baseline temperature levels by more than six degrees. A similar
analysis for the weekly hot-day index reveals heightened values in the
northern regions, indicating not only high daytime temperatures but also
warmer-than-normal nights (see Figure 8). In terms of humidity, Figure 8 shows
lower-than-baseline levels in Central European regions, while southern regions
exhibited higher humidity levels. This is consistent with the typical negative
correlation between temperature and humidity (Fischer Knutti, 2013).
Furthermore, the northern French and Belgian regions exhibited higher ozone
levels compared to baseline concentrations, as Figure 8 suggest.

()

()

()

()

Figure 8: Weekly average of the daily maximum temperature anomalies (a), daily
hot day index (b), daily average relative humidity anomalies (c), and daily
average ozone anomalies (d) in the 30th ISO week of the year 2018 in the NUTS
3 European regions.

###  5.2 Calibration

####  5.2.1 Calibrating the baseline model for weekly death counts

We calibrate the baseline model for weekly death counts, as outlined in
Equations (3.1) and (3.2), by minimising the penalized Poisson negative log-
likelihood, as discussed in Section 4.1. This penalization ensures a smooth
variation of the fitted parameters $\smash{\beta_{p}^{(r)}}$ in the baseline
model across neighbouring regions, which is particularly useful for regions
with a low number of weekly death counts. We consider the calibration period
$\mathcal{T}=\\{2013,2014,...,2019\\}$.

Figure 9 illustrates the fitted parameters in the spatially smoothed baseline
model, showing smooth variations in the parameters between adjacent regions.
The estimated parameter $\smash{\hat{\beta}_{1}^{(r)}}$ for the ISO year
effect exhibits a negative sign in northern European regions and a slightly
positive sign in southern regions. This suggests that, from 2013 to 2019,
mortality rates slightly increased in the south while declining in the north.
Moreover, the larger magnitude of the estimated parameters related to the sine
and cosine Fourier terms in southern regions suggests that the associated
seasonal patterns (with 52-week or 26-week periods) have a more pronounced
effect on weekly death counts in southern regions compared to northern
regions.

Figure 9: The fitted parameters $\hat{\beta}_{0}^{(r)}$,
$\hat{\beta}_{1}^{(r)}$, and $\hat{\beta}_{2}^{(r)}$ (top panels) and
$\hat{\beta}_{3}^{(r)}$, $\hat{\beta}_{4}^{(r)}$, and $\hat{\beta}_{5}^{(r)}$
(bottom panels) in the spatially smoothed Poisson GLM belonging to the
intercept, ISO year $t$ and the sine and cosine Fourier terms with a period of
52 or 26 weeks, see Equation (3.2).

Figure 10 illustrates the estimated baseline mortality rates
$\hat{q}_{t,w}^{(r)}$ in red, as obtained from Equation (2.1) when using the
estimated force of mortality specified in Equation (3.2). In line with the
findings from Figure 9, we observe, for example, that the seasonal variation
in Barcelona is more pronounced than in Milano and Stockholm. Furthermore, the
mortality rates in the southern regions show a relatively stable (Barcelona)
or slightly increasing (Milano) trend, whereas the mortality rates in the
region of Stockholm clearly exhibit a decreasing trend.

####  5.2.2 Tuning and calibrating the machine learning model

After the calibration of the baseline model for the weekly death counts in
Section 5.2.1, we calibrate the XGBoost model as explained in Section 4.2.
This machine learning model incorporates the environmental features engineered
in Section 5.1 and listed in Table 3. Additionally, we incorporate the one-
week lagged values of the features in Table 3 and the season. We furthermore
allow for spatial variations in the impact of specific features on the weekly
deviations from the baseline death counts by incorporating the longitude
(long) and latitude (lat) coordinates of the centre of each NUTS 3 region.
This leads to a total of 56 features used as inputs in the machine learning
model.

We consider the following tuning grid for the six XGBoost parameters of
interest:

$\displaystyle\texttt{nrounds}\in\\{10,20,...,5\ 000\\},\hskip
8.5359pt\texttt{eta}\in\\{0.01,0.05,0.1\\},\hskip
8.5359pt\texttt{min\\_child\\_weight}\in\\{10,100,1000\\}$
$\displaystyle\texttt{subsample}\in\\{0.50,0.75\\},\hskip
8.5359pt\texttt{colsample\\_bytree}\in\\{0.50,0.75\\},\hskip
8.5359pt\texttt{max.depth}\in\\{1,3,5,7,9\\}.$

We tune these six parameters by means of 7-fold cross-validation, following
the approach outlined in Section 4.2.181818We consider a 7-fold cross-
validation as the calibration period $\mathcal{T}$ consists of seven years,
i.e., the years 2013-2019. Throughout this cross-validation process, the
offset, representing the estimated baseline deaths, is considered to be known.
While, in principle, we should recalibrate this offset for each set of
training fold combinations, we refrain from doing so due to computational
constraints. Moreover, the robust construction of the baseline model ensures
that excluding one year from the training set will only have a negligible
impact on the baseline fit and consequently on the offset.

The optimal set of tuning parameters consists of $490$ boosting iterations
(nrounds), a learning rate of $0.01$ (eta), a minimum child weight of $1\ 000$
(min_child_weight), a maximum depth of $7$ for each tree (max.depth), and a
subsample ratio of $75\%$ and $50\%$ for the training data (subsample) and the
columns (colsample_bytree), respectively. Subsequently, we retrain the machine
learning model using this optimal parameter set over the entire calibration
period spanning the years 2013-2019.

###  5.3 Insights in the machine learning model

####  5.3.1 In-sample fit and model performance

Figure 10 displays the observed mortality rates in grey, alongside the
estimated mortality rates from both the baseline model (red) and the machine
learning model integrating environmental features (blue). The comparison
reveals that the machine learning model performs better in terms of an in-
sample fit compared to the baseline model. In Barcelona and Milano, the
proposed model effectively captures the excess mortality during the winters of
2015/2016 and 2017/2018, and the excess mortality observed during the summer
of 2015. In Stockholm, the model identifies the excess mortality observed
during the summer of 2014 and 2018. These findings are confirmed by Figure 10,
which visualizes the residuals of the estimated mortality rates, calculated
as:

$\displaystyle\text{resid}_{t,w}^{(r)}=\hat{q}_{t,w}^{(r)}-q_{t,w}^{(r)},$

where $\smash{\hat{q}_{t,w}^{(r)}}$ and $\smash{q_{t,w}^{(r)}}$ represent the
estimated and observed mortality rate in week $w$ of year $t$ in region $r$,
respectively. The estimated mortality rate is derived from either the baseline
or the machine learning model (XGBoost). We clearly observe less peaks in the
residuals of the machine learning model compared to the residuals of the
baseline model, particularly during the summer and winter weeks.

()

()

Figure 10: Panel (a): the observed and estimated weekly mortality rates during
the period 2013-2019 in the NUTS 3 regions of Barcelona (left), Milano
(middle), and Stockholm (right). The observed weekly mortality rates are shown
in grey and the mortality rates as obtained from the baseline model and the
XGBoost model in red and blue, respectively. Panel (b): residuals of the
estimated weekly mortality rates derived from the baseline model (red) and the
XGBoost model (blue) in Barcelona (left), Milano (middle), and Stockholm
(right).

####  5.3.2 Interpretation tools

##### Feature importance.

We analyse which features contribute most to the machine learning model that
captures deviations in mortality from the baseline model. Hereto, we use the
feature importance measure outlined in Section 4.3. Figure 11 presents the
features for which their feature importance exceeds $1\%$. Moreover, we show
the 95$\%$ non-parametric bootstrap confidence intervals of the feature
importance measures, derived from 1 000 bootstrap samples of the same size as
the training data (DiCiccio Efron, 1996).191919To create a non-parametric
bootstrap sample from the training data, we randomly sample observations with
replacement. For computational efficiency, we fit the machine learning model
on each bootstrap sample using the tuning parameter configuration obtained in
Section 5.2.2. The $x$-axis represents the feature importance in percentages,
indicating the relative contribution of each feature to the predictive model.
The six features with the highest feature importance are all associated with
temperature. Additionally, based on the bootstrap confidence intervals, these
features consistently exhibit the highest feature importance across the
bootstrapped datasets. Among the top three features, we find
w_avg_Tmin_anom_l1, w_avg_T.ind95, and w_avg_T.ind5_l1 which denote the lagged
weekly average of the daily minimum temperature anomalies, the weekly average
of the daily hot day index, and the lagged weekly average of the daily cold-
day index, respectively. In contrast, features related to weather factors such
as Wind, Hum, and Rain demonstrate limited importance in the model, as they do
not rank among the top 10 predictors. This is in line with earlier studies
such as Braga . (2002) and Alberdi . (1998).

Figure 11: The features in the XGBoost model with a feature importance
exceeding $1\%$. The feature importance of the features related to weather are
visualized in red, air pollution in black, the NUTS 3 region in green, and the
season in blue. Furthermore, the black error bars denote the 95$\%$ non-
parametric bootstrap confidence interval for each feature based on $1\ 000$
bootstrap samples of the training data.

##### ALE main effects.

Figure 12 shows the ALE main effects for the features with an importance
surpassing $1\%$. Features with lower importance in the machine learning model
yield less reliable ALE main effects and we, consequently, exclude these from
the illustration. Additionally, we show the $95\%$ pointwise bootstrap
confidence interval of the ALE main effect of each feature based on the fit
constructed on each of the 1 000 bootstrap copies of the training data.

()

()

()

()

()

Figure 12: The ALE main effects for features whose feature importance exceeds
1$\%$, denoted within parentheses in each panel’s legend. To enhance
visualization, the prefix w_avg_ was removed from the feature names in the
legend. We furthermore present the 95$\%$ pointwise confidence interval based
on 1000 bootstrap samples of the training data. The different panels display
the ALE main effects of the extreme temperature indices (a), temperature
anomalies (b), weather anomalies including humidity, wind speed, and rainfall
(c), extreme air pollution indices (d), and air pollution anomalies (e).

In Figure 12, we observe an increasing trend in the ALE main effects
corresponding to the four extreme temperature indices. We estimate that higher
weekly averages of the daily hot or cold day indices, potentially lagged by
one week, correspond to substantially higher excess deaths relative to the
baseline. In accordance with the bootstrap-generated confidence intervals,
these features have a statistically significant impact on the model’s
predictions. In Figure 12, we show the ALE main effect for the four
temperature anomalies. The effect is the strongest and most significant for
the feature w_avg_Tmin_anom_l1. When the lagged weekly average of the daily
minimum temperature anomalies is higher, the ALE main effect is lower. This
effect might vary substantially from season to season. However, such an
interaction effect can not be captured through ALE main effects. In Figure 12,
we illustrate the ALE main effects of the (lagged) weekly average of the daily
humidity, precipitation, and wind speed anomalies. To a smaller extent, lower-
than-baseline humidity levels from the previous week lead to more excess
deaths compared to the effect of higher-than-baseline humidity levels, while
higher precipitation levels from the previous week have the opposite effect.
The other ALE main effects do not show any strong variation across the feature
levels and their pointwise bootstrap confidence intervals contain or are close
to zero. For the weekly averages of the extreme air pollution indices in
Figure 12, we also see, to a lesser extent, an increasing pattern in the ALE
main effect. These ALE main effects are statistically significant as they do
not encompass the value zero. Finally, in Figure 12, most of the bootstrap
confidence intervals of the ALE main effects stay close to zero, indicating
that variations in the air pollution levels do not have a strong impact on
excess deaths. However, for the feature w_avg_no2_anom_l1, we observe that
high lagged nitrogen dioxide levels above its baseline have a slightly, but
significant, higher impact on excess deaths relative to low lagged nitrogen
dioxide levels. We conclude that the (lagged) hot- and cold-week index
features have the most substantial impact on excess mortality.

##### ALE regional effects.

We proceed by examining regional differences in how environmental factors
impact deviations from the mortality baseline. To achieve this, we re-compute
the ALE main effect of a specific feature by restricting the training data to
region $r$. Figures 13 and 13 visualize the results for the hot and lagged
cold-week index respectively, at a value of 1.5, see Equation (5.3). Figure 13
shows the ALE regional effect of the feature w_avg_NO2_anom_l1, i.e., the
lagged weekly average of the daily NO2 anomalies, at a value of 20 $\mu
g/m^{3}$. This corresponds to the situation wherein the daily NO2 baseline
concentrations are surpassed by 20 $\mu g/m^{3}$ on average throughout the
previous week.

()

()

()

Figure 13: ALE regional effect of the features w_avg_Tind95 (a),
w_avg_Tind5_l1 (b), and w_avg_NO2_anom_l1 (c) fitted on the training data
restricted to each region and visualized at the values 1.5, 1.5, and 20 $\mu
g/m^{3}$ respectively.

We observe that the hot-week index exhibits the most pronounced effect in the
southern NUTS 3 regions of Spain, Portugal, Italy, and certain Eastern
European countries. Conversely, northern regions appear less impacted by high
temperatures. The cold-week indicator follows a somewhat comparable pattern.
In the right panel, we see that the Central European countries appear to be
slightly more impacted by elevated air pollution in terms of PM10
concentrations.

Figure 14 illustrates the regional ALE effects for the NUTS 3 regions of
Barcelona, Milano, and Stockholm for the same features shown in Figure 13:
w_avg_Tind95, w_avg_Tind5_l1, and w_avg_NO2_anom_l1. Figure 14 shows a similar
marginal impact of the hot-week index on deviations from the mortality
baseline across the three considered regions, with a slightly lower observed
impact in Stockholm at values for the hot-week index between 1 and 2. This
could explain the behaviour in Figure 13. Figure 14 indicates that the
marginal impact of the lagged cold-week index on mortality deviations from the
baseline is clearly lower in Stockholm, and clearly higher in Barcelona. This
implies that individuals in Stockholm may tolerate extreme cold temperatures
better. For the lagged nitrogen dioxide anomalies in Figure 14, we see a
similar marginal impact across the three considered regions.

()

()

()

Figure 14: ALE regional effect in the NUTS 3 regions of Barcelona (red),
Milano (green), and Stockholm (blue) for the features w_avg_Tind95 (a),
w_avg_Tind5_l1 (b), and w_avg_NO2_anom_l1 (c) fitted on the training data
restricted to region $r$. We show a point-wise confidence interval based on 1
000 bootstrap samples of the training data.

##### ALE interaction effects.

Using the method outlined in Section 4.3, we calculate the ALE interaction
effect of the feature w_avg_T.ind95 with three other features:
w_avg_Tind95_l1, w_avg_PM10_anom, and w_avg_Tmax_anom. Figure 15 illustrates
that elevated levels of the hot-week index for both the current and preceding
week have a larger impact on excess mortality. However, when the previous week
is hot while the current week is not, we observe the opposite effect. This is
an indication of harvesting, as further investigated in Section 5.3.3.
Moreover, Figure 15 shows that heightened PM10 values in combination with a
hot week lead to higher excess mortality. Lastly, Figure 15 indicates that
when the hot-week indicator is high and maximum temperatures significantly
exceed their baseline, a higher impact is observed.

()

()

()

Figure 15: ALE interaction effects between the feature w_avg_T.ind95 and
w_avg_Tind95_l1 (a), w_avg_PM10_anom (b), and w_avg_Tmax_anom (c). The grey
squares indicate that there is no data available.

####  5.3.3 Harvesting effects

Harvesting effects occur when, for example, environmental related excess
mortality in the previous week leads to a mortality deficit in the current
week (Schwartz, 2001). We conduct an exploratory analysis to identify
potential harvesting effects within our training data spanning from 2013 to
2019, across all considered NUTS 3 regions. We then assess if we observe
similar effects using the predictions obtained from the machine learning
model, calibrated in Section 5.2.2.

To quantify the difference in observed or estimated deaths relative to the
mortality baseline for a specific week $w$ in year $t$ and region $r$, we
define the observed and estimated Excess Death Proportion (EDP) as:

$\displaystyle\text{EDP}_{t,w}^{(r)}=\dfrac{d_{t,w}^{(r)}-\hat{b}_{t,w}^{(r)}}{\hat{b}_{t,w}^{(r)}},\hskip
14.22636pt\widehat{\text{EDP}}_{t,w}^{(r)}=\dfrac{\hat{d}_{t,w}^{(r)}-\hat{b}_{t,w}^{(r)}}{\hat{b}_{t,w}^{(r)}},$

for $t\in\mathcal{T}$, $w\in\mathcal{W}_{t}$, and $r\in\mathcal{R}$.
Furthermore, $d_{t,w}^{(r)}$, $\hat{b}_{t,w}^{(r)}$, and $\hat{d}_{t,w}^{(r)}$
represent, respectively, the observed deaths, the deaths estimated by the
baseline mortality model (see Equation (3.3)), and the deaths estimated by the
machine learning model (see Equation (3.5)). Next, we consider a specific
feature $\texttt{f}_{t,w}$ from Table 3 in week $w$ of year $t$ and denote its
lagged version as $\texttt{f\\_l1}_{t,w}:=\texttt{f}_{t,w-1}$.202020Note that
$\texttt{f}_{t,0}$ represent $\texttt{f}_{t-1,53}$ (leap year) or
$\texttt{f}_{t-1,52}$ (non-leap year) We then partition the domain of this
feature into $B$ equal intervals denoted as $[a_{0},a_{1})$, $[a_{1},a_{2})$,
…, $[a_{B-1},a_{B}]$, where $a_{0}$ and $a_{B}$ denote the minimum and maximum
values observed in the training data. This creates a two-dimensional grid for
the interaction f$\times$f_l1. We calculate the average relative change in
observed deaths with respect to the mortality baseline in each grid cell
$[a_{i-1},a_{i})\times[a_{j-1},a_{j}]$, for $i,j\in\\{1,2,...,B\\}$, as:

$\displaystyle\mathcal{D}_{\texttt{f}}^{(i,j)}=\dfrac{1}{n_{i,j}}\displaystyle\sum_{r\in\mathcal{R}}\sum_{t\in\mathcal{T}}\sum_{w\in\mathcal{W}_{t}}\text{EDP}_{t,w}^{(r)}\cdot\mathbbm{1}\left\\{\texttt{f\\_l1}_{t,w}^{(r)}\in[a_{i-1},a_{i})\>\wedge\>\texttt{f}_{t,w}^{(r)}\in[a_{j-1},a_{j})\right\\},$
(5.5)

where $\mathbbm{1}\\{\cdot\\}$ is an indicator function, equalling one when
the argument holds true and zero otherwise, and $n_{i,j}$ denotes the number
of training observations falling inside the grid cell
$[a_{i-1},a_{i})\times[a_{j-1},a_{j})$. We define
$\smash{\hat{\mathcal{D}}_{\texttt{f}}^{(i,j)}}$ as the average relative
change in estimated deaths from the XGBoost model with respect to the
mortality baseline, and calculate it by using
$\smash{\widehat{\text{EDP}}_{t,w}^{(r)}}$ in Equation (5.5). Comparing the
observed and estimated relative differences provides insights into how
effectively the XGBoost model captures interaction effects present in the
training data.

Figure 16 illustrates the results of the aforementioned method for the four
extreme environmental indices w_avg_T.ind95, w_avg_T.ind5, w_avg_NO2.ind95,
and w_avg_O3.ind95. The term “harvesting effect” is employed when, e.g., high
temperatures in the current week lead to excess deaths, yet the subsequent
week with normal temperatures results in a mortality deficit compared to the
baseline. This concept applies to any weather or air pollution factor. Based
on Figure 16, the harvesting effect appears most pronounced for the hot-week
index, see Figure 16. This suggests that following a very hot week succeeded
by a week of normal temperatures, observed death counts are lower on average
by approximately 5$\%$ relative to the baseline (more green colours).
Conversely, Figure 16 reveals that the cold-week index does not exhibit a
harvesting effect in our observations, since we even observe additional excess
deaths one week after a cold week, regardless of the temperature. Moreover, no
harvesting effect is detected for the weekly extreme NO2 index (see Figure 16,
while Figure 16 shows a slight harvesting effect for the weekly extreme ozone
index. We also observe that two consecutive weeks with high air pollution
levels of NO2 or O3 result in higher excess deaths, indicating the presence of
interaction effects in the data. Lastly, we observe a close correspondence
between the observed and estimated relative differences in deaths relative to
the mortality baseline, indicating that the XGBoost model effectively captures
the interaction effects present in the data.

()

()

()

()

Figure 16: The average relative difference in observed and estimated deaths
relative to the weekly mortality baseline model, defined in grid cells of the
feature and its lagged version for w_avg_T.ind95 (a), w_avg_T.ind5 (b),
w_avg_NO2.ind95 (c), and w_avg_O3.ind95 (d).

###  5.4 Temporal and spatial aggregation

Stochastic mortality models in actuarial literature often focus on annual and
country-level population mortality. In this section, we evaluate the
advantages of initially crafting a mortality model that is region-specific,
weekly, and incorporates weather and air pollution factors, followed by
aggregating over the temporal and spatial dimensions of the observations.

##### Temporal aggregation.

We start by introducing the following notations at the NUTS 3 geographical
level and an annual time scale:

$\displaystyle
d_{t}^{(r)}=\displaystyle\sum_{w\in\mathcal{W}_{t}}d_{t,w}^{(r)},\hskip
7.11317pt\hat{b}_{t}^{(r)}=\displaystyle\sum_{w\in\mathcal{W}_{t}}\hat{b}_{t,w}^{(r)},\hskip
7.11317pt\hat{d}_{t}^{(r)}=\displaystyle\sum_{w\in\mathcal{W}_{t}}\hat{d}_{t,w}^{(r)},\hskip
7.11317ptE_{t,w}^{(r)}=\displaystyle\sum_{w\in\mathcal{W}_{t}}E_{t,w}^{(r)},$

for $t\in\mathcal{T}$ and $r\in\mathcal{R}$. Using mortality data at annual
level, a simple approach is to estimate the annual number of deaths using a
Poisson GLM, i.e., we assume:

$\displaystyle
D_{t}^{(r)}\sim\text{Poisson}\left(E_{t}^{(r)}\cdot\mu_{t}^{(r)}\right),$
(5.6)

where the annual force of mortality $\mu_{t}^{(r)}$ for each region $r$
follows:

$\displaystyle\log\mu_{t}^{(r)}=\beta_{0}^{(r)}+\beta_{1}^{(r)}t.$

Figure 17 illustrates the observed deaths (grey dots) alongside the estimated
annual death counts derived from the Poisson GLM in Equation (5.6) (green
line), and the estimated weekly deaths aggregated to an annual time scale as
derived from the weekly Serfling-type baseline model in Equation (3.1) (red
line) and the XGBoost model incorporating environmental features in Equation
(3.4) (blue line). We observe a better in-sample fit for the NUTS 3 regions of
Barcelona and Milano, indicating the usefulness at annual level of adding
region-specific environmental factors into a weekly mortality model. For
Stockholm, the improvement is less pronounced as we see a slightly worse in-
sample fit for the years 2014-2016, but a slightly better in-sample fit for
the years 2018-2019. Appendix D.4 outlines a statistical in-sample comparison
between the three considered models by means of the Poisson deviance. We find
that approximately $80\%$ of the regions favour the XGBoost model.

Figure 17: Observed annual death counts (grey dots), and the estimated annual
death counts from the Poisson GLM (green line), the weekly Serfling-type
baseline model (red line), and the XGBoost model including environmental
features (blue line). We display the results for the NUTS 3 regions Barcelona
(middle), Milano (middle), and Stockholm (right).

##### Spatial aggregation.

We aggregate the region-specific, annual observed and estimated death counts
to country-level. Hereto, let $\mathfrak{C}$ be the set of countries
considered in this case study and $\mathcal{R}_{c}$ the set of NUTS 3 regions
that belong to country $c\in\mathfrak{C}$. We define:

$\displaystyle
d_{t}^{(c)}=\displaystyle\sum_{r\in\mathcal{R}_{c}}d_{t}^{(r)},\hskip
7.11317pt\hat{b}_{t}^{(c)}=\displaystyle\sum_{r\in\mathcal{R}_{c}}\hat{b}_{t}^{(r)},\hskip
7.11317pt\hat{d}_{t}^{(c)}=\displaystyle\sum_{r\in\mathcal{R}_{c}}\hat{d}_{t}^{(r)},\hskip
7.11317ptE_{t}^{(c)}=\displaystyle\sum_{r\in\mathcal{R}_{c}}E_{t}^{(r)},$

for every $c\in\mathfrak{C}$. We aggregate the estimated death counts from the
weekly mortality baseline and the machine learning model and compare these
with the estimated death counts from a GLM using a Poisson assumption for the
annual number of deaths in each country $c\in\mathfrak{C}$:

$\displaystyle
D_{t}^{(c)}\sim\text{Poisson}\left(E_{t}^{(c)}\cdot\mu_{t}^{(c)}\right),$
(5.7)

where the annual force of mortality for each country $c$ follows:

$\displaystyle\log\mu_{t}^{(c)}=\beta_{0}^{(c)}+\beta_{1}^{(c)}t.$

Figure 17 displays the observed death counts (grey dots) alongside the
estimated annual, country-specific death counts obtained from the Poisson GLM
in Equation (5.7)(green line), and the estimated weekly deaths aggregated to a
country-level and annual time scale as derived from the weekly Serfling
baseline model (red line), and the machine learning model (blue line).
Consistent with the regional findings, the XGBoost model is the preferred
model for Spain and Italy, showcasing the usefulness of environmental factors.
When evaluating the in-sample Poisson deviance in each considered country for
each of the three models, we find that the XGBoost model including
environmental features outperforms the Serfling and simple Poisson GLM in 18
out of the 20 countries. The GLM performs best in Liechtenstein, while the
Serfling model performs best in Finland.

Figure 18: Observed annual death counts (grey dots), and the estimated
country-specific annual death counts from the Poisson GLM (green line), the
weekly Serfling-type baseline model (red line), and the XGBoost model
including environmental features (blue line). We display the results for the
countries Spain (left), Italy (middle), and Sweden (right).

###  5.5 Back-testing

We back-test our proposed weekly mortality modelling framework by projecting
the death counts in the European NUTS 3 regions for the year 2019. This
projection relies on our constructed weekly mortality model, which we re-
calibrate on weekly mortality and environmental data for the years 2013 to
2018.

Hereto, we first calibrate the weekly mortality baseline model, as presented
in Equations (3.1) and (3.2), using the weekly death counts for the years
2013-2018. This calibration follows the strategy outlined in Section 4.1. For
computational purposes, we adopt the same penalty parameter values as obtained
during the calibration process of Section 5.2.1 using death counts from 2013
to 2019. Nevertheless, we anticipate that this choice will only have a minor
impact on the calibrated baseline death counts. We denote the resulting fitted
baseline structure as $\smash{\hat{\mu}_{t,w}^{(r)}}$, for week $w$ in year
$t$ and region $r$. Subsequently, we use this calibrated baseline model to
forecast the death counts in the year 2019 as follows:

$\displaystyle\hat{b}_{2019,w}^{(r)}=E_{2019,w}\cdot\text{exp}\left(\left(\hat{\boldsymbol{\beta}}^{(r)}\right)^{T}z_{2019,w}\right),$

with $\hat{\boldsymbol{\beta}}$ the estimated parameter vector obtained from
the calibration over the years 2013-2018 and $z_{2019,w}$ the vector
consisting of the covariates in the baseline model evaluated in the year 2019,
see Section 4.1.

Next, we continue with the environmental data spanning the years 2013 to 2018,
as outlined in Section 2.2. We create environmental anomalies and extreme
environmental indices following the feature engineering process detailed in
Section 5.1. These engineered features serve as input in the machine learning
model, enabling us to estimate the excess or deficit deaths relative to the
calibrated baseline for the aforementioned years. During the calibration phase
of the machine learning model, we solely focus on the tuning parameters
related to the maximum tree depth and the number of trees and maintain the
remaining parameters at the values determined in Section 5.2.2. Following the
calibration, we proceed to construct environmental features for the year 2019,
using the calibrated environmental baseline models established from data
spanning 2013 to 2018. Subsequently, these 2019 environmental features are
used as inputs for the calibrated machine learning model to forecast death
counts across the European NUTS 3 regions for the year 2019.

()

()

()

()

()

()

Figure 19: The observed and projected death counts for the year 2019 in the
NUTS 3 regions of Barcelona (a), Milano (b), Stockholm (c), Lisboa (d), Varese
(e), and Rhône (f). The observed death counts are visualized in grey and the
death counts as obtained from the weekly baseline mortality model and the
XGBoost model in red and blue, respectively.

Figure 19 illustrate the outcomes for the NUTS 3 regions of Barcelona, Milano,
Stockholm, Lisboa, Varese, and Rhône. The results indicate a clear advantage
of the machine learning model (blue line) over the baseline model (red line)
in Milano, see Figure 19. Specifically, the XGBoost model correctly predicts
higher-than-baseline deaths in the first weeks and during the Summer of 2019
in Milano. However, in Barcelona and Stockholm in Figures 19 and 19, we
conclude that, for the year 2019, the calibrated baseline model consistently
overestimates the observed deaths and henceforth fails to accurately capture
the region-specific seasonal mortality trend. Nevertheless, the machine
learning model fulfils its task as the pattern and spikes, present in the
observed and estimated death counts from the XGBoost model, largely coincide.
This is particularly true for Barcelona, while in the region around Stockholm,
it is less pronounced. However, this aligns with our earlier observation that
the northern regions appear to be less impacted by extreme environmental
levels, see Figure 13 for example. It is worth noting that since the back-test
relies on data from only a single year, we refrain from drawing further
conclusions.

To compare the performance of our weekly mortality model, integrating
environmental features, across all considered European NUTS 3 regions, Figure
3a in Appendix D.5 displays the relative change in the Poisson deviance of the
XGBoost model relative to the one of the weekly baseline mortality model. We
find that, based on the one year-prediction results, adding environmental
features is beneficial in about 55$\%$ of the NUTS 3 regions. Northern regions
seem to experience little benefit from the integration of environmental
features, while southern regions generally benefit more. One reason for the
low percentage of regions where an improvement is observed is the inaccurately
estimated baseline trend. This is illustrated with the plots for Barcelona and
Stockholm in the year 2019 in Figures 19 and 19. Figures 19, 19, and 19
illustrate the projected death counts obtained with the weekly baseline model
and the machine learning model for the regions with the biggest relative
improvement (Lisboa and Varese) and with the biggest relative deterioration
(Rhône). In Rhône we estimate high excess deaths in the summer relative to the
baseline. This peak is also visible in the observed death counts, but starts
at a lower level. In the starting weeks of 2019, the model additionally fails
to capture the sudden drop in the death counts.

##  6 Conclusions

This paper proposes a comprehensive framework for integrating environmental
features, related to weather and air pollution factors, into a weekly
mortality model. Building upon the established weekly Serfling-type mortality
model, we enhance this so-called weekly baseline model by incorporating
environmental anomalies and extreme environmental indices through a machine
learning approach to explain deviations from the mortality baseline.

Through a case study spanning over 550 NUTS 3 regions across 20 European
countries, we find that temperature anomalies and extreme indices are
important in explaining deviations from the mortality baseline. While some air
pollution features also contribute significantly, other weather-related
variables such as humidity, rainfall, and wind speed are only of minor
importance. To address the correlation among the environmental features, we
employ Accumulated Local Effect plots to investigate the marginal impact of
each feature on the excess or deficit deaths relative to the baseline. We
hereby find that the hot-week index, indicating the frequency of days
surpassing the 95$\%$ temperature threshold within a week, exhibits the most
pronounced effect on excess deaths. Moreover, regional disparities emerge,
with southern regions estimated to have higher excess deaths associated with
both hot and cold week indices compared to northern regions. Our analysis
further unveils a harvesting effect for the hot-week index, suggesting a
mortality deficit following hot preceding weeks. Next, we demonstrate the
added contribution of environmental factors when aggregating the estimated
weekly mortality rates from the machine learning model on an annual, country-
level basis. Additionally, through back-testing, we validate the model’s
ability to accurately explain excess deaths relative to the baseline in an
out-of-sample year.

Furthermore, the proposed modelling framework could provide a valuable tool
for evaluating the potential impact of specific environmental scenarios on
mortality outcomes. By generating region-specific environmental scenarios for
temperature, weather factors, and air pollution concentrations, we can assess
the associated risks of excess deaths across different regions. Through
feature engineering and utilization of machine learning techniques like
XGBoost, we would be able to estimate risk ratios, representing multipliers
that augment or diminish the baseline number of deaths. Such findings can
inform policymakers about the potential health implications of various
environmental changes, guiding proactive measures to mitigate risks and
protect vulnerable populations.

In future research, it would be valuable to conduct a similar study on cause-
of-death mortality data to investigate how environmental factors influence
excess deaths across various causes. Additionally, we can explore the
robustness of our proposed weekly mortality modelling framework by exploring
alternative weekly mortality models as baselines, as suggested by Schöley
(2021) and by using weather-station data instead of fine-grained gridded data
sets from the Copernicus Climate Data Store and the Copernicus Atmospheric
Monitoring Service.

## References

  * Alberdi . (1998) alberdi1998dailyAlberdi, JC., Díaz, J., Montero, JC. Mirón, I.  1998\.  Daily mortality in Madrid community 1986–1992: Relationship with meteorological variables Daily mortality in Madrid community 1986–1992: Relationship with meteorological variables.  European Journal of Epidemiology146571–578.
  * Almon (1965) almon1965distributedAlmon, S.  1965\.  The distributed lag between capital appropriations and expenditures The distributed lag between capital appropriations and expenditures.  Econometrica: Journal of the Econometric Society331178–196.  10.2307/1911894
  * American Academy of Actuaries . (2016) aciAmerican Academy of Actuaries, Casualty Actuarial Society, Canadian Institute of Actuaries Society of Actuaries.  2016\.  Actuaries Climate Index. Actuaries Climate Index.  https://actuariesclimateindex.org.  Accessed 05-02-2024
  * Apley Zhu (2020) apley2020visualizingApley, DW. Zhu, J.  2020\.  Visualizing the effects of predictor variables in black box supervised learning models Visualizing the effects of predictor variables in black box supervised learning models.  Journal of the Royal Statistical Society Series B: Statistical Methodology8241059–1086.  10.1111/rssb.12377
  * Armstrong (2006) armstrong2006modelsArmstrong, B.  2006\.  Models for the relationship between ambient temperature and daily mortality Models for the relationship between ambient temperature and daily mortality.  Epidemiology176624–631.  10.1097/01.ede.0000239732.50999.8f
  * Balakrishnan . (2019) balakrishnan2019impactBalakrishnan, K., Dey, S., Gupta, T., Dhaliwal, R., Brauer, M., Cohen, AJ.Dandona, L.  2019\.  The impact of air pollution on deaths, disease burden, and life expectancy across the states of India: The Global Burden of Disease Study 2017 The impact of air pollution on deaths, disease burden, and life expectancy across the states of India: The Global Burden of Disease Study 2017.  The Lancet Planetary Health31e26–e39.  10.1016/S2542-5196(18)30261-4
  * Barnett (2010) barnett2010measureBarnett, A.  2010\.  What measure of temperature is the best predictor of mortality? What measure of temperature is the best predictor of mortality?  Environmental Research1106604–611.  10.1097/01.ede.0000362214.90491.53
  * Basu Samet (2002) basu2002relationBasu, R. Samet, JM.  2002\.  Relation between elevated ambient temperature and mortality: A review of the epidemiologic evidence Relation between elevated ambient temperature and mortality: A review of the epidemiologic evidence.  Epidemiologic Reviews242190–202.  10.1093/epirev/mxf007
  * Braga . (2001) braga2001timeBraga, ALF., Zanobetti, A. Schwartz, J.  2001\.  The time course of weather-related deaths The time course of weather-related deaths.  Epidemiology126662–667.  10.1097/00001648-200111000-00014
  * Braga . (2002) braga2002effectBraga, ALF., Zanobetti, A. Schwartz, J.  2002\.  The effect of weather on respiratory and cardiovascular deaths in 12 US cities The effect of weather on respiratory and cardiovascular deaths in 12 US cities.  Environmental Health Perspectives1109859–863.  10.1289/ehp.02110859
  * Breiman Ihaka (1984) breiman1984nonlinearBreiman, L. Ihaka, R.  1984\.  Nonlinear Discriminant Analysis Via Scaling and ACE Nonlinear discriminant analysis via scaling and ACE.  Department of Statistics, University of California.
  * Brunekreef Holgate (2002) brunekreef2002airBrunekreef, B. Holgate, ST.  2002\.  Air pollution and health Air pollution and health.  The Lancet36093411233–1242.  10.1016/S0140-6736(02)11274-8
  * Center for International Earth Science Information Network – CIESIN – Columbia University (2018) nasadatasetCenter for International Earth Science Information Network – CIESIN – Columbia University.  2018\.  Gridded Population of the World, Version 4 (GPWv4): Population Count, Revision 11. Gridded population of the world, version 4 (GPWv4): Population count, revision 11.  Palisades, New York: NASA Socioeconomic Data and Applications Center (SEDAC).  Accessed 01-08-2023  10.7927/H4JW8BX5
  * Chen Guestrin (2016) chen2016xgboostChen, T. Guestrin, C.  2016\.  XGBoost: A scalable tree boosting system XGBoost: A scalable tree boosting system.  Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining ( 785–794).  10.1145/2939672.2939785
  * Chen . (2019) chen2019packageChen, T., He, T., Benesty, M. Khotilovich, V.  2019\.  Package ‘xgboost’ Package ‘xgboost’.  R version901–66.
  * Copernicus Climate Change Service, Climate Data Store (2020) cdsCopernicus Climate Change Service, Climate Data Store.  2020\.  E-OBS daily gridded meteorological data for Europe from 1950 to present derived from in-situ observations. E-OBS daily gridded meteorological data for Europe from 1950 to present derived from in-situ observations.  Copernicus Climate Change Service (C3S) Climate Data Store (CDS).  Accessed 01-09-2023  10.24381/cds.151d3ec6
  * de Schrijver . (2021) de2021comparativede Schrijver, E., Folly, CL., Schneider, R., Royé, D., Franco, OH., Gasparrini, A. Vicedo-Cabrera, AM.  2021\.  A comparative analysis of the temperature-mortality risks using different weather datasets across heterogeneous regions A comparative analysis of the temperature-mortality risks using different weather datasets across heterogeneous regions.  GeoHealth55.  10.1029/2020GH000363
  * Devalia . (1994) devalia1994effectDevalia, JL., Rusznak, C., Herdman, MJ., Trigg, CJ., Davies, RJ. Tarraf, H.  1994\.  Effect of nitrogen dioxide and sulphur dioxide on airway response of mild asthmatic patients to allergen inhalation Effect of nitrogen dioxide and sulphur dioxide on airway response of mild asthmatic patients to allergen inhalation.  The Lancet34489381668–1671.  10.1016/S0140-6736(94)90458-8
  * DiCiccio Efron (1996) diciccio1996bootstrapDiCiccio, TJ. Efron, B.  1996\.  Bootstrap confidence intervals Bootstrap confidence intervals.  Statistical science113189–228.  10.1214/ss/1032280214
  * Dong . (2024) dong2023assessingDong, M., Bruhn, A., Shang, HL. Hui, F.  2024\.  Assessing the Impact of Climate Risk Stresses on Life Insurance Portfolios Assessing the impact of climate risk stresses on life insurance portfolios.  Asia-Pacific Journal of Risk and Insurance18187–114.  10.1515/apjri-2023-0010
  * Dong . (2022) dong2022airDong, ZM., Shang, HL. Bruhn, A.  2022\.  Air Pollution and Mortality Impacts Air pollution and mortality impacts.  Risks106126.  10.3390/risks10060126
  * EIOPA (2021) eiopa2021opinionEIOPA.  2021\.  Opinion on the Supervision of the Use of Climate Change Risk Scenarios in ORSA. Opinion on the supervision of the use of climate change risk scenarios in ORSA.  European Insurance and Occupuations Pensions Authority (EIOPA) Frankfurt.  Accessed 1-03-2024
  * Fischer Knutti (2013) fischer2013robustFischer, EM. Knutti, R.  2013\.  Robust projections of combined humidity and temperature extremes Robust projections of combined humidity and temperature extremes.  Nature Climate Change32126–130.  10.1038/nclimate1682
  * Friedman (2001) friedman2001greedyFriedman, JH.  2001\.  Greedy function approximation: A gradient boosting machine Greedy function approximation: A gradient boosting machine.  Annals of Statistics2951189–1232.  10.1214/aos/1013203451
  * Gasparrini Armstrong (2011) gasparrini2011impactGasparrini, A. Armstrong, B.  2011\.  The impact of heat waves on mortality The impact of heat waves on mortality.  Epidemiology22168.  10.1097/EDE.0b013e3181fdcd99
  * Gasparrini Armstrong (2013) gasparrini2013reducingGasparrini, A. Armstrong, B.  2013\.  Reducing and meta-analysing estimates from distributed lag non-linear models Reducing and meta-analysing estimates from distributed lag non-linear models.  BMC Medical Research Methodology1311–10.  10.1186/1471-2288-13-1
  * Gasparrini . (2010) gasparrini2010distributedGasparrini, A., Armstrong, B. Kenward, MG.  2010\.  Distributed lag non-linear models Distributed lag non-linear models.  Statistics in Medicine29212224–2234.  10.1002/sim.3940
  * Gasparrini . (2012) gasparrini2012multivariateGasparrini, A., Armstrong, B. Kenward, MG.  2012\.  Multivariate meta-analysis for non-linear and other multi-parameter associations Multivariate meta-analysis for non-linear and other multi-parameter associations.  Statistics in Medicine31293821–3839.  10.1002/sim.5471
  * Gasparrini . (2015) gasparrini2015mortalityGasparrini, A., Guo, Y., Hashizume, M., Lavigne, E., Zanobetti, A., Schwartz, J.Armstrong, B.  2015\.  Mortality risk attributable to high and low ambient temperature: A multicountry observational study Mortality risk attributable to high and low ambient temperature: A multicountry observational study.  The Lancet3869991369–375.  10.1016/S0140-6736(14)62114-0
  * Gasparrini Leone (2014) gasparrini2014attributableGasparrini, A. Leone, M.  2014\.  Attributable risk from distributed lag models Attributable risk from distributed lag models.  BMC Medical Research Methodology1411–8.  10.1186/1471-2288-14-55
  * Gervini Yohai (2002) gervini2002classGervini, D. Yohai, VJ.  2002\.  A class of robust and fully efficient regression estimators A class of robust and fully efficient regression estimators.  The Annals of Statistics302583–616.  10.1214/aos/1021379866
  * Guo . (2017) guo2017heatGuo, Y., Gasparrini, A., Armstrong, BG., Tawatsupa, B., Tobias, A., Lavigne, E.Tong, S.  2017\.  Heat wave and mortality: A multicountry, multicommunity study Heat wave and mortality: A multicountry, multicommunity study.  Environmental Health Perspectives1258087006.  10.1289/EHP1026
  * Hastie . (2001) hastie2001elementsHastie, T., Tibshirani, R. Friedman, J.  2001\.  The elements of statistical learning The elements of statistical learning.  Springer New York, NY.  10.1007/978-0-387-84858-7
  * Hulme Jenkins (1998) hulme1998climateHulme, M. Jenkins, G.  1998\.  Climate Change Scenarios for the United Kingdom – Technical Report No 1) Climate change scenarios for the United Kingdom – Technical report No 1) Publisher’s version .
  * Institut national de l’environnement industriel et des risques (Ineris) . (2022) camsdataInstitut national de l’environnement industriel et des risques (Ineris), Aarhus University, Norwegian Meteorological Institute (MET Norway), Jülich Institut für Energie- und Klimaforschung (IEK), Institute of Environmental Protection – National Research Institute (IEP-NRI), Koninklijk Nederlands Meteorologisch Instituut (KNMI)Barcelona Supercomputing Center (BSC)  2022\.  CAMS European air quality forecasts, ENSEMBLE data. CAMS European air quality forecasts, ENSEMBLE data.  https://ads.atmosphere.copernicus.eu/cdsapp#!/dataset/cams-europe-air-quality-reanalyses?tab=overview.  Copernicus Atmosphere Monitoring Service (CAMS) Atmosphere Data Store (ADS).  Accessed 01-08-2023
  * Jdanov . (2021) stmfnoteJdanov, DA., Galarza, AA., Shkolnikov, VM., Jasilionis, D., Németh, L., Leon, DA.Barbieri, M.  2021\.  The short-term mortality fluctuation data series, monitoring mortality shocks across time and space The short-term mortality fluctuation data series, monitoring mortality shocks across time and space.  Scientific Data81235.  10.1038/s41597-021-01019-1
  * Keatinge . (2000) keatinge2000heatKeatinge, WR., Donaldson, GC., Cordioli, E., Martinelli, M., Kunst, AE., Mackenbach, JP.Vuori, I.  2000\.  Heat related mortality in warm and cold regions of Europe: Observational study Heat related mortality in warm and cold regions of Europe: Observational study.  BMJ3217262670–673.  10.1136/bmj.321.7262.670
  * Li Tang (2022) li2022jointLi, H. Tang, Q.  2022\.  Joint extremes in temperature and mortality: A bivariate POT approach Joint extremes in temperature and mortality: A bivariate POT approach.  North American Actuarial Journal26143–63.  10.1080/10920277.2020.1823236
  * Lovett . (1986) lovett1986analysingLovett, AA., Bentham, CG. Flowerdew, R.  1986\.  Analysing geographic variations in mortality using Poisson regression: The example of ischaemic heart disease in England and Wales 1969–1973 Analysing geographic variations in mortality using Poisson regression: The example of ischaemic heart disease in England and Wales 1969–1973.  Social Science & Medicine2310935–943.  10.1016/0277-9536(86)90251-0
  * Meehl Tebaldi (2004) meehl2004moreMeehl, GA. Tebaldi, C.  2004\.  More intense, more frequent, and longer lasting heat waves in the 21st century More intense, more frequent, and longer lasting heat waves in the 21st century.  Science3055686994–997.  10.1126/science.1098704
  * Molnar (2019) molnar2020interpretableMolnar, C.  2019\.  Interpretable machine learning Interpretable machine learning.  [Online]. Available: https://christophm.github.io/interpretable-ml-book/.
  * Mudway Kelly (2000) mudway2000ozoneMudway, IS. Kelly, FJ.  2000\.  Ozone and the lung: A sensitive issue Ozone and the lung: A sensitive issue.  Molecular Aspects of Medicine211–48.  10.1016/s0098-2997(00)00003-0
  * Nelder Wedderburn (1972) nelder1972generalizedNelder, JA. Wedderburn, RWM.  1972\.  Generalized linear models Generalized linear models.  Journal of the Royal Statistical Society: Series A (General)1353370–384.  10.2307/2344614
  * Nielsen . (2018) nielsen2018influenzaNielsen, J., Krause, TG. Mølbak, K.  2018\.  Influenza-associated mortality determined from all-cause mortality, Denmark 2010/11-2016/17:Tthe FluMOMO model Influenza-associated mortality determined from all-cause mortality, Denmark 2010/11-2016/17:tthe FluMOMO model.  Influenza and Other Respiratory Viruses125591–604.  10.1111/irv.12564
  * Nielsen . (2019) nielsen2019europeanNielsen, J., Vestergaard, LS., Richter, L., Schmid, D., Bustos, N., Asikainen, T.Krause, TG.  2019\.  European all-cause excess and influenza-attributable mortality in the 2017/18 season: Should the burden of influenza B be reconsidered? European all-cause excess and influenza-attributable mortality in the 2017/18 season: Should the burden of influenza B be reconsidered?  Clinical Microbiology and Infection25101266–1276.  10.1016/j.cmi.2019.02.011
  * Orellano . (2020) orellano2020shortOrellano, P., Reynoso, J., Quaranta, N., Bardach, A. Ciapponi, A.  2020\.  Short-term exposure to particulate matter (PM10 and PM2.5), nitrogen dioxide (NO2), and ozone (O3) and all-cause and cause-specific mortality: Systematic review and meta-analysis Short-term exposure to particulate matter (PM10 and PM2.5), nitrogen dioxide (NO2), and ozone (O3) and all-cause and cause-specific mortality: Systematic review and meta-analysis.  Environment International142105876.  10.1016/j.envint.2020.105876
  * Orru . (2017) orru2017interplayOrru, H., Ebi, KL. Forsberg, B.  2017\.  The interplay of climate change and air pollution on health The interplay of climate change and air pollution on health.  Current Environmental Health Reports4504–513.  10.1007/s40572-017-0168-6
  * Pascal . (2014) pascal2014shortPascal, M., Falq, G., Wagner, V., Chatignoux, E., Corso, M., Blanchard, M.Larrieu, S.  2014\.  Short-term impacts of particulate matter (PM10, PM10–2.5, PM2.5) on mortality in nine French cities Short-term impacts of particulate matter (PM10, PM10–2.5, PM2.5) on mortality in nine French cities.  Atmospheric Environment95175–184.  10.1016/j.atmosenv.2014.06.030
  * Pattenden . (2003) pattenden2003mortalityPattenden, S., Nikiforov, B. Armstrong, B.  2003\.  Mortality and temperature in Sofia and London Mortality and temperature in Sofia and London.  Journal of Epidemiology and Community Health578628–633.  10.1136/jech.57.8.628
  * Pitacco (2009) pitacco2009modellingPitacco, E.  2009\.  Modelling longevity dynamics for pensions and annuity business Modelling longevity dynamics for pensions and annuity business.  Oxford University Press.
  * Pope Schwartz (1996) pope1996timePope, CA. Schwartz, J.  1996\.  Time series for the analysis of pulmonary health data Time series for the analysis of pulmonary health data.  American Journal of Respiratory and Critical Care medicine1546S229–S233.  10.1164/ajrccm/154.6_Pt_2.S229
  * Rückerl . (2011) ruckerl2011healthRückerl, R., Schneider, A., Breitner, S., Cyrys, J. Peters, A.  2011\.  Health effects of particulate air pollution: A review of epidemiological evidence Health effects of particulate air pollution: A review of epidemiological evidence.  Inhalation Toxicology2310555–592.  10.3109/08958378.2011.593587
  * Schöley (2021) Scholey2021robustnessSchöley, J.  2021\.  Robustness and bias of European excess death estimates in 2020 under varying model specifications Robustness and bias of European excess death estimates in 2020 under varying model specifications.  MedRxiv.  10.1101/2021.06.04.21258353
  * Schwartz (2000) schwartz2000distributedSchwartz, J.  2000\.  The distributed lag between air pollution and daily deaths The distributed lag between air pollution and daily deaths.  Epidemiology113320–326.  10.1097/00001648-200005000-00016
  * Schwartz (2001) schwartz2001thereSchwartz, J.  2001\.  Is there harvesting in the association of airborne particles with daily deaths and hospital admissions? Is there harvesting in the association of airborne particles with daily deaths and hospital admissions?  Epidemiology12155–61.  10.1097/00001648-200101000-00010
  * Serfling (1963) serfling1963methodsSerfling, RE.  1963\.  Methods for current statistical analysis of excess pneumonia-influenza deaths Methods for current statistical analysis of excess pneumonia-influenza deaths.  Public Health Reports786494–506.
  * Sunyer . (1996) sunyer1996airSunyer, J., Castellsagué, J., Sáez, M., Tobias, A. Antó, JM.  1996\.  Air pollution and mortality in Barcelona Air pollution and mortality in Barcelona.  Journal of Epidemiology & Community Health50Suppl 1s76–s80.  10.1136/jech.50.suppl_1.s76
  * The Eurowinter Group (1997) group1997coldThe Eurowinter Group.  1997\.  Cold exposure and winter mortality from ischaemic heart disease, cerebrovascular disease, respiratory disease, and all causes in warm and cold regions of Europe Cold exposure and winter mortality from ischaemic heart disease, cerebrovascular disease, respiratory disease, and all causes in warm and cold regions of Europe.  The Lancet34990621341–1346.  10.1016/S0140-6736(96)12338-2
  * Verhoeff . (1996) verhoeff1996airVerhoeff, AP., Hoek, G., Schwartz, J. van Wijnen, JH.  1996\.  Air pollution and daily mortality in Amsterdam Air pollution and daily mortality in Amsterdam.  Epidemiology73225–230.  10.1097/00001648-199605000-00002
  * Vestergaard . (2020) vestergaard2020excessVestergaard, LS., Nielsen, J., Richter, L., Schmid, D., Bustos, N., Braeye, T.Mølbak, K.  2020\.  Excess all-cause mortality during the COVID-19 pandemic in Europe – Preliminary pooled estimates from the EuroMOMO network, March to April 2020 Excess all-cause mortality during the COVID-19 pandemic in Europe – preliminary pooled estimates from the EuroMOMO network, March to April 2020.  Eurosurveillance25262001214.  10.2807/1560-7917.ES.2020.25.26.2001214
  * Wood (2015) wood2015packageWood, SN.  2015\.  Package “mgcv” Package “mgcv”.  R Package Version1291–7.
  * Wood (2017) wood2017generalizedWood, SN.  2017\.  Generalized additive models: An introduction with R Generalized additive models: An introduction with R.  CRC Press.  10.1201/9781315370279
  * Zanobetti . (2000) zanobetti2000generalizedZanobetti, A., Wand, MP., Schwartz, J. Ryan, LM.  2000\.  Generalized additive distributed lag models: Quantifying mortality displacement Generalized additive distributed lag models: Quantifying mortality displacement.  Biostatistics13279–292.  10.1093/biostatistics/1.3.279

##  Appendix A Methodological approaches in the epidemiological and medical
literature: a literature review

In exploring the relationship between environmental factors and mortality
statistics, various methodologies have been proposed in the epidemiological
and medical literature.

##### Minimum mortality temperature bands.

Keatinge . (2000) conduct an observational study to investigate heat-related
mortalities across Europe. Hereto, they compute the daily death rates of men
and women aged 65-74 years per million population in several regions including
North Finland, South Finland, Baden-Württemberg, Netherlands, London, North
Italy, and Athens from 1988 to 1992. They then retrieve region-specific, daily
mean temperatures from the Royal Meteorological Office and calculate the
average death rate per million population across successive temperature bands
of $3^{\circ}$C, each separated by $0.1^{\circ}$C. Hereto, let $T^{\circ}$ be
any temperature in degrees Celsius and define the set:

$\displaystyle\mathcal{S}_{r}(T^{\circ})=\left\\{(t,w,d)\mid
T^{\circ}\leq\texttt{Tavg}_{t}^{(r)}<T^{\circ}+3^{\circ}\text{C}\right\\},$

with $\texttt{Tavg}_{t}^{(r)}$ the region-specific daily mean temperature at
time $t$. The average region-specific death rate per million population in
temperature band $[T^{\circ},T^{\circ}+3^{\circ}\text{C}]$ is then defined as:

$\displaystyle\dfrac{1}{|\mathcal{S}_{r}(T^{\circ})|}\displaystyle\sum_{(t,w,d)\in\mathcal{S}(T^{\circ}}\dfrac{d_{t}^{(r)}}{E_{t}^{(r)}}\cdot
10^{6},$

where $d_{t}$ are the observed death counts at time $t$ (daily basis) and
$E_{t}$ the corresponding population exposures. Subsequently, they identify
the temperature band with the lowest death rate in each region, also called
the region-specific minimum mortality temperature (band), and find that these
bands are significantly higher in regions with hot summer temperatures, e.g.,
Athens with $[22.7^{\circ}\text{C},25.7^{\circ}\text{C}]$, compared to regions
with cold summer temperatures, e.g., north Finland with
$[14.3^{\circ}\text{C},17.3^{\circ}\text{C}]$. Subsequently, they calculate
the annual heat-related mortality in each region as the number of days the
minimum mortality temperature band was exceeded in a year, multiplied by the
mean difference of the daily death rates per million population at
temperatures above the minimum mortality temperature band with the daily death
rate per million population within this band. They find that regions with hot
summers did not exhibit significantly higher annual heat-related mortality
compared to colder regions. They conclude that European populations have
effectively acclimated to average summer temperatures, and can be expected to
adapt to the predicted global warming, see Hulme Jenkins (1998), with only a
small increase in heat-related mortality. However, such an approach typically
overlooks the non-linear relationship between temperature and the (daily)
mortality statistics.

##### Time series models.

Yet another strand in the literature relies on time series regression models.
Here, researchers typically assume an overdispersed Poisson distribution for
the daily death counts, with as an exploratory variable the observed daily
temperature values or daily measurements of a particular air pollutant,
alongside other time-varying confounding factors (Armstrong, 2006). As an
example, the structure of such a time series regression model, with the daily
average temperature Tavg as the explanatory variable of interest, is typically
represented as:

$\displaystyle\log\mathbb{E}\left[D_{t}^{(r)}\right]=\alpha^{(r)}+f\left(\texttt{Tavg}_{t}^{(r)}\mid\boldsymbol{\beta}^{(r)}\right)$
$\displaystyle+(\texttt{confounding covariates})_{t}^{(r)}+$ (A.1)
$\displaystyle(\texttt{smooth function of time})_{t}^{(r)},$

where $D_{t}^{(r)}$ is the death count random variable from a specific region
$r$ at time $t$ and $\alpha$ denotes the intercept term. Furthermore
$\smash{f^{(r)}\left(\texttt{Tavg}_{t}^{(r)}\mid\boldsymbol{\beta}\right)}$
denotes the function capturing the relationship between the region-specific
daily average temperature and the daily death counts at time $t$ and region
$r$, also known as the temperature-mortality association, parametrized through
a vector $\boldsymbol{\beta}$. Examples of confounding factors include daily
measurements of air pollution or other weather-related covariates, such as
humidity, rain fall, and wind speed, that may influence the daily death
counts. The smooth function of time is incorporated into the model to
accommodate for seasonal effects and demographic shifts (Armstrong, 2006).
While earlier studies such as The Eurowinter Group (1997) and Lovett . (1986)
favoured Poisson regression models with a linear functional form for modelling
the temperature-mortality association, later research by Braga . (2002)
suggests a more nuanced U-, V-, or J-shaped association between daily death
counts and daily temperature values, necessitating a non-linear functional
form. Hereto, Armstrong (2006) mentions that natural cubic splines are
typically used for modelling the temperature-mortality association.

Studies by Braga . (2001) and Pattenden . (2003) reveal that the impact of a
cold spell on mortality may persist for a week or longer, while the effects of
heat are more immediate. Additionally, Verhoeff . (1996) observe a positive
correlation between ozone levels and daily mortality in Amsterdam, with
effects persisting for up to two days. Ozone is a significant contributor to
smog at ground level and can harm the respiratory system (Mudway Kelly, 2000).
Sunyer . (1996) show similar findings in Barcelona for other air pollutants
such as sulphur and nitrogen dioxide. Sulphur dioxide mainly originates from
fossil fuel burnings and irritates the respiratory tract and eyes, posing
heightened risks to asthmatic patients, while nitrogen dioxide is a toxic gas
that induces respiratory issues and can lead to a reduced lung function
(Devalia ., 1994). These studies indicate that the effects of a particular
event, such as a cold spell, heat wave, or elevated air pollution levels, may
extend beyond the period in which it occurs, and possibly exhibit delayed
impacts over time.

##### Distributed lag models.

To adequately address such delayed effects, researchers have proposed
distributed lag models (DLMs). A DLM is a linear, additive regression model
for time series data that predicts the current values of a dependent variable
using both current and past (lagged) values of an explanatory variable. Almon
(1965) initially proposes the DLM to model quarterly capital expenditures in
manufacturing industries, but the framework has been applied in other fields
such as epidemiology, as highlighted by Pope Schwartz (1996). Schwartz (2000)
further refines the DLM and applies it to a case study to estimate a