# Lack of confidence in approximate Bayesian computation model choice

^{a}Université Paris-Dauphine, 75775 Paris cedex 16, France;^{b}Institut Universitaire de France, France;^{c}Centre de Recherche en Économie et Statistique (CREST), 92245 Malakoff cedex, France;^{d}Centre de Biologie pour la Gestion des Populations (CBGP), French National Institute for Agricultural Research (INRA), 34988 Montferrier-sur-Lez cedex, France;^{e}Unité Mixte de Recherche Centre National de la Recherche Scientifique (CNRS) 5149, Université Montpellier 2, 34095 Montpellier, France; and^{f}Department of Statistics, Harvard University, Cambridge, MA 02138-2901

See allHide authors and affiliations

Edited by Stephen E. Fienberg, Carnegie Mellon University, Pittsburgh, PA, and approved July 21, 2011 (received for review February 23, 2011)

## Abstract

Approximate Bayesian computation (ABC) have become an essential tool for the analysis of complex stochastic models. Grelaud et al. [(2009) *Bayesian Anal* 3:427–442] advocated the use of ABC for model choice in the specific case of Gibbs random fields, relying on an intermodel sufficiency property to show that the approximation was legitimate. We implemented ABC model choice in a wide range of phylogenetic models in the Do It Yourself-ABC (DIY-ABC) software [Cornuet et al. (2008) *Bioinformatics* 24:2713–2719]. We now present arguments as to why the theoretical arguments for ABC model choice are missing, because the algorithm involves an unknown loss of information induced by the use of insufficient summary statistics. The approximation error of the posterior probabilities of the models under comparison may thus be unrelated with the computational effort spent in running an ABC algorithm. We then conclude that additional empirical verifications of the performances of the ABC procedure as those available in DIY-ABC are necessary to conduct model choice.

Inference on population genetic models such as coalescent trees is one representative example of cases when statistical analyses such as Bayesian inference cannot easily operate because the likelihood function associated with the data cannot be computed in a manageable time (1⇓–3). The fundamental reason for this impossibility is that the model associated with coalescent data has to integrate over trees of high complexity.

In such settings, traditional approximation tools such as Monte Carlo simulation (4) from the posterior distribution are unavailable for practical purposes. Indeed, due to the complexity of the latent structures defining the likelihood (like the coalescent tree), their simulation is too unstable to bring a reliable approximation in a manageable time. Such complex models call for a practical if cruder approximation method, the approximate Bayesian computation (ABC) methodology (1, 5). This rejection technique bypasses the computation of the likelihood via simulations from the corresponding distribution (see refs. 6 and 7 for recent surveys, and ref. 8 for the wide and successful array of applications based on implementations of ABC in genomics and ecology).

We argue here that ABC is a generally valid approximation method for doing Bayesian inference in complex models. However, without further justification, ABC methods cannot be trusted to discriminate between two competing models when based on insufficient summary statistics. We exhibit simple examples in which the information loss due to insufficiency leads to inconsistency, i.e., when the ABC model selection fails to recover the true model, even with infinite amounts of observation and computation. On the one hand, ABC using the entire data leads to a consistent model-choice decision, but it is clearly infeasible in most settings. On the other hand, too much information loss due to insufficiency leads to a statistically invalid decision procedure. The challenge is in achieving a balance between information loss and consistency. Theoretical results that mathematically validate model choice for insufficient statistics are currently lacking on a general basis.

Our conclusion at this stage is to opt for a cautionary approach in ABC model choice, handling it as an exploratory tool rather than trusting the Bayes factor approximation. The corresponding degree of approximation cannot be evaluated, except via Monte Carlo evaluations of the model selection performances of ABC. More empirical measures such as those proposed in the DIY-ABC software (3) and in ref. 9 thus seem to be the only available solution at the current time for conducting model comparison.

We stress that, although refs. 10 and 11 repeatedly expressed reservations about the formal validity of the ABC approach in statistical testing, those criticisms were rebutted in refs. 12⇓–14 and are not relevant for the current paper.

## Statistical Methods

### The ABC Algorithm.

The setting in which ABC operates is the approximation of a simulation from the posterior distribution *π*(** θ**|

**y**) ∝

*π*(

*θ*)

*f*(

**y**|

**) when distributions associated with both the prior**

*θ**π*and the likelihood

*f*can be simulated (the latter being unavailable in closed form). The first ABC algorithm was introduced by ref. 5 as follows: Given a sample

**y**from a sample space , a sample (

*θ*

_{1},…,

*θ*

_{M}) is produced by

**Algorithm 1: ABC sampler**

**for** *i* = 1 to *N* **do**

**repeat**

Generate *θ*^{′} from the prior distribution *π*(·)

Generate **z** from the likelihood *f*(·|*θ*^{′})

**until** *ρ*{*η*(**z**),*η*(**y**)} ≤ *ϵ*

set *θ*_{i} = *θ*^{′},

**end for**

The parameters of the ABC algorithm are the so-called summary statistic *η*(·), the distance *ρ*{·,·}, and the tolerance level *ϵ* > 0. The approximation of the posterior distribution *π*(** θ**|

**y**) provided by the ABC sampler is to instead sample from the marginal in

**of the joint distribution where denotes the indicator function of**

*θ**B*and The basic justification of the ABC approximation is that, when using a sufficient statistic

*η*and a small (enough) tolerance

*ϵ*, we have

In practice, the statistic *η* is necessarily insufficient (because only exponential families enjoy sufficient statistics with fixed dimension; see ref. 15) and the approximation then converges to the less informative *π*(** θ**|

*η*(

**y**)) when

*ϵ*goes to zero. This loss of information is a necessary price to pay for the access to computable quantities and

*π*(

**|**

*θ**η*(

**y**)) provides a convergent inference on

**when**

*θ***is identifiable in the distribution of**

*θ**η*(

**y**) (16). While acknowledging the gain brought by ABC in handling Bayesian inference in complex models, and the existence of involved summary selection mechanisms (17, 18), we demonstrate here that the loss due to the ABC approximation may be arbitrary in the specific setting of Bayesian model choice via posterior model probabilities.

### ABC Model Choice.

The standard Bayesian tool for model comparison is the marginal likelihood (19) which leads to the Bayes factor for comparing the evidences of models with likelihoods *f*_{1}(**y**|*θ*_{1}) and *f*_{2}(**y**|*θ*_{2}), As detailed in ref. 12, it provides a valid criterion for model comparison that is naturally penalized for model complexity.

Bayesian model choice proceeds by creating a probability structure across *M* models (or likelihoods). It introduces the model index as an extra unknown parameter, associated with its prior distribution, (*m* = 1,…,*M*), whereas the prior distribution on the parameter is conditional on the value *m* of the index, denoted by *π*_{m}(*θ*_{m}) and defined on the parameter space Θ_{m}. The choice between those models is then driven by the posterior distribution of , where *w*_{k}(**y**) denotes the marginal likelihood for model *k*.

Although this posterior distribution is straightforward to interpret, it offers a challenging computational conundrum in Bayesian analysis. When the likelihood is not available, ABC represents the almost unique solution. Ref. 5 describes the use of model choice based on ABC for distinguishing between different mutation models. The justification behind the method is that the average ABC acceptance rate associated with a given model is proportional to the posterior probability corresponding to this approximative model, when identical summary statistics, distance, and tolerance level are used over all models. In practice, an estimate of the ratio of marginal likelihoods is given by the ratio of observed acceptance rates. Using Bayes formula, estimates of the posterior probabilities are straightforward to derive. This approach has been widely implemented in the literature (see, e.g., refs. 20⇓⇓–23).

A representative illustration of the use of an ABC model-choice approach is given by ref. 21, which analyses the European invasion of the Western corn rootworm, North America’s most destructive corn pest. Because this pest was initially introduced in Central Europe, it was believed that subsequent outbreaks in Western Europe originated from this area. Based on an ABC model-choice analysis of the genetic variability of the rootworm, the authors conclude that this belief is false: There have been at least three independent introductions from North America during the past two decades.

The above estimate is improved by regression regularization (24), where model indices are processed as categorical variables in a polychotomous regression. When comparing two models, this involves a standard logistic regression. Rejection-based approaches were lately introduced by refs. 3, 25, and 26, in a Monte Carlo simulation of model indices as well as model parameters. Those recent extensions are already widely used in population genetics, as exemplified by refs. 27⇓⇓⇓⇓⇓⇓⇓⇓–36. Another illustration of the popularity of this approach is given by the availability of four softwares implementing ABC model-choice methodologies:

ABC-SysBio, which relies on a sequential Monte Carlo (SMC)-based ABC for inference in system biology, including model choice (26).

ABCToolbox, which proposes SMC and Markov chain Monte Carlo implementations, as well as Bayes factor approximation (37).

DIY-ABC, which relies on a regularized ABC model choice (ABC-MC) algorithm on population history using molecular markers (3).

PopABC, which relies on a regular ABC-MC algorithm for genealogical simulation (38).

As exposed in, e.g., refs. 25, 39, or 40, once is incorporated within the parameters, the ABC approximation to its posterior follows from the same principles as in regular ABC. The corresponding implementation is as follows, using for the summary statistic a statistic ** η**(

**z**) = {

*η*

_{1}(

**z**),…,

*η*

_{M}(

**z**)} that is the concatenation of the summary statistics used for all models (with an obvious elimination of duplicates):

**Algorithm 2: ABC-MC**

**for** *i* = 1 to *N* **do**

**repeat**

Generate *m* from the prior

Generate *θ*_{m} from the prior *π*_{m}(*θ*_{m})

Generate **z** from the model *f*_{m}(**z**|*θ*_{m})

**until** *ρ*{** η**(

**z**),

**(**

*η***y**)} ≤

*ϵ*

Set *m*^{(i)} = *m* and *θ*^{(i)} = *θ*_{m}

**end for**

The ABC estimate of the posterior probability is then the frequency of acceptances from model *m* in the above simulation . This also corresponds to the frequency of simulated pseudodatasets from model *m* that are closer to the data **y** than the tolerance *ϵ*. In order to improve the estimation by smoothing, ref. 3 follows the rationale that motivated the use of a local linear regression in ref. 2 and relies on a weighted polychotomous regression to estimate based on the ABC output. This modeling is implemented in the DIY-ABC software.

### The Difficulty with ABC-MC.

There is a fundamental discrepancy between the genuine Bayes factors/posterior probabilities and the approximations resulting from ABC-MC.

The ABC approximation to a Bayes factor, *B*_{12} say, resulting from Algorithm 2, is An alternative representation is given by where the pairs (*m*^{t},*z*^{t}) are simulated from the joint prior and *T* is the number of simulations necessary for *N* acceptances in Algorithm 2. In order to study the limit of this approximation, we first let *T* go to infinity. (For simplification purposes and without loss of generality, we choose a uniform prior on the model index.) The limit of is then where and denote the densities of ** η**(

**z**) when

**z**∼

*f*

_{1}(

**z**|

*θ*_{1}) and

**z**∼

*f*

_{2}(

**z**|

*θ*_{2}), respectively. By L’Hospital formula, if

*ϵ*goes to zero, the above converges to namely the Bayes factor for testing model 1 versus model 2 based on the sole observation of

**(**

*η***y**). This result reflects the current perspective on ABC: The inference derived from the ideal ABC output when

*ϵ*= 0 uses only the information contained in

**(**

*η***y**). Thus, in the limiting case, i.e., when the algorithm uses an infinite computational power, the ABC odds ratio does not account for features of the data other than the value of

**(**

*η***y**), which is why the limiting Bayes factor depends only on the distribution of

**under both models.**

*η*When running ABC for point estimation, the use of an insufficient statistic does not usually jeopardize convergence of the method. As shown, e.g., in ref. 16, Theorem 2, the noisy version of ABC as an inference method is convergent under usual regularity conditions for model-based Bayesian inference (41), including identifiability of the parameter for the insufficient statistic ** η**. In contrast, the loss of information induced by

**may seriously impact model-choice Bayesian inference. Indeed, the information contained in**

*η***(**

*η***y**) is lesser than the information contained in

**y**and this even in most cases when

**(**

*η***y**) is a sufficient statistic for

*both models*. In other words,

**(**

*η***y**)

*being sufficient for both*

*f*

_{1}(

**y**|

*θ*_{1})

*and*

*f*

_{2}(

**y**|

*θ*_{2})

*does not usually imply that*

**(**

*η***y**)

*is sufficient for*{

*m*,

*f*

_{m}(

**y**|

*θ*_{m})}. To see why this is the case, consider the most favorable case, namely when

**(**

*η***y**) is a sufficient statistic for both models. We then have by the factorization theorem (15) that (

*i*= 1,2); i.e., [1]Thus, unless

*g*

_{1}(

**y**) =

*g*

_{2}(

**y**), as in the special case of Gibbs random fields detailed below, the two Bayes factors differ by the ratio

*g*

_{1}(

**y**)/

*g*

_{2}(

**y**), which is equal to one only in a very small number of known cases. This decomposition is a straightforward proof that a modelwise sufficient statistic is usually not sufficient across models, hence for model comparison. An immediate corollary is that the ABC-MC approximation does not always converge to the exact Bayes factor.

The discrepancy between limiting ABC and genuine Bayesian inferences does not come as a surprise, because ABC is indeed an approximation method. Users of ABC algorithms are therefore prepared for some degree of imprecision in their final answer, a point stressed by refs. 16 and 42 when they qualify ABC as exact inference on a wrong model. However, the magnitude of the difference between *B*_{12}(**y**) and expressed by Eq. **1** is such that there is no direct connection between both answers. In a general setting, if ** η** has the same dimension as one component of the

*n*components of

**y**, the ratio

*g*

_{1}(

**y**)/

*g*

_{2}(

**y**) is equivalent to a density ratio for a sample of size O(

*n*); hence it can be arbitrarily small or arbitrarily large when

*n*grows. Contrastingly, the Bayes factor is based on an equivalent to a single observation and hence does not necessarily converge with

*n*to the correct limit, as shown by the Poisson and normal examples below and in

*SI Text*. The conclusion derived from the ABC-based Bayes factor may therefore completely differ from the conclusion derived from the exact Bayes factor and there is no possibility of a generic agreement between both, or even of a manageable correction factor. This discrepancy means that a theoretical validation of the ABC-based model-choice procedure is currently missing and that, due to this absence, potentially costly simulation-based assessments are required when calling for this procedure.

Therefore, users must be warned that ABC approximations to Bayes factors do not perform as standard numerical or Monte Carlo approximations, with the exception of Gibbs random fields detailed in the next section. In all cases when *g*_{1}(**y**)/*g*_{2}(**y**) differs from one, no inference on the true Bayes factor can be derived from the ABC-MC approximation without further information on the ratio *g*_{1}(**y**)/*g*_{2}(**y**), most often unavailable in settings where ABC is necessary.

Ref. 40 also derived this relation between both Bayes factors in their formula 18. Although they still advocate the use of ABC model choice in the absence of sufficient statistic, we stress that no theoretical guarantee can be given on the validity of the ABC approximation to the Bayes factor and hence of its use as a model-choice procedure.

Note that the authors of ref. 43 resort to full allelic distributions in an ABC framework, instead of choosing summary statistics. They show how to apply ABC using allele frequencies to draw inferences in cases where selecting suitable summary statistics is difficult (and where the complexity of the model or the size of dataset prohibits the use of full-likelihood methods). In such settings, ABC-MC does not suffer from the divergence exhibited here because the measure of distance does not involve a reduction of the sample. The same comment applies to the ABC-SysBio software of ref. 26, which relies on the whole dataset. The theoretical validation of ABC inference in hidden Markov models by ref. 44 should also extend to the model-choice setting because the approach does not rely on summary statistics but instead on the whole sequence of observations.

## Results

### The Specific Case of Gibbs Random Fields.

In an apparent contradiction with the above, ref. 25 showed that the computation of the posterior probabilities of Gibbs random fields under competition can be done via ABC techniques, which provide a converging approximation to the true Bayes factor. The reason for this result is that, for these models in the above ratio Eq. **1**, *g*_{1}(**y**) = *g*_{2}(**y**). The validation of an ABC comparison of Gibbs random fields is thus that their specific structure allows for a sufficient statistic vector that runs across models and therefore leads to an exact (when *ϵ* = 0) simulation from the posterior probabilities of the models. Each Gibbs random field model has its own sufficient statistic *η*_{m}(·), and ref. 25 exposed the fact that the vector of statistics ** η**(·) = (

*η*

_{1}(·),…,

*η*

_{M}(·)) is also sufficient for the joint parameter .

The authors of ref. 40 point out that this specific property of Gibbs random fields can be extended to any exponential family (hence to any setting with fixed-dimension sufficient statistics; see ref. 15). Their argument is that, by including all sufficient statistics and all dominating measure statistics in an encompassing model, models under comparison are submodels of the encompassing model. The concatenation of those statistics is then jointly sufficient across models. Although this encompassing principle holds in full generality, in particular when comparing models that are already embedded, we think it leads to an overly optimistic perspective about the merits of ABC for model choice: In practice, most complex models do not enjoy sufficient statistics (if only because they are beyond exponential families). The Gibbs case processed by ref. 25 therefore happens to be one of the very few realistic counterexamples. As demonstrated in the next section and in the normal example in *SI Text*, using insufficient statistics is more than a mere loss of information. Looking at what happens in the limiting case when one relies on a common modelwise sufficient statistic is a formal but useful study because it sheds light on the potentially huge discrepancy between the ABC-based and the true Bayes factors. To develop a solution to the problem in the formal case of the exponential families does not help in understanding the discrepancy for nonexponential models.

### Arbitrary Ratios.

The difficulty with the discrepancy between *B*_{12}(**y**) and is that this discrepancy is impossible to evaluate in a general setting, whereas there is no reason to expect a reasonable agreement between both quantities. A first illustration was produced by ref. 45 in the case of MA(*q*) models. A second illustration is detailed in *SI Text* for the normal mode; see Fig. S1.

A simple illustration of the discrepancy due to the use of a modelwise sufficient statistic is a sample **y** = (*y*_{1},…,*y*_{n}) that could come either from a Poisson distribution or from a geometric distribution, already introduced in ref. 25 as a counterexample to Gibbs random fields and later reprocessed in ref. 40 to support their sufficiency argument. In this case, the sum is a sufficient statistic for both models but not across models. The distribution of the sample given *S* is a multinomial distribution when the data are Poisson, whereas it is the uniform distribution over the **y**’s such that in the geometric case, because *S* is then a negative binomial variable. The discrepancy ratio is therefore When simulating *n* Poisson or geometric variables and using prior distributions as exponential and uniform on the parameters of the respective models, the exact Bayes factor is available and the distribution of the discrepancy is therefore available. Fig. S2 gives the range of *B*_{12}(**y**) versus , showing that is then unrelated with *B*_{12}(**y**): The values produced by both approaches have nothing in common. As noted above, the approximation based on the sufficient statistic *S* is producing figures of the magnitude of a *single* observation, whereas the true Bayes factor is of the order of the sample size.

The discrepancy between both Bayes factors is in fact increasing with the sample size, as shown by the following result:

Consider model selection between model 1: with prior distribution *π*_{1}(*λ*) equal to an exponential distribution and model 2: with a uniform prior distribution *π*_{2} when the observed data **y** consists of independent observations with expectation . Then is the minimal sufficient statistic for both models and the Bayes factor based on the sufficient statistic *S*(**y**), , satisfies

Therefore, the Bayes factor based on the statistic *S*(**y**) is *not* consistent; it converges to a nonzero, finite value.

In this specific setting, ref. 40 shows that adding to the *S* creates a statistic (*S*,*P*) that is sufficient across both models. Although this is mathematically correct, it does not provide further understanding of the behavior of ABC model choice in realistic settings: Outside formal examples such as the one above and well-structured if complex exponential families such as Gibbs random fields, it is not possible to devise completion mechanisms that ensure sufficiency across models, or even select well-discriminating statistics. It is therefore more fruitful to study and detect the diverging behavior of the ABC approximation as given, rather than attempting at solving the problem in a specific and formal case.

### Population Genetics.

We recall that ABC has first been introduced by population geneticists (2, 5) for statistical inference about the evolutionary history of species, because no likelihood-based approach existed apart from very simple and hence unrealistic situations. This approach has since been used in an increasing number of biological studies (20, 24, 46), most of them including model choice. It is therefore crucial to get insights into the validity of such studies, particularly when they deal with species of economical or ecological importance (see, e.g., ref. 47). To this end, we need to compare ABC estimates of posterior probabilities to reliable likelihood-based estimates. Combining different modules based on ref. 48, it is possible to approximate the likelihood of population genetic data through importance sampling (IS) even in complex scenarios. In order to evaluate the potential discrepancy between ABC-based and likelihood-based posterior probabilities of evolutionary scenarios, we designed two experiments based on simulated data with limited information content, so that the choice between scenarios is problematic. This setting thus provides a wide enough set of intermediate values of model posterior probabilities, in order to better evaluate the divergence between ABC and likelihood estimates.

In the first experiment, we consider two populations (1 and 2) having diverged at a fixed time in the past and a third population (3) having diverged from one of those two populations (scenarios 1 and 2, respectively). Times are set to 60 generations for the first divergence and to 30 generations for the second divergence. One hundred pseudo observed datasets have been simulated, represented by 15 diploid individuals per population genotyped at five independent microsatellite loci. These loci are assumed to evolve according to the strict stepwise mutation model; i.e., when a mutation occurs, the number of repeats of the mutated gene increases or decreases by one unit with equal probability. The mutation rate, common to all five loci, has been set to 0.005 and effective population sizes to 30. In this experiment, both scenarios have a single parameter: the effective population size, assumed to be identical for all three populations. We chose a uniform prior *U*[2,150] for this parameter (the true value being 30). The IS algorithm was performed using 100 coalescent trees per particle. The marginal likelihood of both scenarios has been computed for the same set of 1,000 particles, and they provide the posterior probability of each scenario. The ABC computations have been performed with DIY-ABC (3). A reference table of 2 million datasets has been simulated using 24 usual summary statistics (provided in Table S1) and the posterior probability of each scenario has been estimated as their proportion in the 500 simulated datasets closest to the pseudo observed one. This population genetic setting does not allow for a choice of sufficient statistics, even at the model level.

The second experiment also opposes two scenarios including three populations, two of them having diverged 100 generations ago and the third one resulting from a recent admixture between the first two populations (scenario 1) or simply diverging from population 1 (scenario 2) at the same time of 5 generations in the past. In scenario 1, the admixture rate is 0.7 from population 1. Pseudo observed datasets (100) of the same size as in experiment 1 (15 diploid individuals per population, 5 independent microsatellite loci) have been generated for an effective population size of 1,000 and mutation rates of 0.0005. In contrast with experiment 1, analyses included the following six parameters (provided with corresponding priors): admixture rate (*U*[0.1,0.9]), three effective population sizes (*U*[200,2000]), the time of admixture/second divergence (*U*[1,10]), and the time of the first divergence (*U*[50,500]). To account for a higher complexity in the scenarios, the IS algorithm was performed with 10,000 coalescent trees per particle. Apart from this change, both ABC and likelihood analyses have been performed in the same way as experiment 1.

Fig. 1 shows a reasonable fit between the exact posterior probability of model 1 (evaluated by IS) and the ABC approximation in the first experiment on most of the 100 simulated datasets, even though the ABC approximation is biased toward 0.5. When using 0.5 as the decision boundary between model 1 and model 2, there is hardly any discrepancy between both approaches, demonstrating that model choice based on ABC can be trusted in this case. Fig. 2 considers the same setting when moving from 24 to 15 summary statistics (given in Table S1): The fit somehow degrades. In particular, the number of opposite conclusions in the model choice moves to 12%. In the more complex setting of the second experiment, the discrepancy worsens, as shown in Fig. 3. The number of opposite conclusions reaches 26% and the fit between both versions of the posterior probabilities is considerably degraded, with a correlation coefficient of 0.643.

The validity of the importance sampling approximation can obviously be questioned in both experiments; however, Figs. S3 and S4 display a strong stability of the posterior probability IS approximation across 10 independent runs for 5 different datasets and give proper confidence in this approach. Increasing the number of loci to 50 and the sample size to 100 individuals per population (see *SI Text*) leads to posterior probabilities of the true scenario overwhelmingly close to one (Fig. S5), thus blurring the distinction between ABC and likelihood-based estimates but also reassuring the ability of ABC to provide the right choice of model with a higher information content of the data. We note that, for this experiment, all ABC-based decisions conclude in favor of the correct model. As shown in Fig. S6, this second experiment requires an increase in the number of importance sampling simulations because of a higher variability in the likelihood.

## Discussion

Since its introduction by refs. 1 and 5, ABC has been extensively used in areas involving complex likelihoods, both for point estimation and testing of hypotheses. In realistic settings, with the exception of models such as Gibbs random fields, which are resilient with respect to their sufficient statistics, the conclusions drawn on model comparison cannot be trusted per se but require further simulation analyses as to the pertinence of the (ABC) Bayes factor based on the summary statistics. This paper has examined in details only the case when the summary statistics are sufficient for both models, while practical situations imply the use of insufficient statistics. The rapidly increasing number of applications estimating posterior probabilities by ABC indicates a clear need for further evaluations of the worth of those estimations, especially because our population genetic experiments showed that those ABC approximations were selecting the proper model.

Further research is needed for producing trustworthy approximations to the posterior probabilities of models. At this stage, unless the whole data are involved in the ABC approximation as in ref. 43, our conclusion on ABC-based model choice is to exploit the approximations in an exploratory manner as measures of discrepancies rather than genuine posterior probabilities. This direction relates with the analyses found in ref. 9. Furthermore, a version of this exploratory analysis is already provided in the DIY-ABC software of ref. 3. An option in this software allows for the computation of a Monte Carlo evaluation of false allocation rates resulting from using the ABC posterior probabilities in selecting a model as the most likely. For instance, in the setting of both our population genetic experiments, DIY-ABC gives false allocation rates equal to 20% (under scenarios 1 and 2) and 14.5% and 12.5% (under scenarios 1 and 2), respectively. This evaluation obviously shifts away from the performances of ABC as an approximation to the posterior probability toward the performances of the whole Bayesian apparatus for selecting a model, but this nonetheless represents a useful and manageable quality assessment for practitioners.

## Acknowledgments

The authors are grateful to the reviewers and to Michael Stumpf for their comments. Computations were performed on the INRA CBGP and MIGALE clusters. C.P.R., J.-M.C., and J-M.M. have been partly supported by Agence Nationale de la Recherche via the 2009–2013 project EMILE (Études de Méthodes Inférentielles et Logiciels pour l’Évolution). N.S.P. gratefully acknowledges National Science Foundation Grant DMS 1107070.

## Footnotes

- ↵
^{1}To whom correspondence should be addressed. E-mail: Christian.Robert{at}ceremade.dauphine.fr.

Author contributions: C.P.R., J.-M.M., and N.S.P. designed research; C.P.R., J.-M.M., and N.S.P. performed research; J.-M.C. and J.-M.M. analyzed data; and C.P.R., J.-M.C., and J.-M.M. wrote the paper.

The authors declare no conflict of interest.

This article is a PNAS Direct Submission.

This article contains supporting information online at www.pnas.org/lookup/suppl/doi:10.1073/pnas.1102900108/-/DCSupplemental.

Freely available online through the PNAS open access option.

## References

- ↵
- ↵
- ↵
- Cornuet JM,
- et al.

- ↵
- Robert C,
- Casella G

- ↵
- ↵
- ↵
- ↵
- ↵
- Ratmann O,
- Andrieu C,
- Wiujf C,
- Richardson S

- ↵
- ↵
- Templeton A

- ↵
- ↵
- ↵
- Berger J,
- Fienberg S,
- Raftery A,
- Robert C

- ↵
- Lehmann E,
- Casella G

- ↵
- Fearnhead P,
- Prangle D

- ↵
- Joyce P,
- Marjoram P

- ↵
- Nunes MA,
- Balding DJ

- ↵
- Jeffreys H

- ↵
- ↵
- Miller N,
- et al.

- ↵
- ↵
- ↵
- Fagundes N,
- et al.

- ↵
- Grelaud A,
- Marin JM,
- Robert C,
- Rodolphe F,
- Tally F

- ↵
- ↵
- ↵
- ↵
- Excoffier C,
- Leuenberger D,
- Wegmann L

- ↵
- Ghirotto S,
- et al.

- ↵
- Guillemaud T,
- Beaumont M,
- Ciosi M,
- Cornuet JM,
- Estoup A

- ↵
- ↵
- ↵
- ↵
- ↵
- Wegmann D,
- Excoffier L

- ↵
- ↵
- Lopes JS,
- Balding D,
- Beaumont MA

- ↵
- Toni T,
- Stumpf M

- ↵
- ↵
- Bernardo J,
- Smith A

- ↵
- Wilkinson RD

- ↵
- ↵
- Dean T,
- Singh S,
- Jasra A,
- Peters G

- ↵
- Marin J,
- Pudlo P,
- Robert C,
- Ryder R

- ↵
- ↵
- ↵
- Stephens D,
- Donnelly P

## Citation Manager Formats

## Article Classifications

- Physical Sciences
- Statistics

- Biological Sciences
- Biophysics and Computational Biology