[R-sig-dyn-mod] nls and dynamic models

[Firstborn B]

Dear all,
I attempt to estimate parameters of a dynamic model using nls (and probably
later nlme). I have problems to use an ode model as a function to be used
in nls. This seems to concern both passing initial values as well as
parameters to the model.


library(deSolve)
seed=2423
dat2<-data.frame(days=(runif(20)+1)*10, X1=runif(20), X2=runif(20))
dat2$y<-0.4*exp(dat2$X1)+0.6*exp(dat2$X2)+rnorm(20, sd=0.3)
# example intentionally simple. I would usually solve it analysically
#***************************************************
#*Model definition
#***************************************************
decomp<-function(t, state, parameters){
  with(as.list(c(state, parameters)), {
    dX1<-a1*X1
    dX2<-a2*X2+a1*X1
    list(c(dX1, dX2))
  } )
}
# Testing the code to demonstrate that it works
parameters<-c(a1=0.05,a2=0.05)
state<-c(X1=0.5, X2=0.5)
times<-seq(0,100, by=1)
out<-ode(y=state, times=times, func=decomp, parms=parameters)
out[100,]
#*****************************
#* Wrapper function to be passed to nls or nlme
#**************************************
calcdecom<-function(a1,a2,t,x1, x2)
{
  state<-c(X1=x1, X2=x2)
  times<-c(0,t)
  parameters<-c(a1=a1,a2=a2)
  out<-ode(y=state, times=times, func=decomp, parms=parameters)
  return(as.numeric(out[2,2]+out[2,3]))
}
# ******************** test
calcdecom(0.1,0.1,5,0.3,0.3)
test<-nls(y~calcdecom(a1, a2, days, 0.5, 0.5 ),
          start=list(a1=0.02, a2=0.4), data=dat2)
test<-nls(y~calcdecom(a1, a2, days, X1, X2 ),
          start=list(a1=0.02, a2=0.4), data=dat2)


Error messages are for the first call to nls:

test<-nls(y~calcdecom(a1, a2, days, 0.5, 0.5 ),
+           start=list(a1=0.02, a2=0.4), data=dat2)
DINTDY-  T (=R1) illegal
In above message, R1 = 12.643

      T not in interval TCUR - HU (= R1) to TCUR (=R2)
In above message, R1 = 15.9793, R2 = 16.5415

DLSODA-  Trouble in DINTDY.  ITASK = I1, TOUT = R1
In above message, I1 = 1

And for the second call to nls
 Error in eval(substitute(expr), data, enclos = parent.frame()) :
  object 'X1' not found

	[[alternative HTML version deleted]]