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

[Thomas Petzoldt]

On 4/5/2012 6:32 PM, mhimanshu wrote:
> 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)
> )

This data set is wrong because time and Z have 8 elements but Y has only 
7. Let's assume the missing Y is 6.0.

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

This means that the step adjustment algorithm was unable to determine an 
appropriate step size (h), possibly because of wrongly negative 
parameter values.

> Can you explain me why this is happening??  and Secondly,  I dont understand
> why i am getting parameters values in negatives after fitting??

You can restrict the parameter ranges by using the optional arguments 
"upper" and "lower", see ?modFit for details.

> Can you
> please help me out in this... :)

The reason for this can be manifold, e.g. wrong model specification, 
unrealistic data, inappropriate accuracy settings, inadequate start 
parameters or just unidentifiable parameters, e.g. if some of them 
depend strongly on each other.

The example below showed clearly that the parameters are strongly 
correlated (0.998 or even 1.000 !!!). This means that not all parameters 
can be identified simultaneously because of high interdependency (see my 
comment in my first post), so you should try to reduce the number of 
fitted parameters. Note however, that Ymax needs to be fitted (or 
adjusted) as well.

Fitting nonlinear models can be an art and requires mathematical 
understanding of the applied techniques (differential equations, 
identifiability analysis, numerical optimization), a good understanding 
of the properties of your (differential equation) model -- and some 
feeling and experience.


Thomas

##---------------------------------------------------------------------

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=0.8)
#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)


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

modelCost <- function(p) {
   out <- ode(y = y, func = derivs, parms = p, times = yobs$time)
   return(modCost(out, yobs, weight="mean"))
}

## start values for the parameters
startpars <- c(a1 = 5, a2 = 0.002, a3 = 0.002, a4 = 0.02, Ymax=7)
plower <- c(a1 = 1, a2 = 0.0001, a3 = 0.0001, a4 = 0.0001, Ymax=0.001)
pupper <- c(a1 = 10, a2 = 2, a3 = 1, a4 = 1, Ymax=10)


## fit the model; nprint = 1 shows intermediate results
fit <- modFit(f = modelCost, p = startpars,
   upper=pupper, lower=plower, control = list(nprint = 1))
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)

summary(fit)

#Parameters:
#      Estimate Std. Error t value Pr(>|t|)
#a1   5.610e+00  2.913e+03   0.002    0.998
#a2   2.000e+00  2.908e+03   0.001    0.999
#a3   1.872e-03  3.162e-03   0.592    0.566
#a4   1.188e-04  1.224e-01   0.001    0.999
#Ymax 8.908e+00  4.615e+03   0.002    0.998
#
#Residual standard error: 0.08001 on 11 degrees of freedom
#
#Parameter correlation:
#           a1       a2       a3      a4     Ymax
#a1    1.00000  1.00000 -0.09916 -0.1031  1.00000
#a2    1.00000  1.00000 -0.09917 -0.1032  1.00000
#a3   -0.09916 -0.09917  1.00000  0.9981 -0.09917
#a4   -0.10315 -0.10315  0.99807  1.0000 -0.10316
#Ymax  1.00000  1.00000 -0.09917 -0.1032  1.00000