Approximate Bayesian computation

A common incarnation of Bayes’ theorem relates the conditional probability (or density) of a particular parameter value ${\displaystyle \theta }$ given data ${\displaystyle D}$ to the probability of ${\displaystyle D}$ given ${\displaystyle \theta }$ by the rule

${\displaystyle p(\theta |D)={\frac {p(D|\theta )p(\theta )}{p(D)}}}$,

where ${\displaystyle p(\theta |D)}$ denotes the posterior, ${\displaystyle p(D|\theta )}$ the likelihood, ${\displaystyle p(\theta )}$ the prior, and ${\displaystyle p(D)}$ the evidence (also referred to as the marginal likelihood or the prior predictive probability of the data). Note that the denominator ${\displaystyle p(D)}$ is normalizing the total probability of the posterior density ${\displaystyle p(\theta |D)}$ to one and can be calculated that way.

The prior represents beliefs or knowledge (such as f.e. physical constraints) about ${\displaystyle \theta }$ before ${\displaystyle D}$ is available. Since the prior narrows down uncertainty, the posterior estimates have less variance, but might be biased. For convenience the prior is often specified by choosing a particular distribution among a set of well-known and tractable families of distributions, such that both the evaluation of prior probabilities and random generation of values of ${\displaystyle \theta }$ are relatively straightforward. For certain kinds of models, it is more pragmatic to specify the prior ${\displaystyle p(\theta )}$ using a factorization of the joint distribution of all the elements of ${\displaystyle \theta }$ in terms of a sequence of their conditional distributions. If one is only interested in the relative posterior plausibilities of different values of ${\displaystyle \theta }$, the evidence ${\displaystyle p(D)}$ can be ignored, as it constitutes a normalising constant, which cancels for any ratio of posterior probabilities. It remains, however, necessary to evaluate the likelihood ${\displaystyle p(D|\theta )}$ and the prior ${\displaystyle p(\theta )}$. For numerous applications, it is computationally expensive, or even completely infeasible, to evaluate the likelihood,[11] which motivates the use of ABC to circumvent this issue.

All ABC-based methods approximate the likelihood function by simulations, the outcomes of which are compared with the observed data.[12][13][14] More specifically, with the ABC rejection algorithm — the most basic form of ABC — a set of parameter points is first sampled from the prior distribution. Given a sampled parameter point ${\displaystyle {\hat {\theta }}}$, a data set ${\displaystyle {\hat {D}}}$ is then simulated under the statistical model ${\displaystyle M}$ specified by ${\displaystyle {\hat {\theta }}}$. If the generated ${\displaystyle {\hat {D}}}$ is too different from the observed data ${\displaystyle D}$, the sampled parameter value is discarded. In precise terms, ${\displaystyle {\hat {D}}}$ is accepted with tolerance ${\displaystyle \epsilon \geq 0}$ if:

${\displaystyle \rho ({\hat {D}},D)\leq \epsilon }$,

where the distance measure ${\displaystyle \rho ({\hat {D}},D)}$ determines the level of discrepancy between ${\displaystyle {\hat {D}}}$ and ${\displaystyle D}$ based on a given metric (e.g. Euclidean distance). A strictly positive tolerance is usually necessary, since the probability that the simulation outcome coincides exactly with the data (event ${\displaystyle {\hat {D}}=D}$) is negligible for all but trivial applications of ABC, which would in practice lead to rejection of nearly all sampled parameter points. The outcome of the ABC rejection algorithm is a sample of parameter values approximately distributed according to the desired posterior distribution, and, crucially, obtained without the need to explicitly evaluate the likelihood function.

Parameter estimation by approximate Bayesian computation: a conceptual overview.

The probability of generating a data set ${\displaystyle {\hat {D}}}$ with a small distance to ${\displaystyle D}$ typically decreases as the dimensionality of the data increases. This leads to a substantial decrease in the computational efficiency of the above basic ABC rejection algorithm. A common approach to lessen this problem is to replace ${\displaystyle D}$ with a set of lower-dimensional summary statistics ${\displaystyle S(D)}$, which are selected to capture the relevant information in ${\displaystyle D}$. The acceptance criterion in ABC rejection algorithm becomes:

${\displaystyle \rho (S({\hat {D}}),S(D))\leq \epsilon }$.

If the summary statistics are sufficient with respect to the model parameters ${\displaystyle \theta }$, the efficiency increase obtained in this way does not introduce any error.[15] Indeed, by definition, sufficiency implies that all information in ${\displaystyle D}$ about ${\displaystyle \theta }$ is captured by ${\displaystyle S(D)}$.

As elaborated below, it is typically impossible, outside the exponential family of distributions, to identify a finite-dimensional set of sufficient statistics. Nevertheless, informative but possibly insufficient summary statistics are often used in applications where inference is performed with ABC methods.

A non-negligible ${\displaystyle \epsilon }$ comes with the price that one samples from ${\displaystyle p(\theta |\rho ({\hat {D}},D)\leq \epsilon )}$ instead of the true posterior ${\displaystyle p(\theta |D)}$. With a sufficiently small tolerance, and a sensible distance measure, the resulting distribution ${\displaystyle p(\theta |\rho ({\hat {D}},D)\leq \epsilon )}$ should often approximate the actual target distribution ${\displaystyle p(\theta |D)}$ reasonably well. On the other hand, a tolerance that is large enough that every point in the parameter space becomes accepted will yield a replica of the prior distribution. There are empirical studies of the difference between ${\displaystyle p(\theta |\rho ({\hat {D}},D)\leq \epsilon )}$ and ${\displaystyle p(\theta |D)}$ as a function of ${\displaystyle \epsilon }$,[26] and theoretical results for an upper ${\displaystyle \epsilon }$-dependent bound for the error in parameter estimates.[27] The accuracy of the posterior (defined as the expected quadratic loss) delivered by ABC as a function of ${\displaystyle \epsilon }$ has also been investigated.[28] However, the convergence of the distributions when ${\displaystyle \epsilon }$ approaches zero, and how it depends on the distance measure used, is an important topic that has yet to be investigated in greater detail. In particular, it remains difficult to disentangle errors introduced by this approximation from errors due to model mis-specification.[12]

As an attempt to correct some of the error due to a non-zero ${\displaystyle \epsilon }$, the usage of local linear weighted regression with ABC to reduce the variance of the posterior estimates has been suggested.[9] The method assigns weights to the parameters according to how well simulated summaries adhere to the observed ones and performs linear regression between the summaries and the weighted parameters in the vicinity of observed summaries. The obtained regression coefficients are used to correct sampled parameters in the direction of observed summaries. An improvement was suggested in the form of nonlinear regression using a feed-forward neural network model.[29] However, it has been shown that the posterior distributions obtained with these approaches are not always consistent with the prior distribution, which did lead to a reformulation of the regression adjustment that respects the prior distribution.[30]

Finally, statistical inference using ABC with a non-zero tolerance ${\displaystyle \epsilon }$ is not inherently flawed: under the assumption of measurement errors, the optimal ${\displaystyle \epsilon }$ can in fact be shown to be not zero.[28][31] Indeed, the bias caused by a non-zero tolerance can be characterized and compensated by introducing a specific form of noise to the summary statistics. Asymptotic consistency for such “noisy ABC”, has been established, together with formulas for the asymptotic variance of the parameter estimates for a fixed tolerance.[28]

Summary statistics may be used to increase the acceptance rate of ABC for high-dimensional data. Low-dimensional sufficient statistics are optimal for this purpose, as they capture all relevant information present in the data in the simplest possible form.[14] However, low-dimensional sufficient statistics are typically unattainable for statistical models where ABC-based inference is most relevant, and consequently, some heuristic is usually necessary to identify useful low-dimensional summary statistics. The use of a set of poorly chosen summary statistics will often lead to inflated credible intervals due to the implied loss of information,[14] which can also bias the discrimination between models. A review of methods for choosing summary statistics is available,[32] which may provide valuable guidance in practice.

One approach to capture most of the information present in data would be to use many statistics, but the accuracy and stability of ABC appears to decrease rapidly with an increasing numbers of summary statistics.[12][14] Instead, a better strategy is to focus on the relevant statistics only—relevancy depending on the whole inference problem, on the model used, and on the data at hand.[33]

An algorithm has been proposed for identifying a representative subset of summary statistics, by iteratively assessing whether an additional statistic introduces a meaningful modification of the posterior.[34] One of the challenges here is that a large ABC approximation error may heavily influence the conclusions about the usefulness of a statistic at any stage of the procedure. Another method[33] decomposes into two main steps. First, a reference approximation of the posterior is constructed by minimizing the entropy. Sets of candidate summaries are then evaluated by comparing the ABC-approximated posteriors with the reference posterior.

With both of these strategies, a subset of statistics is selected from a large set of candidate statistics. Instead, the partial least squares regression approach uses information from all the candidate statistics, each being weighted appropriately.[35] Recently, a method for constructing summaries in a semi-automatic manner has attained a considerable interest.[28] This method is based on the observation that the optimal choice of summary statistics, when minimizing the quadratic loss of the parameter point estimates, can be obtained through the posterior mean of the parameters, which is approximated by performing a linear regression based on the simulated data.

Methods for the identification of summary statistics that could also simultaneously assess the influence on the approximation of the posterior would be of substantial value.[36] This is because the choice of summary statistics and the choice of tolerance constitute two sources of error in the resulting posterior distribution. These errors may corrupt the ranking of models and may also lead to incorrect model predictions. Indeed, none of the methods above assesses the choice of summaries for the purpose of model selection.

It has been shown that the combination of insufficient summary statistics and ABC for model selection can be problematic.[15][22] Indeed, if one lets the Bayes factor based on the summary statistic ${\displaystyle S(D)}$ be denoted by ${\displaystyle B_{1,2}^{s}}$, the relation between ${\displaystyle B_{1,2}}$ and ${\displaystyle B_{1,2}^{s}}$ takes the form:[15]

${\displaystyle B_{1,2}={\frac {p(D|M_{1})}{p(D|M_{2})}}={\frac {p(D|S(D),M_{1})}{p(D|S(D),M_{2})}}{\frac {p(S(D)|M_{1})}{p(S(D)|M_{2})}}={\frac {p(D|S(D),M_{1})}{p(D|S(D),M_{2})}}B_{1,2}^{s}}$.

Thus, a summary statistic ${\displaystyle S(D)}$ is sufficient for comparing two models ${\displaystyle M_{1}}$ and ${\displaystyle M_{2}}$ if and only if:

${\displaystyle p(D|S(D),M_{1})=p(D|S(D),M_{2})}$,

which results in that ${\displaystyle B_{1,2}=B_{1,2}^{s}}$. It is also clear from the equation above that there might be a huge difference between ${\displaystyle B_{1,2}}$ and ${\displaystyle B_{1,2}^{s}}$ if the condition is not satisfied, as can be demonstrated by toy examples.[15][19][22] Crucially, it was shown that sufficiency for ${\displaystyle M_{1}}$ or ${\displaystyle M_{2}}$ alone, or for both models, does not guarantee sufficiency for ranking the models.[15] However, it was also shown that any sufficient summary statistic for a model ${\displaystyle M}$ in which both ${\displaystyle M_{1}}$ and ${\displaystyle M_{2}}$ are nested is valid for ranking the nested models.[15]

The computation of Bayes factors on ${\displaystyle S(D)}$ may therefore be misleading for model selection purposes, unless the ratio between the Bayes factors on ${\displaystyle D}$ and ${\displaystyle S(D)}$ would be available, or at least could be approximated reasonably well. Alternatively, necessary and sufficient conditions on summary statistics for a consistent Bayesian model choice have recently been derived,[37] which can provide useful guidance.

However, this issue is only relevant for model selection when the dimension of the data has been reduced. ABC-based inference, in which the actual data sets are directly compared—as is the case for some systems biology applications (e.g., see [38])—circumvents this problem.

As the above discussion makes clear, any ABC analysis requires choices and trade-offs that can have a considerable impact on its outcomes. Specifically, the choice of competing models/hypotheses, the number of simulations, the choice of summary statistics, or the acceptance threshold cannot currently be based on general rules, but the effect of these choices should be evaluated and tested in each study.[13]

A number of heuristic approaches to the quality control of ABC have been proposed, such as the quantification of the fraction of parameter variance explained by the summary statistics.[13] A common class of methods aims at assessing whether or not the inference yields valid results, regardless of the actually observed data. For instance, given a set of parameter values, which are typically drawn from the prior or the posterior distributions for a model, one can generate a large number of artificial datasets. In this way, the quality and robustness of ABC inference can be assessed in a controlled setting, by gauging how well the chosen ABC inference method recovers the true parameter values, and also models if multiple structurally different models are considered simultaneously.

Another class of methods assesses whether the inference was successful in light of the given observed data, for example, by comparing the posterior predictive distribution of summary statistics to the summary statistics observed.[13] Beyond that, cross-validation techniques[39] and predictive checks[40][41] represent promising future strategies to evaluate the stability and out-of-sample predictive validity of ABC inferences. This is particularly important when modeling large data sets, because then the posterior support of a particular model can appear overwhelmingly conclusive, even if all proposed models in fact are poor representations of the stochastic system underlying the observation data. Out-of-sample predictive checks can reveal potential systematic biases within a model and provide clues on to how to improve its structure or parametrization.

Fundamentally novel approaches for model choice that incorporate quality control as an integral step in the process have recently been proposed. ABC allows, by construction, estimation of the discrepancies between the observed data and the model predictions, with respect to a comprehensive set of statistics. These statistics are not necessarily the same as those used in the acceptance criterion. The resulting discrepancy distributions have been used for selecting models that are in agreement with many aspects of the data simultaneously,[42] and model inconsistency is detected from conflicting and co-dependent summaries. Another quality-control-based method for model selection employs ABC to approximate the effective number of model parameters and the deviance of the posterior predictive distributions of summaries and parameters.[43] The deviance information criterion is then used as measure of model fit. It has also been shown that the models preferred based on this criterion can conflict with those supported by Bayes factors. For this reason, it is useful to combine different methods for model selection to obtain correct conclusions.

Quality controls are achievable and indeed performed in many ABC-based works, but for certain problems, the assessment of the impact of the method-related parameters can be challenging. However, the rapidly increasing use of ABC can be expected to provide a more thorough understanding of the limitations and applicability of the method.

This section reviews risks that are strictly speaking not specific to ABC, but also relevant for other statistical methods as well. However, the flexibility offered by ABC to analyze very complex models makes them highly relevant to discuss here.

The specification of the range and the prior distribution of parameters strongly benefits from previous knowledge about the properties of the system. One criticism has been that in some studies the “parameter ranges and distributions are only guessed based upon the subjective opinion of the investigators”,[44] which is connected to classical objections of Bayesian approaches.[45]

With any computational method, it is typically necessary to constrain the investigated parameter ranges. The parameter ranges should if possible be defined based on known properties of the studied system, but may for practical applications necessitate an educated guess. However, theoretical results regarding objective priors are available, which may for example be based on the principle of indifference or the principle of maximum entropy.[46][47] On the other hand, automated or semi-automated methods for choosing a prior distribution often yield improper densities. As most ABC procedures require generating samples from the prior, improper priors are not directly applicable to ABC.

One should also keep the purpose of the analysis in mind when choosing the prior distribution. In principle, uninformative and flat priors, that exaggerate our subjective ignorance about the parameters, may still yield reasonable parameter estimates. However, Bayes factors are highly sensitive to the prior distribution of parameters. Conclusions on model choice based on Bayes factor can be misleading unless the sensitivity of conclusions to the choice of priors is carefully considered.

Model-based methods have been criticized for not exhaustively covering the hypothesis space.[24] Indeed, model-based studies often revolve around a small number of models, and due to the high computational cost to evaluate a single model in some instances, it may then be difficult to cover a large part of the hypothesis space.

An upper limit to the number of considered candidate models is typically set by the substantial effort required to define the models and to choose between many alternative options.[13] There is no commonly accepted ABC-specific procedure for model construction, so experience and prior knowledge are used instead.[14] Although more robust procedures for a priori model choice and formulation would be beneficial, there is no one-size-fits-all strategy for model development in statistics: sensible characterization of complex systems will always necessitate a great deal of detective work and use of expert knowledge from the problem domain.

Some opponents of ABC contend that since only few models—subjectively chosen and probably all wrong—can be realistically considered, ABC analyses provide only limited insight.[24] However, there is an important distinction between identifying a plausible null hypothesis and assessing the relative fit of alternative hypotheses.[12] Since useful null hypotheses, that potentially hold true, can extremely seldom be put forward in the context of complex models, predictive ability of statistical models as explanations of complex phenomena is far more important than the test of a statistical null hypothesis in this context. It is also common to average over the investigated models, weighted based on their relative plausibility, to infer model features (e.g., parameter values) and to make predictions.

Large data sets may constitute a computational bottleneck for model-based methods. It was, for example, pointed out that in some ABC-based analyses, part of the data have to be omitted.[24] A number of authors have argued that large data sets are not a practical limitation,[13][45] although the severity of this issue depends strongly on the characteristics of the models. Several aspects of a modeling problem can contribute to the computational complexity, such as the sample size, number of observed variables or features, time or spatial resolution, etc. However, with increasing computing power, this issue will potentially be less important.

Instead of sampling parameters for each simulation from the prior, it has been proposed alternatively to combine the Metropolis-Hastings algorithm with ABC, which was reported to result in a higher acceptance rate than for plain ABC.[36] Naturally, such an approach inherits the general burdens of MCMC methods, such as the difficulty to assess convergence, correlation among the samples from the posterior,[26] and relatively poor parallelizability.[13]

Likewise, the ideas of sequential Monte Carlo (SMC) and population Monte Carlo (PMC) methods have been adapted to the ABC setting.[26][48] The general idea is to iteratively approach the posterior from the prior through a sequence of target distributions. An advantage of such methods, compared to ABC-MCMC, is that the samples from the resulting posterior are independent. In addition, with sequential methods the tolerance levels must not be specified prior to the analysis, but are adjusted adaptively.[49]

It is relatively straightforward to parallelize a number of steps in ABC algorithms based on rejection sampling and sequential Monte Carlo methods. It has also been demonstrated that parallel algorithms may yield significant speedups for MCMC-based inference in phylogenetics,[50] which may be a tractable approach also for ABC-based methods. Yet an adequate model for a complex system is very likely to require intensive computation irrespectively of the chosen method of inference, and it is up to the user to select a method that is suitable for the particular application in question.

High-dimensional data sets and high-dimensional parameter spaces can require an extremely large number of parameter points to be simulated in ABC-based studies to obtain a reasonable level of accuracy for the posterior inferences. In such situations, the computational cost is severely increased and may in the worst case render the computational analysis intractable. These are examples of well-known phenomena, which are usually referred to with the umbrella term curse of dimensionality.[51]

To assess how severely the dimensionality of a data set affects the analysis within the context of ABC, analytical formulas have been derived for the error of the ABC estimators as functions of the dimension of the summary statistics.[52][53] In addition, Blum and François have investigated how the dimension of the summary statistics is related to the mean squared error for different correction adjustments to the error of ABC estimators. It was also argued that dimension reduction techniques are useful to avoid the curse-of-dimensionality, due to a potentially lower-dimensional underlying structure of summary statistics.[52] Motivated by minimizing the quadratic loss of ABC estimators, Fearnhead and Prangle have proposed a scheme to project (possibly high-dimensional) data into estimates of the parameter posterior means; these means, now having the same dimension as the parameters, are then used as summary statistics for ABC.[53]

ABC can be used to infer problems in high-dimensional parameter spaces, although one should account for the possibility of overfitting (e.g., see the model selection methods in [42] and [43]). However, the probability of accepting the simulated values for the parameters under a given tolerance with the ABC rejection algorithm typically decreases exponentially with increasing dimensionality of the parameter space (due to the global acceptance criterion).[14] Although no computational method (based on ABC or not) seems to be able to break the curse-of-dimensionality, methods have recently been developed to handle high-dimensional parameter spaces under certain assumptions (e.g., based on polynomial approximation on sparse grids,[54] which could potentially heavily reduce the simulation times for ABC). However, the applicability of such methods is problem dependent, and the difficulty of exploring parameter spaces should in general not be underestimated. For example, the introduction of deterministic global parameter estimation led to reports that the global optima obtained in several previous studies of low-dimensional problems were incorrect.[55] For certain problems, it might therefore be difficult to know whether the model is incorrect or, as discussed above, whether the explored region of the parameter space is inappropriate.[24] More pragmatic approaches are to cut the scope of the problem through model reduction,[14] discretisation of variables and the use of canonical models such as noisy models. Noisy models exploit information on the conditional independence between variables.[56]