Being Bayesian in the 2020s: opportunities and challenges in the practice of modern applied Bayesian statistics

The ability to determine and characterise underlying mechanisms in complex systems is paramount to pioneering research and scientific advancement in the modern era. Over the last decade, the rise of data generation from sensor and internet enabled devices has catalysed the advancement of data collection technologies and analysis methods used to extract meaningful information from complex systems. However, the sheer size of these complex systems (e.g., natural ecosystems like the Great Barrier Reef and river networks) and the expense of data collection means that data cannot be collected throughout the whole system. Further, practical constraints like connectivity, accessibility, and data storage issues reduce our ability to sample frequently through time. This has led to innovation in statistical methods for data collection, promoting an emerging era of “intelligent data collection” where data are collected for a particular purpose such as understanding mechanisms for change, monitoring biodiversity, and identifying threats or vulnerabilities to long-term sustainability. Bayesian optimal experimental design is one such area of recent innovation.

Bayesian design offers a framework for optimising the collection of data specified by a design $\boldsymbol{d}$ for a particular experimental goal, which may be to increase precision of parameter estimates, maximise prediction accuracy and/or distinguish between competing models. More specifically, Bayesian design is concerned with maximising an expected utility, $U(\boldsymbol{d})=\mathbb{E}\,u(\boldsymbol{d},\boldsymbol{\theta},\boldsymbol{y})$ through the choice of design $\boldsymbol{d}$ within a design space $\mathcal{D}$, while accounting for uncertainty about, for example, the parameter $\boldsymbol{\theta}\in\Theta$ and all conceivable data sets we might observe $\boldsymbol{y}\in\mathcal{Y}$. A Bayesian optimal design $\boldsymbol{d}^{*}$ can therefore be expressed as,

where $p(\boldsymbol{\theta},\boldsymbol{y};\boldsymbol{d})$ defines the joint distribution of $\boldsymbol{\theta}$ and $\boldsymbol{y}$ given a design $\boldsymbol{d}$.

Unfortunately, determining $\boldsymbol{d}^{*}$ can be challenging. Firstly, the utility function $u(\boldsymbol{d},\boldsymbol{\theta},\boldsymbol{y})$ typically involves computing some form of expected value with respect to the posterior distribution, which itself is typically analytically intractable. Further, $U(\boldsymbol{d})$ is itself an expectation taken with respect to the prior-predictive distribution, which is also typically intractable. This means numerical or approximate methods are needed, which may impose substantial compute time and/or require a stochastic approximation. For example, Monte Carlo integration has been proposed as an approach to form an approximation to the expected utility as follows:

where $\boldsymbol{\theta}^{(m)}\sim p(\boldsymbol{\theta};\boldsymbol{d})$ and $\boldsymbol{y}^{(m)}\sim p(\boldsymbol{y}|\boldsymbol{\theta}^{(m)};\boldsymbol{d})$, for some large value of $M$. Thus, computations involving $M$ different individual posterior distributions are required just to approximate the expected utility of a design.

Secondly, $\boldsymbol{d}$ may be high-dimensional, meaning that a potentially large optimisation problem needs to be solved for a computationally expensive and noisy objective function. Accordingly, the majority of research in Bayesian design has focused on developing new methods to address one of these two challenges. Below we provide a brief summary of some of the relevant literature from the last 25 years.

Since the conception of statistical decision theory [1, 2] upon which the decision-theoretic framework of Bayesian design is based [3], there have been numerous strategies presented in the literature to address the above challenges. Curve-fitting methods were proposed by [4] to the approximate expected utility in Equation (1). Here, the fitted curve is optimised as a surrogate for the true expected utility to determine the choice of (approximately optimal) design. An alternative simulation-based method proposed by [5] formed the following augmented joint distribution on $\boldsymbol{d}$, $\boldsymbol{\theta}$, and $\boldsymbol{y}$:

where it can be shown that the marginal distribution of $\boldsymbol{d}$ is proportional to $U(\boldsymbol{d})$. Markov chain Monte Carlo (MCMC) methods were then used to sample from this distribution, and subsequently to approximate the marginal mode of $\boldsymbol{d}$. Extensions of this approach were given by [6, 7] which include adopting a sequential Monte Carlo algorithm to more efficiently sample from the augmented distribution as $J$ increases. However, such approaches are limited to low-dimensional design problems (i.e.,  3–4 design points) and simple models due to difficulties in sampling efficiently in high dimensions.

Recently, there has been a shift from sampling-based methods to rapid, approximate posterior inference methods. Combined with a Monte Carlo approximation as given in Equation (1), this has enabled expected utility functions to be efficiently approximated for realistic design problems. This includes those based on complex models (such as non-linear models) and models for data that exhibit complex dependence structures (such as those with different sources of variability including spatially and between groups). Such approximate inference methods include the Laplace approximation [8] and variational Bayes [9], which have been combined with new optimisation algorithms (e.g., the approximate coordinate exchange algorithm; ACE [10]) to solve the most complex and high-dimensional design problems to date.

The most prominent application of Bayesian design methods appears in the clinical trial literature [11]. Recently, this has been exacerbated by the outbreak of COVID-19 where it has been desirable to conduct clinical trial assessments as quickly as possible, with Bayesian (adaptive) designs shown to yield more resource efficient and ethical clinical trials [12, 13]. More recently, Bayesian design methods have been proposed as a basis to efficiently monitor large environmental systems like the Great Barrier Reef [14, 15]. In the following case study, we show how such methods can be used to form sampling designs to monitor a coral reef system, and extend these methods to provide flexible designs that address major practical constraints when sampling real-world ecosystems.

Coral reefs are biodiversity hot-spots for marine species under threat from anthropogenic impacts related to climate change, water pollution and over-exploitation, among other factors [16]. Coral cover is a commonly used indicator to infer the health of coral reef environments [17], where data collection relies on a series of images taken underwater, along a transect (a line across a habitat). Monitoring of coral reef environments is expensive in terms of monetary, human and technological costs, particularly for remote locations. Informative data is critical to support conservation decisions, but with limited resources to invest in monitoring programs, the need is to optimise in-field activities that will result in the intelligent collection of data. Following [18], we consider monitoring submerged shoals which are coral reefs that exist at depths of around 18 to 40m below sea level. Data collection at such depths requires unmanned vehicles to be deployed along a design, i.e., a series of transects which specify where images should be collected. However, spatially precise sampling is known to be difficult in deeper reefs due to unpredictable weather and water currents. Therefore, our aim is to provide Bayesian designs that offer flexibility in where transects will be placed while taking into consideration the complex nature of the systems we are monitoring such as the spatial dependence of natural processes.

In order to define a design, we specify the placement of each transect $k=1,\ldots,q$ on the shoal by its midpoint given in Easting and Northing coordinates, i.e., $E_{k}$ and $N_{k}$, and the angle of the transect, $\alpha_{k}$, in degrees. Each transect line is expressed as a design point $\boldsymbol{d}_{k}=(E_{k},N_{k},\alpha_{k})$. The exact sampling locations (equally spaced along the fixed-length transect) are specified as $\boldsymbol{s}_{i}$. For each transect, we introduce a radius parameter $r_{k}>0$ for the purpose of allowing the sampled image locations to disperse by $\delta_{1},\delta_{2}\stackrel{{\scriptstyle\rm iid}}{{\sim}}\mbox{Unif}(-r_{k},r_{k})$, i.e., sampling at $\boldsymbol{s}+\boldsymbol{\delta}$. For image $i$, a number, $n_{i}$, of randomly selected points on the image are classified as either hard coral or not. Accordingly, the number of points within an image that contain hard coral, $y_{i}$, is modelled as

for regression parameters $\boldsymbol{\beta}=(\beta_{0},...,\beta_{n})^{\top}$, covariance kernel parameters $\boldsymbol{\gamma}=(\gamma_{1},...,\gamma_{m})^{\top}$ for spatially-correlated random effect Z, and covariates $\boldsymbol{x}_{i}$. The priors for $\beta$ and $\gamma$ are based on consideration of historical data (depth and depth squared) collected on the shoal. See [18] for further details.

As a basis for improved monitoring, we consider the amount learned from the data regarding parameters of the above model as our goal of data collection. For this, we specify our utility function as the Kullback-Liebler divergence of the posterior from the prior distribution, where larger values suggest the data is more informative with respect to model parameters.

For a computationally-efficient approximation of the utility for a given design, we employ a Laplace approximation of the posterior distribution, i.e., an approximation of the form

where $\boldsymbol{\theta}^{*}=\mbox{arg max}_{\boldsymbol{\theta}\in\Theta}\mbox{log }p(\boldsymbol{y},\boldsymbol{\theta};\boldsymbol{d})$ and $\mathrm{H}(\boldsymbol{\theta}^{*})$ is the Hessian matrix evaluated at $\boldsymbol{\theta}^{*}$. Here $\boldsymbol{\theta}=(\boldsymbol{\beta},\boldsymbol{\gamma})$, and marginalization of $\boldsymbol{Z}$ is performed approximately using Monte Carlo integration. To obtain an optimal design to for monitoring of the shoal, we propose a two-step approach:

1. 1.

   Firstly, a global search for the Bayesian optimal design $\boldsymbol{d}^{*}=(\boldsymbol{d}_{1}^{*},...,\boldsymbol{d}_{q}^{*})$, where $q=3$ (the total number of transects) is conducted. We consider a discretised design space, and find designs via a discrete version of ACE; and

2. 2.

   Secondly, we form design efficiency windows (illustrating robustness to imprecise sampling) across $r_{k}$ for each transect $k=1,\ldots,q$. To do so, we specify a zero-mean Gaussian process (GP) prior for the approximate expected utility across $\boldsymbol{r}\in\mathbb{R}^{q}$ by $\widehat{U}(\boldsymbol{r};\boldsymbol{d}^{*})$, i.e.,  $\widehat{U}(\boldsymbol{r};\boldsymbol{d}^{*})\sim\mathcal{GP}(\textbf{0},\mathrm{K}(\cdot)+\zeta_{0}\mathrm{I})$, for some kernel matrix $\mathrm{K}(\cdot)$, and $\zeta_{0}>0$. The windows are then obtained as follows:

   1. (a)

      Center the radius on $\boldsymbol{d}^{*}$, i.e., the Bayesian design from (i), and specify a maximum value for $r_{k}$ for $k=1,\ldots,q$;

   2. (b)

      Randomly sample $\delta_{i,1},\delta_{i,2}\stackrel{{\scriptstyle\rm iid}}{{\sim}}\mbox{Unif}(-r_{k},r_{k})$, where $k$ is the transect from which image $i$ is obtained, and evaluate the approximate expected utility of the design at locations $\boldsymbol{s}_{i}+\boldsymbol{\delta}_{i}$;

   3. (c)

      Fit a Gaussian process defined on $\boldsymbol{r}\in\mathbb{R}^{q}$ to the approximate expected utilities;

   4. (d)

      Emulate the expected utility surface across values of $\boldsymbol{r}$ using the posterior predictive mean of the Gaussian process, denoted $\bar{U}(\boldsymbol{r})$;

   5. (e)

      Normalise the predicted expected utility values by that of the original Bayesian design as follows:

      $${\rm eff}(\boldsymbol{r})=\frac{\bar{U}(\boldsymbol{r};\boldsymbol{d}^{*})}{\bar{U}(\boldsymbol{0};\boldsymbol{d}^{*})},$$

      (2)

      and use the above to obtain design efficiency contours (plotted in Figure 1, left). For some design efficiency contour value $c>0$, the corresponding sampling window is the region in space defined by radii $\boldsymbol{r}(c)$ that satisfy ${\rm eff}(\boldsymbol{r}(c))=c$.

Based on the approach, the Bayesian design, $\boldsymbol{d}^{*}$, shown in Figure 1, situates the transects in shallower areas of the reef, but at different depths in the shallow areas, presumably to provide information about the depth effects, $\boldsymbol{\beta}$. Avoiding the deeper regions of the shoal makes sense physiologically, as the corals monitored here are photosynthetic organisms, and therefore rely on light to survive. This design thus avoids the collection of data in areas where there is little chance of observing coral. The design efficiency contours are also shown in Figure 1. If we consider a design efficiency of 0.99 in Equation (2), then possible radius values are $(50,0,44)$ for the three transects, with the flexibility (sampling windows) this provides shown around each transect. As can be seen, transect $\boldsymbol{d}^{*}_{2}$ is more sensitive than the other two, suggesting more effort should be placed in sampling this transect precisely. In practical terms, sampling from shallow areas of the reef, $\boldsymbol{d}^{*}_{1}$ and $\boldsymbol{d}^{*}_{3}$, can be undertaken when the conditions are more unpredictable (e.g., strong currents), and samples from $\boldsymbol{d}^{*}_{2}$ can be obtained when field conditions are more preferable.

In conclusion, Bayesian optimal design addresses a fundamental problem in science: the intelligent collection of data resulting in greater information efficiency, reduced sampling cost and improved estimation. Such benefits have been observed in clinical trials [19, 20, 21, 22] and environmental monitoring [23, 15], and we have shown how they can be used to offer flexible yet efficient sampling in a real-world context. One limitation of the approach is the potential reliance of designs on a number of assumptions, e.g., an assumed model for the data, so we would encourage future research in areas that reduce this reliance and thus provide more robust designs for data collection.

Recent advancements in digital technologies have led to the large-scale collection of more advanced data such as sensor data, satellite imagery, and a host of varied resolution imagery and video including those taken using 360-degree cameras. It is possible, using these advances in technology, to enable location-specific data to remote researchers for analysis. The emergence of virtual reality (VR) technology, for example, acts as a way to connect the public and the scientific community, creating innovative pathways for environmental conservation research by immersing subjects in an otherwise inaccessible vivid virtual scene for elicitation purposes [24, 25]. The opinions and knowledge extracted from this process is itself new data, which can be used for educational purposes [26] or incorporated into statistical models.

Increases in the volume of more complex types of data has led to the development of more effective and efficient analysis methodology. Recently, Bayesian models have seen use as a method to evaluate subject elicitation in the areas of coral reef conservation [27], jaguar and koala habitat suitability assessments [28, 29], and the aesthetic value of sites in the Antarctic Peninsula.

Citizen science (CS) represents one of the most popular emerging data sources in scientific research. CS involves engaging members of the general population in one or more parts of the scientific process. Its applications can be found across almost all disciplines of science, especially in ecology and conservation where scientists are harnessing its power to help solve critical challenges such as climate change and the decline in species abundance. Examples of citizen scientists’ contributions include reporting sightings of species, measuring environmental variables and identifying species on images. Hundreds of citizen science projects can be found in popular online platforms including Zooniverse [33], eButterfly [34], eBird [35], and iNaturalist [36]. A fundamental issue often discussed surrounding citizen science is the quality of the data produced, which is generally error-prone and biased. For example, bias can arise in citizen science datasets due to (1) the unstructured nature of the data, (2) collecting data opportunistically, with more observations from frequently visited locations [37] or at irregular frequencies across time [38], and (3) as a result of differing abilities of the participants to perform tasks such as detecting or identifying species [39, 40]. However, recent advances in statistics, machine learning, and data science are helping realize its full potential and increase trustworthiness [41, 42, 40].

Frequently, citizen science data is elicited via image classification. For example, asking the participants whether images contain a target class or species. In this section, we illustrate two modelling approaches for these types of data.

In the first approach, we consider a binary response variable $y_{ij}$ representing whether the category has been correctly identified by the participant ($i=1,\ldots,m$) in the image ($j=1,\ldots,n$). The probability of obtaining a correct answer can be modeled using an item response model such as the three-parameter logistic model (3PL) [43, 40],

where each $\eta_{j}\in(0,1)$ is a pseudo-guessing parameter accounting for a participants’ chance of answering correctly by guessing, $Z_{i}$ is the latent ability of the $i^{\textrm{th}}$ participant, $\alpha_{j}>0$ is the slope parameter, and $B_{j}$ is the latent difficulty of the $j^{\textrm{th}}$ image. Sometimes, the correct answer for certain images is unknown. In this case, we estimate the latent labels for the images via the estimates of $Z_{i}$, by using the latter as weights in popular methods such as majority or consensus voting. Code to fit these models and exemplar datasets can be found in [44].

The second approach is for the case where we are interested in the proportion of species in elicitation points in images. Here, we compute a statistic $\hat{y}_{ij}\in[0,1]$ giving the apparent proportion of species in a number of elicitation points in image $j$ classified by the participant $i$. The true latent proportion $Y_{j}$ can be estimated based on each participant’s overall performance measures ${\rm se}_{i}$ and ${\rm sp}_{i}$ (which denote the sensitivity and specificity scores of participant $i$, respectively). A Beta prior is placed on the true proportion, yielding the model

where $\alpha_{j}$ and $\beta_{j}$ are the shape and the scale parameters in the beta distribution, respectively. The above model can be parametrized via a specified prior mean $\mu_{j}$ for each $Y_{j}$, and a common precision parameter $\phi$, via ${\alpha_{j}=\mu_{j}\phi}$ and ${\beta_{j}=-\mu_{j}\phi+\phi}$, which in turn implies that ${\mathrm{Var}[Y_{j}]=\frac{\mu_{j}(1-\mu_{j})}{(1+\phi)}}$. Covariates can also be incorporated by defining a beta regression with $\textrm{logit}(\mu_{j})=\boldsymbol{\xi}^{\top}\boldsymbol{x}_{j}+U_{j}+\varepsilon_{j}$, where $\varepsilon_{j}$ are error terms, and $U_{j}$ are spatially-dependent random effects. Both approaches account for spatial variation (captured in $B_{j}$ or $U_{j}$ for the first and second approach, respectively) using different spatial structures (e.g., conditional autoregressive (CAR) priors, covariance matrices, or Gaussian random fields). See more details in [42, 45, 40].

The following case study illustrates the estimation of the latent proportion of hard corals across the Great Barrier Reef in Australia, obtained from underwater images classified by citizen scientists. Figure 2 shows 15 spatially balanced random points in one of the images used in the study. The apparent proportion of hard coral in the image was obtained using the number of points selected by participants containing this category out of 15. Using Equation (4), the (biased) estimates obtained from the citizen scientists can be corrected producing a similar density to the latent unobserved proportions.

The integration of citizen science data with current monitoring efforts from Australian federal agencies and non-governmental organisations is a breakthrough to increase the amount of information about changes along the Great Barrier Reef, learn about climate change impacts and adapt management actions consequently. This model introduced here is the root of a digital platform that estimates the health of the Great Barrier Reef using all available information. Australian Cancer Atlas (atlas.cancer.org.au) is a citizen science program that involves non-experts to collect information online and models it appropriately via a Bayesian approach. This study contributes to increasing the trust of citizen science and produce reliable data for environmental conservation while engaging and arising awareness about coral reefs.

In many areas of study (health, business, and environmental science, for example), the collective data set of interest one wishes to use in modelling is often under the control of different data custodians, i.e., parties responsible for ensuring data is only used or released in instances deemed appropriate to governance requirements. Such requirements often stipulate that the data itself, and information pertaining to it can only be shared in a manner deemed sufficiently private.

Federated learning is the process of fitting a model in the setting where data resides with multiple custodians. Approaches typically place privacy as having the utmost importance, but computational efficiency is important from a practical perspective. There are two broad data settings that occur, each requiring their own style of algorithm. The first is the horizontal setting. Here, multiple data custodians have the same set of variables for different entities. In contrast, the vertical federated learning setting has data custodians who possess different variables for the same entities. An example of a horizontal setting is where different countries possess the data for those primarily residing within. In contrast, an example of vertical federated learning would be where two companies possess their respective sales data for the same collection of customers. The term “federated learning” originated in the deep learning literature with the introduction of the FedAvg algorithm [46]. FedAvg involves updating parameter values of a global model to be the weighted average of parameter values obtained by updating the same model locally (possibly many times) at each iteration. This work led to many related optimization algorithms, e.g., FedProx[47], and FedNova [48] which account for heterogeneous (non i.i.d.) data sources, and the Bayesian nonparametric approach for learning neural networks of [49], where local model parameters are matched to a global model via the posterior distribution of a Beta-Bernoulli process [50]. To date, practical federated analyses appear restricted to the frequentist setting. Examples include the prediction of breast cancer using distributed logistic regression [51] and modelling of the survival of oral cavity cancer through a distributed proportional Cox hazards model [52]. Both these approaches conduct parameter estimation via a Newton-Raphson algorithm [53] and result in equivalent maximum likelihood estimates to those obtained in a standard, non-federated setting. Algorithms for the maximum likelihood estimation of log-linear and logistic regression models in vertical federated learning settings [54, 55, 56, 57, 58] use ideas such as secure multiparty computation [59], and formulating the parameter estimation task as a dual problem [60]. Several overarching software infrastructures such as VANTAGE6 [61] ensure the correct and secure use of the data of each custodian within the specified algorithm, given acceptable (model- and application-specific) rules for information exchange.

Despite the potentially enabling capabilities of federated methods, to our knowledge, Bayesian federated learning methods have yet to impact real-world applications. In the Bayesian inference setting, the “learning” task becomes one of performing posterior inference, e.g., via Markov chain Monte Carlo (MCMC) or variational inference (VI) techniques. Note that Bayesian federated learning approaches may involve multiple communication rounds, though this is only sometimes the case. For example, many distributed MCMC approaches (e.g., [62, 63, 64]), combine individually-fit model posteriors, requiring only a single communication step from each local node. A recent intermediate approach [65] is to construct a surrogate likelihood of the complete data set via an arbitrarily-specified number of communication steps. After constructing the surrogate likelihood, an MCMC algorithm is run on a single device. As the number of communication steps increases, the approximation error introduced by the surrogate likelihood decreases.

In certain cases, carrying out federated Bayesian inference is (at least in principle) relatively straightforward. For example, a naive MCMC algorithm would be trivial to construct for a simple model class, such as any generalized linear model (which assumes the data are independent), provided that one is not concerned with the number of communication steps. To see this, note that the (log-)posterior density function decomposes as

Hence, for the horizontal setting, all that is required is the nodes sharing the sum of their respective log-likelihood terms with the server. However, this approach would require a minimum of two communication steps per iteration of the Markov chain. Recent MCMC methods, similar in style to the FedAvg algorithm (which use Langevin dynamics to update the Markov chain), require only a single communication step per iteration [66, 67]. Such approaches exploit gradient information which decomposes as a sum similarly to (5), though eschew the usual Metropolis–Hastings correction and are hence asymptotically inexact. In some instances, a formally-justified notion of privacy may be required, as opposed to simply an intuitive one given by aggregation of terms. Differential Privacy (DP) (e.g.,[68]) provides such guarantees, and there are variants of MCMC that ensure this, such as DP-MCMC [69], which accomplishes privacy guarantees at the cost of a slight perturbation of stationary distribution of the chain. It is worth noting that all of the above examples mentioned are specific to the horizontal setting, with the vertical setting proving especially challenging as one does not have a beneficial decomposition like that of (5).

As the above alludes to, the development and use of Bayesian federated learning algorithms are complex for several reasons. A method is only suitable for a prescribed application if it satisfies a combination of requirements, such as being able to work with the desired model, computational and communication costs, privacy, and accuracy. For each application, the choice of model and federated method will depend on where the priorities lie, e.g., accuracy, efficiency, or privacy. In some cases, there may be no feasible algorithm (an example is given in the upcoming case study). Thus, inference approaches that improve upon some (or even all) of these aspects are important and warrant future research.

The ultimate goal of federated Bayesian analysis is to circumvent the need for data merging [70, 71] in scenarios where merging is considered infeasible. However, for Bayesian federated learning to reach this point, these approaches must offer custodians and interested parties an accurate inference for complex models while maintaining a level of privacy acceptable to those data custodians. Thus, the methodological development that enables federated inference for more advanced Bayesian models efficiently and/or with additional privacy guarantees is likely to emerge as a critical area of interest in the coming years.

The greatest hindrance to employing federated learning in real-world applications is the lack of possible model types that current algorithms address. Commonly, applied statistical modelling involves incorporating hierarchical structures, and latent variables [72]. To our knowledge, there are no federated Bayesian analysis algorithms at all for such models. To briefly illustrate the unique challenges and the need for developments that account for the nuances of different models, the case study considers spatially-dependent latent variables based on neighbourhood structures. For simplicity, the focus is on the Intrinsic Conditional AutoRegressive (ICAR) prior [73], although variations such as the Besag–York–Mollie [74] and Leroux [75] models are similar in what follows (the latter is used for example, in the Australian Cancer Atlas described in Section 7). The ICAR prior posits a vector of spatially-dependent latent variables, denoted here as $\boldsymbol{Z}$. Each element of $\boldsymbol{Z}$ corresponds to a latent area-level effect of a “site”, which is influenced by neighbouring sites. Writing $i\sim j$ to denote that sites $i$ and $j$ are considered neighbours, and assuming the graph arising from the neighbourhood structure is fully connected, the ICAR prior with precision hyperparameter $\tau$ has log-density

The above may be problematic if the data custodians insist that the latent variables corresponding to their areas must be kept private to themselves. To see why, consider the case that there are two (2) data custodians, with the sets ${\mathcal{C}}_{1}$ and ${\mathcal{C}}_{2}$ containing the indices of data possessed by the first and second custodian, respectively. Then,

where terms in blue are those relevant to the sites under the first custodian, and those in red to the second. When computing the log-posterior density (as required, for example, in Markov chain Monte Carlo algorithms), the first two terms on the right-hand side above can be aggregated and sent to the central server. However, the final term can not as each individual summand requires the individual latent variables to be processed. This is because the latter term considers interactions across custodian boundaries.

Consequently, solutions such as (i) employing judicious reparameterization of the latent variables (possibly compatible with the one that is often also required to enforce identifiability), (ii) changing the model to add additional auxiliary variables, or (iii) otherwise approximating the troublesome term, are required. An additional challenge is that even if inferences on individual latent variables are only available to their respective custodians, they may nevertheless “leak” information across custodian boundaries to neighbouring sites due to the underlying dependency structure.

While the above certainly highlights particular challenges, the first two terms of (6) split nicely across custodians and hint that latent variables need not always be problematic. For more straightforward cases such as the latter, specialized accurate and efficient inference approaches that allow individual custodians to avoid ever sharing their latent variables (either directly or indirectly) or data are the subject of forthcoming work by the authors of this section, who have a longer-term goal of tackling more challenging cases such as ICAR and its relatives in different settings.

The appetite for developing more realistic data-driven models of complex systems is continually rising. Development of such complex models can improve our understanding of the underlying mechanisms driving real phenomena and produce more accurate predictions. However, the calibration of such models remains challenging, as the associated likelihood function with complex models is often too computationally cumbersome to permit timely statistical inferences. Despite this, it is often the case that simulating the model is orders of magnitude faster than evaluating the model’s likelihood function. Such models with intractable likelihoods that nevertheless remain simulable are often referred to as implicit models. Such models are now prevalent across many areas of science (see e.g., various application chapters in [76]).

Currently, the most popular statistical approach amongst practitioners for performing Bayesian inference for implicit models is approximate Bayesian computation (ABC), popularised by [77]. A related method called generalised likelihood uncertainty estimation [78, 79] predates ABC, and [80] explore its connections to ABC. The ABC approach approximates the true posterior as

Here, x denotes simulated data that has the same structure as y, $||\cdot||$ is some norm (i.e., $||\textbf{x}-\textbf{y}||$ measures the closeness of the simulated data to the observed data), and $\epsilon$ stipulates what is considered “close”. Intuitively, values of $\boldsymbol{\theta}$ more likely to produce simulated data x close enough to y have increased (approximate) posterior density. Rather than compare y and x directly, it can be more efficient to compare y and x in a lower dimensional space via a summarisation function that aims to retain as much information from the full data set as possible. For the posterior in (7) to equate to the exact posterior, we require that the observed and simulated datasets are matched perfectly (i.e., as $\epsilon\rightarrow 0$) in terms of some sufficient summarisation. However, in the majority of practical applications, a low-dimensional sufficient statistic does not exist and it is computationally infeasible to take $\epsilon\rightarrow 0$, so we must accept some level of approximation.

Given the wide applicability of the approach, i.e., that only the ability to simulate the model is required to conduct inference, there has been an explosion of research in the past 10–15 years advancing ABC and related methods that lie within the more general class of so-called likelihood-free inference methods. A substantial portion of methodologically-focused ABC research considers aspects including the effective choice of $||\cdot||$ (e.g., [81, 82]), efficient sampling algorithms to explore the approximate posterior in (7) (e.g., [83]) and ABC’s theoretical properties (e.g. [84]). Many of the developments of ABC and some related methods (e.g. Bayesian synthetic likelihood [85, 86]) prior to 2018 are discussed in [76], the first-ever monograph on ABC.

The following case study considers a popular class of sampling algorithms for ABC based on sequential Monte Carlo (SMC). SMC-based ABC algorithms improve efficiency compared to sampling naively from the prior by gradually reducing the ABC tolerance $\epsilon$ where the output produced at iteration $t$ is used to improve the proposal distribution of $\boldsymbol{\theta}$ at iteration $t+1$. The output of the algorithm is $N$ samples, or “particles”, from the ABC posterior in (7) with a final $\epsilon$ that is either pre-specified or determined adaptively by the algorithm. Each particle has attached to it a “distance”, which is the value of $||\textbf{x}-\textbf{y}||$ for x simulated from the model based on the particle’s parameter value. Here, we use the adaptive SMC-ABC algorithm in [87], itself a minor modification of the replenishment algorithm of [88]. The algorithm is summarised below.

1. 1.

   Draw $N$ samples from the prior, and for each sample, simulate the model and compute the corresponding distance. Initialise $\epsilon$ as the largest distance amongst the set of particles.

2. 2.

   Set the next $\epsilon$ as the $\alpha$-quantile of the set of distances. Retain the $N_{\alpha}$ particles with distance less than or equal to $\epsilon$.

3. 3.

   Resample the retained particle set $N-N_{\alpha}$ times so that there are $N$ particles.

4. 4.

   Run MCMC on each of the resampled $N-N_{\alpha}$ particles with stationary distribution (7) with the current $\epsilon$. This step helps to remove duplicate particles created from the previous resampling step. The number of MCMC iterations can be adaptively set based on the MCMC acceptance rate.

5. 5.

   Repeat steps (ii)-(iv) until a desired $\epsilon$ is reached or the MCMC acceptance rate in (iv) is too small (i.e., the number of MCMC steps becomes too large for the computational budget).

A key computational inefficiency of ABC and closely-related methods such as BSL is that many of the model simulations yield MCMC proposals that are rejected. To obtain a suitable quality of approximation, it is not uncommon to require continuing the algorithm past the point where $\epsilon$ is small enough to have average acceptance probabilities of $10^{-2}$ or less. To overcome this issue, there has been significant attention devoted to machine learning based approaches to likelihood-free inference, especially in the past five years. These methods use model simulations (from different parameters) as training data for building a conditional density estimator of the likelihood (e.g., [89]), likelihood ratio (e.g., [90]) or posterior density (e.g., [91]). Following this estimation, standard methods from the Bayesian inference toolkit can be used. Many machine learning approaches to likelihood-free inference can be implemented sequentially, so that samples from the approximate posterior in the previous (or all previous) iterations can comprise an increasingly informed training set that yields a more accurate conditional density estimator in regions of non-negligible posterior probability. For the case study below, we compare the SMC-ABC approach with the sequential neural likelihood (SNL) method of [89], which is outlined below.

1. 1.

   Set the initial proposal distribution of parameter values as the prior, i.e., $q(\boldsymbol{\theta})=p(\boldsymbol{\theta})$.

2. 2.

   Generate a training data set by drawing $M$ parameter/simulated data pairs according to $q(\boldsymbol{\theta})p(\mathbf{x}|\boldsymbol{\theta})$. Fit a conditional normalising flow (a flexible type of regression-density estimator) to the training data to estimate the conditional density of $\mathbf{X}|\boldsymbol{\theta}$.

3. 3.

   Run MCMC to obtain approximate posterior samples, using the learned conditional density of $\mathbf{X}|\boldsymbol{\theta}$ evaluated at the observed data $\mathbf{y}$ as the approximation to the likelihood. Samples from this approximate posterior can also be used to update the proposal distribution $q(\boldsymbol{\theta})$.

4. 4.

   Repeat steps (ii) and (iii) for a desired number of rounds.

In this case study, we apply likelihood-free methods SMC-ABC and SNL for calibrating a complex agent-based model (ABM) of tumour growth. We briefly compare the methods in terms of computationally efficiency and their ability to fit simulated and real tumour growth data.

ABMs have been used in cancer modelling for some time now as they provide a spatial representation of the inherent cellular heterogeneity and stochasticity of tumours [92, 93, 94]. Largely, these models account for the individual cell based behaviours of proliferation, movement and death and aim to predict the impact of stochasticity on spatial tumour growth over time. Previous works have considered this in the context of angiogenesis [95], immune involvement [96, 97] and also treatment [98]. In some cases, data has been used to calibrate or validate aspects of the models [99, 100], although due to the computational cost and their intractable likelihood it is not always easy to infer parameters in an ABM using data.

For this case study, we use a previously published ABM called a Voronoi Cell-Based model (VCBM) [101, 102]. In this model, cancer cells and healthy tissue cells are considered agents, whose centre is modelled by a point on a 2D lattice, and whose boundary is defined by a Voronoi tessellation. To mimic tumour growth and spatial tissue deformation, the model captures cell movement using force-balance equations derived from Hooke’s law. In this way, cell movement is captured off-lattice and is a function of the local cell-neighbourhood pressure, determined using a Delaunay Triangulation.

Tumour growth is captured by introducing a probability of an individual cancer cell proliferating $P$, which is a function of a cell’s distance to the boundary of the tumour: $P=p_{0}\left(1-\frac{d}{d_{\mathrm{max}}}\right)$, where $p_{0}$ is the probability of proliferation, $d$ is the cell’s Euclidean distance to the tumour boundary (measured from the cell centre to the nearest healthy cell centre) and $d_{\mathrm{max}}$ is the maximum radial distance a cell can be from the boundary and still proliferate. In this way, the model evolves stochastically over time with cells either proliferating or moving in a given timestep. The model also uses $g_{\mathrm{age}}$ to define the time taken for a cell to be able to proliferate and uses $p_{\mathrm{psc}}$ as the probability of cancer cell invasion. Hence, the model parameter $\boldsymbol{\theta}$ to be estimated is $(p_{0},p_{\mathrm{psc}},d_{\mathrm{max}},g_{\mathrm{age}})$.

To validate the VCBM, we use published in vivo tumour growth measurements for ovarian cancer [103]. In these experiments, tumour volume was recorded by measuring the tumour width and length as perpendicular axis using calipers and then calculating the approximate tumour volume. We simulate the VCBM in 2D and calculate the corresponding tumour volume measurements equivalently. We also consider one simulated dataset generated with parameter value $\boldsymbol{\theta}=(0.2,{10}^{-5},31,114)$. The datasets are shown as solid black lines in Figure 5. The prior distribution on $\boldsymbol{\theta}$ is given by $p_{0}\sim\mathrm{Beta}(1,1)$, $p_{\mathrm{psc}}\sim\mathrm{Beta}(1,{10}^{4})$, $d_{\mathrm{max}}\sim\mathrm{LogNormal}(\mathrm{log}(30),1)$ and $g_{\mathrm{age}}\sim\mathrm{LogNormal}(\mathrm{log}(160),1)$, with parameters assumed independent a priori [104].

We run SMC-ABC until around 100,000 model simulations have been generated for each dataset. We use the SBI package [105] to implement SNL with five rounds of 10,000 model simulations for each dataset. To compare the performance of SMC-ABC and SNL, we compute the posterior predictive distribution for each dataset. For SNL, we choose the round that visually produces the most accurate posterior predictive distribution. We find that for SNL the performance can degrade with increasing rounds in three ovarian cancer datasets.

The results are shown in Figure 5. It can be seen that SMC-ABC produces posterior predictive distributions that tightly enclose the time series of tumour volumes for three real-world ovarian cancer datasets. It is evident that SNL produces an accurate posterior predictive distribution for the synthetic dataset, with substantially fewer model simulations than that used for SMC-ABC. This result is aligned with other synthetic examples in the literature (e.g., [106]). However, the SNL results for the real data are mixed, and for the three real datasets SMC-ABC produces more accurate posterior predictive distributions. Further, we do not necessarily see an improvement in SNL when increasing the number of rounds (i.e.,  number of model simulations). We suggest a reason for the potential poor performance of SNL is that the real data are more noisy than than what the simulator is able to produce, and this may lead to a poorly estimated likelihood generated by SNL when evaluated at the observed data. In contrast, SMC-ABC produces results that are more robust to this misspecification, albeit at a higher computational cost in terms of model simulations. The potential poor performance of SNL under misspecification, and possible remedies to this problem, require further research. In terms of computational cost, SMC-ABC takes $\sim$3 hours for each dataset. For SNL, it takes $\sim$5 minutes to generate model simulations, $\sim$20 minutes to train the conditional normalising flow and $\sim$3 hours to generate approximate posterior samples (using the slice sampler as in [89]) for each round. The C++ and Python codes used in this study are available at [107].

During the first year of the COVID-19 pandemic in 2020, border closures, lock-downs and other favourable conditions meant that Australia was spared from the high per capita case numbers and COVID-19 related deaths that were experienced in many other countries. When vaccines became available in February, 2021 [128], the low number of COVID-19 related fatalities in Australia was coupled with uncertainty around highly publicised rare adverse-events for the vaccines; Thrombosis and Thrombocytopenia Syndrome (TTS) from AstraZeneca [129] and myocarditis from Pfizer [130]. This led to high levels of vaccination hesitancy in the general public [131]. Although emerging scientific evidence was increasingly available on the risks of the vaccines and their effectiveness against both becoming infected and becoming severely ill once infected [132] [133], compiling and assessing this information from Scientific journals and Government reports is impossible for the majority of the population. Collating this evidence into an easily understood format that could be used by people to make an informed decision on COVID-19 vaccination in the Australian context became crucial.

Bayesian networks [134] are conditional probability models commonly represented as directed-acyclic graphs, with nodes and links representing variables of interest and the interactions between them. Conditional probabilities for the dependent child-nodes are stored in conditional probability tables (Figure 6), which determine the probability of a node being in a given state for each possible combination of parent node states. Using Bayes Theorem [134], the model calculates the probability of a given outcome for any defined scenario. Bayesian networks are widely used in a range of decision support settings including public health [135, 136], environmental conservation [137] and natural resource management [138].

There are several characteristics that make Bayesian networks attractive for an evidence-based COVID-19 risk-benefit calculator. First the conditional probability tables can be populated from different sources, such as data from government reports, results from scientific studies, or recommendations from experts and advisory committees. Second, the probabilistic output means that the model can respond to user-defined scenarios such as “how likely is it that I will get sick” rather than just “will I get sick”. Finally, Bayesian networks are highly interpretable models [139], as they allow exploration of the effect of different observed values (evidence) on the probability of certain outcomes.

The COVID-19 Risk Calculator (CoRiCAL - https://corical.immunisationcoalition.org.au) was developed to help the general public, as well as the doctors advising them, weigh-up the risks and benefits of receiving a COVID-19 vaccination. A Bayesian network model was constructed and parameterised based on the best available evidence from a range of sources that can be used to determine a person’s risk of developing symptomatic COVID-19, dying or other adverse effects from COVID-19, or suffering from adverse effects (including death) from the vaccine itself [140]. The model relied on Australian data to represent the context as accurately as possible, however in cases where local data was lacking, international data was used [141] [142]. Full model information, along with model code is available via the link [143]. A web-based interface (Figure 7) was developed to create a user-friendly tool that considers a person’s age and sex, the brand of the vaccine, how many vaccines they have had already, and the current levels of transmission within the community and displays their chances of an adverse event alongside common relatable risks. As the pandemic landscape changes, it remains crucial that the evidence for making informed choices on COVID-19 vaccination is made accessible. The model is updated in light of new variants, and as new vaccines become available and recommended (for example booster shots).

Recent projections estimate that 99% of the world’s coral reefs will suffer from frequent marine heatwaves under 1.5°C of warming due to climate change [144]. Important ecological and socio-economic changes already occur in tropical oceans because of the decline of key corals that support thousands of species [145]. The latter impacts about one billion people whose income, food supply, coastal protection, and cultural practices depend on coral reef biodiversity. Robust estimation of changes in coral communities at large spatial scales is challenging because there are a lack of observations due to the remoteness of coral reefs and absence of monitoring programs in sea-countries. Also, the fine-scale variability of changes in coral cover result in disparate long-term coral trajectories at reef locations situated only few hundred meters apart [146]. For reef managers, these challenges (among others) contribute to slowdown the development of strategies that aim to reduce impacts of climate change on coral reefs.

A spatio-temporal Bayesian model is developed to estimate the coverage of total coral cover across spatial scales and predict coverage values at unsampled locations. The approach uses outputs from artificial intelligence algorithms trained to classify points on images [147]. For each Marine Ecoregion Of the World (MEOW, [148]), a set of images $j=1,2,\ldots,J$, each composed of $k=1,2,\ldots,50$ elicitation points is used across years of monitoring. Counts, $y_{it}$ for observation $i$ sampled at location $\boldsymbol{s}_{i}$ and time $t$, are modelled using a binomial distribution (with $p$ the probability of positive and $n_{i}$ the total number of positive cases) and controlled by additional components including the fixed effects of environmental disturbances (cyclones and mass coral bleaching events), sampling nested design (depth, transect, site, reef and monitoring program) modelled as independent and identically distributed Gaussian random effects, and spatio-temporal random effects.

The novelty in this model is the incorporation of a spatio-temporal random effects composed of a first-order autoregressive process in time and a Gaussian field that is approximated using a Gaussian Markov random field (GMRF), where the covariance is determined by a Matérn kernel. We employed the GMRF representation as a stochastic partial differential equation, using the method of [149], implemented in the R package INLA [150]. The spatial domain is based on the observed data locations and a buffer with adjacent MEOWs to allow information sharing between units. Spatial predictions are estimated at a grid level of $5\times 5$ km resolution and posterior distributions used to reconstruct coral cover values at coarser spatial scales including MEOWs units and country level. Finally, estimations of coral cover are weighted by the proportion of coral reefs within a MEOW unit following the methodology developed as part of the Global Coral Reef Monitoring Network [151]. We use the default INLA priors for different types of model parameters as discussed in [152]. The model is as follows,

The priors for the autoregressive parameter $\phi$ and independent Gaussian random effects $V_{i}$ used are the INLA defaults. Research efforts focus on developing new technologies to assess the status of coral reefs in rapid and cost-effective ways through automatic image detection [147] and learn about impacts of multiple disturbances and management strategies [153, 154]. ReefCloud is an open-access digital tool that support coral reef monitoring and decision-making by integration of data analyses and reporting (https://reefcloud.ai/). The online collection of worldwide data provides a unique opportunity to model these data together to 1) increase understanding on the impacts of environmental disturbances and 2) reduce uncertainty when estimating coral trajectories at large spatial scales. The pilot product version is developed using the most extensive monitoring program in the world surveying the Great Barrier Reef, Australia. Machine learning outputs from one million of reef images are used to predict values in coral cover across 3,000 coral reefs from 2004 onward. The ReefCloud online dashboard makes knowledge about reef changes accessible to everyone (Figure 8). The project also educates the reef research community and managers on how Bayesian statistical modelling can help to increase our understanding of the impacts of climate change on coral reefs and supporting decision-making from local to global scales.

Cancer is the leading cause of disease burden in Australia, which has comprehensive cancer incidence reporting for all cancers except common skin cancers [155]. Yet because Australia’s population is heavily concentrated in specific coastal areas, cancer rates are commonly reported only for large regions. Difficulties when using sparse data for smaller areas include the reliability of estimates and the risk of identifying individuals. Yet, detailed small-area analyses have immense power to identify and understand inequities in cancer outcomes.

Using Bayesian hierarchical Poisson models incorporating Leroux priors [156] for spatial smoothing, robust and reliable cancer incidence and 5-year survival estimates were generated across Australian statistical areas level 2 (SA2; 2148 areas). These areas represent communities which interact together and while population sizes vary, the median is around 10,000 people [157]. Innovative visualisations helped rapidly identify areas which differed from the national average. Further details on the methods and visualisations are available in [158]. Example code for the Bayesian spatial models is available in Sections 7.3 and 9.8.2 of [159].

In September 2018 the Australian Cancer Atlas (atlas.cancer.org.au) was launched, providing the highest geographic resolution nationwide estimates available (Figure 9). The website is designed to be interactive and engaging, featuring the ability to download all estimates, export pdfs of specific views, filter regions, rapidly compare different cancer types and rates for two areas, and more! There has been strong uptake and positive feedback, and in 2021 estimates were updated and cancer types expanded.

The Atlas has received prestigious spatial industry awards and is currently being replicated internationally. Australian Cancer Atlas 2.0 is underway, which will examine spatio-temporal patterns, and further include cancer risk factors, some types of cancer treatment and selected cancer clinical/stage patterns. Underpinned by Bayesian methods, the Atlas will continue to provide the methods and visualisations necessary for accurate estimation, interpretation and making decisions.