dynamic compartmental models - ex. SIR models

[jason andrews]

Dear PyMCers,

I wanted to see if people had more examples of dynamic compartmental
models, such as difference equation models, monte carlo simulations
using tau-leap, etc. to estimate transition parameters.

I'm new to PyMC and was trying to do a simple SIR model to explore
it. SIR models (Susceptible-Infected-Recovered) are simple,
compartmental infectious disease models. They start with the
population being distributed in one of the three states (SIR), and
then susceptible individuals can move to the infected state at a rate
(beta) times the product of the population of infected and susceptible
individuals. Individuals then recover at a rate gamma. There are a
lot of great examples along with Python code here:
http://www.warwick.ac.uk/~masfz/ModelingInfectiousDiseases/Chapter2/Program_2.1/index.html

I want to take observed data about the population in those three
states over time to estimate the transmission parameter (beta) and
recovery parameter (gamma)

The simple SIR model can be
S[t+1] = S[t] - beta*I[t]*S[t]/N[t]
I[t+1]= I[t] + beta*I[t]*S[t]/N[t] - gamma*I[t]
R[t+1] = R[t] + gamma*I[t]

where N[t]=S[t]+I[t]+R[t]
(there are several ways to write the SIR model depending on
assumptions about transmission).

I'm a beginner to this and haven't had luck implementing it. Below is
my attempt at an SI model (without R) with made-up data and
parameters. Ideally would incorporate stochasticity at each time step,
below the 'actual data' (SIR populations) was deterministic for
simplicity.

I would be grateful for any examples, comments, or suggestions for
where to learn more. I've looked at the Fisheries Surplus model on
the PyMC site.

thank you,
Jason

##SI model
from pymc import*
import numpy as np
from numpy import *
#observed data
susceptible_array= np.array([999,997,996,994,993,992,990,989,986,984])
infected_array= np.array([1,2,5,6,7,8,9,11,13,15])

S = np.zeros(len(susceptible_array))
I = np.zeros(len(infected_array))
S[0] = 999
I[0] = 1

beta = Normal('beta', mu=0, tau=5)
gamma = Normal('gamma', mu=0, tau=5)
for i in range(1,len(susceptible_array)):
S[i] = S[i-1] - 0.05*beta*S[i-1]*I[i-1]/(S[i-1]+I[i-1])
I[i] = I[i-1] + 0.05*beta*S[i-1]*I[i-1]/(S[i-1]+I[i-1]) -
gamma*I[i-1]

A = Poisson('A', mu= S, value=susceptible_array, observed=True)
B = Poisson('B', mu= I, value=infected_array, observed=True)

[Anand Patil]

Hi Jason,

What went wrong with your attempt?

Anand 

[Chris Fonnesbeck]

On Oct 12, 6:45 pm, jason andrews <jasona...@gmail.com> wrote:
> Dear PyMCers,
>
> I wanted to see if people had more examples of dynamic compartmental
> models, such as difference equation models, monte carlo simulations
> using tau-leap, etc. to estimate transition parameters.
>
> I'm new to PyMC and was trying to do a simple SIR model to explore
> it.  SIR models (Susceptible-Infected-Recovered) are simple,
> compartmental infectious disease models.  They start with the
> population being distributed in one of the three states (SIR), and
> then susceptible individuals can move to the infected state at a rate
> (beta) times the product of the population of infected and susceptible
> individuals.  Individuals then recover at a rate gamma. There are a

> lot of great examples along with Python code here:http://www.warwick.ac.uk/~masfz/ModelingInfectiousDiseases/Chapter2/P...

Hi Jason,

I fit some SIR models using PyMC earlier this year. They follow the
methods described in Zipkin et al:

Zipkin, E.F., Jennelle, C.S. and Cooch, E.G. 2010. A primer on the
application of Markov chains to the study of wildlife disease
dynamics. Methods in Ecology and Evolution. 1: 192-198

There are several different versions, but I will post one of them
among the examples on the PyMC project page. Look for it today.

Hope that helps,
cf

[A. Flaxman]

It seems like many people had the same amount of delay in their responses.  I also fit compartmental models with PyMC, and I have often wondered if this is more common than I know of, and if not, why not?

 

The issue with your code is you need to make the compartments a deterministic variable.  I’ve cut-and-pasted your code into a github gist so that it has nice formatting:  http://gist.github.com/626689

 

--Abie

--
You received this message because you are subscribed to the Google Groups "PyMC" group.
To post to this group, send email to py...@googlegroups.com.
To unsubscribe from this group, send email to pymc+uns...@googlegroups.com.
For more options, visit this group at http://groups.google.com/group/pymc?hl=en.

[jason andrews]

Dear All,

Thanks for your quick responses and helpful comments.

Anand: When I ran the model with that code, the resultant parameters
(beta and gamma) were just giving me the distribution that I specified
rather than converging to fit the model. I didn't realize I had to
specify the compartments as deterministic parameters.

Chris: I'd love to see some more examples, that would be very
helpful. Thanks for the reference as well.

Abie: Many thanks for correcting my problem. This works perfectly. I
really appreciate your time and help.

thanks again,

Jason

On Oct 14, 2:21 pm, "A. Flaxman" <a...@uw.edu> wrote:
> It seems like many people had the same amount of delay in their responses.  I also fit compartmental models with PyMC, and I have often wondered if this is more common than I know of, and if not, why not?
>
> The issue with your code is you need to make the compartments a deterministic variable.  I've cut-and-pasted your code into a github gist so that it has nice formatting:  http://gist.github.com/626689
>
> --Abie
>
> From: py...@googlegroups.com [mailto:py...@googlegroups.com] On Behalf Of Anand Patil
> Sent: Thursday, October 14, 2010 1:54 AM
> To: py...@googlegroups.com
> Subject: Re: [pymc] dynamic compartmental models - ex. SIR models
>

> On Wed, Oct 13, 2010 at 12:45 AM, jason andrews <jasona...@gmail.com<mailto:jasona...@gmail.com>> wrote:
> Dear PyMCers,
>
> I wanted to see if people had more examples of dynamic compartmental
> models, such as difference equation models, monte carlo simulations
> using tau-leap, etc. to estimate transition parameters.
>
> I'm new to PyMC and was trying to do a simple SIR model to explore
> it.  SIR models (Susceptible-Infected-Recovered) are simple,
> compartmental infectious disease models.  They start with the
> population being distributed in one of the three states (SIR), and
> then susceptible individuals can move to the infected state at a rate
> (beta) times the product of the population of infected and susceptible
> individuals.  Individuals then recover at a rate gamma. There are a

> lot of great examples along with Python code here:http://www.warwick.ac.uk/~masfz/ModelingInfectiousDiseases/Chapter2/P...

> To post to this group, send email to py...@googlegroups.com<mailto:py...@googlegroups.com>.
> To unsubscribe from this group, send email to pymc+uns...@googlegroups.com<mailto:pymc+uns...@googlegroups.com> .