[R] Using apply() with a function involving ode()

[Adam Zeilinger]

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