[R] How to find best parameter values using deSolve n optim() ?

[Thomas Petzoldt]

Hi Himanshu,

the use of optim is described on its help page.

In addition to this package FME provides additional functionality for 
fitting ODE models. See FME help files, package vignettes and the 
example below for details.

Hope it helps

Thomas Petzoldt

################################################################

library(deSolve)
library(FME)

## function derivs
derivs <- function(time, y, pars) {
   with (as.list(c(y, pars)), {
     dy = a1 * Y * (1 - Y/Ymax) - a2 * (Y + 0.001)
     dz = a3 * Y - a4 * Z;
     return(list(c(dy, dz)))
   })
}

## parameters
pars <- c(a1 = 0.9, a2 = 0.7, a3 = 0.06, a4 = 0.02)
Ymax <- c(0.8)
## initial values
y <- c(Y = 0.2, Z = 0.1)
## time steps
times <- c(seq(0, 10, 0.1))
## numerical solution of ODE
out <- ode(y = y, parms = pars, times = times, func = derivs)

## example observation data
yobs <- data.frame(
   time = 0:7,
   Y = c(0.2, 0.195, 0.19, 0.187, 0.185, 0.183, 0.182, 0.18),
   Z = c(0.1, 0.11, 0.119, 0.128, 0.136, 0.145, 0.153, 0.16)
)

## model cost function, see help file and vignette for details
modelCost <- function(p) {
   out <- ode(y = y, func = derivs, parms = p, times = yobs$time)
   return(modCost(out, yobs))
}

## start values for the parameters
startpars <- c(a1 = 1, a2 = 0.6, a3 = 0.1, a4 = 0.1)

## fit the model; nprint = 1 shows intermediate results
fit <- modFit(f = modelCost, p = startpars, control = list(nprint = 1))

## Note the high correlation between a1 and a2, resp a3 and a4
## that can make parameter identification difficult
summary(fit)

## graphical result
out2 <- ode(y = y, parms = startpars, times = times, func = derivs)
out3 <- ode(y = y, parms = fit$par, times = times, func = derivs)
plot(out, out2, out3, obs = yobs)

legend("topleft", legend=c("original", "startpars", "fitted"),
   col = 1:3, lty = 1:3)