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

Question

I'm using the deSolve package to plot a couple 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 eventual goal is to create an iterative function or process (for loop) to plot how changes in certain parameters (beta and gamma) will affect the solution. The preferred output would be a list that contains all each ode solution for each specified value of beta in the loop. I'm running into issues for integrating a loop into the setup that the deSolve package requires for the ode function.

In the code below, I am trying to plot how the range of values (1 to 2 by increments of 0.1) in parameter beta will affect the plot of the 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'm getting states that the argument parameters in the SIR function is missing

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

I don't understand why this error is being reported when I've assigned parameters in the lines previous.

Answer 1

You might as well define your gradient function (and the 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 hold 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))
}

Now you have a list, each element of which contains the results from one run. You can sapply or lapply over this list, e.g. to get a matrix of the last states from each run:

t(sapply(res,tail,1))

Or if you want the results as one long data frame ...

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,...) would work nearly as well as bind_rows(), but bind_rows's .id argument is convenient for adding the beta values to each data frame. You could also leave the results as a list and loop over them while plotting with separate lines() calls, or (e.g.) bind just the infective columns together and use matplot() to plot them all at the same time. This is just a matter of style and idiom.

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