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.
You might as well define your gradient function (and the other non-changing elements) outside the loop:
Now define the vector of values to try (not necessary but convenient):
and define a list to hold the results:
Now you have a list, each element of which contains the results from one run. You can
sapply
orlapply
over this list, e.g. to get a matrix of the last states from each run:Or if you want the results as one long data frame ...
do.call(rbind,...)
would work nearly as well asbind_rows()
, butbind_rows
's .id
argument is convenient for adding thebeta
values to each data frame. You could also leave the results as a list and loop over them while plotting with separatelines()
calls, or (e.g.) bind just the infective columns together and usematplot()
to plot them all at the same time. This is just a matter of style and idiom.