Robust Approximate Bayesian Inference with Synthetic Likelihood

Bayesian synthetic likelihood (BSL) is now a well-established method for conducting approximate Bayesian inference in complex models where exact Bayesian approaches are either infeasible, or computationally demanding, due to the intractability of likelihood function. Similar to other approximate Bayesian methods, such as the method of approximate Bayesian computation, implicit in the application of BSL is the maintained assumption that the data generating process (DGP) can generate simulated summary statistics that capture the behaviour of the observed summary statistics. This notion of model compatibility with the observed summaries is critical for the performance of BSL and its variants. We demonstrate through several examples that if the assumed DGP differs from the true DGP, model compatibility may no longer be satisfied and BSL can give unreliable inferences. To circumvent the issue of incompatibility between the observed and simulated summary statistics, we propose two robust versions of BSL that can deliver reliable performance regardless of whether or not the assumed DGP can generate simulated summary statistics that mimic the behavior of the observed summaries. Simulation results and two empirical examples demonstrate the performance of this robust approach to BSL, and its superiority over standard BSL when model compatibility is not in evidence.


Introduction
In situations where the likelihood of the underlying model is intractable, approximate Bayesian methods are often the only feasible solution to conduct Bayesian inference. Indeed, approximate Bayesian methods are an increasingly common tool in the arsenal of the practicing statistician and allow users to conduct reliable inference in models where exact Bayesian inference procedures are either infeasible, or too computationally demanding.
We demonstrate through several examples that when the model is misspecified BSL can deliver unreliable results. To circumvent this issue, we propose two novel versions of BSL that can deliver reliable results regardless of whether the model is correctly specified. These robust versions of BSL are based on adjusting the simulated summaries to ensure a form of compatibility between the assumed model and the chosen summary statistics.
The first approach is we propose is based on adjusting the simulated mean of the summaries, while the second approach adjusts the variance. We demonstrate both approaches through a series of Monte Carlo experiments and empirical examples. The results of these experiments yields a striking juxtaposition between our new robust approach and standard BSL: when the model is misspecified, standard BSL inference can fail completely, due to extremely low acceptance rates within the MCMC sampling algorithm, while the robust versions of BSL proposed herein remain well-behaved even under large levels of model misspecification. The improved computational efficiency of our approaches allow practitioners to explore sources of model misspecification to inform further model development when using BSL as the inferential tool.
The remainder of the paper is organized as follows. In Section two we give a brief overview of BSL and model misspecification in the context of BSL. Section three presents our robust approach to BSL. Section four contains a mix of Monte Carlo and empirical examples that demonstrate the performance of this robust BSL approach, and the poor performance of standard BSL. Section five concludes.

Bayesian Synthetic Likelihood Framework
We observe the sequence of data y = (y 1 , . . . , y n ) ⊤ , n ≥ 1, and denote by P n 0 the true distribution of the sample. The true distribution is unknown and instead we consider a class of probability measures {θ ∈ Θ ⊂ R d θ , n ≥ 1 : P n θ }, that we assume, for some value of θ, have generated the data and which admit the corresponding conditional density p n θ (·). Given prior beliefs over the unknown parameters in the model θ, represented by the probability measure Π(θ) and its density π(θ), our aim is to produce draws from the exact posterior density π(θ | y) ∝ p n θ (y)π(θ).
However, in certain examples, the density p n θ (y) is intractable, and sampling from π(θ | y) can be computationally costly or infeasible.
In situations where the likelihood is intractable, so-called likelihood-free methods can be used to conduct inference on the unknown parameters θ. The most common implementations of these methods are approximate Bayesian computation (ABC) and Bayesian synthetic likelihood (BSL). Both ABC and BSL generally degrade the data down to a vector of summary statistics and then perform posterior inference on the unknown θ, now conditional only on this vector of summary statistics, by selecting values of θ that yield a close match between the summary statistics calculated from the observed data and those calculated under the model.
However, given the complexity of the assumed model, P n θ , it is unlikely that the structure of G n (·|θ) is any more tractable than the original likelihood function. Therefore, simulation-based sampling schemes must be applied to generate samples from π[θ | η(y)].
The approximate methods of ABC and BSL differ in how g n (·|θ) is estimated. ABC forms an implicit nonparametric estimator of g n (·|θ), while BSL uses a parametric or semi-parametric  approximation to g n (·|θ). Given the nonparametric nature of ABC, it is not surprising that the curse of dimensionality inherent in ABC, with regards to the dimension of the summary statistics (Blum, 2010), as well as in the dimension of the parameters themselves (Frazier et al., 2018), ensures that ABC based approaches are heavily restricted in terms of the number and type of summaries they can choose. Due to this need to use only a small number of summary statistics within the analysis, standard-ABC-based inference approaches can often lead to imprecise inference.
In contrast to ABC, BSL employs a parametric, or semi-parametric, approximation to G n (·|θ). In particular, BSL replaces -an estimate of -the (intractable) density g n (·|θ) by a multivariate Gaussian approximation: where µ m (θ) and Σ m (θ) are simulated estimates of the mean and variance of the summaries (conditional on θ) and where each simulated data set z i , i = 1, . . . , m, are generated iid from P n θ . The Gaussian approximationḡ n [η(y) | θ] is then directly used within an MCMC sampling scheme to sample from the following approximation to the partial posterior, hereafter referred to as the BSL posterior,π Due to the parametric nature of (1), BSL can often treat summary statistics of larger dimension than ABC and can lead to sharper inference in some cases. While the validity of the Gaussian approximation is often warranted if the underlying summaries satisfy a central limit theorem (Wood, 2010), even in cases where the summary statistics are far from Gaussian, BSL has displayed some insensitivity to violations of this assumption (Price et al., 2018). However, if the statistics are very far from being Gaussian, this can result in less efficient sampling from (1), see An et al. (2018) for a demonstration.

Model Incompatibility and its Consequences
Approximate Bayesian inference methods based on summary statistics implicitly assume that the model being used to generate the simulated summary statistics η(z) can replicate the behavior of the observed summary statistics η(y). That is, BSL is not required to match every aspect of the data, but, rather, only those features of the data that are captured via the summary statistics η(y). This differs from a standard Bayesian framework based on a likelihood, where, under general regularity conditions, the posterior ultimately gives higher probability mass to values of θ ∈ Θ that ensure the Kullback-Liebler (KL) divergence D (P n 0 P n θ ) = log is as close to zero as possible. When the model is misspecified, i.e., when P n 0 = P n θ for any θ ∈ Θ, following Kleijn and Van der Vaart (2012), the posterior eventually places increasing mass on the value θ * ∈ Θ that minimizes the KL-divergence.
In contrast, BSL is not based on trying to find the P n θ that is closest to P n 0 in KL divergence but on trying to match simulated and observed summary statistics. Therefore, the meaningful concept of model misspecification in BSL is not KL divergence but that the choice of the assumed model, allied with our specific choice of summary statistics, are able to replicate the behavior of the observed summary statistic. Using the framework of Marin et al. (2014), and following Frazier et al. (2017), we formalize this notion of misspecification as follows.
Definition 1. We say that the model P n θ × Π and summary statistic map Intuitively, compatibility means that η(y) must be in the support of µ m (θ), at least asymptotically. Compatibility implies that, for some value of θ, η(z) can replicate the behavior of η(y) when z is simulated under P n θ . By focusing inference on summary statistics, approximate inference methods have a type of in-built robustness against model misspecification: the goal of these methods is not to find the region of the parameter space Θ that allow us to match all the features of the DGP; instead, approximate methods only seek to find regions of Θ that allow us to match the features of the data measured by η(·). Therefore, we argue that the researcher should not be concerned about the correct specification of the entire DGP but should instead concern herself with the satisfaction of the above compatibility condition. If this condition is satisfied, even though the model P n θ may be misspecified, we can still, in some sense, adequately capture the features of the data that are described by η(y), and which are likely to be of primary interest to the analysis.
However, as recently discussed by Frazier et al. (2017), in the context of ABC, when approximate methods are based on a model and summary statistic combination that are not compatible, the resulting inference can be suspect. Given that ABC and BSL are based on very similar principles, it is highly likely that BSL will suffer from the same issues as ABC when the above compatibility condition is not satisfied. While we demonstrate this with several examples, both toy and empirical in Section 4, we first consider an artificially simple example where BSL should perform very well, but, due to this incompatibility issue, BSL can deliver unreliable inference.

Consequences of Incompatibility: Toy Normal Model
To demonstrate the performance of BSL under model/summary incompatibility, we consider an artificially simple example: our goal is to conduct inference on the mean parameter θ in the model Assume that we believe with certainty that σ = 1, and so the model we will use to generate data within BSL is given by However, we will expressly consider the scenario where the assumption σ = 1 is false and the true DGP for the observed data maintains σ = 1, which implies that the assumed DGP for z differs from the actual DGP for y. We consider as summary statistics the sample mean η 1 (y) = 1 n n i=1 y i and the sample variance η 2 (y) = 1 n−1 n i=1 (y i − η 1 (y)) 2 . When σ = 1, the model is not compatible with the variance summary, and, hence, is misspecified.
While this example is artificially simple, it will be useful to demonstrate the potential issues with BSL when the model is misspecified. Indeed, regardless of the misspecification, one might think that BSL should perform well: the first summary statistic, η 1 (y), is Gaussian, the statistic η 2 (y) satisfies the central limit theorem, and η 1 (y) is sufficient for θ.
For this experiment, we generate simulated data sets for y such that each data set corresponds to a different value of σ, with σ taking values from σ = 0.2 to σ = 2 with evenly spaced increments of 0.1. Across all the data sets we fix the random numbers used to generate the "observed" data and only change the value of σ to isolate the impact of model misspecification; i.e., for each of the data sets, we use the same random numbers. The sample size across the experiments is taken to be n = 50. Our prior for θ is θ ∼ N (0, 10).
We generate m = 10, 000 simulated data sets to estimate µ m (θ) and Σ m (θ), and carry out BSL using a random-walk Metropolis algorithm. The sampler is initialized at the true value of the parameters, θ = 1, and run for 10,000 iterations. Since the chain is started at the truth, we retain all values for inference.
First, we analyze how the BSL posterior mean responds to varying levels of model misspecification. Panel A of Figure 1 plots the posterior means and corresponding 95% credible intervals for each data set and Panel B gives the acceptance rates across the various data sets. From the results in Panel A, it seems that BSL displays resilience to model misspecification, except perhaps for some instability at large levels of misspecification. However, analyzing the results in Panel B, we see that the stability implied by the results in Panel A is not genuine. In particular, the results in Panel B demonstrate that at high levels of misspecification, σ either very small or very large, BSL has a very difficult time exploring the parameter space. As a consequence, the MCMC chain stays stuck for long periods of time and the acceptance rates in BSL plummet from a peak of about 70%, to a low of 0.85% at σ = 2 despite the large number of model simulations used to estimate the synthetic likelihood. At higher levels of misspecification, the MCMC chain rarely moves and so most values simply reflect the starting value of θ = 1, which gives the results in Panel A of Figure 1 the false sense that BSL is performing well across all levels of misspecification.
In this scenario, the poor performance of BSL can be directly traced to the incompatibility of the second summary statistic η 2 (y): no matter the value of θ simulated, we can never hope to match the value of η 2 (y). This gap between η 2 (y) and η 2 (z) causes the MCMC chain to stick, and results in unreliable inference. Moreover, with real data it is not clear a priori how one could hope to reliably detect such incompatibility issues, since the true DGP is unknown. To circumvent this poor performance of BSL in contexts where the summaries are not compatible, in the following section, we propose a robust version of BSL that delivers reliable inference when the summaries are not compatible and which is capable of detecting precisely which summaries, if any, are incompatible.

Robust Bayesian Synthetic Likelihood
We propose two possible strategies for conducting inference using BSL when the model and summaries are incompatible. The first approach augments the mean of the simulated summaries with additional free parameters, while the second approach augments the variance of the simulated summaries with additional free parameters. Both specifications will allow us to conduct robust inference on the model parameters and determine which of the summaries may be incompatible with the data. Throughout the remainder, we refer to these approaches as "mean" and "variance" robust BSL (R-BSL).

Mean Adjustment and Prior Specification
Incompatibility implies that the observed statistics η(y) are not in the support of the simulated mean µ m (θ), for any θ ∈ Θ, with high probability. Therefore, one approach to create a BSL procedure that is robust to this incompatibility issue is to adjust the vector of simulated means µ m (θ) by adding to µ m (θ) an additional free parameter Γ so that η(y) is always in the support of this new simulated mean. Defining the joint vector of unknown parameters as ζ := (θ ⊤ , Γ ⊤ ) ⊤ ∈ Θ × G ⊂ R d θ × R dη , we define the vector of means for use in BSL as Note that, by considering the scaled adjustment term diag(Σ 1/2 m (θ))Γ, we ensure that these components are measured in the same units as µ m (θ), which allows us to treat Γ as if they were unitless.
Given this linear adjustment, and under weak conditions on the summary statistics and the parameter space Θ × G, it is simple to see that φ m (ζ) will be compatible with η(y) for any prior choice on Γ such that each individual component has support over R.
Denoting our prior on ζ by π(ζ), following Price et al. (2018), the augmented BSL target, which we refer to as the Robust BSL-mean (R-BSL-M) posterior, is to generate samples from (2)

Prior Choice: Laplace Prior
To ensure that the observed summary η(y) is in the support of φ m (ζ), our prior on the components of Γ should allow for diag[Σ 1/2 m (θ)]Γ to escape the support of µ m (θ) with non-negligible probability. However, given that some components of the original µ m (θ) are likely compatible with some components of η(y), we want to make sure that the adjustment term diag[Σ 1/2 m (θ)]Γ does not unduly perturb those components that are compatible. Therefore, we should choose a prior that places the vast majority of its mass near the origin. In this way, our prior choice for Γ should induce sparsity in the components of Γ: only the components of Γ that correspond to incompatible summaries should receive significant posterior probability away from zero, while the components of Γ that correspond to compatible summaries should have the majority of their posterior mass near the origin.
With these dual requirements in mind, and given that each component of Γ has the same prior scale, we propose to follow the Bayesian Lasso literature (Park and Casella, 2008) and use independent Laplace (i.e., double-exponential) priors for each component of Γ, with fixed location 0 and common scale λ > 0: We will often denote this prior by La(0, λ). The Laplace prior for Γ guarantees that the majority of prior mass for γ j is near the origin, but has thick enough tails so that φ m (ζ) is compatible with virtually any η(y) that would be used in practice.
The hyper-parameter λ should be chosen so that the prior support of π(Γ) compliments the support of θ → µ m (θ). That is, λ should be chosen so that the prior support of Γ is larger than the support of µ m (θ), which will allow us to detect deviations from compatibility, but not so large as to cause the statistic φ m (ζ) to have very heavy tails. Indeed, if the tails of φ m (ζ) are too heavy, the implicit normality assumption made in BSL will be strongly violated and can result in inefficient sampling.
Since there is no reason to believe a priori that θ and Γ are related, we take as our overall prior on ζ in R-BSL-M to be π(ζ) := π(θ)La(0, λ).

Variance Compatibility and Prior Specification
Compatibility requires that η(y) be in the support of µ m (θ), for some θ ∈ Θ. While one approach to ensure this compatibility is to adjust the mean of the simulated summaries, an alternative is to inflate the variance of the simulated summaries to ensure that η(y) is always in the support of µ m (θ). Under regularity conditions and for fixed m, it is likely to be the case that the centered , with the variance of Z n,m (θ) decreasing like 1/(mn). Consequently, for n (or m) large enough, if for a given value of θ the statistic {η(y) − b(θ)} is more than a few standard deviations (as measured by Σ 1/2 m (θ)) away from {µ m (θ) − b(θ)}, we can effectively view the summaries as being incompatible. In such cases, the BSL posterior is likely to have low acceptance rates in this region of the parameter space.
Given this characterization, we see that an approach to ensure that µ m (θ) is compatible with η(y) is to artificially inflate the variance Σ m (θ) so that the variance of Z n,m (θ) never completely collapses to zero, and thus {η(y) − b(θ)} can always be found in the support of Z n,m (θ), albeit perhaps with small probability. More specifically, we propose to artificially inflate the variance used within BSL by adding to Σ m (θ) the free parameters Γ = (γ 1 , . . . , γ dη ) ⊤ .
Recalling ζ = (θ ⊤ , Γ ⊤ ) ⊤ , a robust BSL procedure based on adjusting the variance can be implemented by re-defining the variance of the simulated statistics used within BSL to be where [Σ m (θ)] ii denotes the (i, i) element of Σ m (θ). Given the structure of V m (ζ), γ i can be interpreted as an inflation factor operating on the standard deviations of the original BSL variance. An equivalent interpretation is that the i-th BSL variance is multiplied by the factor 1 + γ 2 i . Using V m (ζ) in place of Σ m (θ) within the standard BSL posterior target, (1), and for π(ζ) an appropriate prior on ζ, the Robust BSL-variance (R-BSL-V) posterior is given as:

Prior Choice: Exponential Prior
Note that, by considering the standardization in (4), we ensure that each γ i has the same scale and can be considered as unitless. Moreover, similar to the case of the mean adjustment BSL approach, there is no reason to believe there is any a priori correlation between θ and Γ, so we can consider independent priors, i.e, π(ζ) := π(θ)π(Γ). While several prior choices exist for Γ, following the arguments for the prior choice in the mean adjustment procedure, we need to choose a prior for the components of Γ that have large probability mass near the origin, but that also have enough mass out in the tails to ensure we can detect incompatible summaries in a robust manner.
To this end, we consider independent exponential priors for each component γ i , (i = 1, . . . , d η ), with common scale λ > 0: The hyper-parameter λ should be chosen so that a large amount of prior mass is close to the origin, so as not to over-inflate the variance of the simulated summaries that are compatible. While this choice of prior is not, strictly speaking, a shrinkage prior, it is still the case that we should observe some shrinkage like behavior for summaries that are compatible. That is, for the summaries that are compatible, this additional inflation by Γ is unnecessary and we expect that, for appropriate choices of λ, the addition of this component will not greatly affect the corresponding components in the variance. In contrast, for the summaries that are not compatible, this adjustment term is critical to ensure that the variance of the summaries is large enough to contain the observed summary η(y).
Our robust BSL methods operate on an extended state space over θ and Γ with target distribution To sample these target distributions we use a component-wise MCMC algorithm that updates, in turn, θ conditional on Γ and then Γ conditional on θ. The update for θ is the same as in standard BSL, but where the adjusted mean or inflated variance is computed as appropriate using the current value of Γ. As before, the update for θ involves generating m model simulations, η(z i ), i = 1, . . . , m.
The update for Γ holds the currently accepted model simulations fixed, and thus µ m (θ) and Σ m (θ) are fixed within the update step for Γ. Each component of Γ, γ j , for j = 1, . . . , d η , is updated separately, conditional on the current values of the remaining components (denoted γ /j ). The full conditional distribution for γ j is given by where γ * j is a realisation of γ j and ζ * = [θ, γ * j , γ /j ]. We sample this full conditional distribution using a slice sampler, and, in particular, the "stepping out" and "shrinkage" procedures detailed in Neal (2003). The appeal of the slice sampler is that the acceptance probability is one, and thus there are no tuning parameters that can affect the statistical efficiency. However, a stepping out width needs to be selected, which can impact the speed of the slice sampler. Given that the components of Γ are eventually scaled by the summary statistic standard deviation, and given that our prior choices effectively penalise large values, we expect each component of Γ to be O(1). Thus, we set the stepping out width to be 1, except for the lower bound of γ j in the variance inflation, which is immediately set to 0; and hence the stepping out procedure is not required. We find this choice of width to be suitable, and since updating Γ does not require any model simulations, the slice sampler is very fast. Hence, importantly, our robust BSL methods do not require any additional tuning and the run-time per iteration for non-trivial applications is not noticeably slower. The full MCMC algorithm to sample from the R-BSL posteriors is provided in Algorihtm 1.
It is important to note that under either R-BSL approach, we recover the original BSL target when γ 1 = · · · = γ dη = 0.
Therefore, if the model P n θ can generate summaries η(z) that match η(y), under an appropriate prior specification, the posteriorπ[Γ | η(y)] should not differ substantially from the prior, with most of the posterior mass located near the origin. Indeed, if one wished, this relationship between prior and posterior could be tested to detect precisely which of the components are incompatible.
While obviously related, the two R-BSL adjustment procedures operate on very different principles. The R-BSL-M approach is based on finding joint values (θ, Γ) under which φ m (ζ) − η(y) is "small", which we can interpret as, for fixed θ, finding values of Γ so that µ m (θ) + diag Σ 1/2 (θ) Γ matches η(y). In contrast, the R-BSL-V procedure is not required to find values of θ such that µ m (θ) − η(y) is small. Rather, for a fixed value of θ, the R-BSL-V approach seeks a value of Γ so that the weighted norm is "small". Therefore, if µ m (θ) − η(y) is already small, no variance inflation is needed, and the resulting value of Γ should be close to zero; if instead µ m (θ) − η(y) is "large", the term µ m (θ) − η(y) Vm(ζ) can always be made small simply by choosing a large value of Γ.
1 Estimate µ m (θ 0 ) and Σ m (θ 0 ) via m independent model simulations at θ 0 %%% Obtain Γ i via the following, which does not require any model simulations Update γ i j with a slice sampler with target density Estimate µ m (θ * ) and Σ m (θ * ) via m independent model simulations at θ *

Returning to the Normal Example
In this section, we compare the performance of BSL and R-BSL under model misspecification in the toy normal example. Recall that our goal is inference on θ in the model where we explicitly assume that σ = 1, and generate data in BSL according to We consider as summary statistics the sample mean η 1 (y) = 1 n n i=1 y i and the sample variance η 2 (y) = 1 n−1 n i=1 (y i − η 1 (y)) 2 . The sample size across all experiments is taken to be n = 50, and our prior is θ ∼ N (0, 10). Recall that, under this setting, we observed that the BSL could lead to unreliable inferences when σ 2 was much larger or smaller than unity.
For the experiments in this section, we use precisely the same "observed" data sets in Section 2.3, but also apply the two different versions of the R-BSL approach, based on the mean and variance adjustments, to the data sets. Similar to BSL in the previous experiment, for both R-BSL procedures we use a random-walk Metropolis algorithm that is initialized at the true value of the parameters, θ = 1, and that we run for 10,000 iterations. Since the chain is started at the truth, we retain all values for inference.
For the mean adjustment approach, our prior on Γ = (γ 1 , γ 2 ) ⊤ is independent La(0, λ) with λ = 0.5 across both parameters, which places the vast majority of prior mass between -3 and 3. For the variance adjustment, we give each parameter independent Exp(λ) priors and choose λ so that the prior mean is 0.3. Figure 2 presents the results and mirrors the output in Figure 1 but for the robust versions of BSL.
The results in Panels A and B of Figure 2 demonstrate that the R-BSL approaches displays resilience to model misspecification. Moreover, the results in Panel C demonstrates that the acceptance rate issue that plagued the standard BSL approach is not in evidence for the R-BSL approaches. Indeed, the acceptance rate for the variance adjustment R-BSL approach is nearly unaffected by the level of model misspecification, while the mean adjustment version does display some degradation for σ > 1.5, but still maintains acceptance rates above 5% in all cases. In comparison, when σ = 2, the standard BSL acceptance rate was 0.85%. Figure 3 displays the resulting posterior densities for Γ across the two adjustment procedures, and across all levels of misspecification. Panels A and B give the results for the mean adjustment approach, and correspond to the components γ 1 and γ 2 , respectively, while panels C and D give the results for the variance adjustment approach. For comparison purposes, the black figure in each panel represents the prior densities, and the color-coding in each figure represents the level of misspecification, and where the red colored density encodes correct model specification (i.e., σ = 1).
Focusing on Panel A, we see that the posterior densities for the γ 1 component, which captures our ability to match the first observed statistic (the mean), are indistinguishable from the prior across all the observed data sets, which implies that we can match the mean of this model regardless of model misspecification. In contrast, in Panel B we see that the second component, which captures our ability to capture the second observed statistics (the variance), looks nothing like the prior, except perhaps at low levels of misspecification.  Figure 2: Panels A and B give the posterior means and 95% credible sets for the two different R-BSL procedures. Panel C gives the acceptance rates of the different R-BSL procedures and compares them with the standard BSL approach across the different data sets.
Panels C and D describe precisely the same story as in Panels A and B but correspond to γ 1 and γ 2 in the variance adjustment approach. That is, in the case of the variance adjustment BSL approach (R-BSL-V), under correct specification (σ = 1) the posteriors are indistinguishable from the priors, we are easily able to detect departures from compatibility for the second summary statistics, and the posteriors for the first adjustment term remain indistinguishable from the prior across all data sets. These striking results demonstrate that both R-BSL approaches are capable of producing meaningful inference on the unknown parameter and are capable of determining which features of the model we are not above to replicate.

Moving Average Model
To further demonstrate this methodology, we turn to a common toy examples used to demonstrate approximate inference methodology, the moving average model of order two MA(2). While the likelihood for this model is tractable, and many exact Bayesian routines can be used to conduct inference on the model parameters, it has nonetheless gained widespread use in approximate inference as a toy example to demonstrate approximate inference methodology (see, e.g., Marin et al., 2012). Assume the researcher believes y is generated according to an MA(q) model: where, say, e t ∼ N (0, 1) i.i.d and θ 1 , ..., θ q are such that the roots of the polynomial θ i x i all lie outside the unit circle. Specializing this model to the case where q = 2, we have that z t = e t + θ 1 e t−1 + θ 2 e t−2 , and the unknown parameters θ = (θ 1 , θ 2 ) ⊤ are assumed to obey − 2 < θ 1 < 2, θ 1 + θ 2 > −1, θ 1 − θ 2 < 1.
Since we are interested in examining the ability of BSL and the robust variants proposed herein to deal with model incompatibility, we consider that, while the researcher believes the data is generated according to an MA(2) model, equation (7), the actual DGP for y evolves according to the stochastic volatility (SV) model 0 < ρ < 1, 0 < σ v < 1, u t and v t and both iid standard Gaussian. In this case, if one takes η (y) := (η 0 (y), η 1 (y), η 2 (y)) ⊤ it follows that, under the DGP in (9), For any value of ω, σ v and ρ such that it follows that (P n θ × Π, η) is not compatible. From the definition of b(θ) and b 0 , it also follows that the value that minimizes b(θ) − b 0 Σ(θ) is θ * = (0, 0) ⊤ , and is the value we would expect the BSL posterior to concentrate onto in the limit.
To understand how BSL and the R-BSL alternatives proposed herein perform under this misspecified model, we consider the following Monte Carlo analysis: we generate n = 100 observations from the SV model in (9) and use BSL, R-BSL-M, and R-BSL-V to conduct inference on the unknown parameters θ in the misspecified MA(2) model. Each BSL approach uses m = 10 simulated data sets to estimate the mean and the variance.
Given the form of incompatibility in (10), we would expect that the R-BSL-M and R-BSL-V procedures would detect incompatibility in the first summary statistic, the sample variance, while the corresponding adjustment components for the other summaries would not be indistinguishable from the prior. This analysis follows from noting that, under the true DGP in (9), the first two autocorrelations are zero for all values of (ω, ρ, σ v ) ⊤ , and hence the MA(2) parameters should concentrate posterior mass around zero.
For the R-BSL-M and R-BSL-V procedures, we consider the same priors specification as for the simple normal example, and the exact same choice of hyper-parameters: namely, for R-BSL-M we consider that γ i ∼ La(0, λ = 0.5), and for the R-BSL-V approach, we take as our prior γ i ∼ Exp(λ) and set λ so that the prior mean is 0.3. For each procedure we consider starting values obtained from the maximum likelihood estimators of the MA(2) model, and the values of θ are obtained via a random-walk Metropolis sampler with fixed covariance matrix given by 0.1 0 0 0.1 .
We run the MCMC sampler for 50,000 iterations and discard the first 10,000 for burn-in. Under this Monte Carlo design, the acceptance rates for R-BSL-V and R-BSL-M are 11% and 6.5%, respectively, while the acceptance rate for the standard BSL approach is just 0.58%. Indeed, the resulting posterior samples from standard BSL are unusable and deviate sharply from those obtained by the R-BSL approaches. The posterior densities for θ = (θ 1 , θ 2 ) ⊤ across the R-BSL methods are displayed in Figure 4. 1 Both R-BSL-M and R-BSL-V display significant posterior concentration around 0. This is significant as the value that asymptotically minimizes the weighted norm µ m (θ) − η(y) Σm(θ) , and, thus, which we would expect BSL to concentrate onto, is θ * = (0, 0) ⊤ . To understand the poor performance of standard BSL in this example, we can examine the marginal posteriors for the different adjustment terms across R-BSL-V and R-BSL-M, which are given in Figure 5. The results demonstrate that BSL can not reliably match the first summary statistic, the sample variance, while the posteriors that correspond to the first and second-order autocorrelations do not differ from their priors (respectively given as Laplace and exponential) and signify that we can indeed match the second and third summaries. This finding is perfectly in line with the above suggested theoretical behavior: to minimize µ m (θ) − η(y) BSL must place high posterior mass around θ = (0, 0) ⊤ , but doing so means that we cannot adequately capture the sample variance. While the adjustment procedures can account for this issue, the BSL procedure can not and hence yields low acceptance rates and inaccurate posterior inference. 2

Background
Collective cell spreading models are often used to gain insight into the biological mechanisms governing, for example, wound healing and skin cancer growth (e.g. Vo et al. (2015b,a)). Browning et al. (2018) develop a simulation-based model where cells are able to move freely in continuous space. They calibrate the model to real in vitro data collected from a cell proliferation assay experiment using a rejection-based ABC algorithm. Here, using our new synthetic likelihood methods, we demonstrate that the model is not compatible with the observed summary statistic and provide insight into what aspects of the data that the model is not able to recover.
The model of Browning et al. (2018) is a stochastic individual based model where cells move and interact in a two dimensional space. Here we provide only brief details of the model and refer to Browning et al. (2018) for the full description. Proliferation (cell birth) and motility (movement) for each cell evolves in continuous time according to a Poisson process. The intrinsic rates are given by p and m for proliferation and motility events, respectively. The rates of these processes are also neighbourhood-dependent, with rates decreasing as the amount of crowding around a cell increases. The closeness of cells is governed by a Gaussian kernel that depends on a fixed cell diameter, σ. When a cell proliferates, it places a new cell randomly in its neighbourhood according to an uncorrelated two dimensional Gaussian centered at the cell location with component variances of σ 2 . When a motility events occurs, the cell moves a distance of σ. The direction of the move depends on cell density, biased towards lower cell density. A parameter used to help determine the move direction, γ b , is part of a Gaussian kernel used to measure the closeness of cells. The parameter of interest is θ = (p, m, γ b ) ⊤ .
In the experiments of Browning et al. (2018), images of the cell population are taken every 12 hours starting at 0 hours with the final image taken at 36 hours. Browning et al. (2018) use the number of cells and the pair correlation computed from each image as the summary statistics, resulting in a six dimensional summary statistic here. The pair correlation is the ratio of the number of pairs of agents separated by some pre-specified distance to an expected number of cells separated by the same distance if the cells were uniformly distributed in space.
The prior distribution is set as p ∼ U(0, 10), m ∼ U(0, 0.2) and γ b ∼ U(0, 20) with no dependence amongst parameters, as in Browning et al. (2018). We use MCMC to sample the posterior with 50,000 iterations and no burn-in as we initialise the chain at the point estimate θ = (1, 0.04, 6) ⊤ reported in Browning et al. (2018). We sample over the space of a logit-type transformation of θ so that any proposal is within the prior bounds. We use a multivariate Gaussian random walk proposal on the transformed space with a covariance matrix obtained via some pilot MCMC runs. We use m = 50 model simulations to estimate the synthetic likelihood at each MCMC iteration. We run our methods on both simulated (using the point estimate of Browning et al. (2018)) and real data.

Results
For the priors for each component of γ, we use a Laplace distribution with a scale of 0.5 for the mean adjustment method and an exponential prior with mean of 0.5 for the variance inflation method.
Firstly we present results for the simulated data (model correctly specified). As shown in Figures 6 and 7, the posterior distributions on θ are similar regardless of whether BSL or R-BSL is applied. The MCMC acceptance rates for BSL and R-BSL-M are both 21% and 20%, respectively. The variance inflation seems to allow for a slightly increased acceptance rate (24%) compared to mean adjustment.
The posterior distribution for each component of γ is shown in Figure 8. It can be seen that most posteriors are similar to the prior. For both R-BSL methods, there is no indication that any of the statistics are incompatible with the model, as expected.
For the real data, the MCMC acceptance rate using BSL is only 3% as the variance of the synthetic likelihood is high generating long periods of no acceptance. Applying R-BSL-M and R-BSL-V results in an MCMC acceptance rate of roughly 12% and 18%, respectively, permitting statistical inference. Again, the variance inflation seems to produce an improved acceptance rate. The univariate posterior distributions for γ for R-BSL-M and R-BSL-V are shown in Figures  9 and 10, respectively. Both methods identify that the model is not compatible with the pair correlation statistic at 12 and 36 hours.
In Figure 11, we show R-BSL-M posterior predictive distributions of the summary statistics without (left) and with (right) using estimated mean adjustment parameters with the observed summaries overlaid. The corresponding plots for the variance inflation is shown in Figure 12. From both figures, it is evident from the plots on the left that the model is successful in tracking the number of cells over time. However, the model underestimates the rate of decrease in the pair correlation over time. This is valuable information that might enable mathematical biologists to extend the model so that this data feature can be better captured. Both adjustment methods have allowed us to make this inference. It can be seen in the second column of Figure 11 that the mean adjustment is able to shift the predictions so that observed statistic does not lie so far in the tails. From the second column of Figure 12, the variance adjustment expands the predictions so that the observed statistic does not lie so far in the tails. Finally, Figure 13 compares the posterior distributions for the mean adjustment and variance inflation, together the standard synthetic likelihood results. It can be seen that the posterior distributions are broadly similar, except that the standard synthetic likelihood results suffer from substantial Monte Carlo error due to the small acceptance rate.

Background
We consider an individual-based model of a species called Fowler's Toads (Anaxyrus fowleri ) developed by Marchand et al. (2017), which was also analysed by An et al. (2018). Here we give very brief details, with more information in Marchand et al. (2017) and An et al. (2018). The model assumes that a toad hides in its refuge site in the daytime and moves to a randomly chosen foraging place at night. GPS location data are collected on n t toads for n d days, i.e. the observation matrix Y is of dimension n d ×n t (n t = 66 and n d = 63 here). Then Y is summarised to 4 sets comprising the relative moving distances for time lags of 1, 2, 4, 8 days. For instance, y 1 consists of the displacement information of lag 1 day, Simulating from the model involves two distinct processes. For each toad, we first generate an overnight displacement, ∆y, then mimic the returning behaviour with a simplified model. The overnight displacement is assumed to belong to the Lévy-alpha stable distribution family, with stability parameter α and scale parameter δ. The total returning probability is a constant p 0 , if a return occurs on day i, 1 ≤ i ≤ m, then the return site is the same as the refuge site on day i, where i is selected randomly from 1, 2, . . . , m with equal probability. Here we consider both simulated and real datasets. For the synthetically generated data we take θ = (α, δ, p 0 ) ⊤ = (1.8, 45, 0.6) ⊤ , which is informed by the parameter estimates obtained in Marchand et al. (2017). We use a uniform prior over (1, 2) × (0, 100) × (0, 0.9). Marchand et al. (2017) consider three variations on the model. Here, we consider their 'Model 2' since there is a strong indication from their results that this model is not able to recover some of the chosen summary statistics.
As in Marchand et al. (2017), the dataset of displacements is split into two components. If the absolute value of the displacement is less than 10 metres, it is assumed the toad has returned to its starting location. For the summary statistic, we consider the number of toads that returned (Marchand et al. (2017)). For the non-returns (absolute displacement greater than 10 metres) we consider a larger collection of summaries. We calculate the log difference between adjacent pquantiles with p = 0, 0.1, . . . , 1 and also the median. These statistics are computed separately for the four time lags. This results in 48 statistics in total, which is hard to handle for conventional ABC methods. An et al. (2018) demonstrate that BSL is computationally efficient enough to analyse simulated data for this application with a similar number of summary statistics.
For the γ parameters of the mean adjustment and variation inflation procedures we use the same priors as the previous example.

Results
We first consider the simulated dataset, where we use n = 300 simulations to estimate the synthetic likelihood at each MCMC iteration. Standard BSL, together with the two incompatability extensions, produce approximate posteriors shown in Figure 14. As can be seen, the adjustments produce posteriors remarkably similar to BSL with slightly inflated variances. The MCMC acceptance rates for BSL, R-BSL-M and R-BSL-V are 11%, 9% and 22%, respectively. As consistent with previous results, the variance adjustment improves the computational efficiency, even when the model is correctly specified. The posterior distributions for γ of R-BSL are shown in Figure 15. In all cases, the posteriors are not too dissimilar to the prior.
For the real data, we required n = 2000 simulations for estimating the synthetic likelihood to obtain an acceptance rate of 9% for BSL. However, the chain still suffered from periods of stickiness. In contrast, with only n = 500, the R-BSL-M and R-BSL-V produce acceptance rates of 7% and 15%, respectively, without substantial stickiness. The variance inflation method offers a computational improvement of about one order of magnitude over BSL when accounting for both acceptance rate and number of simulations. The posterior distributions for the components of γ for the R-BSL methods are shown in Figure 16. It is evident from the plots that our methods have identified that there are three or four statistics that the model is not compatible with. The statistic with the largest incompatibility is the number of returns for lag 1. For R-BSL-V, the 95% posterior predictive interval for this statistic is (262, 346) with an observed value of 234. Other statistics showing some incompatibility are the first quantile differences of the non-returns for lags 3 and 4. Figure 17 confirms that the observed data are not consistent with the posterior predictive distribution of the (log) non-return distances for lags 3 and 4, in that the model generally predicts larger non-return distances. The mean adjustment results are similar (not shown). Our adjustment methods permit in-depth analyses such as these and may provide practitioners valuable information for improving the model.
The bivariate posterior distributions for the parameters based on the adjustment methods are shown in Figure 18. It is evident that the estimated posterior distributions are similar, with  the R-BSL-M posteriors slightly more concentrated than R-BSL-V. Univariate posteriors for all BSL approaches are shown in Figure 19. It appears that the adjustments result in a significant shift of the posterior and variance inflation. The standard BSL results are adversely affected by a few long periods of no MCMC acceptance. This indicates that statistical inferences are unlikely to be robust to the model misspecification.

Discussion
This paper has made two significant contributions to the literature on approximate Bayesian methods. Firstly, to our knowledge, this is the first piece of research to demonstrate that, similar to the method of approximate Bayesian computation (ABC), the method of Bayesian synthetic likelihood (BSL) can deliver unreliable inference when the assumed model is misspecified. Secondly, to circumvent the poor behavior of BSL under model misspecification, we have proposed a modification of BSL that is robust to model misspecification. Several Monte Carlo and empirical examples are used to illustrate the performance of this new method, with the results demonstrating the overwhelming performance gains of this robust approach relative to standard BSL when the model is misspecified. In addition to being robust to model misspecification, this new approach also allows the user to detect precisely which summary statistics are incompatible with the assumed data generating process. Incorporating this information within subsequent rounds of model building could lead to better models that can more accurately capture the behavior exhibited by the observed summary statistics. In this sense, the robust BSL approach can be viewed as a BSL based version of the model criticism approach of Ratmann et al. (2009). In the context of ABC, Ratmann et al. (2009) propose an approach to detect aspects of the model that the summary statistics can not adequately capture. Their approach relies on treating the ABC tolerance as an unknown parameter, and augmenting the original ABC inference problem with this additional parameter. The authors argue that posterior realizations for the tolerance parameter that are "large" indicate the possibility of a mismatch between the model and the observed data.
While useful, in the case of multivariate summaries the approach of Ratmann et al. (2009) requires a tolerance parameter for each summary statistic used in the analysis, with posterior inference then required on the full set of model parameters and tolerance parameters. Therefore, even for a moderate number of summaries, this approach can exacerbate the underlying curse-of-dimensionality in ABC, as it pertains to overall number of parameters in the analysis (Frazier et al., 2018). Moreover, since summaries differ in their responsiveness to the values used to simulate the data, the scales of the individual tolerance posteriors can vary greatly and, therefore, it is not entirely clear how to accurately determine when model misspecification is in evidence. In contrast, the approach considered herein has as a direct benchmark with which to gauge the impact of misspecification on the summaries, namely the prior distribution of the adjustment components. If the corresponding posterior for the adjustment component in the robust version of BSL does not resemble the prior, this is strong evidence that this summary can not be matched by the assumed model. If one wished to put a numerical value, or conduct a formal hypothesis test, on the difference between the prior and posterior, any number of techniques could be used; for example, the Kolmogorov-Smirnov test could be used to formally test that the prior and posterior are different.
The examples illustrated a strong indication that, in particular, the variance inflation approach can significantly improve the MCMC acceptance rate, especially under model misspecification. The variance inflation approach bears some resemblance to MCMC ABC approaches that assign a distribution to the ABC tolerance to facilitate MCMC mixing by proposing a relatively large tolerance value (e.g. Bortot et al. (2007)). However, improving mixing is not our primary focus, but is simply a useful by-product. Our R-BSL approaches may also be useful for initial explorations of the parameter space when it is not known where the bulk of the posterior support is since, even for a correctly specified model, a poor parameter value will not be able to recover the observed statistic.
Lastly, we note that, given the similarities between BSL and ABC, a natural question posed during this research was whether or not the mean and variance adjustment approaches discussed in this current paper were applicable in the context of ABC. In concurrent work from the authors, preliminary investigations into a similar type of mean and variance adjusted ABC have revealed that such an approach can mitigate the poor performance of ABC under model misspecification (Frazier et al., 2017).