如何在R中通过多个参数集循环微分方程？

Question

我正在 R 中对传染病的传播进行建模，我需要通过多组参数循环模型。到目前为止，我有一个工作模型（下面是一个虚拟版本），但仅适用于一组变量。此外，我可以通过一个参数的不同值循环模型，但我不知道如何通过参数值的多个向量循环它。非常感谢任何帮助！

我还需要代码来保存每个场景的模型输出，因为我将使用cowplot和ggplot创建一个组合图来比较场景之间的S-I-R动态。

以下是示例场景：

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

这是 SIR 模型：

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

重复但具有多个参数空间

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

我需要代码在每个场景中运行模型并保存模型输出，这样我就可以使用cowplot和ggplot创建一个比较S-I-R动态的组合图。我不确定如何从这里继续，尽管我认为 for 循环是必要的。

Answer 1

可以在循环中完成此操作，但大多数 R 编码员会使用

pmap

和

unnest

。 更优雅！

最简单的方法是编写一个接受所有模拟参数（初始条件和速率常数）的函数，调用 ode 求解器，并返回带有数字列的数据框。

然后定义一个数据框，每个条件一行。

pmap

通过条件进行“循环”的繁重工作。

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