Synthetic Likelihood for Intractable QSP Models

Why Use a Synthetic Likelihood?

A QSP simulator may generate realistic time courses or biomarker measurements while making the probability of the observed data impossible to calculate directly. Calibration still requires a way to compare a parameter value with the observations. One option is to compare selected summaries of simulated and observed data. Synthetic likelihood assigns a probability density to the observed summaries by approximating their simulator-induced distribution with a multivariate normal distribution.

This approach is useful when summaries are informative and their joint distribution is reasonably close to Gaussian at each parameter value of interest. It is not an exact likelihood for the raw data: it is an approximate likelihood for the chosen summaries. The introductory material on likelihood-free calibration established why such an approach may be needed; here the focus is constructing and using the Gaussian approximation.

Definition: Let \(\theta\) denote a vector of model parameters, let \(\mathbf{s}_{obs}\) be a \(d\)-dimensional vector of summaries calculated from observed data, and let \(\mathbf{S}(\theta)\) be the corresponding vector calculated from a data set simulated at \(\theta\). A Gaussian synthetic likelihood assumes \(\mathbf{S}(\theta)\) is approximately multivariate normal, with mean \(\boldsymbol{\mu}(\theta)\) and covariance matrix \(\boldsymbol{\Sigma}(\theta)\). The synthetic likelihood is the density of that normal distribution evaluated at \(\mathbf{s}_{obs}\).

The summaries must be defined identically for observed and simulated data. They might include a biomarker value at a prespecified time, a maximum response, or a measure of between-visit variability. Their scale and scientific meaning matter: changing the summaries changes the target being calibrated.

Estimating and Evaluating the Likelihood

At a fixed parameter value \(\theta\), generate \(m\) independent simulator replicates and calculate a summary vector from each one, denoted \(\mathbf{S}^{(1)}(\theta),\ldots,\mathbf{S}^{(m)}(\theta)\). Estimate the conditional mean and covariance using the sample mean and sample covariance:

\widehat{\boldsymbol{\mu}}_\theta = \frac{1}{m}\sum_{r=1}^{m}\mathbf{S}^{(r)}(\theta), \qquad \widehat{\boldsymbol{\Sigma}}_\theta = \frac{1}{m-1}\sum_{r=1}^{m} \left(\mathbf{S}^{(r)}(\theta)-\widehat{\boldsymbol{\mu}}_\theta\right) \left(\mathbf{S}^{(r)}(\theta)-\widehat{\boldsymbol{\mu}}_\theta\right)^{T}.

Using the sample covariance, whose divisor is \(m-1\), is a conventional choice for estimating the covariance. If this estimate is positive definite, the plug-in Gaussian synthetic likelihood is

\widehat{L}_{SL}(\theta) = \frac{ \exp\!\left[ -\frac{1}{2} (\mathbf{s}_{obs}-\widehat{\boldsymbol{\mu}}_\theta)^T \widehat{\boldsymbol{\Sigma}}_\theta^{-1} (\mathbf{s}_{obs}-\widehat{\boldsymbol{\mu}}_\theta) \right] }{ (2\pi)^{d/2}\left|\widehat{\boldsymbol{\Sigma}}_\theta\right|^{1/2} }.

The quadratic term measures the discrepancy from the simulated mean in units of the estimated joint variability. The determinant term is essential: it accounts for the volume and scale of the simulated distribution. A parameter value is not automatically preferred merely because it has a larger covariance and therefore places the observed summaries close to its mean.

For computation, it is usually more stable to work with the log synthetic likelihood:

\log \widehat{L}_{SL}(\theta) = -\frac{d}{2}\log(2\pi) -\frac{1}{2}\log\left|\widehat{\boldsymbol{\Sigma}}_\theta\right| -\frac{1}{2} (\mathbf{s}_{obs}-\widehat{\boldsymbol{\mu}}_\theta)^T \widehat{\boldsymbol{\Sigma}}_\theta^{-1} (\mathbf{s}_{obs}-\widehat{\boldsymbol{\mu}}_\theta).

Why the Gaussian Density Has This Form

The formula follows by transforming a vector of independent standard normal variables. Suppose \(\mathbf{Z}\) has independent \(N(0,1)\) components and \(\boldsymbol{\Sigma}\) is positive definite. Let \(\boldsymbol{\Sigma}=\mathbf{L}\mathbf{L}^T\), where \(\mathbf{L}\) is a lower triangular Cholesky factor, and define \(\mathbf{S}=\boldsymbol{\mu}+\mathbf{L}\mathbf{Z}\). Then \(\mathbf{S}\) has mean \(\boldsymbol{\mu}\) and covariance \(\boldsymbol{\Sigma}\). The inverse transformation is \(\mathbf{Z}=\mathbf{L}^{-1}(\mathbf{S}-\boldsymbol{\mu})\), whose density transformation contributes the Jacobian factor \(1/|\mathbf{L}|\). Since \(|\boldsymbol{\Sigma}|=|\mathbf{L}|^2\), that factor is \(1/|\boldsymbol{\Sigma}|^{1/2}\), and \(\mathbf{Z}^T\mathbf{Z}=(\mathbf{S}-\boldsymbol{\mu})^T\boldsymbol{\Sigma}^{-1}(\mathbf{S}-\boldsymbol{\mu})\). Substituting the standard normal density gives the stated multivariate normal density.

For an approximate Bayesian calibration, multiply the synthetic likelihood by the prior density. For a continuous parameter vector, the resulting approximate posterior density is proportional to that product; for discrete candidate values, normalize the prior-weighted likelihoods to sum to one. The posterior is approximate because the summary distribution has been approximated and its mean and covariance have been estimated by simulation.

Worked Example: A One-Summary Synthetic Likelihood

Suppose the observed summary is a hypothetical mean biomarker change of \(11\) units. At a candidate parameter value, five independent simulations produce summary values \(8, 9, 10, 11,\) and \(12\) units. Calculate the Gaussian synthetic likelihood.

The simulated mean is \(\widehat{\mu}_\theta=(8+9+10+11+12)/5=10\) units. The sum of squared deviations is \(4+1+0+1+4=10\) units squared, so the sample variance is \(\widehat{\sigma}^2_\theta=10/(5-1)=2.5\) units squared.

For one summary, the density is

\widehat{L}_{SL}(\theta) = \frac{1}{\sqrt{2\pi(2.5)}} \exp\left[-\frac{(11-10)^2}{2(2.5)}\right] = 0.2523\times 0.8187 \approx 0.2066\ {\rm units}^{-1}.

Thus the estimated density at the observed summary is approximately \(0.2066\) per unit. It is a density, not the probability that the summary equals exactly \(11\); a density can exceed one, and probabilities require an interval or region.

Worked Example: A Two-Summary Synthetic Likelihood

Suppose two hypothetical observed summaries are \(110\) units for a biomarker level and \(2.4\) hours for a response-time measure. Simulations at one parameter value give estimated mean \(\widehat{\boldsymbol{\mu}}_\theta=(100,2.0)^T\) and covariance

\widehat{\boldsymbol{\Sigma}}_\theta = \begin{pmatrix} 100 & 4\\ 4 & 0.25 \end{pmatrix}.

The determinant is \(100(0.25)-4(4)=9\), which is positive; the diagonal variances are also positive, so this symmetric \(2\times2\) matrix is positive definite. The inverse is \(\frac{1}{9}\begin{pmatrix}0.25&-4\\-4&100\end{pmatrix}\). The difference between the observed summaries and their estimated mean is \(\mathbf{d}=(10,0.4)^T\). Therefore, the quadratic discrepancy is

\mathbf{d}^T\widehat{\boldsymbol{\Sigma}}_\theta^{-1}\mathbf{d} = \frac{0.25(10)^2-8(10)(0.4)+100(0.4)^2}{9} = \frac{9}{9}=1.

The density is consequently

\widehat{L}_{SL}(\theta) = \frac{\exp(-1/2)}{(2\pi)^{2/2}\sqrt{9}} = \frac{0.6065}{6\pi} \approx 0.03219\ {\rm (unit\cdot hour)}^{-1}.

The off-diagonal covariance is used in the joint discrepancy; treating the summaries as independent would discard their estimated association and produce a different likelihood.

Worked Example: Combining Synthetic Likelihoods With Prior Probabilities

Consider two hypothetical candidate parameter values, \(A\) and \(B\), and a single observed summary of \(1\) unit. Three simulator replicates at \(A\) give \(-1,0,1\); three at \(B\) give \(0.5,1.5,2.5\). Each set has sample mean \(0\) at \(A\) or \(1.5\) at \(B\), and sample variance \(1\) unit squared. Let the prior probabilities be \(P(A)=0.6\) and \(P(B)=0.4\).

The synthetic likelihood at \(A\) is \((2\pi)^{-1/2}\exp(-1/2)\approx0.2420\) per unit. At \(B\), it is \((2\pi)^{-1/2}\exp(-0.5^2/2)\approx0.3521\) per unit. The prior-weighted values are \(0.6(0.2420)=0.1452\) for \(A\) and \(0.4(0.3521)=0.1408\) for \(B\). Normalizing gives

P(A\mid s_{obs})\approx\frac{0.1452}{0.1452+0.1408}=0.5076, \qquad P(B\mid s_{obs})\approx\frac{0.1408}{0.1452+0.1408}=0.4924.

Although \(B\) has the larger synthetic likelihood, the prior favors \(A\). After combining both contributions, the approximate posterior probabilities are close, with \(A\) slightly more probable under this two-candidate model.

Assumptions, Diagnostics, and Practical Choices

The approximation is most defensible when the simulator replicates are generated under the intended model, the summaries are calculated consistently, and their joint conditional distribution at each relevant \(\theta\) is reasonably represented by a normal distribution. If the model includes measurement error or other observation variability, the simulation must represent it when it is part of the data-generating process. Otherwise the estimated covariance may omit real variability.

Several checks are important before relying on the resulting calibration:

Inspect summary distributions. At representative parameter values, check simulated summary histograms or quantile plots for marked skewness, heavy tails, multiple modes, or boundary effects. A Gaussian approximation may be poor when these features are prominent.
Check covariance stability. The sample covariance must be nonsingular for the ordinary density formula. With \(m\) simulated vectors in \(d\) dimensions, its rank cannot exceed \(m-1\); therefore \(m>d\) is necessary, but often not sufficient for a stable estimate.
Increase simulation replicates where needed. Repeated estimates of the mean, covariance, and log synthetic likelihood can reveal substantial Monte Carlo variation. Small \(m\) makes the estimated likelihood noisy, even when the Gaussian assumption itself is reasonable.
Check the summaries’ scales and dependence. Standardization can help numerical conditioning, but the same transformation must be applied to observed and simulated summaries. Correlations should be represented by the covariance, not silently discarded.

When covariance estimates are nearly singular, possible remedies include increasing \(m\), reducing redundant summaries, or using a stated regularization such as a shrinkage estimate. For example, a convex combination of the estimated covariance and its diagonal can improve conditioning. Regularization changes the likelihood being used, so its choice should be reported and checked for sensitivity. In computation, a Cholesky factorization is generally preferable to explicitly forming a matrix inverse.

How It Differs From ABC, and Where It Can Fail

ABC methods compare simulated and observed summaries using a distance and a tolerance or kernel. Gaussian synthetic likelihood instead estimates a mean and covariance for the summaries and evaluates a continuous density. This can avoid selecting a tolerance, while accounting for correlations and different summary scales through the covariance. ABC can be attractive when the summary distribution is strongly non-Gaussian or when only a distance-based comparison is trusted. Synthetic likelihood can be more efficient when a Gaussian approximation is adequate and the covariance can be estimated reliably.

Neither method removes the consequences of choosing summaries. Calibration may be weak if important features of the data are omitted, and conclusions concern agreement with the selected summaries rather than every feature of the full data. A Gaussian fit can also be misleading for bounded summaries, rare events, or a distribution with distinct modes. Finally, likelihood estimates vary because they use a finite number of simulator runs. Treating a noisy plug-in estimate as exact can make the computed posterior appear more precise than the simulation supports.

Common Pitfall

Common mistake: “The observed summaries are close to the simulated mean, so this parameter value must have the largest synthetic likelihood.”

Correction: The likelihood depends on both the discrepancy and the covariance normalization. The quadratic term measures discrepancy relative to covariance, while the determinant term penalizes the density’s spread. Evaluate the complete density, and verify that the covariance estimate is stable and positive definite.

Specify summaries and the simulator.
Define \(\mathbf{s}_{obs}\) and compute the identical summaries from each simulated data set, including relevant observation variability.

Estimate the conditional distribution.
At each proposed \(\theta\), simulate replicates and calculate \(\widehat{\boldsymbol{\mu}}_\theta\) and \(\widehat{\boldsymbol{\Sigma}}_\theta\).

Evaluate and diagnose.
Compute the log Gaussian density, check normality and covariance conditioning, and increase simulation effort or revise the summaries if estimates are unstable.

Calibrate with the prior.
Combine the synthetic likelihood with the prior to obtain an approximate posterior, then interpret it as conditional on the model, summaries, and Gaussian approximation.

Check Your Understanding

Answer the questions, then compare with the solutions below.

Why must the observed and simulated summaries be defined in the same way?
For \(d=4\) summaries and \(m=3\) independent simulated summary vectors, can the usual sample covariance be nonsingular? Explain.
What does the determinant term contribute to a Gaussian synthetic likelihood?
A simulated summary distribution is strongly bimodal. Why is a single Gaussian synthetic likelihood questionable?

Solutions. (1) The synthetic likelihood compares the observed vector with the distribution of the same quantity under simulation; changing the summary definition makes that comparison invalid. (2) No. The covariance rank is at most \(m-1=2\), which is less than \(d=4\), so it is singular. (3) It supplies the density normalization for the covariance scale and volume; omitting it would misstate the probability density and could favor overly dispersed estimates. (4) A single normal distribution represents one approximately bell-shaped joint distribution, not two distinct modes, and may assign misleading density between or around the modes.