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