# Bayesian Design with Sampling Windows for Complex Spatial Processes

K. Buchhorn1,2 K. Mengersen1,2 E. Santos-Fernández1,2 E.E. Peterson1,4 J.M.
McGree1,2,3 1School of Mathematical Sciences. Queensland University of
Technology.  2Centre for Data Science (CDS). Queensland University of
Technology.  3Australian Institute of Marine Science (AIMS). Crawley, Western
Australia.  4EP Consulting. Brisbane, Queensland.

###### Abstract

Optimal design facilitates intelligent data collection. In this paper, we
introduce a fully Bayesian design approach for spatial processes with complex
covariance structures, like those typically exhibited in natural ecosystems.
Coordinate Exchange algorithms are commonly used to find optimal design
points. However, collecting data at specific points is often infeasible in
practice. Currently, there is no provision to allow for flexibility in the
choice of design. We also propose an approach to find Bayesian sampling
windows, rather than points, via Gaussian process emulation to identify
regions of high design efficiency across a multi-dimensional space. These
developments are motivated by two ecological case studies: monitoring water
temperature in a river network system in the northwestern United States and
monitoring submerged coral reefs off the north-west coast of Australia.

##  1 Introduction

Sampling designs that optimise one or more specified utility functions are
valuable in many applied contexts, since they facilitate the intelligent
collection of data. That is, an optimal design results in greater information
efficiency, reduced sampling cost and improved estimation. However, such
schemes remain challenging to develop for spatial processes with complex
covariance structures, such as those found in natural ecosystems. Moreover,
sampling at the exact design point can be difficult in practice within these
systems due to practical challenges associated with the collection of samples
and accessibility of the sites. In this case, a more flexible sampling scheme
that identifies windows with near-optimal design efficiency can be highly
useful.

Data collected by environmental monitoring programs are essential for
understanding patterns, trends and vulnerabilities in ecological and
environmental systems. Here, we focus on river network and coral reef systems,
and develop optimal sampling regions given non-standard spatial relationships.
The covariance across locations on a branching river is particularly complex,
since the network topology is embedded in a 2D terrestrial landscape, with
directional water flow, and disparities between flow volume. Furthermore,
monitoring programs are notoriously expensive in terms of monetary, human and
technological resources (Lindenmayer and Likens, 2010; Roelfsema and Phinn,
2010). Balancing short- and long-term requirements, as well as practical
constraints of data collection can add further complication to the design
process (Kang et al., 2016). Our aim is to develop a more flexible design
approach within a Bayesian framework while achieving two goals: 1) the
determination of a set of optimal sampling locations based on a discrete set
of available locations, and 2) the formation of sampling windows within a
specified region where difficulties due to access or otherwise may be
encountered.

Experimental design problems are commonly viewed as optimisation problems. In
the classical framework, optimal experimental designs are often derived using
utilities based on the expected Fisher information matrix (e.g. Atkinson and
Donev, 1992). Such designs have been shown to depend on assumed values for
model parameters, so approaches have been proposed to remove some of this
dependence. For example, Pronzato and Walter (1985) optimised a design
criterion over a prior probability distribution for the model parameters; the
so-called ‘pseudo-Bayesian’ design approach. While pseudo-Bayesian designs are
more robust for the choice of parameter values, the prior information is not
utilised in subsequent analysis (i.e. computation of the utility). A fully
Bayesian framework provides a unified approach for incorporating prior
information and rigorously handling uncertainty about the model parameters and
the statistical model when forming a design and analysing the resulting data
(see Chaloner and Verdinelli (1995); Ryan et al. (2016)).

Various approaches to finding sampling intervals, so-called windows, are
presented in the statistical literature relating to pharmacokinetics (e.g.
Bogacka et al. (2008); Foo et al. (2012); McGree et al. (2012)). These
intervals provide flexibility in the timing of sample collection, while
assuring a high level of information efficiency. An approach for stratified
random sampling in a geostatistical simulation study by Lin et al. (2011) uses
Latin Hypercube sampling to find the optimal sub-sampling of pre-specified
regions (for a given spatial stratification). However, to our knowledge no
approaches have been developed to form sampling windows for spatial processes
in the statistical literature.

We present a Bayesian approach, using a Gaussian process emulator to
approximate the utility surface in high-dimensions. Our method builds on the
approach to stochastic optimisation employed in the Approximate Coordinate
Exchange (ACE) (Overstall and Woods, 2017) algorithm. We also demonstrate how
to derive design efficiency contours using sampling windows, which can then be
used to guide highly efficient and flexible sampling of complex spatial
processes. The new approach is applied to substantive real-world case studies,
described below.

###  1.1 Motivational Case Studies

Example 1: Water temperature is considered a “master” variable in river
ecosystems because it affects the metabolic rates of aquatic organisms, life
history events related to reproduction, hatching and maturation, and restricts
the distribution and abundance of ectothermic species (Isaak et al., 2017).
Water temperature is often monitored to understand the impacts of climate
change and land management on thermal habitats (Isaak et al., 2016; Jackson et
al., 2018) and is also used as the basis for regulatory actions (Todd et al.,
2008). The recent proliferation of inexpensive _in-situ_ sensors has
dramatically increased the potential for monitoring and data-enabled
management decisions (Isaak et al., 2017). Nevertheless, comprehensive
monitoring on river networks using in-situ sensors remains costly and complex
in terms of balancing short- and long-term priorities. Optimal sensor
placement must be assessed when sensors are initially deployed, new sensors
are added to an established monitoring program or the number of in-situ
sensors must be reduced. In addition, field conditions can be unpredictable;
some sites are inaccessible due to steep terrain, a lack of network reception
for data streaming, or private landholders refusing access. This motivates two
considerations: the design of a global optimal sampling scheme across the
river network, and a more targeted approach to identify finer-scale sampling
regions that provide near-optimal design efficiency.

In river networks, spatial relationships are characterised in part by their
topology (e.g. branching network structure, connectivity). There have been
numerous advances in statistical modelling used to describe spatial covariance
in data collected on branching river networks based on distance travelled
along the river, directionality in water flow and differences in flow volume
(Ver Hoef et al., 2006; Ver Hoef and Peterson, 2010; Peterson et al., 2013;
Santos-Fernandez et al., 2021). The ability to model this spatial component
and more accurately predict conditions throughout a river network often
provides a better understanding of habitat suitability (Isaak et al., 2016)
and leads to more effective management decisions (Sharma et al., 2021a, b).
Frequentist sampling methods for river networks have been developed (Som et
al., 2014; Wang et al., 2020), as well as pseudo-Bayesian approaches (Falk et
al., 2014; Pearse et al., 2020). However, a fully Bayesian design framework
for river networks is currently unexplored.

Example 2: Coral reefs are inextricably linked by the symbiotic relationships
formed between corals and the organisms around them (Richmond, 1993). In
particular, deep reefs have the capacity to act as refugia and to assist in
the recovery of more vulnerable shallow water reefs following a disturbance.
However, major challenges remain in understanding the extent of and managing
the vast marine estate for Australia, as well as for many other maritime
nations. Insights into remote reefs are beginning to emerge due to data
captured using new technologies such as remote and autonomous underwater
vehicles. With limited resources to invest in research and monitoring, there
is a need to optimise in-field activities to collect data that will most
efficiently achieve research and management goals.

Hard coral cover is an indicator commonly used to infer the health and
condition of a coral reef (AIMS, 2021). In this example, we consider deep
reefs, or mesophotic coral ecosystems, with data collected between 12 to 50
metres in depth. The extent to which important ecological processes change
along depth gradients in a reef is not well understood and is a critical
knowledge gap for developing future reef policies and management practices
(Kang et al., 2016). Monitoring coral reef environments often relies on a
series of images taken underwater, along a transect line. As such, collecting
data from submerged shoals can be costly given the specialised skills of those
undertaking the monitoring and the cost of equipment and running a ship.
Motivating our second design scenario, we consider the optimal placement of
transect lines across a submerged shoal taking into account the spatial nature
of the process. In practice, underwater sampling along these lines can be
challenging and unpredictable due to water currents. Thus, a more practical
solution would be to define sampling regions, as proposed here, offering more
flexibility and assurance in sampling.

The following section provides details about the general statistical modelling
used in both case studies, then outlines specific details in each case. In
Section 3, the Bayesian optimal design approach is described. Our proposed
approach to finding sampling windows is given in Section 4, along with the
concept of design efficiency contours. Finally, Section 5 illustrates the
results of these two examples and demonstrates the scope of applicability of
Bayesian design and sampling windows for spatial processes.

##  2 Model

###  2.1 Bayesian Spatial Generalised Linear Models

We start with a general description of spatial linear models and then extend
this to models for data collected on river networks and coral reefs. Let
$\textbf{d}=(d^{1},...,d^{\gamma})$ denote a design, that is, a set of
locations in the physical or parameter space. Consider a linear mixed model
with $n\times 1$ response Y, and the $n\times p$ design matrix X of
explanatory variables spatially indexed at locations, d,

$\textbf{Y}=\textbf{X}\boldsymbol{\beta}+\textbf{z}+\boldsymbol{\epsilon},$
(1)

with spatially correlated random effects
$\textbf{z}\sim\mathcal{N}(\bf{0},\boldsymbol{\Sigma_{z}})$ and independent
residuals $\boldsymbol{\epsilon}\stackrel{{\scriptstyle
iid}}{{\sim}}\mathcal{N}(0,\sigma_{0}^{2})$ yielding a covariance of
$\boldsymbol{\Sigma}=\mbox{Cov}[\textbf{Y}]=\boldsymbol{\Sigma}_{z}+\sigma_{0}^{2}\textbf{I}$
where I is the $n\times n$ identity matrix and $\sigma_{0}^{2}$ is the nugget
effect. Typically, we express $\boldsymbol{\Sigma_{z}}$ as a function based on
distance $h$ between locations d, and covariance parameters,
$\boldsymbol{\theta}_{z}$. We define the set of all parameters as
$\boldsymbol{\theta}=(\boldsymbol{\beta},\boldsymbol{\theta}_{z})^{\prime}$.

There are many covariance models that describe how the dependence between
observed data at two spatial locations decays with distance. A common example
is the exponential function,

$C(h|\boldsymbol{\theta}_{z})=\sigma_{1}^{2}\mbox{exp}\Big{(}\frac{-3h}{\alpha}\Big{)},$
(2)

where the partial sill, $\sigma_{1}^{2}$, represents the magnitude of variance
and the range parameter, $\alpha$, describes how fast the covariance decays
with an appropriate measure of distance between the locations, $h$.

The generalised linear modelling (GLM) framework provides a unified approach
to analyse data for which some function of the mean response (link function)
may vary linearly with the predictors rather than the mean response itself. We
denote the linear predictor by
$\textbf{X}\boldsymbol{\beta}=\boldsymbol{\eta}=\beta_{0}+\sum\limits^{p}_{j=1}\textbf{X}_{j}\beta_{j}$.
A link function $g(\cdot)$ relates the linear predictor to the mean of the
outcome variable $g(\mathbb{E}(\textbf{Y}|\textbf{X}))=\boldsymbol{\eta}$. In
a GLM, Y is assumed to be generated from a distribution in the exponential
family, where the probability density (or mass) function is written as,

$p(y_{i})=\mbox{exp}\Bigg{(}\frac{y_{i}\psi_{i}-b(\psi_{i})}{a_{i}(\phi)}+c(y_{i},\phi)\Bigg{)},$

where for observation $i$, $\psi_{i}$ is the location parameter, $\phi$ is the
scale parameter and $a_{i}(\cdot)$, $b(\cdot)$ and $c(\cdot,\cdot)$ are known
functions.

_Example 1: Model for River Networks_

Consider data collected on a river network exhibiting complex patterns of
spatial dependence due to multi-scale processes occurring within and between
the aquatic and terrestrial environments (Peterson et al., 2013). Two points
$d^{k}$ and $d^{s}$ on the river network are said to be _flow-connected_ if
water flows from an upstream site to a downstream site. They are _flow-
unconnected_ if they reside on the same river network and share a common
junction downstream, but do not share flow. Following Ver Hoef and Peterson
(2010), we estimate autocorrelation between points $d^{k}$ and $d^{s}$ by
constructing random variables as the integration of a moving-average function
over a white-noise random process. The tails of these moving-average functions
are unilateral and spatial autocorrelation only occurs when the tails of the
functions overlap. If the tail points upstream, it is referred to as a _tail-
up_ model and by construction, these models restrict correlation to flow-
connected locations. Stream networks are dendritic and so weights are used to
proportionally allocate (i.e. split) the tail-up moving average function at
river junctions. Tail-up models are particularly useful for modelling
organisms or materials that flow passively downstream (e.g. water pollutants).
Conversely, when the tail of the function points downstream it is referred to
as a _tail-down_ model. These models allow spatial correlation between both
flow-connected and flow unconnected locations and may be more suitable for
modelling organisms such as fish that actively move both up and downstream
(Peterson and Fausch, 2003).

A mixed modelling approach is often used for modelling river network data
since it allows for more complex patterns of spatial dependence to be
described within a single model (Peterson and Hoef, 2010). Variance components
are used to expand the spatially dependent random variable z into several
zero-mean random variables,

$\textbf{Y}=\textbf{X}\boldsymbol{\beta}+\textbf{z}_{{}_{EUC}}+\textbf{z}_{{}_{TU}}+\textbf{z}_{{}_{TD}}+\boldsymbol{\epsilon},$
(3)

with correlation structure described using Euclidean (EUC), tail-up (TU) and
tail-down (TD) covariance models. As the notation suggests, the general
covariance matrix can be formulated by a combination of such covariance
models,

$\displaystyle\boldsymbol{\Sigma}$
$\displaystyle=\mbox{Cov}(\textbf{z}_{{}_{EUC}})+\mbox{Cov}(\textbf{z}_{{}_{TU}})+\mbox{Cov}(\textbf{z}_{{}_{TD}})+\mbox{Cov}(\boldsymbol{\epsilon})$
$\displaystyle=\textbf{C}_{{}_{EUC}}(h_{{}_{EUC}}|\boldsymbol{\theta}_{\textbf{z}_{EUC}})+\textbf{C}_{{}_{TU}}(h_{{}_{H}}|\boldsymbol{\theta}_{\textbf{z}_{TU}})+\textbf{C}_{{}_{TD}}(h_{{}_{H}}|\boldsymbol{\theta}_{\textbf{z}_{TD}})+\sigma_{0}^{2}\textbf{I},$
(4)

where $h_{{}_{EUC}}$ is the Euclidean distance and $h_{{}_{H}}$ represents the
hydrological distance, that is, distance _along_ the river. Note that,
distance considered can depend on flow-connectivity of sites and the
appropriate covariance function but for brevity we denote this as
$h_{{}_{H}}$, see Santos-Fernandez et al. (2021) for further details. See
Appendix A.1 for the specification of the covariance matrices
$\textbf{C}_{{}_{EUC}}$, $\textbf{C}_{{}_{TU}}$ and $\textbf{C}_{{}_{TD}}$, as
well as distance $h_{{}_{H}}$.

The generalised linear mixed model for river networks we consider is of the
form,

Y $\displaystyle\sim\mathcal{N}(\boldsymbol{\mu},\boldsymbol{\Sigma})$ (5)
$\displaystyle
p(\textbf{y}|\boldsymbol{\theta},\textbf{d})=(2\pi)^{-\frac{k}{2}}\mbox{det}(\boldsymbol{\Sigma})^{-\frac{1}{2}}$
$\displaystyle\mbox{exp}\big{(}-\frac{1}{2}(\textbf{y}-\boldsymbol{\mu})^{T}\boldsymbol{\Sigma}^{-1}(\textbf{y}-\boldsymbol{\mu})\big{)},$

with link function as the identity
$g(\boldsymbol{\mu})=\textbf{X}\boldsymbol{\beta}$ and $\boldsymbol{\Sigma}$
as defined by Equation (4) at a collection of locations d on the network.

_Example 2: Coral Reef Modelling_

Estimates of hard coral are generally obtained using a number of images, $N$,
of the reef floor at locations $s_{k}$ collected along a transect. In order to
define transect placement, design parameters are introduced: the midpoint of a
transect using Easting and Northing coordinates, $E^{j}$ and $N^{j}$, the
angle of the transect, $\alpha^{j}$ in degrees, and the length of the
transect, $l^{j}$ in meters (m). Transects can therefore be fully specified by
design parameters $d^{j}=(E^{j},N^{j},\alpha^{j},l^{j})$ for transect $j$ (see
Figure 1). Let $s_{k}=(e_{k},n_{k})$ denote the geographic coordinates of
image $k$ specified by a collection of transects d. Another parameter,
$r^{j}$, is introduced to create a transect _region_. For radius $r^{j}>0$,
points disperse across the geographic plane to
$\hat{s}_{k}=(\hat{e}_{k},\hat{n}_{k})=(e_{k}+\delta_{k}^{1},n_{k}+\delta_{k}^{2})$,
with $\delta_{k}^{1},\delta_{k}^{2}\sim\mbox{Unif}(-r^{j},r^{j})$.

For each image a number, $n_{k}$, of randomly selected points are placed on
image, $k$, and classified as hard coral, or not. The number of points within
an image that contain hard coral, $y_{k}$, are assumed to have a binomial
distribution

Y $\displaystyle\sim\mbox{Binomial}(\textbf{q},\textbf{n})$ $\displaystyle
p(\textbf{y}|\boldsymbol{\theta},\textbf{d},\textbf{z})$
$\displaystyle=\prod^{N}_{k=1}{n_{k}\choose
y_{k}}q_{k}^{y_{k}}(1-q_{k})^{n_{k}-y_{k}}.$ (6)

where
$\textbf{q}=\mathbb{E}(\textbf{Y}|\textbf{X},\boldsymbol{\theta},\textbf{z})$,
with spatially correlated random effects
$\textbf{z}\sim\mathcal{N}(\bf{0},\boldsymbol{\Sigma_{z}}+\sigma_{0}^{2}\textbf{I})$
where the matrix $\boldsymbol{\Sigma_{z}}$ has entries of the form in Equation
(2) or some other covariance matrix. The expected value of coral cover $q_{k}$
on image $k$ leads to the specification of a logistic regression model with
link function denoted by
$g(\textbf{q})=\mbox{log}(\textbf{q}/(1-\textbf{q}))$, so that
$\textbf{q}=\mbox{logit}^{-1}(\boldsymbol{\eta}+\textbf{z})=\mbox{exp}(\boldsymbol{\eta}+\textbf{z})/(1+\mbox{exp}(\boldsymbol{\eta}+\textbf{z}))$
with $\boldsymbol{\eta}=\textbf{X}\boldsymbol{\beta}$.

Figure 1: Data is sampled at each geographic point $s_{k}$ along the transect
defined by $d^{j}=(E^{j},N^{j},\alpha^{j},l^{j})$ (left). A transect region
with radius $r^{j}$ is shown, where data is collected at the stochastic sample
locations $\hat{s}_{k}$ within the region (right).

###  2.2 Bayesian Framework

Let $p(\boldsymbol{\theta})$ denote the prior distribution on
$\boldsymbol{\theta}$, and consider the joint posterior distribution,

$p(\boldsymbol{\theta},\textbf{z}|\textbf{y},\textbf{d})=\frac{p(\textbf{y}|\boldsymbol{\beta},\textbf{d},\textbf{z})p(\boldsymbol{\theta})p(\textbf{z}|\boldsymbol{\theta}_{z})}{p(\textbf{y}|\textbf{d})}\propto
p(\textbf{y}|\boldsymbol{\beta},\textbf{d},\textbf{z})p(\boldsymbol{\theta})p(\textbf{z}|\boldsymbol{\theta}_{z}),$
(7)

where
$p(\textbf{y}|\textbf{d})=\int_{\theta}p(\textbf{y},\boldsymbol{\theta})\
d\boldsymbol{\theta}=\int_{\theta}p(\textbf{y}|\boldsymbol{\theta},\textbf{d})p(\boldsymbol{\theta})\
d\boldsymbol{\theta}$ is the marginal distribution of y. Unless the likelihood
and spatial random effect are analytically tractable, that is,
$\int_{\textbf{z}}p(\textbf{y}|\boldsymbol{\beta},\textbf{d},\textbf{z})p(\textbf{z}|\boldsymbol{\theta}_{z})\
d\textbf{z}=p(\textbf{y}|\boldsymbol{\theta},\textbf{d})$ has a closed form
solution, the spatial random effect must be integrated out numerically. We
describe the approach taken in the next section.

##  3 Bayesian Optimal Design for Spatial Models

Optimal design problems are commonly viewed as an optimisation framework with
respect to an experimental aim. A utility function is defined to quantify the
worth of a design, which may be to control systematic error (bias), reduce
random variation, increase the precision of parameter estimates, make
predictions about future observations or discriminate between competing
models. Formally, we denote the utility function as
$u(\textbf{d},\boldsymbol{\theta},\textbf{y})$ which depends on model
parameters $\boldsymbol{\theta}\sim p(\boldsymbol{\theta})$ and the data we
might observe under that statistical model $\textbf{y}\sim
p(\textbf{y}|\boldsymbol{\theta},\textbf{d})$ at design points d from the
design space $\mathcal{D}$. We establish the notation for a design as follows,
$\textbf{d}_{i}=(d_{i}^{1},...,d_{i}^{\gamma})$, and use $d_{i}^{j}$ to refer
to the $j$th coordinate of design $\textbf{d}_{i}$ throughout the paper.

The aim is to maximise the expected utility,
$U(\textbf{d})=\mathbb{E}[u(\textbf{d},\boldsymbol{\theta},\textbf{y})]$ over
the entire parameter space $\Theta$ and all conceivable data sets
$\mathcal{Y}$. The Bayesian optimal design $\textbf{d}_{*}$ can therefore be
expressed as,

$\displaystyle\textbf{d}_{*}$ $\displaystyle=\mbox{arg
max}_{\textbf{d}\in\mathcal{D}}\mathbb{E}[u(\textbf{d},\boldsymbol{\theta},\textbf{y})]$
$\displaystyle=\mbox{arg
max}_{\textbf{d}\in\mathcal{D}}\int_{\textbf{y}\in\mathcal{Y}}\int_{\boldsymbol{\theta}\in\Theta}u(\textbf{d},\boldsymbol{\theta},\textbf{y})\
p(\boldsymbol{\theta},\textbf{y}|\textbf{d})\ d\boldsymbol{\theta}\
d\textbf{y}$ (8) $\displaystyle=\mbox{arg
max}_{\textbf{d}\in\mathcal{D}}\int_{\textbf{y}\in\mathcal{Y}}\int_{\boldsymbol{\theta}\in\Theta}u(\textbf{d},\boldsymbol{\theta},\textbf{y})\
p(\textbf{y}|\boldsymbol{\theta},\textbf{d})\ p(\boldsymbol{\theta})\
d\boldsymbol{\theta}\ d\textbf{y}.$

Unfortunately, the integral in Equation (8) is analytically intractable for
most applications, and indeed the linear mixed-model in Equation (1). In
practice, a commonly used method is Monte Carlo integration, where the
expected utility can be approximated by the empirical mean, that is,

$\hat{U}(\textbf{d})=\frac{1}{M}\sum^{M}_{m=1}u(\textbf{d},\boldsymbol{\theta}^{(m)},\textbf{y}^{(m)})$
(9)

with draws from the prior $\boldsymbol{\theta}^{(m)}\sim
p(\boldsymbol{\theta})$ and then the likelihood $\textbf{y}^{(m)}\sim
p(\textbf{y}|\boldsymbol{\theta}^{(m)},\textbf{d})$. In order to accurately
estimate $\hat{U}(\textbf{d})$, $M$ needs to be large. In some cases, the
median
$\tilde{U}(\textbf{d})=\mbox{median}(u(\textbf{d},\boldsymbol{\theta}^{(m)},\textbf{y}^{(m)}))$
may be preferred to the arithmetic mean in Equation (9), since it is less
sensitive to tail behaviour when the target distribution is asymmetric. We
explore this comparison further in the following sections.

A commonly used utility function is the Kullback-Liebler (KL) divergence
(Kullback and Leibler, 1951) of the posterior from the prior distribution,

$\displaystyle u(\textbf{y},\textbf{d})$
$\displaystyle=D_{KL}[p(\boldsymbol{\theta}|\textbf{y},\textbf{d})||p(\boldsymbol{\theta})]$
$\displaystyle=\int
p(\boldsymbol{\theta}|\textbf{y},\textbf{d})\log\frac{p(\boldsymbol{\theta}|\textbf{y},\textbf{d})}{p(\boldsymbol{\theta})}\
\mbox{d}\theta.$ (10)

The intuition behind Equation (3) is that a large $D_{KL}$ divergence from
posterior to prior implies that the data in y reduces entropy (randomness or
uncertainty) in $\theta$ and hence the data are more informative at design
points d. As suggested by Lindley (1956), this utility can be used in the
interest of maximising the expected information gain on the model parameters
(or functions of these) when performing experiments at design points d. When
both the prior and the posterior distributions follow multivariate normal
distributions, that is,
$\boldsymbol{\theta}\sim\mathcal{N}(\boldsymbol{\mu}_{0},\boldsymbol{\Sigma_{0}})$
and
$\boldsymbol{\theta}|\textbf{y},\textbf{d}\sim\mathcal{N}(\boldsymbol{\mu}_{1},\boldsymbol{\Sigma_{1}})$,
respectively, then the KL divergence becomes

$u(\textbf{y},\textbf{d})=\frac{1}{2}\Big{[}(\boldsymbol{\mu}_{0}-\boldsymbol{\mu}_{1})^{T}\boldsymbol{\Sigma_{0}}^{-1}(\boldsymbol{\mu}_{0}-\boldsymbol{\mu}_{1})+\text{tr}(\boldsymbol{\Sigma_{0}}^{-1}\boldsymbol{\Sigma_{1}})-\text{ln}\frac{|\boldsymbol{\Sigma_{1}}|}{|\boldsymbol{\Sigma_{0}}|}-{\kappa}\Big{]},$

where $\boldsymbol{\mu}_{0},\ \boldsymbol{\mu}_{1}\in\mathbb{R}^{\kappa}$ and
$\boldsymbol{\Sigma_{0}},\
\boldsymbol{\Sigma_{1}}\in\mathbb{R}^{{\kappa}\times{\kappa}}$ (Duchi, 2007).

Note that in Bayesian design, the utility is always a function of the
posterior distribution, which generally means that inference is performed many
thousands of times since the posterior
$p(\boldsymbol{\theta}|\textbf{y},\textbf{d})$ must be evaluated for each
future data set $\\{\boldsymbol{\theta}^{(m)},\textbf{y}^{(m)}\\}$ that is
drawn from the joint distribution
$p(\boldsymbol{\theta},\textbf{y}|\textbf{d})$. As a computationally efficient
approximation to the posterior distribution, we employ the Laplace
approximation which has the following form:

$\boldsymbol{\theta}|\textbf{y},\textbf{d}\sim\mathcal{N}(\boldsymbol{\theta},\textbf{H}(\boldsymbol{\theta}^{*})^{-1}),$
(11)

where $\boldsymbol{\theta}^{*}=\mbox{arg
max}_{\boldsymbol{\theta}\in\Theta}\\{\mbox{log
}p(\textbf{y}|\boldsymbol{\theta},\textbf{d})+\mbox{log
}p(\boldsymbol{\theta})\\}$ and $\textbf{H}(\boldsymbol{\theta}^{*})$ is the
Hessian matrix defined as:

$\textbf{H}(\boldsymbol{\theta}^{*})=\left.\frac{\displaystyle-\partial^{2}\\{\mbox{log
}p(\textbf{y}|\boldsymbol{\theta},\textbf{d})+\mbox{log
}p(\boldsymbol{\theta})\\}}{\displaystyle\partial\boldsymbol{\theta}\partial\boldsymbol{\theta}^{{}^{\prime}}}\right\rvert_{\boldsymbol{\theta}=\boldsymbol{\theta}^{*}}.$

This approximation requires evaluation of the full data likelihood. Thus, the
full data likelihood used in forming the Laplace approximation is replaced
with a numerical approximation to the posterior for $\boldsymbol{\theta}$,

$p(\boldsymbol{\theta}|\textbf{y},\textbf{d})\propto
p(\boldsymbol{\theta})\int_{\textbf{z}}p(\textbf{y}|\boldsymbol{\beta},\textbf{d},\textbf{z})p(\textbf{z}|\boldsymbol{\theta}_{z})\
d\textbf{z}.$ (12)

Finally, an approach for the maximisation required in finding the Bayesian
optimal design in Equation (8) needs to be specified. It is common to use a
Coordinate Exchange (CE) algorithm to search a discrete space, by exchanging
coordinates in the design, d, with other points in the design space, in a
sequential manner. Formally, denote a proposal design by
$\textbf{d}^{j\curvearrowright}=(d^{1},...,d^{j-1},d^{\curvearrowright},d^{j+1},...,d^{\gamma})$,
obtained by exchanging the $j$th design coordinate for some other design point
$d^{\curvearrowright}\in\mathcal{D}\setminus\textbf{d}$. Given the stochastic
nature of the utility approximation, illustrated by draws of the prior and
likelihood in Equation (8), the search algorithm needs to accommodate
different realisations of $\tilde{U}(\textbf{d})$. For this, an acceptance
probability $p^{*}$ can be used to compare a proposed design
$\textbf{d}^{j\curvearrowright}$ to the current best design d. The Approximate
Coordinate Exchange (ACE) algorithm (Overstall and Woods, 2017) is a commonly
used stochastic CE algorithm with acceptance probability defined by

$p^{*}=1-T_{2B-2}\Big{(}-\dfrac{B\hat{U}(\textbf{d}^{j\curvearrowright})-B\hat{U}(\textbf{d})}{\sqrt{2B\hat{v_{b}}}}\Big{)},$

and

$\hat{v_{b}}=\dfrac{\sum^{B}_{l=1}\big{[}u(\textbf{d}^{j\curvearrowright},\boldsymbol{\theta}^{(l)},\textbf{y}^{(l)})-\hat{U}(\textbf{d}^{j\curvearrowright})\big{]}^{2}+\sum^{B}_{l=1}\big{[}u(\textbf{d},\boldsymbol{\theta}^{(l)},\textbf{y}^{(l)})-\hat{U}(\textbf{d})\big{]}^{2}}{2B-2}$

with large $B$ and random draws
$\\{\boldsymbol{\theta}^{(l)},\textbf{y}^{(l)}\\}$ from
$p(\boldsymbol{\theta},\textbf{y}|\textbf{d})$. However, this acceptance
criteria relies on the assumption of a normally distributed utility (at least
approximately).

We propose an alternative version of the stochastic Coordinate Exchange
algorithm (see Appendix Algorithm 2) with a non-parametric acceptance
criterion defined as

$p_{W}^{*}=1-W(U(\textbf{d}^{j\curvearrowright}),U(\textbf{d}))$

where $W$ is the p-value of the one-sided Wilcoxon rank sum test (equivalent
to the Mann-Whitney test). In this case, the null hypothesis is that the
distributions of the proposed design utility,
$U(\textbf{d}^{j\curvearrowright})$, and the current best design utility,
$U(\textbf{d})$, differ by a location shift of $\mu=0$ and the alternative
hypothesis is that $\mu>0$. The ‘one-sided’ specifies that
$U(\textbf{d}^{j\curvearrowright})$ is shifted to the right of
$U(\textbf{d})$. We explore the proposed non-parametric stochastic Coordinate
Exchange algorithm in comparison to the ACE in the Appendix, both with median
and mean central tendency measure.

##  4 Sampling Windows

As mentioned previously, optimal design methods have been used to determine
sampling windows for experiments in pharmacokinetics (Foo et al., 2012;
Bogacka et al., 2008). For spatial processes, sampling windows would allow the
experimenter some flexibility in the sampling location, while preserving a
minimum required efficiency for the design utility. In order to define
sampling windows, we define the utility in a continuous domain. Consider the
design space $\mathcal{D}\subseteq\mathbb{R}^{q}$ across $q$ windows. We first
introduce Gaussian processes (GP) and illustrate their flexibility and
capability as a multi-dimensional utility smoothing technique. Based on the
smoothing GPs provide, we then show how design efficiency contours can be
derived.

GPs generalise the concept of Gaussian distributions over discrete random
variables to the function space. In the Bayesian context, GPs can be seen as a
prior over the function space, where inference takes place. Let
$f(\textbf{x})$ denote the target value of interest where
$\textbf{x}\in\mathbb{R}^{q}$ is a $q$-dimensional random vector containing
predictor values such as spatial locations. We write
$\mathcal{GP}(m(\textbf{x}),k(\textbf{x},\textbf{x}^{{}^{\prime}}))$ to denote
a GP characterised by a mean function $m(\textbf{x})$ and a covariance
function $k(\textbf{x},\textbf{x}^{{}^{\prime}})$. Consider a GP model with
kernel,

$k(\textbf{x}_{p},\textbf{x}_{s})=\sum\limits^{q}_{j=1}\exp(-\zeta_{j}*\mbox{dist}_{j}(\textbf{x}_{p},\textbf{x}_{s}))$
(13)

with hyperparameters $\boldsymbol{\zeta}$ and distance measure defined as
$\mbox{dist}_{j}(\textbf{x}_{p},\textbf{x}_{s})=|x_{p}^{j}-x_{s}^{j}|$, for
$j=\\{1,...,q\\}$. Using a finite collection of inputs
$\\{\textbf{x}_{1},...,\textbf{x}_{n}\\}$ and the function above, we construct
the kernel matrix $\textbf{K}(\cdot)$ over the parameter space.

Details of the proposed Sampling Windows Algorithm are presented in Algorithm
1, using design notation. We begin by defining a design space for windows
(line 2), this may differ from the design space searched in the discrete
stochastic CE algorithm. We fit a GP (line 4) to the criterion of interest,
here, the median utility $\tilde{U}(\textbf{d}_{i})$ computed at each input
$\textbf{d}_{i}$ (line 3), assuming a zero-mean GP prior,
$f\sim\mathcal{GP}(\textbf{0},\textbf{K}(\textbf{D},\textbf{D})+\zeta_{0}\textbf{I})$.
Inclusion of the nugget $\zeta_{0}>0$ ensures that every $f$ has some
magnitude of variance with itself only (independent noise), so that the Monte
Carlo approximations of the median utility surface will be smoothed, not
interpolated. Hyperparameters
$\boldsymbol{\zeta}=(\zeta_{0},...,\zeta_{\kappa})$ are determined (line 5) by
minimising cross-validation error of the GP,

$\displaystyle\hat{\boldsymbol{\zeta}}$
$\displaystyle=\underset{\boldsymbol{\zeta}}{\mathrm{arg\ min}}\
CV(\boldsymbol{\zeta}|f,\textbf{D})\text{, with }$ (14) $\displaystyle
CV(\boldsymbol{\zeta}|f,\textbf{D})$ $\displaystyle=\frac{\displaystyle
1}{\displaystyle\gamma_{1}}\sum\limits^{\gamma_{1}}_{i=1}\Bigg{(}\frac{\displaystyle\tilde{U}(\textbf{d}_{i}|\textbf{d}_{c})-f(\textbf{d}_{i})}{\displaystyle
1-S_{ii}}\Bigg{)}^{2}$ (15)

where
$\textbf{S}=\textbf{K}(\textbf{D},\textbf{D})[\textbf{K}(\textbf{D},\textbf{D})+\zeta_{0}\textbf{I}]^{-1}$.

The posterior predictive mean of $f$ used as an emulator (in line 6) is,

$\bar{f}=\textbf{K}(\textbf{d}_{\star},\textbf{D})[\textbf{K}(\textbf{D},\textbf{D})+\zeta_{0}\textbf{I}]^{-1}\tilde{U}(\textbf{D}|\textbf{d}_{c})$
(16)

for arbitrary inputs $\textbf{d}_{\star}$ derived using standard results on
the conditional distribution of normal random variables. Finally, computing
the design efficiency (line 7) is achieved by normalising the functional
utility to

$eff(\textbf{d}_{\star j})=\frac{\displaystyle\bar{f}(\textbf{d}_{\star
j}|\textbf{D},\hat{\boldsymbol{\zeta}})}{\displaystyle\bar{f}(\hat{\textbf{d}}_{*j}|\textbf{D},\hat{\boldsymbol{\zeta}})}$
(17)

where $\hat{\textbf{d}}_{\star j}=\underset{\textbf{d}_{\star
j}\in\textbf{D}_{\star}}{\mathrm{arg\ max}}\ \bar{f}(\textbf{d}_{\star
j}|\textbf{D},\hat{\boldsymbol{\zeta}})$. Output design efficiency contours
(line 8) as the set of prediction points above threshold value $t$.

1

2If existing, let
$\textbf{d}_{c}=(\textbf{d}_{c1},\textbf{d}_{c2},...,\textbf{d}_{c\gamma})$ be
the current design. Specify $\tilde{U}(\textbf{d}|\textbf{d}_{c})$ as
appropriate and define prior $p(\boldsymbol{\theta})$ and likelihood
$p(\textbf{y}|\boldsymbol{\theta},\textbf{d})$.

3

_Phase I - Construct $q$-dimensional utility grid_

4

5Define design space $\mathcal{D}\subseteq\mathbb{R}^{q}$ across $q$ windows,
sample $\textbf{D}=(\textbf{d}_{1},...,\textbf{d}_{\gamma_{1}})$ with
$\textbf{d}_{i}\in\mathcal{D}$ for training and
$\textbf{D}_{\star}=(\textbf{d}_{\star 1},...,\textbf{d}_{\star\gamma_{2}})$
with $\textbf{d}_{\star j}\in\mathcal{D}$ for predictions.

6Compute median utility $\tilde{U}(\textbf{d}_{i})$ for each point
$\textbf{d}_{i}$

7

_Phase II - Functional utility surface_

8

9Fit a GP; hyperparameters are determined by minimising cross-validation error

10 $\hat{\boldsymbol{\zeta}}=\underset{\boldsymbol{\zeta}}{\mathrm{arg\ min}}\
CV(\boldsymbol{\zeta}|f,\textbf{d})$ with CV defined in Equation (15)

11Evaluate predicted mean values $\bar{f}$ at each prediction point
$\textbf{d}_{\star j}$ as $\bar{f}(\textbf{d}_{\star
j}|\textbf{D},\hat{\boldsymbol{\zeta}})=\sum\limits^{\gamma_{1}}_{i=1}k(\textbf{d}_{\star
j},\textbf{d}_{i}|\hat{\boldsymbol{\zeta}})\alpha_{i}$ for
$\boldsymbol{\alpha}=[\textbf{K}(\textbf{D},\textbf{D})+\hat{\zeta}_{0}\textbf{I}]^{-1}\tilde{U}(\textbf{D}|\textbf{d}_{c})$

12

_Phase III - Design Efficiency contours_

13

14Normalise utility surface to yield efficiencies $eff(\textbf{d}_{\star j})$
as per Equation (17) for each prediction point $\textbf{d}_{\star j}$

15Output set of points where $eff(\textbf{d}_{\star})>t$ for threshold values
$t=[0,1]$

Algorithm 1 Sampling Windows Algorithm

##  5 Case Studies

###  5.1 Example 1: River Network Design

The design problem in this case study is motivated by a costly and complex
balance of priorities, that is, short-term interest in monitoring particular
river segments (e.g. wildfire impacts on juvenile fish habitat), while
maintaining the ability monitor efficiently at established locations
throughout the region. In this example, we use a water temperature dataset
from the Clearwater River Basin, USA, with observation sites located from 500m
up to 46km apart (Figure 2). Temperature data were collected at 15 spatial
locations using in-situ sensors that recorded measurements at 30 minute
intervals (Isaak et al., 2018). The response variable is mean temperature for
July 1, 2013, with stream slope, elevation, watershed area, and air
temperature used as covariates (Santos-Fernandez et al., 2022). Note that
spatial variation in temperature often constitutes the largest proportion of
total variation, since temporal correlation among sites is strong (Isaak et
al., 2017). A tail-down covariance function was found to be the most suitable
covariance function for describing this spatial dependence.

Figure 2: A total of 15 observation sites located in the Clearwater River
Basin, USA. Average daily water temperature for July 1, 2013. The width of the
stream line is proportional to Shreve’s stream order (Shreve, 1966). The
optimal and non-optimal design locations using the proposed stochastic
Coordinate Exchange search are represented by filled and unfilled dots,
respectively.

There are 15 existing monitoring sites in the Clearwater River Basin and the
design problem is to find the best 5 among them. In addition, wildfires
occurred near two segments that provide cold water refugia for juvenile fish
and sensors are to be deployed to monitor impacts on water temperature. For
convenience, we define search neighbourhoods along the river network, that is,
any continuous line along the branching river network and specify
Neighbourhood 1 (N1) and Neighbourhood 2 (N2) in Figure 4.

Figure 3: Trace plot for the coordinate exchange algorithm using Wilcoxon
acceptance criteria and median utility (top-left); proposal probability using
Wilcoxon acceptance criteria, $p^{*}_{W}$ (bottom-left); every 50th utility
distribution in the proposed CE algorithm (top-right); and the optimal utility
distribution (bottom-right).

The stochastic Coordinate Exchange approach outlined in Algorithm 2 (Appendix)
was computed with $K=5$ random starts and $T=15$ outer iterations. Figure 3
shows a selection of the utility distributions found in the proposed search
algorithm, then the distributions of the optimal design, noting the non-
normality of these distributions. The optimal design was found to be at
locations, $\textbf{d}_{*}=(167,169,172,174,183)$ with a median utility of
3.6805. These optimal locations cover a wide range for each of the covariate
values, in particular, optimal locations include the maximum and minimum value
for air temperature. The set of locations $(167,169,174)$ are clustered at a
river confluence and the most downstream location of all sites, $169$, is
included in the design (Figure 2). These characteristics agree with optimal
designs found using previous pseudo-Bayesian approaches for river networks
Falk et al. (2014). We also benchmark the performance of the proposed
Coordinate Exchange algorithm against other variations; a comparison of
acceptance criteria (Wilcoxon vs. ACE), as well as central tendency measure
(mean vs. median) is provided in the Appendix. The choice of the median as a
preferred measure is supported by the apparent skewness and long tails of the
distributions of the utilities (Figure 3).

Recall the specification of neighbourhoods N1 and N2 within the river network
for sensor deployment. The proposed Sampling Windows Algorithm (Algorithm 1),
is applied. Consider the design space $\mathcal{D}\subseteq\mathbb{R}^{2}$
with design points $d_{i}^{1}\in N1$ and $d_{i}^{2}\in N2$ for each design
$\textbf{d}_{i}=(d_{i}^{1},d_{i}^{2})$. The current design $\textbf{d}_{c}$ is
set as the optimal points $\textbf{d}_{*}$ from the downsizing regime
previously found. In Phase I, the median utility is computed as
$\tilde{U}(\textbf{d}_{i}|\textbf{d}_{c})=\tilde{U}(\textbf{d}_{i}\cup\textbf{d}_{c})$
for each $\textbf{d}_{i}$ in the grid, D, sampled from neighbourhoods 500m
apart. Covariate values at $\textbf{d}_{i}$ for air temperature are
interpolated using 3 nearest neighbours based on Euclidean distance, whereas
the remaining covariate values (stream slope, elevation and watershed area)
are unique to each river segment. In Phase II, a GP is fit assuming a kernel
function defined by Equation (13), where the distance measure
$\mbox{dist}_{j}(\textbf{d}_{s},\textbf{d}_{p})=|d_{s}^{j}-d_{p}^{j}|$ for
$j=1,2$ is equivalent to stream distance (along dimension $j$). Results from
Phase III indicate that N1 is more sensitive to location placement than N2, as
shown in Figure 5.

From an applied standpoint, these design outputs allow decisions to be made in
the field (e.g. placement around specific location $n_{15}$ in N1, but
anywhere convenient on the main branch in N2). This provides support for
greater flexibility in sensor placement when access issues arise, without
having to rerun a full discrete optimal design search. A practical approach,
given the design efficiency contours as a roadmap, would be to visit N1 first;
should $n_{15}$ be inaccessible then move along the river to deploy the in-
situ sensor where possible, say $n_{14}$. Take the line segment corresponding
to $eff(\textbf{D}_{*}|n_{14})$ to inform deployment in N2, in this case,
$n_{24}$, to retain $98\%$ optimal utility.

Figure 4: Existing sites (red dots) and potential sites within search
neighbourhoods N1 and N2 (gray dots) located on the river network. A discrete
set of points were sampled from both neighbourhoods spaced at 500m stream
distance apart. Figure 5: Phases of the proposed Sampling Windows Algorithm;
constructed median utility grid (left), Gaussian process mean predictions
(middle), and normalised contours of design efficiency (right). A sensor close
to $n_{15}\in$N1 will yield high design efficiency, while more flexibility can
be afforded for the choice of sensor placement in N2.

###  5.2 Coral Reef Design

The Australian Institute of Marine Science (AIMS) has monitored the condition
of coral reefs around Australia for decades. We considered data collected in
2013 from the Barracouta East submerged coral reef off the north-west coast of
Australia (Heyward et al., 2013). The design problem is to find the best
regions to sample underwater images in. We propose a two step approach using
the methods described above: first, a global search for the optimal transect
lines with $n=3$, the total number of transects, using the proposed stochastic
CE algorithm in a discrete design space; then, a continuous search for the
optimal radius of each transect region with $q=3$, the number of sampling
windows centred at each transect.

Coral cover is assessed using images of the reef floor taken every 5m along a
transect, where $n_{k}=20$ random points on image $k$ are analysed, leading to
the consideration of the logistic regression model described in Equation (6).
Exploratory analyses showed a strong linear relationship between the response
variable, q, and depth of the reef floor which was used as the only covariate.
Note that depth was interpolated using 3 nearest neighbours based on Euclidean
distance from the survey data. To reduce the design space and thus the
complexity of the search, angles are discretised to
$\alpha^{j}=\\{0,45,90,135\\}$, $E^{j}$ and $N^{j}$ are defined by a discrete
set of points. Length remains fixed at $l^{j}=500$ so that each transect
compromises of 100 images.

The proposed Stochastic Coordinate Exchange algorithm with the Wilcoxon
acceptance criteria and median utility (Algorithm 2) was applied to find the
optimal transect lines across the shoal. When computing the utility, the
Laplace approximation to the posterior requires evaluating the full data
likelihood, as shown in Equation (12). Accordingly, the spatially correlated
random effect, z, needs to be integrated out. For this, we chose to use Monte
Carlo integration with $M_{z}=50$ draws from the marginal likelihood due to
computational restraints. The Monte Carlo integration for the median utility
in Equation (9) is evaluated with $M=350$ and $B=600$ draws from the prior and
full likelihood for the evaluation of
$\hat{U}(\textbf{d}^{j\curvearrowright})$ and $p^{*}_{W}$, respectively.

Figure 6: Trace plot of proposed stochastic coordinate exchange algorithm
(top-left), acceptance probability (bottom-left), and optimal transect design
in geographic space (right).

Spatial covariance is characterised by exponential decay, defined in Equation
(2). To further ease the computational burden, we construct a spatial grid (10
$\times$ 10) over the entire shoal, with $h_{{}_{EUC}}$ the Euclidean distance
between grid centre points (a lower-dimensional surrogate to the $N\times N$
distance matrix using coordinates $\textbf{x}_{k}$). This grid size is a
compromise; not so small that the computational gains are marginal, but not so
large that the spatial relationships are not captured between grid rectangles.
Results in Figure 6 show the optimal placement of transects
$\textbf{d}_{*}=(d_{*}^{1},d_{*}^{2},d_{*}^{3})$ in the shallower areas of the
reef, each with varying depth gradients. Thus, it would seem reasonable to
sample in along these transects to estimate coral cover.

The proposed Sampling Windows Algorithm (Algorithm 1) is applied over the
parameter space $r$, with the kernel function defined by Equation (13) and
distance $\mbox{dist}_{j}(r_{s},r_{p})=|r_{s}^{j}-r_{p}^{j}|$ for $j=1,2,3$.
In the first step, the optimal transect design previously found is fixed,
while the design space $\mathcal{D}\subseteq\mathbb{R}^{3}$ is constructed for
radius parameter, $r$, corresponding to each transect. Since stochasticity is
introduced in the choice of geographic locations, $\hat{s}_{k}$, we average
median utility across 3 location samples for a given radius, each with $M=600$
Monte Carlo draws from the prior and likelihood. The 3-dimensional design
contours are shown in Figure 7. Optimal design efficiency,
$eff(\textbf{D}_{*})=101\%$, is achieved at $r_{*}=(44,0,44)$. This implies
optimality when sampling within a region for the shallower areas, and more
strictly across the transect $d_{*}^{2}$ in the deeper area of the reef.
Generally, these contours can be used to guide decision making on the field,
for example, sample $d_{*}^{2}$ when conditions are good, whereas sampling in
$d_{*}^{1}$ and $d_{*}^{3}$ can be done with stronger underwater currents.

Figure 7: The median utility (left); Gaussian process mean predictions of
median utility (middle), and normalised contours of design efficiency (right)
found at each phase of the proposed Sampling Windows Algorithm. The radius
parameter, $r^{j}$, specifies a radius for transect $d^{j}_{*}$ for
$j=\\{1,2,3\\}$.

##  6 Discussion

Coral reefs and freshwater ecosystems are biodiversity hotspots that are
particularly sensitive to anthropogenic impacts of climate change, water
pollution and over-exploitation, among other factors, resulting in rapid
population declines (Strayer and Dudgeon, 2010; Tickner et al., 2020; Wagner
et al., 2020). Informative data are needed to support critical and timely
conservation decisions, but various practical constraints make precise
sampling difficult or even impossible. In this paper, we have considered
Bayesian optimal design for collecting data where greater flexibility is
needed in the design outputs. To this end, we proposed a spatial Sampling
Windows Algorithm using Gaussian process emulation to estimate design
efficiency contours. While we focused on two important environmental systems,
such an approach is also applicable to a broad range of complex systems that
require rigorous design. For example, optimal design methods for sub-sampling
can be used as a computationally efficient way to analyse Big Data (Drovandi
et al., 2017). Rather than extracting optimal design points from the Big Data
set, which may or may not exist, points could be selected within windows.

There are several novel aspects to the approaches described here. Statistical
models for branching river networks use distance travelled along the river
network, directionality, and flow volume to describe spatial relationships.
There is novelty in the application of Bayesian design for river network
systems, as of yet to be presented in the statistical literature. Since
utility distributions were not normally distributed, we also proposed the use
of a non-parametric acceptance probability in the stochastic Coordinate
Exchange, based on the one-sided Wilcoxon Rank Sum Test. We successfully
demonstrated how Bayesian design works in complex spatial settings and
introduced a method to find sampling windows, facilitating intelligent data
collection.

The novel methodological advances presented here provide practical solutions
to common monitoring challenges. Results from the proposed Sampling Windows
Algorithm highlighted in the river network case study showed that design
efficiency was best achieved around a smaller range in Neighbourhood 1, while
almost agnostic to choice in Neighbourhood 2; recalling that neighbourhoods
are any continuous line along the river. In particular, the design efficiency
contours identify which, and quantify how much, design points are sensitive to
precise sampling. Similarly, results in the coral reef case study showed high
design efficiency around regions in shallow areas of the reef compared to
deeper regions with less hard coral. Providing fine-scale information about
optimality allows practitioners in the field to make informed judgements on
jointly optimal variables based on physical conditions not included in the
utility (e.g. accessibility, safety, etc.), without the need to re-run time
consuming optimisations.

There are also numerous opportunities for future research in this space. One
useful extension would be to incorporate additional monitoring objectives
related to cost into the utility function. For example, the optimal deployment
of a combination of high-cost sensors that produce high-quality data and
lower-cost sensors that produce less reliable data. Sampling and monitoring
costs are a core consideration for any monitoring program, but it is
especially true in remote locations. In such cases, it would also be useful to
consider the travel time between sampling locations in the optimal design. New
Bayesian spatio-temporal models for river networks (Santos-Fernandez et al.,
2021) provide the opportunity to extend the approach developed here for
spatio-temporal sampling windows. This would allow practitioners to make more
effective use of mobile in-situ sensors for spatio-temporal data collection.
Thus, the novel methods presented here (along with future extensions)
facilitate more efficient data collection; providing insights into spatial
(and spatio-temporal) processes which can be used to better inform data-
enabled management decisions in complex and vulnerable ecosystems.

## Acknowledgments

This project was supported by the Australian Research Council (ARC) Linkage
Project (LP180101151) “Revolutionising water-quality monitoring in the
information age”. We thank Dan Isaak and Dona Horan from the US Forest Service
Rocky Mountain Research Station for sharing the data and creating the
SpatialStreamNetwork object used in the river network case study. We also
acknowledge the Seascape Health and Resilience team of the Australian
Institute of Marine Science (AIMS) for the curation and provision of data
towards the coral reef case study.

## Code Accessibility

The computations were performed using the R package bayesDesign, which
contains functions for fitting models, optimisations and visualising design
domains. The code can be found online at github.com/KatieBuc/bayesDesign. Note
that the code relies on high performance computing infrastructure.

## References

  * AIMS (2021) AIMS (2021) Annual Summary Report of Coral Reef Condition 2020/21.  (Available from https://www.aims.gov.au/reef-monitoring/gbr-condition-summary-2020-2021).
  * Atkinson and Donev (1992) Atkinson, A. C. and Donev, A. N. (1992) Optimum experimental designs.  Oxford University Press, New York.
  * Bogacka et al. (2008) Bogacka, B., Johnson, P., Jones, B. and Volkov, O. (2008) D-efficient window experimental designs.  Journal of statistical planning and inference, 138, 160–168.
  * Carlon et al. (2020) Carlon, A. G., Dia, B. M., Espath, L., Lopez, R. H. and Tempone, R. (2020) Nesterov-aided stochastic gradient methods using Laplace approximation for Bayesian design optimization.  Computer Methods in Applied Mechanics and Engineering, 363, 112909.
  * Chaloner and Verdinelli (1995) Chaloner, K. and Verdinelli, I. (1995) Bayesian experimental design: A review.  Statistical Science, 273–304.
  * Drovandi et al. (2017) Drovandi, C. C., Holmes, C., McGree, J. M., Mengersen, K., Richardson, S. and Ryan, E. G. (2017) Principles of experimental design for Big Data analysis.  Statistical Science: A Review Journal of the Institute of Mathematical Statistics, 32, 385.
  * Duchi (2007) Duchi, J. (2007) Derivations for linear algebra and optimization.  Berkeley, California, 3, 2325–5870.
  * Falk et al. (2014) Falk, M. G., McGree, J. M. and Pettitt, A. N. (2014) Sampling designs on stream networks using the pseudo-Bayesian approach.  Environmental and Ecological Statistics, 21, 751–773.
  * Foo et al. (2012) Foo, L. K., McGree, J. and Duffull, S. (2012) A general method to determine sampling windows for nonlinear mixed effects models with an application to population pharmacokinetic studies.  Pharmaceutical statistics, 11, 325–333.
  * Hastings (1970) Hastings, W. K. (1970) Monte Carlo sampling methods using Markov chains and their applications.  Biometrika, 57, 97–109.
  * Heyward et al. (2013) Heyward, A., Speed, C., Meekan, M., Cappo, M., Colquhoun, J., Fisher, R., Meeuwig, J. and Radford, B. (2013) Montara: Barracouta east, goeree and vulcan shoals survey 2013.  Report for PTTEP Australasia (Ashmore Cartier) Pty Ltd. Australian Institute of Marine Science, Perth.
  * Isaak et al. (2018) Isaak, D. J., Luce, C. H., Chandler, G. L., Horan, D. L. and Wollrab, S. P. (2018) Principal components of thermal regimes in mountain river networks.  Hydrology and Earth System Sciences, 22, 6225–6240.
  * Isaak et al. (2017) Isaak, D. J., Wenger, S. J., Peterson, E. E., Ver Hoef, J. M., Nagel, D. E., Luce, C. H., Hostetler, S. W., Dunham, J. B., Roper, B. B., Wollrab, S. P. et al. (2017) The NorWeST summer stream temperature model and scenarios for the western US: A crowd-sourced database and new geospatial tools foster a user community and predict broad climate warming of rivers and streams.  Water Resources Research, 53, 9181–9205.
  * Isaak et al. (2016) Isaak, D. J., Young, M. K., Luce, C. H., Hostetler, S. W., Wenger, S. J., Peterson, E. E., Ver Hoef, J. M., Groce, M. C., Horan, D. L. and Nagel, D. E. (2016) Slow climate velocities of mountain streams portend their role as refugia for cold-water biodiversity.  Proceedings of the National Academy of Sciences, 113, 4374–4379.
  * Jackson et al. (2018) Jackson, F. L., Fryer, R. J., Hannah, D. M., Millar, C. P. and Malcolm, I. A. (2018) A spatio-temporal statistical model of maximum daily river temperatures to inform the management of Scotland’s Atlantic salmon rivers under climate change.  Science of the Total Environment, 612, 1543–1558.
  * Kang et al. (2016) Kang, S. Y., McGree, J. M., Drovandi, C. C., Caley, M. J. and Mengersen, K. L. (2016) Bayesian adaptive design: improving the effectiveness of monitoring of the Great Barrier Reef.  Ecological applications, 26, 2637–2648.
  * Kullback and Leibler (1951) Kullback, S. and Leibler, R. A. (1951) On information and sufficiency.  The Annals of Mathematical Statistics, 22, 79–86.
  * Lin et al. (2011) Lin, Y.-P., Chu, H.-J., Huang, Y.-L., Cheng, B.-Y. and Chang, T.-K. (2011) Modeling spatial uncertainty of heavy metal content in soil by conditional latin hypercube sampling and geostatistical simulation.  Environmental Earth Sciences, 62, 299–311.
  * Lindenmayer and Likens (2010) Lindenmayer, D. and Likens, G. (2010) The science and application of ecological monitoring.  Biological Conservation, 143, 1317–1328.
  * Lindley (1956) Lindley, D. V. (1956) On a measure of the information provided by an experiment.  The Annals of Mathematical Statistics, 27, 986–1005.
  * Long et al. (2013) Long, Q., Scavino, M., Tempone, R. and Wang, S. (2013) Fast estimation of expected information gains for Bayesian experimental designs based on Laplace approximations.  Computer Methods in Applied Mechanics and Engineering, 259, 24–39.
  * McGree et al. (2012) McGree, J., Drovandi, C. and Pettitt, A. (2012) A sequential Monte Carlo approach to derive sampling times and windows for population pharmacokinetic studies.  Journal of pharmacokinetics and pharmacodynamics, 39, 519–526.
  * Overstall et al. (2018) Overstall, A. M., McGree, J. M. and Drovandi, C. C. (2018) An approach for finding fully Bayesian optimal designs using normal-based approximations to loss functions.  Statistics and Computing, 28, 343–358.
  * Overstall and Woods (2017) Overstall, A. M. and Woods, D. C. (2017) Bayesian design of experiments using approximate coordinate exchange.  Technometrics, 59, 458–470.
  * Pearse et al. (2020) Pearse, A. R., McGree, J. M., Som, N. A., Leigh, C., Maxwell, P., Ver Hoef, J. M. and Peterson, E. E. (2020) Ssndesign—an r package for pseudo-Bayesian optimal and adaptive sampling designs on stream networks.  PloS one, 15, e0238422.
  * Peterson and Fausch (2003) Peterson, D. P. and Fausch, K. D. (2003) Upstream movement by nonnative brook trout (salvelinus fontinalis) promotes invasion of native cutthroat trout (oncorhynchus clarki) habitat.  Canadian Journal of Fisheries and Aquatic Sciences, 60, 1502–1516.
  * Peterson and Hoef (2010) Peterson, E. E. and Hoef, J. M. V. (2010) A mixed-model moving-average approach to geostatistical modeling in stream networks.  Ecology, 91, 644–651.
  * Peterson et al. (2013) Peterson, E. E., Ver Hoef, J. M., Isaak, D. J., Falke, J. A., Fortin, M.-J., Jordan, C. E., McNyset, K., Monestiez, P., Ruesch, A. S., Sengupta, A., Som, N., Steel, E. A., Theobald, D. M., Torgersen, C. E. and Wenger, S. J. (2013) Modelling dendritic ecological networks in space: an integrated network perspective.  Ecology Letters, 16, 707–719.
  * Pronzato and Walter (1985) Pronzato, L. and Walter, É. (1985) Robust experiment design via stochastic approximation.  Mathematical Biosciences, 75, 103–120.
  * Richmond (1993) Richmond, R. H. (1993) Coral reefs: present problems and future concerns resulting from anthropogenic disturbance.  American Zoologist, 33, 524–536.
  * Roelfsema and Phinn (2010) Roelfsema, C. and Phinn, S. R. (2010) Calibration and validation of coral reef benthic community maps: integration of field data with high spatial resolution multi spectral satellite imagery.  Journal of Applied Remote Sensing, 4, 043527.
  * Ryan et al. (2016) Ryan, E. G., Drovandi, C. C., McGree, J. M. and Pettitt, A. N. (2016) A review of modern computational algorithms for Bayesian optimal design.  International Statistical Review, 84, 128–154.
  * Santos-Fernandez et al. (2022) Santos-Fernandez, E., Hoef, J. M. V., McGree, J. M., Isaak, D. J., Mengersen, K. and Peterson, E. E. (2022) SSNbayes: An R package for Bayesian spatio-temporal modelling on stream networks.  arXiv preprint arXiv:2202.07166.
  * Santos-Fernandez et al. (2021) Santos-Fernandez, E., Hoef, J. M. V., Peterson, E. E., McGree, J., Isaak, D. and Mengersen, K. (2021) Bayesian spatio-temporal models for stream networks.  arXiv preprint arXiv:2103.03538.
  * Senarathne et al. (2020) Senarathne, S., McGree, J. M. and Müller, W. G. (2020) Bayesian design for minimising uncertainty in spatial processes.  arXiv preprint.  URL: https://arxiv.org/abs/2001.08308.
  * Sharma et al. (2021a) Sharma, A., Dubey, V., Johnson, J., Rawal, Y. and Sivakumar, K. (2021a) Dendritic prioritization through spatial stream network modeling informs targeted management of himalayan riverscapes under brown trout invasion.  Journal of Applied Ecology, 58, 2415–2426.
  * Sharma et al. (2021b) — (2021b) Dendritic prioritization through spatial stream network modeling informs targeted management of himalayan riverscapes under brown trout invasion.  Journal of Applied Ecology, 58, 2415–2426.
  * Shreve (1966) Shreve, R. L. (1966) Statistical law of stream numbers.  The Journal of Geology, 74, 17–37.
  * Shun (1997) Shun, Z. (1997) Another look at the salamander mating data: A modified Laplace approximation approach.  Journal of the American Statistical Association, 92, 341–349.
  * Som et al. (2014) Som, N. A., Monestiez, P., Ver Hoef, J. M., Zimmerman, D. L. and Peterson, E. E. (2014) Spatial sampling on streams: principles for inference on aquatic networks.  Environmetrics, 25, 306–323.
  * Strayer and Dudgeon (2010) Strayer, D. L. and Dudgeon, D. (2010) Freshwater biodiversity conservation: recent progress and future challenges.  Journal of the North American Benthological Society, 29, 344–358.
  * Tickner et al. (2020) Tickner, D., Opperman, J. J., Abell, R., Acreman, M., Arthington, A. H., Bunn, S. E., Cooke, S. J., Dalton, J., Darwall, W., Edwards, G., Harrison, I., Hughes, K., Jones, T., Leclère, D., Lynch, A. J., Leonard, P., McClain, M. E., Muruven, D., Olden, J. D., Ormerod, S. J., Robinson, J., Tharme, R. E., Thieme, M., Tockner, K., Wright, M. and Young, L. (2020) Bending the Curve of Global Freshwater Biodiversity Loss: An Emergency Recovery Plan.  BioScience, 70, 330–342.  URL: https://doi.org/10.1093/biosci/biaa002.
  * Todd et al. (2008) Todd, A. S., Coleman, M. A., Konowal, A. M., May, M. K., Johnson, S., Viera, N. K. and Saunders, J. F. (2008) Development of new water temperature criteria to protect colorado’s fisheries.  Fisheries, 33, 433–443.
  * Ver Hoef et al. (2006) Ver Hoef, J. M., Peterson, E. and Theobald, D. (2006) Spatial statistical models that use flow and stream distance.  Environmental and Ecological Statistics, 13, 449–464.
  * Ver Hoef and Peterson (2010) Ver Hoef, J. M. and Peterson, E. E. (2010) A moving average approach for spatial statistical models of stream networks.  Journal of the American Statistical Association, 105, 6–18.
  * Wagner et al. (2020) Wagner, D., Friedlander, A. M., Pyle, R. L., Brooks, C. M., Gjerde, K. M. and Wilhelm, T. A. (2020) Coral reefs of the high seas: Hidden biodiversity hotspots in need of protection.  Frontiers in Marine Science, 7, 1–13.
  * Wang et al. (2020) Wang, Y., Moe, C. L., Dutta, S., Wadhwa, A., Kanungo, S., Mairinger, W., Zhao, Y., Jiang, Y. and Teunis, P. F. (2020) Designing a typhoid environmental surveillance study: A simulation model for optimum sampling site allocation.  Epidemics, 31, 100391.

##  Appendix A Appendix

###  A.1 Covariance on River Networks

Again, two points $d^{k}$ and $d^{s}$ on a river network are said to be _flow-
connected_ if they share water flow, and _flow-unconnected_ if they reside on
the same network and do not share flow. The tail-up covariance between
locations separated by stream distance $h_{{}_{H}}$ can then be defined as
follows:

$C_{{}_{TU}}(h_{{}_{H}}|\boldsymbol{\theta}_{\textbf{z}_{TU}})=\begin{cases}\textbf{W}_{ks}C(h_{{}_{H}}|\boldsymbol{\theta}_{\textbf{z}_{TU}})&\text{if
flow-connected},\\\ 0&\text{if flow-unconnected},\end{cases}$

where $\textbf{W}_{ks}$ represents the spatial weights between $d^{k}$ and
$d^{s}$ defined by the branching structure of the river network and
$C(h_{{}_{H}}|\boldsymbol{\theta}_{\textbf{z}_{TU}})$ is the (unweighted)
exponential covariance function defined in Equation (2). The tail-down
covariance between locations separated by stream distance $h$ is defined by
the exponential covariance function
$C_{{}_{TD}}(h_{{}_{H}}|\boldsymbol{\theta}_{\textbf{z}_{TD}})=C(h_{{}_{H}}|\boldsymbol{\theta}_{\textbf{z}_{TD}})$
noting that $h_{{}_{H}}$ is the hydrologic distance separating $d^{k}$ and
$d^{s}$. Similarly, the covariance between locations based on Euclidean
distance is
$C_{{}_{EUC}}(h_{{}_{EUC}}|\boldsymbol{\theta}_{\textbf{z}_{EUC}})=C(h_{{}_{EUC}}|\boldsymbol{\theta}_{\textbf{z}_{EUC}})$.
Note that, the exponential model is the only known covariance function that
produces a statistically valid covariance matrix when Euclidean distance is
replaced with total hydrologic distance (Ver Hoef et al., 2006). In matrix
notation, $\textbf{C}_{{}_{TU}}$, $\textbf{C}_{{}_{TD}}$ and
$\textbf{C}_{{}_{EUC}}$ are matrices derived from the above covariance
functions, respectively.

###  A.2 Stochastic Coordinate Exchange Algorithm

In Algorithm 2, an initial random design
$\textbf{d}_{0}=(d_{0}^{1},...,d_{0}^{\gamma})$ is chosen. In the inner-most
for loop, coordinate $j$ in the design d is swapped for a point in the design
space $\mathcal{D}\setminus\textbf{d}$, denoted by
$\textbf{d}^{j\curvearrowright}$. With each swap, the median utility is
computed. Then, the best found utility
$\tilde{U}(\textbf{d}^{j\curvearrowright})$ is compared to the current best
utility $U_{best}$ via the acceptance criteria. If accepted, the best design
is updated. There are $K$ random starts to mitigate convergence to a local
maxima, which can be computed in (embarrassing parallel) fashion.

1 Specify $\tilde{U}(\textbf{d})$ as appropriate. Define prior
$p(\boldsymbol{\theta})$ and likelihood
$p(\textbf{y}|\boldsymbol{\theta},\textbf{d})$. Define design space,
$\mathcal{D}$, as set of points with $\Gamma=|\mathcal{D}|$.

2

3for  _$k=1:K$_ do

4 Initial random design, $\textbf{d}_{0}$

5 Set $\textbf{d}_{k}=\textbf{d}_{0}$ and evaluate
$U_{best}=\tilde{U}(\textbf{d}_{k})$

6 for  _$t=1:T$_ do

7 for  _$j=1:\gamma$_ do

8 for  _$\iota=1:(\Gamma-\gamma)$_ do

9 Set $\textbf{d}_{k}^{j\curvearrowright}$ as $\textbf{d}_{k}$ with the $j$th
coordinate exchanged with the $\iota$th design point not in $\textbf{d}_{k}$

10 Evaluate $\tilde{U}(\textbf{d}_{k}^{j\curvearrowright})$, with large $M$

11 end for

12 Acceptance probability
$p_{W}^{*}=1-W(\mbox{max}(\tilde{U}(\textbf{d}_{k}^{j\curvearrowright})),U_{best})$,
with larger $B$

13 If accepted, $\textbf{d}_{k}=\mbox{arg max
}\tilde{U}(\textbf{d}_{k}^{j\curvearrowright})$ and
$U_{best}=\mbox{max}(\tilde{U}(\textbf{d}_{k}^{j\curvearrowright}))$

14

15 end for

16

17 end for

18

19 end for

Evaluate $\tilde{U}(\textbf{d}_{k})$ with largest $B_{k}$,
$\textbf{d}_{*}=\mbox{arg max }\tilde{U}(\textbf{d}_{k})$

Algorithm 2 Stochastic Coordinate Exchange Algorithm

Figure A.8: Stochastic Coordinate Exchange Algorithm trace plots using ACE acceptance criteria (left) and the proposed Wilcoxon acceptance criteria (right) for the expected KL divergence utility. Table 1: Optimal design scenarios for the river network case study Acceptance criteria | Utility measure | Optimal design
---|---|---
Wilcoxon | Mean | (167, 169, 172, 174, 183)
ACE | Mean | (167, 169, 172, 174, 183)
Wilcoxon | Median | (167, 169, 172, 174, 183)
ACE | Median | (167, 169, 172, 174, 183)

In our example, the Wilcoxon and ACE acceptance criteria resulted in the same
optimal design, with respect to both the mean and median utility CE algorithm
(Table 1). It is noted that with more uncertain priors or multi-modal
utilities, the comparative results between acceptance criteria may be more
pronounced. Both optimal designs, and indeed most of the best designs for each
random start, contain the points at confluence (167, 169, 174).

###  A.3 Median Utility: MCMC vs. Laplace

In general, the form of the posterior is unknown, and therefore difficult to
estimate directly. In such cases, algorithms like Markov Chain Monte Carlo
(MCMC) (e.g. Metropolis–Hastings algorithm) are employed to draw samples which
are used to approximate the posterior distribution (Hastings, 1970).

To efficiently approximate the expected utility in Bayesian design, fast
methods for approximating the posterior distribution are required. For this,
we consider the Laplace approximation to the posterior distribution, see Long
et al. (2013); Overstall et al. (2018); Carlon et al. (2020); Senarathne et
al. (2020). There are some advantages in using the Laplace approximate method
compared to numerical Markov chain Monte Carlo (MCMC) methods. Besides the
huge gains in computational efficiency for low-dimensional design problems,
there is no issue of choosing a suitable proposal density, for which to tune
for a desired acceptance rate, and similarly no need to ensure convergence of
a Markov chain, as well as no need to determine the length of a suitable
“burn-in” phase. However, the Laplace approximation is limited to the
assumption of an approximately normal target distribution, has been shown to
underestimate posterior variance (Shun, 1997).

Therefore, we compared the median utility, $\tilde{U}(\textbf{d})$, using the
Laplace approximation vs. an MCMC (Metropolis–Hastings algorithm)
approximation to the posterior for 50 random designs. Figure A.9 shows the
linear relationship between the Laplace and MCMC inference methods for the KL-
divergence utility in the design problem outlined in Section 5.1.

Figure A.9: The median utility, $\tilde{U}(\textbf{d})$, computed for 50
random designs using inference methods: MCMC (Metropolis–Hastings algorithm)
vs. Laplace approximation. The optimal design found via the proposed
stochastic Coordinate Exchange algorithm (Algorithm 2) is shown in blue.