# Bayesian Methods for Modeling Cumulative Exposure to Extensive Environmental
Health Hazards

Rob<EMAIL_ADDRESS>[ Jesse
<EMAIL_ADDRESS>[ Jon
<EMAIL_ADDRESS>[ Joseph N.S.
<EMAIL_ADDRESS>[ Yang
<EMAIL_ADDRESS>[ Department of Statistics,
Oregon State Universitypresep=, ]e1  Department of Epidemiology , School of
Public Health University of Michigan, Ann Arborpresep=, ]e2,e3,e4  Department
of Statistics, University of Michigan, Ann Arborpresep=, ]e5

(2024)

###### Abstract

Measuring the impact of an environmental point source exposure on the risk of
disease, like cancer or childhood asthma, is well-developed. Modeling how an
environmental health hazard that is extensive in space, like a wastewater
canal, impacts disease risk is not well-developed. We propose a novel Bayesian
generative semiparametric model for characterizing the cumulative spatial
exposure to an environmental health hazard that is not well-represented by a
single point in space. The model couples a dose-response model with a log-
Gaussian Cox process integrated against a distance kernel with an unknown
length-scale. We show that this model is a well-defined Bayesian inverse
model, namely that the posterior exists under a Gaussian process prior for the
log-intensity of exposure, and that a simple integral approximation adequately
controls the computational error. We quantify the finite-sample properties and
the computational tractability of the discretization scheme in a simulation
study. Finally, we apply the model to survey data on household risk of
childhood diarrheal illness from exposure to a system of wastewater canals in
Mezquital Valley, Mexico.

00X00,

00X00,

00X00,

Cumulative Exposure,

Bayesian Model,

Cox Process,

Gaussian Process,

Childhood Diarrhea,

###### keywords:

[class=MSC]

###### keywords:

††volume: TBA††issue: TBA

2010.00000

, and

##  1 Introduction

To measure the impact of environmental exposures, epidemiologists must account
for both the intensity of exposure from the environmental source of interest,
factors impacting individual susceptibility to infection, and other unobserved
potential sources of infection. Epidemiologists often employ regression
modeling to learn the relationship between exposure, susceptibility and
disease risk [2]. When a point source exposure to pathogens or other health
hazard is suspected, methods have been developed to incorporate these sources
of disease into a regression modeling framework [15, 17]. Exposure to the
environmental source is typically operationalized as a function of distance to
the point source [16], and the modeler learns how quickly the magnitude of the
exposure increases as the distance to the point source decreases. For example,
if there is an abnormal clustering of cancer cases near a chemical plant, we
may suspect that this cluster arose in part due to exposure to hazardous
chemicals.

The assumption that legitimizes the use of these methods is that the location
of the potential exposure is certain: the origin of the cancer risk posed by
the chemical plant can reasonably be assumed to originate from a single,
discrete point in space. Thus, while there may be uncertainty in the
parameters governing the exposure as a function of distance to the point
source, the modeler assumes that there is no uncertainty in the distance for
each unit of analysis. This method has been used to quantify the risk of
larynx cancer with respect to distance to an industrial incinerator [16], the
risks of various cancers in relation to petrochemical plant exposure, [9], the
risk of multi-drug-resistant tuberculosis infection for individuals living
near a prison in Peru, [37], and to understand the risks of fast food
restaurant proximity on childhood obesity [27].

When the environmental source of disease is spatially extensive, such as a
river, lake, canal network, or pipe system, the fundamental assumption of the
point-source approach, i.e. that each unit’s exposure can be summarized using
a single distance from the source, no longer holds. There is no longer a
plausible single point to which exposure can be assigned; all points that
comprise the source could pose a hazard to health. In order to use the point-
source method for these types of sources, we need to assume that an
individual’s exposure can still be summarized by their distance to a single
location along the source. The typical assumption is to take the shortest
distance [11, 24]. This assumption will necessarily understate the uncertainty
in exposure because many points along the source may contribute to exposure at
a given community location. These methods also do not allow for variation in
risk concentration along the environmental hazard, which in certain
applications may not accord with reality. For instance when modeling the risk
of diarrheal illness with respect to wastewater runoff, enteric pathogen
concentrations decline with distance to the source of runoff [7]. Thus units
near areas of the waterway that are closer to the wastewater runoff will have
higher exposure compared to units that are far away from runoff. This can
result in biased estimates of unit exposure, whereby average intensity is
correctly estimated, but the risk for high-exposure units is underestimated
and vice-versa for low-exposure units.

Instead, we need to take into account cumulative exposure to the environmental
hazard, where every point contributes to exposure. The approach has several
analogues in public health literature. The first analogue is in determining
exposure to fine particulate matter ($\text{PM}_{2.5}$) and its effect on
birth weight in [3]. One measurement of exposure proposed was a summation over
a fixed window of time of daily $\text{PM}_{2.5}$ minus a predefined
threshold, only on days that $\text{PM}_{2.5}$ exceeded the threshold. This
measure accounts for the cumulative exposure to $\text{PM}_{2.5}$. When
considering exposure to an environmental hazard, instead of summing over time,
we can sum exposure over the spatial extent of the hazard. Another analogous
technique can be seen in joint models of longitudinal outcomes and time-to-
event data, summarized in [23] and further generalized in [6], where we model
the time-dependent hazard of an event dependent on parameters learned from the
model for longitudinal outcomes. The form of the interdependence between the
time-to-event model and the longitudinal outcome is specified by the modeler.
One formulation of the joint model specifies that the log-hazard ratio of the
time-to-event model at time $t$ depends on integral over the interval $[0,t]$
of the latent parameter governing the longitudinal outcome [1], which allows
the event hazard to depend on the cumulative exposure to the latent process.
This corresponds to our problem setting, but substituting space for time.

The problem of assigning risk to a spatially-continuous source of disease is
widespread, and the applications of a coherent model for these scenarios are
myriad. The most direct application is in modeling how infectious disease risk
depends on household and occupational proximity to waterways. In multiple
settings, diarrheal illness risk has been observed to cluster around rivers
[35], canal systems, and other water sources. Legionella is also known to
spread through water and rivers provide one route to infection [11]. In order
to develop an effective response to these public health risks, authorities
would benefit from a detailed understanding of the intensity of disease risk
at different points along waterways. This could aid in determining points at
which to sample water quality, and predictive modeling could suggest groups of
households to inform about the health risks. Beyond infectious disease,
childhood respiratory diseases such as asthma and bronchitis have been linked
to vehicle emissions [26]. Learning how respiratory disease risk changes as a
function of cumulative exposure to freeway traffic, as well as how this
exposure varies along freeway segments would be a boon to public health in
this context as well. Officials could target air pollution mitigation efforts
at areas where emissions concentrate and to which households have high
exposure. Urban planners could use model predictions to build new housing
safely away from the worst highways.

We propose a flexible, generative model for quantifying environmental
exposures that naturally extends the point-source approach to account for
spatially extensive sources of risk. The model addresses two key problems in
the measurement of exposure from spatially extensive point sources. The first
is the problem of uncertainty in the point of exposure, e.g. for a ‘single
hit’ model in which the disease outcome is the result of a single exposure
(e.g. infection), and the second in which it is impacted by the accumulation
of exposure over space and/or time. To accomplish this, we allow each unit’s
exposure to be integrated across the entirety of the environmental hazard. We
use a log-Gaussian process to parameterize the risk at distance zero to the
source of risk, which accounts for differences in risk at distinct points
along the source. We use Bayesian inference implemented in the software
package CmdStan to estimate the unknown model parameters [10], where we’re
able to take advantage of parallel computing of the likelihood. We demonstrate
the model’s ability to infer environmental exposure under different data
generating processes using simulated data, and, finally, we apply our model to
data collected on childhood diarrheal disease in Mezquital Valley, Mexico to
show how the model yields new insights that cannot be obtained using existing
methods.

##  2 Modeling environmental exposure

Suppose we observe outcome data $Y_{it}$, where $i$ indexes the individual, or
stratum of the population, and $t$ designates a time interval,
$t=[t_{1},t_{2}]$. Typically $Y_{it}$ is discrete, either representing counts
of incident cases of a disease when $i$ represents a stratum, or an indicator
variable representing the event that individual $i$ is infected with the
disease of interest. These cases are observed in a spatial domain
$\mathcal{R}$ and each observed unit is associated with a location
$s_{i}\in\mathcal{R}$. The spatial domain also contains a potential
environmental source of disease, defined as a set $\mathcal{C}$, which may be
modeled as a lower-dimensional manifold of $\mathcal{R}$, with associated
coordinate function $\ell:\mathcal{C}\to\mathcal{R}$. Then the spatial domain
of the environmental hazard is
$\ell\left(\mathcal{C}\right)\subseteq\mathcal{R}$.

How can we go about learning how the location of the unit $i$ with respect to
the set of points $\mathcal{C}$ influences the risk of disease? We might fit a
model to the observations, by specifying a distribution $P$ for observations
dependent on a parameter $\mu$ specific to each location $i$ and interval $t$
and nuisance parameters $\delta$:

$\displaystyle Y_{it}$ $\displaystyle\sim P(\mu_{it},\delta).$ (1)

We then model how $\mu_{it}$ depends on unit location $s_{i}$ with respect to
the hazard $\mathcal{C}$, and also potentially on observed covariates, through
a known function $h$ and unknown parameters $\boldsymbol{\theta}$:

$\mu_{it}=h(s_{i},\mathcal{C}\mid\boldsymbol{\theta},X_{it}).$

The parameters are identifiable from the observed data via differential
exposure to the hazard between units that vary in location. If we can learn
these parameters, we’ll know to what extent the environmental hazard endangers
those exposed.

##  3 Existing approaches to modeling environmental exposure

When the risk factor is a point-source the domain $\mathcal{C}$ comprises a
single point, $c$, exposure can be approximated via a function,
$\mathcal{K}(d)$, of the scaled Euclidean distance of location $s_{i}$ to the
location of the point $c$, given by the function $\ell(c)$:
$d_{i}=\mathcal{\|}s_{i}-\ell(c)\|_{2}$. For instance, if $Y_{i}$ are binary
indicators of disease, we would model the outcome as a Bernoulli random
variable with mean parameter $\mu_{it}$:

$Y_{it}\sim\text{Bernoulli}(\mu_{it}),$

and model $h^{-1}(\mu_{it})$ as an additive decomposition of baseline risk of
disease $\lambda$ and the increased risk from the hazard,
$f(t)\mathcal{K}(d_{i}/\rho)$. If we define $\mathcal{K}(d)$ to be a strictly
monotone decreasing function that evaluates to $1$ at $d=0$ and tends to $0$
as $d\to\infty$, then parameter $\rho$ quantifies how quickly the risk of
disease decreases as one moves away from the hazard $c$, and the parameter
$f(t)$ gives the increased risk of disease at distance zero to the source:

$\displaystyle h^{-1}(\mu_{it})=\lambda+f(t)\mathcal{K}(d_{i}/\rho),$ (2)

We impose the constraint that $\lambda,f(t)>0\forall t$, and $h^{-1}$ is a
link function. For a fixed distance, $\mathcal{K}$ is increasing in $\rho$.
Some or all of the parameters that govern $h^{-1}(\mu_{it})$, $\rho$,
$\lambda$, and $f(t)$ will typically be unknown and will need to be learned
from observed data. Popular choices for $\mathcal{K}$ are the kernel of a
Gaussian density and the kernel of an exponential density [37, 17].

[16] introduced the model for noninfectious disease case-control studies,
namely cancer. Two independent nonhomogeneous Poisson processes with intensity
functions $\lambda_{y_{i}=1}(s_{i})$ and $\lambda_{y_{i}=0}(s_{i})$, are
assumed to govern the spatial point pattern of cases and controls
respectively. If we define the intensity function for cases
$\lambda_{y_{i}=1}(s_{i})=r\lambda_{y_{i}=0}(s_{i})(1+\alpha\exp(-\mathcal{\|}s_{i}-\ell(c)\|_{2}^{2}/\rho^{2}))$,
then $\alpha$ represents the incremental intensity for cases at distance 0 to
the environmental hazard compared to the base rate of disease $r$. As
$d_{i}\to\infty$ the risk approaches $r\lambda_{y_{i}=0}(s_{i})$. Inference
proceeds by first estimating $\lambda_{y_{i}=0}(s_{i})$, and subsequently
estimating $r$, $\alpha$ and $\rho$. A later paper [17] shows that inference
of the nonhomogeneous Poisson process intensity function can be avoided by
conditioning on the locations of the cases and controls and fitting a
nonlinear binary regression. The modeled outcome is the binary event that an
observation is a disease case conditional on the location of the observation.
The probability of a case is $\mu_{i}$:

$\displaystyle\mu_{i}=\frac{r(1+\alpha\exp(-\mathcal{\|}s_{i}-\ell(c)\|_{2}^{2}/\rho^{2}))}{1+r(1+\alpha\exp(-\mathcal{\|}s_{i}-\ell(c)\|_{2}^{2}/\rho^{2}))}.$
(3)

The model can be derived by splitting a single Poisson process with intensity

$r\lambda_{y_{i}=0}(s_{i})(1+\alpha\exp(-\mathcal{\|}s_{i}-\ell(c)\|_{2}^{2}/\rho^{2}))+\lambda_{y_{i}=0}(s_{i})$

into cases and controls with a location-dependent thinning probability
$\mu_{i}$. Then $P(Y_{i}=1\mid s_{i})=\mu_{i}$.

The model in [17] is a special case of equation (2) and corresponds to
modeling the odds of disease, $h^{-1}$ as

$h^{-1}(\mu)=\frac{\mu}{1-\mu},$

and

$h^{-1}(\mu_{i})=r+r\alpha\exp(-\mathcal{\|}s_{i}-\ell(c)\|_{2}^{2}/\rho^{2})$
(4)

$\lambda=r$, and $f(t)=r\alpha$. This equation can be rearranged to show that
the environmental exposure multiplies the base odds of disease:

$h^{-1}(\mu_{i})=r(1+\alpha\exp(-\mathcal{\|}s_{i}-\ell(c)\|_{2}^{2}/\rho^{2}))$
(5)

When modeling the effects of multiple environmental hazards on a disease of
interest, the odds model in equation (5) can be extended either
multiplicatively:

$\displaystyle
h^{-1}(\mu_{i})=r\prod_{q=1}^{Q}\,(1+\alpha_{q}\exp(-\mathcal{\|}s_{i}-\ell(c)\|_{2}^{2}/\rho_{q}^{2})),$
(6)

as in [17, 29], or additively:

$\displaystyle
h^{-1}(\mu_{i})=r\,(1+\sum_{q=1}^{Q}\alpha_{q}\exp(-\mathcal{\|}s_{i}-\ell(c)\|_{2}^{2}/\rho_{q}^{2})),$
(7)

as in [5], though the authors note the inferential issues pertaining to
correlated hazards.

###  3.1 Extending the model to infectious disease

The point-source exposure model has been successfully employed in infectious
disease settings too. [37] models the probability of multi-drug-resistant
tuberculosis (MDR-TB) with respect to distance to a prison. In this case, a
probit model is used, corresponding to $h^{-1}(\mu)\equiv\Phi^{-1}(\mu)$, or
the inverse cumulative distribution function of the standard normal
distribution.

The key difference between the model used in [37] and [17] is that in [37] no
assumption is made about the independence of cases and noncases of MDR-TB. In
[17] the cases and controls are modeled to have arisen from independent
nonhomogeneous Poisson processes, which is not a valid assumption when
modeling infectious disease.

###  3.2 Extensive environmental hazards

When the environmental exposure is not a point source, but is instead
extensive in space, such as a river, or a lake, modelers typically resort to
using the shortest distance to the source as a proxy for exposure [11]. If
$\mathcal{C}$ is the set of all points comprising the source, then we can
apply the same model above, defined in equation (2), but with $d_{i}$ defined
as:

$\displaystyle d_{i}$
$\displaystyle=\min_{c\in\mathcal{C}}\mathcal{\|}s_{i}-\ell(c)\|_{2}$ (8)

where $\ell$ maps an element of $\mathcal{C}$ to its spatial coordinates. This
accounting doesn’t capture the true extent of the unit’s exposure, however,
because the unit is exposed to the entirety of the environmental hazard. This
method also imposes the constraint that all points on the canal impart equal
exposure at distance zero, which doesn’t allow the exposure to change as a
function of $c$.

##  4 Cumulative exposure to extensive environmental hazards

As discussed in section 1, not accounting for each unit’s cumulative exposure
to the environmental hazard can introduce bias in estimating the risk of
disease for low- and high-exposure units. In order to ameliorate this
deficiency, we can extend the methodology introduced in section 3 by
partitioning the hazard into mutually exclusive subsets $C_{m}$,
$m\in\\{1,\dots,M\\}$, with spatial centroids $\bar{C}_{m}$, and spatial areas
$\Delta(C_{m})$, such that $\bigcup_{m=1}^{M}C_{m}=\mathcal{C}$. Then we could
treat each section $C_{m}$ of the hazard as a separate point source:

$\displaystyle
h^{-1}(\mu_{i})=r\,\left(1+\sum_{m=1}^{M}f_{m}\Delta(C_{m})\mathcal{K}(\mathcal{\|}s_{i}-\ell(\bar{C}_{m})\|_{2}/\rho)\right).$
(9)

Further, we can jointly model the dependence between different subsets by
specifying a joint distribution for the vector of exposures, $\mathbf{f}$,

$\displaystyle(f_{1},f_{2},\dots,f_{M})\sim P_{f}.$ (10)

This model is nearly a direct extension of model (2). Equation 9 accounts for
cumulative exposure while allowing for the concentration of the disease-
causing agent to change along the hazard.

In infectious disease modeling, we often model how risk depends additively on
exposures [14], rather than multiplicatively, like in (9) and (6). We can
change the model to allow for an additive relationship between the base rate
$\lambda$ and the environmental exposure $\mathcal{C}$:

$\displaystyle h^{-1}(\mu_{i})$
$\displaystyle=\lambda+\sum_{m=1}^{M}f_{m}\Delta(C_{m})\mathcal{K}(\mathcal{\|}s_{i}-\ell(\bar{C}_{m})\|_{2}/\rho).$
(11)

In section 5 we will show how an additive decomposition of risk arises
naturally from a generative model for infection and environmental exposure to
a point source hazard. Then in section 5.3 we will show how equation (11) can
be derived for exposure to an extensive hazard when pathogens along the hazard
are distributed according to a nonhomogeneous Poisson process. Finally, we
will specify a nonparametric model for $P_{f}$ in line (10), which further
extends the model to the generative scenario where pathogens along the source
are distributed according to a log-Gaussian Cox process.

##  5 A new perspective on environmental exposure

Researchers and statisticians often find value in fitting generative models to
observed data because these models can tell coherent probabilistic stories
about how the data arose, and allow researchers to encode scientific
information about the modeled phenomenon through distributions. This then
enables probabilistic model checking, which can yield new models or suggest
new datasets to collect that would refine or enhance the scientific insights.
In the context of environmental health, generative models also allow us to
examine counterfactual scenarios in order to estimate, e.g. the proportion of
observed disease risk attributable to a given exposure. These and other
benefits of generative modeling led us to recast the environmental hazard
model as a natural consequence of a specific probabilistic story about how
environmental disease data came about. This recasting allows us to extend the
point-source model via an expanded generative process that yields a more
detailed picture of environmental exposure.

###  5.1 Dose-response model

Our treatment begins with a model for disease called the exponential dose-
response model. It is a generative model for infectious disease and, as such,
can be used as a building block in a more realistic model of infection from a
point-source exposure. In quantitative microbial risk assessment, the model is
typically used to infer the dose of a pathogen that would lead to a 50%
probability of infection or symptomatic disease.

Let a time interval beginning at time $0$, and ending at time $T$ be denoted
$[0,T]$. Let the nonhomogeneous Poisson process rate $N_{i}(T)$ be the number
of disease-causing pathogens to which a unit $i$ is exposed over the time
interval $[0,T]$. Let the cumulative intensity be defined as
$\Lambda(T)=\int_{0}^{T}\lambda(t)dt$:

$\displaystyle N_{i}(T)\sim\text{Poisson}(\Lambda(T)).$

Subsequently, let the number of pathogens infecting unit $i$ over $[0,T]$,
$K_{i}(T)\mid N_{i}(T)$, be conditionally binomial with probability of success
parameter $r_{i}$:

$\displaystyle K_{i}(T)\mid N_{i}(T)$
$\displaystyle\sim\text{Binomial}(N(T),r_{i}).$

Marginalizing over $N_{i}(T)$ yields

$\displaystyle K_{i}(T)$ $\displaystyle\sim\text{Poisson}(r_{i}\Lambda(T)).$
(12)

Then the probability that unit $i$ becomes infected over the time interval is
$1-\exp(-r_{i}\Lambda(T))$, or the probability that at least one of these
pathogens infects person $i$. The parameter $r_{i}$ represents the
susceptibility of unit $i$ to the disease and may depend on covariates
$X_{i}$. See [8] for more discussion on the merits of the one-hit does-
response model.

Given that $N_{i}(T)$ is a Poisson process on $\R^{+}$, $K_{i}(T)$ is a
thinned Poisson process on $\R^{+}$ with mean measure $r_{i}\lambda(u)$. The
probability of $i$ becoming infected over the interval $[0,T]$ is the
probability of observing the first jump of a Poisson process in the interval,
which is again $1-\exp(-r_{i}\Lambda(T))$.

We may also model the dose of a disease causing agent to which $i$ is exposed
as the result of a probabilistic process governed by the distance to a
hazardous environmental source. Recall that the environmental source of
disease is denoted $c$ and is located at $\ell(c)$. Suppose $W(T)$ is a second
Poisson process with mean measure $f(t)$ that defines the pathogens emitted
from this hazardous source over the period $[0,T]$; let
$F(T)=\int_{0}^{T}f(t)dt$.

$\displaystyle W(T)\mid F(T)\sim\text{Poisson}(\textstyle\int_{0}^{T}f(t)dt).$
(13)

Let $N_{i}(T)\mid W(T)$ be the dose that reaches unit $i$ over the time
interval. We assume that conditional on $W(T)$, $N_{i}(T)$ is binomial
distributed with a parameter dependent on distance to the hazard:

$\displaystyle N_{i}(T)\mid W(T)\sim\text{Binomial}(W(T),p(d_{i})).$

If individual $i$ is located at $s_{i}$, then we can define the probability
that a pathogen reaches individual $i$ as a function of $i$’s distance to the
source: $p(d_{i})=\mathcal{K}(\mathcal{\|}s_{i}-\ell(c)\|_{2}/\rho)$.
Marginalizing over $W(T)$ yields the marginal distribution for $N_{i}(T)$

$\displaystyle N_{i}(T)\mid
F(T)\sim\text{Poisson}(F(T)\mathcal{K}(\mathcal{\|}s_{i}-\ell(c)\|_{2}/\rho)).$
(14)

Using the results from equation (12) allows the derivation of the marginal
distribution of pathogens that infect individual $i$.

$\displaystyle K_{i}(T)\mid F(T)$
$\displaystyle\sim\text{Poisson}(r_{i}F(T)\mathcal{K}(\mathcal{\|}s_{i}-\ell(c)\|_{2}/\rho)).$

Let $Y_{i}(T)$ be the event that $i$ becomes infected from the hazardous
source located at $\ell(x)$ over time interval $[0,T]$. Then $Y_{i}(T)$ is
conditionally Bernoulli distributed:

$\displaystyle Y_{i}(T)\mid f$
$\displaystyle\sim\text{Bernoulli}(1-\exp(-r_{i}F(T)\mathcal{K}(\mathcal{\|}s_{i}-\ell(c)\|_{2}/\rho))).$
(15)

Equation (15) is similar to equation (2) but with a different inverse link
function: $g^{-1}(\theta)=-\log(1-\theta)$,

$\displaystyle Y_{i}(T)\mid\theta_{i}(T)$
$\displaystyle\sim\text{Bernoulli}(\theta_{i}(T))$ (16) $\displaystyle
g^{-1}(\theta_{i}(T))$
$\displaystyle=r_{i}F(T)\mathcal{K}(\mathcal{\|}s_{i}-\ell(c)\|_{2}/\rho)).$
(17)

###  5.2 Expanding the model to include background rate of exposure

The generative model can readily accommodate a term for representing
spatially-invariant background exposure if the baseline exposure is modeled as
an independent homogeneous Poisson process. Let the mean measure of the
background exposure process be $\lambda_{b}dt$ . Then the total dose is the
sum of the two exposures:

$\displaystyle
N_{i}(T)\sim\text{Poisson}(F(T)\mathcal{K}(\mathcal{\|}s_{i}-\ell(c)\|_{2}/\rho)+\lambda_{b}T).$

This then yields a model for binary infection as:

$\displaystyle P(Y_{i}(T)=1\mid f)$
$\displaystyle=1-\exp(-r_{i}(F(T)\mathcal{K}(\mathcal{\|}s_{i}-\ell(c)\|_{2}/\rho)+\lambda_{b}T)).$
(18)

The difference between the model in equation (3) and equation (18) is that the
risk from the environmental hazard adds to the background risk. In equation
(3) the odds of an observation being a disease case is modeled as:

$\displaystyle\text{odds}(P(Y_{i}(T)=1\mid f))$
$\displaystyle=r_{i}(1+\alpha\exp(-\mathcal{\|}s_{i}-c\|_{2}^{2}/\rho^{2})).$
(19)

so we see that risk from the environmental hazard multiplies the background
rate, $r_{i}$.

###  5.3 Modeling exposure to extensive environmental hazards

If the source is extensive in space, we can extend the generative model for
point source hazards. Let $\mathcal{C}$ be the set of points that make up the
environmental hazard and allow the intensity of the Poisson distributed
pathogens, $f(t)$ in equation (13), to depend on $c\in\mathcal{C}$: $f(c,t)$,
and model the pathogens present at the source as a nonhomogeneous Poisson
process with domain $\mathcal{C}\times\R^{+}$. Then the number of pathogens in
a subset $C\subset\mathcal{C}$ with differential volume element $dc$ is :

$\displaystyle W(T)\mid f\sim\text{Poisson}\left(\int_{C\times
T}f(c,t)dcdt\right).$ (20)

If we allow $f(c,t)$ to be a stochastic process itself, we can model $W$ as a
doubly-stochastic Poisson process. Then the event that unit $i$ gets infected
over time interval $[0,T]$ from section $C$ of the environmental hazard is
Bernoulli distributed:

$\displaystyle Y_{i}(T)\mid f$
$\displaystyle\sim\text{Bernoulli}\left(1-\exp\left(-r_{i}\int_{C\times
T}\mathcal{K}(\mathcal{\|}s_{i}-\ell(c)\|_{2}/\rho)f(c,t)dcdt\right)\right).$
(21)

The assumption that $W(T)$ is a nonhomogeneous Poisson process ensures that
$f(c,t)$ is integrable, and thus given the property that $0<\mathcal{K}\leq 1$
the integral in Equation 21 is well-defined.

A question remains, however, as to how to model $f(c,t)$. Our model must
guarantee almost surely integrable functions $f(c,t)$ but be flexible so as to
adapt to many different scenarios. For this reason we model $\log f(c,t)$ as a
Gaussian process. In the next two subsections we describe the properties of
Gaussian processes and doubly-stochastic Poisson processes with log-Gaussian
intensities.

####  5.3.1 Gaussian processes

Gaussian processes are stochastic processes over a domain $\mathcal{V}$ where
every finite-dimensional joint distribution of the stochastic process is
multivariate normally distributed:

$\mathbf{z}\sim\text{Multivariate
Normal}(\boldsymbol{\mu},\boldsymbol{\Sigma}),$

for any points $v_{i}$, $i\in[1,\dots,N]$, with
$\boldsymbol{\mu}_{[i]}=\mu(v_{i})$ and $\Sigma_{[i,j]}=\sigma(v_{i},v_{j})$
and each element of $\mathbf{z}$ associated to a single $v_{i}$. Gaussian
processes are completely characterized by the mean function $\mu(v)$ and
covariance kernel $\sigma(v,v^{\prime})$. The Gaussian process represents a
useful prior for unknown functions because the sample paths of a GP can be
considered random functions from $\mathcal{V}\to\R$. The properties of the
sample paths, such as differentiability and continuity, are defined by
$\mu(v)$ and $\sigma(v,v^{\prime})$. An example of a covariance kernel is the
exponentiated quadratic kernel:

$\displaystyle\sigma(v,v^{\prime}\mid\alpha,\omega)=\alpha^{2}\exp\left(-\frac{1}{2\omega^{2}}\mathcal{\|}v-v^{\prime}\|_{2}^{2}\right),$
(22)

which generates infinitely differentiable sample paths as long as
$\sigma(v,v^{\prime}\mid\alpha,\omega)$ is positive definite within
$\mathcal{V}$[30] . The hyperparameter $\omega$ controls the nonlinearity of
the function, in that large values of $\omega$ lead to almost-linear
functions, while smaller values of $\omega$ lead to functions with many peaks
and troughs over the domain. The hyperparameter $\alpha$ controls the marginal
variance of the Gaussian random variable at a given location $v$. In our case,
the domain of the Gaussian process is $(\mathcal{C}\times\R^{+})$ and
$v=(c,\tau)$. Next we describe how to use Gaussian processes to define the
log-intensity of a Poisson process.

####  5.3.2 Log-Gaussian Cox processes for integrated exposure

A log-Gaussian Cox process (LGCP) is a doubly-stochastic Poisson process where
the intensity function $\lambda(c,t)$ is an exponentiated Gaussian process:
$\lambda(c,t)=\exp Z(c,t)$ [25]. This means that conditional on a draw for
$Z(c,t)$, the LGCP is a nonhomogeneous Poisson process with intensity $\exp
Z(c,t)$.

As such, conditional on the intensity process, the LGCP inherits the complete
randomness property from homogeneous Poisson processes. Namely if
$C_{1},\dots,C_{k}\subset(\mathcal{C}\times\R^{+})$ and $C_{i}\cap
C_{j}=\emptyset$ then $W(C_{1}),\dots,W(C_{k})$ are independent Poisson random
variables with intensities $\int_{C_{i}}\lambda(c,t)dcdt$. Let the collection
of pairs $(C_{m},T_{l}),m\in[1,\dots,M],l\in[1,\dots,L]$ be a partition of
$(\mathcal{C}\times T)$. Then if $Z(c,t)$ is a GP with domain
$(\mathcal{C}\times\R^{+})$, with a valid covariance function $\sigma$ (as
defined in [25, p. 453]), $W(C_{m},T_{l})\mid
Z(c,t)\sim\text{Poisson}(\int_{C_{m}\times T_{l}}\exp Z(c,t)dcdt)$, and
$W(C_{m},T_{l})\perp\\!\\!\\!\perp W(C_{j},T_{l})\mid Z(c,t)\forall j\neq m$.

We can use the LGCP to build our model for extensive environmental hazards. If
we treat section $C_{m}$ of the source as if it were a point source, we may
define the number of pathogens from $C_{m}$ to which individual $i$ is exposed
over the time interval $T_{l}$ as:

$\displaystyle
K(C_{m},T_{l})\mid\mathbf{z}\sim\text{Poisson}(\mathcal{K}\left(\tfrac{\mathcal{\|}\ell(\bar{C}_{m})-s_{i}\|_{2}}{\rho}\right)\exp(\mathbf{z}_{[m,l]})\Delta(C_{m})\Delta(T_{l})).$

Given the independence of $W(C_{m},T_{l})\perp\\!\\!\\!\perp
W(C_{j},T_{l})\mid Z(c,t)\forall j\neq m$, we can express the total number of
particles individual $i$ is exposed to as

$\displaystyle\sum_{l=1}^{L}\sum_{m=1}^{M}K(C_{m},T_{l})\mid\mathbf{z}\sim\text{Poisson}\left(\sum_{l=1}^{L}\sum_{m=1}^{M}\mathcal{K}\left(\tfrac{\mathcal{\|}\ell(\bar{C}_{m})-s_{i}\|_{2}}{\rho}\right)\exp(\mathbf{z}_{[m,l]})\Delta(C_{m})\Delta(T_{l})\right).$

Finally, given a susceptibility parameter $r_{i}$ representing the probability
that a single pathogen infects individual $i$, the probability that individual
$i$ becomes infected over the interval $[0,T]$ is

$\displaystyle Y_{i}(T)$ $\displaystyle\mid\mathbf{z}\sim$ (23)
$\displaystyle\text{Bernoulli}\left(1-\exp\left(-r_{i}\,\sum_{l=1}^{L}\sum_{m=1}^{M}\mathcal{K}\left(\tfrac{\mathcal{\|}\ell(\bar{C}_{m})-s_{i}\|_{2}}{\rho}\right)e^{\mathbf{z}_{[m,l]}}\Delta(C_{m})\Delta(T_{l})\right)\right).$
(24)

Using the same argument as above, we can show that as $M,L\to\infty$ our
approximation converges to the integral in Equation 21, substituting $\exp
Z(c,t)$ for $f(c,t)$, under conditions on the kernel and mean function of the
Gaussian process. The proof is shown in the Supplementary Material.

The full model is then

$\displaystyle Y_{it}\mid
Z(c,t)\sim\text{Bernoulli}\left(1-\exp\left(-r_{i}\,\int_{\mathcal{C}\times[0,T]}\mathcal{K}\left(\tfrac{\mathcal{\|}\ell(c)-s_{i}\|_{2}}{\rho}\right)\exp(Z(c,t))dcdt\right)\right).$
(25)

###  5.4 Modeling susceptibility dependent on covariates

Often there are covariates $\mathbf{x}_{i}\in\R^{K}$ associated with each
individual $i$ that may predict the individual’s susceptibility to infection.
For example, these might be measurements on age, diet, or comorbidities. Then
we could use a log-linear model for $\lambda_{e}r_{i}$ conditional on
covariates $\mathbf{x}_{i}$, with $\boldsymbol{\gamma}\in\R^{K}$:

$\displaystyle\lambda_{e}r(\mathbf{x}_{i})=e^{\boldsymbol{\gamma}^{\prime}\mathbf{x}_{i}}.$

Including the background rate of exposure along with the covariates would
result in the approximate observational model:

$\displaystyle\begin{split}Y_{i}(T)\mid&\mathbf{z}\sim\\\
&\text{Bernoulli}\left(1-e^{-e^{\boldsymbol{\gamma}^{\prime}\mathbf{x}_{i}}\left(\lambda_{b}+\sum_{l=1}^{L}\sum_{m=1}^{M}\mathcal{K}\left(\tfrac{\mathcal{\|}\ell(\bar{C}_{m})-s_{i}\|_{2}}{\rho}\right)e^{\mathbf{z}_{[m,l]}}\Delta(C_{m})\Delta(T_{l})\right)}\right).\end{split}$
(26)

The exact model is thus

$\displaystyle
Y_{i}(T)\mid\mathbf{z}\sim\text{Bernoulli}\left(1-e^{-e^{\boldsymbol{\gamma}^{\prime}\mathbf{x}_{i}}\left(\lambda_{b}+\int_{\mathcal{C}\times[0,T]}\mathcal{K}\left(\tfrac{\mathcal{\|}\ell(c)-s_{i}\|_{2}}{\rho}\right)e^{Z(c,t)}dcdt\right)}\right).$
(27)

The approximate and exact models are well-defined in the sense that they admit
proper posterior measures. Furthermore, the discretizaton error in Equation 26
compared to Equation 27 induces an $\mathcal{O}(M^{-s})$ error in the
respective posterior expectations calculated under each model.

######  Theorem 5.1.

Accuracy of approximate posterior moments
Let $\Phi(u,\theta,w)$ be the exact negative log-likelihood implied by the
generative model Equation 27 and let $\Phi^{M}(u,\theta,w)$ be the approximate
negative log-likelihood corresponding to the generative model Equation 26,
using either an exponential distance kernel, $\mathcal{K}(x)=\exp(-x)$, or a
Gaussian distance kernel, $\exp(-x^{2})$. For $u\in C^{0,s}(\Omega)$, let the
two posteriors be defined in terms of their Radon-Nikodym derivatives with
respect to the Gaussian measure on $C^{0,s}(\Omega)$, $\mu_{0}(\cdot)$:

$\frac{\mathrm{d}\mu}{\mathrm{d}\mu_{0}}(u,\theta)=\frac{1}{\int_{\mathcal{X}\times\Theta}\exp(-\Phi(\upsilon,\vartheta,w))\mathrm{d}\mu_{0}(\upsilon)\mathrm{d}\pi(\vartheta)}\exp(-\Phi(u,\theta,w)),$

and

$\frac{\mathrm{d}\mu^{M}}{\mathrm{d}\mu_{0}}(u,\theta)=\frac{1}{\int_{\mathcal{X}\times\Theta}\exp(-\Phi^{M}(\upsilon,\vartheta,w))\mathrm{d}\mu_{0}(\upsilon)\mathrm{d}\pi(\vartheta)}\exp(-\Phi^{M}(u,\theta,w)).$

Then if the prior for $\rho^{-1}$ has an MGF under an exponential kernel or
the prior for $\rho^{-2}$ has an MGF when the kernel is Gaussian the posterior
moments for functions of $u,\theta$ with finite second moments under each
posterior differ by a factor which is $\mathcal{O}(M^{-s})$.

The proof of the theorem relies on novel extensions of theorems in [13, 34,
32] and is shown in the Supplementary Materials.

The result stems from a finite-sample bound on the Hellinger distance between
the posteriors. These results guide our selection of $M$ and $L$ when setting
up the approximate likelihood Equation 23

While we have presented the model in its full generality, with
$f=\exp(Z(c,t))$ nonparametrically dependent on $t$, we’ll make the
simplifying assumption going forward that the intensity of the pathogen
generation from the environmental hazard is constant in time with rate
$\lambda_{e}$. In other words, $f(c,t)=f(c)\lambda_{e}$, so
$Z(c,t)=Z(c)+\log(\lambda_{e})$, for computational tractability. This leads to
the probability model for $Y_{it}$ being expressed as

$\displaystyle Y_{it}\mid
Z(c)\sim\text{Bernoulli}\left(1-\exp\left(-e^{\boldsymbol{\gamma}^{T}\mathbf{x}_{i}}\,\left(\lambda_{b}+\lambda_{e}T\int_{\mathcal{C}}\mathcal{K}\left(\tfrac{\mathcal{\|}\ell(c)-s_{i}\|_{2}}{\rho}\right)\exp(Z(c))dc\right)\right)\right).$
(28)

####  5.4.1 Computational considerations

Given the approximate exposure term:

$\sum_{l=1}^{L}\sum_{m=1}^{M}\mathcal{K}\left(\tfrac{\mathcal{\|}\ell(\bar{C}_{m})-s_{i}\|_{2}}{\rho}\right)\exp(Z(\bar{C}_{m}),\bar{T}_{l})\Delta(C_{m})\Delta(T_{l}),$

and the true modeled exposure

$\int_{\mathcal{C}\times[0,T]}\mathcal{K}\left(\tfrac{\mathcal{\|}\ell(c)-s_{i}\|_{2}}{\rho}\right)\exp(Z(c,t))dcdt,$

we can bound the approximation error with

$\displaystyle\sum_{l=1}^{L}\sum_{m=1}^{M}\int_{C_{m}\times
T_{l}}\mathcal{K}_{\rho}(\bar{C}_{m})\left|\exp(Z(c,t))-\exp(Z(\bar{C}_{m}),\bar{T}_{l})\right|dcdt$
(29)
$\displaystyle+\sum_{l=1}^{L}\sum_{m=1}^{M}\left(\mathcal{K}\left(\tfrac{\inf_{c}\mathcal{\|}\ell(c)-s_{j}\|_{2}}{\rho}\right)-\mathcal{K}\left(\tfrac{\sup_{c}\mathcal{\|}\ell(c)-s_{j}\|_{2}}{\rho}\right)\right)\int_{C_{m}\times
T_{l}}\exp(Z(c,t))dc\,dt.$ (30)

as shown in the Supplementary Materials. Of most consequence is the term
involving differences of the distance kernel:

$\mathcal{K}\left(\tfrac{\inf_{c}\mathcal{\|}\ell(c)-s_{j}\|_{2}}{\rho}\right)-\mathcal{K}\left(\tfrac{\sup_{c}\mathcal{\|}\ell(c)-s_{j}\|_{2}}{\rho}\right).$

For units that are close to the hazard, this term is approximately
$1-\mathcal{K}\left(\tfrac{\sup_{c}\mathcal{\|}\ell(c)-s_{j}\|_{2}}{\rho}\right)$.
When $\rho$ is also close to zero this term will be near 1, leading to large
integration error. Thus, when there are many units that are near the
environmental hazard, the computational grid should be finer than it would
need to be if few units were near the hazard in order to guarantee small
approximation error.

###  5.5 Model Identifiability

The probability model in Equation 26 is not identifiable as written. We define
the GP prior for intensity to be mean zero and fix the marginal variance
parameter, or $\alpha$ in Equation 22, to be one. Neither of these
restrictions will prevent the Gaussian process from representing unknown
functions [22]. Furthermore, given that the total risk of infection is a
product of the individual hazard of infection and the sum of the instantaneous
exposure from the environmental hazard and background hazard, we cannot infer
the scale of the individual hazard of infection. We fix the intercept term of
$\boldsymbol{\gamma}$ to be $0$. The term
$e^{\boldsymbol{\gamma}^{\prime}\mathbf{x}_{i}}$ now models the relative risk
of infection for two individuals at the same location. The identified model is

$\displaystyle Y_{it}$
$\displaystyle\sim\text{Bernoulli}\left(1-\exp\left(-e^{\boldsymbol{\gamma}^{\prime}\mathbf{x}_{i}}T\left(\lambda_{b}+\sum_{m=1}^{M}\mathcal{K}\left(\tfrac{\mathcal{\|}\ell(\bar{C}_{m})-s_{i}\|_{2}}{\rho}\right)\exp(\boldsymbol{\eta}_{[m]})\Delta(C_{m})\right)\right)\right)$
(31) $\displaystyle\boldsymbol{\eta}$ $\displaystyle\sim\text{Multivariate
Normal}(0,\boldsymbol{\Sigma}_{1,\omega}).$ (32)

$\lambda_{b}$ is identified as
$\mathcal{\|}\ell(\bar{C}_{m})-s_{i}\|_{2}\to\infty$, and
$e^{\boldsymbol{\gamma}^{\prime}\mathbf{x}_{i}}$ is identified by comparing
individuals within households with different values of covariates.

Given the nature of the inverse problem, the intensity of the exposure at
distance zero is not strictly identifiable, but is instead identifiable within
the context of the Gaussian process prior [34]. This is a limitation of our
problem set up and not a limitation of the model.

##  6 Canal system simulation study

Our simulation study set up is similar to our intended application of the
model: we simulate a survey of childhood diarrheal illness in households
located near a system of wastewater canals within a region extending
$10\textrm{km}$ horizontally and $4\textrm{km}$ vertically. This region is
denoted $\mathcal{R}$. The system of canals and diarrheal illness risk is
shown in Figure 1. The left-hand plot shows the geographic location of the
canal in dashed red lines, while the flow of the wastewater is shown in solid
blue arrows. There are three canal segments: the segment $x_{1}$ runs
horizontally along the bottom edge of the region, segment $y$ runs vertically
through the middle of the region, and segment $x_{2}$ runs horizontally and
intersects $y_{1}$ at $2\frac{2}{3}\textrm{km}$. In keeping with the notation
developed in Section 5.3, the extent of the environmental hazard is the set
$\mathcal{C}=\\{x_{1},x_{2},y\\}$.

Sources of wastewater are denoted as $\upsilon$ and are indexed by the canal
segment to which they are associated; sinks are denoted $\delta$ and are
similarly indexed. The diarrhea-causing pathogens along segment $x_{1}$ are
generated according to a nonhomogeneous Poisson process with
$\Lambda_{x_{1}}(c)=0.15+c^{2}/100$, while the pathogens along segment $y$ are
generated with $\Lambda_{y}(c)=\Lambda_{x_{1}}(5)+c^{2}/16$. Canal segment
$x_{2}$ has an intensity of
$\Lambda_{x_{2}}(c)=\Lambda_{y}(8/3)-1/4+c^{2}/100$. These intensities are
such that $\Lambda_{y}(0)=\Lambda_{x_{1}}(5)$ and
$\Lambda_{y}(8/3)=\Lambda_{x_{2}}(5)$, and are respectively indicated as
$\Lambda_{x_{1}\times y}$, and $\Lambda_{x_{2}\times y}$.

We simulate two populations of household locations to investigate our method’s
sensitivity to the distribution of households. One scenario, which we term the
“uniform” scenario, all houses are uniformly distributed within $\mathcal{R}$,
while in the “clustered” scenario, the houses are distributed near the canal
system. We simulate $200{,}000$ household locations, from which we draw simple
random sample of size $J$, where $J\in\\{500,1000,2000\\}$.

For each household in the population, indexed by $j$ and with geographic
location $s_{j}$, we can define the cumulative exposure to the wastewater
pathogens from a canal segment $\nu\in\mathcal{C}$ with endpoints
$\nu_{1},\nu_{2}$ as

$\int_{\nu_{1}}^{\nu_{2}}\exp\left(-\tfrac{\mathcal{\|}\ell_{\nu}(c)-s_{j}\|_{2}}{\rho}\right)\Lambda_{\nu}(c)dc.$

Then the total exposure for household $j$ from the entirety of the canal is

$\displaystyle\mathcal{E}_{j}=\sum_{\nu\in[x_{1},x_{2},y]}\int_{\nu_{1}}^{\nu_{2}}\exp\left(-\tfrac{\mathcal{\|}\ell_{\nu}(c)-s_{j}\|_{2}}{\rho}\right)\Lambda_{\nu}(c)dc.$
(33)

Of note, we have chosen the exponential kernel,
$\exp\left(-\tfrac{\mathcal{\|}\ell_{\nu}(c)-s_{j}\|_{2}}{\rho}\right)$, for
the true measure of exposure at a given distance from a differential element
of the canal.

For the $i$-th observation within the $j$-th household, we observe
$Y_{ij}\in\\{0,1\\}$, where $Y_{ij}$ is the binary indicator for disease. We
simulate $I$ total draws per household, which takes the values $10$ and $100$
for each $J$ within each scenario. The table of simulation scenarios is shown
in the Supplementary Material.

These $Y_{ij}$ are conditionally independent Bernoulli draws:

$\displaystyle
Y_{ij}\mid\Lambda,s_{j}\sim\text{Bernoulli}\bigg{(}1-\exp\bigg{(}-\exp(X_{j}\,\gamma)\bigg{(}\lambda+\sum_{\nu\in[x_{1},x_{2},y]}\text{exposure}(\nu,s_{j},\rho)\bigg{)}\bigg{)}\bigg{)},$

with $\lambda=0.05$, $\rho=0.1$, $\gamma=-0.15$, and
$X_{j}\sim\text{Normal}(0,1)$. The integrals are numerically evaluated using
Gauss-Kronod quadrature implemented in base R language’s integrate [28].

The right-hand graph in figure 1 shows the function $P(Y_{ij}=1\mid s)$ with
$X_{j}=0$. The graph shows that the risk of disease concentrates close to the
canal system and decays as the distance to the canal increases. Figure 1 also
shows that the risk of disease is higher for a fixed $y$ coordinate and an
increasing $x$ coordinate.

Figure 1: Left: Dashed lines indicate the geographic location of the canal
segments $x_{1},x_{2},y$. Blue arrows indicate the flow of wastewater. Crosses
indicate points of interest on the canal network: $\upsilon$s are sources of
wastewater, $\delta$s are sinks of wastewater, and $p$s are canal
intersections. $\Lambda$s denote the intensity function of the canal segment
or point. Right: True probability surface, with arrows depicting the flow of
water through the canal.

###  6.1 Inferential model likelihood

The inferential model is that of Equation 31 applied to the canal system shown
in Figure 1. We define the finite dimensional realization of
$\log\Lambda_{\nu}(c)=Z_{\nu}(c)$ as $\mathbf{z}_{\nu}$, with dimension
$M_{\nu}$, for canal segment $\nu$. This finite-dimensional draw of the
Gaussian process prior is associated with partition
$\\{(C_{m},\bar{C}_{m},\Delta(C_{m}))\mid m=1,\dots,M_{\nu}\\}$ such that
$\bigcup_{m=1}^{M_{\nu}}C_{m}=\nu$. As in Section 5.3.2 the centroid of
partition section $m$ is $\bar{C}_{m}$, and, as such, the $m^{\textrm{th}}$
element of $\mathbf{z}_{\nu}$ is $Z_{\nu}(\bar{C}_{m})$. Then we can define
the approximate modeled exposure as:

$\displaystyle
F(\nu,s_{j},\rho,\mathbf{z}_{\nu})=\sum_{m=1}^{M_{\nu}}\exp\left(-\tfrac{\mathcal{\|}\ell_{\nu}(C_{m})-s_{j}\|_{2}}{\rho}\right)\exp((\mathbf{z}_{\nu})_{m})\Delta(C_{m}).$
(34)

Let $\theta^{\mathrm{environ}}_{j}$ be the total approximate exposure:

$\displaystyle\theta^{\mathrm{environ}}_{j}=\sum_{\nu\in\\{x_{1},x_{2},y\\}}F(\nu,s_{j},\rho,\mathbf{z}_{\nu}).$
(35)

We assume the functional form for the kernel, namely the exponential kernel,
is known.

The full inferential model is:

$\displaystyle Y_{ij}\mid
x_{j},\mathbf{z}_{x_{1}},\mathbf{z}_{x_{2}},\mathbf{z}_{y}\sim\text{Bernoulli}\left(1-\exp\left(-\exp(x_{j}\,\gamma)\left(\lambda+\theta^{\mathrm{environ}}_{j}\right)\right)\right).$

Given the discussion in Section 5.4.1, we would expect larger error the
clustered scenario vs. the uniform scenario. We thus use two grid sizes to
investigate the impact of approximation error on our inferences. We let
$M_{x_{1}}=M_{x_{2}}=M_{y}=M$ for $M=40,160$. Then $\mathbf{z}_{\nu}$ is in
$\R^{M}$ for each $\nu$. The partition associated with $x_{1}$ and $x_{2}$ is

$\displaystyle\\{([10\,\tfrac{n-1}{M},10\,\tfrac{n}{M}],10\,\tfrac{2n-1}{2M},10/M)\mid
n=1,\dots,M\\}.$ (36)

Thus the $n^{\mathrm{th}}$ element of $\mathbf{z}_{x_{1}}$ is
$Z_{x_{1}}(10\,\frac{2n-1}{M})$, while the $n^{\mathrm{th}}$ element of
$\mathbf{z}_{x_{2}}$ is $Z_{x_{2}}(10\,\frac{2n-1}{M})$. The partition
associated with $\mathbf{z}_{y}$ is

$\displaystyle\\{([\tfrac{8}{3}\,\tfrac{n-1}{M/2},\tfrac{8}{3}\,\tfrac{n}{M/2}],\tfrac{8}{3}\tfrac{2n-1}{M},\tfrac{16}{3M})\mid
n=1,\dots,M/2\\},$ (37)
$\displaystyle\\{([\tfrac{8}{3}+\,\tfrac{4}{3}\tfrac{n-21}{M/2},4\,\tfrac{n-20}{M/2}],\tfrac{4}{3}\tfrac{2(n-M/2)-1}{M}+\tfrac{8}{3},\tfrac{8}{3M})\mid
n=M/2+1,\dots,M\\}.$ (38)

####  6.1.1 Log-intensity priors

The Gaussian process prior we use for $\mathbf{z}_{x_{1}}$,
$\mathbf{z}_{x_{2}}$ and $\mathbf{z}_{y}$ reflects that
$\Lambda_{x_{1}}(5)=\Lambda_{y}(0)$ and
$\Lambda_{x_{2}}(5)=\Lambda_{y}(\frac{8}{3})$. We impose the constraint by
conditioning the values of $\mathbf{z}_{x}$, $\mathbf{z}_{x_{2}}$, and
$\mathbf{z}_{y}$ at the intersection points to be equal. This is akin to the
construction of string Gaussian processes introduced in [31], which explores
the formal construction of Gaussian process priors connected at intersection
points such that the Gaussian process defined at the intersection is finitely
differentiable. More details on the prior construction in the simulation study
can be found in the Supplementary Materials.

In order to formulate priors for the parameters $\lambda$, $\gamma$, $\rho$,
and $\alpha$, we sampled from the prior predictive distribution, as advised in
[18]:

$p(y)=\int
p(y\mid\boldsymbol{\theta})p(\boldsymbol{\theta})d\boldsymbol{\theta},$

if the vector $\boldsymbol{\theta}$ represents a concatenation of all of the
model parameters. The goal is to generate plausible observations from our
model with the joint prior distribution $p(\boldsymbol{\theta})$.

For $\lambda$ and $\gamma$ we use independent $\text{Normal}^{+}(0,0.3)$
priors, and for $\alpha$ we use a $\text{Gamma}(4,1)$ prior. For $\rho$, we
use a weakly-informative prior ([20]) of $\text{Normal}^{+}(0,0.5)$. Note that
this prior for $\rho$ does not yield an MGF for $\frac{1}{\rho}$ as required
by Theorem 5.1; despite this we did not observe a marked difference in model
performance compared to simuluation studies with MGFs for $\frac{1}{\rho}$.

###  6.2 Target estimands

Our model’s inferential target is the cumulative exposure from the canal,
defined above in Equation 33 for household $j$ as $\mathcal{E}_{j}$. In order
to measure how well our model predicts this exposure, we measured the bias,
and posterior credible interval coverage for this quantity Recall the model’s
approximate exposure is $\theta^{\mathrm{environ}}_{j}$, defined in eqs. 34 to
35. Let $\rho^{\star}$ be the value of $\rho$ that generated the simulated
data, in our case $0.1$. Then the absolute bias in this estimand is

$\displaystyle\left|\theta^{\text{environ}}_{j}-\mathcal{E}_{j}\right|$
$\displaystyle=\Bigg{|}\sum_{\nu\in[x_{1},x_{2},y]}\sum_{m=1}^{M_{\nu}}\int_{C_{m}}\Big{(}\exp\left(-\tfrac{\mathcal{\|}\ell_{\nu}(\bar{C}_{m})-s_{j}\|_{2}}{\rho}\right)\exp(Z(\bar{C}_{m}))$
$\displaystyle\phantom{\sum_{\nu\in[x_{1},x_{2},y]}\sum_{m=1}^{M_{\nu}}\int_{C_{m}}}-\exp\left(-\tfrac{\mathcal{\|}\ell_{\nu}(c)-s_{j}\|_{2}}{\rho^{\star}}\right)\Lambda_{\nu}(c)\Big{)}dc\Bigg{|}.$

Let
$\mathcal{K}_{\rho}(c)=\exp\left(-\tfrac{\mathcal{\|}\ell_{\nu}(c)-s_{j}\|_{2}}{\rho}\right)$
The error within a partition interval $C_{m}$ is bounded above by

$\left|\theta^{\text{environ}}_{j}-\mathcal{E}_{j}\right|\leq\sum_{\nu\in\\{x_{1},x_{2},y\\}}\int_{\nu}\left|\exp(Z_{\nu}(c))-\Lambda_{\nu}(c)\right|dc+\sup_{c\in\nu}\Lambda_{\nu}(c)\int_{\nu}\left|\mathcal{K}_{\rho}(c)-\mathcal{K}_{\rho^{\star}}(c)\right|dc.$

The upper bound on the absolute bias in the estimand is thus a function of the
integrated absolute error in the intensity approximation, and the integrated
error in our inference for $\rho$ weighted by the true intensity function, and
the resolution of the partition for $\nu$. Thus it is of interest to quantify
the approximate integrated absolute and mean-squared error in
$\Lambda_{x_{1}},\Lambda_{x_{2}},\Lambda_{y}$, as well as the bias for $\rho$.

####  6.2.1 Error estimates

The bias for a point estimator $\hat{\phi}$ with true value $\phi^{\star}$ is
calculated as
$\mathrm{bias}(\hat{\phi},\phi^{\star})=\hat{\phi}-\phi^{\star}$. Our point
estimator for each parameter is the posterior mean so in the results that
follow, the bias for a given dataset $\mathcal{D}$ is
$\mathrm{bias}(\mathbb{E}\left[\phi\mid\mathcal{D}\right],\phi^{\star})$. The
expectation over datasets is approximated by the empirical mean over $S$
simulated datasets $\\{\mathcal{D}_{s},s=1,\dots,S\\}$ is
$\frac{1}{S}\sum_{s}\mathrm{bias}(\mathbb{E}\left[\phi\mid\mathcal{D}_{s}\right],\phi^{\star})$.
Similarly, the mean-squared error is calculated as
$\frac{1}{S}\sum_{s}\mathrm{bias}(\mathbb{E}\left[\phi\mid\mathcal{D}_{s}\right],\phi^{\star})^{2}$.
We also compute the empirical coverage of the equi-tailed $80\%$-credible
intervals for the household-level environmental exposure for a posterior
quantile function for $\theta_{j}$ given dataset $\mathcal{D}$
$Q_{\theta_{j}\mid\mathcal{D}}(p)$ as
$\mathrm{cover}\left(Q_{\theta_{j}\mid\mathcal{D}},\theta_{j}^{\star}\right)=\mathbbm{1}\left(\theta^{\star}_{j}\in\left(Q_{\theta_{j}\mid\mathcal{D}}(0.1),Q_{\theta_{j}\mid\mathcal{D}}(0.9)\right)\right).$
Then the empirical mean coverage across simulations is given as
$\tfrac{1}{S}\textstyle\sum_{s}\tfrac{1}{J}\textstyle\sum_{j}\mathrm{cover}(Q_{\theta_{j}\mid\mathcal{D}_{s}},(\theta_{j})_{s}^{\star}).$

###  6.3 Inference procedure

We run full Bayesian inference inference in CmdStanR, an implementation of the
Stan modeling language and inference algorithms using dynamic Hamiltonian
Monte Carlo [10, 4, 19]. Each model was run with four Markov chain Monte Carlo
chains for $2{,}000$ iterations of warmup and $2{,}000$ iterations post-warmup
samples with a target Metropolis acceptance rate of $0.95$ during warmup.
Convergence was monitored using the Gelman-Rubin diagnostic, $\hat{R}$, [21,
36]. All parameters achieved $\hat{R}$ near 1 ($\max\hat{R}<1.01$), and the
minimum bulk and tail effective sample size divide by the total post-warmup
samples across all parameters and simulations, was $0.07$ and $0.04$. While
these figures are lower than the recommended $10\%$ minimum effective sample
size cutoff, we note that when $M=160$, only $2$ out of $1{,}200$ datasets had
minimum tail effective sample size less than $10\%$ of the total post-warmup
sample size, while only $1$ out of $1{,}200$ dataset had minimum bulk
effective sample size out of total post-warmup sample size of less than
$10\%$. For $M=40$, the number was $4$ and $1$ out of $1{,}200$ for minimum
tail and bulk effective sample sizes of less than $10\%$. There were divergent
transitions during sampling for small minority of models. The total divergent
transitions were small compared to the total post-warmup samples ($\approx
0.2\%$).

Given our reliance on MCMC sampling, our posterior mean and quantile
estimators are Monte-Carlo estimators.

###  6.4 Results

The most salient result from the simulation study is that the distribution of
households with respect to the canal has a large influence on the accuracy of
the inferences. All results that follow use the posterior mean of the
parameter as the estimator. In figure 3, we can see that we estimate $\rho$
more precisely when households are clustered near the canal, but, on the
contrary, we estimate $\lambda$, the spatially-invariant risk of disease, more
precisely when houses are uniformly distributed. This makes sense, as the
units that are most informative about $\lambda$ are those that are far from
the canal, and we have far fewer of those observations when houses are
clustered near the canal. Naturally, we have more households that are far from
the canal when the houses are arrayed uniformly on the $[0,10]\times[0,4]$
plot of land.

Figure 2: Integrated mean absolute error for $\Lambda_{x_{1}}$ and
$\Lambda_{x_{2}}$ with $\pm 2$ standard errors plotted as black bars, $10$
observations per household, grid resolution of $M=160$.

However, there is more information about the intensities $\Lambda_{\nu}(c)$
near the canal, so we see in Figure 2 that the clustered household scenario
allows for smaller integrated mean absolute error compared to the uniform
scenario.

When the model is applied to either clustered or uniformly sampled households,
the 80% intervals achieve the nominal coverage. The uniform scenario yields
negatively biased estimates of the sample average environmental exposure,
$\frac{1}{J}\sum_{j}\theta_{j}^{\mathrm{environ}}$. This is likely due to the
fact that the prior for $\rho$ shrinks the posterior towards zero so with less
information about $\rho$ in the observed data in the uniform scenario, the
prior continues to shrink the integrated risk towards zero.

Figure 3: MSE for $\rho$ and $\lambda$ with $\pm 1.96$ standard errors plotted
as black bars, $x$-jittered for clarity on the plot for $\rho$, $10$
observations per household, grid resolution of $M=160$.

Figures for $100$ observations per household are shown in the Supplementary
Material, and a comparison of results for grid resolution of $M=40$ vs.
$M=160$ is presented. The comparisons show very low rates of coverage of
$\theta_{j}$ for the $80\%$ posterior intervals when $M=40$. This coverage
gets worse as the number of households increase. The results also show that
bias in our posterior-mean estimators for $\theta_{j}$ persists despite
increasing data sizes. These results highlight the importance of using a grid
resolution that is appropriate for the problem setting. Our intuition that
approximation error would be worse for the clustered scenario is borne out in
figures shown in the Supplementary Material. The picture is complicated by the
fact that, on a percentage basis, the bias is worse for the uniform scenario,
as seen in figures contained in the Supplementary Material. This is due to the
fact that many households in the uniform scenario have very small risk from
the canal system because $\rho=100\mathrm{m}$ and the household locations are
uniformly distributed in a $10\mathrm{km}$ by $4\mathrm{km}$ square.

##  7 Application

We apply our integrated risk model to survey data collected from 2017 through
2019 measuring how household proximity to a wastewater canal influences
childhood diarrheal incidence in Mezquital Valley, Mexico. See [12] for more
detail on data collection, and descriptive statistics. The data are
longitudinal measurements of diarrheal disease in children by household. These
households are located along and near the wastewater canals, and grouped into
small localities. GPS coordinates were taken for each household, along with
the GIS data for the canals. Privacy concerns prevent us from sharing the full
map of the households.

###  7.1 Models

We model $Y_{tijk}$, survey responses of diarrheal illness for child $i$ in
household $j$ at survey wave $t$ in locality $k$. The model must account for
changes in susceptibility due to age of the child, wealth of the household,
parental education, and the intra-local correlation of exposure. We fit two
models, the first of which is the model presented in subsection 5.4, the
second of which is the model fitted in [12].

The portion of the Mezquital Valley wastewater canal system on which we are
focused has $43$ segments. We index these segments $\nu$ by $q$, of which
there are $Q$ total segments: $\nu_{q},q=1,\dots,Q$. Let the parameters
accounting for age-related differences in susceptibility be
$\beta_{\texttt{age}}$, the parameter for differences in susceptibility over
time be $\beta_{\texttt{wave}[t]}$, and the wealth and education-related
parameters be $\beta_{\texttt{wealth}},\beta_{\texttt{educ}}$, respectively.
Let $\beta_{k}$ be the increased exposure for locality $k$ with respect to
locality $1$. The parameter $\rho$ is defined as the spatial scale of exposure
to the canal, and $\lambda$ is defined as the spatially-invariant exposure to
diarrheal illness.

As in Equation 34 we define $F(\nu,s_{j},\rho,\mathbf{z}_{\nu})$ to be the
exposure at household location $s_{j}$ to canal segment $\nu$ for a given
bandwidth $\rho$:

$F(\nu,s_{j},\rho,\mathbf{z}_{\nu})=\sum_{m=1}^{M_{\nu}}\exp\left(-\tfrac{\mathcal{\|}\ell_{\nu}(\bar{C}_{m})-s_{j}\|_{2}}{\rho}\right)\exp((\mathbf{z}_{\nu})_{m})\Delta(C_{m}),$

Let $\Sigma_{\nu,\omega}$ be the marginal covariance matrix for multivariate
Gaussian random variable $\mathbf{z}_{\nu}$, and let
$\Sigma_{\nu,\omega}(\nu_{1},\nu_{2})$ be the conditional covariance matrix
conditional on the values of the Gaussian process at points $\nu_{1}$,
$\nu_{2}$. Let $\mu_{\nu,\omega}(\nu_{1},\nu_{2})$ be the conditional mean
function also dependent on the values of the random field at
$\nu_{1},\nu_{2}$. Then we may define the full inferential model as

$\displaystyle\begin{split}Y_{tijk}&\sim\text{Bernoulli}(1-\exp(-\lambda_{tijk}))\\\
\lambda_{tijk}&=\exp(\beta_{\texttt{age}[it]}+\beta_{\texttt{wave}[t]}+\beta_{\texttt{wealth}[j]}+\texttt{educ}_{j}\beta_{\texttt{educ}})\\\
&\quad\times(\exp(\beta_{\texttt{local}[k]})+\sum_{q=1}^{Q}F(\nu_{q},s_{j},\rho,\mathbf{z}_{\nu}))\\\
\mathbf{z}_{\nu}\mid
Z(\nu_{1}),Z(\nu_{2})&\sim\text{GP}(\mu_{\nu,\omega}(\nu_{1},\nu_{2}),\Sigma_{\nu,\omega}(\nu_{1},\nu_{2}))\\\
\rho\sim\text{Normal}^{+}(0,0.5),\quad&\beta_{\texttt{local}{j}}\sim\text{Normal}(0,1),\quad\omega\sim\text{Gamma}(3.7,0.9)\end{split}$
(39)

Note that this prior for $\rho$ does not yield an MGF for $\frac{1}{\rho}$ as
required by Theorem 5.1; despite this we did not observe a substantial
difference in model inferences compared to those where
$\frac{1}{\rho}\sim\text{Normal}(0,b^{2})$. Distance is measured in
kilometers, so we discretized the canal in 50-meter-long segments so
$\Delta(C_{m})=0.05\,\forall\,m$. This simplifies the Gaussian process
construction in that we do not need to scale $\omega$ in order to define
distances between canal segments with different discretization sizes. The
prior for the length-scale of the Gaussian process puts 99% of its mass
between $7\text{km}$ and $130\text{km}$, which enforces the soft constraint
that intensity slowly varies along the canal. The total length of the canal is
about $4{,}400\text{km}$, so a $130\text{km}$ length scale is still relatively
local compared to the total length of the canal.

We compare this model to a version of the minimum distance model presented in
subsection 3.2, which we refer to as the shortest-distance model.
Specifically, this model is a logistic regression with a predictor for
shortest distance to the canal, and is presented in detail in the
Supplementary Materials, as well as in [12]. The most important detail is that
the predictors include a term for the shortest distance to the canal for the
household. This term enters into the log-odds additively, among other
covariate effects, as:

$\beta_{\text{canal}}\log\left(\min_{c\in\mathcal{C}}\mathcal{\|}s_{j}-\ell(c)\|_{2}\right).$

The model will help elucidate the differences between our new method and the
simpler methods currently in use. We will call this model the shortest-
distance model, while we refer to our proposed model as the integrated
exposure model.

###  7.2 Model inferences

The integrated exposure model infers that there is a small increased risk of
diarrheal infection as distance to a point on the canal decreases. The
posterior mean of $\rho$ is $0.01$ with a standard deviation of $0.006$. We
estimate the posterior mean of $\lambda$ to be $0.016$ with a standard
deviation of $0.005$.

This results in nearly zero exposure to wastewater at distances of greater
than 200 meters. On the other hand, the shortest distance model shows a much
slower decline in risk in the right-hand panel of Figure 4. For instance, we
show the change in odds of diarrheal illness as distance to the canal
increases from ten meters to one kilometer compared to the odds of diarrheal
illness at ten meters in figure 4. The odds of diarrhea for the integrated
exposure model, given a household with location $s_{j}$, is given by

$\exp(-(\lambda+\sum_{q=1}^{Q}F(\nu_{q},s_{j},\rho,\mathbf{z}_{\nu_{q}})))-1$

so the change in odds for a household located at $s_{1}$ compared to $s_{2}$
is

$\frac{\exp(-(\lambda+\sum_{q=1}^{Q}F(\nu_{q},s_{2},\rho,\mathbf{z}_{\nu_{q}})))-1}{\exp(-(\lambda+\sum_{q=1}^{Q}F(\nu_{q},s_{1},\rho,\mathbf{z}_{\nu_{q}})))-1}-1.$

The change in odds for the shortest-distance model is given by

$\left(\frac{d_{1}}{d_{2}}\right)^{\beta_{\texttt{canal}}}-1$

where $d_{1}$ and $d_{2}$ are the shortest distances to the canal for
locations $s_{1}$ and $s_{2}$. It is clear from figure 4 that the reduction in
odds as distance to the canal increases is more extreme for the integrated
model than for the shortest-distance model. The difference between the left-
hand and right-hand panels of Figure 4 is likely due to the fact that the
change in odds under the shortest distance model is assumed to follow a power
law while there is no functional form assumed for the integrated exposure
model. This demonstrates the benefits of using a semiparametric Bayesian model
in this circumstance.

Figure 4: Posterior realizations of change in odds of diarrhea versus change
in distance to the canal compared to 10 meters. Odds for the integrated model
show the change in odds for a single household that moves laterally from the
westernmost edge of the canal. Red lines indicate posterior means.

Our integrated exposure model decocts the two processes that contribute to
exposure from environmental hazards: the geometry of the hazard with respect
to the at-risk population, and the variable concentration of enteric pathogens
along the hazard.

##  8 Discussion

We have presented a new model for exposure to environmental sources of disease
that are extensive in space. Our model is unique in that it models the total
exposure as the result of two separate contributions: the orientation of the
units at risk with respect to the environmental hazard and the distribution of
antigens within and along the hazard. This allows us to generate more detailed
policy recommendations compared to models that rely on only the shortest
distance to the hazard as a proxy for exposure. We showed how our was
connected to existing methods for inferring exposure to point-source hazards.
This connection elucidates the connection between the generative model for
disease given exposure to an environmental hazard and the discriminative
models of environmental exposure. We believe that a benefit to our method is
that extensions to the model are more readily apparent than the discriminative
approach to modeling. For example, the exposure at a certain distance to the
hazard could be modified to account for patterns of movement; this suggests
that we could incorporate human mobility data into our model. This would enter
into the model through Equation 14.

We showed in the Supplementary Material that the model satisfies basic
requirements for Bayesian models, namely that the posterior over the intensity
process is well-defined, and, more importantly, that the posterior under the
approximate model is close the posterior under the exact (though intractable)
model. This ensures that the posterior moments calculated under the
approximate model Equation 31 are close to the posterior moments one would
obtain from the exact model Equation 25.

Through extensive simulation studies, we verified that the model was able to
recover known data generating parameters under several sampling scenarios.
These scenarios were chosen to mimic idealized data collection and a more
realistic data collection scheme. We showed that the realistic data collection
scheme had several benefits over the idealized data collection scheme, namely
better inferences, in the mean-squared-error sense, for the parameters of
interest related to exposure from the environmental hazard. Our simulated data
scenarios also allowed us to test the impact of different partition sizes to
investigate how the fidelity of our integral approximation impacted
inferences. We showed that computational grids that allowed for tractable
model fitting still yielded inferences that had close to nominal interval
coverage, and declining mean squared error as the number of observed spatial
units increased.

We also used our model to infer how environmental exposure to wastewater
canals impacted the household risk of childhood diarrhea in Mezquital Valley,
Mexico. This dataset served as the inspiration to develop our model based on
our earlier investigations in [12]. We compared the new model’s inferences to
the inferences from [12], where we were able to show the benefits of our
model.

###  8.1 Limitations

Our model has several limitations. The first limitation is the need to compute
an approximate integral. We opted to use a simple integration scheme, but in
higher dimensions, this integration scheme may converge slowly.

Our implementation of the model in Stan relies using the covariance function
in the density for the Gaussian process. This is a computationally inefficient
method for evaluating the density, and faster approximate schemes are
available [33]. Thus, as it stands, given both the integration scheme and the
density evaluation, our method may not scale well as the dimension of the
environmental hazard increases, or if we include time in the Gaussian process.

The implementation of the model also uses the exponentiated quadratic
covariance function, also known as the squared exponential covariance
function, to parameterize the Gaussian process. This covariance function
results in analytic sample functions from the GP prior, which may be
unrealistically smooth for the underlying physical process that we are
modeling. The downside of this covariance function is that the posterior
concentrates on analytic functions, which means that if the true underlying
function is not analytic it will not be in the support of the posterior. The
proofs in the Supplementary Material do not rely on these smoothness
assumptions, so using a Mátern covariance function would allow the model to
capture a broader space of functions. We did not implement GPs with the Mátern
covariance kernel because Stan has not implemented a Mátern covariance
function with gradients, though it is feasible that this will change in future
versions of Stan.

Another limitation is that it is not clear whether $\rho$ and $\exp(Z(c))$ are
separately identifiable. The results of the simulation study suggest the
answer may be no. More research is needed to understand if the model is
identified, and, if not, how to do so.

The model relies on a parametric model for infection given exposure, namely
the single-hit Poisson model [8]. To the extent that this model does not
describe the relationship between exposure and infection, our method may
suffer bias in the inferences for health risk imposed by the environmental
hazard.

###  8.2 Extensions and future work

Incorporating more realistic models for exposure at the observational-unit
level, and potentially including human mobility data, is an important
direction for future work. With more data on mobility patterns, one could
include a more realistic model for the interactions between at-risk
populations and the environmental hazard.

In this vein, our model currently treats the households as fixed with respect
to the hazard with risk emanating from the hazard and received at the
household. One could instead treat the household inhabitants as interacting
with the environmental hazard according to a radially symmetric Poisson line
process and treat the pathogen exposure as occurring when the line intersects
the environmental hazard.

Our kernel specification depends only on distance to the segment or area of
the environmental hazard, but the point-source literature has investigated the
use of kernels that take into account direction as well as distance. This
could be useful in applications where the built environment can provide
further barriers to exposure. For example, in the Mezquital dataset, some
canal segments are only reachable via fields whereas other segments abut local
roads lined with houses.

Another immediate extension, as presented in section 5.3, is to allow the
concentration of disease-causing agents to be time-varying and to be modeled
using the kernel of the Gaussian process. This extension can be directly
applied to the Mezquital data example given that there are observations
recorded over distinct time points.

In applications where environmental monitoring of health hazards is feasible,
such as in air quality monitoring near roadways, we can augment our models
with direct observations of concentrations of hazardous material at the
source. This is another extension that could be applied to the Mezquital data
because wastewater samples were take at several sites along the wastewater
canal.

[Acknowledgments] The authors would like to thank Horacio Riojas-Rodríguez,
Eunice Felix-Arellano, and Sandra Rodríguez-Dozal for leading the fieldwork
responsible for collecting the Mezquital Valley dataset and Christina Siebe
for supporting data collection. The authors would also like to thank Dan
Simpson for helpful conversations about Bayesian inverse problems, and Tate
Jacobson for helpful conversations about bounds.

Y. Chen is supported by NSF DMS 2113397, NSF PHY 2027555, NASA
22-SWXC22_2-0005 and NASA 22-SWXC22_2-0015.

Supplement with accompanying proofs related to posterior convergence and
approximation error, and further results for simulation studies.

## References

  * Andrinopoulou et al. [2017] Andrinopoulou, E.-R., Rizopoulos, D., Takkenberg, J. J., and Lesaffre, E. (2017).  “Combined Dynamic Predictions Using Joint Models of Two Longitudinal Outcomes and Competing Risk Data.”  Stat Methods Med Res, 26(4): 1787–1801.
URL https://doi.org/10.1177/0962280215588340

  * Bender [2009] Bender, R. (2009).  “Introduction to the Use of Regression Models in Epidemiology.”  In Verma, M. (ed.), Cancer Epidemiology, Methods in Molecular Biology, 179–195. Totowa, NJ: Humana Press.
URL https://doi.org/10.1007/978-1-59745-416-2_9

  * Berrocal et al. [2011] Berrocal, V. J., Gelfand, A. E., Holland, D. M., Burke, J., and Miranda, M. L. (2011).  “On the Use of a PM2.5 Exposure Simulator to Explain Birthweight.”  Environmetrics, 22(4): 553–571.
URL https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3116241/

  * Betancourt [2018] Betancourt, M. (2018).  “A Conceptual Introduction to Hamiltonian Monte Carlo.”  arXiv:1701.02434 [stat].
  * Biggeri and Lagazio [1999] Biggeri, A. and Lagazio, C. (1999).  “Case-control analysis around putative sources.”  Disease Mapping and Risk Assessment for Public Health. Lawson, Bertollini, Biggeri, Böhning, Lesaffre and Viel (eds). Wiley, London, UK, 271–286.
  * Brilleman et al. [2018] Brilleman, S., Crowther, M., Moreno-Betancur, M., Buros Novik, J., and Wolfe, R. (2018).  “Joint longitudinal and time-to-event models via Stan.”  StanCon 2018. 10-12 Jan 2018. Pacific Grove, CA, USA.
URL https://github.com/stan-dev/stancon_talks/

  * Brouwer et al. [2017a] Brouwer, A. F., Eisenberg, M. C., Remais, J. V., Collender, P. A., Meza, R., and Eisenberg, J. N. S. (2017a).  “Modeling Biphasic Environmental Decay of Pathogens and Implications for Risk Analysis.”  Environ Sci Technol, 51(4): 2186–2196.
  * Brouwer et al. [2017b] Brouwer, A. F., Weir, M. H., Eisenberg, M. C., Meza, R., and Eisenberg, J. N. S. (2017b).  “Dose-Response Relationships for Environmentally Mediated Infectious Disease Transmission Models.”  PLOS Computational Biology, 13(4): e1005481.
  * Calculli et al. [2010] Calculli, C., Pollice, A., and Bisceglia, L. (2010).  “Spatial Variation of Multiple Diseases in Relation to an Environmental Risk Source.”  38 WPG2010.
URL https://aisberg.unibg.it/handle/10446/950#.XbtLTEFKg5k

  * Carpenter et al. [2017] Carpenter, B., Gelman, A., Hoffman, M. D., Lee, D., Goodrich, B., Betancourt, M., Brubaker, M., Guo, J., Li, P., and Riddell, A. (2017).  “Stan : A Probabilistic Programming Language.”  Journal of Statistical Software, 76(1).
URL https://www.osti.gov/pages/biblio/1430202-stan-probabilistic-programming-
language

  * Cassell et al. [2018] Cassell, K., Gacek, P., Warren, J. L., Raymond, P. A., Cartter, M., and Weinberger, D. M. (2018).  “Association Between Sporadic Legionellosis and River Systems in Connecticut.”  The Journal of Infectious Diseases, 217(2): 179–187.
  * Contreras et al. [2020] Contreras, J. D., Trangucci, R., Felix-Arellano, E. E., Rodríguez-Dozal, S., Siebe, C., Riojas-Rodríguez, H., Meza, R., Zelner, J., and Eisenberg, J. N. (2020).  “Modeling spatial risk of diarrheal disease associated with household proximity to untreated wastewater used for irrigation in the Mezquital Valley, Mexico.”  Environmental Health Perspectives, 128(7): 077002.
  * Cotter et al. [2010] Cotter, S. L., Dashti, M., and Stuart, A. M. (2010).  “Approximation of Bayesian inverse problems for PDEs.”  SIAM journal on numerical analysis, 48(1): 322–345.
  * Crawford et al. [2019] Crawford, F. W., Marx, F. M., Zelner, J., and Cohen, T. (2019).  “Transmission modeling with regression adjustment for analyzing household-based studies of infectious disease: application to tuberculosis.”  Epidemiology (Cambridge, Mass.).
  * Diggle et al. [1997] Diggle, P., Morris, S., Elliott, P., and Shaddick, G. (1997).  “Regression Modelling of Disease Risk in Relation to Point Sources.”  Journal of the Royal Statistical Society: Series A (Statistics in Society), 160(3): 491–505.
URL
https://rss.onlinelibrary.wiley.com/doi/abs/10.1111/j.1467-985X.1997.00076.x

  * Diggle [1990] Diggle, P. J. (1990).  “A Point Process Modelling Approach to Raised Incidence of a Rare Phenomenon in the Vicinity of a Prespecified Point.”  Journal of the Royal Statistical Society. Series A (Statistics in Society), 153(3): 349–362.
URL http://www.jstor.org/stable/2982977

  * Diggle and Rowlingson [1994] Diggle, P. J. and Rowlingson, B. S. (1994).  “A Conditional Approach to Point Process Modelling of Elevated Risk.”  Journal of the Royal Statistical Society. Series A (Statistics in Society), 157(3): 433–440.
  * Gabry et al. [2019] Gabry, J., Simpson, D., Vehtari, A., Betancourt, M., and Gelman, A. (2019).  “Visualization in Bayesian Workflow.”  Journal of the Royal Statistical Society: Series A (Statistics in Society), 182(2): 389–402.
URL http://rss.onlinelibrary.wiley.com/doi/abs/10.1111/rssa.12378

  * Gabry and Češnovar [2021] Gabry, J. and Češnovar, R. (2021).  cmdstanr: R Interface to ’CmdStan’.  Https://mc-stan.org/cmdstanr, https://discourse.mc-stan.org.
  * Gelman et al. [2008] Gelman, A., Jakulin, A., Pittau, M. G., and Su, Y.-S. (2008).  “A Weakly Informative Default Prior Distribution for Logistic and Other Regression Models.”  Ann. Appl. Stat., 2(4): 1360–1383.
URL https://projecteuclid.org/euclid.aoas/1231424214

  * Gelman and Rubin [1992] Gelman, A. and Rubin, D. B. (1992).  “Inference from Iterative Simulation Using Multiple Sequences.”  Statist. Sci., 7(4): 457–472.
URL https://projecteuclid.org/euclid.ss/1177011136

  * Ghosal and Roy [2006] Ghosal, S. and Roy, A. (2006).  “Posterior consistency of Gaussian process prior for nonparametric binary regression.”  The Annals of Statistics, 34(5).
URL https://projecteuclid.org/journals/annals-of-
statistics/volume-34/issue-5/Posterior-consistency-of-Gaussian-process-prior-
for-nonparametric-binary-regression/10.1214/009053606000000795.full

  * Hickey et al. [2016] Hickey, G. L., Philipson, P., Jorgensen, A., and Kolamunnage-Dona, R. (2016).  “Joint Modelling of Time-to-Event and Multivariate Longitudinal Outcomes: Recent Developments and Issues.”  BMC Medical Research Methodology, 16(1): 117.
URL https://doi.org/10.1186/s12874-016-0212-5

  * Jalava et al. [2014] Jalava, K., Rintala, H., Ollgren, J., Maunula, L., Gomez-Alvarez, V., Revez, J., Palander, M., Antikainen, J., Kauppinen, A., Räsänen, P., Siponen, S., Nyholm, O., Kyyhkynen, A., Hakkarainen, S., Merentie, J., Pärnänen, M., Loginov, R., Ryu, H., Kuusi, M., Siitonen, A., Miettinen, I., Domingo, J. W. S., Hänninen, M.-L., and Pitkänen, T. (2014).  “Novel Microbiological and Spatial Statistical Methods to Improve Strength of Epidemiological Evidence in a Community-Wide Waterborne Outbreak.”  PLOS ONE, 9(8): e104713.
  * Møller et al. [1998] Møller, J., Syversveen, A. R., and Waagepetersen, R. P. (1998).  “Log Gaussian Cox Processes.”  Scandinavian Journal of Statistics, 25(3): 451–482.
  * Perez et al. [2012] Perez, L., Lurmann, F., Wilson, J., Pastor, M., Brandt, S. J., Künzli, N., and McConnell, R. (2012).  “Near-Roadway Pollution and Childhood Asthma: Implications for Developing “Win–Win” Compact Urban Development and Clean Vehicle Strategies.”  Environ Health Perspect, 120(11): 1619–1626.
URL https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3556611/

  * Peterson and Sanchez [2018] Peterson, A. and Sanchez, B. (2018).  “rstap: An R Package for Spatial Temporal Aggregated Predictor Models.”  arXiv:1812.10208 [stat].
URL http://arxiv.org/abs/1812.10208

  * R Core Team [2021] R Core Team (2021).  R: A Language and Environment for Statistical Computing.  R Foundation for Statistical Computing, Vienna, Austria.
URL https://www.R-project.org/

  * Ramis et al. [2011] Ramis, R., Diggle, P., Cambra, K., and López-Abente, G. (2011).  “Prostate Cancer and Industrial Pollution: Risk around Putative Focus in a Multi-Source Scenario.”  Environment International, 37(3): 577–585.
URL http://www.sciencedirect.com/science/article/pii/S0160412010002461

  * Rasmussen and Williams [2006] Rasmussen, C. E. and Williams, C. K. I. (2006).  Gaussian processes for machine learning.  The MIT Press.
  * Samo and Roberts [2015] Samo, Y.-L. K. and Roberts, S. (2015).  “String and Membrane Gaussian Processes.”  arXiv:1507.06977 [stat].
URL http://arxiv.org/abs/1507.06977

  * Simpson et al. [2016] Simpson, D., Illian, J. B., Lindgren, F., Sørbye, S. H., and Rue, H. (2016).  “Going off Grid: Computationally Efficient Inference for Log-Gaussian Cox Processes.”  Biometrika, 103(1).
  * Simpson et al. [2012] Simpson, D., Lindgren, F., and Rue, H. (2012).  “In order to make spatial statistics computationally feasible, we need to forget about the covariance function: SPDES, GMRFS, AND KERNEL METHODS.”  Environmetrics, 23(1): 65–74.
  * Stuart [2010] Stuart, A. M. (2010).  “Inverse problems: A Bayesian perspective.”  Acta Numerica, 19: 451–559.
URL
https://www.cambridge.org/core/product/identifier/S0962492910000061/type/journal_article

  * Thompson et al. [2015] Thompson, C. N., Zelner, J. L., Nhu, T. D. H., Phan, M. V., Hoang Le, P., Nguyen Thanh, H., Vu Thuy, D., Minh Nguyen, N., Ha Manh, T., Van Hoang Minh, T., Lu Lan, V., Nguyen Van Vinh, C., Tran Tinh, H., von Clemm, E., Storch, H., Thwaites, G., Grenfell, B. T., and Baker, S. (2015).  “The Impact of Environmental and Climatic Variation on the Spatiotemporal Trends of Hospitalized Pediatric Diarrhea in Ho Chi Minh City, Vietnam.”  Health Place, 35: 147–154.
URL https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4664115/

  * Vehtari et al. [2020] Vehtari, A., Gelman, A., Simpson, D., Carpenter, B., and Bürkner, P.-C. (2020).  “Rank-normalization, folding, and localization: An improved Rhat for assessing convergence of MCMC.”  Bayesian Analysis.
  * Warren et al. [2018] Warren, J. L., Grandjean, L., Moore, D. A. J., Lithgow, A., Coronel, J., Sheen, P., Zelner, J. L., Andrews, J. R., and Cohen, T. (2018).  “Investigating Spillover of Multidrug-Resistant Tuberculosis from a Prison: A Spatial and Molecular Epidemiological Analysis.”  BMC Medicine, 16(1): 122.
URL https://doi.org/10.1186/s12916-018-1111-x