Gibbs Sampling

Gibbs Sampling in Bayesian Evidence Synthesis and Decision Models

Gibbs sampling simulates a Bayesian posterior distribution by breaking one hard sampling problem into a sequence of easier ones, updating one parameter or block of parameters at a time while the rest stay at their latest values. After enough cycles the stored values behave like correlated draws from the joint posterior. In health economics it sits behind many of the Bayesian evidence syntheses that feed decision models, because the BUGS family of software, used for the meta-analysis code in NICE Decision Support Unit Technical Support Document 2 (TSD 2), was built around it. Its draws are dependent, and strong correlation between parameters can make the chain crawl.

One parameter at a time from its full conditional

Gelman and colleagues, in Bayesian Data Analysis, describe the Gibbs sampler as a Markov chain Monte Carlo method that divides the parameter vector into components and, in each iteration, draws each component conditional on the value of all the others. With $d$ components there are $d$ steps per iteration.

$$\theta_j^{(t)} \sim p\left(\theta_j \mid \theta_1^{(t)}, \ldots, \theta_{j-1}^{(t)}, \theta_{j+1}^{(t-1)}, \ldots, \theta_d^{(t-1)}, y\right), \qquad j = 1, \ldots, d$$

where $\theta_j$ is the $j$th component or block of the parameter vector, $t$ is the iteration number, $y$ is the observed data, and the conditioning set holds components already updated in iteration $t$ at their new values and the rest at their iteration $t-1$ values. The right-hand side is the full conditional, the posterior distribution of $\theta_j$ with everything else treated as known.

After warm-up, posterior quantities are estimated by averaging over the retained draws.

$$\hat{E}\left[h(\theta) \mid y\right] = \frac{1}{M} \sum_{t=1}^{M} h\left(\theta^{(t)}\right)$$

where $M$ is the number of draws kept after warm-up, $\theta^{(t)}$ is the parameter vector at retained iteration $t$, and $h$ is any function of interest, such as a relative treatment effect or an incremental net benefit.

Why no draw is ever rejected

Gelman and colleagues show that Gibbs sampling is a special case of Metropolis-Hastings in which the proposal for each step is the full conditional itself. The acceptance ratio then equals 1 at every step, so every jump is accepted. A Gibbs sampler therefore has no proposal scale to tune, at the price that each full conditional must be available in a form that can be sampled.

Conditional conjugacy and the models that suit it

Many models in health economics are built as a chain of conditional distributions, such as trial-level effects drawn from a population distribution in a hierarchical model. Gelman and colleagues note that the conditional distributions in such models are often conjugate and easy to simulate, and they call the Gibbs sampler their first choice for conditionally conjugate models. Where some conditionals are standard and others are not, the sampler can use exact Gibbs draws where possible and one-dimensional Metropolis updates elsewhere.

Plummer, introducing JAGS, notes that BUGS stands for Bayesian inference Using Gibbs Sampling and that it made models with random effects, measurement error and missing data accessible to a wide audience.

Unobserved quantities can also be added to the parameter vector, which Gelman and colleagues call data augmentation, with indicators for a mixture model as one example. They trace the name of the sampler to Geman and Geman, whose 1984 image restoration paper made local changes by sampling from a local conditional probability distribution.

Worked example with an illustrative mean cost

All numbers in this example are illustrative. A cost model needs the mean annual cost per patient of a new care pathway, in GBP thousands, from $n = 20$ patients with sample mean 4.0 and standard deviation 2.0. Costs are modelled as normal with unknown mean $\mu$ and precision $\phi$, the inverse of the variance. The prior for $\mu$ is normal with mean 3.0 and standard deviation 2.0, and the prior for $\phi$ is gamma with shape 1 and rate 1. With these two independent priors the joint posterior has no standard form, but each full conditional does.

$$\mu \mid \phi, y \sim N\left(\frac{p_0 \mu_0 + n \phi \bar{y}}{p_0 + n \phi}, \frac{1}{p_0 + n \phi}\right)$$

where $\mu_0 = 3.0$ is the prior mean, $p_0 = 1/2.0^2 = 0.25$ is the prior precision, $\bar{y} = 4.0$ is the sample mean, and the second argument is the conditional variance. This is the textbook normal result, in which posterior precision is prior plus data precision and the posterior mean is a precision-weighted average.

$$\phi \mid \mu, y \sim \mathrm{Gamma}\left(a + \frac{n}{2}, b + \frac{S(\mu)}{2}\right), \qquad S(\mu) = (n-1)s^2 + n(\bar{y} - \mu)^2$$

where $a = 1$ and $b = 1$ are the prior shape and rate, $s = 2.0$ is the sample standard deviation, so $(n-1)s^2 = 19 \times 4 = 76$, and $S(\mu)$ is the sum of squared deviations of the costs from the current value of $\mu$.

1. Iteration 1, update the mean. Starting from $\phi = 0.20$, a variance of 5, the conditional precision is 0.25 + 20 × 0.20 = 4.25. The numerator of the conditional mean is 0.25 × 3.0 + 20 × 0.20 × 4.0 = 16.75, so the mean is 16.75 / 4.25 = 3.941, with standard deviation $\sqrt{1/4.25} = 0.485$. A standard normal draw of 0.40 gives $\mu$ = 3.941 + 0.40 × 0.485 = 4.135.

2. Iteration 1, update the precision. With $\mu$ = 4.135 the gap to the sample mean is 0.135, whose square is 0.0182, so $S$ = 76 + 20 × 0.0182 = 76.36. The gamma conditional has shape 1 + 10 = 11 and rate 1 + 76.36 / 2 = 39.18, so its mean is 11 / 39.18 = 0.281. Suppose the draw is 0.270, a variance of 1 / 0.270 = 3.70.

3. Iteration 2, update the mean. The precision is now 0.25 + 20 × 0.270 = 5.65. The numerator is 0.75 + 20 × 0.270 × 4.0 = 22.35, so the conditional mean is 22.35 / 5.65 = 3.956, with standard deviation $\sqrt{1/5.65} = 0.421$. A standard normal draw of minus 1.10 moves the mean down by 1.10 × 0.421 = 0.463, giving $\mu$ = 3.956 minus 0.463, or 3.493.

4. Iteration 2, update the precision. The gap to the sample mean is 0.507, whose square is 0.257, so $S$ = 76 + 20 × 0.257 = 81.14 and the gamma conditional has shape 11 and rate 1 + 81.14 / 2 = 41.57, with mean 11 / 41.57 = 0.265.

5. Run and summarise. The exact marginal posterior, found by numerical integration, has mean 3.95 for $\mu$ and a 95% central credible interval of 3.07 to 4.82, that is GBP 3,070 to GBP 4,820, and one run of 60,000 iterations with the first 10,000 discarded as warm-up matched each figure to within 0.01. A rough check agrees. At a cost precision near 0.25 the data precision is about 20 × 0.25 = 5 against a prior precision of 0.25, so the data carry weight 5 / 5.25 = 0.952 and 0.952 × 4.0 + 0.048 × 3.0 = 3.95.

The mean is pulled slightly below the sample mean of 4.0 towards the prior mean of 3.0. Each retained value of $\mu$ already averages over the uncertainty in $\phi$, and those draws, not a summary mean and interval, are what a probabilistic decision model would take as its cost input. Another run would give slightly different decimals because of Monte Carlo error.

Correlated parameters slow the chain

Updating one component at a time works badly when components are strongly correlated in the posterior. Gelman and colleagues show this with a bivariate normal posterior with unit variances, correlation $\rho$ and data at zero, where the two full conditionals are normal with means $\rho\theta_2$ and $\rho\theta_1$ and variance one minus $\rho^2$. Substituting one update into the next shows how much each sweep remembers.

$$\theta_1^{(t+1)} = \rho\,\theta_2^{(t)} + e_1, \qquad \theta_2^{(t)} = \rho\,\theta_1^{(t)} + e_2 \quad \Rightarrow \quad \theta_1^{(t+1)} = \rho^2\,\theta_1^{(t)} + \rho e_2 + e_1$$

where $e_1$ and $e_2$ are independent normal errors with mean zero and variance one minus $\rho^2$. The chain for $\theta_1$ therefore has autocorrelation $r_k = \rho^{2k}$ at lag $k$. Gelman and colleagues define the effective number of draws in terms of these autocorrelations.

$$n_{\mathrm{eff}} = \frac{N}{1 + 2\sum_{k=1}^{\infty} r_k} = N \, \frac{1 - \rho^2}{1 + \rho^2}$$

where $N$ is the number of retained draws, $r_k$ is the lag-$k$ autocorrelation and $n_{\mathrm{eff}}$ is the number of independent draws carrying the same information about the mean. The second form uses the sum of the geometric series $\rho^2 / (1 - \rho^2)$.

Spiegelhalter notes that baseline levels and treatment differences estimated from the same studies are likely to be strongly correlated, and TSD 2 advises choosing a reference treatment at the centre of a network meta-analysis to reduce correlations between treatment effects that can slow convergence. Gelman and colleagues note that the Gibbs sampler is most efficient when parameterised in terms of independent components. Remedies include linear reparameterisation, updating correlated parameters as one block, parameter expansion, or Hamiltonian Monte Carlo, which they describe as mixing faster in high dimensions. Diagnostics for MCMC output are covered under convergence.

Gibbs sampling in HTA evidence synthesis

TSD 2 notes that Bayesian MCMC has for many years been a mainstay of decision analysis and gives code for WinBUGS 1.4.3. Its beta-blocker results are based on 20,000 iterations on 3 chains after a burn-in of 10,000.

Posterior sampling also keeps the correlation induced by estimating parameters from the same data, and TSD 2 states that where the evidence contains closed loops the posterior samples should be used in the decision model, or the correlations otherwise propagated, because distributions based on the marginal summaries are not adequate. Spiegelhalter calls this the integrated approach, generally implemented with MCMC, and contrasts it with fitting distributions first and then propagating them in a probabilistic sensitivity analysis. The integrated route needs no parametric shapes for the posteriors and carries the dependence between unknowns.

Beyond meta-analysis, Sharples used the Gibbs sampler to estimate posterior distributions of transition rates between grades of coronary disease after cardiac transplantation, modelled as a continuous-time Markov process, the kind of input a Markov model needs.

Choices that change the output

Starting values and warm-up length decide how much of the early chain reflects the start rather than the posterior, and blocking and parameterisation decide how fast it mixes. Priors on variance parameters need care in sparse networks. TSD 2 describes a vague gamma prior on the precision, such as Gamma(0.001, 0.001), as once popular but since fallen out of favour, and notes that it puts more weight on between-trial standard deviations near zero. Where sparse data make MCMC sampling get stuck at zero, it adds, a gamma prior may improve numerical stability.

How it differs from neighbouring methods

Metropolis-Hastings proposes a whole move and accepts or rejects it, while Gibbs moves one coordinate at a time from the exact conditional and always accepts. Hamiltonian Monte Carlo uses gradients and momentum to travel further per iteration, and Gelman and colleagues describe Stan as implementing it, so Stan is not a Gibbs sampler.

The Monte Carlo step of a probabilistic sensitivity analysis usually takes independent draws from distributions fitted to each input, whereas a Gibbs sampler produces dependent draws from a posterior it is still exploring and often supplies those inputs.

Limitations and common misreadings

Most problems with Gibbs output come from treating the draws as better than they are. The result is a decision model input with the wrong precision or the wrong correlation.

1. Treating draws as independent. Consecutive values are correlated, so at a posterior correlation of 0.99, 10,000 iterations carry the information of about 100 independent draws.
2. Assuming conditionals are always easy. Logistic regressions and many other models have conditionals with no standard form, so the sampler has to update them another way, for example with Metropolis steps.
3. Reading convergence from one chain. A chain stuck in one region can look stable, which is why TSD 2 asks for several chains from dispersed starts.
4. Confusing computation with model adequacy. A well-mixed chain describes the posterior of the model as specified, whether or not its likelihood and priors are right.

Sources

• Dias S, Welton NJ, Sutton AJ, Ades AE. NICE DSU Technical Support Document 2: A generalised linear modelling framework for pairwise and network meta-analysis of randomised controlled trials. Decision Support Unit; August 2011, last updated September 2016. Sections 2.1 and 6. https://www.sheffield.ac.uk/sites/default/files/2022-02/TSD2-General-meta-analysis-corrected-2Sep2016v2.pdf
• Gelman A, Carlin JB, Stern HS, Dunson DB, Vehtari A, Rubin DB. Bayesian Data Analysis. 3rd ed. CRC Press; 2013. Sections 2.5, 11.1 to 11.5 and 12.1, the chapter 11 bibliographic note and the chapter 12 introduction. http://www.stat.columbia.edu/~gelman/book/BDA3.pdf
• Geman S, Geman D. Stochastic relaxation, Gibbs distributions, and the Bayesian restoration of images. IEEE Transactions on Pattern Analysis and Machine Intelligence. 1984;PAMI-6(6):721-741. https://doi.org/10.1109/TPAMI.1984.4767596
• Plummer M. JAGS: a program for analysis of Bayesian graphical models using Gibbs sampling. In: Hornik K, Leisch F, Zeileis A, eds. Proceedings of the 3rd International Workshop on Distributed Statistical Computing (DSC 2003). Vienna; 2003. https://www.r-project.org/conferences/DSC-2003/Proceedings/Plummer.pdf
• Sharples LD. Use of the Gibbs sampler to estimate transition rates between grades of coronary disease following cardiac transplantation. Statistics in Medicine. 1993;12(12):1155-1169. https://doi.org/10.1002/sim.4780121205
• Spiegelhalter DJ. Incorporating Bayesian ideas into health-care evaluation. Statistical Science. 2004;19(1):156-174. https://doi.org/10.1214/088342304000000080