## Abstract

Until recently, the use of Bayesian inference was limited to a few cases because for many realistic probability models the likelihood function cannot be calculated analytically. The situation changed with the advent of likelihood-free inference algorithms, often subsumed under the term approximate Bayesian computation (ABC). A key innovation was the use of a postsampling regression adjustment, allowing larger tolerance values and as such shifting computation time to realistic orders of magnitude. Here we propose a reformulation of the regression adjustment in terms of a general linear model (GLM). This allows the integration into the sound theoretical framework of Bayesian statistics and the use of its methods, including model selection via Bayes factors. We then apply the proposed methodology to the question of population subdivision among western chimpanzees, *Pan troglodytes verus*.

WITH the advent of ever more powerful computers and the refinement of algorithms like MCMC or Gibbs sampling, Bayesian statistics have become an important tool for scientific inference during the past two decades. Consider a model creating data (DNA sequence data, for example) determined by parameters from some (bounded) parameter space Π ⊂ ℝ* ^{m}* whose joint prior density we denote by . The quantity of interest is the posterior distribution of the parameters, which can be calculated by Bayes rule aswhere is the likelihood of the data and is a normalizing constant. Direct use of this formula, however, is often prevented by the fact that the likelihood function cannot be calculated analytically for many realistic probability models. In these cases one is obliged to use stochastic simulation. Tavaré

*et al.*(1997) propose a rejection sampling method for simulating a posterior random sample where the full data are replaced by a summary statistic

*s*(like the number of segregating sites in their setting). Even if the statistic does not capture the full information contained in the data , rejection sampling allows for the simulation of approximate posterior distributions of the parameters in question (the scaled mutation rate in their model). This approach was extended to multiple-parameter models with multivariate summary statistics by Weiss and von Haeseler (1998). In their setting a candidate vector of parameters is simulated from a prior distribution and is accepted if its corresponding vector of summary statistics is sufficiently close to the observed summary statistics

**s**

_{obs}with respect to some metric in the space of

**s**,

*i.e.*, if dist(

**s**,

**s**

_{obs}) < ϵ for a fixed tolerance ϵ. We suppose that the likelihood of the full model is continuous and nonzero around

**s**

_{obs}. In practice the summary statistics are often discrete but the range of values is large enough to be approximated by real numbers. The likelihood of the truncated model obtained by this acceptance–rejection process is given by(1)where is the ϵ-ball in the space of summary statistics and Ind(·) is the indicator function. Observe that degenerates to a (Dirac) point measure centered at

**s**

_{obs}as . If the parameters are generated from a prior , then the distribution of the parameters retained after the rejection process outlined above is given by(2)

We call this density the *truncated prior*. Combining (1) and (2) we get(3)Thus the posterior distribution of the parameters under the model for **s** = **s**_{obs} given the prior is exactly equal to the posterior distribution under the truncated model given the truncated prior . If we can estimate the truncated prior and make an educated guess for a parametric statistical model of *M*_{ϵ}(**s**_{obs}), we arrive at a reasonable approximation of the posterior even if the likelihood of the full model is unknown. It is to be expected that due to the localization process the truncated model will exhibit a simpler structure than the full model and thus be easier to estimate.

Estimating is straightforward, at least when the summary statistics can be sampled from in a reasonable amount of time: Sample the parameters from the prior , create their respective statistics **s** from , and save those parameters whose statistics lie in in a list . The empirical distribution of these retained parameters yields an estimate of . If the tolerance ϵ is small, then one can assume that is close to some (unknown) constant over the whole range of . Under that assumption, Equation 3 shows that . However, when the dimension *n* of summary statistics is high (and for more complex models dimensions like *n* = 50 are not unusual), the “curse of dimensionality” implies that the tolerance must be chosen rather large or else the acceptance rate becomes prohibitively low. This, however, distorts the precision of the approximation of the posterior distribution by the truncated prior (see Wegmann *et al.* 2009). This situation can be partially alleviated by speeding up the sampling process; such methods are subsumed under the term *approximate Bayesian computation* (ABC). Marjoram *et al.* (2003) develop a variant of the classical Metropolis–Hastings algorithm (termed ABC–MCMC in Sisson *et al.* 2007), which allows them to sample directly from the truncated prior . In Sisson *et al.* (2007) a sequential Monte Carlo sampler is proposed, requiring substantially less iterations than ABC–MCMC. But even when such methods are applied, the assumption that is constant over the ϵ-ball is a very rough one, indeed.

To take into account the variation of within the ϵ-ball, a postsampling regression adjustment (termed ABC-REG in the following) of the sample of retained parameters is introduced in the important article by Beaumont *et al.* (2002). Basically, they postulate a (locally) linear dependence between the parameters and their associated summary statistics **s**. More precisely, the (local) model they implicitly assume is of the form , where **M** is a matrix of regression coefficients, **m**_{0} a constant vector, and a random vector of zero mean. Computer simulations suggest that for many population models ABC–REG yields posterior marginal densities that have narrower highest posterior density (HPD) regions and are more closely centered around the true parameter values than the empirical posterior densities directly produced by ABC samplers (Wegmann *et al.* 2009). An attractive feature of ABC–REG is that the posterior adjustment is performed directly on the simulated parameters, which makes estimation of the marginal posteriors of individual parameters particularly easy. The method can also be extended to more complex, nonlinear models as demonstrated, *e.g.*, in Blum and Francois (2009). In extreme situations, however, ABC–REG may yield posteriors that are nonzero in parameter regions where the priors actually vanish (see Figure 1B for an illustration of this phenomenon). Moreover, it is not clear how ABC–REG could yield an estimate of the marginal density of model at **s**_{obs}, information that is useful for model comparison.

In contrast to ABC–REG we treat the parameters as exogenous and the summary statistics **s** as endogenous variables and we stipulate for a general linear model (GLM in the literature—not to be confused with the generalized linear models that unfortunately share the same abbreviation). To be precise, we assume the summary statistics **s** created by the truncated model's likelihood to satisfy(4)where **C** is a *n* × *m* matrix of constants, **c**_{0} an *n* × 1 vector, and a random vector with a multivariate normal distribution of zero mean and covariance matrix :

A GLM has the advantage of taking into account not only the (local) linearity, but also the strong correlation normally present between the components of the summary statistics. Of course, the model assumption (4) can never represent the full truth since its statistics are in principle unbounded whereas the likelihood is supported on the ϵ-ball around **s**_{obs}. But since the multivariate Gaussians will fall off rapidly in practice and not reach far out off the boundary of , this is a disadvantage we can live with. In particular, the ordinary least squares (OLS) estimate outlined below implies that for the constant **c**_{0} tends to **s**_{obs} whereas the design matrix **C** and the covariance matrix both vanish. This means that in the limit of zero tolerance ϵ = 0 our model assumption yields the true posterior distribution of .

## THEORY

In this section we describe the above methodology—referred to as ABC–GLM in the following—in more detail. The basic two-step procedure of ABC–GLM may be summarized as follows.

#### GLM1:

Given a model creating summary statistics **s** and given a value of observed summary statistics **s**_{obs}, create a sample of retained parameters , with the aid of some ABC sampler (rejection sampling, ABC–MCMC, or ABC–PRC) based on a prior distribution and some choice of the tolerance ϵ > 0.

#### GLM2:

Estimate the truncated model as a general linear model and determine, on the basis of the sample , from the truncated prior an approximation to the posterior according to Equation 3.

Let us look more closely at these two steps.

##### GLM1: ABC sampling:

We refer the reader to Marjoram *et al.* (2003) and Sisson *et al.* (2007) for details concerning ABC algorithms and to Marjoram and Tavaré (2006) for a comprehensive review of computational methods for genetic data analysis. In practice, the dimension of the summary statistics is often reduced by a principal components analysis (PCA). PCA also has a certain decorrelation effect. A more sophisticated method of reducing the dimension of summary statistics, based on partial least squares (PLS), is described in Wegmann *et al.* (2009). In a recent preprint, Vogl *et al.* (C. Vogl, C. Futschik and C. Schloetterer, unpublished data) propose a Box–Cox-type transformation of the summary statistics that makes the likelihood close to multivariate Gaussian. This transformation might be especially efficient in our context as we assume normality of the error terms in our model assumption.

To fix the notation, let be a sample of vector-valued parameters created by some ABC algorithm simulating from some prior and be the sample of associated summary statistics produced by the model . Each parameter is an *m*-dimensional column vector and each summary statistic is an *n*-dimensional column vector . The samples and can thus be viewed as *m* × *N* and *n* × *N* matrices **P** and **S**, respectively.

The empirical estimate of the truncated prior is given by the discrete distribution that puts a point mass of 1/*N* on each value . We smooth out this empirical distribution by placing a sharp Gaussian peak over each parameter value . More precisely, we set(5)whereandis the covariance matrix of ϕ that determines the width of the Gaussian peaks. The larger the number *N* of sampled parameter values is, the sharper the peaks can be chosen to still get a rather smooth π_{ϵ}. If the parameter domain Π is normalized to [0, 1]* ^{m}*, say, then a reasonable choice is σ

_{k}= 1/

*N*. Otherwise, σ

_{k}should be adapted to the parameter range of the parameter component θ

_{k}. Too small values of σ

_{k}will result in wiggly posterior curves, and too large values might unduly smear out the curves. The best advice is to run the calculations with several choices for . If π

_{ϵ}induces a correlation between parameters, a nondiagonal might be beneficial. In practice, however, the posterior estimates are most sensitive to the diagonal values of .

##### GLM2: general linear model:

As explained in the Introduction, we assume the truncated model to be normal linear; *i.e.*, the random vectors **s** satisfy (4). The covariance matrix encapsulates the strong correlations normally present between the components of the summary statistics. **C**, **c**_{0}, and can be estimated by standard multivariate regression analysis (OLS) from the sample , created in step GLM1. [Strictly speaking, one must redo an ABC sample from uniform priors over Π to get an unbiased estimate of the GLM if the prior is not uniform already. On the other hand, ordinary least-squares estimators are quite insensitive to the prior's influence. In practice, one can as well use the sample to do the estimate. We applied both estimation methods to the toy models presented in the examples from population genetics section and found no significant difference between the estimated posteriors. The same holds true for the so-called feasible generalized least-squares (FGLS) estimator; see Greene (2003). In this two-stage algorithm the covariance matrix is first estimated as in our setting but in a second round the design matrix **C** is newly estimated. When we applied FGLS to our toy models, we found a difference in the estimated matrices only after the eighth significant decimal. FGLS is a more efficient estimator only when the sample sizes are relatively small as is often the case in economical data sets but not in ABC situations. In theory, both OLS and FGLS are consistent estimators but FGLS is more efficient.] To be specific, set **X** = (**1**⋮**P*** ^{t}*), where

**1**is an

*N*× 1 vector of 1's.

**C**and

**c**

_{0}are determined by the usual least-squares estimatorand for we have the estimate(6)where are the residuals. The likelihood for this model—dropping the hats on the matrices to unburden the notation—is given by(7)

An exhaustive treatment of linear models in a Bayesian (econometric) context is given in Zellner's book (Zellner 1971).

Recall from (3) that for a prior and an observed summary statistic **s**_{obs}, the parameter's posterior distribution for our full model is given by(8)where is the likelihood of the truncated model given by (7) and is the estimated (and smoothed) truncated prior given by (5).

Performing some matrix algebra (see appendix a), one can show that the posterior (8) is—up to a multiplicative constant—of the form , whereHere **T**, **t*** ^{j}*, and

**v**

*are given by(9)and*

^{j}**t**

*=*

^{j}**Tv**

*, where(10)*

^{j}From this we get(11)where(12)

When the number of parameters exceeds two, graphical visualization of the posterior distribution becomes impractical and marginal distributions must be calculated. The marginal posterior density of the parameter θ_{k} is defined bywhere integration is performed along all parameters except θ_{k}.

Recall that the marginal distribution of a multivariate normal with respect to the *k*th component is the univariate normal density . Using this fact, it is not hard to show that the marginal posterior of parameter θ_{k} is given by(13)where τ_{k,k} is the *k*th diagonal element of the matrix **T**, *t _{k}^{j}* is the

*k*th component of the vector

**t**

*, and is still determined according to (12). The normalizing constant*

^{j}*a*could, in principle, be determined analytically but is in practice more easily recovered by a numerical integration. Strictly speaking, the integration should be done only over the bounded parameter domain Π and not over the whole of ℝ

*. But this no longer allows for an analytic form of the marginal posterior distribution. For large values of*

^{m}*N*the diagonal elements in the matrix Σ

_{θ}can be chosen so small that the error is in any case negligible.

#### Model selection:

The principal difficulty of model selection methods in nonparametric settings is that it is nearly impossible to estimate the likelihood of at **s**_{obs} due to the high dimension of the summary statistics (curse of dimensionality); see Beaumont (2007) for an approach based on multinomial logit. Parametric models on the other hand lend themselves readily to model selection via Bayes factors. Given the model , one must determine the marginal densityIt is easy to check from (1) and (2) that

Here(14)is the acceptance rate *p* of the rejection process. It can easily be estimated with aid of ABC–REJ: Sample parameters from the prior create the corresponding statistics **s** from and count what fraction of the statistics fall into the ϵ-ball centered at **s**_{obs}.

If we assume the underlying model of to be our GLM, then the marginal density of at **s**_{obs} can be estimated as(15)where the sum runs over the parameter sample ,and

For two models and with prior probabilities π_{A} and π_{B} = 1 – π_{A}, the Bayes factor *B _{AB}* in favor of model over model is(16)where the marginal densities and are calculated according to (15). The posterior probability of model is

## EXAMPLES FROM POPULATION GENETICS

#### Toy models:

In Figure 1 we present the comparison of posteriors obtained with rejection sampling, ABC–REG and ABC–GLM, with those determined analytically (“true posteriors”). As a toy model we inferred the population–mutation parameter θ = 4*N*μ of a panmictic population model from the number of segregating sites *S* of a sample of sequences with 10,000 bp for different observed values and tolerance levels. Estimations are always based on 5000 simulations with dist(*S*, *S*_{obs}) < ϵ, and we report the average of 25 independent replications per data point. Estimation bias of the different approaches was assessed by computing the total variation distance between the inferred posterior and the true one obtained from analytical calculations using the likelihood function introduced by Watterson (1975). Recall that the *L*_{1}-distance of two densities *f*(θ) and *g*(θ) is given byIt is equal to 1 when *f* and *g* have disjoint supports and it vanishes when the functions are identical.

When we used a uniform prior θ ∼ Unif([0.005, 10]) (Figure 1, A–C), both ABC–REG and ABC–GLM give comparable results and improve the posterior estimation compared to the simple rejection algorithm except for very low tolerance values ϵ where the rejection algorithm is expected to be very close to the true posterior. The average total variation distances over all observed data sets and tolerance values ϵ are 0.236, 0.130, and 0.091 for the rejection algorithm, ABC–REG, and ABC–GLM, respectively. Note that perfect matches between the approximate and the true posteriors are difficult to obtain because all approximate posteriors depend on a smoothing step that may not give accurate results close to borders of their supports. However, when we used a discontinuous prior θ ∼ Unif([0.005, 3] ∪ [6, 10]) with an admittedly extremely artificial “gap” in the middle, we observed a quite distinct pattern (Figure 1, D and E). One clearly recognizes that posteriors inferred with ABC–REG are frequently misplaced and often even farther away from the true posterior (in total variation distance) than the prior, especially for cases where the likelihood of the observed data is maximal within the gap. The reason for this is that in the regression step of ABC–REG parameter values may easily be shifted outside the prior support. This behavior of ABC–REG has been observed earlier (Beaumont *et al.* 2002; Estoup *et al.* 2004; Tallmon *et al.* 2004) and as an *ad hoc* solution Hamilton *et al.* (2006) proposed to transform the parameter values prior to the regression step by a transformation of the form , where *a* and *b* are the lower and upper borders of the prior support interval. For more complex priors—like the discontinuous prior used here—this transformation may not work. ABC–GLM is much less affected by the gap prior than ABC–REG. The average total variation distances over all observed data sets and tolerance values ϵ are 0.221, 0.246, and 0.094 for the rejection algorithm, ABC–REG, and ABC–GLM, respectively. Example posteriors with *S*_{obs} = 16 based on 5000 simulations with dist(*S*, *S*_{obs}) < 10 are shown in Figure 2.

The success of ABC–GLM depends on how well a general linear model fits the truncated model . Under the null hypothesis that the fit is perfect the estimated residuals **r**_{j} (see Equation 6) are independently multivariate normally distributed random vectors. Hence the Mahalanobis distances(17)follow a χ^{2}-distribution with *n* degrees of freedom. As a quantification of model assessment we propose to report the Kolmogorov–Smirnov test statistic for the empirical distribution of *d _{j}* and the reference χ

^{2}-distribution. (Reporting

*P*-values will be of little use in practice since the null hypothesis does never hold exactly and hence the

*P*-values will become very small due to the large sample size.)

When the summary statistics are created from a general linear model, the fit should be optimal. This is indeed the case as the simulation results in Table 1 show. We performed 200 simulations of randomly created general linear models with *m* = 3 parameters, *n* = 4 summary statistics, and a multivariate normal prior. The observed statistics were also created from the respective models. For each simulated observed statistic and different acceptance rates *p* = 1.00, 0.50, 0.10, 0.05, and 0.01 we calculated the approximate posterior distributions π_{ϵ}, π_{REG}, and π_{GLM} for the rejection algorithm, ABC–REG, and ABC–GLM, respectively. As the prior is multivariate normal, the true posterior π_{0} can be analytically determined. Table 1 contains the means and standard deviations over the 200 simulations of the total variation distances of the approximate posteriors to the true posterior π_{0} as well as the mean and standard deviations of the Kolmogorov–Smirnov test statistics for the GLM model fit. As is expected, the model fit is perfect [*i.e.*, the Kolmogorov–Smirnov (KS) statistic is close to 0] for acceptance rate *p* = 1. As the acceptance rate becomes lower, the model fit deteriorates since the truncated model of a GLM is no longer exactly a general linear model. The total variation distance to the true posterior increases slightly as *p* gets smaller but the improved rejection posterior π_{ϵ} mostly outbalances the poorer model fit. As is expected in this ideal situation, ABC–GLM and ABC–REG substantially improve the posterior estimation over the pure rejection prior.

To test the other extreme we performed 200 simulations for a nonlinear one-parameter model with uniformly rather than normally distributed error terms; the prior was again a normal distribution. (The details of this toy model are described in appendix b.) As Table 2 shows, the GLM model fit is already poor for an acceptance rate of *p* = 1.00 (KS statistic ∼0.10) and further deteriorates as *p* decreases. Note that the approximate posteriors π_{REG} and π_{GLM} are closer to the true posterior in average than π_{ϵ} and that both adjustment methods perform similarly. As expected, the accuracy of the posteriors increases with smaller acceptance rates, despite the fact that the model fit within the ϵ-ball decreases. This suggests that the rejection step contributes substantially to the estimation accuracy, especially when the true model is nonlinear. We should mention that in ∼30% of the simulations both ABC–GLM and ABC–REG actually increased the distance to the true posterior in comparison to the rejection posterior π_{ϵ}. As a rule of thumb we suggest that posterior adjustments obtained by ABC–GLM or ABC–REG should not be trusted without further validation if the Kolmogorov–Smirnov statistic for the GLM model fit exceeds a value of, say, 0.10. In that case linear models are not sufficiently flexible to account for effects like nonlinearity in the parameters and nonnormality and heteroscedasticity in the error terms. In the setting of ABC–REG a wider class of models is introduced in Blum and Francois (2009), where machine-learning algorithms are applied for the parameter estimations. Whether these extensions can be applied in our context remains to be seen. The advantage of the general linear model is that estimations can be done with ordinary least squares and the important quantities like marginal posteriors and marginal likelihoods can be obtained analytically. For more complex models these quantities will probably be accessible only via numerical integration, Monte Carlo methods, etc.

#### Application to chimpanzees:

In standard taxonomies, chimpanzees, the closest living relatives of humans, are classified into two species: the common chimpanzee (*Pan troglodytes*) and the bonobo (*P. paniscus*). Both species are restricted to Africa and diverged ∼9 MYA (Won and Hey 2005; Becquet and Przeworski 2007). The common chimpanzees are further subdivided into three large populations or subspecies on the basis of their separation by geographic barriers. Among them, the western chimpanzees (*P. troglodytes verus*) form the most remote group. Interestingly, recent multilocus studies found consistent levels of gene flow between the western and the central (*P. t. troglodytes*) chimpanzees (Won and Hey 2005; Becquet and Przeworski 2007). Nonetheless, a recent study of 310 microsatellites in 84 common chimpanzees supports a clear distinction between the previously labeled populations (Becquet *et al.* 2007). Using a PCA analysis, indication for substructure within the western chimpanzees was found in the same study.

To demonstrate the applicability of the model selection given in the theory section we contrast two different models of the western chimpanzee population with this data set: a model of a single panmictic population with constant size and a finite island model of constant size and constant migration among demes. While we estimated θ = 2*N*_{e}μ, priors were set on *N*_{e} and μ separately with log_{10}(*N*_{e}) ∼ Unif([3, 5]) and μ ∼ *N*(5 × 10^{−4}, 2 × 10^{−4}) truncated on μ ∈ [10^{−4}, 10^{−3}]. In the case of the finite island model, we had an additional prior *n*_{pop} ∼ Unif([10, 100]) on the number of islands, and individuals were attributed randomly to the different islands.

We obtained genotypes for all 50 individuals reported to be of western chimpanzee origin from the study of Becquet *et al.* (2007), excluding captive-born hybrids. We checked the proposed (Becquet *et al.* 2007) mutation pattern for each individual locus, and all alleles not matching the assumed stepwise mutation model were set as missing data. A total of 265 loci were used, after removing the loci on the X and the Y chromosome as well as those being monomorphic among the western chimpanzees. All simulations were performed using the software SIMCOAL2 (Laval and Excoffier 2004) and we reproduced the pattern of missing data observed in the data set. Using the software package Arlequin3.0 (Excoffier *et al.* 2005), we calculated two summary statistics on the data set: the average number of alleles per locus, *K*, and *F*_{IS}, the fixation index within the western chimpanzees. We performed a total of 100,000 simulations per model.

In Figure 3 we report the Bayes factor of the island model according to (16) for different acceptance rates *A*_{ϵ}; see (14). While there is a large variation for very small acceptance rates, the Bayes factor stabilizes for *A*_{ϵ} ≥ 0.005. Note that *A*_{ϵ} ≤ 0.005 corresponds to <500 simulations and that the ABC–GLM approach, based on a model estimation and a smoothing step, is expected to produce poor results since the estimation of the model parameters is unreliable due to the small sample size. The good news is that the Bayes factor is stable over a large range of tolerance values. We may therefore safely reject the panmictic population model in favor of population subdivision among western chimpanzees with a Bayes factor of *B* ≈ 10^{5}.

## DISCUSSION

Due to still increasing computational power it is nowadays possible to tackle estimation problems in a Bayesian framework for which analytical calculation of the likelihood is inhibited. In such cases, approximate Bayesian computation is often the choice. A key innovation in speeding up such algorithms was the use of a regression adjustment, termed ABC–REG in this article, which used the frequently present linear relationship between generated summary statistics **s** and parameters of the model in a neighborhood of the observed summary statistics **s**_{obs} (Beaumont *et al.* 2002). The main advantage is that larger tolerance values ϵ still allow us to extract reasonable information about the posterior distribution and hence less simulations are required to estimate the posterior density.

Here we present a new approach to estimate approximate posterior distributions, termed ABC–GLM, similar in spirit to ABC–REG, but with two major advantages: First, by using a GLM to estimate the likelihood function, ABC–GLM is always consistent with the prior distribution. Second, while we do not find the ABC–GLM approach to substantially outperform ABC–REG in standard situations, it is naturally embedded into a standard Bayesian framework, which in turn allows the application of well-known Bayesian methodologies such as model averaging or model selection via Bayes factors. Our simulations show that the rejection step is especially beneficial if the true model is nonlinear for both ABC approaches. ABC–GLM is further compatible with any type of ABC sampler, including likelihood-free MCMC (Marjoram *et al.* 2003) or population Monte Carlo (Beaumont *et al.* 2009). Also, more complicated regression regimes taking nonlinearity or heteroscedacity into account may be envisioned when the GLM is replaced by some more complex model. A great advantage of the current GLM setting is its simplicity, which renders implementation in standard statistical packages feasible.

We showed the applicability of the model selection procedure via Bayes factors by opposing two different models of population structure among the western chimpanzees *P. troglodytes verus*. Our analysis strongly suggests population substructure within the western chimpanzees since an island model is significantly favored over a model of a panmictic population. While none of our simple models is thought to mimic the real setting exactly, we still believe that they capture the main characteristics of the demographic history influencing our summary statistics, namely the number of alleles *K* and the fixation index *F*_{IS}. While the observed *F*_{IS} of 2.6% has been attributed to inbreeding previously (Becquet *et al.* 2007), we propose that such values may easily arise if diploid individuals are sampled in a randomly scattered way over a large, substructured population. While it was almost impossible to simulate the value *F*_{IS} = 2.6% in the model of a panmictic population, it easily falls within the range of values obtained from an island model.

## APPENDIX A: PROOFS OF THE MAIN FORMULAS

To keep this article self-contained, we present a proof of formulas (11) and (15). The argument is an adaptation from the proof of Lemma 1 in Lindley and Smith (1972). By linearity it clearly suffices to show the formulas for one fixed sampled parameter . The results then follow.

Theorem. *Suppose that, given the parameter vector* , *the distribution of the statistics vector* **s** *is multivariate normal*,*and*, *given the fixed parameter vector* , *the distribution of the parameter* *is**Then:*

1. *The distribution of* given **s** *is**where* *and* .

2. *The marginal distribution of* **s** *is**where* *and* .

*Proof*. By Bayes' theoremThe product on the right-hand side is of the form , whereIn the last step we completed the square with respect to and used the fact that **T** is symmetric. Up to a constant that does not depend on we hence getwhere . This proves the first part of the theorem and—by linear superposition—the validity of Equation 11.

To prove the second part of the theorem, observe that with and with . Putting these equalities together, we getThis, being a linear combination of independent multivariate normal variables, is still multivariate normal with mean and its covariance matrix is given byThis proves the second part of the theorem as well as formula (15).▪

## APPENDIX B: NONLINEAR TOY MODELS

In this section we describe a class of toy models that are nonlinear in the parameter θ ∈ ℝ and have nonnormal, possibly heteroscedastic error terms. Still their likelihoods are easy to calculate analytically. We setHere *f _{i}*(θ) are monotonically increasing continuous functions of θ and ϵ

_{i}(θ) are independent, uniformly distributed error terms in the interval [–

*u*(θ),

_{i}*u*(θ)] ⊆ ℝ, where

_{i}*u*(θ) are nondecreasing, continuous functions:

_{i}It is straightforward to check that for a prior π(θ) the posterior distribution of θ given is (up to a normalizing constant)whereandFor the simulations in Table 2 we chose *n* = 5, , and . The prior was .

## Acknowledgments

We are grateful to Laurent Excoffier, David J. Balding, Christian P. Robert, and the anonymous referees for their useful comments on a first draft of this manuscript. This work has been supported by grant no. 3100A0-112072 from the Swiss National Foundation to Laurent Excoffier.

## Footnotes

↵1 These authors contributed equally to this work.

Communicating editor: J. Wakeley

- Received August 27, 2009.
- Accepted September 2, 2009.

- Copyright © 2010 by the Genetics Society of America