[R] ode first step

[Benoit Boulinguiez]

Thanks to Dieter Menne and Spencer Graves I started to get my way through
lsoda()
Now I need to use it in the nls() to assess the parameter.

I've tried 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.

I'm not sure that I'm using nls() and lsoda() correctly.
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.
>
>