> 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=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
>
> ______________________________________________
> 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.
