MCMC for Stochastic Epidemic Models

1 / 68

# MCMC for Stochastic Epidemic Models - PowerPoint PPT Presentation

MCMC for Stochastic Epidemic Models . Philip D. O’Neill School of Mathematical Sciences University of Nottingham. This includes joint work with…. Tom Britton (Stockholm) Niels Becker (ANU, Canberra) Gareth Roberts (Lancaster) Peter Marks (NHS, Derbyshire PCT). Contents.

I am the owner, or an agent authorized to act on behalf of the owner, of the copyrighted work described.

## PowerPoint Slideshow about 'MCMC for Stochastic Epidemic Models' - shelby

Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author.While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server.

- - - - - - - - - - - - - - - - - - - - - - - - - - E N D - - - - - - - - - - - - - - - - - - - - - - - - - -
Presentation Transcript
MCMC for Stochastic Epidemic Models

Philip D. O’Neill

School of Mathematical Sciences

University of Nottingham

This includes joint work with…
• Tom Britton (Stockholm)
• Niels Becker (ANU, Canberra)
• Gareth Roberts (Lancaster)
• Peter Marks (NHS, Derbyshire PCT)
Contents
• 1. MCMC: overview and basics
• 2. Example: Vaccine efficacy
• 3. Data augmentation
• 4. Example: SIR epidemic model
• 5. Model choice
• 6. Example: Norovirus outbreak
• 7. Other topics
Markov chain Monte Carlo (MCMC) Overview and basics
• The key problem is to explore a density function π known up to proportionality.
• The output of an MCMC algorithm is a sequence of samples from the correctly normalised π.
• These samples can be used to estimate summaries of π, e.g. its mean, variance.
How MCMC works
• Key idea is to construct a discrete time Markov chain X1, X2, X3, … on state space S whose stationary distribution is π.
• If P(dy,dx) is the transitional kernel of the chain this means that
How MCMC works (2)
• Subject to some technical conditions,

Distribution of Xn→ πas n →

• Thus to obtain samples from π we simulate the chain and sample from it after a “long time”.
Example: π (x)  x e-2x

X1, X2, …

XN = 1.2662

Example: π (x)  x e-2x

XN = 1.2662

XN+1 = 1.7840

Example: π (x)  x e-2x

XN = 1.2662

XN+1 = 1.7840

XN+2 = 0. 6629

Example: π (x)  x e-2x
• Suppose Markov chain output is

..., XN = 1.2662, XN+1 = 1.7840, XN+2 = 0.6629, …. ,XN+M = 1.0312

(i.e. discard initial N values, burn-in)

How to build the Markov chain
• Surprisingly, there are many ways to construct a Markov chain with stationary distribution π.
• Perhaps the simplest is the Metropolis-Hastings algorithm.
Metropolis-Hastings algorithm
• Set an initial value X1.
• If the chain is currently at Xn = x, randomly propose a new position Xn+1 = y according to a proposal density q(y | x).
• Accept the proposed jump with probability
• If not accepted, Xn+1 = x.
Why the M-H algorithm works
• Let P(dx,dy) denote the transition kernel of the chain.
• Then P(dx,dy) is approximately the probability that the chain jumps from a region dx to a region dy.
• We can calculate P(dx,dy) as follows:
Why M-H works (3)

This last equation shows that πis a stationary distribution for the Markov chain.

• The choice of proposal q(y|x) is fairly arbitrary.
• Popular choices include

q(y|x) = q(y) (Independence sampler)

q(y|x) ~ N(x, 2) (Gaussian proposal)

q(y|x) = q(|y-x|) (Symmetric proposal)

• In practice, MCMC is almost always used for multi-dimensional problems. Given a target density π(x1, x2, … ,xn) it is possible to update each component separately, or even in blocks, using different M-H schemes.
• A popular multi-dimensional scheme is the Gibbs sampler, in which the proposal for a component xi is its full conditional density

π (xi | (x1,…,xi-1, xi+1, … ,xn) )

• The M-H acceptance probability is equal to one in this case.
• How to check convergence? There is no guaranteed way.

Visual inspection of trace plots; diagnostic tools (e.g. looking at autocorrelation).

• Starting values – try a range
• Acceptance rates – not too large/small
• Mixing – how fast does the chain move around?
Contents
• 1. MCMC: overview and basics
• 2. Example: Vaccine efficacy
• 3. Data augmentation
• 4. Example: SIR epidemic model
• 5. Model choice
• 6. Example: Norovirus outbreak
• 7. Other topics
2. Example: Vaccine Efficacy
• Outbreak of Variola Minor, Brazil 1956
• Data on cases in households (size 1 to 12)
• 338 households: 126 had no cases
• 1542 individuals:

809 vaccinated, 85 cases

733 unvaccinated, 425 cases

Objective: estimate vaccine efficacy

Disease transmission model
• Population divided into separate households.
• Divide transmission into community-acquired and within-household.
• q = P( individual avoids outside infection )
•  = P ( one individual fails to infect another

in the same household )

Vaccine response model
• For a vaccinated individual, three responses can occur: complete protection; vaccine failure; or partial protection and infectivity reduction.
• c = P(complete protection)
• f = P(vaccine failure)
• a = proportionate susceptibility reduction
• b = proportionate infectivity reduction
Vaccine response : (A,B)
• A convenient way of summarising the random response is to suppose that an individual’s susceptibility and infectivity reduction is given by a bivariate random variable (A,B). Thus

P[ (A,B) = (0, -)] = c

P[ (A,B) = (1,1) ] = f

P[ (A,B) = (a,b) ] = 1-c-f

Efficacy Measures
• Furthermore it is sensible to define measures of vaccine efficacy using (A,B).
• VES = 1- E[A]

= 1 - f - a(1-f-c) is a protective measure

• VEI = 1 - E[AB] / E[A]

= 1 - [f + ab(1-f-c)] / [f + a(1-f-c)]

is a measure of infectivity reduction

• Note both are functions of basic model parameters
Bayesian inference
• Object of inference is the posterior density

(  | n ) = ( a,b,c,f,q, | n )

where n is the data set. By Bayes’ Theorem

(| n )  (n | )  (),

where (| n ) is the likelihood,

and  () is the prior density for .

MCMC details
• There are six parameters: a,b,c,f,q,
• Each parameter has range [0,1]
• Update each parameter separately using a Metropolis-Hastings step with Gaussian proposal centered on the current value
MCMC pseudocode
• Initialise parameters (e.g. a = 0.5, b=0.5,…)
• User input burn-in (B), sample size (S), thinning gap (T)
• LOOP: counter from –B to (S x T)

Update a, update b, …, update 

IF (counter > 0) AND (counter/T is integer)

THEN store current values

• END LOOP
Updating details for a
• Propose ã~ N(a, 2)
• Accept with probability
• Note that the (symmetric) proposal cancels out
• The other parameters are updated similarly
Results for VES

Posterior mean: VES = 1 – E(A) = 0.84

Posterior Standard Deviation = 0.03

These results are easily obtained using the raw

Markov chain output for the model parameters.

Contents
• 1. MCMC: overview and basics
• 2. Example: Vaccine efficacy
• 3. Data augmentation
• 4. Example: SIR epidemic model
• 5. Model choice
• 6. Example: Norovirus outbreak
• 7. Other topics
4. Data augmentation
• Suppose we have a model with unknown parameter vector  = (1,2,…,n).
• Available data are y = ( y1, y2,…, ym ).
• If the likelihood π(y |  ) is intractable…
• …one solution is to introduce extra parameters (“missing data”)

x = (x1, x2,…, xp)

such that π(y , x |  ) is tractable.

Data augmentation (2)
• The extra parameters x = (x1, x2,…, xp) are simply treated as unknown model parameters as before.
• To obtain samples from π( y |  ), take samples from π(y , x |  ) and ignore x.
• Such a scheme is often easy using MCMC.
Data augmentation (3)
• Can also add parameters to improve the mixing of the Markov chain (auxiliary variables).
• Choosing how to augment data is not always obvious!
Contents
• 1. MCMC: overview and basics
• 2. Example: Vaccine efficacy
• 3. Data augmentation
• 4. Example: SIR epidemic model
• 5. Model choice
• 6. Example: Norovirus outbreak
• 7. Other topics
4. SIR Epidemic Model
• Suppose we observe daily numbers of cases during an epidemic outbreak in some fixed population.
• Objective is to say something about infection rates and infectious period duration of the disease.
Model definition
• Population of N individuals
• At time t there are:

St susceptibles

It infectives

Rt recovered/immune individuals

Thus St + It + Rt = N for all t.

Initially (S0, I0 ,R0 ) = (N-1,1,0).

Model definition (2)
• Each infectious individual remains so for a length of time TI~ Exponential().
• During this time, infectious contacts occur with each susceptible according to a Poisson process of rate /N.
• Thus overall infection rate is  St It/N.
• Two model parameters,  and .
Data, likelihood, augmentation
• Suppose we observe removals at times 0 ≤ r1≤ r2 ≤ …≤ rn ≤ .
• Define r = ( r1,r2 , …,rn ).
• The likelihood of the data, π (r | , ), is practically intractable.
• However, given the (unknown) infection times i = ( i1,i2 , …,in ), π (i ,r | , ) is tractable.
MCMC algorithm
• Specifically,
• It follows that if π() ~ Gamma distribution thenπ( | …) ~ Gamma distribution also.
• Same is true for .
• So can update  and  using a Gibbs step.
MCMC algorithm – infection times
• It remains to update the infection times i = ( i1,i2 , …,in )
• Various ways of doing this.
• A simple way is to use a M-H scheme to randomly move the times.
• For example, propose a new ik by picking a new time uniformly at random in (0,).
Updating infection times

Updating I2 :

Acceptance prob. π (i*,r | ,) / π (i,r | ,)

I6

I2

I4

I6

I4

I2*

Extensions
• Epidemic not known to be finished by 
• Non-exponential infectious periods
• Multi-group models (e.g. age-stratified data)
• More sophisticated updates of infection times
• Inclusion of latent periods
Contents
• 1. MCMC: overview and basics
• 2. Example: Vaccine efficacy
• 3. Data augmentation
• 4. Example: SIR epidemic model
• 5. Model choice
• 6. Example: Norovirus outbreak
• 7. Other topics
5. Model Choice
• Bayesian model choice problems can also be implemented using MCMC.
• So-called “transdimensional MCMC” (alias “Reversible Jump MCMC”) is used.
• The basic idea is to construct the Markov chain on the union of the different sample spaces and (essentially) use M-H.
Simple example
• Model 1 has two parameters: , 
• Model 2 has one parameter: 
• The Markov chain moves between models and within models
• E.g. Xn = (1, , , ) for model 1, ignore 
• Practical question – how to jump between models?
Contents
• 1. MCMC: overview and basics
• 2. Example: Vaccine efficacy
• 3. Data augmentation
• 4. Example: SIR epidemic model
• 5. Model choice
• 6. Example: Norovirus outbreak
• 7. Other topics
6. Example: Norovirus outbreak

Outbreak of gastroenteritis in summer 2001 at school in Derbyshire, England.

A single strain of Norovirus virus found to be the causative agent.

Outbreak data
• 15 classrooms, each child based in one.
• Absence records plus questionnaires.
• 492 children of whom 186 were cases.
• Data include age, period of illness, times of vomiting episodes in classrooms.
Question of interest
• Does vomiting play a significant role in transmission?
• Total of 15 vomiting episodes in classrooms.
Stochastic transmission model
• Assumption:
• A susceptible on weekday t remains so on day t+1 if they avoid infection from each infective child;
• per-infective daily avoidance probabilities are
• classmate : qc
• schoolmate: qs
• in class, vomiters : qv
Transmission model (2)
• At weekends, a susceptible remains so by avoiding infection from all infectives in the community,
• per-infective avoidance probability is q.
Two models
• M1 : Full model: qc, qv, qs, q
• Vomiters treated separately
• M2 : Sub-model: qc, qv=qc, qs, q
• Vomiters classed as normal infectives
MCMC algorithm
• Construct Markov chain on state space

S = { ( qc, qs, qv, q, M) }

where M = 1 or 2 is the current model

• Model-switching step to update M
• Random walk updates for the q’s
Between-model jumps
• Full model  sub-model:

Propose qv = qc

• Sub-model  full model:

Propose qv = qc + N(0, 2 )

• Acceptance probabilities straightforward
Results
• Uniform(0,1) q. priors; P(M1) = 0.5
• P(M1 | data) = 0.0001 (Full model)
• P(M2 | data) = 0.999 (Sub model)
Contents
• 1. MCMC: overview and basics
• 2. Example: Vaccine efficacy
• 3. Data augmentation
• 4. Example: SIR epidemic model
• 5. Model choice
• 6. Example: Norovirus outbreak
• 7. Other topics
7. Other topics
• 1. Improving algorithm performance
• 2. Perfect simulation
• 3. Some conclusions
Improving algorithm performance
• Choose parameters to reduce correlation
• Trade-off between ease of computation and mixing behaviour of chain
• Choice of M-H proposal distributions
• Choice of blocking schemes
Perfect simulation
• Detecting convergence can be a real problem in practice.
• Perfect simulation is a method for constructing a chain that is known to have converged by a certain time.
• Unfortunately it is far less applicable than MCMC.
Some conclusions
• MCMC methods are hugely powerful
• The methods enable analysis of very complicated models
• Sample-based methods easily permit exploration of both parameters and functions of parameters
• Implementation is often relatively easy
• Software available (e.g. BUGS)