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)
non-linear optimisation ODE models
5 messages · Malgorzata Wieteska, Berend Hasselman, Jim Lemon +2 more
On 15 Feb 2017, at 11:32, Malgorzata Wieteska via R-help <r-help at r-project.org> 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 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)
Totally unreadable code because you did not read the Posting Guide.
[[alternative HTML version deleted]]
Do not post in HTML. Berend Hasselman
______________________________________________ R-help at r-project.org 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.
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 <r-help at r-project.org> 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)
R-help at r-project.org 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.
On Feb 15, 2017, at 1:43 PM, Jim Lemon <drjimlemon at gmail.com> 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 <r-help at r-project.org> 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)
>>
>> R-help at r-project.org 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.
>
> ______________________________________________
> R-help at r-project.org 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
9 days later
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 <drjimlemon at gmail.com> 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.