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.gra...@prodsyse.com] Envoyé : vendredi 15 mai 2009 05:28 À : Benoit Boulinguiez Cc : dieter.me...@menne-biomed.de; r-help@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.me...@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.boulingu...@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@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-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.