simulate_stem does not respect initial conditions from simulation_parameters

[damonbayer]

Working with slightly modified code from the SIR example.

library(stemr)
popsize = 1e4 # population size

true_pars =
  c(R0     = 1.5,  # basic reproduction number
    mu_inv = 2,    # infectious period duration = 2 days
    rho    = 0.5,  # case detection rate
    phi    = 10)   # negative binomial overdispersion

# initialize model compartments and rates
strata <- NULL # no strata
compartments <- c("S", "I", "R")

# rates initialized as a list of rate lists
rates <-
  list(rate(rate = "beta * I", # individual level rate (unlumped)
            from = "S",        # source compartment
            to   = "I",        # destination compartment
            incidence = T),    # compute incidence of S2I transitions, required for simulating incidence data
       rate(rate = "mu",       # individual level rate
            from = "I",        # source compartment
            to   = "R",        # destination compartment
            incidence = TRUE)) # compute incidence of I2R transitions (not required for simulating data)

# list used for simulation/inference for the initial state, initial counts fixed.
# state initializer a list of stem_initializer lists.
state_initializer <-
  list(stem_initializer(
    init_states = c(S = popsize-10, I = 10, R = 0), # must match compartment names
    fixed = F)) # initial state fixed for simulation, we'll change this later

# set the parameter values - must be a named vector
parameters =
  c(true_pars["R0"] / popsize / true_pars["mu_inv"], # R0 = beta * P / mu
    1/true_pars["mu_inv"],
    true_pars["rho"],
    true_pars["phi"])
names(parameters) <- c("beta", "mu", "rho", "phi")

# declare the initial time to be constant
constants <- c(t0 = 0)
t0 <- 0; tmax <- 40

# compile the model
dynamics <-
  stem_dynamics(
    rates = rates,
    tmax = tmax,
    parameters = parameters,
    state_initializer = state_initializer,
    compartments = compartments,
    constants = constants,
    compile_ode = T,   # compile ODE functions
    compile_rates = F, # compile MJP functions for Gillespie simulation
    compile_lna = F,   # compile LNA functions
    messages = F       # don't print messages
  )

# list of emission distribution lists (analogous to rate specification)
emissions <-
  list(emission(meas_var = "S2I", # transition or compartment being measured (S->I transitions)
                distribution    = "negbinomial",         # emission distribution
                emission_params = c("phi", "S2I * rho"), # distribution pars, here overdispersion and mean
                incidence       = TRUE,                  # is the data incidence
                obstimes        = seq(1, tmax, by =1)))  # vector of observation times

# compile the measurement process
measurement_process <-
  stem_measure(emissions = emissions,
               dynamics  = dynamics,
               messages  = F)

# put it all together into a stochastic epidemic model object
stem_object <-
  make_stem(
    dynamics = dynamics,
    measurement_process = measurement_process)

Without specifying simulation_parameters in simulate_stem, we randomly sample initial compartment counts:

set.seed(200)
a <- simulate_stem(stem_object = stem_object,
                   nsim = 1,
                   method = "ode", paths = F)


b <- simulate_stem(stem_object = stem_object,
                   nsim = 1,
                   method = "ode", paths = F)

tail(stem_object$dynamics$parameters, 3)
#>  S_0  I_0  R_0 
#> 9990   10    0
a$natural_paths[[1]][1,]
#>       time          S          I          R 
#>    0.00000 9989.73211   10.26789    0.00000
b$natural_paths[[1]][1,]
#>       time          S          I          R 
#>    0.00000 9992.58791    7.41209    0.00000

However, if we do specify simulation_parameters in simulate_stem, including some initial compartment sizes, these initial compartment sizes are not respected:

set.seed(200)

stem_object$dynamics$parameters
#>     beta       mu      rho      phi      S_0      I_0      R_0 
#> 7.50e-05 5.00e-01 5.00e-01 1.00e+01 9.99e+03 1.00e+01 0.00e+00

a <- simulate_stem(stem_object = stem_object,
                   nsim = 1, method = "ode", paths = F,
                   simulation_parameters = list(stem_object$dynamics$parameters))

b <- simulate_stem(stem_object = stem_object,
                   nsim = 1, method = "ode", paths = F,
                   simulation_parameters = list(stem_object$dynamics$parameters))

tail(stem_object$dynamics$parameters, 3)
#>  S_0  I_0  R_0 
#> 9990   10    0
a$natural_paths[[1]][1,]
#>       time          S          I          R 
#>    0.00000 9989.73211   10.26789    0.00000
b$natural_paths[[1]][1,]
#>       time          S          I          R 
#>    0.00000 9992.58791    7.41209    0.00000