Composite pseudo-likelihood functions

For approaches that require a likelihood function, we must select a functional form appropriate to the data. Recall that the likelihood we seek to evaluate is \(\mathcal L(X, \theta)\). This conditional probability, viewed as a function of \(X\) and \(\theta\) describes the loss; essentially, how well the simulation matches the data.

Because we rarely know the true likelihood function, common practice in Bayesian analysis is to use a pseudo-likelihood function, selecting the functional form appropriate to the nature of the observed data \(y^{\text{obs}}\). In doing so, we often make the simplifying assumption that various different data sources are independent, resulting in a composite pseudo-likelihood composed of \(J\) independent components,

\[ \mathcal L(X, \theta) = \prod_{j=1}^{J} \mathcal L_{j}(X, \theta). \]

Applied examples of model calibration, available in other documents, require knowledge of a few common likelihood forms. After briefly describing how to compute likelihoods, we formally introduce these common forms, explaining use cases, applications, and notes for each.

In practice, you may include weights \(w_j>0\) in the composite log-likelihood \(\ell(\theta)=\sum_j w_j\,\ell_j(\theta)\) to balance scales and credibility across data sources.

Note that composite likelihoods can underestimate uncertainty by ignoring cross-component dependence. An advanced technique to address this issue is to use sandwich/Godambe adjustments ¹.

Computing likelihoods

A few different computational approaches are available to estimate pseudo-likelihoods.

The most straightforward approach is Monte Carlo sampling of a selected parameter configuration, \(\theta^{(k)}\). Note, however, that doing so effectively integrates out the latent state trajectory \(X\) by averaging over \(N\) independent realizations of the simulator, \(\tilde x^{(r)} \sim p(\cdot \mid \theta^{(k)})\) for \(r=1,\ldots,N\). The result is an unbiased estimate of the marginal likelihood at \(\theta^{(k)}\):

\[ \mathcal {\widehat L}(\theta^{(k)}) = \frac{1}{N}\sum_{r=1}^N \mathcal L(\tilde x^{(r)}, \theta^{(k)}). \]

The hat over \(\mathcal{\widehat L}\) indicates that this is an estimate to the marginal likelihood \(\mathcal L(\theta^{(k)}) = p(y^{\text{obs}} \mid \theta)\). The estimate is unbiased, but subject to Monte Carlo error that diminishes as the number of simulation replicates increases.

When using Monte Carlo to estimate a composite likelihood, it is important to take the average outside of the product:

\[ \mathcal {\widehat L}(\theta^{(k)}) = \frac{1}{N}\sum_{r=1}^N \prod_{j=1}^J \mathcal L_j(\tilde x^{(r)}, \theta^{(k)}). \]

Alternatively, you can use a particle filter, an SMC method for state-space models. As with Monte Carlo (above), this approach requires running the model many (\(N\)) times. However, here, instead of running simulations to completion, we "resample" when the weights (partial likelihoods) have high variance. In doing so, some simulations will be terminated and others will be duplicated (with random seed perturbation or "rejuvenation"). The result is often a lower-variance estimate than Monte Carlo at the cost of increased complexity, "broken" paths due to resampling, and additional bookkeeping for reproducibility.

The broken-path problem of SMC can be rectified using particle Gibbs methods that always maintain one complete "reference" simulation.

Computationally, it's often more efficient to work with the log of the likelihood function, here denoted \(\ell(X, \theta)\), in which case the composite is formed by summing logs. Please use log-sum-exp to to avoid incorrectly averaging log-likelihoods Monte Carlo replicates.

Common pseudo-likelihood forms

Normal

Use case

The Normal (Gaussian) distribution is a common choice for continuous data, especially when measurement errors are approximately symmetric and data represent averages of many independent events. By the Central Limit Theorem, means of repeated samples tend toward normality even if the underlying process is not normal.

Application

Let the observed data \(y^{\text{obs}}\) be a continuous value (or vector of values). For each replicate \(i\) at input \(\theta^{(i)}\), the simulation produces latent trajectory \(\tilde x^{(i)}\) and output \(\tilde y^{(i)}\). We center the normal distribution at the simulated output and set the variance to reflect the measurement uncertainty in the observed data:

\[ \mu = \tilde y^{(i)}, \qquad \sigma^2 = \sigma^2_{\mathrm{obs}}. \]

If a confidence interval is reported with the observed mean, \(\sigma_{\mathrm{obs}}\) can be back-calculated.

PDF and kernel

\[ \begin{align} \mathcal{N}(x\mid \mu,\sigma^2) &= \frac{1}{\sqrt{2\pi\sigma^2}}\exp\!\left(-\frac{(x-\mu)^2}{2\sigma^2}\right),\\ \mathcal{L}(\tilde x^{(i)},\theta^{(i)}) &= \mathcal{N}\!\big(y^{\text{obs}} \mid \tilde y^{(i)}, \sigma^2\big). \end{align} \]

Notes

In log space,

\[ \ell(\tilde x^{(i)},\theta^{(i)}) = -\tfrac{1}{2}\log(2\pi\sigma^2) - \tfrac{(y^{\text{obs}}-\tilde y^{(i)})^2}{2\sigma^2}. \]

For stochastic models, variation across replicates is integrated out via Monte Carlo, particle filtering, or other methods, not added into the variance, \(\sigma^2\).
When performing parameter calibration using optimization methods to find the best parameters, maximizing the normal likelihood is equivalent to minimizing squared error. To see this:

\[ \begin{align*} \ell(\tilde x^{(i)}, \theta^{(i)}) &= \log \mathcal{N}(y^{\text{obs}} \mid \tilde y^{(i)}, \sigma^2_{\mathrm{obs}}) \\ &= -\tfrac{1}{2}\log(2\pi\sigma^2_{\mathrm{obs}}) - \frac{(\tilde y^{(i)} - y^{\text{obs}})^2}{2\sigma^2_{\mathrm{obs}}}. \end{align*} \]

Since the first term does not depend on \(\theta\), maximizing \(\ell(\tilde x^{(i)}, \theta^{(i)})\) is equivalent to minimizing \((\tilde y^{(i)}-y^{\text{obs}})^2\). Thus, normal likelihood maximization and least-squares fitting yield the same "best" parameter estimates.

Beta–binomial

Use case

Extends the binomial to allow for over-dispersion (e.g., cluster or survey heterogeneity). The beta binomial distribution is appropriate when the data are “successes out of trials.” An example is the number of positive test results out of total tests conducted in a clinical trial or population survey.

Application

Let the observed data be \(y^{\text{obs}} = (s,n)^{\text{obs}}\) representing \(s^{\text{obs}}\) successes out of \(n^{\text{obs}}\) trials. For each simulation replicate \(i\), the simulation provides \(\tilde s^{(i)}\) successes out of \(\tilde n^{(i)}\) trials, from which we form a probability, \(\tilde p^{(i)}\).

No smoothing: \(\tilde p^{(i)} = \tilde s^{(i)}/\tilde n^{(i)}\)
Laplace smoothing: \(\tilde p^{(i)} = (\tilde s^{(i)}+1)/(\tilde n^{(i)}+2)\)

Let \(\kappa>0\) be a per-observation concentration (e.g., derived from reported CIs); then

\[ a_i = \kappa\,\tilde p^{(i)}, \qquad b_i = \kappa\,(1-\tilde p^{(i)}). \]

PMF and kernel

\[ \begin{align} s^{\text{obs}} &\sim \mathrm{BetaBinomial}\!\big(n^{\text{obs}}, a_i, b_i\big)\\ \mathcal{L}(\tilde x^{(i)},\theta^{(i)}) &= \mathrm{BetaBinomial}\!\big(s^{\text{obs}} \mid n^{\text{obs}}, a_i, b_i\big). \end{align} \]

Notes

Larger \(\kappa\) means less over-dispersion; as \(\kappa\to\infty\), this approaches the binomial with probability \(\tilde p^{(i)}\).
If \(\kappa\) is unavailable, use the Binomial kernel.

Gamma–Poisson (negative-binomial)

Use case

The gamma–Poisson (negative binomial) models counts that are over-dispersed relative to the Poisson. The variance exceeds the mean, common in clustered or heterogeneous event data. It arises when the Poisson rate \(\lambda\) has a gamma-distributed prior.

Application

From the data, we observe count \(c^{\text{obs}}\) over "exposure" \(e^\text{obs}\) that could represent person-time. For each replicate \(i\), the simulation produces \(\tilde c^{(i)}\) events over exposure \(\tilde e^{(i)}\). Then compute,

\[ \tilde r^{(i)} = \frac{\tilde c^{(i)}}{\tilde e^{(i)}}, \qquad \mu_i = e^{\text{obs}}_i\,\tilde r^{(i)}. \]

Let \(k>0\) be a per-observation size (concentration). Use the negative-binomial parameterization with mean \(\mu_i\) and variance \(\mu_i + \mu_i^2/k\). In the \((n,p)\) form,

\[ p_i = \frac{k}{k + \mu_i} \]

PMF and kernel

\[ \begin{align} c^{\text{obs}}_i &\sim \mathrm{NegBin}(n=k, p=p_i)\\ \mathcal{L}(\tilde x^{(i)},\theta^{(i)}) &= \mathrm{NegBin}\!\big(c^{\text{obs}}_i \mid n=k, p=p_i\big). \end{align} \]

Notes

As \(k\to\infty\), this converges to the Poisson with mean \(\mu_i\).
If \(k\) is unavailable, use the Poisson kernel.

Dirichlet–multinomial

Use case

Counts across \(K\) categories with over-dispersion relative to the multinomial.

Application

Observed \(\mathbf{c}^{\text{obs}} = (c^{\text{obs}}_1,\ldots,c^{\text{obs}}_K)\) with \(n^{\text{obs}}=\sum_j c^{\text{obs}}_j\). For replicate \(i\), from simulated counts \(\tilde{\mathbf c}^{(i)}\) form smoothed proportions to avoid zero cells:

\[ \hat{\boldsymbol{\pi}}^{(i)} = \frac{\tilde{\mathbf c}^{(i)} + \alpha_s\,\mathbf 1}{\sum_j \tilde c^{(i)}_j + K\alpha_s}, \quad \alpha_s>0. \]

Let \(\kappa>0\) be a per-observation concentration and set

\[ \boldsymbol{\alpha}^{(i)} = \kappa\,\hat{\boldsymbol{\pi}}^{(i)}. \]

PMF and kernel

\[ \begin{align} \mathbf{c}^{\text{obs}} &\sim \mathrm{DirMult}\!\big(n^{\text{obs}},\,\boldsymbol{\alpha}^{(i)}\big)\\ \mathcal{L}(\tilde x^{(i)},\theta^{(i)}) &= \mathrm{DirMult}\!\big(\mathbf{c}^{\text{obs}} \mid n^{\text{obs}}, \boldsymbol{\alpha}^{(i)}\big). \end{align} \]

Notes

Larger \(\kappa\) reduces over-dispersion; as \(\kappa\to\infty\) this approaches the multinomial with probabilities \(\hat{\boldsymbol{\pi}}^{(i)}\).
If \(\kappa\) is unavailable, use the multinomial kernel:

\[ \mathbf{c}^{\text{obs}} \sim \mathrm{Multinomial}\!\big(n^{\text{obs}},\,\hat{\boldsymbol{\pi}}^{(i)}\big). \]

Evaluate in log space with lgamma for stability.

References

Godambe, V. P. (1964). A unified theory of estimating equations. Biometrika, 50(1-2), 177–185. https://doi.org/10.1093/biomet/50.1-2.177 ↩