# SIR modeling in Soss

**URL:** <https://discourse.julialang.org/t/sir-modeling-in-soss/36680>\
**Category:** Probabilistic Programming\
**Created:** [March 29, 2020, 3:57pm UTC](https://discourse.julialang.org/t/sir-modeling-in-soss/36680 "2020-03-29T15:57:28Z")\
**Posts on this page:** 17\
**Page:** 1

<div class="post-metadata">

**Author:** ![cscherrer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cscherrer/32/7631_2.png) [@cscherrer](https://discourse.julialang.org/u/cscherrer)\
**Post date:** [March 29, 2020, 3:57pm UTC](https://discourse.julialang.org/t/sir-modeling-in-soss/36680/1 "2020-03-29T15:57:28Z")

</div>

I spent some time yesterday getting as SIR epidemiology model working in Soss. I see this a lot as a differential equation, but that doesn’t make sense for real data, because what we really observe are counts. So I used a discrete Markov chain. The transitions are “global” relative to the Markov chain, so this was an interesting design challenge.

EDIT: I should say, it doesn’t make sense _to me_ for real data _with small counts_.

I ended up with this for the Markov kernel:

```julia
using Soss

mstep = @model pars,state begin
    # Parameters
    α = pars.α # Transmission rate
    β = pars.β # Recovery rate
    γ = pars.γ # Case fatality rate

    # Starting counts
    s0 = state.s # Susceptible
    i0 = state.i # Infected 
    r0 = state.r # Recovered
    d0 = state.d # Deceased
    n0 = state.n # Population size
    
    # Transitions between states
    si ~ Binomial(s0, α * i0 / n)
    ir ~ Binomial(i0, β)
    id ~ Binomial(i0-ir, γ)
    
    # Updated counts
    s = s0 - si 
    i = i0 + si - ir - id
    r = r0 + ir
    d = d0 + id
    n = n0 - id

    next = (pars=pars, state=(s=s,i=i,r=r,d=d, n=n))
end;

```

This is wrapped in another model that puts priors on the rate parameters:

```julia
m = @model s0 begin
    α ~ Uniform()
    β ~ Uniform()
    γ ~ Uniform()
    pars = (α=α, β=β, γ=γ)
    x ~ MarkovChain(pars, mstep(pars=pars, state=s0))
end

```

_EDIT: Benjamin Vincent caught some errors in the model and its interpretation, see [here](https://discourse.julialang.org/t/sir-modeling-in-soss/36680/16)._

And here’s inference:

```julia
julia> s0 = namedtuple(df[1,:])
(i = 1, r = 0, d = 0, n = 331000000, s = 330999999, si = 0, ir = 0, id = 0)
julia> post = dynamicHMC(m(s0=s0),(x=namedtuple.(eachrow(df)),));
julia> ppost=particles(post)
(γ = 0.00423 ± 0.00011, β = 0.00234 ± 7.9e-5, α = 0.279 ± 0.00084)

```

The next obvious thing is to use this for a forecast. I don’t want to share this yet, as the result is really bleak. I haven’t yet done any model checking, but **I expect the model is much too simple for realistic forecasting**.

To draw from the prior, you can do

```julia
s0 = (i = 1, r = 0, d = 0, n = 331000000, s = 330999999);

r = rand(m(s0=s0));

for (n,s) in enumerate(r.x)
    n>10 && break
    println(s)
end

```

* * *

Note that I had to add a few things to get this to work. [Full code is available here](https://gist.github.com/cscherrer/d649f7d4477bde4f0843f2492442cf20)

_EDIT: I’ve updated the model to work with the current `master` branch of Soss; `MarkovChain` is now included._

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [March 29, 2020, 4:09pm UTC](https://discourse.julialang.org/t/sir-modeling-in-soss/36680/2 "2020-03-29T16:09:21Z")

</div>

> [@cscherrer](#):
>
> I see this a lot as a differential equation, but that doesn’t make sense for real data, because what we really observe are counts.

Everything is an approximation, and discrete time is also an approximation. Dropping that gives a continuous time discrete model, also known as a Gillespie model, is a form of “discrete differential equation” which limits to stochastic differential equation and an ordinary differential equation. This model is fairly widely used in epidemic modeling, and coincidentally the DiffEq documentation has been using an epidemic model in its tutorial for this:

[https://docs.sciml.ai/latest/tutorials/discrete\_stochastic\_example/#Discrete-Stochastic-(Gillespie)-Equations-1](https://docs.sciml.ai/latest/tutorials/discrete_stochastic_example/#Discrete-Stochastic-(Gillespie)-Equations-1)

```julia
using DiffEqBiological
sir_model = @reaction_network SIR begin
    c1, s + i --> 2i
    c2, i --> r
end c1 c2

p = (0.1/1000,0.01)
prob = DiscreteProblem([999,1,0],(0.0,250.0),p)
jump_prob = JumpProblem(prob,Direct(),sir_model)
sol = solve(jump_prob,FunctionMap())
using Plots; plot(sol)

```

![image](https://global.discourse-cdn.com/julialang/original/3X/2/6/268f80b875e3b58e4d3de13e6bc805288e27e884.png)

As the copy number (the number of individuals) approaches infinity, the stochasticity in the model drops (since that grows as sqrt(N) as opposed to the mean of N), and thus it approaches the ODE behavior (which should be clear from the behavior). That’s why the ODE is widely used, because then you can do things like [bifurcation analysis](https://rveltz.github.io/PseudoArcLengthContinuation.jl/dev/) to predict stability of the model and do things like derive an approximating equation for the peak of the distribution.

But there are many ways to discretize the continuous time discrete equations, and it turns out your model is equivalent to the binomial leaping method on this kind of model but with a fixed dt: [Binomial leap methods for simulating stochastic chemical kinetics - PubMed](https://www.ncbi.nlm.nih.gov/pubmed/15549913) . Note that there’s a lot of studies on how you can then improve the model by making it implicit, adapting dt, etc. to stabilize it. So the other thing is just that differential equations (and discrete continuous-time Markov chains) are just a good place to do numerical analysis to figure out how things can go unstable and how to derive more accurate/stable methods.

That model should directly work in Turing.jl with DiffEqBayes.jl, but it would be nice to get an example together with SOSS!

---

<div class="post-metadata">

**Author:** ![DoktorMike](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/doktormike/32/2736_2.png) [@DoktorMike](https://discourse.julialang.org/u/DoktorMike)\
**Post date:** [March 29, 2020, 5:37pm UTC](https://discourse.julialang.org/t/sir-modeling-in-soss/36680/3 "2020-03-29T17:37:37Z")

</div>

I wrote a blog post about the SEIR model applied to the covid19 situation in Denmark. As a next step I would indeed like to have the rates be distributions instead to get a reasonable quantification of the uncertainty involved. Do I need DiffEqBayes for using Turing or Soss to accomplish that? I’ve never tampered with stochasticity in ODEs before. The blog post is here: [mike@home - COVID-19 in Denmark - An epidemic in a small country](https://doktormike.gitlab.io/post/covid-19/) and the code is written in Julia.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [March 29, 2020, 5:41pm UTC](https://discourse.julialang.org/t/sir-modeling-in-soss/36680/4 "2020-03-29T17:41:38Z")

</div>

You don’t need DiffEqBayes to use it with Turing. If you just use `concrete_solve` in there even the Zygote-based autodiff will just work out fine. The DiffEqBayes is just a simplified interface for folks who are not so familiar with probabilistic programming and just want to churn out the standard Bayesian parameter fit, and it’s setup inside to always try to “do the right thing” in that case. (It’s also a very good way to continue testing that DiffEq and Turing work together).

---

<div class="post-metadata">

**Author:** ![DoktorMike](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/doktormike/32/2736_2.png) [@DoktorMike](https://discourse.julialang.org/u/DoktorMike)\
**Post date:** [March 29, 2020, 6:05pm UTC](https://discourse.julialang.org/t/sir-modeling-in-soss/36680/5 "2020-03-29T18:05:41Z")

</div>

Got it.

---

<div class="post-metadata">

**Author:** ![Paul\_Schrimpf](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/paul_schrimpf/32/10134_2.png) [@Paul\_Schrimpf](https://discourse.julialang.org/u/Paul_Schrimpf)\
**Post date:** [March 30, 2020, 2:10am UTC](https://discourse.julialang.org/t/sir-modeling-in-soss/36680/6 "2020-03-30T02:10:20Z")

</div>

That’s a nice blog post. I computed some Bayesian estimates of a similar model.

[https://schrimpf.github.io/CovidSEIR.jl/covid/](https://schrimpf.github.io/CovidSEIR.jl/covid/)

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [March 30, 2020, 5:37am UTC](https://discourse.julialang.org/t/sir-modeling-in-soss/36680/7 "2020-03-30T05:37:24Z")

</div>

Thanks for a very interesting analysis 🙂 I have a question if you don’t mind.

- For [italy](https://schrimpf.github.io/CovidSEIR.jl/italy/), one of the chains (orange) seems very different from the others, does this correspond to a different mode of the posterior? If so, could there be other modes that have not been explored?

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [March 30, 2020, 7:54am UTC](https://discourse.julialang.org/t/sir-modeling-in-soss/36680/8 "2020-03-30T07:54:09Z")

</div>

@Paul_Schrimpf, it appears that you have very bad mixing, this should be fixed before any kind of conclusion is drawn from the model (kudos for including the relevant plots to diagnose this, I wish all researchers did this).

---

<div class="post-metadata">

**Author:** ![Paul\_Schrimpf](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/paul_schrimpf/32/10134_2.png) [@Paul\_Schrimpf](https://discourse.julialang.org/u/Paul_Schrimpf)\
**Post date:** [March 30, 2020, 6:39pm UTC](https://discourse.julialang.org/t/sir-modeling-in-soss/36680/9 "2020-03-30T18:39:18Z")

</div>

Reply, yes there could be other modes. In general, the chains do not look great. I’m running them for more iterations now.

---

<div class="post-metadata">

**Author:** ![Paul\_Schrimpf](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/paul_schrimpf/32/10134_2.png) [@Paul\_Schrimpf](https://discourse.julialang.org/u/Paul_Schrimpf)\
**Post date:** [March 30, 2020, 6:40pm UTC](https://discourse.julialang.org/t/sir-modeling-in-soss/36680/10 "2020-03-30T18:40:43Z")

</div>

Agreed. Thanks for the feedback. I am working on improving it.

---

<div class="post-metadata">

**Author:** ![Paul\_Schrimpf](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/paul_schrimpf/32/10134_2.png) [@Paul\_Schrimpf](https://discourse.julialang.org/u/Paul_Schrimpf)\
**Post date:** [March 31, 2020, 1:55am UTC](https://discourse.julialang.org/t/sir-modeling-in-soss/36680/11 "2020-03-31T01:55:00Z")

</div>

@baggepinnen, @Tamas_Papp After running for many more iterations, the chains still appear different from one another for Italy, Canada, and the US. @baggepinnen, I think you’re right that the posterior is simply multi-modal. In some sense this isn’t surprising — the model has many unobserved state variables, and seeing a short part of the path of the three observed variables is not enough to separately identify the parameters.

For China and Korea, where the number of active cases has begun decreasing, the chains look better and the posterior appears uni-modal. On the other hand, the model does not fit the data as well in these countries. This too is not very surprising. A model with constant transmission rates cannot explain the evolution of cases seen in those countries.

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [March 31, 2020, 7:44am UTC](https://discourse.julialang.org/t/sir-modeling-in-soss/36680/12 "2020-03-31T07:44:40Z")

</div>

> [@Paul\_Schrimpf](#):
>
> the model has many unobserved state variables, and seeing a short part of the path of the three observed variables is not enough to separately identify the parameters.

I wonder if this is the right intuition though. Usually non-identified parameters just result in “ridges” and their higher dimension equivalents, not separate modes. Model misspecification is usually a more common reason for these kind of mixing problems.

Whatever the reason may be, my point was that unless the chains mix well, the inference will not be valid.

Maybe a model with time-varying transmission rates could be interesting. Eg consider a multilevel model where the transmission rate comes from common lognormal (with parameters estimated). Could this pick up the policy changes?

---

<div class="post-metadata">

**Author:** ![pbayer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pbayer/32/11675_2.png) [@pbayer](https://discourse.julialang.org/u/pbayer)\
**Post date:** [March 31, 2020, 8:23am UTC](https://discourse.julialang.org/t/sir-modeling-in-soss/36680/13 "2020-03-31T08:23:38Z")

</div>

You may be interested to look at the [SIR-X model](https://www.medrxiv.org/content/10.1101/2020.02.18.20024414v1). See also the [COVID-19 event horizon](http://rocs.hu-berlin.de/corona/docs/forecast/results_by_country/) of Humboldt university Berlin.

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [March 31, 2020, 8:31am UTC](https://discourse.julialang.org/t/sir-modeling-in-soss/36680/14 "2020-03-31T08:31:18Z")

</div>

It is indeed an interesting model (thanks for mentioning it), but assumes continuous parameter/policy changes. It is my impression that most countries implemented large, discrete changes.

In any case what I am most skeptical about is the measurement side of the models. Even assuming a sophisticated model for the transmission, symptoms, and recovery/death, only serious cases show up in the data with certainty. Testing for the rest of the population is very heterogeneous by country and time, and probably not modeled well with an IID approach (eg testing susceptible contacts).

---

<div class="post-metadata">

**Author:** ![drbenvincent](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/drbenvincent/32/3046_2.png) [@drbenvincent](https://discourse.julialang.org/u/drbenvincent)\
**Post date:** [April 5, 2020, 10:03am UTC](https://discourse.julialang.org/t/sir-modeling-in-soss/36680/15 "2020-04-05T10:03:11Z")

</div>

Very minor point, but I believe γ is not the case fatality rate as that is the fatality rate per time unit. So perhaps the case fatality rate is γ/β ?

---

<div class="post-metadata">

**Author:** ![cscherrer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cscherrer/32/7631_2.png) [@cscherrer](https://discourse.julialang.org/u/cscherrer)\
**Post date:** [April 5, 2020, 3:11pm UTC](https://discourse.julialang.org/t/sir-modeling-in-soss/36680/16 "2020-04-05T15:11:13Z")

</div>

I agree γ is really the fatality rate per time unit, but I don’t yet how the case “eventual fatality” rate is γ/β. I’d think that’s an odds ratio, and the actual rate is γ/(β+γ). But I might well be missing something.

Also, thanks for [your note about changing population size](https://twitter.com/inferencelab/status/1246677342078095361)!

---

<div class="post-metadata">

**Author:** ![drbenvincent](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/drbenvincent/32/3046_2.png) [@drbenvincent](https://discourse.julialang.org/u/drbenvincent)\
**Post date:** [April 5, 2020, 3:31pm UTC](https://discourse.julialang.org/t/sir-modeling-in-soss/36680/17 "2020-04-05T15:31:10Z")

</div>

Yes, I think you’re right.