# Goodness of Fit for Bayesian Generative Models with Applications in Population Genetics

Guillaume Le Mailloux, Paul Bastide, Jean-Michel Marin, Arnaud Estoup

<figure><img src="image_1.png" /></figure>

## To cite this version:

Guillaume Le Mailloux, Paul Bastide, Jean-Michel Marin, Arnaud Estoup. Goodness of Fit for Bayesian Generative Models with Applications in Population Genetics. 2025. hal-05344254

HAL Id: hal-05344254

https://hal.science/hal-05344254v1

Preprint submitted on 3 Nov 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Goodness of Fit for Bayesian Generative Models with Applications in Population Genetics*

Guillaume Le Mailloux¹, Paul Bastide², Jean-Michel Marin†¹, and Arnaud Estoup³

¹IMAG, Univ. Montpellier, CNRS, Montpellier, France

²MAP5, Univ. Paris Cité, CNRS, Paris, France

³CBGP, INRAE, CIRAD, IRD, Montpellier SupAgro, Univ. Montpellier, Montpellier, France

January 28, 2025

**Abstract**

In population genetics and other application fields, models with intractable likelihood are common. Approximate Bayesian Computation (ABC) or more generally Simulation-Based Inference (SBI) methods work by simulating instrumental data sets from the models under study and comparing them with the observed data set $y_{obs}$, using advanced machine learning tools for tasks such as model selection and parameter inference. The present work focuses on model criticism, and more specifically on Goodness of fit (GoF) tests, for intractable likelihood models. We introduce two new GoF tests: the pre-inference GoF tests whether $y_{obs}$ is distributed from the prior predictive distribution, while the post-inference GoF tests whether there is a parameter value such that $y_{obs}$ is distributed from the likelihood with that value. The pre-inference test can be used to prune a large set of models using a limited amount of simulations, while the post-inference test is used to assess the fit of a selected model. Both tests are based on the local outlier factor (LOF, Breunig et al., 2000). This indicator was initially defined for outlier and novelty detection. It is able to quantify local density deviations, capturing subtleties that a more traditional kNN-based approach may miss. We evaluated the performance of our two GoF tests on simulated datasets from three different model settings of varying complexity. We then illustrate the utility of these approaches on a dataset of single nucleotide polymorphism (SNP) markers for the evaluation of complex evolutionary scenarios of modern human populations. Our dual-test GoF approach highlights the flexibility of our method: the pre-inference GoF test provides insight into model validity from a Bayesian perspective, while the post-inference test provides a more general and traditional view of assessing goodness of fit.

**Keywords—** generative models, goodness-of-fit, statistical tests, novelty detection, likelihood-free, population genomics, Single Nucleotide Polymorphism

<footnote>*All authors contributed equally to this work.</footnote>

<footnote>†Corresponding author: jean-michel.marin@umontpellier.fr</footnote>
# 1 Introduction

In numerous applications, the direct calculation of likelihood functions is unfeasible, thereby posing considerable challenges for inference. Two common cases serve to exemplify this issue. The first case pertains to latent variable models, wherein the likelihood function necessitates the evaluation of an intractable integral: $f(\mathbf{y}|\theta) = \int f(\mathbf{y}, \mathbf{u}|\theta)\mu(\mathbf{du})$. This is particularly challenging due to the hidden structure of the model. The second case pertains to likelihoods that are expressed as follows: $f(\mathbf{y}|\theta) = g(\mathbf{y}, \theta)/Z(\theta)$ and $Z(\theta)$ intractable. Sophisticated inferential techniques are required in both cases. In this paper, we will focus on cases where it is possible to simulate from the likelihood, which is typically more straightforward in the case of a latent variable model.

## 1.1 Statistical Inference with Intractable Likelihood

Traditional solutions to these scenarios include expectation-maximisation (EM) algorithms, Gibbs sampling, pseudo-marginal Markov chain Monte Carlo (MCMC) methods and variational approximations. However, these methods rely on the explicit computation of the likelihood $f(\mathbf{y}, \mathbf{u}|\theta)$. In Population Genetics and many other application fields, complex models with intractable likelihoods are common. Approximate Bayesian Computation (ABC; Marin et al., 2012) has emerged as a significant methodology for inference when the likelihood is intractable but data generation from the model is feasible. ABC, also known as Simulation-Based Inference (SBI) in the context of machine learning [Cranmer et al., 2020], circumvents the direct calculation of the likelihood by comparing simulated pseudo-datasets to the observed dataset, thereby facilitating inference without the explicit computation of the likelihood.

Since the seminal contributions of Tavaré et al. [1997] and Pritchard et al. [1999], ABC/SBI has undergone substantial developments. Most popular and effective strategies involve applying Machine Learning tools to a training set generated by a Bayesian generative model (the reference table). We review below some of the recent developments in the field, with a focus on applications to Population Genetics, and refer to Marin et al. [2012] and Cranmer et al. [2020] for comprehensive reviews, and to Lueckmann et al. [2021] for an extensive benchmark of simulation-based inference methods. Recently, Papamakarios and Murray [2016] introduced a method for fast likelihood-free inference using Bayesian Conditional Density Estimation with Mixture Density Networks, which has been shown to be effective in estimating the posterior distribution, particularly in the context of complex models. Pudlo et al. [2016] put forth an approach for dependable ABC model selection utilizing classification random forests, markedly enhancing model choice accuracy by capitalizing on machine learning techniques in a bioinformatics context. Moreover, Sheehan and Song [2016] employed deep learning techniques, specifically unsupervised pre-training via autoencoders, in the context of population genetic inference, thereby illustrating the potential of neural networks for the extraction of informative features in the context of ABC. Building on these advancements, Raynal et al. [2018] developed ABC random forests for Bayesian parameter inference, employing regression random forests to estimate parameters within the ABC framework.

In a more recent development, Papamakarios et al. [2019] introduced Sequential Neural Likelihood Estimation, utilizing autoregressive flows for neural density estimation of the likelihood. This enables fast and accurate inference with normalizing flows, including Masked Autoregressive Flows [Papamakarios et al., 2017]. Greenberg et al. [2019] introduced a method called automatic posterior transformation for likelihood-free inference, which also uses neural density estimation. This method is different from neural likelihood estimation because it uses normalizing flows to directly estimate the posterior. This method is part of the class of Sequential Neural Posterior Estimation, like Lueckmann et al. [2017]. Another notable
contribution was made by Hermans et al. [2020], who introduced Likelihood-free MCMC with approximate likelihood ratios through Sequential Neural Ratio Estimation. This method employs neural networks to estimate likelihood ratios, thereby providing a novel approach to efficient likelihood-free inference.

## 1.2 Model Criticism

The present work focuses on model criticism, which is essential for the assessment of whether an observed dataset, $y_{obs}$, is likely to have been generated by a given model. In the context of likelihood-free inference, Ratmann et al. [2009] discuss the benefit of incorporating model diagnostics within an ABC framework, and demonstrate how this method diagnoses model mismatch and guides model refinement. Wilkinson [2013] shown that ABC gives exact results under the assumption of model error. More recently, Frazier et al. [2020] explore the effects of model misspecification in ABC, showing that different ABC methods produce varying results when the simulation model deviates from the true data-generating process. It finds that accept-reject ABC concentrates on a pseudotrue value but lacks valid frequentist coverage, while local regression ABC targets a different pseudotrue value. The study also proposes diagnostic methods to detect model misspecification, supported by theoretical analysis and practical examples. Finally, Ward et al. [2022] discuss the issue of model misspecification in Neural Posterior Estimation, highlighting how traditional approaches often assume that simulators are perfectly specified representations of reality. It is demonstrated that such assumptions result in unreliable inferences when discrepancies exist between simulated and observed data. To address this issue, they propose Robust Neural Posterior Estimation, a methodology that explicitly models these discrepancies, thereby enabling both model criticism and robust inference. This approach ultimately mitigates the risks associated with misspecification and produces more reliable posterior estimates.

The present work focuses on Bayesian model checking. This type of analysis is important for model validation. It can also be used for model pruning when the number of candidates to be compared is excessive, especially in the context where data simulation is expensive. Bayesian model checking entails the assessment of the degree of fit between a Bayesian model and the observed data. As discussed by Box [1980], the prior predictive check assesses model assumptions without consideration of the data. It involves a comparison of observed data with data simulated from the prior predictive distribution (the likelihood is integrated over the prior, also called the evidence), with the aim of identifying major discrepancies between the simulated data and the observed dataset. In the posterior predictive check, as outlined by Guttman [1967], Rubin [1984], and Gelman et al. [1996], observed data is compared to data simulated from the posterior predictive distribution (the likelihood is integrated over the posterior distribution). More recent approaches, such as the holdout method proposed by Moran et al. [2023], utilize a portion of the data as a validation set to evaluate model performance on unseen data. Lemaire et al. [2016] presented a test for goodness of fit in the prior predictive check framework, that makes use of k-nearest neighbour (k-NN) density estimation. This approach circumvents the necessity for explicit likelihood functions by estimating the density of observed and simulated data to measure misfit. In the present work, we focus on (i) improving the model misfit measure by drawing from the novelty detection literature and (ii) building a correctly calibrated post-inference goodness of fit test using the holdout approach.

## 1.3 Goodness of Fit as a Novelty Detection Problem

In this paper, we link the question of likelihood-free goodness of fit to the field of novelty and outlier detection in order to improve on the results obtained by Lemaire et al. [2016]. Novelty/outlier detection
is used to identify anomalous or infrequent data points within a given dataset. Novelty detection is the
process of identifying previously unseen patterns in new data, whereas outlier detection is the method of
locating data points that deviate significantly from the majority in an existing dataset. These techniques are
of particular importance in fields such as anomaly detection, fraud detection, and fault monitoring. Goldstein
and Uchida [2016] discuss various unsupervised outlier detection techniques, emphasizing their scalability
and adaptability to high-dimensional data. Han et al. [2022] explore more recent advancements in anomaly
detection methods, particularly in the context of machine learning.

One frequently employed method for novelty/outlier detection is the Local Outlier Factor (LOF; Breunig et al. [2000]. The LOF identifies anomalous data points based on their local density relative to their neighbors, thereby enabling the effective detection of outliers in datasets where clusters of varying point densities are present. This method assesses the degree of isolation of a given point in comparison to its surrounding context, rendering it particularly useful for datasets exhibiting diverse density structures, where more global methods may prove ineffective. Figure 1 illustrates this behavior. It shows the first principal component axes of summary statistics of data points simulated according to a simple Population Genetics model (namely, the scenario 2 of the so-called Dep-Indep setting detailed below). The simulated distribution of points (blue circles) has an heterogeneous density, with a densely populated region at the center, and a sparser region on the lower left corner of the plot. The first test point (black triangle) is located at the border of the dense region, so that its distance to its neighbors as well as its local density are typical of the simulated points around. The prior GoF test described in this article does not reject the null hypothesis, whether using the kNN or LOF based score. The second test point (black square) is isolated, so that its local density is lower than that of its closest simulated points: in this case the prior GoF test based on the LOF score identifies this local density discrepancy and rejects the null hypothesis. On the other hand, the GoF test based on the kNN score does not reject the null hypothesis in this case again, because the raw distance of the second test point to its neighbors is not atypical among the simulated points at the level of the dataset as a whole due to the presence of sparse regions, especially in the bottom left area. In what follows, we will show that this behavior extends consistently to large and complex models, which leads to an improvement in the power of the test based on the LOF score.

1.4 Statistical Inference in Population Genetics

A significant number of ABC/SBI methods have been developed by population geneticists [Tavaré et al., 1997, Pritchard et al., 1999, Beaumont et al., 2002]. These methods are particularly valuable in this area due to the high complexity of demo-genetic models and the often intractable or implicit nature of the likelihood functions associated with these models. ABC/SBI permit researchers to test evolutionary hypotheses through model choice and estimate parameters even when traditional likelihood-based approaches are infeasible. For example, applications include the inference of demographic history, selection, migration rates, and recombination patterns. By focusing on summary statistics that capture essential key patterns in the data, ABC/SBI methods allow for scalable and flexible analyses, thus making them essential tools in modern population genomics [Beaumont, 2010, 2019].

To illustrate how powerful our likelihood-free goodness of fit test can be in the field, we worked on a
case study (relying both on simulated and real genomic datsets) involving complex evolutionary scenarios
among human populations Pudlo et al. [2016]. We started with six hypothetical evolutionary scenarios
to model four modern human populations. We chose two of these scenarios, which are quite similar but
really important in terms of real human SNP data, to look at in more detail. This setup is a great way to test
how well modern computational approaches can distinguish between different evolutionary paths, which are
often complex and related, and to assess to what extent do the observed data fit with the scenario considered
to be the best among those proposed.

**1.5 Outline**

The paper is structured as follows. In section 2, we outline our method, detailing the computational frame-
work. In section 3, we test our method through simulation studies, applying it to a range of hypothetical
evolutionary scenarios to evaluate its accuracy and robustness. Section 4 then illustrates our approach using
a real SNP dataset from modern human populations, providing empirical evidence of its practical utility. A
discussion concludes the paper.

<figure><img src="image_3.png" /><figcaption>Figure 1: **Illustrative example of different or similar conclusions when using tests based on the LOF or the kNN score from the Dep-Indep setting** (see Section 3.1.2). Axes 1 and 3 from a PCA on the summary statistics are represented. Simulated points (blue) are from scenario 2, while test points (black) are from scenario 1 of the Dep-Indep setting. The LOF score is computed with the “max-LOF” heuristic with *k* varying between 5 and 20, while the kNN uses *k* = 1 (see Section 2.2 Outlier Score). The square point (top) is correctly rejected by the LOF, but not rejected by the kNN. It is far from the other points locally, so that it has a high LOF score, but its distance to nearby points is not an outlier when compared to typical distances in the whole dataset (particularly in the bottom left region), so its kNN score is not outlying. The triangle point is not rejected by either the LOF or kNN scores: its distance from other points does not differ from that of its neighbors.</figcaption></figure>
# 2 Methods

## 2.1 Model

Let $\Theta_m$ be the space of possible parameter values for a given model $m$, with a prior distribution $\pi_m$ and likelihood function $f_m(\mathbf{y} | \theta)$ that gives a distribution on the observation space $\mathcal{Y}$ given a vector of parameters $\theta \in \Theta_m$. We are typically interested in cases where this likelihood is intractable, and we assume that the model is approximated by simulated data recorded in a reference table $\mathbf{D}_m$:

$$ \mathbf{D}_m = \{(\theta_i, \mathbf{y}_i), i \in [1, N_{\text{ref}}]\}, \text{ with } \theta_i \sim \pi_m \text{ and } \mathbf{y}_i \sim f_m(\cdot | \theta_i). $$

Typically, the observation space $\mathcal{Y}$ has a very high dimension, and a set of summary statistics is used to reduce the data, so that the reference tables has summarized observations $\eta(\mathbf{y}_{\text{obs}})$ instead of full observations. Recent methods use machine learning techniques to learn good summary statistics directly from the data [Fearnhead and Prangle, 2012, Forbes et al., 2022], but this remains a difficult task in a general setting. However, in some specific fields such as population genetics, informative summary statistics based on expert knowledge are available [Pudlo et al., 2016, Raynal et al., 2018]. In all the following, we assume that a set of such summary statistics is already known for the problem at hand. With a slight abuse of notation, we will write $\mathbf{y}_{\text{obs}}$ for $\eta(\mathbf{y}_{\text{obs}})$, so that all observations are in the summary statistics space.

## 2.2 Outlier Score

For any observation $\mathbf{y} \in \mathcal{Y}$, we define an outlier score $T(\mathbf{y}; \mathbf{D}_m)$ that should be designed to be small for points $\mathbf{y}$ that are “inside” the distribution represented by $\mathbf{D}_m$, and large for points that are “outside” of this distribution.

Lemaire et al. [2016] use the mean kNN score, that is defined as the mean distance to the $k$ nearest neighbors of a point. More specifically, for $\mathbf{y} \in \mathcal{Y}$, we denote by $N_k(\mathbf{y}; \mathbf{D}_m)$ the set of the $k$ nearest neighbors of $\mathbf{y}$ in the reference table $\mathbf{D}_m$ (excluding $\mathbf{y}$), and set:

$$ \text{knn}(\mathbf{y}, \mathbf{D}_m) = \frac{1}{k} \sum_{\mathbf{y}' \in N_k(\mathbf{y}; \mathbf{D}_m)} d(\mathbf{y}, \mathbf{y}'), \quad (1) $$

with $d$ the Euclidean distance between $\mathbf{y}$ and $\mathbf{y}'$. Note that in Lemaire et al. [2016], this score is not explicitly linked to the kNN mean distance. They also set $k$ so that a given fraction $\varepsilon$ of the reference table is used in the computation of the score (typically, $\varepsilon = 1\%$). In the following, we use our simulation study to give some guidelines to select the value of $k$ (see Section 3.1.6).

In the present work, inspired from the literature on novelty detection, we use the LOF score [Breunig et al., 2000]. It relies on the k-local reachable density in $\mathbf{y}$ with respect to the reference table $\mathbf{D}_m$:

$$ \text{lrd}_k(\mathbf{y}; \mathbf{D}_m) = \left( \frac{1}{k} \sum_{\mathbf{y}' \in N_k(\mathbf{y}; \mathbf{D}_m)} \text{reach-dist}_k(\mathbf{y}, \mathbf{y}') \right)^{-1}, \quad (2) $$

where $\text{reach-dist}_k(\mathbf{y}, \mathbf{y}') = \max\{d(\mathbf{y}, \mathbf{y}'); \text{k-dist}(\mathbf{y}'; \mathbf{D}_m)\}$, with $\text{k-dist}(\mathbf{y}'; \mathbf{D}_m)$ the $k^{th}$ nearest distance between $\mathbf{y}'$ and points in $\mathbf{D}_m$. Note that if we were to replace the reach-dist by the simple euclidean distance in the definition above, we would have $\text{lrd}_k(\mathbf{y}; \mathbf{D}_m) = \text{knn}(\mathbf{y}, \mathbf{D}_m)^{-1}$. Breunig et al. [2000] use the reach-dist to stabilize the lrd. To take into account distributions with heterogeneous densities, where some regions of the space are much more dense than others (see Fig. 1 for an illustration) the LOF statistic is defined as a ratio:
$$ LOF_k(\mathbf{y}; \mathbf{D}_m) = \frac{\frac{1}{k} \sum_{\mathbf{y}' \in N_k(\mathbf{y}; \mathbf{D}_m)} \text{lrd}_k(\mathbf{y}'; \mathbf{D}_m)}{\text{lrd}_k(\mathbf{y}; \mathbf{D}_m)} \quad (3) $$

This statistic is fully local and can be easily interpreted: for instance, $LOF_k(\mathbf{y}; \mathbf{D}_m) \approx 3$ means that the density in $\mathbf{y}$ is three times smaller than the average density of its neighborhood. As in the kNN case, the value of $k$ must be chosen by the practitioner. Breunig et al. [2000] suggest to define a “max-LOF” statistic: $\max-LOF(\mathbf{y}; \mathbf{D}_m) = \max_{k \in I} LOF_k(\mathbf{y}; \mathbf{D}_m)$ where $I$ is an interval of integers between, e.g., 5 and 20.

## 2.3 Pre-Inference Prior GoF

The goal of the pre-inference goodness of fit test (hereafter pre-inference GoF) is to check whether an observed data $\mathbf{y}_{\text{obs}} \in \mathcal{Y}$ is consistent with the prior predictive distribution of a given model $m$. Let $F$ be the true distribution function that generated the data. The null hypothesis associated with the prior-GoF test can hence be written as:

$$ H_0^{\text{prior}}: F(\cdot) = p_m(\cdot) := \int_{\Theta_m} f_m(\cdot | \theta) \pi_m(d\theta). \quad (4) $$

Following Lemaire et al. [2016], we then define the p-value associated with this test by the probability of the observed score $T(\mathbf{y}_{\text{obs}}; \mathbf{D}_m)$ to be larger than the score $T(\mathbf{y}_{\text{rep}}; \mathbf{D}_m)$ of replicates $\mathbf{y}_{\text{rep}} \sim p_m$ randomly drawn from the prior predictive:

$$ P^{\text{prior}}(\mathbf{y}_{\text{obs}}) = \mathbb{P}_{\mathbf{y}_{\text{rep}} \sim p_m} [T(\mathbf{y}_{\text{rep}}; \mathbf{D}_m) > T(\mathbf{y}_{\text{obs}}; \mathbf{D}_m)]. \quad (5) $$

Note that by definition, this p-value is well calibrated: when $H_0$ is true, $\mathbf{y}_{\text{obs}} \sim p_m$, and the p-value is uniformly distributed between 0 and 1. To estimate this p-value, we use $N_{\text{calib}}$ calibration points drawn from the prior predictive $\mathbf{y}_{\text{j}} \sim p_m$, and compare the observed score with the score of the replicates:

$$ \hat{P}^{\text{prior}}(\mathbf{y}_{\text{obs}}) = \frac{1}{N_{\text{calib}}} \sum_{j=1}^{N_{\text{calib}}} \mathbb{I}\{T(\mathbf{y}_{\text{j}}; \mathbf{D}_{m}) > T(\mathbf{y}_{\text{obs}}; \mathbf{D}_{m})\}. \quad (6) $$

## 2.4 Post-Inference Holdout GoF

The prior GoF as described above is a frequentist test, based on an empirical p-value, of a Bayesian null hypothesis that the data was distributed according to the prior predictive distribution. Taking one step further in the frequentist framework, our goal with the post-inference holdout GoF (hereafter post-inference GoF) is to test whether the true distribution $F$ is equal to the maximized likelihood of the model:

$$ H_0^{\text{holdout}}: F(\cdot) = f_m(\cdot | \theta_m^*) \quad \text{with} \quad \theta_m^* = \underset{\theta \in \Theta_m}{\operatorname{argmin}} \operatorname{KL}[F(\cdot) || f_m(\cdot | \theta)], \quad (7) $$

with KL the Kullback-Leibler divergence. From this null hypothesis and as in the previous section, we can define the following theoretical p-value:

$$ P^{\text{freq}}(\mathbf{y}_{\text{obs}}) = \mathbb{P}_{\mathbf{y}_{\text{rep}} \sim f_m(\cdot | \theta_m^*)} [T(\mathbf{y}_{\text{rep}}; \mathbf{D}_m^*) > T(\mathbf{y}_{\text{obs}}; \mathbf{D}_m^*)], $$

where $\mathbf{D}_m^*$ would be a reference table drawn from $f_m(\cdot | \theta_m^*)$ instead of $p_m(\cdot)$. As in the pre-inference-GoF test, this p-value is, by design, uniformly distributed under the null hypothesis. Note that we could use any reference table in the score $T$ used above, including the original $\mathbf{D}_m$ drawn from the prior predictive, and this
p-value would still be well calibrated. However, using the (theoretical) $D_m^*$ reference table should increase the discrimination power of the test, as, by design, we expect $T(y_{obs}; D_m^*)$ to be large when $y_{obs}$ is not drawn according to $H_0$, and small otherwise.

Unfortunately, as opposed to the prior-GoF test, where generating data from the null hypothesis $p_m(\cdot)$ was straightforward, generating data according to $f_m(\cdot | \theta_m^*)$ is not feasible, as this quantity is unknown. The strategy applied by the holdout predictive check is to approximate this quantity using the posterior distribution learned from the data. Indeed, under regularity conditions, the posterior distribution asymptotically concentrates around the maximum likelihood estimate as the sample size increases, following the Bernstein-von Mises theorem [van der Vaart, 1998]. In order to avoid any double use of the data, Moran et al. [2023] assume that part of the observations was hold-out, so that we have access to $y_{new} \sim F(\cdot)$, an independent observation drawn from the true distribution, and that we draw replicates from the predictive posterior distribution $p_m(\cdot | y_{obs}) = \int_{\Theta_m} f_m(\cdot | \theta) p(d\theta | y_{obs})$ learned from $y_{obs}$:

$$P^{\text{holdout}}(\mathbf{y}_{\text{obs}}, \mathbf{y}_{\text{new}}) = \mathbb{P}_{\mathbf{y}_{\text{rep}} \sim p_m(\cdot | \mathbf{y}_{\text{obs}})} \left[ T(\mathbf{y}_{\text{rep}}; \mathbf{D}_m^{\|\mathbf{y}_{\text{obs}}}) > T(\mathbf{y}_{\text{new}}; \mathbf{D}_m^{\|\mathbf{y}_{\text{obs}}}) \mid \mathbf{y}_{\text{obs}}, \mathbf{y}_{\text{new}} \right], \quad (8)$$

where $D_m^{\|\mathbf{y}_{\text{obs}}}$ would be a reference table drawn from $p_m(\cdot | \mathbf{y}_{\text{obs}})$. Here, the draws from the replicates no longer follow the same distribution as $\mathbf{y}_{\text{obs}}$ under $H_0$, so this p-value is no longer exactly calibrated. However, under some common assumption, Moran et al. [2023] showed that the p-value is asymptotically calibrated when the number of observations grows.

In the context of generative models, we typically do not have access to $p_m(\cdot | \mathbf{y}_{\text{obs}})$, or even to a good approximation of it. However, ABC techniques can allow us to get observations that are approximately drawn from $p_m(\cdot | \mathbf{y}_{\text{obs}})$. Assume that we have access to $N_{\text{post}} = N'_{\text{calib}} + N'_{\text{ref}}$ such draws, that we split between $N'_{\text{calib}}$ calibration points $\mathbf{y}_1, \cdots, \mathbf{y}_{N'_{\text{calib}}}$, and $N'_{\text{ref}}$ reference points stored in a new approximated posterior reference table $\hat{D}_m^{\|\mathbf{y}_{\text{obs}}}$.

$$\hat{P}^{\text{holdout}}(\mathbf{y}_{\text{obs}}, \mathbf{y}_{\text{new}}) = \frac{1}{N'_{\text{calib}}} \sum_{j=1}^{N'_{\text{calib}}} \mathbb{I} \left\{ T(\mathbf{y}_j; \hat{D}_m^{\|\mathbf{y}_{\text{obs}}}) > T(\mathbf{y}_{\text{new}}; \hat{D}_m^{\|\mathbf{y}_{\text{obs}}}) \right\}. \quad (9)$$

This quantity can be computed from the observed dataset, as long as we assume that we can split the data into two independent observations $(\mathbf{y}_{\text{obs}}, \mathbf{y}_{\text{new}})$, and that we can draw points from the approximated posterior. This last step is crucial, as the quality of the estimation of the posterior will control the quality of the approximation of the maximum likelihood that we are aiming at, and hence will impact the calibration of our p-value. In an ABC framework, a number of techniques have been proposed for posterior estimation, including methods based on a simple rejection, a local linear regression [Beaumont et al., 2002], a regression with a Ridge regulation algorithm, or a neural network regression [Blum and François, 2010]; we refer to Beaumont [2010], Marin et al. [2012] for reviews of available techniques.

# 3 Simulation Study

## 3.1 Material and Methods

We considered three datasets, described here from the most simple to the most complex, with the last two taken from established models in population genetics. In each case, we defined a model $m_0$ that was used to simulate the reference table, and a model $m_1$ that we used to simulate “pseudo-observed datasets” (PODs),
i.e., datasets that mimic observations coming from an empirical dataset. Our goal is to apply our goodness of fit procedures (prior-GoF and post-GoF) on such simulated datasets to assess whether model $m_0$ under scrutiny is good enough to describe each POD from $m_1$.

### 3.1.1 Toy “Laplace-Gaussian” Model

We first considered a “toy” dataset, with $m_0$ (Laplace) and $m_1$ (Gaussian) defined as follow:

$$ m_0: \left\{ \begin{array}{l} \theta = (\mu, \sigma) \\ \mu \sim \mathcal{U}(-5,5); \\ \sigma \sim \mathcal{U}(1,4); \\ \mathbf{z}_i \stackrel{iid}{\sim} \mathcal{L}(\mu, \sigma/\sqrt{2}) \text{ for } 1 \le i \le d; \\ \mathbf{y}_{\text{obs}} = \text{L-moments}(\mathbf{z}; m); \end{array} \right. \qquad m_1: \left\{ \begin{array}{l} \theta = (\mu, \sigma) \\ \mu \sim \mathcal{U}(-5,5); \\ \sigma \sim \mathcal{U}(1,4); \\ \mathbf{z}_i \stackrel{iid}{\sim} \mathcal{N}(\mu, \sigma^2) \text{ for } 1 \le i \le d; \\ \mathbf{y}_{\text{obs}} = \text{L-moments}(\mathbf{z}; m). \end{array} \right. $$

for each raw data $\mathbf{z}_i$, $1 \le i \le d$, both models have conditional mean $\mu$ and standard deviations $\sigma$ as variable parameters draw into uniform prior distributions with bounds (-5,5) and (1,4), respectively. As the Laplace model has heavier tails than the Gaussian model, rejecting the Laplace model with a Gaussian observation can be difficult, which is why we selected this configuration. We simulated from these models with the R statistical programming language, using the function rlaplace from the package VGAM [Yee, 2010] for Laplace simulation. To mimic the standard case where summary statistics are used, the simulated vectors $\mathbf{z}$ of dimension $d = 350$ were not observed directly, but instead summarized using the first $m = 20$ sample L-moments ratios [Hosking, 1990] using function salmu from package lmom [Hosking, 2024].

### 3.1.2 Simple “Dep-Indep” population genetics model with four populations

We considered a case study with two evolutionary scenarios of four populations using 1000 single nucleotide polymorphism (SNP) genetic markers generated following the “Hudson” criterion [Cornuet et al., 2014] and genotyped for 10 individuals per population. These two scenarios are depicted in Fig. 2. In scenario 1 (“dep”), populations diverged serially from each other, population 1 being the ancestral population (pop 1 → pop 2 → pop 3 → pop 4). In scenario 2 (“indep”), populations 2, 3 and 4 diverged independently from the ancestral population 1 (pop 1 → pop 2, pop 1 → pop 3, and pop 1 → pop 4).
<figure><img src="image_4.png" /><figcaption>Figure 2: Illustrations of the two Dep-Indep evolutionary scenarios considered in the population genetics example of Section 3.1.2. In Scenario 1 (“dep”), populations diverged serially from each other, population 1 being the ancestral population. In Scenario 2 (“indep”), the populations 3, 4 and 4 diverged independently from the ancestral population 1.</figcaption></figure>

These scenarios are often opposed when analyzing invasion or colonization histories, with scenario 1 corresponding to the case of a cascade of secondary invasion (colonization) events after a primary invasion (colonization) event and scenario 2 corresponding to the case of multiple independent invasion (colonization) events from a single source population. In both models, the parameters are $\theta = (N1, N2, N3, N4, t2, t3, t4)$, with $Ni, 1 \le i \le 4$, the effective population size of population $i$ (in number of diploid individuals), and $ti, 2 \le i \le 4$, the split time (in number of generations since present) between populations $i$ and population $i-1$ for scenario 1, or between population $i$ and population 1 for scenario 2. We used the two following Bayesian models:

$$m_0 (\text{Sc. 2}): \left\{ \begin{array}{l} Ni \sim \mathcal{U}(10^3, 10^4), 1 \le i \le 4; \\ t4 \sim \mathcal{U}(1, 60); \\ t3 \sim \mathcal{U}(61, 120); \\ t2 \sim \mathcal{U}(121, 180); \\ \mathbf{z} \sim \text{Scenario 2}(\theta); \\ \mathbf{y}_{\text{obs}} = \text{SumStat}(\mathbf{z}; m); \end{array} \right. \qquad m_1 (\text{Sc. 1}): \left\{ \begin{array}{l} Ni \sim \mathcal{U}(10^3, 10^4), 1 \le i \le 4; \\ t4 \sim \mathcal{U}(1, 60); \\ t3 \sim \mathcal{U}(61, 120); \\ t2 \sim \mathcal{U}(121, 180); \\ \mathbf{z} \sim \text{Scenario 1}(\theta); \\ \mathbf{y}_{\text{obs}} = \text{SumStat}(\mathbf{z}; m). \end{array} \right.$$

We simulated SNP datasets from each scenario with the diyabc simulation program (Collin et al., 2021¹; release v1.1.51) associated to the package DIYABC-RF v1.0 (Collin et al., 2021²). The observed and simulated SNP datasets were summarized using a total of $m = 130$ summary statistics describing genetic variation within populations (e.g., proportion of monomorphic loci, heterozygosity, population-specific FST) and between pair, triplet or quadruplet of populations (e.g., Nei's distance, Fst-related statistics, Patterson's allele-sharing f-statistics, coefficients of admixture) to describe genetic variation among various population combinations (see Collin et al., 2021 for details).

<footnote n="1">¹Available at: https://github.com/diyabc/diyabc.</footnote>

<footnote n="2">²Available at: https://github.com/diyabc/diyabc.github.io.</footnote>
### 3.1.3 More complex models inspired by evolutionary scenarios for human populations

For this third simulation-based case study, we considered two complex evolutionary scenarios derived from the study of four populations of modern humans initially comprising six evolutionary scenarios and used previously in other methodological studies (e.g., Pudlo et al., 2016). The two scenarios considered here correspond to scenarios 2 and 3 of the six scenarios detailed in the Supplementary Material Fig. S7. We chose to focus on these two scenarios 2 and 3 because they correspond to the scenarios that are both the most difficult to discriminate and the most compatible with the human SNP data available (see Section 4). For computational efficiency, our simulation-based study is limited to 5000 SNP genetic markers (simulated with the Hudson criterion; see Collin et al. 2021) for 10 individuals sampled per population. The parameters for both scenarios are $\theta = (N1, N2, N3, N4, N34, Na, t1, t2, t3, t4, Nbn3, Nbn4, Nbn34, d3, d4, d34, ra)$. In both cases, the ancestral population has effective size (in number of diploid individuals) $Na$, and becomes population 2 (ancestor to YRI) at time $t4$ (in number of generations since present), with effective size $N2$. Population 34 (ancestor to GBR and CHB) splits from population 2 at time $t3$, and goes through a bottleneck with size $Nbn34$ for a duration $d34$, before reaching a population size of $N34$. Populations 3 (ancestor to CHB) and 4 (ancestor to GBR) then split at time $t2$, and both undergo a bottleneck with respective sizes $Nbn3$ and $Nbn4$ for duration $d3$ and $d4$, before reaching their respective population sizes $N3$ and $N4$. Then at time $t1$, population 1 (ancestor to ASW and with effective size $N1$) is the result of an admixture event with population 4 (in scenario 2) or population 3 (in scenario 3) as parental populations, with an admixture rate $ra$ (proportion of genes with a non-African origin). In these scenarios, we imposed the condition $t4 > t3 > t2$. We used the same simulation program diyabc [Collin et al., 2021] and summary statistics as in the previous “Dep-Indep” case study.

We used the same priors for both models $m_0$ and $m_1$:

$$ \theta \sim \mathcal{P} \iff \left\{ \begin{array}{l} Ni \sim \mathcal{U}(10^3, 10^4) \text{ for } i \in \{1, 2, 3, 4, 34\}; \\ Na \sim \mathcal{U}(10^2, 10^3); \\ t1 \sim \mathcal{U}(1, 30); \\ ti \sim \mathcal{U}(10^2, 10^3) \text{ for } i \in \{2, 3, 4\}; \\ Nbni \sim \mathcal{U}(5, 500) \text{ for } i \in \{3, 4, 34\}; \\ di \sim \mathcal{U}(5, 500) \text{ for } i \in \{3, 4, 34\}; \\ ra \sim \mathcal{U}(0.05, 0.95), \end{array} \right. $$

and defined the two Bayesian models:

$$ m_0 (\text{Sc. 2}): \begin{cases} \theta \sim \mathcal{P}; \\ \mathbf{z} \sim \text{Scenario 2}(\theta); \\ \mathbf{y}_{\text{obs}} = \text{SumStat}(\mathbf{z}; m); \end{cases} \quad m_1 (\text{Sc. 3}): \begin{cases} \theta \sim \mathcal{P}; \\ \mathbf{z} \sim \text{Scenario 3}(\theta); \\ \mathbf{y}_{\text{obs}} = \text{SumStat}(\mathbf{z}; m). \end{cases} $$

### 3.1.4 Power Computations

For a given level of test $\alpha$, we define the power $\Pi$ of a GoF test for model $m_0$ by the probability of its p-value to be smaller than $\alpha$ when applied to a model that is different from $H_0$ under the alternative model $m_1$. Power computations for pre- and post-inference GoF are detailed in the Supplementary Material Section B, where the number of pseudo-observed datasets (POD) used a test datasets is also discussed.
### 3.1.5 Test Calibration

To test for the calibration of both tests, we reproduced the analyses described in the Supplementary Material Section B, but this time using PODs from the same prior predictive $p_{m_0}$ used for the reference table, and checked whether the estimated p-values from the $N_{test}$ PODs approximately followed a uniform distribution.

### 3.1.6 Score Computation

We computed the max-LOF score over $k$ varying from 5 to 20, following standard recommendations [Breunig et al., 2000]. For the kNN, we used $k = 1$ (nearest neighbor only), which we found to empirically give better results. Results for other values of $k$ are presented in the Supplementary Material Sections C.2 and C.4.

### 3.1.7 Practical Implementation

All the scores and GoF tests are implemented in the R package `abcgof`, freely available at https://github.com/pbastide/abcgof. We implemented critical functions such as distance computations using C++ and package Rcpp [Eddelbuettel and François, 2011] for improved speed, and parallelized calculations using package `doParallel` [Weston, 2022].

To compute p-values, both the kNN and LOF scores require the computation of the $k$ distances to the nearest neighbors of the test and calibration points to the reference points, for a varying value of $k$. A naive implementation can compute all distances between calibration and reference points, and then sort the distances. Here, we used the kNN implementation from the `dbscan` R package, that relies on a k-d tree search to speed up the computation of the $k$ nearest neighbors Hahsler et al. [2019].

All the code to reproduce the figures and analyses of the paper is available on GitHub: https://github.com/pbastide/gof_sbi_paper.

## 3.2 Results

The main results of our simulation study for the three settings presented above are described here. Further analyses of these simulations are detailed in Supplementary Material Section C.

### 3.2.1 LOF is a good metric for GoF tests

Compared to the previous study of Lemaire et al. [2016], we chose here to use the LOF score instead of the kNN in all p-values computations, with the idea that, the LOF score being specifically adapted for outlier detection, it will be better at discriminating data points that are not from the null distribution. In all of our tests, we saw that the LOF score indeed performed better than the kNN, with higher power (see Fig. 3 and 4), and better calibration (see Fig. 5).
<figure><img src="image_5.png" /><figcaption>Figure 3: **Estimated power for the pre-inference (prior) GoF test**, at a level of 5%, for the three simulation settings (Laplace-Gaussian, Dep-Indep and Human-Like). The LOF score (black) is the “max-LOF” score for *k* between 5 and 20. The kNN score (light orange) is computed for *k* = 1. In all settings, *N*<sub>ref</sub> = *N*<sub>calib</sub>, and the total number of simulated particles is *N*<sub>ref</sub> + *N*<sub>calib</sub> (x axis).</figcaption></figure>

**3.2.2 Pre-inference GoF is Powerful and Well Calibrated**

In particular, the pre-inference (prior) GoF test associated with the LOF metric had very good power, above
or close to 0.9 even for a small number of simulated points in the two “Laplace-Gaussian” and “Dep-Indep”
simple simulation settings (Fig. 3). In the more complex “Human-like” setting, which includes two closely
related scenarios, More simulated particles were needed to achieve a reasonable power of about 0.65 with
the LOF score, which is much higher than the about 0.3 obtained with the kNN score. As expected, the pre-
inference GoF also proved to be well calibrated, with p-values under the null hypothesis being approximately
uniformly distributed (see Supplementary Material Figure S3). We also tested a localized version of the pre-
inference GoF, as in Lemaire et al. [2016]. This test did not exhibit a significative increase in power, despite
being much more computationally intensive (see Supplementary Material Figure S1).

**3.2.3 Post-inference GoF Power is Robust to the Posterior Approximation**

The post-inference GoF test had good power, i.e. above 0.8, even in the most complex “Human-Like”
simulation setting (Fig. 4). The posterior inference method had little impact on the computed power, the
simple rejection method performing as well or better than the local linear or ridge regression (see Supple-
mentary Figure S5). Doubling the number of particles, from 50 to 100 thousands, or changing the level of
localization, from 0.5% to 4%, had little effect.
<figure><img src="image_6.png" /><figcaption>Figure 4: **Estimated power for the post-inference holdout GoF test**, at a level of 5%, for the three simulation settings (Laplace-Gaussian, Dep-Indep and Human-like) and the rejection ABC posterior estimation method. The posterior was localized with $\epsilon$ such that the number of posterior particles is equal to $N_{\text{post}}$ (x axis), starting from $N_{\text{ref}} = 50\,000$ (full line) and $100\,000$ (dashed line) total number of simulated particles. This corresponds to 0.5% to 4% of the simulations, depending on $N_{\text{ref}}$. The LOF score (black) is the “max-LOF” score for $k$ between 5 and 20. The kNN score (light orange) is computed for $k = 1$. In all settings, and $N'_{\text{calib}} = N'_{\text{ref}} = N_{\text{post}}/2$.</figcaption></figure>

### 3.2.4 Post-inference GoF Calibration is Sensitive to the Posterior Approximation

The simple Rejection method with the LOF score seemed to be relatively robust, with the empirical quantiles of the p-values close to the theoretical uniform quantiles (Fig. 5), yielding for approximately calibrated tests. As expected, the quality of the posterior approximation had however a substantial effect on the calibration of the post-inference GoF test (see Supplementary Figure S6). In the simple “Laplace - Gaussian” setting, that has only two parameters, the local linear regression seems enough to correct the posterior estimation, and to provide calibrated p-values. In the “Dep-Indep” setting, the ridge regression seems to be needed to correctly accommodate for the greater number of parameters. Finally, in the more complex “Human Like” setting, none of the simple posterior correction measures seemed to provide for calibrated p-values.
<figure><img src="image_7.png" /><figcaption>Figure 5: **Distribution of p-values of the post-inference GoF test** when the PODs were simulated according to the null model, for the three simulation settings (Laplace-Gaussian, Dep-Indep and Human-like) and the rejection ABC posterior estimation method. The total number $N_{\text{ref}}$ of simulated particles was taken to be equal to 50 000 (full lines) or 100 000 (dashed lines), and 1 000 points were taken in the posterior sample, with $N'_{\text{calib}} = N'_{\text{ref}} = N_{\text{post}}/2$. The LOF score (black) is the “max-LOF” score for $k$ between 5 and 20. The kNN score (light orange) is computed for $k = 1$.</figcaption></figure>

# 4 Illustration using a real SNP dataset of modern Human populations

## 4.1 Material and Methods

### 4.1.1 Data

To illustrate our methodological developments, we analyzed a real SNP dataset obtained from individuals originating from four Human populations (30 unrelated individuals sequenced per population) using the freely accessible public dataset from The 1000 Genomes Project Consortium [2012]<fnref n="3" />. The four Human populations included the Yoruba population (Nigeria) as representative of Africa (encoded YRI in the 1000 genome database), the Han Chinese population (China) as representative of the East Asia (encoded CHB), the British population (England and Scotland) as representative of Europe (encoded GBR), and the population composed of Americans of African Ancestry in SW-USA (encoded ASW). The SNP loci were selected from the 22 autosomal chromosomes using the following criteria: (i) all 30 × 4 analyzed individuals have a genotype characterized by a quality score (GQ) > 10 (on a PHRED scale), (ii) polymorphism is present in at least one of the 30 × 4 individuals in order to fit the SNP simulation algorithm used in DIYABC Random Forest v1.0 [Collin et al., 2021], (iii) the minimum distance between two consecutive SNPs is 1 kb in order to minimize linkage disequilibrium between SNP, and (iv) SNP loci showing significant deviation from Hardy-Weinberg equilibrium at a 1% threshold [Wigginton et al., 2005] in at least one of the four populations has been removed (35 SNP loci concerned). After applying the above criteria, we obtained a dataset including 51 250 SNP loci scattered over the 22 autosomes (with a median distance between two consecutive SNPs equal to 7 kb), among which two subsets each of 12 000 SNP loci (thereafter referred to as dataset 1 and dataset 2) were randomly chosen for applying our GoF methods.

<footnote n="3">³vcf format files including variant calls are available at http://www.1000genomes.org/data.</footnote>
4.1.2 Evolutionary Scenarios

We considered six scenarios (i.e. models) of evolution of the four Human populations which differ from
each other by one ancient and one recent historical event: (i) A single out-of-Africa colonization event
giving an ancestral out-of-Africa population which secondarily split into one European and one East Asian
populational lineage, versus two independent out-of-Africa colonization events, one giving the European
lineage and the other one giving the East Asian lineage. The possibility of a second ancient (i.e. older
than 100000 years) out of Africa colonization event through the Arabian peninsula toward Southern Asia
has been suggested by archaeological studies (e.g., Rose et al., 2011). (ii) The possibility (or not) of a
recent genetic admixture of the Americans of African Ancestry in SW-USA between their African ancestors
and individuals of European or East Asia origins (e.g., Bryc et al., 2015). We detail the six considered
evolutionary scenarios in the Supplementary Material Figure S7. Human population history has been and is
still studied in details in a number of studies using genetic data (e.g., Lohmueller and Nielsen, 2021). Our
aim here is to illustrate the potential of our pre-inference and post-inference GoF methodologies using real
data on a well known case-study that exhibits complex evolutionary histories.

4.1.3 Simulation of SNP Datasets and Summary Statistics

We used the same parameters and priors for simulations as in Section 3.1.3, with conditions on time events
such that t4 > t3 > t2 for scenarios 1, 2 and 3, and t4 > t3 and t4 > t2 for scenarios 4, 5 and 6. We used
the same diyabc and DIYABC-RF simulation tools Collin et al., 2021 and the same 130 summary statistics
as in Section 3.1.2.

4.1.4 Pre-inference GoF for Scenario Pruning

We used the pre-inference (prior) GoF method to quickly, i.e. with a low computational cost, eliminate from
the six proposed scenarios those characterized by a low propensity to generate multidimensional particles
with summary statistics close to the particle corresponding to the observed dataset 1 (or the observed dataset
2). To this aim, we simulated a reference table including 500 to 5000 particles (i.e. SNP datasets sum-
marized by the above mentioned summary statistics) per scenario. We then used the pre-inference GoF
module of our new abcgof R package to compute the probability values of each scenario using the kNN
outlier score (with *k* = 1) or the LOF outlier score (using the “max-LOF” score for *k* between 5 and 20); see
Section 2 for details. The probability values of each scenario were calculated for a number of simulations
equal to 500 (from which 250 particles were used for calibration), 1000 (500 for calibration), 2000 (1000
for calibration), 3000 (1500 for calibration), 4000 (2000 for calibration), and 5000 (2500 for calibration).
95% confidence intervals were computed based on asymptotic standard error estimation or on 500 bootstrap
replicates, as detailed in Section 4.2 below. To avoid inflated Type I error rates due to multiple testing
issues, we applied the Benjamini-Hochberg [Benjamini and Hochberg, 1995] correction to the bootstrap
upper bounds of the LOF p-values computed with 2000 particles.

4.1.5 Scenario Choice Using ABC Random Forest

Following Pudlo et al. [2016], we used Random Forests to select the best scenario between the two scenarios
that were not eliminated after the pre-inference GoF treatment (i.e. scenarios 2 and 3). We simulated a
training set (i.e. reference table) including 10000 particles for each of the scenarios 2 and 3 using the same
set of 130 statistics plus one LDA axis (i.e., the number of scenarios minus 1 as predictors [Pudlo et al.,
[2016]. We then used the “Random Forest analyses” module of DIYABC Random Forest v1.0 [Collin et al., 2021] to process model choice on the training set to predict the best scenario and estimate its posterior probability. We fixed the number of trees in the constructed Random Forests to 1000, as this number turned out to be large enough to ensure a stable estimation of the global error rate (results not shown).

### 4.1.6 Post-Inference GoF Test on the Selected Scenario

We applied the post-inference GoF test on the above-selected scenario (i.e. scenario 2) to evaluate whether the selected scenario missed some important evolutionary features to explain the observed dataset, so that the multi-parameter posterior distribution conditional on the observed dataset has a poor propensity to generate multi-dimensional summary statistics particles close to the particle corresponding to the observed dataset. To this aim, we simulated a reference table including 50000 or 100000 particles under scenario 2. The target (i.e. the particle of the observed dataset) was either the particle of the observed dataset 1, with the particle of the observed dataset 2 used as replicate, or the particle of the observed dataset 2, with the particle of the observed dataset 1 used as replicate. Two overlapping criteria were used to implement our post-inference GoF calculation. The first criterion was to obtain a sufficiently strong localization around the particle corresponding to the observed dataset by retaining only the $\varepsilon$ simulated particles, with $\varepsilon < 5\%$, that are closest. The second criterion was to have a sufficiently large number of retained particles (more than 500) to allow a reasonably robust linear or ridge regression in the local region defined above, as well as a reasonably robust calibration using half, and thus more than 250, of the retained particles for calibration (the other 250 being used to calculate probabilities); see Section 2 for details. We thus computed post-inference GoF probabilities using $\varepsilon$ equal 1%, 2% and 4% for the reference table with 50000 particles or $\varepsilon$ equal 0.5%, 1% and 2% for the reference table with 100000 particles. We used the post-inference GoF module of our new abcgoF R package to compute the post-inference probability value of scenario 2 using the rejection method or the rejection followed by a regression (linear or ridge) method, with the kNN score (for $k=1$) or the LOF score (fusing the “max-LOF” score for $k$ between 5 and 20). As explained in Section 2.4, the GoF post-inference method necessitates to have simulation replicate of each particle included in the subset of the $\varepsilon$ retained particles. We used the {o option of the command-line simulation software diyabc v1.1.51 to simulate replicates of the $\varepsilon$ retained particles from their associated parameter vector values. 95% confidence intervals were computed based on asymptotic standard error estimation or on 500 bootstrap replicates, as detailed in Section 4.2 below.

## 4.2 Quantification of uncertainty on p-values

In the pre-inference GoF case, the estimator of Equation (6) is unbiased ($E[\hat{p}_{\text{prior}}(\mathbf{y}_{\text{obs}})] = P^{\text{prior}}(\mathbf{y}_{\text{obs}})$), and has a variance $\text{Var}[\hat{p}_{\text{prior}}(\mathbf{y}_{\text{obs}})] = P^{\text{prior}}(\mathbf{y}_{\text{obs}})(1 - P^{\text{prior}}(\mathbf{y}_{\text{obs}}))/N_{\text{calib}}$ that decreases linearly with $N_{\text{calib}}$.

Under the assumptions for good posterior estimation, the same is true for the post-inference GoF case. This formula gives us a computationally cheap way of computing asymptotic confidence intervals, using a simple plug-in estimator for the standard error, from a single draw of calibration points among the reference table.

A different approach, more costly in terms of computing time, to uncertainty assessment is to use a bootstrap procedure. In the prior GoF test, given a certain number of simulations, the first step is to choose which points we use as “calibration” points, and which we use as “reference” points, which we do by simply sampling uniformly without replacements points to be in the calibration set. We can repeat this sampling several times, in order to get different sets of calibration points from the simulated points, on which we apply the procedure to compute the p-values. In the post-inference GoF test, the same procedure can be
applied after the localization step, among the selected or inferred points from the posterior. In both case,
this gives us a distribution of p-values, that we can use to get a point estimate (here using the median) and
a confidence interval, here using the highest density interval (HDI, Hyndman, 1996, Kruschke, 2015) at a
95% level.

4.3 Results

4.3.1 Pre-inference GoF for Scenario Pruning

Prior-GoF probability values were high (i.e. larger than 5%) for only two of the six scenarios considered, i.e. for scenarios 2 and 3 (see Figure 6 for results for the SNP dataset 1). This result was well supported for both the LOF and kNN outlier scores, with upper bounds of the 95% confidence intervals (using the bootstrap method) larger than 5%, when the number of simulations was larger than 500 particles per scenario. For a number of simulations of 500 particles per scenario, the probabilities of scenario 5 were close to 5% (i.e. 0.05 and 0.06 for the LOF and kNN scores, respectively), and the upper bounds of the 95% intervals became close to 10% for this scenario. We note that the 95% confidence intervals were significantly narrower when calculated using the asymptotic rather than the bootstrap method (see Supplementary Material Figure S8). Consistent with our simulation-based studies, which show a higher power for the LOF score than for the kNN score, the pre-inference GoF probability values tend to be larger for the kNN than for the LOF score. When applying the Benjamini-Hochberg [Benjamini and Hochberg, 1995] correction to the bootstrap upper bounds of the LOF p-values with N<sub>ref</sub> + N<sub>calib</sub> = 2000, we still rejected all scenarios but scenarios 2 and 3 at the 5% level. We obtained similar results for the real SNP dataset 2 (results not shown).

As a reminder, the purpose of the scenario pruning step is to limit the simulation of unnecessary data
from clearly “out-of-game” models due to the high computational cost involved. We show here on real data
that even for complex scenarios with a large number of parameters, this goal can be achieved with a limited
number of simulations on the order of 1000 particles per scenario. Overall, our pre-inference GoF results
show that the scenarios 1, 4, 5 and 6 can be set aside, leaving only the scenarios 2 and 3 to be kept for
further simulations and inference. The next inference step is to apply statistical treatments of model choice
to select the best between scenarios 2 and 3. To this end, we simulated a total of 10000 particles for each of
the scenarios 2 and 3 and used random forests for the selection of the best scenario. Scenario 2 was clearly
the best, with posterior probabilities of 1.000 and 0.999 for the real datasets 1 and 2, respectively.
<figure><img src="image_8.png" /><figcaption>Figure 6: **Pre-inference GoF results on the real SNP dataset 1 of modern Human populations for different simulation efforts using the LOF or the kNN outlier score.** Estimated p-values for the prior GoF test, for the six scenarios of Figure S7. The LOF score (black) is the “max-LOF” score for *k* between 5 and 20. The kNN score (light orange) is computed for *k* = 1. In all settings, *N*<sub>ref</sub> = *N*<sub>calib</sub>, and the total number of simulated particles is *N*<sub>ref</sub> + *N*<sub>calib</sub> (x axis). Vertical bars represent the 95% Highest Density Intervals (HDI) computed from the bootstrap procedure described in Section 4.2 with 500 replicates. See Supplementary Material Figure S8 for a comparison with asymptotic intervals. The 5% rejection threshold is shown as a dashed line. Only scenarios 2 and 3 are consistently not rejected, even at low computational cost (i.e. 1000 simulations per scenario). The black and orange curves are superimposed for scenarios 1 and 4, so that only the orange curve is visible for these scenarios.</figcaption></figure>

4.3.2 Post-inference GoF on the selected scenario

For the rejection method, post-GoF probability values of scenario 2 were small for both the LOF and kNN
scores, with upper bounds of the 95% CIs below 5%, indicating that scenario 2 poorly fits the observed
datasets (see Figure 7 for results for the SNP dataset 1). This result was slightly less well supported when
the number of simulations was 50000 (compared to 100000), as indicated by upper bounds of the 95% HDI
that were above 5% (but below 10%) for the LOF score (but not for the kNN score). Our simulation-based
studies on the human-like case have shown a similar power for the LOF score than for the kNN score. In
agreement with this, the post-inference GoF probability values obtained on real data were similar for both
outlier scores, but 95% CIs were narrower for the kNN than the LOF scores. We note again that, as for the
pre-inference GoF, the 95% confidence intervals were significantly narrower when calculated from a single
draw of calibration particles than when calculated from 500 bootstrapped iterations of calibration particles
(see Supplementary Material Figure S9 and S10).

For the local linear and ridge regression methods, the post-inference GoF probabilities of scenario 2 are higher than for the rejection method (see Supplementary Material Figure S9). Probability values of scenario 2 using the local linear and ridge regression methods show strong variations depending on the number of simulations and ε parameter, and generally do not indicate that scenario 2 fits the observed data
poorly (contrary to what was found using the rejection method). Overall, these results are consistent with our simulation-based studies on the human-like case, which show a higher power for the rejection method than for the local linear (but not for the ridge regression) methods, as well as a clearly better calibration for the rejection method.

For both with or without-regression post-inference GoF methods, we observe a substantial decrease in the probability values and their variation when we move from 50 000 to 100 000 simulations, as well as when $ε$ values are getting smaller (Figure 7 and Supplementary Material Figure S9). However, even in such more “optimal” conditions, the evidence in favour of the rejection of the scenario 2 remain weak and uncertain for the two post-inference GoF methods using a regression step. For instance for 100 000 simulations and $ε = 0.005$, the post-inference GoF probability of the method using local linear regression is 0.072 (95%CI [0.036, 0.104]) for kNN and 0.056 (95%CI [0.004; 0.216]) for LOF. We obtained similar results for the real dataset 2 (results not shown).

<figure><img src="image_9.png" /><figcaption>Figure 7: **Post-inference GoF results on the real SNP dataset 1 of modern Human populations for different simulation efforts using the LOF or the kNN outlier score.** Estimated p-values for the post-inference GoF test for the selected scenario 2 of Figure S7, using the rejection ABC posterior estimation method. The posterior was localized with $ε$ such that the number of posterior particles is equal to $N_{\text{post}}$ (x axis), starting from 50 000 (left panel) and 100 000 (right panel) total number of simulated particles. In all settings, $N'_{\text{calib}} = N'_{\text{ref}} = N_{\text{post}}/2$. The LOF score (black) is the “max-LOF” score for *k* between 5 and 20. The kNN score (light orange) is computed for *k* = 1. Vertical bars represent the 95% Highest Density Intervals (HDI) computed from the bootstrap procedure described in Section 4.2 with 500 replicates. The 5% rejection threshold is shown as a dashed line. Tested scenarios 2 gets consistently rejected by the post-inference test. See Supplementary Material Figures S9 and S10 for a comparison with other posterior estimation methods and for asymptotic intervals.</figcaption></figure>

# 5 Discussion

In population genetics as well as many other application fields, inferences about models with intractable likelihoods rely on simulation-based methods, such as Approximate Bayesian Computation and Simulation-Based Inference, for tasks such as parameter inference and model selection. Once a model has been selected, the assessment of its Goodness of Fit (GoF) presents a significant challenge, particularly in the absence of a likelihood function. We introduce two novel GoF tests based on Local Outlier Factor (LOF) that more effectively captures local density deviations than a more traditional k-Nearest Neighbors approach. The pre-inference (prior) GoF test serves to prune the less plausible models (with a low computational cost),
while the post-inference GoF test evaluates the fit of the selected model in greater detail. This dual approach provides a comprehensive assessment of model validity, combining Bayesian and general goodness-of-fit insights. The dual-test GoF approach employed by our methodology underscores its flexibility. The pre-inference GoF test offers insight into model validity from a Bayesian perspective, while the post-inference test provides a more general and traditional view of assessing goodness of fit.

## 5.1 Biological Insights for the Human Population Genetics Dataset

Our application of the pre-inference prior GoF methodology for modern human evolutionary scenario prun-
ing demonstrates that it is feasible to avoid the need for extensive, data-intensive simulations for four ap-
parently implausible scenarios among the six complex scenarios examined. This goal can be effectively
achieved with a modest simulation volume, i.e. around 1 000 particles per scenario. The four scenarios
identified as unlikely include events that are inconsistent with current knowledge of the evolution of *Homo
sapiens*. Specifically, they posit two independent colonization events outside Africa, one leading to Eu-
ropean populations and the other to northeast Asian populations (e.g. Chinese populations distinct from
those of South Asia along coastal regions; López et al. [2015]), or they overlook genetic admixture among
American individuals of African descent [Bryc et al., 2015].

The two scenarios that continue to be supported by high probabilities of pre-inference GoF are closely
related in terms of genetic history, as they both propose a single colonization event outside Africa and
an admixture source for African American ancestry, either from Europe (scenario 2) or Asia (scenario 3).
These two regions contain genetically weakly differentiated populations, making it difficult to distinguish
conclusively between scenarios 2 and 3 at this stage. A more definitive decision in favor of scenario 2
requires a formal scenario choice test, which we perform here using a random forest based classification
on 10 000 simulations generated specifically for scenarios 2 and 3. This formal approach provides a more
reliable means of distinguishing between the scenarios, given their nuanced genetic similarities.

Applying safely our post-inference GoF methodology to scenario 2 requires a higher number of simula-
tions for the selected scenario only, of the order of 100 000. The post-GoF treatments, at least those based on
the rejection method, show that (even) scenario 2 does not give a good account of the real observed Human
datasets (i.e. low post-inference GoF probability well below 5%). This result is not surprising as a number
of human population evolutionary modalities highlighted in the last two decades are not modeled in scenario
2. These include more complex demographic fluctuations than those modeled in scenario 2, the mixing of
ancestral European and Asian populations with Neanderthal and Denisovan human lineages respectively,
several successive waves of colonization of the European area, and recurrent episodes of back into Africa
migrations (e.g. reviewed in López et al., 2015).

## 5.2 Calibrated Post-Inference GoF Requires Good Posterior Estimation

The accuracy of the post-inference likelihood test is contingent on the quality of the posterior approxima-
tion. While traditional methods, such as rejection-based Approximate Bayesian Computation (ABC), offer
a straightforward approach to posterior approximation, they frequently exhibit deficiencies related to bias
and inefficiency. These limitations necessitate the employment of more sophisticated posterior approxima-
tion techniques. Recent advancements in machine learning have introduced powerful methods for posterior
approximation, such as normalizing flows, specifically Masked Autoregressive Flows (MAF, Papamakarios
et al. [2017]), and Distributional Random Forests (DRF, Cevid et al. [2022]). MAF offer a flexible way to ap-
proximate complex posterior distributions by learning dependencies between parameters in a tractable man-
ner. Sequential Neural Likelihood [Papamakarios et al., 2019], a technique that has emerged as a promising method for refining posterior approximations through iterative learning, thereby improving the convergence and accuracy of the inference, takes advantage of MAF. Furthermore, importance sampling and adaptive ABC methods present additional opportunities to enhance the posterior by refining simulations in regions of higher posterior density. Also, the use of a non-parametric approach, like the DRF, can adapt to the structure of the data while mitigating the risk of overfitting.

## 5.3 More Useful Models

Finally, as George Box famously stated, “All models are wrong, but some are useful”. This perspective is particularly relevant in the context of Bayesian generative models, where the goal is not to find a perfect representation of reality but rather to develop a model that is useful for making inferences and predictions. Given the inherent simplifications and assumptions made in modeling, it is crucial to interpret the results of GoF tests with an understanding of their limitations. Model misspecification is an unavoidable aspect of statistical inference. While GoF tests are essential for identifying discrepancies between observed data and model predictions, they can be viewed as tools for iterative model refinement rather than definitive tests of validity. When a model fails a GoF test, several actionable steps can be considered: (i) refinement of priors and summary statistics: adjusting priors to better capture known biological constraints or selecting more appropriate summary statistics can significantly improve model fit; and (ii) incorporating additional complexity: introducing features such as population structure, environmental influences, or time-varying parameters can enhance model realism. The iterative process of model criticism and refinement is key to advancing our understanding of complex systems, and GoF tests can serve as a critical component of this process.

## Scripts and Data

Data and scripts to reproduce the analyses presented in this paper. Also available on GitHub: https://github.com/pbastide/gof_sbi_paper

## Acknowledgments

We thank Paul Verdu for useful discussions about recent findings regarding Human populations evolution. We are grateful to the INRAE MIGALE bioinformatics facility (MIGALE, INRAE, 2020. Migale bioinformatics Facility, doi: 10.15454/1.5572390655343293E12) for providing computing resources.

## Funding

This work has been supported by funds from the project AgroStat (reference: ANR-23-EXMA-0002) of the Maths-VivES France 2030 program handled by the French *Agence Nationale de la Recherche*.
# Supplementary Material

## A Localized Prior GoF

### A.1 Localized Prior GoF

Lemaire et al. [2016] developed a localized version of the prior GoF, that tests the prior null hypothesis of Equation (4), but, instead of using the simple reference table $D_m$, use a localized version of it around the observation $y_{obs}$. The estimated p-value for this test is then defined as:

$$ \hat{P}^{\text{loc prior}}(\mathbf{y}_{\text{obs}}) = \frac{1}{N_{\text{calib}}} \sum_{j=1}^{N_{\text{calib}}} \mathbb{I}\{T(\mathbf{y}_j; \hat{\mathbf{D}}_{m}^{\mathbf{y}_{\text{obs}}}) > T(\mathbf{y}_{\text{obs}}; \hat{\mathbf{D}}_{m}^{\mathbf{y}_{\text{obs}}})\}, \quad (10) $$

where $\hat{\mathbf{D}}_{m}^{\mathbf{y}_{\text{obs}}}$ is obtained as in Equation (9) of the main text, by re-simulating the lines of the reference tables that are the closest to the observation (in a simple rejection approach). This test hence has the same computational cost as the post-inference HPC GoF.

Note that this test is not a posterior predictive check, as the null hypothesis corresponds to the prior predictive distribution, which is why we refer to it as a localized prior GoF test.

### A.2 Performance on simulations

We used this localized prior GoF on the same three simulation settings as in the main text, simulating 50 or 100 thousand points for each model, and using a rejection localized with $\epsilon = 1\%$. We compared the results of the localized prior GoF to the results of the standard prior GoF described in Section 2.3 of the main text.

Figure S1 shows the compared power of the two approaches. We can see that, despite being much more computationally intensive (it requires 10 to 200 times more simulations from the model), the localized GoF does not improve the power at all in the Laplace-Gaussian and Dep-Indep examples (Fig. S1, first two lines). In the more complex Human-like example, the power increases only slightly.

We hence advocate for using the light weight prior GoF as described in the main text, and, if computational power is available, to turn to the holdout GoF instead of the localized prior GoF, that does not exhibit clear power gains.
<figure><img src="image_10.png" /><figcaption>Figure S1: Estimated power for the prior (left) and localized prior (right) GoF tests, at a level of 5%, for the three simulation settings. The LOF score (black) is the “max-LOF” score for k between 5 and 20. The kNN score (light orange) is computed for k = 1. In all settings, N<sub>ref</sub> = N<sub>calib</sub>, and the total number of simulated particles is N<sub>ref</sub> + N<sub>calib</sub> (x axis).</figcaption></figure>

B Power Computations

For a given level of test α, we define the power Π of a GoF test for model m₀ by the probability of its p-value to be smaller than α when applied to a model that is different from H₀ under the alternative model m₁.

B.1 Pre-inference GoF

For the pre-inference (prior) GoF of Equations 5, we get:

$$
\begin{align*}
\Pi^{\text{prior}} &= \mathbb{P}_{\mathbf{y}_{\text{obs}} \sim p_{m_1}} [\mathcal{P}^{\text{prior}}(\mathbf{y}_{\text{obs}}) \le \alpha] \\
&= \mathbb{P}_{\mathbf{y}_{\text{obs}} \sim p_{m_1}} \left[ \mathbb{P}_{\mathbf{y}_{\text{rep}} \sim p_{m_0}} \left[ T(\mathbf{y}_{\text{obs}}; \mathbf{D}_{m_0}) > T(\mathbf{y}_{\text{rep}}; \mathbf{D}_{m_0}) \right] \le \alpha \right].
\end{align*}
$$

To estimate this quantity, we draw N<sub>test</sub> observations from the prior predictive of alternative model *m*₁: **y**<sub>*k*</sub><sup>1</sup> ~
*p*<sub>*m*₁</sub>, 1 ≤ *k* ≤ N<sub>test</sub>, and we approximate both probabilities by their frequencies as in Equation 6, assuming
that we have calibration points from the prior predictive of $m_0$: $\mathbf{y}_j^0 \sim p_{m_0}$, $1 \le i \le N_{\text{calib}}$:

$$
\begin{align*}
\hat{\Pi}^{\text{prior}} &= \frac{1}{N_{\text{test}}} \sum_{k=1}^{N_{\text{test}}} \mathbb{I}\{\hat{P}^{\text{prior}}(\mathbf{y}_k^1) \le \alpha\} \\
&= \frac{1}{N_{\text{test}}} \sum_{k=1}^{N_{\text{test}}} \mathbb{I}\left\{\frac{1}{N_{\text{calib}}} \sum_{j=1}^{N_{\text{calib}}} \mathbb{I}\{T(\mathbf{y}_j^0; \mathbf{D}_{m_0}) > T(\mathbf{y}_k^1; \mathbf{D}_{m_0})\} \le \alpha\right\} \\
&= \frac{1}{N_{\text{test}}} |\{k \mid T(\mathbf{y}_k^1; \mathbf{D}_{m_0}) \ge T(\cdot; \mathbf{D}_{m_0})_{1-\alpha}, 1 \le k \le N_{\text{test}}\}|,
\end{align*}
$$

with $T(\cdot; \mathbf{D}_{m_0})_{1-\alpha}$ the $1-\alpha$ quantile of $\{T(\mathbf{y}_j^0; \mathbf{D}_{m_0}) \mid 1 \le j \le N_{\text{calib}}\}$, and $|\cdot|$ the cardinal of the set.

For the computation of pre-inference GoF power, we took $N_{\text{ref}}$ between 250 and 2500 for the reference table, and used $N_{\text{calib}} = N_{\text{ref}}$ in Formula 6, so that the total number of simulations from $m_0$ varied from 500 to 5000.

## B.2 Post-inference GoF

For the power of the post-inference (holdout) GoF test, we use the same definitions for the post-inference GoF, using Equations 8 and 9. We first draw $\mathbf{y}_k^1 \sim p_{m_1}$, $1 \le k \le N_{\text{test}}$, observations from the prior predictive of model $m_1$. Then, for each $k$, we draw a replicate $\mathbf{y}_k^{1,\text{new}}$ using the same model $m_1$ on the parameters $\theta_k$ that were used for $\mathbf{y}_k^1$, so that $\mathbf{y}_k^{1,\text{new}}$ and $\mathbf{y}_k^1$ are drawn from the same model, with the same parameter values. The estimated power in this case is given by:

$$
\begin{align*}
\hat{\Pi}^{\text{holdout}} &= \frac{1}{N_{\text{test}}} \sum_{k=1}^{N_{\text{test}}} \mathbb{I}\{\hat{P}^{\text{holdout}}(\mathbf{y}_k^1, \mathbf{y}_k^{1,\text{new}}) \le \alpha\} \\
&= \frac{1}{N_{\text{test}}} \sum_{k=1}^{N_{\text{test}}} \mathbb{I}\left\{\frac{1}{N'_{\text{calib}}} \sum_{j=1}^{N'_{\text{calib}}} \mathbb{I}\{T(\mathbf{y}_j^0; \hat{\mathbf{D}}_{m_0}^{1,\text{new}}) > T(\mathbf{y}_k^{1,\text{new}}; \hat{\mathbf{D}}_{m_0}^{1,\text{new}})\} \le \alpha\right\}.
\end{align*}
$$

For the computation of post-GoF power, we took $N_{\text{ref}} \in \{5000, 10000\}$, and, for each pseudo-observed dataset (POD) $\mathbf{y}_k$, $1 \le k \le N_{\text{test}}$, we estimated the posterior distribution using a simple rejection-based ABC (see, e.g., Marin et al., 2012) by selecting the lines of the original reference table $\mathbf{D}_{m_0}$ that were within a distance of $\varepsilon$ of $\mathbf{y}_k$:

$$
\mathbf{D}_{m_0}^{1, \varepsilon} = \{d \in \mathbf{D}_{m_0} \mid \|d - \mathbf{y}_k\| \le \varepsilon\}.
$$

We chose $\varepsilon$ so as to keep a total of $N_{\text{post}}$ particles in the localized dataset, with $N_{\text{post}}$ varying between 500 and 2000. This corresponds to 0.5% to 4% of the simulations, depending on $N_{\text{ref}}$. For each line $i$ of $\mathbf{D}_{m_0}^{1, \varepsilon}$, we then extracted the parameters $\theta_i$ from the reference table, and re-simulated a set of observations $\mathbf{y}_i$ using the same model $m_0$ and the same parameters $\theta_i$, to get $\hat{\mathbf{D}}_{m_0}^{1, \varepsilon}$. This re-simulation aimed at approximating draws from the unknown posterior $p_m(\cdot | \mathbf{y}_k)$.

It is clear from Equation (9) that the posterior approximation method should have a strong impact on the calibration of the post-inference GoF test. In addition to the simple Rejection ABC, we tested the local linear regression with an Epanechnikov kernel [Beaumont et al., 2002], taking the median predictions as parameter point estimates, and the ridge regression, taking the median predictions from three penalisation parameter values $\lambda$ equal to 0.0001, 0.001, 0.01 [Blum and François, 2010]. In both cases, we selected the same lines
$\mathbf{D}_{m_0}^{\left|\mathbf{y}_k^1\right|, \varepsilon}$ from the reference table, but then computed $\hat{\theta}_i$ using the regression, and finally re-simulated $\hat{\mathbf{D}}_{m_0}^{\left|\mathbf{y}_k^1\right|, \varepsilon}$ using these estimated parameters.

In all cases, we used again half of the lines for the calibration, that we picked at random in the posterior reference table $\hat{\mathbf{D}}_{m_0}^{\left|\mathbf{y}_k^1\right|, \varepsilon}$, so that $N'_{calib} = N_{post}/2$ and $\hat{\mathbf{D}}_{m_0}^{\left|\mathbf{y}_k^1\right|}$ had $N'_{ref} = N_{post}/2$ lines in Formula 9.

## B.3 Number of test PODs

We took $N_{test} = 1000$ for the first two simulation settings (i.e. the Toy “Laplace-Gaussia” toy Models and “Dep - Indep” population genetics scenarios), and $N_{test} = 500$ for the third one (i.e. the “Human-like” evolutionary scenarios), which was more computationally demanding.

# C Additional Simulation Results

## C.1 Parameter Transform

When using a regression approach (local linear or ridge), the correction term does not take the constraints of the parameters into account, and can produce posterior predictions that lie outside of the set boundaries. Instead of just rejecting these prediction, another approach is to apply the regression correction in a transformed, unconstrained parameter space [Beaumont, 2010].

For the “Laplace-Gaussian” and “Dep-Indep” examples, only independent bound constraints were applied on parameters, so that we used independent logistic transforms on each parameter, as in [Blum and François, 2010]. As the bounds are inclusive, to avoid infinite values, we used slightly inflated bounds, by removing 0.49 to the lower bound, and adding 0.49 to the upper bound. As all parameters are assumed to be integers for simulations, we then rounded the obtained back-transformed values, which ensured that all values lay between the bounds (bounds included).

For the “Human-like” and empirical Human datasets, in addition to bound constraints, relative order constraints were applied on splitting times, with $t4 > t3 > t2$. We dealt with this constraint using the logistic transformation of differences. More specifically, setting all the constraints together, we have:

$$t2 < t3 < t4; 100 \le t2 \le 998; 101 \le t3 \le 999; 103 \le t4 \le 1000.$$

The transform we applied is then:

$$ \begin{cases} u2 = t2 \\ u3 = \max(t3 - t2 - 1; \alpha) \\ u4 = \max(t4 - t3 - 1; \alpha), \end{cases} \quad \text{and} \quad \begin{cases} v2 = \logit(u2; 100,998) \\ v3 = \logit(u3; 0,998 - t2 + 0.49) \\ v4 = \logit(u4; 0,999 - t3 + 0.49), \end{cases} $$

where $\alpha = 0.25$ is chosen so that the lower bound of 0 on the difference is never attained, thus avoiding infinite values, and $\logit(\cdot; a, b)$ is the logit (inverse logistic) transform with upper and lower bounds $a$ and $b$. This transformation sends the constrained parameter values into $\mathbb{R}^3$, and can be inverted as:

$$ \begin{cases} u2 = \logit^{-1}(v2; 100,998) \\ u3 = \logit^{-1}(v3; 0,998 - t2 + 0.49) \\ u4 = \logit^{-1}(v4; 0,999 - t3 + 0.49), \end{cases} \quad \text{and} \quad \begin{cases} t2 = \text{round}(u2) \\ t3 = \text{round}(u3 + t2 + 1) \\ t4 = \text{round}(u4 + t3 + 1), \end{cases} $$
where the rounding operation ensure that we recover integers. Note that this is not a one to one map, as we are mapping a finite set (constrained integer parameters) into $ℝ³$, and it is only injective. Other approaches based on unconstrained parameters as in the “Dep-Indep” example could avoid the use of this ad-hoc transformation, that has known issues [Beaumont, 2010].

## C.2 Impact of *K* value on the prior GoF

Figure S2 shows the impact of the choice of *k* in the LOF and kNN scores on the power of the prior GoF test. For the kNN score (light orange), the horizontal line corresponding to *k* = 1 as in the main text almost always yields the highest power, especially in the more complex example (last column). Similarly, for the LOF score (black), the horizontal line corresponding to the “max-LOF” score for *k* between 5 and 20 as in the main text allows for powers that are close to the optimal over all *k* values. For both scores, the value of *k* has a relatively small impact, which justify the default values we used in the main text.

<figure><img src="image_11.png" /></figure>

Figure S2: **Estimated power for the prior GoF test**, at a level of 5%, for the three simulation settings. The LOF score is computed for *k* between 1 and 20 (black dots), and using the “max-LOF” score for *k* between 5 and 20 (horizontal black line). The kNN score is computed for *k* between 1 and 20 (light orange dots), and using a fixed *k* = 1 (horizontal light orange line) as in the main text for comparison. In all settings, *N*<sub>ref</sub> = *N*<sub>calib</sub>, and the total number of simulated particles is *N*<sub>ref</sub> + *N*<sub>calib</sub> (lines).
### C.3 Prior GoF p-value Calibration

Figure S3 shows the distribution of p-values of the prior GoF test when the PODs are simulated from the null model. It shows that the test is well calibrated, in all settings and for both the LOF and kNN scores, with the p-values being approximately uniformly distributed.

<figure><img src="image_12.png" /></figure>

Figure S3: **Distribution of p-values of the post-inference GoF test** when the PODs are simulated according to the null model, for the three simulation settings and the three ABC posterior estimation methods. The total number of simulated particles varies from 500 to 5000, with $N_{\text{ref}} = N_{\text{calib}}$ in all settings. The LOF score (black) is the “max-LOF” score for *k* between 5 and 20. The kNN score (light orange) is computed for *k* = 1.

### C.4 Impact of *K* value on the Post Inference GoF Power

Figure S4 shows the impact of the choice of *k* in the LOF and kNN scores on the power of the holdout GoF test. The results are similar to the ones we obtained for the prior GoF test, with the default value of *k* we chose having reasonable performances.
<figure><img src="image_13.png" /><figcaption>Figure S4: **Estimated power for the holdout GoF test**, at a level of 5%, for the three simulation settings. The posterior was localized with $\epsilon$ such that the number of posterior particles is equal to $N_{post} = 1000$, starting from 50 000 (top line) and 100 000 (bottom line) total number of simulated particles. The LOF score is computed for $k$ between 1 and 20 (black dots), and using the “max-LOF” score for $k$ between 5 and 20 (horizontal black line). The kNN score is computed for $k$ between 1 and 20 (light orange dots), and using a fixed $k = 1$ (horizontal light orange line) as in the main text for comparison. In all settings, $N'_{ref} = N'_{calib} = 500$.</figcaption></figure>

## C.5 Impact of the Posterior Estimation Method on the Post Inference GoF Power

Figure S5 shows the impact of the posterior estimation procedure on the power of the holdout GoF test. The three methods tested here yield similar powers.
<figure><img src="image_14.png" /><figcaption>Figure S5: **Estimated power for the post-inference holdout GoF test**, at a level of 5%, for the three simulation settings and the three ABC posterior estimation methods. The posterior was localized with $\epsilon$ such that the number of posterior particles is equal to $N_{\text{post}}$ (x axis), starting from $N_{\text{ref}} = 50\,000$ (full line) and $100\,000$ (dashed line) total number of simulated particles. This corresponds to 0.5% to 4% of the simulations, depending on $N_{\text{ref}}$. The LOF score (black) is the “max-LOF” score for $k$ between 5 and 20. The kNN score (light orange) is computed for $k = 1$. In all settings, and $N'_{\text{calib}} = N'_{\text{ref}} = N_{\text{post}}/2$.</figcaption></figure>

### C.6 Impact of the Posterior Estimation Method on the Post Inference GoF Calibration

Figure S6 shows the strong impact of the posterior estimation procedure on the calibration of the holdout GoF test. While simple corrections seem to yield correctly calibrated tests in simple datasets (first two rows), they do not seem to be sufficient in the complex example (last row).
<figure><img src="image_15.png" /><figcaption>Figure S6: Distribution of p-values of the post-inference GoF test when the PODs were simulated according to the null model, for the three simulation settings (Laplace-Gaussian, Dep-Indep and Human-like) and the three ABC posterior estimation methods (rejection, local linear or ridge regression). The total number $N_{\text{ref}}$ of simulated particles was taken to be equal to 50 000 (full lines) or 100 000 (dashed lines), and 1 000 points were taken in the posterior sample, with $N'_{\text{calib}} = N'_{\text{ref}} = N_{\text{post}}/2$. The LOF score (black) is the “max-LOF” score for $k$ between 5 and 20. The kNN score (light orange) is computed for $k=1$.</figcaption></figure>
D Human Population Genetics Example: The six Considered Evolutionary Scenarios

<figure><img src="image_16.png" /><figcaption>Figure S7: **Six scenarios of evolution of four modern Human populations.** The genotyped populations are Yoruba (YRI, Nigeria, Africa), Han (CHB, China, East Asia), British (GBR, England and Scotland, Europe), and Americans of African Ancestry (ASW, USA). The six scenarios differ from each other by one ancient and one recent historical event: (i) a single out-of-Africa colonization event giving an ancestral out-of-Africa population which secondarily split into one European and one East Asian population lineage (scenarios 1, 2 and 3), versus two independent out-of-Africa colonization events, one giving the European lineage and the other one giving the East Asian lineage (scenarios 4, 5 and 6). (ii) the possibility (or not; scenarios 1 and 4) of a recent genetic admixture of ASW individuals with their African ancestors and individuals of European (scenarios 2 and 5) or East Asia (scenarios 3 and 6) origins. The scenarios 2 and 3 are those considered in the third case study of our power simulation-based study (see Section 3.1.3).</figcaption></figure>
# E Human Population Genetics Example: Additional Analyses

## E.1 Human Population Genetics Example: Prior p-values

Figure S8 compares the confidence intervals obtained with the asymptotic and bootstrap methods (see Section 4.2 of the main text). Asymptotic intervals tend to be narrower than bootstrap intervals that, although more computationally intensive, might better represent the uncertainty around the point estimate of the p-value.

<figure><img src="image_17.png" /></figure>

Figure S8: Pre-inference GoF results on the real SNP dataset1 of modern Human populations for different simulation efforts using the LOF or the kNN outlier score. Estimated p-values for the prior GoF test, for the six scenarios of Figure S7. The LOF score (black) is the “max-LOF” score for *k* between 5 and 20. The kNN score (light orange) is computed for *k* = 1. In all settings, *N*<sub>ref</sub> = *N*<sub>calib</sub>, and the total number of simulated particles is *N*<sub>ref</sub> + *N*<sub>calib</sub> (x axis). Vertical bars represent the 95% confidence interval, either using the asymptotic formula (left column) or the Highest Density Intervals (HDI) computed from the bootstrap procedure with 500 replicates (right column). The 5% rejection threshold is shown as a dashed line. Only scenarios 2 and 3 are consistently not rejected, n particular for relatively low simulation costs (i.e. 1000 simulations per scenario).
E.2 Human Population Genetics Example: Post Inference p-values

Figures S9 and S10 show the p-values estimated for the holdout test applied on Scenario 2, using various posterior estimation methods, and either the asymptotic or bootstrap uncertainty estimation procedure. Although Scenario 2 gets consistently rejected when using the rejection method (first line), the result is less clear for the linear and ridge regression procedures. However, as shown in Figure 5 of the main text, using these posterior estimation methods tend to worsen the calibration of the test compared to the more simple rejection method, so that we are less confident in the results of these methods. Better posterior estimation procedures would be needed to circumvent this issue. As previously, asymptotic intervals tend to be narrower than bootstrap intervals

<figure><img src="image_18.png" /></figure>

Figure S9: **Post-inference GoF results on the real SNP dataset1 of modern Human populations for different simulation efforts using the LOF or the kNN outlier score.** Estimated p-values for the holdout GoF test for the selected scenario 2 of Figure S7, using the three ABC posterior estimation methods (rejection, local linear or ridge regression). The posterior was localized with $\epsilon$ such that the number of posterior particles is equal to $N_{post}$ (x axis), starting from 50 000 (left panel) and 100 000 (right panel) total number of simulated particles. In all settings, $N'_{calib} = N'_{ref} = N_{post}/2$. The LOF score (black) is the “max-LOF” score for $k$ between 5 and 20. The kNN score (light orange) is computed for $k = 1$. Vertical bars represent the 95% Highest Density Intervals (HDI) computed from the bootstrap procedure described in Section 4.2 with 500 replicates. The 5% rejection threshold is shown as a dashed line.
<figure><img src="image_19.png" /><figcaption>Figure S10: **Post-inference GoF results on the real SNP dataset1 of modern Human populations for different simulation efforts using the LOF or the kNN outlier score.** Estimated p-values for the holdout GoF test for the selected scenario 2 of Figure S7, using the three ABC posterior estimation methods (rejection, local linear or ridge regression). The posterior was localized with $\epsilon$ such that the number of posterior particles is equal to $N_{post}$ (x axis), starting from 50 000 (left panel) and 100 000 (right panel) total number of simulated particles. In all settings, $N'_{calib} = N'_{ref} = N_{post}/2$. The LOF score (black) is the “max-LOF” score for $k$ between 5 and 20. The kNN score (light orange) is computed for $k = 1$. Vertical bars represent the 95% Highest Density Intervals (HDI) computed from the asymptotic procedure described in Section 4.2. The 5% rejection threshold is shown as a dashed line.</figcaption></figure>
References

Mark A. Beaumont. Approximate Bayesian Computation in Evolution and Ecology. *Annual Review of Ecology, Evolution, and Systematics*, 41(1):379–406, 2010.

Mark A. Beaumont. Approximate Bayesian Computation. *Annual Review of Statistics and Its Application*, 6:379–403, 2019.

Mark A. Beaumont, Wenyang Zhang, and David J. Balding. Approximate Bayesian Computation in Population Genetics. *Genetics*, 162(4):2025–2035, 2002.

Yoav Benjamini and Yosef Hochberg. Controlling the False Discovery Rate: A Practical and Powerful Approach to Multiple Testing. *Journal of the Royal Statistical Society: Series B*, 57(1):289–300, 1995.

Michael G. B. Blum and Olivier François. Non-linear regression models for Approximate Bayesian Computation. *Statistics and Computing*, 20(1):63–73, 2010.

George E. P. Box. Sampling and Bayes' Inference in Scientific Modelling and Robustness. *Journal of the Royal Statistical Society: Series A*, 143(4):383, 1980.

Markus M. Breunig, Hans-Peter Kriegel, Raymond T. Ng, and Jörg Sander. LOF: Identifying Density-Based Local Outliers. *ACM SIGMOD Record*, 29(2):93–104, 2000.

Katarzyna Bryc, Eric Y. Durand, J. Michael Macpherson, David Reich, and Joanna L. Mountain. The Genetic Ancestry of African Americans, Latinos, and European Americans across the United States. *The American Journal of Human Genetics*, 96(1):37–53, 2015.

Domagoj Cevid, Loris Michel, Jeffrey Näf, Peter Bühlmann, and Nicolai Meinshausen. Distributional random forests: Heterogeneity adjustment and multivariate distributional regression. *Journal of Machine Learning Research*, 23(333):1–79, 2022.

François-David Collin, Ghislain Durif, Louis Raynal, Eric Lombaert, Mathieu Gautier, Renaud Vitalis, Jean-Michel Marin, and Arnaud Estoup. Extending approximate bayesian computation with supervised machine learning to infer demographic history from genetic polymorphisms using DIYABC random forest. *Molecular Ecology Resources*, 21(8):2598–2613, 2021.

Jean-Marie Cornuet, Pierre Pudlo, Julien Veyssier, Alexandre Dehne-Garcia, Mathieu Gautier, Raphaël Leblois, Jean-Michel Marin, and Arnaud Estoup. DIYABC v2.0: a software to make approximate Bayesian computation inferences about population history using single nucleotide polymorphism, DNA sequence and microsatellite data. *Bioinformatics*, 30(8):1187–1189, 2014.

Kyle Cranmer, Johann Brehmer, and Gilles Louppe. The frontier of simulation-based inference. *Proceedings of the National Academy of Sciences*, 117(48):30055–30062, 2020.

Dirk Eddelbuettel and Romain François. Rcpp: Seamless random integration. *Journal of Statistical Software*, 40(8), 2011.

Paul Fearnhead and Dennis Prangle. Constructing summary statistics for approximate Bayesian computation: semi-automatic approximate Bayesian computation. *Journal of the Royal Statistical Society: Series B*, 74(3):419–474, 2012.
Florence Forbes, Hien Duy Nguyen, TrungTin Nguyen, and Julyan Arbel. Summary statistics and discrepancy measures for approximate bayesian computation via surrogate posteriors. *Statistics and Computing*, 32(5), 2022.

David T. Frazier, Christian P. Robert, and Judith Rousseau. Model misspecification in approximate bayesian computation: consequences and diagnostics. *Journal of the Royal Statistical Society: Series B*, 82(2): 421–444, 2020.

Andrew Gelman, Xiao-Li Meng, and Hal Stern. Posterior Predictive Assessment of Model Fitness Via Realized Discrepancies. *Statistica Sinica*, 6(4):733–760, 1996.

Markus Goldstein and Seiichi Uchida. A Comparative Evaluation of Unsupervised Anomaly Detection Algorithms for Multivariate Data. *PLOS ONE*, 11(4):e0152173, 2016.

David Greenberg, Marcel Nonnenmacher, and Jakob Macke. Automatic posterior transformation for likelihood-free inference. In Kamalika Chaudhuri and Ruslan Salakhutdinov, editors, *Proceedings of the 36th International Conference on Machine Learning*, volume 97 of *Proceedings of Machine Learning Research*, pages 2404–2414, 2019.

Irwin Guttman. The Use of the Concept of a Future Observation in Goodness-Of-Fit Problems. *Journal of the Royal Statistical Society: Series B*, 29(1):83–100, 1967.

Michael Hahsler, Matthew Piekenbrock, and Derek Doran. dbscan: Fast density-based clustering with R. *Journal of Statistical Software*, 91(1), 2019.

Songqiao Han, Xiyang Hu, Hailiang Huang, Minqi Jiang, and Yue Zhao. Adbench: Anomaly detection benchmark. In S. Koyejo, S. Mohamed, A. Agarwal, D. Belgrave, K. Cho, and A. Oh, editors, *Advances in Neural Information Processing Systems*, pages 32142–32159. Curran Associates, Inc., 2022.

Joeri Hermans, Volodimir Begy, and Gilles Louppe. Likelihood-free MCMC with amortized approximate ratio estimators. In Hal Daumé III and Aarti Singh, editors, *Proceedings of the 37th International Conference on Machine Learning*, volume 119 of *Proceedings of Machine Learning Research*, pages 4239–4248, 2020.

J. R. M. Hosking. L-moments: Analysis and estimation of distributions using linear combinations of order statistics. *Journal of the Royal Statistical Society: Series B*, 52(1):105–124, 1990.

J. R. M. Hosking. *L-Moments*, 2024. R package, version 3.2.

Rob J. Hyndman. Computing and Graphing Highest Density Regions. *The American Statistician*, 50(2): 120–126, 1996.

John K. Kruschke. *Doing Bayesian data analysis*. Elsevier, Academic Press, Amsterdam, 2015.

Louisiane Lemaire, Flora Jay, I-Hung Lee, Katalin Csilléry, and Michael G. B. Blum. Goodness-of-fit statistics for approximate bayesian computation. *arXiv*, 1601.04096, 2016.

Kirk E. Lohmueller and Rasmus Nielsen, editors. *Human Population Genomics: Introduction to Essential Concepts and Applications*. Springer International Publishing, Cham, Switzerland, 2021.
Jan-Matthis Lueckmann, Pedro J Gonçalves, Giacomo Bassetto, Kaan Öcal, Marcel Nonnenmacher, and Jakob H Macke. Flexible statistical inference for mechanistic models of neural dynamics. In I. Guyon, U. Von Luxburg, S. Bengio, H. Wallach, R. Fergus, S. Vishwanathan, and R. Garnett, editors, *Advances in Neural Information Processing Systems*, volume 30. Curran Associates, Inc., 2017.

Jan-Matthis Lueckmann, Jan Boelts, David S. Greenberg, Pedro J. Gonçalves, and Jakob H. Macke. Benchmarking simulation-based inference. In Arindam Banerjee and Kenji Fukumizu, editors, *Proceedings of the 24th International Conference on Artificial Intelligence and Statistics*, volume 130 of *Proceedings of Machine Learning Research*, pages 343–351, 2021.

Saioa López, Lucy Van Dorp, and Garrett Hellenthal. Human dispersal out of Africa: A lasting debate. *Evolutionary Bioinformatics*, pages 57–68, 2015.

Jean-Michel Marin, Pierre Pudlo, Christian P. Robert, and Robin J. Ryder. Approximate bayesian computational methods. *Statistics and Computing*, 22(6):1167–1180, 2012.

Gemma E. Moran, David M. Blei, and Rajesh Ranganath. Holdout predictive checks for bayesian model criticism. *Journal of the Royal Statistical Society: Series B*, 86(1):194–214, 2023.

George Papamakarios and Iain Murray. Fast $\varepsilon$-free inference of simulation models with bayesian conditional density estimation. In D. Lee, M. Sugiyama, U. Luxburg, I. Guyon, and R. Garnett, editors, *Advances in Neural Information Processing Systems 29*, pages 1028–1036. Curran Associates, Inc., 2016.

George Papamakarios, Theo Pavlakou, and Iain Murray. Masked autoregressive flow for density estimation. In I. Guyon, U. Von Luxburg, S. Bengio, H. Wallach, R. Fergus, S. Vishwanathan, and R. Garnett, editors, *Advances in Neural Information Processing Systems*, volume 30. Curran Associates, Inc., 2017.

George Papamakarios, David Street, and Iain Murray. Sequential neural likelihood: Fast likelihood-free inference with autoregressive flows. In Kamalika Chaudhuri and Masashi Sugiyama, editors, *Proceedings of the 22nd International Conference on Artificial Intelligence and Statistics*, volume 89 of *Proceedings of Machine Learning Research*, pages 837–848, 2019.

J K Pritchard, M T Seielstad, A Perez-Lezaun, and M W Feldman. Population growth of human Y chromosomes: a study of Y chromosome microsatellites. *Molecular Biology and Evolution*, 16(12):1791–1798, 1999.

Pierre Pudlo, Jean-Michel Marin, Arnaud Estoup, Jean-Marie Cornuet, Mathieu Gautier, and Christian P. Robert. Reliable ABC model choice via random forests. *Bioinformatics*, 32:859–866, 2016.

Oliver Ratmann, Christophe Andrieu, Carsten Wiuf, and Sylvia Richardson. Model criticism based on likelihood-free inference, with an application to protein network evolution. *Proceedings of the National Academy of Sciences of the United States of America*, 106(26):10576–10581, 2009.

Louis Raynal, Jean-Michel Marin, Pierre Pudlo, Mathieu Ribatet, Christian P. Robert, and Arnaud Estoup. ABC random forests for Bayesian parameter inference. *Bioinformatics*, 35(10):1720–1728, 2018.

Jeffrey I. Rose, Vitaly I. Usik, Anthony E. Marks, Yamandu H. Hilbert, Christopher S. Galletti, Ash Parton, Jean Marie Geiling, Viktor Černý, Mike W. Morley, and Richard G. Roberts. The nubian complex of dhofar, oman: An african middle stone age industry in southern arabia. *PLOS ONE*, 6(11):e28239, 2011.
Donald B. Rubin. Bayesianly Justifiable and Relevant Frequency Calculations for the Applied Statistician. *The Annals of Statistics*, 12(4):1151–1172, 1984.

Sara Sheehan and Yun S. Song. Deep learning for population genetic inference. *PLOS Computational Biology*, 12(3):1–28, 2016.

Simon Tavaré, David J. Balding, R. C. Griffiths, and Peter Donnelly. Inferring coalescence times from dna sequence data. *Genetics*, 145:505–518, 1997.

The 1000 Genomes Project Consortium. An integrated map of genetic variation from 1,092 human genomes. *Nature*, 491:56–65, 2012.

A. W. van der Vaart. *Asymptotic Statistics*. Cambridge University Press, 1998.

Daniel Ward, Patrick Cannon, Mark Beaumont, Matteo Fasiolo, and Sebastian Schmon. Robust Neural Posterior Estimation and Statistical Model Criticism. In S. Koyejo, S. Mohamed, A. Agarwal, D. Belgrave, K. Cho, and A. Oh, editors, *Advances in Neural Information Processing Systems*, volume 35, pages 33845–33859. Curran Associates, Inc., 2022.

Steve Weston. *doParallel: Foreach parallel adaptor for the parallel package*, 2022. R package, version 1.0.17.

Janis E. Wigginton, David J. Cutler, and Gonçalo R. Abecasis. A note on exact tests of hardy-weinberg equilibrium. *The American Journal of Human Genetics*, 76(5):887–893, 2005.

Richard David Wilkinson. Approximate bayesian computation (abc) gives exact results under the assumption of model error. *Statistical Applications in Genetics and Molecular Biology*, 12(2):129–141, 2013.

Thomas W. Yee. The vgam package for categorical data analysis. *Journal of Statistical Software*, 32(10): 1–34, 2010.