2. Multinomial in denim
2.1. Model definition
To specify the proportion that goes to each outgoing compartment, define the transition using the following syntax in model definition:
"proportion * compartment -> out_compartment" = [transition]
Note that the proportion here is the proportion of
compartment that will end up in
out_compartment at equilibrium.
Example: An SIRD model where 90% of infected can recover while the remaining 10% will die.
transitions <- list(
"S -> I" = "beta * S * (I / N) * timeStep",
"0.9 * I -> R" = d_gamma(1/3, 2),
"0.1 * I -> D" = d_exponential(0.1)
)Mathematical formulation for the model
\[ \begin{cases} dS = -\beta S \frac{I}{N} \\ dIR_1 = 0.9 *\beta S \frac{I}{N} -\frac{1}{3}IR_1 \\ dIR_2 = \frac{1}{3}IR_1 - \frac{1}{3}IR_2 \\ dID = 0.1*\beta S \frac{I}{N} - 0.1*ID \\ dR = \frac{1}{3}IR_2 \\ dD = 0.1*ID \end{cases} \]
Proportion normalization
denim will automatically normalize the proportions if
they don’t sum up to 1.
Example 2: The following model definition is equivalent to the one in Example 1
2.2. Example model
To further demonstrate the implementation of multinomial in
denim, we provide the equivalent model implemented in
deSolve and compare the output of 2 implementations.
Model definition in denim
Equivalent model definition in deSolve
# model in deSolve
transition_func <- function(t, state, param){
with(as.list( c(state, param) ), {
dS = - beta * S * (IR1 + IR2 + ID)/N
# apply linear chain trick for I -> R transition
# 0.9 * to specify prop of I that goes to I->R transition
dIR1 = 0.9 * beta * S * (IR1 + IR2 + ID)/N - rate*IR1
dIR2 = rate*IR1 - rate*IR2
dR = rate*IR2
# handle I -> D transition
# 0.1 * to specify prop of I that goes to I->D transition
dID = 0.1 * beta * S * (IR1 + IR2 + ID)/N - exp_rate*ID
dD = exp_rate*ID
list(c(dS, dIR1, dIR2, dID, dR, dD))
})
}
desolveInitialValues <- c(
S = 999,
# note that internally, denim also allocate initial value based on specified proportion
IR1 = 0.9,
IR2 = 0,
ID = 0.1,
R = 0,
D = 0
)Run simulation and compare
Model settings
Run model
# --- run denim model ----
mod <- sim(transitions = transitions,
initialValues = denimInitialValues,
parameters = parameters,
simulationDuration = simulationDuration,
timeStep = timeStep)
# run deSolve model
times <- seq(0, simulationDuration)
ode_mod <- ode(y = desolveInitialValues, times = times, parms = parameters, func = transition_func)
ode_mod <- as.data.frame(ode_mod)
ode_mod$I<- rowSums(ode_mod[, c("IR1", "IR2", "ID")])Output Comparison
