Hello,
I am trying to get the output from the numerical simulation of a system
of ordinary differential equations for a range of values for three
parameters. I am using the ode() function (deSolve package) to run the
numerical simulation and apply() to run the simulation function for each
set of parameter values. I am having trouble getting the apply()
function to work.
Here is an example. Consider an epidemiology model of West Nile Virus:
wnv.model <- function(Time, State, Pars){
with(as.list(c(State, Pars)), {
dSB <- -(alphaB*betaB*SB*IM)/(SB + IB)
dIB <- (alphaB*betaB*SB*IM)/(SB + IB) - deltaB*IB
dSM <- bM*(SM + EM + IM) - (alphaM*betaB*IB*SM)/(SB
+ IB) - dM*SM
dEM <- (alphaM*betaB*SM*IB)/(SB + IB) - kappaM*EM -
dM*EM
dIM <- kappaM*EM - dM*IM
return(list(c(dSB, dIB, dSM, dEM, dIM)))
})
}
# I would like to run a numerical simulation of this model for each
combination of values for the parameters alphaB, betaB, # and kappaM:
av <- seq(0, 1, by = 0.2) # vector of values for alphaB
bv <- seq(0, 1, by = 0.2) # vector of values for betaB
kv <- seq(0, 1, by = 0.2) # vector of values for kappaM
# Here is my function with ode() for the numerical simulation. The
function returns the last value of the simulation for the # "IB" state
variable ("Infected Birds").
library(deSolve)
wnv.sim <- function(x){
Pars <- c(alphaB = x[1], betaB = x[2], deltaB = 0.2, bM =
0.03, dM = 0.03, alphaM = 0.69, kappaM = x[3])
State <- c(SB = 100, IB = 0, SM = 1500, EM = 0, IM = 0.0001)
Time <- seq(0, 60, by = 1)
model.out <- as.data.frame(ode(func = wnv.incubation,
y = State,
parms = Pars,
times = Time))
model.out[nrow(model.out),]$IB
}
# Finally, here is my apply() function:
dat <- apply(expand.grid(av, bv, kv), 1, wnv.sim)
However, the apply() function returns an error:
Error in eval(expr, envir, enclos) : object 'alphaB' not found
I am not sure what is wrong. Any help would be much appreciated.
Sincerely,
Adam
--
Adam Zeilinger
Post Doctoral Scholar
Department of Entomology
University of California Riverside
www.linkedin.com/in/adamzeilinger
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
