How to loop differential equations through multiple parameter sets in R?

Question

I'm modeling the spread of an infectious disease in R, and I need to loop the model through multiple sets of parameters. I have a working model so far (a dummy version is below), but only for a single set of variables. Additionally, I can loop the model through a different values for one parameter, but I don't know how to loop it through multiple vectors of parameter values. Any help is greatly appreciated!

I also need the code to save the model outputs for each scenario, since I will be using cowplot and ggplot to create a combined figure comparing the S-I-R dynamics between scenarios.

Here is the SIR model:

library(deSolve)
library(ggplot2)
library(dplyr)
library(tidyverse) 
parms = c("beta" = 0.00016, "gamma" = 0.12)
CoVode = function(t, x, parms) {
 S = x[1] # Susceptible
 I = x[2] # Infected
 R = x[3] # Recovered
 beta = parms["beta"]
 gamma = parms["gamma"]
 dSdt <- -beta*S*I
 dIdt <- beta*S*I-gamma*I
 dRdt <- gamma*I
 output = c(dSdt,dIdt,dRdt)
 names(output) = c('S', 'I', 'R')
 return(list(output))
} 
# Initial conditions
init = numeric(3)
init[1] = 10000
init[2] = 1
names(init) = c('S','I','R')
# Run the model
odeSim = ode(y = init, 0:50, CoVode, parms)
# Plot results
Simdat <- data.frame(odeSim)
Simdat_long <- Simdat %>%
 pivot_longer(!time, names_to = "groups", values_to = "count")
ggplot(Simdat_long, aes(x= time, y = count, group = groups, colour = groups)) +
 geom_line() 

Repeat but with multiple parameter spaces:

parmdf <- data.frame(scenario=c("A", "B", "C", "D", "E", "F", "G", "H", "I"),
 beta=c(0.0001,0.0001, 0.0001, 0.001, 0.001, 0.001, 0.124, 0.124, 0.124), 
 gamma=c(0.1, 0.2, 0.3, 0.1, 0.2, 0.3, 0.1, 0.2, 0.3)
)

I need the code to run the model through each scenario and save the model outputs so I can then use cowplot and ggplot to create a combined figure comparing the S-I-R dynamics. I am unsure how to proceed from here, although I assume a for-loop will be necessary.

Answer 1

It's possible to accomplish this in a loop, but most R coders would instead use a pmap and an unnest. More elegant!

It's easiest to do this by writing a function that accepts all simulation parameters (both the initial conditions and the rate constants), calls the ode solver, and returns a data frame with numeric columns.

And then defining a data frame with one row per condition.

pmap does the heavy lifting of "looping" through the conditions.

library(tidyverse) 
library(deSolve)
theme_set(theme_bw())
CoVode <- function(t, x, parms) {
 
 S = x[1] # Susceptible
 I = x[2] # Infected
 R = x[3] # Recovered
 
 beta = parms["beta"]
 gamma = parms["gamma"]
 
 dSdt <- -beta*S*I
 dIdt <- beta*S*I-gamma*I
 dRdt <- gamma*I
 
 output = c(S = dSdt, I = dIdt, R = dRdt)
 return(list(output))
 
} 
sim <- function(beta, gamma, S0, I0, R0) {
 # simulate at specified conditions
 p <- c(beta = beta, gamma = gamma)
 y0 <- c(S = S0, I = I0, R = R0)
 ode(y0, 0:50, CoVode, p) |> 
 as_tibble() |> 
 mutate(across(everything(), as.numeric))
}
# check
sim(0.00016, 0.12, 10000, 1, 0)
# grid of conditions
disease <- expand_grid(
 beta = c(0.0001, .0002, .0003),
 gamma = c(0.1, 0.2),
 S0 = 10000, I0 = 1, R0 = 0
)
disease$solution <- disease |> 
 pmap(sim)
# solution is a list; will want to unnest
disease <- disease |> 
 unnest(solution)
# plot
print(
 disease |> 
 pivot_longer(c(S, I, R), names_to = "groups", values_to = "counts") |> 
 ggplot() + 
 aes(time, counts, color = groups) +
 geom_line() + 
 facet_wrap(~paste("beta =", beta, "gamma =", gamma))
)

enter image description here

• Ben Bolker: This is a fine answer, but there's also a lot to be said for the transparency of for loops: see Chapter 4 of the R Inferno.
• Roland: "most R coders would instead use a pmap and an unnest" That's a bold statement. I don't know if it's true. How do you know?
• Kit McLean: @BenBolker I tried to build the follwing for loop: for (i in parmdf) { parms = as.numeric(unlist(parmdf[i,-1])) names(parmdf) = names(parms) odeSim = ode(y = init, 0:50, CoVode, parms) } but I get the error message "Error in as.numeric(parmdf[i, -1]) : 'list' object cannot be coerced to type 'double'".

Answer 2

You could asplit your parmdf or use any method giving you named vectors in a list for each scenario. Then put it in Map.

res <- Map(\(x) ode(y=init, 0:50, CoVode, parms=x), asplit(parmdf[-1], 1) |> 
 setNames(parmdf$scenario))

Gives:

> lapply(res, head, 3)
$A
 time S I R
[1,] 0 10000.000 1.000000 0.0000000
[2,] 1 9998.378 2.459440 0.1621746
[3,] 2 9994.392 6.047242 0.5609795
$B
 time S I R
[1,] 0 1.000000e+04 1.000 0.0000
[2,] 1 3.553107e+03 6240.940 206.9528
[3,] 2 7.295285e-01 8095.132 1905.1388
$C
 time S I R
[1,] 0 1.000000e+04 1.000 0.000
[2,] 1 1.115403e-10 7425.454 2575.546
[3,] 2 -3.325534e-11 5500.910 4500.090