Re: [R] How to find best parameter values using deSolve n optim() ?
Dear Thomas, Thank you for the quick reply. Might be the Bimodal/Biphasic behavior is due to the ODE model. Now I am working on "modMCMC" function of FME. This is strange this function sometimes works fine but other times I am getting following errors and warnings. Error in modCost(out, yobs, weight = "std") : error: cannot estimate weighing for observed variable: Y In addition: Warning messages: 1: In lsoda(y, times, func, parms, ...) : an excessive amount of work (> maxsteps ) was done, but integration was not successful - increase maxsteps 2: In lsoda(y, times, func, parms, ...) : Returning early. Results are accurate, as far as they go I have read through the manual also but couldn't figure out the error. The code works fine if I fix some of the parameter (pars) but it gives above errors if i pass all the parameters(pars). Can you explain why this is happening? -- View this message in context: http://r.789695.n4.nabble.com/How-to-find-best-parameter-values-using-deSolve-n-optim-tp4506042p4633100.html Sent from the R help mailing list archive at Nabble.com. __ 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.
Re: [R] How to find best parameter values using deSolve n optim() ?
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 __ 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.
Re: [R] How to find best parameter values using deSolve n optim() ?
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 -- View this message in context: http://r.789695.n4.nabble.com/How-to-find-best-parameter-values-using-deSolve-n-optim-tp4506042p4632528.html Sent from the R help mailing list archive at Nabble.com. __ 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.
Re: [R] How to find best parameter values using deSolve n optim() ?
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.0020.998
#a2 2.000e+00 2.908e+03 0.0010.999
#a3 1.872e-03 3.162e-03 0.5920.566
#a4 1.188e-04 1.224e-01 0.0010.999
#Ymax 8.908e+00 4.615e+03 0.0020.998
#
#Residual standard error: 0.08001 on 11 degrees of freedom
#
#Parameter correlation:
# a1 a2 a3 a4 Ymax
#a11.0 1.0 -0.09916 -0.1031 1.0
#a21.0 1.0 -0.09917 -0.1032 1.0
#a3 -0.09916 -0.09917 1.0 0.9981 -0.09917
#a4 -0.10315 -0.10315 0.99807 1. -0.10316
#Ymax 1.0 1.0 -0.09917 -0.1032 1.0
__
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.
Re: [R] How to find best parameter values using deSolve n optim() ?
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. __ 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.
Re: [R] How to find best parameter values using deSolve n optim() ?
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)
__
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.
[R] 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.
__
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.
