R help - Apr 2012 - Using apply() with a function involving ode()

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

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