Roever 2020 - DIRECT Algorithm Implementation and Calibration

Summary

bayesmeta() evaluates the 2-D posterior without MCMC: the marginal of $\tau$ is computed by one-dimensional numerical integration, and the marginal of $\mu$ (and of each ${\theta }_i$ and the prediction) is a finite mixture of normals whose support points ${\tau }_j$ and weights $w_j$ are chosen by the DIRECT algorithm (Röver and Friede 2017) to bound the Kullback-Leibler divergence by $\delta$ (default 1%). Correctness is verified by simulation-based calibration using probability integral transform (PIT) values.

Overview

The algorithm is only summarised in this paper (details: Röver and Friede 2017, Discrete approximation of a mixture distribution via restricted divergence). This note records what is specified here.

Main Content

Internals (Sec. 2.11)

• The posterior is 2-D; $p(\tau \mid y)$ known up to a constant, $p(\mu \mid \tau ,y)$ normal (Eq. 24).
• Joint density = likelihood $\times$ priors (analytic); $\tau$-marginal normalised with R’s integrate(); CDF by numerical integration; quantiles by root-finding uniroot() on the CDF.
• $\mu$-marginal $\approx {\sum }_{j}N(\hat{\mu }({\tau }_j),\hat{\sigma }({\tau }_j{)}^2)w_j$ (Eq. 26). The tuning parameters of the DIRECT algorithm are delta ($\delta$, bound on divergence) and epsilon ($ϵ$). In the example, 17 components suffice to keep the Kullback-Leibler divergence below $\delta =1\%$; they are stored in ...$support.
• Mixture density and CDF are linear combinations of component densities/CDFs; random generation and moments are simple (Seidel 2010; Lindsay 1995). Heavy reliance on integrate(), uniroot(), optimize(), optim().

Symmetrised divergence between two normals (App. A.4, Eq. 45)

For $N({\mu }_A,{\sigma }_A^2)$ and $N({\mu }_B,{\sigma }_B^2)$:

$$D_s(p_A\Vert p_B)=\frac{({\mu }_A-{\mu }_B{)}^2}{2}\left(\frac{1}{{\sigma }_A^2}+\frac{1}{{\sigma }_B^2}\right)+\frac{({\sigma }_A^2-{\sigma }_B^2{)}^2}{2{\sigma }_A^2{\sigma }_B^2}.$$

The DIRECT algorithm uses this to decide how close neighbouring ${\tau }_j$ must be. For the shrinkage distributions the maximum divergence across the effect posterior and all $k$ shrinkage distributions is used so a single $\tau$-grid serves all three mixture types.

(The paper’s Eq. 45 is typeset with nested fractions that did not extract cleanly; the closed form above is the standard symmetrised Kullback-Leibler divergence $KL(A\Vert B)+KL(B\Vert A)$ for two normals, which is the quantity the text names.)

Calibration check (App. A.5)

PIT uniformity under correct Bayesian computation (Cook, Gelman and Rubin 2006; Gneiting et al. 2007)

If ${\theta }^{\ast }\sim p(\theta )$, $y^{\ast }\sim p(y\mid {\theta }^{\ast })$ and $p^{\ast }=P(\theta \le {\theta }^{\ast }\mid y^{\ast })$ is the posterior CDF at the true value, then $p^{\ast }\sim U(0,1)$ exactly for any sample size and any proper prior. Hence uniformity of PIT values over many simulations validates software. Bayesian credible intervals are calibrated on average over the prior (versus frequentist coverage at each parameter value).

Simulation protocol (1000 replicates): $\mu \sim N(0,4^2)$, $\tau \sim \text{half-normal}(0.5)$, $k\in \{2,3,5,10,20\}$ drawn at random, ${\sigma }_i\sim U(0.2,1.0)$, $y_i\sim N(\mu ,{\sigma }_i^2+{\tau }^2)$; analyse with matching priors; record pposterior(mu=mu) and pposterior(tau=tau). Fig. 8: empirical CDFs of PIT values for $\mu$ and $\tau$ are in close agreement with the uniform. One-sided 95% credible limits show about 95% empirical coverage. For $\tau$, a curve above the diagonal reads as conservative (tendency to overestimate heterogeneity); for $\mu$, a conservative procedure sits below the diagonal at the low end and above at the high end (wider intervals).

Examples

import numpy as np
from scipy import stats
def calibration_check(posterior_cdf_mu, nsim=1000, seed=1):
    """PIT of mu under the prior-generating process of Röver's App. A.5.
    posterior_cdf_mu(y, sigma) -> callable F(mu) must be supplied by the analysis under test."""
    rng = np.random.default_rng(seed); pit = []
    for _ in range(nsim):
        mu = rng.normal(0, 4); tau = abs(rng.normal(0, 0.5)); k = rng.choice([2,3,5,10,20])
        sigma = rng.uniform(0.2, 1.0, k); y = rng.normal(mu, np.sqrt(sigma**2+tau**2))
        pit.append(posterior_cdf_mu(y, sigma)(mu))
    return np.array(pit)          # should be ~ U(0,1): check with stats.kstest(pit, "uniform")

Connections

• Marketing: any in-house Bayesian aggregation tool (a pooled-incrementality model, an MMM prior builder) should be validated this way (prior-predictive data in, PIT uniformity out); see Posterior Sampling and Summarization and the validation cluster in Model Validation and Interpretability - Overview.

See Also

• Roever 2020 - Pediatric Liver Transplant Worked Example
• Roever 2020 - Bayesian Random-Effects Meta-Analysis with bayesmeta - Overview