The effect of the distribution of transition rate on dynamics

[Fuhan-Yang]

Hi Dr.King,

Recently, I have strange results fitting an SIRS model: my data showed sustained dynamics, while the estimated Reff is below 1. After sanity check using ode(method = 'euler') from deSolve package, I found the euler-multinomial transition may cause this issue.

Using the code at the end, I built a deterministic SIRS model and generated three trajectories of Infected with ode(), pomp() with binomial transition, and euler-multinomial transition. The output is the same between ode() and the pomp() using binomial transition. While the pomp with multinomial transition showed epidemic even the R0 = 0.195. This is strange to me. The differences between binomial rate and euler-multinomial rate are small, but is it possible that the small differences change the dynamic a lot and cause epidemic?

Can you please comment on this? Or please let me know if there is any error in the code (see below). I appreciate your help. Thanks!

library(pomp)
library(dplyr)

### prolong the dynamic time to 100 years ###
inc_df <- data.frame(
  daynum = c(1:36520),
  ili = rpois(36520, 60)
)

sirs_init <- Csnippet("
  S = N-300;
  I1 = 300;
  I2 = 0;
  R1 = R2 = R3 = 0;
  J = 0;
")

### process model: deterministic SIRS ###
sirs_multinom_step <- Csnippet("
  double pR = 1-exp(-3/pi*dt);
  double pI = 1-exp(-2*mu_IR*dt);
  double dN_SI = S*(1-exp(-beta*(I1+I2)/N*dt));
  double dN_I1 = I1*pI;
  double dN_I2 = I2*pI;
  double dN_R1 = R1*pR;
  double dN_R2 = R2*pR;
  double dN_R3 = R3*pR;
  S  += dN_R3 - dN_SI;
  I1 += dN_SI - dN_I1;
  I2 += dN_I1 - dN_I2;
  R1 += dN_I2 - dN_R1;
  R2 += dN_R1 - dN_R2;
  R3 += dN_R2 - dN_R3;
  J  += dN_SI;
")

sirs_binom_step <- Csnippet("
  double pR = 3/pi*dt;
  double pI = 2*mu_IR*dt;
  double dN_SI = S*(beta*(I1+I2)/N*dt);
  double dN_I1 = I1*pI;
  double dN_I2 = I2*pI;
  double dN_R1 = R1*pR;
  double dN_R2 = R2*pR;
  double dN_R3 = R3*pR;
  S  += dN_R3 - dN_SI;
  I1 += dN_SI - dN_I1;
  I2 += dN_I1 - dN_I2;
  R1 += dN_I2 - dN_R1;
  R2 += dN_R1 - dN_R2;
  R3 += dN_R2 - dN_R3;
  J  += dN_SI;
")

### measure model ###
inc_dmeas <- Csnippet(r"{
  if (ISNA(ili))
    lik = (give_log) ? 0 : 1;
  else
    lik = dpois(ili,rho*J,give_log);
}")

inc_rmeas <-  Csnippet(r"{
  ili = rpois(rho*J);}")

## create pomp ##

## to build pomp with binomial rate, simply replace the rprocess component

sirs_det <- pomp(
  data = inc_df,
  times = "daynum",
  t0 = 1,
  rinit = sirs_init,
  rprocess = pomp::euler(sirs_multinom_step,delta.t = 1),
  rmeasure = inc_rmeas,
  dmeasure = inc_dmeas,
  obsnames = "ili",
  statenames = c("S","I1","I2","R1","R2","R3","J"),
  paramnames = c("beta","mu_IR","pi","N","rho"),
  accumvars = "J"
)

## simulation ##
params <- c(beta = 0.19, rho = 0.108, pi = 100, mu_IR = 0.2,N = 1e7)

sirs_det |>
  simulate(
    nsim = 1,
    params = params,
    format = "data.frame"
  ) |>
  mutate(
    I = I1 + I2,
    R = R1 + R2 + R3
  ) -> po_sims

##### deSolve version ####
library(deSolve)

sirs_2I_3R <- function(t, y, params) {

  S <- y[1]
  I1 <- y[2]
  I2 <- y[3]
  R1 <- y[4]
  R2 <- y[5]
  R3 <- y[6]
  J <- y[7]

  beta <- params["beta"]
  gamma <- params["mu_IR"]
  pi <- params["pi"]
  N <- params["N"]

  dS <- -beta*S*(I1+I2)/N + 3/pi*R3
  dI1 <- beta*S*(I1+I2)/N - 2*gamma*I1
  dI2 <- 2*gamma*I1 - 2*gamma*I2
  dR1 <- 2*gamma*I2 - 3/pi*R1
  dR2 <- 3/pi*R1 - 3/pi*R2
  dR3 <- 3/pi*R2 - 3/pi*R3

  ## cumulative infected
  dJ <- beta*S*(I1 + I2)/N

  result <- c(dS, dI1, dI2, dR1, dR2, dR3, dJ)

  list(result)
}

sirs_2I_3R_solve <- function(params, I_init = 300, burn_in = 3650, total = 7300) {

  params <- params
  N <- params["N"]
  start <- c(N-I_init,I_init,0,0,0,0,0)
  names(start) <- c("S","I1","I2","R1","R2","R3","J")

  time <- seq(1,total,1)

  out <- ode(
    y = start,
    times = time,
    func = sirs_2I_3R,
    parms = params,
    method = rkMethod("euler")
  )
  out <- data.frame(out)
  out$inc <- c(0,diff(out$J))

  ## Exclude burn-in period
  out[burn_in:total,]
}

out <- sirs_2I_3R_solve(params, burn_in = 1, total = 36520)

out |>
  mutate(
    I = I1 + I2,
    R = R1 + R2 + R3
  ) -> out

op <- par(mfrow = c(2,1),mar = c(2,3,2,2))
plot(out$I,type = "l",main = "deSolve",ylab = "")
plot(po_sims$I,type = "l",main = "Euler-multinomial",ylab = "")
par(op)

[kingaa]

You seem to suppose that the two computations will yield similar results. On what basis do you suppose this?

The Euler method applied to the differential equation you implement in your "deSolve" version is not equivalent to Euler steps applied to the deterministic map you supply in your "Euler-multinomial" version, is it?

However, the two should become closer with smaller step sizes. For example, in the following, I have replaced the timestep of 1 day with a timestep of 1/100 day.

library(pomp)
library(dplyr)

sirs_init <- Csnippet("
  S = N-300;
  I1 = 300;
  I2 = 0;
  R1 = R2 = R3 = 0;
  J = 0;
")

sirs_det_step <- Csnippet("
  double pR = 1-exp(-3/pi*dt);
  double pI = 1-exp(-2*mu_IR*dt);
  double dN_SI = S*(1-exp(-beta*(I1+I2)/N*dt));
  double dN_I1 = I1*pI;
  double dN_I2 = I2*pI;
  double dN_R1 = R1*pR;
  double dN_R2 = R2*pR;
  double dN_R3 = R3*pR;
  S  += dN_R3 - dN_SI;
  I1 += dN_SI - dN_I1;
  I2 += dN_I1 - dN_I2;
  R1 += dN_I2 - dN_R1;
  R2 += dN_R1 - dN_R2;
  R3 += dN_R2 - dN_R3;
  J  += dN_SI;
")

simulate(
  times=seq(1,36520,by=1),
  t0 = 1,
  rinit = sirs_init,
  rprocess = pomp::euler(sirs_det_step,delta.t = 0.01),
  statenames = c("S","I1","I2","R1","R2","R3","J"),
  paramnames = c("beta","mu_IR","pi","N","rho"),
  accumvars = "J",
  params=c(beta=0.19, rho=0.108, pi=100, mu_IR=0.2, N=1e7),
  format="data.frame"
) |>
  mutate(
    I = I1 + I2,
    R = R1 + R2 + R3
  ) -> po_out

##### deSolve version ####
library(deSolve)

sirs_2I_3R <- function(t, y, params) {

  S <- y[1]
  I1 <- y[2]
  I2 <- y[3]
  R1 <- y[4]
  R2 <- y[5]
  R3 <- y[6]
  J <- y[7]

  beta <- params["beta"]
  gamma <- params["mu_IR"]
  pi <- params["pi"]
  N <- params["N"]

  dS <- -beta*S*(I1+I2)/N + 3/pi*R3
  dI1 <- beta*S*(I1+I2)/N - 2*gamma*I1
  dI2 <- 2*gamma*I1 - 2*gamma*I2
  dR1 <- 2*gamma*I2 - 3/pi*R1
  dR2 <- 3/pi*R1 - 3/pi*R2
  dR3 <- 3/pi*R2 - 3/pi*R3
  ## cumulative infected
  dJ <- beta*S*(I1 + I2)/N

  list(c(dS, dI1, dI2, dR1, dR2, dR3, dJ))
}

sirs_2I_3R_solve <- function (params, I_init = 300, burn_in = 3650, total = 7300) {
  N <- params["N"]
  start <- c(N-I_init,I_init,0,0,0,0,0)
  names(start) <- c("S","I1","I2","R1","R2","R3","J")
  time <- seq(1,total,1)
  out <- ode(
    y = start,
    times = time,
    func = sirs_2I_3R,
    parms = params,
    method = rkMethod("euler")
  )
  out <- data.frame(out)
  out$inc <- c(0,diff(out$J))
  out[burn_in:total,]
}

sirs_2I_3R_solve(
  params = c(beta=0.19, rho=0.108, pi=100, mu_IR=0.2, N=1e7),
  burn_in = 1,
  total = 36520
) |>
  mutate(
    I = I1 + I2,
    R = R1 + R2 + R3
  ) -> ds_out

op <- par(mfrow = c(3,1))
plot(ds_out$I,type = "l",main = "deSolve Euler",ylab = "")
plot(po_out$I,type = "l",main = "pomp Euler",ylab = "")
plot(ds_out$I,po_out$I,type = "p",xlab = "deSolve Euler",ylab="pomp Euler")
abline(a=0,b=1)
par(op)

[Fuhan-Yang]

Thanks for the clarification. I mixed R0 when using different rates. R0 is no longer beta/mu_IR when using euler-multinomial transitions. It should be (1-exp(-beta*dt))/(1-exp(-mu_IR*dt)).
Using the same beta and mu_IR, euler-multinomial transitions yield higher R0 than binomial transitions, while the differences get small when dt is small. Still, we may need to interpret the beta and R0 in a different way when we use euler-multinomial transitions. Thanks a lot!

[kingaa]

I'm not sure what you mean. I didn't see you define anything in terms of Euler-multinomial or binomial distributions in the foregoing code, despite the fact that you use those terms. It appears that your understanding of those terms does not agree with their usage in the package or the literature. For this reason, I don't what the conclusion you draw actually means.

In general, though, you are correct that one's definition of R0 depends strongly on the model structure. It's also worth reflecting on the fact that the meaning of a parameter is determined by its role in the dynamical system. Model parameters don't have a meaning outside of the context of the model. From this point of view, the meaning of "beta" is different in the various formulations you have....

[Fuhan-Yang]

Sorry for being unclear. By default, I calculated R0 from estimated parameters based on my ode in deSolve, which is beta/mu_IR. While this equation is not valid when using the euler-multinomial distribution in pomp. Failing to take account of this leads to misinterpretation of R0, as shown in my question, the epidemic still happens when the "R0" is below 1. When switching to binomial distribution in pomp, the results are consistent with deSolve. That's a relief!