Skip to contents

Migrate to denim

Source: vignettes/deSolve_to_denim.Rmd
deSolve_to_denim.Rmd
## Warning: package 'deSolve' was built under R version 4.3.1

Migrate deSolve code to denim

Original code in deSolve

The model used for demonstrating the process of migrating code from deSolve to denim is as followed

# --- Model definition in deSolve
transition_func <- function(t, state, param){
  with(as.list( c(state, param) ), {
      gamma_rate = 1/scale
      dS = -rate*S
      dI1 = rate*S - gamma_rate*I1
      dI2 = gamma_rate*I1 - gamma_rate*I2
      dI =  dI1 + dI2
      dR = gamma_rate*I2
      list(c(dS, dI, dI1, dI2, dR))
  })
}

# ---- Model configuration 
parameters <- c(rate = 0.2, scale = 3, beta = 0.12) 
initialValues <- c(S = 999, I = 0, I1 = 0, I2=0, R=0)

# ---- Run simulation
times <- seq(0, 20) # simulation duration
ode_mod <- ode(y = initialValues, times = times, parms = parameters, func = transition_func)

# --- show output
ode_mod <- as.data.frame(ode_mod)
head(ode_mod[ode_mod$time %in% 1:20, c("time", "S", "I", "R")])
##   time        S        I          R
## 2    1 817.9120 178.1031   2.984842
## 3    2 669.6497 310.0129  19.337352
## 4    3 548.2628 397.6768  53.060395
## 5    4 448.8796 447.4499 102.670468
## 6    5 367.5116 467.1136 164.374867
## 7    6 300.8930 464.2888 233.818204

Model definition

Unlike deSolve where transitions between compartments are defined by a system of ODEs, transitions in denim must be defined explicitly using the built-in distributions.

User must first identify the distribution type that best describe the transition in their deSolve model.

Identify distribution
# --- Model definition in deSolve
transition_func <- function(t, state, param){
  with(as.list( c(state, param) ), {
      gamma_rate = 1/scale
      
      # For S -> I transition, constant rate is assumed --> exponentially distributed dwell time
      dS = -rate*S
      # For I -> R transition, linear chain trick is applied --> gamma distributed dwell time 
      dI1 = rate*S - gamma_rate*I1
      dI2 = gamma_rate*I1 - gamma_rate*I2
      dI =  dI1 + dI2
      dR = gamma_rate*I2
      list(c(dS, dE, dI, dI1, dI2, dR))
  })
}

With the transitions identified, user can then define the model in denim.

When using denim, the model structure is given as a list of key-value pairs where

  • key is a string showing the transition direction between compartments

  • value is the built-in distribution function that describe the transition

# --- Transition def for denim
transitions <- list(
  "S -> I" = d_exponential(0.2),
  "I -> R" = d_gamma(3, 2) # shape is 2 from number of I sub compartments
)

Model configurations

Similar to deSolve, denim also ask users to provide the initial values and any additional parameters in the form of named vectors.

For the example deSolve code, while users can use the initalValues as is (denim will ignore unused I1, I2 compartments), it is recommended to remove redundant compartments.

For parameters, since rate and scale are already defined in the distribution functions, users only need to keep beta from the initial parameters vector.

# remove I1, I2 compartments
denim_initialValues <- c(S = 999, I = 0, R=0)
denim_parameters <- c(beta = 0.12)

Simulation

Lastly, users need to define the simulation duration and time step for denim to run. Unlike deSolve which takes a time sequence, denim only require the simulation duration and time step.

Since denim is a discrete time model, time step must be set to a small value for the result to closely follow that of deSolve.

mod <- sim(transitions = transitions,
             initialValues = denim_initialValues, 
             parameters = denim_parameters,
             simulationDuration = 20,
             timeStep = 0.01)
plot(mod)

Compare output