Approximate Bayesian Computation: A Nonparametric Perspective

Modelling the demographic history of human North African genomes points to a recent soft split divergence between populations

BACKGROUND

North African human populations present a complex demographic scenario due to the presence of an autochthonous genetic component and population substructure, plus extensive gene flow from the Middle East, Europe, and sub-Saharan Africa.

RESULTS

We conducted a comprehensive analysis of 364 genomes to construct detailed demographic models for the North African region, encompassing its two primary ethnic groups, the Arab and Amazigh populations. This was achieved through an Approximate Bayesian Computation with Deep Learning (ABC-DL) framework and a novel algorithm called Genetic Programming for Population Genetics (GP4PG). This innovative approach enabled us to effectively model intricate demographic scenarios, utilizing a subset of 16 whole genomes at > 30X coverage. The demographic model suggested by GP4PG exhibited a closer alignment with the observed data compared to the ABC-DL model. Both point to a back-to-Africa origin of North African individuals and a close relationship with Eurasian populations. Results support different origins for Amazigh and Arab populations, with Amazigh populations originating back in Epipaleolithic times, while GP4PG supports Arabization as the main source of Middle Eastern ancestry. The GP4PG model includes population substructure in surrounding populations (sub-Saharan Africa and Middle East) with continuous decaying gene flow after population split. Contrary to ABC-DL, the best GP4PG model does not require pulses of admixture from surrounding populations into North Africa pointing to soft splits as drivers of divergence in North Africa.

CONCLUSIONS

We have built a demographic model on North Africa that points to a back-to-Africa expansion and a differential origin between Arab and Amazigh populations.

Flexible Bayesian inference on partially observed epidemics

Individual-based models of contagious processes are useful for predicting epidemic trajectories and informing intervention strategies. In such models, the incorporation of contact network information can capture the non-randomness and heterogeneity of realistic contact dynamics. In this article, we consider Bayesian inference on the spreading parameters of an SIR contagion on a known, static network, where information regarding individual disease status is known only from a series of tests (positive or negative disease status). When the contagion model is complex or information such as infection and removal times is missing, the posterior distribution can be difficult to sample from. Previous work has considered the use of Approximate Bayesian Computation (ABC), which allows for simulation-based Bayesian inference on complex models. However, ABC methods usually require the user to select reasonable summary statistics. Here, we consider an inference scheme based on the Mixture Density Network compressed ABC, which minimizes the expected posterior entropy in order to learn informative summary statistics. This allows us to conduct Bayesian inference on the parameters of a partially observed contagious process while also circumventing the need for manual summary statistic selection. This methodology can be extended to incorporate additional simulation complexities, including behavioural change after positive tests or false test results.

Bat Motion can be Described by Leap Frogging

We present models of bat motion derived from radio-tracking data collected over 14 nights. The data presents an initial dispersal period and a return to roost period. Although a simple diffusion model fits the initial dispersal motion we show that simple convection cannot provide a description of the bats returning to their roost. By extending our model to include non-autonomous parameters, or a leap frogging form of motion, where bats on the exterior move back first, we find we are able to accurately capture the bat's motion. We discuss ways of distinguishing between the two movement descriptions and, finally, consider how the different motion descriptions would impact a bat's hunting strategy.

On predictive inference for intractable models via approximate Bayesian computation

Approximate Bayesian computation (ABC) is commonly used for parameter estimation and model comparison for intractable simulator-based statistical models whose likelihood function cannot be evaluated. In this paper we instead investigate the feasibility of ABC as a generic approximate method for predictive inference, in particular, for computing the posterior predictive distribution of future observations or missing data of interest. We consider three complementary ABC approaches for this goal, each based on different assumptions regarding which predictive density of the intractable model can be sampled from. The case where only simulation from the joint density of the observed and future data given the model parameters can be used for inference is given particular attention and it is shown that the ideal summary statistic in this setting is minimal predictive sufficient instead of merely minimal sufficient (in the ordinary sense). An ABC prediction approach that takes advantage of a certain latent variable representation is also investigated. We additionally show how common ABC sampling algorithms can be used in the predictive settings considered. Our main results are first illustrated by using simple time-series models that facilitate analytical treatment, and later by using two common intractable dynamic models.

SUPPLEMENTARY INFORMATION

The online version contains supplementary material available at 10.1007/s11222-022-10163-6.

Retrospective exposure reconstruction using approximate Bayesian computation: A case study on perfluorooctanoic acid and preeclampsia

BACKGROUND

In environmental epidemiology, measurements of toxicants in biological samples are often used as individual exposure assignments. It is common to obtain only one or a few exposure biomarkers per person and use those measurements to represent each person's relevant toxicant exposure for a given health outcome, even though most exposure biomarkers can fluctuate over time. When the timing of the exposure reflected by the biomarker measurement is misaligned with disease development especially if it occurs after the disease outcome, results could be subject to reverse causality or exposure measurement error.

OBJECTIVE

This study aimed to use an approximate Bayesian computation (ABC) method to improve PFOA exposure estimates and characterize the effects of PFOA on preeclampsia in the C8 Studies.

METHODS

Serum PFOA concentrations were measured in blood samples collected during 2005-2006 in West Virginia and Ohio (the C8 Studies), and residential and water use histories and pregnancy outcomes were obtained from self-reports. Our previous results may have been influenced by the choice of methods for characterizing PFOA exposures. Here we use an ABC method to combine measured PFOA serum concentrations and environmentally modeled PFOA concentrations to reconstruct historical PFOA exposures. We also expanded our previous work by assuming more realistic lognormal distributions for key input parameters in the exposure and pharmacokinetic models.

RESULTS

Compared to using fixed values of model parameters and Monte Carlo simulations, ABC produced similar Spearman correlations between estimated and measured serum PFOA concentrations, yet substantially reduced the mean squared error by over 50%. Based on ABC, compared to previous studies, we found a similar adjusted odds ratio (AOR) for the association between PFOA and preeclampsia.

CONCLUSIONS

Bayesian combination of modeled exposure and measured biomarker concentrations can reduce exposure measurement error compared to modeled exposure.

Optimal Bayesian design for model discrimination via classification

Performing optimal Bayesian design for discriminating between competing models is computationally intensive as it involves estimating posterior model probabilities for thousands of simulated data sets. This issue is compounded further when the likelihood functions for the rival models are computationally expensive. A new approach using supervised classification methods is developed to perform Bayesian optimal model discrimination design. This approach requires considerably fewer simulations from the candidate models than previous approaches using approximate Bayesian computation. Further, it is easy to assess the performance of the optimal design through the misclassification error rate. The approach is particularly useful in the presence of models with intractable likelihoods but can also provide computational advantages when the likelihoods are manageable.

SUPPLEMENTARY INFORMATION

The online version contains supplementary material available at 10.1007/s11222-022-10078-2.

Inference with selection, varying population size, and evolving population structure: application of ABC to a forward-backward coalescent process with interactions

Genetic data are often used to infer demographic history and changes or detect genes under selection. Inferential methods are commonly based on models making various strong assumptions: demography and population structures are supposed a priori known, the evolution of the genetic composition of a population does not affect demography nor population structure, and there is no selection nor interaction between and within genetic strains. In this paper, we present a stochastic birth-death model with competitive interactions and asexual reproduction. We develop an inferential procedure for ecological, demographic, and genetic parameters. We first show how genetic diversity and genealogies are related to birth and death rates, and to how individuals compete within and between strains. This leads us to propose an original model of phylogenies, with trait structure and interactions, that allows multiple merging. Second, we develop an Approximate Bayesian Computation framework to use our model for analyzing genetic data. We apply our procedure to simulated data from a toy model, and to real data by analyzing the genetic diversity of microsatellites on Y-chromosomes sampled from Central Asia human populations in order to test whether different social organizations show significantly different fertilities.

Simulation and inference algorithms for stochastic biochemical reaction networks: from basic concepts to state-of-the-art

Stochasticity is a key characteristic of intracellular processes such as gene regulation and chemical signalling. Therefore, characterizing stochastic effects in biochemical systems is essential to understand the complex dynamics of living things. Mathematical idealizations of biochemically reacting systems must be able to capture stochastic phenomena. While robust theory exists to describe such stochastic models, the computational challenges in exploring these models can be a significant burden in practice since realistic models are analytically intractable. Determining the expected behaviour and variability of a stochastic biochemical reaction network requires many probabilistic simulations of its evolution. Using a biochemical reaction network model to assist in the interpretation of time-course data from a biological experiment is an even greater challenge due to the intractability of the likelihood function for determining observation probabilities. These computational challenges have been subjects of active research for over four decades. In this review, we present an accessible discussion of the major historical developments and state-of-the-art computational techniques relevant to simulation and inference problems for stochastic biochemical reaction network models. Detailed algorithms for particularly important methods are described and complemented with Matlab^® implementations. As a result, this review provides a practical and accessible introduction to computational methods for stochastic models within the life sciences community.

Towards end-to-end likelihood-free inference with convolutional neural networks

Complex simulator-based models with non-standard sampling distributions require sophisticated design choices for reliable approximate parameter inference. We introduce a fast, end-to-end approach for approximate Bayesian computation (ABC) based on fully convolutional neural networks. The method enables users of ABC to derive simultaneously the posterior mean and variance of multidimensional posterior distributions directly from raw simulated data. Once trained on simulated data, the convolutional neural network is able to map real data samples of variable size to the first two posterior moments of the relevant parameter's distributions. Thus, in contrast to other machine learning approaches to ABC, our approach allows us to generate reusable models that can be applied by different researchers employing the same model. We verify the utility of our method on two common statistical models (i.e., a multivariate normal distribution and a multiple regression scenario), for which the posterior parameter distributions can be derived analytically. We then apply our method to recover the parameters of the leaky competing accumulator (LCA) model and we reference our results to the current state-of-the-art technique, which is the probability density estimation (PDA). Results show that our method exhibits a lower approximation error compared with other machine learning approaches to ABC. It also performs similarly to PDA in recovering the parameters of the LCA model.

An ABC Method for Whole-Genome Sequence Data: Inferring Paleolithic and Neolithic Human Expansions

Species generally undergo a complex demographic history consisting, in particular, of multiple changes in population size. Genome-wide sequencing data are potentially highly informative for reconstructing this demographic history. A crucial point is to extract the relevant information from these very large data sets. Here, we design an approach for inferring past demographic events from a moderate number of fully sequenced genomes. Our new approach uses Approximate Bayesian Computation, a simulation-based statistical framework that allows 1) identifying the best demographic scenario among several competing scenarios and 2) estimating the best-fitting parameters under the chosen scenario. Approximate Bayesian Computation relies on the computation of summary statistics. Using a cross-validation approach, we show that statistics such as the lengths of haplotypes shared between individuals, or the decay of linkage disequilibrium with distance, can be combined with classical statistics (e.g., heterozygosity and Tajima's D) to accurately infer complex demographic scenarios including bottlenecks and expansion periods. We also demonstrate the importance of simultaneously estimating the genotyping error rate. Applying our method on genome-wide human-sequence databases, we finally show that a model consisting in a bottleneck followed by a Paleolithic and a Neolithic expansion is the most relevant for Eurasian populations.

Spectral density-based and measure-preserving ABC for partially observed diffusion processes. An illustration on Hamiltonian SDEs

Approximate Bayesian computation (ABC) has become one of the major tools of likelihood-free statistical inference in complex mathematical models. Simultaneously, stochastic differential equations (SDEs) have developed to an established tool for modelling time-dependent, real-world phenomena with underlying random effects. When applying ABC to stochastic models, two major difficulties arise: First, the derivation of effective summary statistics and proper distances is particularly challenging, since simulations from the stochastic process under the same parameter configuration result in different trajectories. Second, exact simulation schemes to generate trajectories from the stochastic model are rarely available, requiring the derivation of suitable numerical methods for the synthetic data generation. To obtain summaries that are less sensitive to the intrinsic stochasticity of the model, we propose to build up the statistical method (e.g. the choice of the summary statistics) on the underlying structural properties of the model. Here, we focus on the existence of an invariant measure and we map the data to their estimated invariant density and invariant spectral density. Then, to ensure that these model properties are kept in the synthetic data generation, we adopt measure-preserving numerical splitting schemes. The derived property-based and measure-preserving ABC method is illustrated on the broad class of partially observed Hamiltonian type SDEs, both with simulated data and with real electroencephalography data. The derived summaries are particularly robust to the model simulation, and this fact, combined with the proposed reliable numerical scheme, yields accurate ABC inference. In contrast, the inference returned using standard numerical methods (Euler-Maruyama discretisation) fails. The proposed ingredients can be incorporated into any type of ABC algorithm and directly applied to all SDEs that are characterised by an invariant distribution and for which a measure-preserving numerical method can be derived.

Resolving outbreak dynamics using approximate Bayesian computation for stochastic birth-death models

Earlier research has suggested that approximate Bayesian computation (ABC) makes it possible to fit simulator-based intractable birth-death models to investigate communicable disease outbreak dynamics with accuracy comparable to that of exact Bayesian methods. However, recent findings have indicated that key parameters, such as the reproductive number R, may remain poorly identifiable with these models. Here we show that this identifiability issue can be resolved by taking into account disease-specific characteristics of the transmission process in closer detail. Using tuberculosis (TB) in the San Francisco Bay area as a case study, we consider a model that generates genotype data from a mixture of three stochastic processes, each with its own distinct dynamics and clear epidemiological interpretation.       We show that our model allows for accurate posterior inferences about outbreak dynamics from aggregated annual case data with genotype information. As a byproduct of the inference, the model provides an estimate of the infectious population size at the time the data were collected. The acquired estimate is approximately two orders of magnitude smaller than assumed in earlier related studies, and it is much better aligned with epidemiological knowledge about active TB prevalence. Similarly, the reproductive number R related to the primary underlying transmission process is estimated to be nearly three times larger than previous estimates, which has a substantial impact on the interpretation of the fitted outbreak model.

Approximate Bayesian inference for discretely observed continuous-time multi-state models

Inference for continuous time multi-state models presents considerable computational difficulties when the process is only observed at discrete time points with no additional information about the state transitions. In fact, for general multi-state Markov model, evaluation of the likelihood function is possible only via intensive numerical approximations. Moreover, in real applications, transitions between states may depend on the time since entry into the current state, and semi-Markov models, where the likelihood function is not available in closed form, should be fitted to the data. Approximate Bayesian Computation (ABC) methods, which make use only of comparisons between simulated and observed summary statistics, represent a solution to intractable likelihood problems and provide alternative algorithms when the likelihood calculation is computationally too costly. In this article we investigate the potentiality of ABC techniques for multi-state models both to obtain the posterior distributions of the model parameters and to compare Markov and semi-Markov models. In addition, we will also exploit ABC methods to estimate and compare hidden Markov and semi-Markov models when observed states are subject to classification errors. We illustrate the performance of the ABC methodology both with simulated data and with a real data example.

Resolving outbreak dynamics using approximate Bayesian computation for stochastic birth-death models

Earlier research has suggested that approximate Bayesian computation (ABC) makes it possible to fit simulator-based intractable birth-death models to investigate communicable disease outbreak dynamics with accuracy comparable to that of exact Bayesian methods. However, recent findings have indicated that key parameters such as the reproductive number R may remain poorly identifiable with these models. Here we show that the identifiability issue can be resolved by taking into account disease-specific characteristics of the transmission process in closer detail. Using tuberculosis (TB) in the San Francisco Bay area as a case-study, we consider a model that generates genotype data from a mixture of three stochastic processes, each with their distinct dynamics and clear epidemiological interpretation.  We show that our model allows for accurate posterior inferences about outbreak dynamics from aggregated annual case data with genotype information.  As a by-product of the inference, the model provides an estimate of the infectious population size at the time the data was collected. The acquired estimate is approximately two orders of magnitude smaller compared to the assumptions made in the earlier related studies, and much better aligned with epidemiological knowledge about active TB prevalence. Similarly, the reproductive number R related to the primary underlying transmission process is estimated to be nearly three-fold compared with the previous estimates, which has a substantial impact on the interpretation of the fitted outbreak model.

Genomic signatures of environmental selection despite near-panmixia in summer flounder

Rapid environmental change is altering the selective pressures experienced by marine species. While adaptation to local environmental conditions depends on a balance between dispersal and natural selection across the seascape, the spatial scale of adaptation and the relative importance of mechanisms maintaining adaptation in the ocean are not well understood. Here, using population assignment tests, Approximate Bayesian Computation (ABC), and genome scans with double-digest restriction-site associated DNA sequencing data, we evaluated population structure and locus-environment associations in a commercially important species, summer flounder (Paralichthys dentatus), along the U.S. east coast. Based on 1,137 single nucleotide polymorphisms across 232 individuals spanning nearly 1,900 km, we found no indication of population structure across Cape Hatteras, North Carolina (F ST = 0.0014) or of isolation by distance along the coast using individual relatedness. ABC estimated the probability of dispersal across the biogeographic break at Cape Hatteras to be high (95% credible interval: 7%-50% migration). However, we found 15 loci whose allele frequencies were associated with at least one of four environmental variables. Of those, 11 were correlated with bottom temperature. For summer flounder, our results suggest continued fisheries management as a single population and identify likely response mechanisms to climate change. Broadly speaking, our findings suggest that spatial balancing selection can manifest in adaptive divergence on regional scales in marine fish despite high dispersal, and that these conditions likely result in the widespread distribution of adaptive alleles and a high potential for future genetic adaptation in response to changing environmental conditions. In the context of a rapidly changing world, a landscape genomics perspective offers a useful approach for understanding the causes and consequences of genetic differentiation.

Fundamentals and Recent Developments in Approximate Bayesian Computation

Bayesian inference plays an important role in phylogenetics, evolutionary biology, and in many other branches of science. It provides a principled framework for dealing with uncertainty and quantifying how it changes in the light of new evidence. For many complex models and inference problems, however, only approximate quantitative answers are obtainable. Approximate Bayesian computation (ABC) refers to a family of algorithms for approximate inference that makes a minimal set of assumptions by only requiring that sampling from a model is possible. We explain here the fundamentals of ABC, review the classical algorithms, and highlight recent developments. [ABC; approximate Bayesian computation; Bayesian inference; likelihood-free inference; phylogenetics; simulator-based models; stochastic simulation models; tree-based models.]

Statistical Inference in the Wright-Fisher Model Using Allele Frequency Data

The Wright–Fisher model provides an elegant mathematical framework for understanding allele frequency data. In particular, the model can be used to infer the demographic history of species and identify loci under selection. A crucial quantity for inference under the Wright–Fisher model is the distribution of allele frequencies (DAF). Despite the apparent simplicity of the model, the calculation of the DAF is challenging. We review and discuss strategies for approximating the DAF, and how these are used in methods that perform inference from allele frequency data. Various evolutionary forces can be incorporated in the Wright–Fisher model, and we consider these in turn. We begin our review with the basic bi-allelic Wright–Fisher model where random genetic drift is the only evolutionary force. We then consider mutation, migration, and selection. In particular, we compare diffusion-based and moment-based methods in terms of accuracy, computational efficiency, and analytical tractability. We conclude with a brief overview of the multi-allelic process with a general mutation model.

Inferring responses to climate dynamics from historical demography in neotropical forest lizards

We apply a comparative framework to test for concerted demographic changes in response to climate shifts in the neotropical lowland forests, learning from the past to inform projections of the future. Using reduced genomic (SNP) data from three lizard species codistributed in Amazonia and the Atlantic Forest (Anolis punctatus, Anolis ortonii, and Polychrus marmoratus), we first reconstruct former population history and test for assemblage-level responses to cycles of moisture transport recently implicated in changes of forest distribution during the Late Quaternary. We find support for population shifts within the time frame of inferred precipitation fluctuations (the last 250,000 y) but detect idiosyncratic responses across species and uniformity of within-species responses across forest regions. These results are incongruent with expectations of concerted population expansion in response to increased rainfall and fail to detect out-of-phase demographic syndromes (expansions vs. contractions) across forest regions. Using reduced genomic data to infer species-specific demographical parameters, we then model the plausible spatial distribution of genetic diversity in the Atlantic Forest into future climates (2080) under a medium carbon emission trajectory. The models forecast very distinct trajectories for the lizard species, reflecting unique estimated population densities and dispersal abilities. Ecological and demographic constraints seemingly lead to distinct and asynchronous responses to climatic regimes in the tropics, even among similarly distributed taxa. Incorporating such constraints is key to improve modeling of the distribution of biodiversity in the past and future.

Demographic model selection using random forests and the site frequency spectrum

Phylogeographic data sets have grown from tens to thousands of loci in recent years, but extant statistical methods do not take full advantage of these large data sets. For example, approximate Bayesian computation (ABC) is a commonly used method for the explicit comparison of alternate demographic histories, but it is limited by the "curse of dimensionality" and issues related to the simulation and summarization of data when applied to next-generation sequencing (NGS) data sets. We implement here several improvements to overcome these difficulties. We use a Random Forest (RF) classifier for model selection to circumvent the curse of dimensionality and apply a binned representation of the multidimensional site frequency spectrum (mSFS) to address issues related to the simulation and summarization of large SNP data sets. We evaluate the performance of these improvements using simulation and find low overall error rates (~7%). We then apply the approach to data from Haplotrema vancouverense, a land snail endemic to the Pacific Northwest of North America. Fifteen demographic models were compared, and our results support a model of recent dispersal from coastal to inland rainforests. Our results demonstrate that binning is an effective strategy for the construction of a mSFS and imply that the statistical power of RF when applied to demographic model selection is at least comparable to traditional ABC algorithms. Importantly, by combining these strategies, large sets of models with differing numbers of populations can be evaluated.

A tutorial introduction to Bayesian inference for stochastic epidemic models using Approximate Bayesian Computation

Likelihood-based inference for disease outbreak data can be very challenging due to the inherent dependence of the data and the fact that they are usually incomplete. In this paper we review recent Approximate Bayesian Computation (ABC) methods for the analysis of such data by fitting to them stochastic epidemic models without having to calculate the likelihood of the observed data. We consider both non-temporal and temporal-data and illustrate the methods with a number of examples featuring different models and datasets. In addition, we present extensions to existing algorithms which are easy to implement and provide an improvement to the existing methodology. Finally, R code to implement the algorithms presented in the paper is available on https://github.com/kypraios/epiABC.

Likelihood-Free Inference in High-Dimensional Models

Methods that bypass analytical evaluations of the likelihood function have become an indispensable tool for statistical inference in many fields of science. These so-called likelihood-free methods rely on accepting and rejecting simulations based on summary statistics, which limits them to low-dimensional models for which the value of the likelihood is large enough to result in manageable acceptance rates. To get around these issues, we introduce a novel, likelihood-free Markov chain Monte Carlo (MCMC) method combining two key innovations: updating only one parameter per iteration and accepting or rejecting this update based on subsets of statistics approximately sufficient for this parameter. This increases acceptance rates dramatically, rendering this approach suitable even for models of very high dimensionality. We further derive that for linear models, a one-dimensional combination of statistics per parameter is sufficient and can be found empirically with simulations. Finally, we demonstrate that our method readily scales to models of very high dimensionality, using toy models as well as by jointly inferring the effective population size, the distribution of fitness effects (DFE) of segregating mutations, and selection coefficients for each locus from data of a recent experiment on the evolution of drug resistance in influenza.