[R] non-linear optimisation ODE models

[Malgorzata Wieteska]

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)
	[[alternative HTML version deleted]]