# A Bayesian framework for incorporating exposure uncertainty into health
analyses with application to air pollution and stillbirth

Saskia Comess1, Howard H. Chang2, Joshua L. Warren3∗

(1Emmett Interdisciplinary Program in Environment and Resources, Stanford
University, Stanford, CA 94305, USA
2Department of Biostatistics and Bioinformatics, Emory University, Atlanta, GA
30329, USA
3Department of Biostatistics, Yale University, New Haven, CT 06510, USA
(∗E-mail: joshua.warren@yale.edu))

###### Abstract

Studies of the relationships between environmental exposures and adverse
health outcomes often rely on a two-stage statistical modeling approach, where
exposure is modeled/predicted in the first stage and used as input to a
separately fit health outcome analysis in the second stage. Uncertainty in
these predictions is frequently ignored, or accounted for in an overly
simplistic manner, when estimating the associations of interest. Working in
the Bayesian setting, we propose a flexible kernel density estimation (KDE)
approach for fully utilizing posterior output from the first stage
modeling/prediction to make accurate inference on the association between
exposure and health in the second stage, derive the full conditional
distributions needed for efficient model fitting, detail its connections with
existing approaches, and compare its performance through simulation. Our KDE
approach is shown to generally have improved performance across several
settings and model comparison metrics. Using competing approaches, we
investigate the association between lagged daily ambient fine particulate
matter levels and stillbirth counts in New Jersey (2011-2015), observing an
increase in risk with elevated exposure three days prior to delivery. The
newly developed methods are available in the R package KDExp.

Keywords: Air pollution; Environmental health; Kernel density estimation;
Stillbirth; Two-stage modeling; Uncertainty propagation.

##  1 Introduction

Environmental health studies require the linking of environmental exposure
information for each observation in the analysis (e.g., individual or time
point) in order to estimate the association with adverse health outcomes.
Because exposure data are typically not available at every spatial location
and time period covered by the study, researchers often rely on predictions
from a first stage statistical model to fill in the spatiotemporal gaps. For
example, several advanced statistical methods have been developed to
interpolate ambient air pollution concentrations using monitoring data
combined with estimates from deterministic models and other data sources
(e.g., Fuentes and Raftery (2005); Berrocal et al. (2010b, a); McMillan et al.
(2010); Berrocal et al. (2012); Reich et al. (2014); Guan et al. (2019);
Warren et al. (2021)).

These methods yield point predictions of exposure at the relevant locations
and times of interest, but also provide measures of uncertainty. Because many
of the methods are fitted within a Bayesian framework using Markov chain Monte
Carlo (MCMC) techniques, samples from posterior predictive distributions (ppd)
are also available. Incorporating this exposure uncertainty into the
subsequent health analysis is important for correctly characterizing
uncertainty in the association between exposure and health. Specifically,
health observations linked with exposure estimates with higher uncertainties
should contribute less to the overall health effect estimate. However,
previous studies often ignore this uncertainty entirely which may impact
inference for these associations, although the full implications of this
approach are not currently clear. Several methods for propagating exposure
uncertainty have been developed and used previously (e.g., Gryparis et al.
(2009); Lee and Shaddick (2010); Peng and Bell (2010); Chang et al. (2011);
Szpiro et al. (2011); Warren et al. (2012); Szpiro and Paciorek (2013);
Blangiardo et al. (2016); Lee et al. (2017); Huang et al. (2018)), and we
provide full details on many of them in Section 2. Although prior work in this
area has recommended an investigation and comparison of their performances
(Lee et al., 2017), to our knowledge such an analysis has yet to be conducted.

In this work, we present a flexible framework for exposure uncertainty
propagation, carry out a simulation study to compare its performance to
existing methods, and apply several of the methods to better understand the
relationship between acute, population-level fine particulate matter (PM2.5)
exposure and daily stillbirth counts using data from three counties in New
Jersey (NJ), 2011-2015. The proposed framework uses kernel density estimation
(KDE) with a Gaussian kernel function to specify prior distributions for the
exposures within the health model. The resulting model fitting derivations
suggest that this represents a hybrid between two existing approaches;
allowing for more flexibility while also avoiding some of the limiting
assumptions of those approaches. It is also shown to maintain computational
efficiency for several common health outcome analysis types. Through
simulation, we show that the new approach is flexible enough to accurately
characterize uncertainty in the predictions, leading to improved estimation of
the association of interest in the health analysis compared to existing
approaches. Differences between the methods are also observed in the NJ
stillbirth case study results, indicating the importance of selecting the
optimal method in future applications.

##  2 Background

We specify that the primary epidemiological health outcome analysis of
interest consists of $n$ data points (e.g., individuals, time periods) where
an exposure level (e.g., ambient air pollution) is assigned to each data point
in order to examine its association with the outcome. For presentation
purposes, we introduce the framework using a single exposure while noting that
it is straightforward to extend the following results/derivations to
accommodate multiple additive exposures. Statistical modeling of the exposures
and health outcome are assumed to take place in the Bayesian setting, as is
common in the environmental health statistical methodology literature (e.g.,
Berrocal et al. (2010b); Chang et al. (2011)).

We assume that $m$ samples from the exposure ppd have been obtained based on
the modeling and prediction of observed exposures in a first stage Bayesian
framework. While many modeling options are available, the end result is the
same across all approaches; an $n$ by $m$ matrix of ppd samples (i.e.,
$\textbf{Z}^{*}$) is obtained such that

$\textbf{Z}^{*}=\begin{bmatrix}\text{z}_{11}^{*}&\text{z}_{12}^{*}&\text{z}_{13}^{*}&\dots&\text{z}_{1m}^{*}\\\
\text{z}_{21}^{*}&\text{z}_{22}^{*}&\text{z}_{23}^{*}&\dots&\text{z}_{2m}^{*}\\\
\vdots&\vdots&\vdots&\ddots&\vdots\\\
\text{z}_{n1}^{*}&\text{z}_{n2}^{*}&\text{z}_{n3}^{*}&\dots&\text{z}_{nm}^{*}\end{bmatrix}.$

In this matrix, $\text{z}^{*}_{ij}$ represents the $j^{th}$ ppd sample of
exposure for health data point $i$. In terms of notation, it is helpful to
define the exposure matrix in terms of row and column vectors such that

$\textbf{Z}^{*}=\begin{bmatrix}\textbf{z}_{1.}^{*\text{T}}\\\ \vdots\\\
\textbf{z}_{n.}^{*\text{T}}\end{bmatrix}=\begin{bmatrix}\textbf{z}_{.1}^{*}&\ldots&\textbf{z}_{.m}^{*}\end{bmatrix}$

where $\textbf{z}^{*}_{i.}$ represents the complete set of exposure ppd
samples for data point $i$ and $\textbf{z}^{*}_{.j}$ is the vector of
exposures for all $n$ data points from the $j^{\text{th}}$ ppd sample.

Within a column of $\textbf{Z}^{*}$ the exposures for the different data
points could be collected independently (i.e., from the marginal ppds) or
jointly (i.e., from the joint ppd). Sampling from the marginal ppds, which may
be necessary due to computational considerations when working with a large
spatial/temporal domain, results in independence across the rows of the
$\textbf{Z}^{*}$, whereas sampling from the joint ppd retains the correlation
across the rows. In either case, trend in exposures across the rows may be
present depending on the structure of the data points (e.g., air pollution
concentrations across time). We note that independence across the columns of
$\textbf{Z}^{*}$ in practice is achieved through Monte Carlo sampling or MCMC
sampling and thinning of the collected ppd samples.

Once the exposure ppd samples are obtained from the first stage modeling, they
are used in a subsequent epidemiological health outcome analysis to determine
their association with the outcome. This second stage health model typically
follows a regression framework of the form

$\displaystyle\begin{split}&Y_{i}|\mu_{i},\boldsymbol{\zeta}\stackrel{{\scriptstyle\text{ind}}}{{\sim}}f\left(y|\mu_{i},\boldsymbol{\zeta}\right),\
i=1,\ldots,n,\\\
&g\left(\mu_{i}\right)=\text{O}_{i}+\textbf{x}_{i}^{\text{T}}\boldsymbol{\beta}+\text{z}_{i}\theta\end{split}$
(1)

where $n$ is the previously defined number of data points; $Y_{i}$ represents
the health outcome for data point $i$; $f(.|.)$ is the probability density
function (pdf) of the outcome; $\mu_{i}$ is the mean of this distribution with
$\boldsymbol{\zeta}$ representing additional parameters that often define
variance/dispersion or auxiliary variables used to improved the efficiency of
posterior sampling; $g\left(.\right)$ is a link function to connect the mean
with a set of covariates; $\text{O}_{i}$ is an offset term sometimes used in
the modeling of count data but will be zero otherwise; and $\textbf{x}_{i}$ is
a vector of covariates unrelated to the primary exposure of interest,
including an intercept term, with $\boldsymbol{\beta}$ the vector of
corresponding regression parameters.

The true but unobserved exposure for data point $i$ is denoted by
$\text{z}_{i}$, where $\theta$ describes the association between exposure and
outcome. In this work, we assume that the set of true exposures,
$\textbf{z}^{\text{T}}=\left(\text{z}_{1},\ldots,\text{z}_{n}\right)$, follows
the ppd derived in the first stage analysis. However, we avoid the likely
unrealistic assumptions made in some previous work that z is included as one
of the columns of $\textbf{Z}^{*}$ (i.e., $\textbf{z}\neq\textbf{z}^{*}_{.j}$
for any $j=1,\ldots,m$) (e.g., Peng and Bell (2010); Chang et al. (2011)).
Instead, we treat entries of z as unknown parameters in (1) with the first
stage ppd representing our current state of knowledge about their values.
Therefore, z is thought to arise from the same process that produced the
columns of $\textbf{Z}^{*}$, but is not actually observed in the finite set of
$m$ samples collected in the first stage. To complete the model specification,
we assign weakly informative prior distributions for each of the introduced
parameters in (1) with specific settings based on the likelihood choice.

###  2.1 Existing approaches

Given $\textbf{Z}^{*}$ and the fully specified health model from (1), the
question becomes how to efficiently utilize the information contained in the
full set of ppd samples to accurately quantify uncertainty in the exposures
when making inference on $\theta$. A number of approaches, ranging in
conceptual and computational complexity, have been proposed and we detail
several of them below. In Figure S1 of the Supplement, we present an overview
of the different approaches.

#### Plug-in exposures

The simplest approach used in previous work replaces $\text{z}_{i}$ from (1)
with $\widehat{\text{z}}_{i}=T\left(\textbf{z}_{i.}^{*}\right)$ where
$T\left(.\right)$ is a function of the input ppd samples (e.g., median, mean)
(e.g., Warren et al. (2021)). We refer to this as the Plug-in approach. Using
only a summary measure of the ppd samples entirely ignores uncertainty in the
exposures which may affect uncertainty estimation for $\theta$.

#### Multiple imputation

The multiple imputation (MI) approach incorporates uncertainty in the
exposures by fitting the health model in (1) separately for each of the $m$
columns of $\textbf{Z}^{*}$, replacing z with $\textbf{z}^{*}_{.j}$
(Blangiardo et al., 2016). During model fit $j$, $s$ posterior samples (post-
convergence and possibly thinned) from $\theta$ are collected and denoted as
$\boldsymbol{\theta}^{(j)\text{T}}=\left(\theta_{1}^{(j)},\ldots,\theta_{s}^{(j)}\right)$.
After fitting the model to all columns of $\textbf{Z}^{*}$, posterior
inference is conducted based on the combined samples across all $m$ model
fits;
$\boldsymbol{\theta}^{\text{T}}=\left(\boldsymbol{\theta}^{(1)\text{T}},\ldots,\boldsymbol{\theta}^{(m)\text{T}}\right)$
(Zhou and Reiter, 2010). Depending on how long it takes to fit the health
model in (1), which is impacted by the likelihood choice and sample size, MI
may be computationally demanding as $m$ increases. This method also assumes
that the columns of $\textbf{Z}^{*}$ resemble z, which is not necessarily true
and depends on the amount of variability in the first stage ppd.

#### Multiple imputation approximation

The multiple imputation approximation (MIA) approach approximates the results
from MI while only requiring a single fit of the health model in (1);
representing a major computational improvement. Specifically, during each
iteration of the MCMC algorithm developed for the health model in (1), MIA
randomly selects a new column of exposures for the $n$ data points with
replacement from $\textbf{Z}^{*}$ (i.e., $\textbf{z}^{*}_{.j}$) and completes
a full sweep of the algorithm (i.e., collecting samples from all of the
introduced model parameters) (Lee et al., 2017). Posterior inference for
$\theta$ is made based on the $s$ MCMC samples collected from this algorithm.

#### Discrete uniform prior distribution

Similar to MIA, the discrete uniform (DU) approach requires only a single fit
of the health model in (1) and uses columns directly from $\textbf{Z}^{*}$
during model fitting. However, instead of randomly selecting a column during
each MCMC iteration, DU incorporates the health data in the decision making,
resulting in the selection of more probable columns during posterior sampling.
It does so by assigning a prior distribution to z using the collected ppd
samples in $\textbf{Z}^{*}$ and carrying out full Bayesian inference for the
health model in (1). Specifically, DU makes the likely unrealistic assumption
that the true vector of exposures is contained in $\textbf{Z}^{*}$ (i.e.,
$\textbf{z}=\textbf{z}^{*}_{.j}$ for some $j=1,\ldots,m$). Based on this
assumption, a prior distribution for z is specified such that

$\text{P}\left(\textbf{z}=\textbf{z}_{.j}^{*}\right)=\frac{1}{m},\
j=1,\ldots,m$

(Peng and Bell, 2010; Chang et al., 2011). Use of DU yields a semi-conjugate
full conditional distribution for z regardless of the choice for
$f\left(.|.\right)$ in (1), allowing for convenient updating during MCMC
sampling. When $m$ becomes large, a more computationally efficient Metropolis
algorithm can be used to propose/evaluate columns from $\textbf{Z}^{*}$ using
a likelihood ratio calculation.

#### Multivariate normal prior distribution

Similar to DU, the multivariate normal (MVN) approach assigns a prior
distribution to z but avoids the assumption that
$\textbf{z}=\textbf{z}^{*}_{.j}$ for some $j=1,\ldots,m$. Specifically, a MVN
prior distribution for z is specified such that

$\textbf{z}\sim\text{MVN}\left(\widehat{\textbf{z}}=\begin{bmatrix}T\left(\textbf{z}_{1.}^{*}\right)\\\
\vdots\\\ T\left(\textbf{z}_{n.}^{*}\right)\end{bmatrix},\
\widehat{\Sigma}=\frac{1}{m-1}\sum_{j=1}^{m}\left(\textbf{z}_{.j}^{*}-\widebar{\textbf{z}}^{*}\right)\left(\textbf{z}_{.j}^{*}-\widebar{\textbf{z}}^{*}\right)^{\text{T}}\right)$

where $\widehat{\textbf{z}}$ and $T\left(.\right)$ have been previously
described, and
$\widebar{\textbf{z}}^{*}=\frac{1}{m}\sum_{j=1}^{m}\textbf{z}_{.j}^{*}$ is the
length $n$ vector of average exposures across all $m$ ppd samples (Warren et
al., 2012; Lee et al., 2017).

Updating z within an MCMC algorithm is straightforward (i.e., the vector has a
standard, closed-form full conditional distribution) for multiple likelihood
choices that cover a number of relevant health outcome data types, including
Gaussian with identity link function (continuous outcome), Bernoulli with
logit link function (binary outcome), and negative binomial with logit link
function (count data). The latter two likelihood/link function results are
made possible by the work of Polson et al. (2013). Details for deriving this
distribution are provided in Section S1 of the Supplement. However, posterior
sampling will be increasingly time consuming as $n$ increases given the large
dimension of the full conditional distribution covariance matrix.
Additionally, MVN may struggle when the shape of the exposure ppd deviates
substantially from normality (e.g., skewness, multiple modes).

##  3 Kernel density estimation prior distributions

We propose using univariate and multivariate KDE with a Gaussian kernel
function to fully leverage the information contained in $\textbf{Z}^{*}$ when
estimating the health model in (1), detail its intuitive connections with
existing approaches, and consider its computational requirements. We show that
assigning prior distributions for z based on KDE results in a more flexible
hybrid between DU and MVN, allowing us to avoid potentially problematic
assumptions made by existing methods without significantly increasing the
computational burden.

For the univariate version of KDE with a Gaussian kernel function (UKDE), the
prior distributions for $\text{z}_{i},\ i=1,\ldots,n,$ are specified
independently as

$f\left(\text{z}_{i}\right)=\frac{1}{m}\sum_{j=1}^{m}\frac{1}{\sqrt{2\pi
h_{i}^{2}}}\exp\left\\{-\frac{1}{2h_{i}^{2}}\left(\text{z}_{i}-\text{z}_{ij}^{*}\right)^{2}\right\\}$

where $h_{i}$ is the bandwidth variable that controls the level of smoothness
of the density function and is estimated using standard approaches (e.g.,
Sheather and Jones (1991)) based on the the ppd samples in
$\textbf{z}_{i.}^{*}$. This prior distribution represents a mixture of $m$
equally-weighted normal distributions centered at the observed samples from
the ppd (i.e., $\text{z}_{ij}^{*}$) with standard deviation $h_{i}$. This
allows information from each individual ppd sample to be utilized when fitting
the health model in (1) and thereby avoids relying on overly simplistic
summaries of the samples used by existing methods (e.g., Plug-in, MVN).

Because we select a Gaussian kernel function, the full conditional
distribution for each $\text{z}_{i}$ has a closed form for the likelihood/link
function combinations mentioned previously in Section 2, allowing for
convenient updates within an MCMC algorithm. Specifically, based on the health
model in (1) the full conditional distribution for $\text{z}_{i}$ is a mixture
of univariate normal distributions such that

$\displaystyle\begin{split}&f\left(\text{z}_{i}|\boldsymbol{Y},\boldsymbol{\beta},\theta,\textbf{z}_{-i},\boldsymbol{\zeta}\right)=\\\
&\sum_{j=1}^{m}\left(\frac{c_{ij}}{\sum_{k=1}^{m}c_{ik}}\right)\frac{\sqrt{\theta^{2}h_{i}^{2}\Omega_{ii}+1}}{h_{i}}\phi\left(\frac{\sqrt{\theta^{2}h_{i}^{2}\Omega_{ii}+1}}{h_{i}}\left[\text{z}_{i}-\left\\{\frac{\left(\widetilde{Y}_{i}-\text{O}_{i}-\textbf{x}_{i}^{\text{T}}\boldsymbol{\beta}\right)\theta
h_{i}^{2}\Omega_{ii}+\text{z}_{ij}^{*}}{\theta^{2}h_{i}^{2}\Omega_{ii}+1}\right\\}\right]\right)\end{split}$
(2)

where $\boldsymbol{Y}$ is the complete vector of $n$ data points;
$\textbf{z}_{-i}$ is the complete vector of true exposures with $\text{z}_{i}$
removed; $\phi\left(.\right)$ is the pdf of the standard univariate normal
distribution; $\Omega$ is an $n$ by $n$ diagonal matrix with

$\Omega_{ii}=\left\\{\begin{array}[]{cl}1/\sigma^{2}_{\epsilon},&\text{Gaussian/identity}\\\
\stackrel{{\scriptstyle\text{ind}}}{{\sim}}\text{P\'{o}lya-
Gamma}\left(1,\text{O}_{i}+\textbf{x}_{i}^{\text{T}}\boldsymbol{\beta}+\text{z}_{i}\theta\right),&\text{Bernoulli/logit}\\\
\stackrel{{\scriptstyle\text{ind}}}{{\sim}}\text{P\'{o}lya-
Gamma}\left(r+Y_{i},\text{O}_{i}+\textbf{x}_{i}^{\text{T}}\boldsymbol{\beta}+\text{z}_{i}\theta\right),&\text{Negative
binomial/logit;}\end{array}\right.$

and

$\widetilde{Y}_{i}=\left\\{\begin{array}[]{cl}Y_{i},&\text{Gaussian/identity}\\\
\left(Y_{i}-0.50\right)/\Omega_{ii},&\text{Bernoulli/logit}\\\
0.50\left(Y_{i}-r\right)/\Omega_{ii},&\text{Negative
binomial/logit.}\end{array}\right.$

In the case of Gaussian distributed health outcome data, $\boldsymbol{\zeta}$
includes the error variance parameter, $\sigma^{2}_{\epsilon}$, and in the
negative binomial case it includes the dispersion parameter, $r$. The mixture
weights in (2) are defined as

$c_{ij}=\exp\left\\{-\frac{1}{2}\left[\frac{\left(\widetilde{Y}_{i}-\text{O}_{i}-\textbf{x}_{i}^{\text{T}}\boldsymbol{\beta}\right)^{2}h_{i}^{2}\Omega_{ii}+\text{z}_{ij}^{*2}}{h_{i}^{2}}-\frac{\left\\{\left(\widetilde{Y}_{i}-\text{O}_{i}-\textbf{x}_{i}^{\text{T}}\boldsymbol{\beta}\right)\theta
h_{i}^{2}\Omega_{ii}+\text{z}_{ij}^{*}\right\\}^{2}}{h_{i}^{2}\left(\theta^{2}h_{i}^{2}\Omega_{ii}+1\right)}\right]\right\\}$
(3)

where all remaining terms have been previously described. Further details for
these derivations are provided in Section S1 of the Supplement.

The form of the full conditional distribution in (2-3) suggests that the use
of UKDE for specifying exposure prior distributions results in a hybrid
approach, combining features of DU with those of MVN. Specifically for
updating $\text{z}_{i}$, we sample from this mixture distribution in two
steps.

  1. 1.

Choosing a probable observed exposure: We compute $c_{ij}$ in (3), which
depends on the health data and ppd exposure value $\text{z}^{*}_{ij}$, for
every $j=1,\ldots,m$. Each $c_{ij}/\sum_{k=1}^{m}c_{ik}$ is computed and a
random $j$ index corresponding to $\text{z}^{*}_{ij}$ is selected based on
these probabilities.

  2. 2.

Updating exposure based on selected value: Based on the selected
$\text{z}^{*}_{ij}$ ppd sample, the true exposure $\text{z}_{i}$ is drawn from
the normal distribution in (2) whose mean depends on $\text{z}^{*}_{ij}$ and
none of the other collected ppd samples.

The first step resembles DU in that the health data are used to inform a
probable ppd sample. However, once this choice is made, DU assigns it as the
true exposure unlike UKDE. The second step resembles MVN since the true
exposure is simulated from a distribution and is therefore, not required to be
one of the observed columns of $\textbf{Z}^{*}$. However, with MVN the prior
mean of the exposures remains the same across all MCMC iterations for the
corresponding full conditional distribution update (i.e.,
$T\left(\textbf{z}^{*}_{i.}\right)$). For UKDE, the mean in (2) changes each
time a new $\text{z}^{*}_{ij}$ is selected.

Therefore, UKDE represents a more flexible alternative to DU and MVN as it
avoids the assumption that the true exposures are observed in $\textbf{Z}^{*}$
and allows the prior mean of the exposures to effectively vary in the full
conditional distribution update across MCMC iterations. UKDE should also be
better able to handle non-symmetric ppd behavior than MVN; although it does
not directly account for correlation between the exposures corresponding to
different data points. In Section S2 of the Supplement, we extend the UKDE
framework to the multivariate setting (i.e., MKDE) by specifying a prior
distribution on the entire vector z, detail similar connections with existing
approaches, and note some limitations of MKDE due to the large dimension of
the exposure data in most environmental health applications.

##  4 Simulation study

We design a simulation study to compare each of the methods detailed in
Sections 2 and 3 with respect to estimating the association between exposure
and health while accounting for exposure uncertainty.

###  4.1 Data generation

We begin by defining the health model based on (1) and using a Gaussian
likelihood with identity link function, no offset term, and no
intercept/covariates such that

$Y_{i}=\theta\text{z}_{i}+\epsilon_{i},\
\epsilon_{i}\stackrel{{\scriptstyle\text{iid}}}{{\sim}}\text{N}\left(0,1\right)$
(4)

for $i=1,\ldots,n=250$, where all terms have been previously described in
Section 2. Two-dimensional spatial locations for the $250$ individuals in the
study are simulated randomly within the unit square. We consider two different
settings for the primary risk parameter, $\theta=0$ and $\theta=1$. These
settings allows us to investigate the type I error rate and power of the
competing approaches, respectively.

Next, we simulate exposure data for analysis. In each setting we define a
unique ppd and use it to create $\textbf{Z}^{*}$, which in this case
represents a $250\times 1,000$ matrix of exposures (i.e., $m=1,000$). We
define the true exposures by simulating an additional vector from the same ppd
and assigning it to z. The true exposures are used to simulate the health data
from (4) and are not included in $\textbf{Z}^{*}$.

When defining the ppds, we vary several factors to test the performances of
the competing methods across different settings. Specifically, we consider
correlation between exposures (uncorrelated, correlated), skewness of the
exposure distribution (symmetric, skewed), and variability of the means of the
marginal exposure distributions (low, high). Each column of $\textbf{Z}^{*}$
is simulated based on these settings such that

$\displaystyle\begin{split}&\textbf{z}_{.j}^{*}=h\left(\textbf{w}_{.j}^{*}\right),\
j=1,\ldots,m=1000,\\\
&\textbf{w}_{.j}^{*}\stackrel{{\scriptstyle\text{iid}}}{{\sim}}\text{MVN}\left(\boldsymbol{\delta},\Sigma\right),\text{
and }\\\
&\delta_{i}|\tau^{2}\stackrel{{\scriptstyle\text{iid}}}{{\sim}}\text{N}\left(0,\tau^{2}\right),\
i=1,\ldots,n=250.\end{split}$ (5)

We standardize $\textbf{Z}^{*}$ by subtracting the mean of the entire matrix
from each element of $\textbf{Z}^{*}$ and dividing these centered values by
the standard deviation of the entire matrix.

Table 1: Simulation study settings for exposure posterior predictive distributions based on (5). Factor | Setting | Form
---|---|---
Associated | No | $\theta=0.00$
| Yes | $\theta=1.00$
Correlated | No | $\Sigma_{ij}=1\left(i=j\right)$
| Yes | $\Sigma_{ij}=\exp\left\\{2\ln\left(0.05\right)||\textbf{s}_{i}-\textbf{s}_{j}||\right\\}$
Skewed | No | $h\left(x\right)=x$
| Yes | $h\left(x\right)=\exp\left\\{x\right\\}$
$\delta_{i}$ Variance | Low | $\tau^{2}=0.10$
| High | $\tau^{2}=1.00$

For uncorrelated exposures, we define $\Sigma=I_{250}$ (i.e., the $250$ by
$250$ identity matrix) while for correlated data we define
$\Sigma_{ij}=\exp\left\\{-\phi||\textbf{s}_{i}-\textbf{s}_{j}||\right\\}$,
where $||.||$ represents the Euclidean distance between individuals $i$ and
$j$, and $\phi=-2\ln\left(0.05\right)$ which allows for the correlation
between two individuals to equal $0.05$ at a distance of $0.50$ (recall that
individuals are simulated within the unit square where the maximum possible
distance between two individuals is $\sqrt{2}$). For symmetric data,
$h\left(.\right)$ is defined as the identity function while
$h\left(.\right)=\exp\left\\{.\right\\}$ for skewed data. Finally,
$\tau^{2}=0.10$ represents low variability in the means of the marginal ppds
while $\tau^{2}=1.00$ represents high variability. As $\tau^{2}$ increases,
these means become further separated and we expect improved performance across
all methods since there will be less relative uncertainty in the exposure
distributions. Full details on all considered settings are given in Table 1
while randomly selected simulated datasets from each setting are presented in
Figure 1 and Figure S2 of the Supplement.

Figure 1: Simulated datasets from each setting in Table 1 where
$\tau^{2}=0.10$ (first column: uncorrelated, not skewed; second column:
correlated, not skewed; third column: uncorrelated, skewed; fourth column:
correlated, skewed). The first row shows the histogram of true exposures
(i.e., z), the second row shows the sample covariance matrix of the ppd
samples (i.e., $\textbf{Z}^{*}$), and the third row shows the histogram of the
median of each row of $\textbf{Z}^{*}$ (i.e., $\widehat{\text{z}}_{i}$).
Examples for $\tau^{2}=1.00$ are shown in Figure S2 of the Supplement.

For each simulated dataset, we generate a new set of spatial locations for the
$250$ individuals and unique set of exposures (i.e., $\textbf{Z}^{*}$ and z).
For every combination of factors in Table 1, we simulate 500 datasets and
analyze each one using the approaches described in Sections 2 and 3.

###  4.2 Data analysis

The health model, which is shared across all methods, is given as

$Y_{i}=\beta_{0}+\theta\text{z}_{i}+\epsilon_{i},\
\epsilon_{i}|\sigma^{2}_{\epsilon}\stackrel{{\scriptstyle\text{iid}}}{{\sim}}\text{N}\left(0,\sigma^{2}_{\epsilon}\right)$
(6)

with prior distributions $\sigma^{2}_{\epsilon}\sim\text{Inverse
Gamma}\left(0.01,0.01\right)$ (all methods),
$\beta_{0},\theta\sim\text{N}\left(0,100^{2}\right)$ (MIA, MVN, DU, UKDE,
MKDE), and $\beta_{0},\theta\sim\text{Flat}$ (Plug-in, MI). The decision to
use different prior distributions for the regression parameters across the
methods is made for computational reasons only, as the flat priors allow for
more efficient Monte Carlo posterior sampling for Plug-in and MI, saving
considerable computing time. We do not anticipate major changes in the results
due to these differences given that neither distribution contains informative
prior information about the parameters. For UKDE we use the method from
Sheather and Jones (1991) to estimate the bandwidth variables $h_{i}$, while
for MKDE we rely on Scott’s rule (Scott, 2015) and the sample covariance of
$\textbf{Z}^{*}$ to estimate the full bandwidth matrix variable $H$ (i.e.,
$H_{ij}=m^{-2/\left(n+4\right)}\widehat{\Sigma}_{ij}$).

From each method we collect $1,000$ posterior samples with which to make
inference. For the approaches where MCMC sampling is required, we run the
algorithms for $11,000$ iterations, removing the first $1,000$ as a burn-in
period, and thinning the remaining samples by a factor of $10$. For the Monte
Carlo sampling approaches, we directly obtain $1,000$ independent samples from
the joint posterior distributions. We estimate $\theta$ using the posterior
mean and quantify uncertainty using the 95% quantile-based, equal-tailed
credible interval.

We compare the ability of each method to estimate $\theta$ under the different
simulation settings by estimating and comparing the bias and mean squared
error (MSE) of the posterior mean, empirical coverage (EC) of the 95% credible
interval, and power/type I error rate. For reference, we also calculate these
quantities for the model in (6) where the true exposures are used - an
analysis not possible in practice.

###  4.3 Results

In Figure 2 we display boxplots of the $\theta$ estimates across all 500
analyses for each method and correlation/skewness setting for the
$\theta=1.00$ and $\tau^{2}=0.10$ scenario. Similar boxplots for the other
$\theta$ and $\tau^{2}$ combinations are shown in Figures S3-S5 of the
Supplement. In Table 2 we present the simulation study results, also from the
$\theta=1.00$ and $\tau^{2}=0.10$ scenario. Similar results are displayed in
Tables S1-S3 for the other $\theta$ and $\tau^{2}$ combinations. The
$\theta=1.00$, $\tau^{2}=0.10$ results suggest that MI and MIA perform very
similarly overall as expected, and that both struggle to estimate $\theta$
well across all settings. While DU outperforms MI and MIA, it tends to have
larger bias and MSE, and lower EC than some of the remaining approaches. These
methods may only be appropriate when the centers of the ppds are well
separated with low uncertainty such that each column of $\textbf{Z}^{*}$
begins to resemble the true exposures z.

Figure 2: Posterior mean estimates of $\theta$ for each method and
correlation/skewness setting across all 500 analyses for the $\theta=1.00$ and
$\tau^{2}=0.01$ scenario (first panel: uncorrelated, not skewed; second panel:
uncorrelated, skewed; third panel: correlated, not skewed; fourth panel:
correlated, skewed). The solid horizontal line represents the true value of
$\theta$.

When the ppds are not skewed, Plug-in generally performs well overall which is
an encouraging sign for past studies that have implemented this approach.
However, Plug-in struggles greatly when the ppds become skewed, leading to
elevated bias and MSE. MVN generally handles skewness better than Plug-in and
has improved performance when the ppds are correlated. However, it tends to
have lower EC, particularly when the ppds are independent and/or skewed.

MKDE performs similarly to DU across most metrics. It likely struggles due to
the difficulty in selecting the bandwidth matrix variable $H$ when the
dimension of the data is large. Overall, UKDE has the best balance of
performance compared to the other methods, especially when the ppds are
skewed. For $\tau^{2}=1.00$, several of the methods perform more similarly to
each other and we anticipate that this trend will continue as the ppds become
further separated and the relative error is effectively reduced.

When $\theta=0$, we see that all methods perform as expected while MI, MIA,
DU, and MDKE appear to be preferred in terms of bias and MSE primarily because
they produce estimates of $\theta$ near zero regardless of its true value.
Therefore, none of the considered methods are falsely identifying associations
at a higher rate than expected; an important conclusion for previous studies
that found evidence of significant associations between health and exposure
using these techniques.

Table 2: Simulation study results for $\theta=1.00$ and $\tau^{2}=0.10$. Estimates are presented with the range of standard errors for a group of estimates provided in parentheses. Bold entries represent the optimal estimate across an entire row (i.e., closest to zero for bias, smallest for mean squared error (MSE), closest to 95 for empirical coverage (EC), and largest for power). All results are multiplied by 100 for presentation purposes. | Settings |  | Methods
---|---|---|---
Metric | Corr | Skew |  | True | Plug-in | MI | MIA | DU | MVN | UKDE | MKDE
Bias | No | No |  | 0.15 | -1.03 | -90.95 | -90.96 | -71.26 | -33.68 | 5.71 | -71.65
| No | Yes |  | 0.30 | 60.77 | -94.07 | -94.07 | -61.23 | -30.73 | 4.59 | -59.06
| Yes | No |  | -0.05 | -2.58 | -89.66 | -89.64 | -50.49 | -2.25 | 2.69 | -61.78
| Yes | Yes |  | 0.00 | 62.14 | -90.83 | -90.83 | -35.10 | -15.02 | -4.46 | -73.23
|  |  |  | (0.27-0.40) | (1.21-3.24) | (0.11-0.21) | (0.12-0.21) | (0.31-1.80) | (0.61-2.33) | (0.99-2.25) | (0.31-3.57)
MSE | No | No |  | 0.36 | 7.29 | 82.79 | 82.80 | 51.25 | 14.17 | 5.22 | 51.80
| No | Yes |  | 0.46 | 89.38 | 88.56 | 88.57 | 40.90 | 24.01 | 7.89 | 38.99
| Yes | No |  | 0.47 | 7.78 | 80.48 | 80.45 | 30.14 | 1.89 | 5.29 | 60.53
| Yes | Yes |  | 0.82 | 116.24 | 82.72 | 82.72 | 28.44 | 29.36 | 25.41 | 117.26
|  |  |  | (0.02-0.07) | (0.47-10.62) | (0.21-0.37) | (0.21-0.37) | (0.44-1.29) | (0.55-2.83) | (0.29-2.28) | (0.44-10.38)
EC | No | No |  | 96.00 | 97.00 | 0.00 | 0.00 | 0.00 | 58.80 | 93.00 | 0.00
| No | Yes |  | 94.80 | 80.40 | 0.00 | 0.00 | 2.40 | 49.40 | 86.00 | 4.00
| Yes | No |  | 94.40 | 95.00 | 0.00 | 0.00 | 1.00 | 92.80 | 93.80 | 3.80
| Yes | Yes |  | 95.00 | 81.60 | 8.00 | 8.20 | 30.60 | 46.00 | 64.80 | 19.80
|  |  |  | (0.88-1.03) | (0.76-1.78) | (0.00-1.21) | (0.00-1.23) | (0.00-2.06) | (1.16-2.24) | (1.08-2.14) | (0.00-1.76)
Power | No | No |  | 100.00 | 92.60 | 0.00 | 0.00 | 54.80 | 88.40 | 94.80 | 52.40
| No | Yes |  | 100.00 | 70.80 | 0.00 | 0.00 | 67.20 | 54.60 | 98.20 | 70.00
| Yes | No |  | 100.00 | 93.60 | 0.00 | 0.00 | 93.00 | 99.80 | 94.60 | 75.20
| Yes | Yes |  | 100.00 | 70.20 | 0.00 | 0.00 | 85.60 | 78.80 | 89.00 | 63.60
|  |  |  | (0.00-0.00) | (1.09-2.05) | (0.00-0.00) | (0.00-0.00) | (1.14-2.23) | (0.20-2.23) | (0.59-1.40) | (1.93-2.23)

##  5 PM2.5 and stillbirth in New Jersey

Stillbirth is generally defined as the loss of a fetus or baby before or
during a delivery that occurs on or after 20 completed weeks of gestation
(Centers for Disease Control and Prevention, 2022). Recent literature reviews
suggest that exposure to ambient air pollution during pregnancy may be
associated with increased risk of stillbirth (Bekkar et al., 2020; Zhang et
al., 2021), though further studies are needed to better understand this
relationship. In this application, we examine the relationship between acute
exposure to PM2.5 in the days prior to delivery and risk of stillbirth using a
population-level time series approach in NJ, 2011-2015. We also investigate
the impact of the different methods for accounting for exposure uncertainty in
the health analysis, detailed in Sections 2 and 3, on the findings.

###  5.1 Data description

We obtain daily counts of fetal deaths and live birth across the three NJ
counties located in the New York (NY)-White Plains-Wayne, NY-NJ Metropolitan
Division (i.e., Bergen, Hudson, and Passaic) between 2011 and 2015 from the
Division of Family Health Services in the NJ Department of Health. Similar to
Warren et al. (2021), we only include singleton deaths/births and those with a
clinically estimated gestational age of $\geq$ 20 weeks. A fetal death
occurring on or after 20 weeks of gestation is defined as a stillbirth. In
Figure 3, we display the study area as well as the proportion of stillbirths
across time for this area.

Figure 3: The three New Jersey counties included in the study (left) and the
daily proportion of births that resulted in stillbirth in those counties in
2011-2015 (right).

Observed air pollution data and model-derived estimates are both obtained from
the United States Environmental Protection Agency (US EPA). Specifically, for
each day in 2002-2015 we access 24-hour average PM2.5 concentrations
(micrograms per cubic meter ($\mu$g/m3)) measured from all active monitors
located in NJ, NY, Delaware, and Pennsylvania from the US EPA’s Air Quality
System (AQS) (US EPA, 2022a). Model-derived daily estimates of 24-hour average
PM2.5 concentrations ($\mu$g/m3) from the Community Multiscale Air Quality
(CMAQ) model, a deterministic numerical air quality model, are obtained on a
12 by 12 kilometer grid across the same study area and time period (US EPA,
2022b). All data and estimates were downloaded from the US EPA’s Remote
Sensing Information Gateway website (US EPA, 2022c). Figure 4 displays the
locations of the AQS and CMAQ data from 2002-2015.

We also obtain daily estimates of the minimum and maximum temperatures at a
one kilometer resolution across the three NJ counties between 2011-2015 from
Daymet (Oak Ridge National Laboratory, 2022). The Daymet framework employs
statistical modeling techniques to produce spatially-temporally interpolated
temperature estimates using observed ground-based data as input. On each study
day, we average the estimates within the three counties to obtain a single
daily average minimum/maximum temperature estimate for the region.

Figure 4: Air pollution modeling/prediction study area where gray dots
represent locations of estimates from the Community Multiscale Air Quality
model and blue dots represent the location of the Air Quality System monitors
that were active at any point between 2002-2015 (left panel); posterior mean
predicted daily maximum 24-hour average PM2.5 from the model in (7-8) (middle
panel); histogram of posterior predictive exposure samples on August 1, 2011.

###  5.2 Stage 1: PM2.5 modeling and prediction

In the first stage of analysis, we use a hierarchical Bayesian framework for
modeling and predicting the daily PM2.5 concentrations collected from the AQS
using the closest CMAQ estimate as a predictor. Many of the air pollution
monitors in the study region are not active on a given day and the network of
monitors are only sparsely located across the study area (see Figure 4). As a
result, we use this model and the spatiotemporal completeness of the CMAQ
estimates to predict the AQS data at unobserved locations and days. We then
take the maximum of these predictions within the three NJ counties on each day
to estimate the daily maximum of the 24-hour PM2.5 concentrations for the
study area. This is used as the primary exposure of interest in the subsequent
stillbirth epidemiological analysis.

The model for the AQS PM2.5 data uses the closest CMAQ estimate as a predictor
within a flexible spatially- and temporally-varying regression coefficient
framework, similar to the original downscaling work of Berrocal et al.
(2010b), such that

$\ln\left\\{Z\left(\textbf{s},t\right)+0.01\right\\}=\eta_{0}\left(\textbf{s},t\right)+\eta_{1}\left(\textbf{s},t\right)\ln\left\\{\text{C}_{B\left(\textbf{s}\right),t}+0.01\right\\}+\epsilon\left(\textbf{s},t\right)$
(7)

where $Z\left(\textbf{s},t\right)$ is the AQS PM2.5 concentration measured at
the monitor located at s on day $t$ (i.e., $t=1$ is January 1, 2002 and
$t=5,113$ is December 31, 2015); $\text{C}_{B\left(\textbf{s}\right),t}$ is
the corresponding CMAQ estimate at the grid cell centroid located closest to s
(i.e., $B\left(\textbf{s}\right)$) on day $t$; and
$\epsilon\left(\textbf{s},t\right)|\sigma^{2}_{\epsilon}\stackrel{{\scriptstyle\text{iid}}}{{\sim}}\text{N}\left(0,\sigma^{2}_{\epsilon}\right)$.
We work on the log scale during modeling given that the PM2.5 concentrations
are $\geq 0$.

The spatially- and temporally-varying intercept and slope parameters are
represented by $\eta_{0}\left(\textbf{s},t\right)$ and
$\eta_{1}\left(\textbf{s},t\right)$, respectively, and allow the association
between the CMAQ estimates and AQS data to flexibly change across space and
time if appropriate. They are modeled as a function of spatial and temporal
covariates such that

$\eta_{j}\left(\textbf{s},t\right)=\mu_{j}+\text{lat}\left(\textbf{s}\right)\gamma_{1j}+\text{lon}\left(\textbf{s}\right)\gamma_{2j}+\left\\{\text{lat}\left(\textbf{s}\right)\text{lon}\left(\textbf{s}\right)\right\\}\gamma_{3j}+\sum_{k=1}^{4}\text{x}_{k}\left(t\right)\delta_{kj},\
j=0,1$ (8)

where $\text{x}_{k}\left(t\right)$ corresponds to the $k^{\text{th}}$ column
of the B-spline basis matrix for a polynomial spline with four degrees of
freedom (df) on day $t$, and
$\text{lat}\left(\textbf{s}\right)/\text{lon}\left(\textbf{s}\right)$ are the
latitude/longitude at spatial location s, respectively.

We complete the model by specifying non/weakly informative prior distributions
for the introduced model parameters. Specifically, we choose flat prior
distributions for all of the regression parameters (i.e.,
$\mu_{j},\gamma_{kj},\delta_{kj}$), and
$\sigma^{2}_{\epsilon}\sim\text{Inverse Gamma}\left(0.01,0.01\right)$;
resulting in an efficient closed-form Monte Carlo sampling algorithm for
obtaining samples from the joint posterior distribution of the model
parameters. We use it to collect $2,000$ independent posterior samples in
total.

Next, we use composition sampling (Tanner, 1996) to generate independent
posterior predictive samples of $Z\left(\textbf{s},t\right)$ at each of the
nine CMAQ grid cell locations within the three NJ counties on each day of the
study. On each day and for every joint ppd sample, we calculate the maximum of
the predictions across the three counties. In total, we obtain a $5,113$
(i.e., number of days) by $m=2,000$ matrix of maximum 24-hour average PM2.5
ppd samples, denoted by $\textbf{Z}^{*}$, and use it as the exposure for the
stillbirth epidemiological analysis.

###  5.3 Stage 2: Modeling stillbirth and PM2.5

Given $\textbf{Z}^{*}$ from Stage 1, we next turn to the stillbirth
epidemiological analysis. We model the total number of stillbirths occurring
across the three NJ counties on a specific day as a function of time-varying
predictors (i.e., day of week, long-term trend), meteorological variables
(i.e., maximum/minimum temperature), and lagged PM2.5 exposure. The model is
given as

$\displaystyle\begin{split}&Y_{t}|r,p_{t}\stackrel{{\scriptstyle\text{ind}}}{{\sim}}\text{Negative
Binomial}\left(r,p_{t}\right),\ t=1,\ldots,n,\\\
&\text{logit}\left(p_{t}\right)=\text{O}_{t}+\sum_{j=1}^{6}1\left\\{\text{dow}\left(t\right)=j\right\\}\alpha_{j}+\sum_{j=1}^{35}\text{x}_{1j}\left(t\right)\nu_{j}+\sum_{j=1}^{4}\text{x}_{2j}\left(w_{0t}\right)\zeta_{j}+\sum_{j=1}^{4}\text{x}_{3j}\left(w_{1t}\right)\kappa_{j}+\text{z}_{t-l}\theta\end{split}$
(9)

where $Y_{t}$ is the number of stillbirths occurring on day $t$ with $t=1$
corresponding to January 1, 2011 and $t=n=1,826$ to December 31, 2015; $p_{t}$
is the probability parameter that controls the magnitude of counts on day $t$;
$r>0$ represents the dispersion parameter with small values indicating
overdispersion in the data; $\text{O}_{t}$ is the offset variable representing
the log of the total number of births occurring on day $t$; $1\left(.\right)$
is an indicator function; $\text{dow}\left(t\right)$ is the day of week that
day $t$ occurred on with Saturday (i.e., $\text{dow}\left(t\right)=7$) serving
as the reference category; $\text{x}_{1j}\left(t\right)$,
$\text{x}_{2j}\left(w_{0t}\right)$, and $\text{x}_{3j}\left(w_{1t}\right)$ are
the $j^{\text{th}}$ columns of the B-spline basis matrices for a natural cubic
spline with 35, 4, and 4 df for study day, minimum temperature (i.e.,
$w_{0t}$), and maximum temperature (i.e., $w_{1t}$), respectively; and
$\text{z}_{t-l}$ is the true but unobserved maximum 24-hour average PM2.5
exposure $l$ days prior to day $t$.

We choose 35 df for the long-term time trend based on selecting 7 df for each
of the 5 study years as in Samet et al. (2000) while noting that Peng et al.
(2006) found reduced bias in effect estimation with more aggressive smoothing
in similar time series modeling. We consider daily lags from two to six days
(i.e., $l=2,\ldots,6$), similar to previous acute stillbirth and air pollution
modeling work (Faiz et al., 2013; Sarovar et al., 2020; Enebish et al., 2022)
and the estimated timing of 48 hours between fetal death and delivery (Gardosi
et al., 1998).

Using the model in (9) and $\textbf{Z}^{*}$, we test several of the existing
methods from Section 2 for propagating exposure uncertainty in the health
analysis along with the newly developed UKDE. MKDE is not considered given its
poor performance in simulation and long computing time for the large analysis
dataset. MI is not applied due to its lengthy run time (i.e., requires fitting
(9) in an MCMC framework $m=2,000$ times) and its overall similarity with MIA
in the simulation study results. We separately fit each method and exposure
lag ($l=2-6$) and make inference on $\theta$, the parameter that describes the
association between maximum PM2.5 exposure and stillbirth risk.

The prior distributions for the parameters in (9) are chosen as
$r\sim\text{Discrete Uniform}\left[1,100\right]$ and
$\alpha_{j},\nu_{j},\zeta_{j},\kappa_{j},\theta\stackrel{{\scriptstyle\text{iid}}}{{\sim}}\text{N}\left(0,100^{2}\right)$.
From all methods, we collect $10,000$ samples from the joint posterior
distributions after discarding the first $20,000$ prior to convergence and
thinning the remaining $200,000$ by a factor of $20$ to reduce posterior
autocorrelation. We assessed convergence by visually inspecting traceplots of
individual parameters and monitoring Geweke’s diagnostic (Geweke, 1991);
neither tool suggested any obvious signs of non-convergence across all model
fits.

###  5.4 Results

In Figure 4, we display daily predictions of the maximum PM2.5 exposures from
the three NJ counties in 2011-2015 along with a histogram of the ppd samples
on a single day (August 1, 2011; other days were similar). The level of
skewness in the ppd samples resembles that from the simulation study due to
the log transformation used in (7). In Figure S6 of the Supplement, we show a
scatterplot of these predictions and the observed AQS data (daily maximum of
the PM2.5 AQS concentrations across all of NJ). The plot shows that the model
is generally predicting well with respect to the observed data.

In Figure 5, we show results from the stillbirth analyses across all
considered methods and lags. Specifically, we present posterior inference
(i.e., posterior means and quantile-based equal-tailed credible intervals) for
$\exp\left\\{\theta\right\\}$ across all analyses, resulting in a relative
risk interpretation. Because we standardize $\textbf{Z}^{*}$ by subtracting
off the median and dividing by the interquartile range (IQR) prior to each
analysis, the estimates represent the relative risk of stillbirth for an IQR
increase in exposure during a given lag (IQRs ranged from 9.37-9.38 across all
lags). We applied a Bonferroni correction due to the separate lag models being
fit and present $99\%$ credible intervals in the figure.

Figure 5: Posterior mean and 99% credible interval plots for
$\exp\left\\{\theta\right\\}$ across the different models and daily lag
periods for the New Jersey three county stillbirth and maximum daily 24-hour
PM2.5 exposure analysis.

Generally the findings are in agreement with the simulation study results. MIA
and, to a lesser extent, DU tend to pull the point estimates towards the null
in comparison to the other approaches. Plug-in, MVN, and UKDE each suggest
that elevated ambient levels of maximum 24-hour average PM2.5 exposure three
days prior to delivery is associated with an increase in stillbirths.
Specifically for UKDE, an IQR increase in exposure of 9.4 $\mu$g/m3 three days
prior to delivery is associated with a $20.36\%$ increase in stillbirths
($99\%$ credible interval: $2.18-38.59$). While the differences between these
three methods is subtle in this application, UKDE has the shortest credible
interval followed by MVN and Plug-in, which may also be in agreement with the
improved MSE performance of UKDE observed in the simulation study.

As a sensitivity analysis, we repeated each of the stillbirth analyses while
randomly shuffling the order of the rows in $\textbf{Z}^{*}$. This mixing
breaks the temporal ordering of the exposures from the original analysis and
we expect to estimate a null signal in $\theta$ unless there are serious
confounding issues that the model fails to capture. The results shown in
Figure S7 of the Supplement show no significant associations across all
methods/lags, with point estimates near the null overall. This finding
provides further evidence that population-level PM2.5 may play an important
role in explaining stillbirth risk.

##  6 Discussion

In this work we developed “UKDE”, a new framework for exposure uncertainty
propagation in subsequent health outcome analyses, detailed its connection
with existing approaches, derived its closed-form MCMC full conditional
distributions, and created an R package for its implementation within several
common epidemiological analyses (KDExp; https://github.com/warrenjl/KDExp).
Existing methods for quantifying this uncertainty were detailed within a
unified framework, making comparing/contrasting the approaches more
accessible. In a simulation study, we showed that UKDE had improved
performance overall and particularly when the ppds were skewed. The
multivariate extension of UKDE, MKDE, was consistently outperformed by the
other methods, likely because of the difficulty associated with estimating
high dimensional densities using multivariate KDE.

Thijssen and Wessels (2020) evaluated several density estimation techniques
for performing sequential Bayesian inference, outside of the environmental
health setting, using posterior samples collected from a first stage analysis.
However, the focus of these analyses was on making inference on a low-
dimensional vector of parameters ($n=10$ was the maximum considered in the
study) that was shared across two or more datasets analyzed sequentially. The
goals and assumptions of this type of analysis differ from those in the
environmental health setting. In environmental health, the parameters shared
across both modeling stages (i.e., exposures) are typically high dimensional
(e.g., number of participants in a study), limiting the usefulness of some of
the presented approaches (e.g., MKDE). Additionally, sequential analysis is
primarily concerned with estimation of the set of parameters included in both
modeling stages/datasets, whereas in environmental health analyses the
emphasis is on correctly characterizing uncertainty in the exposures to
improve inference for parameters only included in the second stage health
outcome model (i.e., associations between exposure and health).

In our stillbirth data analysis in NJ, we found that elevated exposure to
maximum 24-hour average PM2.5 three days prior to delivery was associated with
elevated risk of stillbirth (relative risk: 1.20 (1.02, 1.39) for IQR
increase). In a recent study based in Mongolia, Enebish et al. (2022)
identified multiple critical daily PM2.5 exposure lags with respect to
elevated stillbirth risk, including three days prior to delivery (odds ratio
(OR): 1.28 (1.00, 1.62) for IQR increase). Another NJ-based analysis estimated
an OR with similar magnitude for an IQR increase in carbon monoxide exposure
two days before delivery (OR: 1.20 (1.05, 1.37)) (Faiz et al., 2013). Sarovar
et al. (2020) similarly identified a link between coarse particulate matter
exposure two days before delivery and increased odds of stillbirth, while
Mendola et al. (2017) saw similar links between ozone and stillbirth (with
similar magnitudes to our findings) on multiple days before delivery,
including day three. A recent meta-analysis found that stillbirth was
positively associated with an increase in ozone of 10 $\mu$g/m3 four days
before delivery, but found no significant short-term effects for PM2.5 in
pooled estimates (Zhang et al., 2021).

Future methods work in this area should focus on developing improved
techniques for incorporating correlation between high dimensional ppds. As
methods for investigating the impact of pollution mixtures on health are
becoming more common, extensions of the UKDE framework for more general
application in these settings is needed. In our work we showed that closed-
form MCMC full conditional updates were still available with additive exposure
models, but this does not necessarily hold for more complex interactions and
hierarchical structures used for the exposures. Additionally, extensions of
this work to accommodate critical exposure window identification/estimation
(e.g., Warren et al. (2012)) is also needed. Overall, UKDE is shown to be a
promising framework for characterizing exposure uncertainty in environmental
health studies, representing a more flexible hybrid approach between existing
methods, and can be implemented in the R package KDExp for several common
regression models.

## Acknowledgments

This research was supported by the National Institute of Environmental Health
Sciences (NIEHS) of the National Institutes of Health (NIH) under Award R01
NIEHS ES028346.

## References

  * Bekkar et al. (2020) Bekkar, B., S. Pacheco, R. Basu, R. Basu, and N. Denicola (2020).  Association of air pollution and heat exposure with preterm birth, low birth weight, and stillbirth in the US: a systematic review.  JAMA Network Open 3(6), 1–13.
  * Berrocal et al. (2010a) Berrocal, V. J., A. E. Gelfand, and D. M. Holland (2010a).  A bivariate space-time downscaler under space and time misalignment.  The Annals of Applied Statistics 4(4), 1942.
  * Berrocal et al. (2010b) Berrocal, V. J., A. E. Gelfand, and D. M. Holland (2010b).  A spatio-temporal downscaler for output from numerical models.  Journal of Agricultural, Biological, and Environmental Statistics 15(2), 176–197.
  * Berrocal et al. (2012) Berrocal, V. J., A. E. Gelfand, and D. M. Holland (2012).  Space-time data fusion under error in computer model output: an application to modeling air quality.  Biometrics 68(3), 837–848.
  * Blangiardo et al. (2016) Blangiardo, M., F. Finazzi, and M. Cameletti (2016).  Two-stage bayesian model to evaluate the effect of air pollution on chronic respiratory diseases using drug prescriptions.  Spatial and Spatio-temporal Epidemiology 18, 1–12.
  * Centers for Disease Control and Prevention (2022) Centers for Disease Control and Prevention (2022).  What is Stillbirth?  https://tools.cdc.gov/medialibrary/index.aspx#/media/id/218136.
  * Chang et al. (2011) Chang, H. H., R. D. Peng, and F. Dominici (2011).  Estimating the acute health effects of coarse particulate matter accounting for exposure measurement error.  Biostatistics 12(4), 637–652.
  * Enebish et al. (2022) Enebish, T., D. Warburton, R. Habre, C. Breton, and N. Tuvshindorj (2022).  The acute lag effects of elevated ambient air pollution on stillbirth risk in Ulaanbaatar, Mongolia.  medRxiv.
  * Faiz et al. (2013) Faiz, A. S., G. G. Rhoads, K. Demissie, Y. Lin, L. Kruse, and D. Q. Rich (2013).  Does ambient air pollution trigger stillbirth?  Epidemiology 24(4), 538–544.
  * Fuentes and Raftery (2005) Fuentes, M. and A. E. Raftery (2005).  Model evaluation and spatial interpolation by Bayesian combination of observations with outputs from numerical models.  Biometrics 61(1), 36–45.
  * Gardosi et al. (1998) Gardosi, J., T. Mul, M. Mongelli, and D. Fagan (1998).  Analysis of birthweight and gestational age in anteparturn stillbirths.  BJOG: An International Journal of Obstetrics & Gynaecology 105(5), 524–530.
  * Geweke (1991) Geweke, J. (1991).  Evaluating the Accuracy of Sampling-based Approaches to the Calculation of Posterior Moments, Volume 196.  Federal Reserve Bank of Minneapolis, Research Department Minneapolis, MN, USA.
  * Gryparis et al. (2009) Gryparis, A., C. J. Paciorek, A. Zeka, J. Schwartz, and B. A. Coull (2009).  Measurement error caused by spatial misalignment in environmental epidemiology.  Biostatistics 10(2), 258–274.
  * Guan et al. (2019) Guan, Y., B. J. Reich, J. A. Mulholland, and H. H. Chang (2019).  Multivariate spectral downscaling for pm2.5 species.  arXiv preprint arXiv:1909.03816.
  * Huang et al. (2018) Huang, G., D. Lee, and E. M. Scott (2018).  Multivariate space-time modelling of multiple air pollutants and their health effects accounting for exposure uncertainty.  Statistics in Medicine 37(7), 1134–1148.
  * Lee et al. (2017) Lee, D., S. Mukhopadhyay, A. Rushworth, and S. K. Sahu (2017).  A rigorous statistical framework for spatio-temporal pollution prediction and estimation of its long-term impact on health.  Biostatistics 18(2), 370–385.
  * Lee and Shaddick (2010) Lee, D. and G. Shaddick (2010).  Spatial modeling of air pollution in studies of its short-term health effects.  Biometrics 66(4), 1238–1246.
  * McMillan et al. (2010) McMillan, N. J., D. M. Holland, M. Morara, and J. Feng (2010).  Combining numerical model output and particulate data using bayesian space–time modeling.  Environmetrics 21(1), 48–65.
  * Mendola et al. (2017) Mendola, P., S. Ha, A. Z. Pollack, Y. Zhu, I. Seeni, S. S. Kim, S. Sherman, and D. Liu (2017).  Chronic and acute ozone exposure in the week prior to delivery is associated with the risk of stillbirth.  International Journal of Environmental Research and Public Health 14(7), 731.
  * Oak Ridge National Laboratory (2022) Oak Ridge National Laboratory (2022).  Daymet.  https://daymet.ornl.gov/.
  * Peng and Bell (2010) Peng, R. D. and M. L. Bell (2010).  Spatial misalignment in time series studies of air pollution and health data.  Biostatistics 11(4), 720–740.
  * Peng et al. (2006) Peng, R. D., F. Dominici, and T. A. Louis (2006).  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.
  * Polson et al. (2013) Polson, N. G., J. G. Scott, and J. Windle (2013).  Bayesian inference for logistic models using Pólya–Gamma latent variables.  Journal of the American Statistical Association 108(504), 1339–1349.
  * Reich et al. (2014) Reich, B. J., H. H. Chang, and K. M. Foley (2014).  A spectral method for spatial downscaling.  Biometrics 70(4), 932–942.
  * Samet et al. (2000) Samet, J. M., F. Dominici, F. C. Curriero, I. Coursac, and S. L. Zeger (2000).  Fine particulate air pollution and mortality in 20 us cities, 1987–1994.  New England Journal of Medicine 343(24), 1742–1749.
  * Sarovar et al. (2020) Sarovar, V., B. J. Malig, and R. Basu (2020).  A case-crossover study of short-term air pollution exposure and the risk of stillbirth in California, 1999–2009.  Environmental Research 191, 110103.
  * Scott (2015) Scott, D. W. (2015).  Multivariate Density Estimation: Theory, Practice, and Visualization.  John Wiley & Sons.
  * Sheather and Jones (1991) Sheather, S. J. and M. C. Jones (1991).  A reliable data-based bandwidth selection method for kernel density estimation.  Journal of the Royal Statistical Society: Series B (Methodological) 53(3), 683–690.
  * Szpiro and Paciorek (2013) Szpiro, A. A. and C. J. Paciorek (2013).  Measurement error in two-stage analyses, with application to air pollution epidemiology.  Environmetrics 24(8), 501–517.
  * Szpiro et al. (2011) Szpiro, A. A., L. Sheppard, and T. Lumley (2011).  Efficient measurement error correction with spatially misaligned data.  Biostatistics 12(4), 610–623.
  * Tanner (1996) Tanner, M. A. (1996).  Tools for Statistical Inference: Methods for the Exploration of Posterior Distributions and Likelihood Functions.  New York: Springer-Verlag.
  * Thijssen and Wessels (2020) Thijssen, B. and L. F. Wessels (2020, mar).  Approximating multivariate posterior distribution functions from Monte Carlo samples for sequential Bayesian inference.  PLoS ONE 15(3), e0230101.
  * US EPA (2022a) US EPA (2022a).  Air Quality System (AQS).  https://www.epa.gov/aqs.
  * US EPA (2022b) US EPA (2022b).  CMAQ Models.  https://www.epa.gov/cmaq/cmaq-models-0.
  * US EPA (2022c) US EPA (2022c).  RSIG-Related Downloadable Data Files.  https://www.epa.gov/hesc/rsig-related-downloadable-data-files.
  * Warren et al. (2021) Warren, J. L., H. H. Chang, L. K. Warren, M. J. Strickland, L. A. Darrow, and J. A. Mulholland (2021).  Critical window variable selection for mixtures: estimating the impact of multiple air pollutants on stillbirth.  arXiv preprint arXiv:2104.09730.
  * Warren et al. (2012) Warren, J. L., M. Fuentes, A. H. Herring, and P. H. Langlois (2012).  Spatial-temporal modeling of the association between air pollution exposure and preterm birth: identifying critical windows of exposure.  Biometrics 68(4), 1157–1167.
  * Warren et al. (2021) Warren, J. L., M. L. Miranda, J. L. Tootoo, C. E. Osgood, and M. L. Bell (2021, mar).  Spatial distributed lag data fusion for estimating ambient air pollution.  The Annals of Applied Statistics 15(1), 1–2.
  * Zhang et al. (2021) Zhang, H., X. Zhang, Q. Wang, Y. Xu, Y. Feng, Z. Yu, and C. Huang (2021).  Ambient air pollution and stillbirth: an updated systematic review and meta-analysis of epidemiological studies.  Environmental Pollution 278, 116752.
  * Zhou and Reiter (2010) Zhou, X. and J. P. Reiter (2010).  A note on Bayesian inference after multiple imputation.  The American Statistician 64(2), 159–163.