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)