在 ODE 中使用循环以图形方式比较不同参数 R

Question

我正在使用 deSolve 包来绘制几个微分方程（如果有兴趣请阅读 http://www.maa.org/press/periodicals/loci/joma/the-sir-model-for-spread-of-disease-the-differential-equation-model）。

我的最终目标是创建一个迭代函数或过程（for 循环）来绘制某些参数（beta 和 gamma）的变化将如何影响解决方案。首选输出将是一个列表，其中包含循环中每个指定 beta 值的所有每个 ode 解决方案。我运行遇到了将循环集成到 deSolve 包需要 ode 函数的设置中的问题。

在下面的代码中，我试图绘制参数 beta 中的值范围（1 到 2，增量为 0.1）将如何影响微分方程的绘图。

for(k in seq(1,2,by=0.1)){ #range of values for beta

    init <- c(S=1-1e-6, I=1e-6, R=0) #initial conditions for odes
    time <- seq(0,80,by=1) #time period
    parameters <- c(beta=k, gamma=0.15) #parameters in ode

SIR <- function(time,state,parameters){ #function containing equaations
    with(as.list(c(state,parameters)),{
        dS <- -beta*S*I
        dI <- beta*S*I-gamma*I
        dR <- gamma*I

        return(list(c(dS,dI,dR)))
    })
}

ode(y=init,times=time,func=SIR()[beta],parms=parameters[k])}

}

我遇到的第一个错误指出 SIR 函数中的自变量参数丢失

Error in as.list(c(init, parameters)) : argument "parameters" is missing, with no default

我不明白为什么当我在前面的行中分配 parameters 时会报告此错误。

Answer 1

你不妨在循环外定义你的渐变函数（和其他不变的元素）：

SIR <- function(time,state,parameters) {
  with(as.list(c(state,parameters)),{
    dS <- -beta*S*I
    dI <- beta*S*I-gamma*I
    dR <- gamma*I
    return(list(c(dS,dI,dR)))
  })
}
init <- c(S=1-1e-6, I=1e-6, R=0) #initial conditions for odes
time <- seq(0,80,by=1) #time period

现在定义要尝试的值向量（不是必需但方便）：

betavec <- seq(1,2,by=0.1)

并定义一个列表来保存结果：

res <- vector(length(betavec),mode="list")

library(deSolve)
for (k in seq_along(betavec)){ #range of values for beta
    res[[k]] <- ode(y=init,times=time,func=SIR,
           parms=c(beta=betavec[k], gamma=0.15))
}

现在您有了一个列表，其中的每个元素都包含一个运行的结果。您可以 sapply 或 lapply 覆盖此列表，例如从每个运行:

中获取最后一个状态的矩阵

t(sapply(res,tail,1))

或者，如果您希望将结果作为一个长数据框...

names(res) <- betavec  ## to get beta value incorporated in results
dd <- dplyr::bind_rows(lapply(res,as.data.frame),.id="beta")
dd$beta <- as.numeric(dd$beta)

do.call(rbind,...) 的效果几乎与 bind_rows() 一样好，但是 bind_rows 的 .id 参数便于将 beta 值添加到每个数据框。您还可以将结果保留为列表并在使用单独的 lines() 调用绘制时循环遍历它们，或者（例如）仅将感染列绑定在一起并使用 matplot() 同时绘制它们。这只是风格和习语的问题。

library(ggplot2); theme_set(theme_bw())
library(viridis)
ggplot(dd,aes(x=time,y=I,colour=beta))+
    geom_line(aes(group=beta))+
    scale_color_viridis()+
    scale_y_log10()

在 ODE 中使用循环以图形方式比较不同参数 R

Using a loop in an ODE to graphically compare different parameters R

r

ode