Metropolis-Hastings

Metropolis-Hastings Algorithm in Bayesian Health Economic Models

The Metropolis-Hastings algorithm is a general family of Markov chain Monte Carlo samplers for fitting Bayesian models in health economics, with Gibbs sampling as a special case. At each step it suggests a new value for the unknown parameters and decides by a probability rule whether the chain moves there or stays put, so that over many steps the chain spends time in each region in proportion to its posterior probability, so the draws trace out the posterior distribution. It matters in health technology assessment because evidence syntheses and calibrated disease models often have posteriors with no standard form, and its draws can be passed directly into a probabilistic sensitivity analysis. The output is a correlated sequence, so the iteration count overstates its information, and a badly tuned proposal can waste a long run.

How Hastings generalised the Metropolis algorithm

Gelman and colleagues, in Bayesian Data Analysis, credit Metropolis and colleagues in 1953 and Metropolis and Ulam in 1949 as apparently the first to describe Markov chain simulation of probability distributions, and Hastings with generalising their algorithm in 1970. The original needed a symmetric proposal, and the generalisation allows asymmetric ones with a correction in the acceptance ratio.

Hastings set out the property that makes the method so useful in Bayesian analysis. The computations depend on the target density only through ratios of its values at two points, so its normalising constant need not be known. By Bayes' theorem a posterior is the likelihood times the prior distribution divided by the probability of the data, and that hard-to-compute divisor cancels in a ratio. He also warned that samples from a Markov chain are correlated, so their error needs more care to assess.

The acceptance rule

At iteration $t$ a candidate $\theta^{*}$ is drawn from a proposal distribution, which Gelman and colleagues call the jumping distribution. It is accepted with the following probability.

$$\alpha = \min\left\{1, \frac{p(\theta^{*} \mid y)\, J(\theta^{(t-1)} \mid \theta^{*})}{p(\theta^{(t-1)} \mid y)\, J(\theta^{*} \mid \theta^{(t-1)})}\right\}$$

where $\theta^{(t-1)}$ is the current value of the parameter vector, $\theta^{*}$ is the candidate, $y$ is the observed data, $p(\cdot \mid y)$ is the posterior density, which may be unnormalised, $J(a \mid b)$ is the density of proposing $a$ when the chain is at $b$, and $\alpha$ is the probability of moving to the candidate. A uniform random number is drawn and the move is made when it falls below $\alpha$.

When the proposal is symmetric the two $J$ terms cancel. As Chib and Greenberg put it, a move uphill to a higher density is always accepted, and a move downhill is accepted with probability equal to the density ratio. If the candidate is rejected the current value is recorded again, and Gelman and colleagues stress that this still counts as an iteration. Deleting the repeats would under-represent the regions where the chain lingers because they are most probable.

Chib and Greenberg derive the rule from reversibility, also called detailed balance. For a chain already sampling the target, the flow of probability between any two points must be equal in both directions, and the acceptance step scales down moves in whichever direction would otherwise carry too much flow. This makes the posterior the stationary distribution, which the chain approaches under mild conditions, chiefly that every region can eventually be visited. Those conditions say nothing about speed, which is why convergence is checked empirically with several chains.

Choosing the proposal distribution

The proposal governs efficiency. Chib and Greenberg describe five families of proposal, including these two. In a random walk chain the candidate is the current value plus noise, often normal. In an independence chain the candidate comes from a fixed distribution, usually an approximation to the posterior, whose tails should dominate those of the target, a requirement they compare with the one on the sampling density in importance sampling.

The scale of a random walk trades two failures against each other. Small steps are nearly always accepted but crawl, and large steps mostly land in the tails and are rejected, and Chib and Greenberg note that both show up as high autocorrelation. For a roughly normal posterior Gelman and colleagues give a rule of thumb.

$$J(\theta^{*} \mid \theta^{(t-1)}) = \mathrm{N}\left(\theta^{*} \mid \theta^{(t-1)}, c^{2}\Sigma\right), \qquad c \approx \frac{2.4}{\sqrt{d}}$$

where $\Sigma$ is the posterior variance matrix or an estimate of it, $d$ is the number of parameters updated together, and $c$ is the scale factor. The acceptance rate is then about 0.44 in one dimension, falling to about 0.23 in high dimensions, and efficiency relative to independent sampling is about 0.3 divided by $d$.

Gelman and colleagues warn that a proposal which keeps adapting as the chain runs can converge to the wrong distribution, so they tune in an adaptive phase and use only draws from a later fixed phase for inference.

Worked example with a readmission probability

All numbers in this example are illustrative. A decision model for a new day-surgery pathway needs the probability $p$ of emergency readmission within 30 days, and a pilot series recorded 6 readmissions among 40 patients. The analysis works on the log odds $\theta$, the logarithm of $p/(1-p)$, with the binomial likelihood and logit link that NICE Decision Support Unit Technical Support Document 2 (TSD 2) describes for event data, and a normal prior on $\theta$ with mean minus 2.0 and standard deviation 1.0, a prior median for $p$ of 0.119. With this non-conjugate prior the posterior has no standard form, though its log density is easy to compute.

$$\ell(\theta) = y\,\theta - n \log\left(1 + e^{\theta}\right) - \frac{(\theta - \mu_0)^{2}}{2\sigma_0^{2}}$$

where $\ell(\theta)$ is the natural log of the unnormalised posterior density, $y = 6$ is the number of readmissions, $n = 40$ is the number of patients, $\mu_0 = -2.0$ is the prior mean and $\sigma_0 = 1.0$ is the prior standard deviation.

1. Score the current value. The chain is at $\theta$ = minus 1.80, a probability of 0.142. The first term is 6 × minus 1.80 = minus 10.8000, the second is 40 times the log of one plus $e^{-1.80}$, or 40 times 0.15298, which is 6.1191, and the prior term is 0.20 × 0.20 / 2 = 0.0200. So $\ell$ = minus 10.8000 minus 6.1191 minus 0.0200 = minus 16.9391.

2. A random walk move that is accepted. The proposal is normal with standard deviation 0.5, and a standard normal draw of 0.60 gives a candidate of minus 1.80 + 0.5 × 0.60 = minus 1.50, a probability of 0.182. Its terms are minus 9.0000, then 40 × 0.20141 = 8.0565, then 0.50 × 0.50 / 2 = 0.1250, so $\ell$ = minus 17.1815. The log ratio is minus 17.1815 + 16.9391 = minus 0.2424, a ratio of 0.785, with no correction because the proposal is symmetric. A uniform draw of 0.42 is below 0.785, so the chain moves to minus 1.50.

3. A move that is rejected. A standard normal draw of minus 1.40 gives minus 1.50 minus 0.70 = minus 2.20, with $\ell$ = minus 13.2000 minus 4.2033 minus 0.0200 = minus 17.4233. The log ratio is minus 17.4233 + 17.1815 = minus 0.2418, a ratio of 0.785. A uniform draw of 0.91 exceeds it, so minus 1.50 is recorded a second time.

4. The Hastings correction. A normal approximation at the mode, minus 1.78 with standard deviation 0.41, suggests an independence proposal with mean minus 1.75 and standard deviation 0.60, wider than the approximation so that its tails dominate the target's. From minus 1.50 it offers minus 2.10, with $\ell$ = minus 12.6000 minus 4.6208 minus 0.0050 = minus 17.2258, and a posterior ratio of $e^{-0.0443}$ = 0.9567. The current value lies 0.417 proposal standard deviations above the proposal mean and the candidate 0.583 below, so the correction is $e^{-0.5 \times (0.25/0.60)^2 + 0.5 \times (0.35/0.60)^2}$ = $e^{0.0833}$ = 1.0869. The full ratio is 0.9567 × 1.0869 = 1.040, so the move is accepted outright. The candidate is slightly less probable, but the proposal offers it less often than the current value, and the correction compensates.

5. Tune the scale and run. Numerical integration gives the exact posterior mean of $p$ as 0.1457, with a 95% credible interval of 0.063 to 0.258, and a standard deviation of 0.418 for $\theta$, so the rule of thumb suggests a scale of 2.4 × 0.418 = 1.00. One run at each of three scales kept 20,000 draws after 2,000 warm-up iterations, with effective sample size estimated from the autocorrelations.

The tuned scale accepted 43.9% of moves, close to the 0.44 expected in one dimension, and gave about 70 times the effective draws of the smallest scale. Another run would give slightly different decimals because of Monte Carlo error.

6. Pass the draws to the model. At an illustrative GBP 2,500 per readmission, the expected readmission cost per patient is 2,500 × 0.1457 = GBP 364, with a 95% interval of 2,500 × 0.063 = GBP 158 to 2,500 × 0.258 = GBP 645. The mean sits between the pilot proportion of 0.150 and the prior median of 0.119, nearer the data, which carry a log odds precision of about 40 × 0.15 × 0.85 = 5.1 against the prior's 1.0. A probabilistic analysis would run each retained draw through the model.

Where it is used in health economic evaluation

Gelman and colleagues call the Gibbs sampler their first choice for conditionally conjugate models and use the Metropolis algorithm for models that are not, with logistic regression as their example. They also mix the two, with one-dimensional Metropolis updates where a full conditional is not standard. Plummer, introducing JAGS in 2003, described its samplers as Gibbs samplers acting on single nodes with an Adaptive Rejection Metropolis sampler as the workhorse, so Metropolis-type steps sat inside the Gibbs cycle without the analyst writing any. The binomial logit models of a network meta-analysis are of the non-conjugate kind.

Calibration is another use. Whyte, Walsh and Chilcott used the Metropolis-Hastings algorithm to estimate the natural history parameters of a colorectal cancer state transition model, fitting cancer incidence, autopsy data on polyps and detection rates from faecal occult blood screening in England. Coded in Visual Basic, it ran overnight on a standard desktop PC, and the posterior parameter sets were used in probabilistic sensitivity analysis. Menzies and colleagues, in a calibration tutorial, note that MCMC can run in WinBUGS or Stan for relatively simple models, while an intractable likelihood, common in microsimulation, can be handled by approximate Bayesian computation instead.

How it differs from neighbouring methods

Gibbs sampling is a special case in which each proposal is the full conditional of one parameter, and Gelman and colleagues show that the acceptance ratio is then exactly 1. Metropolis-Hastings needs no conditional in closed form, at the price of a proposal to tune. They describe Hamiltonian Monte Carlo as a generalisation of the Metropolis algorithm with momentum variables, which moves further per iteration and mixes faster in high dimensions.

Importance sampling and plain Monte Carlo integration use independent draws, not a dependent sequence. Gelman and colleagues describe the Metropolis rule as a stochastic version of stepwise mode-finding, and Chib and Greenberg note that it underlies simulated annealing, but the sampler aims at the whole distribution, not its peak.

Limitations and common misreadings

Many problems with Metropolis-Hastings output come from trusting the iteration count. The posterior is only as good as the model.

1. Reading a high acceptance rate as good mixing. In the example the chain that accepted 95.7% of moves was the least efficient, and Chib and Greenberg caution that a chain at the optimal rate can still be highly autocorrelated.
2. Discarding rejected iterations. A repeated value is a valid draw.
3. Tuning throughout the run. Draws from an adaptive phase are not in general draws from the target.
4. Treating draws as independent. Monte Carlo error depends on the effective number of draws, which can be a small fraction of the total.
5. Confusing computation with model adequacy. A well-mixed chain describes the posterior of the model as specified, right or wrong.

Sources

• Chib S, Greenberg E. Understanding the Metropolis-Hastings algorithm. The American Statistician. 1995;49(4):327-335. https://doi.org/10.1080/00031305.1995.10476177
• 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. 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. Chapters 11 and 12. http://www.stat.columbia.edu/~gelman/book/BDA3.pdf
• Hastings WK. Monte Carlo sampling methods using Markov chains and their applications. Biometrika. 1970;57(1):97-109. https://doi.org/10.1093/biomet/57.1.97
• Menzies NA, Soeteman DI, Pandya A, Kim JJ. Bayesian methods for calibrating health policy models: a tutorial. PharmacoEconomics. 2017;35(6):613-624. https://doi.org/10.1007/s40273-017-0494-4
• 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
• Whyte S, Walsh C, Chilcott J. Bayesian calibration of a natural history model with application to a population model for colorectal cancer. Medical Decision Making. 2011;31(4):625-641. https://doi.org/10.1177/0272989X10384738