"Simulation-based Inference", Lesson 2

[erinclancey]

Hi Dr. King,

I have a question about your example from Lesson 2 in the notes you have online from your short course on Simulation-based Inference for Epidemiological Dynamics. In the example below that I copied from you notes you have two transitions, S to I and I to R. Beta is a parameter and it is multiplied by I, but not S. mu_IR is also a parameter but it is not multiplied by I or anything else. Is this because the size in rbinom(size,prob) serves as the S in Beta*S*I or the I in mu_IR*I? I just want to make sure I am coding my model correctly. Thank you very much!

sir_step <- Csnippet("
double dN_SI = rbinom(S,1-exp(-Beta*I/N*dt));
double dN_IR = rbinom(I,1-exp(-mu_IR*dt));
S -= dN_SI;
I += dN_SI - dN_IR;
R += dN_IR;
H += dN_IR;
")

measSIR |>
  pomp(
    rprocess=euler(sir_step,delta.t=1/7),
    rinit=sir_rinit,
    rmeasure=sir_rmeas,
    dmeasure=sir_dmeas,
    accumvars="H",
    statenames=c("S","I","R","H"),
    paramnames=c("Beta","mu_IR","N","eta","rho","k")
  ) -> measSIR

[kingaa]

In a word, yes. You have it right. A call to rbinom(n,p), returns one random draw from a Binomial(n,p) distribution. The expectation of this distribution is n p and it has variance n p (1-p). In the line,

double dN_SI = rbinom(S,1-exp(-Beta*I/N*dt));

S is the number of individuals at risk (of infection) and Beta*I/N is the hazard (risk per unit time). In this context, it is called the force of infection. In the next line,

double dN_IR = rbinom(I,1-exp(-mu_IR*dt));

I is the number at risk (of recovery) and mu_IR is the hazard.
The definitive reference for calls to R's distribution functions is in the R extensions manual. The related eulermultinomial distribution is described in the pomp documentation.