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