Dear Berend,
Yes, the wnv.sim function should have included "wnv.model" instead of
"wnv.incubation." Sorry for the typo.
Although I had inspected the expand.grid() object, I had not connected
the error with a naming problem. Thanks so much for the help! The
function works perfectly.
Sincerely,
Adam
On 4/25/2012 2:02 AM, Berend Hasselman wrote:
See inline comments.
On 25-04-2012, at 08:40, Adam Zeilinger wrote:
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,
There is no function wnv.incubation.
Assuming you mean wnv.model
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
Have you inspected the object returned by expand.grid?
Do head(expand.grid(av,bv,kv))
and you will see that the columns have names which confuse the rest of your
code.
I was able to repair by doing
pars.grid<- expand.grid(av, bv, kv)
names(pars.grid)<- NULL
dat<- apply(pars.grid, 1, wnv.sim)
HTH,
Berend
--
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.
