# Bayesian Generalized Distributed Lag Regression with Variable Selection

Daniel Dempsey and Jason Wyse  School of Computer Science and Statistics,
Trinity College Dublin

(March 2024)

###### Abstract

Distributed Lag Models (DLMs) and similar regression approaches such as MIDAS
have been used for many decades in econometrics, and more recently in the
study of air quality and its impact on human health. They are useful not only
for quantifying accumulating and delayed effects, but also for estimating the
lags that are most susceptible to these effects. Among other things, they have
been used to infer the period of exposure to poor air quality which might
negatively impact child birth weight. The increased attention DLMs have
received in recent years is reflective of their potential to help us
understand a great many issues, particularly in the investigation of how the
environment affects human health. In this paper we describe how to expand the
utility of these models for Bayesian inference by leveraging latent-variables.
In particular we explain how to perform binary regression to better handle
imbalanced data, how to incorporate negative binomial regression, and how to
estimate the probability of predictor inclusion. Extra parameters introduced
through the DLM framework may require calibration for the MCMC algorithm, but
this will not be the case in DLM-based analyses often seen in pollution
exposure literature. In these cases, the parameters are inferred through a
fully automatic Gibbs sampling procedure.

_Keywords:_ Distributed lag regression, Polya-Gamma regression, Functional
data analysis, Air quality, Pollution exposure, Environmental modelling.

##  1 Introduction

It can often be the case that variables are interconnected or correlated in a
manner that accumulates over time. For example, the extent to which exposure
to poor air quality affects human health will depend on both the concentration
of poor quality air and the length of time an individual was exposed to it.
There are likely to be temporal windows of exposure that will more directly
impact current health status than others. A basic assumption along these lines
is that the effect of exposure will asymptotically (in time) decline, though
of course there is no guarantee that the decline is monotonic, at least in the
short term. Understanding these temporal ‘profiles’ of exposure is key to
understanding the effects of these exposures as a whole. _Distributed Lag
Models_ (DLMs) can help us understand these exposure profiles by applying
weights to the temporal lags of the covariates corresponding to the degree to
which they contribute the predictor’s effect on the response variable.

Distributed lags were introduced by Almon, (1965) and were most strongly
associated with econometrics for many of the intervening decades. See for
example Hannan, (1965), Chen et al., (1972), Sims, (1971) and Dhrymes, (1981)
for some early examples. A related class of models known as _Mixed Data
Sampling_ (MIDAS) models were introduced by Ghysels et al., (2004) to resolve
the issue of analysing data of mixed sampling regimes (for example, measuring
the association of a time series with observations sampled yearly and another
with observations sampled monthly). This too was mainly in service of
econometrics. In recent years, DLMs have been used to investigate the impact
of the environment on different aspects of human health. This appears to have
been popularised by Schwartz, (2000) who implemented a DLM into a generalised
additive model (Hastie and Tibshirani,, 1987) to analyse how daily death count
is correlated with exposure to poor air quality. A similar study was performed
later by Zanobetti et al., (2002).

The utility of DLMs (and MIDAS) is evident from the number of analyses
centered around them; Wilson et al., (2017) used a Bayesian DLM to investigate
the relationship between air quality and adverse health outcomes of newborns.
Warren et al., (2020) implemented random effects to account for spatial
autocorrelation into their DLM when they investigated the effect of air
quality exposure during pregnancy and the birth weight of newborns. Clements
and Galvão, (2008) suggest how to safely incorporate auto-regression terms
into MIDAS models. Xu et al., (2019) synthesised the MIDAS model with an
artificial neural network model for very flexible non-linear fits. Li et al.,
(2021) took this a step further by incorporating the effects of temporal
correlation. Antonelli et al., (2021) developed a Bayesian DLM that utilises
spike-and-slab priors for variable selection (Mitchell and Beauchamp,, 1988).
Mogliani and Simoni, (2020) similarly suggests a MIDAS model with spike-and-
slab priors, but with a Laplacian distribution for the slab to incorporate
Bayesian Group-LASSO variable selection (Meier et al.,, 2008; Xu and Ghosh,,
2015). Mork et al., (2021) combines DLMs with Bayesian additive regression
trees (Chipman et al.,, 2010) to analyse the relationship between air quality
exposure and the weight of newly born children.

In this article we demonstrate how to perform Bayesian regression for
generalised DLF models with particular focus paid to the negative binomial and
binary quantile models. Depending on the exact specifications of the DLM (to
be discussed later), Markov Chain Monte Carlo (MCMC) implementation of these
models are fully automatic beyond setting priors and starting values. As has
become common for Bayesian GLMs since Albert and Chib, (1993), these methods
are facilitated by the incorporation of latent variables for their ease of
implementation. The latent variable approaches we adopt can be leveraged to
infer the posterior probability of particular exposures having an impact on
health or other indicators.

The remainder of the article is laid out as follows. In section 2 we present
the general framework for generalised DLMs, focusing on how we can apply
binary quantile regression and negative binomial regression in particular. We
outline a simulation study in section 3 for assessing the accuracy of the
methods at estimating associated parameters and how well the variable
selection procedure performs. There is an extra document supplementary to this
manuscript that contains further simulation study results (via visualisations)
than what is presented here. In section 4 the proposed approach is
demonstrated on a real data example. Both the real and simulated data examples
are based on Chicago air quality data available in the dlnm package for R
(Gasparrini,, 2011). We conclude in section 5 with a description of merits and
possible extensions of the work.

##  2 Methods

Denote the vector response variables as $\mathbf{y}$. We want to model its
association with a collection of predictors, some of which have a ‘standard’
_contemporaneous_ effect on the response (referred to as _static_ predictors),
but others that we believe have a lagged effect on the response which is
modulated by the dosage or duration of exposure to the longitudinal predictor
(referred to as _dynamic_ predictors). For example, measuring how the exposure
to bad air quality affects respiratory function. We need to measure the effect
of bad air quality as in standard regression problems, but now we also need to
infer how that effect is distributed over its longitudinal units, usually
time. The relationship between the lags of the dynamic predictors and the
response will be referred to as the _lag-response_ model. DLMs are regression
methods built specifically to handle these kinds of problems.

We will denote the matrix of dynamic predictors as $\mathbf{Z}$ and the static
predictors (including the intercept, if used) as $\mathbf{Z}_{\text{S}}$.
Under this notation, DLMs can be expressed in the following general form;

$\displaystyle\mathbf{y}\sim
F\left(g^{-1}(\mbox{\boldmath{$\eta$}}),\mbox{\boldmath{$\Phi$}}\right)$ (1)
$\displaystyle\mbox{\boldmath{$\eta$}}=\mathbf{Z}_{\text{S}}\mbox{\boldmath{$\beta$}}_{\text{S}}+\sum_{i=1}^{p}h_{i}\left(\mathbf{z}_{i},\mbox{\boldmath{$\Theta$}}_{i}\right)$
(2)

$F$ is the assumed distribution of the response $\mathbf{y}$, $g^{-1}(\cdot)$
is the inverse link function that allows for generalised regression, and
$\Phi$ is the collection of parameters of $F$ that are required alongside the
transformed linear predictor. $\mathbf{z}_{i}$ is the $i^{\text{th}}$ dynamic
variable, or equivalently, the $i^{\text{th}}$ column of the $\mathbf{Z}$
matrix, $i=1,\cdots,p$. $\mbox{\boldmath{$\beta$}}_{S}$ is the vector of
linear regression parameters corresponding to the static variables. Notice
that if there are no dynamic predictors, this simply collapses to a standard
GLM.

$h_{i}(\cdot)$ encodes the assumed nature of the lag-response of the
$i^{\text{th}}$ dynamic predictor, which will usually be further parameterised
by $\mbox{\boldmath{$\Theta$}}_{i}$. This can be as simple as a crude
aggregate over the lags (in which case there will be no
$\mbox{\boldmath{$\Theta$}}_{i}$). Another option is to use a linear
combination of the last $\tau$ lags (referred to as the _time window_) of the
dynamic variables, effectively treating each lag as its own predictor. This is
a conceptually simple way to model the lag-response but it is heavily
parameterised with a large degree of co-linearity across the lag predictors. A
common compromise in longitudinal health studies is to use a linear
combination of a basis expansion of the lags (Gasparrini,, 2011). Another
option, commonly seen in the MIDAS literature, is to assign weights to the
lags based on a (usually) parameterised curve. The weights are incorporated
into the model through convolution, scaled by a linear coefficient;

$h_{i}\left(\mathbf{z}_{i},\mbox{\boldmath{$\Theta$}}_{i}\right)=\beta_{i}\sum_{k=0}^{\tau_{i}}w_{i}(k;\mbox{\boldmath{$\theta$}}_{i})Z_{t-k,i}$

where $w_{i}(\cdot)$ is the function that assigns the weights, called a
_distributed lag function_ (DLF) with parameters
$\mbox{\boldmath{$\theta$}}_{i}$, and $Z_{i,j}$ denotes the $i,j^{\text{th}}$
cell of matrix $\mathbf{Z}$.

The DLF can be any function that asymptotically declines to zero as the lags
grow larger, but it must be constrained in some way to ensure the coefficient
$\beta_{i}$ is identifiable. A common constraint is
$w_{i}(k,\mbox{\boldmath{$\theta$}}_{i})>0$ for all $k$ and
$\sum_{k=1}^{\tau_{i}}w_{i}(k,\mbox{\boldmath{$\theta$}}_{i})=1$. This
particular approach is useful as it ensures that a dynamic variable can’t swap
between having positive _and_ negative impacts on the value of the response
(insofar that it has any impact at all) for different lags. A list of typical
choices for DLFs are given in Table 1 of Ghysels et al., (2016).

While the DLF approach offers what appears to be a pleasing compromise between
model flexibility and analyst control, the non-linear relationship between the
values of $\mathbf{y}$ and $\theta$ can present challenges for inference and
therefore requires special consideration in MCMC routines, such as extra
Metropolis-within-Gibbs steps. By contrast, the linear combinations of the
lags or basis expansions can easily be incorporated, as they require no
further alterations of the methodologies we will discuss in the forthcoming
sections.

###  2.1 Model Inference

Going back to (1), in this manuscript we are primarily concerned with the
situation where $F$ is a discrete distribution, used to model either count or
binary responses. When estimating the posterior, straightforward application
of Bayes’ Theorem in either case leads to full conditionals that cannot be
sampled directly regardless of the choice of prior. These are recognised as
challenging modelling problems for Bayesian methods. However, in both cases we
can alleviate this issue with the introduction of latent variables. First let
us introduce some extra notation; let $\beta$ (with the subscript omitted),
denote the vector containing $\mbox{\boldmath{$\beta$}}_{S}$ and the subset of
all $\mbox{\boldmath{$\Theta$}}_{i}$ that are linear coefficients. Similarly
let $\mathbf{X}$ be the column-wise concatenation of $\mathbf{Z}_{\text{S}}$
and the vectors relating to the linear predictors of
$\mbox{\boldmath{$\Theta$}}_{i}$. Note that, when using a linear combination
of basis expansions as $h_{i}(\cdot)$ then $\beta$ contains all regression
parameters and $\mathbf{X}$ contains all predictors.

####  2.1.1 Count Data Inference

Automatic Gibbs sampling for the parameters of the negative binomial
regression model was explained by Pillow and Scott, (2012) and Zhou et al.,
(2012), which is briefly summarised here. We start by assuming $y_{t}$ (in
this case a discrete count) follows a negative binomial distribution with
probability parameter equal to the transformed linear predictor and unknown
stopping parameter $\xi$,

$\displaystyle y_{t}$ $\displaystyle\sim$ $\displaystyle\text{NB}(\xi,p_{t})$
$\displaystyle p_{t}$ $\displaystyle=$
$\displaystyle\frac{\exp(\eta_{t})}{1+\exp(\eta_{t})}$

where $\eta_{t}$ are the components of $\eta$ given in (2). Assuming a
Gaussian prior of mean vector $\mathbf{m}$ and covariance matrix $\mathbf{v}$,
exact and automatic Gibbs sampling from the posterior of $\beta$ can be
accomplished by introducing a latent variable $\omega$ that follows a Poyla-
Gamma distribution (Polson et al.,, 2013),

$\displaystyle\omega_{t}\,|\,\mbox{\boldmath{$\beta$}},\xi$
$\displaystyle\sim$
$\displaystyle\text{PG}\left(y_{t}+\xi,\,\mathbf{z}^{\top}_{t,\text{S}}\mbox{\boldmath{$\beta$}}_{\text{S}}+\sum_{i=1}^{p}h_{i}\left(\mathbf{Z}_{t,i},\mbox{\boldmath{$\Theta$}}_{i}\right)\right)$
$\displaystyle\mbox{\boldmath{$\beta$}}\,|\,\mbox{\boldmath{$\omega$}},\xi$
$\displaystyle\sim$
$\displaystyle\mbox{N}\left(\mathbf{M}_{C},\mathbf{V}_{C}\right)$
$\displaystyle\mathbf{V}_{C}$ $\displaystyle=$
$\displaystyle\left(\mathbf{X}^{\top}\mbox{\boldmath{$\Omega$}}_{C}\mathbf{X}+\mathbf{v}^{-1}\right)^{-1}$
$\displaystyle\mathbf{M}_{C}$ $\displaystyle=$
$\displaystyle\mathbf{V}_{C}\left(\mathbf{X}^{\top}\mbox{\boldmath{$\Omega$}}_{C}\mbox{\boldmath{$\lambda$}}_{C}+\mathbf{v}^{-1}\mathbf{m}\right)$
$\displaystyle\mbox{\boldmath{$\Omega$}}_{C}$ $\displaystyle=$
$\displaystyle\text{diag}(\mbox{\boldmath{$\omega$}})$
$\displaystyle\lambda_{t,C}$ $\displaystyle=$
$\displaystyle(y_{t}-\xi)/2\omega_{t}$

where diag($\cdot$) constructs a diagonal matrix based on the components of
the vector input. The subscript $C$ used throughout to highlight that these
are the parameters corresponding to the count data model; this is to avoid
notational conflict in the next section. Regarding the Polya-Gamma
distribution sampling, section 4 of Polson et al., (2013) describes an exact
and efficient rejection sampler which they have implemented in their R package
BayesLogit.

To sample from the posterior of $\xi$, for conjugacy we assume a gamma
distribution prior with shape and scale parameter $a_{0}$ and $b_{0}$. We also
need to include yet another latent variable $\psi$ to the model. $\xi$ is then
updated as follows:

$\displaystyle\pi(\psi_{t}=j\,|\,\xi)=R_{y_{t},j}$
$\displaystyle\xi\sim\text{Gamma}\left(a_{0}+\sum_{i=1}^{N}\psi_{i},b_{0}-\sum_{i=1}^{N}\ln(1-p_{t})\right)$

where $N$ is the length of $\mathbf{y}$ and $R_{i,j}$ is the $i,j\text{th}$
element of the following (recursively constructed) lower triangular matrix:

$R_{i,j}=\frac{i-1}{i}R_{i-1,j}+\frac{\xi}{i}R_{i-1,j-1}$ (3)

with $R_{1,1}=R_{0,1}=1$, and $R_{i,j}=0$ when $j=0$ or $j>i$. This matrix has
square dimensions equal to the maximum value of $\mathbf{y}$. Construction of
this matrix is the biggest potential source of computational bottleneck as it
has to be re-done at every iteration, cell-by-cell. An alternative non-
recursive method exists (Zhou et al.,, 2012) but it becomes numerically
unstable as $y_{\text{max}}$ or $\xi$ grow too large, so in general we
recommend using (3).

Poisson regression can be approximated by simply fixing $\xi$ to some large
number but otherwise following the above approach, circumventing the need to
compute (3). See D’Angelo and Canale, (2022) for further notes on approximate
Bayesian Poisson regression. Zhou et al., (2012) augmented the above algorithm
to capture broader forms of variation, and also derived associated Variational
Bayes estimates.

####  2.1.2 Possibly Imbalanced Binary Data Inference

Now let us assume $y_{t}$ is binary. Exact Gibbs sampling for probit
regression (and an approximation for logistic regression) via latent variables
was first proposed by Albert and Chib, (1993) and has become the blueprint for
binary response regression; see Holmes and Held, (2006); Frühwirth-Schnatter
and Frühwirth, (2010); Polson et al., (2013); Zens et al., (2023) for more
examples of how it is applied for logistic regression. Here, we wish to use
this same approach as a way of handling DLMs for binary responses, where the
response may be imbalanced. This can occur in, say, medical studies where the
focus is on rare diseases, symptoms or other health outcomes.

The logit link has been shown to result in bias of the MLE of the regression
parameters when it misspecifies the ‘true’ link and that this bias is
particularly diminished by skewness (Czado and Santner,, 1992). For this
reason, the approach we recommend here is binary quantile regression (Benoit
and Van den Poel,, 2012; Benoit and den Poel,, 2017) using the 3-parameter
_Asymmetric Laplace Distribution_ (ALD) (Yu and Zhang,, 2005), which uses a
parameter $q\in\\{0,1\\}$ (corresponding to the desired quantile) that
controls the degree of skew. When $q=0.5$ it collapses to the standard
(symmetric) Laplacian distribution. The ALD Bayesian quantile model is

$\displaystyle y_{t}=I(y_{t}^{\star}>0)$
$\displaystyle\mathbf{y}^{\star}=\mathbf{Z}_{\text{S}}\mbox{\boldmath{$\beta$}}_{\text{S}}+\sum_{i=1}^{p}h_{i}\left(\mathbf{z}_{i},\mbox{\boldmath{$\Theta$}}_{i}\right)+\mbox{\boldmath{$\epsilon$}}$
$\displaystyle\epsilon_{t}\sim\text{ALD}(\mu=0,\,\sigma=1,\,q)$

where $I(\cdot)$ is equal to 1 if the condition inside the parenthesis is met
and 0 otherwise. $\mu$ and $\sigma$ are location and scale parameters for the
ALD. Kozumi and Kobayashi, (2011) found that the ALD can be expressed as
scale-mixture of Gaussians, allowing us to form an algorithm for exact Gibbs
sampling, though we exclude the details and derivations here. First we update
the $\mathbf{y}^{\star}$ latent variable,

$y^{\star}_{t}|y_{t},\mbox{\boldmath{$\beta$}},\mbox{\boldmath{$\theta$}}\sim\begin{cases}\text{ALD}(\mathbf{x}_{t}^{\top}\mbox{\boldmath{$\beta$}},\sigma=1,q)&\text{restricted
to the positive axis when $y_{i}=1$,}\\\
\text{ALD}(\mathbf{x}_{t}^{\top}\mbox{\boldmath{$\beta$}},\sigma=1,q)&\text{restricted
to the negative axis when $y_{i}=0$}.\end{cases}$

This truncated ALD is easily sampled from using the inversion method, or via
the approach proposed in section 4.2 of Benoit and den Poel, (2017). Here we
are jointly updating $\mathbf{y}$ with another latent variable $\nu$, whose
full conditional distribution is a Generalised Inverse Gaussian distribution
(Jorgensen,, 2012),

$\nu_{t}\,|\,y^{\star}_{t},\mbox{\boldmath{$\beta$}}\sim\text{GIG}(1/2,\chi_{t}^{2},\delta^{2})$

where

$\chi_{t}^{2}=\frac{(y^{\star}_{t}-\mathbf{x}_{t}^{\top}\mbox{\boldmath{$\beta$}})^{2}}{\phi^{2}}\,,\quad\delta^{2}=2+\frac{\psi^{2}}{\phi^{2}}\,,\quad\psi^{2}=\left(\frac{1-2q}{q(1-q)}\right)^{2}\,,\quad\phi^{2}=\frac{2}{q(1-q)}.$

The full conditional distribution of $\mathbf{y}$ (including $\nu$) is
truncated Gaussian, but the joint update with $\nu$ improves Monte Carlo
efficiency.

As with negative binomial inference, the resulting full conditional
distribution of the linear coefficients is Gaussian,

$\mbox{\boldmath{$\beta$}}\,|\,\mathbf{y},\mbox{\boldmath{$\theta$}},\mbox{\boldmath{$\nu$}}\sim\mbox{N}(\mathbf{M}_{B},\mathbf{V}_{B})$

where

$\displaystyle\mathbf{V}_{B}$ $\displaystyle=$
$\displaystyle\left(\mathbf{X}^{\top}\mbox{\boldmath{$\Omega$}}_{B}\mathbf{X}+\mathbf{v}^{-1}\right)^{-1}$
$\displaystyle\mathbf{M}_{B}$ $\displaystyle=$
$\displaystyle\mathbf{V}_{B}\left(\mathbf{X}^{\top}\mbox{\boldmath{$\Omega$}}_{B}\mbox{\boldmath{$\lambda$}}_{B}+\mathbf{v}^{-1}\mathbf{m}\right)$
$\displaystyle\mbox{\boldmath{$\Omega$}}_{B}$ $\displaystyle=$
$\displaystyle\left(\phi^{2}\text{diag}(\mbox{\boldmath{$\nu$}})\right)^{-1}$
$\displaystyle\mbox{\boldmath{$\lambda$}}_{B}$ $\displaystyle=$
$\displaystyle\mathbf{y}^{\star}-\psi\mbox{\boldmath{$\nu$}}$

and as before, $\mathbf{m}$ and $\mathbf{v}$ are the mean and covariance of
the prior Gaussian distribution for $\beta$.

In the event that $h(\cdot)$ from equation (1) is chosen in such a way that
all the $\mbox{\boldmath{$\Theta$}}_{i}$ are parameters for linear predictors,
then this is all that is necessary for fitting DLMs with count or binary
responses. As stated before, any non-linear parameters that arise from the
lag-response model (as usually seen when applying DLFs) will require further
consideration.

###  2.2 Covariate Inclusion Inference

Let $\gamma_{i}\in\\{0,1\\}$ represent whether or not a corresponding
predictor (or groups of predictors) be included in the model. Regardless of
which of the above models we use, the approach of inferring the posterior of
$\gamma$ is the same. Holmes and Held, (2006) suggest we update $\gamma$
jointly with $\beta$ through a Metropolis update,

$q(\mbox{\boldmath{$\beta$}}^{\star},\mbox{\boldmath{$\gamma$}}^{\star})=\pi(\mbox{\boldmath{$\beta$}}^{\star}\,|\,\mbox{\boldmath{$\gamma$}}^{\star},\mathcal{D})q(\mbox{\boldmath{$\gamma$}}^{\star}\,|\,\mbox{\boldmath{$\gamma$}})$

where $\mathcal{D}$ denotes every parameter in the model (including latent
variables). The first component is simply the full conditional distribution of
$\beta$ already derived for both models above, but only using the columns of
$\mathbf{X}$ supported by $\mbox{\boldmath{$\gamma$}}^{\star}$. The second is
the user-chosen method for proposing new values of $\gamma$. Under this
proposal scheme, the acceptance rate $\alpha$ is

$\alpha=\min\left\\{1,\frac{|\mathbf{v}_{\mbox{\tiny\boldmath{$\gamma$}}^{\star}}|^{-0.5}|\mathbf{V}_{\mbox{\tiny\boldmath{$\gamma$}}^{\star}}|^{0.5}\exp\left[0.5\left(\mathbf{M}^{\top}_{\mbox{\tiny\boldmath{$\gamma$}}^{\star}}\mathbf{V}^{-1}_{\mbox{\tiny\boldmath{$\gamma$}}^{\star}}\mathbf{M}_{\mbox{\tiny\boldmath{$\gamma$}}^{\star}}+\mathbf{m}_{\mbox{\tiny\boldmath{$\gamma$}}^{\star}}^{\top}\mathbf{v}^{-1}_{\mbox{\tiny\boldmath{$\gamma$}}^{\star}}\mathbf{m}_{\mbox{\tiny\boldmath{$\gamma$}}^{\star}}\right)\right]q(\mbox{\boldmath{$\gamma$}}\,|\,\mbox{\boldmath{$\gamma$}}^{\star})\pi_{\mbox{\tiny\boldmath{$\gamma$}}}(\mbox{\boldmath{$\gamma$}}^{\star})}{|\mathbf{v}_{\mbox{\tiny\boldmath{$\gamma$}}}|^{-0.5}|\mathbf{V}_{\mbox{\tiny\boldmath{$\gamma$}}}|^{0.5}\exp\left[0.5\left(\mathbf{M}^{\top}_{\mbox{\tiny\boldmath{$\gamma$}}}\mathbf{V}^{-1}_{\mbox{\tiny\boldmath{$\gamma$}}}\mathbf{M}_{\mbox{\tiny\boldmath{$\gamma$}}}+\mathbf{m}_{\mbox{\tiny\boldmath{$\gamma$}}}^{\top}\mathbf{v}^{-1}_{\mbox{\tiny\boldmath{$\gamma$}}}\mathbf{m}_{\mbox{\tiny\boldmath{$\gamma$}}}\right)\right]q(\mbox{\boldmath{$\gamma$}}^{\star}\,|\,\mbox{\boldmath{$\gamma$}})\pi_{\mbox{\tiny\boldmath{$\gamma$}}}(\mbox{\boldmath{$\gamma$}})}\right\\}$

where $\pi_{\mbox{\tiny\boldmath{$\gamma$}}}(\cdot)$ is the prior of $\gamma$
and $q(\mbox{\boldmath{$\gamma$}})$ is the proposal distribution. Here,
$\mathbf{M}_{\mbox{\tiny\boldmath{$\gamma$}}}$ and
$\mathbf{V}_{\mbox{\tiny\boldmath{$\gamma$}}}$ correspond to the mean and
variance of the full conditional distribution of $\beta$ for the model
described by $\gamma$. Notice that, as written, $\beta$ is marginalised so no
further adjustments to the acceptance probability is required. This makes
sense as ‘removing’ a variable from the above regression methods is the same
as simply fixing its associated slope at zero, so there is actually no shift
in the dimensionality when any element of $\gamma$ is altered. This is
obviously no longer the case if $h(\cdot)$ is a DLF (or similar), in which
case altering $\gamma$ will now result in genuine dimension jumping and the
acceptance rate will need to be adjusted accordingly (Green,, 1995; Green and
Hastie,, 2009).

In regards to the proposal $q(\mbox{\boldmath{$\gamma$}})$, the approach we
assume for the rest of this paper is that, for each MCMC iteration, one random
element from $\gamma$ is chosen to be swapped from 0 to 1 or vice-versa. This
is a symmetric proposal, so the $q(\cdot)$ terms cancel in the acceptance
probability equation. Note this becomes unsuitable as the number of potential
predictors grow larger since the probability that any specific value is
swapped at any iteration shrinks to zero. Adapting a horseshoe prior (Carvalho
et al.,, 2010; Piironen and Vehtari,, 2017) can be beneficial in such
situations.

##  3 Simulation Study

A simulation study was performed to test the above methods. The simulated data
are partially based on a table of air quality metrics collected between 1987
and 2000 in Chicago, available through the R package dlnm (Gasparrini,, 2011).
Three columns from this dataset are used as dynamic predictors: mean relative
humidity, PM10, and ozone (O3). Two further static variables are generated by
taking the absolute value of draws from the standard normal distribution. The
full Chicago data has 5114 observations; we will use the full dataset for one
set of simulations, whereas we will only use the first 250 observations for
another. Some humidity and PM10 values in the dataset are missing; these are
imputed via predictive mean matching using the mice package in R (van Buuren
and Groothuis-Oudshoorn,, 2011).

The lag-response relationship between the real dynamic variables and the
simulated response were as follows: humidity had a unimodal delayed effect,
PM10 had a linearly decreasing effect, and 03 had a bimodal effect. The
delayed and bimodal effects were created using a 2 degree exponential
polynomial (Almon,, 1965):

$w(k,\mbox{\boldmath{$\theta$}})\propto\exp\left(\theta_{1}k+\theta_{2}k^{2}\right)$
(4)

where $w(\cdot,\cdot)$ is scaled across all lags $k$ within each time window
such that $\sum_{k=0}^{\tau}w(k,\mbox{\boldmath{$\theta$}})=1$. This is so the
lag-response can be interpreted as a weighted average over the lags (the same
scaling also applies for the linear decreasing lag-response). For the delayed
effect, we use $\mbox{\boldmath{$\theta$}}=\\{0.32,-0.02\\}$. For the bimodal
effect we combine two Almon functions $w_{1}$ and $w_{2}$ with
$\mbox{\boldmath{$\theta$}}_{1}=\\{0.2,-0.03\\}$ and
$\mbox{\boldmath{$\theta$}}_{2}=\\{6,-0.1\\}$. We then take the effect to be
$0.7w_{1}+0.3w_{2}$. The resultant lag response profiles are shown as the
dashed blue lines in Figure 1.

The overall effect sizes were set as follows: the two static variables have an
overall effect of $-0.5$ and $0.01$. Humidity, PM10 and ozone had overall
effects sizes (of their lag-response relationships) of $0.5$, $0.01$, and
$-0.5$ respectively. This gives us a mix of strong and weak effects to test
the variable selection scheme. Lag-responses were modelled using cubic
b-splines with three evenly spaced knots.

Datasets were created with both a count response and a binary response. The
count response simulations used the negative binomial distribution with
stopping parameter equal to $50$, and intercept equal to $-1$. The binary
response simulations were simulated by drawing random values from the ALD with
the skew parameter equal to $0.9$ and setting the intercept of the regression
model to 0, with the binary response then being set to 1 for positive values
and 0 otherwise.

In summation, there are four distinct types of simulations; Large/small
datasets with count/binary responses. Four of each dataset was created, and
for each of these datasets we fit the appropriate model 30 times to assess the
Monte Carlo consistency. For brevity only one set of the large simulations
(for both count and binary simulations) will be discussed in this manuscript;
see Figures 1–4. Relevant visualisations for _all_ simulation tests are
compiled in the supplementary document.

(a) Count data simulation results.

(b) Binary data simulation results.

Figure 1: Estimated lag-responses for simulated data compared to their true
values. Grey lines represent each of the 30 MCMC runs, the thick black line is
the average across all runs, and the blue dashed lines are the correct values.
The correct values are the same for both groups, they only appear different
here due to the differing scales of the y-axes. PM10 is omitted here as its
small effect size means that its associated lag-response is ill-defined,
though it is included in the supplementary results.

(a) Count data simulation results.

(b) Binary data simulation results.

Figure 2: Covariate inclusion posteriors for simulated data across all 30 MCMC
runs (intercept omitted).

(a) Count data simulation results.

(b) Binary data simulation results.

Figure 3: The difference between the mean estimated effect sizes and their
true values for simulated data across all 30 MCMC runs. Note that the y-axes
of the above graphs are on different scales. Figure 4: Posterior of the
negative binomial stopping parameter for simulated data. Each grey line
represents a posterior from one of the thirty MCMC runs. The blue dashed line
denotes the true value.

From Figure 1 we see that the MCMC routine was able to accurately infer the
lag-response relationship for relative humidity and ozone for the count-
response data. Binary simulations were somewhat less successful; the lag-
response relationship with relative humidity, while overall accurate, is less
consistent, and the algorithm failed to pick up the smaller mode for the ozone
lag-response relationship. PM10 was omitted as the overall effect is so small
that the lag-response is ill-defined and thus expected to fail. However the
attempt at the estimation is shown in the supplementary materials.

Figure 2 shows the variable inclusion parameters. The weaker effects are
largely discarded from the model with a high degree of certainty. More
interestingly, we see that $\gamma$ are discarding individual b-spline
components even from the non-zero effects. This behaviour could be potentially
leveraged to inform about the parameters of the b-spline construction (knots,
component polynomial degree, etc). The binary response simulations, while
performing still well overall, display a lack of sampling consistency,
especially with the ozone variable.

The difference between estimated and true effect sizes are shown in Figure 3.
We see that the effects were mostly accurate, though some effects were
overestimated, in particular the intercept and humidity slope in the negative
binomial simulations. The other simulations we performed, shown in the
supplementary materials, show further varying degrees of over/under estimation
for some dynamic effects for both count and binary response simulations. By
contrast, the static effects (apart from the intercept) are almost always very
accurate with a high degree of Monte Carlo consistency. But these are simply
regression slopes; on the other hand, dynamic effects are intertwined with
lag-response estimation, increasing the statistical complexity which likely
accounts for the poorer performance. Another thing worth pointing out here are
the differences in the count and binary simulations - as before, imbalanced
binary response is harder to infer from (note the difference in y-axes scales
in Figure 3).

Finally, 4 displays the density plot of the MCMC posterior samples of the
negative binomial stopping parameter. The Monte Carlo posterior mode tends to
slightly underestimate, but the true value still sits comfortably within high
posterior density regions.

Overall, we would say the algorithms perform well, especially for the count
simulations. The fact that binary regression performs somewhat worse is
understandable given that binary data carries less information to infer from -
especially when imbalanced, as is the case for these simulations.

##  4 Chicago Air Quality Data

We fit DLM models to the Chicago air quality dataset that was used to
partially simulate response data in the previous section. This data was
originally part of the National Morbidity, Mortality and Air Pollution Study
database, but to our knowledge the full dataset is no longer available. The
data we use here is available through the _dlnm_ package (Gasparrini,, 2011).
In addition to the relative humidity, PM10 levels, and ozone levels that we
spoke of in the simulation section, the full dataset has two more
environmental variables: mean daily temperature and dew point temperature,
both of which are heavily correlated (Pearson $\approx$ 0.95). It also records
the day of the week and the month of the year the air quality metrics were
recorded. Month of the year is obviously useful for coding seasonality, and
day of week is also known to have an impact on daily reported numbers due to
the nature of healthcare surveillance (Buckingham-Jeffery et al.,, 2017). Day
of the week will be dummy-encoded and included in the model (Sunday excluded).
The months were aggregated into their respective seasons, which was dummy
encoded and included in the model (Autumn excluded).

There are three different measures of daily death rates in the data: the
number of cardiovascular-related deaths, the number of respiratory related
deaths, and the overall number of non-accidental deaths. The daily overall
death rate was chosen as the response; the others were discarded. Missing
values are imputed using predictive mean matching via the mice package (van
Buuren and Groothuis-Oudshoorn,, 2011) with default settings.

The remaining air quality metrics are included as dynamic variables. In each
case we assume a lag length of 40 days, and the lag-response are modelled
using cubic b-splines with 3 evenly spaced knots.

We fit a negative binomial DLM to the dataset as described. We also fit a
binary quantile regression DLM to the dataset by dichotimising the response
variable; any day that had in excess of 135 deaths (roughly 10% of the data)
was assigned a value of 1, and 0 otherwise. Please note that we are not
condoning dichotomisation for serious analysis of this data or others - this
was done purely for illustrative purposes.

For the MCMC fits shown, we adopt independent Gaussian distributions centered
at 0 with variance 100 as priors for the regression coefficients. We use a
Bernoulli distribution with probability parameter equal to 0.5 as prior for
each of the variable inclusion parameters. A Gamma distribution with shape and
rate equal to 2 and 1/50 respectively was used for the stopping parameter
prior. We run the MCMC algorithm 4 times each with their own starting values.
One run starts with the intercept only model, another starts at the full
inclusion model, and the other two use random starts; each non-intercept
predictor is included with probability 0.5. They are all ran for $10^{5}$
iterations, half of which is discarded, and only every $10^{\text{th}}$ value
is retained thereafter.

(a) Covariate inclusion posterior across all 4 runs.

(b) Monte Carlo means of the response-lag relationship for all 4 runs.

Figure 5: Results from the negative binomial fit of the Chicago air quality
data.

(a) Covariate inclusion posterior across all 4 runs.

(b) Monte Carlo means of the response-lag relationship for all 4 runs.

Figure 6: Results from the binary quantile regression fit of the Chicago air
quality data.

The posterior of the covariate inclusion parameters and the means of the
estimated lag-response relationships are given in Figures 5 and 6. The
negative binomial DLM shows very consistent results for both the inclusion
indicators and the lag-responses, even with the large degree of
multicolinearity between temperature and dew-point temperature. There was more
uncertainty within the quantile regression results, as is expected due to
binary imbalance. Both models tend to agree on what variables should be
included, with relatively strong levels of inclusion for temperature, PM10,
and ozone. There’s some agreement with the estimated lag-responses too; PM10
has a 7-day delayed effect on death rates, PM10 has an immediate effect
followed by a sharp decline, and ozone has both an immediate effect and a
30-day delayed effect on death rates. They both disagree on the relationship
between death rates and dew point temperature and humidity, though the models
agreed that these predictors were mostly irrelevant regardless.

##  5 Conclusion

Latent variables are frequently leveraged when fitting Bayesian GLMs
(especially binary response models) as they result in fully tractable
conditional distributions and require no further tuning beyond setting priors
and starting values for MCMC. Here we focused on the implementation of
quantile binary regression and negative binomial regression, but any similar
latent variable approach can be incorporated just as easily. A known downside
of these approaches, beyond the fact that some carry a large computational
cost, is that in certain situations (such as when binary data is heavily
imbalanced) there can be severe autocorrelation in the MCMC posterior samples
due to the close relationship between the parameters and the latent variables,
and may be outperformed by a Metropolis-Hastings approach (Johndrow et al.,,
2018). Zens et al., (2023) proposed using parameter expansion (Liu and Wu,,
1999) to improve the efficiency of fitting a logistic regression model via
Polya-Gamma latent variables (Polson et al.,, 2013). Duan et al., (2018)
suggest improvements to the efficiency by scaling the conditional posterior
variance; these methods and others can be incorporated to the DLM framework
used here if desired.

Beyond the concerns around how the model is fit, we must also consider the
model itself and its validity in capturing lag-response dynamics. In this
manuscript we focused purely on b-splines to model the lag-response. These are
extremely easy to fit as they are simply linear predictors like any other; a
similar (though overly-parameterised approach) would be to include the raw
lags as predictors. However the disadvantage of these approaches are the
potentially irreconcilable interpretations of the fitted model; using linear
predictors to model the lag-response, it is possible that a fitted model may
say that a treatment is both beneficial _and_ detrimental, depending on how
far back in time the subject was exposed to it. While there may be some
situations where this is a reasonable outcome, we cannot imagine a scenario
where that would be true for the effect of air quality exposure to human
health, a topic where these models are commonly applied. Functional
constraints on the lag response such as distributed lag functions often
encountered in MIDAS model literature, can rectify this. Even a simple two
parameter model, such as (4) can be enough to model the delayed nature of
accumulating lag-response, _and_ enforce mono-directional lag effects. The
drawback is that they are more difficult to fit, as the full-conditional
distribution of the parameters are sampling-intractable. Also, since the DLF
parameters are linked to each variable and disappear from the model when the
associated predictors are removed, they need to be given special consideration
during variable inclusion update proposals, i.e. the ‘birth’ and ‘death’
process (Green,, 1995; Green and Hastie,, 2009).

All relevant code for this manuscript can be found in the following link:
https://github.com/DanDempsey/Latent-Variable-DLM-Project/tree/main/Code.

## Acknowledgements

Daniel Dempsey’s work was funded by an Irish Research Council Enterprise
Partnership Scheme Award EPSPG/2017/304. Jason Wyse’s work was partially
supported through a Science Foundation Ireland Frontiers for the Future Award
21/FFP$-$P/10123 and the European Union’s Horizon 2020 research and innovation
programme under the Marie Skłodowska-Curie grant agreement No. 813545.

## References

  * Albert and Chib, (1993) Albert, J. H. and Chib, S. (1993).  Bayesian Analysis of Binary and Polychotomous Response Data.  Journal of the American Statistical Association, 88(422):669–679.
  * Almon, (1965) Almon, S. (1965).  The distributed lag between capital appropriations and expenditures.  Econometrica: Journal of the Econometric Society, pages 178–196.
  * Antonelli et al., (2021) Antonelli, J., Wilson, A., and Coull, B. (2021).  Multiple exposure distributed lag models with variable selection.
  * Benoit and den Poel, (2017) Benoit, D. F. and den Poel, D. V. (2017).  bayesQR: A Bayesian Approach to Quantile Regression.  Journal of Statistical Software, 76(1):1–32.
  * Benoit and Van den Poel, (2012) Benoit, D. F. and Van den Poel, D. (2012).  Binary quantile regression: a bayesian approach based on the asymmetric laplace distribution.  Journal of Applied Econometrics, 27(7):1174–1188.
  * Buckingham-Jeffery et al., (2017) Buckingham-Jeffery, E., Morbey, R., House, T., Elliot, A. J., Harcourt, S., and Smith, G. E. (2017).  Correcting for day of the week and public holiday effects: improving a national daily syndromic surveillance service for detecting public health threats.  BMC Public Health, 17(1):1–9.
  * Carvalho et al., (2010) Carvalho, C. M., Polson, N. G., and Scott, J. G. (2010).  The horseshoe estimator for sparse signals.  Biometrika, 97(2):465–480.
  * Chen et al., (1972) Chen, D., Courtney, R., and Schmitz, A. (1972).  A polynomial lag formulation of milk production response.  American Journal of Agricultural Economics, 54(1):77–83.
  * Chipman et al., (2010) Chipman, H. A., George, E. I., and McCulloch, R. E. (2010).  BART: Bayesian additive regression trees.  The Annals of Applied Statistics, 4(1):266 – 298.
  * Clements and Galvão, (2008) Clements, M. P. and Galvão, A. B. (2008).  Macroeconomic forecasting with mixed-frequency data.  Journal of Business & Economic Statistics, 26(4):546–554.
  * Czado and Santner, (1992) Czado, C. and Santner, T. J. (1992).  The effect of link misspecification on binary regression inference.  Journal of statistical planning and inference, 33(2):213–231.
  * D’Angelo and Canale, (2022) D’Angelo, L. and Canale, A. (2022).  Efficient posterior sampling for bayesian poisson regression.  arXiv preprint arXiv:2109.09520.
  * Dhrymes, (1981) Dhrymes, P. J. (1981).  Distributed lags; problems of estimation and formulation, volume 2.  Holden-Day.
  * Duan et al., (2018) Duan, L. L., Johndrow, J. E., and Dunson, D. B. (2018).  Scaling up data augmentation mcmc via calibration.  The Journal of Machine Learning Research, 19(1):2575–2608.
  * Frühwirth-Schnatter and Frühwirth, (2010) Frühwirth-Schnatter, S. and Frühwirth, R. (2010).  Data augmentation and mcmc for binary and multinomial logit models.  Statistical modelling and regression structures: Festschrift in honour of Ludwig Fahrmeir, pages 111–132.
  * Gasparrini, (2011) Gasparrini, A. (2011).  Distributed lag linear and non-linear models in R: the package dlnm.  Journal of Statistical Software, 43(8):1–20.
  * Ghysels et al., (2016) Ghysels, E., Kvedaras, V., and Zemlys, V. (2016).  Mixed Frequency Data Sampling Regression Models: The R Package midasr.  Journal of Statistical Software, 72(1):1–35.
  * Ghysels et al., (2004) Ghysels, E., Santa-Clara, P., and Valkanov, R. (2004).  The midas touch: Mixed data sampling regression models.  UCLA: Finance.
  * Green and Hastie, (2009) Green, P. and Hastie, D. (2009).  Reversible jump MCMC.  Genetics, 155.
  * Green, (1995) Green, P. J. (1995).  Reversible Jump Markov Chain Monte Carlo Computation and Bayesian Model Determination.  Biometrika, 82(4):711–732.
  * Hannan, (1965) Hannan, E. J. (1965).  The estimation of relationships involving distributed lags.  Econometrica: Journal of the Econometric Society, pages 206–224.
  * Hastie and Tibshirani, (1987) Hastie, T. and Tibshirani, R. (1987).  Generalized additive models: some applications.  Journal of the American Statistical Association, 82(398):371–386.
  * Holmes and Held, (2006) Holmes, C. C. and Held, L. (2006).  Bayesian auxiliary variable models for binary and multinomial regression.  Bayesian Analysis, 1(1):145–168.
  * Johndrow et al., (2018) Johndrow, J. E., Smith, A., Pillai, N., and Dunson, D. B. (2018).  Mcmc for imbalanced categorical data.  Journal of the American Statistical Association.
  * Jorgensen, (2012) Jorgensen, B. (2012).  Statistical properties of the generalized inverse Gaussian distribution, volume 9.  Springer Science & Business Media.
  * Kozumi and Kobayashi, (2011) Kozumi, H. and Kobayashi, G. (2011).  Gibbs sampling methods for Bayesian quantile regression.  Journal of Statistical Computation and Simulation, 81(11):1565–1578.
  * Li et al., (2021) Li, X., Yu, H., Xie, Y., and Li, J. (2021).  Attention-based novel neural network for mixed frequency data.  CAAI Transactions on Intelligence Technology, 6(3):301–311.
  * Liu and Wu, (1999) Liu, J. S. and Wu, Y. N. (1999).  Parameter expansion for data augmentation.  Journal of the American Statistical Association, 94(448):1264–1274.
  * Meier et al., (2008) Meier, L., Van De Geer, S., and Bühlmann, P. (2008).  The group lasso for logistic regression.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 70(1):53–71.
  * Mitchell and Beauchamp, (1988) Mitchell, T. J. and Beauchamp, J. J. (1988).  Bayesian variable selection in linear regression.  Journal of the american statistical association, 83(404):1023–1032.
  * Mogliani and Simoni, (2020) Mogliani, M. and Simoni, A. (2020).  Bayesian MIDAS Penalized Regressions: Estimation, Selection, and Prediction.  arXiv:1903.08025 [econ].
  * Mork et al., (2021) Mork, D., Kioumourtzoglou, M.-A., Weisskopf, M., Coull, B. A., and Wilson, A. (2021).  Heterogeneous distributed lag models to estimate personalized effects of maternal exposures to air pollution.  arXiv preprint arXiv:2109.13763.
  * Piironen and Vehtari, (2017) Piironen, J. and Vehtari, A. (2017).  Sparsity information and regularization in the horseshoe and other shrinkage priors.  Electronic Journal of Statistics, 11(2):5018 – 5051.
  * Pillow and Scott, (2012) Pillow, J. and Scott, J. (2012).  Fully bayesian inference for neural models with negative-binomial spiking.  Advances in neural information processing systems, 25.
  * Polson et al., (2013) Polson, N. G., Scott, J. G., and Windle, J. (2013).  Bayesian Inference for Logistic Models Using Pólya–Gamma Latent Variables.  Journal of the American Statistical Association, 108(504):1339–1349.
  * Schwartz, (2000) Schwartz, J. (2000).  The Distributed Lag between Air Pollution and Daily Deaths.  Epidemiology, 11(3):320–326.
  * Sims, (1971) Sims, C. A. (1971).  Discrete approximations to continuous time distributed lags in econometrics.  Econometrica: Journal of the Econometric Society, pages 545–563.
  * van Buuren and Groothuis-Oudshoorn, (2011) van Buuren, S. and Groothuis-Oudshoorn, K. (2011).  mice: Multivariate imputation by chained equations in r.  Journal of Statistical Software, 45(3):1–67.
  * Warren et al., (2020) Warren, J. L., Luben, T. J., and Chang, H. H. (2020).  A spatially varying distributed lag model with application to an air pollution and term low birth weight study.  Journal of the Royal Statistical Society. Series C, Applied statistics, 69(3):681.
  * Wilson et al., (2017) Wilson, A., Chiu, Y.-H. M., Hsu, H.-H. L., Wright, R. O., Wright, R. J., and Coull, B. A. (2017).  Bayesian distributed lag interaction models to identify perinatal windows of vulnerability in children’s health.  Biostatistics, 18(3):537–552.
  * Xu et al., (2019) Xu, Q., Zhuo, X., Jiang, C., and Liu, Y. (2019).  An artificial neural network for mixed frequency data.  Expert Systems with Applications, 118:127–139.
  * Xu and Ghosh, (2015) Xu, X. and Ghosh, M. (2015).  Bayesian variable selection and estimation for group lasso.  Bayesian Analysis, 10(4):909–936.
  * Yu and Zhang, (2005) Yu, K. and Zhang, J. (2005).  A Three-Parameter Asymmetric Laplace Distribution and Its Extension.  Communications in Statistics - Theory and Methods, 34(9-10):1867–1879.
  * Zanobetti et al., (2002) Zanobetti, A., Schwartz, J., Samoli, E., Gryparis, A., Touloumi, G., Atkinson, R., Le Tertre, A., Bobros, J., Celko, M., Goren, A., et al. (2002).  The temporal pattern of mortality responses to air pollution: a multicity assessment of mortality displacement.  Epidemiology, pages 87–93.
  * Zens et al., (2023) Zens, G., Frühwirth-Schnatter, S., and Wagner, H. (2023).  Ultimate p$\backslash$’olya gamma samplers–efficient mcmc for possibly imbalanced binary and categorical data.  arXiv preprint arXiv:2011.06898.
  * Zhou et al., (2012) Zhou, M., Li, L., Dunson, D., and Carin, L. (2012).  Lognormal and gamma mixed negative binomial regression.  In Proceedings of the… International Conference on Machine Learning. International Conference on Machine Learning, volume 2012, page 1343\. NIH Public Access.