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

[Berend Hasselman]

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