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@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@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@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.
