Using a loop in ODE to graphically compare various R parameters

Question

I use the package deSolveto build a pair of differential equations (read if interested http://www.maa.org/press/periodicals/loci/joma/the-sir-model-for-spread-of-disease-the-differential-equation -model ).

My ultimate goal is to create an iterative function or process (for a loop) to determine how changes in certain parameters (beta and gamma) will affect the decision. The preferred output would be a list containing all of each solution odefor each given beta value in the loop. I am having problems integrating the loop into the configuration required by the package deSolvefor the function ode.

In the code below, I'm trying to tell how the range of values ​​(from 1 to 2 in increments of 0.1) in the beta parameter will affect the graph of differential equations.

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

}

The first error I get indicates that there are no argument parameters in the SIR function

Error in as.list (c (init, parameters)): the argument "parameters" is missing, without a default value

I do not understand why this error is reported when I assigned parametersin the previous lines.

Answer 1

You can also define your gradient function (and other non-changing elements) outside the loop:

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

Now define the vector of values ​​to try (not necessary, but convenient):

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

and define a list to store the results:

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