Thanks to Dieter Menne and Spencer Graves I started to get my way through
lsoda()
Now I need to use it in with nls() to assess parameters
I have a go with a basic example
dy/dt = K1*conc
I try to assess the value of K1 from a simulated data set with a K1 close to
2.
Here is (I think) the best code that I've done so far even though it crashes
when I call nls()
--------------------------------------------------------------
x<-seq(0,10,,100)
y<-exp(2*x)
y<-rnorm(y,y,0.3*y)
test.model<-function(t,conc,parms){
dy.dt = parms["K1"]*conc
list(dy.dt)
}
require(deSolve)
foo<-lsoda(c(conc=1),times=seq(0,10,,100),test.model,parms=c(K1=2))
foo
#####################use of nls
func<-function(K1) {
foo<-lsoda(c(conc=1),times=seq(0,10,,100),test.model,parms=c(K1=K1))
foo[,"conc"]
}
nls(foo~func(K1),start=list(K1=1),data=data.frame(foo=y))
##################### have a look on the SSD ##################### y is the
vector of real data
SSD<-function(K1) {
sum((y-func(K1))^2)
}
data<-seq(1.5,2.1,,100)
plot(data,sapply(data,SSD),type="l")
--------------------------------------------------------------
Regards/Cordialement
Benoit Boulinguiez
-----Message d'origine-----
De : spencerg [mailto:spencer.graves at prodsyse.com]
Envoy? : vendredi 15 mai 2009 05:28
? : Benoit Boulinguiez
Cc : dieter.menne at menne-biomed.de; r-help at r-project.org
Objet : Re: [R] ode first step
Have you looked at the vignette in the deSolve package?
(deS <- vignette('compiledCode')) # opens a
"pdf" file
Stangle(deS$file) # writes an R script file to "getwd()"
In spite of the name, this vignette includes an example entirely in R.
By comparing it with your code, I see that you do NOT provide a connection
between y, parms, K1, C0, m, V, K2 and q. Something like the following
might work:
kinetic.model<-function(t,y,parms){
dq.dt = parms['K1']*y['C0'] -
(parms['K1']*y['m']/y['V']+
parms['K2'])*y['q']
list(dq.dt)
}
This may not be correct, but I hope the changes will help you see how
to make it work.
Bon Chance.
Spencer Graves
Benoit Boulinguiez wrote:> As I do not thoroughly understand the way 'lsoda' works, I face
some
> difficulties to 'get' myself into the function(), though I changed
the
code> as follows:
>
> ------------------------------
> require(deSolve)
>
> qm<-0.36
> y0<-c(0)
> parms<-c("K1","K2")
> times<-seq(0,10000,1)
> kinetic.model<-function(t,y,parms){
> dq.dt = K1*C0 - (K1*m/V+ K2)*q
> list(dq.dt)
> }
>
> foo<-lsoda(y0,times,kinetic.model,parms)
> Error in func(time, state, parms, ...) : object 'K1' not found
> ------------------------------
>
> 'K1' and 'K2' are parameters but 'C' is not a
parameter, it's a dependant
> variable of the time.
> I actually express it as a function of q(t) to get this new equation
> dq/dt= K1*C0 - (K1*m/V+ K2)*q(t)
> where K1 and K2 are the unknown but desired parameters and {C0,m,V} are
> constant known values.
>
> Nevertheless, I still get this 'Error about object 'K1' not
found'.
>
>
>
>
>
> Regards/Cordialement
>
>
> Benoit Boulinguiez
>
>
> -----Message d'origine-----
> De : Dieter Menne [mailto:dieter.menne at menne-biomed.de]
> Envoy? : jeudi 14 mai 2009 12:12
> ? : 'Benoit Boulinguiez'
> Objet : RE: [R] ode first step
>
> Try to hide yourself inside the function(). What would you see? No K1, for
> sure, no C, no K2.
> These are passed through parms, so parms["K1"] would work, but
not for C,
> you should add it.
>
> -----Original Message-----
> From: Benoit Boulinguiez [mailto:benoit.boulinguiez at ensc-rennes.fr]
> Sent: Thursday, May 14, 2009 11:53 AM
> To: 'Dieter Menne'
> Subject: RE: [R] ode first step
>
> ------------------------------
> qm<-0.36
> y0<-c(0)
> parms<-c(K1=1,K2=1)
> times<-seq(0,10000,1)
> kinetic.model<-function(t,y,parms){
> dq.dt<- K1*C*(qm-q)-K2*q
> list(dq.dt)
> }
>
> require(deSolve)
> nls(foo<-lsoda(y0,times,kinetic.model,parms)
>
> ______________________________________________
> R-help at 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.
>
>