Sequential Monte Carlo

Sequential Monte Carlo in Health Economic Modelling

Sequential Monte Carlo (SMC) is a family of simulation algorithms that moves a population of weighted values, called particles, through a chain of probability distributions, reweighting, resampling and sometimes moving the particles at each step. In health economics the chain can follow surveillance data week by week, as a particle filter does for the hidden state of an epidemic, or it can lead in small steps towards the posterior of a disease model's unknown parameters, which is how SMC calibrates such models. The output is a set of weighted parameter or state values that can be run through a cost-effectiveness model in the same way as draws in a probabilistic sensitivity analysis. The main practical weakness is weight degeneracy, where a few particles carry nearly all the weight.

One algorithm for any sequence of target distributions

Doucet and Johansen, in their tutorial on particle filtering and smoothing, note that particle filtering and SMC are often presented as the same thing, while SMC covers a broader range of algorithms. Their generic algorithm gives weighted samples from any sequence of distributions, each known only up to a constant.

$$\pi_n(x) = \frac{\gamma_n(x)}{Z_n}, \qquad Z_n = \int \gamma_n(x)\, dx$$

where $\pi_n$ is the target distribution at step $n$, $\gamma_n$ is an unnormalised version of it that can be evaluated at any point $x$, and $Z_n$ is the normalising constant, which may be unknown. In filtering, $\gamma_n$ is the joint density of the hidden states and the observations so far, and $Z_n$ is the marginal likelihood of those observations.

Sequential importance sampling

The particles are built up by importance sampling one step at a time. Each particle is extended by a draw from a proposal distribution and its weight is updated by an incremental factor, so the cost of each step does not grow with the length of the sequence.

$$w_n(x_{1:n}) = w_{n-1}(x_{1:n-1})\, \alpha_n(x_{1:n}), \qquad \alpha_n(x_{1:n}) = \frac{\gamma_n(x_{1:n})}{\gamma_{n-1}(x_{1:n-1})\, q_n(x_n \mid x_{1:n-1})}$$

where $w_n$ is the unnormalised weight of a particle path $x_{1:n}$ at step $n$, $\alpha_n$ is the incremental importance weight, and $q_n(x_n \mid x_{1:n-1})$ is the proposal density used to extend the path. Doucet and Johansen show that, uncorrected, the variance of such estimates typically grows exponentially with $n$, even with a good proposal.

Resampling and the effective sample size

Resampling is the partial remedy. Particles with low weights are removed with high probability and those with high weights are copied, after which all carry equal weight. Because resampling adds noise of its own, Doucet and Johansen recommend resampling only when the weights have become too uneven, judged by the effective sample size.

$$W_n^{(i)} = \frac{w_n^{(i)}}{\sum_{j=1}^{N} w_n^{(j)}}, \qquad \text{ESS}_n = \left( \sum_{i=1}^{N} \left( W_n^{(i)} \right)^2 \right)^{-1}$$

where $N$ is the number of particles, $w_n^{(i)}$ and $W_n^{(i)}$ are the unnormalised and normalised weights of particle $i$, and $\text{ESS}_n$ lies between 1 and $N$. They read ESS as the number of perfect samples giving about the same variance, give $N/2$ as the typical resampling threshold and name systematic resampling as the most widely used scheme.

Three sequences that matter for health models

The generic algorithm becomes a particular method once the sequence of targets is chosen. Three choices are common in health modelling.

The first follows data in time. With targets $p(x_{1:n} \mid y_{1:n})$ for hidden states $x$ and observations $y$ in a state space model, SMC becomes the particle filter, which estimates quantities such as the number of people currently infectious as each week of case or admission data arrives. The Kalman filter solves the same problem exactly for linear Gaussian models, and the particle filter page gives the bootstrap filter.

The second moves a fixed set of model parameters through artificial intermediate distributions. Del Moral, Doucet and Jasra call these SMC samplers, and in their framework the step index is just a counter. One sequence they use is a tempered posterior.

$$\pi_n(\theta) \propto L(y \mid \theta)^{\phi_n}\, f(\theta), \qquad 0 \le \phi_1 < \dots < \phi_p = 1$$

where $\theta$ is the vector of model parameters, $L(y \mid \theta)$ is the likelihood of the calibration data $y$, $f(\theta)$ is the prior density, and $\phi_n$ is the tempering exponent at step $n$. When each step moves the particles with a Markov chain Monte Carlo kernel that leaves $\pi_n$ unchanged, and the matching reversed kernel is used, their incremental weight reduces to $\gamma_n(\theta)/\gamma_{n-1}(\theta)$ evaluated at the particle's previous value, which for this sequence is $L(y \mid \theta)^{\phi_n - \phi_{n-1}}$.

The third is used when the likelihood cannot be computed, as with many stochastic transmission and agent-based models. Toni and colleagues describe approximate Bayesian computation by SMC, in which particles sampled from the prior pass through distributions defined by tolerances $\epsilon_1 > \dots > \epsilon_T \ge 0$ on the distance between simulated and observed data, so that the population approaches the approximate posterior. They applied it to dynamical systems including SIR epidemic models, and to choosing between candidate models.

Calibrating health policy models with weighted samples

Menzies and colleagues, in a tutorial on Bayesian calibration of health policy models, start with sampling importance resampling, in which parameter sets drawn from the prior are run through the model, weighted by their likelihood against the calibration targets and resampled. This is SMC with only one step, which reaches the posterior at once. In their worked example 100,000 prior draws resampled to 100,000 contained only 797 unique parameter sets, with an ESS of 88, which they took as a sign of non-trivial Monte Carlo error. They warn that the method is inefficient when the prior is dispersed relative to the likelihood. Incremental mixture importance sampling, a method developed for mechanistic models, was 150 to 200 times more efficient on the same problem.

Tempered SMC tackles the same weakness with many small steps, resampling when ESS falls and moving particles with MCMC between steps. Whatever the method, Menzies and colleagues then run the calibrated parameter sets through the model as a conventional Monte Carlo simulation, adding draws from the prior for parameters left out of the calibration, such as the treatment cost in their example. In their example, a dynamic transmission model of a hypothetical sexually transmitted disease, this is how calibrated epidemiology reaches the cost-effectiveness of earlier treatment.

SMC inside MCMC for states and parameters together

Doucet and Johansen note that treating parameters as states with no dynamics gives a degenerate transition, which prevents the standard algorithms from working in practice. Andrieu, Doucet and Holenstein proposed particle MCMC, which uses SMC to build efficient high-dimensional proposals inside an MCMC sampler. In the particle marginal Metropolis-Hastings version, the exact marginal likelihood in each acceptance ratio is replaced by its SMC estimate. Because that estimate is unbiased, the sampler still targets the correct posterior for the parameters. Endo, van Leeuwen and Baguelin introduce the method for disease dynamics modellers, with infectious disease examples in R.

Worked example of a tempered calibration step

All numbers are illustrative. A natural history model needs the annual probability $p$ that an untreated precancerous lesion progresses, and the calibration target is a cohort in which 22 of 100 lesions progressed in a year. The prior is represented by five equally weighted particles, $p$ = 0.10, 0.20, 0.30, 0.40 and 0.50, each with weight 0.2, so the prior mean is 0.30.

1. Likelihood of each particle. The binomial likelihood is proportional to $p^{22}(1-p)^{78}$, since the binomial coefficient cancels when weights are normalised. Relative to the best particle, the likelihoods are 0.0023, 1, 0.2242, 0.0008 and under 0.0001.

2. One jump to the posterior. Weighting by the full likelihood, as in sampling importance resampling, gives normalised weights of 0.0019, 0.8148, 0.1827, 0.0006 and 0.0000. The squared weights sum to about 0.6973, so the ESS is 1 / 0.6973 = 1.43, below the threshold of 5 / 2 = 2.5.

3. A tempered first step. With $\phi_1$ = 0.2 the incremental weights are the relative likelihoods raised to 0.2, namely 0.2975, 1, 0.7415, 0.2374 and 0.0369. These sum to 2.3133 and normalise to 0.1286, 0.4323, 0.3206, 0.1026 and 0.0159. The squared weights sum to 0.3170, so the ESS is 1 / 0.3170 = 3.15 and no resampling is needed yet.

4. Why the moves matter. If the run continued to $\phi$ = 1 with no resampling and no moves, the second incremental weight would be the relative likelihood raised to 0.8 and the final weights would equal those in step 2 exactly. Tempering gains only through what happens between steps. If the particles were resampled after step 3, multinomial resampling would give expected copies of 5 times each weight, or 0.64, 2.16, 1.60, 0.51 and 0.08, and an MCMC move then spreads the copies of 0.20 and 0.30 into new nearby values.

5. Carry the result into the model. With the step 2 weights the posterior mean of $p$ is 0.0019 × 0.10 + 0.8148 × 0.20 + 0.1827 × 0.30 + 0.0006 × 0.40 = 0.2182, against the prior mean of 0.30. Per 1,000 people with an untreated lesion that is about 218 progressions instead of 300, and at an illustrative GBP 15,000 per progression the expected cost is 0.2182 × 1,000 × 15,000 = GBP 3,273,000, against GBP 4,500,000 under the prior.

6. Interpretation. Under a continuous uniform prior on 0.05 to 0.55, the range the five particles discretise, the exact posterior is Beta(23, 79) truncated to that range, with mean 0.2255, the same as 23 / 102 to four decimals. The gap between 0.2182 and 0.2255 comes from using five fixed particles, and the ESS of 1.43 warns that the estimate rests almost entirely on one particle.

Design choices that change the answer

Each of the analyst's choices affects accuracy or run time. Doucet and Johansen show that the weight variance is smallest when the proposal is close to the target. The number of particles sets the Monte Carlo error of every weighted estimate. The tempering schedule $\phi_1, \dots, \phi_p$ sets the size of each step, and the MCMC move sets how quickly copies of one particle become distinct again. In ABC SMC the distance function and final tolerance also define which approximate posterior is targeted.

What SMC output supports in an economic evaluation

A weighted particle set approximates a posterior, so it supports the same uses as posterior draws from MCMC. Means, intervals and the probability of cost-effectiveness are computed with the weights, and Doucet and Johansen prefer estimates taken before resampling, because resampling adds variance.

Doucet and Johansen show that the weights also estimate $Z_n / Z_{n-1}$ at each step, and the product of these ratios estimates the evidence for the model. That supports Bayesian model selection and so speaks to structural uncertainty, as in the use Toni and colleagues made of ABC SMC for candidate epidemic models. No SMC output corrects a misspecified model, and Menzies and colleagues warn that tweaking priors, likelihood or model solely to maximise fit risks overfitting.

Boundaries with nearby methods

SMC and Monte Carlo simulation in probabilistic sensitivity analysis both draw random samples, but PSA samples inputs from distributions already specified, while SMC finds those distributions from data. MCMC builds one dependent chain aimed at a single posterior, whereas SMC runs a population through a sequence of targets, and Del Moral, Doucet and Jasra describe MCMC kernels as frequent ingredients of SMC samplers. Sequential analysis of interim trial data is a different idea that shares only the word.

Where SMC fails and how it is misread

Most failures of SMC in practice come from degeneracy, which Doucet and Johansen call the factor that most often determines whether an algorithm works. Other errors come from expecting more of the output than it can give.

1. Reading the number of particles as the sample size. The ESS, not $N$, measures how much the weighted sample says, and the Menzies example kept 797 unique sets out of 100,000.
2. Trusting path estimates over long series. Doucet and Johansen show that repeated resampling reduces the number of distinct values representing early states, so methods relying on full paths fail over long enough series.
3. Expecting it to scale without limit. They note that problems of high enough dimension and complexity cannot be characterised accurately in reasonable time.
4. Confusing the acronyms. Doucet and Johansen record that SMC with resampling is sometimes called sequential importance resampling (SIR), Menzies and colleagues use SIR for sampling importance resampling, and SIR is also the susceptible, infected, recovered epidemic model to which Toni and colleagues applied ABC SMC.

Sources

• Andrieu C, Doucet A, Holenstein R. Particle Markov chain Monte Carlo methods. Journal of the Royal Statistical Society Series B. 2010;72(3):269-342. https://doi.org/10.1111/j.1467-9868.2009.00736.x
• Del Moral P, Doucet A, Jasra A. Sequential Monte Carlo samplers. Journal of the Royal Statistical Society Series B. 2006;68(3):411-436. https://doi.org/10.1111/j.1467-9868.2006.00553.x
• Doucet A, Johansen AM. A tutorial on particle filtering and smoothing: fifteen years later. Version 1.1, December 2008, with typographical corrections March 2012. https://www.stats.ox.ac.uk/~doucet/doucet_johansen_tutorialPF2011.pdf
• Endo A, van Leeuwen E, Baguelin M. Introduction to particle Markov-chain Monte Carlo for disease dynamics modellers. Epidemics. 2019;29:100363. https://doi.org/10.1016/j.epidem.2019.100363
• 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
• Toni T, Welch D, Strelkowa N, Ipsen A, Stumpf MPH. Approximate Bayesian computation scheme for parameter inference and model selection in dynamical systems. Journal of the Royal Society Interface. 2009;6(31):187-202. https://doi.org/10.1098/rsif.2008.0172