non-linear optimisation ODE models

Previous Topic Next Topic
 
classic Classic list List threaded Threaded
5 messages Options
Reply | Threaded
Open this post in threaded view
|

non-linear optimisation ODE models

R help mailing list-2
Hello,
I'm new to R, so sorry for this question. I found a piece of code on stack overflow community, title: r-parameter and initial conditions fitting ODE models with nls.lm.
I've tried to implement a change suggested, but I get an error: Error in unname(myparms[4], B = 0, C = 0) :   unused arguments (B = 0, C = 0)
I'll appreciate any hint.
Malgosia

#set working directorysetwd("~/R/wkspace")#load librarieslibrary(ggplot2)library(reshape2)library(deSolve)library(minpack.lm)time=c(0,0.263,0.526,0.789,1.053,1.316,1.579,1.842,2.105,2.368,2.632,2.895,3.158,3.421,3.684,3.947,4.211,4.474,4.737,5)ca=c(0.957,0.557,0.342,0.224,0.123,0.079,0.035,0.029,0.025,0.017,-0.002,0.009,-0.023,0.006,0.016,0.014,-0.009,-0.03,0.004,-0.024)cb=c(-0.031,0.33,0.512,0.499,0.428,0.396,0.303,0.287,0.221,0.148,0.182,0.116,0.079,0.078,0.059,0.036,0.014,0.036,0.036,0.028)cc=c(-0.015,0.044,0.156,0.31,0.454,0.556,0.651,0.658,0.75,0.854,0.845,0.893,0.942,0.899,0.942,0.991,0.988,0.941,0.971,0.985)df<-data.frame(time,ca,cb,cc)dfnames(df)=c("time","ca","cb","cc")#plot datatmp=melt(df,id.vars=c("time"),variable.name="species",value.name="conc")ggplot(data=tmp,aes(x=time,y=conc,color=species))+geom_point(size=3)#rate functionrxnrate=function(t,c,parms){  #rate constant passed through a list called  k1=parms$k1  k2=parms$k2  k3=parms$k3  #c is the concentration of species  #derivatives dc/dt are computed below  r=rep(0,length(c))  r[1]=-k1*c["A"] #dcA/dt  r[2]=k1*c["A"]-k2*c["B"]+k3*c["C"] #dcB/dt  r[3]=k2*c["B"]-k3*c["C"] #dcC/dt  return(list(r))}# predicted concentration for a given parametercinit=c(A=1,B=0,C=0)t=df$timeparms=list(k1=2, k2=1, k3=3)out=ode(y=cinit,times=t,func=rxnrate,parms=list(k1=k1,k2=k2,k3=k3))head(out)
ssq=function(myparms){  #initial concentration  cinit=c(A=myparms[4],B=0,C=0)  cinit=c(A=unname(myparms[4],B=0,C=0))  print(cinit)  #time points for which conc is reported  #include the points where data is available  t=c(seq(0,5,0.1),df$time)  t=sort(unique(t))  #parameters from the parameters estimation  k1=myparms[1]  k2=myparms[2]  k3=myparms[3]  #solve ODE for a given set of parameters  out=ode(y=cinit,times=t,func=rxnrate,parms=list(k1=k1,k2=k2,k3=k3))  #Filter data that contains time points  outdf=data.frame(out)  outdf=outdf[outdf$time%in% df$time,]  #Evaluate predicted vs experimental residual  preddf=melt(outdf,id.var="time",variable.name="species",value.name="conc")  expdf=melt(df,id.var="time",variable.name="species",value.name="conc")  ssqres=preddf$conc-expdf$conc  return(ssqres)}# parameter fitting using levenberg marquart#initial guess for parametersmyparms=c(k1=0.5,k2=0.5,k3=3,1)cinit=c(A=unname(myparms[4],B=0,C=0))print(cinit)#fittingfitval=nls.lm(par=parms,fn=ssq)#summary of fitsummary(fitval)
        [[alternative HTML version deleted]]

______________________________________________
[hidden email] mailing list -- To UNSUBSCRIBE and more, see
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.
Reply | Threaded
Open this post in threaded view
|

Re: non-linear optimisation ODE models

Berend Hasselman

> On 15 Feb 2017, at 11:32, Malgorzata Wieteska via R-help <[hidden email]> wrote:
>
> Hello,
> I'm new to R, so sorry for this question. I found a piece of code on stack overflow community, title: r-parameter and initial conditions fitting ODE models with nls.lm.
> I've tried to implement a change suggested, but I get an error: Error in unname(myparms[4], B = 0, C = 0) :   unused arguments (B = 0, C = 0)
> I'll appreciate any hint.
> Malgosia
>
> #set working directorysetwd("~/R/wkspace")#load librarieslibrary(ggplot2)library(reshape2)library(deSolve)library(minpack.lm)time=c(0,0.263,0.526,0.789,1.053,1.316,1.579,1.842,2.105,2.368,2.632,2.895,3.158,3.421,3.684,3.947,4.211,4.474,4.737,5)ca=c(0.957,0.557,0.342,0.224,0.123,0.079,0.035,0.029,0.025,0.017,-0.002,0.009,-0.023,0.006,0.016,0.014,-0.009,-0.03,0.004,-0.024)cb=c(-0.031,0.33,0.512,0.499,0.428,0.396,0.303,0.287,0.221,0.148,0.182,0.116,0.079,0.078,0.059,0.036,0.014,0.036,0.036,0.028)cc=c(-0.015,0.044,0.156,0.31,0.454,0.556,0.651,0.658,0.75,0.854,0.845,0.893,0.942,0.899,0.942,0.991,0.988,0.941,0.971,0.985)df<-data.frame(time,ca,cb,cc)dfnames(df)=c("time","ca","cb","cc")#plot datatmp=melt(df,id.vars=c("time"),variable.name="species",value.name="conc")ggplot(data=tmp,aes(x=time,y=conc,color=species))+geom_point(size=3)#rate functionrxnrate=function(t,c,parms){  #rate constant passed through a list called  k1=parms$k1  k2=parms$k2  k3=parms$k3  #c is the concentratio!
 n of species  #derivatives dc/dt are computed below  r=rep(0,length(c))  r[1]=-k1*c["A"] #dcA/dt  r[2]=k1*c["A"]-k2*c["B"]+k3*c["C"] #dcB/dt  r[3]=k2*c["B"]-k3*c["C"] #dcC/dt  return(list(r))}# predicted concentration for a given parametercinit=c(A=1,B=0,C=0)t=df$timeparms=list(k1=2, k2=1, k3=3)out=ode(y=cinit,times=t,func=rxnrate,parms=list(k1=k1,k2=k2,k3=k3))head(out)
> ssq=function(myparms){  #initial concentration  cinit=c(A=myparms[4],B=0,C=0)  cinit=c(A=unname(myparms[4],B=0,C=0))  print(cinit)  #time points for which conc is reported  #include the points where data is available  t=c(seq(0,5,0.1),df$time)  t=sort(unique(t))  #parameters from the parameters estimation  k1=myparms[1]  k2=myparms[2]  k3=myparms[3]  #solve ODE for a given set of parameters  out=ode(y=cinit,times=t,func=rxnrate,parms=list(k1=k1,k2=k2,k3=k3))  #Filter data that contains time points  outdf=data.frame(out)  outdf=outdf[outdf$time%in% df$time,]  #Evaluate predicted vs experimental residual  preddf=melt(outdf,id.var="time",variable.name="species",value.name="conc")  expdf=melt(df,id.var="time",variable.name="species",value.name="conc")  ssqres=preddf$conc-expdf$conc  return(ssqres)}# parameter fitting using levenberg marquart#initial guess for parametersmyparms=c(k1=0.5,k2=0.5,k3=3,1)cinit=c(A=unname(myparms[4],B=0,C=0))print(cinit)#fittingfitval=nls.lm(par=par!
 ms,fn=ssq)#summary of fitsummary(fitval)


Totally unreadable code because you did not read the Posting Guide.

> [[alternative HTML version deleted]]


Do not post in HTML.

Berend Hasselman

>
> ______________________________________________
> [hidden email] mailing list -- To UNSUBSCRIBE and more, see
> 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.

______________________________________________
[hidden email] mailing list -- To UNSUBSCRIBE and more, see
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.
Reply | Threaded
Open this post in threaded view
|

Re: non-linear optimisation ODE models

Jim Lemon-4
Hi Malgorzata,
The function "rxnrate" seems to want three values in a list with the
names "k1", "k2" and "k3". If you are passing something with different
names, it is probably going to complain, so the names "A", "B" and "C"
may be your problem. I can't run the example, so this is a guess.

Jim


> On 15 Feb 2017, at 11:32, Malgorzata Wieteska via R-help <[hidden email]> wrote:
>
> Hello,
> I'm new to R, so sorry for this question. I found a piece of code on stack overflow community, title: r-parameter and initial conditions fitting ODE models with nls.lm.
> I've tried to implement a change suggested, but I get an error: Error in unname(myparms[4], B = 0, C = 0) :   unused arguments (B = 0, C = 0)
> I'll appreciate any hint.
> Malgosia
>
> #set working directorysetwd("~/R/wkspace")#load librarieslibrary(ggplot2)library(reshape2)library(deSolve)library(minpack.lm)time=c(0,0.263,0.526,0.789,1.053,1.316,1.579,1.842,2.105,2.368,2.632,2.895,3.158,3.421,3.684,3.947,4.211,4.474,4.737,5)ca=c(0.957,0.557,0.342,0.224,0.123,0.079,0.035,0.029,0.025,0.017,-0.002,0.009,-0.023,0.006,0.016,0.014,-0.009,-0.03,0.004,-0.024)cb=c(-0.031,0.33,0.512,0.499,0.428,0.396,0.303,0.287,0.221,0.148,0.182,0.116,0.079,0.078,0.059,0.036,0.014,0.036,0.036,0.028)cc=c(-0.015,0.044,0.156,0.31,0.454,0.556,0.651,0.658,0.75,0.854,0.845,0.893,0.942,0.899,0.942,0.991,0.988,0.941,0.971,0.985)df<-data.frame(time,ca,cb,cc)dfnames(df)=c("time","ca","cb","cc")#plot datatmp=melt(df,id.vars=c("time"),variable.name="species",value.name="conc")ggplot(data=tmp,aes(x=time,y=conc,color=species))+geom_point(size=3)#rate functionrxnrate=function(t,c,parms){  #rate constant passed through a list called  k1=parms$k1  k2=parms$k2  k3=parms$k3  #c is the concentratio!
>  n of species  #derivatives dc/dt are computed below  r=rep(0,length(c))  r[1]=-k1*c["A"] #dcA/dt  r[2]=k1*c["A"]-k2*c["B"]+k3*c["C"] #dcB/dt  r[3]=k2*c["B"]-k3*c["C"] #dcC/dt  return(list(r))}# predicted concentration for a given parametercinit=c(A=1,B=0,C=0)t=df$timeparms=list(k1=2, k2=1, k3=3)out=ode(y=cinit,times=t,func=rxnrate,parms=list(k1=k1,k2=k2,k3=k3))head(out)
>> ssq=function(myparms){  #initial concentration  cinit=c(A=myparms[4],B=0,C=0)  cinit=c(A=unname(myparms[4],B=0,C=0))  print(cinit)  #time points for which conc is reported  #include the points where data is available  t=c(seq(0,5,0.1),df$time)  t=sort(unique(t))  #parameters from the parameters estimation  k1=myparms[1]  k2=myparms[2]  k3=myparms[3]  #solve ODE for a given set of parameters  out=ode(y=cinit,times=t,func=rxnrate,parms=list(k1=k1,k2=k2,k3=k3))  #Filter data that contains time points  outdf=data.frame(out)  outdf=outdf[outdf$time%in% df$time,]  #Evaluate predicted vs experimental residual  preddf=melt(outdf,id.var="time",variable.name="species",value.name="conc")  expdf=melt(df,id.var="time",variable.name="species",value.name="conc")  ssqres=preddf$conc-expdf$conc  return(ssqres)}# parameter fitting using levenberg marquart#initial guess for parametersmyparms=c(k1=0.5,k2=0.5,k3=3,1)cinit=c(A=unname(myparms[4],B=0,C=0))print(cinit)#fittingfitval=nls.lm(par=pa!
 r!
>  ms,fn=ssq)#summary of fitsummary(fitval)
>
> [hidden email] mailing list -- To UNSUBSCRIBE and more, see
> 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.

______________________________________________
[hidden email] mailing list -- To UNSUBSCRIBE and more, see
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.
Reply | Threaded
Open this post in threaded view
|

Re: non-linear optimisation ODE models

David Winsemius

> On Feb 15, 2017, at 1:43 PM, Jim Lemon <[hidden email]> wrote:
>
> Hi Malgorzata,
> The function "rxnrate" seems to want three values in a list with the
> names "k1", "k2" and "k3". If you are passing something with different
> names, it is probably going to complain, so the names "A", "B" and "C"
> may be your problem. I can't run the example, so this is a guess.

There's a more readable version at:
http://stackoverflow.com/questions/42256509/how-to-feed-data-into-ide-while-doing-optimisation

It can be run, but does not produce the errors offered when I do so.

--
David.

>
> Jim
>
>
>> On 15 Feb 2017, at 11:32, Malgorzata Wieteska via R-help <[hidden email]> wrote:
>>
>> Hello,
>> I'm new to R, so sorry for this question. I found a piece of code on stack overflow community, title: r-parameter and initial conditions fitting ODE models with nls.lm.
>> I've tried to implement a change suggested, but I get an error: Error in unname(myparms[4], B = 0, C = 0) :   unused arguments (B = 0, C = 0)
>> I'll appreciate any hint.
>> Malgosia
>>
>> #set working directorysetwd("~/R/wkspace")#load librarieslibrary(ggplot2)library(reshape2)library(deSolve)library(minpack.lm)time=c(0,0.263,0.526,0.789,1.053,1.316,1.579,1.842,2.105,2.368,2.632,2.895,3.158,3.421,3.684,3.947,4.211,4.474,4.737,5)ca=c(0.957,0.557,0.342,0.224,0.123,0.079,0.035,0.029,0.025,0.017,-0.002,0.009,-0.023,0.006,0.016,0.014,-0.009,-0.03,0.004,-0.024)cb=c(-0.031,0.33,0.512,0.499,0.428,0.396,0.303,0.287,0.221,0.148,0.182,0.116,0.079,0.078,0.059,0.036,0.014,0.036,0.036,0.028)cc=c(-0.015,0.044,0.156,0.31,0.454,0.556,0.651,0.658,0.75,0.854,0.845,0.893,0.942,0.899,0.942,0.991,0.988,0.941,0.971,0.985)df<-data.frame(time,ca,cb,cc)dfnames(df)=c("time","ca","cb","cc")#plot datatmp=melt(df,id.vars=c("time"),variable.name="species",value.name="conc")ggplot(data=tmp,aes(x=time,y=conc,color=species))+geom_point(size=3)#rate functionrxnrate=function(t,c,parms){  #rate constant passed through a list called  k1=parms$k1  k2=parms$k2  k3=parms$k3  #c is the concentrati!
 o!
>> n of species  #derivatives dc/dt are computed below  r=rep(0,length(c))  r[1]=-k1*c["A"] #dcA/dt  r[2]=k1*c["A"]-k2*c["B"]+k3*c["C"] #dcB/dt  r[3]=k2*c["B"]-k3*c["C"] #dcC/dt  return(list(r))}# predicted concentration for a given parametercinit=c(A=1,B=0,C=0)t=df$timeparms=list(k1=2, k2=1, k3=3)out=ode(y=cinit,times=t,func=rxnrate,parms=list(k1=k1,k2=k2,k3=k3))head(out)
>>> ssq=function(myparms){  #initial concentration  cinit=c(A=myparms[4],B=0,C=0)  cinit=c(A=unname(myparms[4],B=0,C=0))  print(cinit)  #time points for which conc is reported  #include the points where data is available  t=c(seq(0,5,0.1),df$time)  t=sort(unique(t))  #parameters from the parameters estimation  k1=myparms[1]  k2=myparms[2]  k3=myparms[3]  #solve ODE for a given set of parameters  out=ode(y=cinit,times=t,func=rxnrate,parms=list(k1=k1,k2=k2,k3=k3))  #Filter data that contains time points  outdf=data.frame(out)  outdf=outdf[outdf$time%in% df$time,]  #Evaluate predicted vs experimental residual  preddf=melt(outdf,id.var="time",variable.name="species",value.name="conc")  expdf=melt(df,id.var="time",variable.name="species",value.name="conc")  ssqres=preddf$conc-expdf$conc  return(ssqres)}# parameter fitting using levenberg marquart#initial guess for parametersmyparms=c(k1=0.5,k2=0.5,k3=3,1)cinit=c(A=unname(myparms[4],B=0,C=0))print(cinit)#fittingfitval=nls.lm(par=p!
 a!

> r!
>> ms,fn=ssq)#summary of fitsummary(fitval)
>>
>> [hidden email] mailing list -- To UNSUBSCRIBE and more, see
>> 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.
>
> ______________________________________________
> [hidden email] mailing list -- To UNSUBSCRIBE and more, see
> 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.

David Winsemius
Alameda, CA, USA

______________________________________________
[hidden email] mailing list -- To UNSUBSCRIBE and more, see
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.
Reply | Threaded
Open this post in threaded view
|

Re: non-linear optimisation ODE models

Thomas Petzoldt-4
Hi,

fitting ODE models may also be done with package FME, see:

Soetaert K, Petzoldt T. Inverse modelling, sensitivity and Monte Carlo
analysis in R using package FME. Journal of Statistical Software.
2010(33): 1–28. http://dx.doi.org/10.18637/jss.v033.i03

or the (interactive) poster at:

http://desolve.r-forge.r-project.org/user2014/examples/FME/fit_twocomp.svg

Regards,

Thomas P.


Am 16.02.2017 um 19:41 schrieb David Winsemius:

>
>> On Feb 15, 2017, at 1:43 PM, Jim Lemon <[hidden email]> wrote:
>>
>> Hi Malgorzata,
>> The function "rxnrate" seems to want three values in a list with the
>> names "k1", "k2" and "k3". If you are passing something with different
>> names, it is probably going to complain, so the names "A", "B" and "C"
>> may be your problem. I can't run the example, so this is a guess.
>
> There's a more readable version at:
> http://stackoverflow.com/questions/42256509/how-to-feed-data-into-ide-while-doing-optimisation
>
> It can be run, but does not produce the errors offered when I do so.
>

______________________________________________
[hidden email] mailing list -- To UNSUBSCRIBE and more, see
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.