# On the PM2.5 – Mortality Relationship: A Bayesian Dynamic Model for Spatio-
Temporal Confounding

Carlo Zaccardi
Department of Economics
University G. d’Annunzio of Chieti-Pescara
Pescara, Italy
<EMAIL_ADDRESS>
&Pasquale Valentini
Department of Economics
University G. d’Annunzio of Chieti-Pescara
Pescara, Italy
&Luigi Ippoliti
Department of Economics
University G. d’Annunzio of Chieti-Pescara
Pescara, Italy
&Alexandra M. Schmidt
Department of Epidemiology, Biostatistics and Occupational Health
McGill University
Montréal, Canada

###### Abstract

Spatial confounding, often regarded as a major concern in epidemiological
studies, relates to the difficulty of recovering the effect of an exposure on
an outcome when these variables are associated with unobserved factors. This
issue is particularly challenging in spatio-temporal analyses, where it has
been less explored so far. To study the effects of air pollution on mortality
in Italy, we argue that a model that simultaneously accounts for spatio-
temporal confounding and for the non-linear form of the effect of interest is
needed. To this end, we propose a Bayesian dynamic generalized linear model,
which allows for a non-linear association and for a decomposition of the
exposure effect into two components. This decomposition accommodates
associations with the outcome at fine and coarse temporal and spatial scales
of variation. These features, when combined, allow reducing the spatio-
temporal confounding bias and recovering the true shape of the association, as
demonstrated through simulation studies. The results from the real-data
application indicate that the exposure effect seems to have different
magnitudes in different seasons, with peaks in the summer. We hypothesize that
this could be due to possible interactions between the exposure variable with
air temperature and unmeasured confounders.

Keywords: Air pollution; Bayesian; Dynamic generalized linear model;
Mortality; Non-linearity; Spatio-temporal confounding.

##  1 Introduction

Air pollution is a major environmental risk factor and undermines public
health, causing a wide range of health problems, including premature deaths.
To date, there is a very rich literature demonstrating the short- and long-
term consequences of air pollution on health outcomes (see, e.g., Dominici et
al. (2002), Janes et al. (2007), Zeger et al. (2008), Greven et al. (2011),
Dai et al. (2014), Di et al. (2017), Wu et al. (2020), Brunekreef et al.
(2021b), Renzi et al. (2021), Dominici et al. (2022), Gamerman et al. (2022),
Requia et al. (2024)).

Within the ecological regression framework, one of the challenges encountered
is the determination of the shape of the exposure–outcome relationship, which
is commonly assumed to be constant and linear. We say that the relationship is
linear if the exposure enters linearly in the specification of the link
function. However, if this assumption is violated, misspecification bias would
result. We acknowledge that all models are inherently flawed, but some of them
are so misspecified that they can lead to incorrect inferential statements
(Gelman and Shalizi, 2013).

In this paper, the ultimate goal is to correctly evaluate the short-term
association between particulate matter with a diameter less than $2.5\mu m$
(PM2.5) and all-cause mortality. The data span the period between 2018 and
2022 and are available for each health district of two Italian regions (see
Section 2). We argue that the assumption of a constant and linear effect can
be unrealistic for several reasons. The physical properties and chemical
composition of the PM mixture may vary in time and space, as they are
influenced by meteorological conditions, human activities, and the morphology
of the study region (Chen et al., 2013, Masselot et al., 2022, Peng et al.,
2005). Temperature and sunlight, for example, can improve the photochemical
reactions that produce secondary aerosols (Chen et al., 2013, Peng et al.,
2005). Exceptional events such as sand storms (Chen et al., 2013) and
wildfires (Zhou et al., 2021) cause an increase in the concentration of
hazardous particles, so they can modify the association with the outcome.
Human activities, such as traffic, industrial production, domestic heating,
and air conditioning, change throughout the year, so they may also affect the
concentration and the composition of pollutants (Chen et al., 2013, Masselot
et al., 2022, Peng et al., 2005, Kelly and Fussell, 2012). The physiological
response to air pollution may also vary by season, as extreme temperatures may
increase a person’s susceptibility to PM effects (Bergmann et al., 2020, Peng
et al., 2005, Roberts, 2004, Li et al., 2017). Finally, during the summer
people spend more time outdoors or have windows open, and this increased
natural ventilation may enhance people’s exposure to air pollutants and reduce
the exposure measurement error (Bergmann et al., 2020, Peng et al., 2006,
Roberts, 2004, Stafoggia et al., 2008).

To account for the potential non-linearity of the exposure–outcome
relationship, various strategies are available. An exposure-response function
is often estimated (see, among others, Brunekreef et al., 2021b, Cork et al.,
2023, Dominici et al., 2022, Renzi et al., 2021). Alternatively, it would be
appropriate to specify a model that incorporates the inherent spatial or
temporal variations of the exposure effect (Gelfand et al., 2003). Here, we
focus on the temporal evolution of the impact of PM2.5 on human health. By
explicitly modeling the temporal variations in the coefficients, we can
analyze whether the strength of the relationship between PM2.5 and health
outcomes have changed over different time periods. Furthermore, the temporal
dynamic of the exposure effect allows investigating the influence of various
factors that have emerged since 2020, such as lockdowns, the use of facial
masks, physical distancing and social restrictions. These measures have
significantly altered people’s habits and interactions with the environment,
potentially affecting the relationship between PM2.5 exposure and health
outcomes.

Time-varying coefficients models (TVCMs) have been previously considered in
the literature. For example, the exposure effect is modeled either as a
periodic function (Chen et al., 2013, Peng et al., 2005) or as a season-
specific step function (Dai et al., 2014, Klompmaker et al., 2023). Moreover,
Chiogna and Gaetan (2002) and Lee and Shaddick (2008) employ dynamic
generalized linear models (DGLMs; West and Harrison, 1997) when analyzing the
relationship between air pollution and mortality with time-series data.

A second major challenge in ecological regression models relates to the
exposure effect being influenced by several confounders, that is factors
associated with both the exposure and the outcome variables. Although
information on some of these confounders is available, such as meteorological
variables, it is likely that others are unmeasured or difficult to quantify,
such as influenza epidemics or changes in population structure (Peng et al.,
2006). If the model does not correctly account for unmeasured confounding, the
estimator of the exposure effect is biased. We denote this situation as
spatio-temporal confounding. More specifically, we extend the notion of
spatial confounding (Paciorek, 2010) to the spatio-temporal set up.

Despite the complexity of this problem, just few authors have recently
proposed approaches that explicitly account for spatio-temporal confounding.
For instance, Reich et al. (2021) discuss extensions of causal inference
methods. Prates et al. (2022) propose the transformed multivariate Gaussian
Markov random field model (TMGMRF), an extension of the approach by Prates et
al. (2015), where Poisson intensities or Bernoulli rates are modeled directly
and without the need of a link function. A copula structure captures the
spatio-temporal dependence and guarantees the disentanglement between fixed
and random effects, hence their model should be able to control for
confounding bias. A spatio-temporal linear model is considered by Adin et al.
(2023), where the random effects are orthogonalized with respect to the fixed
effects: in other words, this work extends the restricted spatial regression
model (RSR; Hodges and Reich, 2010, Reich et al., 2006). In all these
approaches, however, the exposure–outcome association is always assumed to be
constant and linear.

A comprehensive data-driven approach that can capture the temporal dynamics of
the exposure effect while simultaneously controlling for spatio-temporal
confounding is still lacking. To tackle both problems, we propose a Bayesian
DGLM, which has several advantages. First, it is able to approximate any
smooth function of time (Gamerman and Lopes, 2006, Granger, 2008), hence the
proposed model should not be affected by misspecification bias. Second, the
proposed approach can be seen as a generalization of other methods in the
literature. For example, the approaches by Chen et al. (2013) and Peng et al.
(2005) cannot account for exceptional events such as the outbreak of the SARS-
CoV-2 pandemic. In addition to this issue, the models proposed by Dai et al.
(2014) and Klompmaker et al. (2023), among others, implicitly assume the
possibility of abrupt changes. However, the latter ones are unlikely in
practice, because the exposure effect depends, for example, on meteorological
conditions and on human activity patterns that change gradually over time
(Peng et al., 2005, Bhaskaran et al., 2013). In contrast, the proposed model
is able of overcoming both drawbacks, thanks to the ability of DGLMs to
capture gradual changes and to emphasize both seasonal and non-systematic
variations (Lütkepohl, 2005, Petris et al., 2009, Bitto and Frühwirth-
Schnatter, 2019).

A third major advantage of the proposed approach is that, following Greven et
al. (2011) and Janes et al. (2007), the exposure effect is decomposed into two
orthogonal components, each describing the (non-linear) association with the
health outcome at a different scale of variation. This decomposition should
help to distinguish the effect of interest from that of unmeasured confounding
factors. Finally, the time-varying spatial-specific intercepts in our
specification should allow capturing seasonal and trend components due to the
unmeasured confounders.

This paper is organized as follows. Section 2 offers a brief description of
the data that motivated this study and presents preliminary analyses using
competing approaches. Sections 3 and 4 provide detailed information on the
proposed model and on the inferential procedure, respectively. In Section 5, a
simulation study is conducted to investigate the ability of the proposed
approach in controlling for spatio-temporal confounding when the true
association is known. The results from the real-data application are analyzed
in Section 6. Finally, Section 7 provides a discussion.
ARXIVZaccardi˙et˙al˙ARXIVSupplementary Material is available at
zaccardi2024supplement.

##  2 Motivation

This paper analyzes the association between daily PM2.5 concentrations and
all-cause mortality counts among the elderly (i.e. people aged 65 or more) in
Piemonte and Lombardia regions (Italy) between 2018 and 2022. In this Section,
after describing the dataset, we emphasize the importance of capturing the
temporal dynamics of the regression coefficients, an aspect that is not
adequately addressed by existing approaches.

###  2.1 Data Description

The data used in this study are retrieved from publicly available sources and
full details are given in Section A of the Supplementary Material. Mortality
counts are collected from the Italian National Institute of Statistics (ISTAT,
2023a), and are observed on an irregular lattice comprising ${N}=117$ non-
overlapping areal units termed health districts (HDs). The study period
comprises ${T}=1826$ consecutive days, hence $Y_{it}$ denotes the number of
deaths (by any cause) occurred in district $i=1,\dots,{N}$ during day
$t=1,\dots,{T}$. The PM2.5 concentrations are used as a measure of the
exposure in district $i$ during day $t$, denoted as $X_{it}$, and are
collected from the Copernicus Atmosphere Monitoring Service’s (CAMS)
Atmosphere Data Store (INERIS et al., 2022, CAMS, 2023). We also acquire
information on temperature and relative humidity, which are potential
confounders, from the ERA5 reanalyses dataset available at the Copernicus
Climate Change Service’s (C3S) Climate Data Store (C3S, 2023, Hersbach et al.,
2018). The remote-sensed measurements are spatially and temporally aggregated
to match the resolution of the health data.

A measure usually employed for exploratory analyses is the standardized
mortality ratio (SMR), which is calculated as $SMR_{it}=Y_{it}/E_{it}$, for
$i=1,\dots,{N}$ and $t=1,\dots,{T}$. The denominator, $E_{it}$, denotes the
expected number of deaths, and is computed allowing for the underlying age–sex
stratification of the population within each HD (see also Section A of the
Supplementary Material). Figure 1 (left panel) depicts the average daily SMR
for each HD over the years 2018–2022, whose mean value is $1.25$ over the
whole study domain, and spatial clustering is noticeable. The right panel
shows the long-term average PM2.5 concentrations for each HD, and on average
the Lombardia region is more polluted than Piemonte ($19.41\,\mu g/m^{3}$ and
$15.18\,\mu g/m^{3}$, respectively). Other descriptive statistics are reported
in Section B of the Supplementary Material.

Figure 1: Long-term averages over the whole study period (2018–2022) of the
SMR due to all-cause mortality (left) and of the PM2.5 levels (right) for each
HD. The thicker solid line separates Piemonte (the westernmost region) from
Lombardia (the easternmost region).

Figure 2 shows the average daily SMR (on the logarithmic scale) for Lombardia
and Piemonte. The seasonal pattern is similar in both regions, with the lowest
values in summer, except for the year 2020. Two peaks can be noticed, where
the SMR more than doubled: the first one occurred in April 2020 and the other
one in November 2020. These are due to the first and second “waves” of the
SARS-CoV-2 pandemic, respectively. The first wave was more severe in
Lombardia, while the second one was more severe in Piemonte.

Figure 2: Time series of the average daily SMR for Piemonte (top) and
Lombardia (bottom) regions, on the log scale. The vertical dashed lines
separate each calendar year.

Figure 2 also highlights other two aspects. One is that the highest values of
the SMR in 2018 were reached at the beginning of January, suggesting that the
influenza epidemic reached its peak before the end of 2017. This is different
from the year 2019, where the SMR peaks were reached in February. Another
difference is that there was an increase in mortality during summer 2022.
Possible reasons may be the occurrence of prolonged heat waves (indeed, the
average daily temperature time-series has its absolute maximum in July 2022),
or the outbreak of an influenza epidemic.

###  2.2 Evidence from Competing Approaches

It is assumed that the outcome follows a Poisson distribution, such that, for
$i=1,\dots,{N}$ and $t=1,\dots,{T}$:

$Y_{it}|\mu_{it}\stackrel{{\scriptstyle
ind}}{{\sim}}\mathcal{P}oi(\mu_{it})\,,$ (1)

where $\mu_{it}$ indicates the mean. Then the linear predictor,
$\vartheta_{it}=\ln(\mu_{it})$, is modeled as a function of exposure,
$X_{it}$, and measured confounders, denoted as
${}^{m}\mathbf{Z}_{it}=({{}^{m}Z}_{1it},\dots,{{}^{m}Z}_{{p_{m}}it})^{\prime}$,
comprising basis functions of meteorological variables, and indicator
variables for weekdays, national holidays, and lockdown periods (hence
${p_{m}}=19$; see also Section A of the Supplementary Material). Additional
covariates and basis functions can be included in the linear predictor to
account for unmeasured confounding (see, e.g., Peng et al., 2006, Peng and
Dominici, 2008). Several approaches from the literature are then fitted under
the Bayesian paradigm: (i) the Null model, which does not account for
unmeasured confounding; (ii) the GLMadj model, which adjusts for unmeasured
confounding through a natural cubic spline of time (see, e.g., Peng and
Dominici, 2008, Zhou et al., 2021, Masselot et al., 2022); (iii) the JDZ model
by Janes et al. (2007), which is similar to the previous approach but further
decomposes the exposure effect in two orthogonal parts; (iv) the Dummy model
(see, e.g, Dai et al., 2014, Klompmaker et al., 2023), which allows for
season-specific exposure effects; and (v) the Periodic model by Peng et al.
(2005) and Chen et al. (2013), which assumes that the exposure effect is a
linear combination of trigonometric functions. A hierarchical structure is
assumed for each of them, with area-specific regression coefficients having
the same mean. Further details are provided in Section C of the Supplementary
Material.

Figure 3 shows, for each approach, the (overall) estimated percent change of
risk associated with a $10\ \mu g/m^{3}$ increase in PM2.5 levels, and it can
be noticed that the estimates differ considerably. The Null model yields the
highest percent change on average ($3.92\%$), whereas the JDZ model estimates
the lowest percent change on average ($0.47\%$). The models (i)–(iii) assume
that the association of interest is constant over time, however this
assumption is unlikely for the arguments outlined in Section 1.

Figure 3: Estimated percent change of risk associated with a $10\ \mu g/m^{3}$
increase in the PM2.5 concentrations under several approaches. Null: the model
that does not account for confounding; GLMadj: the model with a smooth
function of time; Dummy: the model similar to GLMadj, but it also captures
season-specific exposure effects; Periodic: the model similar to GLMadj where
the exposure effect is assumed to be a periodic function; JDZ: the model by
Janes et al. (2007). All models assume area-specific intercepts.

The only two approaches able to capture varying strengths of the association
are the Dummy and the Periodic models. They seem to agree on the temporal
evolution of the effect, showing peaks in summer and troughs in winter.
Further discussion on this seemingly unexpected finding is delayed to Section
6. The Dummy and the Periodic models, however, have some limitations, since
both assume that the effect changes within the year, but is the same across
years. This is not consistent with the evidence from the data, as the pandemic
outbreak may have altered the strength of the association. In fact, consider
Figure 4, showing the average daily PM2.5 concentrations for Lombardia and
Piemonte regions. Their temporal evolutions correspond only in part to those
of the SMR. In particular, the two peaks in the SMR in 2020 are not reflected
in the PM2.5 time series. Additionally, it can be noticed from Table B.3 that
the overall 2020 annual mean of PM2.5 levels was the lowest among the five
years considered in our analysis. In other words, the steep increase in
mortality during that year was not caused by air pollution. Instead, it is
expected that people were less exposed to fine particulate matter because of
the lockdown measures implemented in Italy to contain the spread of the SARS-
CoV-2 virus.

Although the Dummy model allows for abrupt changes in the exposure effect,
these are unlikely in practice as pointed out by Bhaskaran et al. (2013) and
Peng et al. (2005). Moreover, the results may be too sensitive to the choice
of the boundaries of each season, and a definition based on the calendar
seasons may result suboptimal.

Figure 4: Time series of the average daily PM2.5 concentrations for Piemonte
and Lombardia regions. The vertical dashed lines separate each calendar year.

##  3 The Proposed Model

For $i=1,\dots,{N}$ and $t=1,\dots,{T}$, let $Z_{it}$ denote an unmeasured
confounder. We assume that ${{}^{m}\mathbf{Z}}_{it}$ and $Z_{it}$ are
spatially and temporally smoother than the exposure, and that the outcome
follows the distribution in Equation (1). The linear predictor is then
specified as:

$\displaystyle\vartheta_{it}$
$\displaystyle=o_{it}+\beta_{0i}+f(X_{it};\bm{\gamma}_{f})+{{}^{m}\mathbf{Z}}_{it}^{\prime}\bm{\alpha}_{m}+Z_{it}+u_{it}\,,$
(2)

where $o_{it}$ is the offset, $\beta_{0i}$ is an area-specific intercept,
$f(\cdot;\cdot)$ is an unknown smooth function parameterized by the vector
$\bm{\gamma}_{f}$, $\bm{\alpha}_{m}$ is a ${p_{m}}$-dimensional vector of
coefficients, and $u_{it}$ is a spatially and temporally uncorrelated Gaussian
random variable, assumed to have zero mean and constant variance, i.e.
$u_{it}\stackrel{{\scriptstyle iid}}{{\sim}}\mathcal{N}(0,\tau_{u}^{2})$.

If all confounders were measured and $f(\cdot;\cdot)$ was modeled using basis
functions, the model specification in Equations (1)–(2) would belong to the
family of generalized additive models (GAMs; see Wood, 2017). Alternatively,
considering that $f(\cdot;\cdot)$ can be adequately approximated by using the
first-order Taylor’s expansion, we have:

$f(X_{it};\bm{\gamma}_{f})=f(\widehat{X_{t}};\bm{\gamma}_{f})+f^{\prime}(\widehat{X_{t}};\bm{\gamma}_{f})(X_{it}-\widehat{X_{t}})+\varpi_{it}\,,$
(3)

where $\widehat{X_{t}}$ denotes the large-scale (or overall) temporal trend,
$f^{\prime}(\widehat{X_{t}};\bm{\gamma}_{f})$ is the first derivative of
$f(\cdot;\bm{\gamma}_{f})$ with respect to its first argument, evaluated at
$\widehat{X_{t}}$, and $\varpi_{it}$ denotes the Taylor’s remainder. The
overall temporal trend can be calculated as a spatial average at each
$t=1,\dots,{T}$. Alternatively, $X_{it}$ can be regressed on a smooth function
of time (see also Janes et al., 2007):

$X_{it}=h(t;df)+\nu_{it}\,,$ (4)

where $df$ are the degrees of freedom of the smooth term, and $\nu_{it}$ is an
uncorrelated Gaussian random variable with zero mean and constant variance.
The overall temporal trend is then defined as the fitted values obtained from
the model (4), i.e. $\widehat{X_{t}}=\widehat{h}(t;df)$.

Equation (3) provides a decomposition of the exposure effect into two
components describing large and small (or local) scales of variation. These
components are functions of $\widehat{X_{t}}$ and $X_{it}-\widehat{X_{t}}$,
respectively. They are orthogonal because $X_{it}-\widehat{X_{t}}$ is the
residual term from the regression in Equation (4). Henceforth, the
decomposition above should help to disentangle the effect of interest from
that of unmeasured confounding factors. To clarify this, consider that
$\widehat{X_{t}}$ captures trend and seasonal components of the exposure,
therefore it is expected to be highly correlated with $Z_{it}$, given that the
latter is smooth, while the correlation between $X_{it}-\widehat{X_{t}}$ and
$Z_{it}$ should be close to zero. Therefore, the unmeasured confounder is more
likely to affect the association between outcome and exposure at large scales
(Greven et al., 2011, Janes et al., 2007), so the first term in the
decomposition of $f(X_{it};\bm{\gamma}_{f})$ in Equation (3) might be biased,
but it is less likely that the second term,
$f^{\prime}(\widehat{X_{t}};\bm{\gamma}_{f})(X_{it}-\widehat{X_{t}})$, is
affected by confounding bias.

By decomposing $f(X_{it};\bm{\gamma}_{f})$ in Equation (2), we thus have:

$\vartheta_{it}=o_{it}+\beta_{0i}+f(\widehat{X_{t}};\bm{\gamma}_{f})+f^{\prime}(\widehat{X_{t}};\bm{\gamma}_{f})(X_{it}-\widehat{X_{t}})+{{}^{m}\mathbf{Z}}_{it}^{\prime}\bm{\alpha}_{m}+Z_{it}+\tilde{u}_{it}\,.$
(5)

where $\tilde{u}_{it}=u_{it}+\varpi_{it}$ is such that
$\tilde{u}_{it}\stackrel{{\scriptstyle
iid}}{{\sim}}\mathcal{N}(0,\tilde{\tau}_{u}^{2})$.

The unknown function evaluated at the overall temporal trend of the exposure,
$f(\widehat{X_{t}};\bm{\gamma}_{f})$, and its first derivative in Equation (5)
can be seen as functions of time, therefore a time-varying coefficients can
adequately represent them (see also West and Harrison, 1997, Section 3.1).
Therefore, we propose to change the linear predictor in Equation (2) for:

$\vartheta_{it}=o_{it}+\beta_{0it}+\beta_{1t}(X_{it}-\widehat{X_{t}})+{{}^{m}\mathbf{Z}}_{it}^{\prime}\bm{\alpha}_{m}+\tilde{u}_{it}\,,$
(6)

where $\beta_{1t}=f^{\prime}(\widehat{X_{t}};\bm{\gamma}_{f})$ is the time-
varying coefficient representing the local association between exposure and
outcome (in the sense of Janes et al. (2007)), while $\beta_{0it}$ is a time-
varying area-specific intercept,

$\beta_{0it}=\beta_{0i}+f(\widehat{X_{t}};\bm{\gamma}_{f})+Z_{it}\,,$

that should account for both the large-scale exposure effect and the
unmeasured confounder.

We now provide two examples of the function $f(\cdot;\cdot)$ and of the
resulting time-varying coefficients in Equation (6). These show that our model
generalizes other methods in the literature.

######  Example 3.1.

If $f(\cdot;\bm{\gamma}_{f})$ is a constant function, i.e.
$f(X_{it};\bm{\gamma}_{f})=\gamma_{f0}$, then $\beta_{1t}=0$ and
$\beta_{0it}=\beta_{0i}+\gamma_{f0}+Z_{it}$.

######  Example 3.2.

If $f(\cdot;\bm{\gamma}_{f})$ is a linear function, i.e.
$f(X_{it};\bm{\gamma}_{f})=\gamma_{f0}+\gamma_{f1}X_{it}$, then
$\beta_{1t}=\gamma_{f1}$ and
$\beta_{0it}=\beta_{0i}+\gamma_{f0}+\gamma_{f1}\widehat{X_{t}}+Z_{it}$. The
model proposed by Janes et al. (2007) can be cast in this case, with
$\gamma_{f0}=0$. The authors, however, assume further that the intercepts are
constant, so the overall trend in the exposure is included as a covariate.
They also include a smooth function of time in their specification to control
for seasonality and long-term trends common to all areas.

###  3.1 Specification of the Time-Varying Coefficients

We implement a DGLM with random walk processes to model the temporal evolution
of the regression coefficients $\beta_{0it}$ and $\beta_{1t}$. We propose the
following state-space formulation for the latent process and the regression
coefficients, for $t=1,\dots,{T}$:

$\displaystyle\bm{\vartheta}_{t}$
$\displaystyle=\mathbf{o}_{t}+\mathbf{\tilde{X}}_{t}\tilde{\bm{\beta}}_{t}+{{}^{m}\mathbf{Z}}_{t}\bm{\alpha}_{m}+\tilde{{\bf
u}}_{t}\,,\quad\tilde{{\bf u}}_{t}\sim\mathcal{N}({\bf
0},\tilde{\tau}_{u}^{2}{\bf I}_{N})$ (7) $\displaystyle\tilde{\bm{\beta}}_{t}$
$\displaystyle=\tilde{\bm{\beta}}_{t-1}+\mathbf{w}_{t},\quad{\bf
w}_{t}\sim\mathcal{N}({\bf 0},\bm{\varSigma}_{w})\,,$ (8)

where $\bm{\vartheta}_{t}=(\vartheta_{1t},\dots,\vartheta_{{N}t})^{\prime}$,
$\mathbf{o}_{t}=(o_{1t},\dots,o_{{N}t})^{\prime}$,

${}^{m}\mathbf{Z}_{t}=\begin{bmatrix}{}^{m}\mathbf{Z}_{1t}^{\prime}\\\
\vdots\\\ {}^{m}\mathbf{Z}_{{N}t}^{\prime}\end{bmatrix}\,,$

and
$\tilde{\bm{\beta}}_{t}=\left(\beta_{01t},\dots,\beta_{0{N}t},\beta_{1t}\right)^{\prime}$
is an $(N+1)$-dimensional vector. The observation matrix is
$\mathbf{\tilde{X}}_{t}=\begin{bmatrix}\mathbf{I}_{N}&\mathbf{\tilde{X}}_{\,\bm{\cdot}t}\end{bmatrix}$,
where

$\mathbf{\tilde{X}}_{\,\bm{\cdot}t}=\begin{bmatrix}X_{1t}-\widehat{X_{t}}\\\
\vdots\\\ X_{{N}t}-\widehat{X_{t}}\end{bmatrix}\,.$

##  4 Bayesian Inference and Computations

The first stages of the Bayesian hierarchical model are defined by Equations
(1) and (7)–(8). To complete it, a prior distribution must be assigned on the
parameter set,
$\bm{\Theta}=\\{\tilde{\bm{\beta}}_{0},\bm{\alpha}_{m},\bm{\varSigma}_{w},\tilde{\tau}_{u}^{2}\\}$.

###  4.1 Prior Specification

Independent prior distributions are imposed on each parameter in
$\bm{\Theta}$. As for the vector of time-varying coefficients,
$\tilde{\bm{\beta}}_{t}$, distributional assumptions are required on the
initial state (i.e., when $t=0$), and on the innovation error covariance
matrix, $\bm{\varSigma}_{w}$. We impose that
$\tilde{\bm{\beta}}_{0}\sim\mathcal{N}(\mathbf{0},\mathbf{V}_{\beta})$, where
$\mathbf{V}_{\beta}=\begin{bmatrix}10\mathbf{I}_{N}&\mathbf{0}\\\
\mathbf{0}&4\times 10^{-6}\end{bmatrix}$.

The prior specification for $\bm{\varSigma}_{w}$ requires more care. We assume
that this matrix is block-diagonal, i.e.
$\bm{\varSigma}_{w}=\begin{bmatrix}\bm{\varSigma}_{w0}&\mathbf{0}\\\
\mathbf{0}&\sigma_{w1}^{2}\end{bmatrix}$. Given that the area-specific
intercepts should capture the effect of the unmeasured confounder, $Z_{it}$,
we expect them to be spatially correlated, so $\bm{\varSigma}_{w0}$ is an
${N}\times{N}$ matrix for which an inverse-Wishart prior is assigned, i.e.
$\bm{\varSigma}_{w0}\sim\mathcal{IW}(a_{w0},\mathbf{S}_{w0})$, with
$a_{w0}={N}+2$ and $\mathbf{S}_{w0}=k_{0}\mathbf{I}_{N}$. This specification
allows for a free-form and data-driven covariance structure, and facilitates
inference, since the posterior full conditional distribution of
$\bm{\varSigma}_{w0}$ is available in closed-form (Brown et al., 1994). An
inverse-gamma prior is imposed on $\sigma_{w1}^{2}$, i.e.
$\sigma_{w1}^{2}\sim\mathcal{IG}(a_{w1},b_{w1})$. Other specifications were
also considered (e.g., $\bm{\varSigma}_{w}$ was assumed to be diagonal), but
they did not provide good mixing of the posterior samples chains. Recall that
capturing gradual changes is a desired feature of the proposed model, and the
matrix $\bm{\varSigma}_{w}$ plays a crucial role to this end. For instance, if
$k_{0}$ is “large”, the short-term variation in the intercepts becomes
unlikely pronounced, whereas if $k_{0}$ is “small”, the “true” variation is
captured only in part. The definition of “large” and “small” mainly depends on
the sample size and on how much the prior dominates the likelihood. Note also
that extreme values of $k_{0}$ may lead to numerical non-positive definiteness
of $\bm{\varSigma}_{w0}$.

Therefore, several values for the hyperparameters $k_{0}$, $a_{w1}$ and
$b_{w1}$ are considered to evaluate the sensitivity of the posterior
distribution to the prior. In the case of both simulated and real data (see
Sections 5 and 6, respectively), we find that $k_{0}=10^{-3}$, $a_{w1}=11$,
and $b_{w1}=10^{-6}$ provide a good balance between the prior and the
likelihood.

Finally, non-informative normal priors are imposed on the coefficients of the
measured confounders,
$\bm{\alpha}_{m}\sim\mathcal{N}(\mathbf{0},10^{6}\mathbf{I}_{{p_{m}}})$.
Following Wakefield (2007), it is assumed that the residual relative risks
$\tilde{u}_{it}$, follow a log-t distribution with two degrees of freedom, and
with $95\%$ of prior mass within the interval $(0.98,1.02)$, hence we obtain
the prior $\tilde{\tau}_{u}^{2}\sim\mathcal{IG}(1,2.1\times 10^{-5})$.

###  4.2 Computations

As first step of the analysis, the overall temporal trend, $\widehat{X_{t}}$,
is obtained as the fitted values from the model in Equation (4), where
$h(t;df)$ represents a natural cubic spline of time with $df$ degrees of
freedom. This two-step strategy should be preferred to a joint modeling of the
exposure and the outcome, as it allows avoiding feedback from the outcome to
the exposure that may lead to biased estimators of the exposure effect (see,
e.g., Stephens et al., 2023).

An MCMC scheme is then implemented to simulate from the posterior distribution
of $\bm{\Theta}$. The algorithm proposed by Chan and Jeliazkov (2009) is used
to sample from the full conditionals of the time-varying coefficients. It is
an efficient simulation algorithm that allows to sample all the regression
coefficients (for all the time instants) in a single step. The remaining full
conditional distributions are composed of inverse-gamma, and inverse-Wishart
distributions. The latent process $\vartheta_{it}$ is an exception and is
sampled using a Metropolis-Hastings step. Further details can be found in
Section D of the Supplementary Material.

##  5 Analysis of Artificial Data

To validate the proposed approach, we conduct two simulation studies, both
inspired by the real data example presented in Section 2. In the first one,
which is described next, we want to assess the ability of the proposed model
to mitigate the confounding bias.

#### Simulation Study 1: Confounding Bias Mitigation

We consider ${N}=100$ contiguous HDs of Piemonte and Lombardia, then we assume
that ${T}=300$, that there are no measured confounders (i.e., ${p_{m}}=0$) and
that the exposure $X_{it}$ and the unmeasured confounder $Z_{it}$ are
correlated vector autoregressive processes of order $1$ (VAR($1$)) with a
linear correlation coefficient
$\rho_{xz}={\rm\mathbb{C}or}(X_{it},Z_{it})=0.5$, for $i=1,\dots,{N}$, and
$t=1,\dots,{T}$, and with spatially-structured, temporally independent
innovations (Cressie and Wikle, 2015, Sahu, 2022). A proper conditional
autoregressive (PCAR) specification (Banerjee et al., 2014, Besag, 1974, Sahu,
2022) is assumed to model the spatial dependence of these innovations. The
marginal spatio-temporal means of exposure and unmeasured confounder are
generated from Fourier basis expansions with period ${T}$. The correlation
matrix of $X_{it}$ is parameterized by spatial and temporal autocorrelation
parameters, i.e. $\phi_{x}^{S}\in(0,1)$ and $\phi_{x}^{T}\in(0,1)$,
respectively. Similarly, $\phi_{z}^{S}$ and $\phi_{z}^{T}$ represent the
autocorrelation parameters of $Z_{it}$, which is sampled conditional on
$X_{it}$. These four parameters represent the strength of the spatial and
temporal dependence, and define three scenarios that are summarized in Table
1. The outcome variable is sampled from the Poisson distribution in Equation
(1) using the following linear predictor:

$\vartheta_{it}=o_{it}+\beta_{0i}+\beta_{1t}X_{it}+Z_{it}+u_{it},\quad
u_{it}\stackrel{{\scriptstyle iid}}{{\sim}}\mathcal{N}(0,\tau_{u}^{2})\,,$

where $\tau_{u}^{2}=10^{-5}$. Additional information on the data generation
mechanism is provided in Section E of the Supplementary Material.

Table 1: Simulation scenarios and corresponding definitions. Scenario | Unmeasured confounder, $Z_{it}$ | $\phi_{x}^{S}$ | $\phi_{z}^{S}$ | $\phi_{x}^{T}$ | $\phi_{z}^{T}$
---|---|---|---|---|---
S1 | Smoother than exposure in space and time | 0.2 | 0.98 | 0.2 | 0.98
S2 | Smoother than exposure only in space | 0.2 | 0.98 | 0.98 | 0.2
S3 | Smoother than exposure only in time | 0.98 | 0.2 | 0.2 | 0.98

The proposed approach is compared to the Null model and to the GLMadj model,
where natural cubic spline of time with $7$ degrees of freedom is used for the
latter. Following Peng et al. (2006), the degrees of freedom for the spline
are chosen by fitting several models with different numbers of bases and
selecting the minimum number from which the estimates appear stabilized.

Figure 5 depicts the results in terms of percent change of risk associated
with a 10-unit increase in the exposure, defined as
$P_{t}=100(\exp\\{10\beta_{1t}\\}-1)$. For each scenario, the black solid line
represents the case where true exposure effect is a constant (top), a periodic
(middle), or a more complex function of time (bottom). The estimate of $P_{t}$
obtained from the proposed DGLM is shown in red, with the solid line
indicating its average, and the shaded area indicating the percentiles $2.5$
and $97.5$ of its distribution over 100 replicates.

The dashed gray and dashed green lines represent the estimated percent changes
of risk from the Null and GLMadj models, averaged over the 100 replicates.
Whereas the Null model is evidently biased in all cases, the GLMadj model and
the proposed approach significantly reduce the confounding bias in Scenario
S1, although without entirely removing it. In Scenario S2, the confounding
bias reduction is still achieved by the DGLM compared to the Null model, but
to a lesser extent. This could be because the unmeasured confounder varies at
high frequencies in time, so it cannot be captured by the time-varying
intercepts that only account for gradual changes. Also, the difference
$X_{it}-\widehat{X_{t}}$ is approximately as rough as $Z_{it}$, and, due to
their correlation, a portion of the variability of the unmeasured confounder
is explained by the local exposure effect. Scenario S3 is very similar to
Scenario S1, suggesting that the lack of smoothness in space of the unmeasured
confounder does not seem to alter the performance (in terms of confounding
bias) of the proposed approach.

Figure 5: Simulation Study 1: results are expressed in terms of percent change
of risk associated with a 10-unit increase in the exposure, $P_{t}$. The black
solid line represents the true exposure effect (whose shape differs in each
row), while the dashed gray line represents the estimated exposure effect from
the Null model. The red solid line and the shaded area represent the average
and the percentiles $2.5$ and $97.5$ of the posterior means of $P_{t}$
obtained from the proposed DGLM over 100 replicated datasets.

These results are confirmed also by Table E.4 in the Supplementary Material,
where the average mean squared error ($MSE_{avg}$), the maximum absolute bias
($MAB$), and the average mean absolute error ($MAE_{avg}$) are reported for
each model and each scenario (see the Supplementary Material, Section E.2, for
definitions). The GLMadj model has the lowest values when the exposure effect
is constant. This is likely because the GLMadj model correctly assumes a
constant exposure–outcome relationship, whereas the proposed model learns its
shape from the data. Nevertheless, in the top subplots of Figure 5, horizontal
straight lines can be drawn through the red shaded areas, indicating that the
DGLM correctly identifies a constant exposure effect. In the remaining cases,
the proposed model seems to perform best, as it is able to capture the
temporal variations of the true exposure effect.

#### Model Reparameterization

Section E.3 of the Supplementary Material also motivates, from a practical
viewpoint, the decomposition of the exposure effect into the two components
described in Equation (3). We compare the proposed DGLM with a reparameterized
model that does not remove the overall trend from the exposure. The results
suggest that the proposed approach is more effective in reducing the
confounding bias and almost always has lower $MSE_{avg}$, $MAB$ and
$MAE_{avg}$ values.

#### Simulation Study 2: Recovery of the Local Exposure Effect ($f^{\prime}$)

In the second simulation study, described in Section F of the Supplementary
Material, we want to assess whether the proposed DGLM is able to recover
$f^{\prime}(\cdot;\bm{\gamma}_{f})$, i.e. the first derivative of the unknown
functional form. The data are generated using Equations (1)–(2), and we show
that the DGLM is able to capture the dynamics of the exposure effect, but
there remains a residual bias which is only due to confounding (i.e., not due
to Taylor’s approximation).

##  6 Application to the Italian Mortality Data

We now apply the proposed DGLM to the Italian mortality dataset (see Section
3). The overall trend, $\widehat{X_{t}}$, is obtained using a natural cubic
spline of time with $4$ degrees of freedom per year. For each HD and time
instant, ${}^{m}{\bf Z}_{it}$ contains bases from a natural cubic spline of
daily temperature and daily relative humidity (with 7 and 4 degrees of
freedom, respectively), indicator variables for weekdays, for holidays, and
for the lockdown periods.

The MCMC algorithm described in Section 4 is run for 100,000 iterations, and
posterior inference is based on the final 80,000 draws. The model’s chains are
monitored for convergence visually through trace plots and using the
R-statistic of Gelman (1996) on two chains that are started from different
initial points simultaneously.

We assess the adequacy of the model fit by analyzing Bayesian $p$-values for
each HD and time instant. The chosen measure evaluates how well the model
aligns with the observed data, and involves calculating the quantity
$\Pr\bigl{(}Y_{it}>Y^{rep}_{it}\bigr{)}+\frac{1}{2}\Pr\bigl{(}Y_{it}=Y^{rep}_{it}\bigr{)}$,
where $Y^{rep}_{it}$ is generated from the predictive distribution in each
MCMC iteration. Bayesian $p$-values close to $0.5$ suggest that the generated
data align well with the model. Conversely, values approaching the boundaries
of the unit interval, i.e. $<0.1$ or $>0.9$, indicate a poor fit (Gelman et
al., 2013). The estimated mean $p$-value is around $0.5$ suggesting an
adequate fit.

Figure 6 depicts the estimated local exposure effect, expressed in terms of
percent change of risk associated with a $10\ \mu g/m^{3}$ increase in the
PM2.5 concentrations. The red solid line represents its posterior mean, while
the shaded area the corresponding 95% credible interval (CI). The exposure
effect is not constant in time, as no horizontal straight line can be drawn
within the CI. This effect, instead, exhibits varying peaks across different
years, and three phases are clearly identifiable: before (2018-19), during
(2020-21) and after (2022) the pandemic outbreak.

During the first two years, the local exposure effect exhibits a seasonal
pattern, with the peaks occurring during the summer months. The estimated
association is bounded away from zero from July to October 2018 from April to
September 2019, with the highest change in mortality risk of $7.47\%$ (95% CI:
$[3.88\%,11.30\%]$) occurring on June 29, 2019.

Figure 6: Time-varying percent change of risk associated with a $10\ \mu
g/m^{3}$ increase in the PM2.5 concentrations estimated by the proposed DGLM.
The red solid line represents its posterior mean, while the shaded area the
corresponding 95% CI.

It is worth noting that a higher percent change of mortality risk occurs when
the fine particulate matter is generally less present in the air, so this may
indicate that a threshold below which the exposure effect is null might not
exist (see, among others, Brunekreef et al., 2021b, Cork et al., 2023,
Dominici et al., 2022, Renzi et al., 2021). The seasonal behavior seems to be
confirmed by other studies (Nawrot et al., 2007, Peng et al., 2006, Roberts,
2004, Stafoggia et al., 2008, Sun et al., 2015, Tian et al., 2018).

These seasonal changes of the risk might be due to potential interaction
effects of PM2.5 concentrations with air temperature. As noted in Section 1,
increased ventilation and physiological susceptibility to air pollution may
explain the higher percent change of risk during summer. Of course,
interaction with other variables, such as other pollutants, are also possible
(see, e.g., Liu et al., 2023, Ji et al., 2022, Siddika et al., 2019).

Previous studies have also reported that particle composition may modify the
effect of PM2.5 on mortality. In fact, fine particulate matter is a complex
mixture of organic (such as black carbon and organic carbon) and inorganic
compounds (such as sulfate, nitrate, and ammonium), and the relative
proportion of these components may vary in time (Franklin et al., 2008,
Masselot et al., 2022, Vinson et al., 2020). Indeed, their precursors are
emitted from both natural (e.g., oceans and volcanoes) and anthropogenic
sources, such as fossil fuel combustion, industrial processes, and residential
heating, and their transportation and transformation to PM2.5 in the
atmosphere are mainly influenced by meteorological conditions (Hao et al.,
2023, Jing et al., 2020, Kelly and Fussell, 2012). Toxicological studies
indicate that the PM constituents have different effects on human health (see,
e.g., Kelly and Fussell, 2012), and some authors have suggested that ammonium
is probably among the most harmful components (Masselot et al., 2022, Hao et
al., 2023).

During the pandemic phase, the CI always contains the zero except for few days
in October 2020. The literature reveals a lack of consensus regarding the
association between air pollution and adverse COVID-19 health outcomes in the
general population (see, for example, Brunekreef et al., 2021a, Walton et al.,
2021). More recently, the review by Bhaskar et al. (2023) suggests a positive
association in the majority of previous studies. Our findings, instead,
indicate that the local exposure effect in 2020 and 2021 is generally not
bounded away from zero. To justify this, consider first of all that our model
uses all-cause mortality in the elderly population as the health outcome, so
it is not directly comparable to the COVID-19 studies mentioned earlier.
Secondly, the issuance of the lockdown measures by the Italian government, as
well as the mandatory use of facial masks, not only effectively curbed the
transmission of the SARS-CoV-2 virus (and of other viruses), but also
potentially reduced exposure to air pollution. During the lockdown periods,
people were required to stay at home, resulting in reduced exposure to outdoor
air pollutants. The ventilation was limited because the lockdowns occurred
during non-summer months. Moreover, the use of facial masks may have reduced
the inhalation of airborne hazardous particles (Kodros et al., 2021, Zhou et
al., 2018).

Although lockdown periods were not issued in the second half of 2021, the
local exposure effect remains not bounded away from zero. This could be due to
the fact that people continued to adopt measures such as wearing masks and
sanitizing hands. Additionally, it is possible that the phenomenon of
mortality displacement or harvesting may have occurred. This is a consequence
of the fact that the most vulnerable people may have died earlier than
expected in the preceding years due to the SARS-CoV-2 virus in 2020-2021,
resulting in an apparently null (or even protective) effect of the exposure
afterwards (Armstrong et al., 2017, Bhaskaran et al., 2013, Zanobetti and
Schwartz, 2008, Zanobetti et al., 2000).

In the third phase, people started to return to their normal activities,
resulting in increased exposure to air pollution. In fact, the local exposure
effect has a shape similar to that in 2018, but is bounded away from zero only
in July 2022 (probably because of mortality displacement), with the highest
change in mortality risk of $4.82\%$ (95% CI: $[0.65\%,9.46\%]$) occurring on
July 26.

#### Comparison with a Restricted Model

In Section G of the Supplementary Material we want to understand the
consequences of removing the measured confounders (i.e., temperature, relative
humidity, and indicators for weekdays, national holidays, and lockdown
periods) from the model specification. Therefore, we impose the restriction
that the elements of $\bm{\alpha}_{m}$ in Equation (6) are zero. The results
are insightful, as the estimated exposure effect greatly overlaps the estimate
obtained when the measured confounders are included, with differences
occurring only in summer. This suggests that the proposed DGLM can provide
accurate estimates even in complex scenarios, and that the confounding effect
may be time-varying and may be due to heat waves.

##  7 Discussion

Motivated by a real case-study, in this paper we have proposed a DGLM as a
solution for two types of bias that may affect the estimation and
interpretation of the exposure–outcome relationship: misspecification and
confounding bias. Misspecification bias arises from assuming that the exposure
effect is constant in time when it is not. When analyzing the health effects
of air pollution, Chiogna and Gaetan (2002) and Peng et al. (2005), among
others, argue that there is uncertainty about the “true” functional form of
the exposure–outcome association. Therefore, assuming a constant effect may be
too restrictive and result in biased estimates. We acknowledge this fact and
assume that there exists a smooth (but unknown) function of the exposure,
$f(\cdot;\cdot)$, describing the association of interest. Confounding bias is
due to the presence of an unmeasured confounder. To control for it, we
decompose the exposure effect into large- and small-scale components following
Janes et al. (2007), but, differently from them, we assume that the local
effect is time-varying.

We have shown that the proposed approach enjoys some advantages. On one hand,
it is able to capture any time-varying exposure effect, thus overcoming the
misspecification bias. On the other hand, thanks to the first-order Taylor’s
expansion, the DGLM can model the local exposure effect within a flexible
framework, thus mitigating the confounding bias. In this sense, the proposed
approach constitutes an extension of the model suggested by Janes et al.
(2007).

The effectiveness of the proposed approach has been demonstrated through
artificial data analyses. We have shown that it performs better than the Null
and the GLMadj models. In particular, the proposed DGLM not only can reduce
the bias, but also can capture complex temporal fluctuations of the exposure
effect and almost recover its true functional form. The bias is less
alleviated when the unmeasured variable is less smooth than the exposure in
time, but it is worth noticing that this situation is rather unrealistic in
practice. We have seen that the GLMadj model could achieve the lowest average
mean squared error in the case of a constant exposure effect (see Table E.4),
however, the shape of the exposure–outcome association is typically unknown.
Unless there is strong prior information about this shape, we thus suggest
using the proposed DGLM. Furthermore, as a proof of concept, in Section F of
the Supplementary Material, we set up a scenario where a complete absence of
confounders is assumed. In this case, the exposure effect results fully
recovered.

Many previous empirical works assume an exposure–outcome association constant
in time. However, the evidence from the analysis of the Italian mortality data
have indicated that this assumption is conceivably incorrect. In fact, the
estimated local exposure effect exhibits a seasonal pattern, with a higher
percent change of risk observed during the summer months, when the overall
PM2.5 level is lower. Similar conclusions were reached in other studies (Peng
et al., 2005, 2006, Roberts, 2004, Stafoggia et al., 2008, Sun et al., 2015,
Tian et al., 2018). Although this fact is apparently counterintuitive, it has
been observed that it might be due to possible interaction effects of PM2.5
concentrations with meteorological variables, such as temperature, and other
air pollutants. We have also noted that our result may be related to possible
effect modification due to changes in the particulate matter mixture (Franklin
et al., 2008, Hao et al., 2023, Masselot et al., 2022). The relative
proportions of the PM2.5 components have not been analyzed in this study, so
they could be the subject of future research. However, it is important to note
that Hao et al. (2023) found a non-linear long-term association between some
PM2.5 constituents and mortality.

An exception to the seasonal behavior are the years 2020 and 2021, when the
exposure effect is not bounded away from zero. This can be justified by the
fact that the lockdown measures, such as stay-at-home orders and the use of
masks, have likely reduced the exposure to air pollution. It is worth
recalling that none of the competing approaches in Section 2 were able to
capture this behavior. Additionally, our findings may contradict the results
of the majority of studies reviewed in Bhaskar et al. (2023). However, there
are differences in the definition of outcomes between our analysis and theirs,
as we focused on all-cause mortality in the elderly population. Moreover, not
all studies related to COVID-19 properly account for potential confounders, so
their findings should be interpreted with caution (Bhaskar et al., 2023,
Brunekreef et al., 2021b, Walton et al., 2021). In contrast, the proposed DGLM
is able to almost completely capture the exposure effect even in the absence
of measured confounders (see Figure G.5 in the Supplementary Material).

To account for the residual spatial autocorrelation in the outcome, possibly
due to unmeasured confounding, the proposed model can potentially provide for
the adjacency of areal units, a feature often ignored in multi-site time-
series studies (see, e.g., Dai et al., 2014, Janes et al., 2007, Peng et al.,
2005, Qiu et al., 2022, Zhou et al., 2021). In fact, this could be achieved by
including into the linear predictor a spatially structured random effect with
a conditional autoregressive specification. On the other hand, our model
cannot consider spatial variations of the local exposure effect, so it would
be interesting to specify a spatially-varying coefficients model, which would
allow discovering differences due to urbanicity and morphology of the study
area, for example.

The decomposition of the exposure effect into two orthogonal components allows
Janes et al. (2007) to test for the presence of unmeasured confounders. The
authors state that if the overall and local exposure effects are equal, then
it is likely that the confounding bias is not present. In our case, this
testing procedure is theoretically possible, but it is not straightforward to
implement in practice, since specific assumptions on the functional form of
$f(\cdot;\bm{\gamma}_{f})$ are required.

It is arguable that, even if the proposed DGLM could recover the true exposure
effect, $f(\cdot;\cdot)$, its functional form still remains unknown. In other
words, the complete list of interaction terms cannot be determined. However,
we believe that this is not a drawback only of the proposed model. In fact,
the list of interaction terms is unknown when fitting any regression model,
and it is up to the researcher to choose the most appropriate form for
$f(\cdot;\cdot)$. This decision may involve some degree of uncertainty. The
proposed approach, instead, allows for more flexibility and should be able to
capture the effect of lockdowns as well as interaction effects between the
exposure and other variables, such as other pollutants, so that inferential
problems due to misspecification are overcome.

In conclusion, the proposed DGLM is a flexible and effective approach to
estimate the exposure–outcome association in the presence of time-varying
effects and unmeasured spatio-temporal confounding. The model is able to
capture the complex dynamics of the exposure effect, and to control for
confounding. Also, in the scenario simulated without confounding, the true
effect is correctly recovered. The proposed approach is a valuable tool for
researchers interested in the health effects of air pollution, and it can be
extended to other fields where the exposure–outcome association is of
interest.

## Acknowledgements

Hersbach et al. (2018) was downloaded from the C3S (2023). INERIS et al.
(2022) was downloaded from the CAMS (2023). The results contain modified C3S
and CAMS information. Neither the European Commission nor ECMWF is responsible
for any use that may be made of the Copernicus information or data it
contains. The research of Zaccardi, Valentini and Ippoliti is funded by the
European Union - NextGenerationEU, research project PRIN2022 CoEnv - Complex
Environmental Data and Modeling (CUP 2022E3RY23). Schmidt is grateful for
financial support from the Natural Sciences and Engineering Research Council
(NSERC) of Canada (Discovery Grant RGPIN-2017-04999).

## References

  * Adin et al. (2023) A. Adin, T. Goicoa, J. S. Hodges, P. M. Schnell, and M. D. Ugarte.  Alleviating confounding in spatio-temporal areal models with an application on crimes against women in india.  _Statistical Modelling_ , 23(1):9–30, 2023.
  * Alduchov and Eskridge (1996) O. A. Alduchov and R. E. Eskridge.  Improved magnus form approximation of saturation vapor pressure.  _Journal of Applied Meteorology and Climatology_ , 35(4):601–609, 1996.
  * Analitis et al. (2008) A. Analitis, K. Katsouyanni, X. Pedeli, U. Kirchmayer, P. Michelozzi, and B. Menne.  Investigating the independent and synergistic effects of heat waves and air pollution on health: The euroheat project.  _Epidemiology_ , 19(1):S214–S215, 2008.
  * Analitis et al. (2014) A. Analitis, P. Michelozzi, D. D’Ippoliti, F. De’Donato, B. Menne, F. Matthies, R. W. Atkinson, C. Iñiguez, X. Basagaña, A. Schneider, et al.  Effects of heat waves on mortality: effect modification and confounding by air pollutants.  _Epidemiology_ , pages 15–22, 2014.
  * Analitis et al. (2018) A. Analitis, F. De’Donato, M. Scortichini, T. Lanki, X. Basagana, F. Ballester, C. Astrom, A. Paldy, M. Pascal, A. Gasparrini, et al.  Synergistic effects of ambient temperature and air pollution on health in europe: results from the phase project.  _International journal of environmental research and public health_ , 15(9):1856, 2018.
  * Armstrong et al. (2017) B. Armstrong, M. L. Bell, M. de Sousa Zanotti Stagliorio Coelho, Y.-L. Leon Guo, Y. Guo, P. Goodman, M. Hashizume, Y. Honda, H. Kim, E. Lavigne, et al.  Longer-term impact of high and low temperature on mortality: an international study to clarify length of mortality displacement.  _Environmental Health Perspectives_ , 125(10):107009, 2017.
  * Banerjee et al. (2014) S. Banerjee, B. P. Carlin, and A. E. Gelfand.  _Hierarchical modeling and analysis for spatial data_.  CRC press, 2014.
  * Bergmann et al. (2020) S. Bergmann, B. Li, E. Pilot, R. Chen, B. Wang, and J. Yang.  Effect modification of the short-term effects of air pollution on morbidity by season: A systematic review and meta-analysis.  _Science of The Total Environment_ , 716:136985, 2020.
  * Besag (1974) J. Besag.  Spatial interaction and the statistical analysis of lattice systems.  _Journal of the Royal Statistical Society: Series B (Methodological)_ , 36(2):192–225, 1974.
  * Bhaskar et al. (2023) A. Bhaskar, J. Chandra, H. Hashemi, K. Butler, L. Bennett, J. Cellini, D. Braun, and F. Dominici.  A literature review of the effects of air pollution on covid-19 health outcomes worldwide: Statistical challenges and data visualization.  _Annual Review of Public Health_ , 44:1–20, 2023.
  * Bhaskaran et al. (2013) K. Bhaskaran, A. Gasparrini, S. Hajat, L. Smeeth, and B. Armstrong.  Time series regression studies in environmental epidemiology.  _International journal of epidemiology_ , 42(4):1187–1195, 2013.
  * Bitto and Frühwirth-Schnatter (2019) A. Bitto and S. Frühwirth-Schnatter.  Achieving shrinkage in a time-varying parameter model framework.  _Journal of Econometrics_ , 210(1):75–97, 2019\.
  * Brown et al. (1994) P. J. Brown, N. D. Le, and J. V. Zidek.  Inference for a covariance matrix.  _Aspects of uncertainty: a tribute to DV Lindley_ , pages 77–92, 1994\.
  * Brunekreef et al. (2021a) B. Brunekreef, G. Downward, F. Forastiere, U. Gehring, D. Heederik, G. Hoek, M. Koopmans, L. Smit, and R. Vermeulen.  Air pollution and covid-19: Including elements of air pollution in rural areas, indoor air pollution and vulnerability and resilience aspects of our society against respiratory disease, social inequality stemming from air pollution.  _Study for the committee on Environment, Public Health and Food Safety, Policy Department for Economic, Scientific and Quality of Life Policies, European Parliament_ , PE 658.216, 2021a.  URL https://www.europarl.europa.eu/thinktank/en/document/IPOL_STU(2021)658216.
  * Brunekreef et al. (2021b) B. Brunekreef, M. Strak, J. Chen, Z. J. Andersen, R. Atkinson, M. Bauwelinck, T. Bellander, M.-C. Boutron, J. Brandt, I. Carey, et al.  Mortality and morbidity effects of long-term exposure to low-level pm2.5, bc, no2, and o3: an analysis of european cohorts in the elapse project.  _Research Reports: Health Effects Institute_ , 2021, 2021b.
  * C3S (2023) C3S.  Era5 hourly data on single levels from 1940 to present.  _Copernicus Climate Change Service (C3S) Climate Data Store (CDS)_ , 2023.  doi: 10.24381/cds.adbb2d47.  Accessed on 22-03-2023.
  * CAMS (2023) CAMS.  European air quality reanalyses.  _Copernicus Atmosphere Monitoring Service (CAMS) Atmosphere Data Store (ADS)_ , 2023.  URL https://ads.atmosphere.copernicus.eu/cdsapp#!/dataset/cams-europe-air-quality-reanalyses?tab=overview.  Accessed on 22-03-2023.
  * Chan and Jeliazkov (2009) J. C. Chan and I. Jeliazkov.  Efficient simulation and integrated likelihood estimation in state space models.  _International Journal of Mathematical Modelling and Numerical Optimisation_ , 1(1-2):101–120, 2009.
  * Chen et al. (2013) R. Chen, R. D. Peng, X. Meng, Z. Zhou, B. Chen, and H. Kan.  Seasonal variation in the acute effect of particulate air pollution on mortality in the china air pollution and health effects study (capes).  _Science of the total environment_ , 450:259–265, 2013\.
  * Chiogna and Gaetan (2002) M. Chiogna and C. Gaetan.  Dynamic generalized linear models with application to environmental epidemiology.  _Journal of the Royal Statistical Society Series C: Applied Statistics_ , 51(4):453–468, 2002.
  * Cork et al. (2023) M. Cork, D. Mork, and F. Dominici.  Methods for estimating the exposure-response curve to inform the new safety standards for fine particulate matter.  _arXiv preprint arXiv:2306.03011_ , 2023.
  * Cressie and Wikle (2015) N. Cressie and C. K. Wikle.  _Statistics for spatio-temporal data_.  John Wiley & Sons, 2015.
  * Dai et al. (2014) L. Dai, A. Zanobetti, P. Koutrakis, and J. D. Schwartz.  Associations of fine particulate matter species with mortality in the united states: a multicity time-series analysis.  _Environmental health perspectives_ , 122(8):837–842, 2014.
  * Di et al. (2017) Q. Di, L. Dai, Y. Wang, A. Zanobetti, C. Choirat, J. D. Schwartz, and F. Dominici.  Association of short-term exposure to air pollution with mortality in older adults.  _Journal of the American Medical Association_ , 318(24):2446–2456, 2017.
  * Dominici et al. (2002) F. Dominici, M. Daniels, S. L. Zeger, and J. M. Samet.  Air pollution and mortality: estimating regional and national dose-response relationships.  _Journal of the American Statistical Association_ , 97(457):100–111, 2002.
  * Dominici et al. (2022) F. Dominici, A. Zanobetti, J. Schwartz, D. Braun, B. Sabath, and X. Wu.  Assessing adverse health effects of long-term exposure to low levels of ambient air pollution: Implementation of causal inference methods.  _Research Reports: Health Effects Institute_ , 2022(211):1–56, 2022.
  * Franklin et al. (2008) M. Franklin, P. Koutrakis, and J. Schwartz.  The role of particle composition on the association between pm2.5 and mortality.  _Epidemiology_ , 19(5):680, 2008.
  * Gamerman and Lopes (2006) D. Gamerman and H. F. Lopes.  _Markov chain Monte Carlo: stochastic simulation for Bayesian inference_.  CRC press, 2006.
  * Gamerman et al. (2022) D. Gamerman, L. Ippoliti, and P. Valentini.  A dynamic structural equation approach to estimate the short-term effects of air pollution on human health.  _Journal of the Royal Statistical Society Series C: Applied Statistics_ , 71(3):739–769, 2022.
  * Gelfand et al. (2003) A. E. Gelfand, H.-J. Kim, C. Sirmans, and S. Banerjee.  Spatial modeling with spatially varying coefficient processes.  _Journal of the American Statistical Association_ , 98(462):387–396, 2003.
  * Gelman (1996) A. Gelman.  Inference and monitoring convergence.  _Chapter in Gilks, Richardson, and Spiegelhalter (1996)_ , 1996.
  * Gelman and Shalizi (2013) A. Gelman and C. R. Shalizi.  Philosophy and the practice of bayesian statistics.  _British Journal of Mathematical and Statistical Psychology_ , 66(1):8–38, 2013.
  * Gelman et al. (2013) A. Gelman, J. B. Carlin, H. S. Stern, D. B. Dunson, A. Vehtari, and D. B. Rubin.  _Bayesian data analysis, third edition_.  CRC Press, Nov 2013.  ISBN 9781439898208.
  * Granger (2008) C. W. Granger.  Non-linear models: Where do we go next-time varying parameter models?  _Studies in Nonlinear Dynamics & Econometrics_, 12(3), 2008.
  * Greven et al. (2011) S. Greven, F. Dominici, and S. Zeger.  An approach to the estimation of chronic air pollution effects using spatio-temporal information.  _Journal of the American Statistical Association_ , 106(494):396–406, 2011.
  * Hao et al. (2023) H. Hao, Y. Wang, Q. Zhu, H. Zhang, A. Rosenberg, J. Schwartz, H. Amini, A. van Donkelaar, R. Martin, P. Liu, et al.  National cohort study of long-term exposure to pm2.5 components and mortality in medicare american older adults.  _Environmental Science & Technology_, 57(17):6835–6843, 2023.
  * Hersbach et al. (2018) H. Hersbach, B. Bell, P. Berrisford, G. Biavati, A. Horányi, J. Muñoz Sabater, J. Nicolas, C. Peubey, R. Radu, I. Rozum, et al.  ERA5 hourly data on single levels from 1940 to present.  _Copernicus Climate Change Service (C3S) Climate Data Store (CDS)_ , 2018.  doi: 10.24381/cds.adbb2d47.  Accessed on 22-03-2023.
  * Hodges and Reich (2010) J. S. Hodges and B. J. Reich.  Adding spatially-correlated errors can mess up the fixed effect you love.  _The American Statistician_ , 64(4):325–334, 2010\.
  * INERIS et al. (2022) INERIS, Aarhus University, MET Norway, IEK, IEP-NRI KNMI, METEO FRANCE, TNO, SMHI, et al.  CAMS European air quality forecasts, ENSEMBLE data.  _Copernicus Atmosphere Monitoring Service (CAMS) Atmosphere Data Store (ADS)_ , 2022.  URL https://ads.atmosphere.copernicus.eu/cdsapp#!/dataset/cams-europe-air-quality-reanalyses?tab=overview.  Accessed on 22-03-2023.
  * ISTAT (2023a) ISTAT.  Decessi e cause di morte: cosa procude l’istat, 2023a.  URL https://www.istat.it/it/archivio/240401.  Istituto Nazionale di Statistica. Accessed on 22-03-2023.
  * ISTAT (2023b) ISTAT.  Tavole di mortalità: Classi di età, 2023b.  URL http://dati.istat.it/Index.aspx?DataSetCode=DCIS_MORTALITA1.  Istituto Nazionale di Statistica. Accessed on 20-03-2023.
  * ISTAT (2023c) ISTAT.  Popolazione residente al 1∘ gennaio, 2023c.  URL http://dati.istat.it/Index.aspx?DataSetCode=DCIS_POPRES1.  Istituto Nazionale di Statistica. Accessed on 22-03-2023.
  * Janes et al. (2007) H. Janes, F. Dominici, and S. L. Zeger.  Trends in air pollution and mortality: an approach to the assessment of unmeasured confounding.  _Epidemiology_ , pages 416–423, 2007.
  * Ji et al. (2022) J. S. Ji, L. Liu, J. Zhang, H. Kan, B. Zhao, K. G. Burkart, and Y. Zeng.  No2 and pm2. 5 air pollution co-exposure and temperature effect modification on pre-mature mortality in advanced age: a longitudinal cohort study in china.  _Environmental Health_ , 21(1):97, 2022.
  * Jing et al. (2020) Z. Jing, P. Liu, T. Wang, H. Song, J. Lee, T. Xu, and Y. Xing.  Effects of meteorological factors and anthropogenic precursors on pm2.5 concentrations in cities in china.  _Sustainability_ , 12(9):3550, 2020.
  * Kelly and Fussell (2012) F. J. Kelly and J. C. Fussell.  Size, source and chemical composition as determinants of toxicity attributable to ambient particulate matter.  _Atmospheric environment_ , 60:504–526, 2012.
  * Klompmaker et al. (2023) J. O. Klompmaker, F. Laden, P. James, M. B. Sabath, X. Wu, J. Schwartz, F. Dominici, A. Zanobetti, and J. E. Hart.  Effects of long-term average temperature on cardiovascular disease hospitalizations in an american elderly population.  _Environmental research_ , 216:114684, 2023.
  * Kodros et al. (2021) J. K. Kodros, K. O’Dell, J. M. Samet, C. L’Orange, J. R. Pierce, and J. Volckens.  Quantifying the health benefits of face masks and respirators to mitigate exposure to severe air pollution.  _GeoHealth_ , 5(9):e2021GH000482, 2021.
  * Lee and Shaddick (2008) D. Lee and G. Shaddick.  Modelling the effects of air pollution on health using bayesian dynamic generalised linear models.  _Environmetrics_ , 19(8):785–804, 2008.
  * Li et al. (2017) J. Li, A. Woodward, X.-Y. Hou, T. Zhu, J. Zhang, H. Brown, J. Yang, R. Qin, J. Gao, S. Gu, et al.  Modification of the effects of air pollutants on mortality by temperature: a systematic review and meta-analysis.  _Science of the total environment_ , 575:1556–1570, 2017\.
  * Liu et al. (2023) C. Liu, R. Chen, F. Sera, A. M. Vicedo-Cabrera, Y. Guo, S. Tong, E. Lavigne, P. M. Correa, N. V. Ortega, S. Achilleos, et al.  Interactive effects of ambient fine particulate matter and ozone on daily mortality in 372 cities: two stage time series analysis.  _The British Medical Journal_ , 383, 2023.
  * Lütkepohl (2005) H. Lütkepohl.  _New introduction to multiple time series analysis_.  Springer Science & Business Media, 2005.
  * Masselot et al. (2022) P. Masselot, F. Sera, R. Schneider, H. Kan, É. Lavigne, M. Stafoggia, A. Tobias, H. Chen, R. T. Burnett, J. Schwartz, et al.  Differential mortality risks associated with pm2.5 components: a multi-country, multi-city study.  _Epidemiology_ , 33(2):167–175, 2022.
  * Nawrot et al. (2007) T. Nawrot, R. Torfs, F. Fierens, S. De Henauw, P. Hoet, G. Van Kersschaever, G. De Backer, and B. Nemery.  Stronger associations between daily mortality and fine particulate air pollution in summer than in winter: evidence from a heavily polluted region in western europe.  _Journal of Epidemiology & Community Health_, 61(2):146–149, 2007.
  * Paciorek (2010) C. J. Paciorek.  The importance of scale for spatial-confounding bias and precision of spatial regression estimators.  _Statistical Science_ , 25(1):107, 2010.
  * Peng and Dominici (2008) R. D. Peng and F. Dominici.  Statistical methods for environmental epidemiology with r.  _R: a case study in air pollution and health_ , 2008.
  * Peng et al. (2005) R. D. Peng, F. Dominici, R. Pastor-Barriuso, S. L. Zeger, and J. M. Samet.  Seasonal analyses of air pollution and mortality in 100 us cities.  _American journal of epidemiology_ , 161(6):585–594, 2005.
  * Peng et al. (2006) R. D. Peng, F. Dominici, and T. A. Louis.  Model choice in time series studies of air pollution and mortality.  _Journal of the Royal Statistical Society: Series A (statistics in society)_ , 169(2):179–203, 2006.
  * Petris et al. (2009) G. Petris, S. Petrone, and P. Campagnoli.  _Dynamic linear models with R_.  Springer Science & Business Media, 2009.
  * Prates et al. (2015) M. O. Prates, D. K. Dey, M. R. Willig, and J. Yan.  Transformed gaussian markov random fields and spatial modeling of species abundance.  _Spatial Statistics_ , 14:382–399, 2015.
  * Prates et al. (2022) M. O. Prates, D. R. Azevedo, Y. C. MacNab, and M. R. Willig.  Non-separable spatio-temporal models via transformed multivariate gaussian markov random fields.  _Journal of the Royal Statistical Society Series C: Applied Statistics_ , 71(5):1116–1136, 2022.
  * Qiu et al. (2022) X. Qiu, M. Danesh-Yazdi, Y. Wei, Q. Di, A. Just, A. Zanobetti, M. Weisskopf, F. Dominici, and J. Schwartz.  Associations of short-term exposure to air pollution and increased ambient temperature with psychiatric hospital admissions in older adults in the usa: a case–crossover study.  _The Lancet Planetary Health_ , 6(4):e331–e341, 2022.
  * Reich et al. (2006) B. J. Reich, J. S. Hodges, and V. Zadnik.  Effects of residual smoothing on the posterior of the fixed effects in disease-mapping models.  _Biometrics_ , 62(4):1197–1206, 2006.
  * Reich et al. (2021) B. J. Reich, S. Yang, Y. Guan, A. B. Giffin, M. J. Miller, and A. Rappold.  A review of spatial causal inference methods for environmental and epidemiological applications.  _International Statistical Review_ , 89(3):605–634, 2021.
  * Renzi et al. (2021) M. Renzi, S. Marchetti, F. De’Donato, M. Pappagallo, M. Scortichini, M. Davoli, L. Frova, P. Michelozzi, and M. Stafoggia.  Acute effects of particulate matter on all-cause mortality in urban, rural, and suburban areas, italy.  _International Journal of Environmental Research and Public Health_ , 18(24):12895, 2021.
  * Requia et al. (2024) W. J. Requia, A. M. Vicedo-Cabrera, H. Amini, and J. D. Schwartz.  Short-term air pollution exposure and mortality in brazil: Investigating the susceptible population groups.  _Environmental Pollution_ , 340:122797, 2024.
  * Roberts (2004) S. Roberts.  Interactions between particulate air pollution and temperature in air pollution mortality time series studies.  _Environmental research_ , 96(3):328–337, 2004\.
  * Sahu (2022) S. K. Sahu.  _Bayesian modeling of spatio-temporal data with R_.  CRC Press, 2022.
  * Siddika et al. (2019) N. Siddika, A. K. Rantala, H. Antikainen, H. Balogun, A. K. Amegah, N. R. Ryti, J. Kukkonen, M. Sofiev, M. S. Jaakkola, and J. J. Jaakkola.  Synergistic effects of prenatal exposure to fine particulate matter (pm2.5) and ozone (o3) on the risk of preterm birth: A population-based cohort study.  _Environmental research_ , 176:108549, 2019.
  * Stafoggia et al. (2008) M. Stafoggia, J. Schwartz, F. Forastiere, and C. Perucci.  Does temperature modify the association between air pollution and mortality? a multicity case-crossover analysis in italy.  _American journal of epidemiology_ , 167(12):1476–1485, 2008.
  * Stephens et al. (2023) D. A. Stephens, W. S. Nobre, E. E. Moodie, and A. M. Schmidt.  Causal inference under mis-specification: Adjustment based on the propensity score (with discussion).  _Bayesian Analysis_ , 18(2):639–694, 2023.
  * Sun et al. (2015) S. Sun, P. Cao, K.-P. Chan, H. Tsang, C.-M. Wong, and T.-Q. Thach.  Temperature as a modifier of the effects of fine particulate matter on acute mortality in hong kong.  _Environmental Pollution_ , 205:357–364, 2015.
  * Tian et al. (2018) L. Tian, F. Liang, Q. Guo, S. Chen, S. Xiao, Z. Wu, X. Jin, and X. Pan.  The effects of interaction between particulate matter and temperature on mortality in beijing, china.  _Environmental Science: Processes & Impacts_, 20(2):395–405, 2018.
  * Vinson et al. (2020) A. Vinson, A. Sidwell, O. Black, and C. Roper.  Seasonal variation in the chemical composition and oxidative potential of pm2.5.  _Atmosphere_ , 11(10):1086, 2020.
  * Wakefield (2007) J. Wakefield.  Disease mapping and spatial regression with count data.  _Biostatistics_ , 8(2):158–183, 2007.
  * Walton et al. (2021) H. Walton, D. Evangelopoulos, M. Kasdagli, L. Selley, D. Dajnak, and K. Katsouyanni.  Investigating links between air pollution, covid-19 and lower respiratory infectious diseases.  _Research Report of the Environmental Research Group, Imperial College London_ , 2021.  URL https://www.imperial.ac.uk/school-public-health/environmental-research-group/research/air-pollution-epidemiology/air-pollution-and-covid-19/.
  * West and Harrison (1997) M. West and J. Harrison.  _Bayesian Forecasting and Dynamic Models_.  Springer Series in Statistics. Springer-Verlag, New York, 1997.  ISBN 9780387947259.  doi: 10.1007/b98971.
  * Wood (2017) S. N. Wood.  _Generalized additive models: an introduction with R_.  CRC press, 2017.
  * Wu et al. (2020) X. Wu, D. Braun, J. Schwartz, M. Kioumourtzoglou, and F. Dominici.  Evaluating the impact of long-term exposure to fine particulate matter on mortality among the elderly.  _Science advances_ , 6(29):eaba5692, 2020.
  * Zanobetti and Schwartz (2008) A. Zanobetti and J. Schwartz.  Mortality displacement in the association of ozone with mortality: an analysis of 48 cities in the united states.  _American Journal of Respiratory and Critical Care Medicine_ , 177(2):184–189, 2008.
  * Zanobetti et al. (2000) A. Zanobetti, M. P. Wand, J. Schwartz, and L. M. Ryan.  Generalized additive distributed lag models: quantifying mortality displacement.  _Biostatistics_ , 1(3):279–292, 2000.
  * Zeger et al. (2008) S. L. Zeger, F. Dominici, A. McDermott, and J. M. Samet.  Mortality in the medicare population and chronic exposure to fine particulate air pollution in urban centers (2000–2005).  _Environmental Health Perspectives_ , 116(12):1614–1619, 2008.
  * Zhou et al. (2018) S. S. Zhou, S. Lukula, C. Chiossone, R. W. Nims, D. B. Suchmann, and M. K. Ijaz.  Assessment of a respiratory face mask for capturing air pollutants and pathogens including human influenza and rhinoviruses.  _Journal of thoracic disease_ , 10(3):2059, 2018\.
  * Zhou et al. (2021) X. Zhou, K. Josey, L. Kamareddine, M. C. Caine, T. Liu, L. J. Mickley, M. Cooper, and F. Dominici.  Excess of covid-19 cases and deaths due to fine particulate matter exposure during the 2020 wildfires in the united states.  _Science advances_ , 7(33):eabi8789, 2021.

## Supplementary Materials

##  Appendix A Data Collection

We gather daily data to analyze the short-term association between PM2.5 and
all-cause mortality in Piemonte and Lombardia regions for people aged 65 or
more. The data span the period between 2018 and 2022, and are collected at the
health district (HD) level.

  * •

Mortality data. All-cause mortality data are collected from the Italian
National Institute of Statistics (ISTAT, 2023a) for each municipality in the
two regions, and are then spatially aggregated to the HD level, which
represents the areal unit. The spatial domain comprises ${N}=117$ HDs where
$84$ belong to the Lombardia region and $33$ to the Piemonte region. The
municipality of Campione d’Italia is excluded from the analysis since it has a
small area and the time series contains many zeroes. Since the data refer to
five consecutive years, the temporal dimension is ${T}=1826$.

For $i=1,\dots,{N}$ and $t=1,\dots,{T}$, the standardized mortality ratio
(SMR) is calculated as the ratio between the observed number of deaths,
$Y_{it}$, and the expected number of deaths, $E_{it}$. The latter is obtained
by indirect standardization (Banerjee et al., 2014, Sahu, 2022) as follows.
The resident population aged 65 or more represents the population at risk and
is obtained from the Italian National Institute of Statistics for each year
(ISTAT, 2023c). The population is stratified by age class (65–69, 70–74,
75–79, 80–84, 85–89, 90–94, or 95+ years) and gender (male or female). For
each stratum $k=1,\dots,K$, the mortality rate for the northwest Italy is
retrieved from the Italian National Institute of Statistics for each year
(ISTAT, 2023b). The expected counts are thus computed as follows:

$E_{it}=\sum_{k=1}^{K}Pop_{itk}r_{tk}\,,$

where $Pop_{itk}$ and $r_{tk}$ are the population at risk and the mortality
rate in the $k$-th stratum, respectively.

  * •

Exposure data. Data on PM2.5 concentrations (expressed in $\mu g/m^{3}$) are
collected from the Copernicus Atmosphere Monitoring Service’s (CAMS)
Atmosphere Data Store (INERIS et al., 2022, CAMS, 2023). Hourly data are
available at a spatial resolution of $0.1^{\circ}\times 0.1^{\circ}$, so we
use zonal statistics to aggregate the estimates to the HD level, and then
average across time to obtain daily data. To evaluate the short-term effect of
interest, a lag of one day is applied to the exposure’s time series.

  * •

Potential confounders. Daily temperature, daily relative humidity, indicator
variables for the days of the week, for national holidays and for lockdown
periods are treated as potential measured confounders. We gather satellite
measurements of air temperature and dew-point temperature (both expressed in
Kelvin degrees) from the ERA5 reanalyses dataset available at the Copernicus
Climate Change Service’s (C3S) Climate Data Store (C3S, 2023, Hersbach et al.,
2018). Relative humidity ($RH$) is calculated from air temperature ($Temp$)
dew-point temperature ($DT$) using the August-Roche-Magnus approximation
formula (Alduchov and Eskridge, 1996):

$RH=100\times\exp\left(\frac{17.625\times DT}{243.04+DT}-\frac{17.625\times
Temp}{243.04+Temp}\right)\,.$

These hourly data are spatially and temporally averaged, similarly to the
exposure.

The lockdowns are defined as follows: the first period is from March 9, 2020,
to May 3, 2020, the second one is from December 21, 2020, to January 15, 2021,
and the third one from March 6, 2021, to April 25, 2021.

##  Appendix B Descriptive Statistics

Descriptive statistics on the SMR data in all HDs are reported in Table B.1.
The average SMR is $1.28$ and $1.18$ in Lombardia and Piemonte, respectively,
and the standard deviation is slightly higher in Lombardia ($0.34$) than in
Piemonte ($0.26$).

Table B.1: Descriptive statistics of daily SMR across all HDs in 2018–2022. Q1, Q2, and Q3 stand for first, second, and third quartiles, respectively. Region | Mean | SD | Min | Q1 | Q2 | Q3 | Max
---|---|---|---|---|---|---|---
Lombardia | 1.28 | 0.29 | 0.75 | 1.14 | 1.25 | 1.39 | 3.52
Piemonte | 1.18 | 0.22 | 0.68 | 1.03 | 1.14 | 1.29 | 2.10
Overall | 1.25 | 0.26 | 0.75 | 1.12 | 1.22 | 1.37 | 3.09

Descriptive statistics on PM2.5 levels in all HDs during the whole study
period are reported in Table B.2. The average daily PM2.5 concentration is
$18.21\,\mu g/m^{3}$ (standard deviation = $10.36$) over the entire spatial
domain, but Lombardia is more polluted than Piemonte. Table B.2 also reports
descriptive statistics of temperature and relative humidity. The average daily
temperature is $284.91\,K$ ($7.60$), while the average daily relative humidity
is $91.76\%$ ($3.59$). These are very similar to the region-specific values
(not shown).

Table B.2: Descriptive statistics of daily PM2.5 concentrations, temperature, and relative humidity across all HDs in 2018–2022. Q1, Q2, and Q3 stand for first, second, and third quartiles, respectively. Variable | Mean | SD | Min | Q1 | Q2 | Q3 | Max
---|---|---|---|---|---|---|---
PM2.5 ($\mu g/m^{3}$) | 18.21 | 10.36 | 2.32 | 10.49 | 15.20 | 24.12 | 63.07
Temperature ($K$) | 284.91 | 7.60 | 265.15 | 278.17 | 284.68 | 292.02 | 301.59
Relative humidity (%) | 91.76 | 3.59 | 79.55 | 89.43 | 92.12 | 94.54 | 98.43

Finally, Table B.3 reports the annual death counts and annual averages of the
daily SMR and daily PM2.5 for each administrative region.

Table B.3: Annual mortality counts, and annual averages of daily SMR and daily PM2.5 concentrations for Piemonte and Lombardia (2018–2022). |  | Year
---|---|---
Variable | Region | 2018 | 2019 | 2020 | 2021 | 2022
Mortality Counts | Lombardia | 89742 | 90613 | 123586 | 97365 | 101069
Piemonte | 48739 | 48479 | 60364 | 51202 | 53581
Overall | 138481 | 139092 | 183950 | 148567 | 154650
SMR | Lombardia | 1.32 | 1.32 | 1.23 | 1.29 | 1.26
Piemonte | 1.20 | 1.20 | 1.16 | 1.18 | 1.17
Overall | 1.29 | 1.29 | 1.21 | 1.26 | 1.24
PM2.5 ($\mu g/m^{3}$) | Lombardia | 18.56 | 20.41 | 18.36 | 19.42 | 20.30
Piemonte | 14.00 | 16.47 | 14.33 | 15.17 | 15.91
Overall | 17.27 | 19.30 | 17.22 | 18.22 | 19.06

##  Appendix C Details on Competing Models

The proposed DGLM is compared to several approaches from the literature. These
are fitted under the Bayesian paradigm, and the general hierarchical structure
is as follows, for $i=1,\dots,{N}$ and $t=1,\dots,{T}$:

$\displaystyle Y_{it}|\vartheta_{it}$
$\displaystyle\sim\mathcal{P}oi(\exp(\vartheta_{it}))$
$\displaystyle\vartheta_{it}$
$\displaystyle=o_{it}+\mathbf{A}_{it}^{\prime}\mathbf{a}_{i}+\varepsilon_{it},\quad\varepsilon_{it}\stackrel{{\scriptstyle
iid}}{{\sim}}\mathcal{N}(0,\tau_{\varepsilon}^{2})$
$\displaystyle\mathbf{a}_{i}$
$\displaystyle\sim\mathcal{N}(\bm{\mu}_{a},\bm{\varSigma}_{a})$
$\displaystyle\bm{\mu}_{a}$ $\displaystyle\sim\mathcal{N}(0,10^{6}{\bf
I}_{p_{a}})$ $\displaystyle\bm{\varSigma}_{a}$
$\displaystyle\sim\mathcal{IW}(p_{a}+1.1,10^{-4}{\bf I}_{p_{a}})$
$\displaystyle\tau_{\varepsilon}^{2}$
$\displaystyle\sim\mathcal{IG}(1,2.1\times 10^{-5})$

where $\mathbf{A}_{it}$ is a $p_{a}$-dimensional vector of covariates (i.e.,
intercept, exposure and measured confounders) and $\mathbf{a}_{i}$ is the
corresponding vector of coefficients. The vector $\bm{\mu}_{a}$ is of main
interest, specifically its element(s) corresponding to the overall exposure
effect.

The specification of $\mathbf{A}_{it}$ depends on the approach under
consideration as listed below.

  * •

Null: the model that does not account for unmeasured confounding, for which
$\mathbf{A}_{it}=(1,X_{it},{{}^{m}{\bf Z}}_{it}^{\prime})^{\prime}$. The
regression coefficients are allowed to vary across areas but not in time.

  * •

GLMadj: the model that includes smooth function of time, $h(t;df)$, modelled
as a natural cubic spline with $df$ degrees of freedom. Let $\mathbf{B}_{t}$
denote the vector of basis functions evaluated at time $t$. Then, we define
$\mathbf{A}_{it}=(1,X_{it},{{}^{m}{\bf
Z}}_{it}^{\prime},\mathbf{B}_{t}^{\prime})^{\prime}$.

  * •

Dummy: Dai et al. (2014), Klompmaker et al. (2023) assume that the long-term
exposure–outcome association vary by season. Following also Chen et al.
(2013), Peng et al. (2005), we may extend this approach to our setting by
including interaction terms. Let $D_{jt}$ indicate an indicator variable for
the $j$-th calendar season, for $j=1,\dots,4$, that is:

$D_{jt}=\begin{cases}1&\mbox{if the day $t$ is in the $j$-th season}\\\
0&\mbox{otherwise}\end{cases}$

Then, $\mathbf{A}_{it}=(1,X_{it}D_{1t},\dots,X_{it}D_{4t},{{}^{m}{\bf
Z}}_{it}^{\prime},\mathbf{B}_{t}^{\prime})^{\prime}$. Henceforth, the Dummy
model is similar to GLMadj, but it also captures season-specific exposure
effects.

  * •

Periodic: the exposure effect is assumed to be a linear combination of
trigonometric functions, as in Peng et al. (2005), Chen et al. (2013). In this
case,
$\mathbf{A}_{it}=(1,X_{it},X_{it}\cos(2\pi t/365),X_{it}\sin(2\pi
t/365),{{}^{m}{\bf Z}}_{it}^{\prime},\mathbf{B}_{t}^{\prime})^{\prime}$.

  * •

JDZ: the model proposed by Janes et al. (2007), where the local exposure
effect is of interest, hence
$\mathbf{A}_{it}=(1,X_{it}-\widehat{X_{t}},\widehat{X_{t}},{{}^{m}{\bf
Z}}_{it}^{\prime},\mathbf{B}_{t}^{*\prime})^{\prime}$. Here, $\widehat{X_{t}}$
denotes the overall temporal trend (see also Section 3), and
$\mathbf{B}_{t}^{*}$ denotes the vector of basis functions from a natural
cubic spline with $df-1$ degrees of freedom, evaluated at time $t$.

The unknown function $h(\cdot;\cdot)$ is modelled as a natural cubic spline of
time with 8 degrees of freedom per year (so $df=40$), and is used to control
for unmeasured confounding. The overall trend, $\widehat{X_{t}}$, is obtained
as the fitted values of a regression of $X_{it}$ on a natural cubic spline of
time with $4$ degrees of freedom per year. Bases from natural cubic splines of
temperature (with 7 degrees of freedom) and of relative humidity (with 4
degrees of freedom) are included as measured confounders. ${{}^{m}{\bf
Z}}_{it}$ also includes the indicator variables for the day of the week, for
holidays, and for the lockdown periods.

##  Appendix D Posterior Full Conditional Distributions for the Proposed DGLM

We provide the computational details to make inference on the model proposed
in Section 3 and following the discussion in Section 4 of the main manuscript.
In particular, stacking Equation (7) of the main text over the ${T}$ time
instants, we have (see also Chan and Jeliazkov, 2009):

$\bm{\vartheta}=\mathbf{o}+\mathbf{G}\bm{\gamma}+\tilde{\mathbf{u}}\,,$

where
$\bm{\vartheta}=(\bm{\vartheta}_{1}^{\prime},\dots,\bm{\vartheta}_{T}^{\prime})^{\prime}$,
$\mathbf{o}=(\mathbf{o}_{1}^{\prime},\dots,\mathbf{o}_{T}^{\prime})^{\prime}$,
$\tilde{\mathbf{u}}=(\tilde{\mathbf{u}}_{1}^{\prime},\dots,\tilde{\mathbf{u}}_{T}^{\prime})^{\prime}$,
$\bm{\gamma}=(\tilde{\bm{\beta}}_{1}^{\prime},\tilde{\bm{\beta}}_{2}^{\prime},\dots,\tilde{\bm{\beta}}_{{T}}^{\prime},\bm{\alpha}_{m}^{\prime})^{\prime}$,
and

$\mathbf{G}=\begin{bmatrix}\mathbf{\tilde{X}}_{\,\bm{\cdot}1}&&&{}^{m}\mathbf{Z}_{1}\\\
&\ddots&&\vdots\\\
&&\mathbf{\tilde{X}}_{\,\bm{\cdot}{T}}&{}^{m}\mathbf{Z}_{T}\end{bmatrix}\,.$

Posterior samples are obtained using a Gibbs sampler with the following full
conditionals:

$\displaystyle\bm{\gamma}|\bullet\sim\mathcal{N}(\overline{\mathbf{V}}\overline{\bm{\mu}},\overline{\mathbf{V}})$
$\displaystyle\sigma_{w1}^{2}|\bullet\sim\mathcal{IG}\left(a_{w1}+\frac{{T}-1}{2},b_{w1}+\frac{1}{2}\sum_{t=2}^{{T}}(\beta_{1t}-\beta_{1,t-1})^{2}\right)$
$\displaystyle\bm{\varSigma}_{w0}|\bullet\sim\mathcal{IW}\left(\mathbf{a}_{w0}+{T}-1,\mathbf{S}_{w0}+\sum_{t=2}^{{T}}(\tilde{\bm{\beta}}_{0t}-\tilde{\bm{\beta}}_{0,t-1})(\tilde{\bm{\beta}}_{0t}-\tilde{\bm{\beta}}_{0,t-1})^{\prime}\right)$
$\displaystyle\vartheta_{it}|\bullet\propto\exp\left\\{Y_{it}\vartheta_{it}-\exp(\vartheta_{it})-\frac{1}{2\tau_{u}^{2}}(\vartheta_{it}-o_{it}-\vartheta_{it}^{*})^{2}\right\\}\quad\mbox{for}\quad
i=1,\dots,{N}\quad\mbox{and}\quad t=1,\dots,{T}$
$\displaystyle\tilde{\tau}_{u}^{2}|\bullet\sim\mathcal{IG}\left(a_{u}+\frac{{N}{T}}{2},b_{u}+\frac{1}{2}\sum_{t=1}^{{T}}(\bm{\vartheta}_{t}-\mathbf{o}_{t}-\bm{\vartheta}_{t}^{*})^{\prime}(\bm{\vartheta}_{t}-\mathbf{o}_{t}-\bm{\vartheta}_{t}^{*})\right)$

where $\vartheta_{it}^{*}$ is a generic element of
$\bm{\vartheta}_{t}^{*}=(\vartheta_{11}^{*},\vartheta_{21}^{*},\dots,\vartheta_{{N}{T}}^{*})^{\prime}=\mathbf{G}\bm{\gamma}$,
and the symbol $\bullet$ denotes data and all the parameters except for the
parameter that is being updated. For $\vartheta_{it}$, $i=1,\dots,{N}$,
$t=1,\dots,{T}$, we use Metropolis-Hastings steps with a Gaussian proposal
distribution.

For the posterior full conditional of $\bm{\gamma}$, two more matrices need to
be defined:

$\mathbf{H}=\begin{bmatrix}\mathbf{I}_{N}&&&&&\\\
-\mathbf{I}_{N}&\mathbf{I}_{N}&&&&\\\ &-\mathbf{I}_{N}&\mathbf{I}_{N}&&&\\\
&&\ddots&\ddots&&\\\ &&&-\mathbf{I}_{N}&\mathbf{I}_{N}&\\\
&&&&&\mathbf{I}_{p_{m}}\end{bmatrix}\quad\mbox{and}\quad\mathbf{S}=\begin{bmatrix}\mathbf{V}_{\beta}&&&&&\\\
&\bm{\varSigma}_{w}&&&&\\\ &&\ddots&&&\\\ &&&\bm{\varSigma}_{w}&&\\\
&&&&\bm{\varSigma}_{w}&\\\ &&&&&10^{6}\mathbf{I}_{{p_{m}}}\end{bmatrix}\,.$

Hence, this posterior full conditional is multivariate normal,
$\bm{\gamma}|\bullet\sim\mathcal{N}(\overline{\mathbf{V}}\overline{\bm{\mu}},\overline{\mathbf{V}})$,
with

$\overline{\mathbf{V}}=(\mathbf{K}+\mathbf{G}^{\prime}\mathbf{G}/\tilde{\tau}_{u}^{2})^{-1}\quad\mbox{and}\quad\overline{\bm{\mu}}=\mathbf{G}^{\prime}\bm{\vartheta}/\tilde{\tau}_{u}^{2}\,,$

where $\mathbf{K}=\mathbf{H}^{\prime}\mathbf{S}^{-1}\mathbf{H}$. We refer to
Chan and Jeliazkov (2009) for further computational details.

##  Appendix E Additional Details on Simulation Study 1

###  E.1 Data Generation Mechanism

To simulate the variables, we use the following data generation mechanism. The
spatial domain comprises ${N}=100$ contiguous HDs, while the temporal
dimension is set to ${T}=300$. Let ${\bf X}=({\bf X}_{1},\dots,{\bf
X}_{{T}})^{\prime}$, with ${\bf X}_{t}=(X_{1t},\dots,X_{{N}t})^{\prime}$,
$t=1,\dots,{T}$, be the exposure, and let ${\bf Z}=({\bf Z}_{1},\dots,{\bf
Z}_{{T}})^{\prime}$, with ${\bf Z}_{t}=(Z_{1t},\dots,Z_{{N}t})^{\prime}$,
$t=1,\dots,{T}$, be the unmeasured confounder. It is thus assumed that

$\begin{bmatrix}\mathbf{X}\\\
\mathbf{Z}\end{bmatrix}\sim\mathcal{N}\left(\begin{bmatrix}\bm{\mu}_{x}\\\
\bm{\mu}_{z}\end{bmatrix},\begin{bmatrix}\bm{\varSigma}_{x}&\bm{\varSigma}_{xz}\\\
\bm{\varSigma}_{zx}&\bm{\varSigma}_{z}\end{bmatrix}\right)\,,$ (9)

where $\bm{\mu}_{x}$ and $\bm{\mu}_{z}$ are the spatio-temporal means
generated from Fourier basis expansions such that their correlation is
approximately $0.5$. It is assumed that the covariance matrices can be
factorized as
$\bm{\varSigma}_{x}=\bm{\varSigma}_{x}^{1/2}\bm{\varSigma}_{x}^{1/2\prime}$,
$\bm{\varSigma}_{z}=\bm{\varSigma}_{z}^{1/2}\bm{\varSigma}_{z}^{1/2\prime}$,
$\bm{\varSigma}_{xz}=\bm{\varSigma}_{zx}^{\prime}=\rho_{xz}\bm{\varSigma}_{x}^{1/2}\bm{\varSigma}_{z}^{1/2\prime}$,
and $\rho_{xz}={\rm\mathbb{C}or}(X_{it},Z_{it})$, for $i=1,\dots,{N}$, and
$t=1,\dots,{T}$.

We then assume that exposure and unmeasured confounder are spatial VAR(1)
processes (Cressie and Wikle, 2015, Sahu, 2022). We now explain how
$\bm{\varSigma}_{x}$ is parameterized, and the same applies to
$\bm{\varSigma}_{z}$. If a stable VAR(1) model is assumed for the exposure
where the temporal autocorrelation is controlled by a scalar,
$\phi_{x}^{T}\in(0,1)$, then it can be shown (see, for example, Lütkepohl,
2005, Section 2.1.4) that $\bm{\varSigma}_{x}=\bm{\Gamma}_{x,ar}\otimes{\bf
C}_{x}(0)$, where

${\bf C}_{x}(0)={\rm\mathbb{E}}\left[{\bf X}_{t}{\bf
X}_{t}^{\prime}\right]=\frac{\tau_{x}^{2}}{1-(\phi_{x}^{T})^{2}}\bm{\Omega}_{x,pcar}^{-1}$

indicates a zero-lag covariance matrix, and $\bm{\Gamma}_{x,ar}$ is the
correlation matrix of a first-order autoregressive process with generic
element
$\left(\bm{\Gamma}_{x,ar}\right)_{jk}=\left({\phi_{x}^{T}}\right)^{|j-k|}$.
The spatial dependence is defined by the matrix $\bm{\Omega}_{x,pcar}$, which
has a proper conditional autoregressive (PCAR) specification (Banerjee et al.,
2014, Besag, 1974, Sahu, 2022):

$\bm{\Omega}_{x,pcar}={\bf D}_{w}-\phi_{x}^{S}{\bf W}\,,$ (10)

where ${\bf D}_{w}$ is a diagonal matrix with the number of neighbors of each
HD on the diagonal, ${\bf W}$ is the binary adjacency matrix, and
$\phi_{x}^{S}\in(0,1)$ is the spatial autocorrelation parameter. ${\bf W}$ is
an ${N}\times{N}$ matrix with entries $w_{ij}=1$ if areal units $i$ and $j$
share a border, and $w_{ij}=0$ otherwise; also, $w_{ii}=0$.

We set $\rho_{xz}=0.5$, $\tau_{x}^{2}=4\times 10^{-4}$, and
$\tau_{z}^{2}=1\times 10^{-3}$, and then use different values for the spatial
and temporal autocorrelation parameters to define the scenarios in Table 1 of
the main manuscript.

Once all the matrices in Equation (9) are defined, we simulate the exposure
from its marginal distribution, then the unmeasured confounder from the
conditional distribution ${\bf Z}|{\bf
X}\sim\mathcal{N}\left(\bm{\mu}_{z|x},\bm{\varSigma}_{z|x}\right)$ where:

$\displaystyle\bm{\mu}_{z|x}$ $\displaystyle={\rm\mathbb{E}}\left[{\bf
Z}|\mathbf{X=x}\right]=\bm{\mu}_{z}+\bm{\varSigma}_{zx}\bm{\varSigma}_{x}^{-1}(\mathbf{x}-\bm{\mu}_{x}),$
$\displaystyle\bm{\varSigma}_{z|x}$
$\displaystyle={\rm\mathbb{V}ar}\left[\mathbf{Z}|\mathbf{X=x}\right]=\bm{\varSigma}_{z}-\bm{\varSigma}_{zx}\bm{\varSigma}_{x}^{-1}\bm{\varSigma}_{xz}\,.$

Finally, the outcome is generated using the following model:

$\displaystyle Y_{it}|\mu_{it}$ $\displaystyle\stackrel{{\scriptstyle
ind}}{{\sim}}\mathcal{P}oi(\mu_{it})$ (11) $\displaystyle\vartheta_{it}$
$\displaystyle=o_{it}+\beta_{0i}+\beta_{1t}X_{it}+Z_{it}+u_{it},\quad
u_{it}\stackrel{{\scriptstyle iid}}{{\sim}}\mathcal{N}(0,\tau_{u}^{2})\,,$
(12)

where $\mu_{it}=\exp(\vartheta_{it})$, $o_{it}=\log(E_{it})$ is the offset
term, and $\tau_{u}^{2}=10^{-5}$. The expected counts, $E_{it}$, are obtained
by using the resident population aged $\geq 65$ years and the age-specific
standardized mortality rates of year 2018. The population and the rates are
assumed to be constant in time.

###  E.2 Performance Assessment

To complete the results of Simulation Study 1, Table E.4 reports, for each
scenario and model, the value for the three indices, namely:

  * •

average mean squared error,
$MSE_{avg}=\frac{1}{{T}}\sum_{t=1}^{{T}}\left[\frac{1}{100}\sum_{i=1}^{100}\left({\rm\mathbb{E}}[P_{t}|\bullet]-P_{t}\right)^{2}\right]$

  * •

maximum absolute bias, $MAB=\max_{t}|{\rm\mathbb{E}}[P_{t}|\bullet]-P_{t}|$

  * •

average mean absolute error,
$MAE_{avg}=\frac{1}{{T}}\sum_{t=1}^{{T}}\left[\frac{1}{100}\sum_{i=1}^{100}|{\rm\mathbb{E}}[P_{t}|\bullet]-P_{t}|\right]$

Here, the differences are computed with respect to $P_{t}$, i.e. the true
simulated percent change of risk.

Table E.4 indicates that, if the exposure effect is constant, the GLMadj model
should be preferred, otherwise the DGLM is the best choice. Note, however,
that the shape of the exposure effect is not known in practice, so the DGLM is
the most robust choice across all scenarios.

Table E.4: Simulation Study 1: $MSE_{avg}$, $MAB$, and $MAE_{avg}$ calculated for each model considered under each Scenario. The lowest values are indicated in italic. |  | $MSE_{avg}$ | $MAB$ | $MAE_{avg}$
---|---|---|---|---
Scenario | Exposure Effect | DGLM | Null | GLMadj | DGLM | Null | GLMadj | DGLM | Null | GLMadj
S1 | Constant | 0.642 | 63.158 | 0.368 | 0.728 | 7.929 | 0.568 | 0.686 | 7.929 | 0.568
Periodic | 0.697 | 115.672 | 3.276 | 0.833 | 12.678 | 3.109 | 0.717 | 10.643 | 1.490
Complex | 0.672 | 66.741 | 1.283 | 1.217 | 9.912 | 2.423 | 0.691 | 8.099 | 0.989
S2 | Constant | 6.886 | 19.741 | 6.622 | 3.660 | 4.435 | 2.547 | 2.571 | 4.435 | 2.547
Periodic | 7.177 | 32.742 | 10.620 | 3.337 | 7.566 | 4.935 | 2.625 | 5.531 | 2.900
Complex | 6.782 | 20.939 | 8.094 | 3.691 | 6.285 | 4.476 | 2.547 | 4.472 | 2.664
S3 | Constant | 0.251 | 45.706 | 0.085 | 0.406 | 6.760 | 0.260 | 0.331 | 6.760 | 0.260
Periodic | 0.276 | 86.516 | 2.706 | 0.480 | 11.224 | 2.823 | 0.360 | 9.189 | 1.396
Complex | 0.271 | 48.538 | 0.981 | 0.911 | 8.717 | 2.126 | 0.325 | 6.904 | 0.852

###  E.3 On The Removal of Large-Scale Temporal Trend

We argue that the term $X_{it}-\widehat{X_{t}}$, for $i=1,\dots,{N}$ and
$t=1,\dots,{T}$, should be used in place of simply $X_{it}$ in the DGLM.
Removing the overall temporal trend from the exposure not only follows from
theoretical arguments (see Equation (3) in the main paper), but also it is
crucial to control for confounding bias. The local variations in the exposure
are in fact less likely to be confounded by the unmeasured confounder, which
is assumed to be spatially and temporally smooth (Janes et al., 2007, Greven
et al., 2011) – see also Section 3 of the manuscript. To illustrate this with
simulated data, we compare the proposed DGLM with a reparameterized model that
does not remove the overall trend from the exposure. Essentially, Equation (6)
in the paper can be rewritten as:

$\vartheta_{it}=o_{it}+\beta_{0it}^{*}+\beta_{1t}X_{it}+{{}^{m}\mathbf{Z}}_{it}^{\prime}\bm{\alpha}_{m}+\tilde{u}_{it}\,,$

where

$\beta_{0it}^{*}=\beta_{0it}-\beta_{1t}\widehat{X_{t}}\,.$

We refer to this model as DGLM_single, since it is a single-step approach that
does not require to estimate the overall temporal trend.

The results, in terms of percent change of risk associated with a 10-unit
increase in the exposure, $P_{t}$, for each scenario considered in the main
paper are depicted in Figure E.1. The estimates tend to differ at the
beginning and at the end of the time window, where the DGLM_single model is
more biased. In addition, note that the DGLM_single model is not able to
recover the constant effect in the top-left subplot. For completeness, Table
E.5 reports the $MSE_{avg}$, $MAB$, and $MAE_{avg}$ for both model considered
under each scenario, where it is confirmed that the proposed DGLM outperforms
the DGLM_single model in all cases.

Figure E.1: Simulation Study 1: results are expressed in terms of percent change of risk associated with a 10-unit increase in the exposure, $P_{t}$. The black line represents the true exposure effect (whose shape differs in each row). The solid red and dashed brown lines represent the average of the posterior means of $P_{t}$ obtained from DGLM and DGLM_no over 100 replicates, respectively. The shaded areas represent the percentiles $2.5$ and $97.5$ of the posterior means of $P_{t}$ under each approach. Table E.5: Simulation Study 1: $MSE_{avg}$, $MAB$, and $MAE_{avg}$ calculated for the DGLM and DGLM_single models under each Scenario. The lowest values for each index are indicated in italic. | Exposure | $MSE_{avg}$ | $MAB$ | $MAE_{avg}$
---|---|---|---|---
Scenario | Effect | DGLM | DGLM_single | DGLM | DGLM_single | DGLM | DGLM_single
S1 | Constant | 0.642 | 0.873 | 0.728 | 2.002 | 0.686 | 0.773
Periodic | 0.697 | 0.816 | 0.833 | 1.762 | 0.717 | 0.763
Complex | 0.672 | 0.890 | 1.217 | 1.955 | 0.691 | 0.799
S2 | Constant | 6.886 | 7.711 | 3.660 | 4.351 | 2.571 | 2.707
Periodic | 7.177 | 7.695 | 3.337 | 4.178 | 2.625 | 2.706
Complex | 6.782 | 7.454 | 3.691 | 4.333 | 2.547 | 2.660
S3 | Constant | 0.251 | 0.254 | 0.406 | 1.033 | 0.406 | 0.418
Periodic | 0.276 | 0.250 | 0.480 | 0.957 | 0.428 | 0.417
Complex | 0.271 | 0.282 | 0.911 | 1.183 | 0.418 | 0.436

##  Appendix F Simulation Study 2

We now illustrate Simulation Study 2, where we want to assess whether the
proposed DGLM is able to recover the first derivative of the unknown
functional form, which represents the marginal effect of the exposure on the
outcome. To this end, we consider Scenarios S1 and S2 again, but we assume
that the linear predictor is as follows:

$\vartheta_{it}=o_{it}+\beta_{0i}+f(X_{it};\bm{\gamma}_{f})+Z_{it}+u_{it}\,,$

where
$f(X_{it};\bm{\gamma}_{f})=\gamma_{f0}+\gamma_{f1}X_{it}+\gamma_{f2}X_{it}^{2}$,
$\gamma_{f0}=\gamma_{f1}=0$, $\gamma_{f2}=2.5\times 10^{-5}$. Henceforth, the
true values of the time-varying coefficients are $\beta_{1t}=5\times
10^{-5}\widehat{X_{t}}$, and $\beta_{0it}=\beta_{0i}+2.5\times
10^{-5}\widehat{X_{t}}^{2}+Z_{it}$. All the other parameters are the same as
in the previous Section.

#### Scenario S1.

Figure F.2 refers to Scenario S1 and consists of two subplots. The top subplot
compares spatial averages of the time-varying intercepts. That is,

$\overline{\beta_{0t}}=\frac{1}{{N}}\sum_{i=1}^{{N}}\beta_{0it}\quad\mbox{for}\quad
t=1,\dots,{T}$

and

$\overline{{\rm\mathbb{E}}[\beta_{0t}|\bullet]}=\frac{1}{{N}}\sum_{i=1}^{{N}}{\rm\mathbb{E}}[\beta_{0it}|\bullet]\quad\mbox{for}\quad
t=1,\dots,{T}$

denote the average true intercept and average posterior mean, respectively.
These quantities are different across the 100 replicates, so the blue (green)
line corresponds to the average $\overline{\beta_{0t}}$
($\overline{{\rm\mathbb{E}}[\beta_{0t}|\bullet]}$), while the shaded area of
the same color defines the percentiles $2.5$ and $97.5$. The blue and green
shaded areas do not completely overlap, hence the unmeasured confounder can be
captured only in part. The bottom subplot shows that the true value of
$\beta_{1t}$ (black line) is almost always within the red shaded area,
indicating that confounding bias is alleviated but not completely ruled out.

Figure F.2: Simulation Study 2: results from Scenario S1 with
$f(X_{it};\bm{\gamma}_{f})=\gamma_{f0}+\gamma_{f1}X_{it}+\gamma_{f2}X_{it}^{2}$.
Top: the blue solid line and the blue shaded area are the average and the
percentiles $2.5$ and $97.5$ of the spatial mean of the simulated intercepts
over 100 replicated datasets; the green solid line and green shaded area
depict the same quantities but for the corresponding spatially-averaged
posterior mean obtained from the proposed DGLM. Bottom: the black solid line
represents the true exposure effect, while the red solid line and the shaded
area represent the average and the percentiles $2.5$ and $97.5$ of the
posterior means of $\beta_{1t}$.

#### Scenario S2.

As for Scenario S2, Figure F.3 indicates that it is not possible to completely
recover the true function, $f(\cdot;\cdot)$. The unmeasured confounder has low
temporal autocorrelation, and the posterior means of the intercepts tend to
have a rougher behavior. However, some variation in the unmeasured confounder
is “absorbed” by the local exposure effect, which remains biased.

Figure F.3: Simulation Study 2: results from Scenario S2 with
$f(X_{it};\bm{\gamma}_{f})=\gamma_{f0}+\gamma_{f1}X_{it}+\gamma_{f2}X_{it}^{2}$.
Top: the blue solid line and the blue shaded area are the average and the
percentiles $2.5$ and $97.5$ of the spatial mean of the simulated intercepts
over 100 replicated datasets; the green solid line and green shaded area
depict the same quantities but for the corresponding spatially-averaged
posterior mean obtained from the proposed DGLM. Bottom: the black solid line
represents the true exposure effect, while the red solid line and the shaded
area represent the average and the percentiles $2.5$ and $97.5$ of the
posterior means of $\beta_{1t}$.

#### Scenario without confounding.

To understand whether the residual bias in Figures F.2 and F.3 is due to the
Taylor’s approximation, we now assume that the outcome is not affected by any
confounder, i.e. ${\bf Z}={\bf 0}$. We only report the results for the
Scenario S1 where $\phi_{x}^{S}=\phi_{x}^{T}=0.2$. The estimated exposure
effect is depicted in Figure F.4, which indicates that the proposed DGLM is
able to recover the true effect, with no visible difference between the black
and red lines for almost all the temporal domain.

Figure F.4: Simulation Study 2: results from Scenario S1 without confounding.
The black solid line represents the true exposure effect. The red solid line
and the shaded area represent the average and the percentiles $2.5$ and $97.5$
of the posterior means of $\beta_{1t}$ obtained from the proposed DGLM over
100 replicated datasets.

##  Appendix G Application: Comparison with Restricted Model

We assess the impact that imposing some constraints on the proposed DGLM has
on the estimated exposure effect. Specifically, we remove the measured
confounders from the model specification. This is equivalent to imposing that
$\bm{\alpha}_{m}={\bf 0}$ in Equation (6) of the main text. Hence, the linear
predictor and the temporal dynamics of the coefficients are specified as
follows:

$\displaystyle\bm{\vartheta}_{t}$ $\displaystyle={\bf
o}_{t}+\mathbf{\tilde{X}}_{t}\tilde{\bm{\beta}}_{t}+\tilde{{\bf u}}_{t}\,,$
$\displaystyle\tilde{\bm{\beta}}_{t}$
$\displaystyle=\tilde{\bm{\beta}}_{t-1}+\mathbf{w}_{t}\,.$

We refer to the implied model as DGLM-X to emphasize that the exposure is the
only covariate. Figure G.5 shows the posterior mean of the percent change of
risk, as defined above, under the DGLM-X (in green) and the proposed DGLM (in
red). It suggests that, in general, the temporal dynamic of the exposure
effect would not be affected much if the measured confounders were not
included in the model. Henceforth, the proposed DGLM can provide accurate
estimates even in complex scenarios.

However, it is worth noting that the DGLM-X provides greater estimates than
the proposed estimates during the summer, while during the winter the
estimates are quite similar. This suggests that the confounding effect may
vary across time. We hypothesize that this is due to air temperature that
exceeds a certain threshold; in other words, confounding may occur due to heat
waves. In fact, heat waves are demonstrated to have an increased effect on
mortality during summer days with relatively high levels of PM2.5
concentrations, and this positive correlation may be the source of the
overestimation (Analitis et al., 2008, 2014). Another study observed a similar
result for the population aged 75 years or more (Analitis et al., 2018).
Another possible explanation for the overestimation could be the modification
of the PM2.5 chemical composition, with high levels of the most toxic
compounds developing in summer (Vinson et al., 2020).

Figure G.5: Estimated time-varying percent change of risk associated with a
$10\ \mu g/m^{3}$ increase in the PM2.5 concentrations. The red line
represents the posterior mean of the proposed DGLM, while the green line the
posterior mean of the DGLM-X.