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.

Arthur · Accepted Answer · 2025-01-22 20:39:28Z

3

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))
)

answered Jan 22 at 20:39

Arthur

2,49210 silver badges25 bronze badges

Sign up to request clarification or add additional context in comments.

3 Comments

Ben Bolker Jan 22 at 20:50

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 Jan 23 at 5:45

"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 Feb 12 at 18:00

@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'".

jay.sf · Accepted Answer · 2025-01-23 11:42:57Z

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

Collectives™ on Stack Overflow

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

2 Answers 2

3 Comments

Comments

Your Answer

Linked

Hot Network Questions

Collectives™ on Stack Overflow

2 Answers 2

3 Comments

Comments

Your Answer

Sign up or log in

Post as a guest

Linked

Related