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 ______________________________________________ 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.