# Estimating a Causal Exposure Response Function with a Continuous Error-Prone
Exposure: A Study of Fine Particulate Matter and All-Cause Mortality

Kevin P. Josey  Department of Biostatistics, Harvard T.H. Chan School of
Public Health, Boston, MA  Corresponding Author<EMAIL_ADDRESS>Priyanka deSouza  Department of Urban and Regional Planning, University of
Colorado, Denver, CO  Xiao Wu  Department of Statistics, Stanford University,
Stanford, CA  Stanford Data Science, Stanford University, Stanford, CA
Danielle Braun  Department of Biostatistics, Harvard T.H. Chan School of
Public Health, Boston, MA  Department of Data Science, Dana-Farber Cancer
Institute, Boston, MA  Rachel Nethery  Department of Biostatistics, Harvard
T.H. Chan School of Public Health, Boston, MA

###### Abstract

Numerous studies have examined the associations between long-term exposure to
fine particulate matter (PM2.5) and adverse health outcomes. Recently, many of
these studies have begun to employ high-resolution predicted PM2.5
concentrations, which are subject to measurement error. Previous approaches
for exposure measurement error correction have either been applied in non-
causal settings or have only considered a categorical exposure. Moreover, most
procedures have failed to account for uncertainty induced by error correction
when fitting an exposure-response function (ERF). To remedy these
deficiencies, we develop a multiple imputation framework that combines
regression calibration and Bayesian techniques to estimate a causal ERF. We
demonstrate how the output of the measurement error correction steps can be
seamlessly integrated into a Bayesian additive regression trees (BART)
estimator of the causal ERF. We also demonstrate how locally-weighted
smoothing of the posterior samples from BART can be used to create a more
accurate ERF estimate. Our proposed approach also properly propagates the
exposure measurement error uncertainty to yield accurate standard error
estimates. We assess the robustness of our proposed approach in an extensive
simulation study. We then apply our methodology to estimate the effects of
PM2.5 on all-cause mortality among Medicare enrollees in New England from
2000-2012.

##  1 Introduction

Methods for conducting causal inference with continuous exposures have gained
significant traction in recent years. The goal of these methods is to estimate
a causal exposure response function (ERF) that is free of confounding bias
(Kennedy et al.,, 2017; Wu et al.,, 2018). Estimators of an ERF typically rely
on estimates of either the generalized propensity score (GPS), a marginalized
model of the outcome process, or some combination of the two to adjust for
confounding. However, most of these methods make the implicit assumption that
the exposure is measured without error, which is implausible in many
observational study settings. For example, in modern environmental health
research, air pollution exposure measurements are often derived from model-
based predictions of pollutant concentrations rather than the exact pollutant
concentrations experienced by an individual. Moreover, because many studies
record only areal measures of residential locations such as ZIP-codes, cities,
or counties, exposure measures often represent aggregated predictions across
these areas. Using aggregated exposures assumes that concentrations are
homogeneous within areas and experienced by all individuals residing in the
area. As such, exposure measurement error is prevalent in many air pollution
epidemiological studies (Kioumourtzoglou et al.,, 2014). Using an error-prone
exposure (EPE) in place of the true exposure violates standard causal
inference assumptions and may result in biased ERF estimates. Due to the
policy-relevance of air pollution ERFs, and particularly those produced using
causal inference techniques, invalid inferences caused by EPEs could have
severe consequences.

While measurement error has been studied extensively outside of causal
inference settings (Carroll et al.,, 2006; Cole et al.,, 2006), accounting for
EPEs in causal inference is a relatively new endeavor and is mainly confined
to scenarios with binary and categorical exposures (Lewbel,, 2007; Braun et
al.,, 2017; Wu et al.,, 2019), or cases of measurement error in a confounder
instead of the exposure (Lenis et al.,, 2017; Webb-Vargas et al.,, 2017).
Beyond the issues encountered when using an EPE in a typical generalized
linear model setting, accommodating an EPE in a causal model presents
additional challenges in conjunction with resolving confounding bias (Braun et
al.,, 2017), many of which remain unaddressed. Moreover, propagating
measurement error-related uncertainty into ERF estimates is paramount for
proper inference in not just causal settings, but more broadly for many
environmental epidemiological studies.

We develop a multiple imputation framework for estimating causal ERFs which
incorporates corrections for various types of exposure measurement error
commonplace in environmental epidemiology. This approach jointly samples from
three sub-models: 1) a model of the EPE, which imputes the true exposures
using information from covariates and validation data when available, 2) a GPS
model, which improves the accuracy of the imputations, and 3) an outcome
model, which estimates the ERF using the imputed exposures while adjusting for
confounders. Markov Chain Monte Carlo (MCMC) methods are used to sample from
the posterior distributions of the unobserved true exposures and the expected
outcome conditioned on the exposures and confounders. A Bayesian additive
regression trees (BART) model is specified for the outcome sub-model to
capture complex, non-linear relationships between the exposure, the
covariates, and the outcome. We also show that implementing a further
smoothing step on the BART-estimated ERF using local regression techniques
provides accurate estimates of the ERF with fewer posterior samples than a
classic Bayesian approach, all while adequately propagating the measurement
error uncertainty. The latter smoothing step transitions our approach from a
Bayesian estimator into one better characterized as a multiple imputation
estimator, while at the same time producing smoother estimates of the ERF than
the notoriously noisy BART output (Nethery et al.,, 2021).

The remainder of this article is structured as follows. We describe the
motivating data example in Section 1.1. In Section 2 we introduce notation,
define the measurement error sources, and identify the causal ERF model.
Section 3 outlines the procedure for correcting the attenuation bias created
by measurement error using multiple imputation methods, and we discuss several
important caveats and limitations to our proposed methodology. Section 4
contains a simulation study. Section 5 applies the proposed method to create a
measurement error-corrected causal ERF for long-term fine particulate matter
(PM2.5) exposure and all-cause mortality in the Medicare population in New
England, using the data described in Section 1.1. We conclude with a
discussion in Section 6.

###  1.1 Motivating Example

PM2.5 is a well-studied air pollutant known to adversely impact numerous
health outcomes (Hajat et al.,, 2002; Dominici et al.,, 2006; Brook et al.,,
2010; Zhu et al.,, 2017), including all-cause mortality (Anderson,, 2009; Di
et al.,, 2017; Wu et al.,, 2020), cardiovascular disease (Dominici et al.,,
2006; Zanobetti et al.,, 2009; Pope III et al.,, 2015; Yazdi et al.,, 2019;
Danesh Yazdi et al.,, 2021), and pulmonary/respiratory diseases (Dominici et
al.,, 2006; Zanobetti et al.,, 2009; Rhee et al.,, 2019). In one of these
studies, using a cohort of Medicare enrollees 2000-2012, Wu et al., (2019)
implemented various causal inference techniques to assess whether long-term
PM2.5 exposure increases the risk of all-cause mortality among older Americans
($>65$ years of age). While Wu et al., (2019) employed a regression
calibration approach to resolve parts of the exposure measurement error
present in this application, there were other error components that were left
uncorrected. Moreover, a simple regression calibration approach may fail to
adequately propagate measurement error-related uncertainty into the causal
effect estimates. Additionally, Wu et al., (2019) considered exposure
categories instead of examining the ERF across a continuum of exposures.
Categorization of exposures implicitly assumes that participants in the same
category are exposed to the same exposure level, which may induce yet another
source of measurement error. We seek to investigate the same scientific
question as Wu et al., (2019) using the same data but employing a novel
analytic approach that: 1) corrects for additional sources of measurement
error in the PM2.5 exposures and 2) models the continuous ERF as opposed to
categorizing the exposure.

We utilize annual average PM2.5 predictions for 1km$\times$1km grid-years
across the US from Di et al., (2016), which were generated via a neural
network that combines information from satellite, ground monitor, and land use
data. These exposure predictions have been widely used in epidemiological
studies (Di et al.,, 2017; Rhee et al.,, 2019; Wu et al.,, 2019; Yazdi et
al.,, 2019; Wu et al.,, 2020; Danesh Yazdi et al.,, 2021). Model assessments
indicate that the accuracy of the predictions vary across regions (Di et al.,,
2016). Because these exposure predictions are error-prone, any inference drawn
from them without correction may be biased. To control for this source of
measurement error, we use the same ground-monitor PM2.5 data that were used to
generate the grid-predicted PM2.5 measurements in the first place. We treat
the measurements from these monitors as the gold standard measurements of
PM2.5 exposure (i.e. error-free measurements). Our model uses the monitored
data as a validation set to re-calibrate the error-prone PM2.5 measurements.
PM2.5 ground monitor locations in New England are overlaid on a ZIP-code map
in Figure 1.

Another limitation of the data is that the residential locations of the
Medicare enrollees can only be mapped to ZIP-codes (Wu et al.,, 2020). Thus,
following common practice in this setting, we use ZIP-code-years as our units
of analysis and all-cause mortality rates within ZIP-code years as our
outcome. This creates a misalignment between the PM2.5 exposure predictions,
which occur on 1km grids, and the ZIP-codes to which we wish to assign PM2.5
exposures.

We also leverage data on demographic and socioeconomic status (SES) factors,
as well as land-use variables which are used to aid in measurement error
correction. The following potential confounder variables are collected from
the US Census Bureau at the ZIP-code-year level: population density, median
household income, percent population and households below poverty level,
racial/ethnic composition, and distribution of educational attainment
(college, some college, high school, not completed high school). Additionally,
individual level information about the Medicare enrollees, including their
sex, age, race, and Medicaid eligibility status (a proxy for low-income) are
obtained from the Medicare record and summarized at the ZIP-code-year level.
The land-use variables, collected at the grid-year level, are: surface
temperature, accumulated precipitation, radiation flux, accumulated total
evaporation, heat flux, precipitation rate, humidity, snow cover, cloud cover,
and wind speed.

##  2 Preliminaries

###  2.1 Notation and Measurement Error

Our data involve measurements at two different spatial scales: ZIP-code-years
(outcome, covariates) and 1km$\times$1km grid-years (exposure, land-use
variables). In general, we refer to grid-years as cells, and cells are assumed
to be nested within ZIP-code-years, which we refer to as clusters. Let
$\mathbf{X}_{i}\in\mathcal{X}$ be a vector of covariates for cluster $i$
($i=1,2,\ldots,n$) and $Y_{i}$ denote the outcome variable in that cluster.
While the methodology we present easily generalizes to other outcome
distributions, in our applied example the $Y_{i}$ are counts (e.g. number of
deaths in ZIP-code-year $i$) with $N_{i}$ representing an offset (e.g. number
of at-risk person-years in ZIP-code-year $i$). The combination of these two
measures allows us to define the rate $\bar{Y}_{i}\equiv Y_{i}/N_{i}$. We
denote the cumulative density functions for the $\mathbf{X}_{i}$ with
$P(\mathbf{x})$ with the associated empirical density function denoted by
$\mathbb{P}(\mathbf{x})$.

We assume there is a true exposure both on the cluster level and on the cell
level. To enable alignment with the outcome data at the cluster level, the
true cluster level exposure, denoted $A_{i}\in\mathcal{A}$, is the latent
variable we are most concerned with accurately predicting. We assume that the
$A_{i}$ are continuous but unobserved for all clusters. The true cell level
exposures within cluster $i$ are indexed by $j=1,2,\ldots,M_{i}$ and are
denoted as $S_{ij}$. The total number of cells is written as
$m\equiv\sum_{i=1}^{n}M_{i}$. We assume, as in our example, that the $S_{ij}$
will only be observed for a subset of cells – the grids containing monitors –
which are indexed by $(i,j)\in\mathcal{S}$. We also have for every cluster and
cell an observed model-predicted value of $S_{ij}$, denoted as
$\tilde{S}_{ij}$. These predictions are the error-prone exposures. The
aggregate value of these EPEs to the cluster level will be denoted with
$\tilde{Z}_{i}=M_{i}^{-1}\sum_{j=1}^{M_{i}}\tilde{S}_{ij}$. Figure 2
illustrates the relationship between $A_{i}$, $S_{ij}$, and $\tilde{S}_{ij}$
in the context of our motivating application. We also allow for the existence
of a set of cell level features, denoted by the vector $\mathbf{W}_{ij}$, that
might inform us about the relationship between $S_{ij}$ and $\tilde{S}_{ij}$.
In our applied example, the $\mathbf{W}_{ij}$ are land-use variables measured
at the grid-year level. $\mathbf{X}_{i}$ may contain aggregated and/or
transformed values of the variables contained in $\mathbf{W}_{ij}$.

Finally, the conditional expected value of any given variable $D$ is denoted
by $\mu_{D}(\cdot)\equiv\mathbb{E}(D|\cdot)$. The conditional probability
density function of $D$ is denoted with $p(D|\cdot)$. The generalized
propensity score (GPS) is the following particular conditional probability
density function: $p(A_{i}|\mathbf{X}_{i})$.

###  2.2 Measurement Error Model

We formulate the measurement error in this context as a special case of
nondifferential classical error. As we noted in Section 2.1, we must contend
with the problem that $S_{ij}$ is rarely observed, and in most areas only the
EPE measurements, $\tilde{S}_{ij}$, are available. Moreover, conducting a
cluster level analysis relating $A_{i}$ to $\bar{Y}_{i}$ requires measurements
of $A_{i}$, not $S_{ij}$. Therefore, we must consider approaches for
summarizing $S_{ij}$ while accounting for the possible confounding influence
of $\mathbf{X}_{i}$ and the varying sizes of $M_{i}$, and ensuring the
uncertainty from this process is propagated into the final estimator.

Throughout this paper, we will assume that $\mathbb{E}(S_{ij}|A_{i})=A_{i}$,
which can be framed as treating the true cell level exposures $S_{ij}$ as
replicate measures of the corresponding cluster level exposure $A_{i}$, each
measured with some amount of error. This ‘replicate measures’
conceptualization is commonly used in the presentation of classical
measurement error methods. An additional dimension of complexity is added by
the need to rely on cell level exposure predictions, $\tilde{S}_{ij}$. This
leads us to formulate the measurement error structure as a variant of the
typical nondifferential classical measurement error, which can be decomposed
into two components as follows:

$U_{ij}=\underbrace{\overbrace{S_{ij}-\tilde{S}_{ij}}^{\text{Prediction}}+\overbrace{\tilde{S}_{ij}-A_{i}}^{\text{Aggregation}}}_{\text{Classical}}.$
(1)

Following conventions in the classical measurement error literature, we make
the assumption that the measurement error $U_{ij}$ is homoscedastic and
conditionally independent such that:

######  Assumption 1 (Conditionally Independent Measurement Error).

We assume $\left(U_{ij}\mathchoice{\mathrel{\hbox
to0.0pt{$\displaystyle\perp$\hss}\mkern
2.0mu{\displaystyle\perp}}}{\mathrel{\hbox
to0.0pt{$\textstyle\perp$\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptstyle\perp$\hss}\mkern
2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptscriptstyle\perp$\hss}\mkern
2.0mu{\scriptscriptstyle\perp}}}U_{ij^{\prime}}\right)|A_{i},\mathbf{W}_{ij}$
for two cells $j\neq j^{\prime}$ in a given cluster $i=1,2,\ldots,n$;

######  Assumption 2 (Homoscedastic Measurement Error).

For all $j=1,2,\ldots,M_{i}$ contained within every $i=1,2,\ldots,n$, we
assume the measurement error conditional on $A_{i}$ has constant variance;
$\mathbb{V}(U_{ij}|A_{i})\equiv\omega^{2}$.

The prediction error component in (1) can be characterized as a Berkson error
term (Haber et al.,, 2020). The aggregation error, on the other hand, is a
type of classical measurement error, implying the composite measurement error
term $U_{ij}$ is also a type of classical error. To help cement this
understanding, ignore for a moment the prediction error term in $U_{ij}$. The
goal of regression calibration is to replace classical measurement error with
Berkson error. It has been noted in the regression calibration literature that
Berkson error resulting from regression calibration often yields less bias
than its classical error counterpart (Carroll et al.,, 2006). Our method,
presented in the following sections, is developed in the same spirit,
replacing $U_{ij}$ with a Berkson error term via multiple imputation. Our
method also allows the user to replace the prediction error component of
$U_{ij}$, if desired, but doing so essentially substitutes the original
Berkson prediction error with another Berkson error term. That said, re-
calibrating the predictions $\tilde{S}_{ij}$ can still be useful if we suspect
$\mathbb{E}[\tilde{S}_{ij}|A_{i}]\neq A_{i}$. In this latter case, if we
assume $\mathbb{E}[S_{ij}|A_{i}]=A_{i}$ for all $(i,j)\in\mathcal{S}$, then we
can use validation data to find an unbiased predictor for $S_{ij}$, denoted by
$\hat{S}_{ij}\equiv\hat{\mu}_{S}(\tilde{S}_{ij},W_{ij})$, which implies
$\mathbb{E}[\hat{S}_{ij}|A_{i}]=A_{i}$.

As we will show in Section 3, corrections for both the classical measurement
error between $S_{ij}$ and $A_{i}$ and any potential bias found in the
predictions $\tilde{S}_{ij}$ can be accommodated under the proposed multiple
imputation framework. The imputations of the true exposures are drawn from an
MCMC sampler using the full conditional likelihood in (2).

###  2.3 Assumptions and Identification

The methods we employ operate within the Neyman-Rubin causal model modified
for continuous exposures Rubin, (1974). Let $\bar{Y}_{i}(a)$ be the outcome
that would occur in cluster $i$ if, possibly contrary to what is observed, it
had received exposure level $A_{i}=a$. We refer to the $\bar{Y}_{i}(a)$ for
any $a\in\mathcal{A}$ as the potential outcomes. For some exposure level $a$
the objective analysis is to estimate the marginal ERF,
$\theta(a)\equiv\mathbb{E}[\bar{Y}_{i}(a)]$. Unlike in the binary or
categorical exposure setting, where there is only a finite number of
unrealized potential outcomes, with a continuous exposure there is an infinite
number of unrealized potential outcomes for every unit in the study. Despite
this perceived challenge, it is relatively straightforward to translate the
assumptions intended for a categorical exposure under the Neyman-Rubin causal
model to a continuous exposure setting necessary to identify and estimate
$\theta(a)$.

To start, we invoke the stable unit treatment value assumption which consists
of two conditions: 1) consistency and 2) no interference. Consistency refers
to the notion that $\bar{Y}_{i}=\bar{Y}_{i}(A_{i})$, i.e. each unit’s observed
outcome is equal to its potential outcome evaluated at the true exposure
experienced. We are careful to define this assumption in the setup to our
problem; recognizing that $A_{i}$ is unobserved. While unobserved, the true
exposure $A_{i}$ still exists, so consistency should hold. On the other hand,
if we naïvely substitute $\tilde{Z}_{i}$ for $A_{i}$, then consistency is
likely violated with respect to the EPE unless $\tilde{Z}_{i}=A_{i}$ implying
$\bar{Y}_{i}(\tilde{Z}_{i})=\bar{Y}_{i}(A_{i})$ for all $i$. As this
assumption is unlikely to hold, it is safe to assume consistency is violated
if measurement error-agnostic approaches are taken. No interference refers to
the assumption that the exposure of one unit does not affect the potential
outcomes of another unit. This condition is difficult to uphold in many
environmental epidemiology studies, including our own applied example. Major
air pollution emissions sources can affect broad areas (e.g., many ZIP-codes),
which may lead to instances where the exposure assignment of one ZIP-code
affects the outcome of another (Papadogeorgou et al.,, 2019). We acknowledge
this limitation, yet addressing it is outside the scope of this work.

The fundamental problem in causal inference with any type of exposure is that
$\theta(a)$ is not intrinsically estimable since $\bar{Y}_{i}(a)$ is not
observed for every $a\in\mathcal{A}$. However, $\theta(a)$ can still be
estimated under the strongly ignorable exposure assumption. With the Neyman-
Rubin causal model, this requires the following two assumptions (Kennedy et
al.,, 2017):

######  Assumption 3 (Overlap/positivity).

$p(A_{i}|\mathbf{X}_{i})>p_{\text{min}}>0$.

######  Assumption 4 (Weak unconfoundedness).

$\left[\bar{Y}_{i}(a)\mathchoice{\mathrel{\hbox
to0.0pt{$\displaystyle\perp$\hss}\mkern
2.0mu{\displaystyle\perp}}}{\mathrel{\hbox
to0.0pt{$\textstyle\perp$\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptstyle\perp$\hss}\mkern
2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptscriptstyle\perp$\hss}\mkern
2.0mu{\scriptscriptstyle\perp}}}A_{i}\right]|\mathbf{X}_{i}$ for all
$a\in\mathcal{A}$.

The overlap assumption requires the exposure assignment mechanism to be non-
deterministic when conditioned on the covariates. In other words, the
probability that a unit is exposed to any level of the exposure along the
support is always greater than zero (or $p_{\text{min}}$). Weak
unconfoundedness states that the potential outcome evaluated at $a$ does not
depend on the true exposure $A_{i}$ when adjusted by the set of observed
confounding variables. This implies there can be no unmeasured confounding. If
all confounders for the relationship between $\bar{Y}_{i}$ and $A_{i}$ are
measured, the latter assumption should hold even when $A_{i}$ exists but is
not observed. When these assumptions hold, it is valid to draw causal
conclusions either with experimental or observational data and a continuous
exposure (Gill and Robins,, 2001), and assuming $A_{i}$ is known, the true
causal ERF can be identified as

$\mathbb{E}\left\\{\mathbb{E}\left[\bar{Y}_{i}\middle|A_{i}=a,\mathbf{X}_{i}\right]\right\\}=\mathbb{E}\left\\{\mathbb{E}\left[\bar{Y}_{i}(a)\middle|\mathbf{X}_{i}\right]\right\\}=\mathbb{E}\left[\bar{Y}_{i}(a)\right]=\theta(a)$

##  3 A Multiple Imputation Approach

###  3.1 Sampling the Exposure and Nuisance Parameters

We introduce our method by first specifying the joint likelihood function for
$\bar{Y}_{i}$, $A_{i}$, and $S_{ij}$, which incorporates the components needed
to address the measurement error described by (1). Supposing that every piece
of data were observed, we have:

$\prod_{i=1}^{n}\left\\{\underbrace{p_{\bar{Y}}(\bar{Y}_{i}|A_{i},\mathbf{X}_{i})}_{\text{Outcome
Model}}\times\underbrace{p_{A}(A_{i}|\mathbf{X}_{i},\phi_{i})\times
p_{\phi}(\phi_{i}|\mathbf{V})}_{\text{GPS
Model}}\times\prod_{j=1}^{M_{i}}\underbrace{p_{S}(S_{ij}|A_{i})}_{\text{EPE
Model}}\right\\}.$ (2)

From a modelling perspective, we assume that the latent exposure variables are
conditionally independent and approximately Normal in distribution with
$A_{i}|\mathbf{X}_{i},\phi_{i}\sim\mathcal{N}\left[\mu_{A}(\mathbf{X}_{i},\phi_{i}),\sigma^{2}\right]$
and $S_{ij}|A_{i}\sim\mathcal{N}(A_{i},\omega^{2})$, in accordance with
Assumptions 1 and 2. A random effect $\phi_{i}$ is included in the model for
$A_{i}$ to control for spatial autocorrelation between the cluster level
exposures (Lee,, 2013). We assume
$\boldsymbol{\phi}\equiv(\phi_{1},\phi_{2},\ldots,\phi_{n})^{T}\sim\mathcal{N}\left[\mathbf{0},\nu^{2}\mathbf{Q}^{-1}(\mathbf{V};\rho)\right]$
where $\mathbf{Q}(\mathbf{V};\rho)$ is a precision matrix controlling the
auto-correlation structure as proposed by Leroux et al., (2000) given a binary
adjacency matrix $\mathbf{V}$ and a hyperparameter $\rho$ controlling the
correlation between adjacent clusters. If we assume the $A_{i}$ are
independent, then we can set $\phi_{i}=0$ and
$p_{\phi}(\phi_{i}|\mathbf{V})=1$ for all $i$. Because $S_{ij}$ is missing, we
can substitute $\tilde{S}_{ij}$ for $S_{ij}$ if we suppose
$\mathbb{E}[\tilde{S}_{ij}|A_{i}]=A_{i}$. However, when given a set of
validation data (i.e. a subset where $S_{ij}$ is observed), a more
conservative approach is to instead assume $\mathbb{E}[S_{ij}|A_{i}]=A_{i}$
and draw a sequence of posterior predictions for $S_{ij}$ for all
$j=1,2,\ldots,M_{i}$ and $i=1,2,\ldots,n$. To do this, we can append the model
$\prod_{(i,j)\in\mathcal{S}}p_{S}(S_{ij}|\tilde{S}_{ij},\mathbf{W}_{ij})$ to
(2) and, supposing that
$S_{ij}|\tilde{S}_{ij},\mathbf{W}_{ij}\sim\mathcal{N}\left[\mu_{S}(\tilde{S}_{ij},\mathbf{W}_{ij}),\tau^{2}\right]$,
sample values of $S_{ij}$. Note that this addition overspecifies the full-
conditional likelihood in (2), so doing so is only possible when validation
data are available.

We use the above likelihood throughout this manuscript to describe how to
sample posterior draws of $A_{i}$ contemporaneously with posterior predictions
for $\mu_{\bar{Y}}(A_{i},\mathbf{X}_{i})$. We will show in the next section
how these posterior samples can be smoothed to create an estimate of the ERF.
Relative to the standard regression calibration approach where only a single
imputation of $A_{i}$ is used in fitting the outcome model, using a set of
posterior samples (i.e. multiple imputations) of the error-corrected exposures
should better propagate uncertainty attributable to the exposure measurement
error into the outcome model.

We use a fully Bayesian joint modeling approach for the measurement error,
GPS, and outcome models. Parameter values are sampled from the posterior
distribution either by Gibbs or Metropolis-Hastings sampling steps, with an
added intermediate step to generate posterior predictive samples of the latent
variables $A_{i}$ and, if necessary, $S_{ij}$. We refer to this intermediate
step as the imputation stage, whereas the steps involving sampling the
parameters are referred to as the analysis stage. See the full details of the
sampling algorithm in Supplemental Section S1.

We must carefully consider the form of $\mu_{A}(\mathbf{X}_{i},\phi_{i})$ and
$\mu_{S}(\tilde{S}_{ij},\mathbf{W}_{ij})$ to properly address bias in the EPE
model while also accounting for confounding. In our simulation study and data
analysis in Sections 4 and 5, we assume linear forms for
$\mu_{A}(\mathbf{X}_{i},\phi_{i})=\phi_{i}+\mathbf{X}^{T}_{i}\boldsymbol{\beta}$
and
$\mu_{S}(\tilde{S}_{ij},\mathbf{W}_{ij})=\left(\tilde{S}_{ij},\mathbf{W}^{T}_{ij}\right)\boldsymbol{\alpha}$.
To better capture nonlinear associations and better avoid issues with model
misspecification, splines or Gaussian processes could alternatively be
specified (Antonelli et al.,, 2020; Ren et al.,, 2021).

If we assume that $\bar{Y}_{i}$ is a linear function of $A_{i}$ and
$\mathbf{X}_{i}$, then the parameters determining
$p_{\bar{Y}}(\bar{Y}_{i}|A_{i}\mathbf{X}_{i})$ can be drawn using a Gibbs
sampling step, assuming the $N_{i}$ are fixed and known. However, we prefer to
use a data-driven model for $p_{\bar{Y}}(\bar{Y}_{i}|A_{i},\mathbf{X}_{i})$
and therefore $\mu_{\bar{Y}}(A_{i},\mathbf{X}_{i})$, because this model is
perhaps the most essential component in (2) for finding accurate estimates of
$\theta(a)$. To this end, we employ an iteratively updated BART model (Chipman
et al.,, 2010). Letting $t=1,2,\ldots,T$ index the trees for a given iteration
of the MCMC, the BART model assumes

$\bar{Y}_{i}=\mu_{\bar{Y}}(A_{i},\mathbf{X}_{i})+\epsilon_{i}\approx\sum_{t=1}^{T}g\left(A_{i},\mathbf{X}_{i};\mathcal{F}_{t},\mathcal{G}_{t}\right)+\epsilon_{i}.$
(3)

Here, $g(\cdot)$ is a function that bins observations into groups with binary
trees formed by the rules contained in $\mathcal{F}_{t}$, and node means
characterized by the set $\mathcal{G}_{t}$. A BART model differs from other
regression tree ensembling methods because of the priors placed on
$\mathcal{F}_{t}$ and $\mathcal{G}_{t}|\mathcal{F}_{t}$, which are used to
sample values from the approximate posterior distribution of
$\mu_{\bar{Y}}(A_{i},\mathbf{X}_{i})$. Posterior samples are denoted by
superscript $(k)$, $k=1,2,\ldots,K$, e.g.,
$\mu^{(k)}_{\bar{Y}}(A^{(k)}_{i},\mathbf{X}_{i})\approx\sum_{t=1}^{T}g\left(A^{(k)}_{i},\mathbf{X}_{i};\mathcal{F}^{(k)}_{t},\mathcal{G}^{(k)}_{t}\right)$.
In each iteration of our MCMC sampler, the posterior samples of the BART
parameters are drawn conditional on the current posterior sample of the
exposures $A_{i}^{(k)}$. The error term is assigned the distribution
$\epsilon_{i}\sim\mathcal{N}(0,\psi^{2}/N_{i})$.

For model-fitting, we assign conjugate inverse-gamma prior distributions to
the variance parameters,
$\omega^{2},\tau^{2},\sigma^{2},\nu^{2},\psi^{2}\sim\mathcal{IG}(0.001,0.001)$,
and set $\rho\sim\mathcal{U}(0,1)$. The regression coefficients for
$\mu_{S}(\cdot)$, $\mu_{A}(\cdot)$, are each assigned independent Gaussian
priors with default values
$\boldsymbol{\alpha},\boldsymbol{\beta}\sim\mathcal{N}(\mathbf{0},10^{6}\mathbf{I})$
where $\mathbf{0}$ is a vector with every entry equal to zero and $\mathbf{I}$
is an identity matrix (with appropriate dimensions). Additional sampling
details are provided in Supplemental Section S1. For an example of how to
implement a generalized linear model for
$\bar{Y}_{i}|A^{(k)}_{i},\mathbf{X}_{i}$ instead of the BART implementation in
(3), please refer to the simulation experiment in Supplemental Section S4.

###  3.2 Smoothing BART Output

While predictions from the BART model (e.g., means of posterior predictive
samples) could be used directly to form an estimate of the ERF, these samples
do not typically form a smooth function of $\theta(a)$ over $a\in\mathcal{A}$.
Smooth ERFs are believed to be most biologically plausible relationship in
many epidemiological and biomedical applications, including the effect of
PM2.5 on all-cause mortality, and are more desirable for identifying causal
relationships (Kim et al.,, 2018). To resolve this problem, we will project a
multiply-imputed pseudo-outcome derived from the BART output onto the support
of $\mathcal{A}$ with local regression techniques. We only need a small number
of imputations, i.e., a small, thinned subset of the BART posterior samples
(indexed by $l=1,2,\ldots,L$, $L\ll K$), to get smooth, reliable estimates of
the ERF. First, we create multiple imputations of the pseudo-outcome as

$\xi^{(l)}\left(A^{(l)}_{i},\mathbf{X}_{i},\bar{Y}_{i}\right)=\left[\bar{Y}_{i}-\mu^{(l)}_{\bar{Y}}\left(A^{(l)}_{i},\mathbf{X}_{i}\right)\right]+\int_{\mathcal{X}}\mu^{(l)}_{\bar{Y}}\left(A^{(l)}_{i},\mathbf{x}\right)d\mathbb{P}(\mathbf{x}).$
(4)

There are two components in this pseudo-outcome. The integral on the right
hand side of the addition symbol is the marginal estimate of the ERF at
$A^{(l)}_{i}$. Under the more typical Bayesian framework, we could find
$\int_{\mathcal{X}}\mu^{(k)}_{\bar{Y}}\left(a,\mathbf{x}\right)d\mathbb{P}(\mathbf{x})$
for each $k=1,2,\ldots,K$. This would be equivalent to a Bayesian version of a
G-computation estimator for the ERF (Keil et al.,, 2018). However, this
process can be time consuming when iterating for all $K$ posterior samples
instead of the $L$ imputations, in addition to yielding a non-smooth ERF. The
term on the left-hand side of the addition symbol in (4) is the residual error
conditioned on $A^{(l)}_{i}$ and $\mathbf{X}_{i}$. Since multiple imputation
couples Bayesian methodology with regression techniques, it is necessary for
each imputation that we approximate the variance for estimates of $\theta(a)$
at each point $a\in\mathcal{A}$ which is facilitated by the addition of this
residual error term. Without the residual error component, the variance
estimates within each imputation would be incorrect (Antonelli et al.,, 2020).

Local regression methods offer perhaps the most flexible means of projecting
the pseudo-outcomes in (4) onto the support $a\in\mathcal{A}$. Given a
bandwidth $h>0$, the ERF estimates we obtain for each imputation
$l=1,2,\ldots,L$ are denoted by $\hat{\theta}^{(l)}_{h}(a)$. A pointwise
estimate of the ERF is summarized by
$\bar{\theta}_{h}(a)=L^{-1}\sum_{l=1}^{L}\hat{\theta}^{(l)}_{h}(a)$. Because
we regress each of the imputed pseudo-outcomes onto the support of the
exposure, $a\in\mathcal{A}$, the values $\hat{\theta}^{(l)}_{h}(a)$ do not
form a proper posterior of $\theta(a)$. This result was noted also by
Antonelli et al., (2020) who found that the posterior distribution of a
regression based estimator like $\hat{\theta}^{(l)}_{h}(a)$ alone did not
adequately characterize the uncertainty. To correct for this, our approach to
estimation combines a kernel-weighted least-squares regression estimator
similar to Kennedy et al., (2017) with the Bayesian approach of Antonelli et
al., (2020). Instead of using a bootstrap estimator to find the empirical MCMC
standard error used by Antonelli et al., (2020), we use an asymptotic standard
error estimator derived by Kennedy et al., (2017). Details for estimating
$\hat{\theta}^{(l)}_{h}(a)$ and the accompanying standard errors for
$\bar{\theta}_{h}(a)$ using multiple imputation combining rules (Rubin,, 2004)
are provided in Supplemental Section S2 along with the details on how to
smooth the BART ERF estimates using kernel-weighted least-squares regression.

A simple regression calibration variation to the above approach would be to
find a single imputation (i.e. $L=1$) of $A_{i}$, let’s say $\hat{A}_{i}$,
then construct an estimator of the outcome mean
$\hat{\mu}_{\bar{Y}}(\hat{A}_{i},\mathbf{X}_{i})$ and substitute that
estimator into (4) replacing
$\mu_{\bar{Y}}^{(1)}\left(A_{i}^{(1)},\mathbf{X}_{i}\right)$. In this single
imputation case, the ERF estimator would equal $\hat{\theta}^{(1)}_{h}(a)$. We
will see in the simulation that for a single imputation, the choice of
$\hat{A}_{i}$ is not so straightforward in the presence of measurement error,
nor does such an approach adequately propagate the uncertainty created by
measurement error without further correction. For more details on the local
regression methods we apply, see Supplemental Section S2.

###  3.3 Issues Stemming from Congeniality

Our model contains two components, an imputation stage and an analysis stage.
The imputation stage generates imputations of the latent exposures, $A_{i}$
and $S_{ij}$, while the analysis stage draws samples from the posterior
distributions of the model parameters. Following the multiple imputation
literature, the imputation component of our model must condition on the
outcome to satisfy assumptions associated with congeniality, explained below.
Congeniality states that the imputation and analysis stage models need to
utilize the same data (Meng,, 1994), and violations of congeniality can bias
parameter estimates. From a Bayesian perspective, non-congeniality can arise
due to cutting feedback or modularizing (Zigler et al.,, 2013) the components
of a joint model. The limiting distribution from an MCMC of $A^{(l)}_{i}$, and
by extension the limiting distribution of the exposure effect, is ill-defined
when cutting feedback between the outcome and exposure model in a Bayesian
setting (Plummer,, 2015). Since we are necessarily using the imputations of
$A_{i}$ to fit the analysis model $\mu_{\bar{Y}}(A_{i},\mathbf{X}_{i})$ and
construct an estimate of the ERF, then to satisfy congeniality and avoid
biasing the ERF estimate, the imputation models must be conditional on the
outcome. This is accomplished in our method by using a fully Bayesian joint
model-fitting scheme for the EPE, GPS, and outcome models.

However, sampling the latent exposures $A_{i}$ conditioned on the outcome
creates bi-directional feedback between the traditional “design stage” (GPS
modeling) and analysis stage that are kept separate in most causal analyses.
This seemingly defies research that emphasizes cutting feedback between the
GPS and outcome model in a Bayesian causal analysis (Zigler et al.,, 2013).
Notice that we include a model for the GPS in (2) that we fit in Section 3.1,
yet we do not utilize the GPS in our outcome model nor our estimator of the
ERF in Section S2. While the outcome data does not appear in the full
conditionals for the GPS model parameters, it is used to generate new
predictions of $A_{i}$ which may indirectly create feedback affecting the GPS
model estimates. To counteract the feedback problems created by congeniality,
we removed the GPS adjustments from the doubly-robust pseudo-outcome supposed
by Kennedy et al., (2017) (which appears in the Supplemental Section S4).
These GPS adjustments appear to provide no utility in our context as
demonstrated by a small simulation study contained within Supplemental Section
S4. In this simulation, we draw particular attention to cases where the
outcome model is misspecified but the GPS is correctly specified. In this
scenario, there is almost no difference between using the GPS adjusted pseudo-
outcome supported by Kennedy et al., (2017), and the pseudo-outcome we suppose
in (4). We contend that this null result is attributable to the feedback
created by the need to use congenial models for the imputation and analysis
stages of the MCMC sampler.

##  4 Numerical Example

###  4.1 Simulation Design

In this simulation study, we examine the performance of the method described
in Section 3 in the presence and absence of exposure measurement error. In
addition, we will evaluate the effects of model misspecification across the
three components of the likelihood in (2), all while examining whether
uncertainty is properly propagated into the final ERF estimate.

We test four different methods for measurement error correction. The first
method naïvely ignores any measurement error and assumes $\tilde{Z}_{i}$ is
the true cluster level exposure. The second method examines an extension to
the regression calibration approach proposed by Wu et al., (2019), which we
adapted to consider a continuous exposure, that corrects for prediction error
but disregards the remaining classical error in (1) created by the remaining
aggregation error. After finding predictions $\hat{S}_{ij}$ for $S_{ij}$ using
least-squares regression conditioned on $\tilde{S}_{ij}$ and
$\mathbf{W}_{ij}$, a single imputation of $A_{i}$ is produced by computing
$\hat{Z}_{i}=M_{i}^{-1}\sum_{j=1}^{M_{i}}\hat{S}_{ij}$. We will refer to these
two implementations that use either $\tilde{Z}_{i}$ or $\hat{Z}_{i}$ as the
“single imputation” approaches. The third and fourth methods that we examine
are two “multiple imputation” approaches that follow the proposal in Section
3, with different outcome model specifications. In one of the multiple
imputation variants, we let
$\mu^{(l)}_{\bar{Y}}\left(A_{i}^{(l)},\mathbf{X}_{i}\right)$ be a BART model.
In another, we let
$\mu^{(l)}_{\bar{Y}}\left(A_{i}^{(l)},\mathbf{X}_{i}\right)$ be a correctly-
specified log-linear Poisson model, assuming unknown coefficient values (see
Supplemental Sections S3 and S4 for details). To highlight the impacts of
correcting measurement error and propagating uncertainty, we utilize
estimation approaches that are as similar as possible across the four
competing methods aside from how they correct for measurement error. For the
single imputation approaches, we construct the pseudo-outcome in (4) in the
exact same manner as in the multiple imputation cases, but for a single
iteration, i.e. setting $L=1$. For the single imputation methods, this means
using the $\tilde{Z}_{i}$ or $\hat{Z}_{i}$ in place of $A^{(l)}_{i}$ and
either $\tilde{\mu}_{\bar{Y}}(\cdot)$ or $\hat{\mu}_{\bar{Y}}(\cdot)$ in place
of $\mu^{(l)}_{\bar{Y}}(\cdot)$, respectively. We also assume that
$\phi_{i}=0$ for all $i=1,2,\ldots,n$ (i.e. $A_{i}$ is independent), both in
the data generation procedure and in the models we fit.

In Supplemental Section S3 we describe the various data generating mechanisms
used to construct the simulation scenarios. For each scenario, we ran $500$
iterations, in each one applying each of the four methods described above. In
each iteration, we obtain pointwise estimates of $\theta(a)$ at $201$ equally
spaced exposure values across the range $a\in[6,14]$. We report the relative
bias, the residual mean square error, and coverage probabilities averaged
across all the $500$ iterations and all the $201$ pointwise ERF estimates. We
also report these same measurements for predictions at a single exposure
level, $a=11$, averaged across the $500$ iterations. Briefly, we vary
$n\in\\{400,800\\}$, $m\in\\{5n,10n\\}$, $\tau^{2}\in\\{0,1,2\\}$, and
$\omega^{2}\in\\{0,1,2\\}$. Note that when $\tau^{2}=0$ the prediction error
is absent, whereas when $\omega^{2}=0$ the remaining classical error is absent
and $U_{ij}$ is limited to the prediction error. We also vary the degree of
misspecification in the component models contained in (2) (Kang and Schafer,,
2007).

###  4.2 Simulation Results

The main results of this simulation illustrate how adding measurement error
can influence estimates of the ERF. Figure 3 shows the average ERF estimated
with each of the four methods in Section 4.1 (relative to the true ERF) under
simulation scenarios with no measurement error, prediction error only,
classical error only, and both sources of error present. These results
demonstrate that generating unbiased imputations of the exposure of interest
is critical, whether using a single or multiple imputations. Note that from
the description of the scenarios, we have
$\mathbb{E}[\tilde{S}_{ij}|A_{i}]=A_{i}$ when the EPE model is correctly
specified. Therefore the naïve approach and the regression calibration
approach, which corrects for prediction error, produce similar results as
shown by the root mean squared errors displayed in Figure 3. This is because
the prediction error is Berkson, so using regression calibration has a
negligible effect. Also, note that the regression calibration and multiple
imputation approaches perform similarly in scenarios without classical error
(Figure 3, top row), but the performance of regression calibration
deteriorates when we simulate measurement error with aggregation error
present. Meanwhile, the multiple imputation approaches retain accuracy (Figure
3, bottom row). Some bias still lingers in each method, particularly at the
extremes of the exposure support (i.e. $a=6$ and $a=14$). This is typical of
local regression methods where there is an implicit bias-variance tradeoff
through the choice of $h$ (which is set to $h=0.2$ when $n=800$ and $h=0.4$
when $n=400$) (Hastie et al.,, 2001). Unsurprisingly, the accuracy of the
different methods we test improves as both $n$ and $m$ increase (Table 1).

We can also see in Table 1 that the coverage probabilities from the 95%
confidence interval estimates among the multiple imputation approaches offer a
marked improvement over the alternatives. When using a correctly specified
log-linear outcome model within a multiple imputation setting, we achieve
coverage probabilities that match the nominal 95% confidence level. However,
the coverage probabilities are imperfect when using a BART outcome model.
BARTs’ tendency to under-report the dispersion of point estimates is well-
documented (Wendling et al.,, 2018; Hahn et al.,, 2020; Nethery et al.,,
2021). The BART output also seems to produce biased estimates at several
points along the curve, contributing to the low coverage probabilities.
Despite this limitation, BART remains the preferred outcome model due to its
strong predictive abilities even when the correct model specification is
unknown (see the following paragraph). Since we average pointwise evaluations
over the support $\mathcal{A}$, we can also see how the coverage probabilities
are affected by using local regression techniques. This can be seen by
examining the single pointwise evaluations at $a=11$, which are closer to the
nominal 95% confidence level than the average coverage probability over the
$201$ pointwise evaluations.

In the second part to our simulation study, where we induce model
misspecification, we observe in Table 1 that the ERF estimates are largely
unaffected by outcome model misspecification alone with a BART model –
remaining relatively unbiased as long as the GPS is still correctly specified.
The same result is obviously not true for a parametric log-linear outcome
model. Even though BART is supplied the untransformed covariates when model
misspecification is present, it is quite adept at identifying complex
nonlinear and interacting effects, making outcome model misspecification a
moot consideration, so long as Assumptions 3 and 4 hold. That said, the
coverage probability decreases considerably in scenarios with a misspecified
outcome model relative to the results of the scenario where all models are
correctly specified.

When the GPS model is misspecified, the level of bias increases substantially
in each of the ERF estimates, both from the single imputation and multiple
imputation models. A perplexing exception to this result is when both the GPS
and the outcome model are misspecified, in which case the bias returns to the
nominal levels observed in cases when the GPS is correctly specified and a
BART outcome model is used. Given these results, and the results of the
simulation in Supplemental Section S4, it is evident that finding an
approximately correct GPS model may be even more important than finding a
correct outcome model as we originally suggested. Doing so might decrease the
precision of the imputations for $A_{i}$. Finally, when the EPE model is
“misspecified”, meaning $\mathbb{E}(\tilde{S}_{ij}|A_{i})\neq A_{i}$ implying
$\mathbb{E}(\tilde{Z}_{i}|A_{i})\neq A_{i}$, then we see the most severe
levels of bias using the naïve (no correction) approach. Since the multiple
imputation and regression calibration approaches correctly model
$S^{(k)}_{ij}$ and $\hat{S}_{ij}$, respectively, such that
$\mathbb{E}(S^{(k)}_{ij}|A_{i})=A_{i}$ and
$\mathbb{E}(\hat{S}_{ij}|A_{i})=A_{i}$, the prediction bias does not affect
the estimates of the exposure response curves using these two approaches.

##  5 Applied Example

Recall the example described in Section 1.1. Wu et al., (2019) used regression
calibration methods on the same data to correct for measurement error for the
grid-year PM2.5 predictions. Subsequently, they aggregated the corrected grid-
year predictions to the ZIP-code-year level, categorized the imputed exposures
based on standard EPA cutoffs, and examined a variety of GPS-based
implementations including matching, sub-classification, and weighting to
estimate the effect of PM2.5 on mortality. Their results suggest that
increasing PM2.5 exposure led to excess mortality events.

In this analysis we re-examine the same data as Wu et al., (2019), applying
the proposed multiple imputation approach to estimate the ERF relating PM2.5
with all-cause mortality. For ZIP-code-year $i$, the exposure $A_{i}$ is the
annual average PM2.5 concentration (in $\mu$g/m3), the outcome $Y_{i}$ is the
number of deaths observed, and $N_{i}$ is the person-years at risk. We have
exposure data on $n=28{,}626$ ZIP-code-years in New England covered by
$m=2{,}395{,}588$ 1km$\times$1km grid-years, each with error-prone PM2.5
exposure predictions $\tilde{S}_{ij}$. A subset of $83$ grids within $75$ ZIP-
codes have error-free PM2.5 measurements from monitoring stations in years
2000-2012 (Figure 1). While it is plausible that the $\tilde{S}_{ij}$ and
$S_{ij}$ within ZIP-code-years are correlated, the majority of this
correlation can be accounted for simply by conditioning on $A_{i}$, thus
enabling Assumptions 1 and 2.

For the multiple imputation approach, we assume the outcome model in (3). As
in the simulation study, we implement two different single imputation
approaches. In the first approach, we ignore prediction error altogether and
suppose $\tilde{Z}_{i}$ is the true exposure (we call this the “no correction”
approach). In the second approach, we use least squares regression on the
validation data to get the regression calibrated $\hat{Z}_{i}$. We then use
these single imputations of the exposure in a BART model to construct an
estimate of the mean function; either
$\tilde{\mu}_{\bar{Y}}(\tilde{Z}_{i},\mathbf{X}_{i})$ or
$\hat{\mu}_{\bar{Y}}(\hat{Z}_{i},\mathbf{X}_{i})$. For both the single ($L=1$)
and multiple imputation approaches ($L=100$), local regression was applied
according to Section S2. We collected $5{,}000$ MCMC samples, thinned by every
$50$ iterations, after a burn-in of another $5{,}000$ iterations. Diffuse
priors specified in Section 3.1 were used for the GPS model parameters, the
EPE model parameters, as well as for the BART model. As opposed to the
simulation study, we account for spatial correlation between the ZIP-code
exposures by sampling $\phi_{i}$ in (2). In this example, the matrix
$\mathbf{V}$ is a binary adjacency matrix with each row and column
representing a different ZIP-code-year.

Figure 1 provides a heat-map displaying the difference between the naïve
exposure imputations of $\tilde{Z}_{i}$ and the posterior means
$L^{-1}\sum_{l=1}^{L}A^{(l)}_{i}$ for each ZIP-code in 2010. Here we can see
substantial differences in the imputed exposure values in some areas of New
England. Figure 4 shows the ERF estimates from each of the three methods we
examined. In the curve estimated with multiple imputation, we can see a modest
increase in the rate of all-cause mortality from $4.7\%$ when annual average
PM2.5 is at $5$ $\mu$g/m3 to $5.0\%$ when annual average PM2.5 is at $15$
$\mu$g/m3. Relative to the uncorrected and regression calibration curves, the
curve estimated with multiple imputation has a steeper gradient at lower
levels of PM2.5. The steepest increase in mortality occurs when PM2.5 changes
from about $5$ $\mu$g/m3to $10$ $\mu$g/m3, which is policy-relevant given the
World Health Organization’s recent decision to lower the recommended limit of
annual average PM2.5 to $5$ $\mu$g/m3 (World Health Organization,, 2021).

We suspect that the minor differences between the three curves are primarily
attributable to the small variation of the classical measurement error,
$\omega^{2}$, which we found to be equal to about 0.0187 ($95\%$ CI:
0.0183-0.0189). The posterior mean for the conditional prediction error
variance, $\tau^{2}$, on the other hand was 1.340 ($95\%$ CI: 1.212, 1.446).
However, as we demonstrated in the simulation study, correcting for the
prediction error is secondary to correcting for the classical measurement
error including the aggregation error, except when
$\mathbb{E}(\tilde{S}_{ij}|A_{i})\neq A_{i}$. This reasoning is further
cemented in this illustrative example – we can see that
$\mathbb{E}(\tilde{S}_{ij}|A_{i})\approx A_{i}$, and so the difference between
the curves with no measurement error correction and regression calibration is
relatively small.

##  6 Discussion

We developed a multiple imputation framework that addresses exposure
measurement error when estimating a causal ERF. We adapted a Bayesian approach
for correcting measurement error to generate imputations of the true exposures
assuming non-differential classical measurement error. We then constructed a
flexible representation of the ERF using kernel smoothed output from a
Bayesian additive regression tree. These two steps interface with one another
over several iterations. Since the final estimator isn’t formally Bayesian due
to post-processing the BART adjusted pseudo-outcome with kernel weighted least
squares regression, we regard it as a multiple imputation-based estimator. We
then described some issues that we encountered with this approach regarding
the seeming conflict that exists between cutting feedback and appeasing rules
about congeniality. We provided several simulated numerical examples that
showcase why correcting for exposure measurement error within a causal-
analysis is important to obtain unbiased estimates and valid inferences.

In the real data example, we estimated the exposure response function
associating all-cause mortality with annual average PM2.5 exposures. At first
glance, the three curves we estimate with varying degrees of measurement error
correction all seem fairly similar for this particular analysis. These results
seem to indicate that there is little benefit to correcting for measurement
error while estimating ERFs, though the attenuation at lower PM2.5 levels is
meaningful. However, the reason for the similarities are mainly due to the low
levels of measurement error. As measurement error increases, so does the
attenuation bias. Indeed, our attempts to correct for measurement error in our
applied example does not seem to have a large impact on the results
associating PM2.5 with mortality. Therefore, this analysis helps validate
other analyses that used similar data but did not correct for measurement
error.

While Medicare data limitations necessitate use of ZIP-code aggregate
exposures in our analysis, there are various simplifications to Sections 2 and
3 that can be made to accommodate scenarios where certain measurement error
sources are absent. If we were provided individual level outcomes and error-
prone exposure data, then the measurement error would be closer to a Berkson-
type rather than classical (Bateson and Wright,, 2010). Accounting for
measurement error in this problem would require a different EPE model than the
one proposed, however, adapting our framework should be straightforward given
the adaptability of the Bayesian framework that we outlined (Carroll et al.,,
2006). Lastly, our approach assumes non-differential measurement error, which
in our context means that $Y_{i}$ is unaffected by $\tilde{S}_{ij}$ given
$A_{i}$. More work is required to generalize this methodology to be useful in
cases where the measurement error is differential.

Another line of future work involves relaxing the independence and
homoscedastic assumptions surrounding the measurement error (Assumptions 1 and
2). Since these models are deployed within an MCMC framework, more
sophisticated spatiotemporal modelling techniques might be employed to better
understand the measurement error residing in the data. Second, alternatives to
the parametric GPS and EPE models were not implemented in this work. However,
more flexible nonparametric methods may help reduce model-dependency of the
ERF. We relaxed some of this model-dependency by using a BART outcome model.
However, the BART implementation we employed assumes the outcome is
approximately Gaussian. There is an option for fitting log-linear BART models
(Murray,, 2021) which might ameliorate some of the problems we encountered
with BART, like the poor coverage probabilities observed in the simulation
study (Hahn et al.,, 2020). However, at the time of writing, code to fit this
alternative BART model could not be obtained. This problem may also be
addressed by choosing a different outcome model as demonstrated in the
simulations.

## Acknowledgements

Funding was provided by the National Institute of Health (NIH) grants (T32
ES007142, R01 ES030616, R01 ES028033, R01 ES026217, K01 ES032458-01) and HEI
grant 4953-RFA14-3/16-4. The contents are solely the responsibility of the
grantee and do not necessarily represent the official views of the funding
agencies. Further, the funding agencies do not endorse the purchase of any
commercial products or services related to this publication.

## Reproducible R Code

The R functions used to fit the model in Section 3 along with the code used to
run the simulation studies in Sections 4 and S4 are available at
https://github.com/kevjosey/causal-me/.

## References

  * Anderson, (2009) Anderson, H. R. (2009).  Air pollution and mortality: A history.  Atmospheric Environment, 43(1):142–152.
  * Antonelli et al., (2020) Antonelli, J., Papadogeorgou, G., and Dominici, F. (In Press 2020).  Causal inference in high dimensions: A marriage between bayesian modeling and good frequentist properties.  Biometrics.
  * Bateson and Wright, (2010) Bateson, T. F. and Wright, J. M. (2010).  Regression calibration for classical exposure measurement error in environmental epidemiology studies using multiple local surrogate exposures.  American journal of epidemiology, 172(3):344–352.
  * Braun et al., (2017) Braun, D., Gorfine, M., Parmigiani, G., Arvold, N. D., Dominici, F., and Zigler, C. (2017).  Propensity scores with misclassified treatment assignment: A likelihood-based adjustment.  Biostatistics, 18(4):695–710.
  * Brook et al., (2010) Brook, R. D., Rajagopalan, S., Pope III, C. A., Brook, J. R., Bhatnagar, A., Diez-Roux, A. V., Holguin, F., Hong, Y., Luepker, R. V., Mittleman, M. A., et al. (2010).  Particulate matter air pollution and cardiovascular disease: an update to the scientific statement from the american heart association.  Circulation, 121(21):2331–2378.
  * Carroll et al., (2006) Carroll, R. J., Ruppert, D., Stefanski, L. A., and Crainiceanu, C. M. (2006).  Measurement error in nonlinear models: a modern perspective.  Chapman and Hall/CRC, Boca Raton, FL, USA.
  * Chipman et al., (1998) Chipman, H. A., George, E. I., and McCulloch, R. E. (1998).  Bayesian cart model search.  Journal of the American Statistical Association, 93(443):935–948.
  * 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.
  * Cole et al., (2006) Cole, S. R., Chu, H., and Greenland, S. (2006).  Multiple-imputation for measurement-error correction.  International journal of epidemiology, 35(4):1074–1081.
  * Danesh Yazdi et al., (2021) Danesh Yazdi, M., Wang, Y., Di, Q., Wei, Y., Requia, W. J., Shi, L., Sabath, M. B., Dominici, F., Coull, B. A., Evans, J. S., et al. (2021).  Long-term association of air pollution and hospital admissions among medicare participants using a doubly robust additive model.  Circulation, 143(16):1584–1596.
  * Di et al., (2016) Di, Q., Koutrakis, P., and Schwartz, J. (2016).  A hybrid prediction model for pm2. 5 mass and components using a chemical transport model and land use regression.  Atmospheric environment, 131:390–399.
  * Di et al., (2017) Di, Q., Wang, Y., Zanobetti, A., Wang, Y., Koutrakis, P., Choirat, C., Dominici, F., and Schwartz, J. D. (2017).  Air pollution and mortality in the medicare population.  New England Journal of Medicine, 376(26):2513–2522.
  * Dominici et al., (2006) Dominici, F., Peng, R. D., Bell, M. L., Pham, L., McDermott, A., Zeger, S. L., and Samet, J. M. (2006).  Fine particulate air pollution and hospital admission for cardiovascular and respiratory diseases.  Jama, 295(10):1127–1134.
  * Gill and Robins, (2001) Gill, R. D. and Robins, J. M. (2001).  Causal inference for complex longitudinal data: the continuous case.  Annals of Statistics, pages 1785–1811.
  * Haber et al., (2020) Haber, G., Sampson, J., and Graubard, B. (2020).  Bias due to berkson error: issues when using predicted values in place of observed covariates.  Biostatistics, 22(4):858–872.
  * Hahn et al., (2020) Hahn, P. R., Murray, J. S., and Carvalho, C. M. (2020).  Bayesian regression tree models for causal inference: Regularization, confounding, and heterogeneous effects (with discussion).  Bayesian Analysis, 15(3):965–1056.
  * Hajat et al., (2002) Hajat, S., Anderson, H., Atkinson, R., and Haines, A. (2002).  Effects of air pollution on general practitioner consultations for upper respiratory diseases in london.  Occupational and environmental medicine, 59(5):294–299.
  * Hastie et al., (2001) Hastie, T., Tibshirani, R., and Friedman, J. (2001).  The Elements of Statistical Learning.  Springer Series in Statistics. Springer New York Inc., New York, NY, USA.
  * Kang and Schafer, (2007) Kang, J. D. Y. and Schafer, J. L. (2007).  Demystifying Double Robustness: A Comparison of Alternative Strategies for Estimating a Population Mean from Incomplete Data.  Statistical Science, 22(4):523 – 539.
  * Keil et al., (2018) Keil, A. P., Daza, E. J., Engel, S. M., Buckley, J. P., and Edwards, J. K. (2018).  A bayesian approach to the g-formula.  Statistical methods in medical research, 27(10):3183–3204.
  * Kennedy et al., (2017) Kennedy, E. H., Ma, Z., McHugh, M. D., and Small, D. S. (2017).  Nonparametric methods for doubly robust estimation of continuous treatment effects.  Journal of the Royal Statistical Society. Series B, Statistical Methodology, 79(4):1229.
  * Kim et al., (2018) Kim, W., Kwon, K., Kwon, S., and Lee, S. (2018).  The identification power of smoothness assumptions in models with counterfactual outcomes.  Quantitative Economics, 9(2):617–642.
  * Kioumourtzoglou et al., (2014) Kioumourtzoglou, M.-A., Spiegelman, D., Szpiro, A. A., Sheppard, L., Kaufman, J. D., Yanosky, J. D., Williams, R., Laden, F., Hong, B., and Suh, H. (2014).  Exposure measurement error in pm 2.5 health effects studies: a pooled analysis of eight personal exposure validation studies.  Environmental Health, 13(1):1–11.
  * Lee, (2013) Lee, D. (2013).  Carbayes: An r package for bayesian spatial modeling with conditional autoregressive priors.  Journal of Statistical Software, 55(13):1–24.
  * Lenis et al., (2017) Lenis, D., Ebnesajjad, C. F., and Stuart, E. A. (2017).  A doubly robust estimator for the average treatment effect in the context of a mean-reverting measurement error.  Biostatistics, 18(2):325–337.
  * Leroux et al., (2000) Leroux, B. G., Lei, X., and Breslow, N. (2000).  Estimation of disease rates in small areas: a new mixed model for spatial dependence.  In Statistical models in epidemiology, the environment, and clinical trials, pages 179–191. Springer.
  * Lewbel, (2007) Lewbel, A. (2007).  Estimation of average treatment effects with misclassification.  Econometrica, 75(2):537–551.
  * Meng, (1994) Meng, X.-L. (1994).  Multiple-Imputation Inferences with Uncongenial Sources of Input.  Statistical Science, 9(4):538 – 558.
  * Murray, (2021) Murray, J. S. (2021).  Log-linear bayesian additive regression trees for multinomial logistic and count regression models.  Journal of the American Statistical Association, 116(534):756–769.
  * Nethery et al., (2021) Nethery, R. C., Mealli, F., Sacks, J. D., and Dominici, F. (2021).  Evaluation of the health impacts of the 1990 clean air act amendments using causal inference and machine learning.  Journal of the American Statistical Association, 116(535):1128–1139.
  * Papadogeorgou et al., (2019) Papadogeorgou, G., Mealli, F., and Zigler, C. M. (2019).  Causal inference with interfering units for cluster and population level treatment allocation programs.  Biometrics, 75(3):778–787.
  * Plummer, (2015) Plummer, M. (2015).  Cuts in bayesian graphical models.  Statistics and Computing, 25(1):37–43.
  * Pope III et al., (2015) Pope III, C. A., Turner, M. C., Burnett, R. T., Jerrett, M., Gapstur, S. M., Diver, W. R., Krewski, D., and Brook, R. D. (2015).  Relationships between fine particulate air pollution, cardiometabolic disorders, and cardiovascular mortality.  Circulation research, 116(1):108–115.
  * Ren et al., (2021) Ren, B., Wu, X., Braun, D., Pillai, N., and Dominici, F. (2021).  Bayesian modeling for exposure response curve via gaussian processes: Causal effects of exposure to air pollution on health outcomes.  arXiv preprint arXiv:2105.03454.
  * Rhee et al., (2019) Rhee, J., Dominici, F., Zanobetti, A., Schwartz, J., Wang, Y., Di, Q., Balmes, J., and Christiani, D. C. (2019).  Impact of long-term exposures to ambient pm2. 5 and ozone on ards risk for older adults in the united states.  Chest, 156(1):71–79.
  * Rubin, (1974) Rubin, D. B. (1974).  Estimating causal effects of treatments in randomized and nonrandomized studies.  Journal of educational Psychology, 66(5):688.
  * Rubin, (2004) Rubin, D. B. (2004).  Multiple imputation for nonresponse in surveys.  John Wiley & Sons, New York, NY, USA.
  * Webb-Vargas et al., (2017) Webb-Vargas, Y., Rudolph, K. E., Lenis, D., Murakami, P., and Stuart, E. A. (2017).  An imputation-based solution to using mismeasured covariates in propensity score analysis.  Statistical methods in medical research, 26(4):1824–1837.
  * Wendling et al., (2018) Wendling, T., Jung, K., Callahan, A., Schuler, A., Shah, N., and Gallego, B. (2018).  Comparing methods for estimation of heterogeneous treatment effects using observational data from health care databases.  Statistics in medicine, 37(23):3309–3324.
  * World Health Organization, (2021) World Health Organization (2021).  New who global air quality guidelines aim to save millions of lives from air pollution.  https://www.who.int/.
  * Wu et al., (2019) Wu, X., Braun, D., Kioumourtzoglou, M.-A., Choirat, C., Di, Q., and Dominici, F. (2019).  Causal inference in the context of an error prone exposure: air pollution and mortality.  The Annals of Applied Statistics, 13(1):520.
  * Wu et al., (2020) Wu, X., Braun, D., Schwartz, J., Kioumourtzoglou, M., and Dominici, F. (2020).  Evaluating the impact of long-term exposure to fine particulate matter on mortality among the elderly.  Science advances, 6(29):eaba5692.
  * Wu et al., (2018) Wu, X., Mealli, F., Kioumourtzoglou, M.-A., Dominici, F., and Braun, D. (2018).  Matching on generalized propensity scores with continuous exposures.  arXiv:1812.06575.
  * Yazdi et al., (2019) Yazdi, M. D., Wang, Y., Di, Q., Zanobetti, A., and Schwartz, J. (2019).  Long-term exposure to pm2. 5 and ozone and hospital admissions of medicare participants in the southeast usa.  Environment international, 130:104879.
  * Zanobetti et al., (2009) Zanobetti, A., Franklin, M., Koutrakis, P., and Schwartz, J. (2009).  Fine particulate air pollution and its components in association with cause-specific emergency admissions.  Environmental Health, 8(1):1–12.
  * Zhu et al., (2017) Zhu, L., Ge, X., Chen, Y., Zeng, X., Pan, W., Zhang, X., Ben, S., Yuan, Q., Xin, J., Shao, W., et al. (2017).  Short-term effects of ambient air pollution and childhood lower respiratory diseases.  Scientific reports, 7(1):1–7.
  * Zigler et al., (2013) Zigler, C. M., Watts, K., Yeh, R. W., Wang, Y., Coull, B. A., and Dominici, F. (2013).  Model feedback in bayesian propensity score estimation.  Biometrics, 69(1):263–273.

## Tables and Figures

Figure 1: ZIP-code and monitor locations (green dots) appearing across New-
England in 2010. The heat map describes the difference between the posterior
mean $L^{-1}\sum_{l=1}^{L}A^{(l)}_{i}$ and $\tilde{Z}_{i}$ in 2010. ZIP-codes
with a grey fill had either missing exposure data or missing outcome data in
2010.

Figure 2: An illustration about the nested structure of the data accompanying
the relationship between the true ZIP-code, the true grid level, and the EPE
grid measurements.

Figure 3: ERF estimates averaged over $500$ simulated iterations when $n=800$
and $M=4000$ and correct model specifications. When classical error is present
then $\omega^{2}=2$, when prediction error is present then $\tau^{2}=1$. The
BART approach to multiple imputation uses a BART outcome model while the GLM
implementation uses a log-linear outcome model. The bar plot overlaying the
ERF estimates represents the root mean squared errors for the corresponding
methods.

$n$ | $m$ | $\omega^{2}$ | $\tau^{2}$ | GPS | Outcome | EPE |  | Relative Bias |  | 95% Coverage Probability
---|---|---|---|---|---|---|---|---|---|---
| No Correction |  | Regression
---
Calibration
| BART Multiple
---
Imputation
| GLM Multiple
---
Imputation
| No Correction |  | Regression
---
Calibration
| BART Multiple
---
Imputation
| GLM Multiple
---
Imputation
400 | 2,000 | 1 | 1 |  |  |  |  | 0.35 (-0.09) | 0.28 (-0.09) | 0.21 (-0.06) | -0.03 (-0.02) |  | 0.44 (0.52) | 0.53 (0.54) | 0.75 (0.78) | 0.83 (0.90)
400 | 2,000 | 1 | 2 |  |  |  |  | 0.42 (-0.11) | 0.30 (-0.09) | 0.28 (-0.06) | -0.05 (-0.02) |  | 0.39 (0.40) | 0.52 (0.46) | 0.74 (0.73) | 0.70 (0.91)
400 | 2,000 | 2 | 1 |  |  |  |  | 0.38 (-0.12) | 0.32 (-0.11) | 0.22 (-0.08) | -0.03 (-0.02) |  | 0.39 (0.45) | 0.45 (0.48) | 0.71 (0.70) | 0.85 (0.90)
400 | 2,000 | 2 | 2 |  |  |  |  | 0.42 (-0.14) | 0.29 (-0.13) | 0.21 (-0.08) | -0.05 (-0.02) |  | 0.36 (0.35) | 0.45 (0.44) | 0.71 (0.62) | 0.75 (0.92)
400 | 4,000 | 1 | 1 |  |  |  |  | 0.30 (-0.06) | 0.21 (-0.06) | 0.20 (-0.04) | -0.03 (-0.02) |  | 0.53 (0.66) | 0.63 (0.69) | 0.78 (0.78) | 0.80 (0.91)
400 | 4,000 | 1 | 2 |  |  |  |  | 0.31 (-0.08) | 0.18 (-0.07) | 0.21 (-0.05) | -0.06 (-0.02) |  | 0.47 (0.58) | 0.63 (0.65) | 0.78 (0.76) | 0.64 (0.92)
400 | 4,000 | 2 | 1 |  |  |  |  | 0.29 (-0.09) | 0.23 (-0.08) | 0.19 (-0.06) | -0.02 (-0.02) |  | 0.47 (0.58) | 0.55 (0.64) | 0.74 (0.72) | 0.86 (0.92)
400 | 4,000 | 2 | 2 |  |  |  |  | 0.30 (-0.11) | 0.19 (-0.10) | 0.15 (-0.07) | -0.04 (-0.02) |  | 0.42 (0.48) | 0.58 (0.52) | 0.76 (0.71) | 0.69 (0.93)
800 | 4,000 | 1 | 1 |  |  |  |  | 0.29 (-0.08) | 0.21 (-0.07) | 0.12 (-0.03) | -0.04 (0.00) |  | 0.42 (0.45) | 0.53 (0.52) | 0.76 (0.75) | 0.77 (0.96)
800 | 4,000 | 1 | 2 |  |  |  |  | 0.35 (-0.11) | 0.22 (-0.09) | 0.12 (-0.04) | -0.07 (0.00) |  | 0.37 (0.38) | 0.50 (0.44) | 0.72 (0.74) | 0.62 (0.96)
800 | 4,000 | 2 | 1 |  |  |  |  | 0.36 (-0.10) | 0.29 (-0.09) | 0.15 (-0.04) | -0.04 (0.00) |  | 0.37 (0.38) | 0.41 (0.39) | 0.73 (0.71) | 0.82 (0.96)
800 | 4,000 | 2 | 2 |  |  |  |  | 0.43 (-0.12) | 0.31 (-0.10) | 0.13 (-0.04) | -0.06 (0.00) |  | 0.35 (0.30) | 0.43 (0.38) | 0.72 (0.72) | 0.65 (0.94)
800 | 8,000 | 1 | 1 |  |  |  |  | 0.20 (-0.06) | 0.15 (-0.05) | 0.10 (-0.03) | -0.04 (-0.01) |  | 0.50 (0.67) | 0.65 (0.69) | 0.78 (0.83) | 0.74 (0.94)
800 | 8,000 | 1 | 2 |  |  |  |  | 0.26 (-0.06) | 0.16 (-0.05) | 0.10 (-0.03) | -0.07 (0.00) |  | 0.44 (0.56) | 0.59 (0.60) | 0.72 (0.76) | 0.59 (0.94)
800 | 8,000 | 2 | 1 |  |  |  |  | 0.22 (-0.08) | 0.17 (-0.07) | 0.10 (-0.03) | -0.04 (-0.01) |  | 0.45 (0.52) | 0.55 (0.58) | 0.78 (0.78) | 0.78 (0.94)
800 | 8,000 | 2 | 2 |  |  |  |  | 0.28 (-0.08) | 0.19 (-0.07) | 0.11 (-0.03) | -0.07 (-0.01) |  | 0.42 (0.45) | 0.56 (0.53) | 0.73 (0.74) | 0.62 (0.94)
800 | 4,000 | 2 | 1 |  |  | ✓ |  | 0.91 (-0.09) | 0.30 (-0.09) | 0.15 (-0.04) | -0.03 (0.00) |  | 0.30 (0.47) | 0.42 (0.40) | 0.71 (0.73) | 0.76 (0.93)
800 | 4,000 | 2 | 1 |  | ✓ |  |  | 0.18 (-0.20) | 0.11 (-0.20) | -0.01 (-0.14) | 0.04 (0.00) |  | 0.25 (0.02) | 0.32 (0.03) | 0.55 (0.14) | 0.32 (0.96)
800 | 4,000 | 2 | 1 |  | ✓ | ✓ |  | 0.60 (-0.22) | 0.11 (-0.20) | -0.01 (-0.14) | 0.05 (0.00) |  | 0.19 (0.02) | 0.30 (0.02) | 0.55 (0.15) | 0.32 (0.97)
800 | 4,000 | 2 | 1 | ✓ |  |  |  | 0.52 (-0.06) | 0.46 (-0.05) | 0.29 (-0.01) | -0.03 (0.00) |  | 0.45 (0.54) | 0.49 (0.58) | 0.75 (0.82) | 0.75 (0.93)
800 | 4,000 | 2 | 1 | ✓ |  | ✓ |  | 1.28 (0.02) | 0.59 (-0.02) | 0.40 (0.02) | -0.03 (0.00) |  | 0.34 (0.66) | 0.49 (0.61) | 0.74 (0.82) | 0.74 (0.93)
800 | 4,000 | 2 | 1 | ✓ | ✓ |  |  | 0.15 (-0.21) | 0.08 (-0.20) | -0.05 (-0.15) | -0.09 (-0.05) |  | 0.18 (0.01) | 0.28 (0.01) | 0.44 (0.08) | 0.37 (0.77)
800 | 4,000 | 2 | 1 | ✓ | ✓ | ✓ |  | 0.52 (-0.25) | 0.02 (-0.23) | -0.09 (-0.16) | -0.09 (-0.05) |  | 0.15 (0.00) | 0.25 (0.00) | 0.41 (0.04) | 0.38 (0.72)
Table 1: Relative bias and coverage probabilities averaged over the simulated
ERF estimates when measurement error is present. The values in parentheses
represent the statistics evaluated at $a=11$. The check-marks indicates
whether the corresponding model labeled in the column header is misspecified.
The GLM approach refers to the multiple imputation implementation using a log-
linear outcome model. The BART approach to multiple imputation uses a BART
outcome model.

Figure 4: ERF estimate of PM2.5 on all-cause mortality in New England between
2000-2012 amongst Medicare recipients under three different approaches to
measurement error correction. The grey ribbon represents the 95% confidence
interval computed from our multiple imputation approach. The histogram
underlying the curves corresponds with the empirical distribution of the EPEs.

Supplement to “Estimating a Causal Exposure Response Function with a
Continuous Error-Prone Exposure: A Study of Fine Particulate Matter and All-
Cause Mortality”

##  S1 MCMC Sampling Algorithm

In this supplement, we describe the MCMC sampling steps and provide additional
details about the multiple imputation measurement error correction method that
we present in Section 3. We will assume a BART model for the outcome, and
Gaussian linear models for the GPS and EPE with
$\mu_{A}(\mathbf{X}_{i},\phi_{i})=\phi_{i}+\mathbf{X}^{T}_{i}\boldsymbol{\beta}$
and
$\mu_{S}(\tilde{S}_{ij},\mathbf{W}_{ij})=\left(\tilde{S}_{ij},\mathbf{W}^{T}_{ij}\right)\boldsymbol{\alpha}$.
Much of this methodology reflects the demonstrations provided in Chapter 9 of
Carroll et al., (2006), if more information is desired.

Later on in the analysis stage we will require prior distributions for the
parameter values in the likelihood model of (2). We suppose that exposure
variance components follow conjugate inverse-gamma distributions, i.e.
$\sigma^{2}\sim\mathcal{IG}(r_{\sigma,1},r_{\sigma,2})$,
$\tau^{2}\sim\mathcal{IG}(r_{\tau,1},r_{\tau,2})$, and
$\omega^{2}\sim\mathcal{IG}(r_{\omega,1},r_{\omega,2})$. Each of the mean
parameters for the GPS and EPE models contained in the vectors
$\boldsymbol{\alpha}$ and $\boldsymbol{\beta}$ are assigned a Gaussian prior
with mean zero and variance $b_{\alpha}^{2}$ or $b_{\beta}^{2}$, respectively.

Before we proceed, there is some additional notation that we will need to
define. Let $\mathbf{I}$ denote a diagonal matrix with dimensions matching
either the length of $\mathbf{X}_{i}$ or
$(\tilde{S}_{ij},\mathbf{W}^{T}_{ij})^{T}$, whichever is appropriate in the
context that it appears. Finally, define the vectors
$\mathbf{S}_{i}\equiv(S_{i1},S_{i2},\ldots,S_{iM_{i}})^{T}$,
$\mathbf{A}^{(k)}\equiv(A_{1}^{(k)},A_{2}^{(k)},\ldots,A_{n}^{(k)})^{T}$,
$\bar{\mathbf{Y}}\equiv(\bar{Y}_{1},\bar{Y}_{2},\ldots,\bar{Y}_{n})^{T}$, and
the matrix
$\mathbf{X}=(\mathbf{X}_{1},\mathbf{X}_{2},\ldots,\mathbf{X}_{n})^{T}$.

To begin, initialize $A^{(1)}_{i}$, $\phi^{(1)}_{i}$, and $S^{(1)}_{ij}$ for
all $j=1,2,\ldots,M_{i}$ and $i=1,2,\ldots,n$. We must also initialize the
parameter values $\boldsymbol{\alpha}^{(1)}$, $\boldsymbol{\beta}^{(1)}$,
$\sigma^{(1)}$, $\tau^{(1)}$, $\omega^{(1)}$, $\rho^{(1)}$, $\nu^{(1)}$,
$\mathcal{F}^{(1)}_{1}$, $\mathcal{G}^{(1)}_{1}$, and $\psi^{(1)}$. For
$k=1,2,\ldots,K$ do:

  * •

Imputation Stage

    1. 1.

Sample
$S^{(k)}_{ij}\sim\mathcal{N}\left[\eta^{(k-1)}\left(A^{(k-1)}_{i},\tilde{S}_{ij},\mathbf{W}_{ij}\right),\left(\Sigma^{(k-1)}_{ij}\right)^{2}\right]$
where

$\begin{split}\Sigma^{(k-1)}_{ij}&=\left[\left(\omega^{(k-1)}\right)^{-2}+\left(\tau^{(k-1)}\right)^{-2}\right]^{-1/2}\quad\text{and}\\\
\eta^{(k-1)}\left(A^{(k-1)}_{i},\tilde{S}_{ij},\mathbf{W}_{ij}\right)&=\left(\Sigma^{(k-1)}_{ij}\right)^{2}\times\left\\{\frac{A^{(k-1)}_{i}}{\left(\omega^{(k-1)}\right)^{2}}+\frac{\mu^{(k-1)}_{S}\left(\tilde{S}_{ij},\mathbf{W}_{ij}\right)}{\left(\tau^{(k-1)}\right)^{2}}\right\\}.\end{split}$

Replace known quantities of $S^{(k)}_{ij}$ when $(i,j)\in\mathcal{S}$;

    2. 2.

To sample $A^{(k)}_{i}$, we will need to consider the full conditional
likelihood:

$p_{A}\left(A_{i}\middle|\mathbf{S}_{i},\mathbf{X}_{i},\bar{Y}_{i},\phi^{(k-1)}_{i}\right)\propto
p_{\bar{Y}}(\bar{Y}_{i}|A_{i},\mathbf{X}_{i})\times\exp\left\\{-\frac{1}{2}\left[\frac{A_{i}-\mu_{A}\left(\mathbf{X}_{i},\phi^{(k-1)}_{i}\right)}{\sigma}\right]^{2}\right\\}\times\exp\left[-\frac{1}{2}\sum_{j=1}^{M_{i}}\left(\frac{S_{ij}-A_{i}}{\omega}\right)^{2}\right].$

For all $i=1,2,\ldots,n$, sample $R_{i}\sim\mathcal{U}(0,1)$ and
$A^{+}_{i}\sim\mathcal{N}(A^{(k-1)}_{i},\Delta^{2})$ for some pre-specified
value $\Delta$ that will achieve a sufficiently high acceptance probability.
Set

$A^{(k)}_{i}=\begin{cases}A^{+}_{i},R_{i}\leq\min\left(1,\delta^{(k)}_{i}\right)\\\
A^{(k-1)}_{i},R_{i}>\min\left(1,\delta^{(k)}_{i}\right)\end{cases}$

where

$\delta^{(k)}_{i}=\frac{p^{(k-1)}_{\bar{Y}}\left(\bar{Y}_{i}\middle|A^{+}_{i},\mathbf{X}_{i}\right)\times\exp\left\\{-\frac{1}{2}\left[\frac{A^{+}_{i}-\mu^{(k-1)}_{A}\left(\mathbf{X}_{i},\phi^{(k-1)}_{i}\right)}{\sigma^{(k-1)}}\right]^{2}\right\\}\times\exp\left[-\frac{1}{2}\sum_{j=1}^{M_{i}}\left(\frac{S^{(k)}_{ij}-A^{+}_{i}}{\omega^{(k-1)}}\right)^{2}\right]}{p^{(k-1)}_{\bar{Y}}\left(\bar{Y}_{i}\middle|A^{(k-1)}_{i},\mathbf{X}_{i}\right)\times\exp\left\\{-\frac{1}{2}\left[\frac{A^{(k-1)}_{i}-\mu^{(k-1)}_{A}\left(\mathbf{X}_{i},\phi^{(k-1)}_{i}\right)}{\sigma^{(k-1)}}\right]^{2}\right\\}\times\exp\left[-\frac{1}{2}\sum_{j=1}^{M_{i}}\left(\frac{S^{(k)}_{ij}-A^{(k-1)}_{i}}{\omega^{(k-1)}}\right)^{2}\right]};$

  * •

Analysis Stage

    1. 3.

Using the subset of individuals from $(i,j)\in\mathcal{S}$, sample the
prediction error mean parameters

$\begin{split}\boldsymbol{\alpha}^{(k)}\sim\mathcal{N}\bBigg@{3}(\left\\{\sum_{(i,j)\in\mathcal{S}}\left[\left(\tilde{S}_{ij},\mathbf{W}^{T}_{ij}\right)^{T}\left(\tilde{S}_{ij},\mathbf{W}^{T}_{ij}\right)\right]+\left[\left(\tau^{(k-1)}\right)^{2}\middle/b^{2}_{\alpha}\right]\mathbf{I}\right\\}^{-1}\times\sum_{(i,j)\in\mathcal{S}}\left[S_{ij}\left(\tilde{S}_{ij},\mathbf{W}^{T}_{ij}\right)^{T}\right],\\\
\left(\tau^{(k-1)}\right)^{2}\left[\sum_{(i,j)\in\mathcal{S}}\left[\left(\tilde{S}_{ij},\mathbf{W}^{T}_{ij}\right)^{T}\left(\tilde{S}_{ij},\mathbf{W}^{T}_{ij}\right)\right]+\left[\left(\tau^{(k-1)}\right)^{2}\middle/b^{2}_{\alpha}\right]\mathbf{I}\right\\}^{-1}\bBigg@{3}).\end{split}$

and set
$\mu^{(k)}_{S}\left(\tilde{S}_{ij},\mathbf{W}_{ij}\right)=\left(\tilde{S}_{ij},\mathbf{W}^{T}_{ij}\right)\boldsymbol{\alpha}^{(k)}$;

    2. 4.

Sample the GPS mean parameters

$\phi^{(k)}_{i}\sim\mathcal{N}\left(\frac{\rho^{(k-1)}\sum_{i^{\prime}=1}^{n}V_{ii^{\prime}}\phi_{i^{\prime}}^{(k-1)}}{\rho^{(k-1)}\sum_{i^{\prime}=1}^{n}V_{ii^{\prime}}+1-\rho^{(k-1)}},\frac{\left(\nu^{(k-1)}\right)^{2}}{\rho^{(k-1)}\sum_{i^{\prime}=1}^{n}V_{ii^{\prime}}+1-\rho^{(k-1)}}\right),$

where $V_{ii^{\prime}}$ is the $i$th row and $i^{\prime}$th column of
$\mathbf{V}$, and

$\begin{split}\boldsymbol{\beta}^{(k)}\sim\mathcal{N}\bBigg@{3}(\left\\{\sum_{i=1}^{n}\left(\mathbf{X}_{i}\mathbf{X}^{T}_{i}\right)+\left[\left(\sigma^{(k-1)}\right)^{2}\middle/b^{2}_{\beta}\right]\mathbf{I}\right\\}^{-1}\times\sum_{i=1}^{n}\left[\mathbf{X}_{i}\left(A^{(k)}_{i}-\phi^{(k)}_{i}\right)\right],\\\
\left(\sigma^{(k-1)}\right)^{2}\left\\{\sum_{i=1}^{n}\left(\mathbf{X}_{i}\mathbf{X}^{T}_{i}\right)+\left[\left(\sigma^{(k-1)}\right)^{2}\middle/b^{2}_{\beta}\right]\mathbf{I}\right\\}^{-1}\bBigg@{3}).\end{split}$

Then set
$\mu^{(k)}_{A}\left(\mathbf{X}_{i},\phi^{(k)}_{i}\right)=\phi_{i}^{(k)}+\mathbf{X}^{T}_{i}\boldsymbol{\beta}^{(k)}$;

    3. 5.

Sample the prediction error variance

$\left(\tau^{(k)}\right)^{2}\sim\mathcal{IG}\left\\{r_{\tau,1}+\frac{\left|\mathcal{S}\right|}{2},r_{\tau,2}+\frac{1}{2}\sum_{(i,j)\in\mathcal{S}}\left[S^{(k)}_{ij}-\mu^{(k)}_{S}\left(\tilde{S}_{ij},\mathbf{W}_{ij}\right)\right]^{2}\right\\};$

    4. 6.

Sample the GPS variance

$\left(\sigma^{(k)}\right)^{2}\sim\mathcal{IG}\left\\{r_{\sigma,1}+\frac{n}{2},r_{\sigma,2}+\frac{1}{2}\sum_{i=1}^{n}\left[A^{(k)}_{i}-\mu^{(k)}_{A}\left(\mathbf{X}_{i},\phi^{(k)}_{i}\right)\right]^{2}\right\\};$

    5. 7.

Sample the classical error variance

$\left(\omega^{(k)}\right)^{2}\sim\mathcal{IG}\left[r_{\omega,1}+\frac{m}{2},r_{\omega,2}+\frac{1}{2}\sum_{i=1}^{n}\sum_{j=1}^{M_{i}}\left(S^{(k)}_{ij}-A^{(k)}_{i}\right)^{2}\right];$

    6. 8.

Sample $\rho^{(k)}$ and $\nu^{(k)}$ according to the suggestions of Leroux et
al., (2000) and Lee, (2013), treating $A^{(k)}_{i}$ as the outcome,
$\phi^{(k)}_{i}$ the random effect, and $\mathbf{X}_{i}$ the covariates;

    7. 9.

For tree $t\in\left\\{1,...,T_{k}\right\\}$:

      1. (a)

Sample $\mathcal{F}_{t}^{(k)}$ from the likelihood
$p_{\mathcal{F}}\left(\mathcal{F}_{t}|\mathbf{A}^{(k)},\mathbf{X},\bar{\mathbf{Y}};\psi^{(k-1)}\right)$
via the Metropolis-Hastings search algorithm given by Chipman et al., (1998);

      2. (b)

Sample $\mathcal{G}_{t}^{(k)}$ from the posterior
$p_{\mathcal{G}}\left(\mathcal{G}_{t}|\mathbf{A}^{(k)},\mathbf{X},\bar{\mathbf{Y}};\psi^{(k-1)},\mathbf{\mathcal{F}}_{t}^{(k)}\right)$
using draws from independent Gaussian distributions as described in Chipman et
al., (2010).

    8. 10.

Sample $\psi^{(k)}$ using a conjugate inverse-gamma prior distribution
according to Chipman et al., (2010), supposing $N_{i}$ as the weights.

With the posterior samples of the outcome model,
$\mu^{(k)}_{\bar{Y}}\left(A^{(k)}_{i},\mathbf{X}_{i}\right)$, and the exposure
imputations $A^{(k)}_{i}$, we can proceed to fit the ERF using the local
regression method described in Section S2.

##  S2 Kernel Weighted Least Squares Regression Details

In this supplement, we provide the details of the smoothing step outlined in
Section 3.2. Given a bandwidth $h>0$, the pseudo-outcomes in (4) are regressed
onto fixed points $a\in\mathcal{A}$ by solving for
$\hat{\theta}^{(l)}_{h}(a)=\mathbf{c}^{T}_{ha}(a)\hat{\boldsymbol{\lambda}}^{(l)}_{ha}$
where

$\hat{\boldsymbol{\lambda}}^{(l)}_{ha}=\operatorname*{arg\,min}_{\boldsymbol{\lambda}\in\Re^{2}}\sum_{i=1}^{n}q_{ha}\left(A^{(l)}_{i}\right)\left\\{\xi^{(l)}\left(A^{(l)}_{i},\mathbf{X}_{i},\bar{Y}_{i}\right)-\mathbf{c}^{T}_{ha}\left(A^{(l)}_{i}\right)\boldsymbol{\lambda}\right\\}^{2}$
(5)

and

$\mathbf{c}_{ha}\left(A^{(l)}_{i}\right)=\left[1,h^{-1}\left(A^{(l)}_{i}-a\right)\right]^{T}.$

For a given bandwidth representing the standard deviation of a Gaussian
probability density function centered around $a\in\mathcal{A}$, the kernel
weights are equivalent to

$q_{ha}\left(A^{(l)}_{i}\right)=\frac{N_{i}}{h\sqrt{2\pi}}\exp\left[-\frac{\left(A^{(l)}_{i}-a\right)^{2}}{2h^{2}}\right].$

As we mentioned already in the manuscript, the values
$\hat{\theta}^{(l)}_{h}(a)$ do not form a proper posterior of $\theta(a)$.
Instead, think of the $\xi^{(l)}(A^{(l)}_{i},\mathbf{X}_{i},Y_{i})$ as one of
$L$ imputations for the ERF evaluated at $A^{(l)}_{i}$, meaning we must use
combining rules intended for multiple imputation to estimate the variance of
$\bar{\theta}_{h}(a)$. For each $l=1,2,\ldots,L$, we require the pointwise
variance estimates for each $\hat{\theta}_{h}^{(l)}(a)$. For this we can use a
subtle modification to the locally-weighted least-squares variance estimator
following the form

$\hat{\boldsymbol{\Omega}}^{(l)}_{h}(a)=\sum_{i=1}^{n}\left(\mathbf{D}^{(l)}_{ha}\right)^{-1}\left[\hat{\boldsymbol{\zeta}}^{(l)}_{ha}\left(A^{(l)}_{i},\mathbf{X}_{i},\bar{Y}_{i}\right)\right]^{\otimes
2}\left(\mathbf{D}^{(l)}_{ha}\right)^{-1},$

where
$\mathbf{D}^{(l)}_{ha}=\sum_{i=1}^{n}q_{ha}\left(A^{(l)}_{i}\right)\mathbf{c}_{ha}\left(A^{(l)}_{i}\right)\mathbf{c}^{T}_{ha}\left(A^{(l)}_{i}\right)$
and

$\begin{split}\hat{\boldsymbol{\zeta}}^{(l)}_{ha}\left(A^{(l)}_{i},\mathbf{X}_{i},\bar{Y}_{i}\right)&=q_{ha}\left(A^{(l)}_{i}\right)\left\\{\xi^{(l)}\left(A^{(l)}_{i},\mathbf{X}_{i},\bar{Y}_{i}\right)-\mathbf{c}^{T}_{ha}\left(A^{(l)}_{i}\right)\hat{\boldsymbol{\lambda}}^{(l)}_{ha}\right\\}\mathbf{c}_{ha}\left(A^{(l)}_{i}\right)\\\
&\qquad+\int_{\mathcal{A}}q_{ha}\left(t\right)\left[\mu^{(l)}_{\bar{Y}}(t,\mathbf{X}_{i})-\int_{\mathcal{X}}\mu^{(l)}_{\bar{Y}}(t,\mathbf{x})d\mathbb{P}(\mathbf{x})\right]\mathbf{c}_{ha}\left(t\right)d\mathbb{Q}(t).\end{split}$
(6)

Like the distribution for $\mathbf{X}_{i}$, the cumulative and empirical
density functions for $A_{i}$ are written as $Q(a)$ and $\mathbb{Q}(a)$,
respectively. At this point we can obtain a variance estimate using the
combining rules of Rubin, (2004) with

$\hat{\mathbb{V}}\left[\bar{\theta}_{h}(a)\right]\approx\left(1+L^{-1}\right)(L-1)^{-1}\sum_{l=1}^{L}\left[\hat{\theta}^{(l)}_{h}(a)-\bar{\theta}_{h}(a)\right]^{2}+L^{-1}\sum_{l=1}^{L}\left[\hat{\Omega}^{(l)}_{h,(1,1)}(a)\right]^{2}$

where $\hat{\Omega}^{(l)}_{h,(1,1)}(a)$ is the (1,1) element of
$\hat{\boldsymbol{\Omega}}^{(l)}_{h}(a)$. This approach is ideologically
similar to Antonelli et al., (2020) who also use combining rules after
regressing posterior samples onto the exposure to estimate the ERF.

The bandwidth $h$ can be chosen by cross-validation. This is accomplished by
partitioning $\left[\xi^{(l)}_{i}(\cdot),A^{(l)}_{i}\right]$,
$i=1,2,\ldots,n$, into one of several folds to find $h>0$ that minimizes the
cross-validated mean squared error via grid search.

##  S3 Simulation Details

Sections S3.1 and S3.2 contains a complete description of the data generating
mechanisms for the simulation study design in Section 4.1. We also provide the
additional simulation results when at least one of the measurement error
sources is set to $0$, $\tau^{2}=0$, $\omega^{2}=0$, or both in Section S3.1.
Contained within this table are results that partially align with the
subfigures in Figure 3 (top-row and bottom-left subfigures when $n=800$ and
$m=4000$).

###  S3.1 Correct Specification

For cluster $i=1,2,\ldots,n$, let
$X_{i1},X_{i2},X_{i3},X_{i4}\sim\mathcal{N}(0,1)$,
$\mathbf{X}_{i}=(X_{i1},X_{i2},X_{i3},X_{i4})^{T}$, which we use to generate

$A_{i}\sim\mathcal{N}\left(10+0.5X_{i1}-0.5X_{i2}-0.5X_{i3}+0.5X_{i4},\sigma^{2}\right)$

with $\sigma^{2}=4$. For $j=1,2,\ldots,M_{i}$, we sample
$S_{ij}\sim\mathcal{N}(A_{i},\omega^{2})$ with $\omega^{2}\in\\{0,1,2\\}$. The
associated predictions for $S_{ij}$ are subsequently generated with the
distribution $\tilde{S}_{ij}\sim\mathcal{N}\left(0,\tau^{2}\right)$ where
$\tau^{2}\in\\{0,1,2\\}$. The outcome counts are generated from a Poisson
distribution $Y_{i}\sim\mathcal{P}[\mu_{Y}(A_{i},\mathbf{X}_{i})]$ with

$\begin{split}\log\left[\mu_{Y}(A_{i},\mathbf{X}_{i})\right]&=-3-0.5X_{i1}-0.25X_{i2}+0.25X_{i3}+0.5X_{i4}+\log(N_{i})\\\
&\qquad+0.25(A_{i}-8)-0.75\cos\left[\pi(A_{i}-6)/4\right]-0.25(A_{i}-10)\times
X_{i1}.\end{split}$

In addition, we vary the total number of cells, $m\in\\{2000,5000\\}$ when
$n=400$ and $m\in\\{4000,8000\\}$ when $n=800$. The number of measurements in
each of the $n$ clusters is distributed uniformly across
$M_{i}\in\\{1,2,3,4,6,7,8,9\\}$ when $n=400$ and $m=2000$ or $n=800$ and
$m=4000$. When $n=400$ and $m=4000$ or $n=800$ and $m=8000$, then the number
of measurements is uniformly distributed across
$M_{i}\in\\{2,4,6,8,12,14,16,18\\}$. The offsets are generated from the
uniform distribution $N_{i}\sim\mathcal{U}(10,1{,}000)$. We fix the proportion
of exposures that are within the validation set $\mathcal{S}$ (i.e. when
$S_{ij}$ is “observed”) to $0.1$.

###  S3.2 Incorrect Specification

We would also like to examine the potential sources of bias and error that
might occur when we misspecify each of the three models in (2) - the outcome
model, the generalized propensity score (GPS), and the error prone exposure
model (EPE). Under correct specification, the data are generated using the the
models described above in Section S3.1. Using the same definitions for
$\mathbf{X}_{i}$ as above, define $\tilde{X}_{i1}=\exp\left(X_{i1}/2\right)$,
$\tilde{X}_{i2}=X_{i2}(1+X_{i1})^{-1}+10$,
$\tilde{X}_{i3}=(X_{i1}X_{i3}/25+0.6)^{3}$, and
$\tilde{X}_{i4}=(X_{i2}+X_{i4}+10)^{2}$ (Kang and Schafer,, 2007). The
transformations
$\tilde{\mathbf{X}}_{i}=(\tilde{X}_{i1},\tilde{X}_{i2},\tilde{X}_{i3},\tilde{X}_{i4})^{T}$
are subsequently scaled and centered to have a mean of zero and a marginal
variance of one. Under a misspecified GPS scenario we generate

$A_{i}\sim\mathcal{N}\left[10+0.5\tilde{X}_{i1}-0.5\tilde{X}_{i2}-0.5\tilde{X}_{i3}+0.5\tilde{X}_{i4},\sigma^{2}\right],$

yet the GPS models are fit using the original covariates $\mathbf{X}_{i}$.
Likewise, scenarios where the outcome model is misspecified means we generate
$Y_{i}\sim\mathcal{P}\left[\tilde{\mu}_{Y}(A_{i},\mathbf{X}_{i})\right]$ with

$\begin{split}\log\left[\tilde{\mu}_{Y}(A_{i},\mathbf{X}_{i})\right]&=-3-0.5\tilde{X}_{i1}-0.25\tilde{X}_{i2}+0.25\tilde{X}_{i3}+0.5\tilde{X}_{i4}+\log(N_{i})\\\
&\qquad+0.25(A_{i}-8)-0.75\cos\left[\pi(A_{i}-6)/4\right]-0.25(A_{i}-10)\times\tilde{X}_{i1}.\end{split}$

However, we continue to fit the outcome models with the original (i.e.
untransformed) covariates similar to the scenarios where the GPS is
misspecified. Finally, when the EPE model is misspecified, then we introduce
prediction bias, generating
$\tilde{S}_{ij}\sim\mathcal{N}\left[S_{ij}-1+0.5W_{ij1}+0.5W_{ij2},\tau^{2}\right]$
where $W_{ij1}\sim\mathcal{N}(1,2)$ and $W_{ij2}\sim\mathcal{N}(X_{i2},1)$. In
this scenario, the naïve approach will produce biased results since there is
no adjustment for the added prediction bias. In the regression calibration and
multiple imputation approaches, however, this bias is corrected for when
addressing the prediction error, and therefore the downstream estimates of the
ERF should be unbiased when the other models are well approximated. For these
misspecification scenarios, we fix $\sigma^{2}=4$, $\omega^{2}=1$,
$\tau^{2}=0.5$, $n=400$, and $m=4000$.

##  S4 Congeniality Simulation

In this supplement, we run a small simulation study to demonstrate how using a
GPS-based estimator has almost no influence over the accuracy of the ERF
estimate while adhering to rules of congeniality (Meng,, 1994). As we
mentioned in Section 3.3, the requirements of congeniality seem to conflict
with the principles of cutting feedback in causal models fit with Bayesian
methods (Zigler et al.,, 2013). To demonstrate how the feedback created by
congeniality can affect estimators that condition on the GPS, we will need to
incorporate the GPS into our pseudo-outcome. This implementation was described
by Kennedy et al., (2017), who specify the following pseudo-outcome which we
have adapted slightly to fit into the multiple imputation framework:

$\tilde{\xi}^{(l)}\left(A^{(l)}_{i},\mathbf{X}_{i},\bar{Y}_{i}\right)=\frac{\left[\bar{Y}_{i}-\mu^{(l)}_{Y}\left(A^{(l)}_{i},\mathbf{X}_{i}\right)\right]}{p^{(l)}_{A}\left(A^{(l)}_{i}\middle|\mathbf{X}_{i},\phi_{i}^{(l)}\right)}\int_{\mathcal{X}}p^{(l)}_{A}\left(A^{(l)}_{i}\middle|\mathbf{x},\phi_{i}^{(l)}\right)d\mathbb{P}(\mathbf{x})+\int_{\mathcal{X}}\mu^{(l)}_{Y}\left(A^{(l)}_{i},\mathbf{x}\right)d\mathbb{P}(\mathbf{x})$
(7)

The first part of this pseudo-outcome, on the left-hand side of the addition
symbol, characterizes an estimator of residual error conditioned on the
stabilized inverse probability weights constructed from the posterior samples
of the GPS. If the outcome model is correctly specified, then the average of
these debiased residuals will approach zero as $n$ approaches infinity. The
second part of (7), on the right-hand side of the addition symbol, is the
marginalized outcome model. When the outcome model is misspecified, then the
estimator of the residual error will counteract the asymptotic bias generated
by the marginalized outcome model, thus allowing for an alternative means to
debias the estimate so long as the GPS is correctly specified. Thus, this
pseudo-outcome facilitates doubly-robust estimation when used in tandem with
locally-weighted regression methods, at least in cases where imputing the
exposure is not necessary.

Consider the simulation scenarios in Section S3, and in particular the
misspecification scenarios in Section S3.2. We will once again fit the model
in (2), but substitute $\xi^{(l)}(\cdot)$ with $\tilde{\xi}^{(l)}(\cdot)$.
Moreover, in an effort to better examine the doubly-robust properties of the
pseudo-outcome in (7), which can be masked by the flexibility of a BART model,
we will fit the likelihood in (2) assuming
$p_{\bar{Y}}(\bar{Y}_{i}|A_{i},\mathbf{X}_{i})$ is a Poisson model,
$Y_{i}\sim\mathcal{P}\left[\mu_{Y}(A_{i},\mathbf{X}_{i})\right]$, where

$\log[\mu_{Y}(A_{i},\mathbf{X}_{i})]=\gamma_{0}+\gamma_{1}(A_{i}-10)+\gamma_{2}\cos\left[\pi(A_{i}-6)/4\right]+\gamma_{3}(A_{i}-10)X_{i1}+\boldsymbol{\gamma}_{4}^{T}\mathbf{X}_{i}+\log(N_{i}).$
(8)

Note that the simulated outcome data without misspecification is generated
using the same model, where $\gamma_{0}$–$\boldsymbol{\gamma}_{4}$ are known.
Otherwise, the smoothing steps outlined in Section S2 remain unchanged other
than to substitute
$\mu^{(l)}_{\bar{Y}}\left(A^{(l)}_{i},\mathbf{X}_{i}\right)$ with the form in
(8). The Gaussian models that we assumed for $p_{A}(A_{i}|\mathbf{X}_{i})$
with $\phi_{i}=0$ for all $i=1,2,\ldots,n$, $p_{S}(S_{ij}|A_{i})$, and
$p_{S}(S_{ij}|\tilde{S}_{ij},\mathbf{W}_{ij})$ in Section 4 are also
implemented in the same fashion here. The algorithm described in Section S1 is
used to generate posterior draws of the latent exposures and unknown parameter
values.

We are most interested in the scenarios where the outcome model is
misspecified and the GPS is correctly specified. For an estimator to be
doubly-robust, it must produce unbiased estimates in these scenarios. There
are two curves for every subfigure in Figure S1, one where we use
$\xi^{(l)}(\cdot)$ and the other where we use $\tilde{\xi}^{(l)}(\cdot)$. Both
curves use Gaussian kernel weighted least squares regression to project the
respective pseudo-outcomes onto $a\in\mathcal{A}$. Examining the scenarios
where the GPS is correctly specified but the outcome model is misspecified, we
can see that the two curves are nearly identical despite one using a doubly-
robust form. Therefore, the utility of a doubly-robust estimator in this
framework yields little to no benefit. The observed bias that we detect when
the outcome model is misspecified and the GPS model is correctly specified, we
contend, is due to the feedback created between the GPS and the outcome model
from the requirements of congeniality to support multiple imputation.

The above analysis guided the development of our framework in the main
manuscript. In an effort to limit any potential deleterious impacts created by
feedback that may or may not exist given our limited testing, we decided to
refrain from incorporating the GPS into our estimator of the ERF, even though
we lose the possibility of achieving the doubly-robust property. However, from
this simulation study, we did not see any benefit to using a doubly-robust
estimator since it would seem that the outcome model needs to be correctly
specified anyway. To make up for this loss, we decided to incorporate BART
into our estimator of the ERF. BART is a nonparametric approach that provides
added flexibility in cases of outcome model misspecification. We concede that
more work is required to better understand how to ameliorate the feedback
issue for GPS-based estimators of the ERF with a latent exposure variable.

## Tables and Figures

Figure S1: ERF estimates averaged over $500$ simulated iterations when $n=800$ and $M=4000$ under various model misspecification conditions. $n$ | $m$ | $\omega^{2}$ | $\tau^{2}$ |  | Relative Bias |  | 95% Coverage Probability
---|---|---|---|---|---|---|---
| No Correction |  | Regression
---
Calibration
| BART Multiple
---
Imputation
| GLM Multiple
---
Imputation
| No Correction |  | Regression
---
Calibration
| BART Multiple
---
Imputation
| GLM Multiple
---
Imputation
400 | 4,000 | 0 | 0 |  | 0.19 (-0.03) | 0.20 (-0.03) | 0.20 (-0.03) | 0.02 (-0.02) |  | 0.65 (0.75) | 0.65 (0.72) | 0.75 (0.80) | 0.90 (0.87)
400 | 4,000 | 0 | 1 |  | 0.31 (-0.04) | 0.23 (-0.03) | 0.25 (-0.02) | 0.00 (-0.02) |  | 0.56 (0.64) | 0.66 (0.67) | 0.79 (0.78) | 0.91 (0.88)
400 | 4,000 | 0 | 2 |  | 0.26 (-0.07) | 0.14 (-0.07) | 0.15 (-0.05) | 0.00 (-0.02) |  | 0.51 (0.64) | 0.63 (0.69) | 0.80 (0.78) | 0.90 (0.88)
400 | 4,000 | 1 | 0 |  | 0.20 (-0.06) | 0.21 (-0.06) | 0.18 (-0.05) | 0.03 (-0.02) |  | 0.57 (0.70) | 0.57 (0.70) | 0.74 (0.76) | 0.93 (0.92)
400 | 4,000 | 2 | 0 |  | 0.24 (-0.07) | 0.24 (-0.08) | 0.20 (-0.06) | 0.02 (-0.02) |  | 0.50 (0.62) | 0.50 (0.62) | 0.69 (0.75) | 0.93 (0.90)
400 | 2,000 | 0 | 0 |  | 0.34 (0.01) | 0.35 (0.01) | 0.39 (0.02) | 0.02 (-0.02) |  | 0.64 (0.75) | 0.64 (0.76) | 0.74 (0.78) | 0.91 (0.90)
400 | 2,000 | 0 | 1 |  | 0.27 (-0.06) | 0.17 (-0.06) | 0.19 (-0.04) | 0.00 (-0.02) |  | 0.53 (0.66) | 0.64 (0.66) | 0.76 (0.75) | 0.88 (0.92)
400 | 2,000 | 0 | 2 |  | 0.31 (-0.10) | 0.18 (-0.08) | 0.19 (-0.05) | -0.01 (-0.02) |  | 0.42 (0.48) | 0.59 (0.54) | 0.76 (0.71) | 0.85 (0.86)
400 | 2,000 | 1 | 0 |  | 0.25 (-0.08) | 0.24 (-0.07) | 0.19 (-0.06) | 0.01 (-0.02) |  | 0.49 (0.58) | 0.49 (0.62) | 0.69 (0.73) | 0.93 (0.90)
400 | 2,000 | 2 | 0 |  | 0.35 (-0.09) | 0.35 (-0.09) | 0.24 (-0.07) | 0.02 (-0.02) |  | 0.44 (0.47) | 0.45 (0.52) | 0.67 (0.72) | 0.92 (0.89)
800 | 8,000 | 0 | 0 |  | 0.09 (-0.02) | 0.10 (-0.02) | 0.10 (-0.02) | 0.00 (-0.01) |  | 0.69 (0.81) | 0.68 (0.80) | 0.78 (0.84) | 0.93 (0.94)
800 | 8,000 | 0 | 1 |  | 0.13 (-0.04) | 0.08 (-0.03) | 0.10 (-0.02) | -0.01 (-0.01) |  | 0.57 (0.76) | 0.68 (0.78) | 0.79 (0.81) | 0.88 (0.94)
800 | 8,000 | 0 | 2 |  | 0.18 (-0.06) | 0.10 (-0.05) | 0.10 (-0.02) | -0.01 (-0.01) |  | 0.49 (0.62) | 0.55 (0.70) | 0.78 (0.84) | 0.87 (0.92)
800 | 8,000 | 1 | 0 |  | 0.16 (-0.04) | 0.15 (-0.04) | 0.12 (-0.02) | 0.01 (-0.01) |  | 0.57 (0.74) | 0.56 (0.74) | 0.75 (0.84) | 0.97 (0.96)
800 | 8,000 | 2 | 0 |  | 0.21 (-0.05) | 0.21 (-0.05) | 0.15 (-0.03) | 0.01 (-0.01) |  | 0.50 (0.62) | 0.52 (0.68) | 0.72 (0.80) | 0.97 (0.96)
800 | 4,000 | 0 | 0 |  | 0.11 (-0.02) | 0.11 (-0.01) | 0.12 (-0.02) | 0.01 (0.00) |  | 0.69 (0.82) | 0.70 (0.84) | 0.79 (0.90) | 0.94 (0.94)
800 | 4,000 | 0 | 1 |  | 0.21 (-0.05) | 0.14 (-0.04) | 0.13 (-0.03) | -0.02 (0.00) |  | 0.50 (0.62) | 0.65 (0.67) | 0.76 (0.76) | 0.87 (0.93)
800 | 4,000 | 0 | 2 |  | 0.28 (-0.09) | 0.18 (-0.07) | 0.11 (-0.03) | -0.03 (-0.01) |  | 0.39 (0.42) | 0.55 (0.56) | 0.74 (0.76) | 0.76 (0.94)
800 | 4,000 | 1 | 0 |  | 0.19 (-0.06) | 0.18 (-0.06) | 0.13 (-0.04) | 0.01 (-0.01) |  | 0.48 (0.57) | 0.48 (0.56) | 0.72 (0.76) | 0.96 (0.94)
800 | 4,000 | 2 | 0 |  | 0.28 (-0.09) | 0.28 (-0.08) | 0.17 (-0.04) | 0.01 (-0.01) |  | 0.42 (0.45) | 0.42 (0.48) | 0.69 (0.76) | 0.96 (0.94)
Table S1: Relative bias and coverage probabilities averaged over the simulated
ERF estimates when either $\tau^{2}=0$, $\omega^{2}=0$, or both. The values in
parentheses represent these same statistics evaluated at $a=11$.