R help - Mar 2012 - How to find best parameter values using deSolve n optim() ?

Hello,
Can someone please help me out with the optim() function.

# parameters 
pars<- c(a1= 0.9, a2= 0.7, a3= 0.06, a4=0.02)

y<- c(Y=0.2, Z=0.1)
Ymax<- c(0.8)
#### fucntion deriv
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)))
 })
 }

 times <- c(seq(0, 10, 0.1))
 out <- ode(y = y, parms = pars, times = times, func = derivs)
 as.data.frame(out)

I have a set of 2 differential equations, and i want to find the best
parameter values for a1 & a2. 
I have an experimental data for dy so that I have to find the best values
for a1 & a2 according to the experimental values so that the observed values
for dy should come near to experimental dy values.

can someone please suggest me how to proceed further or may give their
opinion.

Thanks alot in advance..!!
Himanshu

--
View this message in context:
http://r.789695.n4.nabble.com/How-to-find-best-parameter-values-using-deSolve-n-optim-tp4506042p4506042.html
Sent from the R help mailing list archive at Nabble.com.

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)

Hi Thomas,

Thank you so much for your suggestion. 
I tried your code and it is working fine. Now when I change the values of Y
in yobs I am getting so many warnings.

say, 
yobs <- data.frame(
 time = 0:7,
 Y = c(0.00, 3.40, 4.60 ,5.80, 5.80, 6.00, 6.00 ),
 Z = c(0.1, 0.11, 0.119, 0.128, 0.136, 0.145, 0.153, 0.16)
)

So when i fit the model with the same code that you have written, i got the
following warnings:
DLSODA-  Warning..Internal T (=R1) and H (=R2) are                              
      such that in the machine, T + H = T on the next step                      
     (H = step size). Solver will continue anyway.                              
      In above,  R1 =  0.1484502806322D+01   R2 =  0.2264549048113D-16

and I have got so many such warnings.

Can you explain me why this is happening??  and Secondly,  I dont understand
why i am getting parameters values in negatives after fitting?? Can you
please help me out in this... :)

Thanks
Himanshu

--
View this message in context:
http://r.789695.n4.nabble.com/How-to-find-best-parameter-values-using-deSolve-n-optim-tp4506042p4535368.html
Sent from the R help mailing list archive at Nabble.com.

On 6/6/2012 3:50 PM, mhimanshu wrote:
> Hello Thomas,
>
> This code seems to be fine and its now working well.
>
> I read the about the FME package, but I have one doubt, as in the data set
> given in the paper, it showing a nice kinetics of the viral growth, so my
> question is what if there is a sudden increase in viral growth after some
> interval, say Bimodal growth curve?
>
> How does it fits the bimodal growth curve? I tried with FME but I am not
> getting the desired results. May be you can explain me a little, I would be
> really grateful to you. :)
>
> Thanks a lot,
> Himanshu
Hi and thanks for the feedback,

regarding your problem with a "bimodal growth curve" I am not
completely
sure what you mean. However, I suspect that failing to fit a bimodal 
behaviour may not be caused by parameter fitting with FME, but instead 
would need extension of the underlying ODE model.

Thomas